# Numerical Simulation of Underwater Shock Wave Propagation in the Vicinity of Rigid Wall Based on Ghost Fluid Method

Authors: Ru-Chao Shi; Rui-Yuan Huang; Guang-Yong Wang; You-Kai Wang; Yong-Chi Li
Journal: Shock and Vibration
 (2015)

---

## Abstract

This paper presents numerical simulation of underwater shock wave propagation nearby complex rigid wall. The Ghost Fluid Method (GFM) for the treatment of complex rigid wall is developed. The theoretical analysis on the GFM-based algorithm and relevant numerical tests demonstrate that the GFM-based algorithm is first-order accurate as applied to complex rigid wall. A large number of challenging numerical tests show that the GFM-based algorithm is robust and quite simple in various practical problems. The numerical results on shock wave propagation in the vicinity of rigid wall are verified by comparing to exact solution and the results by body-fitted-grid method.


---

## Body

## 1. Introduction

Numerical simulation of underwater shock wave propagation nearby rigid wall has been a hot area of research in computational physics [1–3], computational biological fluid mechanics, [4] and computer graphics [5]. On the side of application for practical engineering problems, structure somewhile has to be assumed to be rigid wall to simplify computational models [6–10] or to present comparisons for numerical results on deformable solid and elastic-plastic solid [11–14]. There are three kinds of methods by so far to simulate the shock wave propagation in the vicinity of rigid wall. One of them, which is also the most popular, is body-fitted-grid method that need body-fitted grids to be constructed in curvilinear coordinate system according to the shape of rigid wall [15–19]. After the computational regions are discretized by body-fitted grids, numerical methods such as finite difference method, finite volume method, and finite element method can be employed to solve the flow field. Among these three numerical methods, finite difference method may lead to geometrically induced errors [20] while finite volume method and finite element method have shown stability in numerical simulation; finite element method especially is competent to various complicated problems. Another kind of method is Cartesian-grid methods that require the rigid wall to be treated by special algorithm that may be very complicated and can be divided into first-order method [21–24] and second-order method [25]. The application of this kind of method for underwater shock wave propagation can be found in [1] whereas there are very few applications of the Cartesian-grid methods since this kind of method is quite complicated in numerical implementation. The third kind is meshless methods (e.g., Smooth Particle Hydrodynamic method) that have widely applications for complex computational boundary but also have few applications for underwater shock wave propagation because of the large computational resources taken in numerical modelling and simulation.In this work, we employ the idea of Ghost Fluid Method (GFM) in Cartesian grids to treat complex rigid wall to simulate underwater shock wave propagation. Actually, GFMs are not accurate for defining boundary condition, such as solid-liquid interface, solid-gas interface, liquid-gas interface, and liquid-liquid interface. However in recent years GFMs have shown wide range applications [26–29] for shock wave propagation while interacting with material interface because of simplification and robustness. GFMs can be divided into original GFM (OGFM) and several improved versions. OGFM is developed by Fedkiw et al. [30] and the essence of nonoscillation has been verified. Nevertheless, OGFM is incapable in numerical simulation of multimedium that is very stiff at interface or where the density ratio is very large though it can be applied for low pressure gas-liquid flow. Then a modified GFM (MGFM) is developed by Liu et al. [31], which is robust for gas-liquid flow where the density ratio of the medium is very large. Next, real GFM (RGFM) [32] is developed to treat moving interface in some extreme situations. Subsequently some attempts in employing GFMs to treat liquid-solid interface and gas-solid interface have been made. A few of simple applications of GFM for fluid/solid coupling are presented by Fedkiw [33] but GFM is still found to be inaccurate in some situations. To overcome this problem, an essentially MGFM for fluid/elastic solid coupling is proposed and developed by Liu et al. [34] by employing Navier equation to solve solid and constructing an approximate Riemann problem at the interface.As GFM has the advantages of simplification, robustness, and flexibility, we apply the idea and essence of GFM for treatment of complex rigid wall. There are several problems that need to be solved. The dominant problem is how to define Riemann problem for water-rigid wall interface. The second problem is how to evaluate the flow states at ghost grid nodes. The last one is how to analyze the numerical accuracy of the GFM-based algorithm. The text below is arranged as follows. In Section2, we briefly present the methods, governing equations and numerical schemes, which are employed to solve the flow field. In Section 3, the GFM for the treatment of water-rigid wall interface and analysis on conservation errors are presented in detail. In Section 4, the theoretical solution for underwater shock wave propagation in the vicinity of rigid wall is provided to further test the numerical error of GFM-algorithm developed. In Section 5, various challenging tests are carried out by GFM-based algorithm, body-fitted-grids method, and theoretical analysis for further comparisons and discussions. In last section, brief conclusions are summarized and given.


## 2. Methods and Governing Equations

In this work, Euler equations are employed to solve flow field. While we apply GFM to rigid wall, level set method and related techniques are used to keep track of water-rigid wall interface and other material interfaces.


### 2.1. Euler Equations

The flow field is governed by Euler equations:(1)
ρ
ρ
u
ρ
v
ρ
w
E
t
+
ρ
u
ρ
u
2
+
p
ρ
u
v
ρ
u
w
u
E
+
p
x
+
ρ
v
ρ
v
u
ρ
v
2
+
p
ρ
v
w
v
E
+
P
y
+
ρ
w
ρ
w
u
ρ
w
v
ρ
w
2
+
p
w
E
+
P
z
=
0
,where ρ is density and u, v, and w are velocities at x, y, and z direction, respectively. E is total energy which can be expressed as(2)
E
=
ρ
u
2
+
v
2
+
w
2
2
+
ρ
e
.Here, e is internal energy per unit mass. 2D Euler equations are obtained by setting w
=
0. 1D Euler equations are obtained by setting v
=
0 and w
=
0 and can be given as(3)
∂
U
→
∂
t
+
∂
F
→
∂
x
=
0
,where U
→
=
ρ
,
ρ
u
,
E
T and F
→
=
ρ
u
,
ρ
u
2
+
p
,
u
E
+
p
T. The total energy E can be written as(4)
E
=
ρ
u
2
2
+
ρ
e
.In calculations, Euler equations are solved with fifth-order WENO spatial discretization [35] and second-order Runge-Kutta (R-K) time discretization.


### 2.2. Equations of State

For gas, we employγ-law which is expressed as(5)
p
=
ρ
e
γ
-
1
,where p is pressure, ρ is density, e is internal energy per unit mass, and γ is the ratio of specific heats. For high pressure gas, we can also use JWL equation [36]:(6)
p
=
A
1
-
ω
R
1
ρ
0
ρ
exp
⁡
-
R
1
ρ
0
ρ
+
B
1
-
ω
R
2
ρ
0
ρ
exp
⁡
-
R
2
ρ
0
ρ
+
ω
ρ
e
,where p is pressure, ρ is density, e is internal energy per unit mass, and A, B, R
1, R
2, and ω are constants.For water, Tait equation [36] and isentropic one-fluid model [37] are employed here and given as(7)
ρ
=
ρ
0
p
+
B
-
A
B
1
/
N
,
if
p
>
p
sat
⁡
,
ρ
=
k
ρ
g
cav
+
ρ
l
cav
p
+
B
-
A
/
p
cav
+
B
-
A
-
1
/
N
+
k
p
/
p
cav
-
1
/
γ
,
if
p
≤
p
sat
⁡
,where p
sat
⁡ is the saturated pressure of water, N
=
7.15, γ
=
1.33, ρ
0
=
1.0
×
1
0
3
kg
/
m
3, B
=
3.31
×
1
0
8
Pa, A
=
1.0
×
1
0
5
Pa, p
cav
=
p
sat
⁡, and the initial value of k is set to be k
=
0.001. For every iteration step, k would be updated by k
new
=
0.9
k if the result of p is beyond p
sat
⁡. Generally, p
sat
⁡
=
5000.0
Pa. To succinctly express the relationship between pressure p and density ρ, we rewrite the EOS of water that does not include gas phase as(8)
ρ
=
f
p
or
p
=
f
-
1
ρ
.One may find that pressure p of water can be determined by density ρ directly. For water, the internal energy per unit mass e is also the function of density ρ, which can be expressed as [36](9)
e
=
B
ρ
γ
-
1
γ
-
1
ρ
0
γ
+
B
+
A
ρ
.Since the internal energy e of water can be determined by the density ρ directly, it is not necessary to solve corresponding energy conservation equation in calculations.


### 2.3. Level Set Method

We use level set equation [38] to keep track of material interface:(10)
ϕ
t
+
u
ϕ
x
+
v
ϕ
y
+
w
ϕ
z
=
0
,where u, v, and w are velocities in x, y, and z directions, respectively. ϕ is the sign distance function of arbitrary point x
→ in the domain Ω and can be written as(11)
ϕ
=
-
Dis
∂
Ω
,
x
→
,
if
x
→
∈
Ω
-
,
0
,
if
x
→
∈
∂
Ω
,
Dis
∂
Ω
,
x
→
,
if
x
→
∈
Ω
+
.For 2D problems, w in (10) is set to be 0. For 1D problems, both v and w are set to be 0. The velocity extension technique [39] is employed to improve the accuracy of results by level set method. In calculations, level set equation is discretized by fifth-order HJ-WENO [40] and second-order R-K method.


### 2.4. Constant Extrapolation Approach

One of the most important steps of GFM is extrapolation that can extrapolate flow states from real region to ghost region. The partial equations for extrapolations can be expressed as [30, 41](12)
∂
I
∂
t
±
n
→
∇
I
=
0
,where t is artificial time, I is the flow states that need to be extrapolated, and n
→ is the vector normal material interface and will be given as n
→
=
∇
ϕ
/
∇
ϕ. The sign “±” is taken as “+” if the flow states need to be extended from Ω
- to Ω
+. If the flow states are extrapolated from Ω
+ to Ω
-, we take “-.” In calculations, (12) is simply solved by first-order upwind MinMod scheme.


## 2.1. Euler Equations

The flow field is governed by Euler equations:(1)
ρ
ρ
u
ρ
v
ρ
w
E
t
+
ρ
u
ρ
u
2
+
p
ρ
u
v
ρ
u
w
u
E
+
p
x
+
ρ
v
ρ
v
u
ρ
v
2
+
p
ρ
v
w
v
E
+
P
y
+
ρ
w
ρ
w
u
ρ
w
v
ρ
w
2
+
p
w
E
+
P
z
=
0
,where ρ is density and u, v, and w are velocities at x, y, and z direction, respectively. E is total energy which can be expressed as(2)
E
=
ρ
u
2
+
v
2
+
w
2
2
+
ρ
e
.Here, e is internal energy per unit mass. 2D Euler equations are obtained by setting w
=
0. 1D Euler equations are obtained by setting v
=
0 and w
=
0 and can be given as(3)
∂
U
→
∂
t
+
∂
F
→
∂
x
=
0
,where U
→
=
ρ
,
ρ
u
,
E
T and F
→
=
ρ
u
,
ρ
u
2
+
p
,
u
E
+
p
T. The total energy E can be written as(4)
E
=
ρ
u
2
2
+
ρ
e
.In calculations, Euler equations are solved with fifth-order WENO spatial discretization [35] and second-order Runge-Kutta (R-K) time discretization.


## 2.2. Equations of State

For gas, we employγ-law which is expressed as(5)
p
=
ρ
e
γ
-
1
,where p is pressure, ρ is density, e is internal energy per unit mass, and γ is the ratio of specific heats. For high pressure gas, we can also use JWL equation [36]:(6)
p
=
A
1
-
ω
R
1
ρ
0
ρ
exp
⁡
-
R
1
ρ
0
ρ
+
B
1
-
ω
R
2
ρ
0
ρ
exp
⁡
-
R
2
ρ
0
ρ
+
ω
ρ
e
,where p is pressure, ρ is density, e is internal energy per unit mass, and A, B, R
1, R
2, and ω are constants.For water, Tait equation [36] and isentropic one-fluid model [37] are employed here and given as(7)
ρ
=
ρ
0
p
+
B
-
A
B
1
/
N
,
if
p
>
p
sat
⁡
,
ρ
=
k
ρ
g
cav
+
ρ
l
cav
p
+
B
-
A
/
p
cav
+
B
-
A
-
1
/
N
+
k
p
/
p
cav
-
1
/
γ
,
if
p
≤
p
sat
⁡
,where p
sat
⁡ is the saturated pressure of water, N
=
7.15, γ
=
1.33, ρ
0
=
1.0
×
1
0
3
kg
/
m
3, B
=
3.31
×
1
0
8
Pa, A
=
1.0
×
1
0
5
Pa, p
cav
=
p
sat
⁡, and the initial value of k is set to be k
=
0.001. For every iteration step, k would be updated by k
new
=
0.9
k if the result of p is beyond p
sat
⁡. Generally, p
sat
⁡
=
5000.0
Pa. To succinctly express the relationship between pressure p and density ρ, we rewrite the EOS of water that does not include gas phase as(8)
ρ
=
f
p
or
p
=
f
-
1
ρ
.One may find that pressure p of water can be determined by density ρ directly. For water, the internal energy per unit mass e is also the function of density ρ, which can be expressed as [36](9)
e
=
B
ρ
γ
-
1
γ
-
1
ρ
0
γ
+
B
+
A
ρ
.Since the internal energy e of water can be determined by the density ρ directly, it is not necessary to solve corresponding energy conservation equation in calculations.


## 2.3. Level Set Method

We use level set equation [38] to keep track of material interface:(10)
ϕ
t
+
u
ϕ
x
+
v
ϕ
y
+
w
ϕ
z
=
0
,where u, v, and w are velocities in x, y, and z directions, respectively. ϕ is the sign distance function of arbitrary point x
→ in the domain Ω and can be written as(11)
ϕ
=
-
Dis
∂
Ω
,
x
→
,
if
x
→
∈
Ω
-
,
0
,
if
x
→
∈
∂
Ω
,
Dis
∂
Ω
,
x
→
,
if
x
→
∈
Ω
+
.For 2D problems, w in (10) is set to be 0. For 1D problems, both v and w are set to be 0. The velocity extension technique [39] is employed to improve the accuracy of results by level set method. In calculations, level set equation is discretized by fifth-order HJ-WENO [40] and second-order R-K method.


## 2.4. Constant Extrapolation Approach

One of the most important steps of GFM is extrapolation that can extrapolate flow states from real region to ghost region. The partial equations for extrapolations can be expressed as [30, 41](12)
∂
I
∂
t
±
n
→
∇
I
=
0
,where t is artificial time, I is the flow states that need to be extrapolated, and n
→ is the vector normal material interface and will be given as n
→
=
∇
ϕ
/
∇
ϕ. The sign “±” is taken as “+” if the flow states need to be extended from Ω
- to Ω
+. If the flow states are extrapolated from Ω
+ to Ω
-, we take “-.” In calculations, (12) is simply solved by first-order upwind MinMod scheme.


## 3. The GFM for Rigid Wall

GFMs include assuming that medium also exists outside the boundary of real medium, how to set ghost grid nodes for real medium and then how to update the flow states at these ghost grid nodes. In order to develop the GFM for rigid wall, we may first construct meshes on the entire domain neglecting whether or not the meshes are on the side of fluid (water). Then we extend the values at the grid nodes on the side of fluid (water) across fluid-rigid wall interface to the other side where ghost fluid is assumed to exist and grid nodes are correspondingly set as ghost nodes. While we try to complete the two steps mentioned above, the interface between real fluid and ghost fluid needs to be defined and we choose level set function to describe this interface. The algorithm for updating status of ghost fluid nodes should retrain real fluid from escaping when real fluid interacts with the interface and should behave the essential features of rigid wall on the opposite side of real fluid.


### 3.1. GFM-Based Algorithm

Supposing fluid is on left side of interface while rigid wall is on the right side of interface, we follow the ideology of [30, 31] and mathematically described the 1D GFM for treating fluid-rigid wall interface as(13)
∂
U
→
∂
t
+
A
∂
U
→
∂
x
=
0
,
U
→
x
,
0
=
U
→
L
,
if
x
<
x
I
,
U
→
R
∗
,
if
x
>
x
I
,where the subscripts L and R, respectively, represent the status on the left and right side of interface and the sign “∗” refers to the ghost status at the related grid nodes. The density and velocity on the side of rigid wall can be considered as(14)
ρ
R
=
∞
,
u
R
=
0
.The interfacial status at next time step will be predicted. Since GFM needs to behave the features of slip wall boundary for inviscid flow, we have(15)
u
I
pre
=
u
R
,where superscript “pre” indicates that the value of u
I is predicted for next time step. u
I
pre will be extrapolated from interface to ghost region; hence we have(16)
u
R
∗
=
u
I
pre
.By combining (14), (15), and (16), we have(17)
u
R
∗
=
0
.Since density and pressure must satisfy continuous condition, we also extrapolate the values at the real grid nodes that are in the neighborhood of interface to ghost grid nodes, which can be expressed as(18)
ρ
R
∗
=
ρ
L
,
p
R
∗
=
p
L
.The illustration for GFM algorithm mentioned above is shown in Figure 1. For multidimensional problem, velocity resolutions are tangential velocity V
→
T and normal velocity V
→
N. V
→
I
N
pre and V
→
R
N
∗ are obtained the same as u
I
pre and u
R
∗, and V
→
R
T
∗ is treated as ρ
R
∗ and p
R
∗:(19)
V
→
I
N
pre
=
0
,
V
→
R
N
∗
=
0
,
V
→
R
T
∗
=
V
→
L
T
,where V
→
L
T will be determined by(20)
V
→
L
T
=
V
→
L
-
V
→
L
·
n
→
L
.Correspondingly, the schematic for multidimensional GFM algorithm is shown in Figure 2.Figure 1
Illustrations for 1D GFM-based algorithm.Figure 2
Illustrations for multidimensional GFM-based algorithm.


### 3.2. Analysis on Conservation Errors

Suppose the fluid-rigid wall interface is located atx
w
∈
x
A
,
x
B where the rigid wall is on the right side of interface. Since the rigid wall is motionless, the location of interface x
w will be treated as constant. The increase within fluid domain x
A
,
x
w from the time t
n to t
n
+
1 can be given as(21)
E
U
→
=
∫
x
A
x
w
‍
U
→
n
+
1
-
U
→
n
d
x
.The flux that flows out of fluid domain x
A
,
x
w from the time t
n to t
n
+
1 can be expressed as(22)
E
F
→
=
∫
t
n
t
n
+
1
‍
F
→
w
-
F
→
A
d
t
.As U
→
=
ρ
,
ρ
u
,
E
T denotes three conservative variables, density, momentum, and energy, we have(23)
E
U
→
=
0
.By computing double integral of 1D Euler equations (3) with the ranges x
A
,
x
w and t
n
,
t
n
+
1, we have(24)
∫
t
n
t
n
+
1
‍
∫
x
A
x
w
‍
∂
U
→
∂
t
d
x
d
t
+
∫
t
n
t
n
+
1
‍
∫
x
A
x
w
‍
∂
F
→
∂
x
d
x
d
t
=
0
.Furthermore,(25)
∫
x
A
x
w
‍
U
→
n
+
1
-
U
→
n
d
x
+
∫
t
n
t
n
+
1
‍
F
→
w
-
F
→
A
d
t
=
0
.By combining (22), (23), and (25), we have(26)
E
F
→
=
0
.The total conservation error is the sum of E
U
→ and E
F
→ and will be written as(27)
TCE
U
→
,
F
→
=
E
U
→
+
E
F
→
.GFM for rigid wall has no ability to maintain conservation TCE
U
→
,
F
→
=
0 just like other applications of GFM. This is because the computational boundary of real medium is treated as medium-medium interface while this treatment is not absolute accurate. In numerical implementation, errors on mass conservation and shock location may occur. We only need to test the numerical accuracy of Δ
mass due to the linear relation between mass conservation error and shock location error. Let the computational domain be discretized by M meshes. Locate the interface at one mesh x
=
x
C
-
1
/
2
+
α
Δ
x where α
∈
0,1 and x
C
-
1
/
2 is grid node so that interface is between x
C
-
1
/
2
,
x
C
+
1
/
2 (see Figure 3). Then Δ
mass will be given as(28)
Δ
mass
=
∑
i
=
1
C
-
1
‍
ρ
i
Δ
x
+
ρ
C
α
Δ
x
-
∫
x
b
x
w
‍
ρ
exact
d
x
,where x
b is the grid node at the boundary. For WENO scheme and R-K method, the densities of all cells are not solved directly. Thus we use the density of two adjacent nodes to determine the density of cell, which can be expressed as ρ
i
=
ρ
i
-
1
/
2
+
ρ
i
+
1
/
2
/
2 where ρ
i is density of cell i and ρ
i
±
1
/
2 are the densities of grid nodes. The numerical test for conservation errors is presented in Case 1 of Section 5.Figure 3
Schematics for analysis onΔ
mass.


### 3.3. Numerical Implementation

The specific procedures for numerical simulation are summarized as follows.(1)
Construct grid in entire domain and divide the grid nodes into real ones and ghost ones according to fluid-rigid wall interface.
(2)
Define level set function (11) according to the fluid-rigid wall interface.
(3)
Extrapolate the flow states (ρ
L, p
L, and V
→
L
T) at real gird nodes just near interface and interfacial flow states V
→
I
N
pre predicted into a narrow band in ghost region via (12).
(4)
Solve the values at grid nodes in entire domain for next time step.
(5)
Keep track of interface by level set method.
(6)
Update the values of grid nodes according to the numerical results in step(
5
).
(7)
Go to step(
2
) for next time step. In step (
3
), the width of the band is set as 6
Δ
x to satisfy calculations based on fifth-order WENO spatial discretization and second-order R-K time discretization. In general, one can decide the width of band in accordance with the difference schemes for governing equations. The status of the ghost grid nodes outside the above-mentioned band can be set as constants, for example, the initial parameters of fluid. For steps (
2
) and (
5
), it is not necessary to define and solve the location by level set method for each time step because the rigid wall is motionless. However, if there are more than one kind of medium on the side of fluid, steps (
2
) and (
5
) should be repeated every time step.


## 3.1. GFM-Based Algorithm

Supposing fluid is on left side of interface while rigid wall is on the right side of interface, we follow the ideology of [30, 31] and mathematically described the 1D GFM for treating fluid-rigid wall interface as(13)
∂
U
→
∂
t
+
A
∂
U
→
∂
x
=
0
,
U
→
x
,
0
=
U
→
L
,
if
x
<
x
I
,
U
→
R
∗
,
if
x
>
x
I
,where the subscripts L and R, respectively, represent the status on the left and right side of interface and the sign “∗” refers to the ghost status at the related grid nodes. The density and velocity on the side of rigid wall can be considered as(14)
ρ
R
=
∞
,
u
R
=
0
.The interfacial status at next time step will be predicted. Since GFM needs to behave the features of slip wall boundary for inviscid flow, we have(15)
u
I
pre
=
u
R
,where superscript “pre” indicates that the value of u
I is predicted for next time step. u
I
pre will be extrapolated from interface to ghost region; hence we have(16)
u
R
∗
=
u
I
pre
.By combining (14), (15), and (16), we have(17)
u
R
∗
=
0
.Since density and pressure must satisfy continuous condition, we also extrapolate the values at the real grid nodes that are in the neighborhood of interface to ghost grid nodes, which can be expressed as(18)
ρ
R
∗
=
ρ
L
,
p
R
∗
=
p
L
.The illustration for GFM algorithm mentioned above is shown in Figure 1. For multidimensional problem, velocity resolutions are tangential velocity V
→
T and normal velocity V
→
N. V
→
I
N
pre and V
→
R
N
∗ are obtained the same as u
I
pre and u
R
∗, and V
→
R
T
∗ is treated as ρ
R
∗ and p
R
∗:(19)
V
→
I
N
pre
=
0
,
V
→
R
N
∗
=
0
,
V
→
R
T
∗
=
V
→
L
T
,where V
→
L
T will be determined by(20)
V
→
L
T
=
V
→
L
-
V
→
L
·
n
→
L
.Correspondingly, the schematic for multidimensional GFM algorithm is shown in Figure 2.Figure 1
Illustrations for 1D GFM-based algorithm.Figure 2
Illustrations for multidimensional GFM-based algorithm.


## 3.2. Analysis on Conservation Errors

Suppose the fluid-rigid wall interface is located atx
w
∈
x
A
,
x
B where the rigid wall is on the right side of interface. Since the rigid wall is motionless, the location of interface x
w will be treated as constant. The increase within fluid domain x
A
,
x
w from the time t
n to t
n
+
1 can be given as(21)
E
U
→
=
∫
x
A
x
w
‍
U
→
n
+
1
-
U
→
n
d
x
.The flux that flows out of fluid domain x
A
,
x
w from the time t
n to t
n
+
1 can be expressed as(22)
E
F
→
=
∫
t
n
t
n
+
1
‍
F
→
w
-
F
→
A
d
t
.As U
→
=
ρ
,
ρ
u
,
E
T denotes three conservative variables, density, momentum, and energy, we have(23)
E
U
→
=
0
.By computing double integral of 1D Euler equations (3) with the ranges x
A
,
x
w and t
n
,
t
n
+
1, we have(24)
∫
t
n
t
n
+
1
‍
∫
x
A
x
w
‍
∂
U
→
∂
t
d
x
d
t
+
∫
t
n
t
n
+
1
‍
∫
x
A
x
w
‍
∂
F
→
∂
x
d
x
d
t
=
0
.Furthermore,(25)
∫
x
A
x
w
‍
U
→
n
+
1
-
U
→
n
d
x
+
∫
t
n
t
n
+
1
‍
F
→
w
-
F
→
A
d
t
=
0
.By combining (22), (23), and (25), we have(26)
E
F
→
=
0
.The total conservation error is the sum of E
U
→ and E
F
→ and will be written as(27)
TCE
U
→
,
F
→
=
E
U
→
+
E
F
→
.GFM for rigid wall has no ability to maintain conservation TCE
U
→
,
F
→
=
0 just like other applications of GFM. This is because the computational boundary of real medium is treated as medium-medium interface while this treatment is not absolute accurate. In numerical implementation, errors on mass conservation and shock location may occur. We only need to test the numerical accuracy of Δ
mass due to the linear relation between mass conservation error and shock location error. Let the computational domain be discretized by M meshes. Locate the interface at one mesh x
=
x
C
-
1
/
2
+
α
Δ
x where α
∈
0,1 and x
C
-
1
/
2 is grid node so that interface is between x
C
-
1
/
2
,
x
C
+
1
/
2 (see Figure 3). Then Δ
mass will be given as(28)
Δ
mass
=
∑
i
=
1
C
-
1
‍
ρ
i
Δ
x
+
ρ
C
α
Δ
x
-
∫
x
b
x
w
‍
ρ
exact
d
x
,where x
b is the grid node at the boundary. For WENO scheme and R-K method, the densities of all cells are not solved directly. Thus we use the density of two adjacent nodes to determine the density of cell, which can be expressed as ρ
i
=
ρ
i
-
1
/
2
+
ρ
i
+
1
/
2
/
2 where ρ
i is density of cell i and ρ
i
±
1
/
2 are the densities of grid nodes. The numerical test for conservation errors is presented in Case 1 of Section 5.Figure 3
Schematics for analysis onΔ
mass.


## 3.3. Numerical Implementation

The specific procedures for numerical simulation are summarized as follows.(1)
Construct grid in entire domain and divide the grid nodes into real ones and ghost ones according to fluid-rigid wall interface.
(2)
Define level set function (11) according to the fluid-rigid wall interface.
(3)
Extrapolate the flow states (ρ
L, p
L, and V
→
L
T) at real gird nodes just near interface and interfacial flow states V
→
I
N
pre predicted into a narrow band in ghost region via (12).
(4)
Solve the values at grid nodes in entire domain for next time step.
(5)
Keep track of interface by level set method.
(6)
Update the values of grid nodes according to the numerical results in step(
5
).
(7)
Go to step(
2
) for next time step. In step (
3
), the width of the band is set as 6
Δ
x to satisfy calculations based on fifth-order WENO spatial discretization and second-order R-K time discretization. In general, one can decide the width of band in accordance with the difference schemes for governing equations. The status of the ghost grid nodes outside the above-mentioned band can be set as constants, for example, the initial parameters of fluid. For steps (
2
) and (
5
), it is not necessary to define and solve the location by level set method for each time step because the rigid wall is motionless. However, if there are more than one kind of medium on the side of fluid, steps (
2
) and (
5
) should be repeated every time step.


## 4. Exact Solution for Underwater Shock Wave Propagation

Suppose there is an underwater shock wave moving right and its velocity relative to the fluid in front of shock isW
s (W
s
>
0). We build 1D shock wave fixed coordinates as shown in Figure 4. According to equation of continuity, equation of momentum, we have(29)
ρ
1
u
~
1
=
ρ
2
u
~
2
,
p
1
+
ρ
1
u
~
1
2
=
p
2
+
ρ
2
u
~
2
2
,where the subscripts 1 and 2 indicate the status in front of shock and behind shock, respectively. The sign of tilde represents that the flow states are for moving coordinate system. Generally the flow states in front of shock such as pressure, density, and velocity (u
~
1
=
W
s) are already provided, and we can employ EOS for water (8) to close (29). Then the flow states behind shock in moving coordinates including velocity u
~
2, density ρ
2, and pressure p
2 can be solved. Eventually, one may obtain the velocity u
2 behind shock in global coordinates by(30)
u
~
2
=
W
s
+
u
1
-
u
2
,where u
1 and u
2 are the respective velocities in front of shock and behind shock. The velocity of shock wave in global coordinates is u
1
+
W
s.Figure 4
Schematics of moving shock wave (W
s
>
0). Left: global coordinates. Right: shock wave fixed coordinates.Suppose the underwater shock wave mentioned above impacts on rigid wall and results in a reflection wave (see Figure5). Similarly, we build reflection shock fixed coordinates. Due to conservation of mass and momentum, we have(31)
ρ
2
u
¯
2
=
ρ
3
u
¯
3
,
p
2
+
ρ
2
u
¯
2
2
=
p
3
+
ρ
3
u
¯
3
2
,where the subscript 2 indicates the flow states in front of reflection wave which is identical to the flow states behind the incident wave and the subscript 3 denotes the flow states behind reflection wave, and the bars signify that status is in reflection wave fixed coordinates. Since pressure and density in front of reflection wave, p
2 and ρ
2, can be determined by (29) and (30), we only need to use EOS of water and the relationship between shock wave coordinate system and global coordinate system to close (31). This relationship can be given as(32)
u
¯
2
=
W
r
s
,
u
¯
3
=
W
r
s
-
u
2
+
u
3
,where u
2 is the velocity (in global coordinates) in front of reflection wave and will be determined by (29) and (30), u
3 is the velocity (in global coordinates) behind reflection wave and will be set as u
3
=
0, and W
r
s (W
r
s
>
0) is the velocity of reflection shock wave relative to the water in front of shock. By combination of (31), (32), and EOS of water, one can obtain variables u
¯
2, u
¯
3, p
3, ρ
3, and W
r
s. The velocity of reflection wave in global coordinates is u
2
-
W
r
s.Figure 5
Schematics of reflection wave.


## 5. Results and Discussions

In this section, unless otherwise noted all calculations are done with employment of fifth-order WENO scheme to discretize spatial term of Euler equations, fifth-order HJ-WENO to discretize spatial term of level set equations, and second-order R-K method to discretize time term of both Euler equations and level set equations.Case 1.
The purpose of this case is to test the numerical accuracy of GFM-based algorithm. The computational domain is0,1 discretized with M cells (M
=
50,100,200,400,800) and M
+
1 grid nodes. A rigid wall is situated at one cell when M
=
800. The coordinate of rigid wall is set as x
=
0.88
-
α
Δ
x
M
=
800 where α
=
0.25,0.5,0.75 and an underwater shock wave is located at x
=
0.5 as shown in Figure 3. We will test numerical errors when M
=
50,100,200,400 to compare with numerical errors when M
=
800. Since the location of rigid wall is x
r
w
=
0.88
-
α
Δ
x
M
=
800, on the side of rigid wall, there are still six ghost nodal points that can be used to accomplish the numerical calculations when M is equal to 50 and fifth-order WENO and second-order RK scheme are employed. The nondimensional initial parameters in front of incident shock are p
1
=
1.0, ρ
1
=
1.0, u
1
=
0.0, and W
s
=
200.0 (means the velocity of shock wave is W
s
+
u
1
=
200.0). With the employment of theory in Section 4, the exact solution behind incident shock can be solved as p
2
=
4030.089077, ρ
2
=
1.112010, and u
2
=
20.145445. And the exact solution behind reflection wave can be solved as p
3
=
9968.150233, ρ
3
=
1.2034745, u
3
=
0, and W
r
s
=
265.069182 (means the velocity of reflection wave is u
2
-
W
r
s
=
-
244.923736). The computation is proceeded to t
=
2.93
×
1
0
-
3. The conservation errors are presented in Table 1. The numerical results show that GFM-based algorithm is first order. The numerical results on density by GFM depicted in Figure 6 concur very well with the theoretical results.Table 1
Conservation errors on GFM-based algorithm (Case1).
Δ
x
α
=
0.25
α
=
0.5
α
=
0.75
Δ
mass
Order
Δ
mass
Order
Δ
mass
Order
1
/
50
-
3.951216
×
1
0
-
3
—
-
3.886237
×
1
0
-
3
—
-
3.821259
×
1
0
-
3
—
1
/
100
-
2.040762
×
1
0
-
3
0.9532
-
1.975817
×
1
0
-
3
0.9759
-
1.910871
×
1
0
-
3
0.9998
1
/
200
-
1.016700
×
1
0
-
3
1.0052
-
9.517675
×
1
0
-
4
1.0538
-
8.868354
×
1
0
-
4
1.1075
1
/
400
-
5.443664
×
1
0
-
4
0.9012
-
4.794354
×
1
0
-
4
0.9893
-
4.145043
×
1
0
-
4
1.0973
1
/
800
-
2.942627
×
1
0
-
4
0.8875
-
2.293315
×
1
0
-
4
1.0639
-
1.644002
×
1
0
-
4
1.3342Figure 6
Comparisons between numerical results and exact solution, first line:Δ
x
=
1
/
50, second line: Δ
x
=
1
/
100, third line: Δ
x
=
1
/
200, fourth line: Δ
x
=
1
/
400, and fifth line: Δ
x
=
1
/
800 (Case 1).Case 2.
This is a 1D case for practical problem to test the performance of GFM for rigid wall. In this case, shock wave generated by underwater explosion hits rigid wall and then results in a reflection wave. The computational domain is[
0,1
] with 100 meshes uniformly distributed. The interface between explosive and water is at x
=
0.3. The rigid wall is located at x
=
0.89875. The regions of explosive and water are [
0,3.0
] and [
3.0,0.89875
], respectively. The boundary type at x
=
0 and x
=
1 is symmetric and outflow, respectively. The initial flow states for water are p
w
=
1.0
×
1
0
5
Pa, ρ
w
=
1.0
×
1
0
3
kg
/
m
3, and u
w
=
0. The initial status for explosive is p
g
=
7.8039
×
1
0
9
Pa, ρ
g
=
1.63
×
1
0
3
kg
/
m
3, and u
g
=
0 which are the parameters of TNT charge [36] and where JWL EOS needs to be employed. In addition, the grid nodes bordering explosive-water interface is updated by RGFM algorithm [32]. The computation is run to 19.79
ms. Figures 7(a), 7(b), and 7(c), respectively, show pressure, density, and velocity profiles. The theoretical solution is also provided for comparison. There are few discrepancies between the numerical results and exact solution and the shock front are predicted accurately.Figure 7
Numerical results by GFM-based algorithm att
=
19.79
ms (Case 2).
(a)
(b)
(c)Case 3.
This is also a 1D case for practical problem. In this case, the underwater shock wave is generated by underwater explosion whereas explosive gas is solved byγ-law. The computational domain is [
0,1
]. The water-rigid wall interface is located at x
=
0.9. The radius of explosive is 0.5. The left boundary type of entire domain is symmetric. The initial flow states of explosive are γ
g
=
2.0, p
g
=
8.29
×
1
0
3
bar, and ρ
g
=
1.27
×
1
0
3
kg
/
m
3, and the initial flow states of water are p
w
=
1.0
bar, ρ
w
=
1.0
×
1
0
3
kg
/
m
3. The computation is run to t
=
0.251
ms. The numerical results by the GFM presented are depicted in Figure 8 and show that the density coincides very well with the theoretical results.Figure 8
Numerical results by GFM-based algorithm att
=
0.251
ms (Case 3).
(a)
(b)
(c)Case 4.
This is a 2D case to test the performance of GFM for complex rigid wall in multidimensional practical problem. In this case, an underwater explosion of TNT charge occurs near complex rigid wall and then shock waves appear and propagate in the immediate vicinity of rigid wall. While GFM-based algorithm is employed to present numerical results, curvilinear coordinate system (η
,
ξ) is also constructed according to the shape of rigid wall to provide comparisons between numerical results by GFM-based algorithm and body-fitted-grid method. The illustration of computational region is shown in Figure 9. The center of explosive is located at the origin (
0,0
). The initial flow states are p
w
=
1.0
×
1
0
5
Pa, ρ
w
=
1.0
×
1
0
3
kg
/
m
3, u
w
=
0, v
w
=
0, p
g
=
7.8039
×
1
0
9
Pa, ρ
g
=
1.63
×
1
0
3
kg
/
m
3, u
g
=
0, v
g
=
0, and γ
g
=
2.0. The computational domain is [
-
3,3
]
×
[
-
3,3
]. The shape of rigid wall satisfies y
=
(
1
/
4
)
sin
(
x
±
3.0
). For calculations based on GFM, there are 200
×
250 cells distributed uniformly to discretize computational domain with Δ
x
=
Δ
y. For calculations based on body-fitted-grid method, the computational domain on the side of fluid is discretized by 200
×
200 cells where coordinates of nodal points satisfy x
=
ξ and y
=
(
1
/
4
)
sin
(
ξ
)
+
η. Thus the cells in curvilinear coordinate system are smooth and the coordinates at every gird nodes are derivative. For each computational time step, the grid nodes bordering explosive-water interface are updated by RGFM algorithm. In Figures 10(a), 10(b), and 10(c), the differences between the results based on GFM and body-fitted-grid method are shown. Obviously, both of these two methods can be used to capture the underwater shock save reflected from rigid wall. Few discernible discrepancies between the shape and locations of reflection waves, respectively, obtained by GFM and body-fitted-grid method can be perceptible. In Figure 11, we provide pressure distribution at the bottom wall with different grid step. The interface obtained by GFM is implicit; thus, the pressure on water-rigid wall interface is not able to be solved directly. We output the pressure at the grid nodes just bordering the rigid wall. Since there is distance between these grid nodes and rigid wall, and reflection wave propagates with decreasing strength, the pressure distribution output by GFM is below the results by body-fitted-grid method.Figure 9
Schematics of computational domain (Case4).Figure 10
Comparisons between numerical results by GFM-based algorithm (left) and body-fitted-grid method (right). Red dashed line indicates material interface, (a)t
=
78.86
ms, (b) t
=
83.07
ms, and (c) t
=
84.73
ms (Case 4).
(a)
(b)
(c)Figure 11
Pressure distribution in the immediate vicinity of rigid wall (Case4).Case 5.
This is a 3D case which is taken directly form [37]. The related experiment results on this case have been presented in [42]. A sphere PETN charge is located in the center of a container full of water. The diameter and height of the container are 8.889
×
1
0
-
2
m and 2.286
×
1
0
-
1
m, respectively. The diameter of explosive is 3.0
×
1
0
-
2
m. Since we cut the container in half, the computational domain is set to be [
0,0.051816
]
×
[
0,0.10688
]
×
[
0,0.2286
] and discretized by 34
×
70
×
150 uniform cells. The initial status of water and explosive is P
w
=
1.0
×
1
0
5
Pa, ρ
w
=
1.0
×
1
0
3
kg
/
m
3, u
w
=
0, v
w
=
0, w
w
=
0, p
g
=
2.0
×
1
0
9
Pa, ρ
g
=
1.77
×
1
0
3
kg
/
m
3, u
g
=
0, v
g
=
0, w
g
=
0, and γ
g
=
2.0. For each computational time step, the grid nodes bordering explosive-water interface are updated by MGFM algorithm [31]. Figure 12 shows the processes for underwater shock wave propagation in rigid container. The pressure contours on meridian plane by using 2D axis-symmetric Euler equation has been presented in [37]. For ease of comparison we also present the pressure contours and pressure cloud pictures on the meridian plane. Once the explosion starts, an underwater shock wave is generated and propagates to cylindrical wall with decreasing strength resulting in a reflection wave from the wall. The reflection wave hits the explosive-water interface and generates a rarefaction wave as shown in Figure 12(a). The rarefaction wave moves toward column wall and then makes a reflection wave again as shown in Figure 12(b). This new reflection wave is also very strong and creates cavitation flow in the vicinity of rigid wall as shown in Figure 12(c). Meanwhile the rarefaction wave toward basal plan is so strong that cavitation region is created in the immediate vicinity of gas-water interface as shown in Figure 12(d). We may note that numerical results shown in Figure 12 are consistent with numerical results by 2D axis-symmetric Euler equation and experiments. In Figure 13, we compare the results presented, the results by 2D axis-symmetric Euler equation and experimental results. It can be noticed that all the contour lines are almost uniformly distributed. We export the pressure history at the internal grid node 20,14,76 to compare with pressure history at the boundary by 2D axis-symmetric Euler equations. There are few discrepancies between the two profiles. Since elastic-plastic solid is assumed to be rigid body in both the numerical models developed in this paper and the numerical method developed in [37], one may find that the third profile by experiments decreases more rapidly after peaking at about 7000
bar.Figure 12
Pressure cloud pictures and contours by GFM-based algorithm. Blue lines indicate the cavitation region. Red curved surface and red lines represent rigid wall, (a)t
=
30
μ
s, (b) t
=
40
μ
s, (c) t
=
60
μ
s, and (d) t
=
70
μ
s (Case 5).
(a)
(b)
(c)
(d)Figure 13
Comparisons among the results presented, the results by 2D axis-symmetric equations and experimental results (Case5).Case 6.
This challenging case is used to test the performance of GFM on the treatment of quite complex rigid wall. In Figure14, the explosive is located in the center of ellipsoidal vessel whose radii are respectively a
=
1.0
m, b
=
0.5
m and c
=
0.5
m. Constructing Cartesian coordinate system along with the axis of vessel while setting the center of vessel is the origin (
0,0
), we can write the equation of ellipsoid as x
2
/
1
.
0
2
+
y
2
/
0
.
5
2
+
z
2
/
0
.
5
2
=
1.0. The length of explosive is set to be 0.3
m. The angles of inclination between the axis of column and Cartesian coordinate system are 3
π
/
4, 0, and π
/
4 respectively. The computational domain is defined as -
1.25,1.25
×
-
0.75,0.75
×
-
0.75,0.75 with 126
×
76
×
76 grid nodes uniformly distributed. The initial flow states are γ
g
=
2.0, p
g
=
8.29
×
1
0
3
bar, ρ
g
=
1.27
×
1
0
3
kg
/
m
3, p
l
=
1.0
bar, and ρ
l
=
1.0
×
1
0
3
kg
/
m
3. Once the explosion is initiated, a strong shock wave is generated and propagates outwards the wall of vessel with exponentially decreasing strength while high pressure bubble is expanding rapidly. As shown in Figure 15 shock wave hits rigid wall resulting in a strong reflection wave. Figure 16 shows that reflection wave hits the bubble and then generates rarefaction wave while expanding of bubble becomes slow because of the impact of shock wave. Figure 17 shows that the reflection waves from the rigid wall intersect and interact with each other resulting in another new strong reflection shock wave which propagates towards the wall at the end of long axis of vessel. Figure 18 shows that this new reflection wave reflects from the end of vessel and then results in third reflection wave that propagates toward the bubble (shown in Figure 19). This third reflection wave is so strong that cavitation flow (see Figure 20) occurs in large areas nearby rigid wall. One can notice that no contour distortion is found and all the contour lines are nearly regularly distributed. This also demonstrates the efficiency and robustness of GFM for rigid wall.Figure 14
Schematics of computational region (Case6).Figure 15
Numerical results by GFM-based algorithmt
=
0.259
ms. Left: pressure cloud pictures, right: bubble shape (Case 6).Figure 16
Numerical results by GFM-based algorithmt
=
0.437
ms. Left: pressure cloud pictures and pressure contours, right: bubble shape (Case 6).Figure 17
Numerical results by GFM-based algorithmt
=
0.487
ms. Left: pressure cloud pictures, right: bubble shape (Case 6).Figure 18
Numerical results by GFM-based algorithm att
=
0.644
ms. Left: pressure cloud pictures, right: bubble shape (Case 6).Figure 19
Numerical results by GFM-based algorithmt
=
0.816
ms. Left: pressure cloud pictures, right: bubble shape (Case 6).Figure 20
Numerical results by GFM-based algorithmt
=
10.466
ms. Left: pressure cloud pictures, right: bubble shape (Case 6).


## 6. Conclusions

In this work, we propose and develop a GFM-based algorithm for the treatment of complex rigid wall and present numerical simulation of underwater shock wave propagation nearby rigid wall. GFM for complex rigid wall is tested and verified by comparing to exact solution and the results by body-fitted-grid method. Numerical simulation on strong underwater shock wave is more challenging and more difficult than numerical simulation for normal shock wave so most of our tests focus on extreme situations where the pressure is quite high. GFM for complex rigid wall is verified to be first-order accuracy on mass conservation and very simple and robust in various practical problems.



---

*Source: 467376-2015-04-30.xml*