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

Authors: Wei Wang; Xiaozu Su
Journal: Advances in Materials Science and Engineering
 (2017)

---

## Abstract

In order to deal with the divergence and instability due to the ill-posedness of the nonlinear finite element (FE) model of strain-softening 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.


---

## Body

## 1. 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 strain-softening material widely used in civil engineering.Implicit static algorithm (ISA) [2], explicit quasi-static algorithm (EQSA) [2], and dynamic relaxation method (DRM) [3] are the three main methods used for structural static analysis.According to the increment control type, ISA can be divided into load-control, displacement control, and arc-length 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 matrixK 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 displace-control 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 snap-back 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 arc-length 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 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 thea 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–19], and its application scope gradually expanded in recent years [20–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 strain-softening 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.


## 2. 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](1)ε˙Tσ˙≤0,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 t→∞. And at any time point t, 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(2)ψf,a=f-pa=0,where ψ is the residual nodal load vector, f={f1,…,fr} is the external nodal load vector, a={a1,…,ar} is the nodal displacement vector, p={p1,…,pr} is internal nodal load vector, and r is the total number of DOFs in the system.Let the index setI1={1,…,r}, the index set I2={s1,…,sr}⊂I1, and the index set I3={t1,…,tr}⊂I1, where I1 denotes the total DOFs, I2 denote the loading DOFs, I3 denotes the nonloading DOFs, and I1=I2+I3. For an arbitrary vector v={v1,…,vr} with r elements, define v(I2)={vs1,…,vsr} and vI3={vt1,…,vtr}. 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 x={x1,x2,…,xn} is defined as x2=max1≤i≤nxi.During loading increment,a is the unknown and f is the known; in this case, the elements in f(I2) are all not zero, while f(I3)=0. During displacement increment, a(I2) is the known and a(I3) and f(I2) are the unknown, while f(I3)=0. The vector p is computed according to p=∫ΩBTσdΩ [1], where B is nσ×r geometry matrix, σ is nσ×1 stress vector, nσ is the total number of elements in σ, and Ω is the space occupied by the solid.


## 3. 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.


### 3.1. 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 Δfload={Δfs1,…,Δfsr}; when displacement increment is used, the displacement increment Δaload={Δas1,…,Δasr}, where Δaload is a reasonably assumed quantity such as the displacement mode of a structure observed during a strong earthquake.The solution process for increment stepn is as follows.(
1
) At the beginning of the increment n: (1) the known quantities are fn-1=fn-2+Δfn-1=∑i=1n-1Δfi, an-1=an-2+Δan-1=∑i=1n-1Δai, and Cn-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 n-1 (refer to Section 3.3); (3) if n=1, use load increment, fn-1=f0, an-1=a0, and Cn-1=C0, where C0, f0, and a0 represent the status of the discrete system before increment Step  (1).(
2
) When load increment is used, the known quantities are ΔfnI2=Δfload, Δfn(I3)=0, and Cn-1, and the unknown are Δan and Cn. When displacement increment is used, the known quantities are Δan(I2)=Δaload, Δfn(I3)=0, and Cn-1, and the unknown are Δan(I3), Δfn(I2), and Cn.(
3
) The unknown given above are all solved by using DRM, so we get fn=fn-1+Δfn, an=an-1+Δan, and Cn.(
4
) If the convergence criterion is satisfied, go to the solution process for increment step n+1.


### 3.2. 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 intervalh to get the pointsih  (i=0,1,…) along the time axis. During the solution, the discrete time domain is sequentially divided into series {Dk} in which Dk(k=1,2,…) is kth motion segment with application of kinetic damping. In detail, Dk includes an undamped motion process and a final time point at which the kinetic damping is applied. The beginning point of Dk is the ending point jk-1h of Dk-1, and the ending point of Dk is the beginning point jkh of Dk+1. Through the series Dk(k=1,2,…), the dynamic responses fni, ani, Cni (the superscripts represent the time point sequence number) approach the static solution fn, an, Cn.The control process based on dynamic response for solvingDk is as follows:(
1
) Determining the status quantities at the beginning time point of jk-1h of Dk.(
1-1
) If k=1, then j0h=0; the status at this time point is as follows. Firstly, the nodal velocity vector a˙n0=0; secondly, the vector an0=an-1 if load increment is used, or an0=an-1+Δan0 if displacement increment is used, in which Δan0(I2)=Δaload and Δan0(I3)=0; thirdly, the status vector Cn0=Cn-1; fourthly, the vector pn0 is derived from Cn0 according to (7); fifthly, the vector fn0≡fn-1+Δfn if load increment is used, or fn0(I2)=pn0(I2) and fn0(I3)=0 if displacement increment is used.(
1
-
2
) If k>1, then the status at this time point jk-1h is as follows. The vectors a˙njk-1=0, anjk-1, and Cnjk-1 are known; secondly, the vector pnjk-1 is derived from Cnjk-1 according to (7); thirdly, the vector fnjk-1≡fn-1+Δfn if load increment is used, or fnjk-1(I2)=pnjk-1(I2), fnjk-1I3=0 if displacement increment is used.(
2
) Using the status quantities (anjk-1, a˙njk-1, Cnjk-1, pnjk-1, and fnjk-1) at the beginning time point jk-1h of Dk as the initial status, the status quantities at t≥jk-1 can be obtained by solving the undamped motion equation(3)Mna¨nt=fnt-pntt≥jk-1hin which Mn is the diagonal mass matrix, and the method for determining Mn andh is shown in Section 3.3. The status quantities ani, Cni, pni, fni atih time point (i=jk-1,jk-1+1,…) are solved by using explicit time integration algorithm.(
2
-
1
) We compute ani and Cni and the control based on ani.(
2-1-1
) The central difference method formulas a˙ni-3/2=(ani-1-ani-2)/h, a¨ni-1=a˙ni-1/2-a˙ni-3/2/h based on velocity and acceleration are used to construct explicit time integration algorithm to solve (3). The formulas for ani are [30](4)ani=ani-1+h2bni2i=jk-1+1ani=2ani-1-ani-2+h2bnii>jk-1+1,where bni=Mn-1(fni-1-pni-1).(
2-1-2
) We compute Cni from ani, ani-1, Cni-1.(
2-1-3
) If the system kinetic energies Eni-3/2=0.5(a˙ni-3/2)TMna˙ni-3/2, Eni-1/2=0.5(a˙ni-1/2)TMna˙ni-1/2 at time points i-3/2, i-1/2, respectively, obtained by difference calculation from ani, ani-1, and ani-2 that had been calculated satisfy the criterion (5)Eni-1/2<Eni-3/2;then it is taken that the maximum kinetic energy of the system happens at time point (i-1)h, and, thereby, the kinetic damping is applied, and the nodal velocity vector at time point jkh is simultaneously set to zero (a˙njk=0), at the same time, the process Dk ends, which gives the status quantities at time point jkh: anjk=ani-1, a˙njk=0, Cnjk=Cni-1. After that, the process Dk+1 begins.(
2-1-4
) If ani satisfies the divergence criterion when loading increment is used, (6)aniI22≥c1,where c1 is a limiting value, then it is taken that fn≡fn-1+Δfn is near or surpasses the peak point of static equilibrium path, and the computation for increment step n is restarted and the displacement increment is used.(
2-2
) We compute fni and pni and the control based on fni and pni.(2-2-1) If load increment is used,fni≡fn-1+Δfn; if displacement is used, fni(I3)=0 and fni(I2)=pni(I2); the latter guarantees that Δani+1(I2)=Δaload (see (4)).(2-2-2) The vector(7)p=∑e=1neTepe=∑e=1neTe∫ΩeBeTσedΩein which the subscripts n for increment step sequential number and the superscripts i for time point sequential number are both omitted [5]. In (7),e is the element number; ne is the total number of elements; p{e}, B{e}, σ{e}, and Ω{e} are the internal nodal load, geometry matrix, stress vector, and occupied space of element e, respectively; p{e}, B{e}, and σ{e} are nf{e}×1, nσ×nf{e}, and nσ×1 matrix in which nf{e} is the DOFs number of elements e; T{e} is the transformation matrix for element e with the size r×nf{e} where r is the total DOFs of the system.(
2-2-3
) If fni and pni satisfy convergence criterion(8)ψni2fni-fn-12≤c2,where c2 is the limiting value, then the dynamic responses fni, ani, and Cni at the time point i are taken as the solution of increment step n; that is, (9)fn=fni,an=ani,Cn=Cni,and increment step n+1 is started. If displacement increment is used in increment step n and if ΔfnI2TΔan(I2)>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 n+1.(
2-2-4
) If the computed fni and pni do not satisfy the convergence criterion, then i←i+1, and repeat steps “(2-1) and (2-2)”.


### 3.3. Parameter Setting

The explicit time integration algorithm is conditionally stable. For discrete linear system, the condition for stability ish≤2/ωmaxmos, where ωmaxmos is the highest circular eigenfrequency of the system. The computation of ωmaxmos is costly [10]. It is pointed out that [23] if (10)h≤2mkmaxnosthen the algorithm is stable, where m is nodal mass and kmaxnos 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 n through virtual vibration process. Before starting increment step n, the value of h can be chosen so that the nodal mass m can be adjusted for all nodes to satisfy the stability condition. We take(11)m=c3kmaxdonh22≥mmin,where c3 is a coefficient greater than 1 and kmaxdon 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 mmin 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 Mn can be constructed.


## 3.1. 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 Δfload={Δfs1,…,Δfsr}; when displacement increment is used, the displacement increment Δaload={Δas1,…,Δasr}, where Δaload is a reasonably assumed quantity such as the displacement mode of a structure observed during a strong earthquake.The solution process for increment stepn is as follows.(
1
) At the beginning of the increment n: (1) the known quantities are fn-1=fn-2+Δfn-1=∑i=1n-1Δfi, an-1=an-2+Δan-1=∑i=1n-1Δai, and Cn-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 n-1 (refer to Section 3.3); (3) if n=1, use load increment, fn-1=f0, an-1=a0, and Cn-1=C0, where C0, f0, and a0 represent the status of the discrete system before increment Step  (1).(
2
) When load increment is used, the known quantities are ΔfnI2=Δfload, Δfn(I3)=0, and Cn-1, and the unknown are Δan and Cn. When displacement increment is used, the known quantities are Δan(I2)=Δaload, Δfn(I3)=0, and Cn-1, and the unknown are Δan(I3), Δfn(I2), and Cn.(
3
) The unknown given above are all solved by using DRM, so we get fn=fn-1+Δfn, an=an-1+Δan, and Cn.(
4
) If the convergence criterion is satisfied, go to the solution process for increment step n+1.


## 3.2. 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 intervalh to get the pointsih  (i=0,1,…) along the time axis. During the solution, the discrete time domain is sequentially divided into series {Dk} in which Dk(k=1,2,…) is kth motion segment with application of kinetic damping. In detail, Dk includes an undamped motion process and a final time point at which the kinetic damping is applied. The beginning point of Dk is the ending point jk-1h of Dk-1, and the ending point of Dk is the beginning point jkh of Dk+1. Through the series Dk(k=1,2,…), the dynamic responses fni, ani, Cni (the superscripts represent the time point sequence number) approach the static solution fn, an, Cn.The control process based on dynamic response for solvingDk is as follows:(
1
) Determining the status quantities at the beginning time point of jk-1h of Dk.(
1-1
) If k=1, then j0h=0; the status at this time point is as follows. Firstly, the nodal velocity vector a˙n0=0; secondly, the vector an0=an-1 if load increment is used, or an0=an-1+Δan0 if displacement increment is used, in which Δan0(I2)=Δaload and Δan0(I3)=0; thirdly, the status vector Cn0=Cn-1; fourthly, the vector pn0 is derived from Cn0 according to (7); fifthly, the vector fn0≡fn-1+Δfn if load increment is used, or fn0(I2)=pn0(I2) and fn0(I3)=0 if displacement increment is used.(
1
-
2
) If k>1, then the status at this time point jk-1h is as follows. The vectors a˙njk-1=0, anjk-1, and Cnjk-1 are known; secondly, the vector pnjk-1 is derived from Cnjk-1 according to (7); thirdly, the vector fnjk-1≡fn-1+Δfn if load increment is used, or fnjk-1(I2)=pnjk-1(I2), fnjk-1I3=0 if displacement increment is used.(
2
) Using the status quantities (anjk-1, a˙njk-1, Cnjk-1, pnjk-1, and fnjk-1) at the beginning time point jk-1h of Dk as the initial status, the status quantities at t≥jk-1 can be obtained by solving the undamped motion equation(3)Mna¨nt=fnt-pntt≥jk-1hin which Mn is the diagonal mass matrix, and the method for determining Mn andh is shown in Section 3.3. The status quantities ani, Cni, pni, fni atih time point (i=jk-1,jk-1+1,…) are solved by using explicit time integration algorithm.(
2
-
1
) We compute ani and Cni and the control based on ani.(
2-1-1
) The central difference method formulas a˙ni-3/2=(ani-1-ani-2)/h, a¨ni-1=a˙ni-1/2-a˙ni-3/2/h based on velocity and acceleration are used to construct explicit time integration algorithm to solve (3). The formulas for ani are [30](4)ani=ani-1+h2bni2i=jk-1+1ani=2ani-1-ani-2+h2bnii>jk-1+1,where bni=Mn-1(fni-1-pni-1).(
2-1-2
) We compute Cni from ani, ani-1, Cni-1.(
2-1-3
) If the system kinetic energies Eni-3/2=0.5(a˙ni-3/2)TMna˙ni-3/2, Eni-1/2=0.5(a˙ni-1/2)TMna˙ni-1/2 at time points i-3/2, i-1/2, respectively, obtained by difference calculation from ani, ani-1, and ani-2 that had been calculated satisfy the criterion (5)Eni-1/2<Eni-3/2;then it is taken that the maximum kinetic energy of the system happens at time point (i-1)h, and, thereby, the kinetic damping is applied, and the nodal velocity vector at time point jkh is simultaneously set to zero (a˙njk=0), at the same time, the process Dk ends, which gives the status quantities at time point jkh: anjk=ani-1, a˙njk=0, Cnjk=Cni-1. After that, the process Dk+1 begins.(
2-1-4
) If ani satisfies the divergence criterion when loading increment is used, (6)aniI22≥c1,where c1 is a limiting value, then it is taken that fn≡fn-1+Δfn is near or surpasses the peak point of static equilibrium path, and the computation for increment step n is restarted and the displacement increment is used.(
2-2
) We compute fni and pni and the control based on fni and pni.(2-2-1) If load increment is used,fni≡fn-1+Δfn; if displacement is used, fni(I3)=0 and fni(I2)=pni(I2); the latter guarantees that Δani+1(I2)=Δaload (see (4)).(2-2-2) The vector(7)p=∑e=1neTepe=∑e=1neTe∫ΩeBeTσedΩein which the subscripts n for increment step sequential number and the superscripts i for time point sequential number are both omitted [5]. In (7),e is the element number; ne is the total number of elements; p{e}, B{e}, σ{e}, and Ω{e} are the internal nodal load, geometry matrix, stress vector, and occupied space of element e, respectively; p{e}, B{e}, and σ{e} are nf{e}×1, nσ×nf{e}, and nσ×1 matrix in which nf{e} is the DOFs number of elements e; T{e} is the transformation matrix for element e with the size r×nf{e} where r is the total DOFs of the system.(
2-2-3
) If fni and pni satisfy convergence criterion(8)ψni2fni-fn-12≤c2,where c2 is the limiting value, then the dynamic responses fni, ani, and Cni at the time point i are taken as the solution of increment step n; that is, (9)fn=fni,an=ani,Cn=Cni,and increment step n+1 is started. If displacement increment is used in increment step n and if ΔfnI2TΔan(I2)>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 n+1.(
2-2-4
) If the computed fni and pni do not satisfy the convergence criterion, then i←i+1, and repeat steps “(2-1) and (2-2)”.


## 3.3. Parameter Setting

The explicit time integration algorithm is conditionally stable. For discrete linear system, the condition for stability ish≤2/ωmaxmos, where ωmaxmos is the highest circular eigenfrequency of the system. The computation of ωmaxmos is costly [10]. It is pointed out that [23] if (10)h≤2mkmaxnosthen the algorithm is stable, where m is nodal mass and kmaxnos 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 n through virtual vibration process. Before starting increment step n, the value of h can be chosen so that the nodal mass m can be adjusted for all nodes to satisfy the stability condition. We take(11)m=c3kmaxdonh22≥mmin,where c3 is a coefficient greater than 1 and kmaxdon 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 mmin 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 Mn can be constructed.


## 4. Numerical Case Studies

### 4.1. Model I: Plane Truss

Model I is a plane truss constructed by 3 tension bars (Figure1). The section area A of each bar is 1 mm2. Each tension bar is divided into 4 line elements, and linear polynomial is used as the shape function. By calculating the axial stress σl of line element, p{e} can be computed precisely by using direct force-balance method. And the axial stress of the element which is used to calculate σl can be written as εl=(l-lini)/lini, where lini is the undeformed axial length of the element, and l is the deformed axial length.Figure 1
Geometry, meshing, and illustration of stress-strain relation of Model I.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 (Figure1)(12)σl=7.72εle-3.86εl2εl≥εl,unεlσl,unεl,un0<εl<εl,un0εl≤0,where εl,un and σl,un are maximum tension strain and maximum tension stress during the strain history of material, respectively.In Figure1 the loaded nodes and loaded degrees of freedom are also shown. Some of the analysis parameters are valued as h=0.1 s, c2=0.0002, c3=4.0, and mmin=0.3 kg, and the values of other analysis parameters are shown in Table 1.Table 1
Analysis parameters for Models I and II.
Model
Δ
f
l
o
a
d (kN)
Δ
a
l
o
a
d (mm)
c
1 (mm)
I
Δ
f
1
y
=
0.1
×
10
-
3
Δ
a
1
y
=
10
3000
II
Δ
f
1
y
=
Δ
f
5
y
=
5
Δ
a
1
y
=
Δ
a
5
y
=
0.02
20
Δ
f
2
y
=
Δ
f
3
y
=
Δ
f
4
y
=
10
Δ
a
2
y
=
Δ
a
3
y
=
Δ
a
4
y
=
0.02If there is an element whoseεl>1.4 in each of the three bars, then the computation is finished.


### 4.2. Model II: A Concrete Column Subjected to Axial Compression

Model II is a concrete column subjected to axial compression (Figure2), whose thickness t is 150 mm. The meshing is illustrated in Figure 2. The meshing-style in y direction is denoted as IIA, IIB, and IIC whose n×l are 1×400, 2×200, and 4×100, respectively, where n is the number of elements and l is the length of each element in y direction. The element of Model II is rectangular type with 4 nodes, and the quadratic polynomial is chosen as the shape function; p{e} can be computed approximately by using 4 Gauss integration points.Figure 2
Geometry, meshing, and illustration of uniaxial compression stress-strain relation of Model II.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(13)σ=1-ωDini:εin which σ is tress tensor, ω is damage scalar, Dini is initial stiffness tensor, ε is strain tensor, and the symbol “:” indicates the tensorial product contracted on two indices. Dini can be constructed by elasticity modulus E and Poisson ratio ν. The damage scalar ω can be calculated by the following:(14)ω=max⁡ωmax,α+ω++α-ω-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(15)ω±=0ε~≤Kini1-1-A±Kiniε~-A±exp⁡B±ε~-Kiniε~>Kini,where Kini is the initial damage threshold, A+, B+, A-, and B- 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: E=25500 N/mm2, ν=0.2; Kini=9.8×10-5, A+=0.95, B+ = 11,500, A-=1.38, and B- = 2,000. And by following these material parameters, we got the uniaxial compression strength fc=24.4 N/mm2 and corresponding strain εc=1.75×10-3. The illustration of the uniaxial compressive stress-strain relation is shown in Figure 2.In Figure2 the loaded nodes and loaded degrees of freedom are also shown. Some of the analysis parameters are chosen as h=0.1 s, c2=0.0002, c3=4.0, and mmin=0.3 kg, and the values of other analysis parameters are shown in Table 1.When the expressionfn,3y<0.5‖{f1,3y,f2,3y,…,fn,3y}‖ is satisfied, then the computation is finished.


## 4.1. Model I: Plane Truss

Model I is a plane truss constructed by 3 tension bars (Figure1). The section area A of each bar is 1 mm2. Each tension bar is divided into 4 line elements, and linear polynomial is used as the shape function. By calculating the axial stress σl of line element, p{e} can be computed precisely by using direct force-balance method. And the axial stress of the element which is used to calculate σl can be written as εl=(l-lini)/lini, where lini is the undeformed axial length of the element, and l is the deformed axial length.Figure 1
Geometry, meshing, and illustration of stress-strain relation of Model I.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 (Figure1)(12)σl=7.72εle-3.86εl2εl≥εl,unεlσl,unεl,un0<εl<εl,un0εl≤0,where εl,un and σl,un are maximum tension strain and maximum tension stress during the strain history of material, respectively.In Figure1 the loaded nodes and loaded degrees of freedom are also shown. Some of the analysis parameters are valued as h=0.1 s, c2=0.0002, c3=4.0, and mmin=0.3 kg, and the values of other analysis parameters are shown in Table 1.Table 1
Analysis parameters for Models I and II.
Model
Δ
f
l
o
a
d (kN)
Δ
a
l
o
a
d (mm)
c
1 (mm)
I
Δ
f
1
y
=
0.1
×
10
-
3
Δ
a
1
y
=
10
3000
II
Δ
f
1
y
=
Δ
f
5
y
=
5
Δ
a
1
y
=
Δ
a
5
y
=
0.02
20
Δ
f
2
y
=
Δ
f
3
y
=
Δ
f
4
y
=
10
Δ
a
2
y
=
Δ
a
3
y
=
Δ
a
4
y
=
0.02If there is an element whoseεl>1.4 in each of the three bars, then the computation is finished.


## 4.2. Model II: A Concrete Column Subjected to Axial Compression

Model II is a concrete column subjected to axial compression (Figure2), whose thickness t is 150 mm. The meshing is illustrated in Figure 2. The meshing-style in y direction is denoted as IIA, IIB, and IIC whose n×l are 1×400, 2×200, and 4×100, respectively, where n is the number of elements and l is the length of each element in y direction. The element of Model II is rectangular type with 4 nodes, and the quadratic polynomial is chosen as the shape function; p{e} can be computed approximately by using 4 Gauss integration points.Figure 2
Geometry, meshing, and illustration of uniaxial compression stress-strain relation of Model II.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(13)σ=1-ωDini:εin which σ is tress tensor, ω is damage scalar, Dini is initial stiffness tensor, ε is strain tensor, and the symbol “:” indicates the tensorial product contracted on two indices. Dini can be constructed by elasticity modulus E and Poisson ratio ν. The damage scalar ω can be calculated by the following:(14)ω=max⁡ωmax,α+ω++α-ω-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(15)ω±=0ε~≤Kini1-1-A±Kiniε~-A±exp⁡B±ε~-Kiniε~>Kini,where Kini is the initial damage threshold, A+, B+, A-, and B- 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: E=25500 N/mm2, ν=0.2; Kini=9.8×10-5, A+=0.95, B+ = 11,500, A-=1.38, and B- = 2,000. And by following these material parameters, we got the uniaxial compression strength fc=24.4 N/mm2 and corresponding strain εc=1.75×10-3. The illustration of the uniaxial compressive stress-strain relation is shown in Figure 2.In Figure2 the loaded nodes and loaded degrees of freedom are also shown. Some of the analysis parameters are chosen as h=0.1 s, c2=0.0002, c3=4.0, and mmin=0.3 kg, and the values of other analysis parameters are shown in Table 1.When the expressionfn,3y<0.5‖{f1,3y,f2,3y,…,fn,3y}‖ is satisfied, then the computation is finished.


## 5. Results and Discussion

For Model I, the load-deflection curve of node 1 iny direction is given in Figure 3, which shows not only the hardening-stages, softening-stages, and load peaks, but also the load-breaks caused by softening of elements.Figure 3
Load-deflection curve of node 1 iny direction of Model I.The load-displacement curves of node 3 iny 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 equilibrium path are close to their estimated values εcl=1.75×10-3×400=0.7 mm, and fcA=[24.4×(200×150)]/4=183 kN with the assumption that a1y=a2y=a3y=a4y=a5y and 2f1y=f2y=f3y=f4y=2f5y. 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].Figure 4
Load-deflection curve of node 3 iny direction of Model II.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 tension-compression difference of concrete material under complicated condition [34].


## 6. 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 increment 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.



---

*Source: 8060987-2017-07-03.xml*