Advances in Materials Science and Engineering

Volume 2017, Article ID 8060987, 7 pages

https://doi.org/10.1155/2017/8060987

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

College of Civil Engineering, Tongji University, Shanghai 200092, China

Correspondence should be addressed to Wei Wang; moc.qq@0121iewgnaw

Received 31 March 2017; Accepted 8 May 2017; Published 3 July 2017

Academic Editor: Xiao-Yong Wang

Copyright © 2017 Wei Wang and Xiaozu Su. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

#### 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 matrix or the instability of node displacement [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 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 [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 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]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 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 iswhere is the residual nodal load vector, is the external nodal load vector, is the nodal displacement vector, is internal nodal load vector, and is the total number of DOFs in the system.

Let the index set , the index set , and the index set , where denotes the total DOFs, denote the loading DOFs, denotes the nonloading DOFs, and . For an arbitrary vector with elements, define and . Also define the constitutive status of Gaussian points of elements as the set in which is strain vector, is stress vector, and is status vector of stress-strain relation. The norm of the vector is defined as .

During loading increment, is the unknown and is the known; in this case, the elements in are all not zero, while . During displacement increment, is the known and and are the unknown, while . The vector is computed according to [1], where is geometry matrix, is stress vector, 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 ; when displacement increment is used, the displacement increment , where 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.

At the beginning of the increment : the known quantities are , , and , where is the nodal incremental load vector, 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; now either loading increment or displacement increment is used depending on the computation results of increment step (refer to Section 3.3); if , use load increment, , , and , where , , and represent the status of the discrete system before increment Step .

When load increment is used, the known quantities are , , and , and the unknown are and . When displacement increment is used, the known quantities are , , and , and the unknown are , , and .

The unknown given above are all solved by using DRM, so we get , , and .

If the convergence criterion is satisfied, go to the solution process for increment step .

##### 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 interval to get the points* ih* along the time axis. During the solution, the discrete time domain is sequentially divided into series in which 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 *h* of , and the ending point of is the beginning point of . Through the series , the dynamic responses , , (the superscripts represent the time point sequence number) approach the static solution , , .

The control process based on dynamic response for solving is as follows:

Determining the status quantities at the beginning time point of *h* of .

If , then ; the status at this time point is as follows. Firstly, the nodal velocity vector ; secondly, the vector if load increment is used, or if displacement increment is used, in which and ; thirdly, the status vector ; fourthly, the vector is derived from according to (7); fifthly, the vector if load increment is used, or and if displacement increment is used.

If , then the status at this time point is as follows. The vectors , , and are known; secondly, the vector is derived from according to (7); thirdly, the vector if load increment is used, or , if displacement increment is used.

Using the status quantities (, , , , and ) at the beginning time point of as the initial status, the status quantities at can be obtained by solving the undamped motion equationin which is the diagonal mass matrix, and the method for determining and* h* is shown in Section 3.3. The status quantities , , , at* ih* time point are solved by using explicit time integration algorithm.

We compute and and the control based on .

The central difference method formulas , based on velocity and acceleration are used to construct explicit time integration algorithm to solve (3). The formulas for are [30]where .

We compute from , , .

If the system kinetic energies , at time points , , respectively, obtained by difference calculation from , , and that had been calculated satisfy the criterion then it is taken that the maximum kinetic energy of the system happens at time point , and, thereby, the kinetic damping is applied, and the nodal velocity vector at time point is simultaneously set to zero (), at the same time, the process ends, which gives the status quantities at time point : , , . After that, the process begins.

If satisfies the divergence criterion when loading increment is used, where is a limiting value, then it is taken that 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.

We compute and and the control based on and .

(2-2-1) If load increment is used, ; if displacement is used, and ; the latter guarantees that (see (4)).

(2-2-2) The vectorin which the subscripts for increment step sequential number and the superscripts for time point sequential number are both omitted [5]. In (7),* e* is the element number; is the total number of elements; , , , and are the internal nodal load, geometry matrix, stress vector, and occupied space of element , respectively; , , and are , , and matrix in which is the DOFs number of elements ; is the transformation matrix for element with the size where is the total DOFs of the system.

If and satisfy convergence criterionwhere is the limiting value, then the dynamic responses , , and at the time point are taken as the solution of increment step ; that is, and increment step is started. If displacement increment is used in increment step and if 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 .

If the computed and do not satisfy the convergence criterion, then , and repeat steps “ and ”.

##### 3.3. Parameter Setting

The explicit time integration algorithm is conditionally stable. For discrete linear system, the condition for stability is , where is the highest circular eigenfrequency of the system. The computation of is costly [10]. It is pointed out that [23] if then the algorithm is stable, where is nodal mass and 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 takewhere is a coefficient greater than 1 and 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 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 can be constructed.

#### 4. 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, 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 , where is the undeformed axial length of the element, and is the deformed axial length.