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

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

Assuming that the fluid is incompressible, inviscid, and irrotational, a velocity potential exists and would be given by where and are the incident and scattering potential, respectively.

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 of the floating body (the bottom surface) and (the side 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 is Poisson’s ratio and and denote the normal and tangential directions, respectively.

#### 3. Method of Solution

##### 3.1. 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 free-end condition of the beam and their convergence has already been proved by Newman [33]: where and are the symmetric and anti-symmetric modes about , respectively, and and are the positive real roots of equations and may be written in the same form, with replaced by on the right-hand sides of (12).

The modal functions expressed in (11) are orthogonal to each other in the () with the following orthogonality relation: where is Kroenecker’s delta, which is equal to 1 when 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: where represents the solid angle, 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 and memory term in finite water depth [34]: where the instantaneous term and memory term are given in the form, respectively, where 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 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 is defined by The first order time derivative is expressed in the form where

Thus the rate of convergence of (21) depends primarily on . Wehausen and Laitone [34] developed the following expression: where 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 .

If the spherical coordinate is adopted, the nondimensional parameters on are defined by where and lie in the interval and . 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.

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

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

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

##### 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 , and the position of the airplane and its velocity are given by where and are the initial position and velocity, respectively. For simplicity, the load distribution is assumed to be axisymmetric about the center of the moving load . 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 , and denotes the effective radius of the loading.

The total force 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.

#### 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 of and in (28) as follows.

The three variables , , of the function are changed to two arguments and , which are divided into 800 and 200 parts, respectively. The solutions for 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 and , and the time nondimensional parameters lies in the interval . 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).

##### 5.2. Symmetry of Structure

With the discretization of the constant boundary elements, (15) may be expressed as a form of the linear equations

Considering symmetry of the VLFS about - plane and - plane shown in Figure 4, the matrix , vector and , may be divided as follows: The symmetric relationships for the matrix may be formulated as In order to reduce the dimensions of the matrix, the conversions is obtained by taking where the constant coefficient for two planes of symmetry; the transition matrix is given by

Thus (43) can be simplified by substituting (46) where and the matrix 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 and . So the linear equations in region 1 will be In region 2, we have In region 3, we have In region 4, we have

#### 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 functions , , and their spatial derivatives are performed for , 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.

**(a)**

**(b)**

**(c)**

**(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 the -direction and -direction, after referring to the results of Kashiwagi [19], the number of modes in the -direction and in the -direction 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 free-surface Green function and the indirect domain method by using the convolution integral of frequency impulse function.

**(a)**

**(b)**

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 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 s to 0.80 s.

##### 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 Figure 9, where Figure 9(a) shows the deflections in the regular wave condition without the mass impact for s and the deflections in the early stage of the drop test for 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.

**(a)**

**(b)**

##### 6.4. Landing in the Still Water

As mentioned in Section 4.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 s to 60 s but increases again after s. This can be attributed to the radiation of structure waves which impinge the stopped airplane.

**(a)**

**(b)**

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 s are larger than those of the time 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.

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

**(a)**

**(b)**

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.

**(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 time s, it completes the takeoff at position 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.

**(a)**

**(b)**

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 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 s and after time s for troughs and peaks, respectively, and they both stay within 0.9 cm.

**(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 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 s.

**(a)**

**(b)**

Figures 17(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 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)**

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

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

The authors are grateful to the National Science Foundation for Creative Research Groups of China (Grant no. 50921001) for supporting this work.