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 mm2. 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.
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 and 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 s, , , and kg, and the values of other analysis parameters are shown in Table 1.
If there is an element whose 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 (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 , , and , 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; 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 asin which is tress tensor, is damage scalar, is initial stiffness tensor, is strain tensor, and the symbol “” indicates the tensorial product contracted on two indices. can be constructed by elasticity modulus and Poisson ratio . The damage scalar can be calculated by the following:in which 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 arewhere 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: N/mm2, ; , , = 11,500, , and = 2,000. And by following these material parameters, we got the uniaxial compression strength N/mm2 and corresponding strain . The illustration of the uniaxial compressive stress-strain 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 s, , , and kg, and the values of other analysis parameters are shown in Table 1.
When the expression is satisfied, then the computation is finished.
5. 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 hardening-stages, softening-stages, and load peaks, but also the load-breaks 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 equilibrium path are close to their estimated values mm, and kN with the assumption that and . 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 tension-compression difference of concrete material under complicated condition [34].
6. Conclusions
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.
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.
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.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was financially supported by the National Natural Science Foundation of China (51178328).