Seabed Dynamic Response of Offshore Wind Turbine Foundation under Vertical Harmonic Loading : An Analytic Solution

Wave-induced dynamic loads yield displacement and stress propagation on the seabed beneath the offshore wind turbine foundation. Existing attempts to calculate the displacement/stress propagation based on commercial software packages or numerical solvers usually ignore the inertial effects and provide quasi-static solution for the dynamic problemor ignore the damping effect of the soil. Since real time displacement measurements on the seabed are difficult to be carried out, the existing results are not able to be validated and the conclusions are rather arbitrary. Analytic solutions can provide accurate results and are less timeconsuming and more computationally efficient than the numerical approaches but difficult to be achieved. In the present study an exact solution of the wave propagation in the seabed under harmonic excitation of the wind turbine foundation is achieved. The solution is based on integral transforms and stress functions. Verification of the solution is provided and a representative example is solved and commended.


Introduction
Existing attempts to calculate the displacement and stress propagation on the seabed beneath the offshore wind turbine foundation are mainly based on commercial software packages or numerical solvers [1][2][3][4].In the most existing works, the solutions to the wave, structure, and seabed interaction problems are based on a 'one-way coupling' assumption (e.g., [1,3,4]).In [3][4][5], the Biot's equation [6,7] is adopted to model the soil behaviour.The inertial effects are ignored, yielding quasi-static solution for the dynamic displacement propagation in the soil.Moreover, the assumption of the elastic behaviour of the soil does not take into account the damping effect.Since real-time displacement measurements on the seabed beneath the foundation are difficult to be carried out, the existing numerical results and the conclusions are not validated.Also, the usage of the numerical tools can have a high computational cost and is time-consuming for solving the coupling problem in order to achieve a low residual.
Analytic solutions to such problems are more accurate and more efficient compared to the numerical ways but are difficult to achieve.A fundamental contact mechanics solution based on Green's functions and boundary integral equation solution has been achieved in [8] for the stress analysis of the soil-foundation interaction under static torsion.Analytical solutions for dissimilar or homogeneous elastic medium (space and half-space) containing interfacial cracks, under axial, radial, and torsional loading sources are provided in [9][10][11][12][13][14][15][16].Using double integral transforms, the dual integral equations for wave propagation of a cracked medium under transient dynamic loading are originally solved in [17].Analytical solutions for the dynamic response of foundation on elastic and viscoelastic soil are achieved in [18][19][20].Interaction of foundation piles embedded in transversely isotropic soil under dynamic loading is analysed in [21,22].Structural analysis of wind turbine tower under environmental loads is performed in [23], and analytical solution for evaluating the dynamic response of plates under underwater transient loads is presented in [24].In the present study, an exact solution of the wave propagation in the seabed under vertical harmonic excitation of the wind turbine foundation is achieved.The solution is based on integral transforms and stress functions.A representative example is solved and commended in the present work.Verification of the solution is provided by solving a specific case of a sudden point load on the solid surface and the solution is verified against the existing solution in [25].The proposed solution has two benefits: (a) it takes into account the inertial effects in the displacement propagation, and (b) the solution is analytic, yielding better accuracy and negligible computation cost.The limitation of the present version is that it does not include the effect of porous material and damping yet.The analytic solution in the present work can be extended for solving porous media with time-varying dynamic loads taking into account the material damping effect.

Formulation of the Problem
A circular foundation based on an elastic half-space (e.g., rock seabed) is subjected to vertical harmonic load.The coordinate system of the foundation is shown in Figure 1 Since no body-forces exist, the governing equation in terms of displacements (Navier's equation) can be written [26] as where  →  is the vector of the displacement components   ,   ,   in the directions , ,  respectively, ,  are the Lame constants, and  is the material density.Taking into account the axisymmetric geometry and load, the displacement component   is zero.Therefore, the displacement vector can be expressed by the following equation: where  →   ,  →   are the unit vectors in the directions ,  respectively.In terms of scalar and vector potentials ,  → ℎ , (3) can be written [25] as With the aid of ( 4), (2) yields In above equations the parameters  1 ,  2 are given in [25] by the formulae and express the propagation velocity of the dilatational and rotational waves, respectively.The operator ∇ 2 in polar coordinates is given by the following equation: If we use [25] a function , where (6) can further be simplified: With the aid of ( 4), ( 5), (10), and ( 11), the stresses (in polar coordinates) of Appendix A and the displacements can now be written in terms of the scalar functions , .

Solution of the Dynamic Problem
For the solution of the problem, Laplace and Hankel integral transforms will be applied on the governing equations and boundary conditions.The definition and some important properties of the above transforms are given below.

. . Laplace Transform of Any Function
The inverse form of Laplace transform is given by the following formula: . .Hankel Transform of Any Function () The inverse form of Hankel transform is given by the following formula: Taking the Laplace transform of the displacement and stress components given by ( 12)-( 20), formulae (B.1)-(B.7) in Appendix B will be obtained.
. .Boundary Conditions.The surface  = 0 of the half-space is loaded by the stress   (, 0, ) given by (1).Since the surface is free of shear stresses, the following condition is valid: Taking into account the definition of the Hankel and inverse Hankel transform given by ( 23) and ( 24), ( 1) can be written in the form or With the aid of the properties [27] 0 { ⟨−⟩ ;  → } = 0 (29) (27) yields Taking the Laplace transform of above equation, the following formula can be obtained: As a consequence of the definition of the Heaviside function, the following equations can be obtained: Combination of ( 28) and (33) yields With the aid of ( 29), (32), (33), and (34), ( 27) can be written as Using the notation (35) can be written: According to (24), above equation yields Taking the Hankel transform of above equation, the following formula can be obtained: or Then, the Laplace transform of (40) yields Taking the Laplace and Hankel transform to the boundary condition given by ( 25), the following result can be obtained: Applying the Laplace transform in (17) and (18), (, , ) =  {  (, , ) ;  → } Taking the Hankel transform  0 {•;  → } to above equations, the following formulae can be obtained: where Combining ( 41) with ( 45) and ( 42) with (46), the following system in terms of Φ, Ψ can be obtained: . .Determination of the Functions , .In order to solve the system of ( 49) and (50), Laplace and Hankel transform will be applied to the governing equations ( 5) and (11).The Laplace transformation of ( 5) and ( 11) yields Taking the Hankel transform  0 {•;  → } to above equations, the following ordinary differential equations can be obtained: The solution of above equations is where For  → ∞ the functions Φ, Ψ will be unbounded because of the terms  1   and  2   , respectively.Therefore, the following conditions should be valid: Then, With the aid of above two equations, the system of ( 49) and (50), for  = 0, yields Combining ( 47) with ( 57) and ( 48) with (58), the following formulae can be obtained: With the aid of ( 57)-(67), the system of algebraic equations ( 63) and ( 64) can be written as yielding where With the aid of (69) and (70), the functions given by ( 61), (62) can now be obtained: Since all the displacement and stress components given by ( 12)-( 20) are expressed in terms of the functions ,  the inverse Hankel transform  −1 0 {•;  → } of above equation yields the solution of the problem.

. . Inversion of the Hankel Transform. Taking the inverse
Hankel transform  −1 0 {•;  → } in (71), the following result can be obtained: Taking into account (72), the Laplace transformed displacements given in Appendix B can now be written as Taking the inverse Laplace transform of above equations, the final solution for the displacements   (, , ),   (, , ) can be obtained:

Numerical Example and Verification of the Solution
. .Numerical Example.The derived solution given by ( 75) and ( 76) can now be implemented to a representative example.The dynamic load is a uniform pressure acting on a circular surface with radius .The load can be modeled by a harmonic function and can be represented in terms of the variables ,  by (1).The selected data for the proposed loading function are  = 1.0 ,  = 8.2 ,  = 0.94 /, and the material properties for the soil are  1 = 355 /,  2 = 309 /,  = 2.137 × 10 8 / 2 .For the inverse Laplace transforms of ( 75) and (76), the Stehfest [28] method is adopted.According to this method, the functions   (, , ),   (, , ) can be obtained by the following formula: where In above equation, N is a free even number representing the number of terms used in the summation.Values within the range 10 ≤  ≤ 24 yield sufficient accuracy (e.g., [29]).In the present study we adopted the value N=24.According to these findings, the displacement tends to reach the z-axis asymptotically due to the terms  − and  − in (74).
. .Verification of the Solution.The derived solution can be adjusted to provide the surface displacement distribution for the specific case of a suddenly applied point load given by the following equation: where () is the Dirac Delta function of the variable  and ⟨⟩ is the Heaviside step function of time .Solution for this specific case is available in [25]: Following the procedure of [14], the point load  is given by the following equation: Therefore, The above term for the uniform distributed load  should be placed in (73) and (74).On the other hand, the Laplace transform   2 +  2 =  {sin ;  → } (84) should be replaced in ( 73) and (74) by the Laplace transform: For the case of point load , the radius  of the circular surface of the equivalent uniform load  given by ( 83) is  → 0. Therefore, (73) and (74) yield Mathematical Problems in Engineering 7 (86) and (87) yield the known [25] (80) and ( 81).With the aid of the Laplace inversion technique proposed in [25], (87) yields the wave propagation   on the soil surface for the pulse point load of (79): For a unit load P=1, the displacement distributions at time instants t=0.0005s, 0.001s, 0.0015s, and 0.002s are demonstrated in Figures 3(a), 3(b), 3(c), and 3(d).

Conclusions and Discussion
In the present study, an exact solution of the vertical harmonic loading induced dynamic response in the elastic soil has been achieved.The present analytic solution can be utilized for solving the seabed response under vertical harmonic excitation of the wind turbine foundation in waves.
The solution is based on integral transforms and stress functions.Verification to the solution has been conducted and a case study has been provided.The achieved solution for harmonic (with respect to time) pressure has been implemented for suddenly applied point load.The derived displacement functions for this transient point loading source are in full agreement with the existing [25] solution.For t=0.0005 s (Figure 3

Figure 1 :
Figure 1: Circular foundation on elastic seabed: (a) physical model and (b) mechanical model.
. A pressure distribution   (, 0, ) =  [ ⟨⟩ −  ⟨ − ⟩] sin  (1) is assumed to act on the surface of the half-space.The symbol ⟨⟩ denotes the Heaviside function.The term [⟨⟩−⟨− ⟩] represents the uniform pressure distribution in the area 0 ≤  ≤  (Figure 1(b)), the term sin  is the harmonic function of time, and  is the amplitude of the load.The harmonic term sin  causes wave propagation in the domain of the elastic seabed.
(a)) high values for   displacement take place in the vicinity of the point r=0 (loading source location).The displacement peak is surrounded by a valley ring which takes place in the area 0.987<r<0.222m.For t=0.0015 s (Figure3(b)) the peak is flattened and expanded in the area 0.0<r<0.208m, and the valley ring is shallow.For t=0.0015 s (Figure3(c)) the propagating peak is more flat and wide and is expanded in the area 0.0<r<0.388m.For t=0.025 (Figure3(d)) the plotted area 0.0<r<0.

Figure 3 :Figure 4 𝑆 6 )
Figure 3: 3D displacement distribution of   due to transient point load for time instants.Implementation of the proposed solution to transient point load shows full agreement with the existing solution in [25].

Laplace Transform of the Displacement and Stress Components, Given by (12)-(20)
5 m is covered by the flat peak and no valley (in this area) is visible.The proposed method takes into account the inertial effects and does not require much computation power.Since the solution is analytic, the results are more accurate than numerical ones.The present work can be extended for solving viscoelastic porous media response under time-varying dynamic loads. B.