Differential Quadrature Method for Fully Intrinsic Equations of Geometrically Exact Beams

: 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 exa mples 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.


Introduction
Due to the high aspect ratio of a rotor blade, which is usually regarded as a beamlike 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 stressbased 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 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.  Here, ′ ( ) 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. γ ( , )

Fully Intrinsic Beam Equations
x t and κ( , ) x t represent the generalized force and moment strains of the beam, respectively.
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. µ( ) x , ξ( ) 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.

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 µ = 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 Here, is the weight coefficient, and the algebraic expression is as follows: M , and (2) M can be calculated by the respective following formulas: The calculation method of discrete points is obtained by the following formula: By taking the root of the above equation, N discrete points can be obtained. Figure  2 shows the distribution at = 9 N . 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 introduc-  By substituting Equation (9) into Equation (1), the spatial discretization of the intrinsic equations based on DQ−Pade method can be expressed as In Equations (9) and (10), the subscripts containing i or k are variables corresponding to discrete points i x or k x . 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: 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 0 V , Ω 0 , L F , and L M .
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]): k Ck (12) where θ θ θ θθ θ θ The equations after using DQ−Pade discretization can be expressed as γ θ θ θ θ θ κ = = + − + + + = Here, θ is the Rodriguez parameter representing the rotation matrix. 0 u and θ 0 are the geometric boundary conditions. The discretized equations are algebraic equations containing nonlinearity and need to be solved iteratively.

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.
best approximation of the time discrete point calculated by Equations (7) and (8). Additionally, 0 t is the start point of the time domain T . Thus, the time derivative term  q in the discrete Equation (11) can be discretized by the DQ−Pade method: where the weight coefficients W can be calculated using Equation (5)-(8); 0 initial q is the initial condition. Substitute Equation (15) into (11). Then, the following equation is obtained: Equation (16) can also be written in matrix form as 11  1  12  13  1  1   21  22  2  23  2  2   31  32  33  3  3  3   1  2  3 10 0 20 0 The above equation can be written in a simplified form: This is the initial value formula for the time-discrete solution based on the DQ−Pade method. Given the initial condition 0 initial q , 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 0 initial q 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 = 0 T q q (19) Here, T q can be obtained by polynomial interpolation: where k L is the Lagrange interpolation polynomial. Substitute Equation (19) into (18) Similarly, the equation can be simplified as 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.

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 The boundary conditions can be expressed as , when 1 When M elements are divided in the time domain, the initial value Equation (17) is   (26) and the boundary value formula, Equation (22), is transformed into  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.

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.

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 Ta Table 2 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 ( is taken here. The results of the central difference method (CDM) used by the Hodges [4] are also listed in this table.   Table 2 and Figures 3 and 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 −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.   Table 3 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.

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 [32]. A periodic excitation force = 50 sin(20 ) F t 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 Figures 5−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.

tN
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. Figures 9−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, π = 2 / 70 / 72 t (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 −2 . 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 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.  Figure 13. Beam root flapping moment curves with different combinations of discrete elements and points.

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 steadystate 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.