Abstract
This paper reports on the improvement of CraigBampton (CB) method for transient analysis of structures with largescale plastic deformation. As is known, the CB method is effective and accurate in reduced order modeling for linear system. In contrast to this, an improved CB method using tangent modes for nonlinear dynamic analysis has been developed. To do this, the incremental governing equations of nonlinear system are linearized in each time step by using tangent stiffness matrix, and the corresponding tangent modes are proved to be orthogonal with respect to mass matrix as well as with respect to tangent stiffness matrix by incorporating the elasticplastic material behavior. Thus, the tangent modes can be used to assemble the transformation matrix of CB method in nonlinear dynamic analysis. Using the proposed method, two elasticplastic beams loaded impulsively are examined. Simulation results show that the improved CB method is valid and accurate for transient analysis of structures with largescale plastic deformation and has lower computational cost compared with full order model.
1. Introduction
The elasticplastic dynamic analysis of structure is typically performed using finite element method (FEM), and has been remarked as a complicated and timeconsuming work [1]. Because the plastic deformation is path dependent, the external load must be applied onto structure incrementally during the analysis and the time step size should be small enough to ensure the accuracy and reliability of solution. If plastic deformation occurs during the time step, further iterations are required to obtain convergent results. For large or complex structure, the computational cost of elasticplastic dynamic analysis is substantially high. In order to mitigate this problem, model order reduction technique is performed [23].
Modal superposition is a classical method to reduce the model order for liner system. It can offer fairly accurate solution when only a few modes are considered. The first application of modal superposition method into nonlinear system is most probably from Nickell [4] in 1970s. Subsequently, the modal superposition method has been applied for various nonlinear problems. In the literatures, there are basically two strategies to apply this method into nonlinear problems [57]. In the first one, the governing equations of nonlinear problem are written in an incremental form, and are linearized over small time interval by using a tangent stiffness matrix. The eigenvalues and eigenvectors of the nonlinear system are updated at every time step. The works employed this method were proposed by Nickell [4], Morris [8] and Mohraz et al. [9], as well as some local nonlinear and entire nonlinear examples are examined to demonstrate the validity and accuracy of the proposed procedure.
In the second strategy, the nonlinearities of system are approximated as an internal excited load, which is placed on the righthand side of the governing equations. This technique is the socalled pseudoforce method. Thus, the eigenvalues and eigenvectors of the initial state (elastic state) of structure are used throughout the entire analysis. Some nonlinear analyses with such approach were studied and introduced by Bathe [10], Kukreti [11], Muscolino [12], Villaverde [13], and Manoach et al. [14, 15].
In general, the first strategy of applying modal superposition method is accurate, but requires a large computational effort for solving the instantaneous eigenvalue problem in each time step, and the efficiency of this method is questionable in practical analysis, especially for large problems. The advantage of the second approach is that the method avoids the recalculation of the eigenvalue problem in each time step and can be computationally more efficient than a direct integration method. However, this approach is probably best employed in the analysis of only slightly nonlinear systems or systems with only local nonlinearities [6, 7, 10].
Another commonly used reduction method in linear system is component mode synthesis (CMS) technique, which is originally developed for modal analysis of large or complex structure. With this method, first the entire structure is partitioned into a number of substructures, a set of structure modes (e.g. exact eigenmodes, static modes, interface modes, etc.) called component modes were set up to represent the motion of substructures. These component modes then were assembled to get a reduced space. Finally, the entire system is reduced by projecting the governing equations into the reduced modal coordinate space.
In 1960s, Hurty [16] first proposed the CMS technique. In his idea, three types of component modes were used: rigid body modes, constraint or attachment modes, and fixedinterface normal modes. Craig and Bampton [17], subsequently proved that it was not necessary to separate the constraint modes and rigid body modes, and suggested using only constraint modes and fixedinterface normal modes to generate the reduced modal space. The proposed method is referred as CraigBampton (CB) method. After that, various CMS methods have been developed [1822]. However, the CB method is still the most popular and widely used CMS method in structural dynamics because of its simplicity and reliability. Some reviews of CMS method can be found in [2325].
Although, the CB method has been widely used in structural dynamics, the applications for nonlinear system can be cited here are only a small number, and the literatures are mainly concentrated on local nonlinear or weak nonlinear problems. As nonlinear systems do not obey the principle of superposition, the applications of CB method in nonlinear system require certain techniques.
The solutions of contact impact problems with CB method were firstly presented by Wu and Haug [26]. They avoided the contact nonlinearities by using constrain additiondeletion technique. The contact force is addressed as internal force of system. And more deep studies about contact impact problem then were developed [2729] with this idea, including the elastic stress wave propagation, “succession collisions”, and dispersion of stress waves in rod.
In many dynamic problems, the nonlinearity can be limited in a few regions of structure whereas the rest of structure remains linear. For this type of problems, the CB method is typically applied in conjunction with substructure technique. The entire structure is first partitioned into linear component and the other component with nonlinearity, then the linear component is modeled using CB method while the other component is modeled using physical coordinate, these components are finally coupled to dominate the motion of entire structure. This idea has been applied for some practical problems [3032].
If the nature of nonlinearity in systems is weak, it is possible to approximate the nonlinear function by a simpler form, generally a combination of some linear functions. With the help of this idea, the CB method can be applied to certain weak nonlinear problems. Kuether and Allen [33] have calculated the response of a geometrical nonlinear beam with CB method. The nonlinearities are modeled using polynomials in modal coordinate and a series of unite static forces was applied to the full order model to interpolate the coefficients of the polynomial. Bond and Khraishi [34] have proposed an iterative algorithm to approximate the nonlinear material effect as an internal modal force which is placed on the right hand side of governing equations, and used CB method in conjunction with rigid body solution to extend the useful range of rigid body design tools.
Here, we focus on the reduction of transient analysis of structure with largescale plastic deformation. Which is a strong nonlinear dynamic problem, do not obey the principle of superposition, it is not possible to use the original CB method. In order to overcome this difficult, we learn from Nickell’s idea, which was originally proposed to explore the use of mode superposition method for nonlinear problems [4]. That is, the forced motion of nonlinear system over a small time interval can be represented with the nonlinear (tangent stiffness) frequency spectrum. With the help of this idea, the nonlinear frequency spectrum is used to generate the reduced modal space, and the original CB method is improved to be able to reduce the transient analysis of structure with largescale plastic deformation.
The process of this work can be summarized as follows: First, the incremental governing equations of elasticplastic dynamic problem are linearized over a small time interval by using tangent stiffness matrix. The tangent modes of structure based on the nonlinear state at the beginning of the time interval are then proved to be able to control the forced motion of structure, by incorporating the elasticplastic material behavior. Finally, the reduced modal space of CB method is constructed with the tangent modes of substructures.
In the next section, the incremental governing equations of elasticplastic dynamic problem and corresponding solution procedure are described. Following that, the formulations of improved CB method are presented, including the discussion of orthogonality of tangent modes for nonlinear system, basic procedure of the proposed method and some tips to improve the efficiency. Then, some numerical examples are provided to demonstrate the performance of the proposed method in accuracy, validity and efficiency. At the last, the paper is concluded.
2. Equations of elasticplastic dynamic problem
Considering the small strain and deformation type of a motion of an elasticplastic finite body $\mathrm{\Omega}$, by use of the nonlinear finite element theory, the governing equations for the general elasticplastic structure during a small time interval ($t\sim t+\mathrm{\Delta}t$) can be expressed as [1]:
where $\mathbf{M}$ is the constant, symmetric, positivedefinite mass matrix. ${\mathbf{K}}_{t}$ is the tangent stiffness matrix at time $t$, is dependent on the nonlinear constitutive relation of material. ${\mathbf{P}}_{t+\mathrm{\Delta}t}$ is the excited force vector at time $t+\mathrm{\Delta}t$. $\mathbf{F}\left({\mathbf{u}}_{t}\right)$ is the resisting force vector of structure at time $t$, depends on the displacement field. The ${\ddot{\mathbf{u}}}_{t}\text{,}$${\mathbf{u}}_{t}\text{,}$$\mathrm{\Delta}\ddot{\mathbf{u}}$ and $\mathrm{\Delta}\mathbf{u}$ are the nodal acceleration vector, nodal displacement vector, nodal increment acceleration vector and nodal increment displacement vector, respectively. This equations, are the completely incremental form for elasticplastic dynamic problem, and use the “residual force correction” [4]. The damping can be included in Eq. (1), but is not relevant to the following discussion, so it is neglected.
The matrices in Eq. (1) are the sums of the corresponding quantities of every element, and are expressed as:
where ${V}_{e}$ and ${S}_{e}$ denote the volume and surface of element, respectively. $\mathbf{G}$ is a Boolean matrix that represents the transformation from the element coordinate to the physical coordinate. $\mathbf{N}$ is the matrix of shape function of element, $\mathbf{B}$ is straindisplacement matrix of element. $\rho $ is the mass density, ${\mathbf{D}}_{ep}\left(t\right)$ is the elasticplastic matrix of material and ${\mathbf{T}}^{e}$ is the boundary surface traction vector of an element. The elasticplastic matrix ${\mathbf{D}}_{ep}\left(t\right)$ has explicit description in terms of yield function and hardening model, which can be found for more details elsewhere in [1].
The Eq. (1) is usually solved by a stepbystep time integration method, and iterations are needed in the time step once the plastic deformation occurs. There are two basic iterative techniques for determining the plastic deformation in the iterations, which are the NR (NewtonRaphson) and MNR (Modified Newton Raphson) method. Here, we employed the Newmark$\beta $ method in conjunction with MNR iterative technique for solving Eq. (1). By doing so, the Eq. (1) is linearized between time steps. The tangent stiffness matrix of structure at the beginning of each time step is used and kept unchanged in a time interval. The Newmark$\beta $ method has some assumptions on velocity and displacement, which are:
where $\alpha $ and $\beta $ are integration parameter. Considering the relation of $\mathrm{\Delta}\ddot{\mathbf{u}}={\ddot{\mathbf{u}}}_{t+\mathrm{\Delta}t}{\ddot{\mathbf{u}}}_{t}\text{,}$$\mathrm{\Delta}\dot{\mathbf{u}}={\dot{\mathbf{u}}}_{t+\mathrm{\Delta}t}{\dot{\mathbf{u}}}_{t}$, $\mathrm{\Delta}\mathbf{u}={\mathbf{u}}_{t+\mathrm{\Delta}t}{\mathbf{u}}_{t+\mathrm{\Delta}t}$, this yields:
Using Eq. (78) into Eq. (1), the corresponding effective static equilibrium equation of the nonlinear system is given by:
For a simple description of Eq. (9):
where ${\widehat{\mathbf{K}}}_{t}$ denotes the effective stiffness matrix, $\mathrm{\Delta}\widehat{\mathbf{Q}}$ stands for the effective force.
In MNR procedure, the iterations are governed by:
The iterations start from: ${\mathbf{u}}_{t+\mathrm{\Delta}t}^{0}={\mathbf{u}}_{t}$, ${\dot{\mathbf{u}}}_{t+\mathrm{\Delta}t}^{0}={\dot{\mathbf{u}}}_{t}$, ${\ddot{\mathbf{u}}}_{t+\mathrm{\Delta}t}^{0}={\ddot{\mathbf{u}}}_{t}$, ${\mathbf{\sigma}}_{t+\mathrm{\Delta}t}^{0}={\mathbf{\sigma}}_{t}$, and are needed until a predefined convergence criterion is achieved. There are many options in convergence checks that are based on the change in displacement, velocity, acceleration, or force. After the convergent incremental displacements are determined, the state variables can be updated as follows: ${\mathbf{u}}_{t+\mathrm{\Delta}t}^{}={\mathbf{u}}_{t}+\mathrm{\Delta}\mathbf{u}$, ${\dot{\mathbf{u}}}_{t+\mathrm{\Delta}t}^{}={\dot{\mathbf{u}}}_{t}+\mathrm{\Delta}\dot{\mathbf{u}}$, ${\ddot{\mathbf{u}}}_{t+\mathrm{\Delta}t}^{}={\ddot{\mathbf{u}}}_{t}+\mathrm{\Delta}\ddot{\mathbf{u}}$.
For solution of Eq. (1), the tangent stiffness matrix is calculated on the basis of the known nonlinear state at the beginning of a time interval ($t\sim t+\mathrm{\Delta}t$) and the incremental displacements are computed with many iterations. In order to reduce the computational effort of the analysis, we consider using CB method.
3. Formulations of improved CB method for elasticplastic dynamic problem
The CB method is a modalbased reduction technique for linear system. Its essence is to generate a transformation matrix to reduce the size of original system. By partitioning the original system into several fixedinterface substructures, the lower normal modes and constrain modes of substructures are assembled to generate the transformation matrix, and then the equations of original system are projected into a reduced space obtained by transformation matrix. During this process, all the mode shapes used are actual modes of substructure, and the normal modes are orthogonal with respect to mass matrix as well as with respect to stiffness matrix of substructure.
By using MNR iterative procedure, the governing equations expressed in Eq. (1) can be linearized during a time interval ($t\sim t+\mathrm{\Delta}t$), it seems involuntary to apply CB method for elasticplastic dynamic problem, directly. In actual, there is a fundamental issue for the application. The mode shapes deduced form Eq. (1) are not actual modes of structure, they are pseudo and be termed as tangent modes. Although, they are considered to be able to express small forced motion of nonlinear system [4]. The application means to utilize the pseudo modes to make reduction for elasticplastic problem. To this end, it is necessary to deeply understand the property of tangent modes (e.g. orthogonality and linear independence) for elasticplastic problem, and demonstrate the feasibility of application from mathematic view.
3.1. Tangent modes and actual modes of structure with plastic deformation
We will first distinguish the tangent modes and actual modes of structure with plastic deformation. The Eq. (1) dominates the motion of structure with plastic deformation in a small time interval ($t\sim t+\mathrm{\Delta}t$). It also constitutes a free vibration problem as the load is zero. The solution of the free vibration can be taken in the usual form:
where $i=\sqrt{1}$ is the imaginary unit, $\omega $ is a real number and $\stackrel{}{\mathbf{x}}$ represents the amplitudes of generalized displacement. The substitution of Eq. (13) into the free vibration yields an eigenvalue problem:
where ${w}_{i}^{2}$ is eigenvalues and ${\mathrm{\Phi}}_{i}$ is the corresponding eigenvectors, they are termed as tangent (instantaneous) modes. Otherwise, the actual modes of structure with plastic deformation are identical to the frequency spectrum of an elastic structure, which is obtained by completely unloading the structure at time $t$. A clearer explanation is shown in Fig. 1. The tangent modes are completely dependent on the state of nonlinear system at time $t$ (as point D), and actual modes are dependent only on the irreversible deformation (as point B).
Fig. 1The difference of tangent modes and actual modes of structure with plastic deformation
3.2. Properties of tangent modes for elasticplastic dynamic problem
In this section, we will further discuss the property of tangent modes by incorporating the elasticplastic material behavior. As was noted in Eq. (1), $\mathbf{M}$ is real symmetric and positive definite for elasticplastic dynamic problem, the solutions of tangent eigenvalue problem in Eq. (14) are sensitive to ${\mathbf{K}}_{t}$ from matrixalgebra theory. Generally, the eigenvalue problem of elasticplastic system can be classified into two conditions. The first, elastic loading or unloading condition, the eigenvalue problem is identical with linear system. Secondly, plastic loading condition, the system performs nonlinear behavior and the solution of eigenvalue problem is controlled by tangent stiffness matrix ${\mathbf{K}}_{t}$. According to Eq. (3), ${\mathbf{K}}_{t}$ is completely dependent on the elasticplastic matrix ${\mathbf{D}}_{ep}\left(t\right)$, which is a function of plastic law and plastic deformation.
The expression of ${\mathbf{D}}_{ep}\left(t\right)$ can be given explicitly in terms of yield function, hardening model and plastic flow rules. In this paper, the study is subjected to the scope of classical plasticity theory. The plastic deformation is volume incompressible, irreversible, and independent of loading rate. We utilize Von Mises yield function to determine whether the plastic deformation is occurring and PrandtlReuss rule to restrict the flow of plastic deformation. The plastic potential function of the material will be same as subsequent yield function, is known as associative plasticity. The subsequent yield function of material is obtained by adopting the isotropic hardening assumption. Thus, the elasticplastic matrix of material is defined as [1]:
where ${\mathbf{D}}_{e}$ denotes the elastic matrix of material, and is real symmetric. $f$ and $\mathbf{\sigma}\left(t\right)$ are the yield function of material and component of stress vector, respectively. $\mathbf{S}\left(t\right)$ is defined as the deviatoric stress vector. $C$ stands for the hardening modulus of the material as a result of the plastic deformation. Herein, it defines an isotropic hardening modulus. For ideal plastic material, $C$ is equal to zero. $H$ is a scalar term.
The elastic matrix ${\mathbf{D}}_{e}$ is real symmetric, the ${\mathbf{D}}_{e}\mathbf{S}\left(t\right){{\mathbf{S}\left(t\right)}^{{\rm T}}\mathbf{D}}_{e}$ and ${\mathbf{D}}_{ep}\left(t\right)$ are real symmetric, correspondingly. We use ${\mathbf{D}}_{ep}\left(t\right)$ into Eq. (3) and can be proved that the ${\mathbf{K}}_{t}$ is also a real symmetric matrix:
When the structure is upon plastic loading, the following conditions will hold true for any point within the structure:
• For hardening material:
• For ideal plastic material:
We consolidate the two formulas and make a simple presentation with matrix form. Let:
According to plastic theory, $d\mathbf{\sigma}\left(t\right)={\mathbf{D}}_{ep}\left(t\right)d\mathbf{\epsilon}\left(t\right)$, the Eq. (19) can be rewritten as:
In an element, the strain increment of any point can be calculated from $d\mathbf{\epsilon}\left(t\right)=\mathbf{B}d{\mathbf{u}}^{e}\left(t\right)$, this yields:
where $d{\mathbf{u}}^{e}\left(t\right)$ is a vector containing the incremental displacements of all DOFs in an element. We integrate Eq. (21) in the region of an element:
$\mathbf{}\mathbf{}\mathbf{}\mathbf{}\mathbf{}\mathbf{}\mathbf{}={\left(d{\mathbf{u}}^{e}\left(t\right)\right)}^{{\rm T}}{\mathbf{k}}_{t}^{e}d{\mathbf{u}}^{e}\left(t\right)\ge 0.$
And assemble all elements within the structure. The result will take the form:
${=\left(\mathbf{G}d\mathbf{u}\right(t\left)\right)}^{T}{\mathbf{K}}_{t}\left(\mathbf{G}d\mathbf{u}\right(t\left)\right)\ge 0,$
where $d\mathbf{u}\left(t\right)$ is the displacement increment vector of structure. Considering the arbitrariness of vector $d\mathbf{u}\left(t\right)$, the ${\mathbf{K}}_{t}$ is proved to be semipositive.
From matrixalgebra theory, it is easy known that all roots of Eq. (14) are real and nonnegative, the corresponding eigenvectors (tangent modes vectors) hold linear independent. More importantly, these eigenvectors can be proved to be orthogonal with respect to mass matrix as well as with respect to the tangent stiffness matrix. Using this property, the incremental governing equations can be decoupled.
3.3. Procedure of the improved CB method
Referencing the application of CB method in linear system, the basic procedure of improved CB method for elasticplastic dynamic problem can be summarized in Fig. 2.
Fig. 2The basic procedure of CB method for elasticplastic dynamic problem: a) full discretized model, b) substructures Ωii=1,2,…,Ns (, where Ns stands for the number of substructures) and boundary Γ, c) fixed interface boundary treatment, d) reduced order model
a)
b)
c)
d)
3.3.1. Modal analysis of substructure
The original structure (discretized model) is first partitioned into ${N}_{s}$ substructures, which are connected with ﬁxed interface at the boundary ($\mathrm{\Gamma}$) between each other. On the basis of Eq. (1), the undamped incremental governing equations of an elasticplastic substructure during the time interval ($t~t+\mathrm{\Delta}t$) can be written in the partitioned form:
where the subscript $i$ indicates the interior degrees of freedom (DOF), the subscript $b$ indicates the boundary DOF. The superscript s denotes the serial number of the substructure.
The CB method uses two sets of substructure modes to generate transformation matrix. The two sets are named as the fixedinterface normal modes ${\mathrm{\Phi}}^{s}\in {\mathbb{R}}^{{n}_{i}\times {n}_{i}}$ and constrain modes ${\mathrm{\Psi}}^{s}\in {\mathbb{R}}^{{n}_{i}\times {n}_{j}}$ (where ${n}_{i}$ and ${n}_{j}$ stand for the number of interior DOFs and the number of boundary DOFs). The fixedinterface normal modes are defined as the free vibration modes of the constrained substructure, whose interfacial DOFs are set to have zero displacement, i.e., ${\mathrm{\Delta}\mathbf{u}}_{b}^{s}=0$. In Eq. (24), this leads to a tangent eigenvalue problem of substructure:
The constraint modes are the static deflections of substructure which are obtained by imposing successive unit displacement at each boundary DOF:
So, the displacement vector ${\mathrm{\Delta}\mathbf{u}}_{}^{s}$ can be projected into a new space defined by ${\mathrm{\Phi}}^{s}$ and ${\mathrm{\Psi}}^{s}$:
We partition ${\mathrm{\Phi}}^{s}$ into dominate modes ${\mathrm{\Phi}}_{d}^{s}$ and high order modes ${\mathrm{\Phi}}_{h}^{s}$. In order to reduce the size of governing equations, the high normal modes will be neglected in Eq. (27). Then, the physical displacement vector ${\mathrm{\Delta}\mathbf{u}}_{}^{s}$ is approximated by the dominated modes:
This is the transformation of displacement vector from physical coordinate to modal coordinate. The vector ${\mathrm{\Delta}\mathbf{a}}_{}^{s}$ is denoted as the modal displacement, its number equals to $k+{n}_{j}$$(k\ll {n}_{i})$, where $k$ is the number of retained fixedinterface normal modes.
Using ${\mathbf{T}}_{CB}^{s}$ in Eq. (24) and left multiplying ${\mathbf{T}}_{CB}^{s}$ on the two sides of the equations, the reduced governing equations of substructure s can be cast in a modal space:
where ${\mathbf{I}}_{d}^{s}\in {\mathbb{R}}^{k\times k}$ is an identity matrix, ${\left({\mathrm{\Lambda}}_{t}^{s}\right)}_{d}\in {\mathbb{R}}^{k\times k}$ is a diagonal matrix, if ${\mathrm{\Phi}}_{d}^{s}$ is normalized. The coupling item ${\stackrel{}{\mathbf{m}}}_{db}^{s}={\left({\mathrm{\Phi}}_{d}^{s}\right)}^{T}{\mathbf{m}}_{ii}^{s}{\mathrm{\Psi}}^{s}+{\left({\mathrm{\Phi}}_{d}^{s}\right)}^{T}{\mathbf{m}}_{ib}^{s}$.
3.3.2. Assembling equations for whole system
Finally, all the substructures are assembled to form the full system expressed by modal coordinates. The process of assembling is similar to finite element analysis. With the consideration of displacement compatibility and force equilibrium of the interface DOFs at the conjunction region of substructures, and eliminating repeated interfacial DOFs, the dynamic equilibrium equations of full system can be expressed as follows:
where $\mathrm{\Phi}$ is the transformation matrix of full system assembled from all substructures. The $\stackrel{}{\mathbf{M}}={\mathrm{\Phi}}^{T}\mathrm{M}\mathrm{\Phi}$, ${\stackrel{}{\mathbf{K}}}_{t}={\mathrm{\Phi}}^{T}{\mathbf{K}}_{t}\mathrm{\Phi}$ and $\mathrm{\Delta}\stackrel{}{\mathbf{Q}}$ are modal mass matrix, modal tangent stiffness matrix and modal external force vector of full system, respectively.
The motion of original system can be dominated by Eq. (30) in a reduced modal space, the modal mass and stiffness matrices are not diagonal, the equations are coupled yet. Converting Eq. (30) to a decoupled system will decrease the condition number of equations and offer substantial computational saving in the solution of the equations. By using eigenvectors of the eigenvalue problem generated from Eq. (30), the coupled equations will become:
where $\mathbf{I}={\mathrm{\Phi}}_{m}^{T}\stackrel{}{\mathbf{M}}{\mathrm{\Phi}}_{m}^{\mathit{}}$ and ${\stackrel{}{\mathrm{\Lambda}}}_{t}={\mathrm{\Phi}}_{m}^{T}{\stackrel{}{\mathbf{K}}}_{t}{\mathrm{\Phi}}_{m}^{\mathit{}}$ are diagonal matrix. ${\mathrm{\Phi}}_{m}^{\mathit{}}$ and $\mathbf{\eta}$ are projection matrix assembled from all eigenvectors and new displacement vector. It can be proved that the projection will not decrease the accuracy of solutions, as it is an equalspace transformation for Eq. (30). And the computational cost of decoupling is sufficiently small attribute to the application of reduced model.
The focus now is turned to the solution of reduced order dynamic equations expressed in Eq. (31). It can be solved by a stepbystep time integration method discussed in Section 2.2. Comparing with Eq. (1), the Eq. (31) has a much smaller order, which is more efficient in numerical computation.
3.4. Some tips to improve efficiency of improved CB method
From a simple perspective, the application CB method is regarded as coordinate transformation of governing equations. When it is adopted for elasticplastic dynamic problem, the operations of generating transformation matrix will be performed at each time step. The basic premise of the application for elasticplastic dynamic analysis is that the transformation matrix is computationally efficient. For this purpose, we will introduce some tips to improve the efficiency of updating transformation matrix.
• Considering the truncation of fixedinterface normal modes, an algorithm was required that would only need to calculate a few lower normal modes of substructure. It seems that a type of iterative method (e.g. subspace iteration method) is a good choose for solving eigenvalue problem [4].
• When solving Eq. (1), one major purpose of updating tangent stiffness matrix is to speed up the convergence of the iterations. In practical analysis, the tangent stiffness matrix can be updated conditionally. An effective strategy for controlling the tangent stiffness matrix updating is based on maintaining the residual norm of the original system with in a prescribed tolerance all along the time history [6]. Therefore, the transformation matrix of CB method for elasticplastic problem can be updated according to residual norm value conditionally. The weighted residual norm is evaluated as [6]:
where ${R}_{t}={P}_{t}F\left({u}_{t}\right){M\ddot{u}}_{t}$ is residual force, $\epsilon $ is the predefined tolerance and $N$ is the DOFs number of original system.
4. Numerical results and discussion
The procedures of the improved CB method and full order model are available in a numerical program written in Microsoft Visual C++ 6.0. Two examples are simulated to demonstrate the validity, and accuracy of CB method for transient analysis with largescale plastic deformation. The integration parameters of the Newmark method with the values of $\alpha =$0.29, $\beta =$0.56, and were used for all solutions. All computations were carried out on an Intel(R)^{}Core(TM) i32350M @ 3.2 GHz and 4 GB RAM personal computer.
4.1. Model description and response evaluation tool
4.1.1. Beam models
The first example is a simply supported beam with perfect plastic material, the beam is subjected to a uniform distributed load $p\left(t\right)$ applied instantaneously and decayed exponentially. The configuration of this example is shown in Fig. 3. The time variation of the exponentially decaying load is given as $p\left(t\right)=n{p}_{0}\mathrm{e}\mathrm{x}\mathrm{p}(2t/T)$. Where ${p}_{0}$ is static collapse load defined as the load causing a plastic hinge to occur at the midpoint of the beam, $n=$ 1.5 and $T$ is the fundamental time period of the beam.
Fig. 3Problem descriptions of simply supported beam, a) beam element model and substructures, b) material constitutive relation, c) load type
a)
b)
c)
Fig. 4Problem descriptions of fixed beam, a) beam element model and substructures, b) material constitutive relation, c) load type
a)
b)
c)
In second example, considering a fixed beam is subjected to a concentrated pulsing load $f\left(t\right)$ at the midpoint of beam, the material of beam is assumed to be bilinear isotropic plastic. Fig. 4 shows the configuration, material constitutive relation and load type of this example. The time variation of the pulsing load can be given from a piecewise function. Table 1 lists the beam characteristics of the two examples.
Table 1Beams characteristics
Example  $L$ (cm)  $W$ (cm)  $H$ (cm)  $\rho $ (kg/m^{3})  $E$ (GPa)  ${E}_{t}$_{}(GPa)  ${\sigma}_{y}$ (MPa)  Material 
Example 1  76.2  2.54  5.08  7800  206.8  –  344.74  Perfect plastic 
Example 2  300  15  7.5  7800  200  50  235  Bilinear plastic 
4.1.2. Response evaluation tool
In the numerical study, the full order model solutions are used as the reference for all cases. The accuracy of full order model solutions is sensitive to mesh density of structure and time step, it will be poor if the mesh is too course or the time step is too large. In order to get a relatively credible solution, many times calculation with slightly finer mesh or smaller time step, and comparison between the current results and previous results, are needed until the change occurs in results is small enough. Then, the mesh density and time step is probably acceptable for the dynamic problem.
In order to quantitatively compare the solutions of the two methods, a relative approximate ratio is introduced:
where ${\mathbf{x}}^{CB}$ is the variable vector (e.g. deformation, stress, velocity etc.) computed from improved CB method and ${\mathbf{x}}^{Ref}$ is the variable vector computed form full model order. The value of ${R}_{a}$ close to 1.0 indicates strong similarity between the two methods, where value far away from 1 indicates non similarity.
4.2. A simply supported beam
This example is similar to the problem previously investigated by many researchers [3537]. With many times calculation, the beam is meshed using 40 beam elements, the time step equals 1.0e^{5} sec, and the full order model solutions will be used as reference. For the application of the improve CB method, the original finite element model of beam is partitioned into 2 substructures uniformly. The deflection, velocity and stress responses of beam are simulated. All the analysis are carried for 0.02 sec.
4.2.1. Convergence and accuracy
It is known that stress calculation is complicated and needing more DOFs compared with displacement calculation. The stress response of this beam can be adopted to evaluate the convergence of the proposed method.
Fig. 5Convergence study of CB method for example 1
a) Stress response
b) Statistical result of ${R}_{a}$
For this example, the stress response of midpoint at bottom surface is employed. It is obvious that the results from improved CB method and full order model differ very slightly, the differences only occure at the inflection point of the time history plot in Fig. 5(a). The increase of the number of fixedinterface normal modes can reduce the differences to a certain extent. To quantify the influence of normal modes, the relative approximate ratio is calculated in accordance with Eq. (33) and shown in Fig. 5(b). As the number of fixedinterface normal modes is increased, the relative approximate ratio is close to 1.0. This means that the solution of proposed method is convergent to the solution from full order model. And a relatively fast speed of convergence can also be found from Fig. 5(b), the solution from the improved CB method with > 98 % relative approximate ratio is obtained when the fixedinterface normal mode of substructure retains 2. For the succeeding simulations of this example, the division strategy of 2 substructures with 2 fixedinterface normal modes will be used always.
The midpoint deflection response and midpoint velocity response of the beam are shown in Fig. 6. Further, the full order model solutions are also given for comparison. It is obviously that the solutions of proposed method agree well with that of full order model. There is a strong similarity between the two methods for this example.
Fig. 6Responses of the simply supported beam
a) Deflection response (${R}_{a}=$ 99.4 %)
b) Velocity response (${R}_{a}=$ 99.2 %)
4.2.2. Plastic deformation analysis
It is known that the plastic deformation will weaken the stiffness of structure. Accordingly, the tangent frequencies will change, and the more largescale plastic deformation occurs, the more severely the tangent frequencies vary. Fig. 7 shows the variation of the firstorder frequency and the secondorder frequency of the simply supported beam during the forced motion. In this example, the firstorder frequency has once decreased to 41 Hz at time 3.14 ms, only 3 % of fundamental frequency of the beam. This indicates that the beam has undergone largescale plastic deformation at time 3.14 ms.
Fig. 7Variation of tangent frequencies of simply supported beam.
a) Firstorder frequency
b) Secondorder frequency
In order to estimate the plastic deformation at time 3.14 ms, the stress distribution of beam is calculated using the proposed method, shown in Fig. 8. For a clear description, the ratio of height and length of beam is expanded 5 times. The points where the normal stress exceeds 340 MPa are filled with grey grids. The filled regions will be considered as occurrence of plastic deformation, the area of these regions accounts for 19.5 % of the entire beam. Thus, the simply supported beam can be proved with lagerscale plastic deformation at time 3.14 ms, and will show strong nonlinear behaviors during the forced motion.
Fig. 8Plastic deformation of simply supported beam at time 3.14 ms
4.3. A fixed beam
In second example, the time variation of the concentrated pulsing load is given from a piecewise function. The values of parameters in the piecewise function are: ${t}_{1}=$2 ms, ${t}_{2}=$12 ms, ${t}_{3}=$14 ms, ${F}_{m}=$4${f}_{s}$. Where ${f}_{s}$ is the static collapse load defined as the concentrated load causing a plastic hinge to occur at the midpoint of the beam. By some numerical tests, the fixed beam is discretized into 120 beam elements, and time step equals 1.0e^{4} sec. Then, the finite element model is partitioned into 2 substructures uniformly, the response of deflection, velocity and stress are simulated with the presented method. All the solutions are carried out for 0.1 sec and compared with the full order model results.
4.3.1. Convergence and accuracy
We still choose the stress response of midpoint at bottom surface of beam to evaluate the convergence of the improved CB method. As expected, the solution of CB method is convergent quickly to that from full order model, and is shown in Fig. 9. When the fixedinterface normal modes retain 4, the solution of CB method has a high relative approximate ratio over than 98 %. Thus, the division strategy of 2 substructures with 4 fixedinterface normal modes will be used for the succeeding simulation.
Fig. 9Convergence study of CB method for example 2
a) Stress response
b) Statistical result of ${R}_{a}$
Fig. 10 shows the predicted deflection and velocity response of the fixed beam. The overall solutions of proposed method have a quite accurate as full order model as expected, the average approximate ratio between the two method is close to 99 %.
Fig. 10Responses of the fixed beam
a) Deflection response (Ra = 99.2 %)
b) Velocity response (Ra = 99.0 %)
4.3.2. Calculation of plastic deformation
The variations of firstorder and secondorder tangent frequencies are shown in Fig. 11, which are similar in shape to the variation plot of example 1 (in Fig. 7). The Fig. 11 also shows the variation scope, at time 12.1 ms, the firstorder frequency is decreased to the minimum 160 Hz (about 59 % of fundamental frequency of the fixed beam). According, the stress distribution of fixed beam at time 12.1 ms is shown in Fig. 12 (the ratio of height and length is expanded 10 times). The region filled with grey grid stands for where the normal stress exceeds 235 MPa. The area of these regions is 26.7 % of the entire beam.
Fig. 11Variation of tangent frequencies of fixed beam
a) Firstorder frequency
b) Secondorder frequency
Fig. 12Plastic deformation of fixed beam at time 12.1 ms
It can be concluded that by the responses and comparison study, the fixed beam has gone through largescale plastic deformation under the pulsing load, the improved CB method is valid for this problem, and the responses obtained from proposed method have the same accuracy as that obtained from full order model.
4.4. Computational efficiency study
In this section, we aim at the two examples to study the computational efficiency of the proposed method. Although, the examples are simple beam structure, the results still show some intrinsic properties of the proposed method. In the whole simulation, the eigenvalue problems of substructure are solved with subspace iteration method.
4.4.1. The impact of application of reduced model in iterations
For the proposed method (even the original CB method), only internal DOFs of substructure are reduced, all boundary DOFs are kept completely. This is typically a disadvantage of the method if the original structure is partitioned into too many substructures, this will lead to a relative largeorder reduced model. From this point of view, the number of substructures should be minimized. However, fewer divisions in original structure may result in nonsignificant degradation of size for eigenvalue problem, and more normal modes may be retained for the accuracy solution. Such divisions whether will bring computational benefit is uncertain. In practical application of original CB method, the divisions can be in conjunction with some mode selection technique to improve the accuracy and efficiency.
In order to evaluate the efficiency of the improved CB method in the analysis of elasticplastic dynamic problem, many different substructure divisions are designed for the two examples. Four kinds of substructure division designated as case AD are used in the first example. Five other divisions designated as case EI are used in the second example. These division strategies are listed in Table 2. The comparisons of the computational time (the central processing unit (CPU) time) and obtained benefit in different analysis are listed in Table 3 and Table 4.
Table 2Substructure divisions for the two examples
Cases  Case A  Case B  Case C  Case D  Case E  Case F  Case G  Case H  Case I 
Division strategy  ${N}_{s}=$2; $k=$2  ${N}_{s}=$4; $k=$1  ${N}_{s}=$5; $k=$1  ${N}_{s}=$8; $k=$1  ${N}_{s}=$2; $k=$4  ${N}_{s}=$3; $k=$5  ${N}_{s}=$4; $k=$2  ${N}_{s}=$6; $k=$1  ${N}_{s}=$8; $k=$1 
Example  Simply supported beam  Fixed beam 
Table 3Comparisons for example 1 with different division strategies
Cases  Updating cost (ms)  CPU time (s)  Time ratio  ${R}_{a}$ 
Ref (FEM)  4.34  63.5  1.00  100 % 
Case A  7.00  29.1  0.458  98.6 % 
Case B  5.16  23.8  0.375  99.5 % 
Case C  5.16  25.3  0.398  99.5 % 
Case D  5.31  26.7  0.420  99.7 % 
Table 4Comparisons for example 2 with different division strategies
Cases  Updating cost (ms)  CPU time (s)  Time ratio  ${R}_{a}$ 
Ref (FEM)  13.27  153.8  1.00  100 % 
Case E  62.89  70.0  0.455  98.4 % 
Case F  46.6  54.9  0.357  98.2 % 
Case G  27.45  34.5  0.224  98.9 % 
Case H  18.01  26.8  0.174  99.2 % 
Case I  15.85  24.9  0.162  99.6 % 
The updating cost is the average time cost of updating stiffness matrix in simulation with full order model, or of updating transformation matrix in CB method, in each time step. The time ratio is defined as ${T}_{CB}/{T}_{FEM}$, where ${T}_{CB}$ is the CPU time consumed for simulating with improved CB method, and ${T}_{FEM}$ is the CPU time while simulating with the full order model. All the analysis results are ensured to have greater than 98 % relative approximate ratio compared with full order model, the relative approximate ratio is calculated based on stress response. From these two tables we can make the following observation:
• The numerical tests show that proposed method is more efficient than full order model for the elasticplastic dynamic problem, and the solution of both method have a high similarity. The computational time cost of proposed method is only about 0.16~0.45 times as the full order model.
• As expected, the time cost of updating transformation matrix in improved CB method is a little more than the time cost of updating stiffness matrix in full order model, provided that an appropriate division strategy is applied (e.g. case B, C and D for example 1, case I for example 2). This can be attributed to two reasons. First, the eigenvalue problems are all smallscale by dividing the original structure into many substructures. Then, considering the modal truncation, there are only a few lower normal modes are needed to be calculated. After the application of reduced model in iterations in a time step, the computational time cost of solving governing equations has reduced obviously compared with full model analysis. The larger the structural system is, the more likely it is.
• It can also be known that an appropriate substructure division of the original structure will improve computational efficiency significantly on the premise of maintaining same accuracy, e.g., using the division of case I is more than tripled efficient compared with the division of case E for example 2.
Updating the transformation matrix conditionally may be an effective technique to improve the efficiency of the proposed method. To demonstrate this, we choose the most efficient division strategies in Table 3 and Table 4 (they are case B for example 1 and case I for example 2), and update the transformation matrix conditionally in accordance with Eq. (33). The predefined tolerance ε is going to set to 0.05, 0.02 and 0.01.
Table 5 shows the loss of accuracy and saving of computational time for these analyses. With the loose of prescribed tolerance $\epsilon $, the computational efficiency is improved further, as expected. And, the accuracy of solution only degrades slightly. For example 2, the improvement is remarkable especially.
Table 5Validity verification of update strategy
Cases  Update strategy  Update times  CPU time (s)  Time ratio  ${R}_{a}$ 
Case B  Unconditionally  1000  23.8  0.375  99.5 % 
$\epsilon =$0.01  843  22.6  0.356  98.0 %  
$\epsilon =$0.02  769  22.2  0.350  94.0 %  
$\epsilon =$0.05  721  22.0  0.346  93.4 %  
Case I  Unconditionally  1000  24.9  0.162  99.6 % 
$\epsilon =$0.01  864  23.1  0.149  97.8 %  
$\epsilon =$0.02  859  22.7  0.148  97.6 %  
$\epsilon =$0.05  42  13.5  0.088  96.7 % 
4.4.2. The impact of time step
The use of reduced model can have larger time step than use of full order model. On the premise of maintaining an acceptable accuracy of solutions, the larger time step may result in a shorter computational time and high efficiency. Here, we still choose the most efficient division strategies (case B for example 1 and case I for example 2) to calculate stress responses of structure with relative larger time step. And the accuracy of solution will be compared with that obtained from full order model at the same time step.
Table 6 and Table 7 have recorded the computational time for the two methods with relative larger time step. In general, the improved CB method has the same numerical accuracy as full order model when the time step getting greater. Meanwhile, we also have found that the time cost of solving equations increases dramatically in full order moder simulation with the increase of time step. This may be due to the stability of equations in physical coordinate is sensitive to time step and worsen quickly when greater time step is used. This leads to a considerable computational cost for gaining convergent results in each times step. However, the improved CB method can overcome this disadvantage because it converts the equations into modal coordinate to make equations more stable. In our numerical test, the time cost of solving equations for improved CB method almost maintains a constant value with different time step. So, the efficiency of the proposed method for the examples is obvious compared with full order model.
Table 6Comparisons of the computational time for example 1 with different time step
Method  $\mathrm{\Delta}t$(s)  CPU time (s)  Time ratio  ${R}_{a}$ 
Ref (FEM)  1.0e^{5}  63.5  1.0  100 % 
FEM  5.0e^{5}  18.4  0.290  92.7 % 
CB  5.0e^{5}  4.8  0.075  92.1 % 
FEM  1.0e^{4}  21.6  0.340  72.5 % 
CB  1.0e^{4}  2.4  0.038  69.5 % 
Table 7Comparisons of the computational time for example 2 with different time step
Method  $\mathrm{\Delta}t$ (s)  CPU time (s)  Time ratio  ${R}_{a}$ 
Ref (FEM)  1.0e^{4}  153.8  1.00  100 % 
FEM  2.5e^{4}  86.7  0.564  91.1 % 
CB  2.5e^{4}  13.88  0.09  90.1 % 
FEM  5.0e^{4}  134.6  0.875  85.18 % 
CB  5.0e^{4}  7.10  0.046  82.28 % 
In conclusion, the causes of good efficiency of proposed method can be summarized into two. First of all, if the time step required for the proposed method and full order model is same, the computational efforts of iterations are reduced when the presented method is utilized. The second cause is that the reduced order model is more stable for using greater time step compared with the full order finite element model. Generally, both causes of efficiency usually go together.
5. Conclusions
For strong nonlinear problem, the transient analysis of structures with largescale plastic deformation, an improved CB method is proposed. Unlike the original CB method, the transformation matrix of the improved CB method is constructed with a type of nonlinear modes of structure (tangent modes). To do this, the nonlinear problem is linearized over a small time interval, the tangent modes of structure based on the nonlinear state at the beginning of the time interval are proved to be orthogonal with respect to mass matrix as well as with respect to tangent stiffness matrix from mathematic view, by incorporating the elasticplastic material behavior.
Using the improved CB method, the response solutions of structure with largescale plastic deformation are accurate as the solutions obtained from FEM, and the computational time is reduced substantially. The good performance of the improved CB method is demonstrated through two examples of simply supported beam and fixed beam loaded impulsively.
Moreover, the proposed method can be easily extended to the case of large of more complex structure and offer an alternative tool to aid in the dynamic design of structure with the consideration of largescale plastic deformation.
References

Zienkiewicz O. C., et al. The Finite Element Method for Solid and Structural Mechanics. Sixth Edition, Elsevier Ltd., Amsterdam, 2006.

Lulf F. A. Reduced bases for nonlinear structural dynamic systems: a comparative study. Journal of Sound and Vibration, Vol. 332, Issue 15, 2013, p. 38973921.

Corigliano A., Dossi M., Mariani S. Model Order Reduction and domain decomposition strategies for the solution of the dynamic elasticplastic structural problem. Computer Methods in Applied Mechanics and Engineering, Vol. 290, Issue 15, 2015, p. 27155.

Nickell R. E. Nonlinear dynamics by mode superposition. Computer Methods in Applied Mechanics and Engineering, Vol. 7, Issue 7, 1976, p. 107129.

Stricklin J. A., Haisler W. E. Formulations and solution procedures for nonlinear structural analysis. Computers and Structures, Vol. 7, Issue 1, 1977, p. 125136.

Idelsohn S. R., Cardona A. A reduction method for nonlinear structural dynamic analysis. Computer Methods in Applied Mechanics and Engineering, Vol. 49, Issue 3, 1985, p. 253279.

Shah V., Bohm G., Nahavandi A. Modal superposition method for computationally economical nonlinear structural analysis. Journal of Pressure Vessel Technology, Vol. 101, Issue 2, 1979, p. 134141.

Morris N. F. The use of modal superposition in nonlinear dynamics. Computers and Structures, Vol. 7, Issue 1, 1977, p. 6572.

Mohraz B., Elghadamsi F. E., Chang C.J. An incremental modesuperposition for nonlinear dynamic analysis. Earthquake Engineering and Structural Dynamics, Vol. 20, Issue 5, 1991, p. 471481.

Bathe K. J., Gracewski S. On nonlinear dynamic analysis using substructuring and mode superposition. Computers and Structures, Vol. 13, Issues 56, 1981, p. 699707.

Kukreti A. R., Issa H. I. Dynamic analysis of nonlinear structures by pseudonormal mode superposition method. Computers and Structures, Vol. 19, Issue 4, 1984, p. 653663.

Muscolino G. Modesuperposition methods for elastoplastic systems. Journal of Engineering Mechanics, Vol. 115, Issue 10, 1989, p. 21992215.

Villaverde R., Hanna M. M. Efficient mode superposition algorithm for seismic analysis of nonlinear structures. Earthquake Engineering and Structural Dynamics, Vol. 21, Issue 10, 1992, p. 849858.

Manoach E. Dynamic response of elastoplastic mindlin plate by mode superposition method. Journal of Sound and Vibration, Vol. 162, Issue 1, 1993, p. 165175.

Manoach E., Karagiozova D. Dynamic response of thick elasticplastic beams. International Journal of Mechanical Sciences, Vol. 35, Issue 11, 1993, p. 909919.

Hurty W. C. Dynamic analysis of structural systems using component modes. AIAA Journal, Vol. 3, Issue 4, 1965, p. 678685.

Craig R. R. Jr., Bampton M. C. C. Coupling of substructures for dynamic analyses. AIAA Journal, Vol. 6, Issue 7, 1968, p. 13131319.

Benfield W., Hruda R. Vibration analysis of structures by component mode substitution. AIAA Journal, Vol. 9, Issue 7, 1971, p. 12551261.

Macneal R. H. A hybrid method of component mode synthesis. Computers and Structures, Vol. 1, Issue 4, 1971, p. 581601.

Rubin S. Improved componentmode representation for structural dynamic analysis. AIAA Journal, Vol. 13, Issue 8, 1975, p. 9951006.

Suarez L., Singh M. Improved fixed interface method for modal synthesis. AIAA Journal, Vol. 30, Issue 12, 1992, p. 29522958.

Kim J. G., Lee P. S. An enhanced CraigBampton method. International Journal for Numerical Methods in Engineering, Vol. 103, Issue 2, 2015, p. 7993.

Klerk D. D., Rixen D., Voormeeren S. General framework for dynamic substructuring: history, review and classification of techniques. AIAA Journal, Vol. 46, Issue 5, 2008, p. 11691181.

Seshu P. Substructuring and component mode synthesis. Shock and Vibration, Vol. 4, Issue 3, 1997, p. 199210.

Craig J. R. R. Substructure methods in vibration. Journal of Vibration and Acoustics, Vol. 117, 1995, p. 207213.

Wu S. C., Haug E. J. A substructure technique for dynamics of flexible mechanical systems with contactimpact. Journal of Mechanical Design, Vol. 112, Issue 3, 1990, p. 390398.

Guo A., Batzer S. Substructure analysis of a flexible system contactimpact event. Journal of Vibration and Acoustics, Vol. 126, Issue 1, 2004, p. 126131.

Shen Y., Yin X. Dynamic substructure analysis of stress waves generated by impacts on nonuniform rod structures. Mechanism and Machine Theory, Vol. 74, 2014, p. 154172.

Shen Y., Yin X. Analysis of geometric dispersion effect of impactinduced transient waves in composite rod using dynamic substructure method. Applied Mathematical Modeling, Vol. 40, Issue 3, 2016, p. 19721988.

Kim D. K., Bae J. S., Lee I., et al. Dynamic model establishment of a deployable missile control fin with nonlinear hinge. Journal of Spacecraft and Rockets, Vol. 42, Issue 1, 2005, p. 6677.

Krishna I. P., Padmanabhan C. Improved reduced order solution techniques for nonlinear systems with localized nonlinearities. Nonlinear Dynamics, Vol. 63, Issue 4, 2011, p. 561586.

Alamdari M. M., Li J., Samali B., et al. Nonlinear joint model updating in assembled structures. Journal of Engineering Mechanics, Vol. 140, Issue 7, 2013, p. 04014042.

Kuether R. J., Allen M. S. CraigBampton substructuring for geometrically nonlinear subcomponents. Proceedings of the Dynamics of Coupled Structures, Vol. 1, 2014, p. 167178.

Bond J., Khraishi T. Transient nonlinear simulation with component mode synthesis. International Journal of Mechanics and Materials in Design, Vol. 5, Issue 4, 2009, p. 365380.

Nagarajan S., Popov E. P. Elasticplastic dynamic analysis of axisymmetric solids. Computers and Structures, Vol. 4, Issue 6, 1974, p. 1171134.

Yang T., Saigal S. A simple element for static and dynamic response of beams with material and geometric nonlinearities. International Journal for Numerical Methods in Engineering, Vol. 20, Issue 5, 1984, p. 851867.

Liu S. C., Lin T. H. Elasticplastic dynamic analysis of structures using known elastic solutions. Earthquake Engineering and Structural Dynamics, Vol. 7, Issue 2, 1979, p. 147159.
About this article
This research has been supported by Postdoctoral Fund of JiangSu (1402007B), the support is gratefully acknowledged.