# Direct Time Domain Numerical Analysis of Transient Behavior of a VLFS during Unsteady External Loads in Wave Condition

Authors: Yong Cheng; Gang-jun Zhai; Jin-ping Ou
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

The transient response of the VLFS subjected to arbitrary external load is systematically investigated by a direct time domain modal expansion method, in which the BEM solutions based on time domain Kelvin sources are used for hydrodynamic forces. In the analysis, the time domain free-surface Green functions with sufficient accuracy are rapidly evaluated in finite water depth by the interpolation-tabulation method, and the boundary integral equation with a quarter VLFS model is established taking advantage of symmetry of flow field and structure. The validity of the present method is verified by comparing with the time histories of vertical displacements of the VLFS during a mass drop and airplane landing and takeoff in still water conditions, respectively. Then the developed numerical scheme is used in wave conditions to study the combined action taking into account the mass drop/airplane landing/takeoff loads as well as incident wave action. It is found that the elevation of structural waves due to mass drop load can be significantly changed near the impact region, while the vertical motion of runway in wave conditions is dominant as compared with that only generated by airplane.


---

## Body

## 1. Introduction

Very large floating structure (VLFS) is usually regarded as an alternative option of utilizing ocean space because of the scarcity of land to coastal regions, and it has been gradually appearing in various applications such as floating airports, oil storage vessels, floating artificial islands, and floating piers. For a VLFS, the structural length to the vertical length ratio as well as the structural length to wavelength ratio both is larger than unity (see [1]). Thus the deformation will become dominant over the rigid, and the fluid-structure interaction problem should be considered. In the hydroelastic analysis of a VLFS, there are two major methods. One way to tackle this problem is to use an analytical approach [2–6]. If the analytical approach is used, the computational time and memory capacity for a VLFS are not an issue. However, the solved drawback is only applied to simple geometries such as a rectangular plate or a circular plate. Another way to solve the hydroelastic problem of a VLFS is by using a numerical approach. The boundary element method (BEM) (see [7–9]), the finite-element method (FEM) (see [10, 11]), and the hybrid finite element-boundary element (FE-BE) method (see [12–14]) have been presented in previous studies.The element numbers of the wetted surfaces of the VLFS require a large memory capacity, and the time domain model contains the time parameter. Thus, dynamic response of the VLFS is commonly performed in the frequency domain (see [15–17]) when determining the hydroelastic response amplitude operator of the floating body and pertinent response parameters in a steady state condition. However, in real situation, nonharmonic external loads such as a huge mass impact on the structure and landing or taking off of an aircraft can induce the transient behavior of the VLFS and may affect the serviceability of the VLFS. Thus the transient response of the VLFS must be studied by a reliable calculation.Some numerical schemes for transient hydroelastic responses have been treated to date. Watanabe et al. [18] investigated a transient response analysis of a VLFS due to impulsive landing of an airplane by FEM. They applied the wave absorption filter to open boundaries; however, the response analysis was required for a few seconds. Kashiwagi [19, 20] developed an indirect time domain method in which the hydrodynamic effect is evaluated from good performance in the computation of the memory-effect function. Lee and Choi [21] proposed a FE-BE hybrid method to solve the transient responses indirectly by using transient equations, which are derived from the Fourier inverse transform of harmonic equations of motion and the causality condition. Endo [22] presented another time domain method based on FEM to treat the structure and on Wilson’s θ method to solve time-step procedure taking advantage of the memory-effect function for hydrodynamic forces. Shin et al. [23] simulated a transient behavior of a pontoon type VLFS subjected to airplane landing and takeoff, the calculation method for structural deflection is based on a FEM, and the fluid part is based on the BEM. Maeda et al. [24] analyzed the time-series responses without solving the equations of motion in the time domain. Kim and Webster [25] derived the structural motion using a two-dimensional analytical method and also solved the added drag to the aircraft.Though the above-mentioned studies provide enlightening contributions in the research activities related to external loads on the VLFS, some difficulties in carrying out their time domain simulation give restriction in the mathematical model. For example, the integration of memory-effect function is still time-consuming for evaluating at some high-frequency range, and the evaluation of hydrodynamic coefficients such as added mass at infinite frequency commonly neglect some cross-coupling terms in hydrodynamic effects. The author would like to develop a direct time integration method and the method that uses a superposition of modal functions with time-dependent unknown modal amplitudes and solves hydrodynamic diffraction and radiation problems by applying the time domain free-surface Green functions. In this direct time integration method, the present study simulate the transient hydroelastic response of a pontoon type VLFS under the combined action owing to external loads including a huge mass impact on the structure and landing or taking off of an aircraft as well as the incident wave. Numerical results are further addressed, with the time histories of vertical deflections at measured points, the spatial profiles of the VLFS at different times, and the running trajectory of the airplane. There is also a discussion of generated phenomena and the relationship of the vertical movement of the airplane and structural wave propagation.Fast and accurate calculation is necessary for overcoming this difficulty of large CPU time and memory size of computer. Utsunomiya et al. [26] and Teng and Gou [27] have developed the multipole expansion methods for hydroelastic analysis of a VLFS. Kagemoto et al. [28] presented the substructures method that accelerates computation without an appreciable loss of accuracy. Dai [29] has extended the precorrected-FFT method to hydroelastic analysis. However, their calculation models are only valid for the frequency domain studies for the wave-induced hydroelastic response. Huang [30] has put forth a feasible technique to tackle the time-free surface Green functions in infinite water depth; however, the VLFS commonly is placed in finite water depth. This paper derives the expression of the time domain free-surface Green functions and its spatial derivatives in finite water depth, which with sufficient accuracy are rapidly evaluated by using the interpolation-tabulation method. The low number of elements also is an important technique to reduce the memory and CPU time when the pressure distribution is obtained by BEM. According to the symmetry of the VLFS structure and the fluid field (see [31, 32]), this paper is only concerned with numerical simulations of boundary integral equations with a quarter VLFS model.


## 2. Formulation

We consider the time domain transient problem for a pontoon-type VLFS in finite water depth. Figure1 shows the fluid-structure problem and Cartesian coordinate system. The z-axis is pointing upwards, and the x-y plane is on the mean position of the free surface, where h is the water depth and A is the amplitude of the incident wave. The whole fluid domain is defined at Ω which contains the bottom of the VLFS Sb, side of the VLFS Ss, the free surface Sf, the seabed Sd, and the infinite cylindrical surface S∞. The VLFS has a length L, width B, and height hv, and d is the draft of the VLFS in z direction. The problem at hand is to determine the modal deflections under external loads combined action of incident waves.Figure 1
The fluid-structure system and coordinate system.Assuming that the fluid is incompressible, inviscid, and irrotational, a velocity potentialΦ(x,y,z,t) exists and would be given by
(1)Φ(x,y,z,t)=ΦI(x,y,z,t)+ΦS(x,y,z,t)inΩ,
where ΦI(x,y,z,t) and ΦS(x,y,z,t) are the incident and scattering potential, respectively.The velocity potential must satisfy the following Laplace’s equation and boundary conditions on the free surfaceSf, on the sea-bed Sd, the infinity S∞, and on the wetted surface of the floating body Sb (the bottom surface) and Ss (the side surface):
(2)∇2ΦS(x,y,z,t)=0,(3)∂2ΦS(x,y,z,t)∂t2+g∂ΦS(x,y,z,t)∂z=0on⁡Sf,t>0,(4)∂ΦI(x,y,z,t)∂n+∂ΦS(x,y,z,t)∂n=Vnon⁡Sb+Ss,t>0,(5)∂ΦS(x,y,z,t)∂n=0on⁡Sd,t>0,(6)ΦS(x,y,z,t),ΦSt(x,y,z,t),∇ΦS(x,y,z,t)⟶0on⁡S∞,t>0,(7)ΦS(x,y,z,t)=ΦSt(x,y,z,t)=0on⁡Sf,t=0,
where g is the gravitational acceleration, Vn is the normal velocity of the structure, n is a unit normal vector (the positive direction points out of the fluid domain), and ΦSt(x,y,z,t) represent the time derivative of the scattering potential.It is now widely accepted that VLFS response in terms of the vertical deflection can be captured well by modeling the whole VLFS as an elastic plate. In this formulation, assuming the VLFS as an elastic, isotropic, thin plate, the motion of the floating body is governed by the equation of a thin plate resting on a uniform elastic foundation:(8)D∇4W(x,y,z,t)=-msW¨+P(x,y,z,t)-ρgW-PE,
where D=EI is the bending rigidity, E is modulus of elasticity, I is the cross sectional moment of inertia, ms denotes the mass per unit area, ρ denotes density of fluid, PE denotes the external time-dependent loads acting on the VLFS due to a huge mass fall off or an airplane landing or take off, and the dynamic pressure P(x,y,-d,t) relates to the velocity potential on the bottom surface of the VLFS from the linearized Bernoulli’s equation
(9)P(x,y,-d,t)=-ρ∂Φ(x,y,-d,t)∂t.
In the present case, the VLFS is not constrained in the vertical elastic displacement along its edges; the following boundary conditions for a free edge must be satisfied:
(10)∂2W∂n2+v∂2W∂s2=0,∂3W∂n3+(2-v)∂2W∂n∂s2=0,
where v is Poisson’s ratio and n and s denote the normal and tangential directions, respectively.


## 3. Method of Solution

### 3.1. Modal Functions

The vertical elasticity displacementW(x,y,z,t) of the VLFS is the sum of various modes as follows:
(11)W(x,y,t)=∑j=1Nζj(t)fj(x,y)=∑m=0M∑c=0Nζmc(t)fxm(x)fyc(y),
where ζj(t) is the vibration amplitude of the j mode and fj(x,y) is modal function of the j mode. The modal function fj(x,y) of the VLFS can be expressed by the product of modal functions at x and y directions. Here, a case study of the natural modes fxm(x) and fyc(y) of a beam with free end is discussed as follows. The procedures mentioned above are mathematically valid for arbitrarily chosen modal function. However, in order to guarantee the convergence of the solutions, the appropriate choice of the modal functions is essential. In this paper, we employ the following modal functions (dry-modal functions) combining the rigid body motions and the elastic motions because they satisfy the free-end condition of the beam and their convergence has already been proved by Newman [33]:
(12)fx2m(x)={12form=012[cosh⁡(μmSx/(L/2))cosh⁡(μmS)+cos⁡(μmSx/(L/2))cos⁡(μmS)]form=1,2,…,fx2m+1(x)={32xL/2form=012[sinh⁡(μmAx/(L/2))sinh⁡(μmA)+sin(μmAx/(L/2))sin(μmA)]form=1,2,…,
where fx2m(x) and fx2m+1(x) are the symmetric and anti-symmetric modes about x=0, respectively, and μmS and μmA are the positive real roots of equations
(13)tanμmS+tanhμmS=0,tanμmS-tanhμmS=0
and fyc(y) may be written in the same form, with L/2 replaced B/2 by on the right-hand sides of (12).The modal functions expressed in (11) are orthogonal to each other in the (-L/2,L/2) with the following orthogonality relation:
(14)∬Sbfi(x,y)fj(x,y)dxdy=L4·B4δij,
where δij is Kroenecker’s delta, which is equal to 1 when i=j and 0 otherwise.


### 3.2. Fluid Part

The boundary value problems given by (2)–(7) can be solved by using Green’s functions. If free-surface Green’s function satisfying the boundary conditions given by (3), (5), (6), and (7) is considered, the boundary integral equation for the scattering potential can be derived as follows:
(15)αΦS(Q,t)+∬Sb(t)+Ss(t)ΦS(P,t)∂G0(P,Q)∂npdSp=∬Sb(t)+Ss(t)G0(P,Q)∂ΦS(P,t)∂npdSp+∫0tdτ∬Sb(τ)+Ss(τ)[∂ΦS(P,τ)∂npΦS(P,τ)×∂Gτf(P(τ),Q(t),t-τ)∂np-Gτf(P(τ),Q(t),t-τ)×∂ΦS(P,τ)∂np]dSp+1g∫0tdτ∬Cb(τ)[∂ΦS(P,τ)∂npΦS(P,τ)Gττf(P(τ),Q(t),t-τ)-Gτf(P(τ),Q(t),t-τ)×∂ΦS(P,τ)∂np]·Vn(P,τ)dSp,
where α represents the solid angle, Q(x0,y0,z0) and P(x,y,z) represent the source and field point, respectively, and Cb(t) represents the instantaneous waterline of the intersection between the body and the free surface. For Green's function G(P,t,Q,τ), it can be expressed by the superposition of instantaneous term G0 and memory term Gf in finite water depth h [34]:
(16)G(P,t,Q,τ)=G0(P,Q)+Gf(P,t,Q,τ)=1r+1r2-2∫0∞e-khcosh⁡k(z0+h)cosh⁡kh×cosh⁡k(z+h)J0(kR)dk+2∫0∞cosh⁡k(z+h)cosh⁡khsinh⁡khcosh⁡k(z0+h)×{1-cos⁡[(t-τ)gktanhkh]}J0(kR)dk,
where the instantaneous term G0 and memory term Gf are given in the form, respectively,
(17)G0=1r+1r2-2∫0∞e-khcosh⁡k(z0+h)cosh⁡kh×cosh⁡k(z+h)J0(kR)dk,(18)Gf=2∫0∞cosh⁡k(z+h)cosh⁡khsinh⁡khcosh⁡k(z0+h)×{1-cos⁡[(t-τ)gktanhkh]}J0(kR)dk,
where J0 is the Bessel function of the first kind, order zero, R denotes the horizontal distance between field and source point, r denotes the distance between field and source point, and r2 denotes the distance between field and the mirror image of the source field about water surface. In terms of the following nondimensional space and time parameters
(19)X=Rh,Y=-z0h,Z=-zh,T=(t-τ)(gh)1/2,
(17) may be written in the form
(20)G0=1h[F0(X,Y-Z)+F0(X,2-Y-Z)],
where the auxiliary function F0(X,V) is defined by
(21)F0(X,V)=1X2+V2-∫0∞e-usechucosh⁡⁡(uV)J0(uX)du.
The first order time derivative Gτf is expressed in the form
(22)Gτf=-g1/2h-3/2[F(X,Y-Z,T)+F(X,2-Y-Z,T)],
where
(23)F(X,V,T)=∫0∞ktanhkcosh⁡ksinh⁡ksin⁡(Tktanhk)×cosh⁡(kV)J0(kX)dk=2∫0∞ktanhksin⁡(Tktanhk)×ek(V-2)+e-k(V+2)1-e-4kJ0(kX)dk.Thus the rate of convergence of (21) depends primarily on X. Wehausen and Laitone [34] developed the following expression:
(24)F0(X,V)={2∑m=0∞cos⁡[(m+12)πV]K0[(m+12)πX]X>11X2+V2-1X2+(V+2)2-1X2+(V-2)2+1-log⁡2-∑m,namnX2mV2n0≤X≤1,
where K0 is the modified Bessel function of the second kind and the coefficients amn are defined by Newman [35].The rate of convergence of the integral in (23) can be accelerated by adding and subtracting the appropriate function F∞ which may be given in the form
(25)F∞(X,V,T)=lim⁡k→∞F(X,V,T)=2∫0∞ksin(Tk)ek(V-2)J0(kX)dk,
where the vertical coordinate V is restricted to the fluid domain (-1,2).If the spherical coordinate is adopted, the nondimensional parameters onr are defined by
(26)r=[X2+(V-2)2]1/2,τ=Tr1/2,k=kr,V-2=rcos⁡θ,X=rsinθ,
where τ and θ lie in the interval (0,∞) and (0,2π). Thus, we will reduce three arguments to two arguments by substituting (26) into (25). The function F∞ and its partial derivatives are obtained:
(27)F∞=2r-3/2Im⁡{∫0∞ke-kcos⁡θeiτkJ0(ksinθ)dk},F∞X=-2r-5/2Im⁡{∫0∞k32e-kcos⁡θeiτkJ1(ksinθ)dk},F∞V=2r-5/2Im⁡{∫0∞k3/2e-kcos⁡θeiτkJ0(ksinθ)dk},F∞T=2r-2Re{∫0∞ke-kcos⁡θeiτkJ0(ksinθ)dk}.Since the integrands in (27) exhibit slow convergence and high oscillatory, the F∞ and its spatial derivatives can be approximated in terms of the values of parameter τ (see [30, 35]). The function F-F∞ may be solved directly in a straightforward manner due to the oscillatory elimination, and thus F(X,V,T) is obtained:
(28)F(X,V,T)=F∞+F-F∞.Then the boundary surface of (15) is discretized into a number of elements using a standard procedure known as the BEM. Within the boundary elements, physical variables are interpolated by the shape functions, which represent the geometry of each element. In the integration process, the scheme using trapezoidal approximation is applied to the convolution integral. Once (15) is solved, the time history of dynamic pressure (9) can be obtained at any position.


### 3.3. Structure Part

We substitute the fluid pressure and vertical deflection (9) and (11) into (8); we have
(29)∑m=1M∑c=1Nmsfxm(x)fyc(y)ζ¨mc(t)+(D∇4fxm(x)fyc(y)+ρgfxm(x)fyc(y))ζmc(t)=-ρ∂(ΦI+ΦS)∂t-PE,
where m, c are the modal numbers in x and y direction, respectively. Applying Galerkin’s method, we multiply both sides of the above equation by fxi(x)·fyj(y) and integrate over the bottom of the VLFS. Finally, we can obtain a conventional set of equations given by
(30)∑m=1M∑c=1NMmc,ijζ¨mc(t)+∑m=1M∑c=1NKmc,ijζmc(t)=Fij(t)+Eij(t),i=1,…,M,j=1,…,N,
where
(31)Mmc,ij=∬Sbmsfxm(x)fyc(y)fxi(x)fyj(y)ds,(32)Kmc,ij=D∬Sb∇2fxm(x)fyc(y)∇2fxi(x)fyj(y)ds-D(1-v)×∬Sb{∂2fxm(x)fyc(y)∂x2∂2fxi(x)fyj(y)∂y2+∂2fxm(x)fyc(y)∂y2∂2fxi(x)fyj(y)∂x2-2∂2fxm(x)fyc(y)∂x∂y∂2fxi(x)fyj(y)∂x∂y}ds+∬Sbρgfxm(x)fyc(y)fxi(x)fyj(y)ds,(33)Fij=∬Sb-ρ∂(ΦI+ΦS)∂tfxi(x)fyj(y)ds,(34)Eij(t)=-∬SbPE(x,y,t)fxi(x)fyj(y)ds.
It should be noted that Mmc,ij, Kmc,ij, Fij, and Eij are the generalized mass, generalized stiffness, generalized wave force, and generalized external load force, respectively. And the generalized stiffness Kmc,ij shown by (32) has been obtained by taking account of the free-edge boundary conditions, (10), referring to the paper by Kashiwagi [19].In order to solve the deflection of the VLFS, (30) is solved by using the fourth order Runge-Kutta method.


## 3.1. Modal Functions

The vertical elasticity displacementW(x,y,z,t) of the VLFS is the sum of various modes as follows:
(11)W(x,y,t)=∑j=1Nζj(t)fj(x,y)=∑m=0M∑c=0Nζmc(t)fxm(x)fyc(y),
where ζj(t) is the vibration amplitude of the j mode and fj(x,y) is modal function of the j mode. The modal function fj(x,y) of the VLFS can be expressed by the product of modal functions at x and y directions. Here, a case study of the natural modes fxm(x) and fyc(y) of a beam with free end is discussed as follows. The procedures mentioned above are mathematically valid for arbitrarily chosen modal function. However, in order to guarantee the convergence of the solutions, the appropriate choice of the modal functions is essential. In this paper, we employ the following modal functions (dry-modal functions) combining the rigid body motions and the elastic motions because they satisfy the free-end condition of the beam and their convergence has already been proved by Newman [33]:
(12)fx2m(x)={12form=012[cosh⁡(μmSx/(L/2))cosh⁡(μmS)+cos⁡(μmSx/(L/2))cos⁡(μmS)]form=1,2,…,fx2m+1(x)={32xL/2form=012[sinh⁡(μmAx/(L/2))sinh⁡(μmA)+sin(μmAx/(L/2))sin(μmA)]form=1,2,…,
where fx2m(x) and fx2m+1(x) are the symmetric and anti-symmetric modes about x=0, respectively, and μmS and μmA are the positive real roots of equations
(13)tanμmS+tanhμmS=0,tanμmS-tanhμmS=0
and fyc(y) may be written in the same form, with L/2 replaced B/2 by on the right-hand sides of (12).The modal functions expressed in (11) are orthogonal to each other in the (-L/2,L/2) with the following orthogonality relation:
(14)∬Sbfi(x,y)fj(x,y)dxdy=L4·B4δij,
where δij is Kroenecker’s delta, which is equal to 1 when i=j and 0 otherwise.


## 3.2. Fluid Part

The boundary value problems given by (2)–(7) can be solved by using Green’s functions. If free-surface Green’s function satisfying the boundary conditions given by (3), (5), (6), and (7) is considered, the boundary integral equation for the scattering potential can be derived as follows:
(15)αΦS(Q,t)+∬Sb(t)+Ss(t)ΦS(P,t)∂G0(P,Q)∂npdSp=∬Sb(t)+Ss(t)G0(P,Q)∂ΦS(P,t)∂npdSp+∫0tdτ∬Sb(τ)+Ss(τ)[∂ΦS(P,τ)∂npΦS(P,τ)×∂Gτf(P(τ),Q(t),t-τ)∂np-Gτf(P(τ),Q(t),t-τ)×∂ΦS(P,τ)∂np]dSp+1g∫0tdτ∬Cb(τ)[∂ΦS(P,τ)∂npΦS(P,τ)Gττf(P(τ),Q(t),t-τ)-Gτf(P(τ),Q(t),t-τ)×∂ΦS(P,τ)∂np]·Vn(P,τ)dSp,
where α represents the solid angle, Q(x0,y0,z0) and P(x,y,z) represent the source and field point, respectively, and Cb(t) represents the instantaneous waterline of the intersection between the body and the free surface. For Green's function G(P,t,Q,τ), it can be expressed by the superposition of instantaneous term G0 and memory term Gf in finite water depth h [34]:
(16)G(P,t,Q,τ)=G0(P,Q)+Gf(P,t,Q,τ)=1r+1r2-2∫0∞e-khcosh⁡k(z0+h)cosh⁡kh×cosh⁡k(z+h)J0(kR)dk+2∫0∞cosh⁡k(z+h)cosh⁡khsinh⁡khcosh⁡k(z0+h)×{1-cos⁡[(t-τ)gktanhkh]}J0(kR)dk,
where the instantaneous term G0 and memory term Gf are given in the form, respectively,
(17)G0=1r+1r2-2∫0∞e-khcosh⁡k(z0+h)cosh⁡kh×cosh⁡k(z+h)J0(kR)dk,(18)Gf=2∫0∞cosh⁡k(z+h)cosh⁡khsinh⁡khcosh⁡k(z0+h)×{1-cos⁡[(t-τ)gktanhkh]}J0(kR)dk,
where J0 is the Bessel function of the first kind, order zero, R denotes the horizontal distance between field and source point, r denotes the distance between field and source point, and r2 denotes the distance between field and the mirror image of the source field about water surface. In terms of the following nondimensional space and time parameters
(19)X=Rh,Y=-z0h,Z=-zh,T=(t-τ)(gh)1/2,
(17) may be written in the form
(20)G0=1h[F0(X,Y-Z)+F0(X,2-Y-Z)],
where the auxiliary function F0(X,V) is defined by
(21)F0(X,V)=1X2+V2-∫0∞e-usechucosh⁡⁡(uV)J0(uX)du.
The first order time derivative Gτf is expressed in the form
(22)Gτf=-g1/2h-3/2[F(X,Y-Z,T)+F(X,2-Y-Z,T)],
where
(23)F(X,V,T)=∫0∞ktanhkcosh⁡ksinh⁡ksin⁡(Tktanhk)×cosh⁡(kV)J0(kX)dk=2∫0∞ktanhksin⁡(Tktanhk)×ek(V-2)+e-k(V+2)1-e-4kJ0(kX)dk.Thus the rate of convergence of (21) depends primarily on X. Wehausen and Laitone [34] developed the following expression:
(24)F0(X,V)={2∑m=0∞cos⁡[(m+12)πV]K0[(m+12)πX]X>11X2+V2-1X2+(V+2)2-1X2+(V-2)2+1-log⁡2-∑m,namnX2mV2n0≤X≤1,
where K0 is the modified Bessel function of the second kind and the coefficients amn are defined by Newman [35].The rate of convergence of the integral in (23) can be accelerated by adding and subtracting the appropriate function F∞ which may be given in the form
(25)F∞(X,V,T)=lim⁡k→∞F(X,V,T)=2∫0∞ksin(Tk)ek(V-2)J0(kX)dk,
where the vertical coordinate V is restricted to the fluid domain (-1,2).If the spherical coordinate is adopted, the nondimensional parameters onr are defined by
(26)r=[X2+(V-2)2]1/2,τ=Tr1/2,k=kr,V-2=rcos⁡θ,X=rsinθ,
where τ and θ lie in the interval (0,∞) and (0,2π). Thus, we will reduce three arguments to two arguments by substituting (26) into (25). The function F∞ and its partial derivatives are obtained:
(27)F∞=2r-3/2Im⁡{∫0∞ke-kcos⁡θeiτkJ0(ksinθ)dk},F∞X=-2r-5/2Im⁡{∫0∞k32e-kcos⁡θeiτkJ1(ksinθ)dk},F∞V=2r-5/2Im⁡{∫0∞k3/2e-kcos⁡θeiτkJ0(ksinθ)dk},F∞T=2r-2Re{∫0∞ke-kcos⁡θeiτkJ0(ksinθ)dk}.Since the integrands in (27) exhibit slow convergence and high oscillatory, the F∞ and its spatial derivatives can be approximated in terms of the values of parameter τ (see [30, 35]). The function F-F∞ may be solved directly in a straightforward manner due to the oscillatory elimination, and thus F(X,V,T) is obtained:
(28)F(X,V,T)=F∞+F-F∞.Then the boundary surface of (15) is discretized into a number of elements using a standard procedure known as the BEM. Within the boundary elements, physical variables are interpolated by the shape functions, which represent the geometry of each element. In the integration process, the scheme using trapezoidal approximation is applied to the convolution integral. Once (15) is solved, the time history of dynamic pressure (9) can be obtained at any position.


## 3.3. Structure Part

We substitute the fluid pressure and vertical deflection (9) and (11) into (8); we have
(29)∑m=1M∑c=1Nmsfxm(x)fyc(y)ζ¨mc(t)+(D∇4fxm(x)fyc(y)+ρgfxm(x)fyc(y))ζmc(t)=-ρ∂(ΦI+ΦS)∂t-PE,
where m, c are the modal numbers in x and y direction, respectively. Applying Galerkin’s method, we multiply both sides of the above equation by fxi(x)·fyj(y) and integrate over the bottom of the VLFS. Finally, we can obtain a conventional set of equations given by
(30)∑m=1M∑c=1NMmc,ijζ¨mc(t)+∑m=1M∑c=1NKmc,ijζmc(t)=Fij(t)+Eij(t),i=1,…,M,j=1,…,N,
where
(31)Mmc,ij=∬Sbmsfxm(x)fyc(y)fxi(x)fyj(y)ds,(32)Kmc,ij=D∬Sb∇2fxm(x)fyc(y)∇2fxi(x)fyj(y)ds-D(1-v)×∬Sb{∂2fxm(x)fyc(y)∂x2∂2fxi(x)fyj(y)∂y2+∂2fxm(x)fyc(y)∂y2∂2fxi(x)fyj(y)∂x2-2∂2fxm(x)fyc(y)∂x∂y∂2fxi(x)fyj(y)∂x∂y}ds+∬Sbρgfxm(x)fyc(y)fxi(x)fyj(y)ds,(33)Fij=∬Sb-ρ∂(ΦI+ΦS)∂tfxi(x)fyj(y)ds,(34)Eij(t)=-∬SbPE(x,y,t)fxi(x)fyj(y)ds.
It should be noted that Mmc,ij, Kmc,ij, Fij, and Eij are the generalized mass, generalized stiffness, generalized wave force, and generalized external load force, respectively. And the generalized stiffness Kmc,ij shown by (32) has been obtained by taking account of the free-edge boundary conditions, (10), referring to the paper by Kashiwagi [19].In order to solve the deflection of the VLFS, (30) is solved by using the fourth order Runge-Kutta method.


## 4. External Loads

The motion equation (30) can be applied to any time domain transient problem, where the generalized external load force Eij(t) is solved. Here, a huge mass fall off and an airplane landing and takeoff will be considered, for which the external pressure distribution PE in (34) must be discussed as follows.


### 4.1. The Weight Drop Test

In a weight drop test, a weightW was dropped from a height onto the “hit point.” The acceleration of the weight during the impact was a. Therefore, the impact load Fim(t) can be obtained:
(35)Fim(t)=W·a,
and the external pressure distribution PE, appearing in (34) can be expressed as
(36)PE(x,y,t)=Fim(t)·δ(x-xp)·δ(y-yp),
where (xp,yp) is the coordinate of the hit point.Substituting (36) into (34), the generalized external force Eij(t) can be computed as
(37)Eij(t)=-∬SbFim(t)δ(x-xp)δ(y-yp)fxi(x)fyj(y)ds=-Fim(t)fxi(xp)fyj(yp).The numerical parameters for simulation in this paper are given in Table1 by referring to Endo et al. [36], and the vertical displacements at points Z1–Z9 shown in Figure 2 were measured.Table 1
Main parameters of weight drop test.
Total length (m)
Total width (m)
Thickness (m)
Bending rigidity (Nm)
Poisson’s ratio
Weight  (N)
Dropping height (m)
9.75
1.95
0.0163
8985.62
0.3
196
0.12Figure 2
The position of the measured points in the drop test.


### 4.2. The Airplane Landing and Takeoff

A realistic situation is simulated where an airplane landing or takeoff on a VLFS. Here, the time-varying load is assumed to move with a constant initial accelerationα0, and the position ξ(t) of the airplane and its velocity V(t) are given by
(38)ξ(t)=ξ0+V0t+12α0t2,V(t)=V0+α0t,
where ξ0 and V0 are the initial position and velocity, respectively. For simplicity, the load distribution is assumed to be axisymmetric about the center of the moving load (ξ(t),0). In the terms of the relationship between the moving Cartesian coordinate system o--x-yz and the polar coordinate system o--r-θz, the external pressure distribution PE can be expressed as
(39)PE(x,y,t)=Fai(t)·f(r-)=Fai(t)·1RE2e-π(r-/RE)2,
where r-=x-2+y2,  x-=x-ξ(t) and RE denotes the effective radius of the loading.The total forceFai(t) exerted by the landing or takeoff on the VLFS can be given by the difference between the weight W of the airplane and the lift force FL(t):
(40)Fai(t)=W-FL(t),(41)FL(t)=12ρaV2(t)AWaLebLt,
where the parameters aL and bL are given as constants, ρa is the density of air, and AW is the effective wing area of the airplane.Substituting (39) into (34), the generalized external force Eij(t) can be computed as
(42)Eij(t)=-∬SbFai(t)RE2e-π(r-/RE)2fxi(x)fyj(y)ds=-Fai(t)RE2∫-3RE3REe-π(x-/RE)2fxi(x-+ξ(t))ds×∫-3RE3REe-π(y/RE)2fyj(y)ds.The numerical data for simulation in this paper is prepared as listed in Table2 by referring to Kashiwagi [20]. The touch-down position in landing and leave-up position in takeoff are shown in Figure 3, together with measured points (Z1–Z9) for the vertical displacements. It is assumed to land or take off in the following wave direction from the fore-end of the runway.Table 2
Main parameters of landing/takeoff run.
Floating airport
Airplane
Length (L)
5000 m
Weight (W)
3867.08 kN
Width (B)
1000 m
Effective wing area (AW)
511.0 m2
Draft (d)
5.0 m
Effective radius (RE)
10.0 m
Bending rigidity (EI/B)
1.764 × 1011 Nm
Initial position (ξ0)
−1000 m
Initial speed (V0)
Landing
69.35 ms−1
Takeoff
0.0 ms−1
Acceleration (α0)
Landing
−1.263 ms−2
Takeoff
1.026 ms−2
a
L in landing
2.61
a
L in takeoff
1.64 × 10−3
b
L in landing
−0.212
b
L in takeoff
0.125Figure 3
The plan view of VLFS during landing or takeoff.


## 4.1. The Weight Drop Test

In a weight drop test, a weightW was dropped from a height onto the “hit point.” The acceleration of the weight during the impact was a. Therefore, the impact load Fim(t) can be obtained:
(35)Fim(t)=W·a,
and the external pressure distribution PE, appearing in (34) can be expressed as
(36)PE(x,y,t)=Fim(t)·δ(x-xp)·δ(y-yp),
where (xp,yp) is the coordinate of the hit point.Substituting (36) into (34), the generalized external force Eij(t) can be computed as
(37)Eij(t)=-∬SbFim(t)δ(x-xp)δ(y-yp)fxi(x)fyj(y)ds=-Fim(t)fxi(xp)fyj(yp).The numerical parameters for simulation in this paper are given in Table1 by referring to Endo et al. [36], and the vertical displacements at points Z1–Z9 shown in Figure 2 were measured.Table 1
Main parameters of weight drop test.
Total length (m)
Total width (m)
Thickness (m)
Bending rigidity (Nm)
Poisson’s ratio
Weight  (N)
Dropping height (m)
9.75
1.95
0.0163
8985.62
0.3
196
0.12Figure 2
The position of the measured points in the drop test.


## 4.2. The Airplane Landing and Takeoff

A realistic situation is simulated where an airplane landing or takeoff on a VLFS. Here, the time-varying load is assumed to move with a constant initial accelerationα0, and the position ξ(t) of the airplane and its velocity V(t) are given by
(38)ξ(t)=ξ0+V0t+12α0t2,V(t)=V0+α0t,
where ξ0 and V0 are the initial position and velocity, respectively. For simplicity, the load distribution is assumed to be axisymmetric about the center of the moving load (ξ(t),0). In the terms of the relationship between the moving Cartesian coordinate system o--x-yz and the polar coordinate system o--r-θz, the external pressure distribution PE can be expressed as
(39)PE(x,y,t)=Fai(t)·f(r-)=Fai(t)·1RE2e-π(r-/RE)2,
where r-=x-2+y2,  x-=x-ξ(t) and RE denotes the effective radius of the loading.The total forceFai(t) exerted by the landing or takeoff on the VLFS can be given by the difference between the weight W of the airplane and the lift force FL(t):
(40)Fai(t)=W-FL(t),(41)FL(t)=12ρaV2(t)AWaLebLt,
where the parameters aL and bL are given as constants, ρa is the density of air, and AW is the effective wing area of the airplane.Substituting (39) into (34), the generalized external force Eij(t) can be computed as
(42)Eij(t)=-∬SbFai(t)RE2e-π(r-/RE)2fxi(x)fyj(y)ds=-Fai(t)RE2∫-3RE3REe-π(x-/RE)2fxi(x-+ξ(t))ds×∫-3RE3REe-π(y/RE)2fyj(y)ds.The numerical data for simulation in this paper is prepared as listed in Table2 by referring to Kashiwagi [20]. The touch-down position in landing and leave-up position in takeoff are shown in Figure 3, together with measured points (Z1–Z9) for the vertical displacements. It is assumed to land or take off in the following wave direction from the fore-end of the runway.Table 2
Main parameters of landing/takeoff run.
Floating airport
Airplane
Length (L)
5000 m
Weight (W)
3867.08 kN
Width (B)
1000 m
Effective wing area (AW)
511.0 m2
Draft (d)
5.0 m
Effective radius (RE)
10.0 m
Bending rigidity (EI/B)
1.764 × 1011 Nm
Initial position (ξ0)
−1000 m
Initial speed (V0)
Landing
69.35 ms−1
Takeoff
0.0 ms−1
Acceleration (α0)
Landing
−1.263 ms−2
Takeoff
1.026 ms−2
a
L in landing
2.61
a
L in takeoff
1.64 × 10−3
b
L in landing
−0.212
b
L in takeoff
0.125Figure 3
The plan view of VLFS during landing or takeoff.


## 5. Fast Algorithm

### 5.1. Interpolation-Tabulation Method

Accurate and fact computation of the Green function and its derivations is important for saving the CPU time and memory of the computer. The interpolation-tabulation method is applied to the solutions ofF∞ and F-F∞ in (28) as follows.The three variablesX, V, T of the function F∞ are changed to two arguments τ and cos⁡θ, which are divided into 800 and 200 parts, respectively. The solutions for cos⁡θ<0.7 which are efficient in the context of Section 3.2 are described by Beck and Liapis [37]; else the Filon integral scheme is determined to calculate directly. During solving (15), the bilinear interpolation scheme was applied to the effective approximation of F∞ and its derivation.However, theF-F∞ function has three arguments V, X, and T. Here, the space nondimensional parameters V and X are restricted to the region (-1,2) and (0,20), and the time nondimensional parameters T lies in the interval (0,20). First seven X-T planes are adopted in V direction. Next, every X-T plane is divided into 40 parts in the X and T direction, respectively. The slowly varying function F-F∞ and its derivations can be calculated by using Gaussian integration; then the trilinear interpolation scheme is applied to the effective approximation in (15).


### 5.2. Symmetry of Structure

With the discretization of the constant boundary elements, (15) may be expressed as a form of the linear equations
(43)[A]{ΦS}={B}.Considering symmetry of the VLFS aboutx-z plane and y-z plane shown in Figure 4, the matrix [A], vector {ΦS} and {B}, may be divided as follows:
(44)[A]=[A11A12A13A14A21A22A23A24A31A32A33A34A41A42A43A44],{ΦS1ΦS2ΦS3ΦS4},{B1B2B3B4}.
The symmetric relationships for the matrix [A] may be formulated as
(45)A11=A22=A33=A44,A12=A21=A34=A43,A13=A31=A24=A42,A14=A41=A23=A32.
In order to reduce the dimensions of the matrix, the conversions is obtained by taking
(46){Φ^S}=[E]{ΦS},{B^}=[E]{B},[I]=1β[E]2,
where the constant coefficient β=4 for two planes of symmetry; the transition matrix [E] is given by
(47)[E]=[11111-11-111-1-11-1-11].Figure 4
Sketch of the district of symmetry.Thus (43) can be simplified by substituting (46)
(48)[A^]{Φ^S}={B^},
where [A^]=(1/4)[E][A][E] and the matrix [A^] is a block-diagonal matrix whose only nonzero elements occur in the square blocks centered about the principal diagonal, according to the features of the matrices [E] and [A]. So the linear equations in region 1 will be
(49)[A^(1)]{Φ^S(1)}={B^(1)},[A^(1)]=[A11+A12+A13+A14],{Φ^S(1)}={ΦS1+ΦS2+ΦS3+ΦS4},{B^(1)}={B1+B2+B3+B4}.
In region 2, we have
(50)[A^(2)]{Φ^S(2)}={B^(2)},[A^(2)]=[A11-A12+A13-A14],{Φ^S(2)}={ΦS1-ΦS2+ΦS3-ΦS4},{B^(2)}={B1-B2+B3-B4}.
In region 3, we have
(51)[A^(3)]{Φ^S(3)}={B^(3)},[A^(3)]=[A11+A12-A13-A14],{Φ^S(3)}={ΦS1+ΦS2-ΦS3-ΦS4},{B^(3)}={B1+B2-B3-B4}.
In region 4, we have
(52)[A^(4)]{Φ^S(4)}={B^(4)},[A^(4)]=[A11-A12-A13+A14],{Φ^S(4)}={ΦS1-ΦS2-ΦS3+ΦS4},{B^(4)}={B1-B2-B3+B4}.


## 5.1. Interpolation-Tabulation Method

Accurate and fact computation of the Green function and its derivations is important for saving the CPU time and memory of the computer. The interpolation-tabulation method is applied to the solutions ofF∞ and F-F∞ in (28) as follows.The three variablesX, V, T of the function F∞ are changed to two arguments τ and cos⁡θ, which are divided into 800 and 200 parts, respectively. The solutions for cos⁡θ<0.7 which are efficient in the context of Section 3.2 are described by Beck and Liapis [37]; else the Filon integral scheme is determined to calculate directly. During solving (15), the bilinear interpolation scheme was applied to the effective approximation of F∞ and its derivation.However, theF-F∞ function has three arguments V, X, and T. Here, the space nondimensional parameters V and X are restricted to the region (-1,2) and (0,20), and the time nondimensional parameters T lies in the interval (0,20). First seven X-T planes are adopted in V direction. Next, every X-T plane is divided into 40 parts in the X and T direction, respectively. The slowly varying function F-F∞ and its derivations can be calculated by using Gaussian integration; then the trilinear interpolation scheme is applied to the effective approximation in (15).


## 5.2. Symmetry of Structure

With the discretization of the constant boundary elements, (15) may be expressed as a form of the linear equations
(43)[A]{ΦS}={B}.Considering symmetry of the VLFS aboutx-z plane and y-z plane shown in Figure 4, the matrix [A], vector {ΦS} and {B}, may be divided as follows:
(44)[A]=[A11A12A13A14A21A22A23A24A31A32A33A34A41A42A43A44],{ΦS1ΦS2ΦS3ΦS4},{B1B2B3B4}.
The symmetric relationships for the matrix [A] may be formulated as
(45)A11=A22=A33=A44,A12=A21=A34=A43,A13=A31=A24=A42,A14=A41=A23=A32.
In order to reduce the dimensions of the matrix, the conversions is obtained by taking
(46){Φ^S}=[E]{ΦS},{B^}=[E]{B},[I]=1β[E]2,
where the constant coefficient β=4 for two planes of symmetry; the transition matrix [E] is given by
(47)[E]=[11111-11-111-1-11-1-11].Figure 4
Sketch of the district of symmetry.Thus (43) can be simplified by substituting (46)
(48)[A^]{Φ^S}={B^},
where [A^]=(1/4)[E][A][E] and the matrix [A^] is a block-diagonal matrix whose only nonzero elements occur in the square blocks centered about the principal diagonal, according to the features of the matrices [E] and [A]. So the linear equations in region 1 will be
(49)[A^(1)]{Φ^S(1)}={B^(1)},[A^(1)]=[A11+A12+A13+A14],{Φ^S(1)}={ΦS1+ΦS2+ΦS3+ΦS4},{B^(1)}={B1+B2+B3+B4}.
In region 2, we have
(50)[A^(2)]{Φ^S(2)}={B^(2)},[A^(2)]=[A11-A12+A13-A14],{Φ^S(2)}={ΦS1-ΦS2+ΦS3-ΦS4},{B^(2)}={B1-B2+B3-B4}.
In region 3, we have
(51)[A^(3)]{Φ^S(3)}={B^(3)},[A^(3)]=[A11+A12-A13-A14],{Φ^S(3)}={ΦS1+ΦS2-ΦS3-ΦS4},{B^(3)}={B1+B2-B3-B4}.
In region 4, we have
(52)[A^(4)]{Φ^S(4)}={B^(4)},[A^(4)]=[A11-A12-A13+A14],{Φ^S(4)}={ΦS1-ΦS2-ΦS3+ΦS4},{B^(4)}={B1-B2-B3+B4}.


## 6. Results and Discussion

### 6.1. Accuracy in the Interpolation-Tabulation Method

Before starting numerical simulations, it is necessary to confirm good performance in the computation of the time domain free-surface Green functions and its spatial derivatives in finite water depth.In order to examine the validity of the interpolation-tabulation method, computations of the nondimensional functionsF∞, F-F∞, and their spatial derivatives are performed for T=10, V=2 and are compared with corresponding results obtained from Newman. The values are shown in Figures 5 and 6. Obviously, the interpolation-tabulation method can give reliable evaluations which are smooth and in good agreement with Newman’s values.Comparison ofF∞, F∞X, and F∞V by the interpolation-tabulation method with Newman.
(a)
(b)
(c)Comparison ofF-F∞,  (F-F∞)X, and (F-F∞)V by the interpolation-tabulation method with Newman.
(a)
(b)
(c)


### 6.2. Drop Test in Still Water

The numerical simulation of the weight drop test is implemented, corresponding to the experiments conducted by Endo et al. [36] and the numerical results solved by Kashiwagi [19]. The pertinent information for the test model is prepared as listed in Section 4.1.Good convergence is considered for the numbers of modes in thex-direction and y-direction, after referring to the results of Kashiwagi [19], the number of modes in the x-direction M=8 and in the y-direction N=3 is adopted in this case.The comparisons of the vertical deflection time series at measured points are indicated in Figure2 among the present results; the indirect time domain solutions solved by Kashiwagi [19] and experimental tests obtained by Endo et al. [36] are given in Figure 7. It can be seen from this figure that the degree of agreement for these methods is favorable. And the deflections by the present method near the impact point, such as Z1 and Z2, are closer to the measurements than the numerical results solved by Kashiwagi [19]. This may be attributed to the difference in fluid pressure computation between the direct domain method by considering free-surface Green function and the indirect domain method by using the convolution integral of frequency impulse function.Comparison of time histories of the vertical displacement at each measured point.
(a)
(b)The deformed profiles of the VLFS during the mass drop are shown in Figure8. It is seen that the structural wave is transmitted at the longitudinal centerline of the plate, and the shape of the deformation is close to the current static equilibrium configuration at t=1.85 s. The magnitude of the vertical displacements is less than 1.0 cm. The vertical displacements of the plate are very small when x coordinate value is less than 0 but the transient phenomena at the right edge of the VLFS can be obviously seen at times t=0.21 s to 0.80 s.Figure 8
Spatial profiles of the VLFS due to the mass drop.


### 6.3. Drop Test in Regular Wave

In the simulation of the time step procedure, the load force at the right hand side of (30) is divided into two stages. The regular wave comes first from the left side of the VLFS and then the weight drops later by three cycles of the wave period. The wave length is 1.0 m, wave period is 0.8 seconds, wave height is 1 cm, and incident angle is 0 degrees.The deformed profiles of the VLFS during the weight drop are shown in Figure9, where Figure 9(a) shows the deflections in the regular wave condition without the mass impact for t=0 s and the deflections in the early stage of the drop test for t=0.21 and 0.41 s. The figure tells us that the absolute values of vertical displacement at the fore-end of the VLFS in regular wave are about 10 times the ones only generated by the mass drop; however, the magnitudes at the back-end are almost equivalent to the results induced by the mass drop. This means the mass impact load should not be overlooked as compared with the wave load. In other words, the structural wave shape of the VLFS in wave condition is changed when a huge mass falls off on the platform of the VLFS.Spatial profiles of the VLFS during the mass drop in regular wave conditions.
(a)
(b)


### 6.4. Landing in the Still Water

As mentioned in Section4.2, the airplane lands at point Z3 and completes the landing run in 54.9 s. The time histories of the vertical displacements at measurement points Z1, Z4, Z5, Z7, Z8, and Z9 obtained from the present method (direct time domain method) and the indirect time domain method used by Kashiwagi [20] are comparatively shown in Figure 10. The correlation between the two numerical solutions is reasonable. It can be seen that the magnitude of the deflections at measured points is less than about 1.0 cm and is very small as compared with the length value of the runway though the airplane weight is approximately up to 3900 kN. Then the time histories of the vertical displacement are much smoother than the results during weight drop (see Figure 7) and no higher-order deflections are found in time histories of the deflection for the landing on the platform of the VLFS. It is primarily due to the smooth increase of the landing loads (see (40) and (41)). It is also interesting that the vertical displacement of measured points Z7, Z8, and Z9 is not so large at t=55 s to 60 s but increases again after t=60 s. This can be attributed to the radiation of structure waves which impinge the stopped airplane.Time histories of the vertical displacement subjected to landing of airplane.
(a)
(b)Figure11 shows the snapshots of the deflection along the longitudinal centerline of the runway at different times, and the corresponding positions of the airplane are expressed by circles. It is found that the structural waves run after the airplane at time less than 42 s, and then the waves catch up to the airplane at about time 53 s because of the decrease of the airplane speed. After overtaking, structural waves meet the stopped airplane, partial waves are diffracted, and remainder is transmitted. Thus, it is interesting to find that the deformed profiles of the runway at time t=66 s are larger than those of the time t=53 s (see Figure 10). Note that the airplane seems to stay always at the bottom of sunken deflections of the runway during the landing run.Snapshots of the deflection subjected to landing of airplane.
(a)
(b)


### 6.5. Landing in Regular Wave

The interaction of incident wave with the runway during landing is divided into two stages. The regular wave comes first from the fore-end of the VLFS; after three cycles of the wave period, the airplane touches down the runway in the following incident wave direction. The wavelength is 650 m (0.13 times of the runway length), period 23.6 s, height 2.0 m, and incident angle is 0 degrees.The spatial profiles of the runway during the running are shown together with the positions of the airplane in Figure12. It can be seen that the airplane runs faster than the structural waves induced by the regular incident wave in the early stage (at least up 30 s) and then the structure waves overtake the airplane at the final stage of the run (after t=30 s). The maximum vertical displacement in regular wave is 150 cm and is about 150 times the one induced in the still water condition. This means that the wave load is dominant compared with the landing load in the hydroelastic analysis of VLFS.Spatial profiles of the deflection subjected to landing of airplane in regular wave conditions.
(a)
(b)Looking at the history of the vertical displacement of locations in Figure13(a) and the corresponding path during landing in Figure 13(b), it can be seen that the propagating velocity of the structural wave generated by incident wave is slower than the landing speed of airplane in the early stage (at least up to 30 s); however, when the airplane slows down, the deflections of the runway change suddenly in their magnitude and length (20 s–40 s as shown in Figure 13(a)). At the final stage of landing, speed of the airplane decreases to zero and gets left behind by structural waves. During the landing of airplane, the airplane meets two crests within 54.9 s. And thus the vertical motion of the airplane depends mainly on the relative velocity between the structural waves and the airplane.(a) The vertical displacement of locations of airplane. (b) The path of the landing run.
(a)
(b)


### 6.6. Takeoff in the Still Water

As the initial conditions for takeoff, the velocity and acceleration of the airplane are set to zero and constant, respectively. Before starting run, the position of the airplane is assumed to be at the starting point Z3 (−1000 m). Then the airplane suddenly runs from the rest; when timet=60.7 s, it completes the takeoff at position x=890 m. Comparisons between the present method and indirect time domain method used by Kashiwagi [20] are given in Figure 14 for the time histories of vertical deflections at measured points Z2, Z3, Z4, Z5, Z7, and Z9. The figure tells us that the Z3 position has an initial deflection 0.49 cm, which can be attributed to the static weight of the airplane. The two independent solutions of the direct time domain method and the indirect time domain method correlate well with each other.Time histories of the vertical displacement subjected to takeoff of airplane.
(a)
(b)The spatial profiles of the runway and the position indicated with circles during the airplane running are shown in Figure15. The static displacement of the runway by solving (30) can be seen at time t=0 s. Similar to the dynamic behavior of the landing, the position of the airplane stay always the bottom of sunken deflections of the runway, and the airplane always runs faster than the structural waves generated during the takeoff run. When the weight of airplane is equal to the lift force in (41), the airplane completes this takeoff which runs the distance from trough of static deflections to first peak. It can be seen from these deformed profiles that the deflections of the runway have small disturbance in early stage owing to slowly moving of the airplane. As time elapses, the velocity of the airplane and the perturbation of the runway also increase. After completing its takeoff, the structural waves are close to a steady equilibrium configuration, and move to the right at a certain speed. The amplitude of the deflections appears at time t=51 s and after time t=60.7 s for troughs and peaks, respectively, and they both stay within 0.9 cm.Spatial profiles of the runway subjected to takeoff run.
(a)
(b)


### 6.7. Takeoff in Regular Wave

Similar to the simulation of the landing run, the regular wave comes first then the takeoff load arrives later by three cycles of the wave period. The airplane is assumed to takeoff in the following incident direction, and the wave conditions are the same as mentioned in Section6.5.The spatial profiles of the runway during the takeoff and the locations of the airplane along the longitudinal centerline at different times are shown in Figure16. In the takeoff, just like the landing, the interaction of incident wave with the runway is more dominant as compared with the takeoff load. The figure tells us that the magnitude of the deflections which are nearly the same as the elevation of the structural waves is approximately 150 cm and is about 150 times the one induced in the still water condition. At the beginning stage of the run, the propagation velocity of structural waves is faster than the speed of the airplane (at least up 29 s), and then the airplane advances over a crest after t=52 s.Spatial profiles of the deflection subjected to takeoff of airplane in regular wave conditions.
(a)
(b)Figures17(a)-17(b) shows the history of the vertical displacement of position during run and the corresponding path of airplane. Since it takes a considerable amount of time for the airplane to move in the early stage, the fluctuation of the deflections is almost the same as the beginning of the structural wave induced by incident wave; however, when the airplane gains speed, the deflections of the runway change suddenly in their magnitude and length (20 s–60 s as shown in Figure 17(a)). After the speed of the airplane is close to and beyond the velocity of the structural wave, the airplane moves together with the structural wave and then overtakes the second crest at t=43 s. During the takeoff run of airplane, the airplane meets two crests within 60.7 s, and the vertical motion of the airplane depends mainly on the relative velocity between the structural waves and the airplane.(a) The vertical displacement of locations of airplane. (b) The path of the landing run.
(a)
(b)


## 6.1. Accuracy in the Interpolation-Tabulation Method

Before starting numerical simulations, it is necessary to confirm good performance in the computation of the time domain free-surface Green functions and its spatial derivatives in finite water depth.In order to examine the validity of the interpolation-tabulation method, computations of the nondimensional functionsF∞, F-F∞, and their spatial derivatives are performed for T=10, V=2 and are compared with corresponding results obtained from Newman. The values are shown in Figures 5 and 6. Obviously, the interpolation-tabulation method can give reliable evaluations which are smooth and in good agreement with Newman’s values.Comparison ofF∞, F∞X, and F∞V by the interpolation-tabulation method with Newman.
(a)
(b)
(c)Comparison ofF-F∞,  (F-F∞)X, and (F-F∞)V by the interpolation-tabulation method with Newman.
(a)
(b)
(c)


## 6.2. Drop Test in Still Water

The numerical simulation of the weight drop test is implemented, corresponding to the experiments conducted by Endo et al. [36] and the numerical results solved by Kashiwagi [19]. The pertinent information for the test model is prepared as listed in Section 4.1.Good convergence is considered for the numbers of modes in thex-direction and y-direction, after referring to the results of Kashiwagi [19], the number of modes in the x-direction M=8 and in the y-direction N=3 is adopted in this case.The comparisons of the vertical deflection time series at measured points are indicated in Figure2 among the present results; the indirect time domain solutions solved by Kashiwagi [19] and experimental tests obtained by Endo et al. [36] are given in Figure 7. It can be seen from this figure that the degree of agreement for these methods is favorable. And the deflections by the present method near the impact point, such as Z1 and Z2, are closer to the measurements than the numerical results solved by Kashiwagi [19]. This may be attributed to the difference in fluid pressure computation between the direct domain method by considering free-surface Green function and the indirect domain method by using the convolution integral of frequency impulse function.Comparison of time histories of the vertical displacement at each measured point.
(a)
(b)The deformed profiles of the VLFS during the mass drop are shown in Figure8. It is seen that the structural wave is transmitted at the longitudinal centerline of the plate, and the shape of the deformation is close to the current static equilibrium configuration at t=1.85 s. The magnitude of the vertical displacements is less than 1.0 cm. The vertical displacements of the plate are very small when x coordinate value is less than 0 but the transient phenomena at the right edge of the VLFS can be obviously seen at times t=0.21 s to 0.80 s.Figure 8
Spatial profiles of the VLFS due to the mass drop.


## 6.3. Drop Test in Regular Wave

In the simulation of the time step procedure, the load force at the right hand side of (30) is divided into two stages. The regular wave comes first from the left side of the VLFS and then the weight drops later by three cycles of the wave period. The wave length is 1.0 m, wave period is 0.8 seconds, wave height is 1 cm, and incident angle is 0 degrees.The deformed profiles of the VLFS during the weight drop are shown in Figure9, where Figure 9(a) shows the deflections in the regular wave condition without the mass impact for t=0 s and the deflections in the early stage of the drop test for t=0.21 and 0.41 s. The figure tells us that the absolute values of vertical displacement at the fore-end of the VLFS in regular wave are about 10 times the ones only generated by the mass drop; however, the magnitudes at the back-end are almost equivalent to the results induced by the mass drop. This means the mass impact load should not be overlooked as compared with the wave load. In other words, the structural wave shape of the VLFS in wave condition is changed when a huge mass falls off on the platform of the VLFS.Spatial profiles of the VLFS during the mass drop in regular wave conditions.
(a)
(b)


## 6.4. Landing in the Still Water

As mentioned in Section4.2, the airplane lands at point Z3 and completes the landing run in 54.9 s. The time histories of the vertical displacements at measurement points Z1, Z4, Z5, Z7, Z8, and Z9 obtained from the present method (direct time domain method) and the indirect time domain method used by Kashiwagi [20] are comparatively shown in Figure 10. The correlation between the two numerical solutions is reasonable. It can be seen that the magnitude of the deflections at measured points is less than about 1.0 cm and is very small as compared with the length value of the runway though the airplane weight is approximately up to 3900 kN. Then the time histories of the vertical displacement are much smoother than the results during weight drop (see Figure 7) and no higher-order deflections are found in time histories of the deflection for the landing on the platform of the VLFS. It is primarily due to the smooth increase of the landing loads (see (40) and (41)). It is also interesting that the vertical displacement of measured points Z7, Z8, and Z9 is not so large at t=55 s to 60 s but increases again after t=60 s. This can be attributed to the radiation of structure waves which impinge the stopped airplane.Time histories of the vertical displacement subjected to landing of airplane.
(a)
(b)Figure11 shows the snapshots of the deflection along the longitudinal centerline of the runway at different times, and the corresponding positions of the airplane are expressed by circles. It is found that the structural waves run after the airplane at time less than 42 s, and then the waves catch up to the airplane at about time 53 s because of the decrease of the airplane speed. After overtaking, structural waves meet the stopped airplane, partial waves are diffracted, and remainder is transmitted. Thus, it is interesting to find that the deformed profiles of the runway at time t=66 s are larger than those of the time t=53 s (see Figure 10). Note that the airplane seems to stay always at the bottom of sunken deflections of the runway during the landing run.Snapshots of the deflection subjected to landing of airplane.
(a)
(b)


## 6.5. Landing in Regular Wave

The interaction of incident wave with the runway during landing is divided into two stages. The regular wave comes first from the fore-end of the VLFS; after three cycles of the wave period, the airplane touches down the runway in the following incident wave direction. The wavelength is 650 m (0.13 times of the runway length), period 23.6 s, height 2.0 m, and incident angle is 0 degrees.The spatial profiles of the runway during the running are shown together with the positions of the airplane in Figure12. It can be seen that the airplane runs faster than the structural waves induced by the regular incident wave in the early stage (at least up 30 s) and then the structure waves overtake the airplane at the final stage of the run (after t=30 s). The maximum vertical displacement in regular wave is 150 cm and is about 150 times the one induced in the still water condition. This means that the wave load is dominant compared with the landing load in the hydroelastic analysis of VLFS.Spatial profiles of the deflection subjected to landing of airplane in regular wave conditions.
(a)
(b)Looking at the history of the vertical displacement of locations in Figure13(a) and the corresponding path during landing in Figure 13(b), it can be seen that the propagating velocity of the structural wave generated by incident wave is slower than the landing speed of airplane in the early stage (at least up to 30 s); however, when the airplane slows down, the deflections of the runway change suddenly in their magnitude and length (20 s–40 s as shown in Figure 13(a)). At the final stage of landing, speed of the airplane decreases to zero and gets left behind by structural waves. During the landing of airplane, the airplane meets two crests within 54.9 s. And thus the vertical motion of the airplane depends mainly on the relative velocity between the structural waves and the airplane.(a) The vertical displacement of locations of airplane. (b) The path of the landing run.
(a)
(b)


## 6.6. Takeoff in the Still Water

As the initial conditions for takeoff, the velocity and acceleration of the airplane are set to zero and constant, respectively. Before starting run, the position of the airplane is assumed to be at the starting point Z3 (−1000 m). Then the airplane suddenly runs from the rest; when timet=60.7 s, it completes the takeoff at position x=890 m. Comparisons between the present method and indirect time domain method used by Kashiwagi [20] are given in Figure 14 for the time histories of vertical deflections at measured points Z2, Z3, Z4, Z5, Z7, and Z9. The figure tells us that the Z3 position has an initial deflection 0.49 cm, which can be attributed to the static weight of the airplane. The two independent solutions of the direct time domain method and the indirect time domain method correlate well with each other.Time histories of the vertical displacement subjected to takeoff of airplane.
(a)
(b)The spatial profiles of the runway and the position indicated with circles during the airplane running are shown in Figure15. The static displacement of the runway by solving (30) can be seen at time t=0 s. Similar to the dynamic behavior of the landing, the position of the airplane stay always the bottom of sunken deflections of the runway, and the airplane always runs faster than the structural waves generated during the takeoff run. When the weight of airplane is equal to the lift force in (41), the airplane completes this takeoff which runs the distance from trough of static deflections to first peak. It can be seen from these deformed profiles that the deflections of the runway have small disturbance in early stage owing to slowly moving of the airplane. As time elapses, the velocity of the airplane and the perturbation of the runway also increase. After completing its takeoff, the structural waves are close to a steady equilibrium configuration, and move to the right at a certain speed. The amplitude of the deflections appears at time t=51 s and after time t=60.7 s for troughs and peaks, respectively, and they both stay within 0.9 cm.Spatial profiles of the runway subjected to takeoff run.
(a)
(b)


## 6.7. Takeoff in Regular Wave

Similar to the simulation of the landing run, the regular wave comes first then the takeoff load arrives later by three cycles of the wave period. The airplane is assumed to takeoff in the following incident direction, and the wave conditions are the same as mentioned in Section6.5.The spatial profiles of the runway during the takeoff and the locations of the airplane along the longitudinal centerline at different times are shown in Figure16. In the takeoff, just like the landing, the interaction of incident wave with the runway is more dominant as compared with the takeoff load. The figure tells us that the magnitude of the deflections which are nearly the same as the elevation of the structural waves is approximately 150 cm and is about 150 times the one induced in the still water condition. At the beginning stage of the run, the propagation velocity of structural waves is faster than the speed of the airplane (at least up 29 s), and then the airplane advances over a crest after t=52 s.Spatial profiles of the deflection subjected to takeoff of airplane in regular wave conditions.
(a)
(b)Figures17(a)-17(b) shows the history of the vertical displacement of position during run and the corresponding path of airplane. Since it takes a considerable amount of time for the airplane to move in the early stage, the fluctuation of the deflections is almost the same as the beginning of the structural wave induced by incident wave; however, when the airplane gains speed, the deflections of the runway change suddenly in their magnitude and length (20 s–60 s as shown in Figure 17(a)). After the speed of the airplane is close to and beyond the velocity of the structural wave, the airplane moves together with the structural wave and then overtakes the second crest at t=43 s. During the takeoff run of airplane, the airplane meets two crests within 60.7 s, and the vertical motion of the airplane depends mainly on the relative velocity between the structural waves and the airplane.(a) The vertical displacement of locations of airplane. (b) The path of the landing run.
(a)
(b)


## 7. Conclusions

A powerful direct time domain modal expansion method is applied to compute the transient behavior of a VLFS subjected simultaneously to incident wave and external loads including a huge mass drop and landing or takeoff load of an aircraft. The developed time-domain-Kelvin-source-based BEM solutions are carried out for the scatter wave problems, in which the free surface Green functions and its partial derivatives are rapidly and accurately evaluated by using the bilinear and trilinear interpolation-tabulation scheme. The assessed results of the auxiliary functions are generally in good agreement with the rigorous solutions.The computed results of the drop test show that, near the impact region, the displacement histories by the presented method are more consistent with measurements compared with the ones by indirect domain method. For regular wave conditions, the deformed profiles of the VLFS are changed when the mass falls on the VLFS, especially referring to near impact region. In the case of landing, the airplane runs faster than the structural waves in the early stage; as the airplane gradually stops moving, the generated waves catch up the airplane and partial waves are transmitted because of the presence of airplane. For the takeoff case, the runway has an initial deflection due to the static weight of the airplane. After the airplane leaves the runway, the generated waves are simple due to no disturbance on the VLFS. In the still water conditions, no higher order motion exists in the time histories of the vertical displacement and the locations of the airplane stay always at the bottom of sunken deflections of the runway during landing or takeoff run. However, in the following wave conditions, the displacement magnitude of the runway is greater than that only induced by airplane though the airplane weight reaches about 3900 kN, and the deflections of the runway due to the presence of airplane can be ignored. The airplane can surf on a series of progressive waves when the speed of airplane is closer to the one of structural waves.



---

*Source: 628564-2014-07-16.xml*