Model order reduction for structural nonlinear dynamic analysis based on Isogeometric analysis

Model order reduction approach generates lower dimensional approximations to the original system while preserving model’s essential information and computational accuracy. For nonlinear structural dynamic problems, where the stiffness matrix is configuration dependent, an iterative solution procedure is inevitable and a revisit to all the elements is essential for updating the stiffness matrix. In this paper, the nonlinear dynamics of the planar curved beams and 3D cylindrical shells are studied based on the isogeometric analysis and their model order reductions are investigated based on the proper orthogonal decomposition and discrete empirical interpolation method (POD-DEIM). Numerical results show that IGA-based POD-DEIM method significantly improves the computational efficiency of the nonlinear dynamic analysis of the beam and shell structures.


Introduction
Model order reduction approximates the full order system with reduced dimensions while maintaining most of the primary properties of the original system. Proper orthogonal decomposition (POD) is widely used in the reduced order modeling, which generates base vectors and approximates the full order system in a low dimension space. For nonlinear problems where the stiffness matrix is configuration dependent, it is inevitable to revisit all the elements to update the stiffness matrix in each iteration which is not desired for the reduced order modeling [1] . In view of the deficiency of POD, Barrault [2] and Grepl et al. [3] proposed an Empirical Interpolation Method (Empirical Interpolation Method, EIM), it uses a linear combination of base vectors to approximate nonlinear terms where the combination coefficients can be obtained by interpolation. Besides, the base vectors for the ROM are selected based on the greedy algorithm. However, EIM can not be extended to arbitrary nonlinear functions. Chaturanabut et al. [1] proposed a discrete form of EIM, which suits arbitrary ordinary differential equations and is mathematically equivalent to EIM. It uses the selected POD base vectors as input, and DEIM to interpolate the tangent stiffness matrix during the POD model reduction process. Successful example can be found in computational aeroelasticity [4] and structural dynamics [5] . Compared with our previous work [6] , this paper expands the reduced order model of nonlinear structural dynamic analysis of Euler-Bernoulli planar curved beams to that of Kirchhoff-Love shells.
Isogeometric analysis (IGA) proposed by Hughes [7] in 2005, uses Non-Uniform Rational B-Splines (NURBS) to describe the geometry and interpolate the physical field, thus unifies CAD and CAE in a integrated sense. IGA has the advantages of accuracy and high-order continuity which is particularly suitable for the analysis and optimization of thin-wall structures. This paper analyzes the nonlinear

Isogeometric structural analysis
For geometrical nonlinear elastic problems, the governing equations of Euler-Bernoulli planar curved beam and Kirchhoff-Love shell is written as: where n is the membrane force, m is the bending moment, ε is the membrane strain, κ is the curvature, u is the second derivative of displacement u w.r.t. time t, A is the analysis domain,  is the density, dA is the differential area in the reference configuration, 0 T and 0 P are the traction along the Neumann boundary  and the unit surface load, respectively, which are assumed to be dead loads.
According to the concept of IGA, the displacement u and its variation  u can be discretized by the same NURBS basis as the geometry descriptions: where i R is the corresponding NURBS basis, i u and i u  represents sets of three unknown values of each control point, R and u are the metrices of the corresponding quantities. Substituting (3) and (2) into (1), the semi-discretized equilibrium equation is obtained as: Where M is the mass matrix, :

Concept of model order reduction
Neglecting the effect of damping, the discrete form of the system equilibrium equation can be expressed as: where nn   K is the stiffness matrix,  (7), we can obtain: . Multiplying (9) to the left by T V , one can obtain the reduced form of the governing equation:

Proper orthogonal decomposition
The set of orthonormal base vectors V are computed by the Proper Othogonal Decomposition (POD) method. In full order model analysis, the snapshots of the displacement and stiffness matrices are collected and assembled to the matrices If k base vectors are selected to construct the reduced order model, its approximation error to the snapshots in the sense of least square can be expressed as: Then we can choose the number of base vectors according to the following criterial: represents the relative error of the selected base vectors approximating the original system represented by the captured snapshots.

Discrete Empirical Interpolation Method
To further accelerate the computational speed, discrete empirical interpolation method is used to interpolate the stiffness matrix with nonzero-values at selected locations. Suppose that the stiffness column vector Submitting c of (15) into (14) yields the approximated stiffness vectorˆi K :

Numerical examples
In this section, we consider two examples: the Archimedes spiral curved beam and the 3D cylindrical shell. In the HHT- method [10] , a specific set of parameters (   Table 1), which is mainly due to the large support size of the NURBS basis. For comparison, Abaqus solutions of both linear and geometrically nonlinear elements (B21), together with the nonlinear dynamic responses of the FOM and ROM are shown in Fig 3(a). The displacement difference in the x and y directions in Fig 3(b) demonstrates that the calculation results of the full order model and reduced order model agree well, which indicates the accuracy of the proposed reduced order model.    In this example, a simply supported cylindrical shell subjected to a centrally located concentrated force is studied. The geometric descriptions and material properties are given in Fig 1(b). The shell is discretized with 267dofs and 60 elements, with 10 and 6 elements in the circumferential and axial directions, respectively. The polynomial orders in two parametric directions are chosen as p=q=3.

Cylindrical shell
In the offline stage, the time step size is set to be 0.01. Then, the snapshots na In the online stage, the same time span and time step size are kept as the offline stage. The number of base vectors (i_u=22, i_k=9) of displacement and stiffness matrix are determined by (13) and the value of  . The nonzero items of the stiffness matrix of full order model and the reduced order model are shown intuitively in Fig 4(c) and (d), respectively. It can be seen that only 9 positions are need for the ROM, while the number of involved elements is 46, which is mainly due to the tensor product of NURBS basis has larger support size than single NURBS basis.

Summary
In this paper, model order reduction of structural nonlinear dynamics is studied based on the isogeometric analysis. Proper orthogonal decomposition of the displacement field is used for the dimensional reduction. To deal with the nonlinearities of the stiffness matrix, discrete empirical interpolation method combined with POD is adopted to further reduce the computational burden of the ROM. Based on the proposed method, only a few entries of the stiffness matrix are needed for the reconstruction of the full stiffness matrix which significantly improves the computational efficiency of