Time-Domain Three Dimensional BE-FE Method for Transient Response of Floating Structures Under Unsteady Loads

This paper presents a direct time-domain three dimensional (3D) numerical procedure to simulate the transient response of very large floating structures (VLFS) subjected to unsteady external loads as well as moving mass. The proposed procedure employs the Boundary Element and Finite Element methods (FEMBEM).The floating structure and the surrounding fluid are discretized by 4-node isoparametric finite elements (FE) and by 4node constant boundary elements (BE), respectively. Structural analysis is based on Mindlin’s plate theory. The equation of motion is constructed taking into account the effect of inertia loading due to the moving mass. In order to obtain the hydrodynamic forces (added mass and radiation damping), the coupled natural frequencies are first obtained by an iterative method, since hydrodynamic forces become frequency-dependent. Then the Newark integration method is employed to solve the equation of motion for structural system. In order to prove the validity of the present method, a FORTRAN program is developed and numerical examples are carried out to compare its results with those of published experimental results of a scale model of VLFS under a weight drop and airplane landing and takeoff in still water condition. The comparisons show very good agreement.

Latin American Journal of Solids and Structures 13 (2016) 1340-1359 tic deformations of VLFS become dominant over the rigid body motions, and the fluid-structure interaction must be included for accurate analysis.
Many researchers have studied the wave-induced hydrodynamic response of floating structures.The wave response of mat-like structure has been studied by a 1D straight-beam and 2D potential flow theory (Wen and Shinozuka, 1972;Suzuki, 1994).Wen (1974) studied the same problem by 2D rectangular plate and 3D potential flow theory.Wang et al. (1995) presented a hydroelastic response analysis of a box-like floating airport of shallow draft.Hamamoto and Fujita (1995) analyzed the wave response of 3D-module-type floating structure by BE-FE coupled method.Utsunomya et al. (1995) reported the comparison of the results between experimental and theoretical analyses on a mat-type large floating structure.An experimental modal analysis for estimating the dynamic properties of units linked large floating structure model was presented by Endo et al. (1996).Wang and Tay (2010) presented the mathematical formulation for the hydroelastic analysis of pontoon-type VLFS in the frequency domain framework.They employed the coupled finiteelement boundary-element (FE-BE) method to solve the coupled water-plate problem.In order to decouple the water-plate interaction problem, they adopted the modal expansion method.Yoon et al. (2014) proposed a direct coupling method to analyze floating plate structures with multiple hinge connections in regular waves.Their method employed finite element method (FEM) and boundary element method (BEM) to discretize floating plates and surrounding fluid, respectively.They investigated the maximum bending moment and deflection in the plate structures.
External unsteady dynamic loads can act on the very large floating airport during landing and takeoff, for example.Such loads will induce the transient behavior of the structure and may affect the safety of landing and taking off.Therefore, the transient hydrodynamic response of VLFS must be clearly clarified.
Some numerical schemes for transient hydrodynamic responses have been treated to date.Ohmatsu (1998) analyzed the transient response of VLFS indirectly by using impulse response function, which can be determined from the Fourier inverse transform of the frequency response function.Endo and Yago (1999) carried out a series of model tests in which the pontoon-type model is subjected to moving load and vertical impact load which idealizes the airplane landing and takeoff.They also developed a time-domain analysis method in which the structure is modeled by many segmented panels and hydrodynamic effect is evaluated using the memory effect function.In their study, hydrodynamic coefficients were calculated element-wise.Kashiwagi (2000) developed a mode superposition procedure to calculate transient elastic responses of VLFS. Lee et al. (2001) analyzed the transient hydrodynamic response of VLFS based on the mode superposition method.Lee and Choi (2003) developed a (FE-BE) hybrid method to analyze transient hydroelastic response of VLFS.Kashiwagi (2004) used his mode superposition method to numerically simulate the transient responses of a floating airport during landing and takeoff by an airplane using realistic data from a Boeing 747-400 jumbo jet.Qui and Liu (2007) proposed three-dimensional time-domain finite element procedure to analyze the transient hydroelastic responses of VLFS subjected to unsteady external loading.Qui (2009) presented a fully coupled time-domain finite element method for analysis of the dynamic response of free-free flexible beam floating in an unbounded water domain under the effect of moving loads.Endo (2000) applied his time-domain calculation method to a simulation of airplane landing and takeoff taking into consideration the combined action of incident wave.Ismail Latin American Journal of Solids and Structures 13 (2016) 1340-1359 (2000) proposed a direct time-domain method to analyze the dynamic behavior of flexible floating structure under moving loads.In his method, beam finite elements and two dimension constant boundary elements are employed to discretize the floating structure and the surrounding fluid, respectively.Newmark's method was employed to solve the equation of motion.He considered two types of flexible floating structures, namely: flat type and multiple module-linked bridge structures and studied the effect of stiffness flexibility of the former type and the connector flexibility of the latter type on the dynamic behavior of each type.Cheng et al. (2014) proposed a direct time domain modal expansion method to compute the transient behavior of VLFS subjected to simultaneously to incident wave and external loads including weight drop and landing/or takeoff load of an aircraft.
This paper presents a time-domain three dimensional (3D) numerical procedure to simulate the transient response of floating airport subjected to arbitrary unsteady external loads.The proposed method employs the Boundary Element and Finite Element methods (FEM-BEM).The floating structure is discretized by 4-node isoparametric finite element, which is based on Mindlin's plate theory, whereas the 3D surrounding fluid domain boundaries are discretized by four-node constant boundary elements.The equation of motion takes into consideration the effect of inertia loading due to the moving mass.In order to obtain the hydrodynamic forces (added mass and radiation damping), the coupled natural frequencies are first obtained by an iterative method, since hydrodynamic forces become frequency-dependent.Then the Newark integration method is employed to solve the equation of motion for structural system.Based on the presented procedure, a FORTRAN program is developed and its results are verified by comparison with published experimental results of a scale model of VLFS under a weight drop, airplane landing, and takeoff in still water condition.

ANALYTICAL MODEL
The three dimensional floating structure, the coordinate system and the definitions of the fluid domain and boundaries, as shown in Figure 1, are considered in the proposed model.In the present study, the floating structure is assumed to be a plate structure.The plate structure resists exerting forces with bending and torsion moments while the in-plane forces are considered to be negligible.The structure is of constant thickness and linearly elastic.The structure is vertically guided without friction along the sides.The water depth is constant, and is defined as D. The fluid is non-viscous, non-compressible and the fluid motion is irrotational.

Three Dimension BE Model for Fluid Domain
According to the linear wave theory, the fluid motion can be represented by the velocity potential which can be written as where ϕ is the velocity potential of radiated wave, ω is circular frequency, t is the time and i=(- Latin American Journal of Solids and Structures 13 (2016) 1340-1359 Figure 1: The floating structure and coordinate system.
The flow velocity, V  , can be expressed as: where  is the gradient operator.The displacement function of structural motion can be expressed as where s is the vertical displacement amplitude.When a floating structure is vibrating with harmonic motion in the still water, the radiation wave potential can be obtained by solving the Laplace's equation: in which Ω is the fluid domain.The boundary conditions at the sea-bed, at sommerfeld's radiation condition, Dean and Darlrymple (1984), at the free water surface and at the structure-water interface are respectively: Latin American Journal of Solids and Structures 13 (2016) 1340-1359 in which r=(x 2 +y 2 +z 2 ) 1/2 is the distance from the origin, to the point to be considered, g is the gravitational acceleration, and k is the wave number given as follows, where D is the water depth.
The above boundary value problem can be solved by Boundary Element Method (Brebbia, 1989).For this case, weighted residual expression can be written as , Γ=Γ ϕ +Γq, n is the unit normal vector and the terms with bars, f andq , represent known values of the function and its derivatives.
The inverse formula for Equation ( 6) is derived as follows, The fundamental solution ϕ * which is satisfying Equation ( 8) is adopted solution as weighted function, where δ(P,Q) is Dirac delta function.Consequently, by using relation Equation ( 8), the integral equation ( 7) is rewritten as follows, in which fundamental solution ϕ * (P,Q) for three dimensional potential problem is given as: Latin American Journal of Solids and Structures 13 (2016) 1340-1359 where R= [(xP-xQ) 2 +(yP-yQ) 2+ (zP-zQ) 2 ] 1/2 is the distance from the point P(xP,yP,zP) of application of the delta function to any point under consideration, Q (xQ,yQ,zQ).
Incorporating Equation (10) into Equation ( 9), integral equation on boundary may be derived as follows.

{ }
in which c(P) is constant relating boundary Γ.For a smooth boundary, c(P) =1/2.The boundary is assumed to be divided into N constant quadrilateral elements.The values of ϕ and q are assumed to be constant over each element and equal to the value at the mid-element node.Then Equation ( 11) reduces to the following representation before applying any boundary conditions, in which Γj is the boundary of th j element.Now, collocation method is adopted to Equation ( 12) for each element's mid-point Pi ( i=1, 2,... N) Solids and Structures 13 (2016) 1340-1359 The integrals can be calculated by using a 4-points Gauss quadrature role for all elements except the one corresponding to the node under consideration.For this particular case, H ii is zero due to the orthogonality of R and n.The boundary integral G ii can be calculated as where ΔAi and b are the area and aspect ratio of the i-th element, respectively.and The sea bed boundary, Γb, radiation boundary, G ¥ , free surface boundary, Γf , and structural boundaries, Γsh and Γsv , are numbered as 1, 2, 3, 4 and 5, respectively.Substituting boundary conditions from equation (4-a) to equation (4-f), into equation ( 13), equation ( 13) can be written in matrix expressions as follows: where, in which the subscripts indicate the belonging boundaries.The complex matrix, equation ( 15), is solved for the radiation wave potential on each boundary elements, {ϕ}.

FE Model for Structure
The structure is discretized by 4-node isoparametric finite elements on the basis of Mindlin's plate theory.The linear equation of motion of structure is given by in which [Bb] (3x12) and [Bs] (3x12) are the strain-displacement relationship matrices for bending and shear deformation, respectively, [Db] (3x3) and [Ds] (3x3) and are the stress-strain relationship matrices for bending and shear deformation, respectively, [kb] (12x12) , [ks] (12x12), are bending and shear matrices, respectively and J is Jacobian.To avoid the shear locking phenomenon, a selective reduced integration technique is used.The (2x2) and (1x1) Gauss quadrature is used for bending and shear matrices, respectively.Six Gauss points are used through plate thickness.The element mass matrix[m] e (12x12) is given by where, h is the plate thickness, ρst is the mass density of the structure, [N] is the matrix of shape function.The matrix of shape function used here is defined as: Owen and Hinton (1980).
where N1, N2, N3,and N4 have the following expressions: Latin American Journal of Solids and Structures 13 (2016) 1340-1359 The time-dependent external force vector {F(x,y,t)ext} g can be impulsive and/or moving loads.In the case of moving load, there are two models frequently used in the literature (Wu et al., 2000;Esen, 2013): the moving force model and the moving mass model.In the former model, the force is assumed to be constant and equal to the weight of the moving body, where {F(x,y,t)ext} g is the net force applied by the moving mass on x,y point at time t.δ(x-v.t)and δ(y-e) represents the Dirac delta functions in x and y directions, respectively, e is the location of the moving mass in y direction, M is the moving mass and v is the moving mass velocity.This model is called inertia-free loading.In the moving mass model, the inertia forces arising due to combined oscillation of this load and the beam (inertial loading) should be considered.In this case, the global external load vector is expressed as The element nodal forces due to moving mass over the element can be expressed as

Evaluation of Hydrodynamic Force
Since the plate finite element and constant boundary element are used, structure and fluid nodes do not coincide at this stage (Hamamoto and Fujita, 1995).To allocate the fluid force to interface nodes of each finite element, the hydrodynamic pressure on each interface boundary element is related to the plate element nodes by making use of shape function [N] of plate element.The consistent fluid force vector at connecting nodes of each plate element can be obtained as where Pz is the hydrodynamic pressure acting on the central point of each interface boundary element given as where ρω is the mass density of sea water.g can be modified to where α and β are called Rayleigh coefficients and given by Latin American Journal of Solids and Structures 13 (2016) 1340-1359 where ωi is the i-th dry mode frequency and ηi is the damping ratio.In this study, it is assumed that η1= η2= η.By disregarding the damping terms and external force vector in equation 35, the undamped free vibration of floating structure may be governed by The frequency equation of Equation 38 may be given by where w is the wet-mode circular frequency of the structure.As the added mass matrix is dependent on natural frequency, an iterative solution is carried out to obtain the coupled circular natural frequencies and mode shapes.The initial values are assumed by setting [Mw] g = [0] in equation 39.
As a result, the element added mass and radiation damping matrices could be obtained by using w in place of w .The equation of motion, Equation ( 35), can be written at time t+Δt as: , , ( ) The above equation is a second order differential equation with respect to time, and may be solved by the finite difference method.Newmark's δ method is adopted here with choosing the value of δ equal to 1/4.This value corresponds to assume constant variation of acceleration during the time increment Δt.With the aid of this method, the velocity and the acceleration at time t+Δt can be expressed in terms of time t as: Latin American Journal of Solids and Structures 13 (2016) 1340-1359 Equation ( 42) is the simultaneous linear equations, and can be solved easily to obtain the structural displacement vector {uz (t+Δt)} g .Then, the structural velocity and acceleration vectors at time t+Δt are calculated by using equation ( 41).

NUMERICAL SIMULATIONS AND RESULTS
In order to validate the proposed time-domain FE-BE hybrid method for analysis of transient elastic responses of a floating structure under unsteady external loads, a scale model of Mega-Float is considered with length (L) of 9.75 m, breadth (B) of 1.95 m, thickness (h) of 0.053 m, draft (d) of 0.0163 m, bending stiffness (EI) of 17530 N.m 2 , and Poisson ratio of 0.3.Based on this scale model, Endo and Yago (1999) carried out experiments under three loading conditions-weight drop test, weight pull-up test, and weight moving test which idealize the airplane landing and takeoff.The tested model called VL-10 in which the plate is floated on the free surface of a towing tank of length 40 m, width 27.5 m and water depth 1.9 m.In these experiments, the vertical displacement at Z1-Z9 points (equally spaced along the central line) as indicated in Figure 1 were measured.For numerical calculations, the dimensions of fluid domain are taken as 30x3x1.9m, the mass density of water (ρw) is 1x10 3 (kg/m 3 ), and mass density of the structure (ρst) is 315.6 (kg/m 3 ).The floating plate is discretized by 4-node Mindlin's plate elements with selective reduced integration scheme.It should be noted that when a law order Mindlin's plate bending element with selective reduce integration technique is used to discretize the plate structure resting on elastic base, the hourglass mode may occur.To avoid this mode to be occurred, two elements with full integration (2x2) is used.Suzuki and Yoshida (1996) proposed a characteristic length λc, Eq. ( 44), as a rational measure to distinguish VLFS from the conventional ship and floating offshore structures in terms of global response.
Latin American Journal of Solids and Structures 13 (2016) 1340-1359 where kc is the spring constant of hydrostatic restoring force.λc, corresponds to the length of locally deflected region by a static concentrated load.This indicates that the influence of an applied load on the elastic deformation is limited within the region of the length λc.Accordingly, if the length of structure is smaller than the characteristic length, the response is dominated by rigid-body motions, whereas if it is larger than the characteristic length, as typically in VLFS, the response is dominated by elastic deformations.For the tested model, λc = 6.14 m which is less than the plate length, 9.75m.Thus, the response is dominated by elastic deformations.

Simulation for the Weight Drop Test
The weight drop test is considered here as a typical impulsive experiment.In this test, a weight of 196 N was dropped from a height of 0.12 m onto the "Impact point" indicated in Figure 1 (15 cm a head of point Z2).The detected acceleration of the weight during the impact was recorded, and the result is shown in dimensionless form as the ratio of the gravitational acceleration g in Figure 2. Therefore, this dimensionless acceleration multiplied by 196 N is regarded as the impact load.
As the numerical results may be influenced by the finite element mesh, constant damping ration, η, time increment, Δt, the convergence analysis has to be carried out to obtain best comparison with the experimental data.The influence of each factor on the accuracy of the numerical simulation will be studied in the following sub-sections.

Effect of Finite Element Mesh
In order to carry out the convergence analysis of finite element mesh, four finite element meshes are used to discretize the floating structure which are 8x6, 16x6, 32x6, and 64x6 elements.Six elements are used in the width directions while the length is divided by various number of elements.Figure 3 shows the finite element mesh used in discretizing the structure, in the case of 32x6 elements, together with the boundary element mesh used in discretizing the boundaries of water domain.In all considered models, a constant damping ratio, η, is assumed as 0.005 and time step, Δt, is taken as 0.002 s. Figure 4 shows the comparison between the experimental data of the vertical displacement time history at point Z1 and the numerical results obtained for various finite element meshes.It can be noticed that the number of elements per structure length (or per characteristic length) significantly influences the accuracy of the numerical results where the accuracy improves with the increase of number of elements.Further, it is observed that when number of elements increased from 32 to 64, there is no further improvement in the accuracy.Therefore, the mesh 32x6 will be used in the subsequent analysis.

Effect of Structure Damping Ratio
In order to quantify the structure damping ratio, η, three values of damping ratio are considered which are 0.0, 0.005, and 0.05.In all considered models, the plate is discretized by 32x6 (4-node) Mindlin's plate elements with selective reduced integration and time step, Δt, is taken as 0.002 s. be noticed that the damping ratio slightly influences the accuracy of the numerical results though the case with damping ratio of 0.005 gives better numerical results in a comparison with experimental data.Therefore, a value of 0.005 will be assumed in the subsequent simulations.

Effect of Time Step
In order to quantify the time step, Δt, three values of time step are considered which are 0.002s, 0.005s, and 0.01s.In all considered models, the plate is discretized by 32x6 (4-node) Mindlin's plate elements with selective reduced integration and the structure damping ratio, η, is assumed as 0.005.Figure 6 shows the comparison between the experimental data of the vertical displacement time history at point Z1 and the numerical results obtained for various values of time step.It can be noticed that the time step doesn't influence the accuracy of the numerical results as long as time step is taken less than 0.01s.Therefore, a value of 0.002 will be taken in the subsequent simulations.From the above convergence analysis, the plate is discretized by 32x6 (4-node) Mindlin's plate elements with selective reduced integration.A constant damping ration, η, is assumed as 0.005 and time step, Δt, is taken as 0.002 s.The comparisons of vertical displacement time history at five measured points with experimental data are indicated in Figure 7.It can be seen from this figure that the present method can predict the global behavior reasonable well.By comparison with the measurements, one can notice that the displacements near the impact point, Z2 and Z1 are larger than the numerical results.This may be attributed to the difference in the effective area of the impact load between the experiment (the effective area is infinite) and the numerical simulation (the load is applied at one node which cause the deformation to be localized at the impact point).As it can be observed, the maximum displacement is found at the end z1 and at the point Z2 which are near to the impact point.This can be attributed to the local elastic deformations.Also, one can notice that the transient phenomena can be seen even at time t=1.2s.However, the deformation at t=1.8 s is almost the same as that in static equilibrium   above convergence analysis, the plate is discretized by 32x6 4-node Mindlin's plate elements with selective reduced integration, a constant damping ratio, η, is assumed as 0.005 and the time step, Δt, is taken as 0.002 s.The comparisons of vertical displacement time history at five measured points with experimental data are shown in Figure 8.As observed in the figure, the calculated results are found to show a good agreement with the experiment.Also, the figure tells us that at almost same time (t= 0.8s) the points Z1, Z2, and Z3 deflect up while the center point Z5 deflects down.However, the deformation at t=2.5 s is almost the same as that in static equilibrium.The end point Z9 raise up at t=1.s before it sinks at t=1.8s though the magnitude of its deflection is small.These findings may be due to that the response is dominated by the elastic deformations.

Simulation for the Weight Moving Test
In order to simulate airplane motion, the weight moving test is considered.In this test, a small carriage of weight 67.7 N was towed on the rail plotted in figure 1.The rail is shifted 18 cm from the centerline of the model.The towed speed is taken to be 0.61 m/s and the acceleration is ignored.And the starting point is 61 cm aft from z7 point.The plate is discretized by 32x6 4-node Mindlin's plate elements with selective reduced integration.In order to avoid the hourglass mode, two element with full integration are used.A constant damping ration, η, is assumed as 0.005 and the time step, Δt, is taken as 0.002s.The comparisons of vertical deflection time series at five measured points with experimental data are shown in Figure 9.As observed from the figure, the whole features are similar to the experimental data.Further, it can be seen that, while the point Z7 deflects down the end point Z9 deflects up.Also the other end point Z1 deflects up while the point Z3 deflects down.These can be attributed to the fact that the response is dominated by elastic deformations.The maximum deflection is found at Z3 and equal to 1.5 mm which is very small compared to the plate length.1) The number of elements per characteristic length significantly influences the accuracy of the numerical results.Therefore, a convergence analysis for the element mesh should be carried out at a preliminary stage to identify the size of the finite element to be used in the final simulation.
2) The structure damping ratio has slight influence on the response of the studied floating plate.
3) The time step has insignificant influence on the response as long as its value is taken less than 0.01s.4) The numerical results of the proposed procedure are in very good agreement with published experimental results on a rectangular floating structure subjected to external transient loads as well as moving mass.5) The present method can be extended to study the dynamic response of floating structures with various shapes and end conditions under moving vehicles and/or incident waves.Further study can be carried out to investigate the effect of various parameters on the hydrodynamic behavior of such structures.

Figure 2 :
Figure 2: Detected acceleration of weight during impact onto plate.

Figure 3 :
Figure 3: FE Mesh used for discretizing the structure and BE mesh used for discretizing the bounries of fluid domain.

Figure 4 :
Figure 4: Influence of finite element mesh.
Figure 5 shows the comparison between the experimental data of the vertical displacement time history at point Z1 and the numerical results obtained for various values of damping ratio, η.It can Latin American Journal of Solids and Structures 13 (2016) 1340-1359

Figure 6 :
Figure 6: Influence of time step.

Figure 7 :
Figure 7: Time history of vertical displacement in weight drop test.

3. 2
Simulation for the Weight Pull-Up Test The weight pull-up test was carried out by removing a weight of 196 N from the "hit point" indicated in Figure 1.The time delay of the pull-up test is considered in numerical simulations.That is, the external load is supposed to decrease linearly from zero to -196 N during 0.20 s.Based on the Latin American Journal of Solids and Structures 13 (2016) 1340-1359

Figure 8 :
Figure 8: Time history of vertical displacement in weight pull-up test.

Figure 9 :
Figure 9: Time history of vertical displacement in weight moving test.
, a time-domain three dimensional Boundary Element-Finite Element procedure is presented to simulate the transient response of Mega-Float subjected to external transient load as well as moving mass.In the proposed procedure, the floating structure and the surrounding fluid are discretized by 4-node isoparametric finite elements and by 4-nodes constant boundary elements, respectively.Structural analysis is based on Mindlin's plate theory.A FORTRAN program is then developed and its accuracy and efficiency are verified by comparing its results with published experimental results of a scaled model of Mega-Float.Three tests, weight drop test, weight pull-up test, and weight moving test which idealize the airplane landing and takeoff, are numerically simulated.First, the convergence analysis is carried out to study the influence of finite element mesh, structure damping ratio, and the time step.Within this study, conclusions are summarized as follows: Latin American Journal of Solids and Structures 13 (2016) 1340-1359 , [cw] e , and [kw] e are the element matrices of hydrodynamic added mass, hydrodynamic radiation damping, and hydrostatic stiffness, respectively, and are given by Latin American Journal ofSolids and Structures 13 (2016) 1340-1359The radiation wave potential ϕ may be expressed as e