Numerical Implementation of Spatial Elastoplastic Damage Model of Concrete in the Framework of Isogeometric Analysis Approach

This paper is a study of the numerical implementation of the spatial elastoplastic damagemodel 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.Thedamage evolution of concrete driven by the internal variables is expressed by a piecewise function. In the study, the returnmapping 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.


Introduction
Armed with the progress of computer science and the indepth 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, computeraided 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][7][8], strong coupling [9,10], or effective stress approaches [3,[11][12][13][14][15][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 2 Mathematical Problems in Engineering 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.
For  = 0,  ,0 () is defined as For  ≥ 1,  , () can be calculated according to the Coxde 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   ( = 1, 2, . . ., ) 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  , (), where , , and  refer to the number of -spline basis functions of  , (),  , (), and  , () in the directions of , V, and . ,, is the control point of  ,, ,, and  ,, is the weight.

Numerical Formulation for NURBS Solid Model.
The numerical integral of a given function B (i.e., stiffness matrix K) 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 Ωele , the parametric space Ωele , and the physical space Ω ele , 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 are A 1,5
where  pt is the total number of control points for a single NURBS solid and ( pt ,  pt ,  pt ) is the Cartesian coordinates of control point pt.
The integral of function B on Ω can be derived as

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

Basic Functions. The basic functions for elastoplastic damage model are
where ,   , and   are the total, the elastic, and the plastic strains, respectively;  and  are the nominal and the effective stresses, respectively;  is the damage; D  is the initial elastic stiffness matrix; m is the flow vector that can be derived from the plastic potential Φ; h is the plastic modulus;  is the set of internal variables; λ is the plastic multiplier that satisfies the Kuhn-Tucker condition where (, ) is the yield function.Two internal variables of   and   are introduced to record the hardening/softening behaviors of concrete under with where  1 and  2 are the first stress invariant and the second deviatoric stress invariant computed by the effective stress tensor , respectively; σmax is algebraically the maximum principal effective stress with ⟨⟩ = (|| + )/2, the Macaulay bracket function;  =  0 / 0 with  0 and  0 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 0.5 <   < 1.0.

Plastic Potential.
A hyperbolic form of the Drucker-Prager function is used as the plastic potential of concrete: where  is the dilation angle;  0 is the uniaxial tensile yield stress;  is the eccentricity parameter.

Plastic Modulus.
The plastic modulus is written as with ( σ) 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]: with where   =   (  ) 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   = 0 and   = 1.
Based on the experimental data, Yu and Wu [27] assumed that the plastic strain of concrete under uniaxial loading satisfies where ℵ ∈ {, } is the state variable with ℵ =  and ℵ =  referring to uniaxial tension and compression, respectively;  and  are calibrated parameters and are commonly set as  = 1.79 and  = 1.18;  ℵ is the total strain and  0ℵ is the strain corresponding to the elasticity limit.

Mathematical Problems in Engineering
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  0ℵ = 0 is the initial damage;  ℵ and  ℵ are the damage and internal variable of concrete at the peak strain  ℵ ;  2ℵ and  3ℵ are parameters controlling the softening curve of concrete;  1ℵ ,  1ℵ , and  2ℵ 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  2ℵ and  3ℵ on the evolution discipline of  ℵ are shown in Figure 5.

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.

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 +1  =   + Δ.Let the strain increment of concrete within Δ be Δ; the total strain and effective stress of concrete at time +1  then satisfy where +1  trial is the trial stress at time +1 .If ( +1  trial ,  ) ≤ 0, then concrete still remains in elasticity and the trial values are the final values of +1 ; that is, +1  = +1  trial , +1  =  , and +1  = 0; if ( +1  trial ,  ) > 0, then plasticity develops and all trial values should be recalculated.Let all the unknowns to be solved at +1  be +1 x = [ +1  +1  +1  ] T ; we then have where (= 6) and (= 2) are the dimensions of the stress tensor and the number of the internal variables, respectively; I * 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, T , for the next iteration step  + 1, +1 +1 x can be derived as where E( +1  x) and +1  J Local are the residual and Jacobian matrices of (26) with the following expressions: where T 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; when the Euclidean norm of the residuals is less than given tolerance; that is, ‖E( +1  x)‖ ⩽ TOL Local .

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 +1  can then be divided into  parts according to (not necessary equal) where   = Δ  /Δ.

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 dissipationbased 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 Ω is where f int (u) is the equivalent internal force; u is the displacement vector;   is the external load factor; f is a prescribed loading pattern; (u,   ) is the constraint path.

Energy Release Rate.
The energy release rate of Ω can be computed as with where  and V are the exerted power and rate of elastic energy, respectively; f ext =   f is the vector of applied external forces; u is the nodal velocity.Differentiating (40) to time, for the item  1 , we can obtain The rate of the second item  2 can be derived as follows: where Λ=d  /d; C=d/d 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 Δu and Δ  are the displacement increment vector and load increment factor with regard to current load step, respectively; 0 u, 0   , and 0 f * are the converged displacement vector, converged load factor, and f * .The dissipation increment is illustrated as the shaded area in Figure 6.

Solution Scheme.
The solution algorithm can be applied with the help of Sherman-Morrison formula, where the unknowns to be solved at iteration step  are with According to (45), we have In our study, the global solution scheme is carried out in a hybrid scheme, as shown as follows.
(2) Solve a load step from 0   to   (using force control):

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 TOL Local = TOL Global = 0.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.

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

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

Figure 2 :Figure 3 :
Figure 2: Comparison on geometric modeling of a circle by FEM and IGA.

Figure 4 :
Figure 4: Evolution curves of NIGA with different basis orders.

Figure 6 :
Figure 6: Energy dissipation increment for elastoplastic damage model under uniaxial loading.

Figure 7 :
Figure 7: Numerical results compared with experimental data: (a) uniaxial tension and (b) uniaxial compression.

Figure 8 :
Figure 8: Numerical results compared with experimental data: (a) biaxial tension and (b) biaxial compression.

Table 2 :
Material parameters for the yield function and plastic potential.

Table 3 :
Material parameters for uniaxial tension.

Table 4 :
Material parameters for uniaxial compression.