Dependency of Mach number? (same Reynolds number)

Cases
1

Hi everyone,

I am running a 3D simulation of flow around an axisymmetric slender body using PyFR.

The geometry has a sharp pointed nose and overall length-to-diameter ratio of (L/D = 4.0).

The baseline case is:

  • P3
  • Mach 0.1
  • Re_D=200,000
  • Angle of attack = 50 degrees

This case runs without any issues. However, when I increase the Mach number to 0.4, while keeping the Reynolds number, angle of attack, polynomial order, and all other conditions the same, the simulation diverges.

I also tried the following cases:

  • P2, Mach 0.4, AoA = 50 degrees: stable
  • P3, Mach 0.4, AoA = 40 degrees: stable

I also attempted to gradually increase the polynomial order from P1 to P2 and then to P3, but the P3 case at 50 degrees still diverges.

I know that the divergence might be caused by the mesh resolution or quality.

But since the Reynolds number is fixed, I expected the mesh requirements would not change significantly when just increasing the Mach number 0.1 to 0.4.

Of course, the Mach number reaches the transonic or supersonic regime, additional stabilization would be consider, However, in this case, Mach 0.4 is still subsonic.

Does anyone have any suggestions on what I should check?

I have attached the input .ini files for reference.

Best regards, uk

[backend]
precision = double
rank-allocator = linear

[backend-cuda]
device-id = local-rank

[constants]
gamma = 1.4
Pr = 0.72
rhoc = 1.0
M = 0.4
ac = 1.0
uc = %(M)s
Re = 200000.0
mu = %(rhoc)s*%(uc)s/%(Re)s
pi = 3.141592653589793
pinf = %(rhoc)s*%(ac)s*%(ac)s/(%(gamma)s*%(M)s*%(M)s)
qinf = 0.5*%(rhoc)s*%(uc)s*%(uc)s
aoa = 50.0
phi = 0.0
tconv = 20.0
tout_conv = 1.0

[solver]
system = navier-stokes
order = 3
anti-alias = flux, surf-flux

[solver-time-integrator]
scheme = rk45
controller = pi
t0 = 0.0
dt = 1e-6
atol = 1e-5
rtol = 1e-5
tend = 33.33

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

[solver-interfaces-tri]
flux-pts = williams-shunn
quad-deg = 7
quad-pts = williams-shunn

[solver-elements-tet]
soln-pts = shunn-ham
quad-deg = 6
quad-pts = shunn-ham

[soln-bcs-nose]
type = no-slp-adia-wall

[soln-bcs-body]
type = no-slp-adia-wall

[soln-bcs-base]
type = no-slp-adia-wall

[soln-bcs-far]
type = char-riem-inv
rho = rhoc
u = uc*cos(aoa/180*pi)
v = -uc*sin(aoa/180*pi)*sin(phi/180*pi)
w = uc*sin(aoa/180*pi)*cos(phi/180*pi)
p = pinf

[soln-ics]
rho = rhoc
u = uc*cos(aoa/180*pi)
v = -uc*sin(aoa/180*pi)*sin(phi/180*pi)
w = uc*sin(aoa/180*pi)*cos(phi/180*pi)
p = pinf

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

[soln-plugin-dtstats]
flushsteps = 100
file = dtstats.csv
header = true

[soln-plugin-fluidforce-nose]
nsteps = 10
file = nose-forces.csv
header = true 

[soln-plugin-fluidforce-body]
nsteps = 10
file = body-forces.csv
header = true 

[soln-plugin-fluidforce-base]
nsteps = 10
file = base-forces.csv
header = true 

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

[soln-plugin-tavg]
nsteps = 2000
dt-out = 0.1
basedir = .
basename = tavg-{t:.2f}

avg-p = p
avg-u = u
avg-v = v
avg-w = w

avg-uu = u*u
avg-vv = v*v
avg-ww = w*w
avg-uv = u*v
avg-vw = v*w
avg-wu = w*u

fun-avg-Cp = (p - pinf)/qinf
2

Have you looked at what is going on with the flow just before it diverges in one of the higher Mach cases?

Regards, Freddie.

3

Hi Freddie,

Thanks for the suggestion.

Yes, I looked at the flow field, but only briefly so far, since the simulations were run on an isolated server and I have limited access to the data.

In the Mach 0.4 case, I first ran a warm-up simulation at P1 with 1/10 of the Re_D (20,000), and then restarted from that solution using P3 with the standard Re_D (200,000).

The P3 simulation diverged after just about 3.0 time units (D/U_infty, D: diameter).

I then checked the latest output before divergence (dt-out was set to 0.1).

Based on my previous experience, non-physical cells tend to appear first around the base-flow region, so I took a closer look at that region.

However, I did not find any clear indication of negative (or very large) pressure/density there either.

Best regards, uk

4

So it is important to get into a position where you can easily render up the solution as this is by far the best way to determine at an instantaneous level what is causing divergence. With that said, I would try restarting at p = 2 with an intermediate Reynolds number.

Regards, Freddie.

5

Thank you.

Once I can access them again, I will take a look at the instantaneous flow field just before the divergence occurs.

Trying (p = 2) with an intermediate Reynolds number also seems like a useful idea. I will try that, as well as carrying out some other checks, such as comparing the P2 and P3 fields.

Best regards, uk