Algorithm of DRM with Kinetic Damping for Finite Element Static Solution of Strain-Softening Structures

In order to deal with the divergence and instability due to the ill-posedness of the nonlinear finite element (FE) model of strainsoftening structure in implicit static analysis, the dynamic relaxation method (DRM) was used with kinetic damping to solve the static increments in the incremental solution procedure so that the problem becomes well-posed. Moreover, in DRM there is no need to assemble and inverse the stiffness matrix as in implicit static analysis such that the associated computational cost is avoided. The ascending branch of static equilibrium path was solved by load increments, while the peak point and the descending branch were solved by displacement increments. Two numerical examples illustrated the effectiveness of such application of DRM in the FE analysis of static equilibrium path of strain-softening structures.


Introduction
The static equilibrium path analysis is important in structural analysis.The finite element (FE) method is an important means for structural analysis.In nonlinear FE analysis of structures, the incremental method is the key [1].The static equilibrium path analysis of a strain-softening structure is a tough problem due to the existence of descending branch after peak load point.Concrete is a well-known strainsoftening material widely used in civil engineering.
According to the increment control type, ISA can be divided into load-control, displacement control, and arclength control.Firstly, for the load-control method, it was proven mathematically that ISA such as Newton method, modified Newton method, and quasi-Newton method cannot snap the load peak, search the descending branch of static equilibrium path, and handle local instability because of the nonpositive stiffness matrix K or the instability of node displacement a [4].Secondly, comparing to the load-control method, usually a better convergence rate as well as crossing the load peak point can be obtained by adopting the displacecontrol method [5].However, not only do dissymmetry and disbandment of K raise the compute cost [6], but also it is difficult to select the unstable displacement components as displacement control ones [7].Otherwise, it was reported that displacement control method cannot trace the snapback equilibrium path [8] and cannot be fitting for analyzing strain-softening structure [5].Thirdly, the arc-length control method [9], which is usually used in analysis of structure static problem, can snap load peak point and track descending branch caused by geometrical nonlinearity [8] but cannot deal with local instability problem caused by strain-softening or local buckling and cannot snap the bifurcation point [2,10,11].For these reasons, Duan proposed local arclength method which decomposes a structure into damaged zone, failure zone, and dysfunction zone and when doing structure analysis by the local arc-length method, the constraint equation only contains the displacement increment of nodes in failure zone and dysfunction zone [11].From the numerical case studies that have been reported, local arc-length method can conduct analysis of strain-softening structure, but these studies mainly deal with strain-softening under tension [11,12]; moreover, the constraint equation needs to be updated during incremental-iterative procedure 2 Advances in Materials Science and Engineering according to the changing of damaged zone and failure zone.
The literature [2] suggests employing explicit dynamic algorithm (EDA) to do static analysis of structures characterized with local instability; and, in this way, static problem is treated as quasi-static problem in which the inertia effect can be ignored.And we define the EDA as EQSA when it is used to solve quasi-static problem.Some applications of EQSA can be seen in literature [13,14].Being different from implicit algorithm, EQSA uses explicit time integral formula to update displacement vector, which can avoid the construction and inversion of K [5].However, there are difficulties in using EQSA.First of all, EQSA is a conditionally stable algorithm; it is that the stable solution of static problems depends on the minimum size of finite elements under the hypothesis that every element has the same physical properties [2].This condition increases the computational cost of the discrete model with small-dimension elements.Secondly, it is required to reduce the loading rate if a reasonable pseudo static solution is prospected, but it will increase the number of time steps (computational cost).Increasing the loading rate will introduce noise in the loading response (nonnegligible inertia effect) [2] and will even make the computation result completely deviate from the test result [15].In order to reduce the computational cost, the node mass scaling may also introduce noise into the load response [2].The determination of viscous damping coefficient has a strong empirical dependence [2].
The explicit time integration algorithm is used to update the a vector, and the static solution is obtained by the kinetic energy dissipation of virtual dynamic process which includes viscous damping type and kinetic damping type [16].DRM can use the virtual mass and kinetic damping, so as to avoid the situation that computational cost of EQSA depends on the small size of the model, and the viscous damping coefficient depends on the analysis experience.The application of DRM in the static analysis of structures began in the 50th and 60th of the last century [17][18][19], and its application scope gradually expanded in recent years [20][21][22], and recently it was often used in the analysis of tension structures [3,23,24].However, it is less common in the application and analysis of strainsoftening structures.
In the paper, kinetic damping DRM is used to calculate the static equilibrium path of strain-softening structure.Firstly, the paper gives the problem description and then describes the process of structural static problem based on DRM in two levels.And two examples and the corresponding results are described.Finally, the discussion and conclusions are given.

Description of the Problem
The problem to be solved is a boundary problem in static solid analysis.The condition for material softening is given by [25,26] where σ is the stress rate vector and ε is the strain rate vector.The material strain-softening makes the equation of the boundary problem no longer elliptic [27], and accordingly the solution lost well-posedness [28].That introduced time into the problem can make the equation become well-posed [29], but the original boundary problem changes to initial boundary value problem.In order to solve the original boundary problem, we can firstly divide the whole problem solving process into several incremental loads or displacement steps; secondly we attach masses to the nodes to form the diagonal mass matrix M and then use DRM to solve the incremental steps.The essence of DRM is that the original solid problem is turned into the one of well-posed damped vibration problem: the dissipation of mechanical energy caused by damping will make the quasi-static displacement solution of the vibration problem approaches the solution of the original problem as  → ∞.And at any time point , if the displacement vector satisfies the convergence criterion, then it is also the solution of the original static problem.
For the FE discrete model of the solid, all the degrees of freedom (DOFs) can be divided into those with loading and those without loading.Here, we define the DOFs whose nodes are either with load increment or with displace increment as loading DOFs and the remaining DOFs as nonloading DOFs.In the nonlinear FE static analysis of the structure, it is assumed that the boundary of the solid is sufficiently restrained such that it becomes geometrically stable system with no rigid body displacement.The set of nonlinear equilibrium equations to be solved corresponding to all DOFs is where  is the residual nodal load vector, f = { 1 , . . .,   } is the external nodal load vector, a = { 1 , . . .,   } is the nodal displacement vector, p = { 1 , . . .,   } is internal nodal load vector, and  is the total number of DOFs in the system.Let the index set I 1 = {1, . . ., }, the index set I 2 = { 1 , . . .,   } ⊂ I 1 , and the index set I 3 = { 1 , . . .,   } ⊂ I 1 , where I 1 denotes the total DOFs, I 2 denote the loading DOFs, I 3 denotes the nonloading DOFs, and I 1 = I 2 + I 3 .For an arbitrary vector k = {V 1 , . . ., V  } with  elements, define k(I 2 ) = {V 1 , . . ., V  } and k(I 3 ) = {V 1 , . . ., V  }.Also define the constitutive status of Gaussian points of elements as the set C ≜ {, , s} in which  is strain vector,  is stress vector, and s is status vector of stress-strain relation.The norm of the vector During loading increment, a is the unknown and f is the known; in this case, the elements in f(I 2 ) are all not zero, while f(I 3 ) = 0.During displacement increment, a(I 2 ) is the known and a(I 3 ) and f(I 2 ) are the unknown, while f(I 3 ) = 0.The vector p is computed according to p = ∫ Ω B   dΩ [1], where B is   ×  geometry matrix,  is   × 1 stress vector,   is the total number of elements in , and Ω is the space occupied by the solid.

Solution of the Problem
In this paper, the incremental method is used to solve the static problem (2), and DRM is used to solve the increment.In the solution process of DRM, there are several processes in which the kinetic damping is applied.And the method for applying kinetic damping in this paper is to set the nodal velocity vector corresponding to maximum kinetic energy of the system to zero.

Incremental Solution Method and Solution Process for the Static Problem.
Either load increment or displacement increment is used in the global solution of static problem (2) using incremental method.When load increment is used, the load increment Δf load = {Δ 1 , . . ., Δ  }; when displacement increment is used, the displacement increment Δa load = {Δ 1 , . . ., Δ  }, where Δa load is a reasonably assumed quantity such as the displacement mode of a structure observed during a strong earthquake.
The solution process for increment step  is as follows.
(1) At the beginning of the increment : (1) the known quantities are Δa  , and C −1 , where Δf is the nodal incremental load vector, Δa is the nodal incremental displacement vector, and the subscripts denote the sequence number of increment steps; here no superscripts are used because the quantities are the converged ones; (2) now either loading increment or displacement increment is used depending on the computation results of increment step  − 1 (refer to Section 3.3); a 0 , and C −1 = C 0 , where C 0 , f 0 , and a 0 represent the status of the discrete system before increment Step (1).
(3) The unknown given above are all solved by using DRM, so we get f  = f −1 + Δf  , a  = a −1 + Δa  , and C  .
(4) If the convergence criterion is satisfied, go to the solution process for increment step  + 1.

The Increment
Step Solution Method and Solution Process for DRM.Before solving the increment step by using DRM, the time domain is discretized with fixed interval ℎ to get the points ih ( = 0, 1, . ..) along the time axis.During the solution, the discrete time domain is sequentially divided into series {  } in which   ( = 1, 2, . ..) is th motion segment with application of kinetic damping.In detail,   includes an undamped motion process and a final time point at which the kinetic damping is applied.The beginning point of   is the ending point  −1 h of  −1 , and the ending point of   is the beginning point   ℎ of  +1 .Through the series   ( = 1, 2, . ..), the dynamic responses f   , a   , C   (the superscripts represent the time point sequence number) approach the static solution f  , a  , C  .
The control process based on dynamic response for solving   is as follows: (1) Determining the status quantities at the beginning time point of  −1 h of   .

𝑛
according to (7); thirdly, the vector (2) Using the status quantities (a , and f  −1  ) at the beginning time point  −1 ℎ of   as the initial status, the status quantities at  ≥  −1 can be obtained by solving the undamped motion equation in which M  is the diagonal mass matrix, and the method for determining M  and h is shown in Section 3.3.The status quantities a   , C   , p   , f   at ih time point ( =  −1 ,  −1 + 1, . ..) are solved by using explicit time integration algorithm.
(2-1) We compute a   and C   and the control based on a   .
/ℎ based on velocity and acceleration are used to construct explicit time integration algorithm to solve (3).The formulas for a   are [30] where at time points  − 3/2,  − 1/2, respectively, obtained by difference calculation from a   , a −1  , and a −2  that had been calculated satisfy the criterion then it is taken that the maximum kinetic energy of the system happens at time point (−1)ℎ, and, thereby, the kinetic damping is applied, and the nodal velocity vector at time point   ℎ is simultaneously set to zero ( ȧ    = 0), at the same time, the process   ends, which gives the status quantities at time point   ℎ: After that, the process  +1 begins.
(2-1-4) If a   satisfies the divergence criterion when loading increment is used, where  1 is a limiting value, then it is taken that f  ≡ f −1 +Δf  is near or surpasses the peak point of static equilibrium path, and the computation for increment step  is restarted and the displacement increment is used. ( in which the subscripts  for increment step sequential number and the superscripts  for time point sequential number are both omitted [5].In (7) where  2 is the limiting value, then the dynamic responses f   , a   , and C   at the time point  are taken as the solution of increment step ; that is, and increment step  + 1 is started.If displacement increment is used in increment step  and if (Δf  (I 2 ))  Δa  (I 2 ) > 0 at convergence, then the current point is at the hardening segment of the static equilibrium path, and therefore load increment will be used in increment step  + 1.

Parameter Setting.
The explicit time integration algorithm is conditionally stable.For discrete linear system, the condition for stability is ℎ ≤ 2/ mos max , where  mos max is the highest  circular eigenfrequency of the system.The computation of  mos max is costly [10].It is pointed out that [23] if (10) then the algorithm is stable, where  is nodal mass and  nos max is nodal linear displacement stiffness.The nonlinear discrete system can be approximated as linear discrete system within one incremental step.The DRM with kinetic damping solves for the static solution of increment step  through virtual vibration process.Before starting increment step , the value of ℎ can be chosen so that the nodal mass  can be adjusted for all nodes to satisfy the stability condition.We take where  3 is a coefficient greater than 1 and  don max is the maximum linear stiffness among the four displacement directions at each node in the directions of the orthogonal coordinate system in the 2-dimensional case.This stiffness is estimated by applying infinitesimal displacement node by node to calculate the corresponding nodal load.The term  min is lower bound value for the nodal mass, whose function is to avoid the formation of local negative stiffness associated with negative mass.Thereby the global diagonal mass matrix M  can be constructed.

Numerical Case Studies
4.1.Model I: Plane Truss.Model I is a plane truss constructed by 3 tension bars (Figure 1).The section area  of each bar is 1 mm 2 .Each tension bar is divided into 4 line elements, and linear polynomial is used as the shape function.By calculating the axial stress   of line element, p {} can be computed precisely by using direct force-balance method.And the axial stress of the element which is used to calculate   can be written as   = ( −  ini )/ ini , where  ini is the

Model
Δf load (kN) undeformed axial length of the element, and  is the deformed axial length.
The artificial one-dimensional stress-strain constitutive relation used in Model I is so constructed that it is characteristic of strain-softening and history dependence.And it is defined as (Figure 1) where  ,un and  ,un are maximum tension strain and maximum tension stress during the strain history of material, respectively.
In Figure 1 the loaded nodes and loaded degrees of freedom are also shown.Some of the analysis parameters are valued as ℎ = 0.1 s,  2 = 0.0002,  3 = 4.0, and  min = 0.3 kg, and the values of other analysis parameters are shown in Table 1.
If there is an element whose   > 1.4 in each of the three bars, then the computation is finished.

Model II: A Concrete Column Subjected to Axial Compression.
Model II is a concrete column subjected to axial compression (Figure 2), whose thickness  is 150 mm.The meshing is illustrated in Figure 2. The meshing-style in  direction is denoted as IIA, IIB, and IIC whose × are 1×400, 2 × 200, and 4 × 100, respectively, where  is the number of elements and  is the length of each element in  direction.The element of Model II is rectangular type with 4 nodes, and the quadratic polynomial is chosen as the shape function; p {} can be computed approximately by using 4 Gauss integration points.
A scalar damage constitutive model which was put forward by Mazars and Pijaudier-Cabot [31,32] is used in Model II.In the constitutive model, the material is supposed to behave elastically and to remain isotropic.The constitutive equation can be written as  = (1 − ) D ini :  (13) in which  is tress tensor,  is damage scalar, D ini is initial stiffness tensor,  is strain tensor, and the symbol ":" indicates the tensorial product contracted on two indices.D ini can be constructed by elasticity modulus  and Poisson ratio ].The damage scalar  can be calculated by the following:  in which  max is the maximum damage value during the strain history of material and the value of the expression  +  + +  −  − can be regarded as real-time damage value.In (14), the weights  + and  − are functions of the strain state;  + and  − are defined as uniaxial tension and compression damage value [31].The expressions of  + and  − are where  ini is the initial damage threshold,  + ,  + ,  − , and  − are characteristics of the material, and ε is called the notion of equivalent strain, whose definition and expression were given in [31].In Model II, the material parameters are as follows:  = 25500 N/mm 2 , ] = 0.2;  ini = 9.8 × 10 −5 ,  + = 0.95,  + = 11,500,  − = 1.38, and  − = 2,000.And by following these material parameters, we got the uniaxial compression strength   = 24.4N/mm 2 and corresponding strain   = 1.75×10 −3 .The illustration of the uniaxial compressive stressstrain relation is shown in Figure 2.
In Figure 2 the loaded nodes and loaded degrees of freedom are also shown.Some of the analysis parameters are chosen as ℎ = 0.1 s,  2 = 0.0002,  3 = 4.0, and  min = 0.3 kg, and the values of other analysis parameters are shown in Table 1.

Results and Discussion
For Model I, the load-deflection curve of node 1 in  direction is given in Figure 3, which shows not only the hardeningstages, softening-stages, and load peaks, but also the loadbreaks caused by softening of elements.The load-displacement curves of node 3 in  direction of Models IIA, IIB, and IIC are presented in Figure 4.These curves show the hardening-stage and softening-stage of Model II under the axial compression.It is shown in Figure 4 that the different meshes not only present the same load peak value but also give the same equilibrium path before this load peak.The displacement value and load value at peak point of the path are close to their estimated values    = 1.75 × 10 −3 × 400 = 0.7 mm, and    = [24.4 × (200 × 150)]/4 = 183 kN with the assumption that  1 =  2 =  3 =  4 =  5 and 2 1 =  2 =  3 =  4 = 2 5 .After the load peak, different static equilibrium paths are obtained, which is in the accordance with the nonuniqueness of static solution for strain-softening solid [1].
In Model I, an artificial, one-dimensional strain-softening constitutive model is used such that the static equilibrium path for a simple plane truss is solved.In spite of the simple structure, the solution process incorporates the concept associated with the structure of the universal first law of thermodynamics that is basis of the kinetic damping DRM for the solution of the increment steps.That is, it reflects the law of dissipation of kinetic energy and conversion among the external potential energy, internal potential energy, and kinetic energy.Therefore, the application of the method is in essence not limited by the constitutive relation of material and the type of elements.Thereby, further application of the method can be expected in which constitutive model can be calibrated against experimental data so that it is suitable for use in static analysis of engineering structure.
In Model II, the Mazars elastic damage constitutive model with one scalar damage variable is used.Although the resulting stress-strain relation is simple and smooth enough, it is relatively poor in reflecting the constitutive relation under multiaxial state of stress [33].In future studies, the elastic and/or plastic damage constitutive model with 2 scalar damage variables should be used in the static analysis of different concrete members in order to truly reflect the tensioncompression difference of concrete material under complicated condition [34].

Conclusions
(1) The process for the solution of static problems based on nonlinear finite element DRM method is clarified with three depths of global, increment step, and the motion process with application of kinetic damping.The method of DRM makes the original ill-posed static problem become well-posed and solvable by introducing the time domain for the virtual dynamic process such that the static solution is obtained by way of dissipation of mechanical energy.
(2) The global analysis including the ascending branch, the peak point, and the descending branch of the equilibrium path of the strain-softening structure is realized by application of kinetic damping to the motion process, by introducing the rule for conversion between load and displacement increment and by introducing the divergence criterion.
(3) Two instances of structures with strain-softening and the associated ill-posedness, one model with geometrical nonlinearity and the other model with mesh size sensitivity, are successfully solved using DRM.
then the status at this time point  −1 ℎ is as follows.The vectors ȧ  −1  = 0, a  −1  , and C  −1  are known; secondly, the vector p  −1  is derived from C  −1

{i}
Element number i (j) Node number j Load

Figure 1 :
Figure 1: Geometry, meshing, and illustration of stress-strain relation of Model I.

Figure 2 :
Figure 2: Geometry, meshing, and illustration of uniaxial compression stress-strain relation of Model II.

Figure 3 :Figure 4 :
Figure 3: Load-deflection curve of node 1 in  direction of Model I.
We compute f   and p   and the control based on f , e is the element number;   is the total number of elements; p {} , B {} ,  {} , and Ω {} are the internal nodal load, geometry matrix, stress vector, and occupied space of element , respectively; p {} , B {} , and  {} are  {}  × 1,   ×  {}  , and   × 1 matrix in which  {}  is the DOFs number of elements ; T {} is the transformation matrix for element  with the size  ×  {}  where  is the total DOFs of the system.

Table 1 :
Analysis parameters for Models I and II.