## Abstract

The present paper deals with the methods for the evaluation of the hydroelastic interactions that appear during the violent sloshing impacts inside the tanks of liquefied natural gas carriers. The complexity of both the fluid flow and the structural behaviour (containment system and ship structure) does not allow for a fully consistent direct approach according to the present state of the art. Several simplifications are thus necessary in order to isolate the most dominant physical aspects and to treat them properly. In this paper, choice was made of semi-analytical modelling for the hydrodynamic part and finite-element modelling for the structural part. Depending on the impact type, different hydrodynamic models are proposed, and the basic principles of hydroelastic coupling are clearly described and validated with respect to the accuracy and convergence of the numerical results.

## 1. Introduction

Sloshing became a very important practical problem in the last decade due to the increased activities in liquefied natural gas (LNG) transport. Large numbers of LNG carriers were built or are under construction with capacities that are almost double when compared with the classical LNG ships (from 138 000 to 240 000 m^{3}). The most common LNG ships belong to the so-called membrane type with some typical examples of the containment systems (CSs) of the tanks shown in figure 1. Within the membrane-type concept, which is of main concern here, the LNG is kept liquid at very low temperature (−165^{°}C) by a complex insulation system, which is attached to the ship structure.

At the same time as the size of LNG vessels increased, the operational requirements became more and more severe. Indeed, in the past, LNG ships were allowed to operate either in full or empty tank conditions, while today there is a necessity to allow for sailing at any partial filling. This requirement introduces serious difficulties in the design of both the CS and the associated ship structure. Violent sloshing motions may occur (figure 2) and a direct consequence is the occurrence of different impact situations, which may induce extreme structural loadings being devastating for both the CS and ship structure [2]. Correct numerical modelling of the fluid–structure interactions during sloshing impacts is extremely complex, and it is fair to say that, up to now, there is no fully satisfactory numerical model able to treat these situations in a fully consistent manner [1,3]. Even without considering the interaction with the structure (hydroelasticity), the modelling of the pure fluid flow raises serious problems owing to several complex physical phenomena involved (rapid change of the free surface geometry, two (three) phase flow in some situations, gas cushion, low temperature of the LNG (−165^{°}C), important three-dimensional effects [4–6], compressibility [7–9], surface tension, viscous effects and ullage pressure). In addition to these pure fluid mechanics problems, and owing to the flexibility of the CS, another important aspect, which seems to be essential for correct evaluation of the structural responses, is the effect of hydroelasticity. Indeed, owing to the violence of sloshing impacts, the hydrodynamic pressure will often depend on the structural response, that is, fully coupled hydrostructure modelling is necessary [10–15]. In order to better understand the modelling difficulties related to hydroelasticity, in figure 1 two typical CSs that are in use today are shown. The first one is the so-called NO96 system, which is composed of plywood boxes filled with perlite, while the second system, called MARK III, is composed of several levels of foam combined with plywood structure. On the side in contact with LNG, both systems have a membrane made of special metal alloy: Invar for NO96 and stainless steel for MARK III. In the case of NO96 CS, this membrane is flat, while it is corrugated for MARK III CS. Correct structural modelling of such a complex structure is still challenging even for most sophisticated numerical tools based on the well-mastered finite-element method (FEM).

There are several methods in use to solve this problem, each of them having its own drawbacks. The most common method is based on small-scale model tests and its main drawback is the impossibility to derive clear scaling laws for pressure distributions in order to extrapolate them to full scale. In addition, it is not possible to model hydroelasticity within the small-scale model tests and there are also difficulties related to the measurements of the pressure.

The second method could be described as quasi-full-scale measurements. The idea is to use a real CS and impact it with water. This is usually done either simulating drop tests or by impacting the CS with waves in a wave flume. The drawbacks of the method are: impossibility to simulate all possible impact situations, difficulties in measuring the pressures and structural response, impossibility to use LNG as a fluid, and relatively high cost owing to the complex installations.

Both small-scale model tests and quasi-full-scale experiments show very high variability of results, and repeatability of the results is very difficult to achieve. This is because the impact characteristics are extremely sensitive to details of the fluid flow.

The third method, which is used to assess the fluid–structure interaction during sloshing impacts, belongs to mathematical modelling. Numerical modelling of coupled hydrostructure interactions during sloshing impacts is a very challenging problem from both hydrodynamic and structural sides. Indeed, even decoupled, these two problems are very difficult to model properly. Although some attempts have been made to solve the three-dimensional impact problems, the two-dimensional modelling of fluid flow is used most often. A main reason for that are the difficulties associated with the determination of the free surface flow and the exact wetted part of the structure during the impact. In addition, presence of air and possible air/fluid mixing introduce some additional difficulties. The three-dimensional effects of the structural response can be treated by standard FEM codes provided correct characteristics of the CS are available. The determination of the FEM characteristics is far from trivial owing to the complexity of the CS (plywood, foam, steel and mastic ropes) and the associated ship structure. As far as the fluid flow is concerned, the numerical methods, which are used most often in practice, can be subdivided into the pure computational fluid dynamics (CFD) methods and the semi-analytical methods. Owing to strong variation of the free surface during sloshing, the most popular methods belong to the family of the volume of fluid technique and to the so-called smoothed particle hydrodynamics (SPH) method. SPH method has the advantage of being gridless allowing for very strong free surface variations. However, all CFD methods suffer from numerous numerical problems when it comes to the evaluation of the highly localized pressures. The mesh requirements for proper evaluation of the hydrostructure interactions during the impacts become prohibitive, and the stability of different numerical schemes is hard to ensure, especially when hydroelastic analysis needs to be performed. The numerical schemes are usually so time-consuming that it is unfeasible to use them for statistical estimates of tank responses.

Semi-analytical impact models represent another type of methods for sloshing impact problems [1,16]. The idea is to identify the most typical impact situations and then simplify them in order to be able to describe them with simple geometries including a few most important physical parameters only. As far as the fluid flow generated by impact is concerned, both three-dimensional [6,17,18] and two-dimensional [2,16] models and solutions are available for rigid structures. For elastic structures, the fluid flow during impact is difficult to calculate even within two-dimensional semi-analytical models. Nevertheless, hydroelastic behaviour of a conical shell impacting on an incompressible free surface was studied in detail by Scolan [19] within his original semi-analytical model. An approach based on variational inequality formulation of impact problems was used in Gazzola *et al*. [4,5] to solve three-dimensional problems of hydroelastic impact. In view of all the previous assumptions and uncertainties in the estimation of the sloshing impact situations, the two-dimensional approximation of the fluid flow during impact phase appears to be reasonable. However, the two-dimensional assumption of the structural behaviour can hardly be justified because the CS is extremely complex and fundamentally three-dimensional. Within the semi-analytical approach, the three-dimensional effects can be taken into account using the so-called strip technique, which consists of considering the fluid flow as two-dimensional in each strip and considering the structure as three-dimensional. In this way it is possible to couple the semi-analytical two-dimensional methods with a general three-dimensional structural software. These semi-analytical or semi-numerical methods are the main subject of the present paper.

## 2. Methodology

The main idea of the present approach is to ‘take’ good parts of all available hydrodynamic and structural tools and combine them into a well-controlled procedure able to take properly into account the most important physical aspects of the fluid impact with the aim to identify conditions that are the most dangerous from the point of view of the structure integrity.

Global sloshing flow for prescribed ship loading conditions is suggested to be evaluated by combining a seakeeping tool with three-dimensional CFD code and/or small-scale model tests. These simulations can be performed without accounting for flexibility of both the CS and ship structure. Results of global simulations are used to identify the most dangerous places in the tank and impact configurations induced by the violent sloshing motions.

We should keep in mind that the end result of the analysis is not the hydrodynamic pressure but the structural response. Indeed, even if the measured or calculated pressures are extremely high, they are not necessarily dangerous for the structure, because the structural response depends not only on the maximum pressure values but also on spatial and temporal distribution of the pressure and on structural characteristics such as natural periods and the structural damping.

Once the impact conditions have been identified, local hydroelastic analysis is performed using asymptotic models of fluid flow combined with a commercial FEM tool for structural response. The local hydroelastic analysis is applicable only during the impact stage, when the hydrodynamic loads are high and the elastic response of the CS is maximal. Both experiments and numerical computations show that the impact stage is of short duration and impact loads are localized in space. This observation makes it possible to disregard many physical effects, which are of concern in the CFD analysis, such as large dimensions of the tank and its real shape, real profile of the free surface at a distance from the impact region, viscosity of the fluid, its surface tension and gravity effects. However, some effects, which are believed to be of minor importance in the CFD analysis, should be taken into account in the local analysis. These effects are compressibility of the fluid, presence of gas above the fluid surface and in the impact region, aeration of the fluid in the impact region, jetting and fine details of the flow in the jet root region, rapid increase of the wetted surface of the tank wall and the flexibility of the wall. The short duration of the impact stage allows us to simplify the local analysis and to use a combination of analytical and numerical methods instead of direct numerical calculations. The analytical part of the local analysis is very important because it allows us to:

— obtain useful formulae suitable for design needs,

— control numerical results,

— treat properly the coupled problem of fluid–structure interaction during the impact, and

— determine the wetted part of the wall simultaneously with the fluid flow and the pressure distribution.

The last point is crucial in the context of Wagner-type impact [20–22], when the rate of the wetted area increase is greater than the accelerations of the liquid particles and standard schemes of integration in time used in CFD become unappropriate.

In our approach it is suggested to use simplified hydrodynamic models in combination with complex structural models during the impact stage. This idea is based on the experience gained in both theory and applications that semi-analytical models of violent flows during impact stage are comparable in terms of their performance with fully nonlinear calculations with high resolution in space and in time. In many practical cases, both the impact conditions and the level of the fluid aeration in the impact region are not well defined, and small variations of the global conditions may lead to significant variations of the local impact conditions. This is why, in some sense, attempts to reproduce all details of the flow such as the shape of the flow region during impact and all fluid characteristics have no meaning in practice, although they may lead to very interesting mathematical problems.

Three main types of the impact during violent sloshing are distinguished [1,23].

— Wagner-type impact (figure 3

*a*).— Steep wave impact (figure 3

*b*,*d*). In this scenario, the front of the wave is almost vertical. This type of impact occurs when the wave is overturning just before its impact onto the wall. The fluid on and behind the wavefront can be aerated before impact (figure 3*d*) or not (figure 3*b*).— Breaking wave impact, or Bagnold-type impact (figure 3

*c*). In this scenario, the wave is overturning before the wavefront comes into contact with the wall and the cavity filled with gas is closed at the impact instant.

Note that the geometry of the flow region at the impact instant is intentionally simplified in order to reduce the number of parameters of the problem. Parametric analysis of the fluid–structure interaction can be performed for each type of wave impact in order to distinguish the most dangerous impact conditions and to study the effects of these parameters on the structural response.

In this paper, the main concepts of local analysis are demonstrated for each type of wave impact. It will be shown that the suggested combination of semi-analytical models for hydrodynamics and finite-element models of a complex structure leads to a system of ordinary differential equations with respect to the so-called principal coordinates of the structural deflection. This system is truncated and solved by numerical integration. The procedures leading to such systems of ordinary differential equations are detailed for each semi-analytical model in §§4 and 5. In §3, formulation of the semi-analytical model of Wagner-type impact is presented for a deformable wall but analytical results are shown only for a rigid wall. Numerical results are shown in §6 only for the most difficult case of steep wave impact with aeration.

## 3. Wagner-type impact

In this type of wave impact, the wetted area of the vertical wall changes in time during the impact stage [20,21]. We need to find the hydrodynamic loads acting on the wall together with elastic deflection of the wall and the position of the contact region between the fluid and the wall. The variation in time of the contact region makes the problem nonlinear in contrast with other types of wave impact. In this section, we are concerned with the hydrodynamic part of the impact problem assuming that deflection of the wall is given. Only the two-dimensional model is considered.

Two-dimensional unsteady flow in the presence of gravity is governed by the following equations with respect to the velocity potential *Φ*(*x*,*y*,*t*) and the position of the fluid free surface *F*(*x*,*y*,*t*)=0:
3.1
Here *Ω*(*t*) is the flow region at time instant *t*. The fluid boundary ∂*Ω* consists of the following three parts: ∂*Ω*=∂*Ω*_{f}∪∂*Ω*_{v}∪∂*Ω*_{H}. Dynamic and kinematic conditions are satisfied on the free surface ∂*Ω*_{f}(*t*), where *g* is the acceleration owing to gravity. On the wetted part of the wall ∂*Ω*_{v}(*t*), the horizontal velocity of the flow is equal to the normal velocity *w*_{t}(*y*,*t*) of the wall, which can be elastic. Vertical velocity is equal to zero along the horizontal rigid bottom ∂*Ω*_{H}. At *t*=0, the velocity field of the flow ∇*Φ*_{0}(*x*,*y*) and position of the free surface, *F*_{0}(*x*,*y*)=0, are assumed to be given.

When a wavefront with an almost vertical face approaches the vertical wall, the curvature of the free surface at the water line increases. Usually numerical methods fail to simulate properly the liquid flow near such points and analytical study is required. The stage of violent interaction between the wall and the wave is short but this is the stage during which the hydrodynamic pressures take their maximal values.

We shall determine the liquid flow and the pressure distribution during the impact stage under the following assumptions: (i) at *t*=0, the liquid region near the wall can be approximated as a half-strip *D* with depth *H*; (ii) normal deviation of the boundary ∂*Ω*(0) from ∂*D* is of the order of *ϵH*, where *ϵ* is a small non-dimensional parameter; (iii) liquid is ideal and incompressible; (iv) normal velocity ∂*Φ*/∂*n* on the free surface ∂*Ω*_{f} is finite and |∂*Φ*/∂*n*|<*V* along the front of the wave; (v) duration of the impact stage is of the order *ϵH*/*V* ; and (vi) the pressure along the free surface is constant.

We introduce non-dimensional variables and unknown functions as follows:
3.2
where prime stands for non-dimensional variables, *ρ* is the fluid density and *p* is the hydrodynamic pressure. Assumption (i) indicates that the deformation of the wavefront can be described by the equation *x*′=*ϵζ*′(*y*′,*t*′) where *ζ*′(*y*′,*t*′) is a smooth non-dimensional function and both *ζ*′(*y*′,0) and *ζ*′_{t}(*y*′,0) are given. The kinematic condition on the wavefront has the form
3.3
When the wavefront approaches the wall, the water in front of the wave is accelerated upwards. Even for a very short initial stage, which is referred to as the impact stage, the rate of increase of the wetted part of the wall can be significant. This conclusion comes from the assumption that the initial gap between the wavefront and the wall is narrow and its thickness decreases in time. At *t*′=0, the wetted part of the wall corresponds to the interval [0,*b*′(0)]. Water level on the wall moves upwards with the length of the wetted part of the wall being *b*′(*t*′)*H*. The function *b*′(*t*′) is unknown in advance and has to be determined together with the flow and pressure distribution.

Taking formally *ϵ*=0 in the original equations written in the non-dimensional variables, omitting primes and introducing new unknown functions *ϕ*(*x*,*y*,*t*) and *h*(*y*,*t*) as (see [20])
3.4
where *v*(*y*)=−(∂*Φ*_{0}/∂*x*)(*ζ*_{0}(*y*),*y*), we arrive at the following boundary-value problem:
3.5

The boundary-value problem (3.5) can be considered as a water-entry problem for a blunt body, the initial shape of which is given by the equation *x*=−*ζ*_{0}(*y*), penetrating the half-strip of the liquid, *x*>0, 0<*y*<1, at the velocity *v*(*y*)+*w*_{t}(*y*,*t*), where *w*(*y*,*t*) is non-dimensional elastic deflection of the wall. The velocity is different, in the general case, for different points of the body. This means that the fictitious body is deformable, even if *w*≡0, but its shape is prescribed in advance. If the wall is elastic and the function *w*(*y*,*t*) is to be determined together with the pressure distribution, the problem (3.5) has to be supplemented with an equation of the wall deflection and corresponding initial and boundary conditions. The resulting hydroelastic problem can be reduced to a system of ordinary differential equations with respect to the so-called principal coordinates of the structural deflection in a similar way as has been done in Korobkin & Khabakhpasheva [24]. Below some results are shown only for a rigid wall.

For a rigid wall, the function *b*(*t*) is defined by the equation
3.6
where
3.7
This equation follows from the Wagner condition that the free surface deflection *h*(*y*,*t*) is finite at the contact point *y*=*b*(*t*).

The pressure distribution along the wetted part of the rigid wall, 0<*y*<*b*(*t*), reads [20]
3.8
Formula (3.8) shows that the pressure is a square root singular at the contact point *y*=*b*(*t*). The present approach fails to describe the fine flow close to the moving contact point. In particular, the predicted pressure is beyond all bounds as
3.9
where
3.10
Inner asymptotic solution is required to resolve the flow details near the contact point. If the wall is elastic, the problem can be solved by the method of normal modes as is outlined in the following sections.

## 4. Steep wave impact

A main difference of the steep wave impact (figure 4) from the Wagner-type impact is due to the inclination angle of the wavefront from the wall. The wavefront is approximated as vertical within the steep wave impact model. That is, the inclination angle is set to be zero (figure 3*b*,*d*). Steep wave impact can be viewed as a limiting case of the Wagner-type impact [11,23].

There are two models of steep wave impact: without and with mixing of the fluid with the surrounding gas. The first model can be viewed as a model of hydraulic jump propagating towards the wall or a model of a wave that is getting very steep before arriving at the wall. Within the second model, we consider the so-called broken wave with the fluid behind the wavefront being well mixed with the gas. In both cases, the fluid instantly comes in contact with either rigid or elastic wall with a non-zero speed over a region of non-zero area. Compressibility of the fluid should be taken into account to describe this type of impact. The fluid compressibility plays a more significant role in the second model, because aeration reduces the sound speed in the mixture compared with the sound speed of the pure fluid in the first model. In realistic sloshing conditions, the Mach number of the flow, which is the ratio of the wave impact speed to the sound speed in the fluid, is small. This makes it possible to employ the acoustic approximation, where the compressibility is described only by the value of the sound speed in the fluid at rest. The sound speed in the acoustic approximation can be non-uniform but it is independent of time.

Owing to short duration of the impact, both the boundary conditions and equations of motion can be linearized and the boundary conditions can be at the leading order imposed on the initial position of the fluid boundary just before the impact.

Mathematical formulations of the models and the involved effects with and without account for the fluid aeration are rather different. However, the models share the same configuration. In the three-dimensional formulation, the configuration of the flow region is given as follows (figure 4). The fluid of depth *H* occupies the region *x*<0, 0<*y*<*L* and 0<*z*<*H*, where the level *z*=0 corresponds to the flat rigid bottom, and the vertical *z*-axis is directed upward. Before impact, the part of the wall 0<*z*<*H*−*H*_{w} is in contact with the fluid. A wave with vertical front of height *H*_{w} approaches the wall at a constant velocity *U* in the region 0<*y*<*L* and *H*−*H*_{w}<*z*<*H* at the initial time instant *t*=0. Part of the wall *S*_{e} is elastic, and the rest of the wall is rigid.

### (a) Steep wave impact without aeration

Within the acoustic approximation, the flow is described by the velocity potential *φ*(*x*,*y*,*z*,*t*). The problem is formulated in non-dimensional variables that are chosen in such a way that the fluid depth, impact velocity, fluid density and sound velocity in the fluid are equal to unity. The boundary-value problem with respect to the velocity potential has the form [3]
4.1
where *l*=*L*/*H*, *w*(*y*,*z*,*t*) is the non-dimensional deflection of the elastic part of the wall *S*_{e}, *q*(*y*,*z*,*t*) is the hydrodynamic pressure distribution over the wall, and *χ*_{0}(*y*,*z*) and *χ*_{1}(*y*,*z*) are characteristic functions of the impact region and the elastic region of the wall, respectively. In particular, *χ*_{1}(*y*,*z*)=1, where (*y*,*z*)∈*S*_{e}, and *χ*_{1}(*y*,*z*)=0 on the rest of the wall.

Equations (4.1) have been derived using the original equations of hydrodynamics and original boundary conditions, assuming special geometry of the flow region before impact and short duration of the impact stage under consideration [1]. Deflection *w*(*y*,*z*,*t*) is determined numerically using a commercial code for structure dynamics and the modal approach detailed below.

The velocity potential *φ*(*x*,*y*,*z*,*t*) and the deflection *w*(*y*,*z*,*t*) of the flexible surface *S*_{e} are sought in the forms of modal expansions [7]:
4.2
where *V*_{n}(*y*,*z*) and *W*_{k}(*y*,*z*) are the orthonormal eigenmodes for the cross section of the flow domain and for the elastic part of the wall, respectively. The corresponding eigenvalues are *λ*_{n} and *μ*_{k} [3].

Substituting equation (4.2) into equation (4.1) after some algebra and application of the direct and inverse Laplace transforms provides the pressure distribution *q*(*y*,*z*,*t*) over the wall
4.3
Here *J*_{0}(*x*) is the zeroth-order Bessel function of the first kind.

General equations of structural analysis yield that the principal coordinates *a*_{k}(*t*) of the deflection *w*(*y*,*z*,*t*) of the elastic part of the wall *S*_{e} satisfy the following ordinary differential equations:
4.4
where *α*_{k}, non-dimensional eigenfrequencies *ω*_{k}, coefficients of the structural damping *γ*_{k}, and the modes (shape functions) *W*_{k}(*y*,*z*) are obtained using the commercial code ABAQUS. Once equation (4.4) has been integrated in time with zero initial conditions, the functions *a*_{k}(*t*) and *W*_{k}(*y*,*z*) are used in equation (4.2) to calculate the wall deflection *w*(*y*,*z*,*t*).

Equations (4.2)–(4.4) lead to a system of integral and differential equations with respect to the principal coordinates *a*_{k}(*t*) of the deflection:
4.5
where *k*≥1 and
4.6

Note that equations (4.5) serve to compute the principal coordinates *a*_{k}(*t*) of the structural modes *W*_{k}(*y*,*z*) without calculating the hydrodynamic pressure. Pressure can be computed later, if it is necessary, using equation (4.3). The elastic structure in these calculations can be very complex. Stresses at any point of the structure can be calculated once *a*_{k}(*t*) have been determined. If the eigenfrequency *ω*_{1} is relatively small, so that the natural period of the first structural mode 2*π*/*ω*_{1} is much greater than the time scale of the acoustic effects, equations (4.5) can be approximated by differential equations of incompressible fluid impact onto an elastic wall [11]. However, the hydrodynamic pressure during the impact stage cannot be computed within the incompressible fluid model.

### (b) Steep wave impact with air entrainment

In this type of wave impact, the fluid is aerated close to the impact region with air fraction *α*(*x*,*y*,*z*) being a given function of spatial coordinates. The aerated fluid is treated as a fictitious medium with reduced sound speed *c*(*x*,*y*,*z*) and density [14,15,25]. The air fraction is assumed to be zero at some distance *D* from the wall (figure 4*b*). The fluid is treated as homogeneous and incompressible beyond this distance. Compressibility of the aerated fluid close to the wall matters. The sound speed in the aerated fluid can be as small as 25 m s^{−1}.

The liquid flow owing to the impact is irrotational. The velocity potential satisfies the wave equation with variable sound speed:
4.7
where the sound speed *c*(*x*,*y*,*z*) is a given function of spatial coordinates *x*, *y* and *z* only in the region of fluid/air mixture and is set to be where *x*<−*D*/*H*=−*d*. Note that the variation of the fluid density owing to aeration is not taken into account in the present model. Non-dimensional variables are defined in the same way as in §4*a*. The problem (4.7) is similar to equation (4.1) without aeration. The only difference of the present model from the model without aeration is the variable coefficient *c*^{−2}(*x*,*y*,*z*) in the wave equation for the velocity potential. The presence of this coefficient may change the solution significantly.

The boundary-value problem (4.7) is solved by domain-decomposition method. The flow domain is divided into two parts: the region of incompressible fluid (‘external’ domain, *x*<−*d*) and the region of compressible aerated fluid (‘internal’ domain, −*d*<*x*<0). The conditions of continuity of the velocity potential and its normal derivative at the interface, *x*=−*d*, between the aerated and pure fluid and the condition on the wall provide boundary conditions for the wave equation in the internal domain.

The velocity potential in the internal domain *φ*(*x*,*y*,*z*,*t*) and the deflection of the elastic panel *w*(*y*,*z*,*t*) are sought in the forms
4.8
where the functions *V*_{m}(*y*,*z*) and *W*_{k}(*y*,*z*) were introduced in §4*a*. The time-dependent coefficients *ϕ*_{n}(*t*) and *a*_{k}(*t*) are to be determined by using the body boundary condition, conditions at the interface, *x*=−*d*, and the wave equation with non-uniform sound speed. The functions *u*_{n}(*x*) are the non-trivial solutions of the following eigenvalue problem:
4.9
where , and *λ*_{m} is the eigenvalue that corresponds to the eigenfunction *V*_{m}(*y*,*z*) (see §4*a*).

The structural response is described by equation (4.4), which is independent of the aeration effects.

Equations (4.7)–(4.9) and equation (4.4) lead to the following system of ordinary differential equations with respect to the time-dependent coefficients in equation (4.8):
4.10
where ** φ**={

*ϕ*

_{1},

*ϕ*

_{2},…},

**={**

*a**a*

_{1},

*a*

_{2},…},

**={**

*v**v*

_{1},

*v*

_{2},…},

*A*=diag{

*α*

_{1},

*α*

_{2},…}, ,

*D*=diag{

*d*

_{1},

*d*

_{2},…},

*S*={

*S*

_{kn}},

*T*={

*T*

_{kn}}, . The elements of the matrices

*S*are given by 4.11 and

*T*

_{kn}and

*v*

_{k}are given by equation (4.6).

The derived system of differential equations (4.10) can be written as
4.12
with extended unknown vector . The matrix *M* and vector ** C** are given by
4.13
Note that the matrix

*M*and the vector

**are independent of time. This makes it possible to derive the exact solution to the initial value problem (4.12). To do this, the matrix**

*C**M*is represented as a product

*M*=

*R*⋅

*L*⋅

*R*

^{−1}, where

*R*is the complex matrix made by eigenvectors and

*L*is the diagonal matrix formed by the corresponding eigenvalues of the matrix

*M*. Then, the solution is 4.14

Although the initial-value problem (4.12) has the analytical solution (4.14), this solution is only useful for verification of the numerical solver. This is due to the large number of unknowns in equation (4.12) after its truncation and, correspondingly, the large size of the matrices in (4.14). To find the numerical solution to the infinite system of ordinary differential equations (4.12), the system was truncated and integrated by using the fourth-order Runge–Kutta method.

## 5. Bagnold-type impact

This impact type occurs when a wave is overturning in front of the wall and the cavity filled with gas is closed at the impact instant [10,23,26,27]. The presence of the cavity is a main feature of this type of impact. After the time instant when the cavity is closed, the gas in the cavity starts to be compressed and oscillates thereafter. These oscillations are of relatively large period, which is defined by the volume of the cavity and the characteristics of the gas trapped in the cavity.

In general, this problem is very complicated and requires a solver with high resolution. However, during the initial stage of high hydrodynamic loads, the problem can be simplified assuming that the shape of the cavity does not change significantly during this stage and the cavity thickness is small compared with the dimensions of the flow region. This implies, in particular, that the projection of the cavity onto the wall does not vary in time during the stage of interest.

To introduce the model of wave impact with trapped gas cavity, we use the same configuration of the flow region as in §4*a*. The fluid of depth *H* and width *L* occupies the domain *x*<0, 0<*y*<*L* and 0<*z*<*H*, where *z*=0 corresponds to the flat rigid bottom, and the vertical *z*-axis is directed upward. A wave approaches the wall in the positive *x*-direction. The cavity, the projection of which on the wall is shown in figure 5, is closed at *t*=0. In figure 5, the striped area corresponds to the wave impact region, the grey rectangular region corresponds to the flexible part of the wall and the white rectangle corresponds to the cavity *S*_{c}. The part of the wall which is below both the impact region and cavity was wetted before impact. Note that the cavity *S*_{c}, elastic part of the wall *S*_{e} and the impact region are shown as rectangular, which is only for demonstration purposes. These regions can be of any shape and also multi-connected.

The wavefront *H*−*H*_{w}<*z*<*H* approaches the wall at velocity *U* and impacts the wall closing a cavity at *t*=0. The thickness of the cavity *s*_{0}*H* is assumed small and uniform, *s*_{0}≪1. Within the Bagnold model [26], the cavity projection onto the wall does not vary in time and the free surface of the cavity is approximated by a rigid massless plate that moves owing to the hydrodynamic pressure in the fluid and aerodynamic pressure in the cavity. Displacement of this massless plate is *ξ*(*t*)*HU*/*c*_{0} in dimensional variables. Here, *ξ*(0)=0, *ξ*(−0)=*U*_{0}/*U* and *c*_{0} is the sound speed in the fluid at rest. The distance from the wall to the plate is equal to the cavity thickness and is *s*(*t*)*H*=*s*_{0}*H*−*ξ*(*t*)*HU*/*c*_{0}, where functions *s*(*t*) and *ξ*(*t*) are non-dimensional.

The problem is formulated in non-dimensional variables that are chosen in such a way that the fluid depth, impact velocity, fluid density and the sound speed in the fluid at rest are equal to unity. The boundary-value problem with respect to the velocity potential has the form
5.1
where , *v*_{0}=*U*_{0}/*U*, *χ*_{2}(*y*,*z*) equals 1 in *S*_{c}, and zero otherwise, *p*_{c}(*t*) is the adiabatic pressure in the cavity, *γ* is the adiabatic index. Other elements in equation (5.1) were defined for the model of steep wave impact (4.1). Equations (5.1) should be combined with the equations of structural response (4.4). Vibration of the cavity is governed by the equation of the force balance
5.2
where *S* is the area of the cavity projection *S*_{c} onto the vertical wall.

Following the analysis outlined in §4*a* for steep wave impact without aeration, we arrive at the system of integral and differential equations
5.3
where *a*_{0}(*t*)≡*ξ*(*t*), *α*_{0}=0,
5.4

The system (5.3) is similar to the system (4.5) but now we have more unknown functions and more elements to compute such as *Q*_{nk} and *K*_{m0}.

## 6. Numerical results

In this section, results of numerical analysis of steep wave impact with aeration are shown. This type of wave impact is the most difficult in terms of numerical calculations and interesting in terms of physical effects involved. The dimensional variables are used in this section.

We assume that the dimensional sound speed *c*_{dim}(*x*,*y*,*z*) is a function of the vertical coordinate *z* only and is given by [28]
6.1
where the function *α*(*z*) describes the distribution of the gas concentration in the mixture region, *ρ*_{L} and *ρ*_{G} are liquid and gas initial densities, *ρ*_{L}=*ρ*, *γ* is the adiabatic index of polytrophic gas, *γ*=1.4 for air, and *p*_{0} is initial pressure of the gas, *p*_{0}=101 325 N m^{−2} in the present calculations, which corresponds to atmospheric pressure.

In the two-dimensional case, the calculations revealed that aeration can increase strain maximum in elastic structure during impact. The structure was modelled as a steel (*ρ*_{b}=7800 kg m^{−3}) beam of thickness *H*_{b}=0.02 m with Young’s modulus *E*=20.7×10^{10} N m^{−2} and structural damping *γ*_{b}=0.001 s. The fluid density and depth, wave height and impact velocity were *ρ*=1000 kg m^{−3}, *H*=2 m, *H*_{w}=0.25*H* and *U*=1 m s^{−1}, respectively. The beam of 1 m length was simply supported at *z*=1 m and *z*=2 m. Evolutions of bending stresses in the beam for three different levels of aeration are shown in figure 6 in microstrains. Calculations were performed for the aerated region of thickness *D*=0.25 m. The fluid is incompressible in *x*<−*D*. The aerated fluid in −*D*<*x*<0 is treated as a fictitious acoustic medium with uniform sound speed *c*_{dim}=47.14 m s^{−1} (figure 6*a*) and *c*_{dim}=95.46 m s^{−1} (figure 6*b*), and with a non-uniform sound speed *c*_{dim}(*z*) in figure 6*c*. In calculations, *c*_{dim}(*z*)=1500 m s^{−1}, where 0<*z*<1.5 m, with , where 1.5<*z*<2 m and m s^{−1} or 95.46 m s^{−1}.

In figure 6*a*,*b*, the solid lines show evolutions of the bending stresses for the compressible layer of finite thickness, *D*=0.25 m, while the dotted lines are for the compressible fluid of infinite extent, . The values of the sound speed *c*_{dim}=47.14 m s^{−1} and *c*_{dim}=95.46 m s^{−1} in the calculations are such that the corresponding acoustic time scales are close to the ‘wet’ period and ‘dry’ period of the first natural mode of the beam, respectively. It is seen that the stresses for the finite region of aeration may exceed those for the infinite region (figure 6*a*).

Figure 6*c* shows that a variation of the gas concentration along the wall may increase the bending stresses in the wall by almost 50 per cent.

Calculations were performed also for three-dimensional steel boxes shown in figure 7. The boxes are made of steel plates with Young’s modulus *E*=20.7×10^{10} N m^{−2}, density *ρ*_{b}=7800 kg m^{−3}, Poisson ratio *ν*=0.3 and structural damping *γ*=0.0001 s. The boxes were simply supported along the top and bottom edges while the side edges were free. The first box, 1.655 ×1.655 m, is opened at the back with thickness of all plates being 1 cm (figure 7*a*). The second and third boxes, figure 7*b*,*c*, have dimensions 1 ×1 m, both boxes are without side walls but with extra horizontal plates at the middle. The thickness of all plates in the second box is 2 mm. In the third box, only the front plate, which is in contact with the fluid during the impact, is 2 mm thick but other plates are of 2 cm. Figure 8*a*,*b* shows positions of the boxes (grey squares) on the wall (white squares) and the impact regions (striped areas). Figure 8*a* is for the first box, and the figure 8*b* for the second and third ones.

In all cases, the impact velocity was set to *U*=1 m s^{−1}, water depth *H*=2 m, channel width *L*=2 m and the sound speed in the fluid at rest *c*_{0} is equal to 1500 m s^{−1}. In the first case, the centre of the box is at (*y*,*z*)=(1 m,1.0725 m), and the wave height is *H*_{w}=1.2 m, and for the second and third cases, we have (*y*,*z*)=(1 m,1.25 m) and *H*_{w}=0.65 m.

It was revealed that for some conditions of wave impact with aeration, the strains in these structures may exceed those calculated for a homogeneous fluid without the aeration region. In figure 9, the stresses for the second box are shown in microstrains. In this case, the acoustic period of the fluid–gas region was a half of the period of the first mode of the vibration of the elastic box in air. The strains were evaluated at three points on the centre line of the wall: *z*=3/4 (figure 9*a*), *z*=1/2 (figure 9*b*) and *z*=1/4 (figure 9*c*) from the bottom of the box.

Finally, in figure 10, snapshots of the boxes’ deflections are shown, in order to demonstrate that our models predict the three-dimensional displacements of the structure elements throughout the whole structure.

## 7. Conclusions

The semi-analytical (or semi-numerical) approach for fluid–structure interactions which occur in different impact situations has been presented. Within this approach, the fluid flow is treated using the semi-analytical methods while the structural part is solved using the three-dimensional finite-element model. The choice of the simplified semi-analytical approach for the fluid flow was made in order to be able to have a full control of the flow characteristics, which allows for detailed investigations of the influence of different physical parameters. Indeed it seems to be impossible to correctly model all possible impact situations that might occur during the violent sloshing motion and the idea is to identify the most dangerous ones and check the sensitivity of the structural response to different physical quantities involved in the impact dynamics. In that respect, we identified the most important types of fluid impact on the CS and presented the detailed methodology for hydroelastic coupling. Accuracy and different convergence issues are properly investigated so that these models can be safely used in practice provided the validated structural finite-element model is available.

Finally, it is important to mention that, within the overall methodology for sloshing assessment, the above discussed models cannot be used on their own and should be supplemented either by the small-scale model tests or by the CFD calculations. Indeed, in the semi-analytical models, the physical parameters such as impact velocity, wave geometry, aeration and some others are supposed to be known. The determination of these parameters is far from trivial and they can be evaluated only approximately.

## Footnotes

One contribution of 13 to a Theme Issue ‘The mathematical challenges and modelling of hydroelasticity’.

- This journal is © 2011 The Royal Society