IMPLICIT SCHEMES FOR THE SIMULATION OF UNDERWATER IMPLOSION AND EXPLOSION PROBLEMS? Charbel Farhat, Alex Main and Kevin Wang Department of Aeronautics and Astronautics Department of Mechanical Engineering Institute for Computational and Mathematical Engineering Stanford University
Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.
IMPLICIT SCHEMES FOR THE SIMULATION
OF UNDERWATER IMPLOSION
AND EXPLOSION PROBLEMS?
Charbel Farhat, Alex Main and Kevin Wang
Department of Aeronautics and Astronautics
Department of Mechanical Engineering
Institute for Computational and Mathematical Engineering
Stanford University
Stanford, CA 94305
MOTIVATION
Are explicit schemes really required
to perform such simulations?
SHOCK TUBE (GASGAS)
0.5
1.0
P = 1
r = 1
M = 0
g = 1.6
P = 0.1
r = 0.125
M = 0
g = 1.2
POTENTIAL OF IMPLICIT SCHEMES
t1 = 1e1 s
5.0e01
Initial
CFL = 1
4.0e01
CFL= 20
3.0e01
CFL = 40
Density (kg/m3)
2.0e01
1.0e01
0.0e+00
5.0e01
6.0e01
7.0e01
X (m)
OBJECTIVES 2010
 Multifluid and multiphase flow problems
 Multifluid and multiphase fluidstructure interaction problems
1
1
Fj,j+1 = Fj+1/2 (nj,j+1) = (Fj + Fj+1 )  F’ j+1/2 (Wj+1 – Wj)
= Roe (Wj, Wj+1, gs, ps) (stiffened gas)
2
2
@(rf)
(ruf)
@
+ = 0
@t
@x
COMPUTATIONAL FRAMEWORK
j + 1/2
j
j + 1
(conservation form)
FVMERS
Wjn
W*n
W*n
Wj+1n
j  1
j  1/2
j
j + 1/2
j + 1
 Fj,j+1 = Roe (Wjn, W*n, EOSj)
Fj+1,j = Roe (Wj+1n, W*n, EOSj+1)
 W*n and W*n determined from the exact solution of local
twophase Riemann problems
C. Farhat, A. Rallu and S. Shankaran, "A HigherOrder Generalized
Ghost Fluid Method for the Poor for the ThreeDimensional
TwoPhase Flow Computation of Underwater Implosions",
Journal of Computational Physics, Vol. 227, pp. 76747700 (2008)
Wnpj+1
Wnpj
1
1
(RR(pI; pR,rR) 
RL(pI; pL,rL))
uI = (uL + uR) +
2
2
RL(pI; pL,rL) +
RR(pI; pR,rR) + uR – uL = 0
pI, rIL, rIR, uI
 Newton’s method
LOCAL RIEMANN SOLVER
rIL,pI,uI ,rIR
contact discontinuity
rarefaction
shock
t
gas
water
x
j
j + 1/2
j + 1
rLuL pL
rRuR pR
FVMERS (EXPLICIT)
j  1/2
j + 1/2
j  1
j
j + 1
Fj,j+1 = Fj+1,j = Roe (Wjn, Wj+1n, EOSj= EOSj+1)
If fjnfj+1n< 0 then
Fj,j+1 = Roe (Wjn, WjRn(rIL, pI, uI), EOSj)
Fj+1,j = Roe (Wj+1n, W(j+1)Rn(rIR, pI, uI), EOSj+1)
Dt
~
 Wjn+1 = Wjn  (Fj,j+1  Fj,j1) (forward Euler)
Dx
~
 Unpack Wn+1 using fn and solve the levelset equation to get fn+1
 Pack Wpn+1 using fn+1 to get the updated solution Wn+1
FVMERS (IMPLICIT)
j  1/2
j + 1/2
j  1
j
j + 1
Fj,j+1 = Fj+1,j = Roe (Wjn+1, Wj+1n+1,EOSj= EOSj+1)
If fjnfj+1n< 0 then
Fj,j+1 = Roe (Wjn+1, WjRn+1,EOSj)
Fj+1,j = Roe (Wj+1n+1, W(j+1)Rn+1, EOSj+1)
Dt
~
 Wjn+1 = Wjn  (Fj,j+1  Fj,j1) (backward Euler)
Dx
~
 Unpack Wn+1 using fn and solve the levelset equation to get fn+1
 Pack Wpn+1 using fn+1 to get the updated solution Wn+1
IMPLICIT FLUIDFLUID
dFj,j+1 @Fj,j+1 @Fj,[email protected]+1
+
=
[email protected][email protected]+1 @Wjn+1
dFj,j+1 @Fj,[email protected]+1
=
[email protected]+1 @Wj+1n+1
dFj+1,j @Fj+1,j @Fj+1,[email protected](j+1)Rn+1
+
=
[email protected]+1n+1 @W(j+1)Rn+1 @Wj+1n+1
dFj+1,j @Fj+1,[email protected](j+1)Rn+1
=
[email protected](j+1)Rn+1 @Wjn+1
STIFFENED GAS (SG)
uL + FL(rL, pL;pI) = uIL
=
uIR = uR + FR(rR, pR; pI)
@FL @FL @FL
duL +
drL dpL dpI
+
+
@rL @[email protected]
@FR @FR @FR
= duR +
drR dpR dpI
+
+
@rR @[email protected]
dpI
dpI
dpI
dpI
dpI
dpI
, , , , ,
dpL
drL
dpR
drR
duL
duR
STIFFENED GAS
uI = uL + FL(rL, pL;pI) = uR + FR(rR, pR;pI)
duI
duI
duI
duI
duI
duI
, , , , ,
dpL
drL
dpR
drR
duL
duR
rIR, rIL
rIR = RR(rR,uI,pI)
rIL = RL(rL,uI,pI)
OTHER EOSs
 TaitTait
 TaitSG (which also implies TaitPG)
p = A(1  )eR1+ B(1  )eR2 + wre
wr
wr
R1r0
R2r0
r0
r0
r
r
JWL EOS
explosive products of combustion (and in particular
Trinitrotoluene — a.k.a. TNT)
where A, B, R1, R2, w and r0 are material constants
 Highly nonlinear function p(r,e)
 Presence of exponentials
uL + FL(rL, pL;rIL) = uIL
=
uIR = uR + FR(rR, pR; rIR)
GL(rL, pL; rIL) = pIL
pIR= GR(rR, pR; rIR)
=
JWL EOS
(1)
(2)
rIR,uIR ,pIR
t
rarefaction
c(r,p)
rR,uR ,pR
=
r
x
du
+
_
dr
= s
rw+1
p  AeR1+ BeR2
r0
r0
r
r
SGJWL RIEMANN SOLVER
(k)
(1)
(2)
complex Riemann problem
JWL EOS
 JWLJWL
 JWLSG (which also implies JWLPG)
TIMEINTEGRATORS
 BackwardEuler
 Threepoint Backward Differencing Formula (3PBDF)
~
dWi
Win+1  Win
=
dt
Dt
3PBDF
~
dWi
a0Win+1 a1Win + a2Win1
=
dt
Dt
where are constants
a0, a1, a2
3PBDF
 When node i has changed phase between tn1 and tn,
replace Win1 by W(i1)Rn1, the exact solution of the twophase
Riemann problem on the upstream side of the interface
at node i1 and time tn1
n1
n
W(i1)Rn1
i+1
i2
i
i1
2fin+1  2fin
=
Dt
d fi
dt
3PBDF FOR LEVEL SET
1 dfin

2 dt
where the last term can be estimated from the fluxes
at tn
LIMITATION
 Required to address phase change
timestep as necessary
r = 50 (kg/m3)
r = 1000.0 (kg/m3)
u= 0.0 (m/s)
u= 0.0 (m/s)
p= 105 (Pa)
p= 109 (Pa)
SHOCK TUBE PROBLEM
Air
Water
SHOCK TUBE PROBLEM
SHOCK TUBE PROBLEM
TURNER IMPLOSION PROBLEM
(0.5m, 0.5m)
Air (P = 105 Pa)
(0, 0)
Sensor
Water (P = 6.996 MPa)
z
(0.5m, 0.5m)
x
TURNER IMPLOSION PROBLEM
TURNER IMPLOSION PROBLEM
Explicit (RK2), CFL = 0.5
Implicit (3BDF), CFL = 100
TURNER IMPLOSION PROBLEM
TURNER IMPLOSION PROBLEM
TURNER IMPLOSION PROBLEM
 All runs with I/O on were performed with equal amount of I/O
speedup factor ~ 8.74 (FE/BE)
~ 25 (3PBDF/RK2)
pI, rIRus
x = x(t)
contact discontinuity
not involved
rarefaction*
t
structure
fluid
x
i
j
Mij
Wnj
rRuR pR
EMBEDDED FSI FRAMEWORK
n
w(x,0) =W , if x ≥ 0
j
w
F
(w)
= 0
+
x
t
* could also be a shock
u(x(t), t) = u (Mij)∙ nG(Mij)
s
pI, rIRus
x = x(t)
contact discontinuity
not involved
rarefaction*
t
structure
fluid
x
i
j
Mij
Wnj
rRuR pR
ONESIDED RIEMANN PROBLEM
us = uR + R2(pI(2); pR , rR)
 Closed form Jacobians exist as well (SG, Tait)
FLUX COMPUTATION
G
Mij
i
j
fluid 1
fluid 2
Fij= Roe (us, pI(1), Wni , EOS(1), uij)
Fji= Roe (us, pI(2), Wnj , EOS(2), uji)
EMBEDDED FSI (IMPLICIT)
 SGstructure (which also implies PGstructure)
 Taitstructure remains to be done
n1
n
W(i1)Rn1
structure
i+1
i2
i
i1
2D IMP45
air
( p = 14.5 psi )
water
( p = 1500 psi)
X 400
AEROF/DYNA3D
2D IMP45
 Explicit (RK2)
 Implicit (3PBDF)
 Implicit (BE)
2D IMP45
2D IMP45
Explicit (RK2), CFL=0.5
Implicit (3BDF), CFL=100
2D IMP45
 AEROS (nonlinear): implicit midpoint rule timeintegrator
 Same amount of I/O performed in all runs with I/O on
speedup factor = 35.5 (FE/BE)
= 43.0 (RK2/3PBDF)
SUMMARY
 Verification and validation using shock tube and Turner’s
implosion problems
 Achievement of speedup factors of 45
 Verification on a 2D implosion problem
 Achievement of speedup factor of 43 using the midpoint rule
in AEROS