## Abstract

The hydroelastic interaction between an underwater explosion and an elastic plate is investigated num- erically through a domain-decomposition strategy. The three-dimensional features of the problem require a large computational effort, which is reduced through a weak coupling between a one-dimensional radial blast solver, which resolves the blast evolution far from the boundaries, and a three-dimensional compressible flow solver used where the interactions between the compression wave and the boundaries take place and the flow becomes three-dimensional. The three-dimensional flow solver at the boundaries is directly coupled with a modal structural solver that models the response of the solid boundaries like elastic plates. This enables one to simulate the fluid–structure interaction as a strong coupling, in order to capture hydroelastic effects. The method has been applied to the experimental case of Hung *et al.* (2005 *Int. J. Impact Eng.* **31**, 151–168 (doi:10.1016/j.ijimpeng.2003.10.039)) with explosion and structure sufficiently far from other boundaries and successfully validated in terms of the evolution of the acceleration induced on the plate. It was also used to investigate the interaction of an underwater explosion with the bottom of a close-by ship modelled as an orthotropic plate. In the application, the acoustic phase of the fluid–structure interaction is examined, highlighting the need of the fluid–structure coupling to capture correctly the possible inception of cavitation.

## 1. Introduction

Explosions in fluids represent an important topic in many research fields. In medicine, for example, implosions of micro-bubbles can be caused by the ultrasound used to remove kidney stones. For marine applications, two relevant dangerous examples are underwater explosions, which can damage close-by vehicles, and cavitating bubbles on propeller blades possibly leading to erosion. Either when we are talking of micro- or macro-explosions, the common features are: (i) generation of an air bubble, (ii) propagation of a compression wave from the surface of the bubble, and (iii) eventual interaction of the compression wave with the surrounding structures.

Here, this phenomenon is investigated in the framework of underwater explosions, which represents one of the most important causes of loss of navy vessels and it is a danger for close-by submarines, surface ships and offshore structures. This work continues the research proposed in [1]. There, a numerical solution strategy is proposed and described in detail. It is based on a one-dimensional–three-dimensional fluid-dynamic domain-decomposition (DD) coupling. A radial solver is used to study the initial stages of the phenomenon, before the acoustic waves reach a structure, then, locally, a three-dimensional method handles the fluid–structure interaction. The three-dimensional sub-domain is limited by the body boundaries and the one-dimensional–three-dimensional interface across which information is transmitted. In [1], the structure was modelled as rigid and preliminary studies of the fluid–structure interactions accounting for elastic structural behaviour and an approximated coupling with the fluid were performed using only the radial solver over the whole fluid domain and for the whole evolution. On the contrary, in the present contribution the fluid-dynamic problem is described with the same hybrid method but the fluid–elastic structure coupling is fully simulated. This enables one to investigate the hydroelastic behaviour, which destroys any radial symmetry in the flow. Some preliminary results have been documented in [2]. The paper is structured as follows. The proposed method is illustrated in §2, with emphasis on the new aspects. Two selected important cases are analysed in §3. Finally, some conclusions and lines for future research are proposed in §4.

## 2. The numerical solver

Underwater explosions are characterized by several complex phenomena: (i) multi-phase flow, involving at least water and gas cavities resulting from the blast, (ii) flow compression, involving the propagation of acoustic waves during cavity expansion, and (iii) transition from compressible to incompressible flow dynamics during the later cavity oscillations.

The developed solution strategy focuses on underwater explosions and their interaction with structures far from other boundaries, such as the sea floor or the sea surface. The method introduced in [1,3] combines two solvers within a zonal approach to capture the complex and detailed flow features described above. Here, the main aspects of the solver are described, while more details can be found in [1], where, also accuracy and convergence properties of the two involved fluid-dynamic solvers and the structural modelling are assessed. Moreover, the scalability of the parallelized part of the method, i.e. the relationship between required CPU time and number of processors, is provided. Thus, these aspects are not discussed here.

### (a) The compressible solver for two-phase flows: three-dimensional and radial solutions

We examine the flow caused by an underwater explosion assuming that the chemical reactions and the detonation phase have ended and have generated a gas cavity with known initial physical properties. This flow involves at least two phases, gas and water, and it is, in general, three-dimensional and compressible. We neglect viscous effects and solve the problem in a Cartesian Earth-fixed coordinate system (*x*,*y*,*z*). We retain flow nonlinearities since nonlinear features are in general important to properly capture the initiation and propagation of the acoustic waves which, at the first stages, are always rather strong and that can still be important also during a wave–body interaction if the explosion occurs sufficiently close to the wall and/or if the explosion is sufficiently strong. Moreover, nonlinear flow features could be induced by strong hydrodynamic–structural coupling. Under these assumptions, the governing equation can be formally written as
2.1

Here, ** V**=(

*u*,

*v*,

*w*) is the velocity,

*p*the pressure,

*ρ*the density,

*E*=

*ρ*[

*e*+(

*u*

^{2}+

*v*

^{2}+

*w*

^{2})/2] the total energy,

**=[**

*U**u*,

*ρu*,

*ρv*,

*ρw*,

*E*]

^{T}the vector of fluid variables and

**=[**

*F**F*

_{x},

*F*

_{y},

*F*

_{z}]

^{T}the vector of fluxes with

*F*

_{x}=[

*ρu*,

*ρu*

^{2}+

*p*,

*ρuv*,

*ρuw*,(

*E*+

*p*)

*u*]

^{T},

*F*

_{y}=[

*ρv*,

*ρuv*,

*ρv*

^{2}+

*p*,

*ρvw*,(

*E*+

*p*)

*v*]

^{T}and

*F*

_{z}=[

*ρw*,

*ρuw*,

*ρvw*,

*ρw*

^{2}+

*p*,(

*E*+

*p*)

*w*]

^{T}. To complete the problem, an equation of state (EOS) of the form 2.2 is assumed for the specific internal energy

*e*. The functions

*f*

_{f}and

*g*

_{f}depend on the fluid properties, as indicated by the subscript f. In particular, for the gas generated by the explosion, they are set in agreement with the Jones–Wilkins–Lee EOS [4] and for the water they follow an isentropic Tait relation [5]. Equations (2.1) and (2.2) must be complemented by the initial and boundary conditions which apply in the examined case.

Within the three-dimensional framework, equation (2.1) is solved numerically with a second-order finite-difference scheme in space and integrated in time with a total variation diminishing third-order Runge–Kutta scheme (see [6]), with the fluxes in each fluid modelled with a Harten–Lax–van Leer contact Riemann solver (see [7]).

A two-shock approximation to the Riemann problem is enforced at the gas–cavity interface as proposed in [8] solving for the two nonlinear characteristics intersecting at the interface and using mass and momentum jump conditions for the transmitted and reflected shocks. The resulting equation system is nonlinear and it is solved iteratively with a Newton–Raphson method providing values of (** V**⋅

**)**

*n*_{i}=

*u*

_{i},

*p*

_{i},

*ρ*

^{L}

_{i}and

*ρ*

^{R}

_{i}, respectively, the normal velocity and pressure at the interface and the left and right density. Here, ‘left’ means ‘inside the gas cavity’ and ‘right’ means ‘in water’. The estimated values for

*ρ*

^{L}

_{i}and are corrected by enforcing an isobaric condition across the interface so to avoid numerical errors in the pressure. This interface algorithm is based on a ghost fluid method (see [8]) and provides the conditions across the interface to each fluid. Once all interface variables are known, the fluxes

**can be calculated in each fluid and the problem can be stepped forward in time. To follow the evolution of the gas–water interface, the use of the three-dimensional Eulerian solver requires an additional technique. In the present case, a colour-function method is adopted, advecting in time the signed distance (positive in one fluid and negative in the other) from the interface. This is known as level-set function**

*F**ϕ*and its zero value identifies the interface, which can be reconstructed in time by interpolating the

*ϕ*values at the centre of the cells.

*ϕ*is advected in time using the equation 2.3 where

*V*_{i}is the interface velocity calculated like in [8]; in particular the normal velocity comes from the interface solution strategy described above, while the two tangential components are obtained as averages between the solutions in the two fluids. Once the instantaneous interface position is identified as

*ϕ*=0, the Riemann problem described above is solved across the interface to ensure proper reflection and transmission of shock waves.

To make the solution efficient in time, an adaptive mesh refinement is used according to [9]. The grid is subdivided in blocks and each block can be dynamically (i.e. during the simulation) split into eight smaller blocks in case a finer mesh is locally required, or neighbouring blocks can be merged into one single, larger block when the flow variations become sufficiently small. This splitting/merging process is controlled by a criterion which monitors the intensity of the fluid variable ∇*p*. Within such an adaptive algorithm, an interpolation strategy can be required; here, a second-order polynomial expression is chosen consistently with the order of accuracy of the three-dimensional numerical scheme.

The described solver can incorporate the presence of solid and deformable boundaries like, respectively, the sea floor and the sea surface. For both of them, a level-set technique can be applied when the grid is not adapted to the boundary, as may occur for the evolving air–water interface. One must, however, note that special care must be taken when the acoustic wave resulting from an explosion reaches the sea surface, because the interaction involves, among others, partial reflection and cavitation phenomena (e.g. [10]).

Within the radial solver, the governing equation (2.1) of the problem reduces to a one-dimensional Euler equation in the radial direction, *r*, of the form
2.4
with ** U**=[

*ρ*,

*ρu*,

*E*]

^{T},

**=[**

*F**ρu*,

*ρu*

^{2}+

*p*,(

*E*+

*p*)

*u*]

^{T}and

**=2[**

*S**ρu*/

*r*,

*ρu*

^{2}/

*r*,

*u*(

*E*+

*p*)/

*r*]

^{T}. Here,

*u*is the radial velocity,

*p*the pressure,

*ρ*the density and

*E*the total energy

*ρ*(

*e*+

*u*

^{2}/2). The same EOSs of the three-dimensional case are used for the specific internal energy

*e*to close the problem for gas and water. The numerical implementation of the one-dimensional formulation is based on a first-order, finite-difference scheme in space and time, which proved to be also very accurate when compared with other numerical solutions, for fully one-dimensional problems and problems with radial symmetry, even for long-time evolutions [1]. For this reason, no attempt was made to implement a second-order formulation to improve the solution efficiency. The fluxes and the interface conditions are handled as in three dimensions.

For both three-dimensional and radial formulations of the solver, the time step is dynamically chosen on the basis of a Courant–Friedrichs–Lewy stability condition, which in our case reads Δ*t*≤*k*⋅CFL. CFL is min(Δ*r*/max(|*u*−*c*|,|*u*+*c*|)) for the one-dimensional solution and it is similarly defined for the three-dimensional solution. The safety factor *k* is chosen to avoid diverging behaviours of the solution. Here, Δ*r* is the local spatial discretization, and *c* the local speed of sound, so min() and max(), respectively, mean minimum and maximum values in the computational domain.

If the explosion occurs very far from the solid boundaries, a spherical gas cavity is formed, whose evolution is well approximated by a radial symmetry, as gravity does not play a relevant role and the hydrostatic pressure enters the problem only in terms of the surrounding background pressure for the bubble. This has motivated the implementation of a ‘one-dimensional’ radial solution. When the explosion interacts with a structure, the one-dimensional solver can be still combined with a structural model, as long as hydroelastic effects are negligible. This allows a simplified study of the fluid–structure interactions. In particular, this means that a radial solver can be applied until the acoustic wave reflected from the structure reaches the blast bubble and goes back to the body. This alternative has been implemented and assessed in [1] and is also analysed in the physical investigation documented in §3.

Here, a three-solvers coupling is used. The one-dimensional method can be applied at any time in a sub-domain sufficiently far from the region where three-dimensional effects caused by the fluid–body interaction are important, while the three-dimensional method is used to model the near-body flow. In this way, a composite solver, able to describe the whole stages of the fluid–body interaction, can be built and, at the same time, the computational costs (CPU time) are limited as much as possible.

### (b) One-way time–space domain decomposition for compressible two-phase flows

The two compressible solvers described in §2*a* can be combined within a DD strategy to investigate underwater explosions with initial radial symmetry. Here, this hybrid method is proposed as a one-way, time–space coupling, which means that the radial solver is thought to be used over the whole domain for time *t*≤*t**, where *t** represents the time at which the compression wave approaches the one-dimensional–three-dimensional boundary. At *t*=*t**, the DD is switched on and the three-dimensional solver is introduced in a zone containing the structure where three-dimensional effects are expected to be important, with initial conditions provided by the radial solver. The latter is still used in the rest of the fluid domain and provides also the boundary conditions for the three-dimensional solver at the common interface. The position of the one-dimensional–three-dimensional boundary can be stated at the beginning of the simulation or it can be changed to reduce as much as possible the three-dimensional domain.

The one-dimensional–three-dimensional solver–solver interface is not meant to be a sharp surface but it is conceived as an overlapping region between the spherical sub-domain of the radial solution and the prismatic sub-domain of the three-dimensional solution. Within this overlapping layer, the fluid variables are forced to smoothly vary from the values provided by one solver to those estimated by the other. The use of an overlapping region has been found to give a more robust numerical solution than the use of a sharp interface. On the other hand, inside the overlapping region both solvers are applied, so the thickness of this layer must be as small as possible. Here, this thickness is set 4Δ*x*_{loc}, Δ*x*_{loc} being the local size of the three-dimensional computational grid, and a cosine function is used for the smooth variation of the solution across the overlapping region itself. The DD strategy significantly limits the computational costs with respect to using a compressible three-dimensional solver over the whole domain and for the entire time evolution. This can be easily understood if one considers that the computational cost and memory space for the three-dimensional solver increase with *N*^{3}, *N* being the number of points per spatial direction; they increase instead with *N* for a radial solution.

If the explosion occurs sufficiently close to boundaries, such as the sea surface, the sea floor, a marine structure, so that no radial symmetry can be assumed even near the bubble, then the above-proposed DD strategy cannot be applied and the fully three-dimensional solver must be used from the beginning of the simulation and everywhere in the fluid domain.

### (c) Structural modelling

A marine structure exposed to a blast can experience elastic and plastic deformations even leading to rupture. The initial phase, i.e. the acoustic phase, of the incident flow is important for the local forcing on the structure and is analysed by this study. The aim is to investigate the importance of hydroelastic effects and other involved phenomena. This is done assuming only elastic deformations of the body can take place. It means the use of a simplified structural modelling, which has, nevertheless, proved to be able to capture the inception for plastic deformations, as well as the maximum strain rate and dynamic yield stress at hand [1]. Assuming the structure to be the bottom of a vessel, it is modelled as a rectangular orthotropic plate with length *L* and width *B* in the *x* and *y* directions, respectively. The plate is assumed to undergo a linear deformation *w*(*x*,*y*,*t*) governed in time and space by the equation
2.5

Here, *m* is the average plate mass per unit area, *D*_{x} and *D*_{y} are its flexural rigidities in the two main directions and *BB* is its effective torsional rigidity (e.g. [11]). On the right-hand side, *p* is the local hydrodynamic pressure acting on the plate which depends on space, time and *w*, in particular on its time derivatives. This equation is solved with a modal approach by using separation of variables and expressing the deformation *w* as bi-linear combinations of the eigenfunctions of vibrating beams satisfying the plate boundary conditions in both *x* and *y* directions. Each of these combinations represents a plate mode and it is associated with an unknown amplitude depending on time. The number of eigenfunctions in *x* and *y* is taken to be finite, say *N*_{x} and *N*_{y}, leading to *N*=*N*_{x}⋅*N*_{y} plate modes, then the modal expression of *w* is introduced in equation (2.5) and the resulting equation is projected along each of the *N* plate modes. This leads to a linear ordinary differential (in time) equation system for the *N* unknown amplitudes of the modes.

For simplicity purposes, we have here used a ‘thin-plate’ theory (e.g. Kirchhoff–Love plate theory). However, nothing prevents the use of a more complicated model valid for thick plates (e.g. Mindlin–Reissner plate theory). In fact, the DD strategy is independent of the chosen plate theory. Hence, the mentioned strategy can be applied also in combination with higher order plate theories, as long as the plate model gives both the displacement and the displacement velocity back to the fluid solver and receives as input from it the fluid pressure on the plate. When the DD strategy is used, and the three-dimensional solver is used in the fluid region containing the structure, *p* is obtained by solving the fluid–structure coupled problem, i.e. it is affected by the plate elastic behaviour. In more detail, the structure represents a solid deformable interface immersed in the fluid domain. This interface is identified at any instant by a level-set function, similarly to what is done for the gas–water interface. Inside the body, ‘guard’ cells [9] are introduced which can be considered as a sort of ‘ghost cells’. In those cells, the velocity component *v*_{N}(*P*) normal to the body surface at point *P* is made equal to −*v*_{N}(*P*_{mirror})+2*v*_{N wall}, where *v*_{N wall} is the normal velocity of the wall and *P*_{mirror} is the mirror point of *P* inside the fluid. For the tangential components, it is *v*_{τ}(*P*)=*v*_{τ}(*P*_{mirror}). Then, the fluid-dynamic problem (2.1) can be prolonged in time and provides the pressure. The latter can be extrapolated on the structure, i.e. along the iso-surface with zero value of the level-set function associated with the body and used in the structural problem (2.5) to find *w*. Once this is known, the position of the body can be updated in time through use of the level-set technique so that the new body–fluid boundary condition can be enforced again in the fluid dynamic problem. As we assume elastic and small deformations, the errors are very limited if the body-boundary condition is applied on the mean plate configuration.

If the radial solution is used for simulating the incident acoustic wave, the direct full fluid–structure interaction cannot be examined. The coupling is, then, approximated assuming that the linear superposition principle is valid for the fluid–body interaction. Therefore, the pressure *p* is decomposed as
2.6
where *ρ* and *c* are the density and speed of sound in calm water, *p*_{in} is the incident wave pressure and *θ* is the angle between the direction of the incident wave reaching locally the plate and the plate normal vector. The first two terms on the right-hand side of equation (2.6) are important during the acoustic phase and represent the sum of the incident and reflected (from the wall) wave pressure. This leads to twice *p*_{in} for a fixed wall ‘plus’ an acoustic wave-radiation damping in the case of rigidly moving or deformable plate, consistent with the theory developed in [12] for explosion waves interacting with plates. The last term represents an incompressible added-mass contribution, which is important when the water behaves as an incompressible fluid and if the plate oscillates as a consequence of the interaction with the incident acoustic wave. In the present implementation, the acoustic wave-radiation damping term is switched off when *p*_{in} becomes half of the acoustic pressure associated with the incident wave and, then, the added-mass term is switched on.

The added-mass contribution for each plate mode is estimated in an approximate way, studying the radiation problem for each plate mode and assuming that the plate is surrounded by a fixed wall and then by a flat ‘free surface’ with velocity potential *φ*=0. This means that it is a high-frequency added mass estimate and introduces, in general, some inaccuracies. For example, for a ship the structure of interest could be the bottom of the vessel, eventually exposed to an explosion. In this case, the plate would be used to model the grillage on the bottom and this means that it cannot be at the same level of the free surface. The surrounding ‘free surface’ used in the added mass solution is meant to avoid a purely Neumann-type problem, which would lead to a solution in general dependent on an arbitrary constant. The dimensions of the fixed wall are assumed to be large enough that by increasing them there is no significant influence on the added mass. The problem has been solved numerically with a zero-order boundary element method by discretizing the plate, the fixed wall and the flat ‘free surface’ with quadrilateral panels and using an indirect source formulation as in [13]. The added-mass calculations have been successfully assessed against reference solutions for the case of added mass in heave of a rigid plate and for a case of added mass for the elastic modes of a simply-supported uniform plate [1].

## 3. Physical investigations

Here, two cases of underwater explosions are examined, with specific focus on the evolution and characteristics of hydroelastic fluid–structure coupling.

### (a) Underwater explosion: interaction with a uniform plate

The examined case refers to the underwater-explosion experiment of Hung *et al*. [14] and the study is performed using both (i) the radial fluid-dynamic solution and approximate coupling with the structure and (ii) applying the DD strategy and the full fluid–structure coupling. In the tests, a mixture of DP60 detonator and Datasheet causes an underwater explosion interacting with an air-backed aluminium plate mounted on a steel frame (figure 1*a*). Different standoff distances *r*_{p} of the charge centre from the plate were studied but the time histories of the incident-wave pressure and induced plate acceleration are only available for the distance of 0.7 m. Hence, the latter is examined here. To reproduce correctly the incident-wave pressure, the parameters for the EOS (2.2) and, in particular, the initial radius, pressure and density of the gas cavity resulting from the explosion should be identified. The challenge is the unconventionality of the explosive; therefore the following identification strategy for the equivalent TNT parameters was implemented. In [14], Hung and co-authors assumed a radial symmetry of the explosion because of the large distance of its centre from the tank boundaries and expressed the incident-wave pressure at a given radial distance as follows:
3.1

Here, *W* is the charge weight and *r* the standoff distance of the chosen fluid location from the core of the explosion. Moreover, empirical formulae were proposed for the maximum pressure and the pressure rate of decay *θ* as functions of *W* and *r*.

In our model, we take as valid the EOS (2.2) both for water and gas, but, as we do not know the values of the involved coefficients for the used explosive mixture, as well as all required initial conditions, we model the explosion as due to an equivalent TNT charge, for which the coefficients of the related EOS are known. The equivalent explosion is identified enforcing that obtained by using the TNT charge at a given distance *r* is the same as that provided by equation (3.1). In particular, *r*=0.0895 m was chosen as the target location because it is relatively close to the explosion source and, further, is the smallest distance for which tabulated values for and *θ* are provided in [14]. This resulted in a cavity with initial radius *r*_{0}=0.03 m, density *ρ*_{0g}=1630.0 kg m^{−3} and pressure *p*_{0g}=2×10^{10} Pa. Figure 2 shows the comparison between the incident-wave pressure as predicted by the radial numerical solution and by the empirical formula (3.1) at four radial distances from the initial explosion. The results reveal a good description of the pressure decay in time. With reference to the maximum pressure, as expected, the two predictions are identical at *r*=0.0895 m while the present solution tends to underestimate more and more the empirical for increasing distances from the charge. This does not mean that the numerical solution necessarily provides a less accurate pressure prediction moving far from the initial explosion. In fact, the good comparison with the incident-wave pressure observed at *r*=0.7 m, which corresponds to the radial distance of the plate, highlights the good performances of the model. In more detail, figure 3 reveals that, once the rising times of the two solutions are synchronized (this involving a shift of 0.018 ms), the time history of the numerical pressure reproduces very well the experimental one of Hung *et al*. [14]. This means that the numerical explosion has a small error in terms of propagation speed, i.e. it is slightly faster than the physical one, but it can capture better the maximum experimental pressure than the empirical formula proposed in [14] and behaves similarly to such a formula in terms of time decay. Moreover, despite the simplifications used in modelling the explosion, the numerical prediction of 30.5 ms for the first oscillation period of the gas cavity (figure 1*b*) is not far from the value of 26.7 ms observed in the model tests.

The interaction between the underwater explosion and the target plate was studied using, for the fluid dynamic problem, both the full radial solution and the DD strategy. In the latter case, the radial solver extended from the position of the initial explosion up to a radial distance of 0.42 m and the three-dimensional sub-domain overlapped with it within a layer of about four local cells and, thus, included the plate with the steel casing. This means that the gas cavity was always inside the radial sub-domain during the examined evolution (figure 1*b*), therefore the one-dimensional solver simulated a two-phase flow while the three-dimensional method examined a one-phase flow interacting with the structure. In the case of the full radial solution, only the aluminium plate is modelled as a uniform elastic plate; in the DD case, also the steel casing is modelled, but, for numerical convenience (i.e. to make the thickness to be a multiple of the block size), it is assumed slightly thicker than it is in reality, with a thickness of 0.35 m. Since during the examined evolution, the pressure wave does not reach the back side of the real steel casing, the error connected with using a thicker structure is expected to be negligible. In addition, to limit the CPU-time requirements, the symmetry of the problem is accounted for and only one-fourth of the fluid domain is simulated in the three-dimensional sub-domain. In both cases, i.e. using the one-dimensional or the DD fluid-dynamic solver, the simulations were carried out assuming either a restrained or a simply-supported plate with *N*_{x}=*N*_{y}=10 beam modes, leading to *N*=100 plate modes.

In the experiments, the acceleration at the centre of the plate was measured during the initial stage of the wave–body interaction. This is examined in figure 4, which shows a much better comparison of both the radial solution with approximated fluid–structure coupling and the DD solution with full fluid–structure coupling with the experimental acceleration than that given by a recent numerical simulation documented in [15], which employs an analytical technique for the fluid-dynamic problem and an elastic dynamic response of the structure. In particular, both the solvers here proposed: (i) consistently reproduce the rising phase of the measured acceleration and (ii) indicate the inception of oscillations which characterize the experimental data. In addition, the measurements indicate a high-frequency behaviour, which cannot be explained on the basis of the information available from the work in [14]. The one-dimensional and three-dimensional solutions show a similar behaviour during the first stages of the evolution, involving a fast rise of the acceleration up to a maximum value which is followed by a decay. This suggests that at the beginning of the interaction, the radial symmetry is still valid and the coupling between fluid and structure is weak. Data of maximum acceleration at the centre of the plate are summarized in table 1 together with the maximum velocity at the centre of the plate, both referring to the initial phase of the incident-wave pressure interaction with the structure. As in the experiments only the acceleration was measured, and the maximum velocity has been estimated by integrating in time the measured acceleration. In this context, the high-frequency behaviour of the recorded acceleration represents a source of error both for the maximum acceleration and for the maximum velocity. An inspection of the results reveals that both one-dimensional and DD maximum accelerations agree better with the measurements in [14] than the finite element method (FEM) predictions reported in the same paper. The DD results are slightly worse than the one-dimensional predictions. However, such a slightly poorer performance can be explained in view of: (i) the high-frequency oscillatory behaviour in the measurements, (ii) the sensitivity during this first stage of the fluid–body interactions to the physical conditions, and (iii) the numerical choices. The same arguments hold also for the maximum velocity, which is overpredicted by the DD solution. In this case, an additional source of error is related with the time-integration procedure possibly used for the measured accelerations.

The time evolutions of the one-dimensional and DD accelerations start to diverge as the curves approach the first minimum and this leads to clear differences in the later oscillations (figure 4). At this stage, the two DD solutions, assuming either restrained or simply-supported plate, are more consistent with the experiments, better capturing the second oscillation peak and indicating a more complex fluid–structure interaction than that described by the one-dimensional solver.

To better understand the features involved in such a coupling, the DD solution is examined in further detail. The restrained-plate condition is focused upon, as representative of both plate edge conditions (restrained and simply-supported) which behave rather similarly. Figures 5 and 6 illustrate the evolution of the pressure iso-surfaces with *p*=73.6×10^{4} Pa together with, respectively, the pressure and normal velocity on the plate.

As time goes by, the wave pressure interacts with the plate, leading to pressure and velocity which increase on the structure with radial symmetry (see top-left and top-centre panels of figures 5 and 6). At this stage, the acceleration at the plate centre is still limited, but quickly increases (see time labelled as ‘1’ in figure 4). The radial symmetry loosens in time due to the finite dimensions of the plate (see top-right panels of figures 5 and 6). This corresponds, approximately, to the first minimum of the acceleration at the plate centre (see time labelled as ‘2’ in figure 4). When the pressure wave reaches the edges of the plate, the pressure and velocity of the structure are characterized by important contributions from the higher modes (see bottom panels of figures 5 and 6) and this suggests the possible need of more modes in the plate modelling, to properly describe the structural behaviour. An identical need can be found for the simply-supported plate, as confirmed by the similar evolution of the acceleration in figure 4. The need for the increase of modes is confirmed by the solution of figure 7, obtained with *N*=400 plate modes and for two different discretizations. There the acceleration of the plate is compared with the experimental one. The initial behaviour of the acceleration is similarly represented by the two simulations; the shorter duration of the initial peak obtained with the finer mesh is due to the lower numerical dissipation undergone in this case. After *t*=0.47 ms, higher modes play an important role (see also figure 4) and cause a weaker deceleration of the plate. For higher modes, a finer mesh is necessary to capture the hydroelastic coupling.

At this stage, the acceleration at the plate centre reaches a second maximum, which is smaller than the first peak due to the hydroelastic effects (see time labelled as ‘3’ in figure 4).

The reflected compression wave and the propagating one have different intensities because of the effect of the plate motion on the former. Moreover, the transmitted wave changes its curvature, in particular it tends to become cylindrical with axis around the edge of the plate. The combination of these three ingredients causes the interaction between the reflected and transmitted wave to determine the profile illustrated in the bottom–right panels of figures 5 and 6.

At this stage (see the bottom-right panel of figure 6), the plate velocity is very large near the corners as a consequence of the pressure wave reflection from the wall. The last instant shown for the DD prediction of the acceleration at the centre of the plate (see time labelled as ‘4’ in figure 4) corresponds to the last instant of the pressure iso-surfaces evolution discussed here. This corresponds to the time when the pressure wave reflected from the structure has almost reached the interface between the one-dimensional and three-dimensional sub-domains, so the simulation has been stopped. To prolong the solution in time, a dynamic enlargement of the three-dimensional sub-domain at the expenses of the one-dimensional sub-domain should be used, as discussed in §2.

### (b) Underwater explosion: interaction with a navy-ship bottom

Here, the possible consequences of an underwater explosion on the bottom of a navy ship are examined using the DD decomposition and the full fluid–structure coupling.

The bottom grillage arrangement is chosen to be that of configuration 1a among those studied experimentally in [16] (see the sketch of figure 8) and an equivalent orthotropic plate has been identified with main parameters given in table 2 and taken to be restrained at all edges.

The modelled plate has been forced by the underwater explosion caused by a TNT charge documented in [17], which is used as a representative sample of typical explosions, to assess possible ship damages. For this explosion, all parameters needed for the gas and liquid EOSs (2.2) are available and, in particular, the initial conditions of the gas cavity in terms of pressure, density and radius are, respectively, *p*_{0g}=8.381×10^{9} Pa, *ρ*_{0g}=1630.0 kg m^{−3} and *r*_{0}=0.16 m.

A systematic parametric analysis, in terms of the standoff distance *r*_{p}, i.e. the minimum distance of the plate from the initial charge, has been performed in [1] using the radial solver for the fluid-dynamic problem combined with the approximate fluid–structure coupling. For the shortest distance, *r*_{p}=4 m, stresses much larger than the yield value were recorded, suggesting inception of plastic deformations and possible damage of the structure. A decrease of the plate pressure down to the vapour pressure was also documented, suggesting some risk of cavitation. Here, the condition *r*_{p}=4 m is further examined by the DD solver by means of the full fluid–structure coupling. Although this is an academic study, with major simplifications of the geometrical and structural properties, it can still provide useful insights of the phenomena involved when a strong fluid–structure interaction occurs due to the closeness of the blast to the structure.

The left panel of figure 9 illustrates the solution in terms of iso-surfaces of pressure with *p*=10 Pa and contour plots of the plate deformation at time *t*≃2.2 ms, corresponding to the arrival of the acoustic wave on the plate. Taking advantage of the symmetry, only one-fourth of the physical domain is simulated, thus reducing the CPU time and memory space requirements. The plot also indicates the radial extension of the radial (one dimension) sub-domain, which is set within a distance of 1 m from the centre of the initial explosion, while the three-dimensional sub-domain represents the rest of the fluid domain containing the plate. On the right of figure 9, the evolutions of deformation (top sub-panel) and pressure (bottom sub-panel) at the middle of the plate are illustrated. From the results, it is clear that the impact of the compression wave on the wall induces a rapid pressure rise and an upward velocity of the plate. At this stage, the plate displacement is, on the other hand, rather small and tends to increase monotonically. As time goes by, the hydroelastic coupling becomes stronger, the pressure reaches a maximum value and then decreases rather quickly, becoming negative, i.e. smaller than the ambient pressure, at *t*≃2.7 ms and reaching values close to the vapour pressure, which suggests the possible occurrence of cavitation. Then the pressure rises again.

As the incident wave interacts with the structure, the radial symmetry in the pressure iso-surfaces is lost and the iso-surfaces stretch along the main direction of the plate, due to the lower influence of the restrained plate condition with respect to the orthogonal axis (figure 10).

This is clearer from the top view of the pressure iso-surfaces, showing better the stretching of such surfaces and, therefore, the growth of three-dimensional effects in the pressure (compare the middle and right panels of figure 10). Here, the lighter iso-surface corresponds to *p*=0 Pa. Its dimension along the plate main axis is about one-tenth of the plate length, and its presence highlights the inception of a cavitation zone on a smaller surface of the plate inside this iso-surface. Once cavitation occurs, a state change, from liquid to gas, is forced and this affects the structural loads. These phenomena cannot be handled by the present method and require an extension of the solution strategy.

Figure 11 illustrates the same case, by comparing the fully coupled solution with the approximate coupled fluid–structure solution described in §2. During the pressure rise, the two results are very close and predict similar maxima, while later on the approximate solution predicts a pressure decrease with lower rate. As a result, the inception of cavitation is not predicted as it is by the fully coupled solution. Only at a later stage (not shown here) does the pressure from the approximate solution reach values close to the vapour pressure.

## 4. Conclusion

A hybrid numerical solver has been proposed to investigate underwater explosions occurring sufficiently far from the sea floor and the free surface so as to lead to radial symmetric acoustic waves and their interactions with marine structures. The method is based on a one-way time–space DD strategy coupling a radial and a fully three-dimensional solution for compressible inviscid multiphase flows. This is done in a way to achieve a good compromise between capability, accuracy and efficiency. The three-dimensional solver is based on an adaptive mesh refinement technique and a parallelization process which helps to reduce the computational costs (CPU-time and memory-space requirements). The fluid–structure interaction has been implemented both as an approximated coupling and as a full coupling so as to highlight the importance of hydroelastic phenomena. The target structure was modelled as a uniform or orthotropic plate. The latter is a reasonable approximation of the grillage arrangement on the bottom of surface ships. The fluid-dynamic model was used both as a purely radial solver and as a DD-based solver, to investigate two underwater explosion cases. The first example refers to a uniform plate and allowed one to successfully validate the solver for the acoustic (compressible) phase of the fluid–structure interaction. The second case examined an idealized ship bottom and provides a significant assessment of the importance of hydroelastic effects and related consequences. Presently, only the compressible phase has been focused upon, which is relevant for the local consequences on the structure. The later incompressible stage, with eventual gas cavity oscillating and imploding on the structure, is important for global consequences on the vessel and will be examined in a future step of the present research activity. However, excluding the cavity collapse in very tiny bubbles, such a later phase is considered to be less demanding in terms of modelling and implementation, as the authors have already experience in this framework (e.g. [18]). Before moving in this direction, the dynamic DD algorithm will be assessed and the passage of the gas cavity from the radial to the three-dimensional domain will examined, because these two aspects are required to investigate the later incompressible stage.

## Funding statement

This research activity is partially funded by the Research Council of Norway through the Centres of Excellence funding scheme AMOS, project no. 223254, and partially by the Flagship Project RITMARE—the Italian Research for the Sea—coordinated by the Italian National Research Council and funded by the Italian Ministry of Education, University and Research within the National Research Program 2011–2013.

## Footnotes

One contribution of 12 to a Theme Issue ‘Advances in fluid mechanics for offshore engineering: a modelling perspective’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.