comp_flow/
oblique.rs

1//! Weak oblique shock functions
2
3use eqsolver::single_variable::FDNewton;
4use num::Float;
5
6/// Wave angle for weak oblique shock
7///
8/// # Examples
9///
10/// ```
11/// use comp_flow::oblique_beta;
12///
13/// assert_eq!(oblique_beta(2.0_f32, 1.4_f32, 0.1745329_f32), 0.6861576);
14/// assert_eq!(oblique_beta(5.0_f64, 1.4_f64, 0.3490659_f64), 0.5201241529003784);
15///
16///
17/// ```
18pub fn oblique_beta<F: Float>(mach: F, gamma: F, theta: F) -> F {
19    let beta_max: F = oblique_beta_max(mach, gamma);
20    let mut x0 = beta_max;
21    let two = F::from(2.0).unwrap();
22
23    let f = |x: F| {
24        theta.tan()
25            - two / x.tan() * (mach.powi(2) * x.sin().powi(2) - F::one())
26                / (mach.powi(2) * (gamma + (two * x).cos()) + two)
27    };
28    let beta_result = FDNewton::new(f).solve(x0);
29    let mut beta = match beta_result {
30        Ok(x) => x,
31        Err(_) => return F::nan(),
32    };
33
34    while (beta > beta_max) || (beta < F::zero()) {
35        x0 = x0 - F::from(0.1).unwrap();
36        let beta_result = FDNewton::new(f).solve(x0);
37        beta = match beta_result {
38            Ok(x) => x,
39            Err(_) => return F::nan(),
40        };
41    }
42
43    beta
44}
45
46/// Maximum oblique shock angle
47pub fn oblique_beta_max<F: Float>(mach: F, gamma: F) -> F {
48    ((F::one() / (gamma * mach.powi(2))
49        * (((gamma + F::one()) / F::from(4.0).unwrap() * mach.powi(2)) - F::one()
50            + ((gamma + F::one()) * (gamma + F::one()) / F::from(16.0).unwrap() * mach.powi(4)
51                + (gamma - F::one()) / F::from(2.0).unwrap() * mach.powi(2)
52                + F::one())
53            .sqrt()))
54    .sqrt())
55    .asin()
56}
57
58/// Mach number after weak oblique shock
59///
60/// # Examples
61///
62/// ```
63/// use comp_flow::oblique_mach2;
64///
65/// assert_eq!(oblique_mach2(2.0_f32, 1.4_f32, 0.1745329_f32), 1.6405221);
66/// assert_eq!(oblique_mach2(5.0_f64, 1.4_f64, 0.3490659_f64), 3.022151379742916);
67///
68/// ```
69pub fn oblique_mach2<F: Float>(mach: F, gamma: F, theta: F) -> F {
70    let beta = oblique_beta(mach, gamma, theta);
71    let two = F::from(2.).unwrap();
72
73    ((F::one() + (gamma - F::one()) / two * mach.powi(2))
74        / (gamma * mach.powi(2) * beta.sin().powi(2) - (gamma - F::one()) / two)
75        + (mach.powi(2) * beta.cos().powi(2))
76            / (F::one() + (gamma - F::one()) / two * mach.powi(2) * beta.sin().powi(2)))
77    .sqrt()
78}
79
80/// Stagnation pressure ratio across weak oblique shock
81///
82/// # Examples
83///
84/// ```
85/// use comp_flow::oblique_p02_p01;
86///
87/// assert_eq!(oblique_p02_p01(2.0_f32, 1.4_f32, 0.1745329_f32), 0.98464406);
88/// assert_eq!(oblique_p02_p01(5.0_f64, 1.4_f64, 0.3490659_f64), 0.5050701357775494);
89///
90/// ```
91pub fn oblique_p02_p01<F: Float>(mach: F, gamma: F, theta: F) -> F {
92    let beta = oblique_beta(mach, gamma, theta);
93    let mach1n = mach * beta.sin();
94    let two = F::from(2.).unwrap();
95
96    F::one()
97        / ((two * gamma / (gamma + F::one()) * mach1n.powi(2)
98            - (gamma - F::one()) / (gamma + F::one()))
99        .powf(F::one() / (gamma - F::one()))
100            * (two / (gamma + F::one()) / mach1n.powi(2) + (gamma - F::one()) / (gamma + F::one()))
101                .powf(gamma / (gamma - F::one())))
102}
103
104/// Static pressure ratio across weak oblique shock
105///
106/// # Examples
107///
108/// ```
109/// use comp_flow::oblique_p2_p1;
110///
111/// assert_eq!(oblique_p2_p1(2.0_f32, 1.4_f32, 0.1745329_f32), 1.7065787);
112/// assert_eq!(oblique_p2_p1(5.0_f64, 1.4_f64, 0.3490659_f64), 7.037411017501249);
113///
114/// ```
115pub fn oblique_p2_p1<F: Float>(mach: F, gamma: F, theta: F) -> F {
116    let beta = oblique_beta(mach, gamma, theta);
117    let mach1n = mach * beta.sin();
118    F::from(2.).unwrap() * gamma / (gamma + F::one()) * (mach1n.powi(2) - F::one()) + F::one()
119}
120
121/// Static density ratio across weak oblique shock
122///
123/// # Examples
124///
125/// ```
126/// use comp_flow::oblique_rho2_rho1;
127///
128/// assert_eq!(oblique_rho2_rho1(2.0_f32, 1.4_f32, 0.1745329_f32), 1.4584259);
129/// assert_eq!(oblique_rho2_rho1(5.0_f64, 1.4_f64, 0.3490659_f64), 3.3154179190165545);
130///
131/// ```
132pub fn oblique_rho2_rho1<F: Float>(mach: F, gamma: F, theta: F) -> F {
133    let beta = oblique_beta(mach, gamma, theta);
134    let mach1n = mach * beta.sin();
135    (gamma + F::one()) * mach1n.powi(2)
136        / ((gamma - F::one()) * mach1n.powi(2) + F::from(2.).unwrap())
137}
138
139/// Static temperature ratio across weak oblique shock
140///
141/// # Examples
142///
143/// ```
144/// use comp_flow::oblique_t2_t1;
145///
146/// assert_eq!(oblique_t2_t1(2.0_f32, 1.4_f32, 0.1745329_f32), 1.17015128);
147/// assert_eq!(oblique_t2_t1(5.0_f64, 1.4_f64, 0.3490659_f64), 2.122631652901466);
148///
149/// ```
150pub fn oblique_t2_t1<F: Float>(mach: F, gamma: F, theta: F) -> F {
151    let two = F::from(2.).unwrap();
152    let beta = oblique_beta(mach, gamma, theta);
153    let mach1n = mach * beta.sin();
154    (two + (gamma - F::one()) * mach1n.powi(2))
155        * (two * gamma * mach1n.powi(2) - (gamma - F::one()))
156        / ((gamma + F::one()).powi(2) * mach1n.powi(2))
157}
158
159/// Speed of sound ratio across weak oblique shock
160///
161/// # Examples
162///
163/// ```
164/// use comp_flow::oblique_a2_a1;
165///
166/// assert_eq!(oblique_a2_a1(2.0_f32, 1.4_f32, 0.1745329_f32), 1.08173530);
167/// assert_eq!(oblique_a2_a1(5.0_f64, 1.4_f64, 0.3490659_f64), 1.4569254108915342);
168///
169/// ```
170pub fn oblique_a2_a1<F: Float>(mach: F, gamma: F, theta: F) -> F {
171    let two = F::from(2.).unwrap();
172    let beta = oblique_beta(mach, gamma, theta);
173    let mach1n = mach * beta.sin();
174    ((two + (gamma - F::one()) * mach1n.powi(2))
175        * (two * gamma * mach1n.powi(2) - (gamma - F::one()))
176        / ((gamma + F::one()).powi(2) * mach1n.powi(2)))
177    .sqrt()
178}
comp_flow/oblique.rs

comp_flow/
oblique.rs
