INVESTIGATION OF HYDROELASTIC BEHAVIOR OF A PONTOON-TYPE VLFS DURING UNSTEADY EXTERNAL LOADS IN WAVE CONDITION USING A HYBRID FINITE ELEMENT-BOUNDARY ELEMENT ( FE-ME ) METHOD

The hydroelastic behavior of a pontoon-type VLFS subjected to unsteady external loads in wave condition is investigated in the context of the time-domain modal expansion theory, in which the boundary element method (BEM) based on time domain Kelvin sources is used for hydrodynamic forces and the finite element method (FEM) is adopted for solving the deflections of the VLFS. In this analysis, the interpolation-tabulation scheme is applied to assess rapidly and accurately the free-surface Green function in finite water depth, and the boundary integral equation of a quarter VLFS model is further established taking advantage of symmetry of flow field and structure. The VLFS is modeled as an equivalent solid plate based on the Mindlin plate theory. The coupled plate-water model is performed to determine the wave-induced responses and transient behavior under external loads such as a huge mass impact onto the structure and moving loads of an airplane, respectively. These results are verified with existing numerical results and experimental test. Then, the developed numerical tools are used in the study of the combined action taking into account of the mass drop/airplane landing as well as forward or reverse incident wave action. The deflections of the runway, the time history of vertical positions and the trajectory of the airplane are also presented through a systematic time-domain simulation, which illustrates the usefulness of the presently developed numerical solutions.


Introduction
Mat-like very large floating structures (VLFS) are very flexible offshore structures widely regarded as an alternative option of utilizing ocean space for different applications such as floating airports, oil storage vessels, bridges, floating artificial island, etc.When a VLFS moves due to wave or moving load actions, the fluid surrounding the structure changes at the same time.On other hand, the pressure distribution over the body surface also changes in order to satisfy the Bernoulli's equation.If the pressure changes, the motions including rigid body motions and elastic deformations are affected.Thus, the fluid-structure interaction problem is called hydroelasticity and it is necessarily considered in order to obtain the vertical responses of a VLFS.
Numerous contributions have been performed on the research activities related to hydroelastic analysis of VLFS.Watanabe et al. [1] and Eatock Taylor [2] presented a review of these methods on the simulation of hydroelastic responses of VLFS.One way to tackle this problem is to use an analytical approach given by Wu et al. [3], Watanabe et al. [4], Pham et al. [5], and Karmakar and Soares [6].They divided the whole fluid domain into multiple regions and applied the eigenfunction expansion-matching method for obtaining the waveinduced deflections of the VLFS.Based on Mindlin thick plates theory and Wiener-Hopf technique, Zhao et al. [7] presented the dynamical behavior of floating elastic plates acted upon by a localized external load (periodic load).If these analytical methods are used, the computational time and memory capacity for a VLFS are non-issue.However, the solved drawback is only applied to simple geometries such as a rectangular plate and a circular plate.Another way to obtain the hydroelastic responses of VLFS is by using a numerical method.The numerical solutions can be applied to more any shape of VLFS, and the boundary element method (BEM) [8][9][10][11], the finite-element method (FEM) [12][13][14][15][16][17][18], and the hybrid finite element-boundary element (FE-BE) method [19][20][21][22] have been presented in previous studies.
Among various numerical studies, the frequency domain wave-induced behavior has been commonly investigated when determining the hydroelastic response amplitude operator (RAO) of the floating body and pertinent response parameters in a steady state condition.However, in the real situation, we would often like use some time series of the responses, and non-harmonic external loads such as imposed by a huge mass impact onto the structure, landing and taking off of an airplane, can induce the transient behavior of the VLFS and may affect the serviceability of the VLFS.Watanabe et al. [12], Qiu and Liu [13], and Shin et al. [14] applied the finite element methods (FEM) to the transient response analysis of a VLFS due to impulsive landing/takeoff of an airplane, however, the response analysis is required for a few seconds and this simulation model is very much simplified in the treatment of structure or fluid.Kashiwagi [8,9] developed an indirect time-domain method for calculating transient responses of a VLFS, in which the hydrodynamic effect is evaluated from good performance in the computation of the memory-effect function.Lee and Choi [22] proposed a boundary element-finite element (BE-FE) 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.Based on the BE-FE combined method, Endo [23 and 24] simulated the transient behavior of a VLFS subjected simultaneously to takeoff/landing and regular wave load.These mentioned studies provide enlightening contributions in the research activities related to external loads on the VLFS but some difficulties in carrying out their time-domain simulation give restriction on the mathematical model, i.e. the assumption of both small structural motion and wave amplitude, the integration of memory-effect function, and the accurate evaluation of hydrodynamic coefficients.Our approach, instead, is to obtain directly the wave-induced responses and transient phenomena Yong Cheng, Chunyan Ji, Gangjun Zhai, Tianhui Fan behavior of a pontoon-VLFS during unsteady external loads in wave condition 25 by using a direct time-domain modal expansion method.In this analysis, the numerical model and scheme developed by the same authors [21 and 25] for a VLFS edged with dual antimotion plates under wave action (without any unsteady external loads) is extended to incorporate the non-harmonic external loads.The method takes account of applying the timedomain free-surface Green functions in the hydrodynamic diffraction and radiation problems whereas the element used to establish the plate model is the non-conforming quadraticserendipity (NC-QS) Mindlin plate element [20].Such a Mindlin plate element does not suffer from spurious modes and shear locking phenomenon, and considers the effects of shear deformation and rotary inertia, which are neglected in the Classical thin plate theory.The fluid-structure interaction is investigated using a hybrid finite element-boundary element (FE-ME) method, where the FEM is used to solve the coupled plate-water equation whereas the BEM is used to handle the water integral boundary equation.Compared with FEM, the modal expansion method can clearly reduce the number of unknown amplitudes which can be obtained directly from a solution of simultaneous equations at each time step.
Fast and accurate calculation is necessary for overcoming this difficulty of large CPU time and memory size of computer.Utsunomiya et al. [26] has developed the multipole expansion methods for hydroelastic analysis of a VLFS.Kagemoto et al. [27] presented a substructural method that accelerates computation without an appreciable loss of accuracy.Dai [28] has expanded the precorrected-FFT method to hydroelastic analysis.However, their calculation models only satisfy the frequency domain studies for the wave-induced hydroelastic response.Huang [29] 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.Thus the authors derive 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 an 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 using BEM.According to the symmetry of the VLFS structure and the fluid field [30,31], this paper is concerned with numerical simulations of boundary integral equations of a quarter VLFS model.

Mathematical formulation
Fig. 1 The fluid-structure system and coordinate system For the time domain elastic responses of a pontoon-type VLFS in finite water depth, Fig. 1 shows the fluid-structure problem and Cartesian coordinate system.In the model, the zaxis is pointing upwards, and the x-y plane is on the mean position of the free surface, where h is the water depth, 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 A direct time-domain simulation of hydroelastic Yong Cheng, Chunyan Ji, Gangjun Zhai, Tianhui Fan behavior of a pontoon-VLFS during unsteady external loads in wave condition 26 Sf, the seabed Sd and the infinite cylindrical surface S∞.The VLFS has a length L, width B, 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 wave.
Fluid is assumed to be ideal, so that a velocity potential exists and is governed by the Laplace equation: 2 ( , , , ) 0 x y z t    ( 1 ) and the boundary conditions are satisfied as the following: are the incident and scattering potential, respectively.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), By assuming the plate material to be isotropic and obey the Hooke's law, the motion of the floating body is governed by the equation of Mindlin thick plate [32], i.e. () where W is the displacement vector including of the vertical deflection W(x,y,t), the rotation Ψx(x,y,t) about the y-axis and the rotation Ψy(x,y,t) about the x-axis.D and ρs denote the bending rigidity and the structural density, respectively.The differential operators B1 and B2 and the constant matrix B3 are defined in Ref. [18].The loads vector F comprises the hydrostatic pressure -ρgW, hydrodynamic pressure 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) where [Df], [Ds], {χ} and {γ} are the flexural elasticity matrix, the shear elasticity matrix, the flexural strain and the shear strain, respectively, and have been given in Ref. [20].
A direct time-domain simulation of hydroelastic Yong Cheng, Chunyan Ji, Gangjun Zhai, Tianhui Fan behavior of a pontoon-VLFS during unsteady external loads in wave condition 27

Method of solutions
The solution for the hydroelastic responses of the VLFS involves solving the platewater motion equation and the water boundary integral equation by applying the hybrid finite element-boundary element (FE-ME) method.Here, the FEM with non-conforming quadraticserendipity (NC-QS) Mindlin plate element is applied to solving the time series of the plate responses whereas the BEM is used to obtaining the water scattering potential.

Fluid part
The boundary value problems given by Eqs. ( 1)-( 6) can be solved by using the Green's function method.If the free-surface Green's function satisfying the boundary conditions given by Eqs ( 2), ( 4), ( 5), and ( 6) is considered, the boundary integral equation for the scattering potential can be derived as follows: 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 function G(P,t,Q,τ), it can be expressed by the superposition of instantaneous term G 0 and memory term G f [34]: where the instantaneous term G 0 and memory term G f are given in the form, respectively 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 seabed.We non-dimensionalise the spatial parameters using X=R/h, Y=-z0/h, Z=-z/h and the time parameter using ) / )( ( . The instantaneous term G 0 is evaluated using the method of Cheng et al. [25], and the first order time derivative f G  is expressed in the form: A direct time-domain simulation of hydroelastic Yong Cheng, Chunyan Ji, Gangjun Zhai, Tianhui Fan behavior of a pontoon-VLFS during unsteady external loads in wave condition 28 The rate of convergence of the integral in Eq. ( 21) can be accelerated by adding and subtracting the appropriate function  F which may be given in the form ( 2) 0 0 ( , , ) lim ( , , ) where the vertical coordinate V is restricted to the fluid domain (-1, 2).
If the spherical coordinate is adopt, the non-dimensional parameters on r are defined by where τ and θ are lies in the interval (0, ∞) and (0, 2π).Thus, we will reduce three arguments to two arguments by substituting Eq. ( 18) into Eq.( 17).The function  F and its spatial derivatives are obtained: Since the integrands in Eqs. ( 19)-( 22) exhibit slowly convergence and highly oscillation, the F∞ and its spatial derivatives can be approximated using series expansion, asymptotic expansion or Filon quadrature in terms of the values of parameter τ [29 and 34].The function F-F∞ may be solved directly in a straightforward manner due to oscillatory elimination, thus Then the boundary surface of Eq. ( 11) 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 [35], which represent the geometry of each element.In the integration process, the scheme using trapezoidal approximation is applied to the convolution integral.Once Eq. ( 11) is solved, the time history of fluid dynamic pressure in Eq. ( 8) can be obtained at any position.

Structure part
In order to use the FEM for solving the plate equation, Eq. ( 7) is transformed into an equivalent motion equation by using first-order shear deformation plate theory based on Hamilton's principle.The shear correction shear factor is taken as 5/6 according to Ref. [20,[36][37][38].The VLFS model is approximated by a number of the NC-QS Mindlin thick plate elements, in which additional non-conforming basis functions are added to the bending rotations, and thus each element does not exhibit spurious modes and shear locking phenomena.By transforming the coupled plate-water Eq. ( 7) into the equivalent variational equation and minimizing energy functional, we obtain the global form of coupled plate-water equation where e e F f N N P ds where [N] is the NC-QS element interpolation function matrix; superscript T denotes matrix inversion; e means numerical formulation on elements; [Bf] is the flexural strain-displacement matrix whereas [Bs] is the shear strain-displacement matrix as found in [18].Here, the global restoring force matrix Kre is only related to the hydrostatic pressure term in Eq. ( 8), and ignores the effect of rigid body motion [39 and 40].The displacement vector {W} is expressed in terms of modal expansion as follows: Here ζj is the jth generalized modal coordinate; subscript i denotes the node number and j denotes the mode number; fij w ,  are jth natural modes corresponding to vertical displacement, rotation about the y-axis and rotation about the x-axis, respectively, and they can be obtained by solving the eigenvector equation where [λ] is the square diagonal matrix of jth natural circular frequency.By substituting Eq. ( 31) into VLFS-water Eq. ( 24) and premultiplying both sides of the Eq. ( 24) by [f ] T , we can obtain a conventional set of equations given by with A direct time-domain simulation of hydroelastic Yong Cheng, Chunyan Ji, Gangjun Zhai, Tianhui Fan behavior of a pontoon-VLFS during unsteady external loads in wave condition 30 where [GM], [GC], [GK], {GFw} and {GFE} are the generalized mass matrix, the generalized damping matrix, generalized stiffness matrix, generalized wave force vector and generalized external load force vector, respectively.The second order linear differential equations for modal amplitudes shown in Eq. ( 33) are solved at each time step using the fourth order Runge-Kutta method.Then the modal responses are summed up to obtain the total response.

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 of F∞ and F-F∞ in Eq. ( 23) as follows.
Referring to the Eqs.( 19)-( 22), the three variables X, V, T of the function F∞ are changed to two arguments  and  cos , which are divided 800 and 200 parts, respectively.The solutions of F∞ and its partial derivatives for cosθ<0.7 which are efficient in the context of section 3.2 are described by Huang [29], else the Filon integral scheme is determined to calculate directly.During solving the Eq. ( 11), the bilinear interpolation scheme was applied to the effective approximation of F∞ and its' derivation.
However, the F-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), 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 40 parts in X and T direction, respectively.The slowly-varying function F-F∞ and its' derivations can be calculated by using Gaussian integration, then the tri-linear interpolation scheme is applied to the effective approximation in Eq. (11).

Symmetry of structure
With the discretization of the constant boundary elements, the Eq. ( 11) may be expressed as form of the linear equations Considering symmetry of the VLFS about x-z plane and y-z plane shown in Figure 2, the matrix [A], vector   S  and   B , may be divided as follows:

A A A A B A A A A B B A A A A B
It is noted that the subscripts in matrix [A] and vector {B} or superscripts in vector {ΦS} denote the related region of matrix or vector, for example, [A11] is the sub-matrix A direct time-domain simulation of hydroelastic Yong Cheng, Chunyan Ji, Gangjun Zhai, Tianhui Fan behavior of a pontoon-VLFS during unsteady external loads in wave condition 31 corresponding to the region 1, and [A12] is the sub-matrix corresponding to the region 1 and 2, and so on.In addition, the symmetric relationships for the matrix [A] may be formulated: In order to reduce the dimensions of the matrix, the conversions is obtained by taking where the matrix [E] is a transition matrix and the constant coefficient β=4.After substituting Eq. ( 44) into Eq.( 38), the linear system of equations may be obtained for each region from 1 to 4 defined in Fig. 2 and they have been given by Ref. [11].

Accuracy in the interpolation-tabulation
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 non-dimensional functions F∞ , F-F∞ and their time derivatives are performed for and compared with corresponding results obtained from Newman.The values are shown in Figs. 3 and 4. Obviously, the interpolation-tabulation method can give reliable evaluations which are smooth and good agreement with the Newman's values.

Hydroelastic response of VLFS in wave conditions
In order to validate the present formulation and method in the foregoing sections, we first performed the hydroelastic response of VLFS for the loads of regular wave conditions by using the direct time-domain method and the numerical solutions are compared with the analytical results solved by Wu et al. [3] and experimental tests obtained by Utsunomiya et al. [19].The pertinent information for the VLFS design used by Utsunomiya et al. [19] is listed in Table 1.Based on convergence tests, it is found that the ratio of element size to wavelength must be smaller than 0.1, and the vertical deflection obtained is confirmed to be stable and periodic after three cycles of the wave period.In addition, The time step defined as T/60 and modes N=30 are sufficient for the deflections to converge of VLFS.
Figs. 5 and 6 show the distribution of vertical displacement amplitude and bending moment amplitude to incident regular wave with wave periods of T=1.429 s and 2.875 s.In this analysis, the element size to wavelength ratio must be smaller than 0.1 after checking the convergence, and the vertical deflection obtained is confirm to be stable and periodic after three cycles of the wave period.In the two cases, the good agreement is shown between our numerical solutions and Wu et al.'s analytical solutions, and the correlation between the present method and the experimental data is reasonable.The slight discrepancy between the calculated and measured data can be attributed to the possibly un-captured physics (such as vortices and wave breaking), the fluid viscosity, the energy loss induced by the gap between structure and tank wall, inherent nonlinear effects of wave, and the instrument accuracy.The transient phenomena of a weight drop test are then implemented on the VL-10 model, corresponding to the experiment conducted by Endo and Yago [24].In this experiment, the model has a length 9.75m, width 1.95m, thickness 0.0163m and bending rigidity 8985.62Nm.The weight of Wm=196N was dropped from a height of 0.12m onto the "hit point" (see Fig. 7) and then the impact load F0(t) can be obtained: where a denotes the acceleration of the weight during the impact and it has been measured from Kashiwagi [8].The external pressure distribution E P , appearing in Eq. ( 30) can be expressed as where (xp, yp) is the coordinate of the hit point.
Fig. 8 comparatively shows the time histories of the vertical displacement at measured points Z1-Z9 indicated in Fig. 7 among the present method, the numerical results solved by Kashiwagi [8] and experimental tests obtained by Endo and Yago [24].The present numerical results correlate reasonably with the Kashiwagi's numerical solutions and Endo and Yago's experimental results.It is also interesting that the deflections by the present method near the impact point, such as Z1 and Z2, are closer to the measurements than the Kashiwagi's numerical results.This may be due 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.The deformed profiles of the VLFS during the mass drop are shown in Fig. 9.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 current static equilibrium configuration at t=1.85s.The magnitude of the vertical displacements is less than 1.0cm.The vertical displacements of the plate is very small at x coordinate value along the length less than 0, and the transient phenomena at the right edge of the VLFS is obviously seen from t=0.21s to 0.80s.We consider the realistic landing of an airplane on a VLFS with rectangular geometry in plane shown in Fig. 10.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 V0 are given by 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 , and R denotes the effective radius of the loading.

The total force
) ( 0 t F 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  The numerical data for simulation in this paper are prepared as listed in Table 2 by referring to Kashiwagi [9].The airplane lands at touch-down point Z3 indicated in Fig. 10 and it completes the landing run in 54.9s.The time histories of the vertical displacements at measurement point Z1, Z4 ,Z5, Z6, Z7 and Z9 obtained from the present method(direct time domain method), and the indirect time domain method used by Kashiwagi [9], are comparatively shown in Fig. 11.The simulated results generally follow the trend of the numerical curve solved by Kashiwagi [9].It is seen that the magnitude of deflections is less than about 1.2 cm and is very small as compared with the length value of the runway though the airplane weight reaches about 3900 KN.Then the time histories of the vertical displacement is much smoother than the results shown in Fig. 8 during weight drop, so no higher-order deflections are found in time histories of the deflection during the landing on the platform of the VLFS.It is also interesting that the vertical displacement of measured points Z7 and Z9 is not so large at t=55s to 60s but increases again after t=60s.This is mainly due to be the radiation of structure waves which impinge the stopped airplane.12 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 behind the airplane at time less than 42s, and then the waves catch up to the airplane at about time 53s 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 found that the deformed profiles of the runway at time t=66s is larger than that of the time t=53s (see also Fig. 11).Note that the airplane seems to stay always sunken deflections of the runway during the landing run.
A direct time-domain simulation of hydroelastic Yong Cheng, Chunyan Ji, Gangjun Zhai, Tianhui Fan behavior of a pontoon-VLFS during unsteady external loads in wave condition 37

Drop test in regular wave
In the simulation of a weight drop test in regular wave, the load force at the right hand side of Eq. ( 24) is divided into two stages.The regular comes first from the front side of the VLFS then the weight drops later by three cycles of the wave period.The wave length 1.0 m, wave period is 0.8 seconds, wave height is 1 cm and incident angle is 0 degree.
The deformed profiles of the VLFS during the weight drop are shown in Fig. 13, where Figure 13a shows the deflections in the regular wave condition without the mass impact for t=0s, and the deflections in the early stage after impact for t=0.21 and 0.41s.The figures tell 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 magnitude at the back-end is almost equivalent to the one induced by the mass drop.This means the mass impact load should be not overlooked as compared with the wave load.And the structural waves of the VLFS in wave condition are changed when the huge mass falls off on the platform of the VLFS.As mentioned section 5.5, the regular wave comes first from the fore-end of the VLFS, after three cycles of the wave period, the airplane touch down the runway in the following incident direction.The wavelength is 650m (0.13 times of the runway length), period 23.6s, height 2.0m and incident angle 0 degree.
The spatial profiles of the runway during the running are shown together with the positions of the airplane in Fig. 14.From which it can be seen that the maximum vertical displacement in regular wave is 150cm and is about 150 times the one induced in the still water condition.This means the wave load is dominant with compared with the landing load in the hydroelastic analysis of VLFS.
Looking the history of the vertical displacement of locations in Figure 15a and the corresponding path during landing in Figure 15b, 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 20s), however, when the airplane slows down, the deflections of the runway change suddenly in their magnitude and length (20s-40s shown in Fig. 15a).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 structural waves within 54.9s.And thus the vertical motion of the airplane depends mainly on the relative velocity between the structural waves and the airplane.Finally, we consider the run direction of the airplane heading toward the incident wave direction for the purpose.The regular wave comes first from the back-end of the VLFS, after three cycles of the wave period, the airplane touch down the runway heading toward incident direction.The detail of VLFS and wave conditions are given in the previous section.
Fig. 16a shows the history of the vertical location during landing and Fig. 16b shows the deformed profiles of the runway in early stage and later stage of landing together with the path of the airplane.Compared with the curves in Fig. 15, it can be seen that the magnitude of the deflections is nearly the same and it has little to do with the run direction of the airplane owing to the dominant role of wave action.In this particular case, the airplane meets five structural waves within 54.9s of landing run and its vertical motion mainly depends on the structural wave propagation which is raised by incident waves.In Fig. 16a, the periods of the history of vertical displacement during landing are smaller than the ones after the airplane stops run which corresponds to the speed of landing decreasing to zero.

Conclusions
We have presented a time-domain hybrid finite element-boundary element (FE-ME) scheme with a robust tool for investigating the hydroelastic responses of a VLFS subjected to external loads including a weight impact onto the structure, and landing of an aircraft taking into account of the combined action of incident wave.For fluid part, the hydrodynamic forces were directly obtained by using the developed time-domain-Kelvin-source-based BEM solutions but this calculation of Green functions and its partial derivatives requires huge CPU time and memory size of computer.In order to surmount this difficulty, the bilinear and trilinear interpolation-tabulation schemes were present and applied to the fluid-structure interaction.For structural part, the coupled plate-water equation is solved using the FEM with NC-QS Mindlin plate elements.The numerical results were examined for the load cases of regular incident wave condition and of external loads including mass impact and moving load condition, respectively.
After verification by numerical and experimental results, the developed computer program was further used to simulate the hydroelastic motion of the VLFS during a weight mass drop or landing of as airplane in wave condition.From these numerical results, it was found that the vertical motion of the VLFS in regular waves is changed when the mass fall on the VLFS, even though the weight of the mass is 196 N.However, as a result of transient responses during landing, the magnitude of deflections of runway in wave conditions is dominant as compared with the results only generated by airplane.The deflections of the runway due to the presence of airplane after landing run, can be ignored in wave condition.The vertical motion magnitude of the airplane can change abruptly during landing run in the following sea condition due to the reduction of relative velocity between the airplane and structural waves.In the heading toward incident wave direction, the airplane runs up and down on the structural waves induced by incident waves.

A
direct time-domain simulation of hydroelastic Yong Cheng, Chunyan Ji, Gangjun Zhai, Tianhui Fan behavior of a pontoon-VLFS during unsteady external loads in wave condition29

Fig. 2
Fig.2 Sketch of the district of symmetry

Fig. 3 32 Fig. 4
Fig. 3 Comparison of F∞ and F∞T by the interpolation-tabulation method with Newman

Fig. 5 Fig. 6
Fig.5 Comparison of numerical and experimental results for Utsunomiya et al.'s VLFS for regular wave period T=1.429(a) vertical displacement amplitude (b) bending moment amplitude

34 Fig. 7 Fig. 8
Fig. 7 The position of the measured points and hit point in the drop test

Fig. 9 Fig. 10
Fig. 9 Spatial profiles of the VLFS due to the mass drop


is the density of air , and W A is the effective wing area of the airplane.

Fig. 11 Fig. 12
Fig. 11 Time histories of the vertical displacement subjected to landing of airplane

Fig. 13
Fig. 13 Spatial profiles of the VLFS during the mass drop in regular wave conditions 5.6 Landing in the following sea condition

Fig. 14
Fig.14 Spatial profiles of the deflection subjected to landing of airplane in regular wave conditions

Fig. 16 (
Fig. 16 (a) The time history of the vertical location during landing in the heading toward sea condition (b) The path of the Landing run in the heading toward sea condition , [C] is the global viscous damping matrix, which can be neglected in context of potential theory;[M], [Kf], [Ks], [Kre], {Fw} and {FE} are the global mass matrix, the global flexural stiffness matrix, the global shear stiffness matrix, the global restoring force matrix, the global wave force vector and the global external load force vector, respectively.All of these entities can be assembled from corresponding single element matrix [m]e, [kf]e, [ks]e, [kre]e, {fw}e and {fE}e, i.e.

Table 1
Principal details of the VLFS model and sea states Chunyan Ji, Gangjun Zhai, Tianhui Fan behavior of a pontoon-VLFS during unsteady external loads in wave condition

Table 2
Main parameters of landing run