Abstract
This paper is a study of the numerical implementation of the spatial elastoplastic damage model of concrete by isogeometric analysis (IGA) method from three perspectives: the geometric modeling and the numerical formulation via IGA method, the constitutive model of concrete, and the solution algorithms for the local and global problems. The plasticity of concrete is considered on the basis of a nonassociated flow rule, where a three-parameter Barcelona yield surface and a modified Drucker-Prager plastic potential are used. The damage evolution of concrete driven by the internal variables is expressed by a piecewise function. In the study, the return-mapping algorithm and the substepping strategy are used for stress updating, and a new dissipation-based arc-length method with constraint path that considers the combined contribution of plasticity and damage to the energy dissipation is employed to trace the equilibrium path. After comparisons between simulation results and experimental data, the use of the elastoplastic damage model in the framework of IGA approach is proven to be practical in reflecting material properties of concrete.
1. Introduction
Armed with the progress of computer science and the in-depth study of computational theory, the numerical simulation technology of concrete structures has been developing towards highly advanced analysis that requires comprehensive capture of mechanical responses, especially the postcrack performances in three-dimensional (3D) space.
In the traditional finite element modeling, computer-aided design (CAD) and computer-aided engineering (CAE) are constructed on two different platforms based on one-way information transmission which is not only complex but also time-consuming. Therefore, Hughes et al. [1, 2] integrated CAD and CAE and proposed the philosophy of IGA, which directly uses the geometric languages in CAD instead of the commonly used Lagrange functions to construct the analysis domain.
The modeling of concrete mechanism is a crux in the advanced analysis due to high material nonlinearity. In continuum mechanics, concrete behavior can be simulated by elastic, elastic-plastic, or elastoplastic damage models. But the elastic damage model [3, 4] cannot consider the unrecoverable deformation of concrete, and the elastic-plastic model [5, 6] ignores the influence of damage. In comparison, the elastoplastic damage model, a combination of the first two models, is more reasonable because it can truly reflect concrete properties and therefore prevails in ductile failure [6–8], strong coupling [9, 10], or effective stress approaches [3, 11–16], among which the effective stress approach is the most widely used.
The solution algorithm is also a highly nonlinear issue that can be separated into the solutions of local and global problems to be solved at the level of Gaussian point as well as the level of structure, respectively. The difficulty in solving local problem exists in updating concrete stress. To facilitate convergence, Pérez-Foguet et al. [17] presented substepping algorithm combined with return-mapping algorithm. For global problem, the arc-length method is better than the traditional force control or displacement control methods. The classical arc-length method was proposed by Riks [18] and developed by Ramm [19] and Crisfield [20]. But the method was poor in solving material nonlinear problems. Gutiérrez [21], Verhoose et al. [22], and van der Meer et al. [23] established the constraint path of arc-length on the amount of the released energy and proved convergent capability in brittle and ductile analysis.
At present, the use of IGA has been widely explored in geometric modeling and mechanical engineering, but the extension use of it to structural analysis is rare. In this study, the knowledge for the numerical implementation of spatial elastoplastic damage model of concrete via IGA is introduced in detail from three aspects, that is, the geometric modeling and the numerical formulation of IGA, the constitutive model of concrete, and the solution algorithms for the local and global problems.
2. Basic Theory of IGA Solid Modeling
2.1. -Spline Basis Function and NURBS Basis Function
The basic analyzed object of IGA [24] in 1D, 2D, and 3D space is the NURBS (Nonuniform Rational -Spline), NURBS surface, and NURBS solid, respectively. A th-order -spline can be constructed by -spline basis functions defined on the knot sequence and their corresponding control points :
The coordinates of knots should be monotonically nondecreasing (), with two nonoverlapped neighboring knots constructing an element interval.
For , is defined as
For , can be calculated according to the Cox-de Boor recursion formula:
A th-order NURBS basis function of is defined by giving each a proper weight, and a th-order NURBS is constructed by th-order NURBS basis functions and their corresponding control points :where is the weight for . NURBS basis function degenerates to -spline basis function when all the weights take the same value.
The NURBS solid is a tensor product of , , and :where , , and refer to the number of -spline basis functions of , , and in the directions of , , and . is the control point of and is the weight.
2.2. Numerical Formulation for NURBS Solid Model
The numerical integral of a given function (i.e., stiffness matrix ) on a domain consists of NURBS solids that can be achieved via the philosophy of isoparametric transformation. According to the coordinate systems where the NURBS solids locate in the process of transformation, three spaces, that is, the parent space , the parametric space , and the physical space , are defined by the local, knot, and Cartesian coordinates, respectively. Hence,
As shown in Figure 1, we define the mappings of parametric space parent space and physical space parametric space as mapping 1 and mapping 2, respectively. For a NURBS solid, if the coordinates of analyzed point in the parent, parametric, and physical spaces are , , and , respectively, then the Jacobian matrices for mapping 1 and mapping 2 arewhere is the total number of control points for a single NURBS solid and is the Cartesian coordinates of control point .
The integral of function on can be derived aswith determinant .
2.3. Comparison with Traditional FEM
Although the numerical formulations of IGA and FE models are both constructed on the philosophy of isoparametric transformation, the basic ideas of these two are intrinsically different. The FE model, whose basis functions are chosen as Lagrange functions, is only an approximation to the original geometric model. The IGA model, on the other hand, with its element bases directly adopting the NURBS basis functions that exactly describe the geometric model, has no loss of information between the analyzed model and the original one, providing an extremely apparent advantage for curved or complicated shapes, as shown in Figure 2.
(a) FEM modeling
(b) IGA modeling
Considering a cubic as shown in Table 1 and discretizing the cubic to th-order elements, respectively, via FEM and IGA methods, we can calculate the total number of nodes in FE model (NFEM) and the total number of control points in IGA model (NIGA) as functions of and . If linear basis functions are used, then NFEM = NIGA. Once the order of basis functions is elevated, as seen in Figure 3, NIGA is apparently smaller than NFEM for the same value of , indicating that the scale of the equilibrium equation of IGA model is significantly reduced. It is concluded, as for Figure 4, that NIGA is not sensitive to the order changes of basis functions. Even though quartic NURBS basis functions are used, the scale of the equilibrium equation of IGA model is smaller than that of quadratic FE model. It should also be mentioned that the conduction of mesh refinement in IGA model is more convenient than that in FE model, since order elevation and knot insertion can be directly carried out on its linear control model.
(a)
(b)
3. Constitutive Model of Concrete
3.1. Basic Functions
The basic functions for elastoplastic damage model arewhere , , and are the total, the elastic, and the plastic strains, respectively; and are the nominal and the effective stresses, respectively; is the damage; is the initial elastic stiffness matrix; is the flow vector that can be derived from the plastic potential ; is the plastic modulus; is the set of internal variables; is the plastic multiplier that satisfies the Kuhn-Tucker conditionwhere is the yield function.
Two internal variables of and are introduced to record the hardening/softening behaviors of concrete under tension and compression, respectively; that is, , where and are the tensile and compressive equivalent plastic strains, respectively.
3.2. Yield Function
The yield surface of concrete is described by the modified Barcelona functionwithwhere and are the first stress invariant and the second deviatoric stress invariant computed by the effective stress tensor , respectively; is algebraically the maximum principal effective stress with , the Macaulay bracket function; with and being the uniaxial and biaxial yield strengths of concrete in compression, respectively; and are the effective tensile and compressive cohesion stresses of concrete, respectively; is introduced only when concrete is subjected to triaxial compression with .
3.3. Plastic Potential
A hyperbolic form of the Drucker-Prager function is used as the plastic potential of concrete:where is the dilation angle; is the uniaxial tensile yield stress; is the eccentricity parameter.
3.4. Plastic Modulus
The plastic modulus is written aswith being the weight function which can be computed as
3.5. Damage Law
The total damage can be seen as the combination of tensile and compressive damage, which has been discussed in detail by Lubliner et al. [25] and Lee and Fenves [13] and then furtherly developed in ABAQUS [26]:withwhere and are the tensile and compressive damage variables, respectively; and are the stiffness recovery factors reflecting the unilateral effects of concrete in tension and compression, respectively, which are usually set as and .
Based on the experimental data, Yu and Wu [27] assumed that the plastic strain of concrete under uniaxial loading satisfieswhere is the state variable with and referring to uniaxial tension and compression, respectively; and are calibrated parameters and are commonly set as and ; is the total strain and is the strain corresponding to the elasticity limit.
For uniaxial loading, , the effective stress of concrete can be calculated as
Yu and Wu [27] also recommended a piecewise function to describe the damage evolution:where is the initial damage; and are the damage and internal variable of concrete at the peak strain ; and are parameters controlling the softening curve of concrete; , , and can be computed as below according to the continuity requirements:where is the ratio of nominal stress to effective stress at ; is the derivative of effective stress with respect to strain at . The damage law is selected not only because of its applicability in recording both the tensile and the compressive degeneration mechanisms of concrete, but also because of its feasibility in adjusting the softening behaviors of concrete. The influences of and on the evolution discipline of are shown in Figure 5.
(a) Influence of
(b) Influence of
(c) Influence of
(d) Influence of
4. Solution of Local Problem
The solution of local problem, also known as stress updating, can be decomposed into three steps according to the concept of operator split [15, 16, 28], that is, the elastic predictor, the plastic corrector, and the damage corrector. The elastic predictor and the plastic corrector are executed to obtain the effective stresses, the internal variables, and the plastic multiplier that satisfy the updated yield surface, and the damage corrector is used to compute damage and nominal stresses.
4.1. Elastic Predictor and Plastic Corrector
Assuming that the state with variables is already known at time , then the effective stress and damage can be computed as
For time increment , we give . Let the strain increment of concrete within be ; the total strain and effective stress of concrete at time then satisfywhere is the trial stress at time . If , then concrete still remains in elasticity and the trial values are the final values of ; that is, , , and ; if , then plasticity develops and all trial values should be recalculated. Let all the unknowns to be solved at be ; we then havewhere and are the dimensions of the stress tensor and the number of the internal variables, respectively; is the identity matrix of order . A backward Euler scheme is used in the iterative solution of (26). Given a specific iteration step within , that is, , for the next iteration step , can be derived aswhere and are the residual and Jacobian matrices of (26) with the following expressions:where refers to the transposition.
For the first iteration step, the initial solution of all the unknowns of (26) is usually chosen from the elastic trial state; that is, . Convergence is satisfied when the Euclidean norm of the residuals is less than given tolerance; that is, .
4.2. Substepping Strategy
The convergence of (26) may not be easily made in some cases, especially when the strain increment within is large. Therefore, an adaptive substepping strategy [17, 29] is employed if the convergence cannot be made for given tolerance or a maximum number of iterations in the solution. The total strain increment from to can then be divided into parts according to (not necessary equal)where .
4.3. Damage Corrector
Once the solutions of elastic predictor and plastic corrector have been finished, the nominal stress can be updated as
4.4. Consistent Tangent Stiffness
The consistent tangent stiffness at time can be calculated by differentiating both sides of (11) with respect to the strain increment:
If substepping strategy is used, the computations of and should be derived step by step recursively. For substep , (26) can be written as
Let and then differentiating both sides of (32) with respect to we obtain where and .
For substep ,
Combining (34) and (35), we get the following equation according to the recursive relation: where uninterested submatrix in (36) is represented by bold dots . Obviously, and .
5. Solution of Global Problem
In the failure analysis of concrete material, convergence is a challenging problem in both stress update as mentioned above and capture of global responses. The traditional force control method finds it difficult to determine the exact load limit and is unable to trace the softening curve. The displacement control method, performing better than force control, has a strong dependency on selecting the proper loading freedoms. The classical arc-length method, which is able to trace “snap-through” and “snap-back” phenomena, is poor in localization problems of concrete. To help convergence, local arc-length method was proposed, but it requires the failure pattern of concrete to be a priori predictable. Gutiérrez [21], Verhoose et al. [22], and van der Meer et al. [23] developed a new arc-length method whose constraint path is constructed on the energy dissipation of a system. The method was proven to have better applicability than other methods in analyzing brittle/ductile failure. However, all the existing dissipation-based arc-length methods depend on simple constitutive models such as elastic damage model and elastic-plastic model. In this section, we have derived a constraint path, considering the combined contribution of plasticity and damage to the amount of released energy.
The basic function by use of the arc-length method in global solution of a domain iswhere is the equivalent internal force; is the displacement vector; is the external load factor; is a prescribed loading pattern; is the constraint path.
5.1. Energy Release Rate
The energy release rate of can be computed aswithwhere and are the exerted power and rate of elastic energy, respectively; is the vector of applied external forces; is the nodal velocity.
Differentiating (40) to time, for the item , we can obtain
The rate of the second item can be derived as follows:where ; which can be obtained from (31).
Equations (41) and (42) are combined to obtain
Substituting (39) and (43) into (38) yields
A forward Euler discretization is used to obtain the incremental path-following constraint as a function of the path parameter :where and are the displacement increment vector and load increment factor with regard to current load step, respectively; , , and are the converged displacement vector, converged load factor, and . The dissipation increment is illustrated as the shaded area in Figure 6.
5.2. Solution Scheme
The solution algorithm can be applied with the help of Sherman-Morrison formula, where the unknowns to be solved at iteration step arewith
According to (45), we have
In our study, the global solution scheme is carried out in a hybrid scheme, as shown as follows.
The flowchart of the global solution scheme is presented as follows:(1)Initialization: ; set , , , and .(2)Solve a load step from to (using force control):(a)Solve by using Newton-Raphson method and record the iteration numbers .(b)Check . If not, do and ; then go to the next load step ((a)); else, do the following.(c)Compute the released energy according to (44).(d)Compute , with .(e)Check . If not, do or ; then go to step ().(3)Switch to dissipation-based arc-length method:(a)Update , from local computations.(b)Update , , and .(c)Start iteration: Compute and according to (46):(d)Check convergence. If not, go back to (3(a)); else, record iteration numbers .(e)Compute , with .(f)Update ; then go to next load step (3(a)).In the above, and are thresholds that control the step length. The global solution is regarded as converged once the residual of the load satisfies
6. Model Validation
To validate the applicability and effectiveness of the model, a cubic plain concrete specimen with size of 100 × 100 × 100 (unit: mm) is simulated and computed under different loading conditions. The numerical results are compared with experiment data. Material parameters for the yield function and the plastic potential are listed in Table 2. The convergence tolerance values for the local and global solutions are set as %.
6.1. Uniaxial Loading
Two experiments [30, 31] are simulated to investigate the applicability of the model in describing the uniaxial tensile and compressive loading behaviors of concrete. Corresponding parameters are listed in Tables 3 and 4, respectively. As shown in Figure 7, it is seen that the computed stress-strain curve agrees well with experimental data in both ascending and descending branches, indicating that the proposed model is able to capture the uniaxial hardening/softening properties of concrete in both tension and compression.
(a)
(b)
6.2. Biaxial Loading
The same parameters in Tables 3 and 4 are used to observe the performance of the proposed model under biaxial loadings [32], with exception of . In Figure 8, the comparisons of numerical results obtained from different biaxial stress ratios with experimental data are shown. All the values of are regularized by , with taking the values of −29.5 MPa and −32.8 MPa for biaxial tensile and compressive loadings, respectively. Good agreements are seen in these two cases.
(a)
(b)
7. Conclusion
This paper introduces a numerical framework of implementing the spatial elastoplastic damage model of concrete via the recently developed IGA method. For the material model, a nonassociated flow rule is considered to describe the plasticity of concrete, where a three-parameter Barcelona yield surface and a modified Drucker-Prager plastic potential are used. The damage of concrete is considered as a combination of tensile and compressive damage variables, with the evolution of each driven by the internal variables and expressed by a piecewise function. The main procedure is divided into two kinds of solutions: local and global. For the local kind, the decoupling of plasticity and damage is achieved via the philosophy of operator split and both return-mapping algorithm and substepping strategy are used to help convergence. For the global kind, a dissipation-based arc-length method with constraint path is developed, considering contribution of plasticity and damage to the released amount of energy. Through comparisons between numerical simulations and experimental data, the approach thus employed has been proven to be reliable in reflecting the concrete properties and can therefore be further applied to the advanced analysis of concrete structures.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
The research for which this paper is prepared was supported by the Chinese Scholarship Council, whose precious help is gratefully acknowledged.