Acta Polytechnica CTU Proceedings doi:10.14311/APP.2020.26.0133 Acta Polytechnica CTU Proceedings 26:133–138, 2020 © Czech Technical University in Prague, 2020 available online at https://ojs.cvut.cz/ojs/index.php/app COMPARING MORI-TANAKA METHOD AND FIRST-ORDER HOMOGENIZATION SCHEME IN THE VISCOELASTIC MODELING OF UNIDIRECTIONAL FIBROUS COMPOSITES Soňa Valentová, Michal Šejnoha∗, Jan Vorel Czech Technical University in Prague, Faculty of Civil Engineering, Department of Mechanics, Thákurova 7, 166 29 Prague 6, Czech Republic ∗ corresponding author: sejnom@fsv.cvut.cz Abstract. A comparative study of the viscous response of polymer matrix based fibrous composites predicted by the Mori-Tanaka method and finite element simulations based on the 1st order homoge- nization theory is presented. Aligned basalt and carbon fibers embedded into a polymeric matrix are considered to represent a quasi isotropic and transversely isotropic two-phase systems. While differences in the prediction of the macroscopic elastic response are attributed merely to the properties of the fiber phase, the viscoelastic behavior is largely affected by the selected homogenization method. A stiffer response predicted by the Mori-Tanaka method for both creep and relaxation tests is observed for both material systems and supports similar finding found in the literature. Thus suitable modifications of the original formulation of such two-point averaging schemes are needed for them to be applicable in the multi-scale modeling of generally anisotropic yarns in plane weave textile composites. Keywords: Viscoelasticity, homogenization, Mori-Tanaka method, periodic unit cell, transformation field analysis, fibrous composites. 1. Introduction A plane weave textile composites made of woven yarns reinforced by unidirectional fibers bonded to a poly- meric matrix are used in many engineering applica- tions. Since often exposed to rate dependent or long lasting variable loading the viscous behavior of the matrix phase should be properly represented to arrive at reliable predictions of generally time dependent macroscopic response of such composites. It has been shown in [1] that random nature of the yarn microstructure may play a significant role in the prediction of its macroscopic viscoelastic response particularly in case of high contrast in material prop- erties of the fiber and matrix phase. Therein, the concept of statistically equivalent periodic unit cell (SEPUC) [2, 3] was exploited to address adequacy of the often assumed periodic hexagonal arrangement (PHA) of fibers. This issue is also briefly addressed herein. The point is that if sufficiently accurate, the PHA model may be used to validate the applicability of the Mori-Tanaka method because of their corre- spondence in the elastic regime. Following the footsteps of [4] the virtual creep or relaxation tests carried out numerically at the level of yarn could be used to calibrate a suitable macroscopic viscoelastic model to represent the yarn behavior at the level of textile. However, the macroscopically isotropic response of asphalt mixtures assumed in [4] is hardly acceptable for macroscopically orthotropic fi- brous composites. This opens the way to fully coupled multi-scale analysis. Most often the FE2 computa- tional scheme [5] is used, where the SEPUC model at the level of yarns is loaded by increments of the local strain averages derived at the level of yarns for a given increment of the macroscopic load. The SEPUC analysis provides the instantaneous homogenized yarn stiffness and corresponding stress increment. Such a format of stress update may often be very compu- tationally expensive. A potential route appears in replacing detailed finite element simulations by a suit- able micromechanical model such as the Mori-Tanaka averaging scheme [6] assuming a piece-wise uniform distribution of stresses and strains in individual ma- terial phases thus utilizing the so called two-point averaging scheme in case of two-phase materials, see e.g. [7, 8]. However, using the Mori-Tanaka method may raise some doubts because of the method limitations in connection to time dependent or nonlinear analyses. In particular, much stiffer response delivered by the Mori-Tanaka method in comparison to unit cell finite element simulations has often been reported. The pur- pose of this study is to contribute to this subject in the light of purely viscoelastic behavior of the matrix phase and if needed to offer suitable modifications to the original formulation of such micromechanical mod- els while still keeping their computational efficiency. The paper is organized as follows. Section 2 intro- duces the two material systems examined herein. A formulation of the adopted viscoelastic model is pre- sented next in Section 3 followed in Section 4 by short theoretical description of the two homogenization ap- proaches. The principal achievements are summarized in Sections 5 and 6. 133 https://doi.org/10.14311/APP.2020.26.0133 https://ojs.cvut.cz/ojs/index.php/app S. Valentová, M. Šejnoha, J. Vorel Acta Polytechnica CTU Proceedings EA ET GA GT νA a b n cf [GPa] [GPa] [GPa] [GPa] [-] [GPa]−1 [GPa]−1 [-] [-] Carbon fiber 294 13 12 5 0.24 0.67 Basalt fiber 69.68 64.82 28.10 26.14 0.40 0.69 Epoxy matrix 2.03 2.03 0.725 0.725 0.40 0.0474 0.00214 0.3526 Table 1. Material properties of individual phases. 2. Materials Two material systems are examined in the present study to address both the influence of microstructural details and material properties of individual phases. The same epoxy resin was considered for both carbon and basalt fiber composite with their microstructure displayed in Fig. 1(a) and (b), respectively. (a) (b) Figure 1. Example images of local microstructure: a) carbon fibers, b) basalt fibers. Reprinted from [9] with permission from Begell House Publishers. The corresponding material data taken from [9– 11] are listed in Table 1. While basalt fibers are quasi-isotropic the carbon fiber composite represents a highly anisotropic system. On the contrary, the contrast in material properties of the fiber and matrix in the transverse plane is relatively small in case of carbon fibers as oppose to basalt fibers. This particu- lar mismatch in material properties will be reflected in Section 4.2 when constructing the computational model. 3. Local constitutive equations For the selected epoxy resin the last three parameters in Table 1 are the model parameters of a simple power law like formula for the relaxation rate1 R(t,τ) = 1 a + b(t− τ)n , (1) where R(t,τ) is the relaxation function representing the stress at time t due to a unit strain applied at time τ and held constant. To facilitate the numerical solution this function is typically approximated by the Dirichlet series in the form R(t,τ) = M∑ µ=1 Eµ(τ) exp [ τ θµ − t θµ ] , (2) 1Note that these specific values of parameters a,b,n require the time t − τ be given in minutes to yield the viscoelastic modulus in [MPa]. where θµ are the a priory selected relaxation times. Functions Eµ are usually obtained by fitting the re- laxation function, e.g. Eq. (1), via Eq. (2) using the method of least squares. (a) (b) Figure 2. a) Maxwell chain unit, b) Relaxation function. Assuming for example the generalized Maxwell chain model in Fig. 2(a) composed of ten Maxwell elements, we get the spring moduli Eµ listed in Ta- ble 2. Note that sum of all spring moduli provides the elastic modulus of the epoxy matrix in Table 1. For illustration, Fig. 2 plots the relaxation functions provided by Eqns. (1) and (2). µ θµ [MPa·s] Eµ [MPa] 1 2.949089 3.548983×101 2 9.506234 9.965608×101 3 8.905084×101 2.133588×102 4 7.939963×102 3.599087×102 5 6.874448×103 4.325570×102 6 5.709058×104 3.954614×102 7 5.221879×105 2.450281×102 8 4.561310×106 1.384739×102 9 4.654514×107 5.785885×101 10 4.749622×108 5.854287×101 Table 2. Parameters of Maxwell chain model. For the generalized Maxwell chain model it holds σ = M∑ µ=1 σµ, (3) where σµ are stresses in individual Maxwell elements. Integrating the differential form of constitutive equa- tion for each Maxwell element over the time increment 134 vol. 26/2020 Viscoelastic modeling of unidirectional fibrous composites ∆t and using Eq. (3) gives the increment of the total stress ∆σ as ∆σ = Ê (∆ε− ∆µ) , (4) where the instantaneous Young’s modulus Ê and the increment of viscoelastic strain ∆µ read [12] Ê = N∑ µ=1 ∆tEµ θµ ( 1 − exp [ − ∆t θµ ]) , (5) ∆µ = 1 Ê N∑ µ=1 ( 1 − exp [ − ∆t θµ ]) σµ(ti−1), (6) where σµ(ti−1) represents the hidden stress at the beginning of the new time increment. The stress σµ(ti) at the end of the time increment is then provided by σµ(ti) = σµ(ti−1) exp [ − ∆t θµ ] + (7) + Eµ ( 1 − exp [ − ∆t θµ ]) (∆ε− ∆µ). Extension of the above set of one-dimensional equa- tions to a general three-dimensional space is provided, e.g. in [1]. 4. Homogenization schemes In this section we briefly review the Mori-Tanaka mi- cromechanical model for a two phase material system in the framework of transformation field analysis [13] and the first-order homogenization scheme grounding on the existence of periodic fields. 4.1. Mori-Tanaka method Consider a two-phase composite made of elastic trans- versely isotropic fibers of a circular cross-section em- bedded into a viscoelastic isotropic matrix. Assuming piece-wise constant local fields in individual phases r = f,m2 allows us to write the local constitutive equations in the form ∆σf = Lf ∆εf, ∆σm = L̂m(∆εm − ∆µm), (8) where L is the material stiffness matrix and L̂ rep- resents the dependence on the viscoelastic modulus Ê. Suppose that the composite is loaded on its outer boundary either by the prescribed displacements (strain control) or tractions (stress control) compatible with macroscopic uniform strain ∆E or stress ∆Σ increments. The local strains or stresses in terms of their applied macroscopic counterparts follow from Dvorak’s transformation field analysis [13] as ∆εf = Âf ∆E + D̂fm∆µm, (9) ∆εm = Âm∆E + D̂mm∆µm, (10) ∆σf = B̂f ∆Σ − F̂fmL̂m∆µm, (11) ∆σm = B̂m∆Σ − F̂mmL̂m∆µm, (12) 2Subscripts f,m stand for the fiber and matrix phase, re- spectively. where Âr,B̂r and D̂rm,F̂rm,r = f,m, are the me- chanical strain and stress localization factors and strain and stress transformation influence functions, respectively. The Mori-Tanaka method gives the mechanical lo- calization factors in the form Âm = [ cmI + cf T̂f ]−1 , Âf = T̂f Âm, (13) B̂m = [ cmI + cf Ŵf ]−1 , B̂f = Ŵf B̂m,(14) where the so called partial strain and stress localiza- tion factors T̂f, Ŵf can be expressed in terms of the Eshelby tensor Ŝ as T̂f = [ I − P̂ ( L̂m − L̂f )]−1 , (15) Ŵf = [ I − Q̂ ( M̂m − M̂f )]−1 , (16) P̂ = ŜM̂m, P̂ = −L̂m ( Ŝ − I ) , (17) where M = L−1 is the material compliance matrix. Note that in case of the Mori-Tanaka method the Eshelby tensor Ŝ depends on instantaneous proper- ties of the matrix phase. For a cylindrical fiber3 its particular form is available, e.g. in [8]. For a two-phase composite the transformation in- fluence functions are readily provided in terms of the localization factors as D̂fm = ( I − Âf )( L̂m − Lf )−1 L̂m, (18) D̂mm = ( I − Âm )( L̂m − Lf )−1 L̂m, (19) F̂fm = ( I − B̂f )( M̂m − Mf )−1 M̂m, (20) F̂mm = ( I − B̂m )( M̂m − Mf )−1 M̂m. (21) 4.2. 1st order homogenization Assuming again the composite be loaded either by the macroscopic uniform strain ∆E or stress ∆Σ admits the following decomposition ∆u(x) = ∆E x + ∆u∗(x), (22) ∆�(x) = ∆E + ∆�∗(x), (23) where ∆u∗(x) is the fluctuation part of the local dis- placement field ∆u(x). The former one is considered periodic to give the volume average of the fluctuation strain ∆�∗(x) equal to zero. The stepping stone in the derivation of fluctuation displacements ∆u∗(x) is the Hill lemma given by〈 δ�T∆σ 〉 = δET∆Σ. (24) 3Note that for cylindrical fibers the Eshelby tensor depends on the matrix phase Poisson ratio only, which is assumed con- stant in our study, i.e. Ŝ = S . Also note that keeping the time increment ∆t constant yields the viscoelastic modulus Êm also constant. Thus the localization and transformation matrices need to be evaluated only once. 135 S. Valentová, M. Šejnoha, J. Vorel Acta Polytechnica CTU Proceedings where 〈·〉 stands for volume averaging. In the framework of finite element discretization we approximate ∆u∗(x) in terms of their nodal values ∆r to write ∆u∗(x) = N(x)∆r, ∆�∗(x) = B(x)∆r. (25) Substituting from Eq. (25)2 into (24) yields the final system of algebraic equations in the form [ K11 K12 K21 K22 ]{ ∆E ∆r } = { ∆Σ + ∆F 0 ∆f0 } . (26) Individual matrices listed in Eq. (26) are written as K11 = 1 Ω ∫ Ω L̂(x) dΩ, K12 = KT21 = 1 Ω ∫ Ω L̂(x)B(x) dΩ, (27) K22 = 1 Ω ∫ Ω BT(x)L̂(x)B(x) dΩ, and components of the right-hand side vector are ∆F 0 = 1 Ω ∫ Ω L̂(x)µ(x) dΩ, (28) ∆f0 = 1 Ω ∫ Ω BT(x)L̂(x)∆µ(x) dΩ. Recall that the system of equations 26 is adopted when loading the composite by the prescribed macroscopic stress such as in case of creep. Simulating a relaxation test calls for the macroscopic strain ∆E be prescribed. The system of equations 26 then simplifies to (K22)∆r = −(K21)∆E + ∆f0. (29) The increments of macroscopic strains ∆E (Σ pre- scribed) or stress ∆Σ (E prescribed) follow from the volume averaging of their local counterparts. Further details can be found in [1, 8]. 5. Results The first study aimed at comparing the results pro- vided by various unit cell models. Figures 3(b-c) show 5-fiber and 10-fiber periodic unit cells (PUC). Note that these were constructed by simply chang- ing the diameter and slightly also the position of SEPUCs derived originally in [8, Chapter 3] for a carbon-carbon composite to meet the required vol- ume fractions. Such PUC thus should not be termed SEPUC for the present systems. Inadequacy of such a simple approach is confirmed by the results plotted in Fig. 5 for the basalt fiber composite subjected to both creep and relaxation loading conditions displayed in Fig. 4. A highly anisotropic behavior predicted by 5-fiber and 10-fiber PUC in comparison to PHA model in Fig. 3(a) is evident in Figs. 5(a,c) showing a different (a) (b) (c) Figure 3. Periodic unit cells: a) PHA model, b) basalt 5 fiber UC, c) basalt 10 fiber UC. Time [s] Σ (Σ )xy = 100 [MPa]xx Σ Time [s] E xy xxE (E ) = 0.375 [%] (a) (b) Figure 4. Applied loading: a) creep test, b) relax- ation test. response of the two systems when loading the com- posite along X and Y axes, but also in Figs. 5(b,d) showing an evolution of a non-negligible shear strain and stress for the creep and relaxation loading con- ditions, respectively, associated with the two loading cases. While, as expected, the response predicted by the PHA model confirms a material isotropy in the transverse plane. Further study comparing the predictions provided by the FEM simulations and the Mori-Tanaka thus adopted the PHA model only. The results appear in Fig. 6. Figure 6(a) shows the evolution of the macroscopic normal strain Exx caused by the macroscopic normal stress Σxx. For both material systems the initial elas- tic response predicted by the Mori-Tanaka method appears more compliant in comparison to FEM simula- tions, while the subsequent time dependent evolution of strains confirms much stiffer behavior suggested by the Mori-Tanaka method. This is supported by the results derived for the shear loading, Figure 6(b), as well as for the relaxation tests, Figs. 6(c,d) The results clearly show that the predictions depend not only on the material system but also on the applied type of loading conditions. While the differences in the elastic behavior can be associated with the material properties of the fiber phase, the evolution of local creep strain is driven by the selected computational model. A non-uniform distribution of creep strain in the matrix phase caused by the microstructural details, 136 vol. 26/2020 Viscoelastic modeling of unidirectional fibrous composites 0 60 120 180 240 300 360 Time [s] 0,50 0,75 1,00 1,25 E [ % ] PHA 5 fiber UC 10 fiber UC Σ xx applied, E xx plotted Σ yy applied, E yy plotted (a) 0 60 120 180 240 300 360 Time [s] -0,10 -0,08 -0,06 -0,04 -0,02 0,00 E x y [ % ] PHA 5 fiber UC 10 fiber UC Σ xx applied Σ yy applied (b) 0 60 120 180 240 300 360 Time [s] 40 50 60 70 80 Σ [ M P a ] PHA 5 fiber UC 10 fiber UC E xx applied, Σ xx plotted E yy applied, Σ yy plotted (c) 0 60 120 180 240 300 360 Time [s] 0,0 0,5 1,0 1,5 2,0 2,5 3,0 Σ x y [ M P a ] PHA 5 fiber UC 10 fiber UC E xx applied E yy applied (d) . Figure 5. Test of SEPUC applicability: (a-b) creep test - time variation of macroscopic strains due to pre- scribed macroscopic stresses Σxx, Σyy, (c-d) relaxation test - time variation of macroscopic stresses due to prescribed macroscopic strains Exx, Eyy 0 60 120 180 240 300 360 Time [s] 0,8 1,0 1,2 1,4 1,6 1,8 E x x [ % ] FEM PHA Mori-Tanaka Basalt fiber Carbon fiber (a) 0 60 120 180 240 300 360 Time [s] 2,0 3,0 4,0 5,0 6,0 E x y [ % ] FEM PHA Mori-Tanaka Carbon fiber Basalt fiber (b) 0 60 120 180 240 300 360 Time [s] 25 30 35 40 45 50 55 60 Σ x x [ M P a ] FEM PHA Mori-TanakaBasalt fiber Carbon fiber (c) 0 60 120 180 240 300 360 Time [s] 6 8 10 12 14 Σ x y [ M P a ] FEM PHA Mori-Tanaka Basalt fiber Carbon fiber (d) Figure 6. Comparison of macroscopic response pre- dicted by finite element simulations and Mori-Tanaka method for both basalt carbon fiber composites: (a-b) creep test, (c-d) relaxation test - a) Σxx prescribed, b) Σxy prescribed, c) Exx prescribed, d) Exy pre- scribed. 137 S. Valentová, M. Šejnoha, J. Vorel Acta Polytechnica CTU Proceedings recall Fig. 5, and the mismatch in material properties of the fiber and the matrix phase then renders a significantly more compliant macroscopic response in comparison to the Mori-Tanaka method, which builds on the piecewise constant stress and strain fields in individual phases, recall Eqns.(9) - (12). 6. Conclusions The present paper examined a potential application of the Mori-Tanaka method in prediction of effective viscoelastic response of polymer matrix based fibrous composite. The results presented in Fig. 6 confirmed a relatively stiff response delivered by the Mori-Tanaka method when compared to finite element simulations already observed in some previous studies. Thus for the Mori-Tanaka method presented in the framework of Dvorak’s transformation field analysis to be applicable as a very efficient substitute of com- putationally demanding finite element simulations in multi-scale analysis requires some modifications to its original formulation. With regard to generally nonlin- ear material behavior the literature offers a number of potential approaches. Reduction of stresses in the re- inforcing phase, typically assumed linearly elastic, by introducing an artificial damage parameter is one par- ticular example [7, 14]. Slightly different approach has been proposed in [15] for the analysis of quasi-brittle materials where the effect of strain localization was ad- dressed with the aid of interfacial coating to allow for gradual debonding and thus letting the matrix to take all the stresses. Another approach offers more refined transformation field analysis where the concentration and transformation matrices are derived from finite element method using rather crude meshes where each element represents a single material [16]. All these approaches will be examined in our future studies depending on the results provided by an accompanied experimental program. Acknowledgements The financial support provided by the GAČR grant No. 19- 15666S and by the Czech Technical University in Prague within SGS project with the application registered under the No. SGS19/032/OHK1/1T/11 is gratefully acknowl- edged. References [1] M. Šejnoha, J. Zeman. Overall viscoelastic response of random fibrous composites with statistically quasi uniform distribution of reinforcements. Computer Methods in Applied Mechanics and Engineering 191(44):5027–5044, 2002. doi:10.1016/S0045-7825(02)00433-4. [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. doi:10.1016/S0022-5096(00)00027-2. [3] J. Zeman, M. Šejnoha. From random microstructures to representative volume elements. Modelling and Simulation in Materials Science and Engineering 15(4):S325–S335, 2007. doi:10.1088/0965-0393/15/4/S01. [4] 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. doi:10.1615/IntJMultCompEng.v8.i1.100. [5] V. Kouznetsova, M. G. D. Geers, W. A. M. Brekelmans. Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme. International Journal for Numerical Methods in Engineering 54(8):1235–1260, 2002. doi:10.1002/nme.541. [6] Y. Benveniste. A new approach to the application of Mori-Tanaka theory in composite materials. Mechanics of Materials 6:147–157, 1987. doi:10.1016/0167-6636(87)90005-6. [7] 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. [8] M. Šejnoha, J. Zeman. Micromechanics in Practice. WIT Press, Southampton, Boston, 2013. [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. doi:10.1615/IntJMultCompEng.2014011020. [10] F. Fára. Komplexní hodnocení objemových změn kompozitních soustav. Ph.D. thesis, Czech Technical University in Prague, Faculty of Civil Engineerin, Prague, Czech Republic, 1990. [11] S. Valentová, V. Hrbek, J. Vorel, M. Šejnoha. Strength of composite yarn under biaxial loading. Acta Polytechnica CTU Proceedings 15:131–136, 2018. doi:10.14311/APP.2018.15.0131. [12] Z. Bittnar, J. Šejnoha. Numerical methods in structural engineering. ASCE Press, 1996. [13] 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. doi:10.1098/rspa.1992.0062. [14] 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. doi:10.1615/IntJMultCompEng.v2.i4.80. [15] J. Vorel, V. Šmilauer, Z. Bittnar. Multiscale simulations of concrete mechanical tests. Journal of Computational and Applied Mathematics 236(18):4882–4892, 2012. [16] G. Dvorak, Y. Bahei-El-Din, A. Wafa. The modeling of inelastic composite materials with the transformation field analysis. Modelling and Simulation in Materials Science and Engineering 2(3A):571, 1999. doi:10.1088/2F0965-0393.2F2.2F3a.2F011. 138 https://doi.org/10.1016/S0045-7825(02)00433-4 https://doi.org/10.1016/S0022-5096(00)00027-2 https://doi.org/10.1088/0965-0393/15/4/S01 https://doi.org/10.1615/IntJMultCompEng.v8.i1.100 https://doi.org/10.1002/nme.541 https://doi.org/10.1016/0167-6636(87)90005-6 https://doi.org/10.1615/IntJMultCompEng.2014011020 https://doi.org/10.14311/APP.2018.15.0131 https://doi.org/10.1098/rspa.1992.0062 https://doi.org/10.1615/IntJMultCompEng.v2.i4.80 https://doi.org/10.1088/2F0965-0393.2F2.2F3a.2F011 Acta Polytechnica CTU Proceedings 26:134–139, 2020 1 Introduction 2 Materials 3 Local constitutive equations 4 Homogenization schemes 4.1 Mori-Tanaka method 4.2 1st order homogenization 5 Results 6 Conclusions Acknowledgements References