# Application of the Spectral Element Method in a Surface Ship Far-Field UNDEX Problem

Authors: Zhaokuan Lu; Alan Brown
Journal: Shock and Vibration
 (2019)

---

## Abstract

The prediction of surface ship response to a far-field underwater explosion (UNDEX) requires the simulation of shock wave propagation in the fluid, cavitation, fluid-structure interaction, and structural response. Effective approaches to model the fluid include cavitating acoustic finite element (CAFE) and cavitating acoustic spectral element (CASE) methods. Although the spectral element method offers the potential for greater accuracy at lower computational cost, it also generates more spurious oscillations around discontinuities which are difficult to avoid in shock-related problems. Thus, the advantage of CASE remains unproven. In this paper, we present a 3D-partitioned FSI framework and investigate the application of CAFE and CASE to a surface ship early-time far-field UNDEX problem to determine which method has the best computational efficiency for this problem. We also associate the accuracy of the structural response with the modeling of cavitation distribution. A further contribution of this work is the examination of different nonmatching mesh information exchange schemes to demonstrate how they affect the structural response and improve the CAFE/CASE methodologies.


---

## Body

## 1. Introduction

### 1.1. Background

One of the most studied aspects of naval ship survivability is the vulnerability of a ship to an underwater explosion. Shock trials are rarely performed before a naval ship is put into service when valuable information on the ship and onboard equipment response to underwater shock could be obtained. The process is long, extremely expensive, and environmentally unfriendly. Thus, early assessment has been mostly performed with computational prediction.UNDEX events can be categorized as near-field and far-field. In near-field explosions, the explosive is sufficiently close to the ship’s hull and sufficiently large that its effects result in large plastic deformation [1, 2]. Since the energy is mostly absorbed in the ship’s hull by the deformation, the destructive effect is typically local, but severe. In this work, we focus on the early-time far-field UNDEX where the explosive is small or relatively far from the ship’s hull, and the shock wave released by the explosive has decayed to a linear acoustic wave with a discontinuous pressure profile as it reaches the ship [2]. The structural deformation is mostly or entirely elastic. The high-frequency and low-inertia response of the hull plate panels is excited due to the propagation of the energy into the structure. Onboard equipment and people can be exposed to dangerous accelerations. At the point where the hull panels move inward faster than the surrounding water, the water is subject to tension. A cavitation region called local cavitation is developed near the panel due to the inability of water to sustain tension. Eventually, the local cavitation begins to shrink and collapses as the pressure near the panel is restored and the panel rebounds. The collapse creates an impact pressure pulse and panel reloading [3]. At about the same time, the incident shock wave reaches the free surface and the reflection of the shock wave generates a tensile wave. Due to the large impedance difference between water and air, the reflected wave pressure is almost the negative of the incident shock wave, which significantly reduces the total pressure. A large cavitation region is formed when the total pressure drops below vapor pressure and remains cavitated until the pressure is restored. This cavitation region is called bulk cavitation. The closure of bulk cavitation also emits a shock wave, the timing and impact of which depend on the ship’s proximity to the closure point [2]. The combination of incident shock wave loading and subsequent response to cavitation constitutes the early-time response of a surface ship to the far-field UNDEX.


### 1.2. Modeling Techniques for the Fluid Domain in UNDEX

The fluid-structure interaction of the UNDEX problem was first addressed by Taylor [4] who proposed the Taylor plate theory to calculate structural motion under shock impact. The fluid-structure interaction between an exponentially decaying plane shock wave and an infinite air-backed rigid plate was analyzed. The relationship between the wetted surface pressure and plate motion is obtained using Newton’s second law and a linear acoustic assumption. The Taylor plate theory provides a quick and simple way of estimating peak kickoff velocity (the highest upward velocity induced by shock) where hull structural panels are generalized as rigid plates [5]. However, no interaction between cavitated fluid and structure is captured since the fluid domain is not modeled.The boundary element method (BEM) has also been used in far-field UNDEX problems. The advantage of the BEM is that little or no actual fluid modeling is required, thus allowing faster computation. One of the most used BEM methods is the doubly asymptotic approximation (DAA) by Geers [6–8] which applies partial differential equations to relate the pressure and displacement of a structural wet surface to the fluid boundary using the acoustic fluid assumption. The formulation of the first-order DAA includes an early-time approximation (ETA) and late-time approximation (LTA). The name doubly asymptotic suggests that the approximation approaches the exact solution at both early-time/high-frequency and late-time/low-frequency limits. DAA is used in the Underwater Shock Analysis (USA) program [9] which simulates with the UNDEX problems without cavitation. The motion of a submerged cylindrical shell exposed to a far-field underwater explosion, where the high ambient pressure prevents cavitation, was solved using the USA code in [10], and results compare well with experimental data. However, if the cavitation effect cannot be neglected as in the surface ship problem, the fluid domain or a portion of it needs to be explicitly modeled. The response of a floating barge exposed to an UNDEX simulated by the USA code with and without cavitation shows a significant difference in structural response [11]. To improve the pure BEM based on DAA when cavitation presents, several wet-surface approximation (WSA) models have been introduced (e.g., [12]). The investigation of WSA coupled with elastic mass-spring system conducted in [13] revealed the degraded accuracy due to the inability of capturing the fluid accretion (uncavitated fluid between the cavitation upper boundary and the structure).A more rigorous way of simulating the events involving fluid cavitation is to explicitly discretize the fluid domain and couple it with the structural solver. Since the shock wave decays to a constitutively linear acoustic wave in the far-field UNDEX case, finite element-based approaches are typically favored due to their advantages in solving acoustic wave propagation. The foundation of the FEM-type treatment was proposed by Newton [14] who used a scalar displacement potential-based formulation to model the acoustic fluid and a bilinear equation of state for cavitation. Newton also observed the effect of frothing (spurious cavitation region) caused by the discontinuity near the cavitation-fluid boundary [14]. Newton’s model was later used in [15] to study seismic waves and cavitation effects on a dam. Felippa and Deruntz [16] proposed the cavitating acoustic finite element (CAFE) method based on the work of Newton and implemented it in their cavitating fluid analysis (CFA) code. A damping factor was introduced to reduce the frothing. The CAFE approach was later used with the commercial FEM code LS-DYNA and applied in several far-field UNDEX problems (e.g., [17, 18]). The Abaqus explicit acoustic solver has also been used by many researchers [19–23]. Abaqus takes a similar approach to CAFE using the linear acoustic element (AC3D8R), but with a pressure-based formulation [24]. Despite the benefits of using the FEM-based fluid models (e.g., CAFE), the dispersive nature of the FEM elements [25] makes them computationally expensive since high-level mesh refinement is needed to achieve a satisfactory solution.To resolve the dispersive nature of traditional FEM, the spectral element method (SEM) was introduced in the UNDEX problem by Sprague and Geers [26–29] in place of linear FEM. Incorporating the framework proposed by Felippa and Deruntz [16], their new solver was called cavitating acoustic spectral element (CASE) method. Introduced by Patera [30], SEM is a combination of the spectral method and traditional FEM. High-order basis functions with interpolation nodes at the root of orthogonal polynomials are used to reduce the interpolation error which makes SEM less dispersive than linear FEM. Although the high-order SEM suffers from increased spurious oscillation for the discontinuous shock wave (Gibbs phenomenon), the CASE approach exhibited the same accuracy with less computational resources compared to CAFE in a 1D mass-spring plane shock wave problem [26] and a 3D submerged spherical shell problem [27]. It has also been shown that CASE is manageable for a modern PC for the surface ship far-field UNDEX problem [28]. Klenow [31] assessed the superior behavior of CASE in 1D and associated the higher accuracy of the CASE with the better cavitation region modeling. The same submerged spherical shell in [27] with cavitation effect neglected was investigated in [32] using SEM. The results again showed the higher accuracy of the CASE compared to CAFE with the same degrees of freedom. However, the CASE solution has yet to be compared to CAFE in the far-field surface ship UNDEX problem which features complex structure, hull shape, and both local and bulk cavitation regions.Due to the different phenomenon and geometric simplicity of the submerged spherical shell (deep-water UNDEX) compared to the surface ship, the advantage of CASE found in the Sprague and Geers study [26] does not necessarily exist in a complex surface ship application. Thus, it remains to be determined whether high-order SEM possesses a computational advantage over other types of FEM for the surface ship UNDEX problem.


### 1.3. Information Exchange Schemes at the Fluid-Structure Interface

In early CAFE applications, node-to-node matching of the linear FEM-based fluid and the structural mesh was used. Since the fluid element is usually smaller, the structural element size needs to be decreased to meet the requirement on the boundary, which lowers the stable time increment and increases the computational cost. The one-to-one matching scheme also makes it difficult to connect the fluid mesh to a complex structural mesh [33]. This approach is even less applicable to a high-order SEM-based fluid solver coupled to structural linear FEM solver since the nodes are arranged differently in the spectral element. Therefore, pressure and displacement usually have to be interpolated between the nonmatching meshes using various data exchange schemes.The desirable features of effective data exchange schemes are global load (force) conservation, energy conservation, and mapping consistency. Load conservation means the total fluid nodal force does not change after it is received on the structural side. By the same principle, energy conservation requires the balance of total work input/received by the fluid and received/input by the structure. Mapping consistency guarantees the constant pressure (interface force patch test) and a linearly varying displacement field (interface rigid-motion path test) to be exactly mapped [34]. Following the nomenclature and classification of [34], interface treatments can be categorized as prime dual, dual, and direct force-motion transfer (DFMT) methods. Prime dual (also called the mortar method [35]) and dual approaches employ global and localized Lagrange multipliers to enforce the consistency of displacement or pressure during the fluid and structure interaction. The mapping error is minimized by solving a weighted residual-based equation system. The aforementioned features can all be satisfied (e.g., [36]). However, the prime dual and dual methods involve a large amount of computation since all the DOFs on the boundary are coupled and the equation system needs to be solved.In our research, we focus on direct force-motion transfer (DFMT) methods which are local in nature (involving only neighboring nodes for interpolation), thus less computational effort. The implementation of DFMT is also greatly simpler than the prime dual and dual methods. Among the different direct force-motion transfer (DFMT) schemes for nonmatching meshes, the consistent interpolation approach [37–39] and energy-conservative mapping approach [38] are often used. The consistent approach guarantees the mapping consistency of the pressure and the displacement fields. Although the consistent algorithm does not enforce force and energy conservation, it does not necessarily give an unacceptable result as demonstrated in [38]. The conservative approach enforces conservation of energy but loses consistency when the structural mesh is finer than the fluid mesh [40]. It generally guarantees consistency when the fluid mesh is finer than the structural mesh which is the case for UNDEX where the cavitation effect needs to be captured. It has been suggested that the energy conservation is of little importance for a short-duration shock-wave transient simulation [34, 41]. However, no quantification of this is found in the literature.The CASE has previously been coupled to the structure by an average mapping technique which maps the average pressure and displacement fields in [27]. Although this method features straightforward implementation, high-pressure gradients induced by the oblique shock wave are smoothed. Thus, the accuracy of the peak velocity prediction may suffer. It has not been demonstrated whether CASE would preserve superior behavior if a more rigorous information mapping scheme (e.g., consistent scheme) is used.In the research described in this paper, we employ the basic ideas of CAFE and CASE for the fluid modeling and couple them to an Abaqus explicit structural solver. As a preliminary step to applying our approach to a full surface ship problem, we use a floating shock platform (FSP) problem to show the potential relative computational advantage in far-field surface ship UNDEX problems. In addition, consistent, average, and energy-conservative mapping schemes are explored and compared to determine which is the most cost-effective approach.


## 1.1. Background

One of the most studied aspects of naval ship survivability is the vulnerability of a ship to an underwater explosion. Shock trials are rarely performed before a naval ship is put into service when valuable information on the ship and onboard equipment response to underwater shock could be obtained. The process is long, extremely expensive, and environmentally unfriendly. Thus, early assessment has been mostly performed with computational prediction.UNDEX events can be categorized as near-field and far-field. In near-field explosions, the explosive is sufficiently close to the ship’s hull and sufficiently large that its effects result in large plastic deformation [1, 2]. Since the energy is mostly absorbed in the ship’s hull by the deformation, the destructive effect is typically local, but severe. In this work, we focus on the early-time far-field UNDEX where the explosive is small or relatively far from the ship’s hull, and the shock wave released by the explosive has decayed to a linear acoustic wave with a discontinuous pressure profile as it reaches the ship [2]. The structural deformation is mostly or entirely elastic. The high-frequency and low-inertia response of the hull plate panels is excited due to the propagation of the energy into the structure. Onboard equipment and people can be exposed to dangerous accelerations. At the point where the hull panels move inward faster than the surrounding water, the water is subject to tension. A cavitation region called local cavitation is developed near the panel due to the inability of water to sustain tension. Eventually, the local cavitation begins to shrink and collapses as the pressure near the panel is restored and the panel rebounds. The collapse creates an impact pressure pulse and panel reloading [3]. At about the same time, the incident shock wave reaches the free surface and the reflection of the shock wave generates a tensile wave. Due to the large impedance difference between water and air, the reflected wave pressure is almost the negative of the incident shock wave, which significantly reduces the total pressure. A large cavitation region is formed when the total pressure drops below vapor pressure and remains cavitated until the pressure is restored. This cavitation region is called bulk cavitation. The closure of bulk cavitation also emits a shock wave, the timing and impact of which depend on the ship’s proximity to the closure point [2]. The combination of incident shock wave loading and subsequent response to cavitation constitutes the early-time response of a surface ship to the far-field UNDEX.


## 1.2. Modeling Techniques for the Fluid Domain in UNDEX

The fluid-structure interaction of the UNDEX problem was first addressed by Taylor [4] who proposed the Taylor plate theory to calculate structural motion under shock impact. The fluid-structure interaction between an exponentially decaying plane shock wave and an infinite air-backed rigid plate was analyzed. The relationship between the wetted surface pressure and plate motion is obtained using Newton’s second law and a linear acoustic assumption. The Taylor plate theory provides a quick and simple way of estimating peak kickoff velocity (the highest upward velocity induced by shock) where hull structural panels are generalized as rigid plates [5]. However, no interaction between cavitated fluid and structure is captured since the fluid domain is not modeled.The boundary element method (BEM) has also been used in far-field UNDEX problems. The advantage of the BEM is that little or no actual fluid modeling is required, thus allowing faster computation. One of the most used BEM methods is the doubly asymptotic approximation (DAA) by Geers [6–8] which applies partial differential equations to relate the pressure and displacement of a structural wet surface to the fluid boundary using the acoustic fluid assumption. The formulation of the first-order DAA includes an early-time approximation (ETA) and late-time approximation (LTA). The name doubly asymptotic suggests that the approximation approaches the exact solution at both early-time/high-frequency and late-time/low-frequency limits. DAA is used in the Underwater Shock Analysis (USA) program [9] which simulates with the UNDEX problems without cavitation. The motion of a submerged cylindrical shell exposed to a far-field underwater explosion, where the high ambient pressure prevents cavitation, was solved using the USA code in [10], and results compare well with experimental data. However, if the cavitation effect cannot be neglected as in the surface ship problem, the fluid domain or a portion of it needs to be explicitly modeled. The response of a floating barge exposed to an UNDEX simulated by the USA code with and without cavitation shows a significant difference in structural response [11]. To improve the pure BEM based on DAA when cavitation presents, several wet-surface approximation (WSA) models have been introduced (e.g., [12]). The investigation of WSA coupled with elastic mass-spring system conducted in [13] revealed the degraded accuracy due to the inability of capturing the fluid accretion (uncavitated fluid between the cavitation upper boundary and the structure).A more rigorous way of simulating the events involving fluid cavitation is to explicitly discretize the fluid domain and couple it with the structural solver. Since the shock wave decays to a constitutively linear acoustic wave in the far-field UNDEX case, finite element-based approaches are typically favored due to their advantages in solving acoustic wave propagation. The foundation of the FEM-type treatment was proposed by Newton [14] who used a scalar displacement potential-based formulation to model the acoustic fluid and a bilinear equation of state for cavitation. Newton also observed the effect of frothing (spurious cavitation region) caused by the discontinuity near the cavitation-fluid boundary [14]. Newton’s model was later used in [15] to study seismic waves and cavitation effects on a dam. Felippa and Deruntz [16] proposed the cavitating acoustic finite element (CAFE) method based on the work of Newton and implemented it in their cavitating fluid analysis (CFA) code. A damping factor was introduced to reduce the frothing. The CAFE approach was later used with the commercial FEM code LS-DYNA and applied in several far-field UNDEX problems (e.g., [17, 18]). The Abaqus explicit acoustic solver has also been used by many researchers [19–23]. Abaqus takes a similar approach to CAFE using the linear acoustic element (AC3D8R), but with a pressure-based formulation [24]. Despite the benefits of using the FEM-based fluid models (e.g., CAFE), the dispersive nature of the FEM elements [25] makes them computationally expensive since high-level mesh refinement is needed to achieve a satisfactory solution.To resolve the dispersive nature of traditional FEM, the spectral element method (SEM) was introduced in the UNDEX problem by Sprague and Geers [26–29] in place of linear FEM. Incorporating the framework proposed by Felippa and Deruntz [16], their new solver was called cavitating acoustic spectral element (CASE) method. Introduced by Patera [30], SEM is a combination of the spectral method and traditional FEM. High-order basis functions with interpolation nodes at the root of orthogonal polynomials are used to reduce the interpolation error which makes SEM less dispersive than linear FEM. Although the high-order SEM suffers from increased spurious oscillation for the discontinuous shock wave (Gibbs phenomenon), the CASE approach exhibited the same accuracy with less computational resources compared to CAFE in a 1D mass-spring plane shock wave problem [26] and a 3D submerged spherical shell problem [27]. It has also been shown that CASE is manageable for a modern PC for the surface ship far-field UNDEX problem [28]. Klenow [31] assessed the superior behavior of CASE in 1D and associated the higher accuracy of the CASE with the better cavitation region modeling. The same submerged spherical shell in [27] with cavitation effect neglected was investigated in [32] using SEM. The results again showed the higher accuracy of the CASE compared to CAFE with the same degrees of freedom. However, the CASE solution has yet to be compared to CAFE in the far-field surface ship UNDEX problem which features complex structure, hull shape, and both local and bulk cavitation regions.Due to the different phenomenon and geometric simplicity of the submerged spherical shell (deep-water UNDEX) compared to the surface ship, the advantage of CASE found in the Sprague and Geers study [26] does not necessarily exist in a complex surface ship application. Thus, it remains to be determined whether high-order SEM possesses a computational advantage over other types of FEM for the surface ship UNDEX problem.


## 1.3. Information Exchange Schemes at the Fluid-Structure Interface

In early CAFE applications, node-to-node matching of the linear FEM-based fluid and the structural mesh was used. Since the fluid element is usually smaller, the structural element size needs to be decreased to meet the requirement on the boundary, which lowers the stable time increment and increases the computational cost. The one-to-one matching scheme also makes it difficult to connect the fluid mesh to a complex structural mesh [33]. This approach is even less applicable to a high-order SEM-based fluid solver coupled to structural linear FEM solver since the nodes are arranged differently in the spectral element. Therefore, pressure and displacement usually have to be interpolated between the nonmatching meshes using various data exchange schemes.The desirable features of effective data exchange schemes are global load (force) conservation, energy conservation, and mapping consistency. Load conservation means the total fluid nodal force does not change after it is received on the structural side. By the same principle, energy conservation requires the balance of total work input/received by the fluid and received/input by the structure. Mapping consistency guarantees the constant pressure (interface force patch test) and a linearly varying displacement field (interface rigid-motion path test) to be exactly mapped [34]. Following the nomenclature and classification of [34], interface treatments can be categorized as prime dual, dual, and direct force-motion transfer (DFMT) methods. Prime dual (also called the mortar method [35]) and dual approaches employ global and localized Lagrange multipliers to enforce the consistency of displacement or pressure during the fluid and structure interaction. The mapping error is minimized by solving a weighted residual-based equation system. The aforementioned features can all be satisfied (e.g., [36]). However, the prime dual and dual methods involve a large amount of computation since all the DOFs on the boundary are coupled and the equation system needs to be solved.In our research, we focus on direct force-motion transfer (DFMT) methods which are local in nature (involving only neighboring nodes for interpolation), thus less computational effort. The implementation of DFMT is also greatly simpler than the prime dual and dual methods. Among the different direct force-motion transfer (DFMT) schemes for nonmatching meshes, the consistent interpolation approach [37–39] and energy-conservative mapping approach [38] are often used. The consistent approach guarantees the mapping consistency of the pressure and the displacement fields. Although the consistent algorithm does not enforce force and energy conservation, it does not necessarily give an unacceptable result as demonstrated in [38]. The conservative approach enforces conservation of energy but loses consistency when the structural mesh is finer than the fluid mesh [40]. It generally guarantees consistency when the fluid mesh is finer than the structural mesh which is the case for UNDEX where the cavitation effect needs to be captured. It has been suggested that the energy conservation is of little importance for a short-duration shock-wave transient simulation [34, 41]. However, no quantification of this is found in the literature.The CASE has previously been coupled to the structure by an average mapping technique which maps the average pressure and displacement fields in [27]. Although this method features straightforward implementation, high-pressure gradients induced by the oblique shock wave are smoothed. Thus, the accuracy of the peak velocity prediction may suffer. It has not been demonstrated whether CASE would preserve superior behavior if a more rigorous information mapping scheme (e.g., consistent scheme) is used.In the research described in this paper, we employ the basic ideas of CAFE and CASE for the fluid modeling and couple them to an Abaqus explicit structural solver. As a preliminary step to applying our approach to a full surface ship problem, we use a floating shock platform (FSP) problem to show the potential relative computational advantage in far-field surface ship UNDEX problems. In addition, consistent, average, and energy-conservative mapping schemes are explored and compared to determine which is the most cost-effective approach.


## 2. Description of the FSI Framework, Structure, and Fluid Model

### 2.1. Coupled FSI Framework through MpCCI

In this section, we present a partitioned FSI framework for the far-field UNDEX problem. On the structural side, we use the Abaqus Explicit solver which features explicit-time integration. For the fluid solver, the CAFE and CASE methodologies are implemented in our own fluid solver in C++. For the coupling between fluid and structural solvers, MpCCI [42], a multiphysics coupling software is utilized. MpCCI supports most state-of-the-art commercial solvers and our in-house code. Abaqus is coupled through its cosimulation engine, and the in-house code is coupled through MpCCI API functions. The details regarding the resolution of nonmatching meshes and time-step size are introduced later in this section. The skeleton of the coupled framework is shown in Figure 1.Figure 1
Partitioned FSI framework.


### 2.2. Abaqus (Explicit) Structural Solver

Since an important effect of the far-field UNDEX problem is the structural oscillation caused by the stress wave propagation after the shock impact, Abaqus Explicit is chosen due to its advantage in solving problems with transient high-speed events. In comparison, the implicit solver (Abaqus Standard) is more suitable for the low-speed phenomena like structure dynamic motion after the shock wave effect is gone. The semidiscretized structural motion is expressed in(1)Mx¨+Kx=f.The mass matrix is lumped to save memory and facilitate an explicit central difference Newmark scheme detailed in [24].


### 2.3. Cavitating Acoustic Spectral/Finite Element (CASE/CAFE) Method

In this section, we describe the governing equations and spatial discretization of the CASE/CAFE framework which are adapted from [16, 27, 31].The far-field assumption allows the acoustic treatment of the shock wave propagation. Although the acoustic fluid is compressible, we assume the fluctuations between the densityρ and reference density ρ0 are small, thus allowing the continuity equation based on dynamic displacement x→:(2)∂ρ∂t+∇⋅ρ∂x→∂t=0,to be linearized using total condensation S=ρ/ρ0−1 for S≪1:(3)∂S∂t+∇⋅∂x→∂t=0.Neglecting the viscosity and convection term, the Euler momentum equation for the acoustic fluid is written as(4)ρ0∂2x⇀∂t2=ρ0f→−∇pt,where pt is the total fluid pressure and fX→ is the body force field.In this work, we apply the total field formulation (TFM) which comprises the incident, scattered, and static field following the approach from [16]. The incoming shock wave causes the incident field, and its reflection from the boundaries (i.e., structure and free surface) generates the scattered field. The static field is the hydrostatic field of the fluid. The primary variables of TFM only contain the incident and the scattered field, which in combination is called the dynamic field. The static field is not modeled since it is known for all time periods.Using the relationship−∇ψ=ρ0x→ based on irrotational fluid assumption, the governing equations, equations (3) and (4), are converted to the form as follows:(5)s=∇2ψ,(6)∂2ψ∂t2=c2s,s>pcav−peqc2,pcav−peq,otherwise,where s is the dynamic condensation and ψ is the dynamic displacement potential related to the dynamic pressure by p=ψ¨. Equation (6) enforces a cutoff cavitation condition which stops the total pressure going below vapor pressure pcav once cavitation begins. The vapor pressure is small enough to be negligible. peq is the static field pressure.Following the standard Galerkin approach, we multiply equation (5) by test function ϕ and integrate it over the fluid domain Ωf:(7)∫Ωfsϕ−∇2ψϕdΩ=0.Applying Green’s first identity to equation (7) yields the weak form:(8)∫ΩfsϕdΩ+∫Ωf∇ψ⋅∇ϕdΩ=∫Γfϕ∇ψ⋅n→dΓ,where Γf denotes the fluid boundary and n→ is the outward pointing normal vector on the boundary.Due to the fact that the fluid mesh comprises all hexahedral elements, the densified condensation and displacement potential can be written over an element as(9)sξ,η,ζ,t=ϕTst=∑i,j,k=1N+1ϕiξϕjηϕkζsijkt,(10)ψξ,η,ζ,t=ϕTψt=∑i,j,k=1N+1ϕiξϕjηϕkζψijkt,where ϕ is the column vector of the 1D spectral basis function ϕiξ,ϕjη,andϕkζ. For the CASE, the interpolation point and quadrature point are both chosen as Gauss–Lobatto–Legendre (GLL), thus allowing the mass matrix to be fully diagonal. This is the essential difference between the CASE and CAFE which use equally spaced interpolation point and Gauss–Legendre quadrature nodes. As a result, the mass matrix is not diagonal. For the sake of fast computation in the CAFE explicit temporal integration, the mass matrix is diagonalized by the row-sum technique.Substituting equations (9) and (10) into the weak form (equation (8)) gives(11)Qs+Hψ=b,where Q=∫ΩfϕϕTdΩ, H=∫Ωf∇ϕ⋅∇ϕTdΩ, and b=∫Γfϕ∇ψ⋅n→dΓ. Following the nomenclature in [16], we name Q the global capacitance matrix, H the global reactance matrix, and b the global fluid force vector. The governing equation (equation (6)) is evaluated in a node-by-node fashion:(12)ψ¨=c2s∣i,s∣i>−peq∣i/c2,i=1,2,…,nnode,pcav−peq,otherwise.The boundary conditions to be prescribed are at the free surface, fluid-structure interface, and nonreflecting boundary (NRB). At the free surface, the essential boundary condition is imposed as(13)p=0.The fluid-structure interaction and NRB terms are natural boundary conditions governed by the fluid force vectorb which can be expressed as(14)b=bs+bnrb,where(15)bs=∫Γfsϕ∇ψ⋅n→dΓfs=∫Γfsϕ∂ψ∂ndΓfs=−ρ0∫ΓfsϕϕTxsdΓfs,where xs is the displacement vector of the structure along the normal direction of fluid-structure interface taken as positive into the fluid. By the same principle,(16)bnrb=−ρ0∫ΓnrbϕϕTxnrbdΓnrb,where xnrb is the displacement of NRB approximated by curved wave approximation [29]:(17)xnrb−xeq=1ρ0cΘpinc∗+1ρ0ΘR−1pinc∗∗+1ρ0cp∗−peq∗−pinc∗+κρ0p∗∗−peq∗∗−pinc∗∗,where p=pinc+psc, xeq is the displacement at equilibrium, Θ is the cos value of the angle between the incident shock direction and the outward norm n→. κ represents the mean curvature of the NRB. The terms in the first bracket represent the displacement caused by the incident wave, and the terms in the second bracket are the approximation of the displacement caused by the scattered wave. Since the NRB used in later investigation is planar (κ=0), the last term is dropped.


### 2.4. Nonconformal Fluid-Structure Mesh Coupling

Due to the discontinuity of the shock wave and nonlinear cavitation effect, the fluid mesh must typically be finer than the structural mesh. To facilitate the mapping between nonmatching meshes by the methods introduced in this chapter, a general purpose searching and projection algorithm is implemented. The node of interest from the fluid (structure) surface mesh is first paired with the nearest element on the structure (fluid) surface mesh. Then, the node is orthogonally projected onto the paired element. The local coordinate of the projected node is determined to facilitate the later interpolation. If the projection fails, the value of the closest node on the target mesh is directly used. To handle the interpolation of pressure and displacement field, three aforementioned different DFMT mapping algorithms are introduced below.


#### 2.4.1. Consistent and Average Mapping Algorithm

The consistent mapping scheme follows a straightforward implementation by projecting the target nodes to the source mesh and then interpolating the desired properties using the local shape function. Specifically, the structural quadrature nodesχi (i.e., Gauss–Legendre) are mapped to the fluid mesh to obtain the interpolated pressure (equation (18)) and the displacement is obtained by projecting the fluid nodes on the structure mesh and interpolating using the local structural shape function:(18)pχi=∑j=1j=nfNjχipj,(19)uj=∑i=1i=nsNiχjui,where i denotes the nodes on the structure element and j on the fluid element, p is the pressure and u is the nodal displacement, ns and nf are the node number on structure and fluid element, and N is the shape function. In Figure 2, we use a 3rd order CASE coupled to the linear structural element as an illustration of this method.Figure 2
Consistent mapping algorithm.For the average mapping algorithm in [27], the pressure and displacement interpolations are further loosened by taking their averaged values (equations (20) and (21)). Thus, the searching and projection routine is simplified since the local coordinate of the projected node is no longer required. Only the pairing relationship between the node and target element is needed:(20)pχi=p¯j,(21)uj=u¯i.The consistent and average mapping approaches do not attempt to guarantee the force and further energy conservation. The level of conservation depends on the number of Gauss points used.


#### 2.4.2. Energy-Conservative Mapping Algorithm

From the virtual work balance derivation in [38], the nodal force (i.e., pressure flux) mapping matrix must be transposed to the displacement mapping matrix in order to guarantee energy conservation. For instance, if we use the consistent mapping for displacement, then we have to use the same shape function (i.e., the shape function of the structural element) to interpolate the nodal force. Thus, for the fully consistent mapping approach in Section 2.4.1, energy conservation is not guaranteed since the mapping matrices are not related.In the later context, we interpolate the displacement using the consistent approach:(22)uj=∑i=1i=nsNiχjui.Then, the received nodal forcefi on the structure element becomes(23)fi=∑j=1j=nfϕjNiχj,in which χj is the local coordinate of the mapped nf fluid nodes on the structure element and ϕj is the nodal force on the fluid element which can be represented as follows:(24)ϕj=∫ΓfeNjpn→ds,j∈Γfe,where Γfe is the fluid surface.Based on the unit summation property of the shape function (i.e.,∑i=1i=ieNiχj=1), the total received force fts by the structure element becomes(25)fts=∑i=1i=nsfi=∑i=1i=ns∑j=1j=nfϕjNiχj=∑j=1j=nfϕj.Thus, the total force is conservative as well. In Figure3, we use a 3rd-order CASE coupled to the linear structural element as an example for the energy-conservative mapping algorithm derivation.Figure 3
Energy-conservative mapping algorithm.


### 2.5. Temporal Discretization

The governing equations of the CAFE/CASE framework are integrated in time by using a staggered second-order explicit central difference scheme (a.k.a. leapfrog scheme) introduced by Felippa and Deruntz [16]. The approximation of critical time increment is given as(26)Δtcrt=2cλmax1+2β,where λmax is the maximum fluid mesh eigenvalue estimated by Gershgorin’s theorem [29] for CAFE/CASE and β is the artificial damping factor used to mitigate the spurious oscillation. Half of the critical time increment (CFL=0.5) is used in the cases hereafter.The nonmatching time-step interpolation is performed by a monotone cubic Hermite spline implemented in MpCCI to achieve a 3rd-order interpolation accuracy [42]. Different from the subcycling approach, the current approach does not require either solver to cut the time increment to make the ratio integer, thus allowing them to advance at their own optimal step size. Also, this would make the later comparison of CASE and CAFE more consistent by keeping the CFL number constant. This idea is illustrated in Figure 4.Figure 4
Time-step interpolation.


## 2.1. Coupled FSI Framework through MpCCI

In this section, we present a partitioned FSI framework for the far-field UNDEX problem. On the structural side, we use the Abaqus Explicit solver which features explicit-time integration. For the fluid solver, the CAFE and CASE methodologies are implemented in our own fluid solver in C++. For the coupling between fluid and structural solvers, MpCCI [42], a multiphysics coupling software is utilized. MpCCI supports most state-of-the-art commercial solvers and our in-house code. Abaqus is coupled through its cosimulation engine, and the in-house code is coupled through MpCCI API functions. The details regarding the resolution of nonmatching meshes and time-step size are introduced later in this section. The skeleton of the coupled framework is shown in Figure 1.Figure 1
Partitioned FSI framework.


## 2.2. Abaqus (Explicit) Structural Solver

Since an important effect of the far-field UNDEX problem is the structural oscillation caused by the stress wave propagation after the shock impact, Abaqus Explicit is chosen due to its advantage in solving problems with transient high-speed events. In comparison, the implicit solver (Abaqus Standard) is more suitable for the low-speed phenomena like structure dynamic motion after the shock wave effect is gone. The semidiscretized structural motion is expressed in(1)Mx¨+Kx=f.The mass matrix is lumped to save memory and facilitate an explicit central difference Newmark scheme detailed in [24].


## 2.3. Cavitating Acoustic Spectral/Finite Element (CASE/CAFE) Method

In this section, we describe the governing equations and spatial discretization of the CASE/CAFE framework which are adapted from [16, 27, 31].The far-field assumption allows the acoustic treatment of the shock wave propagation. Although the acoustic fluid is compressible, we assume the fluctuations between the densityρ and reference density ρ0 are small, thus allowing the continuity equation based on dynamic displacement x→:(2)∂ρ∂t+∇⋅ρ∂x→∂t=0,to be linearized using total condensation S=ρ/ρ0−1 for S≪1:(3)∂S∂t+∇⋅∂x→∂t=0.Neglecting the viscosity and convection term, the Euler momentum equation for the acoustic fluid is written as(4)ρ0∂2x⇀∂t2=ρ0f→−∇pt,where pt is the total fluid pressure and fX→ is the body force field.In this work, we apply the total field formulation (TFM) which comprises the incident, scattered, and static field following the approach from [16]. The incoming shock wave causes the incident field, and its reflection from the boundaries (i.e., structure and free surface) generates the scattered field. The static field is the hydrostatic field of the fluid. The primary variables of TFM only contain the incident and the scattered field, which in combination is called the dynamic field. The static field is not modeled since it is known for all time periods.Using the relationship−∇ψ=ρ0x→ based on irrotational fluid assumption, the governing equations, equations (3) and (4), are converted to the form as follows:(5)s=∇2ψ,(6)∂2ψ∂t2=c2s,s>pcav−peqc2,pcav−peq,otherwise,where s is the dynamic condensation and ψ is the dynamic displacement potential related to the dynamic pressure by p=ψ¨. Equation (6) enforces a cutoff cavitation condition which stops the total pressure going below vapor pressure pcav once cavitation begins. The vapor pressure is small enough to be negligible. peq is the static field pressure.Following the standard Galerkin approach, we multiply equation (5) by test function ϕ and integrate it over the fluid domain Ωf:(7)∫Ωfsϕ−∇2ψϕdΩ=0.Applying Green’s first identity to equation (7) yields the weak form:(8)∫ΩfsϕdΩ+∫Ωf∇ψ⋅∇ϕdΩ=∫Γfϕ∇ψ⋅n→dΓ,where Γf denotes the fluid boundary and n→ is the outward pointing normal vector on the boundary.Due to the fact that the fluid mesh comprises all hexahedral elements, the densified condensation and displacement potential can be written over an element as(9)sξ,η,ζ,t=ϕTst=∑i,j,k=1N+1ϕiξϕjηϕkζsijkt,(10)ψξ,η,ζ,t=ϕTψt=∑i,j,k=1N+1ϕiξϕjηϕkζψijkt,where ϕ is the column vector of the 1D spectral basis function ϕiξ,ϕjη,andϕkζ. For the CASE, the interpolation point and quadrature point are both chosen as Gauss–Lobatto–Legendre (GLL), thus allowing the mass matrix to be fully diagonal. This is the essential difference between the CASE and CAFE which use equally spaced interpolation point and Gauss–Legendre quadrature nodes. As a result, the mass matrix is not diagonal. For the sake of fast computation in the CAFE explicit temporal integration, the mass matrix is diagonalized by the row-sum technique.Substituting equations (9) and (10) into the weak form (equation (8)) gives(11)Qs+Hψ=b,where Q=∫ΩfϕϕTdΩ, H=∫Ωf∇ϕ⋅∇ϕTdΩ, and b=∫Γfϕ∇ψ⋅n→dΓ. Following the nomenclature in [16], we name Q the global capacitance matrix, H the global reactance matrix, and b the global fluid force vector. The governing equation (equation (6)) is evaluated in a node-by-node fashion:(12)ψ¨=c2s∣i,s∣i>−peq∣i/c2,i=1,2,…,nnode,pcav−peq,otherwise.The boundary conditions to be prescribed are at the free surface, fluid-structure interface, and nonreflecting boundary (NRB). At the free surface, the essential boundary condition is imposed as(13)p=0.The fluid-structure interaction and NRB terms are natural boundary conditions governed by the fluid force vectorb which can be expressed as(14)b=bs+bnrb,where(15)bs=∫Γfsϕ∇ψ⋅n→dΓfs=∫Γfsϕ∂ψ∂ndΓfs=−ρ0∫ΓfsϕϕTxsdΓfs,where xs is the displacement vector of the structure along the normal direction of fluid-structure interface taken as positive into the fluid. By the same principle,(16)bnrb=−ρ0∫ΓnrbϕϕTxnrbdΓnrb,where xnrb is the displacement of NRB approximated by curved wave approximation [29]:(17)xnrb−xeq=1ρ0cΘpinc∗+1ρ0ΘR−1pinc∗∗+1ρ0cp∗−peq∗−pinc∗+κρ0p∗∗−peq∗∗−pinc∗∗,where p=pinc+psc, xeq is the displacement at equilibrium, Θ is the cos value of the angle between the incident shock direction and the outward norm n→. κ represents the mean curvature of the NRB. The terms in the first bracket represent the displacement caused by the incident wave, and the terms in the second bracket are the approximation of the displacement caused by the scattered wave. Since the NRB used in later investigation is planar (κ=0), the last term is dropped.


## 2.4. Nonconformal Fluid-Structure Mesh Coupling

Due to the discontinuity of the shock wave and nonlinear cavitation effect, the fluid mesh must typically be finer than the structural mesh. To facilitate the mapping between nonmatching meshes by the methods introduced in this chapter, a general purpose searching and projection algorithm is implemented. The node of interest from the fluid (structure) surface mesh is first paired with the nearest element on the structure (fluid) surface mesh. Then, the node is orthogonally projected onto the paired element. The local coordinate of the projected node is determined to facilitate the later interpolation. If the projection fails, the value of the closest node on the target mesh is directly used. To handle the interpolation of pressure and displacement field, three aforementioned different DFMT mapping algorithms are introduced below.


### 2.4.1. Consistent and Average Mapping Algorithm

The consistent mapping scheme follows a straightforward implementation by projecting the target nodes to the source mesh and then interpolating the desired properties using the local shape function. Specifically, the structural quadrature nodesχi (i.e., Gauss–Legendre) are mapped to the fluid mesh to obtain the interpolated pressure (equation (18)) and the displacement is obtained by projecting the fluid nodes on the structure mesh and interpolating using the local structural shape function:(18)pχi=∑j=1j=nfNjχipj,(19)uj=∑i=1i=nsNiχjui,where i denotes the nodes on the structure element and j on the fluid element, p is the pressure and u is the nodal displacement, ns and nf are the node number on structure and fluid element, and N is the shape function. In Figure 2, we use a 3rd order CASE coupled to the linear structural element as an illustration of this method.Figure 2
Consistent mapping algorithm.For the average mapping algorithm in [27], the pressure and displacement interpolations are further loosened by taking their averaged values (equations (20) and (21)). Thus, the searching and projection routine is simplified since the local coordinate of the projected node is no longer required. Only the pairing relationship between the node and target element is needed:(20)pχi=p¯j,(21)uj=u¯i.The consistent and average mapping approaches do not attempt to guarantee the force and further energy conservation. The level of conservation depends on the number of Gauss points used.


### 2.4.2. Energy-Conservative Mapping Algorithm

From the virtual work balance derivation in [38], the nodal force (i.e., pressure flux) mapping matrix must be transposed to the displacement mapping matrix in order to guarantee energy conservation. For instance, if we use the consistent mapping for displacement, then we have to use the same shape function (i.e., the shape function of the structural element) to interpolate the nodal force. Thus, for the fully consistent mapping approach in Section 2.4.1, energy conservation is not guaranteed since the mapping matrices are not related.In the later context, we interpolate the displacement using the consistent approach:(22)uj=∑i=1i=nsNiχjui.Then, the received nodal forcefi on the structure element becomes(23)fi=∑j=1j=nfϕjNiχj,in which χj is the local coordinate of the mapped nf fluid nodes on the structure element and ϕj is the nodal force on the fluid element which can be represented as follows:(24)ϕj=∫ΓfeNjpn→ds,j∈Γfe,where Γfe is the fluid surface.Based on the unit summation property of the shape function (i.e.,∑i=1i=ieNiχj=1), the total received force fts by the structure element becomes(25)fts=∑i=1i=nsfi=∑i=1i=ns∑j=1j=nfϕjNiχj=∑j=1j=nfϕj.Thus, the total force is conservative as well. In Figure3, we use a 3rd-order CASE coupled to the linear structural element as an example for the energy-conservative mapping algorithm derivation.Figure 3
Energy-conservative mapping algorithm.


## 2.4.1. Consistent and Average Mapping Algorithm

The consistent mapping scheme follows a straightforward implementation by projecting the target nodes to the source mesh and then interpolating the desired properties using the local shape function. Specifically, the structural quadrature nodesχi (i.e., Gauss–Legendre) are mapped to the fluid mesh to obtain the interpolated pressure (equation (18)) and the displacement is obtained by projecting the fluid nodes on the structure mesh and interpolating using the local structural shape function:(18)pχi=∑j=1j=nfNjχipj,(19)uj=∑i=1i=nsNiχjui,where i denotes the nodes on the structure element and j on the fluid element, p is the pressure and u is the nodal displacement, ns and nf are the node number on structure and fluid element, and N is the shape function. In Figure 2, we use a 3rd order CASE coupled to the linear structural element as an illustration of this method.Figure 2
Consistent mapping algorithm.For the average mapping algorithm in [27], the pressure and displacement interpolations are further loosened by taking their averaged values (equations (20) and (21)). Thus, the searching and projection routine is simplified since the local coordinate of the projected node is no longer required. Only the pairing relationship between the node and target element is needed:(20)pχi=p¯j,(21)uj=u¯i.The consistent and average mapping approaches do not attempt to guarantee the force and further energy conservation. The level of conservation depends on the number of Gauss points used.


## 2.4.2. Energy-Conservative Mapping Algorithm

From the virtual work balance derivation in [38], the nodal force (i.e., pressure flux) mapping matrix must be transposed to the displacement mapping matrix in order to guarantee energy conservation. For instance, if we use the consistent mapping for displacement, then we have to use the same shape function (i.e., the shape function of the structural element) to interpolate the nodal force. Thus, for the fully consistent mapping approach in Section 2.4.1, energy conservation is not guaranteed since the mapping matrices are not related.In the later context, we interpolate the displacement using the consistent approach:(22)uj=∑i=1i=nsNiχjui.Then, the received nodal forcefi on the structure element becomes(23)fi=∑j=1j=nfϕjNiχj,in which χj is the local coordinate of the mapped nf fluid nodes on the structure element and ϕj is the nodal force on the fluid element which can be represented as follows:(24)ϕj=∫ΓfeNjpn→ds,j∈Γfe,where Γfe is the fluid surface.Based on the unit summation property of the shape function (i.e.,∑i=1i=ieNiχj=1), the total received force fts by the structure element becomes(25)fts=∑i=1i=nsfi=∑i=1i=ns∑j=1j=nfϕjNiχj=∑j=1j=nfϕj.Thus, the total force is conservative as well. In Figure3, we use a 3rd-order CASE coupled to the linear structural element as an example for the energy-conservative mapping algorithm derivation.Figure 3
Energy-conservative mapping algorithm.


## 2.5. Temporal Discretization

The governing equations of the CAFE/CASE framework are integrated in time by using a staggered second-order explicit central difference scheme (a.k.a. leapfrog scheme) introduced by Felippa and Deruntz [16]. The approximation of critical time increment is given as(26)Δtcrt=2cλmax1+2β,where λmax is the maximum fluid mesh eigenvalue estimated by Gershgorin’s theorem [29] for CAFE/CASE and β is the artificial damping factor used to mitigate the spurious oscillation. Half of the critical time increment (CFL=0.5) is used in the cases hereafter.The nonmatching time-step interpolation is performed by a monotone cubic Hermite spline implemented in MpCCI to achieve a 3rd-order interpolation accuracy [42]. Different from the subcycling approach, the current approach does not require either solver to cut the time increment to make the ratio integer, thus allowing them to advance at their own optimal step size. Also, this would make the later comparison of CASE and CAFE more consistent by keeping the CFL number constant. This idea is illustrated in Figure 4.Figure 4
Time-step interpolation.


## 3. Numerical Results and Assessment

In this section, we first present a preliminary model assessment and then demonstrate which type of fluid element has the highest computational efficiency for the floating shock platform UNDEX problem.


### 3.1. Model Assessment

In this section, the model and FSI framework introduced in Section2 are assessed using two benchmark problems. The linear CASE is used in the in-house fluid solver, and the consistent mapping algorithm is adopted for the nonmatching mesh interpolation. Excellent agreement with the baseline results provides more confidence in the model setup employed for the numerical experiment described in later sections.


#### 3.1.1. Bleich–Sandler Problem

The Bleich–Sandler plate [43] is a classic UNDEX benchmark problem with cavitation. The problem includes a rigid plate sitting on the free surface of an infinitely deep fluid column (Figure 5). The plate is subject to the load from a step exponential plane shock wave in equation (27) where Ppeak is the peak pressure on the shock front, c is the speed of sound, τ is the decay constant, and xso is the initial location of the shock front at t=0. The plate is kicked off by the incident shock wave, which then causes the cavitation. The plate gravity and the atmospheric pressure acting on top of the plate cause the velocity to decrease to the extent that the motion is reversed. The cavitation eventually closes as the pressure is restored and the plate rebounds. The cavitation collapse emits a pressure impulse which leads to a reloading of the plate. The solution by the method of characteristics is given in [43] and chosen as the benchmark:(27)pincx,t=Ppeakexp−t+x−xso/cτHt+x−xsoc.Figure 5
Bleich–Sandler configuration.Although the problem was originally in 1D, we use a 3D configuration to validate our in-house framework (Figure5). The plate is modeled in Abaqus by a rigid element and coupled with the fluid solver using MpCCI. To limit the size of the fluid domain, the plane nonreflective boundary is prescribed at the bottom of the fluid column. No boundary condition is defined on the sides since they are symmetry boundaries. The detailed parameters are presented in Table 1. We compare our numerical result with the benchmark solution in Figure 6. The result shows excellent agreement between the in-house solver and benchmark solution with an L2 error norm of 0.0322. The L2 error norm used in this case is defined as(28)L2=∫0ttermVbenchmarkt−Vcomputedt2dt∫0ttermVbenchmarkt2dt,where Vbenchmarkt is the benchmark solution time history, Vcomputedt is the computed solution time history, and tterm is the total simulation time.Table 1
Bleich–Sandler parameter configuration.
Parameters
Value
Plate weight
144 kg
Fluid domain size
0.1 m∗ 3.8 m ∗ 0.1 m
Fluid mesh size
0.025 m iny & z 0.0175 m in x
Shock wave peak pressure (Ppeak)
0.712 MPa
Shock wave decay rate (τ)
0.999 ms
Initial location of the shock front (xso)
0 mFigure 6
Comparison of structural response of the in-house solver and benchmark solution (13 ms).


#### 3.1.2. Floating Shock Platform

In this section, we further assess model and FSI framework in a MILS-901D standard floating shock platform (FSP) [44] application with an incident spherical step exponential shock wave (equation (29)) initiated when the shock wave would contact the nearest structure point to the explosive. The incident shock wave pressure is specified as(29)pincX,t=Ppeakexp−t−R−R0/cτHt−R−R0c,where R is the distance from the charge center to the point X and R0 denotes the initial wave radius at the nearest structure point. The peak pressure Ppeak and decay rate τ are calculated from the similitude equation from [45, 46].Since precise experimental data are not available, we assess our results by comparing with a numerical solution calculated by the Abaqus structural-acoustic solver which couples a FEM-based acoustic fluid solver with the Abaqus structural solver. In Abaqus, we employ a CAFE-like approach using the same cutoff cavitation condition and nonreflective boundary condition as used in our fluid solver. The major difference from the CAFE approach is that the linear acoustic fluid element in Abaqus is pressure-based instead of using displacement potential. The nonmatching mesh interpolation is performed by the node-to-surface constraint function in Abaqus which uses the same principle with the energy-conservative mapping algorithm as described in Section2.4.2. The size of the fluid domain is set so that the bulk cavitation caused by the free surface is captured as shown in the Figure 7. The lower bulk cavitation boundary indicated by the red curve is calculated by the algorithm described in [11] assuming the structure is not present. Since the explosive is located directly below the longitudinal center of the FSP, only half of the fluid and structure domain is modeled. A simulation screenshot is shown in Figure 8. The detailed parameters are listed in Table 2.Figure 7
Problem configuration.Figure 8
FSP Abaqus configuration.Table 2
Abaqus parameter configuration.
Parameters
Value
FSP structure size (half)
4.2672 m (14 ft)∗ 1.2192 m (4 ft) ∗ 4.8768 m (16 ft)
Fluid domain size
5.7912 m (19 ft)∗ 3.6576 (12 ft) ∗ 10.9728 m (36 ft)
Structural element type
4-node shell element (S4R)
Structural element size
0.3048 m (1 ft)
Fluid element type
8-node acoustic brick element (AC3D8R)
Fluid element mesh size
0.0762 m (0.25 ft)
Charge stand-off distance
9.144 m (30 ft)
Charge depth
3.048 m (10 ft)
Shock wave peak pressure
16.22 MPa
Shock wave decay rate
0.3208 msThe same FSP structure is modeled in our framework and coupled to the fluid solver using MpCCI. The same base fluid mesh is used as in Abaqus, but with linear CASE. The results are compared at two points. The first point 1646 (−1.8288 m, 0 m, and −0.3048 m) is at the top of an athwartship girder, and the second point 2320 (−0.3048 m, 0 m, and 0.3048 m) is at the center of a stiffened plate. The simulation result of the in-house CASE solver agrees well with the Abaqus acoustic-structure solver solution as shown in Figure9. In order to quantify this comparison, an error norm is used. The error norm is calculated using the approach introduced by Russell [47] which is frequently used for structural vibration problems. The comprehensive error factor εc is obtained by combining the magnitude error factor εm and phase error factor εp together:(30)εc=π4εm2+εp2.Figure 9
Comparison of structural response of in-house FSI framework and Abaqus/Explicit. (a) Point 1646. (b) Point 2320.
(a)
(b)The expression forεm and εp can be found in [47]. Using the Abaqus results as the baseline, the Russell error norm of the in-house solution is 0.0985 at point 1646 and 0.0982 at point 2330, which is satisfactory. The trend at point 1646 is similar to the Bleich–Sandler result which is identified by an initial kickoff caused by the incident wave, panel rebound, and reloading from cavitation collapse. The velocity response at point 2320 in both models shows a relatively low-frequency large-amplitude oscillation because the plate is unsupported at that location.


### 3.2. Comparison of CAFE and CASE

This section compares the behavior of linear finite elements (linear CAFE), linear spectral elements (linear CASE), and high-order spectral elements (high-order CASE) using our in-house solver in the same FSP case configuration used in Section3.1.2. The result at point 1646 is used for plotting and analysis. We add a damping factor β=0.25 to mitigate the spurious oscillation caused by the cavitation. For all methods, we use a base fluid mesh of 7312 (0.3048 m) (1 ft) cubic hexahedral elements. h refinement (linear CAFE and linear CASE) and p refinement (high-order CASE) are performed on the base fluid mesh. The DOF and time increment are shown in the Table 3. Based on our previous convergence study, an N=8 CASE solution is used as the numerical benchmark to quantify the accuracy of the structural response. We first compare the accuracy of structural response to the benchmark using different types of fluid elements and then assess the behavior by plotting the cavitation distribution. The consistent mapping approach described in Section 2.4.1 is used for the cases throughout this section.Table 3
DOF and time increment (s) of different refinement levels.
Refinement level
DOF
CAFE (s)
Linear CASE (s)
CASE (s)
3
210226
2.39473e − 05
1.59649e − 05
8.70518e − 06
4
490637
1.79605e − 05
1.19737e − 05
5.29987e − 06
5
949336
1.43684e − 05
9.57894e − 06
3.56347e − 06
6
1630195
1.19737e − 05
7.98245e − 06
2.55908e − 06
7
2577086
1.02631e − 05
6.8421e − 06
1.92626e − 06


#### 3.2.1. Operation Count Calculation

An operation count is used to quantify computational cost. Due to the large memory usage of storing the global reactance matrixH, the matrix-vector product (Hψ) is typically conducted on the element level (Heψe) and then assembled. The operation count involved in directly multiplying matrix He and vector ψe (direct-matrix product) is(31)2N+16−N+13,where N is the CASE element order. Only the floating point operations (FLOPs) are counted since the assignment operation, and logical operation is subjective based on different implementations.However, for the high-order CASE element, theON6 operation is very inefficient. Following the algorithm detailed in [29], we employ tensor-product factorization (TF) for the matrix-vector multiplication (Heψe) which features ON4 operation (FLOP shown in the following equation)(32)12N+14+15N+13.We plot the FLOP of an element with CASEp refinement and CAFE/CASE h refinement in Figure 10 using the direct-matrix product and tensor-product factorization. As shown in Figure 10, using tensor-product factorization significantly lowers the FLOP of high-order CASE but does not help with linear CAFE/CASE. Thus, we decide to use direct-matrix product for linear CAFE/CASE and tensor-product factorization for high-order CASE.Figure 10
FLOP comparison.The FLOPs in equations (31) and (32) are used to represent the total operation count of one element in a single time step since the other operations (e.g., boundary term calculation) are marginal in comparison. As a result, the total FLOP for the fluid solver can be represented by(33)ndt∗NEL∗2N+16−N+13,for linear CAFE/CASE and(34)ndt∗NEL∗12N+14+15N+13,for high-order CASE where NEL denotes the total element number and ndt the total number of time steps.


#### 3.2.2. Comparison of Structural Response Accuracy

In this section, we compare the accuracy with the different types of elements to the benchmark using the structural velocity response. The Russell error norm (equation (30)) of the structural response is compared to the degree of freedom (Figure 11). This comparison indicates that the error norm for high-order CASE is smaller than the linear CASE which is smaller than linear CAFE with the same DOF. Furthermore, the error decreases faster with DOF in the high-order CASE case. Since the time increment of high-order CASE is significantly smaller than that of the other element types (Table 3), it is necessary to plot the total FLOP against error norm to quantify their relative computational efficiency.Figure 11
Error norm comparison for different element types with respect to DOF.As shown in Figure12, the high-order CASE with tensor-product factorization achieves the highest accuracy compared to the benchmark with the same number of FLOP. The computational efficiency of high-order CASE is superior even without tensor-product factorization (direct-matrix product). Thus, the high-order CASE may still be advantageous if an unstructured tetrahedral mesh without tensor-product factorization is used. Linear CASE outperforms linear CAFE although inferior to the high-order CASE.Figure 12
Operation count comparison for different element types with respect to FLOP.Next, the behavior of high-order CASE of different orders is investigated.h refinement for CASE with different orders is used for this assessment. The base element size is increased to 0.6096 m (2 ft). A case with 8th-order CASE-type p refinement and 2nd-order h refinement is employed for the benchmark. Russell error factor is employed as in the previous cases. The case configuration and the corresponding DOF are listed in Table 4. The time increments for each case are presented in Table 5. The trend in Figure 13 shows the lower computational cost to achieve the same level of accuracy when the element order is higher.Table 4
DOF of different refinement levels of CASE.
p
=
2
p
=
3
p
=
4
h = 1
—
29497
67881
h = 2
64231
210226
490637
h = 3
222325
735310
1725241
h = 4
490637
1630195
3833881
h = 5
1004541
3348961
—
h = 6
1630195
—
—Table 5
Time increment (s) of different refinement levels of CASE.
p
=
2
p
=
3
p
=
4
h = 1
—
1.74104e − 05
1.05997e − 05
h = 2
1.69333e − 05
8.70518e − 06
5.29987e − 06
h = 3
1.12889e − 05
5.80345e − 06
3.53325e − 06
h = 4
8.46667e − 06
4.35259e − 06
2.64993e − 06
h = 5
6.77333e − 06
3.48207e − 06
—
h = 6
5.64444e − 06
—
—Figure 13
Operation count comparison for different high-order SEMs.


#### 3.2.3. Comparison of Dynamic Condensation Distribution and Cavitation Region

In this section, the influence of cavitation region modeling on the accuracy of predicting structural response in the FSP problem is assessed. Our hypothesis is that the interaction between cavitation and structure plays an important role. We compare the dynamic condensation distribution in the depth direction under the center of the FSP bottom at various time periods (10 ms, 20 ms, and 30 ms) following the approach in [31]. A negative dynamic condensation indicates the extent of cavitation based on equation (12). Instead of the Russell error norm definition which is used for quantifying the norm of the oscillatory response history, we use the L2 error norm definition in(35)L2=∫0dsbenchmarkx,t−scomputedx,t2dx∫0dsbenchmarkx,t2dx,where d is the depth and x is the coordinate in the depth direction. Again, the result of the N=8 CASE solution is used as a benchmark.TheL2 norm for the h and p refinement levels from 4 to 6 are tabulated (Tables 6–8), but only the 4th-level h/p-refined linear CASE, linear CAFE, and high-order CASE of the same DOF are plotted for clarity (Figure 14) showing negative dynamic condensation and hence cavitation as a function of depth. The cavitation region initiates and grows as the FSP panel is displaced inward by the initial shock impact, leaving one cavitated region attached to the structure and the other isolated in the fluid. The two cavitation regions keep growing and eventually merge at 30 ms as the panel is still moving inward.Table 6
L
2 norm of the dynamic condensation distribution at 10 ms.
Refinement level
Linear FEM
Linear SEM
High-order SEM
4
0.7840
0.7048
0.3517
5
0.6725
0.5645
0.2369
6
0.5947
0.4634
0.2759Table 7
L
2 norm of the dynamic condensation distribution at 20 ms.
Refinement level
Linear FEM
Linear SEM
High-order SEM
4
0.7301
0.6352
0.3775
5
0.6391
0.5227
0.2406
6
0.5757
0.4446
0.2015Table 8
L
2 norm of the dynamic condensation distribution at 30 ms.
Refinement level
Linear FEM
Linear SEM
High-order SEM
4
0.6757
0.5948
0.5240
5
0.6108
0.5222
0.2900
6
0.5759
0.4981
0.2085Figure 14
Dynamic condensation distribution at (a) 10 ms, (b) 20 ms, and (c) 30 ms.
(a)
(b)
(c)The error norm results show that high-order CASE gives a more accurate cavitation distribution than linear CASE which is better than linear CAFE, all compared to the benchmark. This reflects the strong association between the accuracy of cavitation modeling and structural response. However, much of the error in the dynamic condensation distribution is not observed in the velocity response of the structure. This is an indication that the disturbance below the first cavitation region may not affect the structure as much as the structural response is driven by the fluid directly attached to it. We may conclude that there is a need to accurately capture the shape and magnitude of cavitation regions in the fluid especially when directly in contact with the structure. Thus, high-order spectral elements may be advantageous in modeling the cavitated surface ship problem due to their superior performance in modeling cavitation.


### 3.3. Comparison of Different Mapping Algorithms for CASE

In Section3.2, it is shown that the computational advantage of high-order CASE is preserved with the consistent mapping algorithm. In this section, three different nonmatching mesh mapping approaches (consistent, average, and energy conservation) are compared to determine which one is the most efficient. The position of the sampling point (pt. 5084) is identical to the point 1646 in Section 3.2 which is on the athwartship stringer in the FSP.


#### 3.3.1. Comparison between the Consistent Mapping and Average Mapping Schemes

In this section, a comparison between the consistent and the average mapping algorithms is made. In order to assess the difference, four fluid mesh configurations using CASE of various orders are used, as shown in Figure15. The base structural mesh stays the same at 0.1524 m (0.5 ft) for all the cases. In the first case, fluid mesh uses a 0.1524 m (0.5 ft) base element with 2nd-order h refinement; the second case uses a 0.3048 m (1 ft) base element with 2nd-order h refinement (matching fluid and structure mesh); the third case uses a 0.1524 m (0.5 ft) base element with 2nd-order p refinement; the fourth case uses a 0.3048 m (1 ft) base element with 2nd-order p refinement. Cases 1 and 3, and cases 2 and 4 have the same degree of freedom.Figure 15
Mesh configuration for cases 1 to 4. (a) Case 1. (b) Case 2. (c) Case 3. (d) Case 4.
(a)
(b)
(c)
(d)Instead of directly comparing the response history, we compare the shock response spectrum (SRS) calculated from the acceleration time history. The SRS is important in describing ship shock response especially in terms of ship equipment design. Using different mapping schemes on the fluid-structure wetted surface may alter the oscillation pattern which may then cause large differences in the SRS. The SRS plots in this section are generated by a MATLAB code created by Irvine [48] using the Smallwood algorithm. An amplification factor Q=10 and the 1 Hz sampling rate are applied. Case 3 with consistent mapping result is used as benchmark since we have already shown that higher-order elements with higher DOF have better accuracy. The SRS for cases 1 to 4 is plotted in Figure 16 from 0.1 to 300 Hz which is a typical upper limit for shipboard shock analysis [49]. The L2 error norm of the four cases is presented in Table 9. The SRS result L2 norm difference between consistent and average mapping for each case is presented in Table 10.Figure 16
Shock response spectrum of cases 1 to 4. (a) Case 1. (b) Case 2. (c) Case 3. (d) Case 4.
(a)
(b)
(c)
(d)Table 9
L
2 error norm for cases 1 to 4.
Average mapping
Consistent mapping
Case 1
0.02830
0.01567
Case 2
0.04442
0.03863
Case 3
0.03245
Benchmark
Case 4
0.04577
0.02618Table 10
Difference between consistent and average mapping (L2 norm).
0.1524 m (0.5 ft)
0.3048 m (1 ft)
h = 2
0.01319 (case 1)
0.008317 (case 3)
p
=
2
0.03245 (case 3)
0.06684 (case 4)TheL2 error norm in Table 9 indicates that the consistent mapping algorithm gives a better (lower error norm) result than average mapping in the cases tested. This result makes sense because the displacement/pressure value on quadrature points is more accurate in the consistent mapping case. Average mapping generally overpredicts the peak velocity in the shock response spectrum as shown in Figure 16.If we proceed with the conclusion that consistent mapping gives better results in all the cases, we can treat theL2 norm in Table 10 as the extra mapping error introduced by the average mapping since it is essentially a simplification of the consistent mapping. For the matching mesh configuration (case 2), the difference between consistent and average mapping is the smallest. The difference becomes larger as the fluid element size is increased (case 3 vs case 4). The high-order CASE show a larger difference than the linear element with the same DOF (case 1 vs case 3). Overall, the performance of average mapping becomes poorer for fluid elements with larger size and higher order, which is the case for high-order CASE.


#### 3.3.2. Comparison between Conservative and Consistent Schemes

In this section, the behavior of the consistent mapping scheme is compared to the energy-conservative scheme to examine the importance of energy conservation. The difference between the two algorithms is that energy conservation is not guaranteed by the consistent approach as demonstrated in Section2.4. However, energy conservation may not be important for a short-duration problem such as in far-field UNDEX which usually lasts less than 100 ms. To demonstrate the problem, a 3rd-order CASE fluid mesh coupled with the linear structural mesh of the same element size 0.1524 m (0.5 ft) as presented in Figure 3 is used. For the consistent mapping, we use linear Gauss–Legendre quadrature on the structural element which guarantees total force conservation but is theoretically unable to achieve energy conservation. Although consistency is normally violated for the conservative scheme [50], the nested fluid mesh configuration used in this case guarantees consistency by passing the aforementioned patch test as described in Section 1.3.The SRS comparison in Figure17 shows that there is almost no difference in structural response between the consistent and conservative method. We track the total energy input (positive)/received (negative) by fluid and received/input by the structure in the consistent method. The energy history in Figure 18 indicates that the consistent mapping approach essentially conserves the total energy in this application. This supports the assertion made in [41] that an energy conservation scheme might be of little significance in short-duration problems involving shock. Also, given the knowledge that an energy-conservative method can violate consistency if the fluid elements are fully mismatched (not nested), the consistent mapping approach is sufficient and possibly superior for the far-field UNDEX application.Figure 17
Comparison between consistent and energy-conservative approach in SRS.Figure 18
Energy history on fluid and structure interface.


## 3.1. Model Assessment

In this section, the model and FSI framework introduced in Section2 are assessed using two benchmark problems. The linear CASE is used in the in-house fluid solver, and the consistent mapping algorithm is adopted for the nonmatching mesh interpolation. Excellent agreement with the baseline results provides more confidence in the model setup employed for the numerical experiment described in later sections.


### 3.1.1. Bleich–Sandler Problem

The Bleich–Sandler plate [43] is a classic UNDEX benchmark problem with cavitation. The problem includes a rigid plate sitting on the free surface of an infinitely deep fluid column (Figure 5). The plate is subject to the load from a step exponential plane shock wave in equation (27) where Ppeak is the peak pressure on the shock front, c is the speed of sound, τ is the decay constant, and xso is the initial location of the shock front at t=0. The plate is kicked off by the incident shock wave, which then causes the cavitation. The plate gravity and the atmospheric pressure acting on top of the plate cause the velocity to decrease to the extent that the motion is reversed. The cavitation eventually closes as the pressure is restored and the plate rebounds. The cavitation collapse emits a pressure impulse which leads to a reloading of the plate. The solution by the method of characteristics is given in [43] and chosen as the benchmark:(27)pincx,t=Ppeakexp−t+x−xso/cτHt+x−xsoc.Figure 5
Bleich–Sandler configuration.Although the problem was originally in 1D, we use a 3D configuration to validate our in-house framework (Figure5). The plate is modeled in Abaqus by a rigid element and coupled with the fluid solver using MpCCI. To limit the size of the fluid domain, the plane nonreflective boundary is prescribed at the bottom of the fluid column. No boundary condition is defined on the sides since they are symmetry boundaries. The detailed parameters are presented in Table 1. We compare our numerical result with the benchmark solution in Figure 6. The result shows excellent agreement between the in-house solver and benchmark solution with an L2 error norm of 0.0322. The L2 error norm used in this case is defined as(28)L2=∫0ttermVbenchmarkt−Vcomputedt2dt∫0ttermVbenchmarkt2dt,where Vbenchmarkt is the benchmark solution time history, Vcomputedt is the computed solution time history, and tterm is the total simulation time.Table 1
Bleich–Sandler parameter configuration.
Parameters
Value
Plate weight
144 kg
Fluid domain size
0.1 m∗ 3.8 m ∗ 0.1 m
Fluid mesh size
0.025 m iny & z 0.0175 m in x
Shock wave peak pressure (Ppeak)
0.712 MPa
Shock wave decay rate (τ)
0.999 ms
Initial location of the shock front (xso)
0 mFigure 6
Comparison of structural response of the in-house solver and benchmark solution (13 ms).


### 3.1.2. Floating Shock Platform

In this section, we further assess model and FSI framework in a MILS-901D standard floating shock platform (FSP) [44] application with an incident spherical step exponential shock wave (equation (29)) initiated when the shock wave would contact the nearest structure point to the explosive. The incident shock wave pressure is specified as(29)pincX,t=Ppeakexp−t−R−R0/cτHt−R−R0c,where R is the distance from the charge center to the point X and R0 denotes the initial wave radius at the nearest structure point. The peak pressure Ppeak and decay rate τ are calculated from the similitude equation from [45, 46].Since precise experimental data are not available, we assess our results by comparing with a numerical solution calculated by the Abaqus structural-acoustic solver which couples a FEM-based acoustic fluid solver with the Abaqus structural solver. In Abaqus, we employ a CAFE-like approach using the same cutoff cavitation condition and nonreflective boundary condition as used in our fluid solver. The major difference from the CAFE approach is that the linear acoustic fluid element in Abaqus is pressure-based instead of using displacement potential. The nonmatching mesh interpolation is performed by the node-to-surface constraint function in Abaqus which uses the same principle with the energy-conservative mapping algorithm as described in Section2.4.2. The size of the fluid domain is set so that the bulk cavitation caused by the free surface is captured as shown in the Figure 7. The lower bulk cavitation boundary indicated by the red curve is calculated by the algorithm described in [11] assuming the structure is not present. Since the explosive is located directly below the longitudinal center of the FSP, only half of the fluid and structure domain is modeled. A simulation screenshot is shown in Figure 8. The detailed parameters are listed in Table 2.Figure 7
Problem configuration.Figure 8
FSP Abaqus configuration.Table 2
Abaqus parameter configuration.
Parameters
Value
FSP structure size (half)
4.2672 m (14 ft)∗ 1.2192 m (4 ft) ∗ 4.8768 m (16 ft)
Fluid domain size
5.7912 m (19 ft)∗ 3.6576 (12 ft) ∗ 10.9728 m (36 ft)
Structural element type
4-node shell element (S4R)
Structural element size
0.3048 m (1 ft)
Fluid element type
8-node acoustic brick element (AC3D8R)
Fluid element mesh size
0.0762 m (0.25 ft)
Charge stand-off distance
9.144 m (30 ft)
Charge depth
3.048 m (10 ft)
Shock wave peak pressure
16.22 MPa
Shock wave decay rate
0.3208 msThe same FSP structure is modeled in our framework and coupled to the fluid solver using MpCCI. The same base fluid mesh is used as in Abaqus, but with linear CASE. The results are compared at two points. The first point 1646 (−1.8288 m, 0 m, and −0.3048 m) is at the top of an athwartship girder, and the second point 2320 (−0.3048 m, 0 m, and 0.3048 m) is at the center of a stiffened plate. The simulation result of the in-house CASE solver agrees well with the Abaqus acoustic-structure solver solution as shown in Figure9. In order to quantify this comparison, an error norm is used. The error norm is calculated using the approach introduced by Russell [47] which is frequently used for structural vibration problems. The comprehensive error factor εc is obtained by combining the magnitude error factor εm and phase error factor εp together:(30)εc=π4εm2+εp2.Figure 9
Comparison of structural response of in-house FSI framework and Abaqus/Explicit. (a) Point 1646. (b) Point 2320.
(a)
(b)The expression forεm and εp can be found in [47]. Using the Abaqus results as the baseline, the Russell error norm of the in-house solution is 0.0985 at point 1646 and 0.0982 at point 2330, which is satisfactory. The trend at point 1646 is similar to the Bleich–Sandler result which is identified by an initial kickoff caused by the incident wave, panel rebound, and reloading from cavitation collapse. The velocity response at point 2320 in both models shows a relatively low-frequency large-amplitude oscillation because the plate is unsupported at that location.


## 3.1.1. Bleich–Sandler Problem

The Bleich–Sandler plate [43] is a classic UNDEX benchmark problem with cavitation. The problem includes a rigid plate sitting on the free surface of an infinitely deep fluid column (Figure 5). The plate is subject to the load from a step exponential plane shock wave in equation (27) where Ppeak is the peak pressure on the shock front, c is the speed of sound, τ is the decay constant, and xso is the initial location of the shock front at t=0. The plate is kicked off by the incident shock wave, which then causes the cavitation. The plate gravity and the atmospheric pressure acting on top of the plate cause the velocity to decrease to the extent that the motion is reversed. The cavitation eventually closes as the pressure is restored and the plate rebounds. The cavitation collapse emits a pressure impulse which leads to a reloading of the plate. The solution by the method of characteristics is given in [43] and chosen as the benchmark:(27)pincx,t=Ppeakexp−t+x−xso/cτHt+x−xsoc.Figure 5
Bleich–Sandler configuration.Although the problem was originally in 1D, we use a 3D configuration to validate our in-house framework (Figure5). The plate is modeled in Abaqus by a rigid element and coupled with the fluid solver using MpCCI. To limit the size of the fluid domain, the plane nonreflective boundary is prescribed at the bottom of the fluid column. No boundary condition is defined on the sides since they are symmetry boundaries. The detailed parameters are presented in Table 1. We compare our numerical result with the benchmark solution in Figure 6. The result shows excellent agreement between the in-house solver and benchmark solution with an L2 error norm of 0.0322. The L2 error norm used in this case is defined as(28)L2=∫0ttermVbenchmarkt−Vcomputedt2dt∫0ttermVbenchmarkt2dt,where Vbenchmarkt is the benchmark solution time history, Vcomputedt is the computed solution time history, and tterm is the total simulation time.Table 1
Bleich–Sandler parameter configuration.
Parameters
Value
Plate weight
144 kg
Fluid domain size
0.1 m∗ 3.8 m ∗ 0.1 m
Fluid mesh size
0.025 m iny & z 0.0175 m in x
Shock wave peak pressure (Ppeak)
0.712 MPa
Shock wave decay rate (τ)
0.999 ms
Initial location of the shock front (xso)
0 mFigure 6
Comparison of structural response of the in-house solver and benchmark solution (13 ms).


## 3.1.2. Floating Shock Platform

In this section, we further assess model and FSI framework in a MILS-901D standard floating shock platform (FSP) [44] application with an incident spherical step exponential shock wave (equation (29)) initiated when the shock wave would contact the nearest structure point to the explosive. The incident shock wave pressure is specified as(29)pincX,t=Ppeakexp−t−R−R0/cτHt−R−R0c,where R is the distance from the charge center to the point X and R0 denotes the initial wave radius at the nearest structure point. The peak pressure Ppeak and decay rate τ are calculated from the similitude equation from [45, 46].Since precise experimental data are not available, we assess our results by comparing with a numerical solution calculated by the Abaqus structural-acoustic solver which couples a FEM-based acoustic fluid solver with the Abaqus structural solver. In Abaqus, we employ a CAFE-like approach using the same cutoff cavitation condition and nonreflective boundary condition as used in our fluid solver. The major difference from the CAFE approach is that the linear acoustic fluid element in Abaqus is pressure-based instead of using displacement potential. The nonmatching mesh interpolation is performed by the node-to-surface constraint function in Abaqus which uses the same principle with the energy-conservative mapping algorithm as described in Section2.4.2. The size of the fluid domain is set so that the bulk cavitation caused by the free surface is captured as shown in the Figure 7. The lower bulk cavitation boundary indicated by the red curve is calculated by the algorithm described in [11] assuming the structure is not present. Since the explosive is located directly below the longitudinal center of the FSP, only half of the fluid and structure domain is modeled. A simulation screenshot is shown in Figure 8. The detailed parameters are listed in Table 2.Figure 7
Problem configuration.Figure 8
FSP Abaqus configuration.Table 2
Abaqus parameter configuration.
Parameters
Value
FSP structure size (half)
4.2672 m (14 ft)∗ 1.2192 m (4 ft) ∗ 4.8768 m (16 ft)
Fluid domain size
5.7912 m (19 ft)∗ 3.6576 (12 ft) ∗ 10.9728 m (36 ft)
Structural element type
4-node shell element (S4R)
Structural element size
0.3048 m (1 ft)
Fluid element type
8-node acoustic brick element (AC3D8R)
Fluid element mesh size
0.0762 m (0.25 ft)
Charge stand-off distance
9.144 m (30 ft)
Charge depth
3.048 m (10 ft)
Shock wave peak pressure
16.22 MPa
Shock wave decay rate
0.3208 msThe same FSP structure is modeled in our framework and coupled to the fluid solver using MpCCI. The same base fluid mesh is used as in Abaqus, but with linear CASE. The results are compared at two points. The first point 1646 (−1.8288 m, 0 m, and −0.3048 m) is at the top of an athwartship girder, and the second point 2320 (−0.3048 m, 0 m, and 0.3048 m) is at the center of a stiffened plate. The simulation result of the in-house CASE solver agrees well with the Abaqus acoustic-structure solver solution as shown in Figure9. In order to quantify this comparison, an error norm is used. The error norm is calculated using the approach introduced by Russell [47] which is frequently used for structural vibration problems. The comprehensive error factor εc is obtained by combining the magnitude error factor εm and phase error factor εp together:(30)εc=π4εm2+εp2.Figure 9
Comparison of structural response of in-house FSI framework and Abaqus/Explicit. (a) Point 1646. (b) Point 2320.
(a)
(b)The expression forεm and εp can be found in [47]. Using the Abaqus results as the baseline, the Russell error norm of the in-house solution is 0.0985 at point 1646 and 0.0982 at point 2330, which is satisfactory. The trend at point 1646 is similar to the Bleich–Sandler result which is identified by an initial kickoff caused by the incident wave, panel rebound, and reloading from cavitation collapse. The velocity response at point 2320 in both models shows a relatively low-frequency large-amplitude oscillation because the plate is unsupported at that location.


## 3.2. Comparison of CAFE and CASE

This section compares the behavior of linear finite elements (linear CAFE), linear spectral elements (linear CASE), and high-order spectral elements (high-order CASE) using our in-house solver in the same FSP case configuration used in Section3.1.2. The result at point 1646 is used for plotting and analysis. We add a damping factor β=0.25 to mitigate the spurious oscillation caused by the cavitation. For all methods, we use a base fluid mesh of 7312 (0.3048 m) (1 ft) cubic hexahedral elements. h refinement (linear CAFE and linear CASE) and p refinement (high-order CASE) are performed on the base fluid mesh. The DOF and time increment are shown in the Table 3. Based on our previous convergence study, an N=8 CASE solution is used as the numerical benchmark to quantify the accuracy of the structural response. We first compare the accuracy of structural response to the benchmark using different types of fluid elements and then assess the behavior by plotting the cavitation distribution. The consistent mapping approach described in Section 2.4.1 is used for the cases throughout this section.Table 3
DOF and time increment (s) of different refinement levels.
Refinement level
DOF
CAFE (s)
Linear CASE (s)
CASE (s)
3
210226
2.39473e − 05
1.59649e − 05
8.70518e − 06
4
490637
1.79605e − 05
1.19737e − 05
5.29987e − 06
5
949336
1.43684e − 05
9.57894e − 06
3.56347e − 06
6
1630195
1.19737e − 05
7.98245e − 06
2.55908e − 06
7
2577086
1.02631e − 05
6.8421e − 06
1.92626e − 06


### 3.2.1. Operation Count Calculation

An operation count is used to quantify computational cost. Due to the large memory usage of storing the global reactance matrixH, the matrix-vector product (Hψ) is typically conducted on the element level (Heψe) and then assembled. The operation count involved in directly multiplying matrix He and vector ψe (direct-matrix product) is(31)2N+16−N+13,where N is the CASE element order. Only the floating point operations (FLOPs) are counted since the assignment operation, and logical operation is subjective based on different implementations.However, for the high-order CASE element, theON6 operation is very inefficient. Following the algorithm detailed in [29], we employ tensor-product factorization (TF) for the matrix-vector multiplication (Heψe) which features ON4 operation (FLOP shown in the following equation)(32)12N+14+15N+13.We plot the FLOP of an element with CASEp refinement and CAFE/CASE h refinement in Figure 10 using the direct-matrix product and tensor-product factorization. As shown in Figure 10, using tensor-product factorization significantly lowers the FLOP of high-order CASE but does not help with linear CAFE/CASE. Thus, we decide to use direct-matrix product for linear CAFE/CASE and tensor-product factorization for high-order CASE.Figure 10
FLOP comparison.The FLOPs in equations (31) and (32) are used to represent the total operation count of one element in a single time step since the other operations (e.g., boundary term calculation) are marginal in comparison. As a result, the total FLOP for the fluid solver can be represented by(33)ndt∗NEL∗2N+16−N+13,for linear CAFE/CASE and(34)ndt∗NEL∗12N+14+15N+13,for high-order CASE where NEL denotes the total element number and ndt the total number of time steps.


### 3.2.2. Comparison of Structural Response Accuracy

In this section, we compare the accuracy with the different types of elements to the benchmark using the structural velocity response. The Russell error norm (equation (30)) of the structural response is compared to the degree of freedom (Figure 11). This comparison indicates that the error norm for high-order CASE is smaller than the linear CASE which is smaller than linear CAFE with the same DOF. Furthermore, the error decreases faster with DOF in the high-order CASE case. Since the time increment of high-order CASE is significantly smaller than that of the other element types (Table 3), it is necessary to plot the total FLOP against error norm to quantify their relative computational efficiency.Figure 11
Error norm comparison for different element types with respect to DOF.As shown in Figure12, the high-order CASE with tensor-product factorization achieves the highest accuracy compared to the benchmark with the same number of FLOP. The computational efficiency of high-order CASE is superior even without tensor-product factorization (direct-matrix product). Thus, the high-order CASE may still be advantageous if an unstructured tetrahedral mesh without tensor-product factorization is used. Linear CASE outperforms linear CAFE although inferior to the high-order CASE.Figure 12
Operation count comparison for different element types with respect to FLOP.Next, the behavior of high-order CASE of different orders is investigated.h refinement for CASE with different orders is used for this assessment. The base element size is increased to 0.6096 m (2 ft). A case with 8th-order CASE-type p refinement and 2nd-order h refinement is employed for the benchmark. Russell error factor is employed as in the previous cases. The case configuration and the corresponding DOF are listed in Table 4. The time increments for each case are presented in Table 5. The trend in Figure 13 shows the lower computational cost to achieve the same level of accuracy when the element order is higher.Table 4
DOF of different refinement levels of CASE.
p
=
2
p
=
3
p
=
4
h = 1
—
29497
67881
h = 2
64231
210226
490637
h = 3
222325
735310
1725241
h = 4
490637
1630195
3833881
h = 5
1004541
3348961
—
h = 6
1630195
—
—Table 5
Time increment (s) of different refinement levels of CASE.
p
=
2
p
=
3
p
=
4
h = 1
—
1.74104e − 05
1.05997e − 05
h = 2
1.69333e − 05
8.70518e − 06
5.29987e − 06
h = 3
1.12889e − 05
5.80345e − 06
3.53325e − 06
h = 4
8.46667e − 06
4.35259e − 06
2.64993e − 06
h = 5
6.77333e − 06
3.48207e − 06
—
h = 6
5.64444e − 06
—
—Figure 13
Operation count comparison for different high-order SEMs.


### 3.2.3. Comparison of Dynamic Condensation Distribution and Cavitation Region

In this section, the influence of cavitation region modeling on the accuracy of predicting structural response in the FSP problem is assessed. Our hypothesis is that the interaction between cavitation and structure plays an important role. We compare the dynamic condensation distribution in the depth direction under the center of the FSP bottom at various time periods (10 ms, 20 ms, and 30 ms) following the approach in [31]. A negative dynamic condensation indicates the extent of cavitation based on equation (12). Instead of the Russell error norm definition which is used for quantifying the norm of the oscillatory response history, we use the L2 error norm definition in(35)L2=∫0dsbenchmarkx,t−scomputedx,t2dx∫0dsbenchmarkx,t2dx,where d is the depth and x is the coordinate in the depth direction. Again, the result of the N=8 CASE solution is used as a benchmark.TheL2 norm for the h and p refinement levels from 4 to 6 are tabulated (Tables 6–8), but only the 4th-level h/p-refined linear CASE, linear CAFE, and high-order CASE of the same DOF are plotted for clarity (Figure 14) showing negative dynamic condensation and hence cavitation as a function of depth. The cavitation region initiates and grows as the FSP panel is displaced inward by the initial shock impact, leaving one cavitated region attached to the structure and the other isolated in the fluid. The two cavitation regions keep growing and eventually merge at 30 ms as the panel is still moving inward.Table 6
L
2 norm of the dynamic condensation distribution at 10 ms.
Refinement level
Linear FEM
Linear SEM
High-order SEM
4
0.7840
0.7048
0.3517
5
0.6725
0.5645
0.2369
6
0.5947
0.4634
0.2759Table 7
L
2 norm of the dynamic condensation distribution at 20 ms.
Refinement level
Linear FEM
Linear SEM
High-order SEM
4
0.7301
0.6352
0.3775
5
0.6391
0.5227
0.2406
6
0.5757
0.4446
0.2015Table 8
L
2 norm of the dynamic condensation distribution at 30 ms.
Refinement level
Linear FEM
Linear SEM
High-order SEM
4
0.6757
0.5948
0.5240
5
0.6108
0.5222
0.2900
6
0.5759
0.4981
0.2085Figure 14
Dynamic condensation distribution at (a) 10 ms, (b) 20 ms, and (c) 30 ms.
(a)
(b)
(c)The error norm results show that high-order CASE gives a more accurate cavitation distribution than linear CASE which is better than linear CAFE, all compared to the benchmark. This reflects the strong association between the accuracy of cavitation modeling and structural response. However, much of the error in the dynamic condensation distribution is not observed in the velocity response of the structure. This is an indication that the disturbance below the first cavitation region may not affect the structure as much as the structural response is driven by the fluid directly attached to it. We may conclude that there is a need to accurately capture the shape and magnitude of cavitation regions in the fluid especially when directly in contact with the structure. Thus, high-order spectral elements may be advantageous in modeling the cavitated surface ship problem due to their superior performance in modeling cavitation.


## 3.2.1. Operation Count Calculation

An operation count is used to quantify computational cost. Due to the large memory usage of storing the global reactance matrixH, the matrix-vector product (Hψ) is typically conducted on the element level (Heψe) and then assembled. The operation count involved in directly multiplying matrix He and vector ψe (direct-matrix product) is(31)2N+16−N+13,where N is the CASE element order. Only the floating point operations (FLOPs) are counted since the assignment operation, and logical operation is subjective based on different implementations.However, for the high-order CASE element, theON6 operation is very inefficient. Following the algorithm detailed in [29], we employ tensor-product factorization (TF) for the matrix-vector multiplication (Heψe) which features ON4 operation (FLOP shown in the following equation)(32)12N+14+15N+13.We plot the FLOP of an element with CASEp refinement and CAFE/CASE h refinement in Figure 10 using the direct-matrix product and tensor-product factorization. As shown in Figure 10, using tensor-product factorization significantly lowers the FLOP of high-order CASE but does not help with linear CAFE/CASE. Thus, we decide to use direct-matrix product for linear CAFE/CASE and tensor-product factorization for high-order CASE.Figure 10
FLOP comparison.The FLOPs in equations (31) and (32) are used to represent the total operation count of one element in a single time step since the other operations (e.g., boundary term calculation) are marginal in comparison. As a result, the total FLOP for the fluid solver can be represented by(33)ndt∗NEL∗2N+16−N+13,for linear CAFE/CASE and(34)ndt∗NEL∗12N+14+15N+13,for high-order CASE where NEL denotes the total element number and ndt the total number of time steps.


## 3.2.2. Comparison of Structural Response Accuracy

In this section, we compare the accuracy with the different types of elements to the benchmark using the structural velocity response. The Russell error norm (equation (30)) of the structural response is compared to the degree of freedom (Figure 11). This comparison indicates that the error norm for high-order CASE is smaller than the linear CASE which is smaller than linear CAFE with the same DOF. Furthermore, the error decreases faster with DOF in the high-order CASE case. Since the time increment of high-order CASE is significantly smaller than that of the other element types (Table 3), it is necessary to plot the total FLOP against error norm to quantify their relative computational efficiency.Figure 11
Error norm comparison for different element types with respect to DOF.As shown in Figure12, the high-order CASE with tensor-product factorization achieves the highest accuracy compared to the benchmark with the same number of FLOP. The computational efficiency of high-order CASE is superior even without tensor-product factorization (direct-matrix product). Thus, the high-order CASE may still be advantageous if an unstructured tetrahedral mesh without tensor-product factorization is used. Linear CASE outperforms linear CAFE although inferior to the high-order CASE.Figure 12
Operation count comparison for different element types with respect to FLOP.Next, the behavior of high-order CASE of different orders is investigated.h refinement for CASE with different orders is used for this assessment. The base element size is increased to 0.6096 m (2 ft). A case with 8th-order CASE-type p refinement and 2nd-order h refinement is employed for the benchmark. Russell error factor is employed as in the previous cases. The case configuration and the corresponding DOF are listed in Table 4. The time increments for each case are presented in Table 5. The trend in Figure 13 shows the lower computational cost to achieve the same level of accuracy when the element order is higher.Table 4
DOF of different refinement levels of CASE.
p
=
2
p
=
3
p
=
4
h = 1
—
29497
67881
h = 2
64231
210226
490637
h = 3
222325
735310
1725241
h = 4
490637
1630195
3833881
h = 5
1004541
3348961
—
h = 6
1630195
—
—Table 5
Time increment (s) of different refinement levels of CASE.
p
=
2
p
=
3
p
=
4
h = 1
—
1.74104e − 05
1.05997e − 05
h = 2
1.69333e − 05
8.70518e − 06
5.29987e − 06
h = 3
1.12889e − 05
5.80345e − 06
3.53325e − 06
h = 4
8.46667e − 06
4.35259e − 06
2.64993e − 06
h = 5
6.77333e − 06
3.48207e − 06
—
h = 6
5.64444e − 06
—
—Figure 13
Operation count comparison for different high-order SEMs.


## 3.2.3. Comparison of Dynamic Condensation Distribution and Cavitation Region

In this section, the influence of cavitation region modeling on the accuracy of predicting structural response in the FSP problem is assessed. Our hypothesis is that the interaction between cavitation and structure plays an important role. We compare the dynamic condensation distribution in the depth direction under the center of the FSP bottom at various time periods (10 ms, 20 ms, and 30 ms) following the approach in [31]. A negative dynamic condensation indicates the extent of cavitation based on equation (12). Instead of the Russell error norm definition which is used for quantifying the norm of the oscillatory response history, we use the L2 error norm definition in(35)L2=∫0dsbenchmarkx,t−scomputedx,t2dx∫0dsbenchmarkx,t2dx,where d is the depth and x is the coordinate in the depth direction. Again, the result of the N=8 CASE solution is used as a benchmark.TheL2 norm for the h and p refinement levels from 4 to 6 are tabulated (Tables 6–8), but only the 4th-level h/p-refined linear CASE, linear CAFE, and high-order CASE of the same DOF are plotted for clarity (Figure 14) showing negative dynamic condensation and hence cavitation as a function of depth. The cavitation region initiates and grows as the FSP panel is displaced inward by the initial shock impact, leaving one cavitated region attached to the structure and the other isolated in the fluid. The two cavitation regions keep growing and eventually merge at 30 ms as the panel is still moving inward.Table 6
L
2 norm of the dynamic condensation distribution at 10 ms.
Refinement level
Linear FEM
Linear SEM
High-order SEM
4
0.7840
0.7048
0.3517
5
0.6725
0.5645
0.2369
6
0.5947
0.4634
0.2759Table 7
L
2 norm of the dynamic condensation distribution at 20 ms.
Refinement level
Linear FEM
Linear SEM
High-order SEM
4
0.7301
0.6352
0.3775
5
0.6391
0.5227
0.2406
6
0.5757
0.4446
0.2015Table 8
L
2 norm of the dynamic condensation distribution at 30 ms.
Refinement level
Linear FEM
Linear SEM
High-order SEM
4
0.6757
0.5948
0.5240
5
0.6108
0.5222
0.2900
6
0.5759
0.4981
0.2085Figure 14
Dynamic condensation distribution at (a) 10 ms, (b) 20 ms, and (c) 30 ms.
(a)
(b)
(c)The error norm results show that high-order CASE gives a more accurate cavitation distribution than linear CASE which is better than linear CAFE, all compared to the benchmark. This reflects the strong association between the accuracy of cavitation modeling and structural response. However, much of the error in the dynamic condensation distribution is not observed in the velocity response of the structure. This is an indication that the disturbance below the first cavitation region may not affect the structure as much as the structural response is driven by the fluid directly attached to it. We may conclude that there is a need to accurately capture the shape and magnitude of cavitation regions in the fluid especially when directly in contact with the structure. Thus, high-order spectral elements may be advantageous in modeling the cavitated surface ship problem due to their superior performance in modeling cavitation.


## 3.3. Comparison of Different Mapping Algorithms for CASE

In Section3.2, it is shown that the computational advantage of high-order CASE is preserved with the consistent mapping algorithm. In this section, three different nonmatching mesh mapping approaches (consistent, average, and energy conservation) are compared to determine which one is the most efficient. The position of the sampling point (pt. 5084) is identical to the point 1646 in Section 3.2 which is on the athwartship stringer in the FSP.


### 3.3.1. Comparison between the Consistent Mapping and Average Mapping Schemes

In this section, a comparison between the consistent and the average mapping algorithms is made. In order to assess the difference, four fluid mesh configurations using CASE of various orders are used, as shown in Figure15. The base structural mesh stays the same at 0.1524 m (0.5 ft) for all the cases. In the first case, fluid mesh uses a 0.1524 m (0.5 ft) base element with 2nd-order h refinement; the second case uses a 0.3048 m (1 ft) base element with 2nd-order h refinement (matching fluid and structure mesh); the third case uses a 0.1524 m (0.5 ft) base element with 2nd-order p refinement; the fourth case uses a 0.3048 m (1 ft) base element with 2nd-order p refinement. Cases 1 and 3, and cases 2 and 4 have the same degree of freedom.Figure 15
Mesh configuration for cases 1 to 4. (a) Case 1. (b) Case 2. (c) Case 3. (d) Case 4.
(a)
(b)
(c)
(d)Instead of directly comparing the response history, we compare the shock response spectrum (SRS) calculated from the acceleration time history. The SRS is important in describing ship shock response especially in terms of ship equipment design. Using different mapping schemes on the fluid-structure wetted surface may alter the oscillation pattern which may then cause large differences in the SRS. The SRS plots in this section are generated by a MATLAB code created by Irvine [48] using the Smallwood algorithm. An amplification factor Q=10 and the 1 Hz sampling rate are applied. Case 3 with consistent mapping result is used as benchmark since we have already shown that higher-order elements with higher DOF have better accuracy. The SRS for cases 1 to 4 is plotted in Figure 16 from 0.1 to 300 Hz which is a typical upper limit for shipboard shock analysis [49]. The L2 error norm of the four cases is presented in Table 9. The SRS result L2 norm difference between consistent and average mapping for each case is presented in Table 10.Figure 16
Shock response spectrum of cases 1 to 4. (a) Case 1. (b) Case 2. (c) Case 3. (d) Case 4.
(a)
(b)
(c)
(d)Table 9
L
2 error norm for cases 1 to 4.
Average mapping
Consistent mapping
Case 1
0.02830
0.01567
Case 2
0.04442
0.03863
Case 3
0.03245
Benchmark
Case 4
0.04577
0.02618Table 10
Difference between consistent and average mapping (L2 norm).
0.1524 m (0.5 ft)
0.3048 m (1 ft)
h = 2
0.01319 (case 1)
0.008317 (case 3)
p
=
2
0.03245 (case 3)
0.06684 (case 4)TheL2 error norm in Table 9 indicates that the consistent mapping algorithm gives a better (lower error norm) result than average mapping in the cases tested. This result makes sense because the displacement/pressure value on quadrature points is more accurate in the consistent mapping case. Average mapping generally overpredicts the peak velocity in the shock response spectrum as shown in Figure 16.If we proceed with the conclusion that consistent mapping gives better results in all the cases, we can treat theL2 norm in Table 10 as the extra mapping error introduced by the average mapping since it is essentially a simplification of the consistent mapping. For the matching mesh configuration (case 2), the difference between consistent and average mapping is the smallest. The difference becomes larger as the fluid element size is increased (case 3 vs case 4). The high-order CASE show a larger difference than the linear element with the same DOF (case 1 vs case 3). Overall, the performance of average mapping becomes poorer for fluid elements with larger size and higher order, which is the case for high-order CASE.


### 3.3.2. Comparison between Conservative and Consistent Schemes

In this section, the behavior of the consistent mapping scheme is compared to the energy-conservative scheme to examine the importance of energy conservation. The difference between the two algorithms is that energy conservation is not guaranteed by the consistent approach as demonstrated in Section2.4. However, energy conservation may not be important for a short-duration problem such as in far-field UNDEX which usually lasts less than 100 ms. To demonstrate the problem, a 3rd-order CASE fluid mesh coupled with the linear structural mesh of the same element size 0.1524 m (0.5 ft) as presented in Figure 3 is used. For the consistent mapping, we use linear Gauss–Legendre quadrature on the structural element which guarantees total force conservation but is theoretically unable to achieve energy conservation. Although consistency is normally violated for the conservative scheme [50], the nested fluid mesh configuration used in this case guarantees consistency by passing the aforementioned patch test as described in Section 1.3.The SRS comparison in Figure17 shows that there is almost no difference in structural response between the consistent and conservative method. We track the total energy input (positive)/received (negative) by fluid and received/input by the structure in the consistent method. The energy history in Figure 18 indicates that the consistent mapping approach essentially conserves the total energy in this application. This supports the assertion made in [41] that an energy conservation scheme might be of little significance in short-duration problems involving shock. Also, given the knowledge that an energy-conservative method can violate consistency if the fluid elements are fully mismatched (not nested), the consistent mapping approach is sufficient and possibly superior for the far-field UNDEX application.Figure 17
Comparison between consistent and energy-conservative approach in SRS.Figure 18
Energy history on fluid and structure interface.


## 3.3.1. Comparison between the Consistent Mapping and Average Mapping Schemes

In this section, a comparison between the consistent and the average mapping algorithms is made. In order to assess the difference, four fluid mesh configurations using CASE of various orders are used, as shown in Figure15. The base structural mesh stays the same at 0.1524 m (0.5 ft) for all the cases. In the first case, fluid mesh uses a 0.1524 m (0.5 ft) base element with 2nd-order h refinement; the second case uses a 0.3048 m (1 ft) base element with 2nd-order h refinement (matching fluid and structure mesh); the third case uses a 0.1524 m (0.5 ft) base element with 2nd-order p refinement; the fourth case uses a 0.3048 m (1 ft) base element with 2nd-order p refinement. Cases 1 and 3, and cases 2 and 4 have the same degree of freedom.Figure 15
Mesh configuration for cases 1 to 4. (a) Case 1. (b) Case 2. (c) Case 3. (d) Case 4.
(a)
(b)
(c)
(d)Instead of directly comparing the response history, we compare the shock response spectrum (SRS) calculated from the acceleration time history. The SRS is important in describing ship shock response especially in terms of ship equipment design. Using different mapping schemes on the fluid-structure wetted surface may alter the oscillation pattern which may then cause large differences in the SRS. The SRS plots in this section are generated by a MATLAB code created by Irvine [48] using the Smallwood algorithm. An amplification factor Q=10 and the 1 Hz sampling rate are applied. Case 3 with consistent mapping result is used as benchmark since we have already shown that higher-order elements with higher DOF have better accuracy. The SRS for cases 1 to 4 is plotted in Figure 16 from 0.1 to 300 Hz which is a typical upper limit for shipboard shock analysis [49]. The L2 error norm of the four cases is presented in Table 9. The SRS result L2 norm difference between consistent and average mapping for each case is presented in Table 10.Figure 16
Shock response spectrum of cases 1 to 4. (a) Case 1. (b) Case 2. (c) Case 3. (d) Case 4.
(a)
(b)
(c)
(d)Table 9
L
2 error norm for cases 1 to 4.
Average mapping
Consistent mapping
Case 1
0.02830
0.01567
Case 2
0.04442
0.03863
Case 3
0.03245
Benchmark
Case 4
0.04577
0.02618Table 10
Difference between consistent and average mapping (L2 norm).
0.1524 m (0.5 ft)
0.3048 m (1 ft)
h = 2
0.01319 (case 1)
0.008317 (case 3)
p
=
2
0.03245 (case 3)
0.06684 (case 4)TheL2 error norm in Table 9 indicates that the consistent mapping algorithm gives a better (lower error norm) result than average mapping in the cases tested. This result makes sense because the displacement/pressure value on quadrature points is more accurate in the consistent mapping case. Average mapping generally overpredicts the peak velocity in the shock response spectrum as shown in Figure 16.If we proceed with the conclusion that consistent mapping gives better results in all the cases, we can treat theL2 norm in Table 10 as the extra mapping error introduced by the average mapping since it is essentially a simplification of the consistent mapping. For the matching mesh configuration (case 2), the difference between consistent and average mapping is the smallest. The difference becomes larger as the fluid element size is increased (case 3 vs case 4). The high-order CASE show a larger difference than the linear element with the same DOF (case 1 vs case 3). Overall, the performance of average mapping becomes poorer for fluid elements with larger size and higher order, which is the case for high-order CASE.


## 3.3.2. Comparison between Conservative and Consistent Schemes

In this section, the behavior of the consistent mapping scheme is compared to the energy-conservative scheme to examine the importance of energy conservation. The difference between the two algorithms is that energy conservation is not guaranteed by the consistent approach as demonstrated in Section2.4. However, energy conservation may not be important for a short-duration problem such as in far-field UNDEX which usually lasts less than 100 ms. To demonstrate the problem, a 3rd-order CASE fluid mesh coupled with the linear structural mesh of the same element size 0.1524 m (0.5 ft) as presented in Figure 3 is used. For the consistent mapping, we use linear Gauss–Legendre quadrature on the structural element which guarantees total force conservation but is theoretically unable to achieve energy conservation. Although consistency is normally violated for the conservative scheme [50], the nested fluid mesh configuration used in this case guarantees consistency by passing the aforementioned patch test as described in Section 1.3.The SRS comparison in Figure17 shows that there is almost no difference in structural response between the consistent and conservative method. We track the total energy input (positive)/received (negative) by fluid and received/input by the structure in the consistent method. The energy history in Figure 18 indicates that the consistent mapping approach essentially conserves the total energy in this application. This supports the assertion made in [41] that an energy conservation scheme might be of little significance in short-duration problems involving shock. Also, given the knowledge that an energy-conservative method can violate consistency if the fluid elements are fully mismatched (not nested), the consistent mapping approach is sufficient and possibly superior for the far-field UNDEX application.Figure 17
Comparison between consistent and energy-conservative approach in SRS.Figure 18
Energy history on fluid and structure interface.


## 4. Conclusions

In this paper, research to apply, assess, and improve the CASE algorithm in an early-time far-field UNDEX problem for application to surface ship problems is described. A partitioned fluid-structure interaction framework is assembled and assessed in two benchmark cases. Results are very good. The computational efficiency of linear finite element (linear CAFE), linear spectral element (linear CASE), and high-order spectral element methods (high-order CASE) using a floating shock platform case is compared. Three information exchange schemes (average, consistent, and energy conservative) for the nonmatching meshes are also compared.The results show that superior computational efficiency is achieved with high-order spectral elements in the surface ship FSP far-field UNDEX problem with the consistent mapping approach even though a large cavitation region exists. The high-order CASE possesses the highest computational efficiency even without the use of tensor-product factorization. As the element order increases, the advantage of high-order CASE becomes larger. Although not as effective as high-order CASE, linear CASE still outperforms the fundamental linear CAFE. This result is due to the more accurate modeling of the cavitation region which significantly affects the accuracy of the structural response.The numerical tests of the mapping schemes reveal structural response overestimation by the average approach where the error becomes larger as the coupled mesh nonmatching increases, especially when the fluid element is larger than the structural element. This indicates that the average mapping approach is not suitable for high-order CASE which often have larger element size than the structural elements.We also found that the energy-conservative method is not necessary for this type of problem because the consistent approach can sufficiently ensure energy conservation and generate almost the same result with easier implementation and less computational effort.Because the floating shock platform problem may not represent the full properties of a full ship, our future research and next paper will apply this method to a medium surface combatant ship.



---

*Source: 7463134-2019-07-25.xml*