Acta Polytechnica Vol. 44 No.5–6/2004 On Adequacy of Two-point Averaging Schemes for Composites with Nonlinear Viscoelastic

Finite element simulations on fibrous composites with nonlinear viscoelastic response of the matrix phase are performed to explain why so called two-point averaging schemes may fail to deliver a realistic macroscopic response. Nevertheless, the potential of two-point averaging schemes (the overall response estimated in terms of localized averages of a two-phase composite medium) has been put forward in number of studies either in its original format or modified to overcome the inherited stiffness of classical ”elastic” localization rules. However, when the material model and geometry of the microstructure promote the formation of shear bands, none of the existing two-point averaging schemes will provide an adequate macroscopic response, since they all fail to capture the above phenomenon. Several examples are presented here to support this statement.


Introduction
It is in a nature of mankind to search for simplicity, efficiency and stability.When translated into a computational mechanics language, a new landscape for efficient numerical schemes arises from the introduction of the multi-scale or multilevel solution strategies currently at the forefront of engineering interest when studying complex heterogeneous materials and structures.To provide an illustrative example of such a structure, consider the multi-layered wound composite tube shown in Fig. 1.An accurate prediction of the mechanical response of this particular class of structures inevitably calls for analyses on three widely separated length scales.It is now becoming widely accepted that models constructed on the basis of hierarchical or multi-scale modeling offer a reliable route in numerical investigation of deformation and failure processes taking place at individual scales [1][2][3][4][5].
Since their introduction, the need for a feasible computational framework has challenged the applicability of various numerical techniques at individual scales.The high cost of traditional numerical techniques, e.g., the finite element method, then provided an opportunity for classical averaging schemes as a cost-effective and computationally simple alternative.The motivation for the present paper, in particular, arises from the idea of using the well-known variational principles of Hashin and Shtrikman (HS) [6] when studying the behavior of a composite on the level of individual constituents, Fig. 1c.Assuming that the material response of a two-phase composite system is well described by volume averages of local fields, all these methods including the HS principles can be conveniently referred to as two point averaging schemes.A comprehensive overview of various micromechanical techniques can be found in [7].An extension to the loading conditions that promote inelastic deformation is presented in [8][9][10][11][12][13][14].
A number of studies, however, have revealed the main drawback of the so called "elastic localization rule" for the evaluation of local stress and strain averages with the two-point averaging schemes.In particular, defining the localization rules based on the elastic or tangent behavior of individual phases [8] yields a macroscopic response that is significantly stiffer than that provided by the finite element method or a sufficiently refined transformation field analysis [13].A number of approaches have been proposed to address this task with the main goal of improving the way the macro-

Micro-scale m » m
Fig. 1: An example of three-scale modeling scopic strains or stresses are redistributed into the individual phases.This led to the appearance of alternative methods, based on different linearization of the governing non-linear problem such as the secant, modified secant and affine methods [11,12].Another example are the second-order variational estimates for material systems described by convex potentials [11].A drawback of the previously mentioned methods when presented in the framework of incremental formulation is the need for the solution of a nonlinear system of equations for each load increment.A remedy has been offered in [13] where the classical "total localization rule" presented in the framework of transformation field analysis but enhanced by corrected values for the eigenstrains to soften the localization rule was proposed.Although appealing in the macroscopic response that they deliver, most of the above methods require the existence of instantaneous or asymptotic tangent stiffness, which is not always available, so that an extension to more complex and/or rate-dependent constitutive laws might not be an easy task [14].In this work in particular the authors examined the use of the two-point averaging scheme in conjunction with the extended HS principles for modeling of unidirectional composites with a disordered microstructure.When allowing for a nonlinear viscoelastic response of the matrix phase, a rather surprising phenomenon was observed suggesting severe limitations in the application of two-point averaging schemes.
To introduce the subject we recall the essential conclusions put forward in [14] that illustrate the well-known drawback of the so called "elastic localization rule" for the evaluation of local stress and strain averages in a two-phase medium.To begin with, consider the results presented in Fig. 2 featuring macroscopic stress-strain curves derived for hexagonal arrangement of fibers under transverse shear strain loading, see [14] for more details.The inability of the HS principle, when used in its original form, to correctly capture the stress redistribution is further revealed in Fig. 3, which shows plots of "localized" phase averages.
For the present material system (matrix response described by the generalized Leonov model), the deficiency of the family of elastic localization rules can be attributed to the fact that the significant non-uniformity of local fields, which manifested itself in an evolution of the shear bands as will be shown later in Section 5, cannot be accurately represented by the piecewise uniform variation of the local linearized modules.
To further support this statement, we consider an example of a two-phase laminate having uniform distribution of local fields in individual phases, so that the piecewise uniform approximations fit exactly, thus suggesting a perfect match with the response provided by the standard finite element method.Fig. 4 shows a variation of the matrix shear stress due to overall shear strain rate & E 12 4 1 10 = -- s for several values of the fiber volume fractions c f .The results plotted in Fig. 4 indicate an increase of the matrix strain rate for higher values of the fiber volume fraction manifested here by higher values of the plateau stress.This supposition is confirmed in Fig. 5 showing a time evolution of the matrix shear strain.It is worth noting, however, that a similar conclusion cannot be drawn from Fig. 2, when the FE method is used to derive the volume averages of the local fields.Clearly, the volume average of the matrix phase for the composite plots below the curve is obtained for a pure matrix subjected to the same loading conditions.This particular result is in direct contradiction with the observations gained from Figs. 4-5.A sound explanation of this behavior is therefore needed.
Similar experiments were examined in [14] with reference to the primary HS principle in the framework of the two--point integration scheme augmented with a special choice of the reference medium.The results summarized in Fig. 6 provide comparisons of stress-strain curves obtained for this formulation.Clearly, for material systems with uniform fields the modified incremental elastic localization rule is sufficient for all values of fiber volume fraction c f .This result just demonstrates that the resulting mismatch between the HS principles and the finite element simulations, Fig. 2, is a consequence of non uniformity of the involved fields appearing in a composite.This phenomenon is examined below particularly with reference to geometrically more complex microstructures, Fig. 1.
In achieving this goal the behavior of so-called statistically equivalent unit cells is studied.The derivation of such a unit cell is reviewed in Section 2. The first order homogenization scheme is then addressed in Section 3 in the framework of the finite element method.The basic ingredients of the used constitutive model are presented in Section 4. The paper concludes with a detailed discussion of the studied phenomena presented in Section 6.
In the following text, lowercase boldface letters, e.g., a are used for column vectors while capital boldface letters, e.g.A will be used to denote a matrix.Lightface letters are used for scalar quantities, e.g., a. Specific dimensions of individual quantities follow from the context.The inverse of a non-sin-gular matrix is denoted as A -1 while the superscript T indicates the transpose of a matrix.Finally, only the loading due to constant overall strain rate is considered and the analysis is carried out under the generalized plane strain conditions [24, Appendix B] with x 3 being the axis of the fibers.

Geometrical modeling: construction of statistically optimized unit cells
This section offers a certain statistically optimized periodic unit cell (PUC) consisting of only a small number of particles as a suitable representative of a real microstructure.Such a unit cell is found through a minimization procedure formulated in terms of certain statistical descriptors characterizing the geometrical configuration of a random medium.It is believed that if the two systems, the real microstructure and the corresponding unit cell, are the same in the statistical sense then the mechanical response of both systems will also be the same.This idea has been successfully exploited in a number of the authors' previous works [15][16][17][18][19]. See also [20] for other routes to tackle disordered microstructures.
It has been shown that successful completion of this step requires some measures for a reliable quantification of the random microstructure.To do so it is convenient to introduce the concept of an ensemble -a collection of a large number of systems, which are identical from the macroscopic point of view and different in their microscopical details.The morphology of such a composite system is fully characterized by a random function c a r ( , ) x , which is equal to one when a point x lies in the phase r within the sample a, and equal to zero otherwise.With the aid of function c r , the n-point probability function S r1, …, rn can be defined as [21,22] S r rn n r rn n where • denotes the ensemble average.Thus, S r1, …, rn gives the probability of finding n points ( , , ) x x 1 K n randomly thrown into a medium located in the phases r 1 , …, r n .In the following, we limit our attention to functions of the order of one and two.
Analysis of random composites usually relies on various hypotheses, which simplify the computational effort to a great extent.In particular, under the ergodic hypothesis [21] the volume average of function c a r ( , ) x given by is independent of a and identical to the ensemble average where c r is the volume fraction of the r th phase.The statistical homogeneity assumption means that the value of the ensemble average is translation invariant.Then, for example, the two-point matrix probability function reads where x x x ij j i = -.Note that for an ergodic and periodic microstructure, the two-point probability function S rs has the following form where • now denotes the complex conjugate.When introducing a binary image of the actual microstructure and taking into account the periodicity of the RVE we may approximate Eq. ( 6) by the discrete Fourier transform: ( ) where N x and N y are the horizontal and the vertical resolution of the bitmap.Note that this approach requires only ) operations needed by the direct method.Moreover, the possibility of using highly optimized public software packages makes the DFT-based approach very efficient.See, e.g., [16] for detailed discussion and numerical experiments.
Having the original microstructure characterized in an appropriate sense, the construction of an equivalent PUC is relatively straightforward.In particular, the PUC is derived by matching a selected microstructure describing function of the real microstructure and the PUC.To be more concrete, consider a periodic unit cell consisting of N particles.Dimensions H 1 and H 2 and the x and y coordinates of all particle centers determine the geometry of such a unit cell.The particle locations together with an optimal ratio of cell dimensions H 1 /H 2 are found by minimizing an objective function involving two-point probability function where x y x y = 1 1  , , , , K stores the positions of the particle centers, S r s i i 0 ( , ) is the value of a two-point (matrix) probability function corresponding to the original medium evaluated at a point ( , ) r s i i , and N m is the number of points in which both functions are matched.
A closer inspection of the objective function (8) reveals that it is discontinuous and multimodal with a large number of local plateaus.This is a direct consequence of using bitmap images, where individual entries of the searched vector are integer variables.Based on our previous experience with optimization of these functions [17], the Augmented Simulated Annealing Method [19] is implemented to solve the optimization problem.Examples of the resulting 5-and 10-particle PUCs are displayed in Fig. 7.The hexagonal unit cell with the same volume fraction is shown for comparison.

Numerical modeling: finite element discretization
In this section, the numerical analysis is performed in the spirit of the first order homogenization of periodic fields [1][2]23].To that end, consider a representative volume element Y in terms of one of the statistically optimized unit cells derived in the previous section.Next, let the applied loading conditions produce a uniform distribution of macroscopic strain E or the macroscopic stress S fields.When further assuming a periodicity of microstructure, the local displacement field u(x) and strain field e(x) admit the following decomposition where u * (x) and e * ( ) x represent fluctuations of the local fields due to the presence of heterogeneities [21,23].The local stresses are related to the strains by the constitutive equation where L(x) is the material stiffness matrix and l(x) is the vector of eigenstresses (strain-free stresses).Note that l( Denoting K ij individual blocks of the left-hand side of Eq. ( 12), we can formally eliminate the unknown vector r from ( 12) and obtain the resulting macroscopic homogenized constitutive law where When analyzing the non-linear material system, Eq. ( 10) needs to be replaced by an appropriate linearized constitutive law, and relations ( 13)-( 14) are, rather straightforwardly, replaced by an incremental counterpart.This step is discussed in the next section.
Similar results can also be derived with help of other numerical techniques such as the boundary element method [25].

Constitutive modeling: generalized Leonov model
As suggested in the introductory part, a graphite/epoxy material system is selected as a particular example in this study.Note that the fiber is assumed to remain elastic during deformation so that the inelastic effects are limited to the matrix phase.For the composite structure plotted in Fig. 1, PR100/2+EM100E epoxy is used as a bonding agent.An experimental program carried out on this type of material [27,29] demonstrates that the relevant rate dependent response of the epoxy as well described by the generalized Leonov model. Combining where h is the shear-dependent viscosity.In Eq. ( 15), A and t 0 are material parameters; a s that appears in Eq. (16b) is the stress shift function with respect to the zero shear viscosity h 0 (viscosity corresponding to an elastic response).Clearly, the phenomenological representation of Eq. ( 16a) is the Maxwell model with variable viscosity h.Note that a single Leonov mode is not able to describe a realistic response of real polymers, as it only accounts correctly for the initial stiffness and the strain rate dependent yield stress.To describe the multi-dimensional behavior of the material, the generalized compressible Leonov model, equivalent to the generalized Maxwell chain model, can be used.The viscosity term corresponding to the μ-th unit receives the form where t eq is the equivalent shear stress and s ij is the stress deviator tensor.Admitting only small strains and isotropic materials, a set of constitutive equations defining the generalized compressible Leonov model can be written as where s m is the mean stress, e n is the volumetric strain, K is the bulk modulus, G m is the shear modulus of the m-th unit and d ij is the Kronecker delta.The shear modules G m related to relaxation times q m (h q m m m 0 = ×G ) then provide the time dependency of the model.They are usually found by transforming the experimentally derived creep function into the relaxation function using the Laplace transform combined with the method of partial fractions [27].Coefficients of the creep function were determined from a set of creep experiments performed at various stress levels.Thirteen elements of the generalized Kelvin-Voight chain model were used to obtain an accurate description of the linear compliance function.The resulting coefficients for the PR100/2+EM100E epoxy needed in the generalized Leonov model are stored in Table 1.
Modeling the mechanical behavior of a bonding agent in polymer matrix composites requires a reliable and stable procedure for integration of the set of governing equations ( 18)-( 21).To avoid possible numerical instabilities linked to integration schemes a fully implicit Euler backward integration procedure was developed.Providing the total strain rate is constant during integration a new state of stress in the matrix phase at the end of the current time step assumes the form where t i is the current time at the end of the i-th time incre-ment; s m (t i ) is the elastic mean stress, s(t i ) stores the deviatoric part of the stress vector s(t i ) and De is the deviatoric part of the total strain increment.With reference to the backward Euler integration scheme the time dependent variables at time t i receive the form where s m (t i ), m = 1, 2, …, M, is the deviatoric stress vector in individual units of the Maxwell chain model evaluated at the beginning of a new time increment Dt t t i i = --1 , and M is the assumed number of Mawell units in the chain model.The stress shift factor is given by where the equivalent stress t eq follows from and Clearly, the backward Euler step makes all variables nonlinearly dependent on the stress values found at time t i .Therefore, successful completion of a given integration step requires the solution of a system of nonlinear equations.Here, the solution is established employing the Newton--Raphson method.To that end, define a set of residuals with the primary variables Note that the current increment of the vector Dl(t i ), which appears in Eq. ( 23), is considered as a secondary variable.Then, under the condition that De is constant, the Newton-Raphson iterative scheme reads a a r where Jacobian matrix H is given by The initial values of primary variables at time t i for k = 0 are set to forward Euler estimates.More details about the numerical implementation including comparison of an explicit and implicit integration scheme can be found in [27,28].
The accuracy of the proposed numerical procedure is demonstrated in Fig. 9, which shows a typical uniaxial re-Acta Polytechnica Vol.44 No. 5-6/2004 Table 1: Nonlinear viscoelastic material properties of the PR100/2+EM100E epoxy matrix sponse of the PR100/2+EM100E epoxy resin for three different constant strain rates, where the solid lines were obtained experimentally, while the others follow from the numerical analysis.

Inadequacy of two-point averaging
Recall that, in Section 2, rather encouraging results were obtained for (an appropriately modified) two-point averaging scheme for binary laminates.This scheme, however, when applied to more complex microstructures such as hexagonal packing of fibers delivers rather inadequate results, although having a significantly softer prediction of the macroscopic response than the classical one, see Fig. 3.An explanation is offered by examining, among others, also the response for a pure matrix subjected to the same overall shear strain rate as the composite.It is interesting to observe, particularly in view of the results presented in Figs.4-6, that the stress-strain curve for a pure matrix plots above the finite element estimates of the volume average of matrix stress derived for the composite.
To reconcile this "paradox", namely the fact that the volume average of the matrix yield stress decreases with an increase in the volume average of the matrix strain rate, we can set up a very simple geometrical model with the finite element discretization depicted in Fig. 10.Elements that exhibit the same stress response were marked with the same number.Note that 6-noded elements with the 7-point integration rule were used so that the results displayed in Figs.11-12 correspond to element averages.Fig. 11 indeed confirms large non-homogeneity in the distribution of local fields.It also shows an increase in the rate of deformation of the volume average of matrix shear strain in comparison with the response of the pure matrix under the same loading conditions.This is clearly due to the heterogeneity of the material system and is consequently due to the required compatibility of local and overall fields.This condition inevitably leads to an increase in reference yield stress of the matrix phase provided by the HS procedure.Recall also the laminate response studied in the first example, Figs. 4-6.The expected increase in the volume average of the matrix yield stress, however, is not observed with the finite element calculations even for such a simple geometry.This can be attributed to the highly nonlinear dependency of the yield stress on the applied rate of deformation provided by the Leonov model.Also note that the elements found in a closer  vicinity of the fiber experience a considerably slower response that those located away from the fiber phase.This eventually results in a softer response in comparison with the pure matrix.Clearly, such a behavior cannot be represented by any of the existing localization rules dealing with two-phase material systems and piecewise uniform variation of local fields, at least for the present material model.The conclusions of this case study can, in principle be, extended by analogy to more complex periodic unit cells such as those displayed in Fig. 7. To illustrate the governing mechanism responsible for this fact, the distribution of local equivalent strain in a hexagonal array, 5-and 10-fiber unit cells is shown in Fig. 13.Each microstructure was discretized with 6-noded triangular elements with 7 integration points.In particular, the plotted local fields correspond to the value of overall deformation E 12 = 0,1.
For all the studied microstructures, the final deformation pattern has the character of a localized shear layer that is responsible for the plateau regions clearly visible in the graphs depicted in Figs.3-4.This behavior simply cannot be reliably captured by a two-point averaging scheme and finally leads to erroneous response of the simplified method.Note that a quite similar character of the overall response appears even for discretization of the structure with three-noded constant strain triangular elements, see Fig. 4. In this case, however, the approximation is not rich enough to correctly capture the fact that the formed shear band prohibits interaction of the fibers with the matrix phase, Fig. 15.An Additional explanation comes from the basic assumption of the selected material model.Recall that the model draws on the nonlinear viscoelastic response to be limited purely to the deviatoric part of the stress-strain relationship, while the volumetric response remains elastic.As often observed with J2 plasticity, this may eventually lead to what we call volumetric locking and consequently to substantially stiffer response when compared to reality.This, in particular, suggests that even when using a refined variant of the TFA method -"multi-point incremental homogenization" [10] -the global response will probably not be captured properly.

Conclusions
The paper reviewed essential steps of first order homogenization procedure implemented in the framework of the finite element method.A family of so-called statistical equivalent unit cells was introduced to deal with random microstructures.Unlike the elastic analysis, where a 5-fiber unit cell was found sufficient to arrive at reliable estimates of homogenized properties [16,17] at least a 10-fiber unit cell is needed when highly non-homogeneous evolution of local fields may occur.Recall the formation of shear bands displayed in Figs. 13, 14 and corresponding overall response plotted in Fig. 15.The existence of shear bands, zones of highly localized deformation, on the other hand, may raise a question about mesh objectivity when local formulation is used.It is worth noting that in the present study the size of the shear band is explicitly given by the underlying microstructure.Clearly, when refining the finite element mesh such a shear band will be captured more accurately but will not depend on the size of the elements and will always attend a finite depth.Regularization from the mathematical point of view is thus not needed.Nevertheless, problems might be encountered for relatively coarse meshes or for microstructures with a rather low value of the fiber volume fraction.In this regard, owing to the material model used in this study a phenomenon known as volumetric locking occurred for 3-noded triangular elements, recall Figs.14, 15 and the discussion in Section 6.
As for the main goal of this paper, the analysis clearly uncovered major drawbacks of two point averaging schemes when applied to random microstructures with possibly inelastic phases.As already suggested, the existence of regions of highly localized deformation attributed heterogeneity of a microstructure may considerably influence the overall response (this phenomenon can be further magnified if one of the phases is elastic while the other experiences an inelastic response).Here we refer to the results presented in Figs.11, 12 identifying the mechanism of the actual localization rule and clearly suggesting an inadequacy of the "total elastic" localization rule used with the present form of the HS principles.The two-point averaging schemes, if at hand, should therefore be used with caution.Nevertheless, if the "exact" local response is not of the main interest the parameters of a given material model can be adjusted to meet the desired overall response.In such a case, numerical experiments on refined geometry are usually used in place of laboratory measurements.An example of combining such an approach with a two-point averaging scheme has been presented in [14].

Fig. 5 :Fig. 6 :
Fig. 5: Motivation: Laminate response for strain rate & E 12 4 1 10 = -- s -evolution of the matrix strain rate as a function of the volume fractions

Fig. 7 :
Fig. 7: Examples of statistically optimized periodic unit cells the Eyring flow model for the plastic component of the shear strain rate d shear strain rate de e /dt yields the one-dimensional Leonov model[26] The material parameters required by the generalized compressible Leonov model are parameters A, t 0 and G m (h 0m ).The material parameters A and t 0 describe the stress dependency of the model and can be derived from the Eyring plot, Fig.8(b), assuming that at yielding the overall deformation as Acta Polytechnica Vol.44 No.5-6/2004