# Simulation of a jet of air into a shock tube

I am simulating the resulting turbulent flow generated from injecting air at a supply (absolute) pressure of 202.65 kPa through an inlet of diameter 0.5", into a shock tube that has an initial (absolute) pressure of 101.325 kPa. The dimensions of the rectangular shock tube's cross section are 2.5"x1.5". I am only modelling a portion of the length of the shock tube, shown below. When I hit Calculate, the residuals start off below 1, and then the x, y, z equations shoot up to the order of 10's and 100's. I hope somebody can help me troubleshoot what is supposed to be a simple turbulence simulation.

My flow is compressible so I selected the density-based solver. Transient flow is selected as well. I am using the k-omega model to model the turbulence, and my solution method is implicit. I set 3 boundary conditions:

1) inlet: pressure-inlet --> Gauge pressure total pressure is set to 101325 Pa, and supersonic/initial gauge pressure is set to 0 (i also tried different pressures here)

2) left and right planes indicated in the diagram: both set as pressure-outlet --> gauge pressure is set to 0 as that is the gauge pressure in my shock tube at those two locations

Are the boundary conditions correctly set? Is anything missing?

In the solution initialization, I am computing from inlet, where the Y-velocity is automatically set to -329 m/s2.

The below attached images are of the shock tube, the fluid domain (air) that I am modelling, the mesh, and the set boundary conditions in FLUENT.

## Comments

148how's time step selected?

17I tried:

10 time steps with step size of 0.1 seconds

20 time steps with step size 0.05 seconds

148why these time steps? what is the flow convective CFL (hint: it should be ~ 1)

17My Courant number is set at 5. Why should it be 1? I'll try 1.

I chose those time steps because I want to simulate 1 second of injection - as a start at least.

148Convective CFL has nothing to do with the CFL defined in the solution controls. This CFL controls dual time stepping. Convective CFL ~ U dt/dx. What it is in you case? Basically, you should not drive the flow faster than 1 cell per time step, especially if it has sharp gradients.

17oh. so i'll keep the Courant number in the solution controls as 5, and i'll set the convective CFL to <=1. where can i find the convective CFL though?

17oh you mean i should choose a time step (dt) that will give me a CFL = 1 according to the equation U dt/dx. I get it. For the velocity it initializes from (322 m/s), the required time step that gives me a CFL of 1 is 4.86e-5 s. So in order to simulate 1 second of flow, I need 20,573 time steps? That's huge.

148for reporting convective CFL, you can temporarily switch to the pressure-based solver which has convective CFL under velocity. No need to run the solution in pressure-based to inspect the CFL.

Requirement to maintain convective CFL ~ 1 (and acoustic CFL if you are interested in acoustics) is pretty much universal across CFD solvers. The implicit density-based formulation allows you to run with larger CFLs but at a cost of loosing transient information due to numerical dissipation and potentially arriving to a wrong answer. Very large CFLs can also lead to stability problems, as you are observing.

But as long as you understand the trade off between larger time step and accuracy/fidelity of the solution, you ca try larger steps.

17okay thank you! im running the simulation explicitly with a courant number set at 1. im trying to add a plot in the monitor for courant number in order to monitor it and make sure it doesn't go above 1, but i can't find it anywhere when creating a report definition. I read somewhere that im supposed to find it under "velocity" when creating a volume monitor. Any idea how to add that to my monitor? I'm using ANSYS Fluent 2020.

148As per my previous note, convective CFL is available as a post-processing quantity in the pressure-based solver, but not in the density-based, You can construct an expression or custom field function using available post-processing quantities: velocity magnitude and (cell volume)^(1/3) and monitor it.

17Hello. So I ran this simulation using a pressure based solver, since my flow is subsonic and is not THAT compressible. I first ran it steady-state for 1000 iterations and used the solution to initialize the transient sim. The simulation ran very smoothly with no problems at all, the residuals retained a consistent typical zig-zag form through, and was converging with every time-step. Results looked logical, net mass flow rate was almost 0 through all boundaries, etc.

My settings for the steady-state sim are: pressure-based solver, coupled scheme with pseudo transient activated and a timescale factor of 0.1 (that was the only one that got my solution to start converging).

My settings for the transient model are: Pressure-based solver, PISO scheme, Adaptive time-stepping, CFL-based, Courant number = 0.5. Those settings ran perfectly for the simulation done before I added the pipe inlet.

I then added a pipe to the inlet because the previous simulation did not account for the sudden expansion that occurs when the flow leaves the pipe into the shock tube. I also ran the steady state solution first, but the residuals were slightly different than the first sim, in that the continuity residual did not decrease with the others, but was nevertheless stable:

I used the SS solution to initialize the transient sim, the zig-zag form varied slightly throughout, the solution was converging for every time-step, until it started diverging - then crashed. From the moment i started the transient simulation, the net mass flow rate kept decreasing below 0 until it crashed. Below are snap shots of the residuals before the solution began to diverge, and after it diverged.

Here is part of the error message I got:

Reversed flow on 329 faces (99.7% area) of pressure-inlet 6.

Stabilizing k to enhance linear solver robustness.

Negative k in 54992 cells after linear solve.

Stabilizing temperature to enhance linear solver robustness.

Stabilizing temperature using GMRES to enhance linear solver robustness.

Divergence detected in AMG solver: temperature

absolute pressure limited to 5.000000e+10 in 39988 cells on zone 4

turbulent viscosity limited to viscosity ratio of 1.000000e+05 in 479327 cells

Divergence detected in AMG solver: temperature

Reversed flow on 3176 faces (83.4% area) of pressure-outlet 7.

Reversed flow on 124 faces (1.7% area) of pressure-outlet 8.

Divergence detected in AMG solver: temperature

Divergence detected in AMG solver: temperature

Error at host: floating point exception

Error at Node 1: floating point exception

Error at Node 0: floating point exception

Error at Node 2: floating point exception

Error at Node 3: floating point exception

===============Message from the Cortex Process================================

Compute processes interrupted. Processing can be resumed.

==============================================================================

Error: floating point exception

Error Object: #f

------------------------------------------------------------------------------------------------------

Do you think the mesh is the problem? That's the only thing that changed since I added the pipe to the inlet. Here are my mesh metrics, and some snaps of my mesh.

orthogonal quality:

max: 4.2318e-002

min: 0.99817

average: 0.72103

standard deviation: 0.186

skewness:

max: 4.5386e-005

min: 0.95741

average: 0.27731

standard deviation: 0.18758

aspect ratio:

max: 1.1601

min: 147.36

average: 4.4613

standard deviation: 3.0285

element quality:

max: 1.4197e-003

min: 0.99998

average: 0.5242

standard deviation: 0.25719

2,326Hello,

Firstly, this is not a good mesh to run your simulations on. Your mesh quality is poor and needs improvement.

Secondly, what is your inlet condition?

Also, what is the Mach Number at the nozzle exit? I ask this question because you are talking about the propagation of shock.

What material properties are you using? Specifically, density.

Karthik

148To add to Karthik's comments, I am sure you did due diligence and checked the steady-state simulation to confirm the flow developed as expected. Mesh min ortho quality is low but not too awful. Before investing into redoing the mesh, can you confirm what is happening before the divergence? E. g. save contours of the flow field (velocity, pressure, Mach) on the center plane through the domain (YZ plane) every few steps to identify what happens when the divergence occurs. This will give some pointers as to where the problem maybe coming from. Also, monitor your time step to see if it changes suddenly at the time of the crash.

17Hey guys. Thank you for your replies.

1) Yeah I'm not happy with the cylindrical pipe's mesh at all, but I've been having difficulty generating a good one. The cells on the outer surface look very skewed to me, but if my mesh metrics are within the required ranges as specified by the tutorials, then what exactly defines my mesh as good enough or bad enough? What guidelines do you recommend I follow in order to generate a good enough mesh for this simulation? Is the rectangular shock tube's mesh bad too or is it just the pipe?

2) My boundary conditions are:

inlet is set as pressure-inlet --> Gauge total pressure is set to 101325 Pa, and supersonic/initial gauge pressure is set to 0 (i'm going to set this to 5731 Pa gauge because that is the pressure when the flow reaches M=1, is that what this parameter is defined as?)

left and right planes are both set as pressure-outlet --> gauge pressure is set to 0 as that is the gauge pressure in my shock tube at those two locations

3) My maximum Mach number at the pipe exit according to the generated contours at around 94,700 iterations of the transient sim is around 1.08, so there should be shocks in my flow. i miscalculated that earlier. should i switch to a density-based solver in this case? if so, implicit or explicit?

4) I'm using air as a compressible ideal gas, with density 1.225 kg/m3 at standard conditions.

5) I'm re-running the simulation now to capture contours at different stages as advised. I'll have them uploaded by tomorrow, as well as the variation in time steps. When do I know when to stop increasing the Courant number? At 0.5, the time step remains constant at 1e-8, and at 0.55, it also remains constant at 1e-5. But when I set it at 0.6, the time step starts to increase. Is this a sign that I should leave it at 0.55?

17Below are contours of Mach Number, Velocity Magnitude, and Static Pressure t=0.002079 seconds (before it crashed), and at t=0.0021082 seconds (after it crashed). I also monitored the inlet pressure, the mass flow rate, and the time step.

1) The time step remained constant at 1e-8 for the entire time. It did not change even after it crashed. I had the courant number set to 0.55. I tried increasing it to 0.6 and 0.7 at times and the time step started increasing, then i returned it to 0.55. But at a courant of 0.5 and 0.55, it was constant throughout at 1e-8.

2) The average weighted inlet pressure fluctuated around 6.816-6.82 kPa, before it shot straight up to huge numbers and fluctuating after 0.0021077 seconds, until it crashed. I posted screenshots of the plots of both phases.

3) The mass flow rate was steadily decreasing below 0.0001 before it crashed. Plots of both phases are shown below.

4) Mach Numbers before and after

5) Velocity before and after

5) Static Pressure before and after

148since the flow becomes mildly supersonic and the mesh quality is poor, I'd recommend the pressure-based coupled solver (use Coupled option in p-v coupling). It can hold together well into the supersonic regime, no need to switch to the dbns. If there are still problems, the revisiting the mesh would be a good idea.

17okay i'll do that. are there certain guidelines you recommend i follow in order to achieve a good quality mesh? and is it the pipe only that's bad quality? or is it the shock tube as well? i have a bias near the shock tube walls and near the vicinity of the jet.

148I would sweep-mesh the pipe so it would have good quality hex mesh, and then tet mesh the rest. It's hard to say just by looking at the surface mesh where bad-quality cells are, but definitely the pipe has highly skewed elements. So define the sweep method for the pipe volume, mesh it fist and then mesh the shock tube.

17So i fixed the mesh using your suggested method, and it looks much better. Below are the mesh metrics and screenshots of the mesh:

orthogonal quality:

max: 0.10314

min: 0.9997

average: 0.74944

standard deviation: 0.14951

skewness:

max: 4.4457e-009

min: 0.89686

average: 0.2491

standard deviation: 0.15097

elemental quality:

max: 9.3866e-002

min: 0.99997

average: 0.58812

standard deviation: 0.2668

Aspect ratio:

max: 1.1595

min: 16.521

average: 3.9131

standard deviation: 2.4327

Here are screenshots of the residuals and contours for the steady state solution, they seem logical and fine to me:

And here are the residuals and contours for the transient solution. I used coupled scheme here, and i also tried the piso but i wasnt able to get proper residuals. the contours are very much similar to the steady state because (i assume) i used the steady state solution to initialize my transient one, and the time step that was generated using the adaptive method was 1e-8, hence there wasnt much advancement in the flow time since initialization. but the problem is the residuals, and the net mass flow rate at the inlet deviated heavily from 0.

148actually, looks pretty good. Residuals in transient have some different meaning then in steady-state. All you should be looking is residuals go done ~ 1e-03 within each step. Keep in mind that once you switch to transient, you start recovering transient artifacts (propagating waves) which otherwise don't exist in steady-state. Since this is a closed domain, waves inevitable reflect off the walls, maybe travel upstream along the tube. They should eventually dissipate, My suggestion is to switch to a constant time step, use a larger step at flow CFL ~ 1, and keep coupled. This should allow you to flash out initial transient faster.