The aim of this test case is to validate the following parameters of compressible steadystate turbulent flow through a de Laval Nozzle subsonic and supersonic flow regimes:
The rhoSimpleFoam solver is used for the subsonic case, while for the supersonic case the sonicFoam solver is employed. The k−ω SST
$k\omega \text{}SST$model was used to model turbulence. Simulation results of SimScale were compared to analytical results obtained from methods elucidated in [1]. The mesh was created locally using the blockMesh tool and then imported on to the SimScale platform.
Import validation project into workspace
A typical configuration of the de Laval nozzle with a nonsmooth throat was chosen as the geometry (see Table 1 for coordinates). Since the nozzle is axisymmetric, it was modeled as an angular slice of the complete geometry, with a wedge angle of 18 degrees (see Fig.1.).
Point  x
$x$

y
$y$

z
$z$


A  0  0  0 
B  0  5.5434  35 
C  0  5.5434  35 
D  68.68  0  0 
E  68.68  3.1677  20 
F  68.68  3.1677  20 
G  238.8  6.3354  40 
H  238.8  6.3354  40 
A nonuniformlyspaced hexahedral mesh was generated using the blockMesh tool (see Fig.2.). Flow near the nozzle wall was resolved using inflation of y+=30
${y}^{+}=30$for the subsonic and 300
$300$for the supersonic case. In order to keep the flow twodimensional, the mesh was designed to have only one layer in the y
$y$direction.
Tool Type : OPENFOAM®
Analysis Type : Compressible Steadystate (Turbulent)
Mesh and Element types :
Mesh type  Cells in x  Cells in y  Cells in z  Number of nodes  Type 

blockMesh  150  1  100  15000  2D hex 
Fluid:
Table 3 encapsulates the properties of fluids used in the subsonic and supersonic case simulations. The need for using a different fluid for supersonic case arises from the courant number restriction.
Case  m
$m$
g/mol $g/mol$

cp
${c}_{p}$
J/kgK $J/kgK$

mu
$mu$
N/ms $N/ms$

Pr
$Pr$


Subsonic  28.9
$28.9$

1005
$1005$

1.79×10−5
$1.79\times {10}^{5}$

1
$1$

Supersonic  11640.3
$11640.3$

2.5
$2.5$

1.8×10−5
$1.8\times {10}^{5}$

1
$1$

Boundary Conditions:
Parameter  Inlet (ABC)  Outlet (GHI)  Wall (AEHIFC)  Wedges (AGHEBA + AGIFCA) 

Velocity  7.58 ms−1
$7.58\text{}m{s}^{1}$

Zero Gradient  0.0 ms−1
$0.0\text{}m{s}^{1}$

Wedge 
Pressure  Zero Gradient  1.3×105 Nm−2
$1.3\times {10}^{5}\text{}N{m}^{2}$

Zero Gradient  Wedge 
Temperature  300 K
$300\text{}K$

Zero Gradient  Zero Gradient  Wedge 
k
$k$

0.862 m2/s2
$0.862\text{}{m}^{2}/{s}^{2}$

Zero Gradient  Wall Function  Wedge 
ω
$\omega $

484.269 s−1
$484.269\text{}{s}^{1}$

Zero Gradient  Wall Function  Wedge 
αt
${\alpha}_{t}$

0
$0$
(Calculated) 
0
$0$
(Calculated) 
Wall Function  Wedge 
μt
${\mu}_{t}$

0
$0$
(Calculated) 
0
$0$
(Calculated) 
Wall Function  Wedge 
Parameter  Inlet (ABC)  Outlet (GHI)  Wall (AEHIFC)  Wedges (AGHEBA + AGIFCA) 

Velocity  Zero Gradient  Zero Gradient  0.0 ms−1
$0.0\text{}m{s}^{1}$

Wedge 
Pressure  1.2999 Nm−2
$1.2999\text{}N{m}^{2}$

Wave Transmissive 0.0296 Nm−2
$0.0296\text{}N{m}^{2}$

Zero Gradient  Wedge 
Temperature  1.0388 K
$1.0388\text{}K$

Zero Gradient  Zero Gradient  Wedge 
k
$k$

3.75×10−5 m2/s2
$3.75\times {10}^{5}\text{}{m}^{2}/{s}^{2}$

Zero Gradient  Wall Function  Wedge 
ω
$\omega $

0.144 s−1
$0.144\text{}{s}^{1}$

Zero Gradient  Wall Function  Wedge 
αt
${\alpha}_{t}$

0
$0$
(Calculated) 
0
$0$
(Calculated) 
Wall Function  Wedge 
μt
${\mu}_{t}$

0
$0$
(Calculated) 
0
$0$
(Calculated) 
Wall Function  Wedge 
Results for the subsonic cases are calculated from Bernoulli’s equation and the ideal gas equation as follows:
Subsonic flow does not see a significant temperature rise. So the speed of sound remains almost constant and can be calculated as:
c=γRTm−−−−−√
$$c=\sqrt{\frac{\gamma RT}{m}}$$
Here, γ
$\gamma $, T
$T$, m
$m$and R
$R$represent the specific heat ratio, temperature, and molecular weight of the fluid, and the universal gas constant respectively. The Mach number can then be calculated as:
M=AinUinAc
$$M=\frac{{A}_{in}{U}_{in}}{Ac}$$
Assuming a stagnation pressure of Pstag=0.1301 MPa
${P}_{stag}=0.1301\text{}MPa$, the static pressure can be computed as:
Pstat=Pstag−12ρu2
$${P}_{stat}={P}_{stag}\frac{1}{2}\rho {u}^{2}$$
The relation between nozzle crosssection area A
$A$and Mach number M
$M$is what governs flow characteristics in supersonic flow through a de Laval nozzle [1]:
AAt=1M[2γ+1(1+γ−12M2)]γ+12(γ−1)
$$\frac{A}{{A}_{t}}=\frac{1}{M}{\left[\frac{2}{\gamma +1}\left(1+\frac{\gamma 1}{2}{M}^{2}\right)\right]}^{\frac{\gamma +1}{2(\gamma 1)}}$$
Using the known area ratio, the mach number variation is calculated by solving the above equation. The pressure can be calculated using [1]:
p=p0(1+γ−12M2)γ1−γ
$$p={p}_{0}{\left(1+\frac{\gamma 1}{2}{M}^{2}\right)}^{\frac{\gamma}{1\gamma}}$$
A comparison of the Mach number and pressure variation in the nozzle obtained with SimScale with analytical results is given in Fig.3A and 3B for subsonic and Fig.4A and 4B for supersonic flow.
Fig.3. Visualization of Mach number and pressure (A, B) along the nozzle for subsonic flow
Fig.4. Visualization of Mach number and pressure (A, B) along the nozzle for supersonic flow
The deviation from analytical results exists because the latter is calculated with a onedimensional hypothesis. Thus, all parameters are assumed to be uniform in the radial direction. Fig.5. shows that this is in fact not the case – there exists radial variation in flow variables. This is one reason why some deviation is seen between the two.
Fig.5. Contours of velocity and temperature in the nozzle. Clearly, there exists radial variation in these parameters.
[1]  (1, 2, 3) H. W. Leipmann and A. Roshko: Elements of Gasdynamics. *Courier Corporation (1957)* 
This offering is not approved or endorsed by OpenCFD Limited, producer and distributor of the OpenFOAM software and owner of the OPENFOAM® and OpenCFD® trade marks. OPENFOAM® is a registered trade mark of OpenCFD Limited, producer and distributor of the OpenFOAM software.