-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathheating.py
1776 lines (1526 loc) · 69.3 KB
/
heating.py
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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
Implementation of NASA program NQLDW019. All values are stored in SI units unless otherwise stated.
Reference material:
-------------------
[1] - HEATING ON SOUNDING ROCKET TANGENT OGIVE NOSES (NASA) - https://arc.aiaa.org/doi/pdf/10.2514/3.62081
[2] - TANGENT OGIVE NOSE AERODYNAMIC HEATING PROGRAM - NQLDW019 (NASA) - https://ntrs.nasa.gov/citations/19730063810
- There seems to be a typo in Equation (20) in Reference [1] (there's a missing bracket). The correct version seems to be present in Reference [2].
General Information:
-----------------------------
- Station 1 is the nosecone tip, 11 is at the nosecone base, and 12-15 are beyond the base.
- If the rocket is subsonic, or the post-oblique-shock flow is subsonic, the simulation skips the steps (NaN will usually be stored for the data for these steps, so it doesn't show up in graphs).
Assumptions:
-------------------
- Currently assumes an oblique shock in front of the tangent ogive nosecone - may be better to assume a cone shock.
- If a variable wall temperature is used, the nosecone is modelled as a very simple lumped mass for temperature rises.
- Uniform wall temperature throughout the nose cone.
- Zero angle of attack. (if non-zero angle of attack is important, it can be implemented relatively simply and is explained in https://ntrs.nasa.gov/citations/19730063810)
- Constant Cp and Cv (hence constant gamma)
Known issues:
-------------
- The thermo module seems to have trouble finding the viscosity at some point - I think in post normal-shockwave conditions.
- I have not checked my model for wall temperature rise against any real data or the NASA examples.
Things I wasn't sure of:
-------------------------
- For the lumped mass wall temperature rise model, I improvised a way to deal with the infinite qdot at the nose tip (but am not sure how valid it is)
I simply assumed the heat transfer rates at Station 1 (the nose tip) to be the same as those at Station 2
- Calculating for H*(0) - Equation (18) from https://arc.aiaa.org/doi/pdf/10.2514/3.62081
I think that the rho(x) and mu(x) in Equation (18) are just rho(0) and mu(0), since they're evaluated at 'x=0'. This seemed to give the right values.
- Calculating H*(x) - Equation (17) from https://arc.aiaa.org/doi/pdf/10.2514/3.62081
I think I did it right but the integration was a bit confusing
- Some equations only work in Imperial units. I believe I've taken care of them all, but if problems arise, that might be worth looking into
Nomenclature that is normally used:
-----------------------------------
T - Temperature (K)
rho - Density (kg/m^3)
P - pressure (Pa)
q or qdot - Heat transfer rate per unit area (W/m^2)
Q or Qdot - Heat transfer rate (W)
R - Gas constant (J/kg/K)
Cp or cp - Specific heat capacity at constant pressure (J/kg/K)
gamma - Ratio of specific heats (Cp/Cv)
Pr - Prandtl number
mu - Viscosity (Pa s)
k - Thermal conductivity (W/m/K, I think)
V - Velocity (m/s)
M - Mach number
Re - Reynolds number
nu - Prandtl-Meyer function evaluated at a given Mach and gamma
H - As defined in Equations (17) and (18) of Reference [1]
RN - Hemispherical nosecone radius (m)
dVdx0 - As defined in Equation (19) or Reference [1], units of (1/s) if I recall correctly
Subscripts:
-----------
e or x - "Local" value. This usually means it's taken at the edge of the boundary layer.
inf - Freestream (i.e. atmospheric) value.
ref or star - At the 'reference' enthalpy (usually marked with a star, e.g. T*, in the NASA documents). I'm not sure but I think this might be a sort of average boundary layer temperature.
0 - At the stagnation point for a hemispherical nose cone.
w - At wall temperature and local pressure.
rec - At 'recovery enthalpy', which is the same thing as at the 'adiabatic wall temperature' I believe.
turb - With a turbulent boundary layer
lam - With a laminar boundary layer.
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import thermo.mixture, matplotlib.widgets, json, scipy.integrate, scipy.optimize
import gas_dynamics as gd
from ambiance import Atmosphere
from .transforms import (
pos_l2i,
pos_i2l,
vel_l2i,
vel_i2l,
direction_l2i,
direction_i2l,
i2airspeed,
pos_i2alt,
)
__copyright__ = """
Copyright 2021 Daniel Gibbons
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <https://www.gnu.org/licenses/>.
"""
# Compressible flow functions
def prandtl_meyer(M, gamma=1.4):
"""Prandtl-Meyer function
Parameters
----------
M : float
Mach number
gamma : float, optional
Ratio of specific heats (cp / cv). Defaults to 1.4.
Returns
-------
float
nu (Prandtl-Meyer function evaluated at the given Mach and gamma)
"""
if M < 1:
raise ValueError(
"Cannot calculate the Prandtl-Meyer function of a flow with M < 1"
)
return float(
np.sqrt((gamma + 1) / (gamma - 1))
* np.arctan(np.sqrt((gamma - 1) / (gamma + 1) * (M ** 2 - 1)))
- np.arctan(np.sqrt(M ** 2 - 1))
)
def nu2mach(nu, gamma=1.4):
"""Inverse of the Prandtl-Meyer function
Notes
----------
Calculated using a polynomial approximation, described in http://www.pdas.com/pm.pdf
Parameters
----------
nu : float
Value of the Prandtl-Meyer function
gamma : float, optional
Ratio of specific heats (cp / cv). Defaults to 1.4.
Returns
-------
float
Mach number corresponding to the given value of the Prandtl Meyer function
"""
if gamma != 1.4:
raise ValueError("This function will only work for gamma = 1.4")
nuinf = (6 ** 0.5 - 1) * np.pi / 2
y = (nu / nuinf) ** (2 / 3)
A = 1.3604
B = 0.0962
C = -0.5127
D = -0.6722
E = -0.3278
return (1 + A * y + B * y ** 2 + C * y ** 3) / (1 + D * y + E * y ** 2)
def p2p0(P, M, gamma=1.4):
"""Returns static pressure from stagnation pressure, Mach number, and ratio of specific heats
Parameters
----------
P : float
Static pressure
M : float
Mach number
gamma : float, optional
Ratio of specific heats (cp / cv). Defaults to 1.4.
Returns
-------
float
Stagnation pressure
"""
return P * (1 + (gamma - 1) / 2 * M ** 2) ** (gamma / (gamma - 1))
def p02p(P0, M, gamma=1.4):
"""Returns static pressure from stagnation pressure, Mach number, and ratio of specific heats
Parameters
----------
P0 : float
Stagnation pressure
M : float
Mach number
gamma : float, optional
Ratio of specific heats (cp / cv). Defaults to 1.4.
Returns
-------
float
Static pressure
"""
return P0 * (1 + (gamma - 1) / 2 * M ** 2) ** (-gamma / (gamma - 1))
def T2T0(T, M, gamma=1.4):
"""Returns stagnation temperature from static temperature, Mach number, and ratio of specific heats
Parameters
----------
T : float
Static temperature
M : float
Mach number
gamma : float, optional
Ratio of specific heats (cp / cv). Defaults to 1.4.
Returns
-------
float
Stagnation temperature
"""
return T * (1 + (gamma - 1) / 2 * M ** 2)
def T02T(T0, M, gamma=1.4):
"""Returns static temperature from stagnation temperature, Mach number, and ratio of specific heats
Parameters
----------
T0 : float
Stagnation temperature
M : float
Mach number
gamma : float, optional
Ratio of specific heats (cp / cv). Defaults to 1.4.
Returns
-------
float
Static temperature
"""
return T0 * (1 + (gamma - 1) / 2 * M ** 2) ** (-1)
def rho2rho0(rho, M, gamma=1.4):
"""Returns stagnation density from static density, Mach number, and ratio of specific heats
Parameters
----------
T : float
Static density
M : float
Mach number
gamma : float, optional
Ratio of specific heats (cp / cv). Defaults to 1.4.
Returns
-------
float
Stagnation density
"""
return rho * (1 + (gamma - 1) / 2 * M ** 2) ** (1 / (gamma - 1))
def rho02rho(rho0, M, gamma=1.4):
"""Returns static density from stagnation density, Mach number, and ratio of specific heats
Parameters
----------
rho0 : float
Stagnation density
M : float
Mach number
gamma : float, optional
Ratio of specific heats (cp / cv). Defaults to 1.4.
Returns
-------
float
Static density
"""
return rho0 * (1 + (gamma - 1) / 2 * M ** 2) ** (-1 / (gamma - 1))
def pressure_ratio_to_mach(p_over_p0, gamma=1.4):
"""Get Mach number from the static to stagnation pressure ratio
Parameters
----------
p_over_p0 : float
Static pressure divided by stagnation pressure
gamma : float, optional
Ratio of specific heats (cp / cv). Defaults to 1.4.
Returns
-------
float
Mach number
"""
return ((2 / (gamma - 1)) * (p_over_p0 ** ((gamma - 1) / -gamma) - 1)) ** 0.5
def normal_shock(M, gamma=1.4):
"""Normal shock wave calculator
Parameters
----------
M : float
Mach number
gamma : float, optional
Ratio of specific heats (cp / cv). Defaults to 1.4.
Returns
-------
numpy ndarray
Returns array of floats in the following order: [MS, PS/P, TS/T, rhoS/rho]
"""
MS = (
(1 + 0.5 * (gamma - 1) * M ** 2) / (gamma * M ** 2 - 0.5 * (gamma - 1))
) ** 0.5
PSoverP = 1 + 2 * gamma / (gamma + 1) * (M ** 2 - 1)
TSoverT = (
(gamma - 1)
/ (gamma + 1) ** 2
* 2
/ M ** 2
* (1 + 0.5 * (gamma - 1) * M ** 2)
* (2 * gamma / (gamma - 1) * M ** 2 - 1)
)
rhoSoverrho = (gamma + 1) * M ** 2 / (2 * (1 + 0.5 * (gamma - 1) * M ** 2))
return np.array([MS, PSoverP, TSoverT, rhoSoverrho])
# Properties of air
def cp_air(T=298, P=1e5):
"""Specific heat capacity of air at constant pressure (J/kg/K)"""
# air = thermo.mixture.Mixture('air', T=T, P=P)
# return air.Cp
return 1005
def R_air():
"""Gas constant for air (J/kg/K)"""
return 287
def gamma_air(T=298, P=1e-5):
"""Ratio of specific heats (Cp/Cv) for air"""
return 1.4
def Pr_air(T, P):
"""Prandtl number for air"""
air = thermo.mixture.Mixture("air", T=T, P=P)
return air.Pr
# return 0.71
def k_air(T, P):
"""Thermal conductivity of air (W/m/K)"""
air = thermo.mixture.Mixture("air", T=T, P=P)
return air.k
# return 26.24e-3
def mu_air(T, P):
"""Viscosity of air (Pa s)"""
air = thermo.mixture.Mixture("air", T=T, P=P)
return air.mu
# return 1.81e-5
def oblique_shock(theta, M, T, p, rho, gamma=1.4, R=287):
"""Function for generating oblique shock data. Makes use of the gas_dynamics library.
Args:
theta (float): Deflection angle (rad)
M (float): Mach number before shock
T (float): Temperature before shock (K)
p (float): Pressure before shock (Pa)
rho (float): Density before shock (kg/m3)
gamma (float, optional): Ratios of specific heats (cp/cv). Defaults to 1.4.
R (int, optional): Specific gas constant (J/kg/K). Defaults to 287.
Returns:
MS, TS, PS, rhoS, beta: Post-shock Mach number, temperature, pressure, density and shock angle.
"""
# Get beta using the gas_dynamics library, and convert it to radians
beta = gd.shocks.shocks.shock_angle(M, theta * 180 / np.pi)[0] * np.pi / 180
# From CUED 3A3 Datasheet
PS = p * (1 + 2 * gamma / (gamma + 1) * (M ** 2 * np.sin(beta) ** 2 - 1))
rhoS = (
rho
* ((gamma + 1) * M ** 2 * np.sin(beta) ** 2)
/ (2 * (1 + (gamma - 1) / 2 * M ** 2 * np.sin(beta) ** 2))
)
MS = (
(1 + (gamma - 1) / 2 * M ** 2 * np.sin(beta) ** 2)
/ (gamma * M ** 2 * np.sin(beta) ** 2 - (gamma - 1) / 2)
) ** 0.5 / np.sin(beta - theta)
# Ideal gas - p = rho R T
TS = PS / (rhoS * R)
return MS, TS, PS, rhoS, beta
# Aerodynamic heating analysis
class TangentOgive:
"""
Object used to store the geometry of a tangent ogive nose cone.
Inputs
-------
xprime : float
Longitudinal dimension (tip-to-base distance) (m)
yprime : float
Base radius (m)
"""
def __init__(self, xprime, yprime):
# https://arc.aiaa.org/doi/pdf/10.2514/3.62081 used for nomenclature
self.xprime = xprime # Longitudinal dimension
self.yprime = yprime # Base radius
self.R = (xprime ** 2 + yprime ** 2) / (2 * yprime) # Ogive radius
self.theta = np.arctan2(xprime, self.R - yprime)
self.dtheta = 0.1 * self.theta
# Each point (1 to 15) and its distance along the nose cone surface from the nose tip
self.S_array = np.zeros(15)
for i in range(len(self.S_array)):
self.S_array[i] = self.S(i + 1)
def phi(self, i):
# i = 1 to 15
assert (
i >= 1 and i <= 15
), "i refers to stations 1-15, it cannot the less than 1 or more than 15"
return self.theta - (i - 1) * self.dtheta / 2
def r(self, i):
"""
Local nosecone radius (y-dimension) at station
Inputs
-------
i : int
Station number (1-15)
Returns
-------
float
Local nosecone radius (m)
"""
# i = 1 to 15
assert (
i >= 1 and i <= 15
), "i refers to stations 1-15, it cannot the less than 1 or more than 15"
if i <= 11:
return 2 * self.R * np.sin((i - 1) * self.dtheta / 2) * np.sin(self.phi(i))
else:
return 2 * self.R * np.sin((10) * self.dtheta / 2) * np.sin(self.phi(11))
def S(self, i):
"""
Distance from
Inputs
-------
i : int
Station number (1-15)
Returns
--------
float
Distance along the nosecone surface (m), from the nosecone tip to station i.
"""
# i = 1 to 15
assert (
i >= 1 and i <= 15
), "i refers to stations 1-15, it cannot the less than 1 or more than 15"
return self.R * (i - 1) * self.dtheta
class AeroHeatingAnalysis:
"""
Object used to run aerodynamic heating analyses
Inputs
-------
tangent_ogive : TangentOgive
TangentOgive object specifying the nosecone geometry
trajectory_data : dict or pandas DataFrame
Data on the rocket's trajectory, needs to have "pos_i", "vel_i" and "time".
rocket : Rocket
Rocket object. It's only needed to get LaunchSite data for coordinate transforms.
fixed_wall_temperature : bool
If True, the wall temperature is fixed to its starting value. Otherwise a simple model is used to model its temperature change.
starting_temperature : float, optional
Temperature that the nose cone starts with (K). Defaults to None, in which case the rocket starts with the local atmospheric temperature.
nosecone_mass : float, optional
Mass of the nosecone (kg) - used to find its heat capacity. Only needed if you're modelling variable temperatures
specific_heat_capacity : float, optional
Specific heat capacity of the nosecone (J/kg/K). Defaults to an approximate value for aluminium.
turbulent_transition_Rex : float, optional
Local Reynolds number at which the boundary layer transition from laminar to turbulent. Defaults to 7.5e6.
Attributes
----------
i : int
Current index in the timestep array.
tangent_ogive : TangentOgive
A TangentOgive object containing nosecone geometry data
trajectory_data : dict or pandas DataFrame
Contains trajectory data.
trajectory_dict : dict
Same data trajectory_data, but converted to a dictionary if it wasn't already one.
rocket : Rocket
Rocket object used to run the simulation. Its only purpose is to get LaunchSite data for coordinate transforms.
fixed_wall_temperature : bool
If True, the wall temperature is fixed to its starting value. Otherwise a simple model is used to model its temperature change.
turbulent_transition_Rex : float, optional
Local Reynolds number at which the boundary layer transition from laminar to turbulent. Defaults to 7.5e6.
heat_capacity : float
Heat capacity the nosecone (J/K). Is caclulated by doing specific_heat_capacity*nosecone_mass.
M : numpy ndarray
Local Mach number at each station and timestep. Has dimensions (15, N) where N is the length of the "time" array in trajectory_dict.
P : numpy ndarray
Local static pressure at each station and timestep. Has dimensions (15, N) where N is the length of the "time" array in trajectory_dict.
Te : numpy ndarray
Temperature at the edge of the boundary layer, at each station and timestep. Has dimensions (15, N) where N is the length of the "time" array in trajectory_dict.
Tw : numpy ndarray
Wall temperature at each timestep. Local static pressure at each station and timestep. Has dimensions (N) where N is the length of the "time" array in trajectory_dict.
Trec_lam : numpy ndarray
Adiabatic wall tempearture (also known as the 'recovery temperature') for a laminar boundary layer, at each station and timestep. Has dimensions (15, N) where N is the length of the "time" array in trajectory_dict.
Trec_turb : numpy ndarray
Adiabatic wall tempearture (also known as the 'recovery temperature') for a turbulent boundary layer, at each station and timestep. Has dimensions (15, N) where N is the length of the "time" array in trajectory_dict.
Hstar_function : numpy ndarray
Function that needs to be integrated to get H*(x), as defined in the NASA documents. This is not equal to H*(x) itself. This is stored because it's needed for integration. Has dimensions (15, N) where N is the length of the "time" array in trajectory_dict.
Rex : numpy ndarray
Local Reynolds number at each station and timestep. Has dimensions (15, N) where N is the length of the "time" array in trajectory_dict.
q_lam : numpy ndarray
Heat transfer rate (W/m^2) with a laminar boundary layer, at each station and timestep. Has dimensions (15, N) where N is the length of the "time" array in trajectory_dict.
q_turb : numpy ndarray
Heat transfer rate (W/m^2) with a turbulent boundary layer, at each station and timestep. Has dimensions (15, N) where N is the length of the "time" array in trajectory_dict.
q0_hemispherical_nose : numpy ndarray
Heat transfer rate (W/m^2) at the stagnation point at each timestep, if we were using a hemispherical nosecone. Has dimensions (N) where N is the length of the "time" array in trajectory_dict.
"""
def __init__(
self,
tangent_ogive,
trajectory_data,
rocket,
fixed_wall_temperature=True,
starting_temperature=None,
nosecone_mass=None,
specific_heat_capacity=900,
turbulent_transition_Rex=7.5e6,
):
self.tangent_ogive = tangent_ogive
self.trajectory_data = trajectory_data
self.rocket = rocket
self.turbulent_transition_Rex = turbulent_transition_Rex
self.fixed_wall_temperature = fixed_wall_temperature
# Convert the data into a dictionary if it isn't already one:
if type(self.trajectory_data) is dict:
self.trajectory_dict = self.trajectory_data
else:
self.trajectory_dict = self.trajectory_data.to_dict(orient="list")
# If we want a variable wall temperature:
if self.fixed_wall_temperature == False:
assert (
nosecone_mass != None
), "You need to input a value for the nosecone mass if you want to model a variable wall temperature"
self.heat_capacity = specific_heat_capacity * nosecone_mass
# Timestep index:
self.i = 0
# Arrays to store the fluid properties at each discretised point on the nose cone (1 to 15), and at each timestep
self.M = np.full(
[15, len(self.trajectory_dict["time"])], float("NaN")
) # Local Mach number
self.P = np.full(
[15, len(self.trajectory_dict["time"])], float("NaN")
) # Local pressure
# Initialise the wall temperature:
if starting_temperature == None:
starting_temperature = Atmosphere(
pos_i2alt(
self.trajectory_dict["pos_i"][0], self.trajectory_dict["time"][0]
)
).temperature[
0
] # Assume the nose cone starts with ambient temperature
if self.fixed_wall_temperature == True:
self.Tw = np.full(len(self.trajectory_dict["time"]), starting_temperature)
else:
self.Tw = np.full(len(self.trajectory_dict["time"]), float("NaN"))
self.Tw[0] = starting_temperature
self.Te = np.full(
[15, len(self.trajectory_dict["time"])], float("NaN")
) # Temperature at the edge of the boundary layer
self.Tstar = np.full(
[15, len(self.trajectory_dict["time"])], float("NaN")
) # T* as defined in the paper
self.Trec_lam = np.full(
[15, len(self.trajectory_dict["time"])], float("NaN")
) # Temperature corresponding to hrec_lam
self.Trec_turb = np.full(
[15, len(self.trajectory_dict["time"])], float("NaN")
) # Temperature corresponding to hrec_turb
self.Hstar_function = np.full(
[15, len(self.trajectory_dict["time"])], float("NaN")
) # This array is used to minimise number of calculations for the integration needed in H*(x)
self.Rex = np.full(
[15, len(self.trajectory_dict["time"])], float("NaN")
) # Local Reynolds number
# Arrays to store the heat transfer rates
self.q_lam = np.full(
[15, len(self.trajectory_dict["time"])], float("NaN")
) # Laminar boundary layer
self.q_turb = np.full(
[15, len(self.trajectory_dict["time"])], float("NaN")
) # Turbunulent boundary layer
self.q0_hemispherical_nose = np.full(
len(self.trajectory_dict["time"]), float("NaN")
) # At the stagnation point for a rocket with a hemispherical nose cone - used as a reference point
def step(self, print_style=None):
"""
Performs one step of the aerodynamic analysis, starting from the current value of self.i.
Inputs:
-------
print_style : str
Options for print style:
None - nothing is printed
"FORTRAN" - same output as the examples in https://ntrs.nasa.gov/citations/19730063810, printing in the Imperial units listed
"metric" - outputs useful data in metric units
"""
# Get altitude:
alt = pos_i2alt(
self.trajectory_dict["pos_i"][self.i], self.trajectory_dict["time"][self.i]
)
if alt < -5004: # hacky way to fix out of bound altitudes for ambience
alt = 5004
elif alt > 81020:
alt = 81020
# Get ambient conditions:
Pinf = Atmosphere(alt).pressure[0]
Tinf = Atmosphere(alt).temperature[0]
rhoinf = Atmosphere(alt).density[0]
# Get the freestream velocity and Mach number
Vinf = np.linalg.norm(
i2airspeed(
self.trajectory_dict["pos_i"][self.i],
self.trajectory_dict["vel_i"][self.i],
self.rocket.launch_site,
self.trajectory_dict["time"][self.i],
)
)
Minf = Vinf / Atmosphere(alt).speed_of_sound[0]
if print_style == "FORTRAN":
print("")
print("FREE STREAM CONDITIONS")
print(
"XMINF={:<10} VINFY={:.4e} GAMINF={:.4e} RHOINF={:.4e}".format(
0, 3.28084 * Vinf, gamma_air(), 0.00194032 * rhoinf
)
)
print(
"HINFY={:.4e} PINF ={:.4e} (ATMOS) PINFY ={:.4e} (PSF)".format(
0.000429923 * cp_air() * Tinf, Pinf / 101325, 0.0208854 * Pinf
)
)
print("TINFY={:.4e}".format(Tinf))
print("")
if print_style == "metric":
print("")
print("SUBCRIPTS:")
print("0 or STAG : At the stagnation point for a hemispherical nose")
print(
"REF : At 'reference' enthalpy and local pressure - I think this is like an average-ish boundary layer enthalpy"
)
print(
"REC : At 'recovery' enthalpy and local pressure - I believe this is the wall enthalpy at which no heat transfer takes place"
)
print("W : At the wall temperature and local pressure")
print("INF : Freestream (i.e. atmospheric) property")
print("LAM : With a laminar boundary layer")
print("TURB : With a turbulent boundary layer")
print("")
print("FREE STREAM CONDITIONS")
print(
"ALT ={:06.2f} km TINF={:06.2f} K PINF={:06.2f} kPa RHOINF={:06.2f} kg/m^3".format(
alt / 1000, Tinf, Pinf / 1000, rhoinf
)
)
print("VINF={:06.2f} m/s MINF={:06.2f}".format(Vinf, Minf))
print("")
# Check if we're supersonic - if so we'll have a shock wave
if Minf > 1:
# For an oblique shock (tangent ogive nose cone)
oblique_shock_data = oblique_shock(
self.tangent_ogive.theta, Minf, Tinf, Pinf, rhoinf
) # MS, TS, PS, rhoS, beta
oblique_MS = oblique_shock_data[0]
oblique_TS = oblique_shock_data[1]
oblique_PS = oblique_shock_data[2]
oblique_rhoS = oblique_shock_data[3]
oblique_P0S = p2p0(oblique_PS, oblique_MS)
oblique_T0S = T2T0(oblique_TS, oblique_MS)
oblique_rho0S = rho2rho0(oblique_rhoS, oblique_MS)
# For a normal shock (hemispherical nosecone)
normal_shock_data = normal_shock(Minf)
normal_MS = normal_shock_data[0]
normal_PS = normal_shock_data[1] * Pinf
normal_TS = normal_shock_data[2] * Tinf
normal_rhoS = normal_shock_data[3] * rhoinf
normal_P0S = p2p0(normal_PS, normal_MS)
normal_T0S = T2T0(normal_TS, normal_MS)
normal_rho0S = rho2rho0(normal_rhoS, normal_MS)
# Stagnation point heat transfer rate for a hemispherical nosecone
Pr0 = Pr_air(normal_T0S, normal_P0S)
h0 = cp_air() * normal_T0S
hw = cp_air() * self.Tw[self.i]
mu0 = mu_air(normal_T0S, normal_P0S)
rhow0 = normal_P0S / (
R_air() * self.Tw[self.i]
) # p = rho * R * T (ideal gas)
muw0 = mu_air(self.Tw[self.i], normal_P0S)
RN = 0.3048 # Let RN = 1 ft = 0.3048m, as it recommends using that as a reference value (although apparently it shouldn't matter?)
dVdx0 = (2 ** 0.5) / RN * ((normal_P0S - Pinf) / normal_rho0S) ** 0.5
# Equation (29) from https://arc.aiaa.org/doi/pdf/10.2514/3.62081
# Note that the equation only works in Imperial units, and requires you to specify density in slugs/ft^3, which is NOT lbm/ft^3
# Metric density (kg/m^3) --> Imperial density (slugs/ft^3): Multiply by 0.00194032
# Metric viscosity (Pa s) --> Imperial viscosity (lbf sec/ft^2): Divide by 47.880259
# Metric enthalpy (J/kg/s) --> Imperial enthalpy (Btu/lbm): Multiply by 0.000429923
# Note that 'g', the acceleration of gravity, is equal to 32.1740 ft/s^2
self.q0_hemispherical_nose[self.i] = (
0.76
* 32.1740
* Pr0 ** (-0.6)
* (0.00194032 * rhow0 * muw0 / 47.880259) ** 0.1
* (0.00194032 * normal_rho0S * mu0 / 47.880259) ** 0.4
* (0.000429923 * h0 - 0.000429923 * hw)
* dVdx0 ** 0.5
)
# Now convert from Imperial heat transfer rate (Btu/ft^2/s) --> Metric heat transfer rate (W/m^2): Divide by 0.000088055
self.q0_hemispherical_nose[self.i] = (
self.q0_hemispherical_nose[self.i] / 0.000088055
)
if print_style == "FORTRAN":
print("")
print("STAGNATION POINT DATA FOR SPHERICAL NOSE")
print(
"HREF0 ={:<10} TREF0 ={:<10} VISCR0={:<10} TKREF0={:<10}".format(
0, 0, 0, 0
)
)
print(
"ZREF0 ={:<10} PRREF0={:<10} CPREF0={:<10} RHOR0 ={:<10}".format(
0, 0, 0, 0
)
)
print(
"CPCVR0={:.4e} RN ={:.4e} T0 ={:.4e}".format(
gamma_air(), RN / 0.3048, normal_T0S
)
)
print(
"P0 ={:.4e} RHO0 ={:.4e} SR0 ={:<10} TK0 ={:<10}".format(
normal_P0S / 101325, 0.00194032 * normal_rho0S, 0, 0
)
)
print(
"VISC0 ={:.4e} DVDX0 ={:.4e} Z0 ={:<10} CP0 ={:.4e}".format(
mu0 / 47.880259, dVdx0, 0, 0.000429923 * cp_air()
)
)
print(
"A0 ={:<10} TW0 ={:.4e} VISCW0={:.4e} HW0 ={:.4e}".format(
0, self.Tw[self.i], muw0 / 47.880259, 0.000429923 * hw
)
)
print("")
print(
"CPW0 ={:.4e} PR0 ={:.4e}".format(
0.000429923 * cp_air(), Pr0
)
)
print(
"QSTPT ={:.4e} = NOSE STAGNATION POINT HEAT RATE".format(
0.000088055 * self.q0_hemispherical_nose[self.i]
)
)
print(
"H0 ={:.4e} HT ={:<10} RHOW0={:.4e}".format(
0.000429923 * h0, 0, 0.00194032 * rhow0
)
)
print("")
if print_style == "metric":
print("")
print("STAGNATION POINT DATA FOR SPHERICAL NOSE")
print(
"P0 ={:06.2f} kPa T0 ={:06.2f} K RHO0={:06.2f} kg/m^3".format(
normal_P0S / 1000, normal_T0S, normal_rho0S
)
)
print(
"TW ={:06.2f} K RHOW0={:06.2f} kg/m^3".format(
self.Tw[self.i], rhow0
)
)
print(
"QSTAG={:06.2f} kW/m^2".format(
self.q0_hemispherical_nose[self.i] / 1000
)
)
print("")
# Prandtl-Meyer expansion (only possible for supersonic flow):
if oblique_MS > 1:
# Get values at the nose cone tip:
nu1 = prandtl_meyer(oblique_MS)
theta1 = self.tangent_ogive.theta
# Prandtl-Meyer expansion from post-shockwave to each discretised point
for j in range(10):
# Angle between the flow and the horizontal:
theta = self.tangent_ogive.theta - self.tangent_ogive.dtheta * j
# Across a +mu characteristic: nu1 + theta1 = nu2 + theta2
nu = nu1 + theta1 - theta
# Check if we've exceeded nu_max, in which case we can't turn the flow any further
if nu > (np.pi / 2) * (
np.sqrt((gamma_air() + 1) / (gamma_air() - 1)) - 1
):
raise ValueError(
"Cannot turn flow any further at nosecone position {}, exceeded nu_max. Flow will have seperated (which is not yet implemented). Stopping simulation.".format(
j + 1
)
)
# Record the local Mach number and pressure
self.M[j, self.i] = nu2mach(nu)
self.P[j, self.i] = p02p(oblique_P0S, self.M[j, self.i])
# Expand for the last few points using Equations (1) - (6) from https://arc.aiaa.org/doi/pdf/10.2514/3.62081
for j in [10, 11, 12, 13, 14]:
if j >= 10 and j <= 13:
self.P[j, self.i] = (Pinf + self.P[j - 1, self.i]) / 2
elif j == 14:
self.P[j, self.i] = Pinf
self.M[j, self.i] = pressure_ratio_to_mach(
self.P[j, self.i] / oblique_P0S
)
# Now deal with the heat transfer itself
for j in range(15):
# Edge of boundary layer temperature - i.e. flow temperature post-shock and after Prandtl-Meyer expansion
self.Te[j, self.i] = T02T(oblique_T0S, self.M[j, self.i])
# Enthalpies
he = cp_air() * self.Te[j, self.i]
# Prandtl numbers and specific heat capacities
Pre = Pr_air(self.Te[j, self.i], self.P[j, self.i])
#'Reference' values, as defined in https://arc.aiaa.org/doi/pdf/10.2514/3.62081 page 3
hstar = (he + hw) / 2 + 0.22 * (Pre ** 0.5) * (h0 - hw)
self.Tstar[j, self.i] = hstar / cp_air()
Prstar = Pr_air(self.Tstar[j, self.i], self.P[j, self.i])
#'Recovery' values, as defined in https://arc.aiaa.org/doi/pdf/10.2514/3.62081 page 3 - I think these are the wall enthalpies for zero heat transfer
hrec_lam = he * (1 - Prstar ** (1 / 2)) + h0 * (Prstar ** (1 / 2))
hrec_turb = he * (1 - Prstar ** (1 / 3)) + h0 * (Prstar ** (1 / 3))
self.Trec_lam[j, self.i] = hrec_lam / cp_air()
self.Trec_turb[j, self.i] = hrec_turb / cp_air()
# Get H*(x) - I'm not sure about if I did the integral bit right
rhostar0 = normal_P0S / (R_air() * self.Tstar[j, self.i])
mustar0 = mu_air(T=self.Tstar[j, self.i], P=normal_P0S)
rhostar = self.P[j, self.i] / (R_air() * self.Tstar[j, self.i])
mustar = mu_air(T=self.Tstar[j, self.i], P=self.P[j, self.i])
r = self.tangent_ogive.r(j + 1)
V = (
gamma_air() * R_air() * T02T(oblique_T0S, self.M[j, self.i])
) ** 0.5 * self.M[j, self.i]
self.Hstar_function[j, self.i] = (rhostar * mustar * V * r ** 2) / (
rhostar0 * mustar0 * Vinf
)
# Get the integral bit of H*(x) using trapezium rule
integral = np.trapz(
self.Hstar_function[0 : j + 1, self.i],
self.tangent_ogive.S_array[0 : j + 1],
)
# Equation (17) from https://arc.aiaa.org/doi/pdf/10.2514/3.62081
if j == 0:
Hstar = np.inf
else:
Hstar = (
(rhostar * V * r) / (rhostar0 * Vinf) / (integral ** 0.5)
)
# Get H*(0) - Equation (18) - it seems like the (x) values in Equation (18) are actually (0) values
# It seems weird that they still included them though, since they end up cancelling out
# Hstar0 = ( ((2*rhostar/rhostar0)*dVdx0 )/(Vinf * mustar/mustar0) )**0.5 * (2)**0.5
Hstar0 = (2 * dVdx0 / Vinf) ** 0.5 * (2 ** 0.5)
# Laminar heat transfer rate, normalised by that for a hemispherical nosecone
kstar = k_air(T=self.Tstar[j, self.i], P=self.P[j, self.i])
kstar0 = k_air(T=self.Tstar[j, self.i], P=normal_P0S)
Cpw = cp_air()
Cpw0 = cp_air()
# Equation (13) from https://arc.aiaa.org/doi/pdf/10.2514/3.62081 - wasn't sure which 'hrec' to use here but I think it's the laminar one
qxq0_lam = (kstar * Hstar * (hrec_lam - hw) * Cpw0) / (
kstar0 * Hstar0 * (h0 - hw) * Cpw
)
# Now we can find the absolute laminar heat transfer rates, in W/m^2
self.q_lam[j, self.i] = (
qxq0_lam * self.q0_hemispherical_nose[self.i]
)
# Turbulent heat transfer rate - using Equation (20) from https://arc.aiaa.org/doi/pdf/10.2514/3.62081
# THERE LOOKS LIKE THERE'S A TYPO IN THE EQUATION! It should be {1 - (Pr*)^(1/3)}*he, not 1 - {(Pr*)^(1/3)}*he
# For the correct version see page Eq (20) on page 14 of https://ntrs.nasa.gov/citations/19730063810
# Note that the equation only works in Imperial units, and requires you to specify density in slugs/ft^3, which is NOT lbm/ft^3
# Density (kg/m^3) --> Density (slugs/ft^3): Multiply by 0.00194032
# Viscosity (Pa s) --> Viscosity (lbf sec/ft^2): Divide by 47.880259
# Enthalpy (J/kg/s) --> Enthalpy (Btu/lbm): Multiply by 0.000429923
# Thermal conductivity (W/m/K) --> Thermal conductivity (Btu/ft/s/K): Multiply by 0.000288894658
# Velocity (m/s) ---> Velocity (ft/s): Multiply by 3.28084
# Note that 'g', the acceleration of gravity, is equal to 32.1740 ft/s^2