Geometrically Nonlinear Analysis of Composite Beam Based on Space-Time Finite Element Method

1. Introduction

Many engineering components can be idealized as beams that have one of their dimensions larger than the other two, e.g., helicopter rotor blades. In the field of aircraft design, beam structure is widely used as the main load-bearing structure of which dynamic problem has always been one of the important topics.

In recent years, the increasing complexity of rotor structure and analysis conditions has put forward higher requirements for blade dynamics models, among which the fully intrinsic geometrically exact beam theory [1] stands out with its unique form of first-order nonlinear partial differential equations that independent of displacement and rotation variables. Ghorashi [2] applied the general nonlinear intrinsic equations of the beam to solve the elastodynamic response of the hingeless composite blade under the condition of accelerating speed, and successfully analyzed the nonlinear response of the active airfoil blade and its influence on the steady-state load. Amoozgar [3] modeled the rotation dynamics of morphing composite blade by using full intrinsic geometrically exact beam equations, and solved them by finite element method. The convergence results were in good agreement with the experimental results. In [4], modeling methods of fully intrinsic geometrically exact beam under different boundary conditions are explored, and used for static equilibrium, steady-state motion, and linearized dynamic analyses. The results show the full intrinsic method is particularly suitable for blade dynamics analysis.

The derivation of full intrinsic form starts from the study of geometrically exact beams, which first appeared in the work of Kirchhoff and Clebsch [5] in 1944. Subsequently, according to different solving variables, the beam theories can be divided into three types: displacement-based scheme [6] that reduces the complexity of the formula by using high-order truncation, mixed variational scheme [7] that introduces generalized velocity and generalized strain as unknown variables, and fully intrinsic scheme [8] which completely eliminate displacement and rotation variables. The concept of intrinsic was first introduced by Green and Laws [9] in 1966 in the process of studying elastic rods. In 1977, Hegemier and Nair [10] proposed a fully intrinsic elastic rod theory with large deformation and small strain, that governing equations were composed of load and velocity variables, and were used to solve the steady-state motion of helicopter rotor blades. In 1990, Hodges [7] has developed an exact intrinsic formulation, which can be written in a simple matrix form with only second-degree nonlinearities and ideally suited for finite element method. The feasibility of finite element method solve mixed variables equation is verified in [11]. With the canonical form of the equations of motion, in 2003 Hodges [8] derived a set of generalized strain-velocity compatibility relations. Then make up a complete set of fully intrinsic and geometrically accurate partial differential equations of beam dynamics, which are energy conserved.

The fully intrinsic geometrically exact beam equation is a kind of nonlinear time-domain dependent partial differential equation, which is usually solved by a discrete process of first space then time, and the partial differential equations are transformed into time-domain recursive equations. Since the finite difference method is a simple and convenient spatial discreting method, Patil and Hodges [12] applied the fully intrinsic model to the dynamic analysis of highly flexible wing with the central difference scheme. This spatial discreting scheme was later used by Khouli [13] to analyze the hinged swept blade. In 2011, Patil and Althoff [14] used the continuous Galerkin method to discrete the fully intrinsic equation in space, and fitted the variables with Legendre polynomials. The results show that the Galerkin method can obtain more accurate solutions with fewer trial functions. In the same year, Patil and Hodges [15] proposed the concept of variable order finite element by changing the order of Legendre function, which compared the accuracy, computational efficiency and applicability to general problems. The results presented show that one can obtain approximately third-order, fifth-order, seventh-order, and ninth-order convergence for the linear, quadratic, cubic, and quartic finite elements. And one can use higher-order finite elements for a better approximation of mode shapes. In recent years, there have also been many spatial solving formats with better efficiency and accuracy, such as GDQ method [16] and neural network model [17]. However, these numerical algorithms always need to recurse gradually over time, often require a lot of time cost to iterate, have cumulative errors and greatly depend on the initial value.

The space-time finite element method is a high-precision and effective method to solve the problem of the development of time partial differential. Its main feature is to treat the time domain as the space domain for unified discretization, forming an unstructured grid in the sense of space and time, which is consistent for any order of discretization. Higher-order methods can be naturally embedded in these schemes, and have excellent characteristics of high calculation accuracy and moderate calculation amount in the field of nonlinear dynamics problems. Currently, there are relatively few time finite element methods for solving fully intrinsic geometrically exact beams. In [8], Hodges adopted the central difference scheme to discrete the space and time of equations simultaneously. Chen Lidao [18] adopted a DQ-Pade method with high order accuracy to construct the space-time finite element equation of the fully intrinsic model. In 2023, Chen Lidao [19] established a space-time finite element steady-state scheme based on Galerkin weighted residual method with Legendre functions, which has been successfully applied to non-conservative load issues. This discrete scheme is simple in structure and easy to assemble, but only the influence of linear and second-order schemes on the results was compared.

The present work is based on the general steady-state format derived in Ref. [19], The basic process of the Galerkin space-time continuous finite element method of the fully intrinsic geometrically exact beam is: Firstly, the space-time element control equation is obtained by discretization of the global elements along the spanwise and azimuth Angle of the blade. Then, assembling the element equations according to certain laws to get the total algebraic matrix equation of beam. Finally, solve it by Newton method and reverse displacement. Besides, the discontinuity treatment of blade structure, energy conservation of equation, initial value format, conservative load treatment and the influence of variable order are investigated. And a series of static deformation and dynamic response issues are given to verify the effectiveness of proposed algorithm.

2. Fully Intrinsic Beam Equations

Ignoring the warping deformation in the profile, the full intrinsic equations of the geometric exact beam theory can be expressed in a compact form as

$$\begin{array}{l}{F^{\prime }}_B+(\tilde{k}+{\tilde{\kappa }}_B)F_B+{\mathrm{f}}_B={\dot{P}}_B+{\tilde{\mathrm{\Omega }}}_B{P}_B\\ {M^{\prime }}_B+(\tilde{k}+{\tilde{\kappa }}_B)M_B+({\tilde{e}}_1+{\tilde{\gamma }}_B)F_B+{\mathrm{m}}_B={\dot{H}}_B+{\tilde{\mathrm{\Omega }}}_B{H}_B+{\tilde{V}}_B{P}_B\\ {V^{\prime }}_B+(\tilde{k}+{\tilde{\kappa }}_B)V_B+({\tilde{e}}_1+{\tilde{\gamma }}_B){\mathrm{\Omega }}_B={\dot{\gamma }}_B\\ {{\mathrm{\Omega }}^{\prime }}_B+(\tilde{k}+{\tilde{\kappa }}_B){\mathrm{\Omega }}_B={\dot{\kappa }}_B\end{array}$$

(1)

Here, $( {)}^{\prime }$ is the first partial derivative of space, and ${(}^{·})$ is the first partial derivative of time, $k$ is the initial twist/curvature of the beam, $e_1={\left[\begin{array}{ccc}1& 0& 0\end{array}\right]}^T$, $\mathrm{f}$ and $\mathrm{m}$ respectively represent the element distribution force and element distribution moment loaded on the beam that the direction is perpendicular to the beam axis. $F$, $M$, $V$and $\mathrm{\Omega }$ are respectively representing the force, moment, linear velocity and angular velocity of the beam section at the location $x$ and time $t$. $\gamma$, $\kappa$, $P$ and $H$ respectively representing the force strain, moment strain, linear momentum and angular momentum of the corresponding beam section. Note that all the intrinsic variables are expressed in the deformed base B as ( )_B. $\widetilde{( )}$ is the cross product operator, taking $k$ as an example:

$$k=\left[\begin{array}{c}{k}_1\\ k_2\\ k_3\end{array}\right]，\tilde{k}=\left[\begin{array}{ccc}0& -k_3& k_2\\ k_3& 0& -k_1\\ -k_2& k_1& 0\end{array}\right]$$

(2)

The constitutive equations relate the momenta to the velocities and the strain measures to the sectional loads as

$$\left\{\begin{array}{c}\gamma \\ \kappa \end{array}\right\}=\left[\begin{array}{cc}R& Z\\ Z^T& T\end{array}\right]\left\{\begin{array}{c}F\\ M\end{array}\right\}\begin{array}{cccc}& & & \end{array}\left\{\begin{array}{c}P\\ H\end{array}\right\}=\left[\begin{array}{cc}G& K\\ K^T& I\end{array}\right]\left\{\begin{array}{c}V\\ \mathrm{\Omega }\end{array}\right\}$$

(3)

Where, the matrix $R$, $Z$ and $T$ constitute the flexibility characteristics of the beam section, and the matrix $G$, $K$ and $I$ respectively represent the linear density, mass eccentricity and moment of inertia. Equations (1) and (3) constitute a complete set of differential equations for first-order partial derivatives in time and space.

For a hingeless beam with length L, the spatial boundary and periodic continuity conditions of the steady-state problem can be summarized as follows:

$$\begin{array}{l}F(L,t)=F_L,M(L,t)=M_L,V(0,t)=V_0,\mathrm{\Omega }(0,t)={\mathrm{\Omega }}_0\\ \gamma (x,0)=\gamma (x,T),\kappa (x,0)=\kappa (x,T),P(x,0)=P(x,T),H(x,0)=H(x,T)\end{array}$$

(4)

Here, $F_L$ and $M_L$ represents the load condition at the free end, $V_0$ and ${\mathrm{\Omega }}_0$ are the velocity condition at the root of the blade. Let ( )⁰ represent the initial time state, thus the initial value boundary conditions of the transient problem can be similarly summarized as:

$$\begin{array}{l}F(L,t)=F_L,M(L,t)=M_L,V(0,t)=V_0,\mathrm{\Omega }(0,t)={\mathrm{\Omega }}_0\\ \gamma (x,0)={\gamma }^0,\kappa (x,0)={\kappa }^0,P(x,0)=P^0,H(x,0)=H^0\end{array}$$

(5)

3. Galerkin Space-Time Finite Element Method

According to the work in [19], the fully intrinsic equations is naturally adapted to the Galerkin space-time continuous finite element method (STFEM). But only the analysis of low order scheme is given, lacking the comparison of high order algorithm. This paper will construct a variable-order scheme of the STFEM.

3.1. Conservation Laws of the Galerkin Space-Time Equations

For a hingeless rotor blade, consider the following weak form of all the equations (fully intrinsic equation, space boundary condition and time boundary condition), premultiply $V^T$, ${\mathrm{\Omega }}^T$, $F^T$ and $M^T$ by the Galerkin weighted method, and integrate in the spanwise and time domain:

$$\begin{array}{l}{\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left\{V^T\left[\dot{P}+\tilde{\mathrm{\Omega }}P-F^{\prime }-(\tilde{k}+\tilde{\kappa })F-\mathrm{f}\right]\right.}}\\ \begin{array}{cc}& \end{array}+{\mathrm{\Omega }}^T\left[\dot{H}+\tilde{\mathrm{\Omega }}H+\tilde{V}P-M^{\prime }-(\tilde{k}+\tilde{\kappa })M-({\tilde{e}}_1+\tilde{\gamma })F-\mathrm{m}\right]\\ \begin{array}{cc}& \end{array}+F^T\left[\dot{\gamma }-V^{\prime }-(\tilde{k}+\tilde{\kappa })V-({\tilde{e}}_1+\tilde{\gamma })\mathrm{\Omega }\right]\left.+M^T\left[\dot{\kappa }-{\mathrm{\Omega }}^{\prime }-(\tilde{k}+\tilde{\kappa })\mathrm{\Omega }\right]\right\}dxdt\\ +{\displaystyle {\int }_0^T\left\{V^T(L,t)\left[F(L,t)-F^L\right]+{\mathrm{\Omega }}^T(L,t)\left[M(L,t)-M^L\right]\right.}\\ \begin{array}{cc}& \end{array}\left.+F^T(0,t)\left[V^0-V(0,t)\right]+M^T(0,t)\left[{\mathrm{\Omega }}^0-\mathrm{\Omega }(0,t)\right]\right\}dt\\ +{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^L\left\{V^T(x,0)\left[P(x,0)-P^0\right]+{\mathrm{\Omega }}^T(x,0)\left[H(x,0)-H^0\right]\right.}\\ \begin{array}{cc}& \end{array}\left.+F^T(x,0)\left[\gamma (x,0)-{\gamma }^0\right]+M^T(x,0)\left[\kappa (x,0)-{\kappa }^0\right]\right\}dx=0\end{array}$$

(6)

Note that the constitutive equations are already satisfied exactly. With the rules of tilde notation $\widetilde{( )}$, simplify the Eq. (6) as

$$\begin{array}{l}{\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left\{V^T\left[\dot{P}-F^{\prime }-\mathrm{f}\right]\right.}}+{\mathrm{\Omega }}^T\left[\dot{H}-M^{\prime }-\mathrm{m}\right]+F^T\left[\dot{\gamma }-V^{\prime }\right]\left.+M^T\left[\dot{\kappa }-{\mathrm{\Omega }}^{\prime }\right]\right\}dxdt\\ +{\displaystyle {\int }_0^T\left\{V^T(L,t)\left[F(L,t)-F^L\right]+{\mathrm{\Omega }}^T(L,t)\left[M(L,t)-M^L\right]\right.}\\ \begin{array}{cc}& \end{array}\left.+F^T(0,t)\left[V^0-V(0,t)\right]+M^T(0,t)\left[{\mathrm{\Omega }}^0-\mathrm{\Omega }(0,t)\right]\right\}dt\\ +{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^L\left\{V^T(x,0)\left[P(x,0)-P^0\right]+{\mathrm{\Omega }}^T(x,0)\left[H(x,0)-H^0\right]\right.}\\ \begin{array}{cc}& \end{array}\left.+F^T(x,0)\left[\gamma (x,0)-{\gamma }^0\right]+M^T(x,0)\left[\kappa (x,0)-{\kappa }^0\right]\right\}dx=0\end{array}$$

(7)

Integrating by parts in space,

$$\begin{array}{l}{\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left(V^T{F}^{\prime }+F^T{V}^{\prime }\right)dxdt}}={{\displaystyle {\int }_0^T\left.V^{T}F\right|}}_0^{L}dt={\displaystyle {\int }_0^T\left\{V^T(L,t)F(L,t)-F^T(0,t)V(0,t)\right\}dt}\\ {\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left({\mathrm{\Omega }}^T{M}^{\prime }+M^T{\mathrm{\Omega }}^{\prime }\right)dxdt}}={{\displaystyle {\int }_0^T\left.{\mathrm{\Omega }}^{T}M\right|}}_0^{L}dt={\displaystyle {\int }_0^T\left\{{\mathrm{\Omega }}^T(L,t)M(L,t)-M^T(0,t)\mathrm{\Omega }(0,t)\right\}dt}\end{array}$$

(8)

By the constitutive laws, integrating by parts in time,

$$\begin{array}{l}{\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left(V^T\dot{P}+{\mathrm{\Omega }}^T\dot{H}\right)dxdt}}={{\displaystyle {\int }_0^L\left.\frac{1}{2}\left\{\begin{array}{cc}{V}^T& {\mathrm{\Omega }}^T\end{array}\right\}\left\{\begin{array}{c}P\\ H\end{array}\right\}\right|}}_0^{T}dx={{\displaystyle {\int }_0^L\left.\mathrm{T}\right|}}_0^{T}dx\\ {\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left(F^T\dot{\gamma }+M^T\dot{\kappa }\right)dxdt}}={{\displaystyle {\int }_0^L\left.\frac{1}{2}\left\{\begin{array}{cc}{F}^T& M^T\end{array}\right\}\left\{\begin{array}{c}\gamma \\ \kappa \end{array}\right\}\right|}}_0^{T}dx={{\displaystyle {\int }_0^L\left.\mathrm{U}\right|}}_0^{T}dx\end{array}$$

(9)

Here, bring the kinetic energy and strain energy into evidence, which leads to:

$$\begin{array}{l}{\displaystyle {\int }_0^L{\left(\mathrm{T}+\mathrm{U}\right)}^{T}dx}-{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^L\left[V^T(x,0)P^0+{\mathrm{\Omega }}^T(x,0)H^0\right]dx}-{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^L\left[F^T(x,0){\gamma }^0+M^T(x,0){\kappa }^0\right]dx}\\ ={\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left[V^T\mathrm{f}+{\mathrm{\Omega }}^T\mathrm{m}\right]dxdt}}+{\displaystyle {\int }_0^T\left[V^T(L,t)F^L+{\mathrm{\Omega }}^T(L,t)M^L-F^T(0,t)V^0-M^T(0,t){\mathrm{\Omega }}^0\right]dt}\end{array}$$

(10)

The first term on the left side of the equation is the final energy, the second term is the initial kinetic energy, and the third term is the initial potential energy. The first term on the right of the equation is the work done by the distributed load, and the second term is the work done by the concentrated load at both ends of the beam. The equations states that the change of energy of the beam is equal to the work done on the beam. Thus, the above equations is an energy balance equations. The above convidence is valid in every space-time element. When dealing with periodic steady state problem, replace the time initial condition by periodic continuous condition, without breaking the conservation of equation.

3.2. Discretization in the Space-Time Domain

The integral domain of the above equation is then discretized as shown on the left in Figure 1. In the process of rotation, the blade with length L is dispersed into Ns segments along the spanwise. For any motion duration T which is multiple of the rotation period for the steady-state format, is dispersed into Nt segments along the azimuth Angle. Ultimately, form a global space-time cell grid on the scale of $Ns\times Nt$. The Eq. (6) must be satisfied in each element, and simultaneously convert the boundary conditions into continuity between adjacent elements.

Taking the element E_ij, which is ith(i=1~Ns) in space and jth (j=1~Nt) in time, shown on the right in Figure 2. The spacial length of this element is denoted as L_i, and the temporal length is denoted as T_j, while 0 represents the start position or the initial time. The variables within the element are denoted as Q_ij(x,t) that substitute for intrinsic variables $V_{ij}\left(x,t\right),{\mathrm{\Omega }}_{ij}\left(x,t\right),F_{ij}\left(x,t\right),M_{ij}\left(x,t\right),{\gamma }_{ij}\left(x,t\right),{\kappa }_{ij}\left(x,t\right),P_{ij}\left(x,t\right),H_{ij}\left(x,t\right)$among which (x,t) are coordinates in the element. The lines around the element are shared with the four elements around time (E_i,j-1, E_i,j+1) and space (E_i-1,j, E_i+1,j), that is, the continuous variables should be continuous with the adjacent elements or meet the boundary conditions. For a homogeneous blade, the continuity condition can be summarized as follows:

$$\begin{array}{c}{F}_{i,j}\left(L_i,t\right)=F_{i+1,j}\left(0,t\right),M_{i,j}\left(L_i,t\right)=M_{i+1,j}\left(0,t\right),V_{i,j}\left(L_i,t\right)=V_{i+1,j}\left(0,t\right),{\mathrm{\Omega }}_{i,j}\left(L_i,t\right)={\mathrm{\Omega }}_{i+1,j}\left(0,t\right)\\ {\gamma }_{i,j}\left(x,T_j\right)={\gamma }_{i,j+1}\left(x,0\right),{\kappa }_{i,j}\left(x,T_j\right)={\kappa }_{i,j+1}\left(x,0\right),P_{i,j}\left(x,T_j\right)=P_{i,j+1}\left(x,0\right),H_{i,j}\left(x,T_j\right)=H_{i,j+1}\left(x,0\right)\end{array}$$

(11)

The above continuity conditions are applicable to any element with concentrated mass and initial deformation, thus the space-time finite element equation is as following:

$$\begin{array}{l}{\displaystyle {\int }_0^{Tj}{\displaystyle {\int }_0^{Li}\left\{V^T\left[\dot{P}+\tilde{\mathrm{\Omega }}P-F^{\prime }-(\tilde{k}+\tilde{\kappa })F-\mathrm{f}\right]\right.}}\\ \begin{array}{cc}& \end{array}+{\mathrm{\Omega }}^T\left[\dot{H}+\tilde{\mathrm{\Omega }}H+\tilde{V}P-M^{\prime }-(\tilde{k}+\tilde{\kappa })M-({\tilde{e}}_1+\tilde{\gamma })F-\mathrm{m}\right]\\ \begin{array}{cc}& \end{array}+F^T\left[\dot{\gamma }-V^{\prime }-(\tilde{k}+\tilde{\kappa })V-({\tilde{e}}_1+\tilde{\gamma })\mathrm{\Omega }\right]\left.+M^T\left[\dot{\kappa }-{\mathrm{\Omega }}^{\prime }-(\tilde{k}+\tilde{\kappa })\mathrm{\Omega }\right]\right\}dxdt\\ +{\displaystyle {\int }_0^{Tj}\left\{V^T(L_i,t)\left[F(L_i,t)-F_{i+1}(0,t)\right]+{\mathrm{\Omega }}^T(L_i,t)\left[M(L_i,t)-M_{i+1}(0,t)\right]\right.}\\ \left.\begin{array}{cc}& \end{array}+F^T(0,t)\left[V_{i-1}(L_{i-1},t)-V(0,t)\right]+M^T(0,t)\left[{\mathrm{\Omega }}_{i-1}(L_{i-1},t)-\mathrm{\Omega }(0,t)\right]\right\}dt\\ +{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^{Li}\left\{V^T(x,0)\left[P(x,0)-P_{j-1}(x,T_{j-1})\right]+{\mathrm{\Omega }}^T(x,0)\left[H(x,0)-H_{j-1}(x,T_{j-1})\right]\right.}\\ \begin{array}{cc}& \end{array}\left.+F^T(x,0)\left[\gamma (x,0)-{\gamma }_{j-1}(x,T_{j-1})\right]+M^T(x,0)\left[\kappa (x,0)-{\kappa }_{j-1}(x,T_{j-1})\right]\right\}dx=0\end{array}$$

(12)

Divide the above equation into four equations by the preceding multiplicative term:

$$\begin{array}{l}{\displaystyle {\int }_0^{Tj}{\displaystyle {\int }_0^{Li}\left\{V^T\left[\dot{P}+\tilde{\mathrm{\Omega }}P-F^{\prime }-(\tilde{k}+\tilde{\kappa })F-\mathrm{f}\right]\right\}}}dxdt\\ +{\displaystyle {\int }_0^{Tj}\left\{V^T(L_i,t)\left[F(L_i,t)-F_{i+1}(0,t)\right]\right\}}dt+{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^{Li}\left\{V^T(x,0)\left[P(x,0)-P_{j-1}(x,T_{j-1})\right]\right\}}dx=0\end{array}$$

(13)

$$\begin{array}{l}{\displaystyle {\int }_0^{Tj}{\displaystyle {\int }_0^{Li}\left\{{\mathrm{\Omega }}^T\left[\dot{H}+\tilde{\mathrm{\Omega }}H+\tilde{V}P-M^{\prime }-(\tilde{k}+\tilde{\kappa })M-({\tilde{e}}_1+\tilde{\gamma })F-\mathrm{m}\right]\right\}}}dxdt\\ +{\displaystyle {\int }_0^{Tj}\left\{{\mathrm{\Omega }}^T(L_i,t)\left[M(L_i,t)-M_{i+1}(0,t)\right]\right\}}dt+{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^{Li}\left\{{\mathrm{\Omega }}^T(x,0)\left[H(x,0)-H_{j-1}(x,T_{j-1})\right]\right\}}dx=0\end{array}$$

(14)

$$\begin{array}{l}{\displaystyle {\int }_0^{Tj}{\displaystyle {\int }_0^{Li}\left\{F^T\left[\dot{\gamma }-V^{\prime }-(\tilde{k}+\tilde{\kappa })V-({\tilde{e}}_1+\tilde{\gamma })\mathrm{\Omega }\right]\right\}}}dxdt\\ +{\displaystyle {\int }_0^{Tj}\left\{F^T(0,t)\left[V_{i-1}(L_{i-1},t)-V(0,t)\right]\right\}}dt+{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^{Li}\left\{F^T(x,0)\left[\gamma (x,0)-{\gamma }_{j-1}(x,T_{j-1})\right]\right\}}dx=0\end{array}$$

(15)

$$\begin{array}{l}{\displaystyle {\int }_0^{Tj}{\displaystyle {\int }_0^{Li}\left\{M^T\left[\dot{\kappa }-{\mathrm{\Omega }}^{\prime }-(\tilde{k}+\tilde{\kappa })\mathrm{\Omega }\right]\right\}}}dxdt\\ +{\displaystyle {\int }_0^{Tj}\left\{M^T(0,t)\left[{\mathrm{\Omega }}_{i-1}(L_{i-1},t)-\mathrm{\Omega }(0,t)\right]\right\}}dt+{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^{Li}\left\{M^T(x,0)\left[\kappa (x,0)-{\kappa }_{j-1}(x,T_{j-1})\right]\right\}}dx=0\end{array}$$

(16)

The Eqs. (13), (14) are the dynamic equations of the element, and the Eqs. (15), (16) are the kinematic equations. When Eqs. (13)-(16) are both satisfied, Eq. (12) is established.

3.3. Separation of Variables in Space and Time

As for the dimensions of time and space are independent of each other, assume the intrinsic variables can be written as a linear combination of the products of a series of trial functions:

$$Q_{i,j}\left(x,t\right)=p\left({\overline{x}}_i\right)p\left({\overline{t}}_j\right)\overline{q}$$

(17)

Here, $\overline{q}$is the column vector of distinct polynomial weights, and ${\overline{x}}_i=x/L_i$, ${\overline{t}}_j=t/T_j$. The independent trial functions used are the shifted Legendre function [20] that are orthogonal with the domain [0,1], and can well fit the boundary and periodic conditions. The polynomials can be iteratively derived as:

$$p^0\left(\overline{x}\right)=1,p^1\left(\overline{x}\right)=2\overline{x}-1,p^{k+1}\left(\overline{x}\right)={\displaystyle \frac{(2k+1)(2\overline{x}-1)p^k(\overline{x})-i{p}^{k-1}(\overline{x})}{k+1}}$$

(18)

Therefore, the derivatives of variables can be expressed as the derivatives of spatial or temporal interpolation polynomials, respectively:

$${Q^{\prime }}_{i,j}\left(x,t\right)=p^{\prime }\left({\overline{x}}_i\right)p\left({\overline{t}}_j\right)\overline{q},{\dot{Q}}_{i,j}\left(x,t\right)=p\left({\overline{x}}_i\right)\dot{p}\left({\overline{t}}_j\right)\overline{q}$$

(19)

The continuity condition can be expressed as the continuity of the value of interpolating polynomial at ends (0 or 1) of the elements:

$$\begin{array}{l}{Q}_{i,j}\left(L_i,t\right)-Q_{i+1,j}\left(0,t\right)=p\left(1\right)p\left({\overline{t}}^j\right)\overline{q}-p\left(0\right)p\left({\overline{t}}^j\right){\overline{q}}_{i+1}\\ Q_{i,j}\left(x,T_j\right)-Q_{i,j+1}\left(x,0\right)=p\left({\overline{x}}^i\right)p\left(1\right)\overline{q}-p\left({\overline{x}}^i\right)p\left(0\right){\overline{q}}_{j+1}\end{array}$$

(20)

Let $f$, $m$, $v$and $\omega$ respectively represent the weights of corresponding variables $F$, $M$, $V$and $\mathrm{\Omega }$. According to the Eqs. (3), (19), (20), transform the element fully intrinsic equations (13) - (16) into an algebraic system about the weight coefficients:

$$\begin{array}{l}{\displaystyle {\int }_0^{Tj}{\displaystyle {\int }_0^{Li}\left\{p^a{p}^b\left[p^k{\dot{p}}^l\left(Gv+K\omega \right)+\left(p^p{p}^q\tilde{\omega }\right)p^k{p}^l\left(Gv+K\omega \right)\right.\right.}}\\ \left.\left.\begin{array}{cc}& \end{array}-\left(\tilde{k}+p^p{p}^q\left(\widetilde{Z^{T}f}+\widetilde{Tm}\right)\right)p^k{p}^{l}f-{p^{\prime }}^k{p}^{l}f-\mathrm{f}\right]\right\}dxdt\\ +{\displaystyle {\int }_0^{Tj}\left\{p^a\left(1\right)p^b\left[p^k\left(1\right)p^{l}f-p^k\left(0\right)p^l{f}_{i+1}\right]\right\}}dt\\ +{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^{Li}\left\{p^a{p}^b(0)\left[p^k{p}^l\left(0\right)\left(Gv+K\omega \right)-p^k{p}^l\left(1\right)\left(G{v}_{j-1}+K{\omega }_{j-1}\right)\right]\right\}}dx=0\end{array}$$

(21)

$$\begin{array}{l}{\displaystyle {\int }_0^{Tj}{\displaystyle {\int }_0^{Li}\left\{p^a{p}^b\left[p^k{\dot{p}}^l\left(K^{T}v+I\omega \right)+\left(p^p{p}^q\tilde{\omega }\right)p^k{p}^l\left(K^{T}v+I\omega \right)\right.\right.}}\\ \begin{array}{cc}& \end{array}+\left(p^p{p}^q\tilde{v}\right)p^k{p}^l\left(Gv+K\omega \right)-\left(\tilde{k}+p^p{p}^q\left(\widetilde{Z^{T}f}+\widetilde{Tm}\right)\right){p_x}^k{p_t}^{l}m\\ \begin{array}{cc}& \end{array}\left.\left.-\left({\tilde{e}}_1+p^p{p}^q\left(\widetilde{Rf}+\widetilde{Zm}\right)\right)p^k{p}^{l}f-{p^{\prime }}^k{p}^{l}m-\mathrm{m}\right]\right\}dxdt\\ +{\displaystyle {\int }_0^{Tj}\left\{p^a\left(1\right)p^b\left[p^k\left(1\right)p^{l}m-p^k\left(0\right)p^l{m}_{i+1}\right]\right\}}dt\\ +{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^{Li}\left\{p^a{p}^b(0)\left[p^k{p}^l\left(0\right)\left(K^{T}v+I\omega \right)-p^k{p}^l\left(1\right)\left(K^T{v}_{j-1}+I{\omega }_{j-1}\right)\right]\right\}}dx=0\end{array}$$

(22)

$$\begin{array}{l}{\displaystyle {\int }_0^{Tj}{\displaystyle {\int }_0^{Li}\left\{p^a{p}^b\left[p^k{\dot{p}}^l\left(Rf+Zm\right)\right.\right.}}-\left(\tilde{k}+p^p{p}^q\left(\widetilde{Z^{T}f}+\widetilde{Tm}\right)\right){p_x}^k{p_t}^{l}v\\ \left.\left.\begin{array}{cc}& \end{array}-\left({\tilde{e}}_1+p^p{p}^q\left(\widetilde{Rf}+\widetilde{Zm}\right)\right)p^k{p}^l\omega -{p^{\prime }}^k{p}^{l}v\right]\right\}dxdt\\ +{\displaystyle {\int }_0^{Tj}\left\{p^a\left(0\right)p^b\left[p^k\left(1\right)p^l{v}_{i-1}-p^k\left(0\right)p^{l}v\right]\right\}}dt\\ +{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^{Li}\left\{p^a{p}^b(0)\left[p^k{p}^l\left(0\right)\left(Rf+Zm\right)-p^k{p}^l\left(1\right)\left(R{f}_{j-1}+Z{m}_{j-1}\right)\right]\right\}}dx=0\end{array}$$

(23)

$$\begin{array}{l}{\displaystyle {\int }_0^{Tj}{\displaystyle {\int }_0^{Li}\left\{p^a{p}^b\left[p^k{\dot{p}}^l\left(Z^{T}f+Tm\right)-{p^{\prime }}^k{p}^l\omega -\left(\tilde{k}+p^p{p}^q\left(\widetilde{Z^{T}f}+\widetilde{Tm}\right)\right)p^k{p}^l\omega \right]\right\}}}dxdt\\ +{\displaystyle {\int }_0^{Tj}\left\{p^a\left(0\right)p^b\left[p^k\left(1\right)p^l{\omega }_{i-1}-p^k\left(0\right)p^l\omega \right]\right\}}dt\\ +{\displaystyle \frac{1}{2}}{\displaystyle {\int }_0^{Li}\left\{p^a{p}^b(0)\left[p^k{p}^l\left(0\right)\left(Z^{T}f+Tm\right)-p^k{p}^l\left(1\right)\left(Z^T{f}_{j-1}+T{m}_{j-1}\right)\right]\right\}}dx=0\end{array}$$

(24)

Here, a lot of subscripts (ab, kl, pq) are used to distinguish the order of polynomial. Especially, each group of subscripts corresponds to a different variable. And within each group, one is of temporal interpolation order between 0~Rs and the other is of spatial interpolation order between 0~Rt. Furthermore, during the processing of the premultiplication term, each equation is required to satisfy by any $f^T$, $m^T$, $v^T$and ${\omega }^T$, so it can be mentioned outside the equations.

3.4. Assemble Element Algebraic Equations

Take one weight coefficient out of the integral, and sort the integrals in Eqs. (21)~(24): the terms without variables in the integral are classified as linear matrices $S_{ij}$, the cross terms with variables in the integral are classified as nonlinear matrices $C_{ij}(q_{ij})$, the terms without variables outside the integral are classified as constant matrices $D_{ij}$, and the variables outside the integral belong to adjacent elements $q_{i-1,j}$, $q_{i+1,j}$ and $q_{i,j-1}$ are classified as continuity matrices that contain spatial continuity matrices $H_1^{i-1,j}$ and $H_2^{i+1,j}$, and time continuity matrices $h_t^{i,j-1}$. Thus, the element equation can be summarized in the following matrix form:

$$S_{ij}{q}_{ij}+C_{ij}(q_{ij})\cdot q_{ij}+D_{ij}+H_1^{i-1,j}{q}_{i-1,j}+H_2^{i+1,j}{q}_{i+1,j}+h_t^{i,j-1}{q}_{i,j-1}=0$$

(25)

Here, $q_{ij}={\left\{\begin{array}{cccc}v& \omega & f& m\end{array}\right\}}^T$ is the column vector of $12Rs\cdot Rt$ dimension.

The system computing model of the structure is composed of all the space-time elements, and the structural characteristic matrix of the system is composed of superposition of each element coefficient matrix. On the basis of obtaining the above element matrix equation, then assemble the beam system equation according to the order of first space next time, as shown in Figure 2.

First of all, the equation variables are integrated: the jth time node variables $q_j$ are obtained by arranging the space-time element variables $q_{ij}$ according to the spatial sequence of 1~Nt. Then the system unknowns $q$ are obtained by arranging different time node variables according to time sequence of 1~Ns.

Then, Assemble the coefficient matrix of each element according to the arrangement logic of variables:

(1)

Total stiffness matrix$B$: The spatial continuity matrix $H_1^{i-1,j}$, $H_2^{i+1,j}$ and the linear matrix $S_{ij}$ are arranged side by side in the order of space formed as $\left[\begin{array}{ccc}{H}_1^{i-1,j}& S_{ij}& H_2^{i+1,j}\end{array}\right]$, then put it at row $i$ and colume $\left[\begin{array}{ccc}i-1& i& i+1\end{array}\right]$ to form the time node matrix $S_j$. The element time continuity matrix $h_t^{i,j-1}$ is arranged on the diagonal in spatial order to form the time node continuity matrix $H_t^{j-1}$. The time correlation matrix $H_t^{j-1}$ and $S_j$ are arranged in the order of time formed as $\left[\begin{array}{cc}{H}_t^{j-1}& S_j\end{array}\right]$, then put it at row $j$ and colume $\left[\begin{array}{cc}j-1& j\end{array}\right]$ to form the total stiffness matrix.

(2)

Total nonlinear matrix $C(q)$: The element nonlinear matrix $C_{ij}(q_{ij})$ is arranged diagonally in space order formed as $C_j(q_j)$, and then diagonally in time order.

Finally, the constant matrix $D$ is assembled according to the assembly logic of the variables: first, the element constant matrix $D_{ij}$ is placed in the row $i$ to form the time constant matrix $D_j$, and then it is placed in the row $j$ to form the total constant matrix that includes external loads and boundary conditions.

So the assembled algebraic equations can be written as follows:

$$Bq+C\left(q\right)\cdot q+D=0$$

(26)

It should be noted that the time continuity matrix $H_t^m$ on the corner of the stiffness matrix $B$ in the initial value scheme is zero and the initial value condition will be contained in the constant matrix. When the number of elements is 1, the whole domain is taken as a element, eliminating the process of discretization and assembly, but the variables in the domain are required to be consistent, that is, the profile characteristics are continuous, and in the periodic scheme, the matrix $H_t^m$ is implicit in the matrix $S_1$.

Newton-Raphson method is used to solve the nonlinear system equations with the Jacobian matrix that can be expressed accurately. Here's a simple demonstration, supposing the Eq. (26) is

$$B\left[\begin{array}{c}F\\ V\end{array}\right]+\left[\begin{array}{cc}\tilde{V}& \\ & 0\end{array}\right]\left[\begin{array}{c}F\\ V\end{array}\right]+D=0$$

(27)

Then the Jacobian matrix is

$$Jac=B+\left[\begin{array}{cc}\tilde{V}& -\tilde{F}\\ & 0\end{array}\right]$$

(28)

The iterative initial value is taken by dividing the blade stiffness matrix with the constant matrix, i.e $q_0=-D/B$. According to the result of each step, the Jacobian matrix is updated and the variables are iterated until convergence. The result is substituted back to Eq. (17) to get the space and time distribution function of the intrinsic variables.

After obtaining the global distribution of the intrinsic variables and rotation matrix, it is still difficult to imagine the displacement and rotation angle of the beam, so it is necessary to convert the results to the configuration space. For this purpose, the displacement can be reversed by the following displacement-strain equation and displacement-velocity equation:

$$\begin{array}{l}{u}^{\prime }=C^T\left(e_1+\gamma \right)-e_1-\tilde{k}u\\ \dot{u}=C^{T}V-\mathrm{v}-\tilde{\omega }u\end{array}$$

(29)

Here, $u$ is the displacement in the undeformed coordinate system and $\mathrm{v}$ is the linear velocity of the undeformed blade. Similarly, the theta of rotation can be reversed by rotation kinematical equations which is:

$$\begin{array}{l}{C}^{\prime }=C\tilde{k}-(\tilde{k}+\tilde{\kappa })C\\ \dot{C}=C \mathrm{w} \text{˜}-\tilde{\mathrm{\Omega }}C\end{array}$$

(30)

Here, $\mathrm{w}$is the angular velocity of the undeformed blade, and $C$is a 3x3-dimensional tensor matrix, which represents the rotation of the section from the undeformed coordinate to the deformed coordinate with $C·C^T=\mathrm{\Delta }$. The initial value conditions of displacement and rotation are:

$$\begin{array}{l}u(0,t)=U_0,u(x,0)=U^0,\\ C(0,t)=C_0,C(x,0)=C^0\end{array}$$

(31)

To represent the angle $\theta ={\left[\begin{array}{ccc}{\theta }_1& {\theta }_2& {\theta }_3\end{array}\right]}^T$ in the undeformed coordinate system, according to the specific expression of the coordinate transformation matrix which is no singularity, it is obtained:

$$C=\left[\begin{array}{ccc}1& & \\ & \cos {\theta }_1& \sin {\theta }_1\\ & -\sin {\theta }_1& \cos {\theta }_1\end{array}\right]\left[\begin{array}{ccc}\cos {\theta }_2& & -\sin {\theta }_2\\ & 1& \\ \sin {\theta }_2& & \cos {\theta }_2\end{array}\right]\left[\begin{array}{ccc}\cos {\theta }_3& \sin {\theta }_3& \\ -\sin {\theta }_3& \cos {\theta }_3& \\ & & 1\end{array}\right]$$

(32)

3.5. Flow of Space-Time Finite Element Method

The basic flow of space-time finite element method (STFEM) for fully intrinsic geometrically exact beam equations can be summarized as follows:

• Galerkin space-time equation for a hingless blade is Eq. (6). Discrete the space along the beam, and discrete the time along the azimuth angle in each movement period, then obtain the global discrete space-time grid as Figure 1;
• Establish the continuous relationships between units as Eqs. (11). Thus the space-time finite element equation can be obtained by the weighted method as Eq. (12), divide it into four equations according to the different weight functions like Eqs. (13) - (16).
• With the Legendre functions Eq. (18), approximate the variables to a linear combination of a set of spatio-temporal polynomials as Eq. (19). Then variables are only the weight coefficients of polynomials and transform the derivative terms and continuity into the derivative and continuity of polynomials as Eq. (20). Thus the finite element algebraic equations are obtained as Eqs. (21)~(24).
• Take one weight coefficient out of integral, classify the integral functions according to the number of variables, then the matrix equation Eq. (25) of element is obtained.
• Assemble the above element matrix equation according to a certain law as show in Figure 2. Then the global discrete algebraic equation Eq. (26) of beam are obtained whose coefficient matrix only includes stiffness matrix, nonlinear matrix and constant matrix;
• Gain the linear initial value as $q_0=-D/B$ and the Jacobian matrix like Eqs. (27), (28). The Newton method is used to solve the nonlinear equation, and the global distribution of the intrinsic values is sorted out.
• Given the initial conditions of displacement and rotation angle like Eq. (31), obtain the displacement distribution $u$ in the entire domain by Eq. (29) and the transformation matrix distribution $C$ by Eq. (30), then get the theta distribution $\theta$ by Eq. (32).

3.6. Treatment of Blade Structure Discontinuity

It should be noted that in the actual design process, there is a discontinuity in the blade structure, that is, the section structural characteristics are discontinuous which lead to a discontinuity of load. According to the constitutive law, the fact that the strain is continuous in the spatial direction does not mean that the load is uniformly continuous:

$$\left\{\begin{array}{c}F\\ M\end{array}\right\}=\left[S\right]\left\{\begin{array}{c}\gamma \\ \kappa \end{array}\right\}\Rightarrow {\left\{\begin{array}{c}F\\ M\end{array}\right\}}^{-}-{\left\{\begin{array}{c}F\\ M\end{array}\right\}}^{+}={\left[S\right]}^{-}{\left\{\begin{array}{c}\gamma \\ \kappa \end{array}\right\}}^{-}-{\left[S\right]}^{+}{\left\{\begin{array}{c}\gamma \\ \kappa \end{array}\right\}}^{+}=\left({\left[S\right]}^{-}-{\left[S\right]}^{+}\right){\left\{\begin{array}{c}\gamma \\ \kappa \end{array}\right\}}^{-}$$

(33)

Where, $\mathrm{（}\begin{array}{c}\end{array}{\mathrm{）}}^{-}$and $\mathrm{（}\begin{array}{c}\end{array}{\mathrm{）}}^{+}$ represents the left and right parameters of the discontinuous section, [S] is the profile stiffness matrix. Therefore, the spatial continuity of the load at the discontinuity section can be established by the section flexibility matrix [R] as:

$$\left\{\begin{array}{c}\gamma \\ \kappa \end{array}\right\}=\left[R\right]\left\{\begin{array}{c}F\\ M\end{array}\right\}\Rightarrow {\left\{\begin{array}{c}F\\ M\end{array}\right\}}^{+}={\left[S\right]}^{+}{\left\{\begin{array}{c}\gamma \\ \kappa \end{array}\right\}}^{+}={\left[S\right]}^{+}{\left\{\begin{array}{c}\gamma \\ \kappa \end{array}\right\}}^{-}={\left[S\right]}^{+}{\left[R\right]}^{-}{\left\{\begin{array}{c}F\\ M\end{array}\right\}}^{-}$$

(34)

When the left and right characteristics of the section are consistent, Eq. (32) naturally degenerates into Eq. (3).

3.7. Handling of Conservative Loads

In the process of rotor movement, the main external force is the aerodynamic load following the deformation, which completely fits to the coordinate system of the fully intrinsic model. In addition, the blade is also affected by gravity, although the value is small compared with the aerodynamic force, but in many cases can’t be ignored.

Refer to Ref. [21], a coupling scheme of geometrically exact equations combined kinematical equation was proposed by Lagrange multipliers. So that all the variables can vary independently in the calculus. This paper introduces the conservative loads through the coupling of rotation kinematical equations. Thus, by coupling Eq. (30) with Eq. (1) through Lagrange multipliers, it can be obtained:

$$\begin{array}{l}{\displaystyle {\int }_0^T{\displaystyle {\int }_0^L\left\{V^T\left[\dot{P}+\tilde{\mathrm{\Omega }}P-F^{\prime }-(\tilde{k}+\tilde{\kappa })F-C \mathrm{f} \text{¯}\right]\right.}}\\ \begin{array}{cc}& \end{array}+{\mathrm{\Omega }}^T\left[\dot{H}+\tilde{\mathrm{\Omega }}H+\tilde{V}P-M^{\prime }-(\tilde{k}+\tilde{\kappa })M-({\tilde{e}}_1+\tilde{\gamma })F-C \mathrm{m} \text{¯}\right]\\ \begin{array}{cc}& \end{array}+F^T\left[\dot{\gamma }-V^{\prime }-(\tilde{k}+\tilde{\kappa })V-({\tilde{e}}_1+\tilde{\gamma })\mathrm{\Omega }\right]+M^T\left[\dot{\kappa }-{\mathrm{\Omega }}^{\prime }-(\tilde{k}+\tilde{\kappa })\mathrm{\Omega }\right]\\ \begin{array}{cc}& \end{array}\left.+{\lambda }_1^T\left[C^{\prime }-C\tilde{k}+(\tilde{k}+\tilde{\kappa })C\right]+{\lambda }_2^T\left[\dot{C}-C \mathrm{w} \text{˜}+\tilde{\mathrm{\Omega }}C\right]\right\}dxdt=0\end{array}$$

(35)

Here, $\mathrm{f} \text{¯}$ and $\mathrm{m} \text{¯}$ are the conservative force and moment. The Lagrange multipliers ${\lambda }_i$ always can be identified by conservation of energy and it is a three-dimensional column vector that can be deal with Eq. (17). Since ${\lambda }_i$ can be independently and arbitrarily varied, its specific form does not affect the process of the method proposed in this paper. Similarly, make modifications to the load boundary conditions with conservative tip force ${\overline{F}}_L$ and tip moment ${\overline{M}}_L$:

$$F(L,t)-C{\overline{F}}_L=0,M(L,t)-C{\overline{M}}_L=0$$

(36)

4. Numerical Results

In order to verify the effectiveness of the algorithm presented in this paper, the influence of higher order scheme and the ability to solve large deformation problem under conservative load, four cases are presented. Firstly, investigated by applying the concentrated tip moment with different order or elements. Secondly, non-conservative tip force and conservative tip force are applied on the cantilever beam. Thirdly, Minguet beam experiment were conducted to verify the ability under conservative load with common structural coupling. Finally, the dynamic response problem of rotating or non-rotating beam proves the validity of the initial value scheme.

4.1. Large Deformation Simulation of Bending Moment on Elastic Beam Tip

The first case is a benchmark problem for the geometrically nonlinear analysis of beams [22]. A flapping moment $M_L$ is applied to the tip of an elastic beam clamped at root. The parameter of the beam is $L=20ft$, $E{I}_2=9000lb\cdot f{t}^2$. Each point on the beam will be located on the arc with curvature $\rho =M_L/E{I}_2$, and its exact solution can be expressed as:

$$r_1=\rho \sin {\displaystyle \frac{s}{\rho }},r_3=\rho \left(1-\cos {\displaystyle \frac{s}{\rho }}\right)$$

(37)

Under the tip bending moment, only the bending moment exists at each point of the fixed beam, and the value is equal to the tip bending moment. This means the variables are only constants $m_0$ that the Legendre interpolation order is not sensitive to this problem and is only relevant to the algorithm error. The displacement is obtained by solving the approximate solution of the differential equation, whose error is closely related to the deformation degree and mesh density.

Table 1 list the tip position r3 under different tip bending moments with single one element, and Table 2 list the tip position r3 under different number of elements with a tip bending moments of 2500lb-ft2. All spatial elements adopt 9th-order interpolation. Through comparison, it is obvious that with the increase of deformation degree, the precision decreases obviously, and increasing the number of mesh can also improve the accuracy of the result. In Ref. [23], use the 20-element RCAS model, the magnitude of error is 6.00×10^-4 when the tip bending moment is 500lb-ft2, which is far less than the calculation accuracy of a single element. While the magnitude of error is 1.98×10^-4 when the tip bending moment is 2500lb-ft2, which can be achieved with only 3 elements.

Table 2 compares the differences between the interpolation functions order with the same number of elements. We can observe that as the element number increases, the two schemes tend to unify the solution, but when the number of units is less than 10 units, the error of the higher order scheme is significantly smaller.

4.2. Static Deflections of Cantilever with Conservative or Non-Conservative Load

Case 2 is the static deformation problem of a fixed beam under dynamic force [24,25]. The dynamic load is divided into two cases: the non-conservative and conservative force at the free end of the blade. The so-called non-conservative force is the direction of the force always perpendicular to the beam axis, while conservative force keeps in the original directions. The basic parameters of the beam are: length L=100cm, moment of inertia I₂₂=1.666667cm⁴, material Young's modulus E=2.1x10⁷N/cm².

The displacement-force curve of the free end of the beam under different non-conservative tip forces is shown in Figure 3, in which an obvious nonlinear large deformation of the beam can be seen. The simulation results of 5th-order 1-unit, 5th-order 2-units and 9th-order 1-unit are compared with the calculated data in Ref. 24 which is obtained by 10 elements and 30unknows. It can be noted that both the 5th-order 2-elements and the 9th-order 1-element are basically consistent with the reference results, except that the 5th-order 1 element deviates far under large deformation.

At the tip load of about 30kN, u1 deviates from the reference. According to the trend analysis, at this time, the beam is developing from linear to nonlinear deformation, the tip is gradually moving upward, the angle growth trend is gradually slowing down. So the tip should still be extending negatively, that is, the tip has passed the origin along the X-axis, and -u1/L is greater than 1. In addition, comparing the 5th-order 2-element and 9th-order 1-element, it can be found that when the beam is under large deformation, increasing the number of elements is closer to the result than increasing the order. This is because the displacement nonlinearity is more prominent than the load nonlinearity, and the mesh encryption can be more effective in fitting.

The displacement-force curve of the free end of the beam under different conservative tip loads is shown in Figure 4. It can be seen from the figure that both the 9th-order 1-unit and the 5th-order 2-unit are basically consistent with the reference results. Although there is a certain gap between the u3 of the 5th-order 1-unit and the compared data, the time error is not obvious when PL²/EI<6.

4.3. Static Behavior of Composite Blades Under Large Deflections

Composite materials have emerged as a primary material for helicopter rotor blades that leads to the presence of structural couplings such as extension-twisting and bending-twisting. These couplings can be used to improve the dynamic and aero-elastic performance of helicopter rotors. Case 3 is the large deformation problem of composite blade under static load [26,27], and examines the performance of the proposed algorithm for structurally coupled composite blade under large load conditions.

During the experiment, an aluminum chuck was used to fix the composite beam on the rotating platform to control the installation angle. The blade exposed outside the chuck was 550mm length and 30mm width. The strain gauge was attached at 50mm outside the chuck, and the displacement was observed at 300mm and 500mm. When the installation angle is 0°, the blade is placed horizontally. The STFEM with the 9-order 1-element format is used, compared with the finite difference method (FDM) in Ref. [28,29] which discretizes length into 56 nodes. In the experiments, the beam's own weight is included and all displacements are measured from initial position.

The first specimen has a [0/90]3s layup, without any coupling to serve as a reference, as shown in Figure 5. Among them, the displacement w is almost linear, and the displacement u is negative in practice. When the installation Angle is 0°, the displacement v is 0, and the tip deformation is about 40% of its length under maximum load. When the installation angle is turned to 45°, due to the weak rigidity in the normal plane, the deformation tends to be the initial v and w are basically the same. The calculation results of the method in this paper are basically consistent with those in the literature, which accords with the nature law and proves the reliability of the algorithm in the case of no coupling.

The next one is a [45/0]3s, to illustrate the effect of bending-twisting coupling, as shown in Figure 6. Due to the existence of flexural and torsional coupling, even at 0° installation Angle, there is a certain coupling deformation v, and it appears as v at 45°, and w is basically the same, but it is completely different at -45°. In general, although the algorithm in this paper has some shortcomings in u calculation, it is more close to the coupling change of v and w.

Meanwhile, a [20/-70/ -70/20]2a layup was used to illustrate the effect of extension-twist, as shown in Figure 7. Compared with uncoupled layup, the torsional stiffness was significantly enhanced. The results show that although there is a certain gap between the proposed algorithm and the experimental data, it is completely consistent with the results of the FDM.

4.4. Dynamic Response of Cantilever Beam

Two cases for dynamic response [30] are performed for hingeless boundary conditions with nonrotating or rotating to simulate in time domain of helicopter blades. The first one is a nonrotating beam with a tip dynamic forces along vertical orientation F₃=10sin(20t). And the second one is a rotating beam under a tip force along transversal orientation F₂=10sin(20t) with the rotating speed is 70 rad/s. The material properties of this beam are shown in Table 3.

Figure 8, Figure 9, Figure 10 and Figure 11 shows the time history curve of tip displacement, tip rotation, root force and root bending moment of a non-rotating cantilever beam under a given dynamic load F₃=10sin(20t). The algorithm in this paper (STFEM) adopts a 5x5 cell, 5x5 order space-time format, and computes in a step of 0.3s. The calculation takes 4 steps, totaling 1.2s. Results Compared with the Dymore program, the waveform is basically the same, though there is a certain time delay with the development. On the whole, the calculated results in this paper are completely consistent with Ansys, which proves the accuracy and reliability of the proposed method.

Figure 12, Figure 13, Figure 14 and Figure 15 show the time history curve of tip displacement, tip rotation and root bending moment of a rotating clamped beam under a given tip dynamic load F2=10sin(20t). The format used is the same as above, the time step is 0.045s, and the total of 8 steps is 0.36s. According to the results, each peak and valley position can be clearly captured. Although the deviation in the pendulum torque is slightly large compared with other terms, the whole is consistent with the simulation code Dymore, which verifies the feasibility of the proposed algorithm in the aerodynamic analysis of blades.

5. Conclusions

In this paper, a space-time continuous finite element method for solving the fully intrinsic equations of geometrically exact beam is established successfully. The process of the space-time discretization, establishment of element equation, discretization of variables, assembly of equations and solution are described in detail. This space-time finite element scheme is energy conservation and variable order that has the potential of adaptive algorithm. This paper also provides solutions to common discontinuity problem, conservative load problem and initial value problem. The results show that this algorithm has good precision, stable convergence and ability to analyze the problem of blade coupling. By increasing the interpolation order of variables, the calculation accuracy is improved effectively, and the minimum number requirement of elements is significantly reduced that can even be developed into a meshless algorithm. In the following studies, the proposed method can naturally introduce aerodynamic equations to achieve aero-elastic adjoint solution, or as a high-precision dynamics solver, use loose coupling method to analyze the aero-elastic response of the rotor.