Research Article  Open Access
Application of the Spectral Element Method in a Surface Ship FarField UNDEX Problem
Abstract
The prediction of surface ship response to a farfield underwater explosion (UNDEX) requires the simulation of shock wave propagation in the fluid, cavitation, fluidstructure 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 shockrelated problems. Thus, the advantage of CASE remains unproven. In this paper, we present a 3Dpartitioned FSI framework and investigate the application of CAFE and CASE to a surface ship earlytime farfield 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.
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 nearfield and farfield. In nearfield 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 earlytime farfield 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 highfrequency and lowinertia 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 earlytime response of a surface ship to the farfield UNDEX.
1.2. Modeling Techniques for the Fluid Domain in UNDEX
The fluidstructure 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 fluidstructure interaction between an exponentially decaying plane shock wave and an infinite airbacked 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 farfield 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 firstorder DAA includes an earlytime approximation (ETA) and latetime approximation (LTA). The name doubly asymptotic suggests that the approximation approaches the exact solution at both earlytime/highfrequency and latetime/lowfrequency 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 farfield 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 wetsurface approximation (WSA) models have been introduced (e.g., [12]). The investigation of WSA coupled with elastic massspring 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 farfield UNDEX case, finite elementbased approaches are typically favored due to their advantages in solving acoustic wave propagation. The foundation of the FEMtype treatment was proposed by Newton [14] who used a scalar displacement potentialbased 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 cavitationfluid 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 LSDYNA and applied in several farfield 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 pressurebased formulation [24]. Despite the benefits of using the FEMbased fluid models (e.g., CAFE), the dispersive nature of the FEM elements [25] makes them computationally expensive since highlevel 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. Highorder 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 highorder 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 massspring 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 farfield 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 farfield 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 (deepwater 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 highorder SEM possesses a computational advantage over other types of FEM for the surface ship UNDEX problem.
1.3. Information Exchange Schemes at the FluidStructure Interface
In early CAFE applications, nodetonode matching of the linear FEMbased 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 onetoone 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 highorder SEMbased 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 rigidmotion 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 forcemotion 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 residualbased 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 forcemotion 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 forcemotion transfer (DFMT) schemes for nonmatching meshes, the consistent interpolation approach [37–39] and energyconservative 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 shortduration shockwave 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, highpressure 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 farfield surface ship UNDEX problems. In addition, consistent, average, and energyconservative mapping schemes are explored and compared to determine which is the most costeffective 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 farfield UNDEX problem. On the structural side, we use the Abaqus Explicit solver which features explicittime 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 stateoftheart commercial solvers and our inhouse code. Abaqus is coupled through its cosimulation engine, and the inhouse code is coupled through MpCCI API functions. The details regarding the resolution of nonmatching meshes and timestep size are introduced later in this section. The skeleton of the coupled framework is shown in Figure 1.
2.2. Abaqus (Explicit) Structural Solver
Since an important effect of the farfield 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 highspeed events. In comparison, the implicit solver (Abaqus Standard) is more suitable for the lowspeed phenomena like structure dynamic motion after the shock wave effect is gone. The semidiscretized structural motion is expressed in
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 farfield 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 are small, thus allowing the continuity equation based on dynamic displacement :to be linearized using total condensation for :
Neglecting the viscosity and convection term, the Euler momentum equation for the acoustic fluid is written aswhere is the total fluid pressure and 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 based on irrotational fluid assumption, the governing equations, equations (3) and (4), are converted to the form as follows:where is the dynamic condensation and is the dynamic displacement potential related to the dynamic pressure by . Equation (6) enforces a cutoff cavitation condition which stops the total pressure going below vapor pressure once cavitation begins. The vapor pressure is small enough to be negligible. is the static field pressure.
Following the standard Galerkin approach, we multiply equation (5) by test function and integrate it over the fluid domain :
Applying Green’s first identity to equation (7) yields the weak form:where denotes the fluid boundary and 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 aswhere is the column vector of the 1D spectral basis function . 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 rowsum technique.
Substituting equations (9) and (10) into the weak form (equation (8)) giveswhere , , and . Following the nomenclature in [16], we name the global capacitance matrix, the global reactance matrix, and the global fluid force vector. The governing equation (equation (6)) is evaluated in a nodebynode fashion:
The boundary conditions to be prescribed are at the free surface, fluidstructure interface, and nonreflecting boundary (NRB). At the free surface, the essential boundary condition is imposed as
The fluidstructure interaction and NRB terms are natural boundary conditions governed by the fluid force vector b which can be expressed aswherewhere is the displacement vector of the structure along the normal direction of fluidstructure interface taken as positive into the fluid. By the same principle,where is the displacement of NRB approximated by curved wave approximation [29]:where , is the displacement at equilibrium, is the value of the angle between the incident shock direction and the outward norm . 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 (), the last term is dropped.
2.4. Nonconformal FluidStructure 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.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:where denotes the nodes on the structure element and on the fluid element, is the pressure and is the nodal displacement, and are the node number on structure and fluid element, and is the shape function. In Figure 2, we use a 3^{rd} order CASE coupled to the linear structural element as an illustration of this method.
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:
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. EnergyConservative 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:
Then, the received nodal force on the structure element becomesin which is the local coordinate of the mapped fluid nodes on the structure element and is the nodal force on the fluid element which can be represented as follows:where is the fluid surface.
Based on the unit summation property of the shape function (i.e., ), the total received force by the structure element becomes
Thus, the total force is conservative as well. In Figure 3, we use a 3rdorder CASE coupled to the linear structural element as an example for the energyconservative mapping algorithm derivation.
2.5. Temporal Discretization
The governing equations of the CAFE/CASE framework are integrated in time by using a staggered secondorder explicit central difference scheme (a.k.a. leapfrog scheme) introduced by Felippa and Deruntz [16]. The approximation of critical time increment is given aswhere 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 () is used in the cases hereafter.
The nonmatching timestep interpolation is performed by a monotone cubic Hermite spline implemented in MpCCI to achieve a 3^{rd}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.
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 Section 2 are assessed using two benchmark problems. The linear CASE is used in the inhouse 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 is the peak pressure on the shock front, is the speed of sound, is the decay constant, and is the initial location of the shock front at . 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:
Although the problem was originally in 1D, we use a 3D configuration to validate our inhouse framework (Figure 5). 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 inhouse solver and benchmark solution with an L_{2} error norm of 0.0322. The L_{2} error norm used in this case is defined aswhere is the benchmark solution time history, is the computed solution time history, and is the total simulation time.

3.1.2. Floating Shock Platform
In this section, we further assess model and FSI framework in a MILS901D 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 aswhere is the distance from the charge center to the point and denotes the initial wave radius at the nearest structure point. The peak pressure 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 structuralacoustic solver which couples a FEMbased acoustic fluid solver with the Abaqus structural solver. In Abaqus, we employ a CAFElike 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 pressurebased instead of using displacement potential. The nonmatching mesh interpolation is performed by the nodetosurface constraint function in Abaqus which uses the same principle with the energyconservative mapping algorithm as described in Section 2.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.

The 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 inhouse CASE solver agrees well with the Abaqus acousticstructure solver solution as shown in Figure 9. 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 is obtained by combining the magnitude error factor and phase error factor together:
(a)
(b)
The expression for and can be found in [47]. Using the Abaqus results as the baseline, the Russell error norm of the inhouse 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 lowfrequency largeamplitude 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 highorder spectral elements (highorder CASE) using our inhouse solver in the same FSP case configuration used in Section 3.1.2. The result at point 1646 is used for plotting and analysis. We add a damping factor 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 refinement (highorder 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 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.

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 matrix , the matrixvector product () is typically conducted on the element level () and then assembled. The operation count involved in directly multiplying matrix and vector (directmatrix product) iswhere 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 highorder CASE element, the operation is very inefficient. Following the algorithm detailed in [29], we employ tensorproduct factorization (TF) for the matrixvector multiplication () which features operation (FLOP shown in the following equation)
We plot the FLOP of an element with CASE refinement and CAFE/CASE h refinement in Figure 10 using the directmatrix product and tensorproduct factorization. As shown in Figure 10, using tensorproduct factorization significantly lowers the FLOP of highorder CASE but does not help with linear CAFE/CASE. Thus, we decide to use directmatrix product for linear CAFE/CASE and tensorproduct factorization for highorder CASE.
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 byfor linear CAFE/CASE andfor highorder CASE where denotes the total element number and 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 highorder 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 highorder CASE case. Since the time increment of highorder 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.
As shown in Figure 12, the highorder CASE with tensorproduct factorization achieves the highest accuracy compared to the benchmark with the same number of FLOP. The computational efficiency of highorder CASE is superior even without tensorproduct factorization (directmatrix product). Thus, the highorder CASE may still be advantageous if an unstructured tetrahedral mesh without tensorproduct factorization is used. Linear CASE outperforms linear CAFE although inferior to the highorder CASE.
Next, the behavior of highorder 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 8^{th}order CASEtype refinement and 2^{nd}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.


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 error norm definition inwhere is the depth and is the coordinate in the depth direction. Again, the result of the CASE solution is used as a benchmark.
The norm for the h and refinement levels from 4 to 6 are tabulated (Tables 6–8), but only the 4^{th}level h/prefined linear CASE, linear CAFE, and highorder 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.



(a)
(b)
(c)
The error norm results show that highorder 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, highorder 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 Section 3.2, it is shown that the computational advantage of highorder 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 Figure 15. 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 2^{nd}order h refinement; the second case uses a 0.3048 m (1 ft) base element with 2^{nd}order h refinement (matching fluid and structure mesh); the third case uses a 0.1524 m (0.5 ft) base element with 2^{nd}order p refinement; the fourth case uses a 0.3048 m (1 ft) base element with 2^{nd}order refinement. Cases 1 and 3, and cases 2 and 4 have the same degree of freedom.
(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 fluidstructure 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 and the 1 Hz sampling rate are applied. Case 3 with consistent mapping result is used as benchmark since we have already shown that higherorder 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 L_{2} error norm of the four cases is presented in Table 9. The SRS result L_{2} norm difference between consistent and average mapping for each case is presented in Table 10.
(a)
(b)
(c)
(d)


The L_{2} 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 the L_{2} 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 highorder 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 highorder CASE.
3.3.2. Comparison between Conservative and Consistent Schemes
In this section, the behavior of the consistent mapping scheme is compared to the energyconservative 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 Section 2.4. However, energy conservation may not be important for a shortduration problem such as in farfield UNDEX which usually lasts less than 100 ms. To demonstrate the problem, a 3^{rd}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 Figure 17 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 shortduration problems involving shock. Also, given the knowledge that an energyconservative method can violate consistency if the fluid elements are fully mismatched (not nested), the consistent mapping approach is sufficient and possibly superior for the farfield UNDEX application.
4. Conclusions
In this paper, research to apply, assess, and improve the CASE algorithm in an earlytime farfield UNDEX problem for application to surface ship problems is described. A partitioned fluidstructure 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 highorder spectral element methods (highorder 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 highorder spectral elements in the surface ship FSP farfield UNDEX problem with the consistent mapping approach even though a large cavitation region exists. The highorder CASE possesses the highest computational efficiency even without the use of tensorproduct factorization. As the element order increases, the advantage of highorder CASE becomes larger. Although not as effective as highorder 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 highorder CASE which often have larger element size than the structural elements.
We also found that the energyconservative 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.
Data Availability
The data used to support the findings of this study are included in the repository: https://github.com/lzhaok6/Research_data_Application_of_the_Spectral_Element_Method_in_a_Surface_Ship_Farfield_UNDEX_Problem.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the Office of Naval Research for the project on “Naval Ship Distributed System Vulnerability and Battle Damage Recovery in EarlyStage Ship Design” (Grant no. N000141512476).
References
 A. H. Keil, The Response of Ships to Underwater Explosions, Transactions of Society of Naval Architects and Marine Engineers, New York, NY, USA, 1961.
 W. D. Reid, The Response of Surface Ships to Underwater Explosions, DSTO Aeronautical and Maritime Research Laboratory, Melbourne, Victoria, Australia, 1996.
 K. Mäkinen, “Cavitation models for structures excited by a plane shock wave,” Journal of Fluids and Structures, vol. 12, no. 1, pp. 85–101, 1998. View at: Publisher Site  Google Scholar
 G. I. Taylor, “The pressure and impulse of submarine explosion waves on plates,” in The Scientific Papers of Sir Geoffrey Ingram Taylor, pp. 287–303, Cambridge University Press, Cambridge, UK, 1941. View at: Google Scholar
 A. K. Mathew, Modeling Underwater Explosion (UNDEX) Shock Effects for Vulnerability Assessment in Early Stage Ship Design, Virginia Tech, Blacksburg, VA, USA, 2018.
 T. L. Geers, “Doubly asymptotic approximations for vibration analysis of submerged structures,” Journal of the Acoustical Society of America, vol. 64, no. 5, pp. 1500–1508, 1978. View at: Publisher Site  Google Scholar
 T. L. Geers, “Residual potential and approximate methods for threedimensional fluidstructure interaction problems,” Journal of the Acoustical Society of America, vol. 49, no. 5B, pp. 1505–1510, 1970. View at: Publisher Site  Google Scholar
 T. L. Geers and B. J. Toothaker, “Thirdorder doubly asymptotic approximations for computational acoustics,” Journal of Computational Acoustics, vol. 8, no. 1, pp. 101–120, 2000. View at: Publisher Site  Google Scholar
 J. A. DeRuntz, T. L. Geers, and C. A. Felippa, The Underwater Shock Analysis (USA) Code: A Reference Manual, Defense Technical Information Center, Fort Belvoir, VA, USA, 1978.
 P. K. Fox, Y. W. Kwon, and Y. S. Shin, Nonlinear Response of Cylindrical Shells to Underwater Explosion: Testings and Numerical Prediction Using USA/DYNA3D, Naval Postgraduate School, Monterey, CA, USA, 1992.
 S. L. Wood and Y. S. Shin, Cavitation Effects on a Shiplike Box Structure Subjected to an Underwater Explosion, Naval Postgraduate School, Monterey, CA, USA, 1998.
 F. L. DiMaggio, I. S. Sandler, and D. Rubin, “Uncoupling approximations in fluidstructure interaction problems with cavitation,” Journal of Applied Mechanics, vol. 48, no. 4, pp. 753–756, 1981. View at: Publisher Site  Google Scholar
 M. A. Sprague and T. L. Geers, “Computational treatments of cavitation effects in nearfreesurface underwater shock analysis,” Shock and Vibration, vol. 8, no. 2, pp. 105–122, 2001. View at: Publisher Site  Google Scholar
 R. E. Newton, “Effects of cavitation on underwater shock loadingplane problem part 1,” Naval Postgraduate School, Monterey, CA, USA, 1979, NPS 6979007PR. View at: Google Scholar
 O. C. Zienkiewicz, D. K. Paul, and E. Hinton, “Cavitation in fluidstructure response (with particular reference to dams under earthquake loading),” Earthquake Engineering & Structural Dynamics, vol. 11, no. 4, pp. 463–481, 1983. View at: Publisher Site  Google Scholar
 C. A. Felippa and J. A. Deruntz, “Finite element analysis of shockinduced hull cavitation,” Computer Methods in Applied Mechanics and Engineering, vol. 44, no. 3, pp. 297–337, 1984. View at: Publisher Site  Google Scholar
 Y. S. Shin, “Ship shock modeling and simulation for farfield underwater explosion,” Computers & Structures, vol. 82, no. 23–26, pp. 2211–2219, 2004. View at: Publisher Site  Google Scholar
 S. W. Gong and K. Y. Lam, “On attenuation of floating structure response to underwater shock,” International Journal of Impact Engineering, vol. 32, no. 11, pp. 1857–1877, 2006. View at: Publisher Site  Google Scholar
 D. B. Woyak, “Modeling submerged structures loaded by underwater explosions with ABAQUS/explicit,” in Proceedings of the ABAQUS Users Conference, pp. 1–15, Newport, RI, USA, June 2002. View at: Google Scholar
 X. Li, X. Wei, P. Zheng, and W. G. Wu, “The response of a floating structure due to underwater explosion with cavitation effect,” International Society of Offshore and Polar Engineers, vol. 7, pp. 633–638, 2010. View at: Google Scholar
 Z. Aman, Z. Weixing, W. Shiping, and F. Linhan, “Dynamic response of the noncontact underwater explosions on naval equipment,” Marine Structures, vol. 24, no. 4, pp. 396–411, 2011. View at: Publisher Site  Google Scholar
 Z. Zong, Y. Zhao, and H. Li, “A numerical study of whole ship structural damage resulting from closein underwater explosion shock,” Marine Structures, vol. 31, pp. 24–43, 2013. View at: Publisher Site  Google Scholar
 G. Wang, S. Zhang, M. Yu, H. Li, and Y. Kong, “Investigation of the shock wave propagation characteristics and cavitation effects of underwater explosion near boundaries,” Applied Ocean Research, vol. 46, pp. 40–53, 2014. View at: Publisher Site  Google Scholar
 Dassault Systèmes Simulia, Abaqus 6.14 Theory Guide, Dassault Systèmes Simulia, Johnston, RI, USA, 2014.
 J. D. de Basabe and M. K. Sen, “Grid dispersion and stability criteria of some common finiteelement methods for acoustic and elastic wave equations,” Geophysics, vol. 72, no. 6, pp. T81–T95, 2007. View at: Publisher Site  Google Scholar
 M. A. Sprague and T. L. Geers, “Spectral elements and field separation for an acoustic fluid subject to cavitation,” Journal of Computational Physics, vol. 184, no. 1, pp. 149–162, 2003. View at: Publisher Site  Google Scholar
 M. A. Sprague and T. L. Geers, “A spectralelement method for modelling cavitation in transient fluidstructure interaction,” International Journal for Numerical Methods in Engineering, vol. 60, no. 15, pp. 2467–2499, 2004. View at: Publisher Site  Google Scholar
 M. A. Sprague and T. L. Geers, “A spectralelement/finiteelement analysis of a shiplike structure subjected to an underwater explosion,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 1718, pp. 2149–2167, 2006. View at: Publisher Site  Google Scholar
 M. A. Sprague, Advanced Computational Techniques for the Analysis of 3D FluidStructure Interaction with Cavitation, University of Colorado at Boulder, Boulder, CO, USA, 2002.
 A. T. Patera, “A spectral element method for fluid dynamics: laminar flow in a channel expansion,” Journal of Computational Physics, vol. 54, no. 3, pp. 468–488, 1984. View at: Publisher Site  Google Scholar
 B. Klenow, Finite and Spectral Element Methods for Modeling FarField Underwater Explosion Effects on Ships, Virginia Tech, Blacksburg, VA, USA, 2009.
 A.M. Zhang, S.F. Ren, Q. Li, and J. Li, “3D numerical simulation on fluidstructure interaction of structure subjected to underwater explosion with cavitation,” Applied Mathematics and Mechanics, vol. 33, no. 9, pp. 1191–1206, 2012. View at: Publisher Site  Google Scholar
 Y. Shin and N. Schneider, “Ship shock trial simulation of USS Winston S. Churchill (DDG 81): modeling and simulation strategy and surrounding fluid volume effects,” in Proceedings of the 74th Shock and Vibration Symposium, San Diego, CA, USA, October 2003. View at: Google Scholar
 C. A. Felippa, K. C. Park, and M. R. Ross, “A classification of interface treatments for FSI,” in Fluid Structure Interaction II: Modelling, Simulation, Optimization, H.J. Bungartz, M. Mehl, and M. Schäfer, Eds., pp. 27–51, Springer, Berlin, Germany, 2010. View at: Google Scholar
 C. Bernardi, Y. Maday, and A. T. Patera, “Domain decomposition by the mortar element method,” in Asymptotic And Numerical Methods For Partial Differential Equations With Critical Parameters, H. G. Kaper, M. Garbey, and G. W. Pieper, Eds., pp. 269–286, Springer, Dordrecht, Netherlands, 1993. View at: Google Scholar
 M. R. Ross, C. A. Felippa, K. C. Park, and M. A. Sprague, “Treatment of acoustic fluidstructure interaction by localized lagrange multipliers: formulation,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 33–40, pp. 3057–3079, 2008. View at: Publisher Site  Google Scholar
 C. Farhat, M. Lesoinne, and N. Maman, “Mixed explicit/implicit time integration of coupled aeroelastic problems: threefield formulation, geometric conservation and distributed solution,” International Journal for Numerical Methods in Fluids, vol. 21, no. 10, pp. 807–835, 1995. View at: Publisher Site  Google Scholar
 C. Farhat, M. Lesoinne, and P. Le Tallec, “Load and motion transfer algorithms for fluid/structure interaction problems with nonmatching discrete interfaces: momentum and energy conservation, optimal discretization and application to aeroelasticity,” Computer Methods in Applied Mechanics and Engineering, vol. 157, no. 12, pp. 95–114, 1998. View at: Publisher Site  Google Scholar
 R. Lohner, C. Yang, J. Cerbal et al., “Fluidstructurethermal interaction using a loose coupling algorithm and adaptive unstructured grids,” in Proceedings of the 29th American Institute of Aeronautics and Astronautics, Fluid Dynamics Conference, Albuquerque, NM, USA, June 1998. View at: Publisher Site  Google Scholar
 R. K. Jaiman, X. Jiao, P. H. Geubelle, and E. Loth, “Assessment of conservative load transfer for fluid–solid interface with nonmatching meshes,” International Journal for Numerical Methods in Engineering, vol. 64, no. 15, pp. 2014–2038, 2005. View at: Publisher Site  Google Scholar
 M. R. Ross, M. A. Sprague, C. A. Felippa, and K. C. Park, “Treatment of acoustic fluidstructure interaction by localized lagrange multipliers and comparison to alternative interfacecoupling methods,” Computer Methods in Applied Mechanics and Engineering, vol. 198, no. 9–12, pp. 986–1005, 2009. View at: Publisher Site  Google Scholar
 MpCCI 4.5 User Manual, Fraunhofer Institute for Algorithms and Scientific Computing SCAI, Sankt Augustin, Germany, 2017.
 H. H. Bleich and I. S. Sandler, “Interaction between structures and bilinear fluids,” International Journal of Solids and Structures, vol. 6, no. 5, pp. 617–639, 1970. View at: Publisher Site  Google Scholar
 MILS901D, Military Specification: Shock Tests. H.I. (HighImpact) Shipboard Machinery, Equipment, and Systems, Requirements for (17Mar1989) [S/S BY MILDTL901E], Department of the Navy, Washington, DC, USA, 1989.
 R. H. Cole, Underwater Explosions, Princeton University Press, Princeton, NJ, USA, 1948.
 M. M. Swisdak, Explosion Effects and Properties: Part II—Explosion Effects in Water, Naval Surface Warfare Center, Washington, DC, USA, 1978.
 D. M. Russell, “Error measures for comparing transient data: part I: development of a comprehensive error measure Part II: error measures case study,” in Proceedings of the 68th Shock and Vibration Symposium, pp. 3–6, Baltimore, MD, USA, November 1997. View at: Google Scholar
 T. Irvine, An Introduction to the Shock Response Spectrum, Vibrationdata, Madison, AL, USA, 2012.
 E. W. Clements, Shipboard Shock and Navy Devices for its Simulation, Defense Technical Information Center, Fort Belvoir, VA, USA, 1972.
 A. de Boer, A. H. van Zuijlen, and H. Bijl, “Comparison of conservative and consistent approaches for the coupling of nonmatching meshes,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 4950, pp. 4284–4297, 2008. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Zhaokuan Lu and Alan Brown. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.