1. Introduction
Aircraft design is moving to increasingly slender and lightweight wings for higher efficiency and lower CO
2 emissions. These structures bring new design challenges due to their high flexibility and high deflections under normal loads. The geometric nonlinearities introduced with large deformations require models to predict the effects on the aeroelastic phenomena. The importance of aerodynamic and structural geometrical nonlinearities in the aeroelastic behavior of high-aspect-ratio wings has been established by Patil and Hodges in [
1]. Patil et al. [
2] have looked at the effect of structural geometric nonlinearities on the flutter behavior of high-aspect-ratio wings, and they presented the changes in the structural and aeroelastic characteristics of a steady-state deflection of a wing. Their study revealed a significant change in the structural frequencies and a significant reduction in the flutter speed.
Detailed coupled computational fluid dynamics and finite element method formulation for aeroelastic analysis have been widely studied [
3]. These models can be very sophisticated and generally require a large number of calculations which are not efficient for design and optimization. For this reason, the research moves to low-order aeroelastic models, which consent to reduce the computational cost with similar prediction capabilities if compared to high-order models. A popular approach for nonlinear elasticity consists in geometrically exact beam formulation; an example is given by Hodges [
4,
5], who presented an intrinsic formulation for nonlinear dynamics of initially curved and twisted anisotropic beams. These models use equivalent beam properties derived from finite element models (FEMs) [
6,
7]. These formulations have been used in several works for studying very flexible wing structures: Drela [
8] used nonlinear beams to develop an integrated model for aerodynamic and structural simulation of flexible aircraft, while Patil [
9] presented a theory for flight dynamic analysis of highly flexible wing configuration, accounting for geometrically nonlinear structural deformations. A more recent class of low-order structural models relies on high-order modal expansions [
10,
11]. These models require nonlinear static responses of an FEM to identify modal expansion terms. Another model has been presented by Bruni et al. [
12], who expanded the partial differential equations for the beam dynamics up to the third order to explore the effects of static deflection, external trim, gust loads, and aerodynamic stall. The solution was obtained with Galerkin’s method and with a multi-modal approach. Other models for static and dynamic nonlinear analysis of beam structures exist in the literature and employ different solutions to simulate specific conditions. Some models consider only one-dimensional finite elements as in [
13], while other developed models consider all the possible degrees of freedom. Yang et al. [
14] developed a six degree of freedom beam element, including material nonlinearity; they described a procedure for nonlinear static analysis which involves piecewise linearization of the response equations and iterations at each incremental step to achieve equilibrium. Surana et al. [
15] presented a geometrically nonlinear formulation for a 3D curved beam element using a Lagrangian approach and verified the accuracy and efficiency of the formulation with literature results of nonlinear static analysis. More recently, Jin and Yun [
16] developed a three-dimensional beam element for geometrically nonlinear dynamic analysis; the derivation is based on the co-rotational formulation. The model showed good results and can undergo large deflection and rotations, but small strains are assumed.
Low-order beam structural models can be further developed to consider also material couplings and expand the aeroelastic design domain. Anisotropic materials can be used to improve wing box structural performance according to the concept of aeroelastic tailoring [
17,
18]. The advantages of aeroelastic tailoring are enhanced by orthotropic materials. Bakthavatsalam [
19] demonstrated the effect on the flutter speed of aeroelastically tailoring wing and tail surfaces of a closely coupled wing–tail flutter model, and it was shown that tailoring the wing surface produced the largest increase in flutter speed, but tailoring the tail and reducing its stiffness could also produce an increase in flutter speed. Weisshaar [
20,
21,
22] focused on the use of laminated composites to increase the divergence speeds of swept forward wings. Weisshaar included bending–torsion coupling, defining a stiffness parameter that describes the amount of interaction between the bending curvature and twist rate. This parameter is a function of the orientation and stacking sequence of symmetrical laminate plies with respect to a reference axis along the wing. Composite panels layups can be studied to achieve specific aeroelastic performance and considering also functionally graded materials [
23,
24,
25,
26] and variable angle tow [
27,
28,
29,
30]. However, aeroelastic tailoring is not limited to composite materials, as several studies [
31,
32,
33] have shown that the arrangement of stiffeners can be used to control directional stiffness and bending–torsion coupling.
Beam finite elements are particularly suitable for high-aspect-ratio wing analysis. However, traditional beam elements do not consider the bending–torsion couplings granted by the use of oriented orthotropic materials. Recently, Patuelli et al. [
34] presented a beam finite element with bending–torsion coupling formulation (BTCE). The linear finite element was derived with Galerkin’s method, while the bending–torsion coupling was obtained with specific shape functions and the hypothesis of constant torsional moment along the beam element. The resulting model was validated with experimental and numerical results [
34,
35], showing good accuracy for static and dynamic analysis. The scope of the present research work is to develop a procedure to perform dynamic analysis in the presence of geometric nonlinearities. Cestino et al. [
36] studied the flutter instability of high-aspect-ratio wings and considered the phenomenon as the sum of two effects, the geometrical effect (GE) given by the deformed geometry and the stiffness effect (SE), which is the effect caused by the loads at the equilibrium condition on the differential stiffness matrix. They demonstrated that the GE is the main contribution in the nonlinear dynamic analysis of slender structures and that the results of the flutter analysis are verified by experimental evidence either when considering only GE or when accounting for SE.
This research work presents the derivation of a nonlinear beam finite element with bending–torsion coupling formulation (BTCE-NL), which considers both GE and SE. The element stiffness matrix is derived considering the nonlinear terms through the perturbation of a known equilibrium configuration. Moreover, an approach that accounts only for the GE (BTCE-GE) has been developed considering a deformed equilibrium dependent transformation matrix to orient the BTCE. The derived models have been validated with several experimental modal analyses performed with Laser Doppler Vibrometer (LDV). The experimental data gathered have been used to validate the characteristic frequencies and the mode shape predicted with the new models, but they also consented to understand at which level of deflection the dynamic linear analysis is no longer suitable for mode shape prediction, and advanced models are needed. The experimental tests considered a beam structure with a coupling coefficient equal to zero to avoid flapwise vibrations. The results revealed that both models can predict characteristic frequencies and modes with good precision for moderate-to-large deflections and that the BTCE-GE can be sufficient to analyze slender structures with moderate initial deflections. The models has been tested also for a composite beam-box structure with a circumferentially asymmetric stiffness (CAS) configuration described in [
37]. The experimental modal analysis of a beam structure with bending–torsion couplings involves both edgewise and flapwise coupled modes and would require more sophisticated equipment; for this reason, the BTCE models results have been compared with the results of a SHELL FE model solved with NASTRAN SOL106. The methods adopted showed good correlation with the SHELL FE model for moderate deformations. The models presented in this research can be used to perform optimization cycles with low computational costs and find layup configurations able to mitigate the deformation-induced nonlinear effects.
3. Experimental Validation for Isotropic Beam
The models here derived have been validated through experimental modal analysis. The tests were performed on a rectangular section aluminum 6060 beam with dimensions
mm,
mm, and
mm (
Figure 3) and the mechanical properties listed in
Table 1. The beam was clamped in four different positions to gather data of four cases, respectively, with useful length
mm,
mm,
mm, and
mm (
Figure 4). We defined the ratio
with
equal to the tip deflection; one of the scopes of the experimental test is to understand at which level of
the geometric nonlinear effects have an influence on the beam mode shapes and characteristic frequencies, determining the need of a nonlinear modal analysis. The other objective for the experimental testing is to verify the accuracy of the presented models in predicting characteristic frequencies and mode shapes.
The beam was investigated with four experimental tests with similar equipment. The beam was clamped between two steel blocks to guarantee a rigid constraint (
Figure 5) at the first section of the beam. Ten polylactic acid (PLA) targets (
Figure 5) were placed along the beam to gather data at ten equidistant stations. The number of targets was limited to ten units to keep the additional weight on the beam negligible. More targets can be added to improve the acquisition resolution, but the mass and inertia must be considered and can alter the nonlinear effects observed. Each target presented two vertical parts for acquisition, where a squared piece of reflective tape was positioned to improve the surface reflectiveness.
The acquisition was performed with a Polytec PSV-500 Laser Doppler Vibrometer (LDV) system, while the excitation was obtained with an electrodynamic shaker K200xE01. The shaker was placed at
mm from the constraint, perpendicular to the beam axis as represented in
Figure 5. The objective was to excite only the edgewise degree of freedom because edgewise and torsional characteristic frequencies are the most affected by flapwise deflection; moreover, the torsional modes should be visible only when the nonlinear coupling effect becomes important according to [
1].
The experimental validation considers the case where the coupling term K is equal to 0. This allows to reduce the number of variables and keep the interpretation of the results straightforward. Moreover, the term K induces bending–torsion deformations, which means that a mode shape which involves the torsional degree of freedom determines also the flapwise displacements, which need a 3D LDV system to be detected.
Four numerical models were defined for experimental result comparison. Two models were implemented with MATLAB and used the BTCE finite element. One used the linear formulation derived in [
34], accounting for the geometrical effect (BTCE-GE), while the second used the nonlinear BTCE (BTCE-NL). Two additional models were defined in PATRAN and solved in NASTRAN starting with an undeformed configuration; both used SHELL elements to describe the beam geometry, but one was solved with SOL103 and the second one was solved with the nonlinear solution SOL106. The linear modal analysis was performed on the undeformed configuration to obtain linear mode shapes and frequencies for the result comparison. The choice behind the use of SHELL finite elements is the possibility to add bending–torsion coupling terms, which is possible for the BTCE but not for conventional beam elements. The BTCE models were obtained by assembling 10 elements which represent the 10 segments described by the targets positioned on the experimental beam. The first node was constrained, imposing the translations and the rotations equal to 0. For the BTCE-GE model, the modal analysis was performed using the linear stiffness matrix rotated with the equilibrium configuration-dependent transformation matrix
described in the second section of this work. The BTCE-NL model used the stiffness matrix derived in the second section, which depends on the equilibrium static deformation. The mass of each element was lumped at the nodes, and a linear static analysis determined the equilibrium deformation used to complete the element stiffness matrix and perform the nonlinear modal analysis. Alternatively, the deformed configuration can be obtained with a nonlinear static analysis performed with NASTRAN. The linear static analysis for a vertical load does not present edgewise displacements or rotation, while the nonlinear static analysis have a small in-plane component
and
, which can be considered negligible. The SHELL elements was created with 10 mm QUAD4 elements and then solved with SOL103 for the linear modal analysis in the undeformed configuration as reference. The model was completed with an inertial load to perform also the nonlinear modal analysis SOL106 that accounts for the preload. The numerical models results were compared with the experimental results in terms of mode shapes and characteristic frequencies. The linear analysis was performed to understand at which level of deformation it becomes unreliable and a nonlinear formulation becomes needed.
4. Numerical Models Comparison for Composite Beam
A numerical comparison was performed for a case with bending–torsion coupling; the reference structure is a box beam structure with a circumferentially asymmetric stiffness (CAS) laminated composite configuration. The structure was described in [
37]; the section is represented in
Figure 6. The beam was obtained with a unidirectional T700 carbon–epoxy layer bonded onto wooden spars with fibers oriented at
. The structural box had the following dimensions: length,
L = 522 mm; width, w = 20 mm; height, h = 2.8 mm; upper and lower panel thicknesses, t = 0.2 mm; mass per unit length, m = 1.095 ×
kg/mm; and torsional unit inertia, Ip = 4.75 ×
kg/mm. The mechanical properties of the material are listed in
Table 2.
The reference model was defined in PATRAN with SHELL elements (
Figure 7), while a beam model with the formulation presented in this work was obtained assembling 10 elements. The load condition chosen for the numerical comparison is a concentrated tip load. The load was incremented to reach different deformation levels and observe the limits of validity of the presented model. The deformation was evaluated with a nonlinear static analysis, then a nonlinear modal analysis was performed for each load case, and the numerical results were compared in terms of characteristic frequencies. The deformed configuration used to orient and compute the nonlinear beam finite element was retrieved from the nonlinear static analysis performed with NASTRAN. The first eight characteristic frequencies were computed for the two FE models and normalized with the value obtained with a linear modal analysis of the undeformed configuration. The normalized frequencies were compared for each mode at different deformation levels.
5. Results and Discussion
This section presents the experimental evidences collected and compared with the numerical models. The dynamic behavior of the nonlinear BTCE was compared in terms of predicted natural frequencies and mode shapes. The accuracy of the natural frequencies was evaluated in terms of relative difference. The similarity between the FE models and the experimental mode shapes was evaluated with the Modal Assurance Criterion (MAC). MAC is a statistical indicator used to quantifying the similarity between two sets of mode shapes, where a value equal to 1 indicates complete similarity, while 0 indicates no correlation between the modes investigated [
39,
40]. Equation (
25) was applied to the experimental and numerical mode shapes computed to obtain the MAC matrices. The mode sets of the experimental mode shapes were compared with themselves computing the Auto MAC, allowing to verify the existence of similarities between different mode shapes and thus the presence of couplings between the degrees of freedom. The couplings, if present, should show the same pattern for experimental and nonlinear numerical modes. When the structure does not present couplings, the expected matrices for the experimental and numerical linear and nonlinear modes should be diagonal:
5.1. Static Analysis Results
The nonlinear finite element derived depends on the equilibrium deformation under static load. The deformation can be obtained through linear or nonlinear static analysis. In this research, a linear static analysis is used to determine the initial equilibrium deformation for the load cases considered during the experimental tests on the isotropic beam. The results of the deflection at the tip were compared with the result of a SHELL model of the beam solved with SOL106 and experimental results, and the accuracy was evaluated computing the relative difference between numerical and experimental results. The comparison is reported in
Table 3 with the relative difference for each result.
5.2. Experimental Modal Analysis Results
The Frequency Response Functions (FRFs) obtained through experimental modal analysis are reported in
Figure 8. The experimental and numerical results for the characteristic frequencies are reported in
Table 4. For a beam length equal to 1000 mm, the first torsional mode was not detected, while for a length of 1500 mm, the torsional mode was detected, but the peak was significantly smaller than the others. This confirms that the coupling between edgewise bending and torsion is weak for deformations below
. On the other hand, for bigger deformations, the excitation of the edgewise degree of freedom provoked also the detection of the first torsional mode coupled with the edgewise bending mode.
The frequencies reported in
Table 4 show good accordance with the predicted values and the experimental results. In this case, both linear and nonlinear models can be use to determine the characteristic frequencies of the structures. The relative difference between the predicted and observed frequencies, reported within parentheses in
Table 4, is generally below
with some exception compatible with the approximations introduced with the derivation of the BTCE models. Moreover, the difference between the BTCE-GE model and the BTCE-NL model are minimal, confirming the findings reported in [
36,
41] concerning the major contribution of the geometrical effect in this class of analysis.
The experimental Auto MAC matrices are reported in
Figure 9A,
Figure 10A,
Figure 11A and
Figure 12A. The mode order is based on the frequency value, from the mode with the lowest frequency to the one with the highest. With this convention, for
L = 1000 mm and
L = 1500 mm, the torsional mode occupies the third position, while for the other cases, it is placed in the fourth position. The remaining modes represent the edgewise modes. In
Figure 9, it is possible to notice the absence of the torsional mode. As already stated, when nonlinear effects are not present, edgewise displacement and torsion are not coupled, thus exciting the edgewise displacement; the torsional mode cannot be observed. On the other hand, the numerical modes predict five uncoupled modes as expected.
For
L = 1500 mm, the experimental Auto MAC matrix reported in
Figure 10A reveals a certain level of coupling between the torsional mode and the second and third edgewise modes. This coupling is not detected by the linear FE models, while it is present in the nonlinear FE model MAC matrices. The experimental and nonlinear FE model MAC matrices present some differences in the out of the diagonal values. In this case, the torsional mode presents a small peak in the FRF because the nonlinear effect is present but not very relevant, with
deflection. Moreover, the number of targets is relatively low and can cause some discrepancies. However, it is possible to conclude that for
L = 1500 mm, the nonlinear effect is present and can be predicted with the nonlinear BTCE models, but linear modal analysis can be a reasonable approximation for this level of deflection. The Auto MAC matrices for
L = 2000 mm and
L = 2500 mm are reported in
Figure 11A and
Figure 12A. The pattern given by the experimental results is correctly predicted by the nonlinear models; moreover, these cases highlight the lack of accuracy obtained when linear models are considered for the modal analysis of structures with moderate deformations.
Figure 9,
Figure 10,
Figure 11 and
Figure 12 present also the comparison between the experimental modes and the numerical modes calculated with linear and nonlinear FE models. Ideally, if the numerical modes are coincident with the experimental modes, the MAC matrices should be identical to the Auto MAC experimental matrices. In general, it is possible to affirm that the mode shape predicted with the nonlinear BTCE model is in good accordance with the experimental modes; moreover, they are confirmed by the results of the SHELL FE model solved with NASTRAN SOL106. In the first case (
L = 1000 mm), the torsional mode was not detected, and for this reason, the comparison with numerical counterpart is not reported in
Figure 9. The MAC matrices for the beam with
L = 1500 mm reveal a high similarity with the experimental results when nonlinear modal analysis is used, while the similarity is considerably lower when the nonlinear effects are not considered. This is even more evident for
L = 2000 mm and
L = 2500 mm. The fourth case presents a relatively low similarity for the fourth mode (
Figure 12), which corresponds to the third edgewise mode coupled with the torsional mode. This is probably connected to the resolution obtained with the chosen number of targets and can be improved by considering more acquisition points. However, the objective was to keep the mass of the targets negligible for all the cases considered, and for this reason, the number of acquisition points were kept constant throughout the experimental activity.
5.3. Numerical Modal Analysis Results
The results of the numerical modal analysis of the composite box beam structure described by [
37] are represented in
Figure 13,
Figure 14,
Figure 15 and
Figure 16. Eight load cases were considered for a maximum deflection
; six of them correspond to a deflection below
and can help to observe more precisely at which point the nonlinear effects cause the deviation from the linear results. The results of the simulation performed with the BTCE models were compared to the frequencies obtained with a SHELL model solved with NASTRAN SOL106. For this comparison, the first eight modes were investigated. In this case, the comparison was performed on the frequencies computed with the nonlinear models normalized with their linear counterparts computed for the undeformed configuration. With this method, the variation of the characteristic frequency is highlighted. The material orientation causes the flapwise bending–torsion coupling, while the deflection causes the edgewise bending–torsion coupling; for this reason, all the modes involve three degrees of freedom. However, one component of the eigenvector has a considerably higher value than the other; for this reason, the modes where flapwise bending is the major effect are denoted with the letter F, while the modes where the edgewise bending component is predominant are denoted with the letter E, and the mainly torsional modes are denoted with the letter T, as done previously.
The results shows a good correlation between the BTCE models and the SHELL FE model. The first, second, fifth and seventh modes present very similar results, even for large displacements. The third, fourth, sixth and eighth modes present some discrepancies when the deformations are bigger than . A less accurate prediction of the characteristic frequencies can be attributed to many factors. First of all, the number and the nature of the finite element used bring approximations that are necessary to lower the computational costs but can influence the results. Secondly, the hypothesis of inextensibility adopted for the BTCE could be unverified for large nonlinear deformations.Moreover, the curvatures and the rotation matrix are obtained under the hypothesis of moderate-to-large displacements. This comparison shows that the BTCE models could be used for the nonlinear analysis of predeformed structures with deflection below with results comparable to the characteristic modes of a SHELL FE model of the same structure solved with NASTRAN SOL106. Moreover, the results show that the differences between BTCE-GE and BTCE-NL are minimal, up to a deflection of and increase for larger deformations. The models here presented can be further improved with an experimental test involving coupled structures to assess the performance and correctly evaluate the influence of geometrical and stiffness effects.
6. Conclusions
This research work presented the derivation of two models for the dynamic analysis of beam structures with bending–torsion coupling in the presence of geometric nonlinearities. The first model accounts for the geometric effect through the orientation of the beam finite element according to a known equilibrium deformation. The second model accounts for nonlinear terms in the stiffness matrix derivation through the hypothesis of small perturbations of an equilibrium configuration under a static load. The stiffness matrix was derived with Hamilton’s Principle. An experimental activity was carried out with the scope of verifying the level of deflection sufficient to have appreciable nonlinear effects and to assess the accuracy of the nonlinear analysis with the BTCE models. The experimental tests have been performed with a LDV system on an aluminum beam constrained at four different length, this allowed to study the nonlinear effects with four different level of deformation. The experiment showed that for the cases considered, the geometric nonlinearities have minor effects on the characteristic frequencies of the structure. Linear and nonlinear numerical models predicted frequencies generally within an error lower than . Concerning the mode shapes, this research work revealed that for a deformation below , the mode shapes present a low level of coupling, and linear numerical models can be used to study structure under these conditions. On the other hand, for the cases where the deformation was or , the presence of nonlinear couplings determined relevant differences in the mode shapes, which were correctly predicted by the BTCE models here derived. Some of the results presented small differences between the observed modes and the ones predicted with the derived numerical models. These minor discrepancies are connected to the relatively low number of scanning point, which lowered the resolution. Moreover, the BTCE models rely on the equilibrium solution computed with a linear static analysis, which can be less accurate for higher deformations. The experimental activity showed that the stiffness effect plays a minor role for the analysis considered and that the BTCE-GE model can be sufficient for characteristic modes and frequencies prediction. The BTCE models were also tested for a composite structure with bending–torsion coupling. A numerical comparison revealed good accordance with the results obtained with a SHELL FE model solved with NASTRAN SOL106. The use of the presented model can be extended to the study of the aeroelastic performance of wing structures. Moreover, the bending–torsion coupling formulation allows to perform optimization on the material orientation to achieve desired dynamic properties also in the presence of geometrical nonlinearities.