Earthquake dynamic response of large flexible multibody systems

Abstract. In the paper dynamics of large flexible structures imposed on earthquakes and high amplitude vibrations is regarded. Precise dynamic equations of flexible systems are the basis for reliable motion simulation and analysis of loading of the design scheme elements. Generalized Newton–Euler dynamic equations for rigid and flexible bodies are applied. The basement compulsory motion realized because of earthquake or wave propagation is presented in the dynamic equations as reonomic constraints. The dynamic equations, algebraic equations and reonomic constraints compile a system of differential algebraic equations which are transformed to a system of ordinary differential equations with respect to the generalized coordinates and the reactions due to the reonomic constraints. Examples of large flexible structures and wind power generator dynamic analysis are presented.


Introduction
Earthquake shaking affects forced motion of the structure basement.The ground motion is registered by so called strong-motion accelerographs and normally consists in three orthogonal components of the ground acceleration.The velocity and displacement of the ground are obtained integrating the data of the accelerogram.They also can be analyzed to obtain direct estimation of peak ground motion, duration shaking and frequency.To specific seismic zones of the earth the so called seismic Zone factor Z is assigned which correspond to a fraction of the earth acceleration.For example for seismic zone 1 Z = 0.075G; for zone 4 Z = 0.4G (G = 9.81 is the earth acceleration).In Fig. 1 the displacements (all measures are in SI UNITS) of North-South component of El Centro earthquake, California, 1940, as well as, the polynomial approximation (in red) of the numerical data is depicted.To provide guidelines and formulas that constitute minimum legal requirements for design of structures building codes are developed and in use within particular regions.These requirements are intended to achieve satisfactory performance of the structures subject of seismic excitation.Several organizations involved in earthquake-resistant design have published recommendations based on combination of theory, ex-periments and practical observations, and form the foundation of the official codes (Chopra, 2005).
Two methods for resistant design of buildings and structures are used in practice (Chopra, 2005): the lateral force method, and the dynamic model method.The first method is applicable for seismic zones low seismic risk and lowrise buildings.The buildings are modeled as single degree of freedom systems subject to specific excitation.The dynamic method could be applied to any structure but is normally used for high risk zones and structures over five stories.The structure is modeled as a many degree of freedom discrete system with concentrated masses at each level of the building.In earthquake resistant design of buildings, the maximal response include, displacements, accelerations, shear forces, overturning moments, and torques.Because of the complexity of the dynamics based earthquake resistant design of buildings Chopra and Goel (2002) developed a modal pushover analysis procedure for linearly elastic buildings.The procedure is extended to inelastic buildings and the errors in the procedure relative to a nonlinear response history analysis are documented.
In recent years innovative means of enhancing structure functionality and safety against natural hazards have been in various stages of research and development (Nishitani and Published by Copernicus Publications.Inoue, 2001).They can be grouped into three broad areas: base isolation; passive energy dissipation; applications as active control.Base isolation can now be considered a more mature technology with wider compared with the other two.
Passive energy dissipation systems encompass a range of materials and devices for enhancing damping, stiffness and strength.In comparison with passive control, active control of structural response could be characterized by the following two features: a certain amount of external power or energy added as a result of a decision making process based on real time measured data (Soong and Dargush, 1997).
The inventive methods for earthquake structure response analysis and design are based on precise dynamic modeling and simulation.The up-to-date methods for multibody system dynamics simulation taking into account nonlinear behavior of the large flexible structures could be classify in three basic groups: floating frame of references (Shabana, 2005), finite elements in relative coordinates (Zahariev, 2006); and absolute nodal coordinate formulation (Shabana, 1996;Shabana and Mikkola, 2003).The large structures and mechanical devices like skyscrapers, antenna, wind power generators (WPG), etc., are typical multibody systems compiled of rigid and flexible parts with nonlinear behavior.Shaking of their foundations because of ground wave propagation, as well as, the nonlinear interaction ground -structure basement could be effectively analyzed using the recent achievements in the multibody system dynamics simulation methodology.The forward dynamic analysis provides the designers with preliminary information for structure behavior and the effect of the design parameters.
In the paper finite elements in relative coordinates and generalized Newton-Euler dynamic equations (Zahariev, 2006) are applied for dynamic simulation of large flexible structures.Simulation of the ground oscillation is implemented using statistical data for earthquake ground acceleration during the disasters.The structure basement oscillations are regarded as reonomic constraints in the system of Differential Algebraic Equations (DAE).The structure response is analyzed integrating the dynamic equations so derived.Several virtual examples of motion simulation and flexible deflections of structures and WPG imposed on earthquakes are presented in the paper.

Ground-structure basement interaction
The earthquake wave propagation is well studied and sufficient data collection is at the disposal of the scientists and engineers.But, according to the author opinion, information about the structure basement response because of ground motion and wave propagation is still not well studied and additional investigations are needed.Laboratory experiments together with the simulation results could help for revealing the whole picture of the ground-basement interaction during an earthquake.In some cases the ground, imposed on wave propagation, could behave as a viscous material and the interaction with the basement are to be analyzed by means of hydrodynamics methodology.Obviously, further investigations in this direction are to be conducted.
In the paper, for simulation purposes only, the ground is considered as environment with specific elastic and damping properties, which in certain conditions could experience plastic deflections.In Fig. 2 the two dimensional (2-D) kinematic model of ground-basement interaction of a structure imposed of earthquake is shown.The ground is presented by finite element (FE) discretization (Fig. 2a).The structure is compiled of rigid body basement and a flexible tower of FEs.The precise dynamic model should consider the wave propagation of the ground and the interaction with the basement.Using the methodology of the hydrodynamics the reduced forces (F 1 and F 2 , Fig. 2b) loading the basement could be calculated and included in the dynamic equations of the structure.Implementation of this method requires enormous statistical data for the ground properties and local earthquake wave propagation.The simulation of the basement motion is much simpler when it is built on a rigid ground, for example rocks.Then the basement motion could be considered coincident with the ground motion read from the statistical data.
A simplified model for ground-basement interaction is presented in Fig. 2c.Virtual springs and dampers represent the elasticity and the damping properties of the ground, respectively.The spring properties are presented as a resistant force that depends on the deflections.Since the ground could exhibit residual deflections (virtual plastic deflections) the structure could change its location.The linear ground motion is 3-D and additional 2-D rotation (variation of the basement orientation) is possible.( ) Inertia forces in a rigid body and nodes of finite elements discretization.

Generalized quasi-static dynamics of rigid and flexible systems
The quasi-static dynamics assumes the inertia forces as external loading.The inertia forces in the dynamic equations of finite element discretization are not precisely defined in the finite element theory that is discussed by Zahariev (2006).In Fig. 3 a free rigid body, as well as, a node i of a free flexible element and its neighbor nodes i − 1 and i + 1 are presented.
The linear and angular velocities of the rigid body mass center (v C , ω C ) and of the nodes (v i , ω i for a node i) are also shown.The body inertia forces Φ C and torques T C depend on its own velocities and accelerations, while the inertia forces Φ i and torques T i loading the node i depend not only on its own velocities and accelerations but on these of the neighbor nodes i − 1 and i + 1 too.A node of a flexible element is a coordinate system which motion does not depend on the coordinates of the other nodes.The only thing that keeps the flexible element as a common set of nodes it is the stiffness of the material and the stiffness forces between the neighbor nodes.
In the literature of the finite element theory the inertia forces loading the nodes of a flexible element are defined by the expression (Doyle, 2001) is the matrix vector of the inertia forces and torques in all nodes, M is the global mass matrix of the element, . The inertia forces so computed do not include the velocity dependent terms as it is for the rigid body.The generalized Newton-Euler dynamic equations derived by Zahariev (2006) define the inertia forces in a flexible element and the velocity dependent terms are taken into account, i.e.: where larly to definition of a skew-symmetric matrix of a vector angular velocity ω × i , is a generalized skew-symmetric matrix of the linear and angular velocities of node i; As it could be seen from Eq. ( 1) the inertia forces are expressed with respect to the quasi-velocities and accelerations and do not depend on the kind of coordinates.That is significant advantage for application of the numerical algorithm for computer code generation of the dynamic equations and expansion of the numerical method of the Newton-Euler equations in flexible systems.Absolute velocity values are invariant relative to the kind of coordinate this approach, actually, could be applied to every kind of the coordinates (ω i = ω i (q) ; v i = v i (q), Where q is n × 1 matrix-vector of the coordinates.Using the principle of the virtual work all external forces, including the inertia, stiffness and damping forces are reduced to the system coordinates to derive the dynamic equation.Applying the general partitioning methods or the augmented Lagrangian approach the system of differential equations subject to algebraic constraints could be reduced to a system of n × n ODE with respect to the generalized coordinates q, i.e.: M is the mass matrix, q are the generalized coordinates, S are the generalized forces, B (q, q) is the velocity depend term.

Dynamics of multibody systems subject to reonomic constraints
Reonomic are the kinematic constraints that depend on time.For a multibody system with degree of freedom n which motion is described by ODE, Eq. ( 2), subject to m reonomic constraints qi = qi (t) the dynamic equations are presented as follows: In Eq. ( 4) the coordinates qi = qi (t) , i = k+1, k+2, .., k+m are known, while the generalized forces S i = S i (t) , i = k + 1, k + 2..., k + m are unknown.In other words, the solution of equations system, Eq. ( 4), results in solution of mixed direct and inverse dynamic problem to find the values of the coordinates qi = qi (t) , i = 1, ..., k, k +m+1, ..., n and the generalized forces S i = S i (t) , i = k+1, k+2..., k+m.The dynamic equations, subject to reonomic constraints, Eq. ( 4), are transformed with respect to the unknown parameters as follows: where the matrices S, M, M and are as follows: .

Earthquake dynamic response of a WPG
The WPG is compiled of large flexible bodies (the blades) connected to a shaft (rigid or flexible) that transmits the torque to a gear with elastic clutch and electric generator machine.This structure together with the large flexible tower over which it is placed implements complex motion in space.
It is imposed of wind load with transient behavior and is very sensitive of the changing working conditions of the environment.
Mainly two wind turbine concepts with their control strategies have been applied in practice, namely: active stall control wind turbine with induction generator; variable speed, variable pitch wind turbine with doubly-fed induction generator.The first group operates with almost constant speed.In this case, the generator directly couples the grid to drive train.The second one operates with variable speed.In this case, the generator does not directly couple the grid to drive train.Thereby, the rotor is permitted to rotate at any speed by introducing power electronic converters between the generator and the grid.The constant speed configuration is characterized by stiff power train dynamics due to the fact that electrical generator is locked to the grid; as a result, just a small variation of the rotor shaft speed is allowed.The construction and performance of this system are very much dependent on the mechanical characteristic of the mechanical subsystems.In addition, the turbulence and tower shadow induces rapidly fluctuation loads that appear as variations in the power.These variations are undesired for grid-connected wind turbine, since they result in mechanical stresses that decrease the lifetime of wind turbine and decrease the power quality.
Alternatively, variable speed configurations provide the ability to control the rotor speed.This allows the wind turbine system to operate constantly near to its optimum tip-speed ratio.Although the main disadvantage of the variable-speed configuration are the additional cost and the complexity of power converters required to interface the generator and the grid, its use has been increased due the above mentioned advantages.
There are not enough studies on the influence of the seismic loading to the WPG and in many cases the influence of such loading is not taken into account in the design process.The main purpose of the structural model of a wind turbine is to be able to determine temporal variation of the loads in various components in order to estimate fatigue damage.To calculate the deflections and velocities of various components in the wind turbine in the time domain, a precise dynamic model including the inertia terms is needed.The wind turbines are large flexible structures (more than 150 m in diameter) which blades implement complex motion in space with high velocities and rapidly changing accelerations of the blade deflections.The vibrations cause unstable working conditions, noise and random loads of the units.Because of the changing external loading over the blades the turbine shaft is imposed on bending and torsion including also impacts in the bearings.
Structures of WPGs imposed on earthquakes are investigated here in order to prove the reliability of the numerical algorithms for dynamic analysis subject to reonomic constraints.An example of the dynamic model of a WPG and its response because of prescribed ground motion is presented.In Fig. 4 a simplified presentation of a WPG design scheme, as well as, its finite element discretization are shown.The gear increases the blade angular velocity ω w to an optimal value ω r for the electric power generator, ω e .The elasticity of the gear shaft and the clutch causes different values of ω r and ω e .The reducer is assumed a rigid block with stiffness k r and damping c r , which values are provided by the manufacturer.The power unit is adjusted to an elastic tower and rigid body foundation.2-D ground motion, respectively, two coordinates q 1 and q 2 are regarded.The coordinate q 3 is the shaft rotation of the blades.Friction forces F 1 , F 2 appear between the ground and the foundation.The elastic and damping coefficients of the ground are k g c g .The flexible nodes are pointed out by arrows.The dense circles present the rigid bodies.Each flexible node adds six degree of freedom (dof).Some of them could be neglected (for example the longitudinal deflections).The translational joints numbered with 1 and 2 present the displacements of the ground.The mass of the foundation is m f .The flexible blades are connected to a common hub with mass m h and inertia J h .Joints 3 represents the flexible shaft gear and the hub.The mass properties of the gear are m r , J r .Joint 4 represents the flexibility of the gear.Jount 5 represnts the flexibility of the clutch between the electric generator and the gear (stiffness k e and damping c e ).The inertia proprties of the generator are m e , J e .Details for discretization of a rigid and flexible multibody system, as well as, for deriving the dynamic equations could be found in the paper of Zahariev (2006).The dynamic equations of the WPG (Fig. 4) with n dof, for which the kinematic constraints are transformed or the generalized partitioning approach is applied, are presented by second order ODE, Eq. ( 2).The reonomic constraints presenting 2-D ground motion are as follows: q1 = q1 (t) ; q2 = q2 (t) . (6) Since the values of q 1 and q 2 are defined at every instant the Eqs.( 2) and ( 6) are transformed to linear equation system with respect to qi , i = 3, 4, ..., n and the generalized forces in gondola foundation hub 0 joints 1 and 2 (S 1 , S 2 ), i.e.: E and 0 are 2 × 2 and n − 2 × n − 2 unity and zero matrices, respectively, M is the matrix compiled of the first two colunms of M, S = 0 0 S 3 . . .S n \ .Two types of WPG are used in practice, i.e.: constant speed and variable speed (constant torque) WPG.For the first case there is an additional reonomic constraint of the type ω w = q3 = const.For the second type WPG the torque of the blade shaft is approximately constant that is force constraint imposed on the dynamic equations.Integrating numerically Eq. ( 7) one obtains the values of the generalized coordinates, the forces between the ground and the foundation and, as a result, the global motion of the structure with respect to time.In Fig. 5 and in Table 1 the discretization and mass distribution scheme of a virtual WPG are presented, as well as, its design and mass parameters.The structure is compiled of three rigid bodies, foundation, hub and gondola.The tower and the three blades are assumed flexible.Each of them is divided into three substructures with decreasing size from the basement and the hub to the tips.It is assumed two reonomic constraints.The first one q 1 = q 1 (t) is a ground motion perpendicular to the plain of the blades and the second one q 2 = q 2 (t) is prescribed constant rotation of the blade shaft.It is shown the orientation of the coordinate systems of the corresponding substructures and bodies.
Two numerical experiments are conducted.The first one (Fig. 6) is simulation of the WPG earthquake response and the resulting tower and blade tips displacements in case of no rotation of the blades (still position of the turbine).The second one (Fig. 7) is in case of turbine rotation with the prescribed angular velocity ω w = 3 rad s −1 .As it could be seen from the examples the deflections of the tower and the blades are several times higher when the turbine rotates.This is as a result of the inertia forces caused by the velocity dependent  terms and the mutual influence of the ground motion and the rotation of the turbine.The examples demonstrate the effectiveness of the numerical algorithm, as well as, the necessity of the preliminary analysis of large flexible structures in various working conditions and scenarios.

Conclusions
Numerical algorithm for dynamic simulation of large flexible structures imposed on earthquakes is developed.Precise dynamic model of the flexible system is derived using generalized Newton-Euler dynamic equations for rigid and flexible bodies.A procedure for transformation of DAE subject to reonomic constraints to a system of ordinary ODE is pro- posed.Examples of motion simulation of WPG imposed on earthquake are depicted.

Figure 2 .
Figure 2. Modeling the ground-basement motion and force interaction.

Figure 4 .
Figure 4.A design scheme and the discretized structure of a wind power generator.

Figure 5 .
Figure 5. Example of motion simulation of constant speed WPG subject to ground shaking.

Figure 6 .
Figure 6.Displacements of the WPG imposed on earthquake without rotation of the turbine.

Figure 7 .
Figure 7. Displacements of the rotating WPG imposed on earthquake.

Table 1 .
Design and mass parameters of the WPG virtual structure (all measures are in SI).