Nonlinear Dynamic Analysis Framework for Slender Structures Using the Modal Rotation Method

High-aspect-ratio wings have been extensively researched in recent years owing to their vast applications in aircraft. In particular, these wings are known for their ability to suppress the induced drag that occurs when creating lift, which can be useful for high-altitude long-endurance (HALE) aircraft [1]. HALE aircraft need to handle the thin air of the atmosphere. This warrants wings with a lift-to-drag ratio larger than that of the wings adopted in conventional passenger aircraft that fly in the troposphere. High-aspect-ratio wings can be advantageous not only for HALE but also for passenger aircraft [2] as they can ensure high fuel efficiency.

Although high-aspect-ratio wings can improve future aircraft performance, their practical implementation can be challenging. The analyses of high-aspect-ratio wings require nonlinear structural models because the wings experience large deformations, which cannot be described by the linear superposition of displacement modes. Therefore, nonlinear finite element methods are widely used for analyzing large deformations. After obtaining a detailed design for the strength of the wing via static analysis performed using a high-fidelity nonlinear shell/plate or solid element model, dynamic analysis, such as flight simulation, control system design, and gust response analysis, is required. However, the high computational cost of the nonlinear shell/plate and solid model serve as a disadvantage when investigating global motion based on dynamic analysis. Therefore, computationally efficient dynamic analysis frameworks based on low-fidelity nonlinear beam models have been developed for dynamic analyses [3].

Palacios et al. [3] reviewed three nonlinear beam formulations that are widely used for the dynamic analysis of the high-aspect-ratio wings. The first is a displacement-based geometrically exact beam formulation (GEBF) [4] that uses displacements and rotations as variables; the second is an intrinsic beam formulation [5] that uses the internal forces and moments as variables; and the third is a strain-based beam formulation [6] that uses strains and curvatures as variables. Classifications and detailed explanations of these three methods are presented in Refs. [3,7], and modal approaches have been introduced to the nonlinear beam formulations [8–10].

The multibody dynamics approach has also been used for dynamic analysis of high-aspect-ratio wings [11]. As reported in Ref. [12], the floating frame of reference formulation (FFRF) used in multibody dynamics [13–15] is applicable only to small deformations. However, there is a way to express the large deformation of the high-aspect-ratio wings based on FFRF. First, a high-aspect-ratio wing is divided into substructures with small deformation. Second, these substructures are modeled by FFRF. Lastly, these substructures are elastically jointed considering their relative rotation.

Geometrically exact beam formulation is also used in multibody dynamics. GEBF does not require substructure division and is the primary option for analyzing high-aspect-ratio wings [16]. GEBF [17] was implemented in open-source software, namely, the simulation of high-aspect-ratio airplanes in Python (SHARPy) [18,19]. GEBF was practically applied to analyses of aircraft with high-aspect-ratio wings [20]. The development of GEBF has continued over the years, for instance, Fan et al. [21] proposed a GEBF that interpolates only the rotation parameter, whereas the basic GEBF interpolates both the position and rotation independently. Because of this rotation-only interpolation, the beam comprised fewer degrees-of-freedom per finite element node.

The absolute nodal coordinate formulation (ANCF) developed by the multibody dynamics community [22] is another method used for the analysis of high-aspect-ratio wings [23]. ANCF interpolates only the position vector. Because the rotation of the element is described by the position vector gradient, ANCF comprises a constant mass matrix and no singularity [24,25]. ANCF beam elements with various position vector interpolations have been developed [26–28]. The dynamic analysis results of the ANCF beam element with standard polynomial interpolation, namely, the fully parameterized element, are in good agreement with the wind tunnel experimental results [29]. Model reduction methods for ANCF beam elements were developed based on modal transformation [30,31] and vector-strain transformation [32–34] to increase the computational speed. Dynamic analyses using the reduced-order model were validated via wind tunnel experiments [35]. Furthermore, an efficient parallel computation technique has been investigated for ANCF [36]. Owing to these advantages, developments, and experimental validations, the ANCF beam element has become a reliable tool for the dynamic analysis of high-aspect-ratio wings.

Corotational formulations are also employed for wings exhibiting large deformations [37,38]. In corotational formulations, a large deformation is expressed by the large translation and rotation of each substructure or element under the small-strain assumption.

The dynamic analysis of high-aspect-ratio wings can be performed at a low computational cost using the aforementioned nonlinear beam models. However, a gap exists between the static analysis model used for the detailed design of wing strength (i.e., a high-fidelity shell/plate or solid model) and the dynamic analysis model used for flight simulation, controller design, and gust response analysis (i.e., a low-fidelity beam model). If the modal information obtained from the high-fidelity model is used for the dynamic analysis, the computational cost can be reduced and the cross-sectional beam parameter calculation can be avoided. However, the linear superposition of displacement modes cannot be used to describe the geometrically nonlinear deformation. Moreover, a nonlinear normal mode approach [39] can be used only for a relatively simple structure [40]. Although the substructure modal approaches [41,42] can describe large deformations, they require additional tasks, such as the substructure division, boundary condition selection to calculate the substructure displacement modes, and substructure connection using an appropriate technique (e.g., fictitious mass technique [43]). To address these challenges, Drachinsky and Raveh [44] proposed a sophisticated modal approach, namely, the modal rotation method (MRM), which does not require the aforementioned tasks.

Modal rotation method describes a large deformation based on the linear superposition of curvature modes. The curvature modes are calculated using the second-order spatial derivatives of the displacement modes. Figure 1 depicts a conceptual diagram describing large deformations in MRM. A reference line is set in the longitudinal direction of the structure. The reference line serves as a reference along which the structure is divided into segments. The segments serve to represent large deformations based on their curvature. A fine division is required to ensure the change in deflection angle on each segment is small. The MRM does not require the beam cross-sectional parameters because the stiffness matrix of the linear/linearized model can be used directly, and the geometrical nonlinear stiffness is considered by correcting the internal bending moment in the iterative calculation. Furthermore, the computational cost of the MRM is low because the analysis is performed using only linear modes.

The first step in MRM is to obtain the displacement modes from a linear/linearized finite element model. The second step involves calculating the curvature modes. The reference line is selected for the structure in the third step. The reference line is then divided into several parts in the fourth step. These parts are referred to as “segments” rather than “substructures.” The fifth step involves rotating the segments recursively from the beam root to the tip based on the relative rotation angles calculated via the linear superposition of the curvature modes. The entire structure undergoes a large deformation, whereas the strain is small in certain cases. MRM assumes one such case and utilizes the fact that the curvature proportional to the strain is small even for a large deformation. Throughout this process, modal analysis is required only once for the entire structure. MRM does not require the selection of substructure boundary conditions or a substructure connection technique. Three models were proposed in the original study on MRM [44], including two-dimensional (2D), 2D+torsion, and three-dimensional models. The 2D model can be used for slender cantilever structures subjected to bending forces, whereas the 2D+torsion model is suitable for the analysis of slender cantilever wings because the torsion angles determine the aerodynamic forces acting on the wings.

Originally, MRM was proposed as a static analysis method [44] before being applied to frequency-domain flutter analyses [45]; however, MRM has not yet been extended to time-domain dynamic analysis. The time-domain dynamic analysis is essential for designing aircraft with high-aspect-ratio wings. In 2003, an experimental airplane HELIOS with a high-aspect-ratio wing crashed during a flight test because of an unexpected dynamic behavior caused by a large deformation, as reported by the National Aeronautics and Space Administration [46]. According to this report [46], one of the key recommendations to prevent such crashes is to develop more advanced, multidisciplinary (structures, aeroelastic, aerodynamics, atmospheric, materials, propulsion, controls,etc.) “time-domain” analysis methods. Therefore, extending MRM to a time-domain dynamic analysis method can be valuable. Two problems must be resolved to realize this extension.

The first problem is that MRM has been formulated only for static analysis, and a discretized dynamic equation and its solution procedure have not yet been developed. The second problem is that the formulation of the energy stored in a structure and the work performed by external forces have not been established in MRM. In other words, it is unclear whether MRM satisfies the energy conservation law. Therefore, the work performed by the follower external force should be formulated by considering the change in the force direction caused by a large deformation. The follower force constantly changes direction to follow the deformation of the structure. The antonym of “follower force” is “dead force.”

In this study, we develop a novel time-domain dynamic analysis framework referred to as dynamic MRM (DMRM) to analyze slender cantilever structures, such as high-aspect-ratio wings, by solving the two aforementioned problems. The proposed DMRM seamlessly connects the high-fidelity static analysis model of strength design with the low-fidelity dynamic analysis model of flight simulation, controller design, and gust response analysis. The proposed method targets slender structures with small strains and large displacements and considers geometrical nonlinearity, but not material nonlinearity. Studies on the dynamic analysis of the high-aspect-ratio wings, which is the subject of the proposed method, often assume slender structures with small strains and large displacements [42,47]. Therefore, geometrical nonlinearity is more dominant and important in the deformation of high-aspect-ratio wings, and it can be solely considered in this study. We consider two research objectives for the development of DMRM.

The first objective is to propose a discretized dynamic equation and its time integration process. A major advantage of MRM is that geometrical nonlinearity can be considered by recursively rotating the segments based on curvature modes. We utilize this to calculate the modal acceleration based on the modal coordinates and velocities obtained at certain times.

The second objective is to derive an expression for the energy of a structure and to propose a formulation method for the work performed by follower forces. Geometrical nonlinearity is considered as a force term rather than a stiffness term in MRM. In this study, geometrical nonlinearity is considered in the calculation of the work performed by external forces in the same manner. Additionally, we suggest that the work performed by follower forces can be calculated accurately by averaging the directions of the forces with respect to time. These formulations demonstrate that DMRM satisfies the energy conservation law.

The remainder of this paper is organized as follows: In Sec. 2, the theory of the proposed DMRM is presented. In Sec. 3, the simulation results are presented to demonstrate that DMRM derives an efficient dynamic analysis model from linear models and reduces the calculation time. The investigation of appropriate reduction methods required for the stiffness and mass matrices is also presented in this section. Finally, the conclusions of the study and the future scope are presented in Sec. 4.

In this section, the 2D and 2D+torsion models of MRM are extended to dynamic analyses. Before describing the proposed DMRM, the original MRM developed for static analysis [44] is explained using a relatively simple 2D scenario (Fig. 2). This figure depicts an example in which two finite elements and four segments are considered. The structure is divided into finite elements to obtain the displacement modes, and a reference line is selected along the longitudinal direction of the structure. The structure is then divided into segments. While finite elements are used only to obtain the mode shapes, segments are used to represent large deformations. The modal shape at each node of the finite element model is interpolated to each segment.

Here, ΔF_mode(ξ) denotes the force term that corrects the generalized force according to the deformed shape. Equation (4) can be solved for the modal displacement ξ.

where Φ_dθ is composed of selected curvature modes on the reference line and dθ denotes a vector, whose entry is the relative rotation angles of each segment. Just as the widely used general displacement mode is a basis for calculating displacement, the curvature mode is a basis for calculating curvature. By utilizing the curvature technique, large deformations are analyzed based on modes obtained from a linear finite element model, enabling to perform the analysis without using the complex concept of nonlinear modes. The large deformation shape is described by rotating the segments based on the rotation angles in dθ, as explained in Ref. [44] and Appendix  A. The generalized force term ΔF_mode(ξ) can be computed based on the geometrical relation in the deformed shape, as presented in Ref. [44] and Appendix  B.

where M denotes the mass matrix of the linear/linearized model used to compute the displacement modes. The inertial nonlinearity caused by the relative segment rotations is assumed to be negligible in Eq. (8). We demonstrate that this is negligible when a high-aspect-ratio wing exhibits nearly 30% tip deflection against the wingspan in Sec. 3. A 30% tip deflection may occur in relatively severe flight conditions where an aircraft with high-aspect-ratio wings encounters gusts.

Figure 4 depicts the dynamic analysis procedure performed using this equation. The subscript i represents the discrete i-th time point. When the modal displacement ξ_i is given and the modal acceleration $ξ¨i$ is calculated, the modal displacement ξ_i+1 at the subsequent time-step can be determined using a time integration scheme, such as the Euler's method or Runge–Kutta method. However, ΔF_mode(ξ_i) must be calculated to obtain the modal acceleration using Eq. (11).

The acceleration $ξ¨i$ can be calculated based on ξ_i. We capitalized on the advantage of MRM that geometrical nonlinearity can be considered based on deformed shapes. After dθ is calculated using the linear superposition of the curvature modes based on ξ_i, the deformed shape at the i-th time point is calculated based on dθ using the method described in Appendix  A. Then, ΔF_mode(ξ) is calculated based on the deformed shape using the method described in Appendix  B. Therefore, the modal acceleration can be calculated using Eq. (11).

The 2D+torsion model assumes a small torsion deformation independent of the bending deformation. Even when a 50% wingspan deflection occurs, this assumption provided accurate results [44]. Torsion deformations were analyzed based on the linear superposition of torsion modes. The torsion deformations are expressed by rotating the segments based on the obtained torsion angles after locating the segments using the obtained curvatures dθ.

The displacement of each segment that constitutes the shape of the structure is not directly obtained using time-integration schemes. In DMRM, the modal displacement ξ is obtained via a time-integration scheme for each discretized time using Eq. (10). The displacement of each segment is obtained based on the curvature dθ calculated using the modal displacement ξ.

The proposed DMRM can solve the equation of motion on modal coordinates, unlike general nonlinear finite element methods, wherein the equation of motion is solved on generalized coordinates. The three advantages of using modal coordinates are as follows.

First, dynamic analysis can be performed with a few degrees-of-freedom by employing dominant modes. Even if the complexity of the structure increases, the first advantage is retained because the analysis uses the mode shapes that exist only on the reference line. Conversely, in general nonlinear finite element methods, the degrees-of-freedom of the equation of motion increase when the structure needs to be modeled with numerous elements owing to structural complexity.

Second, a time-step is systematically determined based on the periods of the modes used in the analysis. The periods of the modes are calculated using a linear finite element model. Dynamic analysis can be performed with the largest possible time-step within a value that ensures the accuracy of the analysis. A large time-step can be set because DMRM can exclude high-order modes from the analysis, which reduces computational costs. In contrast, the time-step cannot be systematically determined in the case of general nonlinear finite element methods.

Third, the characteristics of the oscillations included in the large dynamic deformation can be determined based on the ratio of each mode. The time history of each modal displacement can be obtained because DMRM is performed on the modal coordinates. Therefore, the dominant modes in the dynamic deformation can be determined, and the information on the vibration modes is useful in the controller design. Conversely, nonlinear finite element methods cannot obtain this information directly.

In this section, the equation for the energy stored in a structure at each time point is derived, and a formulation method is proposed for the work performed by the follower forces at a certain time period. We suggest that the work performed by follower forces can be calculated by averaging the directions of the forces over time. This method cannot be used when the directions of the forces change drastically owing to the deformation occurring in a short time caused by an impact load.

The work performed by the external forces during a certain time period agrees with the energy that is newly stored in the structure during that time period. If an agreement is confirmed in the analysis, the formulation of the DMRM can be validated from the viewpoint of energy conservation.

In the case of elastic energy (Eq. (13)), we use the linear part of the stiffness matrix $Kmodelinear$ obtained from Eq. (2). We assume that the local strain at each position is small and that the stiffness matrix $Kmodelinear$⁠, which is inherent to the structure, does not change according to the deformations. Therefore, the stiffness matrix $Kmodelinear$ is used to calculate elastic energy. In the case of large deformations, the work performed by external forces during a certain time period is affected by geometrical nonlinearity because the positions and directions of the external forces change significantly. Therefore, geometrical nonlinearity must be considered in the calculation of the work performed by external forces.

The calculation of the position vector for each segment in the global coordinates is presented in Appendix  A. The work W_i performed by the dead force by the i-th time point is computed using the position of the segment where the force is acting at that time point.

Three simulation test cases of cantilever structures were performed to confirm the validity and advantages of DMRM. Here, the energy and work performed by the external forces as well as the deformation were analyzed. The calculation time of DMRM was compared with that of nonlinear finite element methods based on ANCF. As reported in Ref. [48], some ANCF beam elements cannot perform accurate analysis when subjected to torsional load or off-axis local concentrated forces. This is because when the number of element divisions is insufficient, points where the displacement is not smooth occur between elements. Therefore, in this study, ANCF beam elements were only used for comparative verification when the element is not subjected to torsional load or off-axis local concentrated forces. In the first two test cases, we used the ANCF beam element model proposed in Ref. [49], where the longitudinal elastic forces were simplified using the L1 model. This element is a thin beam element that has no interpolation in the y and z directions but has interpolation in the x direction based on a cubic polynomial. This element is suitable for analyzing deformations where bending is dominant. In the third test case, the ANCF plate element model proposed in Ref. [50] was used. The elastic forces of this element were formulated under the assumption that the in-plane deformation is smaller than the out-of-plane deformation. It has been shown in Ref. [50] that this element can analyze large out-of-plane deformations. All analyses of DMRM and ANCF were performed using the fourth-order explicit Runge–Kutta method to exclude the effects of the tuning parameters in the implicit methods. The generalized mass and stiffness matrices used in the DMRM analysis were obtained from the linearized modes, mass matrix, and stiffness matrix obtained based on a linear finite element model. The stiffness matrix obtained from a linear finite element model can be used to calculate $Kmodelinear$ within the assumption of small strains. In addition, a mass matrix obtained from a linear finite element model can be used when the inertial nonlinearity caused by the relative rotations of the segments is negligible. The analyses were performed using a single core of the AMD Ryzen 9 5900X Processor at 3.70 GHz and 128 GB RAM and MATLAB R2023a. The calculation time was measured three times in each analysis and averaged. The initial conditions of deformation and velocity of the cantilever structures were set to zero in all three test cases. Reference lines for the cantilever structures were set on the neutral axis.

In the first test case, a 2D model of DMRM was used to analyze the beam deformation caused by a tip dead force, as depicted in Fig. 7. This beam with an I-shaped cross section was used in the original study on MRM [44]. Table 1 lists the beam parameters. The dead force applied was 1000 N. In this test case, 400 linear beam elements were used to obtain the modes, mass matrix, and stiffness matrix for the DMRM analysis. The number of segments in DMRM was set to 800. The first to fifth modes calculated using the linear beam model were employed in DMRM. The period of the fifth mode of the beam was 0.13 s. Therefore, the time-step was set to 1.0 × 10⁻³ s, which was nearly 1/100 of the order of the period. Figure 8 shows the time history of the modal displacement of each mode obtained from the DMRM on a logarithmic scale. The modes normalized by the mass matrix were used in the analysis of modal displacements, where Φ^TM_modeΦ = I. Figure 8 indicates that the lower-order modes are more dominant than the higher-order modes. Therefore, the first and second bending modes were employed in the subsequent DMRM analyses.

The second bending mode period calculated by the linear model was 1.2 s. Therefore, the time-step was set to 1.0 × 10⁻² s for DMRM. During the ANCF analysis, the beam was divided into 10 beam elements and the time-step was set to 1.0 × 10⁻⁴ s. This is because numerical divergence occurred when a time-step of more than 1.0 × 10⁻³ s was employed for ANCF. Figure 9 shows the time histories of the tip positions analyzed using DMRM and ANCF, wherein the results of DMRM agree with those of ANCF. To investigate the effects of geometrical nonlinearity, the results of the linear finite element method (LFEM) are also presented in Fig. 9. The results of LFEM differ from those of DMRM and ANCF, indicating that both DMRM and ANCF can consider geometrical nonlinearity. The result of the LFEM was always x = 30 m because the dead force was applied in the z-direction.

Figure 10 shows the time history of the deformed shapes obtained via DMRM and ANCF. The DMRM results are in good agreement with the ANCF results. The calculation time required for DMRM was 13 s, whereas that required for ANCF was 35 s. DMRM achieved a 63% reduction in the calculation time compared with ANCF because DMRM had fewer variables and a larger time-step than ANCF.

The beam does not contain any energy in its initial state. Therefore, the validity of the analysis from the standpoint of energy is demonstrated if the work performed by the tip force at a certain time period agrees with the energy stored in the structure in that period. Figure 11 shows the time histories of the energies and work calculated using DMRM. The energy stored in the structure during a certain time period is plotted using positive values, and the work performed by the tip force during that period is plotted using negative values. The sum of the energies and work is zero at each time point. This confirms the validity of the proposed calculation method for the work performed by the dead force.

The RMS error changes little above 800 divisions, which indicates the convergence of the analysis with respect to the number of segments.

In the second test case, the deformation of a relatively complex-shaped plate caused by the dead force was analyzed (Fig. 13). The material of the plate was vinyl chloride with a density of 1425 kg/m³, Young's modulus of 2.9 GPa, and a Poisson's ratio of 0.39. The tip force applied was 0.27 N.

A linear plate element with 12 degrees-of-freedom [51] was used to obtain the modes, mass matrix, and stiffness matrix required for the DMRM analysis. Each node of the linear plate element comprises three variables, namely, the displacement in the out-of-plane direction, the gradient of the out-of-plane displacement along the longitudinal direction, and the gradient of the out-of-plane displacement along the direction of the width. The gradients of the out-of-plane displacements along the longitudinal direction can be treated as deflection angles on a reference line. The gradients of the out-of-plane displacements along the width can be treated as torsion angles. In this test case, only the bending deformation was considered because the loading condition in Fig. 13 did not cause torsion deformation.

The complex plate was modeled by assigning plate elements with extremely low stiffness to the perforated areas to obtain the modes of the plate easily (Fig. 14). The figure shows the top view of the plate and the number of divisions of the linear plate elements in each direction. The plate was evenly divided into 100 and 4 elements in the longitudinal and width directions, respectively, including the perforated areas. The Young's modulus, density, and Poisson's ratio in the perforated areas were set to 1.0 × 10⁻³ times the corresponding parameters of vinyl chloride.

Before performing the MRM analysis, the mass and stiffness matrices should be obtained from the linear plate model. However, the number of degrees-of-freedom of the matrices does not agree with those of the segments in MRM. For instance, when a structure is modeled using linear plate elements (Fig. 15), the nodes are distributed in the plane of the structure. The number of degrees-of-freedom in the mass and stiffness matrices corresponds to the number of nodes. However, MRM uses mass and stiffness matrices and the number of degrees-of-freedom corresponding to the number of segments divided only along the reference line. Therefore, the matrix obtained from the linear finite element model should be reduced to a matrix of size corresponding to the number of segments.

For comparison, two reduction methods were used for the DMRM analysis, namely, the improved reduced system (IRS) [52] and the Guyan reduction [53]. In the Guyan reduction method, both stiffness and mass matrices can be reduced by extending the Guyan reduction theory [54]. The number of segments was set to 400, and the first and second bending modes were used in the analysis. The period of the second bending mode of the plate was approximately 5.4 × 10⁻² s. Therefore, the time-step was set to 5.0 × 10⁻⁴ s for DMRM. During the ANCF analysis, the time-step was set to 1.0 × 10⁻⁵ s to prevent numerical divergence, and the plate was divided into eight beam elements (Fig. 16). Although this plate can be analyzed using a beam model, we aim to demonstrate that DMRM analysis can be performed using the modes and the mass and stiffness matrices obtained from the linear plate model rather than the linear beam model. When the structure exhibits a more complex geometry, such nonlinear beam modeling becomes difficult. In contrast, DMRM can analyze a complex geometry after the mode shapes and the reduced mass and stiffness matrices are obtained from a linear shell/plate or solid model.

Figure 17 shows the results of the tip positions analyzed using the four methods: DMRM based on the matrices reduced by IRS, DMRM based on the matrices reduced by Guyan reduction, ANCF beam model, and linear plate element with 12 degrees-of-freedom [51]. The results of DMRM (IRS) and DMRM (Guyan) are consistent with the results of ANCF, whereas the results of the linear plate element are not consistent with them. Therefore, both IRS and Guyan reduction can be used as reduction methods for DMRM. A minor difference exists between the two DMRM results and the ANCF result at the x position. This is because, unlike ANCF, DMRM does not consider axial deformation. The calculation time of DMRM (IRS) using the first two modes was 6 s, whereas that of the ANCF beam element was 19 s. DMRM (IRS) achieved a 68% reduction in the calculation time compared with the ANCF beam model.

Figure 18 shows the time histories of the energy and work calculated by DMRM (IRS). The sum of both values is zero at each time point, validating the proposed calculation method for the work performed by the dead force.

In the third test case, the 2D+torsion model of the DMRM was used to analyze the deformation of a simple plate (Fig. 19). The plate had a span of 0.4 m, width of 0.04 m, and thickness of 0.002 m. The material parameters of the plate were identical to those of the plate described in Sec. 3.2. The applied force was modeled as a distributed follower force. A distributed force of 1.6 N/m was applied at a quarter of the width away from the neutral axis. This loading condition simulated the aerodynamic forces. The bending and torsion deformations of the simple plate were analyzed using the DMRM 2D+torsion model. The loading conditions may have caused bending in the width direction; however, this deformation can be neglected because the bending stiffness in the width direction is larger than that in the thickness direction.

The plate was modeled using the linear plate model [51] to provide the modes and the mass and stiffness matrices for DMRM. The plate was divided into 100 elements in the longitudinal direction and four elements in the width direction. The number of segments was set to 400. The mass and stiffness matrices obtained by the linear plate model were reduced by IRS. The initial two modes of both bending and torsion were used for the analysis. The period of the second bending mode of the plate was approximately 5.5 × 10⁻² s. The period of the second torsion mode of the plate was approximately 6.0 × 10⁻³ s. The time-step of DMRM was set to 5.0 × 10⁻⁵ s, which was nearly 1/100 of the torsional period. The first and second bending modes were the first and second modes of the plate, respectively.

For comparison, the plate was divided into 15 ANCF plate elements [50] along the longitudinal direction. The time-step was set to 1.0 × 10⁻⁶ s.

The positions of the tip on the leading edge (LE) and the trailing edge (TE) calculated using DMRM were compared with those calculated using the ANCF plate elements (Fig. 20). The difference between the LE and TE positions is shown in Fig. 20 to compare the torsion deformations. Although all results are in good agreement, a minor discrepancy is observed between the DMRM and ANCF plates. This can be attributed to the ANCF plate model expressing the width-direction deformation, which is not observed in DMRM.

Figure 21 shows the time histories of the energy and work calculated using DMRM. As their sum is always zero, the proposed calculation method for the work performed by the follower force is valid.

The calculation time of the DMRM 2D+torsion model was 737 s, whereas that of the ANCF plate element was 13,820 s. DMRM achieved a 95% reduction in the calculation time in comparison with the ANCF plate model in the case of a simple plate. This indicates that DMRM can reduce the calculation time compared with high-fidelity nonlinear shell/plate or solid models in dynamic analysis when applied to a practical wing shape that requires the discretization of more elements.

In this study, we developed a novel time-domain dynamic analysis framework, DMRM, for slender cantilever structures. The proposed method targeted slender structures with small strains and large displacements and considered geometrical nonlinearity, but not material nonlinearity. In this framework, we proposed a discretized dynamic equation of MRM and its solution by capitalizing on the fact that the geometrical nonlinearity could be considered based on the deformed shape of the structures. Additionally, we proposed a calculation method for the energy stored in the structure and the work performed by external forces. By averaging the direction of the force vectors in terms of time, the work performed was calculated even in the case of follower forces.

Three simulation cases were presented. First, the beam was analyzed under a dead force, wherein the results of the DMRM 2D model agreed with those of the ANCF beam model. DMRM achieved a 63% reduction in the calculation time compared with the ANCF beam model owing to the effect of mode reduction. Moreover, the energy stored in the structure agreed with the work performed by the dead force.

Second, a complex-shaped plate under a tip-dead force was analyzed. DMRM succeeded in the dynamic large-deformation analysis of a complex structure using only a few modes on the reference line. The agreement between the work performed by the dead force and energy stored in the plate was also confirmed in this complex structure. We further demonstrated that no difference existed in the analysis accuracy between the two reduction methods of the mass and stiffness matrices, namely, Guyan reduction and IRS, used for DMRM.

Third, the large bending deformation and the small torsion deformation caused by the distributed follower force simulating the aerodynamic force were analyzed using the DMRM 2D+torsion model. The bending and torsion deformations agreed with those of the ANCF plate model. The agreement between the work performed by the follower force and energy stored in the structure was also confirmed. Therefore, the proposed averaging method for the direction of the follower force with respect to time was validated. DMRM achieved a 95% reduction in the calculation time compared with the ANCF plate model. This indicates that DMRM is more effective than the high-fidelity nonlinear finite element model for dynamic analysis.

The proposed method has restrictions on the strain magnitude for the stiffness term and on the magnitude of the effect of the relative rotation of the segments for the inertial term. The three test cases reported in Sec. 3 revealed that even with these two restrictions, the proposed method can accurately analyze large deformations. Owing to these two restrictions, the proposed method may fail to analyze structures such as extremely flexible structures undergoing folding back. Therefore, the analysis of such structures will be the subject of future research. In the future, we intend to introduce the nonlinear inertia effects in DMRM for a more accurate analysis. Furthermore, the three-dimensional model of MRM should be extended to dynamic analysis.

• Japan Society for the Promotion of Science (JSPS) KAKENHI (Grant Nos. 21K14341, 22K18853, and 24K17440; Funder ID: 10.13039/501100001691).

• Nohmura Memorial Foundation.

• Suzuki Foundation.

• Kurita Water and Environment Foundation (Funder ID: 10.13039/501100001698).

The datasets generated and supporting the findings of this article are obtainable from the corresponding authors upon reasonable request.

dl =

vector representing geometries of segments before deformations

dR =

rotation matrix for adjacent segments

dθ =

vector of change in the deflection angle at each segment

EM =

external bending moment acting on segments

$F¯$ =

external force vector averaged over time

F =

external force vector

F^l =

external force acting on the l-th segment

F_mode =

generalized external force vector

ΔF =

vector to correct the difference between external forces before and after deformations

ΔF^k =

vector that corrects the difference between external forces before and after deformations in the k-th segment

ΔF_mode(ξ) =

generalized external force vector correcting moments

K =

stiffness matrix

K =

kinetic energy

K_mode(ξ) =

nonlinear generalized stiffness matrix

$Kmodelinear$ =

linear component of generalized stiffness matrix

$ΔKmodenonlinear(ξ)$ =

nonlinear component of generalized stiffness matrix

M =

mass matrix

M =

internal bending moment acting on segments after deformations

M₀ =

internal bending moment acting on segments before deformations

M_mode =

generalized mass matrix

ΔM =

moment added to segments to correct bending moments

R =

rotation matrix of segments with respect to global coordinates

RMS =

root-mean-square of the z displacement

RMS error =

error of the root-mean-square

U =

elastic energy

U_j =

position vector at the j-th discretized time of points where external force acts

U^k =

position vector after deformations of the k-th segment

W =

work performed by external force

X =

position vector of segments before deformations

Δx =

segment length

z_q =

displacement in z direction at the q-th discretized time of point

θ =

deflection angle

ξ =

modal displacement

Φ =

displacement modes matrix

Φ_dθ =

curvature modes matrix

Subscripts and Superscripts

()_j =

j-th discretized time

()_x, ()_y, ()_z =

x, y, z directions

()^dead =

scalar or vector of dead force

()^follower =

vector of follower force

()^k =

k-th segment