Investigation of underwater explosion near composite structures using a combined RKDG-FEM approach

Journal Pre-proof Investigation of underwater explosion near composite structures using a combined RKDG-FEM approach

Liang Ge, A-Man Zhang, Shi-Ping Wang

PII:

S0021-9991(19)30818-6

DOI:

https://doi.org/10.1016/j.jcp.2019.109113

Reference:

YJCPH 109113

To appear in:

Journal of Computational Physics

Received date:

15 March 2017

Revised date:

30 June 2019

Accepted date:

9 November 2019

Please cite this article as: L. Ge et al., Investigation of underwater explosion near composite structures using a combined RKDG-FEM approach, J. Comput. Phys. (2019), 109113, doi: https://doi.org/10.1016/j.jcp.2019.109113.

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Highlights • A combined RKDG-FEM approach is developed for simulating UNDEX problems. • UNDEX involving shock-bubble-structure interactions and flow cavitation is simulated. • Comparison study of UNDEX near steel and composite sandwich structures is conducted.

Investigation of underwater explosion near composite structures using a combined RKDG-FEM approach Liang Ge, A-Man Zhang*, Shi-Ping Wang College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China

Abstract A combined RKDG-FEM approach, using the Runge-Kutta discontinuous Galerkin (RKDG) method [1, 2] for fluid simulating and the finite element method (FEM) [3, 4] for structure simulating, is developed to study the complex interactions between compressible multiphase flows and deformable steel and composite structures in near-field early-time underwater explosion (UNDEX). The level set method (LSM) [5, 6] and the real ghost fluid method (GFM) [7, 8] are applied to capture and treat the moving gas-water interface, respectively. A GFM-based interface treating scheme [9, 10] is adopted to treat the fluid-structure interface, achieving the coupling by solving simultaneously the fluid characteristic equation and the structure equation of motion at the material interface to satisfy the kinematic compatibility and pressure equilibrium conditions. The isentropic one-fluid cavitation model in [11] is employed to capture the transient flow cavitation. Accuracy and effectiveness of the combined RKDG-FEM approach are validated by five benchmark cases, and then it is applied to simulate early-time UNDEX problems near a submerged aluminum foam-cored sandwich plate and a monolithic steel plate of equivalent mass. Numerical results show that the sandwich plate has a higher underwater explosion load resistance through absorbing plastic

*Corresponding

author.

E-mail address: [email protected] (A.-M. Zhang), [email protected] (L. Ge) 1

energy by the foam core than the monolithic steel plate of equivalent mass, and that the foam core plays an important role in the fluid-structure interaction (FSI) and flow cavitation effects, which can be exploited to improve the UNDEX mitigation capability of composite sandwich structures. Key words͵Combined RKDG-FEM approach; UNDEX; Multiphase flows; FSI; Flow cavitation; Composite structure.

1 Introduction Naval structure suffers from severe marine environments during its operation involving corrosive sea water, temperature extremes, and even transient impulsive loads resulting from weapon impact, air explosion, or UNDEX [12, 13]. With the development of precision-guided weapons, the naval structures are far more likely to be attacked by nearly-field UNDEX. The physical process of near-field UNDEX is extremely complex, involving various complicated issues such as shock wave propagation, compressible multiphase flows with strong shock impacting material interface, and transient FSI, etc [14-18]. Composite structures are attractive for marine structures as lightweight construction, environmental benefits, thermal and sound insulations, and corrosive resistance [19]; they can also provide good UNDEX resistance due to high shear and bending resistances, high strength-to-weight ratios, and etc [20]. Many experimental and numerical investigations have been being conducted to research the deformation, damage, and failure modes of the sandwich structure subjected to high-intensity impulsive loads of air explosion and UNDEX [21-24]. The dynamic response and damage characteristic of sandwich structures are extremely complicated due to the material heterogeneity, competing damage mechanisms, complex failure modes, and FSI effects [21, 25-27]. It has been found that the deflection and transmitted impulse experienced by sandwich plates 2

are significant smaller than by equivalent monolithic plates made exclusively of either the facesheet material or the core material [28-30]. Resent researches indicate that the FSI effect is crucial in the compressive response and can be exploited to enhance the blast mitigation capability of sandwich composite structures [24, 31-33]. The simulations of compressible multiphase flows and transient FSI are still challenging topics, especially if the material interfaces are under strong shock impacting, as in near-field UNDEX [34, 35]. Currently available numerical methods founded in the literature for transient FSI problems can be generally classified into three classical categories: Eulerian-Eulerian (E-E) methods, Lagrangian-Lagrangian (L-L) methods, and Eulerian-Lagrangian (E-L) methods [9, 35]. E-E methods treat both the fluid and the structure by a unified Eulerian framework, and they can be seamlessly coupled by the volume of fluid method (VOF) or some other interface capturing techniques [35, 36]. As suffering from severe numerical diffusion and hard to track the time-dependent material properties which are rather important to model material response and damage to dynamic loading, Eulerian numerical methods perform poorly for many solid dynamic calculations, and they are also extremely difficult to construct an appropriate equation of state for complex anisotropic materials such as composite structures [10]. On the contrary, L-L methods treat both the fluid and the structure with pure Lagrangian numerical methods [22, 37]; therefore they encounter difficulties in treating flow with strong shock waves, large deformation, or vorticity, since these can cause large deformation or even distortion of Lagrangian meshes and then stop the calculation [35, 37, 38]. By using Eulerian numerical methods on fluid domain and Lagrangian numerical methods on structure domain, E-L methods combine the advantages of the Eulerian and

3

Lagrangian methods [35, 39]. One of the most popular E-L methods is the arbitrary Lagrangian Eulerian method (ALE), coupling the fluid and structure domains by a “mushy” mesh region with a smoothly velocity varying from the velocity of the stationary Eulerian mesh to that of the moving Lagrangian mesh [39]. Many other E-L methods such as the immersed boundary method (IBM) [40] and the embedded boundary method (EBM) [41] are also developed to simulate transmit FSI problems. Recently, the GFM has also been applied to study the transmit FSI problems [9, 10, 35, 42, 43]. The key point of this approach is to treat the material interface as internal boundary, and then the internal boundary condition is imposed by the

variable extrapolation from one medium to another.

In the present work, we use the RGFM [8] to treat the gas-water interface, the GFM-based interface treating scheme developed in [9, 10] to treat the fluid-structure interface, and the RKDG method and the FEM to simulate the fluid and the structure, respectively. The objective of this work is to develop an efficient combined RKDG-FEM approach to investigate UNDEX near deformable steel and sandwich structures involving complex interactions among compressible multiphase flows, FSI, and flow cavitation. Section 2 presents the governing equations of the Eulerian flow and the Lagrangian structural model. Section 3 presents the detailed numerical implementations of the combined RKDG-FEM approach. Section 4 presents the validation studies for the combined approach. In section 5, simulations are carried out for two configurations: (1) a configuration of UNDEX near a submerged clamped sandwich plate, and (2) a configuration of UNDEX near a submerged steel plate of same areal mass. A brief discussion of current and future works is given in Sect. 6.

4

2 Governing equations 2.1 Eulerian flow model For near-field early-time UNDEX, as the characteristic time of molecular diffusion induced by viscosity is far smaller than that of explosion shockwave and gas dynamics [44], the viscous effects can be neglecting in Eulerian flow model. Hence, the inviscid gas (explosive gas), water, and cavitating flow can all be described by the following Euler equation [45, 46]: U t + Fx ( U ) + G y ( U ) = S ( U ) ,

(1)

where U is the conservative variable vector, F and G is invisicid flux vectors, and S is source term vector, which are defined as

ª ρu º ª ρv º ª ρu º ªu º « ρu 2 + p » « ρ uv » « » 2 « ρu » n − 1 « ρu « » « » ». « » U= , F (U) = , G (U) = , S (U) = − « ρ uv » « ρv2 + p » « ρv » x « ρ uv » « » « » « » « » ¬E¼ ¬( E + p ) u ¼ ¬( E + p ) v ¼ ¬( E + p ) u ¼ Here ρ and p denote the fluid density and pressure, respectively. u and v denote the flow velocities in the x (or radial) and y (or vertical) directions, respectively. E denotes the unit volume total energy, given by

E = ρe +

1 ρ (u 2 + v2 ) , 2

(2)

where e is defined as the specific internal energy, and n is a system parameter with n = 1 and

v = 0 for one-dimensional problems, n = 1 for two-dimensional problems, and n= 2 for axial-symmetric. In order to close the governing system (1), equation of state (EOS) are obtained. The ideal gas and the Tait EOS are employed to model gas and water, respectively. The ideal gas EOS is given by

5

p = ( γ g − 1) ρ e,

(3)

where γ g is the ratio of specific heats, which is set to 1.4 for air and 2.0 for explosive gas [47, 48]. The Tait EOS is given by γ p = B ª( ρ ρ0 ) w − 1º + A, ¬ ¼

(4)

where ρ 0 = 1025 kg/m 3 is the initial density of water, and γ w = 7.15 , A = 1×105 Pa , and

B = 3.31× 108 Pa are material constants [49, 50]. For near-filed UNDEX, as the complex interactions of shock wave, bubble, and material interface, the pressure of the flow region near the material interface may reach toward the limit of vapor pressure, and thus cavitation may occur at this low pressure region. Therefore, an special EOS is required for the cavitation regions. In current work, the isentropic one-fluid cavitation model [11] is used to capture the transient cavitation regions. The averaged density of the cavitating flow is defined as

ρ = α ⋅ ρ g + (1 − α ) ⋅ ρl ,

(5)

where α ∈[ 0,1] is the void fraction, which is set to 0 for pure water and 1 for pure gas; ρ g and

ρ l are the densities of the vapor and water, respectively. In this cavitation model, under the assumption that the cavitating flow is a homogenous and barotropic mixture of gas and liquid [51, 52], the theoretical sound is defined as −

1

­° ª α 1 − α º ½° 2 . + a = ®ª¬αρg + (1 − α ) ρl º¼ « 2 2 »¾ ρ ρ a a « » l l ¼¿ ° ¬ g g ¯°

(6)

Under the assumption that the cavitation evolution is only driven by pressure drop [11], the void fraction α is governed and derived by

6

§ 1 dα 1 = α (1 − α ) ¨ − ¨ ρ a2 ρ a2 dp g g © l l

· ¸¸ . ¹

(7)

As the components of the cavitating mixture are assumed to maintain isentropic property [11], by integrating (6) and (7), resulting in 1

( p pcav ) γ , α =k 1 1−α p p γ w

(

ρ= where

p = p + B − A and

cav

)

(8)

g

k ρ gcav + ρ lcav

(p

pcav )

−

1

γw

+ k ( p pcav )

−

1

,

γg

(9)

pcav = pcav + B − A . k = α 0 (1 − α 0 ) is calculated by a procedure

cav cav developed in [11]. ρ g and ρl are respectively the gas and liquid densities when their pressures

are equal to the cavitation pressure pcav .

2.2 Lagrangian structure model The explicit dynamic FEM [3, 4, 10] along with appropriate material models is used to model the deformable structures. The Lagrangian governing equations for structure are s = ∇ ⋅ ı + f , Ω s × [ 0, T ] , ρsu

(10)

ı ⋅ n = t , Γ t × [ 0, T ] ,

(11)

u s = u , Γ u × [ 0, T ] ,

(12)

u s ( x, 0 ) = u 0 ( x ) ,

x ∈ Ωs ,

(13)

u s ( x , 0 ) = u 0 ( x ) ,

x ∈ Ωs .

(14)

s are the structure displacement, velocity, and acceleration vectors, respectively; Her us , u s , and u

ρ s , ı , and f are the structure density, body force vector, and Cauchy stress tensor, respectively; x and Ωs are the structure spatial coordinates and computational domain, respectively; Γt and 7

Γu denote the boundaries of surface traction and displacement, respectively; t is the predefined surface traction boundary condition, and u is the predefined displacement boundary condition; u0 is the initial displacement condition, and u 0 is the initial velocity condition. With initial conditions, boundary conditions, kinematic relations, and appropriate constitutive models, the finite element discrete equation [10] for system (10)-(14) is given as

s = Fint + F ext , Mu

(15)

where

M = ³ ρ s NT NdV , Fint = − ³ ∇NsdV , F ext = ³ NfdV + ³ Ntd Γ. Ωs

Ωs

Ωs

Γt

Here M is the diagonal lumped mass matrix, F int denotes the vector of internal force, F ext denote the vector of external force, and N denotes the displacement interpolation matrix. Eq. (15) is integrated using the explicit central difference integration rule [10], given as 1

1

(i + ) (i − ) u s 2 = u s 2 +

u

(i +1) s

1  (i) Δt (i +1) + Δt (i) ) u ( s , 2

(i) s

= u + Δt

(i +1)

u

1 (i + ) 2 s

(16)

,

where the superscript ( i ) refers to the increment number;

(i− ) 1 2

and

(i+ ) 1 2

refer to the middle

increment values. In this paper, the commercial software ABAQUS/Explicit [53] is adopted as the FEM solver. The composite sandwich plate comprises two steel facesheets and an aluminum foam core. The steel facesheets are modeled by the Johnson-Cook constitutive model [54, 55]:

(

σ d = A + B (ε p )

n

) (1+ C ln ε )(1− T *

*m

),

(17)

where σ d is the dynamic flow stress, ε p the plastic strain in equivalent form, and A , B , n , C , and m are constants for material. ε* = ε ε0 is the plastic strain rate in dimensionless form, where

8

ε0 =1.0 s −1 is a reference strain rate. T * = (T − Tr ) (Tm − Tr ) ∈ [ 0,1] is the temperature in dimensionless form, where Tr and Tm are the room and the material melting temperature, respectively. The aluminum foam core is modelled as a compressible continuum by an isotropic metal foam constitutive model developed by Deshpande and Fleck in [56]. The material yield stress function for metal foam is defined as

Φ = σˆ − σ y ,

(18)

where Φ is the isotropic yield surface, σ y is the yield stress, and σˆ is equivalent stress, given by

σˆ 2 =

σ e2 + α 2σ m2 , 2 1 + (α 3)

(19)

where σ e denotes the von Mises effective stress, σm denotes the mean stress, and α is a material parameter denoting the yield surface’s shape. This paper employs an over-stress model [56], and the yield stress σ y is defined as

σ y =σ c + ηεˆ p ,

(20)

p where σ c ( εˆ ) denotes the static uniaxial stress versus plastic relation, η is the viscosity, and εˆ p

is the plastic strain rate.

3 Numerical methods 3.1 Review of RKDG method The discontinuous Galerkin method [1, 2] is used to discretize Eq. (1) in space. To define the DG method, we give a partition of the computational domain Ω in a collection of nonoverlapping

9

elements K , and the boundary of element K is defined as ∂K . For any time t ∈ [ 0, T ] , the test function Φ ( x ) and approximate solution U ( x, t ) are sought in the finite element space of discontinuous functions, and the space is defined as

Vh = { ph ∈ L∞ ( Ω) : ph K ∈ Pk , ∀K ∈,h } ,

(21)

where ,h is a triangulation of the domain Ω , and P k is a set of polynomials of degree up to k, where k is set to 2 for all cases in current work. The Eq. (1)’s weak formulation can be obtained by multiplying a test function Φ ( x ) ∈Vh on both sides, integrating over an single cell K , and integrating by parts, resulting in

d U ( x, t ) Φ ( x )dA + ¦ ³ H ( U ( x, t ) ) ⋅ ne, K Φ ( x ) d Γ e dt ³K e∈∂K

(22)

−³ H ( U ( x, t ) ) ⋅∇Φ ( x ) dA = ³ S ( U ( x, t ) )Φ ( x ) dA, K

K

where H ( U ) = ( F ( U ) , G ( U ) ) , and ne, K is the unit normal of edge e in outward direction. The line and area integrals in Eq. (22) are typically discretized by a Gaussian quadrature with sufficient high accuracy as

³

e

L

H ( U ( x, t ) ) ⋅ n e, K Φ ( x ) d Γ ≈ e ¦ ª¬ωi H ( U ( xei , t ) ) ⋅ n e, K Φ ( xei ) º¼,

(23)

i =1

³

K

H ( U ( x, t ) ) ⋅∇Φ ( x ) dA ≈ K

M

¦ ª¬ω H ( U ( x , t ) ) ⋅∇Φ ( x )º¼, j

Kj

Kj

(24)

j =1

where xei and x K j are the Gaussian points in edge e and element K , respectively; e is the length of edge e , and K

is the area of element K . The flux H ( U ( x ei , t ) ) ⋅ n e, K in Eq. (23) is

l ( x , t ) in [57, 58], the exact solution U ( x,t ) is by the replaced by the HLLC numerical flux H ei  ( x, t ) , and the test function Φ ( x ) by Φ  ( x ) ∈ V k , resulting in approximate solution U h

10

L d l (x ,t)Φ  ( x, t ) Φ  ( x )dA + ¦ ¦ ω H  (x ) e − U i ei ei dt ³K e∈∂K i =1 M

(

M

)

(

(25)

)

¦ω j H U ( xKj , t ) ⋅∇Φ ( xKj ) K = ¦ω jS U ( xKj , t ) Φ ( xKj ) K . j =1

j =1

 ( x, t ) of element K is defined as The approximation solution U Nk

 ( x, t ) = ¦ U ( m ) ( t ) Φ  ( K ) ( x ), U K K m

(26)

m=0

K ( m) where U K ( t ) ( m = 0," , N k ) are the degrees of freedom, and Φ (m ) ( x ) ( m = 0, " , N k ) denote

the basis functions. The basis functions which are local orthogonal in [2] are adopted, which are  ( K ) ( x ) = 1, Φ 0

 ( K ) ( x ) = 2 ( x − xK ) , Φ 1 Δx K

 ( K ) ( x ) = 2 ( y − yK ) , Φ 2 Δy K

 ( K ) ( x ) = 4 ( x − xK )( y − yK ) , Φ  ( K ) ( x ) = 4 ( x − xK ) − 1 , Φ  ( K ) ( x ) = 4 ( y − yK ) − 1 , Φ 3 4 5 2 2 3 3 Δx K Δ y K ( Δx K ) ( Δy K ) 2

where

( xK , y K )

2

(27)

is the center of element K ; ΔxK and ΔyK are the lengths of element K ’s

sides in the x and y directions, respectively. The semi-discrete scheme of Eq. (25) can be written as U K ,t = R ( U K ) ,

(28)

where R ( U K ) is a spatial discretization operator. We discretize Eq. (28) in time by third-order total variation diminishing Runge-Kutta scheme [59, 60]:

U(K) = U nK + ΔtR ( U nK ) , 1

3 1 1 1 2 1 U(K ) = U nK + U(K) + ΔtR U(K) , 4 4 4 1 2 2 2 2 U nK+1 = U nK + U(K ) + ΔtR U(K ) . 3 3 3

( ) ( )

(29)

If stronger discontinuities appear in the approximate solution, significant spurious oscillations or even nonlinear instability may occurs. In order to eliminate the oscillations, enhance the stability of the scheme, and suppress the failure of preserving the density or pressure positivity resulting from 11

strong shock impacting material interfaces or complex geometries, a robust high order weighted essentially non-oscillatory (WENO) limiter [61], carried out with the decomposition in a local characteristic space to improve the effectiveness, is applied after each Runge-Kutta inner stage in current works.

3.2 Gas-water interface treatment The level set function [5, 6] is adopted to capture the moving gas-water interface, given by ∂φ f ∂t

+u

∂φ f ∂x

+v

∂φ f ∂y

= 0,

(30)

where φ f is a signed distance function. In order to capture the accurate position of the moving material interface, The Eq. (30) is discretized by a fifth-order WENO scheme [62, 63] and a three-order TVD Runge-Kutta scheme [59] in space and time, respectively. When the velocity gradient at the vicinity of the interface is too large, distortion of the level set contours may occur or even the calculation would break down [8]. To improve this situation, an extension velocity field [64] is constructed for solving Eq. (30). In current work, the extension velocity is obtained by solving a convection equation to extrapolate the material interface velocity, acquired by solving interface Riemann problem, to the whole flow region as in [8]. After obtaining the location of the gas-water interface, the RGFM is used to treat the interface Riemann problem. For the sake of clarity, a one-dimensional illustration is used to illustrate the gas-water interface treating scheme. As shown in Fig. 3.1, for the real fluid cells i − 2 , i − 1 , and i , the cells i + 1 , i + 2 , and i + 3 in the other side of the interface are defined as the ghost fluid elements, and then a gas-liquid Riemann problem is defined and solved using two nonlinear

12

characteristic equations to predict the interfacial states along the interface normal direction. The obtained predicted interfacial states are then used to redefine the flow states of the ghost fluid elements i + 1 , i + 2 , i + 3 and the real fluid element i . For Euler system (1), the two nonlinear characteristic equations [10, 65] along the gas-water interface normal direction are dpI du 1 − n ρ IL u IL a IL2 = 0, along + ρ IL aIL I + dt dt xI dpI du 1 − n ρ IR u IR aIR2 = 0, along − ρ IR aIR I + dt dt xI

dxI = u I + a IL , dt dxI = u I − aIR . dt

(31) (32)

Here subscripts "I" , "IR " , and "IL" denote the interface, right and left side of the interface, respectively. pI and uI are the pressure and velocity at the interface, respectively. ρ IL , uIL , and

aIL are the density, velocity, and sound speed at the left side of the interface, respectively. ρ IR , u IR and aIR are the density, velocity, and sound speed at the right side of the interface, respectively. n is the system parameter defined in Eq. (1). An approximate Riemann solver developed in [34] and

based on double shock approximation is adopted to solve Eqs. (31) and (32). For two-dimensional cases, the interface Riemann problem is built along the normal direction of the gas-water interface. After obtaining the interfacial states and redefining the flow states of real elements adjacent to the material interface, a constant extrapolation technique in [66] is employed to extrapolate the interfacial flow status into the ghost fluid region. The constant extrapolation equation is given by

It ± n f ⋅∇I = 0,

(33)

where I is interface flow variables containing density, pressure, normal velocity, and tangential velocity; n f = ∇φ f ∇φ f

is the unit normal vector of the gas-water interface, calculated by a

central differential scheme. Once the flow status of the ghost fluid elements is defined for each

13

medium, the RKDG solver is used to solve both the gas and the liquid.

Fig. 3.1 One-dimensional schematic diagram for the gas-water interface treating scheme.

3.3 Fluid-structure interface treatment Under the assumptions that the fluid-structure interface is adiabatic, nonreactive, impermeable, and unable to support surface tension, the pressure and normal velocity across the material interface are continuous, and all coupling between the fluid and the structure occurs only at the material interface [46]. An interface treating scheme based on the modified ghost fluid method (MGFM) [34], the fluid characteristic equation [67], and the structure equation of motion is used to couple the RKDG fluid solver and the FEM structure solver at the fluid-structure interface. The structure motion equation and the fluid characteristic equation are solved simultaneously at the fluid-structure interface to obtain the interface states. The obtained interfacial states are then used to define the flow states are then used to define the flow states of the ghost fluid elements and the normal surface traction of the deformable structure. For current scheme, the kinematic compatibility and pressure equilibrium conditions are satisfied at the material interface. As shown in Fig. 3.2, Eulerian fluid and Lagrangian structural elements coexist in the computation domain. We define the region Ωf , covered by structured grids, as the Eulerian fluid domain, and the region Ωs , covered by the unstructured grids, as the Lagrangian structure domain. 14

Γfs is the fluid-structure interface, Ωr = Ω f − Ωs is the real fluid domain, and Ω g = Ω f  Ωs is the ghost fluid domain. At initial, a signed distance function φs is constructed at the center of each fluid element with φs > 0 for the ghost region and φs < 0 for the real region, and the minimum distance from the element center to material interface Γfs is defined as the magnitude of φs . In order to treat the inner fluid boundary condition and let the real fluid elements adjacent to Γfs having a complete stencil, a narrow band of ghost fluid elements with 0 < φs < ε need to be solved. The value of ε is dependent on the numerical scheme stencil which is used in fluid calculation, and it is set to 3 × max ( Δx, Δ y ) for current RKDG solver. The Eq. (30) with the material interface velocity obtained from the FEM solver is used to capture the moving fluid-structure interface in each time step. As shown in Fig. 3.2, in order to define the fluid-structure interfacial status, assuming that

A is an element on boundary located in real fluid domain, we need to search another element on boundary B located in the structure domain, and the angle made by the respective normal directions of elements A and B is minimum. The unit normal vector of each structure element is defined as n s = ∇ φ s ∇ φ s .

Fig. 3.2 Two-dimensional schematic diagram for the fluid-structure interface treating scheme.

15

As shown in Fig. 3.3, an one-dimensional illustration is used to illustrate the fluid-structure interface treating technique. We assume that the fluid domain is located on the left side of the material interface and the structure domain is on the right side. The fixed Eulerian fluid cells are distributed over both the fluid domain and the structure domain, while the moving cells of the Lagrangian structure are only located over the structure domain. The fluid cells located on the fluid and structure domains are defined as real fluid and ghost fluid elements, respectively. In each time step, only a small number of ghost fluid cells near the interface needs to be solved, and the other ghost fluid cells are existence to allow that the fluid-structure interface is always immersed in the Eulerian fluid cells even after structure deformation, so the interface can always be treated as an internal boundary. The non-linear fluid characteristic equation for the fluid on the left side of the interface [10] and the equation of motion for the right structure are solved simultaneously at the fluid-structure interface, given by dpI du 1 − n ρ IL u IL a IL2 = 0, along + ρ IL aIL I + dt dt xI

³

Ωs

dxI = u I + a IL , dt

s = − ³ ∇NsdV + ³ NfdV + ³ Ntd Γ, ρ s NT NdV u Ωs

Ωs

Γt

(34) (35)

where ρIL , uIL , and aIL are the pressure, velocity, and density of the fluid on the interface’s left side, respectively; ρ I , uI , and xI denote the pressure, velocity, and location at the interface, respectively. In each time step, the acoustic impedance ρILaIL can be assumed to be constant [9, 68], and thus Eq. (34) can be discretized as

pI = pIL − ρIL aIL ( uI − uIL ) −

(1 − n) Δt ρ xI

2 IL IL IL

u a ,

(36)

where the interface velocity uI and pressure pI are obtained through solving Eqs. (35) and (36)

16

simultaneously, leading to a half-way coupled scheme [9, 10]. Then pI is used to defined the surface traction boundary condition for the structure region, and pI , uI are the used to define the variables of the flow for the ghost fluid region, as shown in Fig. 3.3. For two-dimensional cases, this interface treating scheme operates on the material interface’s normal direction, and the Eikonal Eq. (33) is used to extrapolate the flow status from the interface to whole ghost fluid region.

Fig. 3.3 One-dimensional schematic diagram for the fluid-structure interface treating scheme.

3.4 Outline of the combined approach To make sure the stability of current combined RKDG-FEM approach, the global time step size should be determined by [9, 10]

Δt = min ( Δt f , Δts ) ,

(37)

where

§ Δxi , j § Δyi, j · ρs · ¸ , Δts = min ¨ Le Δt f = CFL × min ¨ , . ¨ ˆ + 2μˆ ¸¸ i, j ¨ u ¸ a v a + + λ , , , , i j i j i j i j © ¹ © ¹ Here Δt f is the time step size for the fluid calculation, ai, j is the sound speed of fluid, and CFL is the Courant number, which is set to 0.2 for all calculations in current work. Δts is the time step size for the structure calculation, ρ s is the structure density, Le is the characteristic structure

17

element dimension, Ee is Yong’s modulus, ν is Poisson’s ratio, and λˆ = E eν ª¬ (1 + ν )(1 − 2ν ) º¼ and μˆ = Ee 2 (1 + ν ) are the Effective Lame’s constants. Now we briefly summarize the algorithm implementation of the combined RKDG-FEM approach. As illustrated in Fig. 3.4, it contains two computational module: the multiphase module and the FSI module. The core of the multiphase flow module is a gas-water interface treating scheme based on the LSM and the RGFM, and that of the FSI module is a fluid-structure interface treating scheme based on the MGFM, the fluid characteristic equation, and the structure’s motion equation. We assume the solutions at t = t n are known, and the following steps are taken to obtain the solutions at next time step: Step 1: Calculate the global time step size by Eq. (37), which satisfies the stability conditions over the fluid and structure computational domains. Step 2: Advance the fluid level set variable φ f and the structure level set variable φs to the next time step via Eq. (30). Then the new locations of the gas-water and fluid-structure interfaces can be defined. Step 3: Define the ghost fluid elements adjacent to the gas-water interface as shown in Fig. 3.1, calculate the gas-water interface status using Eqs. (31) and (32), and then project the variables of the interface into the ghost fluid elements via Eq. (33). Step 4: Define the ghost fluid elements adjacent to the fluid-structure interface as shown in Fig. 3.3, calculate the fluid-structure interface status using Eqs. (34) and (35), and project the interface variables into ghost fluid elements by the same way as step 3. Step 5: Use the RKDG solver to calculate the gas and water domains, and the FEM solver to

18

calculate the structure domain. Step 6: Check whether the stop criterion is satisfied. If not, repeat Step 1 to Step 5.

Fig. 3.4 Outline of the combined RKDG-FEM approach for simulating UNDEX near composite structures.

4 Validation studies 4.1 One-dimensional underwater explosion problem In order to validate the multiphase flow module of current combined RKDG-FEM approach, a one-dimensional UNDEX case with extremely high pressure in the gas medium is studied, which has also been investigated in [69, 70]. The initial conditions are

­(1630, 0, 7.81×109 ,1.4 ) , for 0 ≤ x ≤ 0.5, ° ( ρ , u , p, γ ) = ® 5 °¯ (1000, 0,1×10 , 7.15 ) , for 0.5 < x ≤ 1.

(38)

As the initial high-pressure region is located on the gas side, the exact solution includes a rarefaction wave in left direction, a contact discontinuities in right direction, and a strong shock wave in right direction. Simulations are calculated using 400 cells. The computed density, pressure, and

19

velocity at t = 0.0001 are shown in Fig. 4.1, which dovetail nicely with the respective analytical results in [69, 70].

20

Fig. 4.1 Comparisons of the density, pressure, and velocity profiles for the one-dimensional UNDEX problem at

t = 0.0001 . Black line: the analytical solution in [69, 70]. Red square: the numerical solution computed by the present combined RKDG-FEM approach.

4.2 One-dimensional cavitation problem in a closed tube For the accuracy and stability validation of the isentropic cavitation model, an one-dimensional cavitation problem in a closed tube is modelled, which has also been studied by Jun Zhu in [71]. The initial conditions are ­(1000, − 100,10 5 , 7.15 ) , for 0 ≤ x ≤ 0.5, ( ρ , u , p , γ ) = °® 5 °¯(1000, + 100,10 , 7.15 ) , for 0.5 < x ≤ 1.0.

(39)

The perfectly reflective boundary condition is imposed on both ends of the tube, and the two shocks reflected from the ends move towards the center, resulting in the interactions of shocks and cavitation. Calculations are performed using 400 cells. The computed density, pressure, and velocity at t = 0.0001 are shown in Fig. 4.2, which agree well with the respective numerical solutions calculated by Jun Zhu in [71].

21

22

Fig. 4.2 Comparisons of the density, pressure, and velocity profiles for the one-dimensional cavitation problem in a closed tube at t = 0.0001 . Black line: the numerical solution in [71]. Red square: the numerical solution computed by the present combined RKDG-FEM approach.

4.3 One-dimensional underwater shockwave impact on a composite plate Following, in order to verify the effectiveness of the RKDG-FEM approach in solving the problems involving composite sandwich structures, an example of one-dimensional underwater shock impact on a three-layers sandwich plate is computed, which has also been studied by experimental and numerical approaches in [9, 22]. As shown in Fig. 4.3, the water domain has a length 1.4 m , density 1000 kg/m3 , and sound speed 1400 m/s , and is discretized by 2000 uniform −t θ is imposed on the grids. A pressure wave satisfying exponentially decaying function p ( t ) = pm e

water domain left side; pm is the peak value of the pressure wave, and θ is the decay time. A sandwich plate is placed on the right end of the water region, which comprises two identical face plates of thickness 12 mm and a foam core plate of thickness 25 mm . The foam core is discretized by 30 uniform grids, and the face sheets are both discretized by 6 uniform grids The steel face sheets are modelled as elastic solids with density 8000 kg/m 3 and Young’s modulus 210 GPa . An

23

elastic-plastic rate dependent solid model is used to simulate the foam core plate with density 236 kg/m3 and Young’s modulus 1 GPa . Free-sliding (unsupported) or rigid fixed (supported) boundary condition is imposed on the back of the sandwich plate. Fig. 4.4 shows the plastic core compression strain of the sandwich plate at different shock impulse with θ = 0.24 ms of both unsupported and supported boundary conditions, which shows good agreement with the numerical results in [9, 22] and the experimental result in [22].

Fig. 4.3 Schematic diagram for the one-dimensional underwater shock impact on a sandwich plate.

Fig. 4.4 Plastic core compressions of sandwich plates with unsupported and supported boundary conditions ( θ = 0.24 ms ). Red square: the experimental solutions in [22]. Black line: the numerical solution in [22]. Green circle: the numerical solution in [9]. Blue triangle: the numerical solution computed by the present combined RKDG-FEM approach.

24

4.4 Two-dimensional underwater explosion in a free field Next, we compute an example of two-dimensional UNDEX in a free field, which has also been studied in [72-74] to verify the effectiveness of different compressible multi-phase flow algorithms. As shown in Fig. 4.5, the size of the computational domain is

( x, y ) ∈ [ −1,1] × [ −1,1] ,

which is

discretized by 400 × 400 uniform structured grids. Outflow boundary condition is imposed on all four sides of the fluid domain. At initial, the bubble with a radius of 0.5 is located at center of the fluid domain. The initial states are ­ (1000, 0, 0,1 × 10 5 , 7.15 ) , for water, ( ρ , u , v , p , γ ) = °® 9 °¯(1630, 0, 0, 7.81 × 10 ,1.4 ) , for bubble.

(40)

Fig. 4.6 shows the density (left) and pressure (right) contours of the flow region at t = 0.0001 , where the dotted lines present the interface of the gas bubble. Fig. 4.7 shows the comparison of the density, pressure, and velocity curves along x axis at t = 0.0001 between two-dimensional numerical results (red square) and the “exact” solutions computed by an one-dimensional axisymmetric RKDG-GFM model (black line) with 5000 cells as treating in [73, 74]. The present numerical solutions coincide with the “exact” solutions, which verifies the validity of the multiphase flow module of current RKDG-FEM approach in solving two-dimensional UNDEX problems.

25

Fig. 4.5 Schematic diagram for the two-dimensional UNDEX in a free field.

Fig. 4.6 The density (left) and pressure (right) contours for the two-dimensional UNDEX in a free field at

t = 0.0001 .

26

Fig. 4.7 A comparison of the density, velocity, and pressure profiles for the two-dimensional UNDEX in a free field at t = 0.0001 . Black line: the numerical solution calculated by an axisymmetric RKDG model in one-dimension with 5000 cells. Red square: the numerical solution computed by the present combined RKDG-FEM approach.

4.5 Two-dimensional underwater shock impact on a composite plate Finally, in order to verify current RKDG-FEM approach in solving more realistic problems involving composite structures, a two-dimensional case of underwater shock impact on a composite plate with three-layers is studied, which has also been studied by a L-L and an E-L approach in [10]. As shown in Fig. 4.8, the sandwich plate consists of two steel faces and an Aluminum foam core.

27

The thickness of the steel faces is 0.08 m , and that of the Aluminum foam is 0.14 m . We assume perfect bonding is performed between different material layers. An exponentially decaying pressure wave satisfying p = p0 e( pressure

− y yd )

is imposed on the fluid domain lower boundary, and the peak

p0 and the characteristic decay length yd are set to 2.0 × 106 Pa and 1.4 m ,

respectively. P1 is a predefined gauging point located at the center of the lower steel face. The steel 3 plate’s density, Young’s modulus, Poisson’s rate and yield strength are chosen as ρ s = 8000 kg/m ,

Es = 2.1×1011 Pa , μs = 0.27 , and σ s = 9.0 ×108 Pa , respectively; the Aluminum foam core plate’s density, Young’s modulus, Poisson’s rate and yield strength are chosen as ρ a = 2700 kg/m3 ,

Ea = 1.0 ×1010 Pa , μa = 0.0 , and σ a = 1.5 ×106 Pa , respectively. After performing a mesh convergence study, Δx = Δy = 0.01 m and Δx = Δy = 0.02 m are determined for the fluid and the structure grids, respectively. For the fluid domain, out-flow boundary condition is applied on the left, right and bottom boundaries. For the deformable structure, the left and right sides are assumed to be rigid fixed. The time step size is set to Δt = 5×10−7 s . Fig. 4.9 shows the pressure and displacement histories of the gauging point P1, comparing with the results calculated by the L-L and E-L approaches in [10]. The numerical results coincide well with the results computed by E-L approach in [10], which verifies the effectiveness of the combined RKDG-FEM approach in solving fluid-structure interaction problems involving composite structures.

28

Fig. 4.8 Schematic diagram for the two-dimensional underwater shock impact on a water-backed, end-clamped composite plate.

Fig. 4.9 Pressure and displacement histories of point P1 for the shock impact on a water-backed, end-clamed sandwich plate. Blue line: the numerical solution calculated by a L-L approach in [10]. Black line: the numerical solution calculated by an E-L approach in [10]. Red square: the numerical solution computed by the present combined RKDG-FEM approach.

5 Application 5.1 Computational set-up In order to investigate the flow characteristics and structure response of the steel and sandwich structures subjected to UNDEX, two distinct configurations described in Fig. 5.1 are simulated in

29

this work: Case 1: UNDEX near a submerged, end-clamped circular steel plate. Case 2: UNDEX near a submerged, end-clamped circular sandwich plate. The steel plate in case 1 has the same mass as the sandwich plate in case 2. Axial-symmetric numerical model is adopted in both cases. For simplicity, we assume the bonding between different material layers of the sandwich plate is perfect, and the interface delamination between the material layers is ignored. The steel plate and the steel facesheets of the sandwich plate’s density, Young’s 5 modulus and Poisson’s ratio are ρ s = 8000 kg/m 3 , Es = 2.1×10 MPa , and ν s = 0.28 [75], and

other material parameters are A=950 MPa , B =560 MPa , n =0.26 , C =0.014 , and m = 1 [75]. The aluminum foam core of the sandwich plate’s density, Young’s modulus, elastic Poisson’s ratio, 3 3 plastic Poisson’s ratio are ρ f = 430 kg/m , E f = 1×10 MPa , ν e = 0.3 and ν p = 0 [19, 26]. We

follow the experimental data from the reference [26] for the relationship between static yield strength

σ c and equivalent plastic strain εˆ p . As shown in Fig. 5.1, WSI1 and WSI2 are the lower and upper water-structure interfaces, respectively, and P1 and P2 are the predicted Lagrangian gauging points located at the center of WSI1 and WSI2, respectively. For both cases, the horizontal dimension of the fluid domain is set to 5.0 m , and vertical dimension is set to 12.0 m . The explosive is replaced by a high-pressure gas

bubble as in [10, 76], and the initial radius and pressure of the bubble are set to 0.5 m and 6.0 GPa , respectively. At initial, the distance between the bubble center and WSI1 is 2.0 m . The radii of the steel and the sandwich circular plate are both taken to be 5.0 m . The thickness of the steel plate is 0.068 m , that of the sandwich plate’s facesheets is 0.03 m , and that of the foam core

30

is 0.15 m . For the fluid domain, the right, down, and up boundaries is applied with the non-reflective boundary condition, and the left boundary is applied with the axial-symmetric boundary condition. For the deformable structure, the left and right sides is applied with the axial-symmetric and rigid fixed boundary conditions, respectively. After performing a grid convergence analysis, Δx = Δy = 0.005 m and Δx = Δy = 0.01 m are used for fluid and structure cells, respectively, and the global time step is set to Δt = 8 ×10−7 s .

Fig. 5.1 Schematic diagrams for UNDEXs near a steel plate (left: case 1) and a sandwich plate (right: case 2).

5.2 Flow characteristics With the computational set-up, simulations are carried out through the combined RKDG-FEM approach developed in h3. In this subsection, we focus on the analysis of the flow characteristics of the both cases. Fig. 5.2 shows the numerical schlieren (left) and pressure distribution (right) contours of the flow domain of case 1 at different time instances. Shock propagation, bubble expansion, cavitation evolution, and complex interaction among shock, bubble, cavitation region, and structures are presented clearly. After detonation, an incident shock wave (SW1) is outdiffusion from the 31

high-pressure bubble, and a rarefaction wave is propagation into the bubble simultaneously, as shown in Fig. 5.2(a). With SW1 reaching and then interacting with WSI1, as the acoustic impedance of water is lower than that of steel, another shock wave (SW2) is reflected back to the flow region below the steel plate, as shown in Fig. 5.2(b). Under the combined impact of SW1 and SW2, the steel plate obtains a transient velocity and immediately generates an obvious upward displacement. Therefore, a shock wave is transmitted to the fluid region above WSI2, and a rarefaction wave (RW1) is reflected back to the fluid region beneath WSI1, as shown in Fig. 5.2(c). After SW2 reaches the top surface of the bubble, due to the lower acoustic impedance than that of water, another rarefaction wave (RW2) is reflected back to water domain above the bubble, and a shock wave (SW3) transmits into the bubble simultaneously, as shown in Fig. 5.2(d). When RW2 encounters RW1 at water region between the plate and the bubble, the fluid pressure of this overlap region is lower than the saturated pressure, and then a flow cavitation bubble occurs, as shown in Fig. 5.2(e), (f), and (g). As illustrated in Fig. 5.2(g), under the combined impact of the complex wave structures between the original gas bubble and the plate, the cavitation bubble collapses from the middle position, and then a cavitation ring is presented. As the propagation speed of shock wave in water is larger than that in air, SW2 propagate faster than SW3, which can be observed in Fig. 5.2(g) and (h). Thus a diffraction wave generates at the bottom of the bubble, as plotted in Fig. 5.2(i).

32

Fig. 5.2 Numerical schlieren (left half) and pressure (right half) contours of case 1: (a) t = 0.08 ms ; (b)

t = 0.72 ms ; (c) t = 0.92 ms ; (d) t = 1.08 ms ; (e) t = 1.16 ms ; (f) t = 1.48 ms ; (g) t = 1.76 ms ; (h) t = 2.16 ms , (i) t = 3.16 ms .

Fig. 5.3 shows the numerical schlieren (left) and pressure distribution (right) contours of the flow domain of case 2. As illustrated in Fig. 5.3(b), the sandwich plate is under the combined compact of the incident shock wave SW1 and the reflected shock wave SW2 as in case 1. The velocity and displacement of the lower facesheet are obviously larger than those of the steel plate in case 1 respectively, and the upper facesheet is still static as the momentum translated from the lower facesheet is entirely absorbed by the aluminum foam core through plastic dissipation. As the large transient velocity and displacement of the lower facesheet, a strong rarefaction wave (RW1) is reflected back to the flow region, resulting in the inception of a cavitation bubble beneath the sandwich plate immediately, as presented in Fig. 5.3(c). As shown in Fig. 5.3(d), another tiny cavitation bubble appears right above the original gas bubble, resulting from the rarefaction wave (RW2) reflected from the bubble surface. As shown in Fig. 5.3(d), after meeting, the two cavitation

33

bubbles merge into a single one that is larger and more lasting than the cavitation bubble in case 1. Comparisons between the two cases show that the foam core has an obvious effect in the FSI and flow cavitation effects, and this effects can be utilized to improve the anti-blast capability of composite sandwich structures.

Fig. 5.3 Numerical schlieren (left half) and pressure (right half) contours of case 2: (a) t = 0.08 ms ; (b)

t = 0.72 ms ; (c) t = 0.92 ms ; (d) t = 1.08 ms ; (e) t = 1.16 ms ; (f) t = 1.48 ms ; (g) t = 1.76 ms ; (h) t = 2.16 ms , (i) t = 3.16 ms .

Fig. 5.4 shows the pressure-time histories of the predefined gauging points in both cases: P1 of case 1 (C1P1), P2 of case 1 (C1P2), P1 of case 2 (C2P1), and P2 of case 2 (C2P2). The first peak pressures of C1P1 and C2P1 appear simultaneously at t ≈ 0.54 ms , caused by the impact of the

34

primary incident shock wave SW1. As the velocity and displacement of the sandwich plate’s lower facesheet are larger than those of the steel plate respectively, the peak pressure of C2P1 is slightly smaller than that of C1P1, and attenuates more quickly. At t ≈ 0.95 ms , the pressure of C2P1 reduces to the valley value, and starts to reload again. The second peak pressure of C2P1 appears at t ≈ 1.2 ms , which is derived from the collapse of the cavitation bubble as shown in Fig. 5.3(d), (e),

and (f). There is only one pressure peak for C1P1, and it decays to the minimum value at t ≈ 1.5 ms . For the upper facesheet of the sandwich plate, as part of the momentum transmitted from the lower facesheet is dissipation by the central foam core, the peak pressure of C2P2 appears obviously later than that of C1P2, and it is slightly smaller. Comparison of the fluid pressure profiles along WSI1 between case 1 and case 2 at different time instances is shown in Fig. 5.5. For both cases, the first fluid peak pressure results from the combined impact of the initial incident shock wave (SW1) and the reflected shock wave (SW2). In each time instance, the first peak pressure of case 2 is clearly lower than that of case 1, and it decays more quickly, which indicate the sandwich plate is excellent for underwater explosion loading resistance. At t ≈ 1.2 ms , another pressure peak arises at middle position of the sandwich plate for case 2, and it stems from the collapse of the cavitation bubble as shown in Fig. 5.3(d), (e), and (f). The second peak pressure is non-negligible comparing to the first one. At t ≈ 1.48 ms and t ≈ 1.8 ms , the second pressure peak moves towards the ends of structure as the expansion of the

cavitation bubble ring, as shown in Fig. 5.3(d), (f), (g), and (h). For case 1, however, the pressure distribution on WSI1 is scarcely influence by the cavitation effect. This analysis shows the sandwich plate changes the load characteristic of the flow region through affecting the FSI and flow cavitation

35

effects; these effects can to be utilized to improve blast mitigation of composite sandwich structures.

Fig. 5.4 Comparison of the pressure histories of the predefined gaging points between cases 1 and 2.

Fig. 5.5 Comparison of the fluid pressure profiles along WSI1 at different time instances between cases 1 and 2.

5.3 Structural response In this subsection, we focus on the analysis of the dynamic responses of the monolithic steel plate in case 1 and the aluminum foam-cored sandwich plate in case 2. Fig. 5.6 shows the Mises stress and PEEQ contours of the steel plate (left) and the sandwich plate (right) at different time instances. Notice that the vertical dimensions of all nephograms are magnified by three times for 36

clarity. For case 1, the steel plate is deformed continuously along with the propagation of stress wave to both sides. At each time instance, the Mises stress of the upper surface of the steel plate is obviously larger than that of the lower surface. As shown in Fig. 5.6(e) and (f), the PEEQ wave appears at the central position of the steel plate. For case 2, the deformation process of every cross-section along z-axis direction of the sandwich plate can be divided into three distinctive phases [77-79]: Phase 1: The lower facesheet obtains an initial momentum from the combined impact of SW1 and SW2 as shown in Fig. 5.3(b) and (c), and starts to compress the foam core, while the upper facesheet is still stationary. Phase 2: The upper facesheet starts to deform till all components of the sandwich plate acquire an identical velocity when the foam core is fully compacted. Phase 3: The lower facesheet, the foam core, and the upper facesheet of the sandwich plate bend and stretch almost perfectly in step. As shown in Fig. 5.6, the PEEQ shock in the foam core is obvious at every time instance, while that is nearly invisible in the lower and the upper facesheets.

37

Fig. 5.6 Mises stress (up) and PEEQ (down) contours of the steel plate (left: case 1) and the sandwich plate (right: case 2) at different time instances: (a) t = 0.8 ms ; (b) t = 1.0 ms ; (c) t = 1.5 ms ; (d) t = 2.0 ms ; (e)

t = 2.5 ms ; (f) t = 3.0 ms . (Notice that the vertical dimensions are magnified by three times.)

Comparisons of the displacement and velocity histories of the predefined gauging points between cases 1 and 2 are shown in Figs. 5.7 and 5.8, respectively. For the sandwich plate’s centre cross-section (P1-P2), the time intervals of three deformation phases are approximately 0.55  0.88 ms , 0.88  1.63 ms and 1.5  3.0 ms , respectively. As illustrated in Fig. 5.7, for case 1,

as the compressive deformation of the steel plate is extremely small, the displacement and velocity of C1P1 coincide closely with those of C1P2, respectively. For case 2, however, neither the displacements nor the velocities of C2P1 and C2P2 are coincidence as the deformation processes of the sandwich plate’s components are inconsistent, and the displacement of the front facesheet (C2P1) is far larger than that of the back facesheet (C2P2) in all three phases. The displacement curves of the

38

steel plate (C1P1 and C1P2) are always located between those of the front and back facesheets. As depicted in Fig. 5.7, in phase 1, the monolithic steel plate (C1P1 and C1P2) and the lower facesheet of the sandwich plate (C2P1) both obtain a transient velocity resulting from the joint impact of the initial incident shock wave (SW1) and the reflected shock wave (SW2), and the peak velocity of the monolithic steel plate is about half of that of the front facesheet. However, the upper facesheet of the sandwich plate (C2P2) is always motionless in this phase. In phase 2, the upper facesheet (C2P2) starts to move, while the speed of the lower facesheet (C2P1) is decreasing uniformly. The lower and upper facesheets basically reach an identical moving velocity at the end of phase 2, and move together during the whole process of phase 3. Comparison of the displacements between C1P2 and C2P2 indicates that composite sandwich structures provide significant benefits for blast resistance.

Fig. 5.7 Comparison of the displacement histories of the predefined gaging points between cases 1 and 2.

39

Fig. 5.8 Comparison of the velocity histories of the predefined gaging points between cases 1 and 2.

Comparison of the PEEQ histories of the predefined gauging points between cases 1 and 2 is shown in Fig. 5.9, and the PEEQ is only presented in phase 3 for both cases. The steel plate is damaged more severely than the sandwich plate as the PEEQs of C1P1 and C1P2 are always larger than those of C2P1 and C2P2, respectively. Comparison of the energy histories between the monolithic steel in case 1 and the sandwich plate in case 2 is shown in Fig. 5.10. For the steel plate, almost all external work induced by the explosion loads is transformed into the kinetic energy of the structure, while the internal energy and plastic dissipation are negligible. However, for the sandwich plate, the aluminum foam core’s plastic deformation absorbs a large amount of external work, and only a small fraction of external work is stored during the deformation of the form core as the elastic energy. As can be observed in Fig. 5.10, the external work of case 2 is much larger than that of case 1 during all three deformation phases. This may due to the fact that the cavitation bubble in case 2 is much larger and more durable than that in case 1, which significantly changes the loading

40

characteristic of the flow domain. This study shows the sandwich plate has a higher underwater explosion load resistance through absorb plastic energy by foam cored than the monolithic steel plate of equivalent mass.

Fig. 5.9 Comparison of the PEEQ histories of the predefined gaging points between cases 1 and 2.

Fig. 5.10 Energy balance for the monolithic steel plate in case 1 and the foam-cored sandwich plate in case 2.

6 Conclusions In this paper, a combined RKDG-FEM approach is developed to simulating near-field

41

early-time UNDEX involving compressible multiphase flows, transient fluid-structure interaction, and flow cavitation. This approach is composed of a RKDG solver for fluid simulating, a FEM solver for structure simulating, the isentropic one-fluid cavitation model for transient cavitation capturing, the LSM for moving gas-water and fluid-structure interfaces capturing, the RGFM for gas-water interface coupling, and an GFM-based interface treating scheme for fluid-structure interface coupling. The RKDG solver is applied to solve the Euler equation for inviscid and compressible flow on an Eulerian framework which uses a fixed, Cartesian grid, and the FEM solver to the equations of motion and material constitution of deformation structure on a Lagrangian framework using a moving grid. The fluid-structure interface treating scheme is accomplished by solving simultaneously the fluid characteristic equation and the structure equation of motion at the material interface to satisfy the kinematic compatibility and pressure equilibrium conditions. Five benchmark cases are studied to check the accuracy and efficiency of the coupling RKDG-FEM approach. Then the flow characteristics and structure response are investigated for UNDEXs near a submerged aluminum foam-cored sandwich plate and a monolithic steel plate of equivalent mass. The simulating results show the aluminum foam core of the sandwich plate has an effective effect in the FSI and flow cavitation effects. Flow characteristics between the cases with sandwich plate (case 2) and with a single steel plate of equivalent mass (case 1) are distinguishing. Only one cavitation bubble is generated in case 1, but two in case 2, and the fused cavitation bubble in case 2 is larger and more durable than the single cavitation bubble in case 1. Phenomena of reflection, refraction, transmission, and diffraction of shock wave are present in both cases, while the wave system between the bubble and the sandwich plate in case 2 is much more complex than that

42

between the bubble and the monolithic steel plate in case 1. Dynamic responses of the sandwich and the steel plate are also difference. The deflection and PEEQ of the sandwich plate’s upper facesheet are significantly lower than those of the monolithic steel plate of equivalent mass, respectively, and the impulse transmitting through the sandwich plate is smaller than through the monolithic steel plate, which indicate that the sandwich plate has a higher UNDEX load resistance through absorbing plastic energy by the foam cored than the monolithic steel plate. Current work mainly focus on the algorithm, verification, and simple application of the combined RKDG-FEM approach. Next, it will be use to investigate more complex problems involving material heterogeneity, damage mechanisms, and complex failure modes of composite sandwich structure when subjected to underwater explosion loads, and to exploit the FSI and flow cavitation effects of UNDEX near sandwich structure so as to improve the UNDEX mitigation capability of composite structures.

Acknowledgements This work was supported by the National Natural Science Foundation of China (No: 51879052) and the Industrial Technology Development Programs (JCKY2017604C002, JCKY2018604C010). The first author also wants to thank Dr. Zhongyu Zhang and Dr. Xiao Huang for helpful discussions.

References [1] B. Cockburn, C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework, Math. Comp., 52 (1989) 411--435. [2] B. Cockburn, C.-W. Shu, The Runge–Kutta Discontinuous Galerkin Method for Conservation Laws V: Multidimensional Systems, Journal of Computational Physics, 141 (1998) 199-224. [3] O.C. Zienkiewicz, R.L. Taylor, O.C. Zienkiewicz, R.L. Taylor, The finite element method, McGraw-hill London, 1977. [4] T. Belytschko, W.K. Liu, B. Moran, K. Elkhodary, Nonlinear finite elements for continua and structures, John 43

wiley & sons, 2013. [5] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations, Journal of Computational Physics, 79 (1988) 12-49. [6] M. Sussman, P. Smereka, S. Osher, A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow, Journal of Computational Physics, 114 (1994) 146-159. [7] R.P. Fedkiw, T. Aslam, B. Merriman, S. Osher, A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method), Journal of Computational Physics, 152 (1999) 457-492. [8] C.W. Wang, T.G. Liu, B.C. Khoo, A Real Ghost Fluid Method for the Simulation of Multimedium Compressible Flow, SIAM Journal on Scientific Computing, 28 (2006) 278-302. [9] Z. Liu, W. Xie, Y.L. Young, Numerical modeling of complex interactions between underwater shocks and composite structures, Computational Mechanics, 43 (2008) 239-251. [10] W. Xie, Z. Liu, Y.L. Young, Application of a coupled Eulerian-Lagrangian method to simulate interactions between deformable composite structures and compressible multiphase flow, International Journal for Numerical Methods in Engineering, 80 (2009) 1497-1519. [11] T.G. Liu, B.C. Khoo, W.F. Xie, Isentropic one-fluid modelling of unsteady cavitating flow, Journal of Computational Physics, 201 (2004) 80-108. [12] S. Avachat, M. Zhou, Response of submerged metallic sandwich structures to underwater impulsive loads, Journal of Mechanics of Materials and Structures, 10 (2015) 17-41. [13] S. Avachat, M. Zhou, Compressive response of sandwich plates to water-based impulsive loading, International Journal of Impact Engineering, 93 (2016) 196-210. [14] R.H. Cole, R. Weller, Underwater explosions, Physics Today, 1 (1948) 35. [15] A. Keil, The response of ships to underwater explosions, in, DTIC Document, 1961. [16] Z. Zhang, L. Wang, V.V. Silberschmidt, Damage response of steel plate to underwater explosion: Effect of shaped charge liner, International Journal of Impact Engineering, 103 (2017) 38-49. [17] F.R. Ming, A.M. Zhang, Y.Z. Xue, S.P. Wang, Damage characteristics of ship structures subjected to shockwaves of underwater contact explosions, Ocean Engineering, 117 (2016) 359-382. [18] P. Wang, A.M. Zhang, F. Ming, P. Sun, H. Cheng, A novel non-reflecting boundary condition for fluid dynamics solved by smoothed particle hydrodynamics, Journal of Fluid Mechanics, 860 (2018) 81-114. [19] M.F. Ashby, T. Evans, N.A. Fleck, J. Hutchinson, H. Wadley, L. Gibson, Metal foams: a design guide, Elsevier, 2000. [20] Z. Zhang, L. Wang, V.V. Silberschmidt, S. Wang, SPH-FEM simulation of shaped-charge jet penetration into double hull: A comparison study for steel and SPS, Composite Structures, 155 (2016) 135-144. [21] V. Deshpande, N. Fleck, One-dimensional response of sandwich plates to underwater shock loading, Journal of the Mechanics and Physics of Solids, 53 (2005) 2347-2383. [22] N.A. Fleck, A. Heaver, V.S. Deshpande, An underwater shock simulator, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 462 (2006) 1021-1041. [23] M.T. Tilbrook, V.S. Deshpande, N.A. Fleck, Underwater blast loading of sandwich beams: Regimes of behaviour, International Journal of Solids and Structures, 46 (2009) 3209-3221. [24] X. Wei, P. Tran, A. de Vaucorbeil, R.B. Ramaswamy, F. Latourte, H.D. Espinosa, Three-dimensional numerical modeling of composite panels subjected to underwater blast, Journal of the Mechanics and Physics of Solids, 61 (2013) 1319-1336. [25] X. Qiu, V.S. Deshpande, N.A. Fleck, Dynamic Response of a Clamped Circular Sandwich Plate Subject to

44

Shock Loading, Journal of Applied Mechanics, 71 (2004) 637. [26] D.D. Radford, G.J. McShane, V.S. Deshpande, N.A. Fleck, The response of clamped sandwich plates with metallic foam cores to simulated blast loading, International Journal of Solids and Structures, 43 (2006) 2243-2259. [27] M.T. Tilbrook, V.S. Deshpande, N.A. Fleck, The impulsive response of sandwich beams: Analytical and numerical investigation of regimes of behaviour, Journal of the Mechanics and Physics of Solids, 54 (2006) 2242-2280. [28] Y. Liang, A.V. Spuskanyuk, S.E. Flores, D.R. Hayhurst, J.W. Hutchinson, R.M. McMeeking, A.G. Evans, The Response of Metallic Sandwich Panels to Water Blast, Journal of Applied Mechanics, 74 (2007) 81. [29] Z. Xue, J.W. Hutchinson, Preliminary assessment of sandwich plates subject to blast loads, International Journal of Mechanical Sciences, 45 (2003) 687-705. [30] Z. Xue, J.W. Hutchinson, A comparative study of impulse-resistant metal sandwich plates, International Journal of Impact Engineering, 30 (2004) 1283-1305. [31] G. Taylor, The pressure and impulse of submarine explosion waves on plates, The scientific papers of GI Taylor, 3 (1963) 287-303. [32] F. Latourte, D. Grégoire, D. Zenkert, X. Wei, H.D. Espinosa, Failure mechanisms in composite panels subjected to underwater impulsive loads, Journal of the Mechanics and Physics of Solids, 59 (2011) 1623-1646. [33] N. Kambouchev, R. Radovitzky, L. Noels, Fluid–Structure Interaction Effects in the Dynamic Response of Free-Standing Plates to Uniform Shock Loading, Journal of Applied Mechanics, 74 (2006) 1042-1045. [34] T.G. Liu, B.C. Khoo, K.S. Yeo, Ghost fluid method for strong shock impacting on material interface, Journal of Computational Physics, 190 (2003) 651-681. [35] R.P. Fedkiw, Coupling an Eulerian Fluid Calculation to a Lagrangian Solid Calculation with the Ghost Fluid Method, Journal of Computational Physics, 175 (2002) 200-224. [36] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics, 39 (1981) 201-225. [37] D.J. Benson, Computational methods in Lagrangian and Eulerian hydrocodes, Computer Methods in Applied Mechanics and Engineering, 99 (1992) 235-394. [38] D.J. Benson, A new two-dimensional flux-limited shock viscosity for impact calculations, Computer Methods in Applied Mechanics and Engineering, 93 (1991) 39-95. [39] C.W. Hirt, A.A. Amsden, J.L. Cook, An arbitrary Lagrangian-Eulerian computing method for all flow speeds, Journal of Computational Physics, 14 (1974) 227-253. [40] C.S. Peskin, The immersed boundary method, Acta numerica, 11 (2002) 479-517. [41] J. Yang, S. Preidikman, E. Balaras, A strongly coupled, embedded-boundary method for fluid–structure interactions of elastically mounted rigid bodies, Journal of Fluids and Structures, 24 (2008) 167-182. [42] X. Ye, Study on transient noise and near-field load characteristics induced by bubble motion [D], in, Harbin Engineering University, 2014. [43] Z. Jin, C. Yin, Y. Chen, H. Hua, Coupling Runge-Kutta discontinuous Galerkin method to finite element method for compressible multi-phase flow interacting with a deformable sandwich structure, Ocean Engineering, 130 (2017) 597-610. [44] A. Pishevar, R. Amirifar, An adaptive ALE method for underwater explosion simulations including cavitation, Shock waves, 20 (2010) 425-439.

45

[45] P.A. Thompson, Compressible-fluid dynamic, McGraw-Hill, 1971. [46] M. Arienti, P. Hung, E. Morano, J.E. Shepherd, A level set approach to Eulerian–Lagrangian coupling, Journal of Computational Physics, 185 (2003) 213-251. [47] F.H. Harlow, A.A. Amsden, Fluid dynamics-an introductory test, in, DTIC Document, 1970. [48] T.G. Liu, B.C. Khoo, C.W. Wang, The ghost fluid method for compressible gas–water simulation, Journal of Computational Physics, 204 (2005) 193-221. [49] C.H. Cooke, T.-J. Chen, Continuous front tracking with subcell resolution, Journal of Scientific Computing, 6 (1991) 269-282. [50] A.B. Wardlaw Jr, H.U. Mair, Spherical solutions of an underwater explosion bubble, Shock and Vibration, 5 (1998) 89-102. [51] G.B. Wallis, One-dimensional two-phase flow, McGraw-Hill Companies, 1969. [52] C.E. Brennen, Cavitation and bubble dynamics, Cambridge University Press, 2013. [53] Hibbitt, Karlsson, Sorensen, ABAQUS/Explicit: User's Manual, Hibbitt, Karlsson and Sorenson Incorporated, 2001. [54] G.R. Johnson, W.H. Cook, Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures, Engineering Fracture Mechanics, 21 (1985) 31-48. [55] G.R. Johnson, W.H. Cook, A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures, in:

Proceedings of the 7th International Symposium on Ballistics, The Hague, The

Netherlands, 1983, pp. 541-547. [56] V.S. Deshpande, N.A. Fleck, Isotropic constitutive models for metallic foams, Journal of the Mechanics and Physics of Solids, 48 (2000) 1253-1283. [57] E.F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL-Riemann solver, Shock Waves, 4 (1994) 25-34. [58] E.F. Toro, Riemann solvers and numerical methods for fluid dynamics: a practical introduction, Springer Science & Business Media, 2013. [59] S. Gottlieb, C.-W. Shu, Total variation diminishing Runge-Kutta schemes, Math. Comp., 67 (1998) 73--85. [60] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, Journal of Computational Physics, 77 (1988) 439-471. [61] X. Zhong, C.-W. Shu, A simple weighted essentially nonoscillatory limiter for Runge–Kutta discontinuous Galerkin methods, Journal of Computational Physics, 232 (2013) 397-415. [62] G.-S. Jiang, C.-W. Shu, Efficient Implementation of Weighted ENO Schemes, J. Comput. Phys., 126 (1996) 202-228. [63] G.-S. Jiang, D. Peng, Weighted ENO Schemes for Hamilton--Jacobi Equations, SIAM Journal on Scientific Computing, 21 (2000) 2126-2143. [64] D. Adalsteinsson, J.A. Sethian, The Fast Construction of Extension Velocities in Level Set Methods, Journal of Computational Physics, 148 (1999) 2-22. [65] W.F. Xie, T.G. Liu, B.C. Khoo, The simulation of cavitating flows induced by underwater shock and free surface interaction, Applied Numerical Mathematics, 57 (2007) 734-745. [66] T.D. Aslam, A partial differential equation approach to multidimensional extrapolation, Journal of Computational Physics, 193 (2004) 349-355. [67] P. Roe, Characteristic-based schemes for the Euler equations, Annual review of fluid mechanics, 18 (1986) 337-365.

46

[68] W.F. Xie, Y.L. Young, T.G. Liu, B.C. Khoo, Dynamic response of deformable structures subjected to shock load and cavitation reload, Computational Mechanics, 40 (2006) 667-681. [69] J. Qiu, T. Liu, B.C. Khoo, Runge–Kutta discontinuous Galerkin methods for compressible two-medium flow simulations: One-dimensional case, Journal of Computational Physics, 222 (2007) 353-373. [70] J. Qiu, T.G. Liu, B.C. Khoo, Simulations of compressible two-medium flow by Runge-Kutta discontinuous Galerkin methods with the ghost fluid method, Communications in Computational Physics, 3 (2008) 479-504. [71] J. Zhu, T. Liu, J. Qiu, B.C. Khoo, RKDG methods with WENO limiters for unsteady cavitating flow, Computers & Fluids, 57 (2012) 52-65. [72] H. Terashima, G. Tryggvason, A front-tracking method with projected interface conditions for compressible multi-fluid flows, Computers & Fluids, 39 (2010) 1804-1814. [73] W. Bo, J.W. Grove, A volume of fluid method based ghost fluid method for compressible multi-fluid flows, Computers & Fluids, 90 (2014) 113-122. [74] K.-M. Shyue, An Efficient Shock-Capturing Algorithm for Compressible Multicomponent Problems, Journal of Computational Physics, 142 (1998) 208-242. [75] S. Oxelosund, Armox Protection Plate, AM International Specification, Edition, 2 (2001). [76] W.F. Xie, Y.L. Young, T.G. Liu, Multiphase modeling of dynamic fluid–structure interaction during close-in explosion, International Journal for Numerical Methods in Engineering, 74 (2008) 1019-1043. [77] F. Zhu, Z. Wang, G. Lu, L. Zhao, Analytical investigation and optimal design of sandwich panels subjected to shock loading, Materials & Design, 30 (2009) 91-100. [78] K.P. Dharmasena, H.N.G. Wadley, Z. Xue, J.W. Hutchinson, Mechanical response of metallic honeycomb sandwich panel structures to high-intensity dynamic loading, International Journal of Impact Engineering, 35 (2008) 1063-1074. [79] C. Qi, S. Yang, L.-J. Yang, Z.-Y. Wei, Z.-H. Lu, Blast resistance and multi-objective optimization of aluminum foam-cored sandwich panels, Composite Structures, 105 (2013) 45-57.

47

Declaration of interests ‫ ܈‬The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ‫܆‬The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: