AP04-Bittnar2.vp 1 Introduction It is in a nature of mankind to search for simplicity, effi- ciency 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 mechani- cal 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 re- liable route in numerical investigation of deformation and failure processes taking place at individual scales [1–5]. Since their introduction, the need for a feasible com- putational framework has challenged the applicability of various numerical techniques at individual scales. The high cost of traditional numerical techniques, e.g., the finite ele- ment 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 individ- ual 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 aver- aging schemes. A comprehensive overview of various micro- mechanical techniques can be found in [7]. An extension to the loading conditions that promote inelastic deformation is presented in [8–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 lo- calization rules based on the elastic or tangent behavior of individual phases [8] yields a macroscopic response that is sig- nificantly 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- 200 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No.5–6/2004 On Adequacy of Two-point Averaging Schemes for Composites with Nonlinear Viscoelastic Phases J. Zeman, R. Valenta, M. Šejnoha 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. Keywords: Fiber-reinforced composite materials, microstructure, nonlinear viscoelastic behavior, Leonov model, energy methods, finite element modeling. Macro-scale m» Meso-scale mm» 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-lin- ear 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” pre- sented 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 ex- tension 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 phenome- non was observed suggesting severe limitations in the applica- tion 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 ar- rangement 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 redis- tribution is further revealed in Fig. 3, which shows plots of “lo- calized” phase averages. For the present material system (matrix response de- scribed 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 exam- ple 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 �E12 4 110� � �s for several values of the fiber volume fractions cf. 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 ob- tained for a pure matrix subjected to the same loading condi- tions. 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 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 201 Acta Polytechnica Vol. 44 No.5–6/2004 Fig. 2: Motivation: macroscopic response for �E12 4 110� � �s – hexagonal packing Fig. 3: Motivation: localized response for �E12 4 110� � �s – hexagonal packing Fig. 4: Motivation: Laminate response for strain rate �E12 4 110� � �s the reference medium. The results summarized in Fig. 6 pro- vide 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 cf. 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 appear- ing 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 con- cludes 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 indi- cates 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 x3 being the axis of the fibers. 2 Geometrical modeling: construction of statistically optimized unit cells This section offers a certain statistically optimized peri- odic 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 charac- terizing 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–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 mor- phology of such a composite system is fully characterized by a random function � �r ( , )x , which is equal to one when a point x lies in the phase r within the sample �� and equal to zero otherwise. With the aid of function � r , the n-point probability function Sr1, …, rn can be defined as [21, 22] Sr rn n r rn n1 1 1 1, , ( , , ) ( , ) ( , ) ,� � �x x x x� � � � � (1) where � denotes the ensemble average. Thus, Sr1, …, rn gives the probability of finding n points ( , , )x x1 � n randomly thrown into a medium located in the phases r1, …, rn. 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 vol- ume average of function � �r ( , )x given by � � �r V r V V ( , ) lim ( )x x y y� � �� � 1 d , (2) is independent of � and identical to the ensemble average � �r r r rS c( ) ( ) ( )x x x� � � , (3) where cr is the volume fraction of the r th phase. The statistical homogeneity assumption means that the value of the ensem- ble average is translation invariant. Then, for example, the two-point matrix probability function reads S Smm mm( , ) ( )x x x1 2 12� , (4) where x x xij j i� � . Note that for an ergodic and periodic microstructure, the two-point probability function Srs has the following form Srs r s( ) ( ) ( )x x x y y� �� 1 � � � � d , (5) 202 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No.5–6/2004 Fig. 5: Motivation: Laminate response for strain rate �E12 4 110� � �s – evolution of the matrix strain rate as a func- tion of the volume fractions Fig. 6: Motivation: Comparison of FEM and HS-based estimates for strain rate �E12 4 110� � �s where � represents the size of the analyzed domain. The Fourier transform of function Srs is given by F F FSrs r s( ) ( ) ( )� � � � �� 1 � , (6) where � now denotes the complex conjugate. When introduc- ing 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: � �S N N mm x y m m� 1 IDFT DFT DFT( ) ( )� � , (7) where Nx and Ny are the horizontal and the vertical resolu- tion of the bitmap. Note that this approach requires only O (NxNy) log (NxNy) operations, instead of O(N Nx y 2 2) opera- tions 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. Dimen- sions H1 and H2 and the x and y coordinates of all particle centers determine the geometry of such a unit cell. The parti- cle locations together with an optimal ratio of cell dimensions H1/H2 are found by minimizing an objective function involv- ing two-point probability function F H H S r s S r sN i i i i i Nm ( , , ) ( ( , ) ( , ))x 1 2 0 2 1 � � � (8) where � �x N N Nx y x y� 1 1, , , ,� stores the positions of the particle centers, S r si i0( , ) is the value of a two-point (matrix) probability function corresponding to the original medium evaluated at a point ( , )r si i , and Nm 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 opti- mization of these functions [17], the Augmented Simulated An- nealing 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. 3 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 condi- tions produce a uniform distribution of macroscopic strain E or the macroscopic stress � fields. When further assuming a periodicity of microstructure, the local displacement field u(x) and strain field �(x) admit the following decomposition u x E x u x x E x( ) ( ), ( ) ( )* *� � � �� � , (9) where u*(x) and �*( )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 � �( ) ( ) ( ) ( )x x x x� �L � , (10) where L(x) is the material stiffness matrix and �(x) is the vec- tor of eigenstresses (strain-free stresses). Note that �(x) can present a variety of physical phenomena like temperature or moisture effects, creep stresses, etc. Employing the Principle of Virtual Work (the Hill lemma) in the form � � � � � � � � � T *T T T ( ) ( ) ( ( ) )( ( ) ( ) ( ) , x x x E x x x E � � � � � L (11) and introducing the standard finite element approximation u x x r x x r* *( ) ( ) , ( ) ( )� �N B� [24], we finally obtain the discre- tized equilibrium equations L L B B L B L B ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) x x x x x x x x x x x x d d d dT T Y Y Y Y � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � E r x x x x x � � � � ( ) ( ) ( ) d dT Y Y B � � � � � � � � � � � . (12) © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 203 Acta Polytechnica Vol. 44 No.5–6/2004 Fig. 7: Examples of statistically optimized periodic unit cells Denoting Kij 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 � �� �LFE FEE (13) where L K K K K K K B FE T FE Td � � � � � � � 1 1 11 12 22 1 12 12 22 1 Y Y Y ( ), ( )� � x x ( ) ( ) .x x x� d Y � � � � �� ! " "" (14) 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 dis- cussed in the next section. Similar results can also be derived with help of other numerical techniques such as the boundary element method [25]. 4 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 ma- trix 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 mate- rial [27, 29] demonstrates that the relevant rate dependent response of the epoxy as well described by the generalized Leonov model. Combining the Eyring flow model for the plastic compo- nent of the shear strain rate d d e t A p � 1 2 0 sinh � � , (15) with the elastic shear strain rate dee / dt yields the one-dimen- sional Leonov model [26] d d d d d d d d d d e t e t e t e t e t e p e p � � � � � 2 , (16a) � � � � � � d d ( e t a p � �0 0 0 0 sinh , (16b) where is the shear-dependent viscosity. In Eq. (15), A and � � are material parameters; a� that appears in Eq. (16b) is the stress shift function with respect to the zero shear viscosity 0 (viscosity corresponding to an elastic response). Clearly, the phenomenological representation of Eq. (16a) is the Maxwell model with variable viscosity . Note that a single Leonov mode is not able to describe a realistic response of real polymers, as it only accounts cor- rectly 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 re- ceives the form � �� � �� �0 1 2, ( ),a eq eq ij ijs s , (17) where �eq is the equivalent shear stress and sij is the stress de- viator tensor. Admitting only small strains and isotropic ma- terials, a set of constitutive equations defining the generalized compressible Leonov model can be written as � �m K� , (18) � ( � � ), ,sij ij ij pG e e� � �� �2 , (19) � ( ) ( ), , , , e s s ij p ij eq ij eqa � � � � � � � � � � 1 2 1 2 0 , (20) � � � � � ij m ij ij ij ij M � � � � s s s, , 1 , (21) where �m is the mean stress, � is the volumetric strain, K is the bulk modulus, G� is the shear modulus of the �-th unit and �ij is the Kronecker delta. The material parameters required by the generalized compressible Leonov model are parameters A, �0 and G�( 0�). The material parameters A and �0 describe the stress depend- ency of the model and can be derived from the Eyring plot, Fig. 8(b), assuming that at yielding the overall deformation as 204 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No.5–6/2004 Fig. 8: (a) uniaxial tensile experiments, (b) the Eyring plot equal to the plastic deformation. Experimental measure- ments needed to calibrate the Eyring flow model appear in Fig. 8(a). The shear modules G� related to relaxation times �� � �� � �0 � # G ) then provide the time dependency of the model. They are usually found by transforming the experi- mentally 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 explicit 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 � � �m i m it t K( ) ( )� ��1 � , (22) s s e( ) ( ) �( ) ( )t t G t ti i i i� � ��1 2 Q� �� , (23) where ti is the current time at the end of the i-th time incre- ment; �m(ti) is the elastic mean stress, s(ti) stores the deviatoric part of the stress vector �(ti) and �e is the deviatoric part of the total strain increment. With reference to the backward Euler integration scheme the time dependent variables at time ti receive the form �( ) ( ) exp ( ) G t G a t t t a t i i i � � � � � � � ! " " � � � � � � � � � � � � � �� � 1 �� � 1 M , (24) � � �( ) exp ( ) ( )t t a t ti i i� � � � � � � � ! " " � � � � � � � � � � 1 1 � � � � � s 1 M , (25) where s�(ti), � � 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 �t t ti i� � �1, and M is the assumed number of Mawell units in the chain model. The stress shift factor is given by a t t ti eq i eq i � � � � � ( ) ( ) sinh ( ) � � � �� ! "" 0 0 , (26) where the equivalent stress �eq follows from �eq i i it t t( ) ( ) ( )� �1 2 1s sT Q , (27) and Q � � �� � �� diag 1 1 1 1 2 1 2 1 2 . (28) Clearly, the backward Euler step makes all variables non- linearly dependent on the stress values found at time ti. 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 � �r � T G A Tas T t t teq i i i� � �� ( ) ( ) ( ) 1 2 1s sT Q , (29) G G t G a t t t a t i i i � � � � � � � � ! " " � � � � � �( ) ( ) exp ( ) � � � � � � �� � 1 � � �� � 1 M , (30) A a t t ti eq i eq i � � � � �� ! "" � � � � � ( ) ( ) sinh ( ) 0 0 , (31) with the primary variables � �a � � �eq i i it G t a t( ) �( ) ( ) T . (32) Note that the current increment of the vector ��(ti), which appears in Eq. (23), is considered as a secondary variable. Then, under the condition that � e is constant, the New- ton-Raphson iterative scheme reads a a rk i k i k kt t� �� �1 1( ) ( ) H , (33) where Jacobian matrix H is given by H � d d r a . (34) The initial values of primary variables at time ti for k � 0 are set to forward Euler estimates. More details about the nu- merical 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- © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 205 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 differ- ent constant strain rates, where the solid lines were obtained experimentally, while the others follow from the numerical analysis. 5 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 esti- mates 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 ele- ment 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 corre- spond to element averages. Fig. 11 indeed confirms large non-homogeneity in the dis- tribution 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 heteroge- neity of the material system and is consequently due to the required compatibility of local and overall fields. This condi- tion 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 206 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No.5–6/2004 Fig 9: Experiments vs. numerical simulations Fig. 10: A simple finite element mesh with groups of elements Fig. 11: Element averages of matrix strains for �E12 4 110� � �s Fig. 12: Element averages of matrix stresses for �E12 4 110� � �s 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 ma- trix. 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, ex- tended by analogy to more complex periodic unit cells such as those displayed in Fig. 7. To illustrate the governing mecha- nism responsible for this fact, the distribution of local equiva- lent 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 E12 � 0,1. For all the studied microstructures, the final deformation pattern has the character of a localized shear layer that is re- sponsible 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 explana- tion comes from the basic assumption of the selected material model. Recall that the model draws on the nonlinear visco- elastic response to be limited purely to the deviatoric part of © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 207 Acta Polytechnica Vol. 44 No.5–6/2004 Fig. 13: Distribution of local equivalent strain fields for discre- tization with 6-noded elements: (a) Hexagonal array, (b) 5-fiber PUC, (c) 10-fiber PUC Fig. 14: Distribution of local equivalent strain fields for discre- tization with 3-noded elements: (a) Hexagonal array, (b) 5-fiber PUC, (c) 10-fiber PUC Fig. 15: Overall response of optimized unit cells – effect of the or- der of discretization 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 conse- quently 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. 6 Conclusions The paper reviewed essential steps of first order homoge- nization procedure implemented in the framework of the finite element method. A family of so-called statistical equiva- lent unit cells was introduced to deal with random microstruc- tures. Unlike the elastic analysis, where a 5-fiber unit cell was found sufficient to arrive at reliable estimates of homo- genized 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 dis- played 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 micro- structure. Clearly, when refining the finite element mesh such a shear band will be captured more accurately but will not de- pend 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 encoun- tered 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 in- elastic phases. As already suggested, the existence of regions of highly localized deformation attributed heterogeneity of a microstructure may considerably influence the overall re- sponse (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 princi- ples. 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 re- fined geometry are usually used in place of laboratory mea- surements. An example of combining such an approach with a two-point averaging scheme has been presented in [14]. Acknowledgments Financial support for this work was provided by GA R grant No. 103/01/D052. References [1] Fish J., Shek K., Pandheeradi M., Shephard M. S.: “Computational plasticity for composite structures based on mathematical homogenization: Theory and practice.” Computer Methods in Applied Mechanics and En- gineering, Vol. 148 (1997) No. 1–2, p. 53–73. [2] Fish J., Shek K.: “Multiscale analysis of large-scale non- linear structures and materials.” International Journal for Computational Civil and Structural Engineering, Vol. 1 (2000), No. 1, p. 79–90. [3] Kouznetsova V. G., Geers M. G. D., Brekelmans W. A. M.: “Multi-scale second-order computational homogeniza- tion of multi-phase materials: A nested finite element so- lution strategy.” Computer Methods in Applied Mechanics and Engineering, (2004), accepted for publication. [4] Massart T. J.: “Multi-scale modeling of damage in ma- sonry structures”, Ph.D. thesis, Technische Universiteit Eindhoven, 2003. [5] Kabele P.: “Assesment of Structural Performance of En- gineered Cementitious Composites by Computer Simu- lation.” CTU Reports, Vol. 5 (2001), No. 4. [6] Hashin Z., Shtrikman S.: “On some variational prin- ciples in anisotropic and nonhomogeneous elasticity.” Journal of the Mechanics ans Physics of Solids, Vol. 10 (1962), p. 335–342. [7] Willis J. R.: “Variational and related methods for the overall properties of composites.” Advances in Applied Mechanics, Vol. 21 (1981), p. 2–4. [8] Lagoudas D. C., Gavazzi A. C., Nigam H.: “Elastoplas- tic behavior of metal matrix composites based on incremental plasticity and the Mori-Tanaka averag- ing scheme.” Computational Mechanics, Vol. 8 (1991), p. 193–203. [9] Dvorak G. J.: “Transformation field analysis of inelastic composite materials.” Proceedings of the Royal Society of London Series A - Mathematical, Physical and Engineering Sciences, Vol. 437 (1992), p. 311–327. [10] Dvorak G. J., Bahei-El-Din Y. A., Wafa A. M.: “Im- plementation of the transformation field analysis for inelastic composite-materials.” Computational Mechanics, Vol. 14 (1994), No. 3, p. 201–228. [11] Castaneda P. P., Suquet P.: ”Nonlinear composites.” Ad- vances in Applied Mechanics, Vol. 34, No. 998, p. 171–302. [12] Masson R., Bornert M., Suquet P., Zaoui A.: “An affine formulation for the prediction of the effective proper- ties of nonlinear composites and polycrystals.” Journal of the Mechanics and Physics of Solids, Vol. 48 (2000), p. 1203–1227. [13] Chaboche J. L., Kruch S, Maire J. F., Pottier T.: “To- wards a micromechanics based inelastic and damage modeling of composites.” International Journal of Plastic- ity, Vol. 17 (2001), No. 4, p. 411–439. [14] Šejnoha M., Valenta R., Zeman J.: “Nonlinear visco- elastic analysis of statistically homogeneous random composites.” International Journal of Multiscale Computa- tional Engineering, (2004), submitted for publication. [15] Šejnoha M., Zeman J.: “Overall viscoelastic response of random fibrous composites with statistically quasi uni- form distribution of reinforcements.” Computer Methods 208 © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ Acta Polytechnica Vol. 44 No.5–6/2004 in Applied Mechanics and Engineering, Vol. 191 (2002), No. 44, p. 5027–5044. [16] Šejnoha M., Zeman J.: “Micromechanical analysis of random composites.” CTU Reports, 6 (1), Czech Tech- nical University in Prague, 2002. [17] Zeman J., Šejnoha M.: “Numerical evaluation of effec- tive properties of graphite fiber tow impregnated by polymer matrix.” Journal of the Mechanics and Physics of Solids, Vol. 49 (2001), No. 1, p. 69–90. [18] Zeman J., Šejnoha M.: “Homogenization of plain weave composites with imperfect microstructure: Part I - Theo- retical formulation.” International Journal of Solids and Structures, (2004), in print. [19] Matouš K., Lepš M., Zeman J. Šejnoha M.: “Apply- ing genetic algorithms to selected topics commonly en- countered in engineering practice.” Computer Methods in Applied Mechanics and Engineering, Vol. 190 (2000), No. 13–14, p. 1629–1650. [20] Shan Z., Gokhale A. M.: “Representative volume ele- ment for non-uniform microstructures.” Computational Materials Science, Vol. 24 (2002), p. 361–379. [21] Beran M. J.: “Statistical Continuum Theories.” Mono- graphs in Statistical Physics, Interscience Publishers, 1968. [22] Torquato S.: “Random heterogeneous materials: Micro- structure and macroscopic properties.” Springer-Verlag, 2002. [23] Michel J. C., Moulinec H., Suquet P.: “Effective proper- ties of composite materials with periodic microstructure: A computational approach.” Computer Methods in Applied Mechanics and Engineering, Vol. 172 (1999), p. 109–143. [24] Bittnar Z., Šejnoha J.: “Numerical methods in structural engineering.” ASCE Press, 1996. [25] Procházka P., Šejnoha J.: “A BEM formulation for homogenization of composites with randomly distrib- uted fibers.” Engineering Analysis with Boundary Elements, Vol. 27 (2) (2003), p. 137–144. [26] Leonov A. I.: “Non-equilibrium thermodynamics and rheology of viscoelastic polymer media.” Rheologica Acta, Vol. 15, (1976), p. 85–98. [27] Valenta R.: „Numerické modelování polymerů“ [Nu- merical modeling of polymers], Master thesis, Czech Technical University in Prague, (2003) (in Czech). [28] Valenta R.: “Numerical implementation of the Leonov model”, In Engineering Mechanics 2004 (I. Zolotarev – A. Poživilová, editors), Institute of Thermomechanics, Czech Academy of Science, (2004), p. 303–304. [29] Valenta R., Šejnoha M.: ”Epoxy resin as a bonding agent in polymer matrix composites: material properties and numerical implementation.” To be presented at ICCES’04 Madeira, Portugal, July 26–29, 2004. Ing. Jan Zeman, Ph.D. phone: +420 224 354 482 Fax: +420 224 310 775 e-mail: zemanj@cml.fsv.cvut.cz Ing. Richard Valenta phone: +420 224 354 472 Fax: +420 224 310 775 e-mail: richard.valenta@fsv.cvut.cz Doc. Ing. Michal Šejnoha, Ph.D. phone: +420 224 354 494 Fax: +420 224 310 775 e-mail: sejnom@fsv.cvut.cz Department of Structural Mechanics Czech Technical University in Prague Faculty of Civil Engineering Thákurova 7 166 29 Prague 6, Czech Republic © Czech Technical University Publishing House http://ctn.cvut.cz/ap/ 209 Acta Polytechnica Vol. 44 No.5–6/2004