Periodic homogenization and damage evolution in RVE composite material with inclusion

This work deals with the coupling between a periodic homogenization procedure and a damage process occurring in a 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, and 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.


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 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][3][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 straingradient (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 multicoated 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 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]: with: ij j x : Represents a linear displacement field, ' () ii ux : is a periodic function representing a modification of the linear displacement field, because to the heterogeneity for the composite structure. 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: The indices to identify the th k pair of two opposite parallel boundary surfaces of a representative volume element (RVE). The periodic function ' () ii ux is considered to be the same for two parallel borders, the difference between the Eqn. 2 and Eqn. 3 is given as follows: 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: The strain energy of a homogeneous composite material over the RVE, is as follows: The stored strain energy of a heterogeneous composite material over the RVE is as follows: with:  ij : The strain tensor. From Eqn. 7 and Eqn. 10, we can therefore write the energy difference, as follows: with: i u : The i th average displacement. and using the equilibrium equation: By replacing Eqn. 12 with Eqn. 11, we obtain: Using Gauss's theorem, we can transform the integration on the volume into the integration on the surface, as follows: Then, at the level of the surface j S , we have: We can conclude that: 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 ; ; ; ; ; 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: Equations for applying periodic boundary conditions at the faces: The mean strain in Eqn. 6 can be transformed as follows: (21) In order to evaluate the effective properties ij C , 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: with: There are only five independent material constants ( ) A material is isotropic when its property is the same in all three directions. In this case, ( )  2  12 33  13  21  12  2  2  2  2  22  33 11  11 13  12 33  12 13 C .C -C S =-=-E C .C -2.C .C -C .C +2.C .C (26) ( ) For an orthotropic/isotropic case transversely and a uniaxial stress-state along the direction X-X, the effective elastic properties are determined:

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.

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: The damage is controlled by the equivalent strain  eq , which allows to translate a tri-axial state by an equivalence to a uniaxial state. The strain tensor in the principal coordinate system: Knowing that the positive part + is defined, such that if  i is the main strain in direction i:

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 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: The maximum ultimate stress of the composite (function of the characteristics of inclusions): The reference length is linked to the height h of the section: The initial modulus of the composite in traction is given by: The ultimate strain corresponding to the total mobilization of the inclusions-matrix adhesion, is given by: When there is tearing of inclusions, the breaking strain of the composite, is given by: 330 This strain is a consequence of opening the cracks too large. This fracture strain of the inclusions is limited as follows: The damage is progressive in the field of strains. For a uni-axial state of tension, the equivalent strain the nonlinear field is such that:

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):

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 (   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 [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 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.

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 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.
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). 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][13][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.   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 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 eff E 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 t A and made the parameter t B varied.       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.
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 t A and we varied the parameter t B under different volume fraction, shows us that the variation of the volume fraction of shaped inclusions (spheroidal, ellipsoidal), and of the parameter t B 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 t B corresponding to this percentage. Shows us that the variation of the parameter t A 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). n j : The unit normal, u i : The i th displacement, S j : The outer limit of the RVE, E 11 , E 22 and E 33 : The moduli of elasticity of extensions along directions 1, 2 and 3, respectively, ν ij (i, j = 1, 2, 3) : Poisson's ratios, G 12 , G 13 , and G 23 : Shear moduli, U i P : The nodal variable at the node P, of degree of freedom i, E : Hooke's matrix, D : The damage variable, ε e : Elastic strain, E b0 : Modulus of longitudinal elasticity of concrete, ϖ : Percentage by volume of fibers, θ 0 : Fiber orientation coefficient, l r : Reference length, β: Model constant, h : Composite cross section height, l f : Fiber length, E f : Elastic modulus of the fiber, n : Fiber-concrete equivalence coefficient, τ u : Ultimate fiber-concrete bond stress, ϕ : Diameter of a fiber, f ft : Composite tensile strength, ε ft : Composite cracking strain, ε rf : Fiber breaking strain, ε rt : Breaking strain of the composite in tension.