{'fraction_non_alphanumeric': 0.02736448272861742, 'fraction_numerical': 0.009484033759679805, 'mean_word_length': 4.606585365853658, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Seismic Transient Simulation of an Operating Wind Turbine Considering the Soil-Structure Interaction

In this contribution we present a strategy to investigate the vibrations of onshore wind turbines subjected to simultaneous aerodynamic loads and seismic loads. The latter result from the propagation of seismic waves through the underlying soil and are governed by the Lamé’s equations. The structure and its foundation as well as parts of the soil are modeled by Finite Element Method (FEM). The infinite half-space is discretized by Scaled Boundary Finite Element Method (SBFEM). Both methods are coupled at the common interface. The seismic excitation is expressed as a 3D seismic wave field and transformed into boundary tractions, which are then applied at the interface between the near and far fields. The aerodynamic actions are generated with an additional aerodynamic tool and are applied at the tower head. The proposed method can be used also to investigate the 3D nonlinear response of the near field, where nonlinear material properties can be assigned to any element of FEM part.


Introduction
Recent standard codes (e.g. [7]) introduced the seismic load combination for the design of wind turbines, which considers a reasonable likelihood of earthquake occurrence during the operational state or an emergency shutdown. A large part of the south European coastal areas presents high seismic hazard and wind conditions, that are sufficiently suitable for financial returns from modern wind turbines. In these areas, the seismic load combination may govern the tower and the foundation design. It is important to consider the interaction between the structure and the underlying soil for achieving a realistic estimation of the seismic response of wind turbines [19,21]. The simulation of the Soil-Structure-Interaction (SSI) issues for operating wind turbines involves two different mechanical problems: on the one hand, the vibrations due to aerodynamic loads excite the structure of the wind turbine and its foundation and, on the other hand, the vibrations are transferred into the surrounding infinite half-space, a process which is governed by the Lamé's equations. For such a complex problem, no analytical or semi-analytical solution is available, so that numerical models are used. For seismic analysis, an important aspect is the simulation of the loading generated by the earthquake. The seismic waves originate in the far field and travel towards the soil-structure interface, where a first part of the energy passes to the structure and a second part is reflected backwards. The wave pattern is further complicated by the presence of layers (or other inhomogeneities), which cause scattering and diffraction.
The aim of this study is to present a numerical model for the simulation of turbine-soil systems, subjected to simultaneous earthquake and aerodynamic loading. The wind turbine can be modeled by using standard methods such as the Finite Element Method (FEM). In contrast, the infinite soil cannot be represented by means of discrete models, as they do not satisfy the radiation condition for which the wave energy must dissipate to infinity. Therefore, we proposed a method based on the coupling between the FEM and the Scaled Boundary Finite Element Method (SBFEM). A detailed derivation of the SBFEM and its coupling to the FEM can be found in Schauer [13].
Most of the studies about the FEM/SBFEM considered only external loads applied to the near field. Few research works presented a method for including both external loads and the earthquake actions caused by the propagation of scattered seismic waves.
Seiphoori et al. [16] studied the nonlinear seismic analysis of concrete faced rockfill dam, obtaining the scattered motion and interaction forces along the canyon using the SBFEM. The dam was subjected to spatially variable P, SV, and SH waves, in order to evaluate the effect of dam-foundation interaction and the reservoir water effects.
Also Bazyar and Song [2] transformed seismic waves into surface tractions at the interface between the near and far field system. They considered incident fields of obliquely plane waves coming from the far field, where the seismic wave inputs are formulated as boundary tractions applied to the near field.
Syed and Maheshwari [17] presented a rigorous algorithm for the computation of ground interaction forces at the interface between near and far field. These are obtained solving a convolution integral between the acceleration unit-impulse matrix and the relative accelerations, which are obtained as the difference between the seismic scattering accelerations and the absolute accelerations at the near-far field interaction nodes. Unfortunately, this rigorous approach with non-locality in both time and space (due to the convolution) results in significant computation effort. For this reason, Syed and Maheshwari proposed an improvement of the seismic FEM/SBFEM [18], where two different approximation techniques -one in time and one in space -helped to reduce the computational time to only 5% of that required using the conventional method.
In the present work, the FEM/SBFEM coupling after Schauer [13] is enhanced with a function for computing the seismic boundary tractions, which simulate the three-dimensional scattered earthquake loading. Efficiency and feasibility of the method are key aspects for the seismic simulation of turbine-soil systems. Therefore, the new feature for the computation of the seismic forces was conceived as independent as possible from the existing "'non-seismic"' algorithm and it is an elegant approximation, based on simple principles of wave propagation.

Numerical Method
For the seismic simulation of the operating wind turbine we use a sub-structuring method, where the problem is subdivided into two sub-systems (cf. fig.1). The structure and its foundation as well as parts of the soil (the near field) are modeled by FEM. The infinite half-space (the far field) is simulated with the SBFEM. The sub-systems are coupled at the common interface Γ [22].

FEM-SBFEM Coupling
The FEM-SBFEM coupling is described briefly, more information about the theoretical background and its efficient implementation can be found in Schauer [13] [14]. The displacementbased FEM in time domain can be derived by the energy theorem here the vectors u, du dt =u and d 2 u dt 2 =ü represent displacement, velocity and acceleration, respectively. κ and ρ denote damping ratio and density. Applied tractions and forces are given by t and f . After inserting the strain displacement relation = Du and the stress strain relation σ = E in eq. (1) and introducing shape functions the given energy theorem can be rewritten in the matrix-vector form: where M , C and K denotes mass, damping and stiffness matrix, respectively. p ΩΩ are the forces applied to the structure (such as the aerodynamic loads) and to the near field; the p ΓΓ are the forces applied at the interface between near and far field (such as the seismic loads). The damping properties are introduced according to the Rayleigh representation, so that the damping matrix is a combination of weighted mass matrix and stiffness matrix The influence of the infinite half-space is described by SBFEM. The interacting forces at the interface Γ are given by the vector here γ is a parameter introduced by the time integration scheme. The unit acceleration impulse matrices M ∞ are assumed to be constant within one time step ∆t and they are computed in a pre-process.
The vector of eq. (5) can be added to the right hand side of eq. (3) in order to get a direct and bidirectional coupled FEM-SBFEM formulation The computation in time domain is conducted by executing the generalized-α time stepping algorithm [5].

Seismic loads
In seismic analyses, the rigorous procedure to compute the interaction forces at the interface Γ would require to use the relative quantities in Eq. 5 in place ofü j andu j . However, this requires a modification of the existing algorithm and a higher computation effort. As an efficient approximation, we consider the seismic forces as external loads contained in the vector p ΓΓ in Eq. 6. We assume that the seismic wave propagates through the far field at a constant speed horizontally or inclined with a certain angle. A first function computes the distribution of the accelerations due to the propagating wave front at each time step in the far field, using the SBFEM on its own. It is possible to input arbitrary time-dependent accelerograms for the wave front. The wave front reaches the near field and accelerates the interactions nodes at different time steps depending on the position of the nodes, which are located on the coupling surface shaped as an half sphere (see Fig. 2). Only the accelerations applied to the outer nodes are considered and the wave front is transferred homogeneously to the near field. In Fig. 2, the nodes in the shadow areas are not loaded during the wave propagation process to avoid inhomogeneous effects.
A second function translates the accelerations at the outer nodes at each time step n into loads, which are compatible with the FEM/SBFEM coupling and are saved into the vector p ΓΓ . The loads can be computed as: where M is the mass matrix, as defined before. 1 is a vector that contains '1' at the degrees of freedom accelerated by the wave front and a(t) contains the scalar information of acceleration.

Demonstration based on the propagation of an Heaviside-like wave front in a free field
A horizontal wave front propagates through a free field (a homogeneous half space without the structure) and the time-dependent amplitude of the accelerations varies as an Heaviside function with a maximum acceleration of 1 ms −2 (see Fig. 3a). The front propagates from the negative towards the positive x-coordinates (from −200 m to 200 m). The free field is modeled as a homogeneous isotropic half space with the material parameters: Young's modulus E = 1.173 ·10 8 N/m 2 , Poison's ratio ν = 0.28 and density ρ = 1937 kg/m 3 , and therefore a P-wave velocity of c p = 246.12 m/s. Fig. 3b shows the discretization of the free field, made up of the near field modeled as a half sphere of radius 200 m with 3438 nodes and 324 elements and the far field modeled an an excavated half space. The far field is evaluated only at the nodes of the interaction nodes on the surface of the half sphere, where the two sub-systems are coupled.
At the time 0 s, the wave front reaches the outermost node at x = −200 m and gradually it accelerates all the other outer nodes. These boundary accelerations are transformed into boundary loads for the near field and the displacements at the free field surface are evaluated at different x-coordinates. Fig. 3c and Fig. 3d shows the displacements at x = 0 m and x = 75 m respectively. Fig. 4 shows the three dimensional visualization of the wave propagation in the near field.

Aerodynamic loads
The estimation of the aerodynamic loads is carried out according to [20]. It is a common approach to uncouple the tower-foundation structure and the rotor-nacelle system. From a structural point of view, for a specific wind speed and turbulence class, the rotor-nacelle system transforms the turbulent wind (the input) into tower top loads (the output). The output of the rotor-nacelle system analysis, that is the tower top aerodynamic loads, become the input for the second sub-system, which is the tower-foundation-soil system. Deriving the aerodynamic loads independently of the seismic motion implies two main assumptions: • the aerodynamic loads are not noticeably influenced by the SSI and can be computed for a fixed base tower; • the system remains linear and the superposition principle is applicable.
The validity of the first assumption depends on the focus of the investigation and the magnitude of the seismic load. Here, the focus of the investigations lies on the demonstration of the strategy for 3D seismic excitations rather than on the aerodynamics. The aerodynamic loads are computed with FAST [8], and then applied to the beam head.

Implementation
To analyze the seismic behavior of wind turbines different third-party software packages as well as different in-house computer codes are combined. The general work flow is shown in Fig. 5. Finite element meshes and scaled boundary finite element meshes are generated by using the geometry and mesh generation toolkit CUBIT [11]. Here the geometry is created as well as the final meshes. To conduct the analysis of the full model of wind turbine, its foundation and the surrounding soil, preliminary steps are required. At first the eigenvalue problem is solved using GNU Octave [6] to compute the eigenvalues and eigenvectors of the wind turbine and its foundation. From the eigenvalues of the system the damping coefficients for the coupled problem are derived. Secondly the response of the infinite half-space is computed, the M ∞ matrices are stored in compressed binary format. The aerodynamic loads are generated by using FAST [10] from National Renewable Energy Laboratory. The process, which solved the coupled FEM-SBFEM computation, has accesses to all previously generated files. The FE-model is build up and the information of the SBFE-model is used to determine the coupling nodes in order to apply the influence of the infinite half-space to the boundary of the FE-mesh. The aerodynamic loads acting on the rotor and nacelle are applied to the tower top. The earthquake loads are generated and applied to the near field according to Sec. 2.2.
Even if the seismic forces are computed with an efficient approximation technique, the SSI interaction forces p b in Eq. 5 are computed rigorously with a convolution integral, which is memory and time consuming. To handle the high memory consumption and the computational effort, model reduction techniques and parallel computing techniques are used [14]. Also with the same objectives in mind, a projection method for using non-matching meshes at the nearfield/far-field interface [15] will be introduced in future works.
Different third party libraries are used to achieve reasonable performance. The program which computes the infinite half-spaces response has to handle dense and sparse matrices as well, that is why ScaLAPACK [3] builds the main core of it. The finite element matrices are sparse and handled with PETSc [1]. At first all input files are parsed to set up the virtual model. This includes in case of the FEM: nodes, elements, material parameters, damping coefficients, boundary conditions, nodal load and nodal accelerations as well as initial conditions. With this information the description of the FEM part is complete. For the SBFEM, nodes and elements are processed and M ∞ matrices are read. Afterward, the interface is initialized and coupling nodes are defined. Then the system matrices M , K, and D are assembled. Here these matrices stay constant during the simulation, since linear elastic material behavior and only small displacements of the structure are assumed. The simulation is conducted with a time step length ∆t until the total time T is reached. Within this time loop, three major steps are processed. At first the response of the infinite half space is computed. The resulting forces are added to the force vector, which contains also the aerodynamic forces and the seismic forces. After solving the equation of motion, by using a generalized-α time stepping scheme (cf. [5]), the velocities at the interface are stored to use them to evaluate the future time steps.

Case Study
For the considered case study, the stiffness and mass properties distribution of the wind turbine properties are summarized in Tab. 1. The rotor blades,the nacelle and the gear box were idealized as a stubby rigid element of 1 m at the top of the tower. The tower is made of steel using the following material parameters: Young's modulus E = 2.1 · 10 11 N/m 2 , Poisson's ratio ν = 0.3 and density ρ = 7850 kg/m 3 . The wind turbine is founded on a monopile embedded in the soil for a length of 38 m and with the same properties of the tower.
The soil is assumed to be homogeneous and consists of sand with a shear modulus of G s = 7.64 · 10 7 N/m 2 , a Poisson's ratio of ν s = 0.28 and a density of ρ s = 1937 kg/m 3 . Control and regulation issues ensure an operational service within the desired and safe range of rotational speeds, inclinations of the blades to the wind and power outputs. These aspects were not considered as analysis parameters of this study and were not included in the model.

Model
The tower with the monopile is modeled as a cantilever Euler-Bernoulli beam. 11 elements, with bending and axial flexibility as well as torsion are used. Each nodes has six degree of freedom  Aus diesem Grund werden zwei Lagerungen betrachtet, um d im Boden ausreichend genau abzubilden und zu vergleiche stellen die zwei Lagerungsarten dar. � x � � y � � z � � Abbildung 6.11: Lagerung ϕ x = ϕ y = ϕ z = 0

Numerische Lösung der Windkraftanlage (Typ2)
Wie bereits beschrieben, wird der Pfahl zunächst mit den belastet (vgl. Abbildung 6.4), um Unterschiede aus der nac statischen Lagerung zu erfassen. In diesem Zusammenhang verläufe aus der Ersatzlastbeanspruchung infolge des horizon in 38,1 m Tiefe und bei 0 m dargestellt. Es werden beide ang betrachtet (vgl. Abbildungen 6.13 bis 6.14).   Fig. 6 shows the FEM near field including the monopile foundation. The far field is represented using scaled boundary finite elements. The nodes and elements are located at the common interface between the finite element (near field) and the scaled boundary finite element (far field). Since only the surface is discretized 97 nodes, with 3 DOF each, are needed. The scaling center is located in the domains center at the tower base with the global coordinates SC = {0, 0, 0}.

Loads
For the computation of the turbulent aerodynamic loads with FAST, a reference speed of 12 m/s was used. The turbulent wind field model used corresponds to the spectral and exponential coherence model of Kaimal, which also meets the requirements of the standards [4]. The resulting loads F x and F y act in fore-aft direction and side-side direction, respectively. F z acts in axial direction. The moments M x , M y , and M z are roll, pitch and yaw moments, respectively.
As explained in Sec. 2.2, a pre-step for the estimation of the seismic loads p ΓΓ at the interaction nodes is necessary. Here a strong earthquake event has been chosen to conduct the numerical analysis. Fig. 8 shows the acceleration-time-plot of the Kobe, Japan earthquake in 1995. It has a moment magnitude scale of M w = 7.2.

Damping
The presented model uses Rayleigh damping and the coefficients c M and c K are needed to build the damping matrix (cf. Eq. (4)). It is possible to compute c M and c K by evaluating the damping coefficient The terms c M 4πf i and c K πf i prescribe mass proportional damping and stiffness proportional damping, respectively. We considered the first two natural frequencies are f 1 = 0.714 Hz and   Fig. 9 shows the response of the tower subjected to aerodynamic loads with and without damping. In the further paragraphs, all the calculations are carried out for a damped system.

Results
The case study is conducted with a time step length of ∆t = 0.001 s for a total period of T = 360 s. First, the tower is subjected only to the seismic actions. Then, the tower is subjected to operating aerodynamic loads (as in Fig. 9) and, additionally, at the time t = 100 s the seismic action is triggered. All the plot shows the different components of the displacements of the tower.   Fig. 10 depicts the results of the towers top and base displacements in the x-, y-and z-direction, when the tower is subjected to a seismic excitation. The horizontal wave front propagates in the near field in x-direction and excites the tower mainly in x-direction. However, the propagation generates also a non-negligible vertical component at the tower base. The amplitude of the base displacements reaches an absolute maximum of 4 cm. The double integral of the accelerations in Fig. 8 can be interpreted as an approximation of the seismic root point excitation without considering the spatial propagation of the wave and it leads to a value of ca. 20 cm, which is 5 times bigger than the computed maximum base displacement. This is due to the three-dimensional effects considered when generating the seismic loads. The maximum top displacements in x-direction is about 15 cm.

Turbine under aerodynamic and seismic loads
At t = 100 s the seismic excitation is applied in form of tracions at the near field boundary. Fig. 11 shows the tower displacements. The additional seismic loads results, as expected, into an increase of the towers top displacement during the seismic event and leads to a maximum top displacement of ca. 37 cm. After the seismic event is subsided, the towers motions returns to values comparable to the aerodynamic response without seismic event, shown is Fig. 9 for the damped case.

Conclusions
In this article we presented a strategy to simulate operating wind turbine subjected to seismic excitation, considering the soil-structure-interaction, with the aid of a coupled FEM-SBFEM approach in time domain. The proposed method is able to represent the salient aspects of the SSI interaction and allows arbitrary transient loading conditions, such as aerodynamic loads, seismic loads and the combinations of them. The seismic loading is expressed as a 3D seismic wave field and transformed into boundary tractions, which are then applied at the interface between the near and far fields. This allows a scattered representation of the incoming seismic wave field and lead to more realistic simulations. Each simulation lasted few hours running on a personal computers. This strategy helps to compute more realistic design loads for the case of a simultaneous interplay of aerodynamic and seismic loads. It can be used to verify if the design loads computed with a standard combination of the sum of the maxima based on independent aerodynamic and seimsic simulations are reasonable or too conservative. This would likely have important implications on the turbine design. In future works, nonlinear effects shall be taken into account, thanks to the suitability of this strategy for transient time-domain analysis. The nonlinearities can be accounted for in the structure, in the near field and at the soil-structure interface. Concerning the soil (near field) nonlinearity, most of the models used in practice consider a simple one-dimensional behavior. The proposed method is able to investigate the 3D nonlinear response of the near field, as the nonlinearity can be assigned to any element of the model. More sophisticated FEM tower models, with openings or non-axial-symmetrical shapes including SSI can be easily investigated with the proposed strategy.