Computational Fluid Dynamics - Advanced (a)

de Laval Nozzle (Converging-Diverging Nozzle)

The converging-diverging nozzle is perhaps the most important piece of engineering hardware associated with propulsion and the high speed flow of gases. It was invented by a Swedish engineer Gustaf de Laval in 1888 and is thus often referred to as the 'de Laval' nozzle. This principle was first used in a rocket engine by Robert Goddard. Very nearly all modern rocket engines that employ hot gas combustion use de Laval nozzles.

The de Laval nozzle operation relies on the different properties of gases flowing at subsonic and supersonic speeds. The gas flow through a de Laval nozzle is isentropic.

In this section results of numerical simulations (viscous and inviscid cases) by means of Fluent and OpenFOAM for two different geometries are compared with the 1-D model results. The first geometry is a converging-diverging nozzle with a sharp kink at the throat and in the second case, this sharp kink is replaced by a curved edge (Fig. a.1). All the comparisons are done for two cases: subsonic and supersonic. In this way, the results of two different solvers are compared, the effects of smoothening of the geometry is shown, and the validity of the 1-D model is examined.

Fig. a.1: Geometries; sharp-shaped (up) and curve-shaped (down)

PART I: 1-D Model

In this section, a 1-D model for a straight axis nozzle will be presented. The results of such a model has to be interpreted as mean-flow results for each cross-section of the nozzle and will be used to compare with CFD simulations of Fluent and OpenFOAM.

Because of the straightness of the axis, we can neglect mean acceleration effect in radial velocity. Furthermore, in order to simplify the problem (but also keeping its supposed main features), an isentropic flow is assumed all over the nozzle. This assumption can be punctually denied by the presence of a shock-wave without making the 1-D model invalid. A further assumption for a reasonable 1-D model is that the changing of section has to be smooth enough to avoid importance of local fluid dynamics effects (like boundary layer separations). According to it, one should expect a certain discrepancy between the 1-D model and the sharp-shaped nozzle, at least close to the kink.

For a steady-state, inviscid ideal gas flow in an adiabatic nozzle, it holds:
\begin{equation} \begin{cases} \frac{A}{A_t} = \frac{1}{M} \left[\frac{1}{1+\gamma}\left(1+\frac{\gamma-1}{2}M^2\right) \right]^{\frac{\gamma+1}{2(\gamma-1)}} \\ \frac{p}{p_0} = \frac{1}{1+\gamma}\left(1+\frac{\gamma-1}{2}M^2\right)^{-\frac{\gamma}{\gamma-1}} \\ \frac{T}{T_0} = \frac{1}{1+\gamma}\left(1+\frac{\gamma-1}{2}M^2\right)^{-1} \end{cases} \end{equation} where $ A $ is a general cross-section of the nozzle, $ A_t $ is the throat cross-section, $ M $ is the average Mach number in $ A $, $ \gamma = 1.4 $, $ p $ is the average static pressure in $ A $, $ p_0 $ is the average total pressure, $ T $ is the average static temperature in $ A $, and $ T_0 $ is the average total temperature.

Specifying the geometry of the nozzle, it is possible to obtain an estimation of two critical Mach numbers such that a shock-wave occurs in the interval between these two critical values. The cases to investigate in the present section do not consider the simulation of a shock-wave. The 1-D model will be then used to define boundary conditions for Fluent and OpenFOAM such that the simulated flows places outside the critical Mach number range.

The chosen geometry for the nozzle considers a throat cross-section diameter equal to $ 0.02 $ m and an outflow one equal to $ 0.04 $ m. The critical Mach number range at the outflow cross-section is then:
\begin{equation} M \in \left[M_{cr1}, M_{cr2}\right] = \left[0.145, 2.945\right] \end{equation} Setting the total pressure and temperature to be respectively $ 1.3 $ atm and $ 300 $ K, the average critical conditions for static pressure and temperature at the outflow cross-section are:
\begin{equation} \begin{cases} p_{cr1} = 129822 Pa \\ T_{cr1} = 298.7 K \\ p_{cr2} = 3895 Pa \\ T_{cr2} = 109.7 K \end{cases} \end{equation} Fig. a.2 shows the section ratio vs. Mach number for 1D isentropic flow in critical conditions.

Fig. a.2: $ A/A_t $ vs. Mach number in choking condition

PART II: Fluent

In this section, the created geometry and mesh using ICEM CFD are presented, then the Fluent setup will be discussed for different simulation cases.

Assuming an axisymmetric flow in the nozzle, the problem is reduced to a 2-dimensional problem as shown in Fig. a.1. Then the mesh is generated using ICEM CFD. The mesh is coarser close to the symmetry line and finer adjacent to the wall since the boundary layer should be resolved accurately in this region (Fig. a.3). After considering the mesh independency test, the number of nodes are set to 150 in axial and 100 in radial direction.

Fig. a.3: Mesh produced by ICEM CFD

After importing the mesh into Fluent, the following procedure is done to setup the software. It should be noted that the setup is explained for one of the cases (subsonic) and the differences are explained for the other case (supersonic).

Solution setup

Before setting the solver, it is necessary to check the scales of the geometry in “General > Mesh > Scale...” and check the quality of the mesh. For the subsonic case, the Pressure-based steady solver, and for the supersonic case, the Density-based steady solver is employed. Also, the solver has to be set to “Axisymmetric” mode for the 2D space.

In the “Models” menu, the energy equation is turned on and the standard $ k $-$ \epsilon $ model with enhanced wall function is selected. Moreover, the material is set to air as an ideal gas and the operating pressure is kept as $ 1 $ atm.

Boundary conditions

For the inlet boundary condition, the mass-flow-inlet with the mass flow rate of $ 0.05 $ kg/s is used (subsonic case). For the supersonic case, the pressure-inlet boundary condition is selected with the gauge total pressure of $ 30400 $ Pa. The turbulence specification method is set to “Intensity and Hydraulic Diameter” for both cases with turbulent intensity of $ 1% $ and hydraulic diameter of $ 0.07 $ m. Also, the temperature is set to $ 300 $ K.

The outlet boundary condition is set to pressure-outlet with the gauge pressure of $28675 $ Pa ($ -98325 $ Pa for the supersonic case) in order to be higher than the critical pressure (see previous section for the details). Consequently, the temperature is set to $ 298.9 $ K ($ 109.8 $ K for the supersonic case).

Then the type of the symmetry boundary condition is set to “axis” and no-slip (adiabatic) boundary condition is used for the wall.

Solution methods

For both the Pressure-based and Density-based steady solvers, the SIMPLE scheme is employed with linear discretization method and the under-relaxation factors are kept as default. The convergence criteria for all the residuals are set to $ 10^{-7} $. Then the standard initialization computing from the inlet is selected. Finally, the calculation is started by the number of iterations of $ 1000 $ with the reporting and update interval of $ 10 $.


Geometry and meshing

In this part, only the sharp-shaped geometry is investigated using OpenFOAM. To build the mesh, two ways were used: converting the ICEM CFD mesh: using the command "fluentMeshToFoam", and building the grid with OpenFOAM using the command "buildMesh". Both meshes are tested and the second one is employed for the simulations.

OpenFOAM solvers

In both subsonic and supersonic cases, the compressibility effects are supposed to play the main role in the evolution of the flow. Therefore, it is very important to use the solvers which are implemented to solve the compressible flows. Furthermore, in order to have a fair comparison with the Fluent simulations, the steady-state solvers are preferred.

Subsonic flow solver

For the subsonic case, only density-based solvers are suitable. In this section, "rhoSimpleFoam", a steady-state density-based solver is employed. This solver implements SIMPLE algorithm for the simulations.

In order to investigate the role of turbulence, "laminar", "$ k $-$ \epsilon $" and "$ k $-$ \omega $ SST" models are implemented and the results are presented in the successive section.

More details are presented in "Subsonic laminar/turbulent flow in de Laval nozzle".

Supersonic flow solver

The only steady-state solvers available for compressible simulations are "rhoPorousSimpleFoam", "rhoSimpleFoam" and "rhoSimplecFoam", but none of them is suggested for trans/supersonic flows by the UserGuide of OpenFOAM. Therefore, for the supersonic case, "sonicFoam", a transient, pressure-based solver, is used. This solver implements PIMPLE algorithm for the simulations.

Similar to the previous case, the "laminar", "$ k $-$ \epsilon $" and "$ k $-$ \omega $ SST" models are implemented.

More details are presented in "Supersonic laminar/turbulent flow in de Laval nozzle".

The codes used to set up the $ k $-$ \omega $ SST model for both cases are given in Appendix B.

"ParaView" and "foamLog"

To post-process the simulation results, mainly two tools were used: physical variables field visualization via "ParaView" using the command "paraFoam", and exporting residuals data using the commands -solver > log and foamLog log.

The outcome from "paraView" is directly used to present "OpenFOAM" results in the next section, whereas exported data from "foamLog" required a further elaboration in Matlab.

Results and discussions

Fluent results

In this section, the results of Fluent simulations for subsonic and supersonic cases for each geometry are presented. For easier comparison all the plots are collected in Appendix A.

Subsonic flow in a sharp-shaped nozzle

In Fig. A1 and A2, the viscous and inviscid models are compared and the contour results in terms of pressure, temperature and Mach number are plotted consequently. In such a case, because of the conditions which are chosen for the simulations, two different simulation results are observed for viscous and inviscid cases.

As it is clear in Fig. A1 (c) and Fig. A2 (a), the boundary layer is developed along the nozzle wall because of the fluid viscosity. In contrast, this boundary layer cannot be observed in the inviscid case. However, in the convergent part of the nozzle, the flow contours are very similar for both the models and it means that the viscosity effects do not play an important role in the convergent part of the nozzle.

The sharp kink of the nozzle causes an infinite pressure gradient at the nozzle throat which results in totally different flow behavior for viscous and inviscid models. In the viscous case, a boundary layer separation is evident which is clearly shown in Fig. A2 (c). This phenomenon results in higher output velocity along the axis in the viscous case (Fig. A2 (a)) since the effective outflow cross-section is reduced by the boundary layer separation (the dark blue region close to the wall). On the other hand, for an inviscid flow, the sharp kink produces a jump which is convected along the streamlines until the end of the nozzle as shown in Fig. A1 (d) as well as Fig. A2 (b).

As it is described, the viscosity affects the evolution of temperature and Mach number in the nozzle. So the temperature and Mach number fields calculated by the inviscid model cannot represent the real fluid behavior at least in a local sense. Nevertheless, the pressure field is well predicted by the inviscid model (Fig. A1 (a) and (b)).

The similar structure for contour fields comparison is used for all the following cases.

Supersonic flow in a sharp-shaped nozzle

In contrast with the subsonic sharp-shaped nozzle, a discontinuity is observed for the inviscid model (Fig. A3 (b) and (d), Fig. A4 (b)). The flow behavior describes a sudden Prandtl-Meyer expansion in the absence of viscosity because of the sharp kink. As shown in Fig. A3 (a) and (c) and also Fig. A4 (a), the main viscosity effect is to smooth such a discontinuity.

However, in the converging part of the nozzle, the viscosity effects are negligible and the two models result in very similar flow fields. Moreover, comparing to the subsonic boundary layer behavior, the supersonic boundary layer is much thinner and remains attached along the diverging part of the nozzle (except for a small region in the neighbor of the kink). According to the viscous model results, a good approximation would be neglecting the boundary layer effects in the whole flow field since the compressibility effects dominate.

Subsonic flow in a curve-shaped nozzle

For the curve-shaped geometry (Fig. A5 and A6), more or less the same differences can be observed between the viscous and inviscid models. In the viscous model, the considered curvature at the throat does not play an important role for a subsonic flow. In the inviscid case, because of the geometry smoothness, no convected jump is observed in the diverging part of the nozzle. This could be observed by comparing Fig. A1 (d) and Fig. A5 (d).

Furthermore, considering Fig. A5 (a) and comparing it with Fig. A1 (a), it could be seen that the curvature produces smoother flow expansion at the throat neighbor cross-sections.

Supersonic flow in a curve-shaped nozzle

Considering contour plots of Fig. A7 and A8, the only relevant effect of the curvature is observed for the inviscid model where the Prandtl-Meyer expansion is spread along the curved kink. This avoids the presence of discontinuities observed in the sharp-shaped supersonic nozzle.

Other effects of the throat curvature will be discussed in the next section where the viscous, inviscid and 1-D model results are compared along the nozzle axis.

Fluent vs. 1-D model

The main aim of this comparison is to validate the viscous model implemented in Fluent. In this section, the flow properties such as pressure, temperature and Mach number for 1-D model calculation and viscous and inviscid simulations (along the axis of the nozzle) are compared. The graphs are shown in Fig. a.4 and Fig. a.5.

Fig. a.4: Fluent vs. 1-D model, subsonic case

Fig. a.5: Fluent vs. 1-D model, supersonic case

In order to have a fair comparison between 1-D model and viscous simulation, static and total pressure and temperature are matched at the inlet. In addition, the Mach number is set at the inlet. For supersonic case these conditions are enough for fully describe the flow because just one isentropic solution is possible inside the channel independent of the outflow conditions. In contrast, for the subsonic case, an outflow condition has to be set as well and for that the pressure matching is used.

In the supersonic case, the nozzle is choked, so we can define the Mach number just in dependence on the geometry. Because the choking relation between $ A/A_t $ and Mach number (Fig. a.2) is not valid in the considered subsonic case, this relation has to be substituted by the continuity of the mass flow. More detailed calculations can be found in the Matlab scripts in Appendix B.

Considering Mach number graphs in Fig. a.4, higher velocity at the outflow can be observed for the viscous case as a result of the boundary layer separation (see also Fig. A2 (a) and Fig. A6 (a)).

Furthermore, the flow properties are compared for the curve- and sharp-shaped nozzles (supersonic case) as shown in Fig. a.5. For the inviscid model the smoothing effect of the curvature is evident.

From these comparisons it is concluded that the 1-D model results are consistent with the flow simulations.

OpenFOAM results

As described before, in OpenFOAM just the sharp-shaped nozzle is simulated. To understand how the turbulence models affect the results, for both sub- and supersonic cases, three setups ("laminar", "$ k $-$ \epsilon $" and "$ k $-$ \omega $ SST") are compared and the results are shown in Fig. A9 to A14. Pressure, temperature and velocity contours are considered for both the regimes, it is concluded that the two different turbulence models do not differ so much. Although there are some differences in temperature field between the laminar and turbulence models.

Since the wall treatment of the $ k $-$ \omega $ SST model is more consistent with the one used in the Fluent simulations, this model is selected for further comparisons with the Fluent results.

OpenFOAM vs. Fluent

In Fig. A15 to A18 pressure, temperature and velocity contours are compared between Fluent and OpenFOAM simulations for the sharp-shaped nozzle case. In the supersonic case, the results are presented in a non-dimensional form rescaling pressure with 1 atm, temperature with $ 288.15 $ K, and velocity with $ 340 $ m/s in order to have a higher efficiency in the supersonic OpenFOAM transient solver.

Apart from small quantitative differences, especially close to the walls, the results are in a good agreement for both the regimes.

As already discussed before, a fair comparison between Fluent and OpenFOAM computational times is only possible for the steady-state subsonic solvers (in supersonic, transient solver is used in OpenFOAM), and the residuals trends are shown in Fig. A19. For comparing the two simulation software, it is necessary to solve one case with both of them with the same settings and compare the computational time. As shown in Fig. A19, OpenFOAM is much more efficient than Fluent.

OpenFOAM vs. Fluent

The following items describe the most important conclusions of this section:
  1. Results of the inviscid model are not valid for the divergent part of the nozzle in the subsonic case since the boundary layer separation is not modeled.
  2. The curvature at the nozzle throat do not have a considerable effect on the flow pattern in the viscous model. But for the inviscid model, the curvature smooths the Prandtl-Meyer expansion downstream of the throat.
  3. According to the comparison along the nozzle axis, the Fluent simulations and the 1-D model results are in good agreement.
  4. Different turbulent models are not playing a dominant role in the flow field evolution.
  5. OpenFOAM is much quicker (about 3 times) than Fluent in the steady-state solvers according to the comparison for subsonic sharp-shaped case.
In this section, only two nozzle shapes are considered; but for a more comprehensive investigation, it is recommended to examine different nozzle shapes. But since one of the main aims of this section is to compare the results of Fluent and OpenFOAM considering the same geometries and conditions, it was sufficient to employ two geometries (considering the limited time frame). Further investigations can be done using different nozzle shapes in future works.

1 comment:

  1. Hello!...
    Thank you so much for this tutorial, it is really nice. I was looking to the Apendix B and I would like to know how did you define the boundary conditions in the OpenFoam Supersonic case. I see that for temperature your are using: uniform 1.038781163 (Kelvins?); and for pressure internalField uniform 1.299995065 (Pa); Can you help me understand this case setup?...