Stratoballoon/vertical.py at main · MateJozsaPhys/Stratoballoon · 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
import numpy as np
import matplotlib.pyplot as plt

sim_length = 1000

z_0 = 0. #initial altitude
v_0 = 0. #initial velocity

#masses
m_gross = 1.
m_gas = 1.
m_added = 1.

m_tot = m_gross + m_gas + m_added

rho_air = 1. #air density

g = 10. #gravitational acceleration
C_d = 1. #drag coefficient
A_b = 1. #balloon reference area
A_p = 1. #parachute reference area

dt = 1. #time step

M_gas = 4*10**(-3) #molar mass (kg/mol) I have to express it this way to obtain the same dimensions everywhere
R = 8.314 #universal gas constant (J/(K*mol))

#pressure parameters
c_p = 1004.68506 #constant pressure specific heat (J/(kg/K))
T_0 = 288.16 #sea level standard temperature (kg)
M = 0.02896968 #molar mass of dry air (kg/mol)
p_0 = 101325 #sea level standard atmospheric pressure (Pa)  N/m^2

def pressure(z_n):
	return p_0*(1 - (g*z_n)/(c_p*T_0))**(c_p*M/R)


V_0 = 2.34 #initial volume of the balloon (m^3)
gamma = 5/3 #adiabatic coefficient for helium (ideal monoatomic gas!!!)

const = p_0*V_0**gamma #p*V^gamma = const

def volume(z_n):
	return (const/pressure(z_n))**(1/gamma)

#initial naive presumption (this is actually valid only up to 10 kms)
L = 0.0065 #temperature lapse rate (K/m)
def Temperature(z_n):
	return T_0 - L*z_n
def rho_air(z_n):
	return M/(R*Temperature(z_n))*pressure(z_n)

################################################################ Ascent ########################################################################

def velocity_ascending(z_n, v_n):
	#print(v_n)
	return (1/m_tot)*(m_tot*v_n + g*rho_air(z_n)*volume(z_n)*dt - g*(m_gross + m_gas)*dt - 1/2*(C_d*rho_air(z_n)*v_n**2*A_b*dt))

################################################################ Descent #######################################################################

def velocity_descending(z_n, v_n):
	return -np.sqrt(2*m_gross*g/(C_d*rho_air(z_n)*A_p))

def altitude(z_n, v_n):
	return z_n + v_n*dt

z_n = z_0
v_n = v_0

V = np.array([])
Z = np.array([])
for i in range(sim_length):
	V = np.append(V, v_n)
	Z = np.append(Z, z_n)
	v_np1 = velocity_ascending(z_n, v_n)
	z_np1 = altitude(z_n, v_n)
	v_n = v_np1
	z_n = z_np1

while(z_n > 0):
	V = np.append(V, v_n)
	Z = np.append(Z, z_n)
	v_np1 = velocity_descending(z_n, v_n)
	z_np1 = altitude(z_n, v_n)
	v_n = v_np1
	z_n = z_np1

print(V)
print(Z)

plt.plot(Z)
plt.grid()
plt.xlabel('time...sort of...')
plt.ylabel('altitude...kindof...')
plt.title('Altitude in time')
plt.show()