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

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 post-crack 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, Perez-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. Gutierrez [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. B-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 B-Spline), NURBS surface, and NURBS solid, respectively. A pth-order B-spline [C.sub.B]([xi]) can be constructed by n B-spline basis functions [N.sub.i,p] defined on the knot sequence [THETA] = {[[xi].sub.1], [[xi].sub.2], ..., [[xi].sub.n+p+1]} and their corresponding control points [A.sub.i]:

[C.sub.B] ([xi]) = [n.summation over (i=1)] [N.sub.i,p] ([xi]) [A.sub.i]

The coordinates of knots should be monotonically non-decreasing [mathematical expression not reproducible], with two nonoverlapped neighboring knots constructing an element interval.

For p = 0, [N.sub.i,0] ([xi]) is defined as

[mathematical expression not reproducible]. (2)

For p [greater than or equal to] 1, [N.sub.i,p]([xi]) can be calculated according to the Coxde Boor recursion formula:

[mathematical expression not reproducible]. (3)

A pth-order NURBS basis function of is [R.sub.i,p] ([xi]) defined by giving each [N.sub.i,p] ([xi]) a proper weight, and a pth-order NURBS is constructed by n pth-order NURBS basis functions and their corresponding control points [A.sub.i]:

[mathematical expression not reproducible], (4)

where [[omega].sub.i] (i = 1,2, ..., n) is the weight for [N.sub.i,p] ([xi]). NURBS basis function degenerates to B-spline basis function when all the weights take the same value.

The NURBS solid is a tensor product of [N.sub.i,p] ([xi]), [M.sub.j,q]([eta]), and [P.sub.k,r] ([zeta]):

[mathematical expression not reproducible], (5)

where n, m, and I refer to the number of B-spline basis functions of [N.sub.i,p] ([xi]), [M.sub.j,q] ([eta]), and [P.sub.k,r]([zeta]) in the directions of u, v, and w. [A.sub.i,j,k] is the control point of [R.sup.p,q,r.sub.i,j,k] and [[omega].sub.i,j,k] is the weight.

2.2. Numerical Formulation for NURBS Solid Model. The numerical integral of a given function B (i.e., stiffness matrix K) on a domain [OMEGA] consists of ne 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 [[??].sup.ele], the parametric space [[??].sup.ele], and the physical space [[OMEGA].sup.ele], are defined by the local, knot, and Cartesian coordinates, respectively. Hence,

[mathematical expression not reproducible], (6)

[mathematical expression not reproducible], (7)

[mathematical expression not reproducible]. (8)

As shown in Figure 1, we define the mappings of parametric space [right arrow] parent space and physical space [right arrow] 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 [mathematical expression not reproducible], and (%, y, z), respectively, then the Jacobian matrices for mapping 1 and mapping 2 are

[mathematical expression not reproducible]. (9)

where [n.sub.pt] is the total number of control points for a single NURBS solid and ([x.sub.ipt], [y.sub.ipt], [z.sub.ipt]) is the Cartesian coordinates of control point pt.

The integral of function B on [OMEGA] can be derived as

[mathematical expression not reproducible] (10)

with determinant [absolute value of (J)] = [absolute value of ([J.sub.2])] x [absolute value of ([J.sub.1])]

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.

Considering a cubic as shown in Table 1 and discretizing the cubic to n x n x n pth-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 n and p. 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 n, 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.

3. Constitutive Model of Concrete

3.1. Basic Functions. The basic functions for elastoplastic damage model are

[mathematical expression not reproducible], (11)

where [epsilon], [[epsilon].sup.[epsilon]], and [[epsilon].sup.p] are the total, the elastic, and the plastic strains, respectively; [sigma] and [bar.[sigma]] are the nominal and the effective stresses, respectively; d is the damage; [D.sup.e] is the initial elastic stiffness matrix; m is the flow vector that can be derived from the plastic potential [PHI]; h is the plastic modulus; k is the set of internal variables; [??] is the plastic multiplier that satisfies the Kuhn-Tucker condition

[??] [greater than or equal to] 0, F([bar.[sigma]], k) [less than or equal to] 0 (12) with [??]F = 0,

where F([bar.[sigma]], k) is the yield function.

Two internal variables of [k.sub.t] and [k.sub.c] are introduced to record the hardening/softening behaviors of concrete under tension and compression, respectively; that is, [mathematical expression not reproducible] are the tensile and compressive equivalent plastic strains, respectively.

3.2. Yield Function. The yield surface of concrete is described by the modified Barcelona function

[mathematical expression not reproducible] (13)

with

[mathematical expression not reproducible], (14)

where [[bar.I].sub.1] and [[bar.J].sub.2] are the first stress invariant and the second deviatoric stress invariant computed by the effective stress tensor [bar.[sigma]], respectively; [[??].sub.max] is algebraically the maximum principal effective stress with <x> = ([absolute value of (x)] + x)/2, the Macaulay bracket function; [theta] = [f.sub.b0]/[f.sub.c0] with [f.sub.b0] and [f.sub.c0] being the uniaxial and biaxial yield strengths of concrete in compression, respectively; [[bar.c].sub.t] (Kt) and [[bar.c].sub.c]([k.sub.c]) are the effective tensile and compressive cohesion stresses of concrete, respectively; [[gamma].sub.F] is introduced only when concrete is subjected to triaxial compression with 0.5 < [K.sub.c] < 1.0.

3.3. Plastic Potential. A hyperbolic form of the Drucker-Prager function is used as the plastic potential of concrete:

[mathematical expression not reproducible], (15)

where [phi] is the dilation angle; [f.sub.t0] is the uniaxial tensile yield stress; [rho] is the eccentricity parameter.

3.4. Plastic Modulus. The plastic modulus is written as

[mathematical expression not reproducible] (16)

with r([??]) being the weight function which can be computed as

[mathematical expression not reproducible]. (17)

3.5. Damage Law. The total damage d 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]:

d = 1 - (1 - [s.sub.t][d.sub.c]([k.sub.c]))(1 - [s.sub.c][d.sub.t]([k.sub.t])) (18)

with

[mathematical expression not reproducible] (19),

where [d.sub.t] = [d.sub.t]([k.sub.t]) and [d.sub.c] = [d.sub.c]([k.sub.c]) are the tensile and compressive damage variables, respectively; [w.sub.t] and [w.sub.c] are the stiffness recovery factors reflecting the unilateral effects of concrete in tension and compression, respectively, which are usually set as [w.sub.t] = 0 and [w.sub.c] = 1.

Based on the experimental data, Yu and Wu [27] assumed that the plastic strain of concrete under uniaxial loading satisfies

[[epsilon].sup.p.sub.N] = [alpha][([[epsilon].sub.N] - [[epsilon].sub.0N]).sup.[beta]], (20)

where N [member of] {t, c} is the state variable with R = t and R = c referring to uniaxial tension and compression, respectively; [alpha] and [beta] are calibrated parameters and are commonly set as [alpha] = 1.79 and [beta] = 1.18; [[epsilon].sub.N] is the total strain and [[epsilon].sub.0R] is the strain corresponding to the elasticity limit.

For uniaxial loading, [k.sub.N] = [[epsilon].sup.p.sub.N], the effective stress of concrete can be calculated as

[mathematical expression not reproducible]. (21)

Yu and Wu [27] also recommended a piecewise function to describe the damage evolution:

[mathematical expression not reproducible], (22)

where [d.sub.0N] = 0 is the initial damage; d/N and k/n are the damage and internal variable of concrete at the peak strain [[epsilon].sub.fN]; [A.sub.2N] and [B.sub.3N] are parameters controlling the softening curve of concrete; [A.sub.1N], [B.sub.1N], and [B.sub.2N] can be computed as below according to the continuity requirements:

[mathematical expression not reproducible], (23)

where [[lambda].sub.N] is the ratio of nominal stress to effective stress at [[epsilon].sub.FN]; [[bar.[sigma]]'.sub.N]([[epsilon].sub.fN] is the derivative of effective stress with respect to strain at e/N. 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 [A.sub.2N] and [S.sub.3N] on the evolution discipline of [d.sub.N] are shown in Figure 5.

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 {[sup.n][epsilon], [sup.n][[epsilon].sup.p], [sup.n]k}| is already known at time [sup.n]f, then the effective stress [sup.n][bar.[sigma]] and damage [sup.n]d can be computed as

[mathematical expression not reproducible]. (24)

For time increment [DELTA]f, we give [sup.n+1]f = [sup.n]f + [DELTA]f. Let the strain increment of concrete within [DELTA]t be [DELTA][epsilon]; the total strain and effective stress of concrete at time [sup.n+1]f then satisfy

[mathematical expression not reproducible], (25)

where [sup.n+1][[bar.[sigma]].sup.trial] is the trial stress at time [supn+1]t. If F([sup.n+1][[bar.[sigma]].sup.trial], [sup.n]k) [less than or equal to] 0, then concrete still remains in elasticity and the trial values are the final values of [sup.n+1]t; that is, [mathematical expression not reproducible], then plasticity develops and all trial values should be recalculated. Let all the unknowns to be solved at [mathematical expression not reproducible]; we then have

[mathematical expression not reproducible], (26)

where n[bar.[sigma]](= 6) and nk(= 2) are the dimensions of the stress tensor and the number of the internal variables, respectively; [I.sub.*] is the identity matrix of order *. A backward Euler scheme is used in the iterative solution of (26). Given a specific iteration step k within [DELTA]t, that is, [mathematical expression not reproducible], for the next iteration step k + 1, [sup.n+1.sub.k+1]x can be derived as

[mathematical expression not reproducible], (27)

where E([sup.n+1.sub.k]x) and [sup.n+1.sub.k][J.sub.Local] are the residual and Jacobian matrices of (26) with the following expressions:

[mathematical expression not reproducible], (28)

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; that is, [sup.n+1.sub.0]x = [[[sup.n+1][[bar.[sigma]].sub.trial] [sup.k] 0].sup.T] Convergence is satisfied when the Euclidean norm of the residuals is less than given tolerance; that is, [parallel]E([sup.n+1.sub.k]x)[parallel] [less than or equal to] [TOL.sub.Local].

4.2. Substepping Strategy. The convergence of (26) may not be easily made in some cases, especially when the strain increment [DELTA][epsilon] within [DELTA]t 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 [DELTA]e from [sup.n]t to [sup.n+1]t can then be divided into m parts according to (not necessary equal)

[mathematical expression not reproducible], (29)

where [[alpha].sub.k] = [DELTA][[epsilon].sub.k]/[[DELTA].sub.[epsilon]].

4.3. Damage Corrector. Once the solutions of elastic predictor and plastic corrector have been finished, the nominal stress can be updated as

[mathematical expression not reproducible]. (30)

4.4. Consistent Tangent Stiffness. The consistent tangent stiffness at time [sup.n+1]t can be calculated by differentiating both sides of (11) with respect to the strain increment:

[mathematical expression not reproducible]. (31)

If substepping strategy is used, the computations of [mathematical expression not reproducible] should be derived step by step recursively. For substep k (k [not equal to] 1), (26) can be written as

[mathematical expression not reproducible]. (32)

Let

[mathematical expression not reproducible]; (33)

and then differentiating both sides of (32) with respect to [DELTA][epsilon] we obtain

[mathematical expression not reproducible], (34)

where [mathematical expression not reproducible].

For substep k = 1,

[mathematical expression not reproducible]. (35)

Combining (34) and (35), we get the following equation according to the recursive relation:

[mathematical expression not reproducible], (36)

where uninterested submatrix in (36) is represented by bold dots *. Obviously, [mathematical expression not reproducible].

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. Gutierrez [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 [OMEGA] is

[mathematical expression not reproducible]. (37)

where [f.sub.int](u) is the equivalent internal force; u is the displacement vector; [[lambda].sub.f] is the external load factor; [??] is a prescribed loading pattern; g(u, [[lambda].sub.F]) is the constraint path.

5.1. Energy Release Rate. The energy release rate of Q can be computed as

G = P - [??] (38)

with

[mathematical expression not reproducible] (39)

[mathematical expression not reproducible], (40)

where P and [??] are the exerted power and rate of elastic energy, respectively; [f.sub.ext] = [[lambda].sub.F][??] is the vector of applied external forces; [??] is the nodal velocity.

Differentiating (40) to time, for the item [V.sub.1], we can obtain

[mathematical expression not reproducible]. (41)

The rate of the second item [V.sub.2] can be derived as follows:

[mathematical expression not reproducible], (42)

where [LAMBDA] = d[[epsilon].sup.p]/d[epsilon]; C=d[sigma]/d[epsilon] which can be obtained from (31).

Equations (41) and (42) are combined to obtain

[mathematical expression not reproducible]. (43)

Substituting (39) and (43) into (38) yields

[mathematical expression not reproducible]. (44)

A forward Euler discretization is used to obtain the incremental path-following constraint as a function of the path parameter [DELTA][tau]:

[mathematical expression not reproducible], (45)

where [DELTA]u and [DELTA][[lambda].sub.F] are the displacement increment vector and load increment factor with regard to current load step, respectively; [sup.0]u, [sup.0][[lambda].sub.F], and [sup.0][f.sup.*] are the converged displacement vector, converged load factor, and [f.sup.*]. 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 k are

[mathematical expression not reproducible] (46)

with

[mathematical expression not reproducible]. (47)

According to (45), we have

[mathematical expression not reproducible]. (48)

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: [[lambda].sub.F] = 0; set [n.sub.d], [DELTA][[lambda].sub.F], [DELTA][[tau].sub.min], and [DELTA][[tau].sub.max].

(2) Solve a load step from [sup.0][[lambda].sub.F] to [[lambda].sub.F] (using force control):

(a) Solve [f.sub.int] = [[lambda].sub.F][??] by using Newton-Raphson method and record the iteration numbers n.

(b) Check n > [n.sub.d]. If not, do [sup.0][[lambda].sub.F] = [[lambda].sub.F] and [[lambda].sub.F] = [[lambda].sub.F] + [DELTA][[lambda].sub.F]; then go to the next load step (2(a)); else, do the following.

(c) Compute the released energy G according to (44).

(d) Compute [DELTA][tau] = [0.5.sup.r] x G, with r = 0.25 x (n - [n.sub.d]).

(e) Check [DELTA][[tau].sub.min] [less than or equal to] [DELTA][tau] [less than or equal to] [DELTA][[tau].sub.max]. If not, do [DELTA][tau] = [DELTA][[tau].sub.min] or [DELTA][tau] = [DELTA]-[[tau].sub.max]; then go to step (3).

(3) Switch to dissipation-based arc-length method:

(a) Update K, [f.sub.int] from local computations.

(b) Update [f.sup.*], v, and w.

(c) Start iteration:

r = [[lambda].sub.F][??] - [f.sub.int] [d.sup.I] = [K.sup.-1] x r, [d.sub.II] = -[K.sup.-1] x [??]. (49)

Compute du and d[[lambda].sub.F] according to (46):

u = u + du, [[lambda].sub.F] = [[lambda].sub.F] + d[[lambda].sub.F].

(d) Check convergence. If not, go back to (3(a)); else, record iteration numbers n.

(e) Compute [DELTA][tau] = [0.5.sup.r] x G, with r = 0.25 x (n - [n.sub.d]).

(f) Update [sup.0][[lambda].sub.F] = [[lambda].sub.F]; then go to next load step (3(a)).

In the above [DELTA][[tau].sub.min] and [DELTA][[tau].sub.max] are thresholds that control the step length. The global solution is regarded as converged once the residual of the load satisfies

[mathematical expression not reproducible]. (51)

6. Model Validation

To validate the applicability and effectiveness of the model, a cubic plain concrete specimen with size of 100 x 100 x 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.sub.LoCal] = [TOL.sub.Global] = 0.1%.

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.

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 [f.sub.fc]. In Figure 8, the comparisons of numerical results obtained from different biaxial stress ratios with experimental data are shown. All the values of [[sigma].sub.1] are regularized by [[beta].sub.p] = [f.sub.fc], with [f.sub.fc] 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.

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.

http://dx.doi.org/10.1155/2016/4273024

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.

References

[1] J. A. Cottrell, T. J. R. Hughes, and Y. Bazilevs, Isogeometric Analysis: Toward Integration of CAD and FEA, John Wiley & Sons, 2009.

[2] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs, "Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement," Computer Methods in Applied Mechanics and Engineering, vol. 194, no. 39-41, pp. 4135-4195, 2005.

[3] J. Mazars and G. Pyaudier-Cabot, "Continuum damage theory--application to concrete," Journal of Engineering Mechanics, vol. 115, no. 2, pp. 345-365,1989.

[4] T. Rabczuk, J. Akkermann, and J. Eibl, "A numerical model for reinforced concrete structures," International Journal of Solids and Structures, vol. 42, no. 5-6, pp. 1327-1354, 2005.

[5] P. Grassl, K. Lundgren, and K. Gylltoft, "Concrete in compression: a plasticity theory with a novel hardening law," International Journal of Solids and Structures, vol. 39, no. 20, pp. 5205-5223, 2002.

[6] W. B. Kratzig and R. Polling, "An elasto-plastic damage model for reinforced concrete with minimum number of material parameters," Computers & Structures, vol. 82, no. 15-16, pp. 1201-1215, 2004.

[7] G. Meschke and R. Lackner, "Anisotropic modelling of cracked concrete based on plasticity-damage theory Computational plasticity: fundamentals and applications," in Proceedings of the International Center for Numerical Methods in Engineering (CIMNE '97), pp. 1543-1550, Barcelona, Spain, 1997.

[8] S. Oller, E. Onate, J. Oliver, and J. Lubliner, "Finite element nonlinear analysis of concrete structures using a 'plastic-damage model"' Engineering Fracture Mechanics, vol. 35, no. 1-3, pp. 219-231,1990.

[9] F. Armero and S. Oller, "A general framework for continuum damage models. I. INFinitesimal plastic damage models in stress space," International Journal of Solids and Structures, vol. 37, no. 48-50, pp. 7409-7436, 2000.

[10] B. Luccioni, S. Oller, and R. Danesi, "Coupled plastic-damaged model," Computer Methods in Applied Mechanics and Engineering, vol. 129, no. 1-2, pp. 81-89,1996.

[11] R. Faria, J. Oliver, and M. Cervera, "A strain-based plastic viscous-damage model for massive concrete structures," International Journal of Solids and Structures, vol. 35, no. 14, pp. 1533-1558, 1998.

[12] J. W. Ju, "On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects," International Journal of Solids and Structures, vol. 25, no. 7, pp. 803-833, 1989.

[13] J. LeeandG. L. Fenves, "Plastic-damage model for cyclic loading of concrete structures," Journal of Engineering Mechanics, vol. 124, no. 8, pp. 892-900, 1998.

[14] J. Lee and G. L. Fenves, "A return-mapping algorithm for plastic-damage models: 3-D and plane stress formulation,"

International Journal for Numerical Methods in Engineering, vol. 50, no. 2, pp. 487-506, 2001.

[15] J. C. Simo and J. W. Ju, "Strain- and stress-based continuum damage models-I. Formulation," International Journal of Solids and Structures, vol. 23, no. 7, pp. 821-840,1987.

[16] J. C. Simo and J. W. Ju, "Strain- and stress-based continuum damage models-II. Computational aspects," International Journal of Solids and Structures, vol. 23, no. 7, pp. 841-869,1987.

[17] A. Perez-Foguet, A. Rodriguez-Ferran, and A. Huerta, "Consistent tangent matrices for substepping schemes," Computer Methods in Applied Mechanics and Engineering, vol. 190, no. 3536, pp. 4627-4647, 2001.

[18] E. Riks, "An incremental approach to the solution of snapping and buckling problems," International Journal of Solids and Structures, vol. 15, no. 7, pp. 529-551,1979.

[19] E. Ramm, "Strategies for tracing the nonlinear response near limit points," in Nonlinear Finite Element Analysis in Structural Mechanics, W. Wunderlich, E. Stein, and K. J. Bathe, Eds., pp. 63-89, Springer, Berlin, Germany, 1981.

[20] M. A. Crisfield, "A fast incremental/iterative solution procedure that handles 'snap-through'," Computers and Structures, vol. 13, no. 1-3, pp. 55-62, 1981.

[21] M. A. Gutierrez, "Energy release control for numerical simulations of failure in quasi-brittle solids," Communications in Numerical Methods in Engineering, vol. 20, no. 1, pp. 19-29, 2004.

[22] C. V. Verhoose, J. J. C. Remmers, and M. A. Gutierrez, "A dissipation-based arc-length method for robust simulation of brittle and ductile failure," International Journal for Numerical Methods in Engineering, vol. 77, no. 9, pp. 1290-1321, 2009.

[23] F. P. van der Meer, C. Oliver, and L. J. Sluys, "Computational analysis of progressive failure in a notched laminate including shear nonlinearity and fiber failure," Composites Science and Technology, vol. 70, no. 4, pp. 692-700, 2010.

[24] V. P. Nguyen, C. Anitescu, S. P. A. Bordas, and T. Rabczuk, "Isogeometric analysis: an overview and computer implementation aspects," Mathematics and Computers in Simulation, vol. 117, pp. 89-116, 2015.

[25] J. Lubliner, J. Oliver, S. Oller, and E. Onate, "A plastic-damage model for concrete," International Journal of Solids and Structures, vol. 25, no. 3, pp. 299-326,1989.

[26] H. Hibbitt, B. Karlsson, and P. Sorensen, Abaqus v6. 12 Documentation-ABAQUS Analysis User's Manual, Dassault Systemes Simulia, Providence, RI, USA, 2012.

[27] H.-X. Yu and J.-H. Wu, "An elastoplastic damage constitutive model for concrete based on undamaged state," Engineering Mechanics, vol. 26, no. 10, pp. 79-86, 2009 (Chinese).

[28] J. C. Simo and T. J. R. Hughes, Computational Inelasticity, vol. 7 of Interdisciplinary Applied Mathematics, Springer, New York, NY, USA, 1998.

[29] A. Saritas and F. C. Filippou, "Numerical integration of a class of 3d plastic-damage concrete models and condensation of 3d stress-strain relations for use in beam finite elements," Engineering Structures, vol. 31, no. 10, pp. 2327-2336, 2009.

[30] V. S. Gopalaratnam and P. S. Surendra, "Softening response of plain concrete in direct tension," ACI Journal Proceedings, vol. 82, pp. 310-323, 1985.

[31] I. D. Karsan and J. O. Jirsa, "Behavior of concrete under compressive loadings," Journal of the Structural Division, vol. 95, pp. 2543-2564,1969.

[32] H. Kupfer, H. K. Hilsdorf, and H. Rusch, "Behavior of concrete under biaxial stresses," ACI Journal Proceedings, vol. 66, no. 8, pp. 656-666,1969.

Cheng Ma, Wei-zhen Chen, and Jian-yuan Sun

Department of Bridge Engineering, Tongji University, Shanghai 200092, China

Correspondence should be addressed to Jian-yuan Sun; sunjy@tongji.edu.cn

Received 16 January 2016; Revised 3 May 2016; Accepted 25 May 2016

Academic Editor: Pierfrancesco Cacciola

Caption: FIGURE 1: Schematic and mappings of [[??].sup.ele], [[??].sup.ele], and [[OMEGA].sup.ele] in two-dimensional space.

Caption: FIGURE 2: Comparison on geometric modeling of a circle by FEM and IGA.

Caption: FIGURE 3: Evolution curves of NFEM and NIGA with respect to n.

Caption: FIGURE 4: Evolution curves of NIGA with different basis orders.

Caption: FIGURE 5: Influences of [A.sub.2N] and [B.sub.3N] on the evolution discipline of [d.sub.N].

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

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

Caption: FIGURE 8: Numerical results compared with experimental data: (a) biaxial tension and (b) biaxial compression.
```TABLE 1: Comparison between NFEM and NIGA.

Analyzing       Order of             NFEM                   NIGA
model            basis
function

[formula         p = 1          [(n + 1).sup.3]        [(n + 1).sup.3]
expression       p = 2     [(n + 1).sup.2] (4n + 1)    [(n + 2).sup.3]
not              p = 3     [(n + 1).sup.2] (7n + 1)    [(n + 3).sup.3]
reproducible]     ...                 ...                    ...
p        [(n + 1).sup.2] [(3p -     [(n + p).sup.3]
2)n + 1]

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

[theta]   [gamma]   [rho]   [phi]

1.15      3.0       0.1     n!6

Table 3: Material parameters for uniaxial tension.

[E.sub.0t] (MPa)    v     [[epsilon].sub.0t]   [[epsilon].sub.ft]
[DELTA][[tau].sub.min] (N x m)   [DELTA][[tau].sub.max] (N x m)

3.1 x [10.sup.4]   0.18        0.000105             0.00012
0.0001                            1.0

[E.sub.0t] (MPa)   [f.sub.ft] (MPa)   [d.sub.0t]

3.1 x [10.sup.4]         3.48            0.0

[E.sub.0t] (MPa)        [A.sub.2t]        [B.sub.3t]

3.1 x [10.sup.4]   0.7 [f.sub.ft.sup.2]      1.5

[E.sub.0t] (MPa)   [DELTA][[lambda].sub.f]   [n.sub.d]

3.1 x [10.sup.4]             1.0                 3

Table 4: Material parameters for uniaxial compression.

[E.sub.0c] (MPa)    v     [[epsilon].sub.0c]   [[epsilon].sub.fc]

3.1 x [10.sup.4]   0.18        0.00067               0.002

[E.sub.0c] (MPa)   [f.sub.fc] (MPa)   [d.sub.0c]   [A.sub.2c]

3.1 x [10.sup.4]        -27.6            0.0          0.75

[E.sub.0c] (MPa)   [B.sub.3c]   [DELTA][[lambda].sub.f]   [n.sub.d]

3.1 x [10.sup.4]      2.0                 1.0                 3

[E.sub.0c] (MPa)   [DELTA][[tau].sub.min] (N x m)

3.1 x [10.sup.4]               0.0001

[E.sub.0c] (MPa)   [DELTA][[tau].sub.max] (N x m)

3.1 x [10.sup.4]                5.0
```
COPYRIGHT 2016 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.