# Tension-Compression Damage Model with Consistent Crack Bandwidths for Concrete Materials

Authors: Wei He; Ying Xu; Yu Cheng; Peng-Fei Jia; Ting-Ting Fu
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

This paper proposes a tension-compression damage model for concrete materials, formulated within the framework of thermodynamics of irreversible processes. The aim of this work is to solve the following problems: the premature divergence of numerical solutions under general loading conditions due to the conflict of tensile and compressive damage bounding surfaces, which is a result of the application of the spectral decomposition method to distinguish tension and compression, and the unsatisfactory reproduction of distinct tension-compression behaviors of concrete by strain-driven damage models. The former is solved by the sign of the volumetric deformation, while the latter is solved via two separated dissipation mechanisms. Moreover, of specific interest is an improved solution to the problem of mesh-size dependency using consistent crack bandwidths, which takes into account situations with irregular meshes and arbitrary crack directions in the context of the crack band approach. The performance of the model is validated by the well-documented experimental data. The simplicity and the explicit integration of the constitutive equations render the model well suitable for large-scale computations.


---

## Body

## 1. Introduction

Concrete is still one of the most widely used construction materials in civil engineering applications [1]. Therefore, developing an efficient constitutive model for concrete is of great importance since it contributes to engineering designs, failure analyses, and life-cycle assessments by reducing costs and time [2]. Microscopically, the mechanical behavior of concrete materials is mainly governed by the nucleation and growth of microcracks and their coalescence into macrocracks [3], which may result in the failure of a structure. Physically, it is an irreversible thermodynamic process accompanying energy dissipation and changes of the microstructure [4]. Such a progressive degradation of material properties can be phenomenologically described by models [5–17] based on continuum damage mechanics (CDM), which is generally accepted by the scientific and engineering communities [18].Despite the achievements of CDM over a wide range of materials (e.g., [19–22]), the rational and accurate CDM modeling of concrete materials still represents a challenging task. Three modeling challenges have been identified in finite element (FE) analyses. The first challenge is the criterion distinguishing tension and compression, which is used for activating/inactivating damage evolutions in tension and compression. The existing models usually utilize the spectral decomposition of a stress or strain tensor into tensile (positive) and compressive (negative) eigenvalues (e.g., [23, 24]), in order to describe different behaviors of concrete under general loading conditions. However, in the use of this method, the following problem would arise: divergence of numerical solutions will occur under biaxial or reverse loading conditions due to the conflict of tensile and compressive damage bounding surfaces induced by Poisson’s ratio effect [25]. A better solution to this problem is the introduction of a weighted damage parameter (e.g., [5, 25]). Moreover, the sign of the volumetric strain ([26]) is used to distinguish what is tension and what is compression. However, it should be noted that the aforementioned method is not applicable to incompressible materials which have a value of Poisson’s ratio around 0.5 (e.g., fresh concrete and concrete under confinement) since the volumetric strain is null. Therefore, the present model is limited to the compressible materials characterized by Poisson’s ratio lower than 0.5.The second challenge is the modeling of tension-compression behaviors of concrete materials which exhibit large differences in their peak strengths and postpeak stages, as well as the unilateral effect characterized by a total or partial recovery of the degraded stiffness during the tension-compression load reversal. The classical strain-driven damage models incorporating one damage criterion alone may be incapable of predicting such distinct behaviors. Among them (e.g., [27–32]), one case can be mentioned as an example: the predicted strength of concrete under biaxial compression is seriously lower than that under uniaxial compression. To solve the problem, this paper proposes two separate damage criteria describing the state of damage under tension and compression, respectively, whereby the distinct behaviors can be reproduced in a physically realistic fashion. Due to the scope of this paper, the unilateral effect will not be discussed.The third challenge is to avoid the lack of mesh objectivity of numerical results when the modeling of strain-softening is taken into account. At the stage of strain-softening, the concrete does not lose its entire load-carrying capacity but exhibits a progressively decreasing residual strength. Therefore, capturing the strain-softening is important because it is responsible for predicting the behavior of overloaded structures. However, damage models based on conventional continuum theories usually produce mesh-dependent results in FE analyses [6]. The mesh dependency can be further classified into mesh-size [33] and mesh-bias [34] dependencies. The mesh-size dependency, which originates from the loss of ellipticity of the governing partial differential equations for static problems [35], which means that both the computed energy dissipation due to strain-softening and the corresponding force-displacement response are governed by the fineness of the FE mesh [36]. It also reflects a physically unrealistic fact that the infinitesimal element size leads to a vanishing energy dissipation during the structural failure [37]. The mesh-bias dependency, which is a result of the inadequacy of the spatial discretization procedure [38, 39], means that the numerical result depends on the orientation of the FE mesh in the sense that the computed direction of the strain localization band tends to propagate along continuous mesh lines [40], which may give rise to spurious crack trajectories. Although the problem of mesh-bias dependency is important for damage modeling of concrete materials, this study focuses only on the solution to the mesh-size dependency.To solve the problem of mesh-size dependency in the framework of continuum mechanics, various types of regularization methods have been developed that can be distinguished by the crack band approach (or fracture energy regularization) [33, 41, 42], nonlocal models including the integral-type [6, 43, 44] and the gradient-type [13, 45], micropolar method [46], and viscous or rate-dependent [11] methods. On one hand, the nonlocal method seems to be more recent (e.g., [28, 30, 47]). However, the nonlocal models have their own difficulties such as (i) incorrect crack initiation, ahead of the crack tip [48], and (ii) incorrect shielding effect with nonzero nonlocal interactions across a crack surface [49]. Furthermore, nonlocal methods have no potential for large-scale computations of concrete structures since they require small enough elements, e.g., their sizes are at least one-third of the characteristic length of the material [50], to subdivide strain localization bands of the structure, which are usually a priori unknown. On the other hand, the crack band approach is widely used in most of the commercial FE packages and is still an important analytical tool for practical applications (e.g., [51, 52]). One goal of the present work is to propose a constitutive model which needs to ensure a high algorithmic efficiency and an easy implementation; therefore, a fracture energy regularization based on the crack band approach [33] is adopted.One key parameter in the crack band approach is the crack bandwidth, through which the fracture energy (regarded to be a material property) is related to the volumetric (or specific) fracture energy, and it therefore ensures an objectivity of the energy dissipation regardless of the mesh refinement. The existing estimations of the crack bandwidth can be divided into three categories: (i) methods based on the square root of the element area (for two-dimensional elements) or the cubic root of the element volume (for three-dimensional elements) (e.g., [53–55]); (ii) the notable methods proposed by Oliver [56]; and (iii) projection methods (e.g., [57–60]). The methods in the first category, which estimate the crack bandwidth as a constant during the numerical analysis process, are less applicable to situations with irregular meshes or arbitrary crack directions [56]. The methods in the last two categories, in which crack bandwidths are varied with respect to crack directions and spatial positions, solve the aforementioned problems. However, they may encounter difficulties such as excessively large estimation [59] or misprediction of the crack bandwidth. This paper therefore develops an improved estimation of the crack bandwidth to circumvent the above difficulties. Moreover, two different crack bandwidths (one for the tensile strain-softening and the other for the compressive strain-softening) are taken into account in the present study, because of which concrete exhibits quite different localization behaviors in tension and compression [61]. It should also be noted that the crack band approach can only be regarded as a partial regularization method for mesh-size dependency in the sense that the objectivity of the energy dissipation is ensured [62], and thus the global structural response can be captured correctly [56], but the mathematical description of the rate boundary value problem (BVP) still remains ill-posed [63].The main purpose of this study is not to describe all the mechanical behaviors of concrete but rather to give methodological efforts to solve the aforementioned difficulties. Therefore, the plastic dissipation is not considered and the isotropic damage is assumed for simplification. Accordingly, the objectives of this paper are fourfold: (i) to define the criterion distinguishing tension and compression by the sign of the volumetric deformation; (ii) to propose a model to describe distinct tension-compression behaviors of concrete materials under a framework similar to that of Cervera et al. [11] and Comi and Perego [26]; (iii) to develop an improved estimation of the crack bandwidth using the projection method for solving the problem of mesh-size dependency; and (iv) to validate the proposed model by the selected benchmark tests.


## 2. Tension-Compression Damage Model with Consistent Crack Bandwidths

### 2.1. Helmholtz Free Energy Potential

The thermodynamics of irreversible processes provides a general framework for formulating constitutive equations. In other words, the state of the material can be described by the thermodynamic potential. Consequently, the stresses and other associated variables can be derived from this potential [64]. Because of the scope of this paper, the present model is developed based on the following assumptions: (i) damage evolutions are assumed to be isotropic and to occur in small strains and (ii) the rate-dependent behavior and the thermal effect are not considered. To take into account distinct tension-compression behaviors of concrete, the postulated thermodynamic potential in terms of Helmholtz free energy ψ is decomposed into volumetric and deviatoric parts ψv and ψd in a similar way as in the study of Cervera et al. [11] and Comi and Perego [26] as follows:(1)ρψε,dt,dc=ρψvεv,dt,dc+ρψdεd,dt,dc=ρψv+εv+,dt,dc+ρψv−εv−,dt,dc+ρψdεd,dt,dc,(2)ρψv+εv+,dt,dc=321−dt1−dcK0εv+:εv+,ρψv−εv−,dt,dc=321−dcK0εv−:εv−,(3)ρψdεd,dt,dc=1−Htrεdt1−H−trεdcG0εd:εd,where ρ is the mass density of the material, ε is the second-order strain tensor, dt and dc are the damage internal variables that represent the level of material degradation in tension and compression, K0 and G0 are the initial bulk and shear moduli of the undamaged material, respectively, H• is the Heaviside function defined by Hx=1ifx≥0andHx=0ifx<0, and tr• is the trace of a tensor. The tensile and compressive damage internal variables are monotonously increasing scalars, both of which grow from zero (undamaged material with the intact stiffness) to one (completely damaged material). εv and εd are second-order volumetric and deviatoric strain tensors, defined by the following equation:(4)εv=Iv:ε,(5)εd=Id:ε,(6)Iv=13I⊗I,(7)Id=Is−Iv,where Iv and Id are the fourth-order volumetric and deviatoric projections, respectively, and I and Is are the second-order identity tensor and the fourth-order symmetric identity tensor, respectively.ε
v
+ and εv− are positive and negative volumetric strain tensors, defined according to(8)εv+=Htrεεv,(9)εv−=H−trεεv.


### 2.2. Damage Criteria

To describe the damage due to tension or compression, it first necessitates a criterion that defines tensile/compressive states. As mentioned in Introduction, the method of spectral decomposition may lead to the premature convergence problem in numerical solutions [25]. Therefore, the sign of the volumetric deformation (relative volume change), which is calculated by trε, is used to distinguish tension and compression in this paper. The similar application can be found in the work of Comi and Perego [26]. In the following, the tensile and compressive strain states are thus defined by trε≥0 and trε<0, respectively.The damage criterion characterizes the state of damage and represents a bounding surface whose boundary corresponds to the states at which the damage will grow. To clearly specify the damage due to tension or compression, two damage criteria are introduced as follows:(10)Fiεieq,κi=εieqε−κi≤0,i=torc,where Ft and Fc are the tensile and compressive damage criteria, εteq and εceq are the equivalent tensile and compressive strains, and κt and κc are the current tensile and compressive damage thresholds, respectively. It is noteworthy that the damage does not develop if Fi<0 and grows if Fi=0.The damage threshold represents the largest value of the equivalent strain ever reached by the material during the loading history [65–67]. Therefore, κt and κc can be further expressed as follows:(11)κit=maxτ≤tεieqτ,i=torc.Similar to the work of Fichant et al. [68], the equivalent strains in equation (10) are assumed to be as follows:(12)εieq=1E0αiI1σ+βiJ2σ,i=tiftrε≥0,i=cotherwise,with(13)I1σ=trσe,J2σ=12σe:σe−16I1σ2,σe=ℂ0:ε,(14)αt=121−ftfc,βt=321+ftfc,αc=1−fcfbc,βc=32−fcfbc,where E0 is the initial Young’s modulus of the undamaged material, I1σ is the first invariant of the effective stress tensor, J2σ is the second invariant of the deviatoric effective stress tensor, σe is the effective stress tensor, ℂ0 is the fourth-order elastic stiffness tensor of the undamaged material, αi and βi are the material constants, ft is the uniaxial tensile strength, and fc and fbc are the uniaxial and biaxial compressive strengths of the material, respectively. Following the experimental result from the study of Kupfer et al. [69], the typical value of fbc/fc=1.16 for concrete is adopted throughout this paper.


### 2.3. Evolution Laws for Internal Variables

From the macroscopic viewpoint within the framework of thermodynamics, the damage variable can be introduced as an internal state variable to quantify the average material degradation [22]. In order to capture distinct tension-compression behaviors of concrete, two types of damage evolution laws, which relate the corresponding damage internal variable to the associated damage threshold, are defined in the following general form:(15)di=giκi,i=torc.In this paper, explicit damage evolution laws are further postulated as follows:(16)dt=gtκt=0,ifκt≤κt0,1−κt0κtκtu−κtκtu−κt0At,ifκt0<κt<κtu,1,ifκt≥κtu,(17)dc=gcκc=1−11+γ−1κc/κc0γ/γ−1,ifκc≤κc0,1−κc0γκc1−κc−κc0κcu−κc02,ifκc0<κc<κcu,1,ifκc≥κcu,where κt0=ft/E0 and κc0=γfc/E0 are the initial tensile and compressive damage thresholds, respectively, with γ being the material parameter. At is the nonnegative material parameter that controls the shape of the tensile stress-strain relation, as shown in Figure 1. It is noteworthy that if At=1, then the damage evolution law in tension represents a linear softening. κtu and κcu are the ultimate tensile and compressive damage thresholds, respectively. According to the crack band approach [33], the ultimate damage threshold depends on the fracture energy and the crack bandwidth. It will be further discussed later.Figure 1
Influence ofAt on uniaxial stress-strain relation in tension.At any material point, damage can evolve only if the current state reaches the boundary of the elastic domain, characterized byFiεieq,κi=0. To reflect the loading, unloading, and reloading conditions more generally, the evolution of damage internal variables is governed by the following Kuhn–Tucker relations:(18)Fiεieq,κi≤0,di•≥0,Fiεieq,κidi•=0,i=torc.


### 2.4. Mechanical Dissipation and Constitutive Laws

The failure of concrete is related to the degradation of physical properties of the material. During this damage process, the energy dissipation is always nonnegative and accompanied by an increase of entropy [66]. It hence implies an irreversible process of thermodynamics. Therefore, the mechanical dissipation ϕd under the isothermal condition needs to satisfy the second principle of thermodynamics, which can be expressed by the Clausius–Duhem inequality as follows:(19)ϕd•=σ:ε•−ρψ•≥0,where σ is the second-order stress tensor.Substitution of equation (1) into equation (19) yields the following equation:(20)ϕd•=σ−ρ∂ψ∂ε:ε•+Ytdt•+Ycdc•≥0,with the tensile and compressive damage energy release rates defined as follows:(21)Yi=−ρ∂ψ∂di=−12ε:∂ℂ∂di:ε,i=torc,where ℂ is the fourth-order secant stiffness tensor of the damaged material.Because of the fact that the inequality of equation (20) must be satisfied by the arbitrary change of the strain tensor ε, the coefficient of the first term of the inequality, therefore, must be zero. This requirement yields the following constitutive laws:(22)σ=ρ∂ψ∂ε,(23)ℂdt,dc=∂σ∂ε=ρ∂ψ2∂ε∂ε.Using equation (22), the Clausius–Duhem inequality expressed by equation (20) is eventually written in the following form:(24)ϕd•=Ytdt•+Ycdc•≥0.Therefore, the second principle of thermodynamics is satisfied as long as the damage variablesdt and dc are monotonically increasing since Yt and Yc are strictly nonnegative because of their quadratic forms, which involve ∂ℂ/∂di≤0.


### 2.5. Fracture Energy Regularization with Consistent Crack Bandwidths

To solve the problem of mesh-size dependency, a fracture energy regularization based on the crack band approach [33] is adopted in the present study. The general idea of the crack band approach is to rescale the softening branch of the stress-strain curve, in order to impose the same energy dissipation per unit area for different element sizes [59]. This adjustment consists in relating the fracture energy (regarded to be a material property) to the volumetric (or specific) fracture energy (Figures 2(c) and 2(d)), which depends on the crack bandwidth. To achieve this goal of regularization, it is finally required to relate the ultimate damage threshold to the fracture energy and the crack bandwidth.Figure 2
Damage zone and fracture energy regularization: (a) damage zone in the concrete specimen under tensile loading; (b) stress-strain response outside the damage zone; (c) traction-separation law in the damage zone; (d) general stress-strain relation for the concrete specimen.The crack band approach assumes that the failure of material occurs within a damage zone having a limited width, whereas the other parts are unloaded [59] during the process of failure (Figure 2). According to Bažant and Oh [33] and Jirásek and Marfia [43], the same energy dissipation per unit area can be ensured using following relations:(25)Gtf=ltPgtf=ltP∫0κtuσtnnκtdκt,(26)Gcf=lcPgcf=lcP∫κc0κcuσcnnκcdκc,where Gtf and Gcf are the Mode-I and compressive fracture energies defined as the energy dissipation per unit area of a surface due to cracking or crushing, ltP and lcP are the tensile and compressive crack bandwidths, gtf and gcf are the tensile and compressive volumetric fracture energies defined as the energy dissipation per unit volume, and σtnn and σcnn are the tensile and compressive stresses normal and tangential to the surface of cracking, respectively. It is noteworthy that Gtf and gtf can be represented by the area under the traction-separation (tensile stress normal to the crack versus crack opening displacement) curve and under the stress-strain curve, respectively, as shown in Figures 2(c) and 2(d).Using equation (16), substitution of equation (22) into equation (25) yields the following equation:(27)κtu=1+AtGtfltPft+1−At2κt0.It is noteworthy that equation (27) is derived under the condition of uniaxial loading. Nevertheless, the previous studies (e.g., [20, 70, 71]) have revealed that it is also an acceptable approximation under general loading conditions.Similarly, the compressive damage thresholdκcu can be expressed as follows:(28)κcu=32GcflcPfc+κc0.Generally, the crack bandwidth can be regarded as an FE discretization parameter, which depends on the chosen element shape (e.g., quadrilaterals or triangles), element size, element orientation (or mesh alignment) with respect to the direction of the strain localization band, interpolation function or order of the displacement approximation (e.g., linear or quadratic), and numerical integration scheme (number of integration points) [53, 59, 60]. For the sake of simplicity, only four-node quadrilateral and three-node triangular elements are used throughout the paper. It should also be noted that only low-order elements (i.e., elements with linear or multi-linear interpolation functions) are adopted in this study since the stress fields obtained from high-order elements (i.e., elements with quadratic interpolation functions) are highly disturbed for crack band simulations [59].Inspired by the works of Červenka et al. [57] and Jirásek and Bauer [59], this paper develops an improved estimation of crack bandwidths in tension and compression by projecting the element onto the directions normal and tangential to the crack direction, respectively. Mathematically, they can be written as follows:(29)liP=limax−limini=torc,with(30)ltmax=maxN=1,Nm−nxN⋅nx,ltmin=minN=1,Nm−nxN⋅nx,lcmax=maxN=1,Nm−nxN⋅tx,(31)lcmin=minN=1,Nm−nxN⋅tx,where ltP and lcP are the tensile and compressive crack bandwidths, limax and limin (i=torc) are the maximum and minimum projection values, xN is the midpoint of the element edge (for four-node quadrilateral elements) or the nodal point (for three-node triangular elements), Nm−n is the number of element midpoints (for four-node quadrilateral elements) or nodes (for three-node triangular elements), x represents the point (e.g., Gaussian point) in the structure where the crack bandwidth is to be evaluated, and nx and tx are the unit vectors with the direction normal and tangential to the crack direction at the point x, respectively. The geometrical meaning of equation (29) is illustrated in Figure 3.Figure 3
Proposed estimation of crack bandwidths.For comparison, the existing estimation formulas of the crack bandwidth mentioned in Introduction are recalled as follows.The tensile crack bandwidthltR, which is estimated by Rots [53] belonging to the first category, can be written as follows:(32)ltR=A,for  quadratic  two‐dimensional  elements,2A,for  linear  two‐dimensional  elements,with A being the area of the element.The tensile crack bandwidthltO, which is estimated according to the first method of Oliver [56], can be expressed as follows:(33)ltO=∑i=1Ncϕi∂Nix∂x⋅nx−1,with(34)ϕi=1,ifxi−xc⋅nx≥0,0,ifxi−xc⋅nx<0,where Nc is the number of corner nodes of the element, Ni is the standard shape function for the corner node i belonging to an element with Nc corner nodes, ϕi is the nodal value at the corner node i, and xi and xc are the corner node and the center of the element, respectively.The tensile crack bandwidthltG, which is estimated by Govindjee et al. [58] belonging to the third category, can be given as follows:(35)ltG=∑i=1Ncτi∂Nix∂x⋅nx−1,with(36)τi=xi−xc⋅nx−τminτmax−τmin,τmax=maxi=1,Ncxi−xc⋅nx,(37)τmin=mini=1,Ncxi−xc⋅nx.The geometrical meanings of the aforementioned tensile crack bandwidths and their differences with respect to the direction ofnx are illustrated in Figures 4 and 5, where θ is the angle between the normal to the crack and the global axis x at the point x. The crack bandwidths denoted by ltO (according to equation (33)) become very large (Figures 4(b) and 5(b)) under certain circumstances, and they also show a discontinuity (Figures 4(c) and 5(c)), while the crack bandwidths denoted by ltP and ltG (according to equations (29) and (35), respectively) are free from such problems. The discontinuous crack bandwidth obtained from equation (33) is caused by the discontinuous nature of ϕi, whose value is either 1 or 0 as defined by equation (34). It is also noteworthy that the present method (labelled as ltP) and the method of Govindjee et al. [58] (labelled as ltG) give exactly the same estimation of the tensile crack bandwidth for the linear triangular element as shown in Figures 5(a) and 5(c) since both of them use the method of projecting corner nodes of the element onto the direction normal to the crack direction. The further comparative study of these estimation methods will be given in the next section.Figure 4
Comparison of estimations of crack bandwidth for a rectangular quadrilateral element: (a) geometrical meanings of crack bandwidths estimated by the present model, Oliver [56], and Govindjee et al. [58]; (b) large crack bandwidth estimated by the first method of Oliver [56]; (c) relations of tensile crack bandwidth and angle θ.
(a)
(b)
(c)Figure 5
Comparison of estimations of crack bandwidth for a right triangular element: (a) geometrical meanings of crack bandwidths estimated by the present model, Oliver [56], and Govindjee et al. [58]; (b) large crack bandwidth estimated by the first method of Oliver [56]; (c) relations of tensile crack bandwidth and angle θ.
(a)
(b)
(c)


## 2.1. Helmholtz Free Energy Potential

The thermodynamics of irreversible processes provides a general framework for formulating constitutive equations. In other words, the state of the material can be described by the thermodynamic potential. Consequently, the stresses and other associated variables can be derived from this potential [64]. Because of the scope of this paper, the present model is developed based on the following assumptions: (i) damage evolutions are assumed to be isotropic and to occur in small strains and (ii) the rate-dependent behavior and the thermal effect are not considered. To take into account distinct tension-compression behaviors of concrete, the postulated thermodynamic potential in terms of Helmholtz free energy ψ is decomposed into volumetric and deviatoric parts ψv and ψd in a similar way as in the study of Cervera et al. [11] and Comi and Perego [26] as follows:(1)ρψε,dt,dc=ρψvεv,dt,dc+ρψdεd,dt,dc=ρψv+εv+,dt,dc+ρψv−εv−,dt,dc+ρψdεd,dt,dc,(2)ρψv+εv+,dt,dc=321−dt1−dcK0εv+:εv+,ρψv−εv−,dt,dc=321−dcK0εv−:εv−,(3)ρψdεd,dt,dc=1−Htrεdt1−H−trεdcG0εd:εd,where ρ is the mass density of the material, ε is the second-order strain tensor, dt and dc are the damage internal variables that represent the level of material degradation in tension and compression, K0 and G0 are the initial bulk and shear moduli of the undamaged material, respectively, H• is the Heaviside function defined by Hx=1ifx≥0andHx=0ifx<0, and tr• is the trace of a tensor. The tensile and compressive damage internal variables are monotonously increasing scalars, both of which grow from zero (undamaged material with the intact stiffness) to one (completely damaged material). εv and εd are second-order volumetric and deviatoric strain tensors, defined by the following equation:(4)εv=Iv:ε,(5)εd=Id:ε,(6)Iv=13I⊗I,(7)Id=Is−Iv,where Iv and Id are the fourth-order volumetric and deviatoric projections, respectively, and I and Is are the second-order identity tensor and the fourth-order symmetric identity tensor, respectively.ε
v
+ and εv− are positive and negative volumetric strain tensors, defined according to(8)εv+=Htrεεv,(9)εv−=H−trεεv.


## 2.2. Damage Criteria

To describe the damage due to tension or compression, it first necessitates a criterion that defines tensile/compressive states. As mentioned in Introduction, the method of spectral decomposition may lead to the premature convergence problem in numerical solutions [25]. Therefore, the sign of the volumetric deformation (relative volume change), which is calculated by trε, is used to distinguish tension and compression in this paper. The similar application can be found in the work of Comi and Perego [26]. In the following, the tensile and compressive strain states are thus defined by trε≥0 and trε<0, respectively.The damage criterion characterizes the state of damage and represents a bounding surface whose boundary corresponds to the states at which the damage will grow. To clearly specify the damage due to tension or compression, two damage criteria are introduced as follows:(10)Fiεieq,κi=εieqε−κi≤0,i=torc,where Ft and Fc are the tensile and compressive damage criteria, εteq and εceq are the equivalent tensile and compressive strains, and κt and κc are the current tensile and compressive damage thresholds, respectively. It is noteworthy that the damage does not develop if Fi<0 and grows if Fi=0.The damage threshold represents the largest value of the equivalent strain ever reached by the material during the loading history [65–67]. Therefore, κt and κc can be further expressed as follows:(11)κit=maxτ≤tεieqτ,i=torc.Similar to the work of Fichant et al. [68], the equivalent strains in equation (10) are assumed to be as follows:(12)εieq=1E0αiI1σ+βiJ2σ,i=tiftrε≥0,i=cotherwise,with(13)I1σ=trσe,J2σ=12σe:σe−16I1σ2,σe=ℂ0:ε,(14)αt=121−ftfc,βt=321+ftfc,αc=1−fcfbc,βc=32−fcfbc,where E0 is the initial Young’s modulus of the undamaged material, I1σ is the first invariant of the effective stress tensor, J2σ is the second invariant of the deviatoric effective stress tensor, σe is the effective stress tensor, ℂ0 is the fourth-order elastic stiffness tensor of the undamaged material, αi and βi are the material constants, ft is the uniaxial tensile strength, and fc and fbc are the uniaxial and biaxial compressive strengths of the material, respectively. Following the experimental result from the study of Kupfer et al. [69], the typical value of fbc/fc=1.16 for concrete is adopted throughout this paper.


## 2.3. Evolution Laws for Internal Variables

From the macroscopic viewpoint within the framework of thermodynamics, the damage variable can be introduced as an internal state variable to quantify the average material degradation [22]. In order to capture distinct tension-compression behaviors of concrete, two types of damage evolution laws, which relate the corresponding damage internal variable to the associated damage threshold, are defined in the following general form:(15)di=giκi,i=torc.In this paper, explicit damage evolution laws are further postulated as follows:(16)dt=gtκt=0,ifκt≤κt0,1−κt0κtκtu−κtκtu−κt0At,ifκt0<κt<κtu,1,ifκt≥κtu,(17)dc=gcκc=1−11+γ−1κc/κc0γ/γ−1,ifκc≤κc0,1−κc0γκc1−κc−κc0κcu−κc02,ifκc0<κc<κcu,1,ifκc≥κcu,where κt0=ft/E0 and κc0=γfc/E0 are the initial tensile and compressive damage thresholds, respectively, with γ being the material parameter. At is the nonnegative material parameter that controls the shape of the tensile stress-strain relation, as shown in Figure 1. It is noteworthy that if At=1, then the damage evolution law in tension represents a linear softening. κtu and κcu are the ultimate tensile and compressive damage thresholds, respectively. According to the crack band approach [33], the ultimate damage threshold depends on the fracture energy and the crack bandwidth. It will be further discussed later.Figure 1
Influence ofAt on uniaxial stress-strain relation in tension.At any material point, damage can evolve only if the current state reaches the boundary of the elastic domain, characterized byFiεieq,κi=0. To reflect the loading, unloading, and reloading conditions more generally, the evolution of damage internal variables is governed by the following Kuhn–Tucker relations:(18)Fiεieq,κi≤0,di•≥0,Fiεieq,κidi•=0,i=torc.


## 2.4. Mechanical Dissipation and Constitutive Laws

The failure of concrete is related to the degradation of physical properties of the material. During this damage process, the energy dissipation is always nonnegative and accompanied by an increase of entropy [66]. It hence implies an irreversible process of thermodynamics. Therefore, the mechanical dissipation ϕd under the isothermal condition needs to satisfy the second principle of thermodynamics, which can be expressed by the Clausius–Duhem inequality as follows:(19)ϕd•=σ:ε•−ρψ•≥0,where σ is the second-order stress tensor.Substitution of equation (1) into equation (19) yields the following equation:(20)ϕd•=σ−ρ∂ψ∂ε:ε•+Ytdt•+Ycdc•≥0,with the tensile and compressive damage energy release rates defined as follows:(21)Yi=−ρ∂ψ∂di=−12ε:∂ℂ∂di:ε,i=torc,where ℂ is the fourth-order secant stiffness tensor of the damaged material.Because of the fact that the inequality of equation (20) must be satisfied by the arbitrary change of the strain tensor ε, the coefficient of the first term of the inequality, therefore, must be zero. This requirement yields the following constitutive laws:(22)σ=ρ∂ψ∂ε,(23)ℂdt,dc=∂σ∂ε=ρ∂ψ2∂ε∂ε.Using equation (22), the Clausius–Duhem inequality expressed by equation (20) is eventually written in the following form:(24)ϕd•=Ytdt•+Ycdc•≥0.Therefore, the second principle of thermodynamics is satisfied as long as the damage variablesdt and dc are monotonically increasing since Yt and Yc are strictly nonnegative because of their quadratic forms, which involve ∂ℂ/∂di≤0.


## 2.5. Fracture Energy Regularization with Consistent Crack Bandwidths

To solve the problem of mesh-size dependency, a fracture energy regularization based on the crack band approach [33] is adopted in the present study. The general idea of the crack band approach is to rescale the softening branch of the stress-strain curve, in order to impose the same energy dissipation per unit area for different element sizes [59]. This adjustment consists in relating the fracture energy (regarded to be a material property) to the volumetric (or specific) fracture energy (Figures 2(c) and 2(d)), which depends on the crack bandwidth. To achieve this goal of regularization, it is finally required to relate the ultimate damage threshold to the fracture energy and the crack bandwidth.Figure 2
Damage zone and fracture energy regularization: (a) damage zone in the concrete specimen under tensile loading; (b) stress-strain response outside the damage zone; (c) traction-separation law in the damage zone; (d) general stress-strain relation for the concrete specimen.The crack band approach assumes that the failure of material occurs within a damage zone having a limited width, whereas the other parts are unloaded [59] during the process of failure (Figure 2). According to Bažant and Oh [33] and Jirásek and Marfia [43], the same energy dissipation per unit area can be ensured using following relations:(25)Gtf=ltPgtf=ltP∫0κtuσtnnκtdκt,(26)Gcf=lcPgcf=lcP∫κc0κcuσcnnκcdκc,where Gtf and Gcf are the Mode-I and compressive fracture energies defined as the energy dissipation per unit area of a surface due to cracking or crushing, ltP and lcP are the tensile and compressive crack bandwidths, gtf and gcf are the tensile and compressive volumetric fracture energies defined as the energy dissipation per unit volume, and σtnn and σcnn are the tensile and compressive stresses normal and tangential to the surface of cracking, respectively. It is noteworthy that Gtf and gtf can be represented by the area under the traction-separation (tensile stress normal to the crack versus crack opening displacement) curve and under the stress-strain curve, respectively, as shown in Figures 2(c) and 2(d).Using equation (16), substitution of equation (22) into equation (25) yields the following equation:(27)κtu=1+AtGtfltPft+1−At2κt0.It is noteworthy that equation (27) is derived under the condition of uniaxial loading. Nevertheless, the previous studies (e.g., [20, 70, 71]) have revealed that it is also an acceptable approximation under general loading conditions.Similarly, the compressive damage thresholdκcu can be expressed as follows:(28)κcu=32GcflcPfc+κc0.Generally, the crack bandwidth can be regarded as an FE discretization parameter, which depends on the chosen element shape (e.g., quadrilaterals or triangles), element size, element orientation (or mesh alignment) with respect to the direction of the strain localization band, interpolation function or order of the displacement approximation (e.g., linear or quadratic), and numerical integration scheme (number of integration points) [53, 59, 60]. For the sake of simplicity, only four-node quadrilateral and three-node triangular elements are used throughout the paper. It should also be noted that only low-order elements (i.e., elements with linear or multi-linear interpolation functions) are adopted in this study since the stress fields obtained from high-order elements (i.e., elements with quadratic interpolation functions) are highly disturbed for crack band simulations [59].Inspired by the works of Červenka et al. [57] and Jirásek and Bauer [59], this paper develops an improved estimation of crack bandwidths in tension and compression by projecting the element onto the directions normal and tangential to the crack direction, respectively. Mathematically, they can be written as follows:(29)liP=limax−limini=torc,with(30)ltmax=maxN=1,Nm−nxN⋅nx,ltmin=minN=1,Nm−nxN⋅nx,lcmax=maxN=1,Nm−nxN⋅tx,(31)lcmin=minN=1,Nm−nxN⋅tx,where ltP and lcP are the tensile and compressive crack bandwidths, limax and limin (i=torc) are the maximum and minimum projection values, xN is the midpoint of the element edge (for four-node quadrilateral elements) or the nodal point (for three-node triangular elements), Nm−n is the number of element midpoints (for four-node quadrilateral elements) or nodes (for three-node triangular elements), x represents the point (e.g., Gaussian point) in the structure where the crack bandwidth is to be evaluated, and nx and tx are the unit vectors with the direction normal and tangential to the crack direction at the point x, respectively. The geometrical meaning of equation (29) is illustrated in Figure 3.Figure 3
Proposed estimation of crack bandwidths.For comparison, the existing estimation formulas of the crack bandwidth mentioned in Introduction are recalled as follows.The tensile crack bandwidthltR, which is estimated by Rots [53] belonging to the first category, can be written as follows:(32)ltR=A,for  quadratic  two‐dimensional  elements,2A,for  linear  two‐dimensional  elements,with A being the area of the element.The tensile crack bandwidthltO, which is estimated according to the first method of Oliver [56], can be expressed as follows:(33)ltO=∑i=1Ncϕi∂Nix∂x⋅nx−1,with(34)ϕi=1,ifxi−xc⋅nx≥0,0,ifxi−xc⋅nx<0,where Nc is the number of corner nodes of the element, Ni is the standard shape function for the corner node i belonging to an element with Nc corner nodes, ϕi is the nodal value at the corner node i, and xi and xc are the corner node and the center of the element, respectively.The tensile crack bandwidthltG, which is estimated by Govindjee et al. [58] belonging to the third category, can be given as follows:(35)ltG=∑i=1Ncτi∂Nix∂x⋅nx−1,with(36)τi=xi−xc⋅nx−τminτmax−τmin,τmax=maxi=1,Ncxi−xc⋅nx,(37)τmin=mini=1,Ncxi−xc⋅nx.The geometrical meanings of the aforementioned tensile crack bandwidths and their differences with respect to the direction ofnx are illustrated in Figures 4 and 5, where θ is the angle between the normal to the crack and the global axis x at the point x. The crack bandwidths denoted by ltO (according to equation (33)) become very large (Figures 4(b) and 5(b)) under certain circumstances, and they also show a discontinuity (Figures 4(c) and 5(c)), while the crack bandwidths denoted by ltP and ltG (according to equations (29) and (35), respectively) are free from such problems. The discontinuous crack bandwidth obtained from equation (33) is caused by the discontinuous nature of ϕi, whose value is either 1 or 0 as defined by equation (34). It is also noteworthy that the present method (labelled as ltP) and the method of Govindjee et al. [58] (labelled as ltG) give exactly the same estimation of the tensile crack bandwidth for the linear triangular element as shown in Figures 5(a) and 5(c) since both of them use the method of projecting corner nodes of the element onto the direction normal to the crack direction. The further comparative study of these estimation methods will be given in the next section.Figure 4
Comparison of estimations of crack bandwidth for a rectangular quadrilateral element: (a) geometrical meanings of crack bandwidths estimated by the present model, Oliver [56], and Govindjee et al. [58]; (b) large crack bandwidth estimated by the first method of Oliver [56]; (c) relations of tensile crack bandwidth and angle θ.
(a)
(b)
(c)Figure 5
Comparison of estimations of crack bandwidth for a right triangular element: (a) geometrical meanings of crack bandwidths estimated by the present model, Oliver [56], and Govindjee et al. [58]; (b) large crack bandwidth estimated by the first method of Oliver [56]; (c) relations of tensile crack bandwidth and angle θ.
(a)
(b)
(c)


## 3. Model Verification

To establish the validity of the proposed constitutive model, the computed results were either verified theoretically or compared with measured ones from experiments. These are presented in the following sections.


### 3.1. Analysis of Biaxial Tests of Concrete

To examine the rationality of the adopted criterion distinguishing tension and compression, and to validate the present model in reproducing the distinct strength features of concrete under biaxial stresses, a series of plain concrete specimens subjected to different stress ratios (denoted byα), which were tested by Kupfer et al. [69], were chosen for the numerical analysis. The material properties used for the computation were Young’s modulus E0=3.2×104 MPa, Poisson’s ratio ν0=0.2, concrete uniaxial tensile strength ft=3.28 MPa, At=1.0, uniaxial compressive strength fc=32.8 MPa, γ=2.3, Mode-I fracture energy Gtf=0.1 N/mm, and compressive fracture energy Gcf=15.0 N/mm. A one-element mesh (Mesh A was 200.0 mm long × 200.0 mm high) with a four-node plane stress quadrilateral element using a 2 × 2 Gaussian integration scheme [44] was used to represent the plain concrete specimen (Figure 6).Figure 6
Biaxial strength envelope of concrete.Figure6 shows the comparison between the experimental and computed results of the strength envelope, which corresponds to the envelope of the maximum stress ever reached in the material during the loading, in a biaxial principal stress space. The results indicate that the predicted biaxial strengths of concrete are very close to the experimental ones, either in the biaxial tensile domain or in the biaxial compressive domain, or even in the tensile-compressive domain. The numerical results in Figure 6 also demonstrate the rationality of the adopted criterion distinguishing tension and compression. Figure 7 shows the comparison of stress-strain curves of numerical and experimental results [69] under biaxial compression (σ3=0), according to the following stress ratios (α=σ2/σ1): α=0, α=0.52, and α=1. It shows that the proposed model can predict well both the strength enhancement of concrete and the corresponding peak strain in biaxial compressive states.Figure 7
Computed stress-strain curves of concrete under biaxial compression.


### 3.2. Analysis of a Double-Edge Notched Concrete Specimen under Direct Tension

To validate the proposed model in capturing the strain-softening, a double-edge notched (DEN) concrete specimen under direct tension, which was investigated experimentally by Hordijk [72], was used for the numerical analysis. The test specimen shown in Figure 8, which had dimensions 125 mm height × 60 mm width × 50 mm thickness with two symmetrical notches of 5 × 5 mm2 in the middle cross-section, was glued to the loading platens and subjected to an axial tensile load under the deformation control. The average deformation δ, which denotes the elongation of the specimen, was measured experimentally by using four pairs of extensometers located at the front and rear faces of the specimen. The gauge length of extensometers was 35 mm (Figure 8). The average stress σavg in the experiment was defined as the vertical load divided by the net cross-sectional area (50 × 50 mm2).Figure 8
DEN concrete specimen: geometry, boundary, and loading conditions (dimensions in mm).The material properties of the test specimen used for the numerical analysis were Young’s modulusE0=1.8×104 MPa, Poisson’s ratio ν0=0.2, concrete uniaxial tensile strength ft=3.5 MPa, At=6.0, uniaxial compressive strength fc=34.0 MPa, and Mode-I fracture energy Gtf=0.07 N/mm.The DEN concrete specimen was modeled by a 298-element mesh with four-node plane stress quadrilateral elements using a 2 × 2 Gaussian integration scheme, as shown in Figure9. According to the boundary conditions of the experiment, all FE nodes at the lower boundary were assumed to be fixed (Figure 9). Figure 10 presents the computed average stress-deformation curve in tension. The numerical result shows a satisfactory agreement with the experimental one, including both the load-carrying capacity and the postpeak behavior.Figure 9
FE mesh for the DEN concrete specimen.Figure 10
Average stress-deformation curves for the DEN concrete specimen.


### 3.3. Analysis of a Bar Subjected to Uniaxial Tension Using Elements with a Constant Strain Field

To prove the efficiency of the proposed estimation of the crack bandwidth with respect to the FE mesh irregularity, a bar subjected to uniaxial tension was used for the numerical analysis by following the work of Oliver [56]. Figure 11 shows the geometry, boundary, and loading conditions of the bar with a thickness of 50 mm.Figure 11
A bar subjected to uniaxial tension.The material properties used for the numerical analysis were Young’s modulusE0=2.0×104 MPa, Poisson’s ratio ν0=0.0, concrete uniaxial tensile strength ft=2.0 MPa, At=1.0, uniaxial compressive strength fc=20.0 MPa, and Mode-I fracture energy Gtf=0.1 N/mm.According to Oliver [56], the localization of the bar due to strain-softening was induced artificially by the prescribed imperfection (Figure 11) with a small reduction of the uniaxial tensile strength in shaded elements, as shown in Figure 12.Figure 12
FE meshes used for the study of mesh irregularity (after the study [56]): (a) Mesh A; (b) Mesh B; (c) Mesh C; (d) Mesh D.
(a)
(b)
(c)
(d)The bar was discretized by the four-node quadrilateral or three-node triangular plane stress elements, both of which adopted a linear interpolation function (implying a constant strain field) and a 1-point Gaussian integration scheme. Figure12 shows the meshes used for the analysis with different mesh irregularities.Due to the fact that the projection methods (e.g., the present method and that of [58]) have solved the difficulty of the excessively large estimation of the crack bandwidth arising in the first method of Oliver [56] shown in the previous section, only the element area-based, Govindjee et al. [58]-based, and the present methods of estimating the crack bandwidth are compared in this section. Figure 13(c) shows that the load-displacement responses obtained by the present method are almost coincident with the theoretical solution and are hardly detectable in the plots, while the corresponding responses, which were obtained by the crack bandwidth estimated by 2A and by Govindjee et al. [58] using equation (35), show an unacceptable mesh dependency with respect to the FE mesh irregularity. Consequently, the present model provides an accurate estimation of the crack bandwidth that yields mesh-size independent results with respect to the mesh irregularity.Figure 13
Force-displacement responses obtained by different estimations of crack bandwidth with different mesh irregularities: crack bandwidth estimated by (a)2A, (b) Govindjee et al. [58], and (c) the present model.
(a)
(b)
(c)


### 3.4. Analysis of a Single-Edge Notched Concrete Beam Subjected to Four-Point Shear Loading

To demonstrate the mesh-size objectivity and the capability of the proposed model in predicting the structural response in shear loading, a single-edge notched (SEN) concrete beam of Arrea and Ingraffea [73] was used for the numerical verification. The geometry of the SEN concrete beam, boundary, and four-point shear loading conditions, which were used to induce a mixed tensile/shear failure, are shown in Figure 14. The load was applied using a steel beam, through which the load was transmitted to the SEN concrete beam at points A and B (Figure 14).Figure 14
SEN concrete beam: geometry, boundary, and loading conditions (dimensions in mm).The material properties of the SEN beam used in the analysis were Young’s modulusE0=2.48×104 MPa, Poisson’s ratio ν0=0.18, concrete uniaxial tensile strength ft=3.5 MPa, At=1.0, uniaxial compressive strength fc=35.0 MPa, γ=2.2, Mode-I fracture energy Gtf=0.075 N/mm, and compressive fracture energy Gcf=15.0 N/mm. The steel plates used as load distributors were assumed to be a linear elastic material with Young’s modulus Es=210 GPa and Poisson’s ratio νs=0.3.The SEN concrete beam and the steel plates were discretized by four-node quadrilateral (with 2 × 2 Gaussian integration) and three-node triangular plane stress elements in the FE analysis, respectively. Three different mesh sizes were used for the mesh-size dependency study: coarse, medium, and fine meshes, which had a total of 480, 810, and 1176 elements, respectively (Figure15). The corresponding minimum element sizes of these three meshes were 19.1 mm, 12.8 mm, and 9.6 mm, respectively.Figure 15
Different FE mesh sizes for the SEN concrete beam in the mesh-size dependency study: (a) coarse mesh; (b) medium mesh; (c) fine mesh.
(a)
(b)
(c)Figure16 presents the predicted structural responses in terms of the load P at point B versus the crack mouth sliding displacement (CMSD). Figure 17 shows the simulated final tensile damage distribution, which is denoted by the dark color corresponding to a high damage value close to one. Figure 18 shows the simulated crack patterns at the final computation stage.Figure 16
Load-CMSD curves for SEN concrete beam obtained with different FE mesh sizes.Figure 17
Final tensile damage profiles in the central part of the SEN concrete beam: (a) coarse mesh; (b) medium mesh; (c) fine mesh.Figure 18
Computed final crack patterns in the central part of the SEN concrete beam: (a) coarse mesh; (b) medium mesh; (c) fine mesh.
(a)
(b)
(c)According to the numerical results, it is clear that (i) the global structural responses are well predicted compared to the experimental results (Figure16); (ii) the proposed model provides an acceptable prediction of the curved macroscopic crack, as compared with the experimentally observed crack trajectory (Figures 17 and 18); and (iii) the simulated crack trajectories obtained with the medium and fine meshes are identical, which indicates that the predicted paths of the crack propagation are almost mesh-size independent.The predicted thicknesses of the macroscopic crack are inconsistent when using different mesh sizes, as shown in Figures17 and 18. As more clearly shown in Figure 18, the simulated localization band tends to occur in one row of elements. This is due to the inherent deficiency of the crack band approach since the rate BVP still remains ill-posed [63]. The similar conclusion to simulations with low-order elements can be found in the study of Jirásek and Bauer [59], who also pointed out that the minimum size of the localization band is the size of one row of integration points (i.e., subelement domains) but not the size of one row of elements in simulations with high-order elements.As also can be observed in Figures17 and 18, the predicted crack trajectories are less curved than its experimental counterpart due to the approximation error of the spatial discretization in the standard finite element method. It can be explained by the fact that the localized damage zone cannot be resolved with sufficient accuracy because of the minimum element size (19.1 mm) of the coarse mesh being too large, as compared with the thickness of the notch (5 mm).


## 3.1. Analysis of Biaxial Tests of Concrete

To examine the rationality of the adopted criterion distinguishing tension and compression, and to validate the present model in reproducing the distinct strength features of concrete under biaxial stresses, a series of plain concrete specimens subjected to different stress ratios (denoted byα), which were tested by Kupfer et al. [69], were chosen for the numerical analysis. The material properties used for the computation were Young’s modulus E0=3.2×104 MPa, Poisson’s ratio ν0=0.2, concrete uniaxial tensile strength ft=3.28 MPa, At=1.0, uniaxial compressive strength fc=32.8 MPa, γ=2.3, Mode-I fracture energy Gtf=0.1 N/mm, and compressive fracture energy Gcf=15.0 N/mm. A one-element mesh (Mesh A was 200.0 mm long × 200.0 mm high) with a four-node plane stress quadrilateral element using a 2 × 2 Gaussian integration scheme [44] was used to represent the plain concrete specimen (Figure 6).Figure 6
Biaxial strength envelope of concrete.Figure6 shows the comparison between the experimental and computed results of the strength envelope, which corresponds to the envelope of the maximum stress ever reached in the material during the loading, in a biaxial principal stress space. The results indicate that the predicted biaxial strengths of concrete are very close to the experimental ones, either in the biaxial tensile domain or in the biaxial compressive domain, or even in the tensile-compressive domain. The numerical results in Figure 6 also demonstrate the rationality of the adopted criterion distinguishing tension and compression. Figure 7 shows the comparison of stress-strain curves of numerical and experimental results [69] under biaxial compression (σ3=0), according to the following stress ratios (α=σ2/σ1): α=0, α=0.52, and α=1. It shows that the proposed model can predict well both the strength enhancement of concrete and the corresponding peak strain in biaxial compressive states.Figure 7
Computed stress-strain curves of concrete under biaxial compression.


## 3.2. Analysis of a Double-Edge Notched Concrete Specimen under Direct Tension

To validate the proposed model in capturing the strain-softening, a double-edge notched (DEN) concrete specimen under direct tension, which was investigated experimentally by Hordijk [72], was used for the numerical analysis. The test specimen shown in Figure 8, which had dimensions 125 mm height × 60 mm width × 50 mm thickness with two symmetrical notches of 5 × 5 mm2 in the middle cross-section, was glued to the loading platens and subjected to an axial tensile load under the deformation control. The average deformation δ, which denotes the elongation of the specimen, was measured experimentally by using four pairs of extensometers located at the front and rear faces of the specimen. The gauge length of extensometers was 35 mm (Figure 8). The average stress σavg in the experiment was defined as the vertical load divided by the net cross-sectional area (50 × 50 mm2).Figure 8
DEN concrete specimen: geometry, boundary, and loading conditions (dimensions in mm).The material properties of the test specimen used for the numerical analysis were Young’s modulusE0=1.8×104 MPa, Poisson’s ratio ν0=0.2, concrete uniaxial tensile strength ft=3.5 MPa, At=6.0, uniaxial compressive strength fc=34.0 MPa, and Mode-I fracture energy Gtf=0.07 N/mm.The DEN concrete specimen was modeled by a 298-element mesh with four-node plane stress quadrilateral elements using a 2 × 2 Gaussian integration scheme, as shown in Figure9. According to the boundary conditions of the experiment, all FE nodes at the lower boundary were assumed to be fixed (Figure 9). Figure 10 presents the computed average stress-deformation curve in tension. The numerical result shows a satisfactory agreement with the experimental one, including both the load-carrying capacity and the postpeak behavior.Figure 9
FE mesh for the DEN concrete specimen.Figure 10
Average stress-deformation curves for the DEN concrete specimen.


## 3.3. Analysis of a Bar Subjected to Uniaxial Tension Using Elements with a Constant Strain Field

To prove the efficiency of the proposed estimation of the crack bandwidth with respect to the FE mesh irregularity, a bar subjected to uniaxial tension was used for the numerical analysis by following the work of Oliver [56]. Figure 11 shows the geometry, boundary, and loading conditions of the bar with a thickness of 50 mm.Figure 11
A bar subjected to uniaxial tension.The material properties used for the numerical analysis were Young’s modulusE0=2.0×104 MPa, Poisson’s ratio ν0=0.0, concrete uniaxial tensile strength ft=2.0 MPa, At=1.0, uniaxial compressive strength fc=20.0 MPa, and Mode-I fracture energy Gtf=0.1 N/mm.According to Oliver [56], the localization of the bar due to strain-softening was induced artificially by the prescribed imperfection (Figure 11) with a small reduction of the uniaxial tensile strength in shaded elements, as shown in Figure 12.Figure 12
FE meshes used for the study of mesh irregularity (after the study [56]): (a) Mesh A; (b) Mesh B; (c) Mesh C; (d) Mesh D.
(a)
(b)
(c)
(d)The bar was discretized by the four-node quadrilateral or three-node triangular plane stress elements, both of which adopted a linear interpolation function (implying a constant strain field) and a 1-point Gaussian integration scheme. Figure12 shows the meshes used for the analysis with different mesh irregularities.Due to the fact that the projection methods (e.g., the present method and that of [58]) have solved the difficulty of the excessively large estimation of the crack bandwidth arising in the first method of Oliver [56] shown in the previous section, only the element area-based, Govindjee et al. [58]-based, and the present methods of estimating the crack bandwidth are compared in this section. Figure 13(c) shows that the load-displacement responses obtained by the present method are almost coincident with the theoretical solution and are hardly detectable in the plots, while the corresponding responses, which were obtained by the crack bandwidth estimated by 2A and by Govindjee et al. [58] using equation (35), show an unacceptable mesh dependency with respect to the FE mesh irregularity. Consequently, the present model provides an accurate estimation of the crack bandwidth that yields mesh-size independent results with respect to the mesh irregularity.Figure 13
Force-displacement responses obtained by different estimations of crack bandwidth with different mesh irregularities: crack bandwidth estimated by (a)2A, (b) Govindjee et al. [58], and (c) the present model.
(a)
(b)
(c)


## 3.4. Analysis of a Single-Edge Notched Concrete Beam Subjected to Four-Point Shear Loading

To demonstrate the mesh-size objectivity and the capability of the proposed model in predicting the structural response in shear loading, a single-edge notched (SEN) concrete beam of Arrea and Ingraffea [73] was used for the numerical verification. The geometry of the SEN concrete beam, boundary, and four-point shear loading conditions, which were used to induce a mixed tensile/shear failure, are shown in Figure 14. The load was applied using a steel beam, through which the load was transmitted to the SEN concrete beam at points A and B (Figure 14).Figure 14
SEN concrete beam: geometry, boundary, and loading conditions (dimensions in mm).The material properties of the SEN beam used in the analysis were Young’s modulusE0=2.48×104 MPa, Poisson’s ratio ν0=0.18, concrete uniaxial tensile strength ft=3.5 MPa, At=1.0, uniaxial compressive strength fc=35.0 MPa, γ=2.2, Mode-I fracture energy Gtf=0.075 N/mm, and compressive fracture energy Gcf=15.0 N/mm. The steel plates used as load distributors were assumed to be a linear elastic material with Young’s modulus Es=210 GPa and Poisson’s ratio νs=0.3.The SEN concrete beam and the steel plates were discretized by four-node quadrilateral (with 2 × 2 Gaussian integration) and three-node triangular plane stress elements in the FE analysis, respectively. Three different mesh sizes were used for the mesh-size dependency study: coarse, medium, and fine meshes, which had a total of 480, 810, and 1176 elements, respectively (Figure15). The corresponding minimum element sizes of these three meshes were 19.1 mm, 12.8 mm, and 9.6 mm, respectively.Figure 15
Different FE mesh sizes for the SEN concrete beam in the mesh-size dependency study: (a) coarse mesh; (b) medium mesh; (c) fine mesh.
(a)
(b)
(c)Figure16 presents the predicted structural responses in terms of the load P at point B versus the crack mouth sliding displacement (CMSD). Figure 17 shows the simulated final tensile damage distribution, which is denoted by the dark color corresponding to a high damage value close to one. Figure 18 shows the simulated crack patterns at the final computation stage.Figure 16
Load-CMSD curves for SEN concrete beam obtained with different FE mesh sizes.Figure 17
Final tensile damage profiles in the central part of the SEN concrete beam: (a) coarse mesh; (b) medium mesh; (c) fine mesh.Figure 18
Computed final crack patterns in the central part of the SEN concrete beam: (a) coarse mesh; (b) medium mesh; (c) fine mesh.
(a)
(b)
(c)According to the numerical results, it is clear that (i) the global structural responses are well predicted compared to the experimental results (Figure16); (ii) the proposed model provides an acceptable prediction of the curved macroscopic crack, as compared with the experimentally observed crack trajectory (Figures 17 and 18); and (iii) the simulated crack trajectories obtained with the medium and fine meshes are identical, which indicates that the predicted paths of the crack propagation are almost mesh-size independent.The predicted thicknesses of the macroscopic crack are inconsistent when using different mesh sizes, as shown in Figures17 and 18. As more clearly shown in Figure 18, the simulated localization band tends to occur in one row of elements. This is due to the inherent deficiency of the crack band approach since the rate BVP still remains ill-posed [63]. The similar conclusion to simulations with low-order elements can be found in the study of Jirásek and Bauer [59], who also pointed out that the minimum size of the localization band is the size of one row of integration points (i.e., subelement domains) but not the size of one row of elements in simulations with high-order elements.As also can be observed in Figures17 and 18, the predicted crack trajectories are less curved than its experimental counterpart due to the approximation error of the spatial discretization in the standard finite element method. It can be explained by the fact that the localized damage zone cannot be resolved with sufficient accuracy because of the minimum element size (19.1 mm) of the coarse mesh being too large, as compared with the thickness of the notch (5 mm).


## 4. Conclusions

This paper presents a tension-compression damage model with consistent crack bandwidths developed for the macroscopic description of concrete behaviors. The model is formulated within the framework of the internal variable theory of thermodynamics so that the constitutive relations are derived without ad hoc assumptions. Based on the numerical and experimental results, the following conclusions can be drawn:(i)
The essential features of the proposed model include the identification of tensile/compressive strain states by the sign of the volumetric deformation, the rational definition of two equivalent strains, two internal variables, and two separated damage criteria describing distinct damage mechanics under tension and compression. Such features enable the present model to provide a better prediction of the distinct tension-compression behaviors of concrete materials. Furthermore, the adopted method using the sign of the volumetric deformation avoids the premature convergence difficult in the numerical analysis, which is often encountered in existing models using the spectral decomposition of a stress or strain tensor.(ii)
The main contribution of this paper is that it provides an improved estimation of the crack bandwidth, which is developed for two-dimensional linear elements by taking into account the element shape, element size, and element orientation with respect to crack directions and spatial positions. Its capability in numerical analyses of irregular meshes and arbitrary crack directions has been validated by the comparative study of existing estimations of the crack bandwidth at the element level. The proposed model incorporating this enhanced description of the crack bandwidth has also been proven to be mesh-size independent at the structural level. The present model is able to provide objective simulation results with respect to the finite element mesh refinements that are very important for the computation of the real engineering structures.(iii)
The proposed model has been validated by the well-documented experimental data of concrete. The good agreement between the predicted and experimental results consequently demonstrates that the present model provides an effective tool for modeling concrete behaviors in tension and compression.(iv)
The present model has a limited number of material parameters, each of which has a clear physical meaning and can be identified by standard laboratory tests. Moreover, the simplicity and the explicit integration of the constitutive equations make the model well suited for large-scale computations.



---

*Source: 2810108-2019-02-07.xml*