Difficulty with starting a subsonic nozzle simulation

Cases
1

Hello,
I was able to setup PyFR with OpenMP on my HPC and run one of the test cases. I’m now trying to simulate a simple axisymmetric subsonic nozzle meshed with hexahedra (degree Q3). I’m having difficulty getting the case to run for more than ~3e-6 seconds, even when running with a very small time step of 1e-9 seconds. I believe the mesh is fine, it has run well in a similar high-order Flux Reconstruction solver. As this is the first case I’ve run with PyFR other than the provided test cases, it’s likely I’m overlooking something simple.

I tried to simplify my original approach for initialization based on recommendations I found here in the forum; such as order=1 rather than 3, scheme=euler rather than rk45, and testing out different approaches for initializing the flowfield, but have not had any success. I have not found any clues from looking at the flowfield prior to crash, just a pressure wave coming from the nozzle inflow.

When running with the euler scheme, it runs for ~3e-6 seconds and crashes with NaNs. When running with rk45, I get an error: “minimum sized time step rejected”.

Does anyone have any ideas for initializing this case?

# Set Point 3
# NPR = 1.197, NTR = 1.0

[backend]
precision = single
rank-allocator = linear

[constants]
gamma = 1.4
mu = 0.0000185
Pr = 0.71
# Custom
NPR = 1.197
NTR = 1.0
Ps = 101325
Ts = 288
Rhos = 1.225
Cp = 1006

[solver]
system = navier-stokes
order = 1
#order = 3
#viscosity-correction = sutherland
viscosity-correction = none
#shock-capturing = entropy-filter

[solver-time-integrator]
formulation = std
scheme = euler
controller = none
#controller = pi
#atol = 1e-5
#rtol= 1e-5
tstart = 0
tend = 1.00
dt = 1e-09

#[solver-entropy-filter]
#d-min = 1e-6
#p-min = 1e-6
#e-tol = 1e-6

[solver-interfaces]
riemann-solver = rusanov
ldg-beta = 0.5
ldg-tau = 0.1

[solver-interfaces-quad]
flux-pts = gauss-legendre-lobatto
quad-deg = 10
quad-pts = gauss-legendre-lobatto

[solver-elements-hex]
flux-pts = gauss-legendre-lobatto
soln-pts = gauss-legendre-lobatto
quad-deg = 10

[soln-ics]
rho = Rhos
u = 3.0
v = 0.0
w = 0.0
p = Ps

[soln-bcs-extinflow]
type = char-riem-inv
rho = Rhos
u = 3.0
v = 0.0
w = 0.0
p = Ps

[soln-bcs-farfield]
type = char-riem-inv
rho = Rhos
u = 3.0
v = 0.0
w = 0.0
p = Ps

[soln-bcs-nozzleextwall]
type = no-slp-adia-wall
#type = slp-adia-wall

[soln-bcs-nozzleinflow]
type = sub-in-ftpttang
pt = NPR*Ps
cpTt = NTR*Ts*Cp
theta = 0.0
phi = 0.0

[soln-bcs-nozzleintwall]
type = no-slp-adia-wall
#type = slp-adia-wall

[soln-bcs-nozzlelipwall]
type = no-slp-adia-wall
#type = slp-adia-wall

[soln-bcs-outflow]
type = sub-out-fp
p = Ps

[soln-plugin-writer]
dt-out = 0.0001
basedir = .
basename = soln-{t:.4f}

[soln-plugin-nancheck]
nsteps = 10

[soln-plugin-residual]
nsteps = 10
file = residual.csv
header = true

2

It looks as if your farfield conditions are trying to impose zero velocity. This seems suspect.

Regards, Freddie.

3

You are right, that does not make sense. That was a relict of trying various things with the boundary conditions to see if I could extend the simulation time.
I edited the .ini to make it more reasonable - using a sub-out-fp boundary condition on the outflow. I also added some coflow, which can sometimes aid in convergence but doesn’t seem to help here. These boundary conditions still result in a quick crash.

4

I was able to get the case to run much longer by simply switching riemann-solver from “rusanov” to “hllc”. However, that uncovered an issue that is probably the cause of my earlier difficulties getting the simulation to run for more than a small fraction of a second: my nozzle is an inlet! It is ingesting flow rather than expelling.

Am I somehow using the sub-in-ftpttang BC incorrectly? I’m initializing the flowfield with ambient pressure, and setting ambient as the downstream and farfield boundary condition; while setting the nozzle Pt to 1.2 * ambient - I don’t understand how it could be ingesting flow.

[constants]
gamma = 1.4
mu = 0.0000185
Pr = 0.71
# Custom
NPR = 1.197
NTR = 1.0
Ps = 101325
Ts = 288
Rhos = 1.225
Cp = 1006

[soln-ics]
rho = Rhos
u = 3.0
v = 0.0
w = 0.0
p = Ps

# Nozzle Inflow
[soln-bcs-nozzleinflow]
type = sub-in-ftpttang
pt = NPR*Ps
cpTt = NTR*Ts*Cp
theta = 0.0
phi = 0.0

# External inflow
[soln-bcs-extinflow]
type = char-riem-inv
rho = Rhos
u = 3.0
v = 0.0
w = 0.0
p = Ps

[soln-bcs-farfield]
type = char-riem-inv
rho = Rhos
u = 3.0
v = 0.0
w = 0.0
p = Ps

# Downstream Outflow
[soln-bcs-outflow]
type = sub-out-fp
p = Ps

5

Do you have the angles correct for the orientation of the flow you want?

6

Also, double check the units you are using to define your total pressure and enthalpy. If you check the sub-in-ftpttang.mako kernel, you’ll see that your values will give you inlet velocities on the order of 170 whereas your ICs/outlet BCs have a velocity of 3. If this is what you want (it does seem consistent with a subsonic nozzle Mach number), then you may need to start this with P0 and let the flow develop and stabilize before restarting with a higher approximation order.

7

I believe my units are correct, I’ve been sticking with SI for simplicity.

I’m starting with order=1 to let the flow develop, but still no luck - each simulation crashes with NaNs once the nozzle starts expelling flow. Example below with a time step just before crash. The velocities are much smaller in magnitude than I’d expect, I’d expect an exit velocity of around 166 m/s, but maybe this is just due to the initial start-up transients.

arn2
arn21165×1036 25.7 KB

Here’s my attempt at an .ini file for the initial flow development. Other than order=1 and scheme=euler, are there any other settings that may be more robust for startup?

# Set Point 3
# NPR = 1.197, NTR = 1.0

[backend]
precision = single
rank-allocator = linear

[constants]
gamma = 1.4
mu = 0.0000185
Pr = 0.71
# Custom
NPR = 1.197
NTR = 1.0
Ps = 101325
Ts = 288
Rhos = 1.225
Cp = 1006

[solver]
system = navier-stokes
order = 1
#order = 3
#viscosity-correction = sutherland
viscosity-correction = none
#shock-capturing = entropy-filter

[solver-time-integrator]
formulation = std
scheme = euler
#scheme = rk4
#scheme = rk45
controller = none
#controller = pi
#atol = 1e-5
#rtol= 1e-5
tstart = 0
tend = 1.00
dt = 1e-10

#[solver-entropy-filter]
#d-min = 1e-6
#p-min = 1e-6
#e-tol = 1e-6

[solver-interfaces]
#riemann-solver = rusanov
riemann-solver = hllc
ldg-beta = 0.5
ldg-tau = 0.1

[solver-interfaces-quad]
flux-pts = gauss-legendre-lobatto
soln-pts = gauss-legendre-lobatto
#flux-pts = gauss-legendre
quad-deg = 10
#quad-pts = gauss-legendre
#quad-pts = gauss-legendre-lobatto

[solver-elements-hex]
#flux-pts = gauss-legendre
#soln-pts = gauss-legendre
flux-pts = gauss-legendre-lobatto
soln-pts = gauss-legendre-lobatto
quad-deg = 10

[soln-ics]
rho = Rhos
u = 3.0
v = 0.0
w = 0.0
p = Ps

[soln-bcs-extinflow]
type = char-riem-inv
rho = Rhos
u = 3.0
v = 0.0
w = 0.0
p = Ps

[soln-bcs-farfield]
type = char-riem-inv
rho = Rhos
u = 3.0
v = 0.0
w = 0.0
p = Ps

[soln-bcs-nozzleextwall]
#type = no-slp-adia-wall
type = slp-adia-wall

[soln-bcs-nozzleinflow]
type = sub-in-ftpttang
pt = NPR*Ps
cpTt = NTR*Ts*Cp
theta = 0.0
phi = 90.0

[soln-bcs-nozzleintwall]
#type = no-slp-adia-wall
type = slp-adia-wall

[soln-bcs-nozzlelipwall]
#type = no-slp-adia-wall
type = slp-adia-wall

[soln-bcs-outflow]
type = sub-out-fp
p = Ps

[soln-plugin-writer]
dt-out = 0.0001
basedir = .
basename = soln-{t:.4f}

[soln-plugin-nancheck]
nsteps = 10

[soln-plugin-residual]
nsteps = 10
file = residual.csv
header = true

8

Order here refers to polynomial order, so I would recommend using order = 0 until the flow has developed (you’ll have to change your solution/flux points to gauss-legendre for that). Feel free to use adaptive time stepping there. Since this reduces to a first-order Godunov scheme, if it fails that indicates that something else in the problem setup is the issue (BCs, Riemann solver, time step, etc.).

9

When running with order=0 and adaptive time stepping, I get an immediate crash: “Minimum sized time step rejected”. I experimented with a few different schemes, toggling adaptive time stepping off and on, and smaller time steps, and all resulted in an immediate crash with either the above error, NaNs detected, or no error message at all.
Interestingly, when running with “riemann-solver = rusanov” I received an error of “Negative mesh jacobians detected”. Does this make sense? If there were negative Jacobians, I would have expected that error the first time I tried running the mesh - and it’s strange that I only receive the error when running order=0, riemann-solver=rusanov.

Pointwise says there are no negative Jacobians, and this mesh has been running fine in a similar FR solver - so I have assumed there is nothing wrong with the mesh. But perhaps something is going wrong on export?

10

That’s an odd error. Can you share your mesh (in GMSH format)?

11

Sure - can I send it to your @tamu.edu email provided on the PyFR Team page?

12

That’s fine with me.