Multi-Scale Damage Model for Quasi-Brittle Composite Materials

: In the present paper, a hierarchical multi-scale method is developed for the nonlinear analysis of composite materials undergoing heterogeneity and damage. Starting from the homogenization theory, the energy equivalence between scales is developed. Then accompanied with the energy based damage model, the multi-scale damage evolutions are resolved by homogenizing the energy scalar over the meso-cell. The macroscopic behaviors described by the multi-scale damage evolutions represent the mesoscopic heterogeneity and damage of the composites. A rather simple structure made from particle reinforced composite materials is developed as a numerical example. The agreement between the full-scale simulating results and the multi-scale simulating results demonstrates the capacity of the proposed model to simulate nonlinear behaviors of quasi-brittle composite materials within the multi-scale framework.


Introduction
Composite materials are widely studied and used due to their excellent mechanical and physical performances. However, the thorough analytical model that represents the salient nonlinear features of composite materials is still an open topic due to the complexity induced by the heterogeneity and its nonlinear evolution in meso-scale. The damage and failure of composite materials and structures are multi-scale in their nature. The full scale numerical simulation which considers the mesoscopic imperfections explicitly during the overall analysis of structure is commonly time-consuming and sometimes unstable. Furthermore, most of the information obtained for the fine scale is often unnecessary to teach us about the structural damage and failure. Thus in recent 20 years, the multiscale modeling of composite materials has drawn multitude attentions following the developments of computational mechanics (see Refs [Miehe, Schotte and Schroder (1999); Trykozko and Zijl (2002) ;Zimmermann, Kleinman and Hordijk (2003); To and Li (2005); Li, Liu, Agrawal et al. (2006); Li and Tong (2015)] and others). The purpose of multi-scale method is to simulate the overall behaviors of structure with appropriate considerations of mesoscopic/mesoscopic structures so that the cost is much less than the full scale analysis. The inter-scale connection developed based on the multi-scale method may bridge the gap between the mesoscopic imperfections and the macroscopic damage evolution, by which the multi-scale damage model could be established. Multi-scale methods have been studied for a long time in many scientific communities. Numbers of theories have been proposed to analyze various of systems. Relying on the degree to the inter-scale coupling, they could be classified into two groups [Rudd and Broughton (2000); Ghosh, Bai and Raghavan (2007)]. The hierarchical/sequential methods [Bakhvalov and Panasenko (1989); Guedes and Kikuchi (1990) ;Fish, Yu and Shek (1999); Chen and Mehraeen (2005); Ren, Chen, Li et al. (2011)] were proposed for the weak coupling systems. By condensing the inter-scale coupling into several parameters, the scales involved in the systems could be simulated in separate calculations conducted in sequence. For the strong coupled systems, the scales must be simulated concurrently. By combining different scales described by different models in different resolutions, the concurrent methods concentrate on developing smooth and effective transitions across the scales [Tadmor, Ortiz and Phillips (1996) ;Shilkrot, Miller and Curtin (2002); Wagner and Liu (2003); Xiao and Belytschko (2004); Ghosh, Bai and Raghavan (2007); Feng and Ren (2017); Ren and Li (2013)]. For the linear elastic problems, the hierarchical method is often implemented by the bottom-up homogenization which provides an extremely efficient way to evaluate the homogenized properties of composite materials. Then the structure could be simulated relying on the continuum model with the material parameters informed from the mesoscale. For the structures undergoing nonlinearity, localization and discontinuity, it is believed that the strong coupling among scales dominates the performance of systems. Thus the concurrent methods with real-time interactive communication among scales are considered to be suitable for these problems. However, the numerical simulation of large scale structures based on concurrent methods is computationally prohibitive at present. Thus the present paper aims to propose a feasible method for the multi-scale nonlinear analysis of composites structures. As the standard analytical tool for the degrading structures, the continuum damage mechanics is adopted for the analysis of structure in macro-level [Feng, Ren and Li (2018) ;Feng, Wang and Wu (2019) ;Feng, Xie, Deng et al. (2019)]. And the mesoscopic cell is developed to describe the mesoscopic nonlinearities and discontinuities. An energy based homogenization approach is proposed based on asymptotic method to construct the relationship between mesoscopic cells and macroscopic damaged continua. The damage initiation and evolution can be explained and quantified by the mesoscopic cells, and there-after can be used for the macroscopic continua. Although certain detailed information in meso-level may be lost due to the homogenization, the proposed hierarchical approach is able to simulate the multi-scale nonlinear behaviors of structures with minimum computational expenses. The proposed model is based on the concept of damage, which is feasible for the analysis of quasi-brittle composite materials. The present paper is organized as follows. Section 2 describes the model problem considered in the present theory and developes the equation of energy homogenization, which considers cracks and nonlinearities in the meso-scale. Based on the energy homogenization equation, the multi-scale damage evolutions are developed in Section 3. In Section 4, a rather simple structure made from particle reinforced composite materials is simulated as a numerical example. In the end, the conclusions are given in Section 5.
2 Model problem and homogenization of energy Figure 1: Two-scale structure of model problem As shown in Fig. 1, two scales are considered for the nonlinear structural simulations. In the macro-scale, the overall structure occupying domain Ω with boundary Γ is described as a continua. In the meso-scale, the heterogeneities and discontinuities of mesoscopic cell are considered. Besides mesoscopic cracksand aggregates in meso-scale are also considered to define the domain with different mechanical properties from the main body. Thus the damage evolution in macro-scale could be informed by the analytical results at meso-scale. And the inter-scale link could be developed based on the homogenization method. To develop the multi-scale analytical model ( Fig. 1), we define the overall coordinate by X = (X 1 , X 2 , X 3 ) in structural level and the local material coordinate by Y = (Y 1 , Y 2 , Y 3 ) for mesoscopic cell respectively. To address the geometric nonlinearities of meso-cell, two coordinates, i.e., Y and y, are introduced in meso-level. In the original coordinate Y , the meso-cell is denoted by Ω Y and Γ Y ; and in the current coordinate y, we have the denotations Ω y and Γ y . The mapping function between the original and the current coordinates is defines as (1) and the displacement for the material point is Consider the equilibrium in the current coordinate y expressed by Cauchy stress τ , we have div y τ = 0 with the natural boundary condition on the inner boundaries On the outer boundary of the meso-cell, we define the following essential boundary condition where ε is the homogenized strain of the meso-cell. Multiplying Eq. (3) by displacement u and integrating over Ω y yield By integration by parts, one obtains Consider the first term on the LHS of Eq. (7). By substituting Eq. (5) one obtains where the homogenized Cauchy stress Thus Eq. (7) further gives where the volume of the deformed meso-cell is In nonlinear continuum mechanics, we also adopt other energy conjugated stress-strain definitions for the convenience of analysis and simulation. Thus the volumetric integration in Eq. (10) could be also expressed by other stress-strain pairs.
where S and F · S are the first and second Piola-Kirchhoff stress tensors; E is Green strain tensor. And the integration domain is the original cell Ω Y or the deformed cell Ω y as appropriate. Define the generalized energy conjugated stress-strain pairs by Σ and Ξ, we have the generalized Helmholtz free energy as follows: And the homogenized Helmholtz free energy of the meso-cell is expressed as The energy integration in Eq. (10) finally converts to Eq. (15) indicates that the energy of homogenized material equals to the averaging energy within the mesoscopic cell, which is the extension of classic Hill's condition [Hill (1963)] by considering the continuum nonlinearities and discontinuities.
3 Energy based multi-scale damage evolution After years of development, the modern framework of damage theory is established based on irreversible thermodynamics with energy based damage definitions [Mazars and Pijaudier-Cabot (1989); Ju (1989); Faria, Oliver and Cervera (1998)].
The single scalar damage theory is often defined by the following equation: where d denotes the damage scalar; ψ is the homogenized HFE of damaged heterogenous material and ψ 0 is the energy density of undamaged material which could be easily expressed by the undamaged material undergoing the same strain with the damaged material. We have where ε is the scenod-order strain tensor and C is the fourth-order elastic modulus tensor. Following the standard procedure of conventional damage theory, Eq. (16) should be substituted into Clausius-Duhem inequality as follows: Then we obtain the constitutive law and the evolutional inequality where the damage energy release rate According to Refs [Ju (1989); Lemaitre and Desmorat (2005); Wu, Li and Faria (2006)] and other works, damage evolution should be driven by the damage energy release rate. Thus we have where r t is given by and the damage evolution function g(·) should be determined by experimentation according to the conventional damage mechanics. The present paper proposed that the damage evolution could be easily resolved by Eq. (16) as follows: And we have proven that the homogenized HFE ψ could be determined by Eq. (15) based on the simulated results of the mesoscopic cell. Hence the damage evolution function g(·) could be resolved by Eq. (24) with the information from mesoscopic cracking and nonlinearity.
To describe the damage modes of quasi-brittle material undergoing different loadings conditions, the bi-scalar damage models [Faria, Oliver and Cervera (1998); Wu, Li and Faria (2006); Feng, Wu, Sun et al. (2017)] were proposed with two damage variables, i.e., d + and d − , corresponding to tensile and compressive degradations respectively. Consider the definition of effective stress as follows: where ε, ε e and ε p are the total strain, elastic strain and plastic strain respectively. Note that here we employed an empirical model for the platicity evolution, which is expressed aṡ where b p denotes the empirical plastic flow parameter in which ξ p ≥ 0 denotes a plastic coefficient to control the plastic strain rate, and usually locates between 0.1 to 0.3; E 0 denotes the initial Young's modulus. Also note that the tensile plastic strain is neglected since it is relatively small. The effective stress represents the mesoscopic stress of the undamaged material. To represent the tensile and compressive loading conditions, Mazars et al. [Mazars and Pijaudier-Cabot (1989)] proposed the decomposition of effective stress tensor as follow: And according to the eigen based decomposition procedure proposed by Faria et al. [Faria, Oliver and Cervera (1998)], the positive and negative components of effective stress could be expressed as whereσ i and p i are the i-th eigen value and eigen vector of effective stress respectively; and the Macaulay brackets are defined as Further decompose HFE into elastic and plastic components, we have By noting the effective stress split in Eq. (28), the elastic HFE could be further decomposed into tensile and compressive components. Thus we have Consider the direct degradation in the energy level [Wu, Li and Faria (2006)], we obtain where the initial elastic HFE is For the plastic component of HFE, consider the degradation of plastic energy induced by tensile and compressive damages, we obtain where ψ p± 0 (κ) denotes the initial plastic HFEs and κ denotes the plastic variable. Substituting HFE decomposition (31) into Clausius-Duhem inequality (18) yields Thus we have the following three expressions Substituting Eq. (37a) into inequality (37b), we obtain the following two inequalities Further noting the Drucker's postulate and integrating over time domain yields In the present paper, we neglect the plastic evolution induced by tensile loading and damage. One obtains Inequality (37c) defines the damage energy release rates as the conjugated force of damages as follows By introducing the D-P plastic evolution potential to Eq. (41), Wu and coworkers [Wu, Li and Faria (2006)] solved the expression of damage energy release rates as follows whereĨ ± 1 are the first invariants ofσ ± ;J ± 2 are the second invariants of the deviatoric components ofσ ± ; E and ν denote the Young's module and Poisson's ratio; b 0 and α are material parameters which could be determined by multiaxial tests. Based on Eq. (42) and the damage consistent condition, Li et al. [Li and Ren (2009)] resolved the energy equivalent strain as follows: Then the damage evolutions undergoing multiaxial loading could be expressed as: where g ± (·) are the uniaxial tensile and compressive damage evolution functions which are often determined by uniaxial experimentations. In the present work, we consider the following damage evolution functions under tension and compression respectively [Feng and Li (2015); Feng, Kolay, Ricles et al. (2016)].
And the definitions of symbols are where E denotes the Young's modulus; f + and f − are the tensile and compressive strengths; ε + and ε − are strains corresponding to f + and f − on the stress-strain curves; and α + and α − are parameters governing the descending parts of tensile and compressive stress-strain curves.
To develop the multi-scale damage evolutions for bi-scalar damage model, we rearrange the HFE (31) as follows: As mentioned above, HFE ψ could be resolved by Eq. (15) based on the numerical simulations of the mesoscopic cell, and DREEs Y + and Y − could be calculated by Eq. (43). Here we propose the following two numerical cases corresponding to the standard mesoscopic cells to determine the tensile and compressive damage evolutions respectively.
• Tensile case As the numerical specimen undergoing tensile loading, we havẽ Then we obtain By rearranging Eq. (47), the tensile damage scalar could be expressed as • Compressive case The compressive loading condition gives Then the compressive damage could be resolved by Eq. (47) as follows Finally, the damage evolution functions g ± (·) could be resolved by Eqs. (50) and (52) respectively with the information of mesoscopic crack propagations.

Numerical example
As an illustrating example,a composite structure undergoing direct tension is studied (Fig. 2). The in-plane dimensions of the plate are l × b = 0.6 m × 0.3 m, and the thickness of the plate is 0.1 m. Totally 6 × 3 = 18 aggregates are uniformly distributed in the matrix. Each aggregate is with a diameter of 0.06 m. In the present example, we consider very strong aggregates that could remain linear elastic throughout the loading process, and the matrix experiences cracking and crushing at the same time. This example is specially designed to compare the proposed hierarchical multi-scale model with the fullscale simulation. Firstly, we perform the full-scale simulation of the composite structure as the benchmark results, although it is rather expensive in computation. The finite element mesh of the full-scale model is developed as Fig. 3. And the stress chains also experience more serious damages due to the stress concentration, as shown by Fig. 6(a). Fig. 5(b) suggests that the mechanism of internal force chains fails due to the strain localization induced by crack. And the major crack is clearly sketched by the concentration of damage in Fig. 6(b).  Figs. 9-10 also find the chains of internal forces and damages in the beginning and strong localization of damage in the end. By using the energy integration Eq. (15) and multi-scale damage characterization Eq. (50), the damage evolution could be resolved as Fig. 11(b). And the homogenized stress-strain curve could be also characterized based on the theory micro-mechanics, as shown in Fig. 11(a). As the present problem is governed by the tensile damage and fracture, we suppress the compressive damage evolution in the multi-scale damage analysis so that the compressive numerical test for the meso-cell is skipped. Note that the tensile damage at point A, where the model reaches its maximum tensile strength, is about 90%. This is just because here we simulated a plain quasi-brittle example. For this kind of material, once it reaches the tensile strength, which means the material cracks and will soon lose its capacity, so the damage of material will be large. This is not the case with some reinforced concrete material [Kakavand, Neuner, Schreter et al. (2018)] which will still has some capacity even after cracking. Finally, we could simulate the overall performances of the composite structure shown in Fig. 2 by using the coarse mesh (6 × 3) shown in Fig. 2. At each integration point of the coarse mesh, the stress-strain relations are represented by the damage model with damage evolutions characterized by Fig. 11(b). The load-displacement curves in Fig. 12 indicate good agreement between the full-scale damage model and the proposed multi-scale damage model. The peak load by full-scale model is 94.9 kN at displacement 0.0001431 m, while the peak load by multi-scale model is 96.7 kN at displacement 0.0001415 m. The differences between the full-scale and multi-scale model are only 1.9% for peak load and 1.1% for corresponding displacement.

Conclusions
For composite materials, the damage phenomenon is rather complex. In the subscale, the heterogeneous material structure often experiences cracks and material nonlinearities in rather complicated patterns. To develop the continuum damage model with considerations of meso-structure and cracks, the form of homogenised Helmholtz free energy is extended to consider both the effects of cracking and nonlinearity. Thereafter the damage evolution could resolved by the homogenised Helmholtz free energy and applied to the simulation of structures. Based on the proposed method, the simulation of structure of composite materials could be rather effective. In the meanwhile, the damage evolution laws informed by the cracking and nonlinearity in meso-level could enhance the capacity of the numerical simulation for the failure of composite material structures.