Mori–Tanaka Formalism-Based Method Used to Estimate the Viscoelastic Parameters of Laminated Composites †

The paper establishes the mechanical properties of a viscoelastic composite material reinforced with fibers, where the fiber is transverse isotropic and the matrix is isotropic (a common case met in engineering practice). A computation method using the Mori–Tanaka mean field method has been developed in order to apply on viscoelastic materials. Using this procedure, the time-dependent response of a viscoelastic composite material can be determined. Schapery’s nonlinear constitutive equation is also used in the compliance matrix determination of the composite material under investigation. Nonlinearity factors were determined by creep tests at different values of stresses and temperatures and for different materials, based on the least squares method. The results obtained experimentally and their comparison with the theoretically obtained values show a good agreement between experiment and calculation.


Introduction
The viscoelastic response of a composite material, consisting of two or more components, is determined by the elastic/viscoelastic properties of the components. The calculation of these properties, based on various mathematical models, using fundamental theorems from the theory of elasticity, variational methods, averaging methods, homogenization theory, etc., was made in a series of papers for biphasic materials [1][2][3][4][5][6]. However, we must mention that there is little work in the literature dedicated to the study of viscoelastic materials subjected to high temperatures. In most of these works, the basic hypothesis considered was that, the constituent materials of the composite are isotropic. In practice, however, there are many cases in which one of the materials has anisotropic properties. In this case, the classical model established for predicting the property of the composite no longer works. However, there are few works in the literature when the reinforcement material has anisotropic or transverse isotropic properties. The problem is, therefore, to establish new relations valid for determining the global properties in these new circumstances. As an example, treated in the paper, a composite will be considered that is reinforced with graphite/epoxy fibers, where the fibers represent a material with anisotropic properties. The reinforcing parallel fibers are considered long enough that the end effects can be neglected.
The calculation method will consider a representative cylindrical element volume (RVE), having the fibers arranged along the axis of the cylinder. The cross section of this cylinder is large enough in Polymers 2020, 12, 2481 2 of 16 relation to the dimensions of the fibers, so that averaging operations can be performed unaffected by the small non-homogeneities in the structure and by the arrangement of the fibers. The fibers will be considered continuous along the entire length RVE [7,8]. It will be hypothesized that the composite is considered transversely isotropic and homogeneous.
In these hypotheses, we use a mathematical model by which the elastic/viscoelastic properties of the material, will be determined by upper and lower bounds of the properties of the components, which are considered to be known.
Using Hill's theory [9][10][11], one may obtain the bounds on bulk modulus, the longitudinal Young's modulus and the Poisson's ratio. These bounds are determined in terms of the volume ratios of fiber v f and matrix v m .
For a reinforced composite with cylindrical fibers, the upper and lower bounds on the effective elastic moduli have been calculated for a geometry of hexagonal fibers array by Hashin [12][13][14][15].
A new direction of research represents the gradual change of the fiber orientation and the study of material at high stress level and temperature. A better model is presented in the following, where exact estimations for some engineering constants are determined.
In the last few years a series of papers studied the viscoelastic creep response of different material subjected to different loads and temperatures.
The creep behavior of linear viscoelastic materials subjected to a combined tension and torsion loading is presented in [16]. The creep response is expressed using fractional-exponential functions. Experimental tests have validated the model for HDPE and organic glass.
All the phases of creep, primary, secondary, tertiary and rupture, are analyzed in [17], based on a nonlinear viscoelastic creep constitutive formula. The parameters of this model are analyzed and evaluated. Other papers studied different aspects of the creep response of different materials [18][19][20][21][22][23][24].

Mathematical Model
To obtain the elastic/viscoelastic coefficients of the biphasic composite, reinforced with long, anisotropic fibers, the Mori-Tanaka theory will be used [25]. For this, RVE as defined in the Introduction section of this paper, will be used, consisting of an epoxy matrix with viscoelastic behavior, reinforced with monotonously aligned and parallel fibers, uniformly distributed within the matrix. Considering the above, the composite obtained will be regarded as orthotropic. If, however, the elliptical cylinders fibers are randomly oriented, the composite body under consideration will be viewed as transverse isotropic. The elliptical fibers are defined by the ratio α = t/w (see Figure 1). If the fibers are circular, this ratio is 1.
The calculation method will consider a representative cylindrical element volume (RVE), having the fibers arranged along the axis of the cylinder. The cross section of this cylinder is large enough in relation to the dimensions of the fibers, so that averaging operations can be performed unaffected by the small non-homogeneities in the structure and by the arrangement of the fibers. The fibers will be considered continuous along the entire length RVE [7,8]. It will be hypothesized that the composite is considered transversely isotropic and homogeneous.
In these hypotheses, we use a mathematical model by which the elastic/viscoelastic properties of the material, will be determined by upper and lower bounds of the properties of the components, which are considered to be known.
Using Hill's theory [9][10][11], one may obtain the bounds on bulk modulus, the longitudinal Young's modulus and the Poisson's ratio. These bounds are determined in terms of the volume ratios of fiber vf and matrix vm.
For a reinforced composite with cylindrical fibers, the upper and lower bounds on the effective elastic moduli have been calculated for a geometry of hexagonal fibers array by Hashin [12][13][14][15].
A new direction of research represents the gradual change of the fiber orientation and the study of material at high stress level and temperature. A better model is presented in the following, where exact estimations for some engineering constants are determined.
In the last few years a series of papers studied the viscoelastic creep response of different material subjected to different loads and temperatures.
The creep behavior of linear viscoelastic materials subjected to a combined tension and torsion loading is presented in [16]. The creep response is expressed using fractional-exponential functions. Experimental tests have validated the model for HDPE and organic glass.
All the phases of creep, primary, secondary, tertiary and rupture, are analyzed in [17], based on a nonlinear viscoelastic creep constitutive formula. The parameters of this model are analyzed and evaluated. Other papers studied different aspects of the creep response of different materials [18][19][20][21][22][23][24].

Mathematical Model
To obtain the elastic/viscoelastic coefficients of the biphasic composite, reinforced with long, anisotropic fibers, the Mori-Tanaka theory will be used [25]. For this, RVE as defined in the Introduction section of this paper, will be used, consisting of an epoxy matrix with viscoelastic behavior, reinforced with monotonously aligned and parallel fibers, uniformly distributed within the matrix. Considering the above, the composite obtained will be regarded as orthotropic. If, however, the elliptical cylinders fibers are randomly oriented, the composite body under consideration will be viewed as transverse isotropic. The elliptical fibers are defined by the ratio w t / = α (see Figure 1). If the fibers are circular, this ratio is 1. Eshelby's [26] approach is used by Zhao and Weng [27] for a cylindrical fiber with an elliptic section in conjunction with Mori-Tanaka's [25] mean-field theory.
We mention that in the paper [28], the two phases are considered to be isotropic materials. In our work it is considered that the fibers have an anisotropic behavior. This approach represents the contribution of the current investigation to the study of these materials for determining the elastic/viscoelastic constants for this type of composite. Obtaining engineering constants as close as possible to reality is of particular importance in practical applications. Applying the Mori-Tanaka method, it is possible to calculate all the necessary modules for characterizing an orthotropic material, in particular for characterizing the transversely isotropic materials. The compliance matrix can be obtained and allows an experimental verification of the results, by obtaining the creep curves.
Let us now consider a two-phase composite reinforced with fibers uniformly dispersed in the matrix material. We denote a comparison material with (CM). Consider both of the RVE, the real composite and the CM, subjected to the same boundary traction σ. Boundary traction represents a traction applied at the end (boundary surface) of the RVE. The loading condition for both, the RVE and the real composite is the same, a traction noted σ. In the following, we will denote by C m the elastic coefficients matrix and by C f the elastic coefficients of the fiber. It is obvious that under the traction σ, the average strain field in the matrix differs from the average strain in (CM), where the mean stress is σ. Let ε represent the difference between the average values of the strain. The same thing happens with the stress field. Between the two kinds of material exists a different mean stress σ. We can conclude: In CM, the relation between the mean strain field ε • and the mean stress field σ is: The mean strain fields in the RVE of the composite are ε m = ε 0 + ε and the mean stress field is The mean strain field in the fiber differs from that in the matrix through an additional term ε pt and hence ε f = ε m + ε pt = ε 0 + ε + ε pt . A similar mean stress field differs by the term σ pt and therefore, σ f = σ + σ + σ pt . The stress-strain relation becomes: or: where the following relation holds: In Appendix A, the Eshelby's transformation tensor P for our case is presented, where the symmetry property P ikjl = P jikl = P ijtk is valid.
To compute the average stress in the whole RVE, the well-known relation as shown below is used: which, in our case, reduces to: For the strain, using a similar procedure, it results in: Hereby I is denoted as the unit tensor. Introducing (7) into (4) yields: which after preliminary calculation reduces to: or: and : Finally, it reduces to: from where: The coefficients A ij can be found in Appendix A [11]. Using the definitions for the shear strain, we have [11]: Equations (13)- (17) shown above will be utilized to determine the elastic/viscoelastic coefficients of an RVE considered as an orthotropic body. The results will be applied to a composite reinforced with graphite/epoxy fibers, while this material is considered as transversely isotropic.
To determine the longitudinal Young's modulus E m of an orthotropic body, we subject to a pure traction σ 11 the studied composite and the comparison material. Then, it follows that σ 11 = E 11 ε 11 for the composite and σ 11 = E m ε 0 11 ; ε 0 22 = ε 0 33 = −ν m ε 0 11 if the comparison material is considered. Using rel. (13) we have: where: a ij = A ij /A, A ij and A are defined in Appendix A; see rel. (A6). It follows that: Similarly, the Young moduli are obtained for the other directions: and To determine the shear moduli, the following relations are used: Recall that: Comparing rel. (22) with rel. (23) it is obtained for G 12 : The other shear moduli are obtained using similar procedures: and To obtain the Poisson's ratio, the following relations are used: Note that: Also or: Now, substitution of (29) and (30) offers: Polymers 2020, 12, 2481 6 of 16 which can be rearranged to: Similarly, it can be shown that: and The method presented in this paper offers the elastic/viscoelastic constants of the composite through uniquely determined values. This is in contrast to methods mostly used for determining these values, which are the bounded methods, the current approach provides the lower and the upper bounds of those values. The originality of the method is dependent on the fact that the fibers are not isotropic but have an anisotropic or transverse isotropic behavior, a situation that can often occur in technique and applications. The time response of the material is evaluated using Schapery's constitutive equation [28]. In this way, nonlinear viscoelastic materials can be studied. So, experimental creep curves can be determined and, in this way, the relations obtained in the paper can be verified.

Testing Facilities and Results
Within the testing program designed to verify the results obtained in the paper, we aimed to design an experimental system that would allow the simultaneous testing of a large number of samples, in different loading conditions and different temperatures.
Most tests have been realized at 23 • C considering the ambient temperature. The relative humidity was controlled in the room where the measurements were performed.
The test period in which the significant results were obtained was 10 h. As a result, tests were performed for all samples for this interval. The experimental setup for creep testing consists of fifteen "single lever" arrangement traction test stations and seven "double lever" arrangement stations with a 10:1 and 25:1 multiplication ratio, respectively [29]. The system is computer controlled and tests on all existing test stations can be performed simultaneously.
The load level for each station is electronically controlled. The creep stress that occurs in the tested specimen is also measured electronically. The system is based on measurements with foreign gages. The installation makes up for the temperature variations that can appear and distort the results.
The disadvantages of this method consist in the high costs of gages taking into account the great number of tests. A strain gage can be used only once. Other disadvantages consist in the great number of fix-up operations, difficulties in adhering the gages, calibration of gages, etc. An alternative is to use extensometers. Figures 2-12 show creep curves for the investigated two-phase composite with an isotropic viscoelastic matrix and transversally isotropic fibers. It is easily seen that fair to good prediction of the experimental data is possible with the theoretical procedure introduced in this section. Further discussion of the results in terms of accuracy and comparison with experimental data is presented in the Conclusion section. The nonlinear viscoelastic characterization was established by testing composites at different applied stresses and test temperature. Five stress levels ranging between 10 and 70% were applied. Each of these stresses was tested at four temperatures, while lowering the level of applied stress with increasing temperature. Some of the obtained results are presented in the following figures. The creep behavior of polyether-ether-ketone (PEEK) and epoxy resin laminates are presented. In these cases, the time response obtained with the proposed relations agrees very well with the experimental results. A study of the diagram of the transverse strain ε 22 shows a very accurate Polymers 2020, 12, 2481 7 of 16 prediction of mechanical constants obtained at a relatively high temperature of 120 • C. Figure 5 shows that the deviation between the theoretical and experimental results is less than 8%.  The creep curves of the epoxy resin; more values of load at 80 °C are presented in Figure 4 and at 120 °C in Figure 5.  The creep curves of the epoxy resin; more values of load at 80 °C are presented in Figure 4 and at 120 °C in Figure 5. The creep curves of the epoxy resin; more values of load at 80 • C are presented in Figure 4 and at 120 • C in Figure 5.

Conclusions
In the first part of our study we have established theoretical values of some mechanical characteristics of a composite material with two components. In this paper it was considered that the reinforcing fibers are homogeneous and transversely isotropic. The matrix is isotropic. With these hypotheses, a prediction was made of the values of elastic/viscoelastic constants of a biphasic composite material depending on the values of elastic constants of the constituent materials [30]. The method applied in the paper used the field method of Mori-Tanaka, extended into the visco-elastic domain to determine the time response of the composite and creep curves [31][32][33][34][35][36][37][38][39][40]. Thus, experimental results can be obtained to validate the theoretical results obtained. A good agreement was obtained between experimental and predicted results for pure resin and 90 degree test specimens.
The same result was not obtained for the 45 degree specimen, where the differences are slightly larger, although they are also very close to the theoretical results. In order to increase the precision, it is proposed to use more precise models, which should take into account the change of the orientation of the fibers, especially at high values of applied tension and temperature.
Creep tests were performed in the current study to determine the nonlinear behavior of neat and carbon fiber-reinforced polyether-ether-ketone (PEEK) and epoxy resin, on {90}]4s and {+45}4s laminates. The same laminates were used to study the transverse properties and shear properties.
The number of experiments in which the temperature, the applied stress and the combination of materials varied was large enough, after comparing the experimental results with the theoretical ones, to provide confidence in the calculation formulas obtained [4].

Conclusions
In the first part of our study we have established theoretical values of some mechanical characteristics of a composite material with two components. In this paper it was considered that the reinforcing fibers are homogeneous and transversely isotropic. The matrix is isotropic. With these hypotheses, a prediction was made of the values of elastic/viscoelastic constants of a biphasic composite material depending on the values of elastic constants of the constituent materials [30]. The method applied in the paper used the field method of Mori-Tanaka, extended into the visco-elastic domain to determine the time response of the composite and creep curves [31][32][33][34][35][36][37][38][39][40]. Thus, experimental results can be obtained to validate the theoretical results obtained. A good agreement was obtained between experimental and predicted results for pure resin and 90 degree test specimens.
The same result was not obtained for the 45 degree specimen, where the differences are slightly larger, although they are also very close to the theoretical results. In order to increase the precision, it is proposed to use more precise models, which should take into account the change of the orientation of the fibers, especially at high values of applied tension and temperature.
Creep tests were performed in the current study to determine the nonlinear behavior of neat and carbon fiber-reinforced polyether-ether-ketone (PEEK) and epoxy resin, on {90}] 4s and {±45} 4s laminates. The same laminates were used to study the transverse properties and shear properties.
The number of experiments in which the temperature, the applied stress and the combination of materials varied was large enough, after comparing the experimental results with the theoretical ones, to provide confidence in the calculation formulas obtained [4].

Eshelby's Tensor for an Elliptic Cylinder
The Eshelby's P ijkl tensor written for fibers with an elliptical cross section has the components [28] and P 1111 = 0 ; P 1122 = 0 ; P 1133 = 0 (A2) where all other P ijkl = 0. Relation (9) offers the liaison between ε 0 ij and ε * ij . For this it is necessary to know that: (see Reference [25]). with the following expressions for N i : where λ f and λ m are the Lame's constants, respectively, for the fiber matrix. From (A3) it results: