DOI:10.14311/APP.2021.30.0114 Acta Polytechnica CTU Proceedings 30:114–120, 2021 © Czech Technical University in Prague, 2021 available online at http://ojs.cvut.cz/ojs/index.php/app APPLICATION OF THE MORI-TANAKA METHOD TO DESCRIBE THE RATE-DEPENDENT BEHAVIOR OF UNIDIRECTIONAL FIBROUS COMPOSITES Soňa Valentováa, Michal Šejnohaa, ∗, Jan Vorela, Radek Sedláčekb, Pavel Padevěta a Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague 6, Czech Republic b Czech Technical University in Prague, Division of biomechanics, Faculty of Mechanical Engineering, Technická 4, 160 00 Prague, Czech Republic ∗ corresponding author: sejnom@fsv.cvut.cz Abstract. 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. Keywords: Fibrous composites, generalized Leonov model, homogenization, Maxwell chain model, Mori-Tanaka method, periodic unit cell, viscoelasticity. 1. Introduction Unidirectional fiber-reinforced composites have al- ready found a way into a variety of industrial fields. For such reasons as low weight, good thermal proper- ties, high corrosion and chemical resistance, high lon- gitudinal stiffness and strength, lower economic and environmental cost compared to other widely used materials, their application appears in aerospace, au- tomotive, 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 char- acter of the matrix phase is necessary. Despite a random nature of the distribution of fibers in the matrix, easily visible from images of ma- terial cross sections, a periodic hexagonal arrange- ment (PHA) model [1] is adopted, because for high volume fraction of fibers it has been suggested as suf- ficiently accurate representation of a statistically ho- mogeneous microstructure [2]. Moreover, when lim- iting attention to elasticity the Mori-Tanaka (MT) method [3] delivers the macroscopic response compa- rable to finite element (FEM) simulations assuming the PHA model [4]. Unfortunately, when moving be- yond elasticity the original Mori-Tanaka scheme ap- pears 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 experimen- tal program to calibrate the Leonov model [7, 8, 10] can be found in [11]. The paper is organized as follows. Section 2 intro- duces 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 dif- ferent 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 Sec- tion 5. 2. Examined composite system Following our recent study on the Mori-Tanaka method [5] we again consider a unidirectional fibrous 114 http://dx.doi.org/10.14311/APP.2021.30.0114 http://ojs.cvut.cz/ojs/index.php/app vol. 30/2021 Viscoelastic modeling of unidirectional fibrous composites Figure 1. Cross section of examined fibrous composite. (a) (b) Figure 2. a) MTS Alliance 30kN, b) polymer dog bone specimens. EA ET GA GT νA cf [GPa] [GPa] [GPa] [GPa] [-] [-] 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 representa- tive cross-section of the yarn and its binary form is displayed in the Fig. 1. An image analysis was car- ried 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 elas- ticity, the PHA model may yield the overall behav- ior orthotropic [5]. Thus the application of the MT methods calls for caution when applied to real mate- rial systems. 3. Generalized Leonov Model The generalized Leonov model is adopted here to de- scribe the nonlinear viscoelastic behavior of the epoxy matrix. The formulation assumes the bulk response be linearly elastic σm = Kεv . (1) The viscoelastic response is thus driven solely by the deviatoric stress and strain components. For the Maxwell chain model the corresponding constitutive equations read dsij dt = M! µ=1 2Gµ( deij dt − de p,µ ij dt ), (2) sij = M! µ=1 s µ ij . (3) 115 S. Valentová, M. Šejnoha, J. Vorel et al. Acta Polytechnica CTU Proceedings Figure 3. Strain rate dependent uniaxial tensile tests. Figure 4. Eyring plot. The creep strain increment follows from the Eyring flow model as dep dt = 1 2A sinh(τ /τ0) , (4) where A and τ0 are model parameters determined ex- perimentally from either uniaxial [10, 12] or torsional tests [8] performed at different strain rates. For ex- ample, rewriting Eq. (4) for a uniaxial tension we get σy = τ0 √ 3 ln(2A √ 3) + τ0 √ 3 ln ε, (5) which upon substituting the experimentally observed yield, or the maximum, tensile stress σy allows us to determine the necessary material parameters, see Figures 3 and 4 and the next section. 3.1. 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 ten- sion in the displacement controlled regime at a spe- cific 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 con- firm 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). Figure 5. Creep experiments at various stress levels. Figure 6. Approximation of master curve. 3.2. 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 J (t) = M! µ=1 Jµ " 1 − exp # − t τµaσ (t) $% , (6) R(t) = M! µ=1 Eµ exp # − t θµaσ (t) $ . (7) where τµ and θµ are the selected retardation and relaxation times, respectively, and aσ is the stress shift factor to define a nonlinear viscosity of a dash- pot element in, e.g. the Kelvin unit ηµ(τeq ) = τµJµaσ (τeq ), τeq = & 1/2sij sij . The compliances Jµ can be derived by minimizing the least square dif- ference of J (t) given by Eq. (6) and experimentally constructed master curve for the selected set of retar- dation times, while exploiting the time-stress super- position principle. Note that τ1 should be sufficiently small to represent in the limit τ → 0 a solid mate- rial. 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 are 116 vol. 30/2021 Viscoelastic modeling of unidirectional fibrous composites (a) (b) Figure 7. a) Periodic unit cell (PUC) with hexagonal arrangement of fibers in yarns, b) mesh of PUC. µ τµ [s] Jµ [MPa−1] θµ [MPa·s] Eµ [MPa] 1 0.001 2.606512×10−4 9.927397×10−3 2.787166×101 2 0.01 1.905071×10−6 9.966502×10−2 1.278184×101 3 0.1 8.808431×10−7 9.815126×10−1 7.056602×101 4 1 4.934025×10−6 9.543319×10+0 1.711529×102 5 10 1.276165×10−5 9.344254×10+1 2.334448×101 6 100 1.969419×10−5 9.580883×10+2 1.418353×102 7 1000 1.290521×10−5 8.275395×10+3 5.659977×102 8 10000 6.291266×10−5 9.647045×10+4 1.586346×102 9 100000 7.887707×10−6 2.005373×10+5 1.944645×103 10 1000000 1.577867×10−3 4.168654×10+5 5.096147×102 Table 2. Parameters of Maxwell chain model. 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 Ta- ble 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]. 4. 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]. 4.1. 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 ∆σf = Lf ∆εf , ∆σm = 'Lm(∆εm − ∆µm), (8) where Lr , r = f, m1, is the phase material stiffness matrix and 'Lm represents the dependence on the cur- rent viscoelastic modulus. The local strain incre- ments written in terms of the prescribed macroscopic strain increment ∆E and the increment of the creep strain ∆µm developed in the matrix follow from Dvo- rak’s transformation field analysis, see [14], in the form ∆εf = 'Af ∆E + 'Df m∆µm, (9) ∆εm = 'Am∆E + 'Dmm∆µm. (10) where 'Ar and 'Drm are the mechanical strain localiza- tion factors and strain transformation influence func- tions, respectively. In the light of the MT method [3] they can be expressed as 'Am = ( cmI + cf 'Tf )−1 , 'Af = 'Tf 'Am, (11) 'Drm = * I − 'Ar + * 'Lm − Lf +−1 'Lm. (12) The partial strain concentration factor 'Tf 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 evi- dent from Fig. 8 comparing the FEM simulations with 1Subscripts f, m stand for the fiber and matrix phase, re- spectively. 117 S. Valentová, M. Šejnoha, J. Vorel et al. Acta Polytechnica CTU Proceedings Figure 8. Stress-strain diagrams for various uniaxial tensile strain rate loading. 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. 4.2. 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 ∆σf = (1 − ω)Lf ∆εf . (13) Two particular formulations describing the evolution of ω are investigated. The one suggested in [12] reads ω = 1 − " τ teq N τ0 /sinh( τ teq N τ0 ) %M , (14) where M, N are the model parameters and the equiv- alent stress τeq is given by τeq = , 1 2 sij sij . (15) The influence of parameters N and M is illustrated in Fig. 9 suggesting a rather similar effect of the two pa- rameters 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 strain Ėxy = 1 × 10−4 s−1 shows that such a formulation does not provide satisfactory results. Therefore, we propose a new formulation written as ω = N " 1 − # τ teq − T τ0 /sinh( τ teq − T τ0 ) $%M , (16) (a) (b) Figure 9. Influence of parameters M, N in describ- ing the evolution of ω according to Eq. (14): (a) influ- ence of parameter M , N =1, (b) influence of parame- ter N , M =1. Figure 10. Comparison of FEM and MT predictions with ω given by Eq. (14). 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 118 vol. 30/2021 Viscoelastic modeling of unidirectional fibrous composites (a) (b) (c) Figure 11. Influence of parameters M, N in describ- ing the evolution of ω according to Eq. (16): (a) in- fluence of parameter M , N =1, T =0; (b) influence of parameter N , M =1, T =0; (c) influence of parameter T , M =N =1. 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 N = 〈0, 1〉, (17) and if τeq < T then ω = 0, (18) Figure 12. Comparison of FEM and MT predictions with ω given by Eq. (16). if τeq > T then ω follows from Eq. (16). (19) Figure 12, showing again the macroscopic stress- strain 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). 5. Conclusions The main goal of this study was to confirm suitabil- ity of the Mori-Tanaka method to replace demand- ing finite element simulations in the prediction of the macroscopic response of fibrous composites. This would prove advantageous particularly in a multi- scale analysis of textile composites where the re- sponse at the level of yarn must be solved indepen- dently from the textile ply. To that end, two modifi- cations 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 consti- tutive 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 re- sponse of two-phase composites made of aligned fibers embedded into a polymer matrix. Acknowledgements The financial support provided by the GAČR grant No. 19-15666S and by the Czech Technical University in 119 S. Valentová, M. Šejnoha, J. Vorel et al. Acta Polytechnica CTU Proceedings Prague within SGS project with the application registered under the No. SGS20/038/OHK1/1T/11 is gratefully ac- knowledged. References [1] J. L. Teplý, G. J. Dvorak. Bound on overall instantaneous properties of elastic-plastic composites. Journal of the Mechanics and Physics of Solids 36(1):29–58, 1988. [2] J. Zeman, M. Šejnoha. Numerical evaluation of effective properties of graphite fiber tow impregnated by polymer matrix. Journal of the Mechanics and Physics of Solids 49(1):69–90, 2001. [3] Y. Benveniste. A new approach to the application of Mori-Tanaka theory in composite materials. Mechanics of Materials 6:147–157, 1987. [4] M. Šejnoha, J. Zeman. Micromechanics in Practice. WIT Press, Southampton, Boston, 2013. [5] S. Valentová, M. Šejnoha, J. Vorel. Comparing Mori- Tanaka method and first-order homogenization scheme in the viscoelastic modeling of unidirectional fibrous composites. Acta Polytechnica CTU Proceedings 26:133–138, 2020. doi:10.14311/APP.2020.26.0133. [6] R. Masson, M. Bornert, P. Suquet, A. Zaoui. An affine formulation for the prediction of the effective properties of nonlinear composites and polycrystals. Journal of the Mechanics and Physics of Solids 48:1203–1227, 2000. [7] M. Šejnoha, R. Valenta, J. Zeman. Nonlinear viscoelastic analysis of statistically homogeneous random composites. International Journal for Multiscale Computational Engineering 2(4):645–673, 2004. [8] R. Valenta, M. Šejnoha, J. Zeman. Macroscopic constitutive law for mastic asphalt mixtures from multiscale modeling. International Journal for Multiscale Computational Engineering 8(1):131–149, 2010. [9] J. Vorel, E. Grippon, M. Šejnoha. Effective thermoelastic properties of polysiloxane matrix based plain weave textile composites. International Journal for Multiscale Computational Engineering 13(3):181–200, 2015. [10] T. A. Tervoort. Constitutive modeling of polymer glasses: Finite, nonlinear visocelastic behaviour of polycarbonate. Ph.D. thesis, Eindhoven University of Technology, Eindhoven, 1996. [11] M. Šejnoha, S. Valentová, J. Vorel, et al. Modeling of viscoelastic response of unidirectional fibrous composites made of basalt and carbon fibers. In W. P. de Wilde, S. Hernandes (eds.), 10th International Conference on High Performance and Optimum Design of Structures and Materials. WIT Press eLibrary, https://www.witpress.com/elibrary, 2020. [12] R. Valenta, M. Šejnoha. Hierarchical modeling of mastic asphalt in layered road structures based on the Mori-Tanaka method. Acta Polytechnica 52(6):48–58, 2012. [13] S. Valentová, V. Hrbek, J. Vorel, M. Šejnoha. Strength of composite yarn under biaxial loading. Acta Polytechnica CTU Proceedings 15:131–136, 2018. [14] G. J. Dvorak, Y. Benveniste. On transformation strains and uniform fields in multiphase elastic media. Proceedings of the Royal Society of London Series A - Mathematical, Physical and Engineering Sciences 437(1900):291–310, 1992. 120 https://doi.org/10.14311/APP.2020.26.0133