APPLICATION OF THE MORI-TANAKA METHOD TO DESCRIBE THE RATE-DEPENDENT BEHAVIOR OF UNIDIRECTIONAL FIBROUS COMPOSITES

This paper examines the possibility of using the Mori-Tanaka micromechanical model describe the rate dependent behavior of the polymer matrix based fibrous composites. The generalized Leonov model is adopted to capture the time and rate dependent character of the selected matrix, while fibers are assumed elastic. The performance of the Mori-Tanaka method is tested against the finite element simulations carried out in the framework of first-order homogenization. For simplicity, the periodic hexagonal array model is chosen to represent the microstructural arrangement of fibers in the yarn cross-section. To match the predictions provided by the two approaches a suitable modification to the original Mori-Tanaka method is proposed. An extensive parametric study is presented to illustrate a considerable improvement of the predictive capability of the modified Mori-Tanaka method.


Introduction
Unidirectional fiber-reinforced composites have already found a way into a variety of industrial fields. For such reasons as low weight, good thermal properties, high corrosion and chemical resistance, high longitudinal stiffness and strength, lower economic and environmental cost compared to other widely used materials, their application appears in aerospace, automotive, marine, and various construction sectors. In order to take full advantage of the benefits the composite material offers, a suitable computational model predicting its macroscopic response including the rate dependent behavior due to the viscous character of the matrix phase is necessary.
Despite a random nature of the distribution of fibers in the matrix, easily visible from images of material cross sections, a periodic hexagonal arrangement (PHA) model [1] is adopted, because for high volume fraction of fibers it has been suggested as sufficiently accurate representation of a statistically homogeneous microstructure [2]. Moreover, when limiting attention to elasticity the Mori-Tanaka (MT) method [3] delivers the macroscopic response comparable to finite element (FEM) simulations assuming the PHA model [4]. Unfortunately, when moving beyond elasticity the original Mori-Tanaka scheme appears too stiff [5] and some action is needed to bring the predictions closer to those provided by FEM [6,7]. In the present study, we build upon the approach proposed in [8] for asphalts. However, their original formulations requires some modification to suit the present material system.
While attention is primarily accorded to testing the Mori-Tanaka method as a very efficient substitute for a highly expensive finite element method in large scale simulation of textile composites [9], the material model representing the epoxy matrix is mentioned only briefly. Further details on extensive experimental program to calibrate the Leonov model [7,8,10] can be found in [11].
The paper is organized as follows. Section 2 introduces the material system examined herein. Next, section 3 shortly describes the Leonov viscoelastic model adopted for the epoxy matrix. For the sake of convenience we briefly summarize the main results of the experimental program needed to calibrate the constitutive model. In particular, creep tests at different stress levels and rate dependent uniaxial tensile strain controlled experiments are described. Finally, Section 4 outlines two types of modifications of the Mori-Tanaka method accompanied by an extensive parametric study to examine their potential. The principal achievements are then summarized in Section 5.

Examined composite system
Following our recent study on the Mori-Tanaka method [5] we again consider a unidirectional fibrous vol. 30/2021 Viscoelastic modeling of unidirectional fibrous composites [-] [-] Carbon fibers 294 13 12 5 0.24 0.56 Table 1. Material properties of carbon fibers and its volume fraction in yarn.
composite system made of elastic carbon fibers and nonlinear viscoelastic polymer matrix. A representative cross-section of the yarn and its binary form is displayed in the Fig. 1. An image analysis was carried out to acquire the volume fraction of the fibers, see Table 1 together with material properties taken from [9].
Because fibers are transversely isotropic and the matrix is isotropic, the overall macroscopic response is expected to be also transversely isotropic. This is precisely what both the MT method and the PHA model would predict in case of elastic response. While the MT method keeps this property even beyond elasticity, the PHA model may yield the overall behavior orthotropic [5]. Thus the application of the MT methods calls for caution when applied to real material systems.

Generalized Leonov Model
The generalized Leonov model is adopted here to describe the nonlinear viscoelastic behavior of the epoxy matrix. The formulation assumes the bulk response be linearly elastic The viscoelastic response is thus driven solely by the deviatoric stress and strain components. For the Maxwell chain model the corresponding constitutive equations read  The creep strain increment follows from the Eyring flow model as where A and τ 0 are model parameters determined experimentally from either uniaxial [10,12] or torsional tests [8] performed at different strain rates. For example, rewriting Eq. (4) for a uniaxial tension we get which upon substituting the experimentally observed yield, or the maximum, tensile stress σ y allows us to determine the necessary material parameters, see

Strain rate dependent tensile tests
Six dog bone specimens in Fig. 2(b) made of 285/500 aero Havel epoxy resin were loaded in simple tension in the displacement controlled regime at a specific strain rate until failure using the MTS Alliance 30kN electromechanical testing machine equipped with 30kN load cell, see Fig. 2(a).
The resulting stress strain diagrams in Fig. 3 confirm the rate dependent response of the epoxy matrix. Replotting the yield stress as a function of the applied strain rate in a logarithmic scale and applying simple linear regression gives the material parameters of the Eyring flow model in Eq. (5).

Creep tests
The Maxwell chain model requires calibration of the creep relaxation function typically found using the Laplace transform of the creep compliance function. One-dimensional format of these functions can be written as where τ µ and θ µ are the selected retardation and relaxation times, respectively, and a σ is the stress shift factor to define a nonlinear viscosity of a dashpot element in, e.g. the Kelvin unit η µ (τ eq ) = τ µ J µ a σ (τ eq ), τ eq = & 1/2s ij s ij . The compliances J µ can be derived by minimizing the least square difference of J (t) given by Eq. (6) and experimentally constructed master curve for the selected set of retardation times, while exploiting the time-stress superposition principle. Note that τ 1 should be sufficiently small to represent in the limit τ → 0 a solid material. Theoretically, τ 1 = 0.001 might not satisfy this requirement but the resulting series validated in [11] proved satisfactory.
To that end, a series of creep tests at different stress levels is performed first. The creep curves in Fig. 5   then horizontally shifted to get the master curve for the selected reference stress. The experimental points as well as the resulting polynomial representation are plotted in Fig. 6.
Ten Maxwell units in parallel were considered in the present study to represent the Maxwell chain model. The resulting parameters are listed in Table 2. The shear moduli in Eq. (2) are obtained from E µ for the selected Poisson ratio set equal to 0.39. Further details including validation of the proposed constitutive model are available in [11].

Homogenization
As mentioned in the introductory part, we consider the predictions provided by FEM as a benchmark to test the performance of the MT method. Owing to a space limitation we present, with reference to FEM simulations, the PHA computational model only, see Fig. 7, and refer the reader for details on theoretical grounds of the first-order homogenization to, e.g. [4,13, and references therein].

Mori-Tanaka method: standard formulation
The MT method, however, deserves more attention.
To begin with, we write the increments of local stress averages in individual phases as where L r , r = f, m 1 , is the phase material stiffness matrix and ' L m represents the dependence on the current viscoelastic modulus. The local strain increments written in terms of the prescribed macroscopic strain increment ∆E and the increment of the creep strain ∆µ m developed in the matrix follow from Dvorak's transformation field analysis, see [14], in the form where ' A r and ' D rm are the mechanical strain localization factors and strain transformation influence functions, respectively. In the light of the MT method [3] they can be expressed as The partial strain concentration factor ' T f follows from the Eshelby transformation field analysis and in general it is a function the shape and orientation of the fiber and instantaneous properties of the matrix.
The principal drawback of the MT method is evident from Fig. 8 comparing the FEM simulations with the MT results for various strain rates of the applied in-plane shear strain. It is reasonable to attribute the observed stiff response of the MT method to its inability to represent localization of inelastic strain since dealing with phase volume averages only.

Mori-Tanaka method: modified formulation
Owing to the previous results we propose a simple modification to match the predictions provided by both MT method and FEM. Being inspired by [8,12] we attempt to describe the stress redistribution due to formation of shear bands by the reduction of stress taken by the fibers through a damage like parameter ω and write the fiber stress increment as Two particular formulations describing the evolution of ω are investigated. The one suggested in [12] reads where M, N are the model parameters and the equivalent stress τ eq is given by The influence of parameters N and M is illustrated in Fig. 9 suggesting a rather similar effect of the two parameters both in terms of stiffness and overall stress reduction. Thus tuning only one of them appears sufficient in adjusting the macroscopic response. However, as seen in Fig. 10, plotting the FEM and MT estimates of the macroscopic response for the applied rate of macroscopic in-plane shear straiṅ E xy = 1 × 10 −4 s −1 shows that such a formulation does not provide satisfactory results. Therefore, we propose a new formulation written as  where M, N, T are again the model parameters, which, unlike to Eq. (14), provide more freedom in the model calibration. The way they control the evolution of ω is evident from Fig. 11. While parameter M adjusts the rate of evolution of ω with τ eq , the other two parameters N, T merely shift the damage curve horizontally and vertically. Thus N controls the minimum stress the fibers may take and T serves to delay the onset of damage evolution. Note that these two parameters are subject to some limitations, particularly and if τ eq < T then ω = 0, if τ eq > T then ω follows from Eq. (16). (19) Figure 12, showing again the macroscopic stressstrain curve for the applied shear strain rateĖ xy = 1 × 10 −4 , promotes this new formulation as almost a perfect match between FEM and MT predictions can be obtained with properly adjusted parameters M, N, T in Eq. (16).

Conclusions
The main goal of this study was to confirm suitability of the Mori-Tanaka method to replace demanding finite element simulations in the prediction of the macroscopic response of fibrous composites. This would prove advantageous particularly in a multiscale analysis of textile composites where the response at the level of yarn must be solved independently from the textile ply. To that end, two modifications to the original MT formulation were tested with particular attention devoted to the nonlinear viscoelastic response of the matrix phase governed by the generalized Leonov model. Both modifications attempt to reduce the stress taken by the fibers by introducing a damage like parameter ω in the constitutive equation for the fiber phase.
The first modification originated from the work of Valenta et al. [8,12]. While successful for asphalt mixtures it proved insufficient in the present study. Therefore, a new modification was proposed to allow for more freedom in the calibration of the evolution law for damage parameter ω. The results clearly show a superiority of this new formulation allowing us to receive a very good agreement between FEM and MT predictions. This also promotes the MT method as a perfect tool for the prediction of macroscopic response of two-phase composites made of aligned fibers embedded into a polymer matrix.
Prague within SGS project with the application registered under the No. SGS20/038/OHK1/1T/11 is gratefully acknowledged.