COMSOL 6.4 - Supersonic Air-to-Air Ejector

This application models compressible turbulent gas flow in a supersonic air ejector using the High Mach Number Flow interface in COMSOL Multiphysics. Ejectors are simple mechanical components used to induce a secondary flow by momentum and energy transfer from a high-velocity primary jet. The high-energy fluid (primary flow) passes through a convergent-divergent nozzle and reaches supersonic conditions. After exiting the nozzle, it interacts with the secondary flow which is accelerated through an entrainment-induced effect. The mixing between both flows takes place along a constant-area duct called the mixing chamber where complex interactions between the mixing layer and shocks can be observed. A diffuser is usually placed before the outlet to recover pressure and bring the flow back to stagnation.

Ejectors are used for a wide range of applications, including industrial refrigeration, vacuum generation, gas recirculation, and thrust augmentation in aircraft propulsion systems. Great efforts have been made to determine their optimum design and operating conditions, as well as how to describe the flow within them (Ref. 1).

This application models an ejector working with air in both the primary and secondary streams. The geometry and boundary conditions are based on Ref. 2, and Ref. 3. The items of interest are the primary and secondary mass flows, the static pressure distribution along the centerline of the ejector, and the resolution of the flow in the mixing region.

Model Definition

Figure 1 shows the geometry of the ejector. Its dimensions can be found in Table 1. A two-dimensional axisymmetric geometry is used to approximate the 3D geometry of the device and to reduce the size of the problem.

The flow velocity in the ejector is large enough to introduce significant variations in the density and temperature of the fluid, and the flow is governed by the fully compressible Navier–Stokes equations. Moreover, the Mach number is expected to be larger than one in the divergent section of the primary nozzle, as well as in the mixing chamber. Interaction between the boundary layers and mixing layers cause the deceleration from supersonic to subsonic flow to take place through a complex succession of shocks called shock train or pseudo-shock wave phenomenon (Ref. 4). Thus, the mesh has to be fine enough to accurately capture this phenomenon.

Figure 1: The geometry of the ejector.

Table 1: Dimensions of the ejector in mm.
d1	d2	dt	dnd	dm	dd
16	160	8	12	24	51
L1	L2	L3	Lm	Ld	NXP
7	23	90	240	70	15

The problem is modeled using the Favre-averaged Navier–Stokes equations and the standard k-ε turbulence model.

Both primary and secondary flows are air with a specific gas constant of 287 J/(kg⋅K) and a ratio of specific heats of 1.4. The dynamic viscosity and thermal conductivity of the air are computed from Sutherland’s law.

Boundary Conditions

Inlet

The flow at the inlets is specified in terms of its total properties: T0 = 300 K and P0 = 5 atm for the primary flow, and T0 = 300 K and P0 = 0.55 atm for the secondary flow.

The inlet conditions are applied using an Inlet feature, where the Flow condition is specified to be . This provides a numerically consistent boundary condition by evaluating the flow characteristics at the inlet.

The velocities at both inlets are unknown. However, they are expected to be very small compared to the velocities inside the nozzle and mixing chamber. The values that must be prescribed at the inlet are the total values of temperature and pressure, which define the energy of the flow. The Mach number can be set to 0 and will be determined by the characteristics based boundary condition at the inlet. This provides a good initial solution, but the total values of pressure and temperature may differ slightly. The solution can be improved if the problem is solved again setting the Mach number to the values computed by the characteristics based boundary condition at the inlets, which are approximately 0.14 and 0.01 for the primary and secondary inlets, respectively.

The inlet values for the turbulent kinetic energy, k, and the turbulent dissipation rate, ε, are approximated from the turbulent intensity, IT, and turbulence length scale, LT. The turbulent intensity is set by default to 0.05 (5%). The length scale can be approximated as 7% of the pipe diameter or hydraulic diameter. This is done automatically when is set to .

For more background on this boundary condition, see the CFD Module User’s Guide.

Outlet

The flow reaches supersonic conditions inside the ejector. However, it is expanded and decelerated along the mixing chamber and the diffuser, reaching subsonic conditions before being discharged into the atmosphere. The outlet is then subsonic with atmospheric static pressure. This is modeled using an Outlet node with the set to subsonic.

Results and Discussion

The mass flows obtained are depicted in Table 2. The distributions of Mach number and velocity inside the ejector are depicted in Figure 2 and Figure 3. The primary flow is accelerated in the convergent section of the nozzle, reaching sonic conditions at the throat, and is expanded further in the divergent section. At the outlet of the primary nozzle, the secondary flow acts as an artificial wall for the primary flow. This leads to the formation of virtual nozzle throats, and a succession of expansion and compression waves can be observed in the region upstream of the mixing zone. Then, the flow decelerates along the constant-area duct and is brought back to stagnation in the diffuser. The region where both flows mix can be visualized by plotting the turbulent kinetic energy, see Figure 4.

Figure 5 plots the distribution of pressure along the centerline and walls of the mixing chamber. At the centerline of the duct, the flow successively changes from supersonic to subsonic flow via multiple shocks. However, this cannot be detected by wall pressure measurements because the surface pressures tend to be smeared out due to the dissipation in the boundary layer (see Ref. 4). The distribution of temperature is shown in Figure 6. Very low temperatures can be observed inside the device. This must be taken into account when designing an ejector, specially when working with two-phase flows. The results obtained correlate well with Ref. 2 and Ref. 3.

Table 2: Mass flows.

0.057 kg/s	0.037 kg/s	0.094 kg/s

Figure 2: Mach number distribution.

Figure 3: The distribution of velocity in the nozzle and the mixing chamber.

Figure 4: The distribution of turbulent kinetic energy. The mixing region is clearly identified.

Figure 5: The distribution of pressure along the mixing chamber.

Figure 6: The distribution of temperature. Very low temperatures can be observed at the outlet of the nozzle. The temperature at the outlet of the ejector is lower than at both inlets.

Notes About the COMSOL Implementation

The present application is highly nonlinear and sensitive to the solution procedure. The mesh needed to capture the interaction between the shocks and the mixing layer, and to resolve the solution near the walls, is extremely fine. However, convergence may be hard to achieve with such a fine mesh unless a good enough initial solution is used. A way to overcome this is to first solve the problem on a coarse mesh and then refine it. The solution on the coarse mesh provides good initial values, but lacks accuracy in three important areas: wall resolution, capture of shocks, and resolution of the mixing layer. The adaptive mesh refinement feature can be used to overcome this. However, in order to fully resolve the mesh, a high element growth rate would be needed, potentially leading to convergence problems. An alternative option is to first refine manually both on the boundary layers and in the nozzle, and then to use the adaptive mesh refinement feature to resolve the mixing layer.

References

1. S. He, Y. Li, and R.Z. Wang, “Progress of Mathematical Modeling on Ejectors,” Renew. Sustain. Energy Rev., vol. 13, pp 1760–1780, 2009.

2. Y. Bartosiewicz, Zine Aidoun, P. Desevaux, and Yves Mercadier, “Numerical and Experimental Investigations on Supersonic Ejectors,” Int. J. of Heat and Fluid Flow, vol. 26, pp 56–70, 2005.

3. P. Desevaux,A. Bouhangel, and E. Gavignet, “Flow Visualization in Supersonic Ejectors Using Laser Tomography Techniques,” Int. J. of Refrigeration, vol. 34, pp 1633–1640, 2010.

4. F. Gnani, H. Zare-Behtash, and K. Kontis, “Pseudo-shock Waves and their Interactions in High-speed Intakes,” Progress in Aerospace Sciences, vol. 82, pp 36–56, 2016.

CFD_Module/High_Mach_Number_Flow/supersonic_ejector

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > > ε (hmnf).

Click .

Click

In the tree, select > .

Click

Geometry 1

The model geometry is available as a parameterized geometry sequence in a separate MPH file.

In the toolbar, click and choose .

Browse to the model’s Application Libraries folder and double-click the file supersonic_ejector_geom_sequence.mph.

In the toolbar, click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
P1	5[atm]	5.0663E5 Pa	Total pressure, primary flow
P2	0.55[atm]	55729 Pa	Total pressure, secondary flow
Pout	1[atm]	1.0133E5 Pa	Pressure, outlet
T1	300[K]	300 K	Total temperature, primary flow
T2	T1	300 K	Total temperature, secondary flow
Rs	287[J/kg/K]	287 J/(kg·K)	Specific gas constant
gamma	1.41	1.41	Ratio of specific heats
iso_diff	0	0	Isotropic diffusion

Add isotropic diffusion to improve convergence when computing the initial solution with a coarse mesh.

Click the

button in the toolbar.

In the dialog, in the tree, select the checkbox for the node > .

Click .

High Mach Number Flow, k-ε (hmnf)

In the window, under click ε.

In the window for ε, click to expand the section.

Find the subsection. Select the checkbox.

In the δid text field, type iso_diff.

Find the subsection. Select the checkbox.

In the δid text field, type iso_diff.

Fluid 1

In the window, under > ε (hmnf) click .

In the window for , locate the section.

From the Rs list, choose . In the associated text field, type Rs.

From the γ list, choose .

From the γ list, choose . In the associated text field, type gamma.

Initial Values 1

In the window, click .

In the window for , locate the section.

Specify the u vector as


0	r
100	z

In the p text field, type 2[atm].

Inlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the p0,tot text field, type P1.

In the T0,tot text field, type T1.

In the Ma0 text field, type 0.14.

Inlet 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the p0,tot text field, type P2.

In the T0,tot text field, type T2.

In the Ma0 text field, type 0.01.

Outlet 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the p0 text field, type Pout.

The next step is to generate a coarse mesh to compute an initial solution.

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Click

Component 1 (comp1)

Generate a second mesh. Refine the boundary layer mesh and increase the mesh resolution on the walls and in the nozzle.

Mesh 2

In the toolbar, click and choose .

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Boundaries 4 and 5 only (select the walls of the nozzle).

Locate the section. From the list, choose .

From the list, choose .

Size 2

Right-click and choose .

In the window for , locate the section.

From the list, choose .

In the list, choose and .

Click

Select Boundaries 6–9 and 11–15 only.

Locate the section. From the list, choose .

From the list, choose .

Size 3

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domain 2 only.

Locate the section. From the list, choose .

From the list, choose .

Size 4

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domain 3 only.

Locate the section. From the list, choose .

From the list, choose .

Size

In the window, click .

In the window for , locate the section.

From the list, choose .

Corner Refinement 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click to select the

toggle button.

From the list, choose .

Free Triangular 1

In the toolbar, click

Boundary Layers 1

In the toolbar, click

Boundary Layer Properties

In the window, click .

In the window for , locate the section.

From the list, choose .

In the list, choose , , and .

Click

Select Boundaries 6, 7, 9, and 11–15 only.

Locate the section. In the text field, type 5.

From the list, choose .

In the text field, type 5e-5.

Boundary Layer Properties 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 10.

From the list, choose .

In the text field, type 1e-5.

Select Boundaries 4, 5, and 8 only (select the walls of both the primary inlet and the nozzle).

Click

Study 1

Step 1: Stationary

In the window, under click .

In the window for , click to expand the section.

In the table, enter the following settings:


Component	Mesh
Geometry 1	Mesh 1

Edit the study step to solve the model in two steps in order to increase the robustness of the solution procedure. First, solve for P2=Pout (no adverse gradient of pressure), and then solve the problem again adding the adverse gradient of pressure (P2 smaller than Pout). Disable the continuation solver since it is mainly suitable for linear applications.

Click to expand the section. Select the checkbox.

Click