A co-rotational formulation for quasi-steady aerodynamic nonlinear analysis of frame structures

The design of structures submitted to aerodynamic loads usually requires the development of specific computational models considering fluid-structure interactions. Models using structural frame elements are developed in several relevant applications such as, the design of advanced aircraft wings, wind turbine blades or power transmission lines. In the case of flexible frame structures submitted to fluid flows, the computation of inertial and aerodynamic forces for large displacements and rotations is a challenging task. In this article, we present a novel formulation for the efficient computation of aerodynamic forces in frame structures, coupling the co-rotational framework with the quasi-steady theory. A numerical procedure is provided considering a tangent matrix for the aerodynamic forces. This formulation is implemented in the open-source library ONSAS, allowing users to reproduce the results or solve other frame nonlinear dynamic problems. The proposed formulation and its implementation are validated through the resolution of four numerical examples. The formulation and the numerical procedure proposed efficiently provide accurate solutions for these challenging problems with large displacements and rotations.


Introduction
Nonlinear structural dynamic problems are formulated in a vast and diverse set of applications such as: developing new wind turbines systems [2,39], designing suspended bridges or aircraft wings [47,6], predicting failures in power transmission lines [41], reducing fruit production losses [8] or even studying the movement of aquatic plants [20].In all of these applications, structures can be modeled using frame elements, and are also submitted to loads caused by the interaction with fluid flows.The development of an efficient and accurate numerical method for the resolution of this type of problems is the main motivation of this article.
Structural design standards have a limited range of application, and are not applicable to most of problems mentioned above [14,38].Given this limitation, alternative approaches are mainly based on experimental tests [8] or numerical simulations [42].Experimental tests might be expensive and/or challenging to design, therefore, new numerical methods for accurate structural dynamics simulations are actively developed [17].
The Finite Element Method (FEM) [48] has become the gold-standard for computational modeling in structural analysis in numerous disciplines.For frame structures, the co-rotational approach has shown several advantages, including a more versatile and less intricate mathematical formulation [10].This approach is based on splitting the element deformation in: one rigid movement and one local deformation [5].Different co-rotational formulations were developed for solving static [31], stability [4] or dynamic [25] structural analysis problems.In [30] a co-rotational formulation using a projector matrix is presented and applied to nonlinear solid analysis.A consistent formulation for three-dimensional nonlinear dynamic analysis of frame structures was presented [26], allowing to accurately simulate deformations with large displacements and rotations using a reduced number of elements.In [44] it is shown that, for structures submitted to large rotations, the consistent formulation is considerably more accurate and efficient than the lumped mass approach.
In the last decades different frame analysis formulations were used for the mentioned applications of interest.In [13], a three-dimensional nonlinear three-node isoparametric element is used for modeling the movement of overhead transmission lines, considering a consistent mass matrix for linear inertial terms.With the same purpose, in [40] a three-dimensional linear frame element was used to simulate cable elements.In [28], a formulation considering nonlinear internal forces with a lumped mass matrix for linear inertial terms was used for modeling wind turbine blades.In [39] a vorticity wind turbine was modeled using the discrete element method concluding that, large displacements must be considered to emulate states of maximum output power.Regarding the nonlinear geometric analysis of wind turbine blades, the linearized equations of motion were solved in [16], obtaining a good level of agreement between the exact beam theory and formulations using shell elements.In [11,29], static co-rotational formulations were used to simulate morphing or highly flexible wings, highlighting the computational efficiency to validate experimental results.However, the performance of a formulation considering consistent inertial terms and aerodynamic forces using the co-rotational framework, has not been reported.
Regarding the availability of software for the numerical resolution of these problems, three specific tools can be mentioned: RIFLEX, FAST and HAWC2.RIFLEX is a proprietary software developed for fluid-structure interaction problems.It uses a co-rotational approach for modeling frame elements, a linear consistent mass matrix, a Rayleigh damping matrix and allows to compute mass, shear and elastic centers [9,12].FAST is a modular open-source framework for fluid structure numerical simulations.This software uses beam elements based on the exact beam theory and a Timoshenko mass matrix [45,32].HAWC2 is a proprietary software for aeroelastic simulation of wind turbines, developed using linear anisotropic Timoshenko beam elements [24].
In [19] and [21], Euler-Bernoulli and Kirchhoff nonlinear beam formulations were used to simulate the deformation of flexible structures submitted to drag and lift forces.The deformation of thin elastic blades was studied in [27], concluding that inertial effects have a significant impact on the numerical results of the model.Furthermore, in [33] a nonlinear aeroelastic study of wings was conducted concluding that deformation is the dominant nonlinearity for long slender wings.
In [35] the quasi-steady theory and an equivalent beam model with a lumped mass matrix were employed to analyze galloping effects on buildings.In [18] the co-rotational framework is introduced for the computation of aerodynamic forces, neglecting the pitch torsional moment, and a linear lumped mass matrix approach is used for the inertial effects.The tangent matrices of the aerodynamic forces vector were also neglected in the iterative numerical scheme applied.Therefore, there is still a research gap in quantifying the benefits of using a consistent inertial formulation and the pitch moment.
In this work, we present a novel unified formulation for consistent co-rotational analysis of frame structures submitted to nonlinear aerodynamic forces.For the first time, the fluid interaction effect is included by considering the quasi-steady theory, a consistent corotational formulation and the Principle of Virtual Work.In particular, in contrast with [18], in the proposed formulation the aerodynamic pitch moment is not neglected.In addition, our results showed that considering the aerodynamic stiffness tangent matrix, which is not usually computed in the literature, improves the accuracy of the solution in static analysis problems.
The proposed formulation is implemented in the open-source structural analysis solver ONSAS [34], providing a new tool for the scientific community.We perform numerical analyses for different flow conditions, cross-sections and magnitudes of displacements and rotations, studying changes in mesh sizes and number of Gauss numerical integration points.The formulation and its implementation are tested through the resolution of four numerical examples.In the first example the implementation is validated using the open-source tool from [19].All the scripts used in the numerical examples, are publicly available allowing any user to automatically reproduce the results presented.
This article is organized as follows.In Section 2 the basic concepts of the co-rotational framework are described.In Section 3, the proposed formulation is presented, with a corresponding numerical procedure for the resolution of the balance equations.In Section 4 the numerical results obtained are presented, and in Section 5, the conclusions obtained are described.

Preliminaries
In this section, the fundamental concepts of the co-rotational frame analysis approach are described.The main kinematic identities and the internal and inertial forces are briefly presented considering [4,26].

Co-rotational kinematics
The main concepts behind the co-rotational approach are the use of different systems of coordinates and the application of the Principle of Virtual Work.Given a two-node frame element and two systems of coordinates (global and local), a vector of generalized nodal displacements  can be represented in the global system of coordinates as   , and in the local system as   .The virtual work of a set of nodal forces  is the same in the local or the global system of coordinates, as it is presented in Equation ( 1): for any vector of virtual displacements .In the co-rotational approach three configurations are defined as shown in Fig. 1: a reference configuration, a rigid-rotation configuration and the total-deformation configuration.As it is shown, four systems of coordinates are defined: , corresponding to the canonical, reference, rigid-rotation and total-deformed configurations, respectively.Orthogonal matrices  0 ,   ,   and  can also be defined as shown in Fig. 1, to rotate the base vectors of these systems of coordinates.
The column vector of nodal displacements written in the canonical system {   } is denoted as   = [ ( 1 )  , ( 1 )  , ( 2 )  , ( 2 )  ]  , where   and   are the column vectors of linear displacements and rotations, respectively, of node .In the co-rotational approach, the displacements of the element are also written considering the system of coordinates ]  .The extension is given by  =   −  0 where   and  0 are the deformed and reference lengths of the element, and the local rotations are given by the vectors   as it is shown in Fig. 2. In order to apply the Principle of Virtual Work, the vectors of variations of the generalized displacements  need to be written in the same system of coordinates.The variation of the local extension verifies Equation (2): and for the vectors of rotations Equation (3): where  × represents a matrix of zeros with  rows and  columns (the sub-indexes are omitted for the 3 × 3 case),  is a matrix given by Equation ( 4): and  is a matrix given by Equation ( 5): with   being the -th entry of the vector   is defined in Equation ( 6): and   being the -th entry of the vector  defined by  = 1 2 ( 1 +  2 ).For a cross-section located at the position , as shown in Fig. 2, with centroid  and deformed base , the variations of the displacements and rotations can also be written in local and global systems, using Equations ( 7) a) and b): and Equations ( 8) a) and b): where  1 and  2 are the linear interpolation functions (for axial displacement) and  3 ,  4 ,  5 and  6 are Hermite interpolation functions (for bending).
The position of  in canonical coordinates can be written as: where   are the local transverse displacements.Considering Equation ( 9), the variations of the displacement and rotation of the point  can be written as: respectively, where  2 =  2  +   and  1 =  +  1  − ũ   , with ũ being the skew operator associated with the vector   .Finally, velocities and accelerations can be obtained using Equation (11): where

Internal and inertial forces
The expressions of the elemental internal and inertial forces in global coordinates can be obtained using the Principle of Virtual Work.Considering Equation (1) for the internal forces, and substituting the relations presented in Equations ( 2) and (3), we obtain Equation ( 12): This identity is valid for any virtual displacement   , thus we obtain the Equation ( 13): where    is the known vector of internal forces , with normal force and bending moments, given by a linear constitutive behavior.
For the inertial term, the kinematic energy  of the element is written as: where  is the area of the cross-section,  is the density of the material and  is the geometric inertia tensor.Considering the variation in both members of Equation ( 14) it can be obtained Equation ( 15): The inertial force vector of the element in global coordinates    is then defined consistently by Equation ( 16):

Methodology
In this section we present the proposed formulation for the computation of the aerodynamic forces, and describe a numerical procedure for the resolution of the governing equations.

Co-rotational quasi-steady aerodynamic forces
Let us consider a frame element, with uniform cross-section, submitted to a fluid flow as shown in Fig. 3.For a section located at  0 with centroid , the deformed position at time  is given by  =   ( 0 ).The element is submitted to forces induced by a fluid with absolute velocities given by the field   (, ) ∶ ℝ 3 × ℝ → ℝ 3 .The velocity of the centroid is ̇( 0 , ) and the relative velocity in the deformed position is defined by: In this definition a fundamental assumption was considered: the movement of the structure does not affect the absolute velocities of the fluid [7].
The interaction between the fluid flow and the frame element produces normal and shear stresses, that are represented by moments and forces (generalized forces) applied at the deformed position of the centroid.These forces are assumed to be uniquely defined in terms of the instantaneous position and velocity of the deformed section [18].In particular, in this formulation, the forces are assumed to depend only on   (the projection of the relative velocity onto the plane Π 23 defined by  2 and  3 ) as it is shown in Fig. 3.It is assumed that the density of the fluid is considerably lower than the density of the structure, therefore the added-mass effect is neglected.The quasi-steady aerodynamic distributed forces for drag, lift and torsional moment are given by the expressions: respectively, where   is the density of the fluid,   is the given characteristic dimension of the cross-section and   ,   and   are the drag, lift and moment coefficients, determined by wind tunnel tests for different Reynolds numbers  and angles of incidence .The angle  is defined by   and the unitary chord vector of the section   , as shown in Fig. 4. It is remarked that drag and lift force vectors are included in the plane Π 23 .
The vector   written in the total-deformed system of coordinates is denoted as (  )  .In the same manner, the notation (•)  is used for any vector in this system and this sub-index is omitted for vectors in the canonical system.The expression of (  )  is given by Equation ( 19): where ( 1 )  = [1, 0, 0]  .Using the rotation matrices of the co-rotational framework as change of basis operators:   =  ()  ,  =  ()  , we can write: Substituting Equation ( 20) in (19) and defining a projection operator  2 we obtain Equation ( 21): Using this we can define the unitary vectors (  )  , (  )  and (  )  in Equations ( 22), ( 23) and ( 24): with  3 = exp([∕2, 0, 0]  ).The angle of incidence  verifies Equation ( 25): and considering that   and   are unitary we obtain the expression: where a convention was considered as shown in Fig. 4. Using Equation (26) the angle of incidence  can be determined for any deformed configuration.Substituting the identities obtained above in Equation ( 18) we obtain: M.C.Vanzulli and J.M. Pérez Zerpa The virtual work corresponding to the aerodynamic forces of the element is given by Equation ( 28): where  + is the sum of   and   .Considering the vector of nodal aerodynamic generalized forces in global coordinates    , the virtual work can also be written as: and substituting Equations (10) a) and b) in ( 28) and using Equation ( 29) we obtain Equation ( 30): Operating we obtain Equation (31): and substituting Equation ( 27) the complete expression of the aerodynamic forces vector is obtained in Equation ( 32): It is important to highlight that, in contrast with [18], in the proposed formulation the pitch moment is not neglected.

Balance equations and numerical resolution procedure
The governing equations are obtained by considering the virtual work for all the elements of the structure for the forces in Equations ( 13), ( 16) and (32).Additionally a vector with external forces not induced by the fluid interaction    can be added.The nonlinear system of governing equations is written in Equation (33): where the arguments of the residual forces   were omitted.The numerical resolution procedure proposed consists in solving the system of nonlinear governing equations using iterative methods.For the static analysis cases the Newton-Raphson method is used, while for dynamic analysis cases the Newmark method with   = 1∕4 and   = 1∕2 [3], and the -HHT method with   = −0.05[26] are used.The computation of the aerodynamic and inertial force vectors is done using numerical integration.Three different formulations are considered and used to generate the numerical results: F1: In this formulation a lumped mass approach is considered for the inertial terms and the aerodynamic forces are computed using the equations presented in Section 3.1, with   = 0 (neglecting torsional moment).Also, the tangent matrices of the aerodynamic force vector are neglected.This formulation is considered to be equivalent to the one presented in [18].F2: In this formulation the consistent approach [26] is considered for the inertial terms, and the aerodynamic force vector is computed using the equations presented in Section 3.1 (without neglecting the torsional moment).The tangent matrices of the aerodynamic force vector are neglected.F3: In this formulation the forces are computed as in F2, however for the tangent matrix the contribution of the aerodynamic forces is added.The aerodynamic stiffness matrix is computed for each frame element using a finite difference approach given by Equation ( 34): where    , is the -th column of the aerodynamic stiffness matrix, ℎ = 1 × 10 −10 m and   is a canonical vector.

Numerical results
In this section the numerical results obtained for four problems are presented.Unless it is specified the fluid considered is air with density   = 1.225 kg/m 3 , kinematic viscosity   = 1.5 × 10 −5 m 2 /s, at 20 • C and atmospheric pressure.Regarding the elastic properties, Poisson's ratio  = 0.3 is considered for the first four examples.For all the problems, homogeneous initial conditions are considered.
All the numerical results presented can be reproduced by running scripts publicly available. 1he results shown were produced using a computer with a Linux OS, a 64-bit architecture, an Intel i7-11370H CPU and 16 Gb of RAM, running the implementation of the formulations in ONSAS on GNU-Octave [15].The visualization is done using Paraview [1] and GNU-Octave.
The stopping criteria considered in all the examples are given by Equation ( 35): where  is the number of iteration and   ,   are scalars to be defined.

Example 1: drag reconfiguration of a cylindrical cantilever beam
In this example a cantilever beam submitted to a flow producing drag forces is considered.The main goal is to obtain results about the performance of the three formulations described in Section 3.2, for a problem with large displacements.The example is based on one of the problems considered in [19], where a reference solution is presented and validated with experimental data.Finally, a brief numerical study on the variation of the results for different numbers of Gauss integration points is presented.

Problem definition
The problem consists in a cantilever beam submitted to a fluid flow with uniform velocity   (, ) =   () 2 , as shown in Fig. 5.The beam is clamped on the boundary at  = 0 m, and the span length is  = 1 m.The cross-section of the beam is circular with diameter  = 1 cm, and the chord length used to compute the aerodynamic forces is   = .For the material of the beam a linear elastic isotropic model is considered, with Young modulus  = 30 MPa and density  = 7000 kg/m 3 .
The fluid considered is water with density   = 1020 kg/m 3 and kinematic viscosity   = 10 −6 m 2 /s.The values of the aerodynamic coefficients are considered from [37] as:   = 1.2 and   =   = 0.

Numerical results: static case
The goal of this analysis case is to study the results obtained with formulations F1, F2 and F3 for fluid velocities   ∈ {0.005, 0.012, 0.029, 0.072, 0.176, 0.432, 1.057, 2.588, 6.336, 15.512} all in m/s.This is a static analysis case, thus the velocity of any point of the beam is u = 0. Substituting this in Equation ( 17) we obtain   (, ) =   (, ).In Fig. 6 the absolute velocity of the fluid and its relative projected transversal component are shown, and for this analysis case, we obtain the identity:  Equation (36) indicates that once the beam is deformed, the norm of the projected velocity decreases, and therefore the drag force decreases.This geometric nonlinearity effect is called reconfiguration.Due to this reconfiguration mechanism the drag load does not increase with the square of the fluid velocity [19].The problem can be studied through the dimensionless Equation (37): where  is the global drag force towards  2 ,   is the Cauchy number that describes the ratio between the stiffness of the beam and the flow load and the reconfiguration number  reflects the geometric nonlinear effect by dividing the drag of the flexible beam to that of a rigid one of the same geometry.
For each velocity   , the numerical formulations are used to obtain the solution and the drag force  is computed.The N-R method is used with only force stopping criteria considering   = 10 −8 .The results obtained for  (using 2 and 20 elements) are shown in Fig. 7a, and a reference solution is also obtained using the source code 2 developed for [19].
It is observed that the three formulations match the reference solution for   ⩽ 0.072 m/s.Moreover, when the aerodynamic load exhibits high geometrical nonlinearities (  > 0.072 m/s) F3 converges while F2 and F1 do not.This demonstrates that in this case, the stiffness aerodynamic matrix is necessary to accurately reproduce the reference results.Furthermore, a larger number of elements is required to verify the reference solution accurately in cases with   ⩾ 1.057 m/s.
It is reported that, for   = 0.072 m/s and using 20 elements, the formulation F3 requires 6 times the execution time required by F1 or F2 formulations.

Numerical results: dynamic case
In this case a nonlinear dynamic analysis is performed and the flow velocity   is given by Equation (38): 2 https://github .com/lm2 -poly /Reconfiguration -Beam.where   = 7 s and two values are considered for   : 6.3355 m/s (Case 1), and 1.0568 m/s (Case 2).The trapezoidal Newmark numerical method is used, with time step Δ = 0.05 s and the stopping criteria are   = 1 × 10 −10 and   = 1 × 10 −4 .The solutions obtained for the displacement   of point A are shown in Fig. 8.The solutions provided by the three formulations F1, F2, and F3, using 20 elements for both cases, converge to the reference solution.Moreover the formulation F1, using 4 elements in Case 1, does not converge after  ≈ 1.5 s.
In order to analyze the numerical behavior of the formulations described, a study considering different number of integration Gauss points for the formulation F2 is presented.The displacement functions obtained at time at  = 7 s for different number of integration points are compared.The difference between the functions is computed considering Equation (39): where the reference solution is obtained using 20 elements and 10 Gauss integration points.The results obtained are shown in Fig. 9.
It is noted that at  = 7 s a considerable curvature is present in the displacement field (  of node A is ≈ 1 m), thus, a high nonlinearity is also present in the aerodynamic forces term.Given the results shown in Fig. 9, 4 Gauss points are considered for integrating the aerodynamic forces.
The results obtained let us conclude that the formulations considered can be used to solve the problem.For static cases, including the aerodynamic stiffness matrix is necessary to accurately solve problems with large bending deformations.For dynamic cases, the three formulations considered are able to provide accurate results when an appropriate discretization is used.Formulation F3 requires more computational time than F1 and F2, therefore it is recommended only for static analyses of structures, submitted to large deformations.

Example 2: simple propeller model
In this example a simple propeller submitted to lift forces is considered.This problem is used to validate the results provided by the proposed formulation in a dynamic case with large displacements and rotations.

Problem definition
The problem consists in a three-blade propeller submitted to a flow with uniform velocity.Each blade has a length  = 3 m and a circular cross-section with diameter  = 0.1 m, as shown in Fig. 10a.Regarding the stiffness of the blades, two cases are considered: a rigid case (with analytic solution) and a flexible case, providing considerably large rotations and bending of the blades.The density  = 6000 kg/m 3 is considered.Regarding boundary conditions, the node O has five degrees of freedom fixed: the three displacements and rotations   and   .The rotation  ,O is free.
A uniform flow   = 1  1 m/s is applied with synthetic aerodynamic coefficients   = 0.2 and   =   = 0. Given this, a uniform lift distributed force   contained in the plane  2 - 3 is induced, as shown in Fig. 10b.These specific settings allow us to obtain an analytic solution for the rigid case.
The -HHT numerical integration method is used, with a time step increment set to Δ = 1 s.

Numerical results: rigid case
For this case, and using the value selected for   , it can be assumed that ̇ ≪   , thus   ≈   .For the considered properties and boundary conditions, and for a Young modulus  = 210 GPa, the bending deformation of the blades can be neglected, allowing to obtain an analytic solution.
Considering  ,O ≈  ,A =   , the angular momentum balance equation is written in Equation ( 40): and using the homogeneous initial conditions, Equation ( 41) is obtained: For the numerical resolution, the tolerances   = 10 −6 and   = 10 −12 are set, and the final simulation time is   = 450 s.
The results obtained for  , are shown in Fig. 11a, where it can be observed that the analytic solution is verified with the results provided by the proposed F2 formulation, even using one element per blade.In contrast, as expected, the formulation F1 requires the use of five elements per blade to match the analytic solution.In Fig. 11b the deformed configurations are shown at  = 100 s.

Numerical results: flexible case
The goal of this case is to test the proposed formulation for a highly-flexible propeller.To do so, a Young modulus  = 2.1 kPa is considered.This local flexible behavior plays a key role in applications such as the design of morphing wings [46,43,11], and represents a challenge for numerical methods.In this case, due to the bending deformation, the rotations of points O and A, shown in Fig. 10a, are considerably different.The numerical results obtained for the point O and point A using formulation F2 are presented in Figs.12a and 12b, respectively.Formulation F1, using 5 and 20 elements, was not able to provide a numerical solution, hence the solutions were discarded.The deformed configurations obtained using F1 and F2 formulations at  = 285 s are shown in Fig. 13.
The results let us conclude that, in this problem, for flexible elements with rotations larger than 2 and 20 elements per blade, formulation F1 is not able to provide a solution and formulation F2 converges for all time steps.
The results obtained let us conclude that the F2 formulation provides accurate results for large rotations and considerable bending deformations.We consider that the differences in performance between the formulations are due to the approach used for computing the inertial terms.In the following example, a non-zero aerodynamic pitch moment is considered, providing a more complete comparison.

Example 3: simplified cantilever blade
In this example a cantilever beam with an airfoil cross-section is considered.Realistic drag, lift and moment aerodynamic coefficients are used.The goal of this example is to compare the solutions obtained with F1 and F2 formulations, considering realistic aerodynamic coefficients.

Problem definition
In this example a cantilever beam submitted to a fluid flow with varying direction, as illustrated in Fig. 14a.The cross-section of the beam is given by a NERL S809 wind turbine airfoil, with  = 10 m and chord length of   = 1 m [23].The aerodynamic coefficient functions, obtained from [36], are shown in Fig. 15.The geometry of the problem is inspired on an specimen presented in [16] and for the material properties, equivalent Young modulus   = 14 GPa and shear modulus   = 5.6 GPa are adopted.The flow velocity is uniform and its expression is given by Equation (42): with   = 30 m/s and  ∈ [0 • , 40 • ].The problem is solved considering a static analysis for each value of  considered.The change in  can be associated with a slow change in the pitch angle during the operation of a wind turbine.

Numerical results
For the numerical resolution 10 co-rotational frame elements with F1 and F2 were used.The N-R method is used, considering   = 5 × 10 −7 and   = 10 −15 .For the computation of the numerical solution, 4 Gauss integration points were considered for computing the aerodynamic forces.
The results obtained for the bending moments   ,   , are presented in Fig. 16a.Additionally, the resultant shear forces   and   at node O, are shown in Fig. 16b.The torsional moments at point 0 (  ()), for both formulations, are illustrated in Fig. 17 It can be observed the formulation F1 largely underestimates the torsional moment at point O.The results obtained let us conclude that F1 should not be used for the resolution of problems in which the torsional moment is relevant.On the other hand, formulation F2 provided appropriate results.

Example 4: simplified wind turbine 4.4.1. Problem definition
In this example a flexible frame structure undergoing significantly large rotations is considered.For this, a simplified wind turbine model is developed, where the fundamental features of the real problem are present.
The problem consists in an idealized wind turbine as shown in Fig. 18.Each blade has a uniform NERL airfoil with the geometry and material properties presented in Section 4.3.The aerodynamic coefficients were extended to obtain values between -30 • and  90 • based on [22].A uniform constant wind velocity   27  1 m/s is considered.The initial conditions are homogeneous both for displacements and velocities of the blades.

Numerical results
The -HHT numerical method is used with a time step Δ  = 0.01 s and a final time   = 30 s.The residual force and displacement tolerances are:   = 10 −5 and   = 10 −10 .The spatial discretization of each blade is done using 30 aerodynamic co-rotational elements.The formulation F1 was not able to provide a numerical solution, while F2 provided the results shown.
The numerical results obtained for the rotation   of point O, are shown in Fig. 19a, and the results of the angular velocity θ  () of point O, are shown in Fig. 19b.The horizontal   and vertical   displacements of node A are shown in Fig. 20.
The results let us conclude that formulation F2 is able to provide a numerical solution of the problem.Moreover it can be observed that, as expected, the angular velocity does not diverge and a quasi-stationary regime is reached.

Conclusions
In this article a new formulation for the numerical analysis of frame structures submitted to aerodynamic forces is presented.The methodology extends the application of the co-rotational approach, for computing the quasi-steady aerodynamic forces in the deformed configuration for large displacements and rotations.The co-rotational approach is used to compute aerodynamic, internal and inertial forces, providing a set of nonlinear governing equations.An aerodynamic stiffness matrix is added to the tangent matrix in the numerical procedure using a finite difference approach.Three formulations were considered: F1, which is considered equivalent to a previous work of the literature, and F2/F3 two variants of the proposed co-rotational formulation.

Fig. 5 .
Fig. 5. Example 1: Diagram of cantilever beam with boundary conditions and fluid flow.

Fig. 9 .
Fig. 9. Example 1: Relative error   for different numbers of Gauss integration points at  = 7 s.
. The torsional moment at point 0 provided by F1 with  = 40 • is 0.03 N m, while F2 provided a moment −1819 N m.