Published: 31 March 2017

Improved Craig-Bampton method for transient analysis of structures with large-scale plastic deformation

Peng-bo Qian1
Lin-fang Qian2
Xiao-chun Yin3
1Huatian Engineering and Technology Corporation, MCC, Nanjing, Jiangsu, 210094, China
2, 1College of Mechanical Engineering, Nanjing University of Science and Technology, Nanjing, Jiangsu, 210094, China
3Department of Mechanics and Engineering Science, Nanjing University of Science and Technology, Nanjing, Jiangsu, 210094, China
Corresponding Author:
Xiao-chun Yin
Views 478
Reads 246
Downloads 1534

Abstract

This paper reports on the improvement of Craig-Bampton (CB) method for transient analysis of structures with large-scale 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 elastic-plastic 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 elastic-plastic beams loaded impulsively are examined. Simulation results show that the improved CB method is valid and accurate for transient analysis of structures with large-scale plastic deformation and has lower computational cost compared with full order model.

1. Introduction

The elastic-plastic dynamic analysis of structure is typically performed using finite element method (FEM), and has been remarked as a complicated and time-consuming 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 elastic-plastic dynamic analysis is substantially high. In order to mitigate this problem, model order reduction technique is performed [2-3].

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 [5-7]. 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 right-hand side of the governing equations. This technique is the so-called pseudo-force 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 fixed-interface 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 fixed-interface normal modes to generate the reduced modal space. The proposed method is referred as Craig-Bampton (CB) method. After that, various CMS methods have been developed [18-22]. 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 [23-25].

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 addition-deletion technique. The contact force is addressed as internal force of system. And more deep studies about contact impact problem then were developed [27-29] 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 [30-32].

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 large-scale 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 large-scale plastic deformation.

The process of this work can be summarized as follows: First, the incremental governing equations of elastic-plastic 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 elastic-plastic 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 elastic-plastic 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 elastic-plastic dynamic problem

Considering the small strain and deformation type of a motion of an elastic-plastic finite body Ω, by use of the nonlinear finite element theory, the governing equations for the general elastic-plastic structure during a small time interval (tt+Δt) can be expressed as [1]:

1
MΔu¨+KtΔu=Pt+Δt-Mu¨t+Fut=ΔQ,

where M is the constant, symmetric, positive-definite mass matrix. Kt is the tangent stiffness matrix at time t, is dependent on the nonlinear constitutive relation of material. Pt+Δt is the excited force vector at time t+Δt. F(ut) is the resisting force vector of structure at time t, depends on the displacement field. The u¨t,ut,Δu¨ and Δ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 elastic-plastic 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:

2
M=eGΤmeG=eGΤVe NΤρNdv,
3
Kt=eGΤkteG=eGΤVe BΤDeptBdvG,
4
Pt=eGΤpteG=eGΤSe NΤTedv,

where Ve and Se denote the volume and surface of element, respectively. G is a Boolean matrix that represents the transformation from the element coordinate to the physical coordinate. N is the matrix of shape function of element, B is strain-displacement matrix of element. ρ is the mass density, Dep(t) is the elastic-plastic matrix of material and Te is the boundary surface traction vector of an element. The elastic-plastic matrix Dep(t) 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 step-by-step 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 (Newton-Raphson) and MNR (Modified Newton Raphson) method. Here, we employed the Newmark-β 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-β method has some assumptions on velocity and displacement, which are:

5
u˙t+Δt=u˙t+[1-αu¨t+αu¨t+Δt]Δt,
6
ut+Δt=ut+u˙tΔt+1/2-βu¨tΔt2+βu¨t+ΔtΔt2,

where α and β are integration parameter. Considering the relation of Δu¨=u¨t+Δt-u¨t,Δu˙=u˙t+Δt-u˙t, Δu=ut+Δt-ut+Δt, this yields:

7
Δu¨=ΔuβΔt2-u˙tβΔt-u¨t2β,
8
Δu˙=u¨tΔt+αΔuβΔt-αu˙tβ-αΔtu¨t2β.

Using Eq. (7-8) into Eq. (1), the corresponding effective static equilibrium equation of the nonlinear system is given by:

9
M/βΔt2+KtΔu=Pt+Δt-Mu¨t+Fut+Mu˙tβΔt+Mu¨t2β.

For a simple description of Eq. (9):

10
K^tΔu=ΔQ^ ,

where K^t denotes the effective stiffness matrix, ΔQ^ stands for the effective force.

In MNR procedure, the iterations are governed by:

11
K^tΔui=ΔQ^i-1 ,
12
ΔQ^i-1=Pt+Δt-Mu¨t+Δti-1+eVe BΤσt+Δti-1dv+Mu˙tβΔt+Mu¨t2β.

The iterations start from: ut+Δt0=ut, u˙t+Δt0=u˙t, u¨t+Δt0=u¨t, σt+Δt0=σ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: ut+Δt =ut+Δu, u˙t+Δt =u˙t+Δu˙, u¨t+Δt =u¨t+Δ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 (tt+Δt) and the incremental displacements are computed with many iterations. In order to reduce the computational effort of the analysis, we consider using C-B method.

3. Formulations of improved C-B method for elastic-plastic dynamic problem

The C-B method is a modal-based 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 fixed-interface 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 (tt+Δt), it seems involuntary to apply C-B method for elastic-plastic 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 elastic-plastic problem. To this end, it is necessary to deeply understand the property of tangent modes (e.g. orthogonality and linear independence) for elastic-plastic 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 (tt+Δ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:

13
x=x- eiωt,

where i=-1 is the imaginary unit, ω is a real number and x- represents the amplitudes of generalized displacement. The substitution of Eq. (13) into the free vibration yields an eigenvalue problem:

14
Kt-wi2M=Φi,

where wi2 is eigenvalues and Φ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

The difference of tangent modes and actual modes of structure with plastic deformation

3.2. Properties of tangent modes for elastic-plastic dynamic problem

In this section, we will further discuss the property of tangent modes by incorporating the elastic-plastic material behavior. As was noted in Eq. (1), M is real symmetric and positive definite for elastic-plastic dynamic problem, the solutions of tangent eigenvalue problem in Eq. (14) are sensitive to Kt from matrix-algebra theory. Generally, the eigenvalue problem of elastic-plastic 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 Kt. According to Eq. (3), Kt is completely dependent on the elastic-plastic matrix Dep(t), which is a function of plastic law and plastic deformation.

The expression of Dep(t) 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 Prandtl-Reuss 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 elastic-plastic matrix of material is defined as [1]:

15
Dept=De-DefσtfσtΤDefσtΤDefσt+C=De-1HDeStStΤDe,

where De denotes the elastic matrix of material, and is real symmetric. f and σ(t) are the yield function of material and component of stress vector, respectively. S(t) 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 De is real symmetric, the DeS(t)S(t)ΤDe and Dept are real symmetric, correspondingly. We use Dept into Eq. (3) and can be proved that the Kt is also a real symmetric matrix:

16
KtT=eGΤVe BΤDep(t)BdvGΤ=eGΤVe BΤDeptΤBdvG=Kt.

When the structure is upon plastic loading, the following conditions will hold true for any point within the structure:

• For hardening material:

17
dσijtdεijt>0,

• For ideal plastic material:

18
dσijtdεijt=0.

We consolidate the two formulas and make a simple presentation with matrix form. Let:

19
dσtΤdεt0.

According to plastic theory, dσt=Deptdεt, the Eq. (19) can be rewritten as:

20
dσtΤdεt=dεtΤDeptdεt0.

In an element, the strain increment of any point can be calculated from dεt=Bdue(t), this yields:

21
dσtΤdεt=dεtΤDeptdεt=duetΤBΤDeptBdue(t)0,

where due(t) is a vector containing the incremental displacements of all DOFs in an element. We integrate Eq. (21) in the region of an element:

22
Ve dσtΤdεtdv=duetΤVe BΤDeptBdvduet
=duetΤkteduet0.

And assemble all elements within the structure. The result will take the form:

23
e Ve dσtΤdεtdv=e du(t)TGTkteGdu(t)
=(Gdu(t))TKt(Gdu(t))0,

where du(t) is the displacement increment vector of structure. Considering the arbitrariness of vector du(t), the Kt is proved to be semi-positive.

From matrix-algebra theory, it is easy known that all roots of Eq. (14) are real and non-negative, 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 elastic-plastic dynamic problem can be summarized in Fig. 2.

Fig. 2The basic procedure of C-B method for elastic-plastic 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

The basic procedure of C-B method for elastic-plastic 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)

The basic procedure of C-B method for elastic-plastic 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

b)

The basic procedure of C-B method for elastic-plastic 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

c)

The basic procedure of C-B method for elastic-plastic 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

d)

3.3.1. Modal analysis of substructure

The original structure (discretized model) is first partitioned into Ns substructures, which are connected with fixed interface at the boundary (Γ) between each other. On the basis of Eq. (1), the un-damped incremental governing equations of an elastic-plastic substructure during the time interval (t~t+Δt) can be written in the partitioned form:

24
miismibsmbismbbsΔu¨isΔu¨bs+(kts)ii(kts)ib(kts)bi(kts)bbΔuisΔubs=0Δqbs,

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 fixed-interface normal modes ΦsRni×ni and constrain modes ΨsRni×nj (where ni and nj stand for the number of interior DOFs and the number of boundary DOFs). The fixed-interface normal modes are defined as the free vibration modes of the constrained substructure, whose interfacial DOFs are set to have zero displacement, i.e., Δubs=0. In Eq. (24), this leads to a tangent eigenvalue problem of substructure:

25
(kts)ii-wj2miisψjs=0, j=1,2ni.

The constraint modes are the static deflections of substructure which are obtained by imposing successive unit displacement at each boundary DOF:

26
Ψs=-ktsii-1ktsib.

So, the displacement vector Δu s can be projected into a new space defined by Φs and Ψs:

27
Δu s=ΔuisΔubs=ΦsΨs0IΔaisΔubs.

We partition Φs into dominate modes Φds and high order modes Φhs. In order to reduce the size of governing equations, the high normal modes will be neglected in Eq. (27). Then, the physical displacement vector Δu s is approximated by the dominated modes:

28
Δu s=ΦdsΦhsΨs00IΔadsΔahsΔubsΦdsΨs0IΔadsΔubs=TCBsΔa s.

This is the transformation of displacement vector from physical coordinate to modal coordinate. The vector Δa s is denoted as the modal displacement, its number equals to k+nj(kni), where k is the number of retained fixed-interface normal modes.

Using TCBs in Eq. (24) and left multiplying TCBs on the two sides of the equations, the reduced governing equations of substructure s can be cast in a modal space:

29
Idsm-dbsm-bdsm-bbsΔa¨dsΔu¨bs+(Λts)d00(k-ts)bbΔadsΔubs=TCBsT0Δqbs,

where IdsRk×k is an identity matrix, (Λts)dRk×k is a diagonal matrix, if Φds is normalized. The coupling item m-dbs=(Φds)TmiisΨs+(Φds)Tmibs.

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:

30
M-Δa¨+K-tΔa=Pt+Δt-Mu¨t+Fut=ΦTΔQ=ΔQ-,

where Φ is the transformation matrix of full system assembled from all substructures. The M-=ΦTMΦ, K-t=ΦTKtΦ and Δ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:

31
IΔη¨+Λ-tΔη=ΦmTΔQ-,

where I=ΦmTM-Φm and Λ-t=ΦmTK-tΦm are diagonal matrix. Φm and η 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 equal-space 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 step-by-step 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 C-B method

From a simple perspective, the application CB method is regarded as coordinate transformation of governing equations. When it is adopted for elastic-plastic dynamic problem, the operations of generating transformation matrix will be performed at each time step. The basic premise of the application for elastic-plastic 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 fixed-interface 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 C-B method for elastic-plastic problem can be updated according to residual norm value conditionally. The weighted residual norm is evaluated as [6]:

32
ε=1NRt2(Pt2+Mu¨t2),

where Rt=Pt-F(ut)-Mu¨t is residual force, ε is the predefined tolerance and N is the DOFs number of original system.

4. Numerical results and discussion

The procedures of the improved C-B 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 large-scale plastic deformation. The integration parameters of the Newmark method with the values of α= 0.29, β= 0.56, and were used for all solutions. All computations were carried out on an Intel(R)Core(TM) i3-2350M @ 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(t) 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 pt=np0exp(-2t/T). Where p0 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

Problem descriptions of simply supported beam, a) beam element model  and substructures, b) material constitutive relation, c) load type

a)

Problem descriptions of simply supported beam, a) beam element model  and substructures, b) material constitutive relation, c) load type

b)

Problem descriptions of simply supported beam, a) beam element model  and substructures, b) material constitutive relation, c) load type

c)

Fig. 4Problem descriptions of fixed beam, a) beam element model and substructures, b) material constitutive relation, c) load type

Problem descriptions of fixed beam, a) beam element model and substructures,  b) material constitutive relation, c) load type

a)

Problem descriptions of fixed beam, a) beam element model and substructures,  b) material constitutive relation, c) load type

b)

Problem descriptions of fixed beam, a) beam element model and substructures,  b) material constitutive relation, c) load type

c)

In second example, considering a fixed beam is subjected to a concentrated pulsing load ft 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)
ρ (kg/m3)
E (GPa)
Et(GPa)
σ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:

33
Ra=1-xCBt-xReft2xReft2×100 %,

where xCB is the variable vector (e.g. deformation, stress, velocity etc.) computed from improved C-B method and xRef is the variable vector computed form full model order. The value of Ra 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 [35-37]. 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

Convergence study of CB method for example 1

a) Stress response

Convergence study of CB method for example 1

b) Statistical result of Ra

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 fixed-interface 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 fixed-interface 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 fixed-interface normal mode of substructure retains 2. For the succeeding simulations of this example, the division strategy of 2 substructures with 2 fixed-interface 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

Responses of the simply supported beam

a) Deflection response (Ra= 99.4 %)

Responses of the simply supported beam

b) Velocity response (Ra= 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 large-scale plastic deformation occurs, the more severely the tangent frequencies vary. Fig. 7 shows the variation of the first-order frequency and the second-order frequency of the simply supported beam during the forced motion. In this example, the first-order 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 large-scale plastic deformation at time 3.14 ms.

Fig. 7Variation of tangent frequencies of simply supported beam.

Variation of tangent frequencies of simply supported beam.

a) First-order frequency

Variation of tangent frequencies of simply supported beam.

b) Second-order 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 lager-scale 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

Plastic 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: t1= 2 ms, t2= 12 ms, t3= 14 ms, Fm= 4fs. Where fs 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 fixed-interface 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 fixed-interface normal modes will be used for the succeeding simulation.

Fig. 9Convergence study of CB method for example 2

Convergence study of CB method for example 2

a) Stress response

Convergence study of CB method for example 2

b) Statistical result of Ra

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

Responses of the fixed beam

a) Deflection response (Ra = 99.2 %)

Responses of the fixed beam

b) Velocity response (Ra = 99.0 %)

4.3.2. Calculation of plastic deformation

The variations of first-order and second-order 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 first-order 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

Variation of tangent frequencies of fixed beam

a) First-order frequency

Variation of tangent frequencies of fixed beam

b) Second-order frequency

Fig. 12Plastic deformation of fixed beam at time 12.1 ms

Plastic 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 large-scale 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 large-order reduced model. From this point of view, the number of substructures should be minimized. However, fewer divisions in original structure may result in non-significant 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 elastic-plastic dynamic problem, many different substructure divisions are designed for the two examples. Four kinds of substructure division designated as case A-D are used in the first example. Five other divisions designated as case E-I 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
Ns= 2;
k= 2
Ns= 4;
k= 1
Ns= 5;
k= 1
Ns= 8;
k= 1
Ns= 2;
k= 4
Ns= 3;
k= 5
Ns= 4;
k= 2
Ns= 6;
k= 1
Ns= 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
Ra
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
Ra
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 TCB/TFEM, where TCB is the CPU time consumed for simulating with improved CB method, and TFEM 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 elastic-plastic 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 small-scale 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 ε, 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
Ra
Case B
Unconditionally
1000
23.8
0.375
99.5 %
ε= 0.01
843
22.6
0.356
98.0 %
ε= 0.02
769
22.2
0.350
94.0 %
ε= 0.05
721
22.0
0.346
93.4 %
Case I
Unconditionally
1000
24.9
0.162
99.6 %
ε= 0.01
864
23.1
0.149
97.8 %
ε= 0.02
859
22.7
0.148
97.6 %
ε= 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
Δt (s)
CPU time (s)
Time ratio
Ra
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
Δt (s)
CPU time (s)
Time ratio
Ra
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 large-scale 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 elastic-plastic material behavior.

Using the improved CB method, the response solutions of structure with large-scale 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 large-scale 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. 3897-3921.
  • Corigliano A., Dossi M., Mariani S. Model Order Reduction and domain decomposition strategies for the solution of the dynamic elastic-plastic structural problem. Computer Methods in Applied Mechanics and Engineering, Vol. 290, Issue 15, 2015, p. 27-155.
  • Nickell R. E. Nonlinear dynamics by mode superposition. Computer Methods in Applied Mechanics and Engineering, Vol. 7, Issue 7, 1976, p. 107-129.
  • Stricklin J. A., Haisler W. E. Formulations and solution procedures for nonlinear structural analysis. Computers and Structures, Vol. 7, Issue 1, 1977, p. 125-136.
  • 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. 253-279.
  • 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. 134-141.
  • Morris N. F. The use of modal superposition in nonlinear dynamics. Computers and Structures, Vol. 7, Issue 1, 1977, p. 65-72.
  • Mohraz B., Elghadamsi F. E., Chang C.-J. An incremental mode-superposition for non-linear dynamic analysis. Earthquake Engineering and Structural Dynamics, Vol. 20, Issue 5, 1991, p. 471-481.
  • Bathe K. J., Gracewski S. On nonlinear dynamic analysis using substructuring and mode superposition. Computers and Structures, Vol. 13, Issues 5-6, 1981, p. 699-707.
  • Kukreti A. R., Issa H. I. Dynamic analysis of nonlinear structures by pseudo-normal mode superposition method. Computers and Structures, Vol. 19, Issue 4, 1984, p. 653-663.
  • Muscolino G. Mode-superposition methods for elastoplastic systems. Journal of Engineering Mechanics, Vol. 115, Issue 10, 1989, p. 2199-2215.
  • Villaverde R., Hanna M. M. Efficient mode superposition algorithm for seismic analysis of non-linear structures. Earthquake Engineering and Structural Dynamics, Vol. 21, Issue 10, 1992, p. 849-858.
  • Manoach E. Dynamic response of elastoplastic mindlin plate by mode superposition method. Journal of Sound and Vibration, Vol. 162, Issue 1, 1993, p. 165-175.
  • Manoach E., Karagiozova D. Dynamic response of thick elastic-plastic beams. International Journal of Mechanical Sciences, Vol. 35, Issue 11, 1993, p. 909-919.
  • Hurty W. C. Dynamic analysis of structural systems using component modes. AIAA Journal, Vol. 3, Issue 4, 1965, p. 678-685.
  • Craig R. R. Jr., Bampton M. C. C. Coupling of substructures for dynamic analyses. AIAA Journal, Vol. 6, Issue 7, 1968, p. 1313-1319.
  • Benfield W., Hruda R. Vibration analysis of structures by component mode substitution. AIAA Journal, Vol. 9, Issue 7, 1971, p. 1255-1261.
  • Macneal R. H. A hybrid method of component mode synthesis. Computers and Structures, Vol. 1, Issue 4, 1971, p. 581-601.
  • Rubin S. Improved component-mode representation for structural dynamic analysis. AIAA Journal, Vol. 13, Issue 8, 1975, p. 995-1006.
  • Suarez L., Singh M. Improved fixed interface method for modal synthesis. AIAA Journal, Vol. 30, Issue 12, 1992, p. 2952-2958.
  • Kim J. G., Lee P. S. An enhanced Craig-Bampton method. International Journal for Numerical Methods in Engineering, Vol. 103, Issue 2, 2015, p. 79-93.
  • 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. 1169-1181.
  • Seshu P. Substructuring and component mode synthesis. Shock and Vibration, Vol. 4, Issue 3, 1997, p. 199-210.
  • Craig J. R. R. Substructure methods in vibration. Journal of Vibration and Acoustics, Vol. 117, 1995, p. 207-213.
  • Wu S. C., Haug E. J. A substructure technique for dynamics of flexible mechanical systems with contact-impact. Journal of Mechanical Design, Vol. 112, Issue 3, 1990, p. 390-398.
  • Guo A., Batzer S. Substructure analysis of a flexible system contact-impact event. Journal of Vibration and Acoustics, Vol. 126, Issue 1, 2004, p. 126-131.
  • Shen Y., Yin X. Dynamic substructure analysis of stress waves generated by impacts on non-uniform rod structures. Mechanism and Machine Theory, Vol. 74, 2014, p. 154-172.
  • Shen Y., Yin X. Analysis of geometric dispersion effect of impact-induced transient waves in composite rod using dynamic substructure method. Applied Mathematical Modeling, Vol. 40, Issue 3, 2016, p. 1972-1988.
  • 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. 66-77.
  • Krishna I. P., Padmanabhan C. Improved reduced order solution techniques for nonlinear systems with localized nonlinearities. Nonlinear Dynamics, Vol. 63, Issue 4, 2011, p. 561-586.
  • 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. Craig-Bampton substructuring for geometrically nonlinear subcomponents. Proceedings of the Dynamics of Coupled Structures, Vol. 1, 2014, p. 167-178.
  • Bond J., Khraishi T. Transient non-linear simulation with component mode synthesis. International Journal of Mechanics and Materials in Design, Vol. 5, Issue 4, 2009, p. 365-380.
  • Nagarajan S., Popov E. P. Elastic-plastic dynamic analysis of axisymmetric solids. Computers and Structures, Vol. 4, Issue 6, 1974, p. 117-1134.
  • 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. 851-867.
  • Liu S. C., Lin T. H. Elastic-plastic dynamic analysis of structures using known elastic solutions. Earthquake Engineering and Structural Dynamics, Vol. 7, Issue 2, 1979, p. 147-159.

About this article

Received
15 June 2016
Accepted
12 November 2016
Published
31 March 2017
SUBJECTS
Mechanical vibrations and applications
Keywords
Craig-Bampton method
component mode synthesis
plastic deformation
transient response
reduction method
Acknowledgements

This research has been supported by Postdoctoral Fund of JiangSu (1402007B), the support is gratefully acknowledged.