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

and Applied Analysis 3 of the floating body Sb (the bottom surface) and Ss (the side surface): ∇ 2 ΦS (x, y, z, t) = 0, (2) ∂ 2 ΦS (x, y, z, t) ∂t + g ∂ΦS (x, y, z, t) ∂z = 0 on Sf, t > 0, (3) ∂ΦI (x, y, z, t) ∂n + ∂ΦS (x, y, z, t) ∂n = Vn on Sb + Ss, t > 0, (4) ∂ΦS (x, y, z, t) ∂n = 0 on Sd, t > 0, (5) ΦS (x, y, z, t) , ΦSt (x, y, z, t) , ∇ΦS (x, y, z, t) 󳨀→ 0 on S∞, t > 0, (6) ΦS (x, y, z, t) = ΦSt (x, y, z, t) = 0 on Sf, t = 0, (7) 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: D∇ 4 W(x, y, z, t) = −ms?̈? + P (x, y, z, t) − ρgW − PE, (8) 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 P (x, y, −d, t) = −ρ ∂Φ (x, y, −d, t) ∂t . (9) 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: ∂ 2 W ∂n + V ∂ 2 W ∂s = 0, ∂ 3 W ∂n + (2 − V) ∂ 2 W ∂n∂s = 0, (10) 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 displacement W(x, y, z, t) of the VLFS is the sum of various modes as follows: W(x, y, t) = N ∑ j=1 ζj (t) fj (x, y) = M ∑ m=0 N ∑ c=0 ζmc (t) fxm (x) fyc (y) , (11) where ζj(t) is the vibration amplitude of the j mode and fj(x, y) is modal function of the jmode.Themodal 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 freeend condition of the beam and their convergence has already been proved by Newman [33]: fx2m (x) = {{{{{ {{{{{ { 1 2 for m = 0 1 2 [ cosh (μ m x/ (L/2)) cosh (μ m ) + cos (μ m x/ (L/2)) cos (μ m ) ] for m = 1, 2, . . . , fx2m+1 (x) = {{{{{ {{{{{ { √3 2 x L/2 for m = 0 1 2 [ sinh (μ m x/ (L/2)) sinh (μ m ) + sin (μ m x/ (L/2)) sin (μ m ) ] for m = 1, 2, . . . , (12) where fx2m(x) and fx2m+1(x) are the symmetric and antisymmetric modes about x = 0, respectively, and μ m and μ m are the positive real roots of equations tanμ m + tanhμ m = 0, tanμ m − tanhμ m = 0 (13) 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: ∫∫ Sb fi (x, y) fj (x, y) dx dy = L 4 ⋅ B 4 δij, (14) 4 Abstract and Applied Analysis 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: αΦS (Q, t) + ∫∫ Sb(t)+Ss(t) ΦS (P, t) ∂G 0 (P, Q) ∂np dSp = ∫∫ Sb(t)+Ss(t) G 0 (P, Q) ∂ΦS (P, t) ∂np dSp


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][3][4][5][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][8][9]), the finite-element method (FEM) (see [10,11]), and the hybrid finite element-boundary element (FE-BE) method (see [12][13][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][16][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 highfrequency 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.

Runway
Incident angel Impact point Landing run Incident wave Dropping Figure 1: The fluid-structure system and coordinate system. This paper derives the expression of the time domain freesurface 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.

Formulation
We consider the time domain transient problem for a pontoon-type VLFS in finite water depth. Figure 1 shows the fluid-structure problem and Cartesian coordinate system. The -axis is pointing upwards, and the -plane is on the mean position of the free surface, where ℎ is the water depth and is the amplitude of the incident wave. The whole fluid domain is defined at Ω which contains the bottom of the VLFS , side of the VLFS , the free surface , the seabed , and the infinite cylindrical surface ∞ . The VLFS has a length , width , and height ℎV, and is the draft of the VLFS in direction. The problem at hand is to determine the modal deflections under external loads combined action of incident waves.
The velocity potential must satisfy the following Laplace's equation and boundary conditions on the free surface , on the sea-bed , the infinity ∞ , and on the wetted surface where is the gravitational acceleration, is the normal velocity of the structure, is a unit normal vector (the positive direction points out of the fluid domain), and Φ ( , , , ) 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: where = is the bending rigidity, is modulus of elasticity, is the cross sectional moment of inertia, denotes the mass per unit area, denotes density of fluid, 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 ( , , − , ) relates to the velocity potential on the bottom surface of the VLFS from the linearized Bernoulli's equation 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: where V is Poisson's ratio and and denote the normal and tangential directions, respectively.

Modal Functions.
The vertical elasticity displacement ( , , , ) of the VLFS is the sum of various modes as follows: where ( ) is the vibration amplitude of the mode and ( , ) is modal function of the mode. The modal function ( , ) of the VLFS can be expressed by the product of modal functions at and directions. Here, a case study of the natural modes ( ) and ( ) 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 freeend condition of the beam and their convergence has already been proved by Newman [33]: where 2 ( ) and 2 +1 ( ) are the symmetric and antisymmetric modes about = 0, respectively, and and are the positive real roots of equations and ( ) may be written in the same form, with /2 replaced /2 by on the right-hand sides of (12).
The modal functions expressed in (11) are orthogonal to each other in the (− /2, /2) with the following orthogonality relation: where is Kroenecker's delta, which is equal to 1 when = and 0 otherwise.

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: where represents the solid angle, ( 0 , 0 , 0 ) and ( , , ) represent the source and field point, respectively, and ( ) represents the instantaneous waterline of the intersection between the body and the free surface. For Green's function ( , , , ), it can be expressed by the superposition of instantaneous term 0 and memory term in finite water depth ℎ [34]: where the instantaneous term 0 and memory term are given in the form, respectively, where 0 is the Bessel function of the first kind, order zero, denotes the horizontal distance between field and source point, denotes the distance between field and source point, and 2 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 (17) may be written in the form where the auxiliary function 0 ( , ) is defined by The first order time derivative is expressed in the form where Abstract and Applied Analysis 5 Thus the rate of convergence of (21) depends primarily on . Wehausen and Laitone [34] developed the following expression: where 0 is the modified Bessel function of the second kind and the coefficients are defined by Newman [35]. The rate of convergence of the integral in (23) can be accelerated by adding and subtracting the appropriate function ∞ which may be given in the form where the vertical coordinate is restricted to the fluid domain (−1, 2).
If the spherical coordinate is adopted, the nondimensional parameters on are defined by 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 ∞ and its partial derivatives are obtained: Since the integrands in (27) exhibit slow convergence and high oscillatory, the ∞ and its spatial derivatives can be approximated in terms of the values of parameter (see [30,35]). The function − ∞ may be solved directly in a straightforward manner due to the oscillatory elimination, and thus ( , , ) is obtained: 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.

Structure
Part. We substitute the fluid pressure and vertical deflection (9) and (11) into (8); we have where , are the modal numbers in and direction, respectively. Applying Galerkin's method, we multiply both sides of the above equation by ( )⋅ ( ) and integrate over the bottom of the VLFS. Finally, we can obtain a conventional set of equations given by where It should be noted that , , , , , and are the generalized mass, generalized stiffness, generalized wave force, and generalized external load force, respectively. And the generalized stiffness , 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.

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

The Weight Drop Test.
In a weight drop test, a weight was dropped from a height onto the "hit point. " The acceleration of the weight during the impact was . Therefore, the impact load im ( ) can be obtained: and the external pressure distribution , appearing in (34) can be expressed as where ( , ) is the coordinate of the hit point. Substituting (36) into (34), the generalized external force ( ) can be computed as The numerical parameters for simulation in this paper are given in Table 1 by referring to Endo et al. [36], and the vertical displacements at points Z1-Z9 shown in Figure 2 were measured.

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 ( ) of the airplane and its velocity ( ) are given by where 0 and 0 are the initial position and velocity, respectively. For simplicity, the load distribution is assumed to be  axisymmetric about the center of the moving load ( ( ), 0). In the terms of the relationship between the moving Cartesian coordinate system − and the polar coordinate system − , the external pressure distribution can be expressed as where = √ 2 + 2 , = − ( ) and denotes the effective radius of the loading. The total force ai ( ) exerted by the landing or takeoff on the VLFS can be given by the difference between the weight of the airplane and the lift force ( ): where the parameters and are given as constants, is the density of air, and is the effective wing area of the airplane.
Substituting (39) into (34), the generalized external force ( ) can be computed as The numerical data for simulation in this paper is prepared as listed in Table 2 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.

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 ∞ and − ∞ in (28) as follows.   The three variables , , of the function ∞ 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 ∞ and its derivation.
However, the − ∞ function has three arguments , , and . Here, the space nondimensional parameters and are restricted to the region (−1, 2) and (0, 20), and the time nondimensional parameters lies in the interval (0, 20). First seven -planes are adopted in direction. Next, every -plane is divided into 40 parts in the and direction, respectively. The slowly varying function − ∞ and its derivations can be calculated by using Gaussian integration; then the trilinear interpolation scheme is applied to the effective approximation in (15).

Symmetry of
In order to reduce the dimensions of the matrix, the conversions is obtained by taking In region 2, we have In region 3, we have In region 4, we have

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 interpolationtabulation method, computations of the nondimensional functions ∞ , − ∞ , and their spatial derivatives are performed for = 10, = 2 and are compared with corresponding results obtained from Newman. The values are shown in Figures 5 and 6. Obviously, the interpolationtabulation method can give reliable evaluations which are smooth and in good agreement with Newman's values.

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 the -direction and -direction, after referring to the results of Kashiwagi [19], the number of modes in thedirection = 8 and in the -direction = 3 is adopted in this case.
The comparisons of the vertical deflection time series at measured points are indicated in Figure 2 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 freesurface 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 Figure 8. 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 = 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 coordinate value is less than 0 but the transient phenomena at the right edge of the VLFS can be obviously seen at times = 0.21 s to 0.80 s.

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 Abstract and Applied Analysis   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 = 55 s to 60 s but increases again after = 60 s. This can be attributed to the radiation of structure waves which impinge the stopped airplane. Figure 11 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 = 66 s are larger than those of the time = 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.  Figure 12. 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 = 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.
Looking at the history of the vertical displacement of locations in Figure 13(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.
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 time = 60.7 s, it completes the takeoff at position = 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 spatial profiles of the runway and the position indicated with circles during the airplane running are shown in Figure 15. The static displacement of the runway by solving (30) can be seen at time = 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 = 51 s and after time = 60.7 s for troughs and peaks, respectively, and they both stay within 0.9 cm.

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 Section 6.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 Figure 16. 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 = 52 s.   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 = 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.

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