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

Tóm tắt A mesoscale numerical approach to predict damage behavior in concrete basing on phase field method: ...cedures. 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 irre...ve 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 equatio...ening 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 macro...

**A mesoscale numerical approach to predict damage behavior in concrete basing on phase field method**, để tải tài liệu về máy bạn click vào nút DOWNLOAD ở trên

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.

File đính kèm:

- a_mesoscale_numerical_approach_to_predict_damage_behavior_in.pdf