A Substructure Synthesis Method with Nonlinear ROM Including Geometric Nonlinearities

: Large ﬂexible aircraft are often accompanied by large deformations during ﬂight leading to obvious geometric nonlinearities in response. Geometric nonlinear dynamic response simulations based on full-order models often carry unbearable computing burden. Meanwhile, geometric nonlinearities are caused by large ﬂexible wings in most cases and the deformation of fuselages is small. Analyzing the whole aircraft as a nonlinear structure will greatly increase the analysis complexity and cost. The analysis of complicated aircraft structures can be more efﬁcient and simpliﬁed if subcomponents can be divided and treated. This paper aims to develop a hybrid interface substructure synthesis method by expanding the nonlinear reduced-order model (ROM) with the implicit condensation and expansion (ICE) approach, to estimate the dynamic transient response for aircraft structures including geometric nonlinearities. A small number of linear modes are used to construct a nonlinear ROM for substructures with large deformation, and linear substructures with small deformation can also be assembled comprehensively. The method proposed is compatible with ﬁnite element method (FEM), allowing for realistic engineering model analysis. Numerical examples with large ﬂexible aircraft models are calculated to validate the accuracy and efﬁciency of this method contrasted with nonlinear FEM.


Introduction
Geometric nonlinearities must be considered in the design process of modern, highperformance aircraft, especially high-altitude long-endurance (HALE) aircrafts. Because of light-weight design and substantial flexibility, the wings of HALE aircraft always have high slenderness and produce large deformation during flight. Such large deformation leads to vital changes in structural stiffness characteristics and aerodynamic configuration, which obviously affect dynamic response, stability and flight performance of aircraft [1][2][3]. Accurate calculation of structural dynamic response has great significance in the design of high-performance aircrafts [4].
Direct numerical simulation approach of geometric nonlinear dynamic response with physical degrees of freedom (DoFs), such as commercial finite element method (FEM) packages, often holds prohibitive computational burden. Especially for complex aircraft structures, detailed FEM models of structures often tend to possess large numbers of DoFs in order to account for extremely detailed geometric features and material distributions. High computational burden leads to analytical difficulties on aeroelasticity or optimization, etc. Reduced-order model (ROM) is an attractive theory for economical geometric nonlinear dynamic response simulation. Parametric ROM is especially receiving attention, in which the equations of motion are expressed explicitly in the form of ordinary differential equations (ODEs) including polynomial terms of coordinates. These can be directly integrated in the time domain, or be used for further analysis by well-established approaches [5]. distributed nonlinear problems, Apiwattanalunggran extended the fixed interface CMS method by using fixed interface nonlinear normal modes to analyze structures with weak nonlinearities [31]. Karpel added virtual mass elements to the substructures to simulate the effects of other substructures components [32]. They analyzed the static and dynamic responses of large flexible structures. This method divided the wing component into segments rather than as a whole structure in order to connect with the fuselage, so it was more like a combination of the nonlinear CMS method and the finite segment method. Kantor extended this method and applied it into a simple large flexible aircraft [33]. With strip theory, geometrical nonlinear aeroelastic analysis was implemented. Joannin and Laxalde extended the classic CMS method by nonlinear complex modes for the study of nonlinear and mistuned cyclic structures [34,35]. Kuether presented two modal substructure methods with nonlinear normal modes and a nonlinear reduced-order model under linearization [36]. He also studied a nonlinear substructure method combined ICE method and C-B method with a simple beam model and gave the comparison of calculation results based on back-bone curves [37,38]. Wu et al. developed interface reduction techniques in the C-B method considering geometric nonlinearities and modal derivatives were used for nonlinear displacement recovery [39]. However, existing substructure methods for structure modeling with geometric nonlinear subcomponents is often independently based on the fixed interface substructure techniques. Under the analysis requirement of complex structures, especially aircraft structures, key points of each component, such as wing and fuselage modeling, are different. Development of a hybrid interface substructure technique considering geometric nonlinearities is essential.
This paper aims to develop a nonlinear hybrid interface substructure method in order to solve nonlinear dynamic response analysis problems of large flexible aircraft structures. Complex structure models can be divided into nonlinear components and linear components. The nonlinear components are modeled with ROM by ICE method using selected linear modes as reduced basis vectors in modal space. A series of static test cases are applied to generate nonlinear stiffness coefficients by a regression procedure. Nonlinear stiffness terms of ROM are approximated with polynomials. The linear components are modeled naturally by linear modes. Each component is assembled to satisfy force equilibrium and compatibility by CMS techniques. A set of nonlinear equations of motion can be integrated into simulation dynamic transient responses under various inputs. The application of the method described can be employed with a relatively detailed FEM and largely reduces the computational burden.
The content of this paper is arranged as follows. Nonlinear ROM construction for nonlinear components is presented in Section 2.1. Nonlinear substructure method based on hybrid interface substructure method is introduced in Section 2.2. Numerical examples of a large flexible aircraft are given in Section 3. Considering the characteristics of aircraft, wing sections are modeled as nonlinear components and the fuselage section is simulated as a linear component. The results validate the validity and efficiency contrasted with nonlinear FEM solutions. In Section 4, we finally draw useful conclusions.

Nonlinear ROM Method
Construction of ROM considering geometric nonlinearities is illustrated in this section. It can be applied not only to structure modeling with global geometric nonlinearities, but also later for geometric nonlinear component modeling in a substructure method.

Structure Equations
The development of the structural ROM is based on equations derived from a Galerkin approach to solve geometric nonlinear dynamics in a weak form [6]. The equation of motion of the structure may often be given in a dynamic equation as: where tensor S is the second Piola-Kirchhoff (P-K) stress tensor, tensor F is the deformation gradient tensor, ρ 0 is the mass density of the structure and b 0 is the vector of the body force. X denotes the position vector of the structure in the reference configuration and u i denotes the deformed position vector. With the Galerkin approach, a truncated basis of linear modes Φ= [Φ 1 ,Φ 2 , · · · ,Φ n ] are chosen as the basis functions. A quadratic and cubic polynomial form describing the nonlinear relationship and the dynamic equation in terms of the generalized coordinates can be expressed as [18]: ijl p q j q l q p = F i (2) where M ij are the terms of the reduced mass matrix, F i are the modal components of the external force, q 1 , q 2 , · · · , q n are general modal coordinates, E ij , E ijl and E ijl p are the components of the tensors of the reduced stiffness. Einstein summation convention is applied to the formulation.
When a truncated basis of the linear modes is chosen as the basis functions, M ij and ij can be expressed in the formulation as: where M i represents the modal mass term of the i th basis function, and E i is the modal stiffness term of the i th basis function. The formulation of the nonlinear dynamic equations corresponding to the i th basis function can be written as: ijl p q j q l q p = F i The nonlinear ROM equations are obtained in modal space with far fewer degrees of freedom than the full-order FE model. Nonlinear modal restoring forces Q i = E (2) ijl q j q l + E (3) ijl p q j q l q p become functions of general modal coordinates q = (q 1 , q 2 , · · · , q n ).
In another way, the equations of motion can be discretized by FEM in physical space: where M is the linear mass matrix and K is the linear stiffness matrix, f NL (u) indicates the nonlinear restoring force function, u and F are displacement vector and external force vector in physical space. The geometric nonlinear theory demonstrates that the quadratic and cubic polynomials of displacement are accurate enough to describe the nonlinear restoring force [40]. Introduced linear mode Φ and reduced system by: The low-order structure dynamic equations of Equation (6) are in the same formulation as Equation (5). When given specified external forces, the time-dependent modal coordinates can be calculated from Equation (5), and physical displacement vector u is solved by Equation (7). ijl p , the static formulation of Equation (5) is: Evidently, if there is a set of specified static loads and corresponding structural deformations, the unknown nonlinear stiffness terms related to the applied loads and the structural displacement resultant can be found by using regression analysis. Corresponding deformation can be calculated by commercial FEM software packages, such as MSC.Nastran [12,13].
Considering that there are NT sets of static test load cases and corresponding deformations, the loads and deformations are transformed into modal space. To solve the unknown nonlinear modal stiffness terms, the left side of Equation (7) can be fitted in a regression progress. The regression problem can be presented as: It should be noted that the superscript without bracket denotes the serial number of test cases from 1 to NT instead of the power value.

Strategy for Generating Test Load Cases
Regression analysis is performed using the commercial software package on the actual deformation and load testing after FEM analysis. The accuracy of the nonlinear stiffness coefficients depends directly on the rationality of the selected test load cases, which is related to the success of recovery of the nonlinear structure equations. Two factors are important in generation of test load cases. The first is the spatial distribution of the load over the FEM model. The second is the overall magnitude of the FEM loading [12]. With the absence of a priori information about the exact nature of the nonlinear stiffness coefficients, a sum of a number of weighted mode shapes are often used as test load cases in ROM [12,13,16]. For a given test load cases: Here, F is a vector of discrete loads in physical nodal space, a r are scalar weighting factors for each mode shape and NR is the number of modes selected for test load case generation. The selected test cases cover all the possible nonlinear characteristics of the structure investigated. The desired displacement at the points of interest can be reached by selecting appropriate weighting factors a r . The scalar weighting factors a r can be obtained as: where W C is the desired transverse displacement at location "C" in the structural model, r is the transverse displacement of Φ r at location "C".

Nonlinear Substructure Method with ROM
To simulate the dynamic response of complex structure by using the substructure synthesis method, the whole structure should be divided into a limited number of substructures. The substructure synthesis method extended to consider geometric nonlinearities is illustrated in this section. The structure can be divided into nonlinear substructure components with fixed interface and linear substructure components with free interface. Nonlinear components can be modeled by ROM proposed before.

Nonlinear Substructure Component
Considering a nonlinear substructure (nl = 1, 2, · · · , p) system, each dynamic equation of substructures in physics coordinates can be expressed as: Subscript t, b represents interior and boundary coordinates, g(u) is the nonlinear function, G b is the boundary load vector and f is the external force vector.
It should be noted that, only the geometrical nonlinear problem is considered here, and the nonlinear problem of interface connection of aircraft structure is not considered. The applied ICE method, g t (u) and g b (u) can be obtained from regression analysis described before. When given u (nl) b = 0 as fixed interface boundary condition, the structure dynamic equation of interior freedoms is: Using linear modes to reduce order with u t (nl) =Φ t (nl) q t (nl) , the low order equation is: 1jl p q j q l q p , E 2jl q j q l +E (3) 2jl p q j q l q p , · · · , E (2) njl p q j q l q p ) (nl)T (17) where M i , E i and Φ t represent the modal mass term, modal stiffness term and mode shape matrix corresponding to substructure with fixed interface. Static constraint modes (SCMs) Ψ b are introduced to translate displacement: The structure dynamic equation can be reformulated as:

Linear Substructure Component
Considering a linear substructure (l = 1, 2, · · · , q) system, each substructure dynamic equation in physics coordinates can be expressed as: According to the basic assumption in the free interface CMS method [25], the displacement of each substructure can be represented by combination of linear modes Φ (l) under free-interface boundary and residual modes Ψ d (l) as follows: Modal coordinates q (l) can represent all freedom, and not only interior freedom, because of the free interface boundary condition. The structure dynamic equation can be reformulated as:

Substructure Synthesis Techniques
For simplicity, it is assumed that the whole structure consists of two nonlinear substructures (Sub-A and Sub-B) components and one linear substructure (Sub-C) component. The connected relationship is shown in Figure 1.  Three components dynamic equations can be presented based on Equations (19) and (21) as: Three components dynamic equations can be presented based on Equations (19) and (21) as: The compatibility equations of the interface coordinates and forces can be written as follows: Here are Bohr matrices for determining position of the connection freedoms between substructures. Substituting Equations (28)-(31) into Equations (25)-(27) and moving nonlinear terms to the right side of equations, coupling of these equations can be given as: Given the specified external force, the modal coordinate response can be obtained from Equation (32) and the physical coordinate response can be solved by Equations (18) and (22). In this manner, the dynamic transient response of a large flexible structure with substructures is successfully simulated using the proposed method.

Numerical Examples
A large flexible aircraft model is presented here to demonstrate the accuracy and efficiency of the proposed method for nonlinear dynamic transient analysis. Static analysis with lateral loads on wing-tip and dynamic transient analysis with harmonic lateral loads on wing-tip are implemented. The response analysis results are compared to nonlinear FEM.

FEM Model
A large flexible aircraft model with flying wing configuration is presented in Figure 2. Main design parameters are given in Table 1  Low order linear mode analysis results of the aircraft model are presented in Figure 3 and Table 2. The first order of elastic mode is 1st symmetric vertical bending mode and the frequency is 3.67 Hz which indicated that the flexibility of the model is high.  Low order linear mode analysis results of the aircraft model are presented in Figure 3 and Table 2. The first order of elastic mode is 1st symmetric vertical bending mode and the frequency is 3.67 Hz which indicated that the flexibility of the model is high.
Width of fuselage/mm 1200.0 Weight/kg 20.0 Low order linear mode analysis results of the aircraft model are presented in Figure 3 and Table 2. The first order of elastic mode is 1st symmetric vertical bending mode and the frequency is 3.67 Hz which indicated that the flexibility of the model is high.   Figure 4. Low order linear modes analysis results of substructures are given in Figures 5 and 6 and Tables 3 and 4. The boundary condition of the fuselage substructure model is free condition and the boundary condition of the wing substructure models is root fixed condition. It is consistent with the synthesis method described before. Due to the symmetricity of the   Figure 4.
Low order linear modes analysis results of substructures are given in Figures 5 and 6 and Tables 3 and 4. The boundary condition of the fuselage substructure model is free condition and the boundary condition of the wing substructure models is root fixed condition. It is consistent with the synthesis method described before. Due to the symmetricity of the aircraft structure, only the right wing substructure analysis results are shown.    Figure 4. Low order linear modes analysis results of substructures are given in Figures 5 and 6 and Tables 3 and 4. The boundary condition of the fuselage substructure model is free condition and the boundary condition of the wing substructure models is root fixed condition. It is consistent with the synthesis method described before. Due to the symmetricity of the aircraft structure, only the right wing substructure analysis results are shown.

ROM Model of Wing Components
The ICE technique is used for ROM conduction of the wing component including geometric nonlinearities. The accuracy of the nonlinear component model should be validated before complex structure modeling. The sample datasets used in these numerical examples are generated in MSC.Nastran. The first vertical bending mode and the first torsion mode are selected for generating test load cases with Equation (10). Applied load on wing model structure, corresponding deformation in physical space can be calculated directly. Choosing a series of appropriate scalar weighting factors to satisfy demand of large deformation situations, the test datasets used include 216 samples and the distribution of vertical deflection of the wing tip by sample numbers is shown in Figure 7.

ROM Model of Wing Components
The ICE technique is used for ROM conduction of the wing component including geometric nonlinearities. The accuracy of the nonlinear component model should be validated before complex structure modeling. The sample datasets used in these numerical examples are generated in MSC.Nastran. The first vertical bending mode and the first torsion mode are selected for generating test load cases with Equation (10). Applied load on wing model structure, corresponding deformation in physical space can be calculated directly. Choosing a series of appropriate scalar weighting factors to satisfy demand of large deformation situations, the test datasets used include 216 samples and the distribution of vertical deflection of the wing tip by sample numbers is shown in Figure 7. The maximum wingtip deflection in the samples is about 20% of the span length of the wing model. It satisfies the requirement of nonlinear analysis of aircraft structures. The ROM model is validated by predicting the static deformation of the wing model. In this part, 70% (151 samples) of the 216 samples are used for training the ROM model and the remaining 30% (65 samples) are used for validation of the model performance.
The least-square technique is used as regression method for finding nonlinear stiffness coefficients. Figure 8

Response Analysis Results of Aircraft
With the ROM of nonlinear wing components established, the aircraft structure model can be assembled. Nonlinear static and dynamic responses are simulated for validations. In nonlinear static response analysis, lateral loads along the z direction are applied on the tip of the left and right wings. The aircraft model is fixed at the center of gravity. The response analysis results of wing tip deflection are given in Figure 9, using the nonlinear substructure method proposed, nonlinear FEM (full order method, Nastran SOL106) and linear analysis method (full order, Nastran SOL101), up to the load level. It can be seen that the analysis results of the nonlinear substructure method can reach good agreement with the nonlinear FEM. The nonlinear substructure method has obvious advantages in computational cost compared to the full order analysis method. In comparison, the deflection results of the linear analysis are obviously higher than the nonlinear analysis method and the difference gradually increases with increasing load. The linear analysis method is not suitable for modeling the large flexible structures.
In nonlinear transient dynamic response analysis, a harmonic lateral load is applied on the tip of the left and right wings. The boundary condition is the same as static response analysis. The harmonic load F is given as: where F 0 is the load amplitude and f is the load frequency. The time step of simulation is set as 0.001 s and the total time is 3 s. Nonlinear FEM (full order method, Nastran SOL400) and linear analysis method (full order, Nastran SOL109) are also implemented with the same time step parameters for contrast. Figures 10 and 11 present dynamic transient analysis results of wing tip with load amplitude 7 N and 20 N. The load frequencies are both 2 Hz. With lower load amplitude, the linear and nonlinear analysis results show no obvious difference. The analysis results using the nonlinear substructure method proposed match well with results from nonlinear FEM.
When the load amplitude increases to 20 N, the geometric nonlinearities lead to an important impact in response analysis. The linear analysis and nonlinear analysis results are obviously different. The analysis results of nonlinear substructure method can still reach a good agreement with nonlinear FEM. In the case of large deformation, the method proposed in this paper shows high accuracy for dynamic response simulation. Figures 12 and 13 give the dynamic transient response analysis results with a large deformation when changing the load frequencies to 1 Hz and 3 Hz. The load amplitude remains 20 N. Load frequencies are lower than the first order elastic mode frequency and the response deformations increase with increasing load frequency. Figure 14 gives the dynamic transient response analysis results with the load frequency changing to 7 Hz. Load frequency is much higher than the first order elastic mode frequency and the response deformation is small. The analysis results from the nonlinear substructure method still maintain good agreement with nonlinear FEM. It is worth noting that the difference between two nonlinear analysis methods does not change significantly as the response deformation increases. With load frequencies of 1 Hz and 3 Hz, the linear analysis method also leads to obviously different response results when considering geometric nonlinearities as before.  SOL106) and linear analysis method (full order, Nastran SOL101), up to the load level. It can be seen that the analysis results of the nonlinear substructure method can reach good agreement with the nonlinear FEM. The nonlinear substructure method has obvious advantages in computational cost compared to the full order analysis method. In comparison, the deflection results of the linear analysis are obviously higher than the nonlinear analysis method and the difference gradually increases with increasing load. The linear analysis method is not suitable for modeling the large flexible structures. In nonlinear transient dynamic response analysis, a harmonic lateral load is applied on the tip of the left and right wings. The boundary condition is the same as static response analysis. The harmonic load F is given as: where 0 F is the load amplitude and f is the load frequency. The time step of simulation is set as 0.001 s and the total time is 3 s. Nonlinear FEM (full order method, Nastran SOL400) and linear analysis method (full order, Nastran SOL109) are also implemented with the same time step parameters for contrast.    High computing efficiency is an important value for low order models. The computational costs of different methods are presented in Table 5. The nonlinear substructure method spends much less time compared to nonlinear FEM with the same computing High computing efficiency is an important value for low order models. The computational costs of different methods are presented in Table 5. The nonlinear substructure method spends much less time compared to nonlinear FEM with the same computing hardware (3.40 Ghz, Intel core i7-3770, 8.00 GB RAM). It will be a huge advantage of the application of this method in optimization or aeroelasticity, etc. which needs massive iterative processing.

Conclusions
In this paper, a novel nonlinear hybrid interface substructure method with ROM has been developed for dynamic transient response analysis for aircraft structure models with large flexibility. Nonlinear ROM based on ICE techniques is an efficient method for large flexible structure modeling. However, for complex aircraft models, only wing components have large deformation. The analysis can be more efficient and practical if the substructure method can be applied. The method proposed divides the structure into nonlinear substructure components and linear substructure components. The nonlinear substructure components are modeled with ICE techniques and the nonlinear stiffness coefficients are approximated with polynomials. Regression analysis is used for finding the nonlinear stiffness coefficients. The linear components are modeled naturally by linear modes. Each component is assembled by CMS techniques and the boundary conditions of force equilibrium and displacement compatibility are satisfied.
Numerical examples with large flexible aircraft models are presented to validate the accuracy and efficiency compared to nonlinear FEM. Under static loads and dynamic harmonic loads, the analysis results reveal high consistency with nonlinear FEM in both small and large deformation conditions. Meanwhile, the computational time is much lower than nonlinear FEM. The nonlinear substructure method proposed in this paper has high accuracy and efficiency at the same time. It can be used well for large flexible structure response analysis.