Micromechanics Based Inelastic and Damage Modeling of Composites

There are many different methods and tools that can be used to deliver the macroscopic constitutive response of heterogeneous materials from a local description of the microstructure behavior. Here we are concerned with non-linear behavior caused by inelasticity of the constituents or with the initiation and growth of damage. In developing the homogenization procedures for non-linear materials we have to define both the homogenization step itself (from local variables to overall variables) and the often more complicated localization step (from overall controlled quantities to the corresponding local quantities. As practical tools in the homogenization framework we have mainly three categories of methods involving the linearly elastic behavior of the composite aggregate:


Introduction
There are many different methods and tools that can be used to deliver the macroscopic constitutive response of heterogeneous materials from a local description of the microstructure behavior.Here we are concerned with non-linear behavior caused by inelasticity of the constituents or with the initiation and growth of damage.In developing the homogenization procedures for non-linear materials we have to define both the homogenization step itself (from local variables to overall variables) and the often more complicated localization step (from overall controlled quantities to the corresponding local quantities.As practical tools in the homogenization framework we have mainly three categories of methods involving the linearly elastic behavior of the composite aggregate: l First, we consider exact solutions and theorems that deliver variational bounds on the overall constitutive parameters. At that level, we can mention Voigt and Reuss bounds, Hashin-Shtrikman bounds in the linear context, and more recent works (Ponte Castaneda, 1991; Ponte Castaneda and Willis, 1995; Ponte Castaneda and Suquet, 1998), whicht give pertinent results for inelastic behavior.
l Third, numerical techniques, most often based on the assumption of microstructure periodicity.In these cases, which may be quite particular, periodic homogenization techniques (Suquet, 1998) are able to describe the local stress (and strain) fields and their evolutions very correctly.They also deliver the overall stress strain behavior of the considered representative volume element of the material.
Use is made of finite element methods or fast Fourier transform solutions.These numerical methods are limited to quasi-periodic situations.Moreover their cost is very high and their direct use in a true structural analysis is at present limited to very special cases.
Nonlinear problems of localization and homogenization are of great importance today.Not only classical composites suffer from deterioration of the material due to hereditary problems (aging, viscoelasticity).On the other hand, composite materials prepared in a special way can improve the properties of another materials, and the resulting effect can be much better than before.In this case, nonlinear and time dependent behavior has to be taken into account.
The present paper develops constitutive equations for inelasticity and damage of heterogeneous materials that benefit from some specificities of a special boundary element method.On the one hand, we need to obtain better approximations of the local stress and strain fields than are provided by Eshelby based approaches, especially when considering damage and failure conditions.On the other band, we want to simplify the numerical techniques of overall homogenization in order to obtain a treatable system of equations that can recover the status of a constitutive equation.
In this paper we restrict ourselves to the application and improvement of a method proposed initially by Laws (Laws, 1973, Dvorak, 1992).The most elaborated version, called transformation field analysis (TFA) (Dvorak and Benveniste, 1992), incorporates in the same framework thermo-elasticity, plasticity and viscoelasticity (Dvorak, 1992), or even piezoelectric-elasticity couplings (Benveniste, 1993).It can be used either with a low number of sub-volumes (or subphases), typically one subvolume per actual phase, then recovering the context of Eshelby type approaches, or with a larger number of sub-volumes, corresponding to some simplifications (or averaging) of the finite element based methods.
We note the main lines of the TFA method, and some of its properties for two-phase exploitations (Dvorak and Benveniste, 1992).Since we know the relative "stiffness" of the method when applied to a two-phase system (uniform fields on each phase), we propose and exploit a correction method that takes advantage of "tangent formulation" under asymptotic conditions (Section 3).Introducing damage effects through a continuum damage mechanism (CDM) formulation within each sub-volume, the models will be based on numerical results obtained via the overall strain over the external boundary of representative volume element (RVA) and the localization process.The boundary element method is used as a very powerful tool for problems of this kind.The method has proved to be appropriate particularly for solving the contact problem of debonding fibers from the matrix, and also for "softening" of the matrix due to damage, cracking, and other continuum phenomena.

Local and overall constitutive equations
We begin by the expression for the local elastic constitutive equations, assuming a uniform elastic stiffness L r over each sub-volume V r .The local stresses s r (x) and strains e r Prescribed distribution of local eigenstrains is denoted m r ( ) x and l m r r r is the corresponding eigenstress fieId.These eigenstrains can be thermal strains, pIastic strains, transformation strains.In the macroscale, the overall (uniform) strain E and stress S are also reIated by: where L is the overall elastic stiffness matrix and l and m are eigenstresses and eigenstrains.
In the case of pure elasticity, eigenstrains and eigenstresses may be given by a change of temperature, swelling, watering, prestressing, for example.They are not unknowns, but are given in advance.Then the following procedures are not of interest to us, as the expressions for internal states in the composite body can be calculated in a much simpler way.
In pure elasticity, concentration factors A and B are sought, such as: e r ( ) ( ) It is immediately clear that the relationship between the two concentration factors follows from (1), ( 2) and (3) as: Note that according to Dvorak and Benveniste, (1992), it holds: In this way the overall stiffness can be calculated, knowing the stiffnesses and strain concentration factors and volume fractions of both phases: c V V r z = .
Similarly, the eigenparameters can be determined from the stress and strain concentration factors and volume fractions as: Note that the Lippman-Schwinger equation is used in this text and will be derived in a slightly different way than usual.The Lippman-Schwinger equation relates strains, strains of a comparative medium and the eigenstrain (eigentress) field.Some particular implementations of this equation enable us to derive very important results, using the theory of generalized functions (distributions).The Lippman-Schwinger equation has the form: [ ] in which e 0 denotes the strain field (in our case considered to be uniform) that would exist in a comparable homogeneous medium L 0 under the same boundary conditions.The kernel G is defined by: [ ] where G ijkl is the Green function of the homogeneous medium L 0 obeying: where d ip is the Kronecker symbol and d( ) x -x is the Dirac continuous functional.

Localization and homogenization
In this section eigenstrain and eigenstress are dealt with and involved into the computation.No body forces are present.
Let us consider the coordinate system 0y 1 y 2 in 2D (for the sake of simplicity our restriction is to 2D, while the generalization to 3D is straightforward), bounded domain W (unit cell) with the boundary We define average quantities < .> by: where meas W stands for the volume (area) of W.
The homogenization can start by introducing the overall (average) strain E, or stress tensor S, or periodic conditions can be prescribed.The latter conditions are the most useful in applications.There are plenty of other boundary conditions which are less important than those we have mentioned.
From (10) we get, for example: Without loss of generality, we focus only on given E. Then we have to solve the problem: The real displacement u and the real strain e may be written in the form of the sum of E and the fluctuating terms ¢ u and ¢ e as: In the case of elastic behavior of both fiber and matrix and the fiber-matrix interface, it holds: The procedure for solving of e is split into two steps: First, let u 0 , e 0 , and s 0 be the known displacement, strain, and stress fields, respectively, defined on a comparative medium L 0 .The linear Hooke law relates the stresses and strains: s e Matrix L 0 is not yet fixed.An obvious option is to put e 0 = E and in the sense of the above definitions L 0 can be expected to be the average sought stiffness.It will be seen that this assumption is not so easy.
In the second step, a geometrically identical body is considered, which is heterogeneous, anisotropic, and may exhibit nonlinear behavior.The displacements u, the strains e, and the stresses s are unknown, and the generalized Hooke law including the eigenstresses l, or eigenstrains m, holds valid as: Similarly to [19] we define the symmetric polarization tensor t by: We also define e e e = -0 , s s s Our aim is to obtain the relations strain, or stresses and eigenstrain, or eigenstresses.Since both s and s 0 are statically admissible, the following equations have to be satisfied in the sense of distributions: where Subtracting ( 17) from ( 16) yields 3 Numerical derivation of quantities Solution of problems involving the linear or nonlinear behavior of composite bodies is mostly formulated in terms of integral equations.Consequently, a natural way to solve these problems is to describe the behavior of such bodies by the boundary integral equation method.
The integral equation equivalent to (21) to ( 24) can be expressed as: u u y py y p y uy y y where the starred quantities are the given kernels.Due to the validity of the boundary equation ( 22 ( The convected term X arises at the internal point of W by the interchange of integration and differentiation when deriving (27) from (26).Note that this unpleasant term may be avoided by introducing Eshelby's trick, which is very famous in the theory of composite materials.Levin in his work on termal bounds uses a similar trick leading to omitting the convected term.

Numerical interpretation of TFA
Now, our goal is to derive the strain -eigenstress relation of the form where A and G are the influence function tensors (A is mostly referred to as the mechanical concentration function tensor).Note that once computed, these matrices do not change their value during the incremental process for nonlinear solution of plasticity and damage.After discretization of the boundary and the domain W into M internal cells, we get the discretized form of the integral equations: where U is a square matrix (2N×2N) and 2N is a number of degrees of freedom on the boundary in 2D (when using, e.g., a linear approximation of both tractions and displacements, N is the number of nodal points), p is the vector (2N) of discretized unknowns tractions at the nodal points on ¶W , S, ¢ S are the matrices (2N×3M) of influences of the strains and eigenstrains in the discretized domain (in 2D three components of the strain tensor are independent), H is a (3M×2N) matrix, and, finally, Y and ¢ Y are square matrices (3M×3M).
Since the system is well-posed, the regular matrix U may be inverted.Elimination of p from both latter types of equations provides where Obviously, W is a regular (6M×6M) matrix, since for the given component of l a unique response e may be expected.The sought influence function tensor G is equal to W T -1 , while W -1 is the mechanical concentration function tensor A. Now we turn to concentration factors for the phases.There still is a certain freedom in selecting matrix L 0 .Let put l = 0 and then successively In the same way t in (21) changes, and the influence on the integral (25) to ( 27 In ( 36) and (37) it is important that the integration is led over only one phase and the second is not explicitly obtained in the formulas.This fact can principally speed up the computation, particularly if the fiber volume ratio is very small (fiber reinforced concrete), or is very large (classical composites).Because of well-known identity one concentration factor can be calculated from either (36) or (37), and the remaining follows from (38).Note that for an unlimited domain W the Lippman-Schwinger equation is fully fulfilled.The last relation makes clear why we concentrated our attention on splitting of the integrals into integration over the fiber and integration over the matrix.Because of (38) we can simply calculate the integrals either over the matrix, or fiber, exceptionally.Then, using (38), we get the second concentration factor needed.
Hence, in the case of elastic behavior, the homogenization is straightforward: and the elastic stiffness matrix appears to be: It is worth noting that the stress concentration factor can be derived in the same way: Note that in the general case relation ( 41) is calculated from the given S.This leads to similar integral equations, but the solution is not the same.Since the external forces are equilibrated, the rigid motion of the unit cell is disregarded, and the solution is then unique.Here we use relation (41) and the possibility to invert L * .
Including the eigenparameters into the calculation, we can derive from (35) to (37) the relations e m The relations (42) are the starting points for the theories established by Dvorak.He assumes that the concentration functions are estimated by approximate formulas following Mori-Tanaka, or the Self-consistent method.The calculations presented in this paper are considered to be very accurate and fast.

Viscoelasticity
In linear viscoelasticity it is always possible to write the stress-strain relation in a similar form to that of elasticity with the terms of the L * matrix now representing suitable differential or integral operators, rather than elastic constants.Thus in an isotropic continuum a pair of operators corresponding to an appropriate pair of elastic constants will appear -while for anisotropic material up to 9 separate operators may be necessary.
Typically, the viscoelastic part of the strain may thus be described by where each term of the viscoelastic matrix, L -1 , may take up a form of Kelvin chains.This, as is well known, can be interpreted as a response of the form: The increments of each such term in a time interval may now be found from the above expression from the knowledge of the current value of appropriate stress component s s and the current value of e n .Thus it becomes necessary to store only a finite number of such terms as e n at their current value to represent the full history effect.The viscoelastic strains can be treated as eigenstrains, value A representing spring and value B representing dashpot.

Damage properties
Damage effects are split into two parts: the first part is connected with exclusion of tensile stress in the principal direction, which overcomes some prescribed value (tensile strength s + ).Note that if the tensile strength is different from zero, the mathematical formulations lead to problems, the correct solution of which is not yet fully proved.Nevertheless, the mechanical results are reasonable and for practical applications they are treated as correct, at least from the mechanical standpoint.If the principal stress is too low, the stiffness is weaker.The tests are carried out in each cell of the discretization into the boundary elements and internal cells.The second criterion of damage is violation of the Mohr-Coulomb hypothesis in the principal shear direction.Because the outcome of such a violation is shear cracks -disconnection of displacements -the shear modulus is calculated for the new stage.

Example
Since the procedure leads to a linear relations at each stage, for given (step by step increasing) values of components of either stress tensor (elastic and relaxation part) or the strain tensor (elastic and viscoelastic part) of "unit impulses" we do not need to compute the influence tensors.Now the procedure fully described in [15] can be used.Note that in [15], the values of the concentration tensors, which are the most important quantities for numerical computation, are computed by very approximate methods (Mori-Tanaka, etc.).
A quarter of a unit cell is considered with a fiber volume ratio equal to 0.6, according to Fig. 2.
We used the following elastic material properties of the phases: Young's modulus of fiber E f = 210 MPa, Poisson's ratio n f = 0.16; for the matrix E m = 17 MPa, and n m = 0.3.For a fiber volume ratio 0.6, the radius of the fiber is r = 0. From the above matrix we can conclude that the responses on the normal unit strains are computed with high accuracy (comparing the symmetry), while the results from the shearing strains are less accurate, but still very precise.
The angle of internal friction was f = 25°, and cohesion c = 270 kPa.The results are depicted in Fig. 3.In the upper and lower part of the ilustrations damage is obvious (debonding occurred) and on both sides viscoelastic behavior prevails.
The computation was run on Pentium IV PC, 2.6 GHz in FORTRAN.The program for generating the meshes of the internal cells and also the boundary nodes had been prepared, as is clearly shown in Fig. 2. According to the wish of the user, the meshing can be improved.The time needed for computing even a large system of equations (150×150), which can be stored into memory without the use of a hard disc or an extended/expanded memory, was negligible in each step.Our illustration does not reach such dimensions of computation.It is also not necessary, in such problems, to increase the precision of the meshing, which would reduce the efficiency.The iteration at each step of loading was also very fast.It is worth noting that similar computations were carried out using FEM, but finer meshing had to be imposed to get comparable results with BEM in the procedure presented here.The comparison was tested in such a way that the sum of the concentration factors would be the unit tensor.

Conclusions
This paper has presented the fundamental idea of a numerical procedure leading to the overall viscoelastic and damage behavior of the matrix in a composite aggregate on a unit cell.Based on Transformation Field Analysis it is possible to obtain strain and stress/strain influence matrices, that relate the strains and the eigenstrains and eigenstresses.Since this problem leads to integral equations, the most suitable numerical tools appear to be BEM.The obvious advantage of this procedure is found in a priori computed influence matrices (concentration tensors).These may be stored into a computer and, hence, the iteration process for solving the nonlinear material behavior of the structures is very efficient.Moreover, attention may be focused on integrals over either the matrix or the fiber.
A very important property of the above procedure is the linearity of the problem at each stage.The accuracy of the overall stiffness (compliance) matrix is not dependent on the size of the step, providing there is no "unwanted" unloading in any internal cell at the current stage.When this is the case, the time is slightly extended, as the standard iterative process has to be carried out.This was not the case for our computations.
Although it is not our intention to discuss cases of eigenstrain fields, it is appropriate to mention the connection between the present formulas and those obtained in studies of the thermoelastic response in composite materials subjected to a uniform change of temperature.
the second integral on the right hand side of (25) disappears.Therefore, we get: equation successively by one of the coordinate x , we arrive at the expression 714.The homogenized elastic matrix L * in this case possessed the following values:

Fig. 2 :
Fig. 2: Geometry of boundary elements and internal cells