PyroShaper/solver.py at master · arda-guler/PyroShaper · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
import math
import tkinter as tk
import time
import matplotlib.pyplot as plt
import numpy as np

from mesh import *

def generate_tubular_grain(msh):
    for c in msh.cells:
        if (msh.Nx/2)**2 >= (c.x - msh.Nx/2)**2 + (c.y - msh.Ny/2)**2 >= (msh.Nx*0.3)**2:
            c.fuel = 100

def ignite_tubular_grain(msh):
    for c in msh.cells:
        if c.fuel > 0:
            c.active = True
            if c.get_burn_rate(msh) > 0 and (msh.Nx*0.45)**2 >= (c.x - msh.Nx/2)**2 + (c.y - msh.Ny/2)**2:
                c.active = True
            else:
                c.active = False

def generate_cocentric_grain(msh):
    for c in msh.cells:
        if (msh.Nx*0.5)**2 > (c.x - msh.Nx/2)**2 + (c.y - msh.Ny/2)**2:
            c.fuel = 100
            if (msh.Nx*0.3)**2 >= (c.x - msh.Nx/2)**2 + (c.y - msh.Ny/2)**2 >= (msh.Nx*0.2)**2:
                c.fuel = 0

def ignite_tubular_grain(msh):
    for c in msh.cells:
        if c.fuel > 0:
            c.active = True
            if c.get_burn_rate(msh) > 0 and (msh.Nx*0.45)**2 >= (c.x - msh.Nx/2)**2 + (c.y - msh.Ny/2)**2:
                c.active = True
            else:
                c.active = False

def low_pass_filter(y, cutoff):
    n = len(y)
    fft_y = np.fft.fft(y)
    fft_y[int(cutoff):int(n-cutoff)] = 0
    filtered_y = np.fft.ifft(fft_y)

    return np.real(filtered_y)

def simulate(msh, dt):
    root = tk.Tk()
    root.title("Solid ")
    cvas = tk.Canvas(root, width=900, height=590, bg="black")
    cvas.grid(row=0, column=0, rowspan=15, columnspan=5)

    thrusts = []
    times = []
    fuels = []

    running = True
    time = 0
    cycle = 0
    while running:
        # PHYSICS
        for c in msh.cells:
            c.update(msh, dt)

        fuel = 0
        thrust = 0
        for c in msh.cells:
            fuel += c.fuel
            thrust += c.get_burn_rate(msh)

        if fuel == 0:
            running = False

        times.append(time)
        thrusts.append(thrust)
        fuels.append(fuel)

        # GRAPHICS
        for c in msh.cells:
            if c.fuel > 0:
                if not c.active:
                    fill_color = "gray"
                #elif c.get_burn_rate(msh) > 0:
                else:
                    fill_color = "red"

                cvas.create_oval(c.x*3-2 + 200, c.y*3-2 + 150, c.x*3+2 + 200, c.y*3+2 + 150, fill=fill_color)

        cvas.create_text(250, 50, text="Time: "+str(time), fill="red")
        cvas.create_text(250, 70, text="Fuel: "+str(fuel), fill="red")
        cvas.create_text(250, 90, text="Thrust: "+str(thrust), fill="red")

        root.update()
        cvas.delete("all")
        # time.sleep(dt)
        cycle += 1
        time = cycle * dt

    root.destroy()

    thrusts_f = low_pass_filter(thrusts, len(thrusts)//20)

    plt.plot(times, thrusts_f)
    plt.title("Thrust Profile")
    plt.xlabel("Time")
    plt.ylabel("Thrust")
    plt.grid()
    plt.show()

    plt.plot(times, fuels)
    plt.title("Fuel Mass Profile")
    plt.xlabel("Time")
    plt.ylabel("Fuel")
    plt.grid()
    plt.show()