A mesoscale numerical approach to predict damage behavior in concrete basing on phase field method

he damage level.
d = 1 represents the fully broken state of the material, d = 0 describes the unbroken state.
As a result, the second term in Eq. (1) can be expressed approximately by∫
Γ
gcdS '
∫
Ω
gcγd(d,∇d)dΩ, (2)
where γd(d,∇d) denotes the crack density function. According to [4, 7], it is defined by
γd(d,∇d) = 12l d
2 +
l
2
∇d · ∇d, (3)
l is a regularization parameter related to the width of the smeared crack. A sharp crack
is revealed as l approaches the value of zero. It is worth mentioning that, l is also an
internal parameter that affects the critical load of crack initiation [8,9]. In order to couple
the fracture phase field with the deformation problem, the energy density function can
be rewritten as
ψ(ε, d) = g(d)ψ(ε), (4)
where g(d) is a parabolic degradation function: g(d) = (1 − d)2 + k and k is a small
numerical parameter to avoid loss of stability in case of fully damaged elements. As a
consequence, the total energy in Eq. (1) reads as
E =
∫
Ω
ψ(ε, d)dΩ+
∫
Ω
gcγd(d,∇d)dΩ. (5)
In order to account for a stress degradation only in tension, the elastic strain density is
decomposed into positive part due to tension and negative part due to compression
ψ(ε, d) = g(d)ψ+(ε) + ψ−(ε). (6)
There are several methods which provide the split decomposition of strain energy [7, 9].
In the present work, the strain decomposition proposed by Miehe et al. [7] is adopted.
The elastic strain is decomposed into extensive ε+ and compressive ε− parts, with
ε = ε+ + ε−, (7)
and
ψ±(ε) =
λ
2
[〈Tr(ε)〉±]2 + µTr{(ε±)2}, (8)
where
ε+ =
D
∑
i=1
〈εi〉+ni ⊗ ni ε− =
D
∑
i=1
〈εi〉−ni ⊗ ni, (9)
where εi and ni are the eigenvalues and eigenvector of ε, i.e satisfying εni = εini. In
Eq. (9), 〈x〉+ = (x + |x|)/2 and 〈x〉− = (x− |x|)/2.
46 Nguyen Hoang Quan, Tran Bao Viet, Nguyen Thanh Tung
Besides, the external energy can be formulated as follows
Wext =
∫
Ω
f · udΩ+
∫
δΩF
F · u dΓ, (10)
where f and F are respectively the prescribed volume and boundary force.
The displacement field u(x) and the phase field d(x) can be determined by applying
the principle of maximum dissipation and energy minimization [5]. The problem can
be split into two quasi independent minimization procedures. First, with a fixed u, the
phase field problem is derived by minimizing the total energy with respect to the phase
field
d(x) = Arg
{
inf
d∈Sd
E(u, d)
}
, (11)
where Sd =
{
∇d · n = 0 on ∂Ω, d|d˙(x) ≥ 0, 0 ≤ d(x) ≤ 1
}
. The evolution law of the
phase field to ensure the irreversibility of the process is derived through a thermody-
namically consistent framework, see e.g for more details [4,10,11]. The weak form of the
phase field problem can be written as∫
Ω
{(
2ψ+ +
gc
l
)
dδd + gcl∇d · ∇(δd)
}
dΩ =
∫
Ω
2ψ+δddΩ. (12)
The history strain energy density functionH(x, t) is introduced to describe a depen-
dence on history and possible loading-unloading
H(x, t) = max
τ∈[0,t]
{
ψ+(x, τ)
}
. (13)
The weak form of the phase field problem is finally rewritten as∫
Ω
{(
2H+ gc
l
)
dδd + gcl∇d · ∇(δd)
}
dΩ =
∫
Ω
2HδddΩ. (14)
Then with a fixed d, the mechanical problem consists in minimizing the total energy
with respect to displacements:
u(x) = Arg
{
inf
u∈Su
(
E(u, d)−Wext
)}
, (15)
where Su = {u|u(x) = u on ∂Ωu, u ∈ H1(Ω)}. We obtain the weak form for u(x) ∈ Su
as follows ∫
Ω
σ : ε(δu)dΩ =
∫
Ω
f · δudΩ+
∫
∂ΩF
FδudΓ ∀δu ∈ H10(Ω), (16)
where the Cauchy stress σ =
∂ψ
∂ε
is given by
σ =
(
(1− d)2 + k
){
λ〈Trε〉+1+ 2µε+
}
+ λ〈Trε〉−1+ 2µε−. (17)
A mesoscale numerical approach to predict damage behavior in concrete basing on phase field method 47
2.2. Finite element discretization
The problem described in Eqs. (14), (16) are solved by a standard FE procedure
in a staggered scheme at each time step, i.e. the phase field problem and the mechanical
problem are solved alternatively. For more theoretical and practical details, the interested
reader can refer to Miehe et al. [7] and Nguyen et al. [4].
In 2D, the vector form for second order tensor can be expressed as: [ε] = {ε11, ε22, 2ε12}T
and [σ] = {σ11, σ22, σ12}T, and [1] = {1, 1, 0}T. The discretization of the system of the
governing equations at element level using the FEM for displacement and phase field
variables can be expressed as follows
u = Nuui, u = Nuδui,
[ε(u)] = Buui, [ε(δu)] = Buδui,
d(u) = Nd(x)di, ∇d(x) = Bd(x)di,
δd(u) = Nd(x)δdi, ∇δd(x) = Bd(x)δdi,
(18)
where ui, di denoting nodal displacements and nodal phase field at time tn+1, respec-
tively. Nu,d, Bu,d are vector of shape function and matrix of shape functions derivatives
for displacement and phase field, respectively.
The finite element equation of phase field problem is given by
Kddn+1 = Fd, (19)
where
Kd =
∫
Ω
{(
gc
l
+ 2Hn
)
NTd Nd + gclB
T
d Bd
}
dΩ, (20)
and
Fd =
∫
Ω
2NTdHndΩ, (21)
whereHn = H(un) is computed from the previous load increment{ Hn+1(x) = ψ+n+1(x) if ψ+n+1(x) > ψ+n (x)
Hn+1(x) = ψ+n (x) if ψ+n+1(x) ≤ ψ+n (x)
(22)
For the mechanical problem, the spectral decomposition of the strain field
(Eqs. (7), (8)) cause a strongly nonlinear mechanical problem. To avoid this nonlinear-
ity, the shifted strain tensor split algorithms proposed by the present authors in Ref [4]
is adopted. Within the context of incremental scheme, the projection tensors defined at
time n + 1, will be evaluated based on the result from previous time step n as follows
ε±n+1 ' P±(εn) : εn+1,
〈Trεn+1〉+ ' R+(εn)Trεn+1, 〈Trεn+1〉− ' R−(εn)Trεn+1,
(23)
48 Nguyen Hoang Quan, Tran Bao Viet, Nguyen Thanh Tung
with R+(εn) = 12
(
sign(Trεn) + 1
)
; R−(εn) = 12
(
sign(−Trεn) + 1
)
. Setting R±(εn) ≡
R±n , P±(εn) ≡ P±n , where P± are the matrix forms associated with the fourth-order ten-
sors P±. Then the stress at time tn+1 can be expressed as
[σn+1] =
(
(1− dn+1)2 + k
){
λR+n ([εn+1] · [1])[1] + 2µP+n [εn+1]
}
+ λR−n ([εn+1] · [1])[1] + 2µP−n [εn+1].
(24)
The finite element equation for the mechanical problem can be written as follows
{K1(dn+1, un) +K2(un))}un+1 = Fn+1, (25)
where
K1(dn+1) =
∫
Ω
BTu
{(
(1− dn+1)2 + k
)(
λR+n [1]T[1] + 2µP+n
)}
BudΩ, (26)
K2 =
∫
Ω
BTu
{
λR+n [1]T[1] + 2µP+n
}
BudΩ, (27)
Fn+1 =
∫
Ω
NTu fdΩ+
∫
∂ΩF
NTu FdΓ. (28)
3. ARRANGEMENT OF AGGREGATE PARTICLE BASED ON
MONTE CARLO SIMULATION
At mesoscopic level, concrete could be represented as biphasic material: coarse ag-
gregates and mortar matrix and an interfacial transition zone (ITZ) between them. The
evaluation of the composite behavior of concrete at mesoscopic level requires the gener-
ation of a random aggregate structure in which the shape, size and distribution of coarse
aggregate closely resemble real concrete in the statistical sense. The shape of aggregate
particles depends on the aggregate types. In general, gravel aggregates have a rounded
shape while crushed stone aggregates have an angular shape. In 2D numerical simula-
tion, the aggregate shape could be simulated by a polygonal shape [12] and elliptical or
circular shape [13].
The size distribution of concrete may be constructed based on an experimental siev-
ing process. Alternatively, the grading of aggregate particle is designed by the Fuller
curve which give an optimal density and strength of concrete mixture. The Fuller curve
can be expressed as follows
P(d) = 100
(
d
dmax
)n
, (29)
where P(d) is the cumulative percentage passing a sieve with aperture diameter d. dmax
is the and maximum size of aggregate and n is the exponent of the equation. Thus, for an
interval [di, di+1] defined by two sequential sieve opening sizes, di and di+1, then the area
of aggregates within a grading segment [di, di+1] can be calculated as:
Aagg[di+1 − di] = P(di+1)− P(di)P(dmax)− P(dmin)Pagg A, (30)
A mesoscale numerical approach to predict damage behavior in concrete basing on phase field method 49
where Aagg[di+1 − di] is the area of aggregate within the grading segment [di, di+1]. dmin
is the minimum size of aggregate, Pagg is the area fraction of all aggregates and A is the
total size of the concrete specimen.
Regarding the simulation of the aggregate spatial distribution, the random sampling
principle of Monte Carlo’s simulation method is used. This method is commonly called
the take-and-place method. The random principle is applied by taking the aggregate
sizes from a grading curve and placing each particle in the mortar matrix randomly so
that intersection between aggregate is avoided. This method has been used by most
researchers including Bazant et al. [14], Schlangen and Van Mier [15]. A different method
can be used such as: the divide and fill method [16], the random particle drop method
[17].
4. APPLICATION
The main purpose of this numerical example to demonstrate the potential of the
phase field method to simulate the crack propagation in highly complex microstructure
of concrete. For this purpose, the tensile test in [18] is numerically analyzed and the
results are compared with the experimental one.
Fig. 2. Geometry of the specimen and
boundary condition
Fig. 2 shows the geometry and bound-
ary conditions for uniaxial tests. It consists of
50 mm × 50 mm square numerical specimens.
The model is fixed at the bottom boundary
and is subjected to a uniformly distributed dis-
placement at the top boundary. In this study,
the aggregate size distribution is generated by
using the Fuller curve. n is chosen equal to 0.5.
The aggregate particle whose size is greater
than 2.36 mm is considered as coarse aggregate
while fine aggregate together with cement ma-
trix is treated as mortar. The interfacial transi-
tion zone between coarse aggregates and mor-
tar matrix is neglected. Here, for the sake of
simplicity, the coarse aggregate particles are
geometrically represented by a circular shape.
For normal strength concrete, the coarse aggre-
gate represents around 40-50% the concrete volume. In this study, the area fraction of
coarse aggregate is assumed to be equal to 40%. Coarse aggregates and mortar are
described by linear elastic behavior. Similar material properties as in [12] are used in
this study. Young’s modulus is 70000 MPa for coarse aggregates and is 25000 MPa for
mortar. Poisson’s ratio of both coarse aggregates and mortar is 0.2. Fracture energy is
gc = 0.06 N/mm for coarse aggregates and is gc = 0.05 N/mm for mortar.
The analyses are performed in plane stress condition and the out of plane thick-
ness was unit. All analyses is ended at a displacement 0.08 mm. The computation is
performed with monotonic displacement increments of u = 10−4 mm. during 800 load
50 Nguyen Hoang Quan, Tran Bao Viet, Nguyen Thanh Tung
increments. The length scale parameter is chosen as l = 0.35 mm. In this study, the do-
main does not contain pre-existing cracks, we can not predict the crack nucleation and
the crack pattern. Thus, in order to detect the crack nucleation, the domain is meshed by
a regular triangular grid element whose characteristic size is about h ≈ 0.15 mm. The to-
tal number of element is 249560. In [19], it showed that an element size h < l/2 is needed
in order to resolve the regularized crack surface Γl(d), such that we have Γl(d) ≈ Γ in the
finite element approximation.
(a) Concrete sample 1 (b) Final crack pattern in concrete sample 1
(c) Concrete sample 2 (d) Final crack pattern in concrete sample 2
Fig. 3. Predicted final crack pattern with different aggregate distribution (sample 1 and sample 2)
Figs. 3(a), 3(c) represent two random generations of aggregates, called concrete sam-
ple 1 and concrete sample 2. The results obtained in terms of final crack patterns are
depicted in Figs. 3(b), 3(d). The cracks are represented by red color. Their correspond-
ing stress-displacement curves are plotted in Fig. 4 and are compared with the result
A mesoscale numerical approach to predict damage behavior in concrete basing on phase field method 51
Fig. 4. Comparison of stress-displacement curves in uniaxial tension test
obtained experimentally by Hordijk [18]. It can be seen that the numerical and experi-
mental results are in good agreement at the elastic stage. There is slightly difference in
the peak stress. However, the numerical results underestimate the softening branch. It
can be explained by the fact that some simplification are adopted in this study: the ITZ is
neglected, the aggregate shape is represented by circular shape. The concrete sample fails
either with one macrocracks (concrete sample 1) or with two crack (concrete sample 2).
All macrocracks are predominantly perpendicular to the load direction. It can be seen
that the post-peak stress obtained from concrete sample 1 drops more quickly than the
one from concrete sample 2. Thus, there is a lower disspated energy in concrete sample 1.
This behaviour may be attributed to smaller fracture area in single crack than two cracks.
(a) Initial stage of crack propagation processes (b) Final stage of crack propagation processes
Fig. 5. Crack propagation processes in uniaxial tension test with pre-existing notch
52 Nguyen Hoang Quan, Tran Bao Viet, Nguyen Thanh Tung
To futher verify the performance of the proposed model for fracture of concrete at
mesoscale, a uniaxial test with pre-existing notch with the length of 20 mm and the depth
of 1 mm is simulated. The same boundary condition, material properties in the aforemen-
tioned test is used. The crack propagation processes of concrete for pre-existing notch are
shown in Fig. 5. At the initial stage, the crack are started from the notch Fig. 5(a). The
crack orientation is consistent with the crack propagation orientation of mode I fracture.
In the subsequent propagation processes, relatively large aggregates encountered and
consequently, the crack path is changed Fig. 5(b). The crack runs around than propa-
gating through aggregates, which is in good agreement with the realistic behaviour of
normal concrete.
5. CONCLUSION
In this paper, we deal with the traditional ”but complex” problem of modeling the
damage and fracture behavior of concrete material. To do this, we construct an numer-
ical approximation based on the phase field thermodynamic framework. Then, we are
interested in only simple numerical tests composed of 2D configuration, circular aggre-
gate, and no ITZ phase. In spite of this simplicity, some numerical results show good
agreement with experimental observations, this is an interesting one and permit us to in-
vestigate this approach for further applications concerning the complex micro-structure
of cement-based composite material.
ACKNOWLEDGMENT
This research is supported by Ministry of Education and Training under the grant
number B2020-GHA-07.
REFERENCES
[1] G. Constantinides and F.-J. Ulm. The effect of two types of CSH on the elasticity of cement-
based materials: Results from nanoindentation and micromechanical modeling. Cement and
Concrete Research, 34, (1), (2004), pp. 67–80. https://doi.org/10.1016/s0008-8846(03)00230-8.
[2] P. Wriggers and S. O. Moftah. Mesoscale models for concrete: Homogenisation and
damage behaviour. Finite Elements in Analysis and Design, 42, (7), (2006), pp. 623–636.
https://doi.org/10.1016/j.finel.2005.11.008.
[3] N. Ile, X.-H. Nguyen, P. Kotronis, J. Mazars, and J. M. Reynouard. Shaking table tests of
lightly rc walls: Numerical simulations. Journal of Earthquake Engineering, 12, (6), (2008),
pp. 849–878. https://doi.org/10.1080/13632460801890430.
[4] T. T. Nguyen, J. Yvonnet, Q. Z. Zhu, M. Bornert, and C. Chateau. A phase field method
to simulate crack nucleation and propagation in strongly heterogeneous materials from di-
rect imaging of their microstructure. Engineering Fracture Mechanics, 139, (2015), pp. 18–39.
https://doi.org/10.1016/j.engfracmech.2015.03.045.
[5] G. A. Francfort and J. J. Marigo. Revisiting brittle fracture as an energy minimization
problem. Journal of the Mechanics and Physics of Solids, 46, (8), (1998), pp. 1319–1342.
https://doi.org/10.1016/s0022-5096(98)00034-9.
[6] X. Li, D. Chu, Y. Gao, and Z. Liu. Numerical study on crack propagation in linear elastic
multiphase composite materials using phase field method. Engineering Computations, 36, (1),
(2019), pp. 307–333. https://doi.org/10.1108/EC-03-2018-0116.
A mesoscale numerical approach to predict damage behavior in concrete basing on phase field method 53
[7] C. Miehe, M. Hofacker, and F. Welschinger. A phase field model for rate-independent
crack propagation: Robust algorithmic implementation based on operator splits. Com-
puter Methods in Applied Mechanics and Engineering, (199), (2010), pp. 2765–2778.
https://doi.org/10.1016/j.cma.2010.04.011.
[8] T. T. Nguyen, Y. Yvonnet, M. Bornert, C. C. Chateau, K. Sab, R. Romani, and B. Le Roy.
On the choice of parameters in the phase field method for simulating crack initiation
with experimental validation. International Journal of Fracture, 197, (2), (2016), pp. 213–226.
https://doi.org/10.1007/s10704-016-0082-1.
[9] H. Amor, J. J. Marigo, and C. Maurini. Regularized formulation of the variational brittle
fracture with unilateral contact: numerical experiments. Journal of the Mechanics and Physics
of Solids, 57, (8), (2009), pp. 1209–1229. https://doi.org/10.1016/j.jmps.2009.04.011.
[10] T. T. Nguyen, J. Yvonnet, Q.-Z. Zhu, M. Bornert, and C. Chateau. A phase field method for
computational modeling of interfacial damage interacting with crack propagation in realistic
microstructures obtained by microtomography. Computer Methods in Applied Mechanics and
Engineering, 312, (2016), pp. 567–595. https://doi.org/10.1016/j.cma.2015.10.007.
[11] T. T. Nguyen, J. Yvonnet, D. Waldmann, and Q. C. He. Phase field modeling of interfacial
damage in heterogeneous media with stiff and soft interphases. Engineering Fracture Mechan-
ics, (2019), pp. 106–547. https://doi.org/10.1016/j.engfracmech.2019.106574.
[12] C. M. Lo´pez, I. Carol, and A. Aguado. Meso-structural study of concrete fracture using in-
terface element I: numerical model and tensile behavior. Materials and Structures, 41, (2007),
pp. 583–599. https://doi.org/10.1617/s11527-007-9314-1.
[13] X. F. Wang, Z. J. Yang, J. R. Yates, A. P. Jivkov, and C. Zhang. Monte carlo
simulation of mesoscale fracture modelling of concrete with random aggre-
gates and pores. Construction and Building Materials, 15, (2015), pp. 35–45.
https://doi.org/10.1016/j.conbuildmat.2014.09.069.
[14] Z. P. Bazant, M. R. Tabbara, M. T. Kazemi, and G. Pijaudier-Cabot. Random particle model for
facture of aggregate or fiber composites. Journal of Engineering Mechanic, 116, (1990), pp. 1686–
1705. https://doi.org/10.1061/(asce)0733-9399(1990)116:8(1686).
[15] E. Schlangen and J. G. M. van Mier. Simple lattice model for numerical simulation of fracture
of concrete materials and structures. Materials and Structures, 25, (156), (1992), pp. 534–542.
https://doi.org/10.1007/bf02472449.
[16] G. D. Schutter and L. Taerwe. Random particle model for concrete based on De-
launay triangulation. Journal of Engineering Mechanic, 26, (156), (1993), pp. 1686–1705.
https://doi.org/10.1007/bf02472853.
[17] J. G. M. van Mier and M. R. A. V. Vliet. Influence of microstructure of concrete on size/s-
cale effects in tensile fracture. Engineering Fracture Mechanics, 70, (16), (2003), pp. 2281–2306.
https://doi.org/10.1016/s0013-7944(02)00222-9.
[18] D. A. Hordijk. Tensile and tensile fatigue behaviour of concrete: experiments, modelling and
analyses. Heron, 37, (1992), pp. 1–79.
[19] C. Miehe, M. Hofacker, and F. Welschinger. Thermodynamically consistent phase-
field models of fractures: variational principles and multi-field FE implementations.
International Journal for Numerical Method in Engineering, (83), (2010), pp. 1273–1311.
https://doi.org/10.1002/nme.2861.