Shot Peening Processes to obtain Nanocrystalline surfaces in metal alloys: K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 319 Periodic homogenization and damage evolution in RVE composite material with inclusion Karim Benyahi, Youcef Bouafia, Mohand Said Kachi, Amel Hamri Department of civil engineering, Mouloud Mammeri University, 15000 Tizi Ouzou, Algeria. karim.benyahi@ummto.dz, https://orcid.org/0000-0003-2033-6217 youcef.bouafia2012@yahoo.com, kachi_ms@yahoo.fr, hamri.a@hotmail.com Sarah Benakli Laboratory LMSSC, National Conservatory of Arts and Crafts (CNAM Paris), France. benakli.s@gmail.com ABSTRACT. This work deals with the coupling between a periodic homogenization procedure, and a damage process occurring in a Representative Volume Element (RVE) of inclusion composite materials. We mainly seek, on the one hand to determine the effective mechanical properties according to the different volume fractions, and forms of inclusions for a composite with inclusions at the macroscopic level. On the other hand, to explore the rupture mechanisms that can take place at the microstructure level. To do this; the first step is to propose a periodic homogenization procedure, to predict the homogenized mechanical characteristics of an inclusion composite. This homogenization procedure is applied to the theory based on finite element analysis, by the Abaqus calculation code. The inclusions are modeled by a random object modeler, and the periodic homogenization method is implemented by python scripts. It is then a matter of introducing the damage into the problem of homogenization, that is to say; once the homogenized characteristics are assessed in the absence of the damage initiated, by microcracks and micro cavitations. It is then possible to introduce damage models, by a subroutine (Umat) in the Abaqus calculation code. The verifications carried out focused on RVE of composite materials with inclusions. KEYWORDS. Periodic homogenization; Elastic; Modeling; Inclusion; Damage. Citation: Benyahi, K., Bouafia, Y., Kachi, MS., Benakli, S., Hamri, A., Periodic homogenization and damage evolution in RVE composite material with inclusion, Frattura ed Integrità Strutturale, 58 (2021) 319-343. Received: 18.07.2021 Accepted: 11.09.2021 Published: 01.10.2021 Copyright: © 2021 This is an open access article under the terms of the CC-BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. INTRODUCTION he term "homogenization" is a process of calculating effective mechanical properties, and the term "localization" is used to determine stress and local strain (Suquet [1]). The first idea when using a homogenization procedure, it is necessary to apply the law of mixtures. Which postulates that the homogenized elastic modulus are the means of T K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 320 the elastic modulus of different constituents, weighted by their volume fractions. This method is rather easy to use but not very rigorous, in order, to make an approach allowing to appreciate the effective mechanical properties. Many analytical homogenization models, making it possible to predict the effective mechanical elasticity properties of a two-phase heterogeneous material. From those of its different constituents are proposed in the literature, those that have received the most attention are the diagrams diluted, Mori-Tanaka, self-consistent and Eshelby's solution [2-4]. However, these analytical or semi-analytical homogenization models, did not make it possible to go back to the local properties of the desired solution. As the complexity of the various microscopic phenomena, gave rise to more recent resolution methods. Nowadays, there are several micromechanical methods. The research having been carried out for the analysis, and prediction of the overall behavior of composite materials, and can be summarized as follows. Michel et al. [5] proposed two different families of numerical methods to solve the problem. The first method is based on the finite element method, by implementing periodicity conditions either by a control of strains or of stresses. The second numerical method is based on fast Fourier transforms (FFT). Which could be an alternative to the finite element method, for the micromechanical analyzes of representative solid elements with periodicity conditions. Sun et al. [6] proposed a procedure for predicting the elastic constant of the composite from the RVE, using the principles of strain energy equivalence. The average strain and stress for RVE are defined using Gauss's theorem and energy equivalence principles, and the full set of elastic constants for a unidirectional composite has been obtained. The appropriate stresses on the RVE under various loads were determined from the conditions of symmetry and periodicity. Xia et al. [7] proposed a method of micromechanical Finite Element Method (FEM) analysis applied to unidirectional, and right-angle laminates subjected to multiaxial loading conditions. On the basis of the general conditions of periodicity stated by Suquet [1]. They presented an explicit form of boundary conditions suitable for Finite Element Method (FEM) analyzes, of parallelepipedic RVE models subjected to multiaxial loads. Several theoretical models have been proposed for determining the effective properties of composite materials. And the research that has been conducted to improve our understanding, of periodic homogenization is summarized below. Yvonnet et al. [08] proposed a method for evaluating the higher-order tensors of an efficient general anisotropic strain- gradient (Mindlin) model, by a homogenization technique to determine the effective parameters. This proposed method uses finite element calculations on RVEs based on the principle of superposition. Jakabcin and Seppecher [09] studied the applicability of these formulas for highly contrasted structures. The study is carried out on structures of which the limiting energy is already known and compares the energies given by the convergence results, the corrective formulas and by a direct numerical simulation of the complete structure. Guinovart-Díaz et al. [10] proposed an analytical method of asymptotic homogenization (AHM) to the calculation of the effective elastic stiffnesses of a composite reinforced with fibers, with an imperfect contact between the matrix, and the fibers for parallelogram-like arrangement of fibers. Savvas et al. [11] proposed a simulation of the homogenization of random heterogeneous media with arbitrarily shaped inclusions. This study is carried out by modeling with the extended finite element method (XFEM), coupled with Monte Carlo simulation (MCS). And where, the influence of the inclusion shape on the effective properties of random media was studied. Tsalis et al. [12] proposed a method for the introduction of the admissible strain fields of materials, with generalized periodicity and present their properties. Tsalis, et al. [13] presented an analytical homogenization technique of elastic composites, with generalized periodicity to explain the effect of microstructure nonlinearity on effective properties. Wu et al. [14] used a method to reduce computational cost, and obtain more accurate approximations to predict the effective thermos-mechanical properties of random heterogeneous materials. By Richardson's extrapolation technique using smaller cell domains, and without the need for calculations in a larger cell domain. Godin [15] proposed an explicit formula to calculate the effective properties tensor of a periodic lattice of two-phase dielectric tubes embedded in a host matrix. Dhimole et al. [16] used a method based on the minimum energy loss of the structural genome, to predict the mechanical properties of three-dimensional (3D) four-directional braided composites. Beicha et al. [17] conducted a comparative study between the effective elastic properties of fiber-reinforced composites respectively built with a hexagonal, and random distribution of non-overlapping fibers. The results show that the observed constraints are higher for random distribution than for periodic distribution. Bonfoh et al. [18] proposed a general formulation of the multi- coated inclusion problem in the general case of ellipsoidal inclusions, and anisotropic elasticity. Rodríguez-Ramos et al. [19] used the Asymptotic Homogenization Method (AHM), and eigenfunction expansion-variational method (EEVM) to obtain the effective elastic moduli of two-phase fibrous periodic composites, for different types of parallelogram cells. Nevertheless, the latter are more and more used for structures used at their limit, which prompts us to look for a method to predict the onset of damage and the resulting behavior. Some analytical investigations directly related to this study are briefly reviewed. Yang and Misra [20] proposed a second gradient stress–strain theory for materials following damage elasticity based on the method of virtual power. This approach made it possible to develop the equations governing cohesive materials undergoing damage. Khoroshun and Shikula [21] described several mathematical models on the K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 321 different coupled deformation, and long-term microdamage processes applied to linear elastic composites of stochastic structure. Berthier et al. [22] developed a time-dependent nonlocal continuous damage model, which takes into account the transfer of elastic energy stored in fibers in breaking energy and viscous dissipation. Dorhmi et al. [23] developed micromechanical model to describe the progressive loss of stiffness observed during plastic straining based on microstructure changes. The study is being carried out on ductile metal-matrix composites subjected to mechanical loadings. Sorić et al. [24] proposed a numerical modelling of the responses to ductile damage in heterogeneous materials, by a second order homogenization approach. The evolution of ductile damage at the microlevel is taken into account, by using the gradient enhanced elastoplasticity. Wu et al. [25] proposed an improved gradient homogenization procedure for fiber reinforced materials. In this model, the fiber is considered as transversely homogeneous isotropic and assumed to remain linear elastic. While the material of the matrix can be considered as homogeneous isotropic. It is modeled as elastoplastic coupled to a damage law described by a non-local constitutive model. The method has been validated by the simulation of a damage process, in unidirectional carbon fiber reinforced epoxy composites, subjected to different load conditions. Devries et al. [26] developed the main results of the homogenization theory of periodic media, as well as the assessment and simulation of damage for composite materials. Their applications focused on the simulation of the evolution of damage by fiber failure in a unidirectional composite, involving parameters defined at the microscopic scale. And the prediction of takeoff near a free edge in a layered structure, using asymptotic extensions of the boundary layer. Xia et al. [27] developed a three-dimensional meso/micro-mechanical finite element multilayer model for the prediction of the overall mechanical behavior of glass fibers [0,903,0] T/epoxy laminate, and for the study of damage mechanisms in reinforced polymer laminates. The epoxy matrix is represented by a nonlinear viscoelastic constitutive model, and a criterion of damage to the epoxy matrix is introduced into the finite element model. The model prediction was in good agreement with the observation of experience, not only in the overall stress-strain response, but also in the initiation and propagation of damage to the transverse matrix. In this article, periodic boundary equations have been established and applied to the theory based on finite element analysis to predict macroscopic elasticity. While considering a microstructure composed by several inclusions with different shapes and orientations, and a random distribution of the latter. The use of the homogenization technique allows us to study various damage mechanisms, by introducing damaged medium models from a microscopic approach. While still utilizing the efficient mechanical properties of inclusion composites. This allowing us to take into account the damaging behavior that occurs in the material, specifically the matrix at the micro scale. This document is organized as follows. Section 2 discusses a method of solving a periodic homogenization problem using the principles of strain energy equivalence, in a coordinated way with finite element analysis. It gives the modifications required in the Abaqus calculation code by python srcipts, in order to introduce periodic boundary conditions at the micro scale. This procedure is performed by the introduction of additional degrees of freedom supporting the components of macroscopic strains. In section 3 are presented several damage models controlled by an equivalent strain tensor in the principal coordinate system. The inclusions are supposed to remain elastic, while the matrix is supposed to obey the actual behavior of the material constituting, it with damage at the micro level. Section 4 presents detailed results obtained by several tests of this selected periodic homogenization method, which will be compared with an estimated homogenization method (Mori-Tanaka model). These results can be used to verify the implementation of the simplified approach proposed in the Abaqus calculation code. Finally, Section 5 provides applications regarding the introduction of damage into the homogenization problem. It is a question of using a local formulation for the evolution of the damage at the level of the microstructure, using the subroutine (Umat) implemented in the Fortran language on the Abaqus calculation code. This proposed method makes it possible to simulate the response of the matrix behavior, on several RVEs requested in uniaxial traction. PERIODIC BOUNDARY CONDITION onsider a periodic structure made up of a periodic network of repeated unit cells. In the case of a composite material with periodic microstructure, we can define a "basic cell" and also periodicity vectors (Fig. 1). The field of displacement for a periodic structure, can be expressed as indicated by Suquet [1]: = + ' ' ( ) . ( ) , iji i j i i iU x x u x u periodic (1) with: C K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 322  . ij jx : Represents a linear displacement field, ' ( )i iu x : is a periodic function representing a modification of the linear displacement field, because to the heterogeneity for the composite structure. Figure 1: Description of a different basic cells in periodic microstructure. For each unit base cell, its boundary surfaces should appear in parallel pairs, and displacements on a pair of opposite parallel boundary surfaces can be expressed as:  + + = + ( ) ( ) ' ij k k i j iu x u (2)  − − = + ( ) ( ) ' ij k k i j iu x u (3) + k and − k : The indices to identify the th k pair of two opposite parallel boundary surfaces of a representative volume element (RVE). The periodic function ' ( )i iu x is considered to be the same for two parallel borders, the difference between the Eqn. 2 and Eqn. 3 is given as follows:   + − + − − = − =  ( ) ( ) ( ) ( ) ( ) k k k k k i i j ijij j ju u x x x (4) with:  k jx : The constants for each pair of parallel boundary surfaces, with  ij known,  ij : The tensor of the macroscopic strains of the periodic structure. The average tensors of stress and strain are obtained by an integration on the RVE, as follows: ( ) ij ijV 1 = x, y, z dV V (5) ( ) ij ijV 1 = x, y, z dV V (6) The strain energy of a homogeneous composite material over the RVE, is as follows:  = ij ij 1 2 U V (7) The stored strain energy of a heterogeneous composite material over the RVE is as follows: K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 323  =  * ij ij V 1 dV 2 U (8) ( )   += − * ij ijij ij V 1 dV 2 U (9) ( )    −= +  * ij ijij ij ij V V 1 1 dV dV 2 2 U (10) with:  ij : The strain tensor. From Eqn. 7 and Eqn. 10, we can therefore write the energy difference, as follows: ( ) ( )( )    − =      − − iij iij* ii V 1 + dV 2 j j U u u uU u x x (11) with: iu : The ith average displacement. and using the equilibrium equation:    = ij 0 jx (12) By replacing Eqn. 12 with Eqn. 11, we obtain: ( )( ) − =   −  iij i * V 1 dV 2 j u x u U U (13) Using Gauss's theorem, we can transform the integration on the volume into the integration on the surface, as follows: ( ) −− =  * iij i j S 1 d 2 j juU SuU (14) Then, at the level of the surface jS , we have: = iiu u (15) We can conclude that: − = * 0U U (16) The above derivation shows that the homogenization, we have used ensures the equivalence between the heterogeneous and homogenized RVE. To apply the periodic boundary conditions at the micro scale, there are three types of set of nodes: faces, edges and corners in the RVE having been modeled in the form of a parallelepiped in 3D. The regions on the boundaries of the RVE are selected and numbered as shown in Fig. 2. Additional equations should be introduced for periodic boundary conditions at the micro scale, where the macroscopic strain matrix has components  = ij with i, j = x, y, z. In order to include the macroscopic deformation of component K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 324  ij , it is necessary to create six reference points (PR-1, PR-2, PR-3, PR-4, PR-5, PR-6). Where their degree of freedom of displacement has a relation with      x xy xz y yz z, , , , , respectively as follows:      − − − − − −= = = = = = 1 1 1 1 1 1 x xy xz y yz z( 1) ( 2 ) ( 3) ( 4 ) ( 5 ) ( 6 ); ; ; ; ;PR PR PR PR PR PRU U U U U U (17) Figure 2: Representation and periodic boundary conditions of an RVE. We create the additional equations of all faces, edges and corners with the python command in order to apply the periodic boundary conditions in the Abaqus calculation code. And we apply these periodic boundary conditions at the micro scale. Finally, to solve the homogenization equations, we use the finite element method. Equations for applying periodic boundary conditions at the corners:                    x yz xzG A x yz yxC E x yz xzD F x yz yxH B u -u = + + u -u = + - u -u = - - u -u = - + (18) Equations for applying periodic boundary conditions at the faces:      xBCGF ADHE FEHG BADC CDHG BAEF u -u = u -u =0 u -u =0 (19) Equations for applying periodic boundary conditions at the edges:                        x xzFG AD x xzHE CB yx yzCG AE yx yzDH BF yx yzEF DC yxHG AB yz u -u = + u -u = - u -u = + u -u = - u -u = + u -u = - (20) K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 325 The mean strain in Eqn. 6 can be transformed as follows:  = −ij ( ) 1 ( ) j i j j i j S u n u n dS V (21) In order to evaluate the effective properties ijC , it is first necessary to evaluate the mean of the stress and the strain of Eqns. 5 and 6 respectively, by the application of different outline conditions. Then to introduce them into the constitutive relation as follows:   = ij ij ij C (22) with: ijC : The stiffness modulus corresponding to the application of the strain mode. From the relation between the flexibility matrix S, and the elastic property, we can calculate the effective elastic properties from the stiffness matrix of the RVE.                                            3121 3311 22 12 32 11 22 33 13 23 33 11 22 12 13 23 1 - 0- 0 0 EE E 1 - - 0 0 0 E E E 1 0 0 0 - - E E ES = 1 00 0 0 0 G 1 0 0 0 0 0 G 0 0 0 0 1 0 G (23) There are only five independent material constants ( ) 11 33 12 31 32, , , ,E E G for a transversely orthotropic / isotropic material. Where 13 depends on ( ) 11 33 12 31, , ,E E and can be calculated by the following relation:  31 13 33 11 = E E (24) A material is isotropic when its property is the same in all three directions. In this case, ( )= = =11 22 33E E E E , ( )   = = =12 31 32 and ( )= = =12 31 32G G G G . There are only two independent material constants ( ),E for an isotropic material. However, anisotropic materials have non-zero values in the upper right and lower left parts of their flexibility and stiffness matrices. There exists a relation between the components of the flexibility matrix, with the components of the rigidity matrix, which is as continuation: ( ) ( ) 2 11 33 13 11 2 2 2 2 1133 11 11 13 12 33 12 13 C .C -C 1 S = = EC .C -2.C .C -C .C +2.C .C (25) K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 326 ( ) ( )  2 12 33 13 21 12 2 2 2 2 2233 11 11 13 12 33 12 13 C .C -C S =- =- EC .C -2.C .C -C .C +2.C .C (26) ( ) 13 31 13 2 3311 33 13 12 33 C S =- =- EC .C -2.C +C .C (27) ( ) ( ) 11 12 33 2 3311 33 13 12 33 C +C 1 S = = EC .C -2.C +C .C (28) 44 44 12 1 1 S = = C G (29) 55 55 13 1 1 S = = C G (30) The determination of the effective elastic properties depends on the stress-state, and also on the property of the composite material (isotropic or orthotropic case) of the RVE. For an isotropic case and a uniaxial stress-state according to the direction X-X, the effective elastic properties are determined:   eff eff 11 11 11 E =E = (31)    eff eff eff 12 13+= 2 (32) with:             eff 22 12 11 eff 33 13 11 =- =- (33) For an orthotropic/isotropic case transversely and a uniaxial stress-state along the direction X-X, the effective elastic properties are determined:   eff 11 11 11 E = (34)    eff 22 12 11 =- (35) GENERAL IMPLEMENTATION FLOWCHART ON THE ABAQUS CALCULATION CODE he proposed periodic homogenization method has been implemented numerically using a Python script, to invoke the solution in the Abaqus finite element computational code. The iterative algorithm is presented in Fig. 3. At each incremental step of the analysis, the material properties subroutine (Umat), where the constitutive law was used to carry out the communication between the Python code, and the Abaqus solver. The Python code modifies the boundary conditions to an Abaqus input file, applied to the RVE. T K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 327 Figure 3: Periodic homogenization procedure implementation flowchart. DAMAGE PROBLEM nce the homogenized characteristics are appreciated in the absence of the damage initiated by microcracks and microcavitations. It is then possible to introduce models of damage by a subroutine (Umat) in the Abaqus calculation code. The materials are supposed to obey the criterion of Von Mises [28], which is based on the second invariant of the deviatoric stress tensor. It is then supposed that the plasticity of the composite, being governed by that of the matrix then the total plasticity criterion will be obtained from that of the matrix. To translate the real behavior of materials into nonlinear elasticity, relatively developed laws are introduced within the framework of damage mechanics. Mazars damage model The Mazars model ([29], [30]) was developed within the framework of damage mechanics (Fig. 4). The stress is given by the following relation: ( )  e= 1-D E (36) The damage is controlled by the equivalent strain eq , which allows to translate a tri-axial state by an equivalence to a uni- axial state. The strain tensor in the principal coordinate system:     22 2 eq 1 2 3+ + + = + + (37) Knowing that the positive part + is defined, such that if  i is the main strain in direction i:  i + =       i i i if 0 0 if 0 (38) eq is an indicator of the state of tension in the material, which generates the damage. This quantity defines the load surface f, such that: O Import the RVE generated in the format (.step) into the Abaqus calculation code Generate the microstructure (geometric generator of random objects) Periodic homogenization by python scripts (imposition of periodic boundary conditions, introduction of elimination equations): - Read the RVE file (.odb) - Calculate 𝜎 from the mean of the volume - Calculate the current Jacobian 𝐶 Start the simulation of the model under the Abaqus calculation code Determination of effective mechanical characteristics Introduce the constitutive laws of the constituents of composite by the subroutine (UMAT) on the Abaqus calculation code: - Read the macroscopic constraints - Read the current Jacobian matrix - Submit to Abaqus K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 328 ( )eqf= -K D =0 (39) with: ( ) eqK D = D=0if (40) Figure 4: Behavior of concrete according to the Mazars model ([29], [30]). The relation between these variables is given as follows:    t t c cD= D + D (41) Usually the value of  is fixed at 1.06. The coefficients t and c carry out a link between the damage, and the state of traction or compression. When traction is activated  =t 1 while  =c 0 , and vice-versa in compression. The damage evolution laws tD and cD are expressed only from the equivalent strain eq . ( ) ( )( )   t t t t eq d0 eq 1-A D =1- -A exp -B - (42) ( ) ( )( )   c c c c eq d0 eq 1-A D =1- -A exp -B - (43) with: d0 : The threshold strain of damage,  : Coefficient which was introduced later to improve the shear behavior, t c t cA A B, , , Band : Material parameters to identify. Damage model from Bouafia et al. To describe the behavior of composites with cylindrical inclusions in tension. The relationships proposed by Bouafia et al. ([31], [32]) within the framework of the theory of beams, with taking into account of the damage which was developed for K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 329 the behavior of the concrete of fibers in tension (Fig. 5). The fibers (cylindrical inclusions) are dispersed in the concrete at random, and the modeling is carried out considering a uniform distribution. Figure 5: Constitutive law (σ-ε) in traction of steel fiber concrete [31]. The variable of the damage in uni-axial traction, is given by: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )                       −     → = −             → = − −     6 u ft u 1 6 t t ft u 6 u u r 2 6 t r u - 1 - - - 1 . 1 - uc ftuc t ct ct uc t ct f D E E D E (44) The maximum ultimate stress of the composite (function of the characteristics of inclusions):      0 =c f u ul (45) The reference length is linked to the height h of the section: = hrl (46) The initial modulus of the composite in traction is given by: ( ) +0 01=ct bE nE (47) The ultimate strain corresponding to the total mobilization of the inclusions-matrix adhesion, is given by:      + 2 = 3 u f u ft f l E h (48) When there is tearing of inclusions, the breaking strain of the composite, is given by:       + + 2 = 3 4 u f f rt ft f l l E h h (49) K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 330 This strain is a consequence of opening the cracks too large. This fracture strain of the inclusions is limited as follows:  =rt ft (50) The damage is progressive in the field of strains. For a uni-axial state of tension, the equivalent strain ( )( )  1 2,t tD D , in the nonlinear field is such that: ( ) ( ) ( ) ( ) ( ) ( )                → = =        → = =  ft u 1 1 1 u r 2 2 2 = 1- 1- = 1- 1- t ct t t t ct t t D E E D D D E E D D (51) and ( ) ( ) ( ) ( ) ( )         =       → =  1 1 ft ft 1 r 2 2 r f , - 0 f , - 0 t t t t t D D D D D (52) Figure 6: Damage model algorithm implementation diagram. The threshold function defines two limit states: - An elastic state ( )ft . - A breaking limit given by the value ( ) r . Such as: ( )      =  = →   = → → ft 1 1 r 2 =0 f , 0 =0 t t t D D D (53) K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 331 IMPLEMENTATION OF DAMAGE MODELS IN ABAQUS CALCULATION CODE he model of Mazars ([29], [30]), the model of Bouafia et al. ([31], [32]) simulating an isotropic elastic-damageable behavior in tension, they are implemented in the Abaqus calculation code. Their algorithms thus obtain are presented in Figs. 6a-6b respectively. DAMAGE SEARCH AND INTRODUCTION FLOWCHART he flowchart for finding, and introducing damage into a periodic homogenization problem, is shown below (Fig. 7): Figure 7: Damage search and introduction flowchart. NUMERICAL COMPARISON OF HOMOGENIZATION MODELS TO EXPERIENCE Lightweight concrete composite n this example, the concrete used is hydraulic lightweight aggregates, the matrix is cementitious while the aggregates are expanded clay [33]. The elastic properties of the constituents of composite are presented in Tab. 1. In this example, it is a question of comparing our numerical results obtained (by the semi-analytical method of homogenization (Mori-Tanaka Model), and by the method of periodic homogenization (Numerical Model)), with the experimental results (Experimental) [33]. These latter are obtained from compression tests carried out on cylindrical specimens (16 x 32 mm), for different volume fractions of light aggregates: 0.125, 0.25, 0.375 and 0.45. The width of the RVE is taken on the order of the diameter of median grain of sand (0.08 mm) [33]. T T I Yes Breaking No Yes ε̃ ≥ ε̃Crit No 𝐢𝐧 𝐜 𝐫 = 𝐢𝐧 𝐜 𝐫 + 𝟏 Periodic homogenization 𝐶0 𝑛, 𝐶 𝑛, 𝜎 𝑛, 𝐸𝑛 Effective properties E, ν, G, K, incrmax ∆𝜀 Imposed Mazars model parameters At , Ac , Bt , Bc , ε̃ Crit Bouafia et al. model parameters ϖ, θ0, lr , β, h, lf, Ef, τu, ϕ, εrt , ε̃Crit incr ≥ incrmax End of simulation K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 332 Lightweight concrete Young's modulus (MPa) Poisson coefficient Matrix 23630 0.20 Aggregates (expanded clay) 5679 0.15 Table 1: The elastic properties of the constituents of lightweight concrete composite [33]. Figure 8: Variation of the equivalent Young's modulus as a function of the percentage of aggregates. We notice in Fig. 8 that the equivalent Young's modulus values found by the semi-analytical homogenization method (Mori-Tanaka Model) are close to the experimental values (Experimental) [33], and that for very low volume fractions (less than 15%). On the other hand, the periodic homogenization model (Numerical Model) gives good results, which approach the experiment allowing a better prediction of the elastic characteristics, and that for different volume fractions. Cementitious composite with grains of sand In this example, the heterogeneous material consists of a matrix and the spheroidal inclusions shown in Fig. 9. The elastic properties of the constituents of composite are presented in Tab. 2. Lightweight concrete Young's modulus (MPa) Poisson coefficient Matrix (cement paste) 7000 0.25 Inclusion (sand) 107000 0.15 Table 2: The elastic properties of the constituents of sand mortar composite [34]. The elastic modulus of sand grains [34] was taken to be equal to that of the siliceous aggregates, ie 107000 MPa, the Poisson's ratio was taken to be 0.15. Measurement of the elastic modulus of cement paste gave 6000 MPa. To improve the calibration of the curve between modulus and volume of sand [34], it was taken equal to 7000 MPa. The Poisson's ratio was taken equal to 0.25. The width of the RVE is taken on the order of the diameter of median grain of sand (0.02 mm) [34]. The RVE used in this example is generated by a random object modeler, and this with a random distribution of spheroidal inclusions (with a number of 10 inclusions in the RVE), and with different volume fractions. The inclusions of the RVE are generated randomly, based on a sphere pattern. Once the random generation of the inclusions of RVE is carried out, we will integrate it into the Abaqus calculation code. The size of the RVE does not change, only the volume fraction K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 333 which evolves. The objective of this study is to determine the effective behavior of a composite (matrix - inclusion of spheroidal shape). A comparison of the mechanical characteristics between the semi-analytical homogenization model (Mori-Tanaka Model), the periodic homogenization model (Numerical Model) and the experimental measurements [34] is shown in Fig. 10. Figure 9: Modeling of the material with spheroidal inclusions on the Abaqus calculation code. Figure 10: Elastic modulus of sand mortars. This involves comparing the results obtained by the semi-analytical homogenization model (Mori-Tanaka Model), and by the periodic homogenization method (Numerical Model) with the experimental results (Experimental) [34]. The elastic properties of isotropic materials (of the cement matrix and those of the sand) are average values, and they are obtained for different volume fractions of spheroidal shaped aggregates. We note that the estimations of the semi-analytical homogenization model (Mori-Tanaka Model), the results obtained by the periodic homogenization method (Numerical Model), and the results from the experiment (Experimental) are quite similar, and this for different volume fractions. ESTIMATION OF THE MECHANICAL CHARACTERISTICS AND TAKING INTO ACCOUNT THE DAMAGE OF AN INCLUSION COMPOSITE n this paragraph, we will deal with several examples in order to validate the periodic homogenization model, and this by varying the percentage and shape of inclusions for RVE of composite materials. In this study, we will compare the results of our modeling by the periodic homogenization method which are represented by (Numerical I K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 334 Homogenization), with the results of a semi-analytical method (Mori-Tanaka model) from the work of Al Kassem [35], which are represented by (Mori-Tanaka Model). Once the effective mechanical characteristics are assessed in the absence of damage, then we introduce damage models for the RVE in composite materials, by a subroutine (Umat) on the Abaqus calculation code. Determination of effective mechanical characteristics In this example the heterogeneous material consists of a matrix and inclusions. The elastic properties of the constituents are presented in Tab. 3. Materials Young's modulus E (MPa) Poisson coefficient ν Matrix 2800 0.35 Inclusion 72000 0.172 Table 3: The elastic properties of the components of composite (matrix-inclusions) [35]. The Representative Volume Element (RVE) of the inclusion composite consists of two phases: the matrix and the inclusion which follow an isotropic behavior, this is a three-dimensional (3D) case in small perturbations. The matrix is modeled as an elementary cube and the inclusion is modeled in several forms (spheroidal, ellipsoidal and cylindrical), and as rigid bodies. The homogeneous boundary conditions introduced in our modeling do not give other stresses than in the direction of perturbation; therefore, they result in a uni-axial stress state under tensile loading. The inclusions of RVE are generated at random (with a number of 10 inclusions in a RVE), based on a pattern of inclusions in the form of a sphere, ellipse and cylinder. Once the random generation of the inclusions of RVE is carried out, we will integrate it into the Abaqus calculation code. The size of RVE does not change (size_rve_x= size_rve_y= size_rve_z=1) [35], only the volume fraction which evolves. The objective of this study is to determine the effective mechanical characteristics of a composite (matrix - inclusions). Figure 11: Modeling of an inclusion composite. The results obtained in Fig. 11, show that all the principal stresses of the microstructure are accumulated in the x direction, or at the level of the perturbation direction using periodic boundary conditions. A comparison of the effective mechanical characteristics between the Mori-Tanaka estimation method (Mori-Tanaka Model) [35], and the periodic homogenization model developed in this study (Numerical Homogenization) is shown below (Figs. 12-14). And this for different volume fractions of inclusions (spheroidal, ellipsoidal, cylindrical). This involves comparing the numerical results obtained by the periodic homogenization method (finite element model with periodic boundary conditions), with the results obtained by the Mori-Tanaka model [35]. The elastic properties of the materials (of the matrix and those of the inclusion) are average values, and they are obtained for different volume fractions of inclusions (spheroidal, ellipsoidal, cylindrical). We note that at low volume fractions, we observe a very small difference between our model (Numerical Homogenization), and the Mori-Tanaka model (Mori-Tanaka Model) [35], while for large volume fractions, we observe a bigger difference. K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 335 Figure 12: Homogenized mechanical characteristics for different volume fractions of spheroidal inclusions. Figure 13: Homogenized mechanical characteristics for different volume fractions of ellipsoidal inclusions. K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 336 Figure 14: Homogenized mechanical characteristics for different volume fraction of cylindrical inclusions. It is also noted that the estimations of the semi-analytical model of Mori-Tanaka [35], and the results of the periodic homogenization method developed in this study are quite close. And that for very low volume fractions (less than 15%) caused by a weak interaction between the phases. On the other hand, for volume fractions greater than 15%, the effective mechanical properties are overestimated, a deviation is observed which can be caused by a phenomenon of interaction between the inclusions. We conclude, that the periodic homogenization method (PBC) gives mechanical characteristics more rigid, than those of the semi-analytical model (Mori-Tanaka model). Introduction of damage into the problem of homogenization We implemented a Mazars damage model ([29], [30]) for spheroidal and ellipsoidal shaped inclusions, and the model of Bouafia et al. ([31], [32]) for cylindrical shaped inclusions in the Abaqus calculation code by the use of a subroutine (Umat). Calculations on a RVE containing 10 inclusions with type localization conditions (PBC), allowed us to study the K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 337 influence of the variation of the volume fraction of inclusions on the mechanical characteristics. Then, once the effective mechanical properties of the composite were known, a damage model was introduced to study the response of the microstructure. The effective mechanical properties of the representative volume element ( ),eff effE are those obtained from the periodic homogenization method, corresponding to each volume fraction. The volume element is stressed in pure tension. The loading is done with displacement imposed in three-dimensional strain. The characteristics of composite are defined in Tab. 4: Table 4: The effective properties and damage parameter of a composite. In a first case, we represent in Figs. 15-16 the local stress-strain response and the evolution of the damage respectively, depending on the strain of a volume element loaded in pure tension, and under various volume fractions. We set the parameter tA and made the parameter tB varied. Figure 15: Influence of the variation of the volume fraction inclusions on the local stress-strain behavior. Matrix- Spheroidal inclusions Volume Fraction (%) Eh (MPa) νh AC BC At Bt εD0 10 3440.43 0.333 1.4 240.83 0.8 2752.35 0.001 15 3764.60 0.335 1.4 263.52 0.8 3011.68 0.001 20 4264.90 0.326 1.4 298.54 0.8 3411.92 0.001 25 4779.03 0.316 1.4 334.53 0.8 3823.23 0.001 30 5324.70 0.312 1.4 372.73 0.8 4259.76 0.001 Matrix- Ellipsoidal inclusions Volume Fraction (%) Eh (MPa) νh AC BC At Bt εD0 10 3391.73 0.340 1.4 237.42 0.8 2713.39 0.001 15 3787.73 0.333 1.4 265.14 0.8 3030.19 0.001 20 4240.37 0.328 1.4 296.82 0.8 3392.30 0.001 25 4784.43 0.321 1.4 334.91 0.8 3827.55 0.001 30 5320.80 0.312 1.4 372.46 0.8 4256.64 0.001 Matrix- Cylindrical inclusions 𝐅𝐜𝐣 (MPa) Ftj (MPa) Eh (MPa) νh ε0 (%) ϕ (mm) lf (mm) τu (MPa) ω (%) θ0 47.76 3.46 5247.08 0.342 0.21 0.04 0.8 7 5 0.405 47.76 3.46 7780.74 0.335 0.21 0.04 0.8 7 10 0.405 47.76 3.46 10360.73 0.327 0.21 0.04 0.8 7 15 0.405 47.76 3.46 13033.63 0.320 0.21 0.04 0.8 7 20 0.405 K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 338 Figure 16: Influence of the variation of the volume fraction inclusions on the evolution of damage according to the strain. Figs. 15-16 show us that the increase in the volume fraction of the inclusions (spheroidal, ellipsoidal), and of the parameter tB allows an improvement in the stiffness of composite. And it makes a change the shape of damage as a function of strain, or the composite becomes less ductile. In a second case, the influence of the parameter tA is given in Figs. 17-18 with respect to the local stress-strain response, and the evolution of damage respectively, of a volume element loaded in pure tension. With an average volume fraction (taken equal to 20%), and a fixed value of the parameter tB corresponding to this percentage. Figs. 17-18 show us that the increase in the volume fraction of shaped inclusions (spheroidal, ellipsoidal), and of the parameter tA results in a sudden drop in the strength of composite (very low residual stress). And it makes a change the shape of damage as a function of the strain, where the composite becomes more fragile. Figure 17: Influence of the parameter tA on the local behavior stress-strain. We represent the local stress-strain response and the evolution of damage respectively, as a function of the strain of a volume element stressed in pure tension, and under different volume fraction of cylindrical inclusions. Figs. 19-20 show us that increasing the percentages of cylindrical inclusions increases elastic properties (the material becomes more rigid). And it causes a sudden drop in the strength of composite (low residual stress), and a ductility bearing having greater strain of the composite. What changes the rate of damage depending on the strain, where the resistance of the softening part of composite is very low. K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 339 Figure 18: Influence of the parameter tA on the evolution of damage according to the strain. Figure 19: Local stress-strain behavior. Figure 20: Evolution of the damage according to the strain. The behavior of the composite follows an elastic phase until the appearance of the first cracks, then a sudden drop in stress results to a residual stress which activates the contribution of cylindrical inclusions to take up the forces developed. K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 340 Finally, this results in a level of ductility allowing greater strain of the composite. Note that the stress corresponding to this level is a function of the percentage of inclusions. Also, we notice, the more we increase the percentage of fibers, the shape of the local stress-strain behavior curve decreases, which corresponds to a less ductile behavior. CONCLUSION he effective properties of heterogeneous materials (Young's elastic modulus (E), Poisson ratio (ν), shear modulus (G) and compressibility modulus (K)) were evaluated for different volume fractions (5 at 30%), taking into account the periodic displacement boundary conditions. For very low volume fractions (<15%), the comparison of the various results obtained by the numerical model of periodic homogenization (PBC), with the results of estimations of the Mori-Tanaka model [35] is satisfactory, caused by a weak interaction between phases. Based analysis of the results, the following conclusions were drawn: - On the other hand, for volume fractions greater than 15%, the effective mechanical properties obtained by our simulation are overestimated, they are more rigid than those of the semi-analytical model of Mori-Tanaka [35]. - It is concluded that the percentage of inclusions as well as their shape, and orientation have an influence on the mechanical characteristics of the composite. The Mazars damage model ([29], [30]) allowed us to introduce a non-local formulation for the evolution of damage at the microstructure level. And also makes it possible to reproduce the softening behavior of the material, which is gradually reduced as the microcracks develop until reaching zero. Which corresponds to a visible macro-crack. Regarding the parametric study, the following ascertainment were drawn: - We set the parameter tA and we varied the parameter tB under different volume fraction, shows us that the variation of the volume fraction of shaped inclusions (spheroidal, ellipsoidal), and of the parameter tB increases the elastic properties of the composite (improvement of the stiffness). And, it slightly changes the shape of damage curve (the composite becomes a little less ductile). - On the other hand for an average volume fraction (taken equal to 20%), and a fixed value of the parameter tB corresponding to this percentage. Shows us that the variation of the parameter tA makes it possible to have a sudden drop in the resistance of the composite (no residual stress), as well as a very low contribution of inclusions in the softening part of the composite. And, it changes the shape of the composite damage curve (the composite becomes more fragile). The damage model of Bouafia et al. ([31], [32]), allowed us to introduce the evolution of damage through a ductility plateau, which is a function of the characteristics of cylindrical inclusions (percentage, diameters, orientation and bond stress). ACKNOWLEDGEMENTS he authors wish to thank the Algerian Ministry of higher education and scientific research for funding the University education research project (PRFU – N° A01L02UN150120180004) and Tassili Project (PHC – 43940NJ). DECLARATION OF INTEREST STATEMENT n behalf of all authors, the corresponding author states that there is no conflict of interest. NOTATION Sj : The jth area, T T O K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 341 nj : The unit normal, ui : The ith displacement, Sj : The outer limit of the RVE, E11, E22 and E33 : The moduli of elasticity of extensions along directions 1, 2 and 3, respectively, νij (i, j = 1, 2, 3) : Poisson's ratios, G12, G13, and G23 : Shear moduli, Ui P : The nodal variable at the node P, of degree of freedom i, E : Hooke's matrix, D : The damage variable, εe : Elastic strain, Eb0 : Modulus of longitudinal elasticity of concrete, ϖ : Percentage by volume of fibers, θ0 : Fiber orientation coefficient, lr : Reference length, β: Model constant, h : Composite cross section height, lf : Fiber length, Ef : Elastic modulus of the fiber, n : Fiber-concrete equivalence coefficient, τu : Ultimate fiber-concrete bond stress, ϕ : Diameter of a fiber, fft : Composite tensile strength, εft : Composite cracking strain, εrf : Fiber breaking strain, εrt : Breaking strain of the composite in tension. REFERENCES [1] Suquet, P. (1987). Elements of homogenization theory for inelastic solid mechanics. Sanchez-Palencia, E., Zaoui, A. (Eds.), Homogenization Techniques for Composite Media. Springer-Verlag, Berlin, pp. 194–275. DOI: 10.1007/3-540-17616-0 [2] Mori, T., Tanaka, K. (1973). Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta metallurgica, 21(5), pp. 571-574. DOI: 10.1016/0001-6160(73)90064-3. [3] Hill, R. (1965). A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids 13, pp. 213-222. DOI: 10.1016/0022-5096(65)90010-4. [4] Eshelby, J.D. (1957). The determination of the elastic field of an ellipsoidal inclusion and related problems. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 241(1226), pp. 376-396. DOI: 10.1098/rspa.1957.0133. [5] Michel, JC., Moulinec, H., Suquet, P. Effective properties of composite materials with periodic microstructure: a computational approach. Computer methods in applied mechanics and engineering 172 (1-4), pp. 109-143. DOI: 10.1016/s0045-7825(98)00227-8. [6] Sun, C T., Vaidya, R S. (1996). Prediction of composite properties from a representative volume element. Composites Science and Technology, 56(2), pp. 171-179. DOI: 10.1016/0266-3538(95)00141-7. [7] Xia, Z., Zhang, Y., Ellyin F. (2003). A unified periodical boundary conditions for representative volume elements of composites and applications. International Journal of Solids and Structures 40, pp. 1907–1921. DOI: 10.1016/s0020-7683(03)00024-6. [8] Yvonnet, J., Auffray, N., and Monchiet, V. (2020). Computational second-order homogenization of materials with effective anisotropic strain-gradient behavior. International Journal of Solids and Structures, 191, pp. 434-448. DOI: 10.1016/j.ijsolstr.2020.01.006. https://doi.org/10.1007/3-540-17616-0 https://doi.org/10.1016/0001-6160(73)90064-3 https://doi.org/10.1016/0022-5096(65)90010-4 https://doi.org/10.1098/rspa.1957.0133 https://doi.org/10.1016/s0045-7825(98)00227-8 https://doi.org/10.1016/0266-3538(95)00141-7 https://doi.org/10.1016/0266-3538(95)00141-7 K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 342 [9] Jakabcin, L. and Seppecher, P. (2020). On periodic homogenization of highly contrasted elastic structures. Journal of the Mechanics and Physics of Solids, 144:104104. DOI: 10.1016/j.jmps.2020.104104. [10] Guinovart-Díaz, R., Rodríguez-Ramos, R., Bravo-Castillero, J., López-Realpozo, J., Sabina, F., Sevostianov, I. (2013). Effective elastic properties of a periodic fiber reinforced composite with parallelogram-like arrangement of fibers and imperfect contact between matrix and fibers. International Journal of Solids and Structures, 50(13), pp. 2022-2032. DOI: 10.1016/j.ijsolstr.2013.02.019. [11] Savvas, D., Stefanou, G., Papadrakakis, M. and Deodatis, G. (2014). Homogenization of random heterogeneous media with inclusions of arbitrary shape modeled by XFEM. Computational Mechanics 54(5), pp. 1221–1235. DOI: 10.1007/s00466-014-1053-x. [12] Tsalis, D., Baxevanis, T., Chatzigeorgiou, G., Charalambakis, N. (2013). Homogenization of elastoplastic composites with generalized periodicity in the microstructure. International Journal of Plasticity, 51, pp. 161–187. DOI: 10.1016/j.ijplas.2013.05.006. [13] Tsalis, D., Chatzigeorgiou, G., Charalambakis, N. (2012). Homogenization of structures with generalized periodicity. Composites: Part B, 43(6), pp. 2495–2512. DOI: 10.1016/j.compositesb.2012.01.054. [14] Wu, Y., Nie Y., and Yang Z., (2014). Prediction of effective properties for random heterogeneous materials with extrapolation. Archive of Applied Mechanics 84(2), pp. 247–261. DOI: 10.1007/s00419-013-0797-7. [15] Godin, Y.A. (2016). Effective properties of periodic tubular structures. Q. J. Mech. Appl. Math., 69(2), pp. 181–193. DOI: 10.1093/qjmam/hbw003. [16] Dhimole, V.K., Chen, Y., Cho, C. (2020). Modeling and Two-Step Homogenization of Aperiodic Heterogenous 3D Four-Directional Braided Composites. J. Compos. Sci., 4(4), 179. DOI: 10.3390/jcs4040179. [17] Beicha, D., Kanit, T., Brunet, Y., Imad, A., El Moumen, A., Khelfaoui, Y. (2016). Effective transverse elastic properties of unidirectional fiber reinforced composites. Mech. Mater. 102, pp. 47–53. DOI: 10.1016/j.mechmat.2016.08.010. [18] Bonfoh, N., Coulibaly, M., Sabar, H. (2014). Effective properties of elastic composite materials with multi-coated reinforcements: A new micromechanical modelling and applications. Compos. Struct. 115, pp. 111–119. DOI: 10.1016/j.compstruct.2014.04.011. [19] Rodríguez-Ramos, R., Yan, P., López-Realpozo, J.C., Guinovart-Díaz, R., Bravo-Castillero, J., Sabina, F.J., Jiang, C.P. (2011). Two analytical models for the study of periodic fibrous elastic composite with different unit cells. Compos. Struct., 93(2), pp. 709–714. DOI: 10.1016/j.compstruct.2010.08.008. [20] Yang, Y. and Misra, A. (2012). Micromechanics based second gradient continuum theory for shear band modeling in cohesive granular materials following damage elasticity. Internat. J. Solids Structures 49(18), pp. 2500-2514. DOI: 10.1016/j.ijsolstr.2012.05.024. [21] Khoroshun, L. P. and Shikula, E. N. (2012). Deformation and damage of linear elastic homogeneous and composite materials (review). International Applied Mechanics, 48, pp. 131–175. DOI: 10.1007/s10778-012-0512-3. [22] Berthier, E., Ponson, L., Dascalu, C. (2014). Quasi-brittle Fracture of Heterogeneous Materials: A Nonlocal Damage Model. Procedia Materials Science, 3, pp. 1878-1883. DOI: 10.1016/j.mspro.2014.06.303. [23] Dorhmi, K., Morin, L., Derrien, K., Hadjem-Hamouche, Z., Chevalier, J.P. (2020). A homogenization-based damage model for stiffness loss in ductile metal-matrix composites. Journal of the Mechanics and Physics of Solids, 137, 103812. DOI: 10.1016/j.jmps.2019.103812. [24] Sorić, J., Lesičar, T., Tonković, Z. (2021). On Ductile Damage Modelling of Heterogeneous Material Using Second- Order Homogenization Approach. Computer Modeling in Engineering and Sciences,126 (3), pp. 915-934. DOI: 10.32604/cmes.2021.014142. [25] Wu, L., Noels, L., Adam, L., Doghri, I. (2012). A multiscale mean-field homogenization method for fiber-reinforced composites with gradient-enhanced damage models. Computer Methods in Applied Mechanics and Engineering 233, pp. 164-179. DOI: 10.1016/j.cma.2012.04.011. [26] Devries, F., Dumontet, H., Duvaut, G., and Lene, F. (1989). Homogenization and damage for composite structures. Int. J. Numer. Meth. Engrg., 27, pp. 285–298. DOI: 10.1002/nme.1620270206. [27] Xia, Z., Chen, Y., Ellyin, F. (2000). A meso/micro-mechanical model for damage progression in glass–fiber/epoxy cross-ply laminates by finite-element analysis. Composite Science and Technology 60, pp. 1171–1179. DOI: 10.1016/s0266-3538(00)00022-1. [28] Mises, R V. (1913). Mechanik der festen Körper im plastisch-deformablen Zustand. Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse, pp. 582-592. [29] Mazars, J. (1986). A description of micro-and macroscale damage of concrete structures. Engineering Fracture Mechanics, 25(5-6), pp. 729-737. DOI: 10.1016/0013-7944(86)90036-6. https://doi.org/10.1016/j.cma.2012.04.011 https://doi.org/10.1016/j.cma.2012.04.011 https://doi.org/10.1002/nme.1620270206 https://doi.org/10.1002/nme.1620270206 https://doi.org/10.1016/0013-7944(86)90036-6 https://doi.org/10.1016/0013-7944(86)90036-6 K. Benyahi et alii, Frattura ed Integrità Strutturale, 58 (2021) 319-343; DOI: 10.3221/IGF-ESIS.58.24 343 [30] Mazars, J., Hamon, F., Grange, S. (2015). A new 3D damage model for concrete under monotonic, cyclic and dynamic loadings. Materials and Structures, 48(11), pp. 3779-3793. DOI: 10.1617/s11527-014-0439-8. [31] Bouafia, Y., Kachi, M S., Fouré, B. (2002). Stress-strain relationship in the case of reinforced steel fiber concrete. Annals of BTP. [32] Belhadj, N., Bouafia, Y., Smahi, R. (2015). Modeling the Behavior of the Fiber Reinforced Concrete by Damage Mechanics. In: Applied Mechanics and Materials. Trans Tech Publications, pp. 386-390. DOI: 10.4028/www.scientific.net/amm.749.386. [33] Nguyen, HG. (2008). Micromechanical approach for modeling the elastoplastic behavior of composites: application to resin mortars. Doctoral thesis. University of Cergy-Pontoise, France. [34] Duplan, F. (2014). Cementitious composites with controlled elastic modulus: design, characterization and micromechanical modeling. Doctoral thesis. University of Toulouse III - Paul Sabatier, France. [35] AL Kassem, G., Weichert, D. (2009). Micromechanical material models for polymer composites through advanced numerical simulation techniques. In: PAMM: Proceedings in Applied Mathematics and Mechanics. Berlin: WILEY‐VCH Verlag, pp. 413-414. DOI: 10.1002/pamm.200910180. https://doi.org/10.1617/s11527-014-0439-8 https://doi.org/10.1617/s11527-014-0439-8 https://doi.org/10.4028/www.scientific.net/amm.749.386 https://doi.org/10.1002/pamm.200910180 https://doi.org/10.1002/pamm.200910180