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

Lu, Zhaokuan; Brown, Alan

doi:https://doi.org/10.1155/2019/7463134

Shock and Vibration

On this page

Abstract Introduction Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2019 | Article ID 7463134 | https://doi.org/10.1155/2019/7463134

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

Zhaokuan Lu¹and Alan Brown¹

Academic Editor: Jean-Mathieu Mencik

Received06 Apr 2019

Accepted11 Jul 2019

Published25 Jul 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.

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.

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.

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

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 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 row-sum 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 node-by-node fashion:

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

The fluid-structure 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 fluid-structure 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 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.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. 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:

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 3rd-order CASE coupled to the linear structural element as an example for the energy-conservative mapping algorithm derivation.

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 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 time-step 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 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 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 in-house 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 in-house solver and benchmark solution with an L₂ error norm of 0.0322. The L₂ 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 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 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 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 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 in-house CASE solver agrees well with the Abaqus acoustic-structure 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 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 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 (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 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 matrix-vector product () is typically conducted on the element level () and then assembled. The operation count involved in directly multiplying matrix and vector (direct-matrix 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 high-order CASE element, the operation is very inefficient. Following the algorithm detailed in [29], we employ tensor-product factorization (TF) for the matrix-vector 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 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.

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 high-order 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 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.

As shown in Figure 12, 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.

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 8^th-order CASE-type 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/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.

(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 Section 3.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 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 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 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 L₂ error norm of the four cases is presented in Table 9. The SRS result L₂ norm difference between consistent and average mapping for each case is presented in Table 10.

(a)

(b)

(c)

(d)

The L₂ 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₂ 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 Section 2.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 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 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.

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.

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_Far-field_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 Early-Stage Ship Design” (Grant no. N00014-15-1-2476).

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 three-dimensional fluid-structure 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, “Third-order 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 Ship-like 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 fluid-structure 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 near-free-surface 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 loading-plane problem part 1,” Naval Postgraduate School, Monterey, CA, USA, 1979, NPS 69-79-007PR.
View at: Google Scholar
O. C. Zienkiewicz, D. K. Paul, and E. Hinton, “Cavitation in fluid-structure 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 shock-induced 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 far-field 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 non-contact 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 close-in 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 finite-element 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 spectral-element method for modelling cavitation in transient fluid-structure 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 spectral-element/finite-element analysis of a ship-like structure subjected to an underwater explosion,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 17-18, pp. 2149–2167, 2006.
View at: Publisher Site | Google Scholar
M. A. Sprague, Advanced Computational Techniques for the Analysis of 3-D Fluid-Structure 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 Far-Field 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 fluid-structure 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 fluid-structure 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: three-field 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 non-matching discrete interfaces: momentum and energy conservation, optimal discretization and application to aeroelasticity,” Computer Methods in Applied Mechanics and Engineering, vol. 157, no. 1-2, pp. 95–114, 1998.
View at: Publisher Site | Google Scholar
R. Lohner, C. Yang, J. Cerbal et al., “Fluid-structure-thermal 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 non-matching 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 fluid-structure interaction by localized lagrange multipliers and comparison to alternative interface-coupling 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
MIL-S-901D, Military Specification: Shock Tests. H.I. (High-Impact) Shipboard Machinery, Equipment, and Systems, Requirements for (17-Mar-1989) [S/S BY MIL-DTL-901E], 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 non-matching meshes,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 49-50, 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.

PDF Download Citation

Download other formats

Order printed copies

Views

909

Downloads

1066

Citations