Differential Quadrature Method for Fully Intrinsic Equations of Geometrically Exact Beams

Abstract

In this paper, a differential quadrature method of high-order precision (DQ−Pade), which is equivalent to the generalized Pade approximation for approximating the end of a time or spatial interval, is used to solve nonlinear fully intrinsic equations of beams. The equations are a set of first-order differential equations with respect to time and space, and the explicit unknowns of the equations involve only forces, moments, velocity and angular velocity, without displacements and rotations. Based on the DQ−Pade method, the spatial and temporal discrete forms of fully intrinsic equations were derived. To verify the effectiveness and applicability of the proposed method for discretizing the fully intrinsic equations, different examples available in the literatures were considered. It was found that when using the DQ−Pade method, the solutions of the intrinsic beam equations are obviously superior to those found by some other usual algorithms in efficiency and computational accuracy.

1. Introduction

Due to the high aspect ratio of a rotor blade, which is usually regarded as a beam-like structure, beam theory has been a very important part of helicopter rotor dynamics modeling and analysis. The blade will have medium or large deflection in complex working environments, which leads to the nonlinear characteristics of the aeroelastic response analysis of the helicopter rotor. Among the many beam theories describing these nonlinear behaviors, the large deformation beam theory has attracted the most attention. It has only a small strain assumption and no restriction on displacements or rotations.

At present, the large deformation beam theory has been applied in many fields, and lots of beam models have been developed. They can be classified according to their solution methods. It was found that there are three categories: displacement-based methods [1], mixed-form methods [2,3] and stress-based methods [4,5]. The main advantage of displacement-based methods is the ease of applying geometric boundary conditions with them. However, they are computationally expensive due to the higher-order nonlinear terms. Even with high-order truncation, the complete formula has to be written over several pages. Mixed-form methods establish the kinematics equation by mixing the rotation and displacement variables with the generalized velocity and strain measures. Lagrange multipliers are used to apply the constitutive and kinematic relations to the governing equations in the mixed formula. Differently from the displacement-based theory and the mixed formula theory, there is no displacement, nor any rotation variables, in the stress-based formula. Green and Laws [6] coined the term “fully intrinsic” for this concept. The fully intrinsic theory was originally proposed by Hegemier and Nair [7]. In 2003, Hodges [4] reconsidered it on the basis of [2] and developed the fully intrinsic equations of geometrically exact beams, which only include forces, moments, velocity and angular velocity. The displacements and rotations can be obtained by recovering them in the post-processing. Due to no truncation in the derivation process, these compact equations are complete and accurate, and the maximum nonlinearity (the highest order of the nonlinear terms in the equations) is only two. Palacios and Cesnik [8] compared the numerical efficiency and ease of use of these three different methods for the analysis of aircraft structure, aeroelastic characteristics, and flight dynamics.

Patil and Hodges [9] analyzed the nonlinear aeroelastic characteristic of wings with a high aspect ratio using the fully intrinsic method. Patil and Althoff [10] presented an energy-continuous Galerkin method for solving the fully intrinsic equations, which can obtain a high precision solution with a low computational cost. However, the Galerkin method was also found to face some difficulties in practical applications, such as time-consuming integration and an inability to deal with spatial discontinuity. Patil and Hodges [11] proposed the concept of the variable-order finite-element Galerkin method. They found that the third-order finite element offers a particularly good balance of accuracy, computational efficiency, and applicability to general problems. Sotoudeh and Hodges [12] developed an incremental method for solving the fully intrinsic equations of statically indeterminate structures. They [13] also extended the application of fully intrinsic theory to more general configurations by establishing different types of fully intrinsic boundary conditions. Khaneh, Masjedi, and Ovesy [14,15] investigated the static analysis method by discretizing the fully intrinsic equations using the Chebyshev collocation method. Tashaorian et al. [16] established non-local, fully intrinsic equations for non-beam-like structures whose one-dimensional structure length is close to the two-dimensional cross-section size, which breaks through the limitation that the fully intrinsic equations can only be used for modeling slender structures. Byeonguk et al. [17] used the exponential function to replace the Rodrigues parameters to represent the finite rotation, and realized the absolute non-singularity expression in the equations, and therefore, the numerical realization of solving the equations became more compact.

The fully intrinsic equations developed by Hodges are the differential equations of first-order partial derivatives in time and space. Usually, two steps are needed to discretize them: the first is to discretize the equations of motion in the spatial domain to obtain the first-order partial differential equations about time; the second is to discretize the equations of motion in the time domain to obtain the algebraic equations, which can be solved directly. There are four commonly used methods for spatial discretization in the literature which the reader can consult. The first is the finite difference method. Hodges [2] used the central difference scheme to discretize the equations in the spatial domain, and mentioned that this method has second-order accuracy. The second is the Galerkin method, which was presented by Patil [10] using Legendre polynomials as weighting functions to establish integral equations. The third method is the Chebyshev collocation method. Khaneh, Masjedi, and Ovesy [14,15] used this method to discretize static intrinsic equations. The last method is the differential quadrature (DQ) method that has been used more recently.

It was proposed by Bellman and Casti [18] in 1971. The basic principle is that the derivative of a function at a specific point can be approximated as a weighted linear summation of the function values at all discrete points in the entire domain. Bert [19] first applied the method to the mechanical analysis of structural elements, and found that the method has the advantages of a simple formula, convenient use, less computation, and high precision. However, with the deepening of research and applications, researchers found that when the number of discrete points is large, the numerical ill-condition will appear when calculating the weighting coefficients. This disadvantage limited the application of the DQ method. Shu and Richards [20,21] overcame the shortcoming by introducing the generalized differential quadrature (GDQ) method based on Lagrange polynomial vector space. By using the GDQ method, the weighting coefficients of arbitrary-order derivatives are calculated using simple algebraic formulas. A notable feature of the GDQ method is that there are no limits to the type and number of discrete points. Oleg et al. [22] constructed a multi-layer wavelet collocation method which uses the Gaussian wavelet function as the interpolation basis function, and the collocation points are selected by dichotomy. An adaptive multi-layer wavelet configuration method [23] and a fast, adaptive wavelet configuration method [24] were presented by the same group. Chen [25] proposed the concept of the differential quadrature element method, which discretizes the domain into several regular sub-domains to deal with irregular boundaries. Wei et al. [26] found that the differential quadrature element method is more effective at dealing with higher-order derivative problems. Fung [27,28] pointed out that the accuracy and stability of the DQ method in the time domain depend on the locations of the sampling discrete points. By using appropriately distributed discrete points, the accuracy of the end value of the time interval can be improved to 2N-1 or 2N (N is the maximum order of polynomials), which is better than the commonly used uniformly spaced point and Chebyshev–Gauss–Lobatto point distribution. They proved that this terminal approximation is equivalent to the generalized Pade approximation.

Amoozgar and Shahverdi [29] firstly applied the differential quadrature method to solve the fully intrinsic equations of a geometrically exact beam. They found that the GDQ method can obtain more accurate results with a lesser computational cost than the traditional finite element method. Based on GDQ method, they also analyzed the dynamic stability of the beam under following forces [30] and the aeroelastic stability of the hingeless rotor [31]. Tashaorian et al. [16] used the GDQ method for the spatial discretization of non-local fully intrinsic equations.

For the first-order partial differential equations with respect to time obtained after the space discretization of the fully intrinsic equations, the difference algorithms [2,10,32] are commonly used to solve them recursively. These algorithms often need to go through a time-consuming recursion process when solving the aeroelastic response of the rotor, especially when solving the steady-state periodic response of the rotor. Borri [33] also pointed out the shortcomings of these algorithms in solving the steady-state response of the rotor’s aeroelasticity: since these algorithms are solved step by step, recursively, an initial value is required. For the nonlinear equations of aeroelasticity, it may be difficult to find an initial value related to a specific periodic solution. Therefore, he proposed that in the time-discrete process of the equations of motion, a set of algebraic equations should be derived by substituting periodic boundary value conditions, so that the nonlinear steady-state periodic solution of the rotor aeroelasticity can be obtained directly.

In this paper, the space-time discretization of the fully intrinsic beam equations is carried out on the basis of a differential quadrature method of high-order precision, and the new calculation formulas for solving the responses of the fully intrinsic dynamic equations of the beam are presented. The validity and applicability of the formulas proposed in this paper are verified by static analysis and modal calculation of the cantilever beam, and the calculation of the beam’s dynamic response, including the transient response and periodic steady-state response.

2. Fully Intrinsic Beam Equations

Figure 1 shows the undeformed and deformed states of the beam. At each point on the axis of the undeformed beam, a reference coordinate system $b(x_1)$ is introduced. At each point along the axis of the deformed beam, a reference coordinate system $B(x_1)$ is introduced. All variables involved in the fully intrinsic equations are expressed in the coordinate system $b$ and $B$. The equations can be expressed in a compact form as

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

(1)

Here, $({)}^{\prime }$ represents the partial derivative with respect to the reference axis of the undeformed beam coordinate system, and ${(}^{·})$ represents the partial derivative with respect to time. $F(x,t)$ and $M(x,t)$ represent the force and moment of the beam section, respectively. $P(x,t)$ and $H(x,t)$ represent the linear momentum and angular momentum of the beam section, respectively. $\gamma (x,t)$ and $\kappa (x,t)$ represent the generalized force and moment strains of the beam, respectively. $V(x,t)$ and $\mathsf{\Omega }(x,t)$ represent the linear and angular velocities of the beam section, respectively. $f(x,t)$ and $m(x,t)$ represent the external force and moment acting on the beam, respectively. $k$ is the initial twist/curvature of the beam that expressed in the undeformed coordinate system $b$. $\widetilde{()}$ is the cross-product operator, and for $k={\left\{\begin{array}{ccc}{k}_1& k_2& k_3\end{array}\right\}}^T$, the operator can be expressed as

$$\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 first two equations in Equations (1) are partial differential equations of linear and angular momentum equilibrium. The latter two are partial differential equations for the kinematic equations, which were derived from the generalized strain-displacement and generalized velocity-displacement equations. This process essentially eliminates the displacement and rotation variables, leaving the equation with only intrinsic quantities.

$$\begin{array}{l}\left\{\begin{array}{c}\gamma \\ \kappa \end{array}\right\}=\left[\begin{array}{cc}R& S\\ S^T& T\end{array}\right]\left\{\begin{array}{c}F\\ M\end{array}\right\}\\ \left\{\begin{array}{c}P\\ H\end{array}\right\}=\left[\begin{array}{cc}\mu \Delta & -\mu \tilde{\xi }\\ \mu \tilde{\xi }& I\end{array}\right]\left\{\begin{array}{c}F\\ M\end{array}\right\}\end{array}$$

(3)

Equations (3) are the equations for the beam constitutive relation and the generalized velocity-momentum relation. $R(x)$, $T(x)$, and $S(x)$ are the flexibility characteristic constants of the beam profile. This linear constitutive relation is valid for small-strain cases. $\mu (x)$, $\xi (x)$, and $I(x)$ represent the beam’s linear density, centroid offset (the centroid of the profile relative to the reference axis), and the moment of inertia per unit length, respectively. Equations (1) and (3) constitute a complete set of differential equations for first-order partial derivatives in time and space.

3. Space-Time Discretization of Fully Intrinsic Equations

3.1. Discretization in the Spatial Domain

To solve the above equations, space discretization must be carried out first. The essence of spatial discretization is to eliminate the derivative terms. The DQ method, which has received much attention in recent years, approximates the derivative term as the weighted linear summation of the function values at all discrete points along the domain. Therefore, the computation of the weight coefficients is very important, which depends on the distribution of discrete points. The approximate solution at the end of time interval obtained by the discrete point distribution used in the [27,28] is equivalent to the approximate solution obtained by using the Pade approximation. Its characteristic is that the accuracy of the approximate solution at the inner points of the domain is no different from it is other points on the distribution, both of which are $p$-order ($p$ is the maximum order of the polynomials). However, the accuracy at the terminal node can be improved to $2p$ where the additional parameter $\mu =1$ (the format is non-dissipative at the high-frequency regime).

This paper firstly adopts this discrete points’ distribution to discretize the spatial domain equations, and this differential discretization method is called the DQ−Pade method in this paper. The $n$th derivative of the function $f(x)$ with respect to the spatial direction $x$ can be expressed as

$$\frac{{\partial }^{(n)}f}{\partial x^{(n)}}=f_x^{(n)}(x_i)={\displaystyle \sum _{k=1}^N{W}_{ik}^{(n)}f(x_k)},\hspace{1em}n=1,\dots ,N$$

(4)

Here, $W_{ik}^{(n)}$ is the weight coefficient, and the algebraic expression is as follows:

$$\begin{array}{l}{W}_{ik}^{(1)}=\frac{M^{(1)}(x_i)}{(x_i-x_k)M^{(1)}(x_k)},\hspace{1em}i,k=0,1,\dots .,N\hspace{1em}but\hspace{1em}i\ne k\\ W_{ii}^{(1)}=\frac{M^{(2)}(x_i)}{2M^{(1)}(x_k)},i=0,1,\dots .,N\end{array}$$

(5)

Among them, $M$, $M^{(1)}$, and $M^{(2)}$ can be calculated by the respective following formulas:

$$\begin{array}{l}M(\zeta )={\zeta }^{N+1}-W_1\zeta -W_2{\zeta }^2-W_3{\zeta }^3-\dots -W_N{\zeta }^N\\ M^{(1)}(\zeta )=(N+1){\zeta }^N-W_1-2W_2\zeta -3W_3{\zeta }^2-\dots -N{W}_N{\zeta }^{N-1}\\ M^{(2)}(\zeta )=N(N+1){\zeta }^{N-1}-2W_2-6W_3\zeta -\dots -N(N-1)W_N{\zeta }^{N-2}\end{array}$$

(6)

where

$$W_k=\frac{{(-1)}^{N-k}N!N!(N+k-2)!}{(k-1)!(k-1)!(N+1-k)!(2N)!}·\frac{2(N+\mu (k-1))}{1+\mu }\hspace{1em},0\le \mu \le 1$$

(7)

The calculation method of discrete points is obtained by the following formula:

$${\zeta }^N-W_1-W_2\zeta -W_3{\zeta }^2-\dots -W_N{\zeta }^{N-1}=0$$

(8)

By taking the root of the above equation, $N$ discrete points can be obtained. Figure 2 shows the distribution at $N=9$. As can be seen in the figure, the discrete points are not evenly distributed, and nodal points (0 or 1) are not included. Therefore, when introducing boundary conditions $V_0$, ${\mathsf{\Omega }}_0$, $F_L$, and $M_L^{}$, Equation (5) should be used to directly calculate the weight coefficient $W_{i0}$ of the start point value or to calculate $W_{iL}^{}$ of end point after some transformations. Thus, the spatial derivative terms in Equation (1) are written as

$$\begin{array}{l}{V}_i^{\prime }={\displaystyle \sum _{k=1}^N{W}_{ik}^{}{V}_k}+W_{i0}{V}_0\\ {\mathsf{\Omega }}_i^{\prime }={\displaystyle \sum _{k=1}^N{W}_{ik}^{}{\mathsf{\Omega }}_k}+W_{i0}{\mathsf{\Omega }}_0\\ F_i^{\prime }={\displaystyle \sum _{k=1}^N{W}_{ik}^{}{F}_k}+W_{iL}{F}_L\\ M_i^{\prime }={\displaystyle \sum _{k=1}^N{W}_{ik}^{}{M}_k}+W_{iL}{M}_L\end{array}$$

(9)

By substituting Equation (9) into Equation (1), the spatial discretization of the intrinsic equations based on DQ−Pade method can be expressed as

$$\begin{array}{l}{\displaystyle \sum _{k=1}^N{W}_{ik}{F}_k+({\tilde{k}}_i+{\tilde{\kappa }}_i)}{F}_i+f_i+W_{iL}{F}_L={\dot{P}}_i+{\tilde{\mathsf{\Omega }}}_i{P}_i\\ {\displaystyle \sum _{k=1}^N{W}_{ik}{M}_k+({\tilde{k}}_i+{\tilde{\kappa }}_i)}{M}_i+({\tilde{e}}_1+{\tilde{\gamma }}_i)F_i+m_i+W_{iL}{M}_L={\dot{H}}_i+{\tilde{\mathsf{\Omega }}}_i{H}_i+{\tilde{V}}_i{P}_i\\ {\displaystyle \sum _{k=1}^N{W}_{ik}{V}_k+({\tilde{k}}_i+{\tilde{\kappa }}_i)}{V}_i+({\tilde{e}}_1+{\tilde{\gamma }}_i){\mathsf{\Omega }}_i+W_{i0}{V}_0={\dot{\gamma }}_i\\ {\displaystyle \sum _{k=1}^N{W}_{ik}{\mathsf{\Omega }}_k+({\tilde{k}}_i+{\tilde{\kappa }}_i)}{\mathsf{\Omega }}_i+W_{i0}{\mathsf{\Omega }}_0={\dot{\kappa }}_i\end{array}$$

(10)

In Equations (9) and (10), the subscripts containing $i$ or $k$ are variables corresponding to discrete points $x_i^{}$ or $x_k$. Arrange the equations with respect to the time term, and then a set of differential equations with respect to time, which are the same in form as other spatial discretization methods, is obtained:

$$A_{ji}{\dot{q}}_i+B_{ji}{q}_i+C_{jik}{q}_i{q}_k+D_j=0$$

(11)

where $q$ is the vector of unknown intrinsic variables, $A$ and $B$ are linear coefficient matrices. $C$ is a nonlinear coefficient matrix. $D$ is the external load vector, which also contains intrinsic dynamic boundary conditions $V_0$, ${\mathsf{\Omega }}_0$, $F_L$, and $M_L^{}$.

It is worth mentioning that the DQ−Pade method can also be used to discretize the equations of displacements and rotations recovering from the intrinsic variables (Hodges [5]):

$$\begin{array}{l}{u}^{\prime }=C^T(e_1+\gamma )-e_1-\tilde{k}u\\ {\mathsf{\theta }}^{\prime }=(\Delta +\frac{1}{2}\tilde{\theta }+\frac{1}{4}\theta {\theta }^T)(\kappa +k-Ck)\end{array}$$

(12)

where

$$C=\frac{[1-(1/4){\theta }^T\theta ]\Delta -\tilde{\theta }+(1/2)\theta {\theta }^T}{1+(1/4){\theta }^T\theta }$$

(13)

The equations after using DQ−Pade discretization can be expressed as

$$\begin{array}{l}{\displaystyle \sum _{k=1}^N{W}_{ik}{u}_k}+W_{i0}{u}_0-C_i^T(e_1+{\gamma }_i)+e_1+{\tilde{k}}_i{u}_i=0\\ {\displaystyle \sum _{k=1}^N{W}_{ik}{\theta }_k}+W_{i0}{\theta }_0-(\Delta +\frac{1}{2}{\tilde{\theta }}_i+\frac{1}{4}{\theta }_i{\theta }_i^T)({\kappa }_i+k_i-C_i{k}_i)=0\end{array}$$

(14)

Here, $\theta$ is the Rodriguez parameter representing the rotation matrix. $u_0$ and ${\theta }_0$ are the geometric boundary conditions. The discretized equations are algebraic equations containing nonlinearity and need to be solved iteratively.

3.2. Discretization in the Time Domain

DQ−Pade was originally derived for the solution of time-domain responses and has been successfully applied in time-domain problems [20,21]. Therefore, it is naturally used in this paper to discretize the first-order partial differential equation with respect to time obtained after space discretization. Compared with the traditional recursive method, this paper presents the global format of the DQ−Pade method, which can get the responses of all discrete time points at once.

If there are M + 1 time points ($t_0^{},t_1,\dots ,t_M$) in a time domain $T$, $t_i(i=1,\dots ,M)$ is the best approximation of the time discrete point calculated by Equations (7) and (8). Additionally, $t_0^{}$ is the start point of the time domain $T$. Thus, the time derivative term $\dot{q}$ in the discrete Equation (11) can be discretized by the DQ−Pade method:

$$\dot{q}(t)={\displaystyle \sum _{k=1}^M{W}_{ik}^{}{q}_k}+W_{i0}{q}_0^{initial}$$

(15)

where the weight coefficients $W$ can be calculated using Equation (5)–(8); $q_0^{initial}$ is the initial condition. Substitute Equation (15) into (11). Then, the following equation is obtained:

$$A{\displaystyle \sum _{k=1}^M{W}_{ik}^{}{q}_k}+A{W}_{i0}{q}_0^{initial}+B{q}_i+C(q_i)q_i+D_i=0,\hspace{1em}i=1,\dots ,M$$

(16)

Equation (16) can also be written in matrix form as

$$\begin{array}{l}\left[\begin{array}{ccccc}A{W}_{11}+B+C(q_1)& A{W}_{12}& A{W}_{13}& \dots & A{W}_{1M}\\ A{W}_{21}& A{W}_{22}+B+C(q_2)& A{W}_{23}& \dots & A{W}_{2M}\\ A{W}_{31}& A{W}_{32}& A{W}_{33}+B+C(q_3)& \dots & A{W}_{3M}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ A{W}_{M1}& A{W}_{M2}& A{W}_{M3}& \dots & A{W}_{MM}+B+C(q_M)\end{array}\right]\left\{\begin{array}{c}{q}_1\\ q_2\\ q_3\\ ⋮\\ q_M\end{array}\right\}\\ +\left\{\begin{array}{c}A{W}_{10}\cdot q_0^{initial}\\ A{W}_{20}\cdot q_0^{initial}\\ A{W}_{30}\cdot q_0^{initial}\\ ⋮\\ A{W}_{M0}\cdot q_0^{initial}\end{array}\right\}+\left\{\begin{array}{c}{D}_1\\ D_2\\ D_3\\ ⋮\\ D_M\end{array}\right\}=0\end{array}$$

(17)

The above equation can be written in a simplified form:

$$A_{total}\cdot q_{total}+A_0\cdot q_0^{initial}+D_{total}=0$$

(18)

This is the initial value formula for the time-discrete solution based on the DQ−Pade method. Given the initial condition $q_0^{initial}$, the responses at all discrete time points can be solved at one time.

If the time domain $T$ is a period, the steady-state periodic response of the intrinsic equations is required. For this case, it is difficult to give the initial value $q_0^{initial}$ satisfying the final steady-state periodic response when Equation (18) is still used. Therefore, it is necessary to deduce the discrete formula with consideration of the periodic boundary condition

$$q_0^{}=q_T$$

(19)

Here, $q_T^{}$ can be obtained by polynomial interpolation:

$$q_T=L_0(1)q_0+L_1(1)q_1+\dots +L_M(1)q_M$$

(20)

where $L_k$ is the Lagrange interpolation polynomial. Substitute Equation (19) into (18) to get

$$q_0^{boundary}=l_1{q}_1+\dots +l_M{q}_M,l_i=L_i/(1-L_0)$$

(21)

$q_0^{boundary}$ is the constraint condition of periodic boundary value problems. Bring it into Equation (16) and write it in matrix form:

$$\begin{array}{l}\left[\begin{array}{ccccc}A{W}_{11}+B+C(q_1)& A{W}_{12}& A{W}_{13}& \dots & A{W}_{1M}\\ A{W}_{21}& A{W}_{22}+B+C(q_2)& A{W}_{23}& \dots & A{W}_{2M}\\ A{W}_{31}& A{W}_{32}& A{W}_{33}+B+C(q_3)& \dots & A{W}_{3M}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ A{W}_{M1}& A{W}_{M2}& A{W}_{M3}& \dots & A{W}_{MM}+B+C(q_M)\end{array}\right]\left\{\begin{array}{c}{q}_1\\ q_2\\ q_3\\ ⋮\\ q_M\end{array}\right\}\\ +\left[\begin{array}{ccccc}A{W}_{10}\cdot l_1& A{W}_{10}\cdot l_2& A{W}_{10}\cdot l_3& \dots & A{W}_{10}\cdot l_M\\ A{W}_{20}\cdot l_1& A{W}_{20}\cdot l_2& A{W}_{20}\cdot l_3& \dots & A{W}_{20}\cdot l_M\\ A{W}_{30}\cdot l_1& A{W}_{30}\cdot l_2& A{W}_{30}\cdot l_3& \dots & A{W}_{30}\cdot l_M\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ A{W}_{M0}\cdot l_1& A{W}_{M0}\cdot l_2& A{W}_{M0}\cdot l_3& \dots & A{W}_{M0}\cdot l_M\end{array}\right]\left\{\begin{array}{c}{q}_1\\ q_2\\ q_3\\ ⋮\\ q_M\end{array}\right\}+\left\{\begin{array}{c}{D}_1\\ D_2\\ D_3\\ ⋮\\ D_M\end{array}\right\}=0\end{array}$$

(22)

Similarly, the equation can be simplified as

$$A_{total}\cdot q_{total}+A_0\cdot q_0^{boundary}+D_{total}=0$$

(23)

Equation (23) is the periodic boundary value formula based on the DQ−Pade method for a time discrete solution. Like the initial value formula of Equation (17) and (18), the derived algebraic equations contain nonlinear terms, $C$. Therefore, it is necessary to solve the algebraic equations via nonlinear iteration, e.g., a Newton iteration method. To get the function’s values at other points along the domain, Lagrange polynomial interpolation can be used.

3.3. DQ−Pade Element Method

It can be seen in the space discrete equation, Equation (10), and the time discrete equation, Equations (18) and (23), that the coefficient matrixes of the equations are fully matrixes (each element of a matrix has a nonzero value). When dealing with large-scale problems, such as rapid changes or discontinuities in spatial or time domains, the quantity of discrete points needs to be increased, and then the cost of storing and operating the coefficient matrix will increase rapidly, reducing the computational efficiency. Another reason that can be seen in Equation (8) is that when the discrete number $N$ is large, it is difficult to find the root of the formula, which leads to a decrease in the calculation accuracy of the weight coefficient. Meanwhile, there are more nodal points in the domain, which is beneficial to the accuracy of the DQ−Pade method to some extent. Therefore, it is necessary to divide the entire time/space domain into multiple sub-domains, similarly to the finite element method. Additionally, the sub-domains are connected through the weight coefficients of the public points.

If $M$ elements are divided in the spatial domain, the spatial discrete Equation (10) is transformed into

$$\begin{array}{l}{\displaystyle \sum _{k=1}^N{W}_{ik}{F}_k^m+({\tilde{k}}_i^m+{\tilde{\kappa }}_i^m)}{F}_i^m+f_i^m+W_{iL}({\displaystyle \sum _{k=1}^N{L}_k(0)}{F}_k^{m+1})={\dot{P}}_i^m+{\tilde{\mathsf{\Omega }}}_i^m{P}_i^m\\ {\displaystyle \sum _{k=1}^N{W}_{ik}{M}_k^m+({\tilde{k}}_i^m+{\tilde{\kappa }}_i^m)}{M}_i^m+({\tilde{e}}_1+{\tilde{\gamma }}_i^m)F_i^m+m_i^m+W_{iL}({\displaystyle \sum _{k=1}^N{L}_k(0)}{M}_k^{m+1})={\dot{H}}_i^m+{\tilde{\mathsf{\Omega }}}_i^m{H}_i^m+{\tilde{V}}_i^m{P}_i^m\\ {\displaystyle \sum _{k=1}^N{W}_{ik}{V}_k^m+({\tilde{k}}_i^m+{\tilde{\kappa }}_i^m)}{V}_i^m+({\tilde{e}}_1+{\tilde{\gamma }}_i^m){\mathsf{\Omega }}_i^m+W_{i0}({\displaystyle \sum _{k=0}^N{L}_k(1)}{V}_k^{m-1})={\dot{\gamma }}_i^m\\ {\displaystyle \sum _{k=1}^N{W}_{ik}{\mathsf{\Omega }}_k^m+({\tilde{k}}_i^m+{\tilde{\kappa }}_i^m)}{\mathsf{\Omega }}_i^m+W_{i0}({\displaystyle \sum _{k=1}^N{L}_k(1)}{\mathsf{\Omega }}_k^{m-1})={\dot{\kappa }}_i^m\end{array}$$

(24)

The boundary conditions can be expressed as

$$\left\{\begin{array}{c}{\displaystyle \sum _{k=0}^N{L}_k(1)}{V}_k^{m-1}=V_0^{},{\displaystyle \sum _{k=1}^N{L}_k(1)}{\mathsf{\Omega }}_k^{m-1}={\mathsf{\Omega }}_0^{},\mathrm{when}m=1\\ {\displaystyle \sum _{k=1}^N{L}_k(0)}{F}_k^{m+1}=F_L,{\displaystyle \sum _{k=1}^N{L}_k(0)}{M}_k^{m+1}=M_L^{},\mathrm{when}m=M\end{array}\right.$$

(25)

When $M$ elements are divided in the time domain, the initial value Equation (17) is transformed into

$$\begin{array}{l}\left[\begin{array}{cccccccccc}{\mathrm{AW}}_{11}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_1^1)& \cdots & {\mathrm{AW}}_{1N}& & & & & & & \\ ⋮& \ddots & ⋮& & & & & & & \\ {\mathrm{AW}}_{N1}& \cdots & {\mathrm{AW}}_{NN}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_N^1)& & & & & & & \\ {\mathrm{AW}}_{10}\cdot L_1(1)& \cdots & {\mathrm{AW}}_{10}\cdot L_N(1)& {\mathrm{AW}}_{11}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_1^2)& \cdots & {\mathrm{AW}}_{1N}& & & & \\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮& & & & \\ {\mathrm{AW}}_{N0}\cdot L_1(1)& \cdots & {\mathrm{AW}}_{N0}\cdot L_N(1)& {\mathrm{AW}}_{N1}& \cdots & {\mathrm{AW}}_{NN}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_N^2)& & & & \\ & & & & & & \ddots & & & \\ & & & \ddots & {\mathrm{AW}}_{10}\cdot L_1(1)& \cdots & {\mathrm{AW}}_{10}\cdot L_N(1)& {\mathrm{AW}}_{11}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_1^M)& \cdots & {\mathrm{AW}}_{1N}\\ & & & & ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ & & & & {\mathrm{AW}}_{N0}\cdot L_1(1)& \cdots & {\mathrm{AW}}_{N0}\cdot L_N(1)& {\mathrm{AW}}_{N1}& \cdots & {\mathrm{AW}}_{NN}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_N^M)\end{array}\right]\\ \cdot \left\{\begin{array}{c}{\mathrm{q}}_1^1\\ ⋮\\ {\mathrm{q}}_N^1\\ {\mathrm{q}}_1^2\\ ⋮\\ {\mathrm{q}}_N^2\\ ⋮\\ {\mathrm{q}}_1^M\\ ⋮\\ {\mathrm{q}}_N^M\end{array}\right\}+\left\{\begin{array}{c}{\mathrm{D}}_1^1+{\mathrm{AW}}_{10}\cdot q_0^{initial}\\ ⋮\\ {\mathrm{D}}_N^1+{\mathrm{AW}}_{N0}\cdot q_0^{initial}\\ {\mathrm{D}}_1^2\\ ⋮\\ {\mathrm{D}}_N^2\\ ⋮\\ {\mathrm{D}}_1^M\\ ⋮\\ {\mathrm{D}}_N^M\end{array}\right\}=0\end{array}$$

(26)

and the boundary value formula, Equation (22), is transformed into

$$\begin{array}{l}\left[\begin{array}{cccccccccc}{\mathrm{AW}}_{11}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_1^1)& \cdots & {\mathrm{AW}}_{1N}& & & & & {\mathrm{AW}}_{10}\cdot L_1(1)& \cdots & {\mathrm{AW}}_{10}\cdot L_N(1)\\ ⋮& \ddots & ⋮& & & & & ⋮& \ddots & ⋮\\ {\mathrm{AW}}_{N1}& \cdots & {\mathrm{AW}}_{NN}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_N^1)& & & & & {\mathrm{AW}}_{N0}\cdot L_1(1)& \cdots & {\mathrm{AW}}_{N0}\cdot L_N(1)\\ {\mathrm{AW}}_{10}\cdot L_1(1)& \cdots & {\mathrm{AW}}_{10}\cdot L_N(1)& {\mathrm{AW}}_{11}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_1^2)& \cdots & {\mathrm{AW}}_{1N}& & & & \\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮& & & & \\ {\mathrm{AW}}_{N0}\cdot L_1(1)& \cdots & {\mathrm{AW}}_{N0}\cdot L_N(1)& {\mathrm{AW}}_{N1}& \cdots & {\mathrm{AW}}_{NN}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_N^2)& & & & \\ & & & & & & \ddots & & & \\ & & & \ddots & {\mathrm{AW}}_{10}\cdot L_1(1)& \cdots & {\mathrm{AW}}_{10}\cdot L_N(1)& {\mathrm{AW}}_{11}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_1^M)& \cdots & {\mathrm{AW}}_{1N}\\ & & & & ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ & & & & {\mathrm{AW}}_{N0}\cdot L_1(1)& \cdots & {\mathrm{AW}}_{N0}\cdot L_N(1)& {\mathrm{AW}}_{N1}& \cdots & {\mathrm{AW}}_{NN}+\mathrm{B}+\mathrm{C}({\mathrm{q}}_N^M)\end{array}\right]\\ \cdot \left\{\begin{array}{c}{\mathrm{q}}_1^1\\ ⋮\\ {\mathrm{q}}_N^1\\ {\mathrm{q}}_1^2\\ ⋮\\ {\mathrm{q}}_N^2\\ ⋮\\ {\mathrm{q}}_1^M\\ ⋮\\ {\mathrm{q}}_N^M\end{array}\right\}+\left\{\begin{array}{c}{\mathrm{D}}_1^1\\ ⋮\\ {\mathrm{D}}_N^1\\ {\mathrm{D}}_1^2\\ ⋮\\ {\mathrm{D}}_N^2\\ ⋮\\ {\mathrm{D}}_1^M\\ ⋮\\ {\mathrm{D}}_N^M\end{array}\right\}=0\end{array}$$

(27)

The main diagonal matrix block in Equation (26) and (27) is the coefficient matrix of each element, and the sub-diagonal matrix is the connection coefficient matrix formed by each element through the public point (the connecting point of two elements). The matrix block in the upper right corner of the Equation (27) is a connected coefficient matrix formed by the periodic boundary condition. The initial value condition is reflected in the first sub-vector in the total load vector of Equation (26).

It can be seen in the matrix form that the coefficient matrix of the equations after dividing is characterized by a sparse banded distribution. Compared with the method of forming full matrix coefficients without dividing elements, the DQ−Pade element method can improve the calculation speed of the matrix and reduce the storage cost, thereby improving the calculation efficiency. The accuracy and efficiency of the method are verified by the following examples.

4. Numerical Results

To verify the applicability and validity of the proposed time/spatial discretization formula for intrinsic beam equations, two examples are considered. Firstly, the spatial discrete formula is verified by static analysis and modal calculation of the cantilever beam. Another example is calculating the dynamic response of a rotating cantilever beam to verify the accuracy and efficiency of the proposed time-discrete formula.

4.1. Static Analysis and Modal Calculation of the Cantilever Beam

Firstly, the static analysis of the following force applied to the tip of the cantilever beam is considered. Then, the velocity condition is applied to its root to calculate its natural frequency [10]. The structural characteristics of the cantilever beam are shown in

Table 1. The structural characteristics of the cantilever beam.

Parameter	Value
Span	16 m
Chord	1 m
Mass density	0.75 kg/m
Moment of inertia (50% chord)	0.1 kg m
Spanwise elastic axis	50% chord
Center of gravity	50% chord
Bending rigidity	2 × 104 Nm2
Torsional rigidity	1 × 104 Nm2
Bending rigidity (chordwise)	4 × 106 Nm2
Shear/extensional rigidity	∞

. For static analysis, apply a following force $F_{tip}=3EI/L^2$ at the tip. Simitses and Hodges [34] gave the exact solution of the dimensionless equation (made by $F_{tip}L$) for the bending moment at the root ${\overline{M}}_{root}=0.8104403623$. Due to static analysis, time-dependent terms are eliminated in Equation (9), and it can be simplified to

$$\begin{array}{l}{\displaystyle \sum _{k=1}^N{W}_{ik}{F}_k+{\tilde{\kappa }}_i}{F}_i+W_{iL}{F}_L=0\\ {\displaystyle \sum _{k=1}^N{W}_{ik}{M}_k+{\tilde{\kappa }}_i}{M}_i+{\tilde{e}}_1{F}_i+W_{iL}{M}_L=0\end{array}$$

Table 2. Results of the root bending moment calculated by different DQ methods.

Method	n = 3	n = 4	n = 5	n = 6	n = 7
CDM	0.795943708	0.804111347	0.806942329	0.808227074	0.808914563
GDQ	0.628913129	0.848504651	0.803100361	0.809702945	0.810900600
WDQ	-	-	-	-	0.811791100
Present	0.806762171	0.809924405	0.810414363	0.810437905	0.810440225
	n = 8	n = 9	n = 10	n = 13	n = 21
	0.809324779	0.809589080	0.809769336	0.810064442	0.810305505
	0.810379136	0.810430350	0.810445453	0.810440413	0.810440363
	-	0.810440438	-	0.810440087	0.810440342
	0.810440354	0.810440362	0.810440363	0.810440356	0.812366715

lists the results of the root bending moment calculated by different DQ methods. The GDQ method is the discretization method based on Lagrange polynomial space that was used to calculate weight coefficients in [29]. The wavelet-differential quadrature (WDQ) method is a DQ method based on the adaptive multi-layer wavelet collocation method in [27,28]. The number of discrete points is determined by the wavelet resolution layer $j$ and the translation factor $N$ ($n=2^{j+1}+2N+1$), and $N=2$, $j=0\sim 4$ is taken here. The results of the central difference method (CDM) used by the Hodges [4] are also listed in this table.

From Table 2 and Figure 3 and Figure 4, it can be clearly seen that the DQ methods has better accuracy than the CDM. Compared with the other two DQ methods, the DQ−Pade method has obvious advantages in convergence speed. With the same computational cost, the accuracy of this method is at least three orders of magnitude higher than those of the other two methods. To achieve the minimum error magnitude of 10⁻¹⁰, only nine discrete points are needed in this method. The computational cost is much lower than 16 points for the GDQ method and 37 points for the WDQ method.

However, it can be seen in Figure 3 that when the number of discrete points exceeds 9, the calculation accuracy of this method decreases as the number of discrete points increases. The reason is that with the increase in the number of discrete points, it will be difficult to calculate the distributed discrete points, resulting in a decrease in the accuracy of calculating discrete weight coefficients. Figure 4 shows the error curve of calculating the bending moment at the root by dividing different discrete element numbers. It can be seen in the figure that dividing discrete elements can effectively improve the calculation accuracy and restrain the divergence of error. At the same time, it is worth noting that as the number of element (M) increases, the rate of error convergence decreases. Additionally, it is observed that when the number of discrete points in single element exceeds 10, the accuracy usually decreases. Therefore, it is recommended to divide the sub-domains as little as possible in the entire domain, and discretize points in any sub-domain no more than 10 times.

If the intrinsic dynamic boundary conditions of the rotating cantilever beam [10] are

$$V_0^T=[0,51.03,0],{\mathsf{\Omega }}_0^T=[0,0,3.189],F_L^T=[0,0,0],M_L^T=[0,0,0]$$

Wright [35] gave an exact solution for the bending frequency in this case:

$${\omega }_{1st,bending}=5.703,{\omega }_{2nd,bending}=18.72,{\omega }_{3rd,bending}=44.50$$

Table 3. Results of a beam rotating at natural frequency.

Frequency	n = 5	n = 6	n = 7	n = 8	n = 9	n = 10
1st bending/Present	5.7024	5.7024	5.7024	5.7024	5.7024	5.7024
1st bending/GDQM	5.6983	5.7027	5.7023	5.7024	5.7024	5.7024
1st bending/CDM	6.5402	6.3841	6.2779	6.2005	6.1417	6.0953
2nd bending/Present	18.7226	18.7229	18.7229	18.7229	18.7229	18.7229
2nd bending/GDQM	18.7054	18.7208	18.7226	18.7228	18.7229	18.7229
2nd bending/ CDM	19.7283	19.7155	19.6650	19.6021	19.5377	19.4790
3rd bending/Present	44.7008	44.5122	44.4997	44.4985	44.4985	44.4985
3rd bending/GDQM	42.3069	42.2699	44.4587	44.5012	44.4983	44.4984
3rd bending/CDM	35.3993	39.3702	41.2891	42.3684	43.0291	43.4587

lists the frequency results of CDM, DQ−Pade, and the GDQ method. The CDM is sensitive to the calculation step, so it needs more discrete points. It was found that the frequency results had considerable accuracy when the number of discrete points exceeds 150. Compared with the GDQ method, our method still has the prominent characteristics of low cost and high precision, so it is an ideal choice for the spatial discretization of the fully intrinsic equations.

4.2. Dynamic Response of the Rotating Beam

In the second example, the dynamic response analysis of rotating cantilever beam is considered, including transient and steady-state responses. The length of the beam is 1 m, and the rotational speed is 70 rad/s. The material properties of the beam are shown in

Table 4. Material properties of the beam.

Parameter	Value
Mass per unit length	0.2 kg/m
$\mathrm{Moment}\mathrm{of}\mathrm{inertia}\mathrm{per}\mathrm{unit}\mathrm{length}{I}_{xx}$	10−4 kg·m
$\mathrm{Moment}\mathrm{of}\mathrm{inertia}\mathrm{per}\mathrm{unit}\mathrm{length}{I}_{yy}$	10−6 kg·m
$\mathrm{Moment}\mathrm{of}\mathrm{inertia}\mathrm{per}\mathrm{unit}\mathrm{length}{I}_{zz}$	10−4 kg·m
$\mathrm{Extensional}\mathrm{rigidity}{K}_{11}$	106 N
$\mathrm{Shear}\mathrm{rigidity}{K}_{22}$	1020 N
$\mathrm{Shear}\mathrm{rigidity}{K}_{33}$	1020 N
$\mathrm{Torsional}\mathrm{rigidity}{K}_{44}$	50 N·m2
$\mathrm{Bending}\mathrm{rigidity}{K}_{55}$	50 N·m2
$\mathrm{Bending}\mathrm{rigidity}\left(\mathrm{chordwise}\right)K_{66}$	1000 N·m2

[32].

A periodic excitation force $F=50\sin (20t)N$ is applied in the flapping direction at the tip of the beam. The transient response results of DYMORE with a calculation step size of 0.001 s is given in [32]. This present method had the same calculation cost: 125 time elements and 8 time points in each element, denoted as M125−N8. Tip displacement, tip rotation angle, root force, and bending moment are shown in Figure 5, Figure 6, Figure 7 and Figure 8, respectively. From these figures, it can be seen that the present method is in good agreement with the calculation results of DYMORE. Even for axial forces with insufficient calculation accuracy in [32], this method can still achieve satisfactory accuracy.

With the same rotational speed, a periodic harmonic force is applied in the direction of tip flapping:

$$F=10\sin (70t)+5\sin (2\cdot 70t)+2\sin (3\cdot 70t)N$$

The authors of [32] used the backward second-order Euler method to calculate the time domain response. We took the Euler method as a reference for the steady-state response. Figure 9, Figure 10, Figure 11 and Figure 12 shows the steady-state periodic response curves of tip displacements, rotations, root forces, and moments calculated by the two methods. In the figures, the steady-state period formula of the DQ−Pade element method uses a combination of eight time elements divided in one period and nine distribution points used in each element (M8−N9), both of which use the same time element step size, $t=2\pi /70/72$ (one time element per 5°). It can be seen in the figures that the two curves basically coincide, and both can well meet the periodic condition. When the calculation step size is 5°, the Euler method required six periods, and it took 114 s, to reach the so-so periodic boundary condition (absolute error < 1). When calculating 15 cycles, which took 286 s, the periodic boundary value condition converged to the order of 10⁻². The results of Euler method in the figures are given for this condition.

The steady-state periodic boundary value formula of the DQ−Pade method itself satisfies the periodic boundary conditions. At the same calculation cost, i.e., calculating 72 time discrete points in a cycle, the calculation accuracy and calculation time were compared, with different combinations of discrete element numbers and discrete point numbers (M × N = 72). Taking M36−N10 (one time element per 1°) as a reference, and taking the bending moment of the beam root as an example, Figure 13 shows the convergence of different combinations. It can be seen that the combination of fewer discrete elements and more discrete points distributed in each element can obtain the best accuracy.

Table 5. The calculation times of different combinations of discrete elements and points.

M/N	36/2	24/3	18/4	12/6	9/8	8/9
CPU Time	23.01 s	25.17 s	28.84 s	38.88 s	53.66 s	63.12 s

gives the calculation times of different combinations. The combination with the best accuracy took the longest, but this was still far less than the calculation time of the Euler method. The reason is that each time discrete point has a state quantity that contains all the structural degrees of freedom, and the increase in discrete points in the time element leads to the multiplication of the matrix operation cost. Therefore, it is needed to make a balance between accuracy and computational efficiency.

5. Conclusions

In this paper, a new method based on the DQ−Pade method was proposed to solve the nonlinear fully intrinsic equations of geometrically exact beams. Starting from the spatial and temporal discretization of the fully intrinsic equations, the DQ−Pade method was used to re-deduce the spatial discretization equations of the fully intrinsic beam equations, and the static analysis and modal calculation of the beam were investigated. It was found that the accuracy and efficiency of this method are better than those of the traditional GDQ method. In terms of time discretization, a calculation formula of initial value and steady-state periodic boundary value is proposed, based on the DQ−Pade method. With the consideration of the efficiency and accuracy of the DQ−Pade method, the DQ−Pade element method was further proposed. The results show that the method is very effective, and it has outstanding performance in terms of accuracy and efficiency, which makes a good foundation for the subsequent rotor vibration load prediction and aeroelastic stability analysis.