Reversed Flow in sub-out-fp Boundary Conditions

Hi! I am hoping to simulate a converging nozzle expelling from atmospheric conditions into an approximately 47kPa chamber. To do this, I’ve set up a sub-in-ftpttang inlet and sub-out-fp outlet condition. My simulation runs for long enough that the flow develops, but always seems to unexpectedly diverge near the outlet (outlets are the right, top, and bottom walls of the outlet chamber). My theory is that reversed flow near the outlet is causing this divergence, and I wanted to ask whether:

  • This is indeed a possibility
  • There may be a way to prevent reversed flow. In my case, I don’t believe it would be problematic to impose zero velocity across the domain should reverse flow occur.

Alternatively, would you recommend I impose a sup-in-fa boundary condition with some velocity? (this doesn’t feel entirely physical).

I have also tried extending my domain and could definitely do this further but it feels there should be a more elegant solution.

I am attaching my .ini file and an image of my solution near divergence.



[backend]
precision = double

[backend-openmp]
cc=gcc
cflags = -O3
schedule = guided, 6

[backend-cuda]
mpi-type = cuda-aware

[constants]
; Dyn Viscosity kg m^-1 s^-1
mu = 0.0000181 
; Gamma dimensionless
gamma = 1.4
; Prandtl number dimensionless
Pr = 0.72
; Static temperature K
Ts = 273
; Isobaric Specific heat J Kg^-1 K^-1
Cp = 1005
; Density kg m^-3
rhoIn = 1.225
; Pressure Pa
Pin = 101325
Pout = 47015
; Gas Constant R J Kg^-1 K^-1
R = 287.058

[solver]
system = navier-stokes
order = 2
shock-capturing = entropy-filter

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

[solver-time-integrator]
scheme = rk45
tstart = 0.0
tend = 0.002
dt = 6e-9
controller = none

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

[solver-interfaces-line]
flux-pts = gauss-legendre

[solver-elements-tri]
soln-pts = williams-shunn

[solver-elements-quad]
soln-pts = gauss-legendre

[soln-plugin-nancheck]
nsteps = 10

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

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

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

[soln-bcs-wall-slip]
type = slp-adia-wall

[soln-bcs-inlet]
type = sub-in-ftpttang

pt = Pin
cpTt = Ts*Cp
theta = 0.0
phi = 0.0

[soln-bcs-outlet]
type = sub-out-fp
; Set up to ramp pressure in 0.0005s. Maybe could be done more smoothly?
;p = Pout + (Pin+11000 - Pout) * (1 - 1/(1 + exp(-4000*(t-0.00040))))
;p = max(Pout, (Pin - (2000*t)*(Pin-Pout)))
p = Pout

[soln-ics]
u = 0.0
v = 0.0
;p = min(Pin,Pin + (Pin-Pout)*(max(-100*x,-1)))
p = min(Pin,Pin + (Pin-Pout)*(max(-200*(x-0.0),-1)))
;p = Pin
;rho = (min(Pin,Pin + (Pin-Pout)*(max(-100*x,-1))))/(R*Ts)
rho = (min(Pin,Pin + (Pin-Pout)*(max(-200*(x-0.0),-1))))/(R*Ts)

Cheers,
Tom

Have you tried switching to the char-riem-inv condition?

Regards, Freddie.

Hi Freddie,

Thanks for the quick reply! I was reluctant to try the char-riem-inv BC as it required a specific velocity (as opposed to just u>0). As I’m studying viscous choking I don’t know the exact mass flow through the nozzle (which would depend on BL thickness) or the velocity distribution downstream.

Seeing your comment I decided to give it a shot with an inviscid guess for u, and this appears to have stabilized the solution while still giving reasonable results. Thanks a lot for the help :slight_smile:

Best wishes,
Tom