Skip to content

aero.ShockWave

@package ShockWave local Rankine Hugoniot equations for shock waves

Functions

Mn_Ps_ratio(Pratio, gamma=defg._gamma)

Computes normal Mach number from static pressure ratio across a shockwave

Parameters:

Name Type Description Default
Pratio float

static pressure ratio

 required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Type Description

normal Mach number in shock reference frame
Source code in aerokit/aero/ShockWave.py 
30
31
32
33
34
35
36
37
38
39
40
41
def Mn_Ps_ratio(Pratio: float, gamma=defg._gamma):
    """Computes normal Mach number from static pressure ratio across a shockwave

    Args:
      Pratio: static pressure ratio
      gamma:  (Default value = defg._gamma)

    Returns:
        normal Mach number in shock reference frame

    """
    return np.sqrt(1 + (Pratio - 1.0) * (gamma + 1.0) / (2.0 * gamma))

Mn_Pt_ratio(ptratio, gamma=defg._gamma)

Parameters:

Name Type Description Default
ptratio
required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Source code in aerokit/aero/ShockWave.py 
 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
def Mn_Pt_ratio(ptratio, gamma=defg._gamma):
    """

    Args:
      ptratio:
      gamma:  (Default value = defg._gamma)

    Returns:

    """

    def ptratio_of_mach(m):
        """

        Args:
          m:

        Returns:

        """
        return Pt_ratio(m, gamma)

    if ptratio > 1:
        print("!!! cannot find Mn for Pt_ratio > 1")
        return 1
    return ITS.secant_solve(ptratio_of_mach, ptratio, 1.5)

Ps_ratio(Mn, gamma=defg._gamma)

Computes static pressure ratio across a shockwave

Parameters:

Name Type Description Default
Mn

upstream normal Mach number in shock reference frame

 required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Type Description

static pressure ratio
Source code in aerokit/aero/ShockWave.py 
16
17
18
19
20
21
22
23
24
25
26
27
def Ps_ratio(Mn, gamma=defg._gamma):
    """Computes static pressure ratio across a shockwave

    Args:
      Mn: upstream normal Mach number in shock reference frame
      gamma:  (Default value = defg._gamma)

    Returns:
        static pressure ratio

    """
    return 1.0 + (2.0 * gamma / (gamma + 1.0)) * (Mn ** 2 - 1.0)

Pt_ratio(Mn, gamma=defg._gamma)

Total pressure ration through shock wave

Parameters:

Name Type Description Default
Mn

upstream relative normal Mach number to param

 required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Source code in aerokit/aero/ShockWave.py 
83
84
85
86
87
88
89
90
91
92
93
def Pt_ratio(Mn, gamma=defg._gamma):
    """Total pressure ration through shock wave

    Args:
      Mn: upstream relative normal Mach number to param
      gamma:  (Default value = defg._gamma)

    Returns:

    """
    return Ps_ratio(Mn, gamma) * Is.PtPs_Mach(downstream_Mn(Mn, gamma)) / Is.PtPs_Mach(Mn, gamma)

Rho_ratio(Mn, gamma=defg._gamma)

Parameters:

Name Type Description Default
Mn

normal Mach number

 required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Source code in aerokit/aero/ShockWave.py 
44
45
46
47
48
49
50
51
52
53
54
def Rho_ratio(Mn, gamma=defg._gamma):
    """

    Args:
      Mn: normal Mach number
      gamma:  (Default value = defg._gamma)

    Returns:

    """
    return ((gamma + 1.0) * Mn ** 2) / (2.0 + (gamma - 1.0) * Mn ** 2)

Ts_ratio(Mn, gamma=defg._gamma)

Parameters:

Name Type Description Default
Mn

normal Mach number

 required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Source code in aerokit/aero/ShockWave.py 
57
58
59
60
61
62
63
64
65
66
67
def Ts_ratio(Mn, gamma=defg._gamma):
    """

    Args:
      Mn: normal Mach number
      gamma:  (Default value = defg._gamma)

    Returns:

    """
    return Ps_ratio(Mn, gamma) / Rho_ratio(Mn, gamma)

conical_Mach_walldeflection_sigma(deflection, sigma, gamma=defg._gamma)

computes upstream Mach number from wall deflection and shock angle sigma

Parameters:

Name Type Description Default
Mach

param deflection:

 required
gamma

Default value = defg._gamma)

 _gamma
deflection
required

Returns:

Source code in aerokit/aero/ShockWave.py 
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
def conical_Mach_walldeflection_sigma(deflection, sigma, gamma=defg._gamma):
    """computes upstream Mach number from wall deflection and shock angle sigma

    Args:
      Mach: param deflection:
      gamma: Default value = defg._gamma)
      deflection:

    Returns:

    """
    assert sigma > deflection

    def local_def(mach):
        """internal wrapping function to iterative solve"""
        # print("local mach sigma", mach, sigma)
        return conical_deflection_Mach_sigma(mach, sigma, gamma)

    return ITS.secant_solve(local_def, deflection, 1.0 / degree.sin(sigma - 0.5 * deflection))

conical_deflection_Mach_sigma(Mach, sigma, gamma=defg._gamma, tol=1e-06)

Parameters:

Name Type Description Default
Mach

param sigma:

 required
gamma

Default value = defg._gamma)

 _gamma
tol

Default value = 1.0e-6)

 1e-06
sigma
required

Returns:

Source code in aerokit/aero/ShockWave.py 
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
def conical_deflection_Mach_sigma(Mach, sigma, gamma=defg._gamma, tol=1.0e-6):
    """

    Args:
      Mach: param sigma:
      gamma: Default value = defg._gamma)
      tol: Default value = 1.0e-6)
      sigma:

    Returns:

    """

    def rhs(phi, data):
        """RHS for conical equations (2 equations)

        Args:
          phi: angle (coordinate of the ODE)
          data: theta (flow angle) and Mach number

        Returns: theta/mach variation rate (RHS)
        """
        th = data[0]
        ma = data[1]
        k = 1.0 - (ma * degree.sin(phi - th)) ** 2
        rhs = np.zeros(2)
        rhs[0] = -degree.sin(th) * degree.cos(phi - th) / degree.sin(phi) / k
        rhs[1] = (
            np.pi
            / 180.0
            * degree.sin(th)
            * degree.sin(phi - th)
            / degree.sin(phi)
            / k
            * ma
            * (1.0 + 0.5 * (gamma - 1) * ma * ma)
        )
        return rhs

    # initial data for integration from downstream shock to wall
    th = deflection_Mach_sigma(Mach, sigma, gamma)
    ma = downstream_Mn(Mach * degree.sin(sigma), gamma) / degree.sin(sigma - th)
    phi = sigma
    thma = np.array([th, ma])
    h = -phi / 10.0
    conv = phi
    while abs(conv) >= tol:
        dthma, err = _rkf45(rhs, phi, thma, h)
        # Accept integration step if error e is within tolerance
        if err <= tol:
            conv = thma[0] - phi
            if conv / (h - dthma[0]) < 1:
                h = h * conv / (h - dthma[0])
            else:
                thma = thma + dthma
                phi = phi + h
                # print(Mach, sigma, phi, thma)
        else:
            h = 0.9 * h * (tol / err) ** 0.2
            # print("new h ",h)

    deflection = 0.5 * (thma[0] + phi)
    return deflection

conical_sigma_Mach_walldeflection(Mach, deflection, gamma=defg._gamma)

computes shock angle sigma from upstream Mach number and wall deflection

Parameters:

Name Type Description Default
Mach

param deflection:

 required
gamma

Default value = defg._gamma)

 _gamma
deflection
required

Returns:

Source code in aerokit/aero/ShockWave.py 
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
def conical_sigma_Mach_walldeflection(Mach, deflection, gamma=defg._gamma):
    """computes shock angle sigma from upstream Mach number and wall deflection

    Args:
      Mach: param deflection:
      gamma: Default value = defg._gamma)
      deflection:

    Returns:

    """

    def local_def(sig):
        """internal wrapping function to iterative solve"""
        return conical_deflection_Mach_sigma(Mach, sig, gamma)

    return ITS.secant_solve(local_def, deflection, degree.asin(1.0 / Mach) + deflection)

deflection_Mach_ShockPsratio(Mach, Pratio, gamma=defg._gamma)

Parameters:

Name Type Description Default
Mach

param Pratio:

 required
gamma

Default value = defg._gamma)

 _gamma
Pratio
required

Returns:

Source code in aerokit/aero/ShockWave.py 
143
144
145
146
147
148
149
150
151
152
153
154
def deflection_Mach_ShockPsratio(Mach, Pratio, gamma=defg._gamma):
    """

    Args:
      Mach: param Pratio:
      gamma: Default value = defg._gamma)
      Pratio:

    Returns:

    """
    return deflection_Mach_sigma(Mach, degree.asin(Mn_Ps_ratio(Pratio, gamma) / Mach), gamma)

deflection_Mach_sigma(Mach, sigma, gamma=defg._gamma)

Parameters:

Name Type Description Default
Mach

param sigma:

 required
gamma

Default value = defg._gamma)

 _gamma
sigma
required

Returns:

Source code in aerokit/aero/ShockWave.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
def deflection_Mach_sigma(Mach, sigma, gamma=defg._gamma):
    """

    Args:
      Mach: param sigma:
      gamma: Default value = defg._gamma)
      sigma:

    Returns:

    """
    sigrad = np.radians(sigma)
    Mn = Mach * np.sin(sigrad)
    return sigma - np.degrees(np.arctan2(np.tan(sigrad), Rho_ratio(Mn, gamma)))

dev_Max(Mach, gamma=defg._gamma)

computes the maximum deviation (always subsonic downstream flow), separation of weak/strong shock

Parameters:

Name Type Description Default
Mach

param gamma: (Default value = defg._gamma)

 required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Source code in aerokit/aero/ShockWave.py 
256
257
258
259
260
261
262
263
264
265
266
def dev_Max(Mach, gamma=defg._gamma):
    """computes the maximum deviation (always subsonic downstream flow), separation of weak/strong shock

    Args:
      Mach: param gamma:  (Default value = defg._gamma)
      gamma:  (Default value = defg._gamma)

    Returns:

    """
    return deflection_Mach_sigma(Mach, sigma_DevMax(Mach, gamma=gamma), gamma=gamma)

dev_Sonic(Mach, gamma=defg._gamma)

computes the deviation angle for a downstream SONIC Mach number

Parameters:

Name Type Description Default
Mach

param gamma: (Default value = defg._gamma)

 required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Source code in aerokit/aero/ShockWave.py 
269
270
271
272
273
274
275
276
277
278
279
def dev_Sonic(Mach, gamma=defg._gamma):
    """computes the deviation angle for a downstream SONIC Mach number

    Args:
      Mach: param gamma:  (Default value = defg._gamma)
      gamma:  (Default value = defg._gamma)

    Returns:

    """
    return deflection_Mach_sigma(Mach, sigma_Sonic(Mach, gamma=gamma), gamma=gamma)

downstreamMach_Mach_ShockPsratio(Mach, Pratio, gamma=defg._gamma)

Parameters:

Name Type Description Default
Mach

param Pratio:

 required
gamma

Default value = defg._gamma)

 _gamma
Pratio
required

Returns:

Source code in aerokit/aero/ShockWave.py 
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
def downstreamMach_Mach_ShockPsratio(Mach, Pratio, gamma=defg._gamma):
    """

    Args:
      Mach: param Pratio:
      gamma: Default value = defg._gamma)
      Pratio:

    Returns:

    """
    Mn0 = Mn_Ps_ratio(Pratio, gamma)
    sig = degree.asin(Mn0 / Mach)
    dev = deflection_Mach_sigma(Mach, sig, gamma)
    Mn1 = downstream_Mn(Mn0, gamma)
    return Mn1 / degree.sin(sig - dev)

downstream_Mn(Mn, gamma=defg._gamma)

Parameters:

Name Type Description Default
Mn

param gamma: (Default value = defg._gamma)

 required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Source code in aerokit/aero/ShockWave.py 
70
71
72
73
74
75
76
77
78
79
80
def downstream_Mn(Mn, gamma=defg._gamma):
    """

    Args:
      Mn: param gamma:  (Default value = defg._gamma)
      gamma:  (Default value = defg._gamma)

    Returns:

    """
    return np.sqrt((1.0 + 0.5 * (gamma - 1.0) * Mn ** 2) / (gamma * Mn ** 2 - 0.5 * (gamma - 1.0)))

sigma_DevMax(Mach, gamma=defg._gamma)

computes the shock angle at maximum deviation (always subsonic downstream flow), separation of weak/strong shock

Parameters:

Name Type Description Default
Mach

param gamma: (Default value = defg._gamma)

 required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Source code in aerokit/aero/ShockWave.py 
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
def sigma_DevMax(Mach, gamma=defg._gamma):
    """computes the shock angle at maximum deviation (always subsonic downstream flow), separation of weak/strong shock

    Args:
      Mach: param gamma:  (Default value = defg._gamma)
      gamma:  (Default value = defg._gamma)

    Returns:

    """
    fogpu = 4.0 / (gamma + 1.0)
    M2 = np.square(Mach)
    # ka = (M2-1.)*(1.+.5*(gamma-1.)*M2)
    # kb = .25*((gamma+1.)*M2-(3.-gamma))*M2 + 1.
    # return degree.atan(np.sqrt((kb+np.sqrt(kb**2+ka))/ka))
    return degree.asin(
        np.sqrt(
            1.0
            / fogpu
            / gamma
            / M2
            * (M2 - fogpu + np.sqrt(M2 * M2 + 2 * (gamma - 1.0) * fogpu * M2 + 4.0 * fogpu))
        )
    )

sigma_Mach_deflection(Mach, deflection, init=None, gamma=defg._gamma)

computes the shock angle from upstream Mach number and deflection angle

Parameters:

Name Type Description Default
Mach

param deflection:

 required
init

Default value = None)

 None
gamma

Default value = defg._gamma)

 _gamma
deflection
required

Returns:

Source code in aerokit/aero/ShockWave.py 
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
def sigma_Mach_deflection(Mach, deflection, init=None, gamma=defg._gamma):
    """computes the shock angle from upstream Mach number and deflection angle

    Args:
      Mach: param deflection:
      init: Default value = None)
      gamma: Default value = defg._gamma)
      deflection:

    Returns:

    """
    if init is None:
        sig0 = degree.asin(1.0 / Mach) + deflection
    else:
        sig0 = init

    def local_f(sig):
        """

        Args:
          sig:

        Returns:

        """
        return deflection_Mach_sigma(Mach, sig, gamma)

    return ITS.secant_solve(local_f, np.abs(deflection), sig0) * np.sign(deflection)

sigma_Sonic(Mach, gamma=defg._gamma)

computes the shock angle for a downstream SONIC Mach number

Parameters:

Name Type Description Default
Mach

param gamma: (Default value = defg._gamma)

 required
gamma

(Default value = defg._gamma)

 _gamma

Returns:

Source code in aerokit/aero/ShockWave.py 
308
309
310
311
312
313
314
315
316
317
318
319
320
321
def sigma_Sonic(Mach, gamma=defg._gamma):
    """computes the shock angle for a downstream SONIC Mach number

    Args:
      Mach: param gamma:  (Default value = defg._gamma)
      gamma:  (Default value = defg._gamma)

    Returns:

    """
    M2 = np.square(Mach)
    ka = gamma - 3 + M2 * (gamma + 1)
    kb = (gamma + 1) * (np.square(M2 - 3) + gamma * np.square(M2 + 1))
    return degree.asin(np.sqrt((ka + np.sqrt(kb)) / (4 * gamma * M2)))

strongsigma_Mach_deflection(Mach, deflection, gamma=defg._gamma)

Parameters:

Name Type Description Default
Mach

param deflection:

 required
gamma

Default value = defg._gamma)

 _gamma
deflection
required

Returns:

Source code in aerokit/aero/ShockWave.py 
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
def strongsigma_Mach_deflection(Mach, deflection, gamma=defg._gamma):
    """

    Args:
      Mach: param deflection:
      gamma: Default value = defg._gamma)
      deflection:

    Returns:

    """
    ka = (1.0 + 0.5 * (gamma + 1.0) * Mach ** 2) * degree.tan(deflection)
    kb = 1.0 - Mach ** 2
    kc = (1.0 + 0.5 * (gamma - 1.0) * Mach ** 2) * degree.tan(deflection)
    kd = (ka ** 2 / 3.0 - kb) / 3.0
    ke = 2.0 * ka ** 3 / 27.0 - ka * kb / 3.0 + kc
    if ke ** 2 - 4.0 * kd ** 3 > 0:
        print("no strong shock wave solution")
        return 90.0
    else:
        phi = np.arccos(-0.5 * ke / np.sqrt(kd ** 3)) + 4 * math.pi
        kf = 2.0 * np.sqrt(kd) * np.cos(phi / 3.0) - ka / 3.0
        return degree.atan(1.0 / kf)

weaksigma_Mach_deflection(Mach, deflection, gamma=defg._gamma)

Parameters:

Name Type Description Default
Mach

param deflection:

 required
gamma

Default value = defg._gamma)

 _gamma
deflection
required

Returns:

Source code in aerokit/aero/ShockWave.py 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
def weaksigma_Mach_deflection(Mach, deflection, gamma=defg._gamma):
    """

    Args:
      Mach: param deflection:
      gamma: Default value = defg._gamma)
      deflection:

    Returns:

    """
    ka = (1.0 + 0.5 * (gamma + 1.0) * Mach ** 2) * degree.tan(deflection)
    kb = 1.0 - Mach ** 2
    kc = (1.0 + 0.5 * (gamma - 1.0) * Mach ** 2) * degree.tan(deflection)
    kd = (ka ** 2 / 3.0 - kb) / 3.0
    ke = 2.0 * ka ** 3 / 27.0 - ka * kb / 3.0 + kc
    if ke ** 2 - 4.0 * kd ** 3 > 0:
        print("no weak shock wave solution")
        return degree.asin(1.0 / Mach)
    else:
        phi = np.arccos(-0.5 * ke / np.sqrt(kd ** 3))
        kf = 2.0 * np.sqrt(kd) * np.cos(phi / 3.0) - ka / 3.0
        return degree.atan(1.0 / kf)