Transonic Flow in a Sajben Diffuser

In the present example the high speed turbulent gas flow in a converging and diverging nozzle is modeled using the High Mach Number Flow interface. The diffuser is transonic in the sense that the flow at the inlet is subsonic, but due to the contraction and the low outlet pressure, the flow accelerates and becomes sonic (Ma = 1) in the throat of the nozzle. After a short region of supersonic flow, a normal shock wave brings the flow back to subsonic flow. This setup has been studied in a number of experiments and numerical simulations by M. Sajben and co-workers (see for example Ref. 1, Ref. 2, Ref. 3, and Ref. 4), in an effort to study unsteady fluctuations in supersonic inlets with applications in supersonic aircraft propulsion systems. The geometry and setup is often referred to as a Sajben diffuser and constitutes a common test case for the simulation of high Mach number internal flows. In this example, the time-averaged transonic flow through a Sajben diffuser is solved for using two different exit pressures. The flow in the diffuser is fully turbulent with a inlet Reynolds number of 7 × 105 based on the inlet fluid properties and the channel height. The model uses the Spalart-Allmaras turbulence model to compute the turbulent viscosity. For the first outlet pressure value a normal shock is present, but the flow remains attached throughout the diverging part. For the second, lower, outlet pressure value, the shock is strong enough to cause a shock-induced separation in the diverging part. Based on the ability to induce flow separation, the shock in the first case is referred to as weak, while in the second case it is termed strong following the definition in Ref. 2.

Model Definition

Figure 1 shows the physical geometry of the converging and diverging nozzle model. The nozzle dimensions correspond to those used in the experiments in Ref. 2 and in the benchmarks simulation of Ref. 5 and Ref. 6. In the central contraction part, the minimum vertical height separating the lower and upper walls, the throat height hth, is 1.7322 in (44 mm). The channel height at the inlet is 1.4hth, and the outlet height is 1.5hth.

Figure 1: Geometry and dimensions of the Sajben diffuser model.

Fluid Properties

The fluid occupying the channel is assumed to be air, by specifying a specific gas constant of 287 J/(kg·K) and a ratio of specific heats of 1.4. The dynamic viscosity is computed using Sutherland’s Law:

where the Tref = 500 R (about 278 K) corresponds to the inlet total temperature, and the Sutherland constant is Sμ = 111 K. The reference viscosity μref is defined from the inlet Reynolds number:

which in turn is used in a parametric sweep, where the model is solved for increasing Reynolds number. The final inlet Reynolds number to be computed is 7 × 105, corresponding to the one used in the simulation in Ref. 5. The thermal conductivity of the gas is defined using the definition of the Prandtl number

which is assumed to be 0.71 in this case. This is a typical number for air around 293 K.

Boundary Conditions

Inlet Condition

The flow at the inlet is subsonic with a flow speed corresponding to a Mach number of 0.46. The inlet conditions are specified in terms of total properties, where the total pressure is defined as 19.58 psi and the total temperature is 500 R. The inlet conditions are applied using an Inlet feature, where the Flow condition is specified to be Characteristics based. This provides a numerically consistent boundary condition by evaluating the flow characteristics at the inlet (for more background on this boundary condition, see the CFD Module User’s Guide).

Outlet Condition

At the outlet the static pressure is specified. The model is solved for the two outlet pressure values specified in Table 1 below.

Table 1: Outlet Pressure.
Pressure	Fraction of inlet total pressure	description
16.05 psi	0.82	Case 1: weak shock
14.10 psi	0.72	Case 2: strong shock

These outlet pressure values are known from experiments and simulations to be low enough to produce sonic conditions at the throat of the nozzle. However, they are not low enough for the flow to stay supersonic throughout the diverging part. The supersonic flow in the divergent part is terminated by a normal shock wave, so that the flow in the following part including the outlet becomes subsonic. The pressure is specified in the model using an Outlet node with the Flow condition set to Subsonic.

Results and Discussion

Below, some of the results from the transonic diffuser model computed in COMSOL Multiphysics are shown and discussed. The results are compared to experimental data from Ref. 2 (strong shock case) and Ref. 4 (weak shock case). The experimental data was extracted as tabulated data from Ref. 5 and Ref. 6 and plotted in COMSOL using interpolation functions.

Figure 2 shows the Mach numbers and velocity streamlines resulting from applying the first outlet pressure, 16.05 psi. It can be seen that the flow accelerates in the converging part, reaches sonic conditions at the throat, after which a region of supersonic flow follows in the diverging part. The supersonic region is terminated by a normal shock wave, which brings the flow back to subsonic conditions. In the remaining part of the channel the flow decelerates subsonically toward the outlet. The zero contour of the x-component velocity is also plotted in the figure, but this is only present on the walls and not visible inside the domain. Hence no separation zone is present, and the flow remains attached throughout the divergent part of the channel. The shock is not able to cause flow separation and is therefore termed weak. These results correlate well with those in Ref. 2 and Ref. 5.

Figure 2: Mach number, flow streamlines, and zero x-component velocity contour resulting from the weak shock case.

Figure 3 shows the same quantities as Figure 2, but uses the results from the second outlet pressure case, pout = 14.10 psi. Due to the lower outlet pressure, the normal shock wave is positioned further downstream in the divergent channel part. More importantly, a flow separation zone can be seen behind the shock, as indicated by the zero x-component velocity contour. The shock wave in this case is apparently strong enough to induce flow separation. This result is in accordance with those presented in Ref. 2 and Ref. 6.

Figure 3: Mach number, flow streamlines, and zero x-component velocity contour (in red) resulting from the strong shock case.

Figure 4 shows the development of the static pressure on the upper wall normalized by the inlet total pressure. Results from both the weak and strong cases are plotted and compared with the experimental data of Ref. 1 and Ref. 2. The results from both outlet cases are in general in very good agreement with the experimental results in the diffuser. Note however that the shock positions in the model are slightly shifted in the downstream direction in comparison with the experiments.

For analysis of the results in the interior of the channel, Figure 5 plots the streamwise velocity profiles from the strong shock case at two different positions in the divergent part of the channel together with experimental results. The velocity profile at the first position, x = 4.611hth, compares very well with the experimental results. Both the velocity magnitude and the size of the separation zone, including reversed flow, are accurately reproduced. Further downstream, at the x = 6.340hth position, the velocity magnitude in the central part of the channel is also in good agreement with the experimental results. Closer to the upper wall, the model results include flow reversal at this position. This is not found in the experimental result, indicating that the separation zone in the model extends further downstream than that in the experiment.

Figure 4: Top wall static pressure normalized by the inlet total pressure. Model results (lines) and experimental results (diamonds) are shown.

Figure 5: Mean x-component velocity at two positions downstream of the strong shock. Model results (lines) and experimental results (diamonds) are shown.

Notes About the COMSOL Implementation

The present model is highly nonlinear and sensitive to the solution procedure. The sensitivity is accentuated by the fact that the model seeks to determine the equilibrium position of a normal shock wave, positioned in a channel with smoothly varying channel height.

You solve the model in two steps. First, apply the higher outlet pressure to simulate the weak shock case. To solve this case, use a parametric sweep where the inlet Reynolds number increases stepwise by decreasing the dynamic viscosity. When you have obtained a converged result for the highest Reynolds number, use this solution as the initial condition for the second case with the lower outlet pressure value. In both cases, use pseudo-time stepping with manually defined CFL number expressions to compute stationary solutions.

References

1. M. Sajben, J.C. Kroutil, and C.P. Chen, “A High-Speed Schlieren Investigation of Diffuser Flows with Dynamic Distortion”, AIAA Paper 77-875, 1977.

2. T.J. Bogar, M. Sajben, and J.C. Kroutil, “Characteristic Frequencies of Transonic Diffuser Flow Oscillations,” AIAA Journal, vol. 21, no. 9, pp. 1232–1240, 1983.

3. J.T. Salmon, T.J. Bogar, and M. Sajben, “Laser Doppler Velocimetry in Unsteady, Separated, Transonic Flow,” AIAA Journal, vol. 21, no. 12, pp. 1690–1697, 1983.

4. T. Hsieh, A.B. Wardlaw Jr., T.J. Bogar, P. Collins, and T. Coakley, “Numerical Investigation of Unsteady Inlet Flowfields,” AIAA Journal, vol. 25, no. 1, pp. 75–81, 1987.

5. NPARC Alliance Validation Archive, “Sajben Transonic Diffuser: Study #1,”, 2008, http://www.grc.nasa.gov/WWW/wind/valid/transdif/transdif01/transdif01.html

6. NPARC Alliance Validation Archive, “Sajben Transonic Diffuser: Study #2,”, 2008, http://www.grc.nasa.gov/WWW/wind/valid/transdif/transdif02/transdif02.html

CFD_Module/High_Mach_Number_Flow/sajben_diffuser

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

Geometry parameters

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
x0	-6.99809[in]	-0.17775 m	Inlet x-position
xEnd	14.98353[in]	0.38058 m	Outlet x-position
h_in	2.44483[in]	0.062099 m	Diffuser inlet height
h_out	2.59830[in]	0.065997 m	Diffuser outlet height
h_th	1.732[in]	0.043993 m	Throat height

Right-click and choose .

In the dialog box, type Geometry parameters in the text field.

Click .

Fluid flow parameters

In the toolbar, click

and choose .

In the window for , type Fluid flow parameters in the text field.

Locate the section. In the table, enter the following settings:


Name	Expression	Value	Description
Rein	7e5	7E5	Inlet Reynolds number
case	1	1	Case number; 1 = weak shock, 2 = strong shock
Min	0.46	0.46	Inlet Mach number
gamma	1.4	1.4	Ratio of specific heats
Pr	0.72	0.72	Prandtl number
Rs	287[J/kg/K]	287 J/(kg·K)	Specific gas constant
Tin_tot	500[R]	277.78 K	Inlet total temperature
Tin_stat	Tin_tot/(1+0.5Min^2(-1+gamma))	266.5 K	Inlet static temperature
pin_tot	19.58[psi]	1.35E5 Pa	Inlet total pressure
pin_stat	pin_tot/(1+0.5Min^2(-1+gamma))^(gamma/(-1+gamma))	1.1677E5 Pa	Inlet static pressure
rhoin	pin_stat/Rs/Tin_stat	1.5267 kg/m³	Inlet density
mu_ref	rhoinu_inh_in/Rein	2.0387E-5 kg/(m·s)	Reference dynamic viscosity
u_in	Minsqrt(gammaRs*Tin_stat+eps)	150.53 m/s	Inlet velocity
pOut	if(case==1,16.05,0)[psi]+if(case==2,14.1,0)[psi]	1.1066E5 Pa	Outlet pressure

Interpolation 1 (int1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Find the subsection. In the table, enter the following settings:


Function name	Position in file
top_pos	1

Click .

Browse to the model’s Application Libraries folder and double-click the file sajben_diffuser_upper_wall.txt.

Click .

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

Interpolation 2 (int2)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Find the subsection. In the table, enter the following settings:


Function name	Position in file
ptop_weak	1

Click .

Browse to the model’s Application Libraries folder and double-click the file sajben_diffuser_ptop_weak.txt.

Click .

Interpolation 3 (int3)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Find the subsection. In the table, enter the following settings:


Function name	Position in file
ptop_strong	1

Click .

Browse to the model’s Application Libraries folder and double-click the file sajben_diffuser_ptop_strong.txt.

Click .

Interpolation 4 (int4)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Find the subsection. In the table, enter the following settings:


Function name	Position in file
u_at4611	1

Click .

Browse to the model’s Application Libraries folder and double-click the file sajben_diffuser_u-xh4611.txt.

Click .

Interpolation 5 (int5)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Find the subsection. In the table, enter the following settings:


Function name	Position in file
u_at6340	1

Click .

Browse to the model’s Application Libraries folder and double-click the file sajben_diffuser_u-xh6340.txt.

Click .

Definitions

Variables 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
CFLnum	if(case==1,CFLweak,0)+if(case==2,CFLstrong,0)		CFL number for pseudo time stepping
CFLweak	1.3^min(niterCMP-1,9)+if(niterCMP>25,5*1.2^min(niterCMP-26,12),0)		CFL number, weak case
CFLstrong	1+if(niterCMP>10,1.2^min(niterCMP-10,12),0)+if(niterCMP>120,1.3^min(niterCMP-120,9),0)+if(niterCMP>220,1.3^min(niterCMP-220,9),0)		CFL number, strong case

The manual CFL number expression for the strong shock corresponds to the implemented automatic expression for turbulent flows. In this case the solution already contains a shock that will move due to the change in the outlet pressure, and a cautious increase of the CFL number is needed. In the weak case simulation, a shock is not yet formed, and the simulation time can be reduced by using a more aggressive ramping of the CFL number.

Geometry 1

Parametric Curve 1 (pc1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type x0[1/in].

In the text field, type xEnd[1/in].

Locate the section. In the text field, type s[in].

In the text field, type top_pos(s).

Click

Polygon 1 (pol1)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the text field, type x0 x0 x0 xEnd xEnd xEnd.

In the text field, type h_in 0 0 0 0 h_out.

Click

Convert to Solid 1 (csol1)

In the toolbar, click

and choose .

Click in the window and then press Ctrl+A to select both objects.

In the window for , click

Add a rectangular domain in the divergent part of the nozzle. This will be used to increase the resolution in the region where the shock is located.

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 0.16.

In the text field, type 0.1.

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

Click

Partition Objects 1 (par1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Find the subsection. Select the

toggle button.

Select the object only.

Click

Add a feature to specify the interior boundaries as mesh control entities. In this manner these entities can be used to control the mesh, but at the same time they will automatically be omitted when defining the physics and when postprocessing results.

Mesh Control Edges 1 (mce1)

In the toolbar, click

and choose .

On the object , select Boundaries 3 and 5 only.

It might be easier to select the boundaries by using the window. To open this window, in the toolbar click and choose . (If you are running the cross-platform desktop, you find in the main menu).