Modeling Study of the Creep Behavior of Carbon-Fiber-Reinforced Composites: A Review

Abstract

The aim of this paper is to present some important practical cases in the analysis of the creep response of unidirectional fiber-reinforced composites. Some of the currently used models are described: the micromechanical model, homogenization technics, the Mori–Tanaka method, and the finite element method (FEM). Each method was analyzed to determine its advantages and disadvantages. Regarding the accuracy of the obtained results, comparisons are made with experimental tests. The methods presented here are applied to carbon-fiber-reinforced composites, but these considerations can also be applied to other types of composite materials.

1. Introduction

The creep phenomenon that can occur in viscoelastic materials is defined as manifesting in three hypostases: primary, secondary, and tertiary. The creep phenomenon is defined as a deformation in time of the studied material, if it is loaded with a known force [1] (Figure 1). Creep phenomena usually manifest at high temperatures. However, there are situations in which the creep can appear at lower temperatures, for example, at room temperature, for some types of materials.

Of course, this phenomenon, which manifests in the elongation of the material over time, can become dangerous in the operation of a machine. Figure 1 shows the three intervals of creep behavior. Current applications refer mostly to the first two stages of creep, when the deformation rate is relatively high. In the primary creep stage, a high rate is observed at the beginning, which slows down over time. The aspect of the creep curve depends on the material, load, and time. In the secondary creep stage, there is a relatively constant rate. A high rate of deformation characterizes the third creep stage. The time interval in which this high increase is observed is short and is associated with the destruction of the material. In engineering practice, it is not necessary to reach this stage; as a result, the study of behavior in this area has not attracted much attention. Designers must know the rate of deformation. This can be determined using measurements or by using a verified calculus model. The paper presents such models, which are useful for design activities [2,3]. Creep behavior is interesting for engineers and studies on this phenomenon are numerous [4,5,6].

The technology of advanced composites has developed to the point where these materials are being increasingly utilized in the commercial, military, and aerospace industries, among others. Composite materials are ideal for structural applications where high strength-to-weight and stiffness-to-weight ratios, improved fatigue resistance, and improved dimensional stability are required. Reinforced fiber polymers date back to the early years of the last century. There are two major steps in the manufacturing of polymer-based laminated composites, namely layup and curing. In the layup stage, continuous filaments are arranged in unidirectional laminae or are interwoven. The fibers are often impregnated with resinous material, such as polyester resin, which later serves as the matrix material. The next step, thermal curing, involves the drying or polymerization of the resinous matrix material and is accomplished in suitable autoclaves. The aim is to form a permanent bond between the fibers and the matrix, as well as between the laminae, in order to obtain lightweight, stiff panels [2].

The materials used in engineering have different purposes and are manufactured according to different technologies; as a result, they have a variety of properties. The creep diagrams of these materials can be very different, even under the same loading and temperature conditions. The simplest way to construct a creep diagram is to perform experimental measurements. However, such an approach is expensive and time consuming. Loads with different constant loads must be considered, and tests must be performed at different temperatures.

In [7], a scheme for accelerated characterization is proposed to analyze the viscoelastic response of general laminated composites. The use of this scheme allows a small number of experimental measurements to be performed. The measurements allow for short-term tests at high temperature, to predict the long-term response [8,9,10].

It would be much more advantageous for designers to have useful creep models which could be used to obtain creep diagrams by calculation.

To study the nonlinear viscoelastic behavior of a unidirectional composite, the well-known FEM method is applied. The symmetry properties of the composite allow for the simplification of such an analysis. A good correlation with the FEM micromechanics models developed in [11] is obtained. The method can also be used to study a composite with a complex topology [12,13,14]. Such a description also offers the possibility of studying the material in a wide range of boundary conditions. Thus, the thermal effects and the expansion due to humidity were included through the initial conditions. In [15], the above equations were used for unidirectional composites reinforced with graphite and glass.

The works [16,17,18] improve the classic models used in the case of nonlinear behavior. An empirical model was developed to achieve this. A method that can be easily implemented using a numerical procedure was thus obtained.

Based on the previously presented studies [16,17,18], a nonlinear viscoelastic model was developed in [19,20]. The developed model and the experimental measurements taken for test specimens allowed for an orthotropic material. The presented procedure can also be applied to study the long-term nonlinear viscoelastic response of laminates.

Other research [21] has shown that a law moisture concentration (at about 1%) can be a critical limit for carbon epoxy laminates. When this limit is exceeded, the viscoelastic rate of deformation occurs faster.

The study of a material made of an epoxy resin reinforced with unidirectional aramid fibers by tests and measurements at high temperatures is presented in [22]. An appropriate mathematical model for this study proved to be the “power law” which can describe behavior in both the linear and nonlinear domains, so that it can model viscoelastic behavior. To study the behavior in the nonlinear field, some nonlinear viscoelastic coefficients are introduced (these coefficients depend both on the stresses to which the materials are subjected and on the temperature). This method of analysis was proven to concur with the nonlinear model presented in [12].

In [23,24,25,26], a variational principle is used in which the time variable also appears, using a relatively simple mathematical description. In [27], the heat-induced stress field in the components of a polymer composite at low temperatures is studied (one application is considered for spacecraft). The geometry of the composite microstructure proves to be important in terms of the field of stresses and the deformation of this type of material under the conditions described above.

In [28], all the engineering constants that define one orthotropic and one transverse isotropic composite are determined. For a transverse isotropic material, the results [29,30,31,32] provide us with the upper and lower limits of engineering constants. In [33], the Mori–Tanaka method presented in [34] is extended.

Paper [35] shows the non-linear viscoelastic/viscoplastic behaviors of graphite/bismaleimide. Paper [12] presents a nonlinear formulation used to study materials at temperatures above 93 °C. A micromechanical analysis for the study of the behavior of a fiber-reinforced composite is described in [36,37]. These studies show good concordance with the findings presented in [11]. Biphasic materials and their mechanical properties have been extensively studied in numerous papers published in recent years [38,39,40,41,42,43,44,45,46,47,48,49]. New results are presented in [50,51,52,53,54].

In this review, the authors present more factors related to the analysis of the creep behavior of a composite material reinforced with fibers. The model’s proposed offer results and a creep curve in the case of different loads. The results presented in this review are mainly based on the results obtained in [55,56,57,58,59].

The creep calculation of composite materials represents an important step in the process of designing a new material. A series of methods are therefore developed to achieve this objective. The problem remains an important one in the context of unprecedented advances in the development of new materials, with increasing numbers of properties that are useful in various applications. To the knowledge of the authors, the systematization and unitary presentation of these methods has not yet been achieved. This study thus makes a significant contribution to the field. The methods based on the homogenization theory are presented in Section 2, Section 3 and Section 4, and those based on the FEM theory are presented in Section 5.

2. The Micromechanical Model in Homogenization Theory

2.1. Model and Constitutive Law

The method of the micromechanical model aims to obtain the overall mechanical parameters of a composite based on models that use the parameters of the individual constituents of the composite and the interaction that exists between them. Consider a unidirectional composite with randomly distributed fibers in the matrix material. In the models used, it is normal to take into account a periodicity in the distribution of fibers. In this way, the existing periodicities allow us to simplify the analysis. Figure 2 presents such a model. The following assumptions can be formulated:

• The fibers are continuous and circular, and oriented in the X₁ direction. They are positioned regularly in a rectangular array in the transversal X₂–X₃ plane;
• The fibers are linearly elastic and anisotropic. The matrix is isotropic and nonlinearly viscoelastic;
• No cracks or holes appear or develop, and the contact fiber matrix is mechanical.

Using the proposed model, it is possible to determine the response of the material if only a single repeating unit cell (RUC) is studied—see Figure 3. As such, the complexity of the problem can be significantly reduced. For this analysis, it is sufficient to study only a quarter of a fiber, as in Figure 3b. The main hypothesis of the theory is that the RUCs are very small, reported to be equal to the dimensions of the studied material.

The RUC refers to a local coordinate system (X₁, $x_2^{(\lambda )}$, $x_3^{(\lambda )}$) (Figure 4). The displacement in each subcell is defined through the following formulas [27,28]:

$$u_i^{(\lambda )}=u_i^{o(\lambda )}+x_2^{(\lambda )}{\xi }_i^{(\lambda )}+x_3^{(\lambda )}{\zeta }_i^{(\lambda )}\hspace{1em};\hspace{1em}i=1,2,3\hspace{0.17em}.$$

(1)

Here, $u_i^o$ is the displacement component of the origin and “$\lambda$” represents both the fiber (when $\lambda =f$) and the matrix (when $\lambda =m$).

Considering the material to be linear, the strain–displacement relations are as follows:

$${\epsilon }_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\hspace{1em};\hspace{1em}i,j=1,2,3 ,$$

(2)

Equation (2) can be written for fiber and matrix in the unified form:

$${\epsilon }_{ij}^{(\lambda )}=\frac{1}{2}\left(\frac{\partial u_i^{(\lambda )}}{\partial x_j}+\frac{\partial u_j^{(\lambda )}}{\partial x_i}\right)\hspace{1em};\hspace{1em}i,j=1,2,3$$

(3)

If $i\ne j$, it can be written as

$${\epsilon }_{ij}^{(\lambda )}=\frac{1}{2}\left(\frac{\partial u_i^{(\lambda )}}{\partial x_j}+\frac{\partial u_j^{(\lambda )}}{\partial x_i}\right)=\frac{{\gamma }_{ij}^{(\lambda )}}{2}\hspace{1em};\hspace{1em}i,j=1,2,3 ;\hspace{1em}i\ne j$$

(4)

The engineering shear strain is denoted as ${\gamma }_{ij}^{(\lambda )}=2{\epsilon }_{ij}^{(\lambda )}\hspace{1em};\hspace{1em}i,j=1,2,3$ $;\hspace{1em}i\ne j$.

Using Equation (1) into Equation (2) and considering Equations (3) and (4), the following equations can be obtained:

$${\epsilon }_{11}^{(\lambda )}=\frac{\partial u_i^{o(\lambda )}}{\partial X_1}$$

(5)

$${\epsilon }_{11}^{(\lambda )}=\frac{\partial {\mathit{u}}_{\mathit{i}}^{\mathit{o}(\lambda )}}{\partial {\mathit{X}}_1}$$

(6)

$${\epsilon }_{33}^{(\lambda )}={\zeta }_3^{(\lambda )}$$

(7)

$${\gamma }_{23}^{(\lambda )}=\left[{\xi }_3^{(\lambda )}+{\zeta }_2^{(\lambda )}\right]$$

(8)

$${\gamma }_{31}^{(\lambda )}=\left[\frac{\partial u_3^{o(\lambda )}}{\partial X_1}+{\zeta }_1^{(\lambda )}\right]$$

(9)

$${\gamma }_{12}^{(\lambda )}=\left[\frac{\partial u_2^{o(\lambda )}}{\partial X_1}+{\xi }_1^{(\lambda )}\right]$$

(10)

Considering a linear and transversely isotropic composite, the constitutive equation can be written as

$${\left\{\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ {\gamma }_{23}\\ {\gamma }_{31}\\ {\gamma }_{12}\end{array}\right\}}^{(\lambda )}={\left[\begin{array}{cccccc}{S}_{11}& S_{12}& S_{12}& 0& 0& 0\\ S_{12}& S_{22}& S_{23}& 0& 0& 0\\ S_{12}& S_{23}& S_{22}& 0& 0& 0\\ 0& 0& 0& S_{44}& 0& 0\\ 0& 0& 0& 0& S_{66}& 0\\ 0& 0& 0& 0& 0& S_{66}\end{array}\right]}^{(\lambda )}{\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{11}\\ {\sigma }_{11}\\ {\tau }_{23}\\ {\tau }_{31}\\ {\tau }_{12}\end{array}\right\}}^{(\lambda )}$$

(11)

or, considering the expression of the engineering constant for this type of material:

$${\left\{\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ {\gamma }_{23}\\ {\gamma }_{31}\\ {\gamma }_{12}\end{array}\right\}}^{(\lambda )}={\left[\begin{array}{cccccc}\frac{1}{E_{11}}& -\frac{{\nu }_{12}}{E_{22}}& -\frac{{\nu }_{12}}{E_{22}}& 0& 0& 0\\ -\frac{{\nu }_{12}}{E_{22}}& \frac{1}{E_{22}}& -\frac{{\nu }_{23}}{E_{22}}& 0& 0& 0\\ -\frac{{\nu }_{12}}{E_{22}}& -\frac{{\nu }_{23}}{E_{22}}& \frac{1}{E_{22}}& 0& 0& 0\\ 0& 0& 0& \frac{1}{G_{23}}& 0& 0\\ 0& 0& 0& 0& \frac{1}{G_{12}}& 0\\ 0& 0& 0& 0& 0& \frac{1}{G_{12}}\end{array}\right]}^{(\lambda )}{\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\tau }_{23}\\ {\tau }_{31}\\ {\tau }_{12}\end{array}\right\}}^{(\lambda )}$$

(12)

Here, $E_{11}$ and $E_{22}=E_{33}$ are Young’s moduli, $G_{23}$ and $G_{12}=G_{13}$ are the shear moduli, and ${\nu }_{23}$ and ${\nu }_{12}={\nu }_{13}$ are the Poisson ratios. The direction of anisotropy is, in our model, X₁, and the plane of isotropy is X₂–X₃. From Equation (12), the following formula is obtained:

$${\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\tau }_{23}\\ {\tau }_{31}\\ {\tau }_{12}\end{array}\right\}}^{(f)}={\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{12}& 0& 0& 0\\ C_{12}& C_{22}& C_{23}& 0& 0& 0\\ C_{12}& C_{23}& C_{22}& 0& 0& 0\\ 0& 0& 0& C_{66}& 0& 0\\ 0& 0& 0& 0& C_{44}& 0\\ 0& 0& 0& 0& 0& C_{44}\end{array}\right]}^{(f)}{\left\{\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ {\gamma }_{23}\\ {\gamma }_{31}\\ {\gamma }_{12}\end{array}\right\}}^{(f)}$$

(13)

with

$$C_{66}=\frac{C_{22}-C_{23}}{2}$$

(14)

Equation (13) can be written as

$${\sigma }^{(f)}=C^{(f)}{\epsilon }^{(f)}$$

(15)

The behavior of a viscoelastic material can be described using Boltzmann’s superposition principle and the results presented in [12]:

$$\epsilon (t)=D_n\hspace{0.17em}{\sigma }_n$$

(16)

with

$$D_n=g_o{D}_o+g{\displaystyle \sum _{j=1}^m{D}_j\left(1-e^{-t/r_j}\right)}$$

(17)

The paper [12] presents us with the possibility of writing the constitutive equations as

$${\epsilon }_{ij}^{(m)}=D_n\left[1+\nu (t)\right]{\sigma }_{ij}^{(m)}-D_n\hspace{0.17em}\nu (t)\hspace{0.17em}{\sigma }_{kk}^{(m)}{\delta }_{ij}$$

(18)

In Equation (18), $D_n$ is obtained using Equation (17), $\nu (t)$ is the Poisson ratio (in our study, this is considered to be independent of time), and ${\delta }_{ij}$ is Kronecker’s delta.

The first step must be to determine the average stresses, and then the average strains. Thus, the general behavior of the material is obtained based on the average stresses and average strains in a RUC.

2.2. Average Stress

In the proposed model, the RUC is considered to be a rectangular parallelepiped with parallel edges. The reference frame axes are (X₁, X₂, X₃) of the volume V. This will be determined as the average stress ${\overline{\sigma }}_{ij}$ in V. This can be obtained via the following relation:

$${\overline{\sigma }}_{ij}=\frac{1}{V}{\displaystyle \underset{V}{\int }{S}_{ij}dV}$$

(19)

Considering one-quarter of a cell, this relation can be written as

$${\overline{\sigma }}_{ij}=\frac{1}{A}\left({\overline{S}}_{ij}^{(f)}{A}_f+{\overline{S}}_{ij}^{(m)}{A}_m\right)$$

(20)

where ${\overline{S}}_{ij}^{(\lambda )}$ are the average stresses. Now, consider a unit depth of the RUC, i.e., V = A × 1. Using the notation presented in Figure 3b, the following is obtained:

$$A={(R+h/2)}^2; A_f^{}=\pi R^2/4; A_m={(R+h/2)}^2-\pi R^2/4$$

(21)

$${\overline{\sigma }}_{ij}=\frac{1}{{\left(R+\frac{h}{2}\right)}^2}\left\{\frac{\pi R^2}{4}{\overline{S}}_{ij}^{(f)}+\left[{\left(R+\frac{h}{2}\right)}^2-\frac{\pi R^2}{4}\right]{\overline{S}}_{ij}^{(m)}\right\}$$

(22)

The partial average stress ${\overline{\mathit{S}}}_{ij}^{(\lambda )}$ is obtained with

$${\overline{S}}_{ij}^{(\lambda )}=\frac{1}{A_{\nu }}{\displaystyle \underset{A}{\int }{\sigma }_{ij}^{(\lambda )}}dA=\frac{1}{A_{\nu }}{\displaystyle \iint {\sigma }_{ij}^{(\lambda )}d{x}_{2}d{x}_3}$$

(23)

Using polar coordinates, the Jacobian can be obtained:

$$J=\frac{\partial (x_2,x_3)}{\partial (r,\theta )}=\left|\begin{array}{cc}\cos \theta & \sin \theta \\ -r\sin \theta & r\cos \theta \end{array}\right|=r$$

(24)

Moreover, Equation (23) for fiber (subcell “f”) becomes

$$\mathrm{e}{\overline{S}}_{ij}^{(f)}=\frac{4}{\pi R^2}{\displaystyle {\int }_0^{\pi /2}{\displaystyle {\int }_0^R{\sigma }_{ij}^{(f)}}}rdrd\theta$$

(25)

where ${\sigma }_{ij}^{(f)}$ is given by Equation (13). Introducing Equations (5)–(10) and (11) into Equation (25) leads to

$${\overline{S}}_{ij}^{(f)}=\left\{\begin{array}{c}{C}_{11}{\epsilon }_{11}^{(f)}+C_{12}\left({\epsilon }_{22}^{(f)}+{\epsilon }_{33}^{(f)}\right)\\ C_{12}{\epsilon }_{11}^{(f)}+C_{22}{\epsilon }_{22}^{(f)}+C_{23}{\epsilon }_{33}^{(f)}\\ C_{12}{\epsilon }_{11}^{(f)}+C_{23}{\epsilon }_{22}^{(f)}+C_{22}{\epsilon }_{33}^{(f)}\\ C_{66}{\gamma }_{23}^{(f)}\\ C_{44}{\gamma }_{31}^{(f)}\\ C_{44}{\gamma }_{12}^{(f)}\end{array}\right\}$$

(26)

or

$${\overline{S}}_{ij}^{(f)}=\left\{\begin{array}{c}{C}_{11}\frac{\partial u_1^{o(f)}}{\partial X_1}+C_{12}\left({\xi }_2^{(f)}+{\zeta }_3^{(f)}\right)\\ C_{12}\frac{\partial u_1^{o(f)}}{\partial X_1}+C_{22}{\xi }_2^{(f)}+C_{23}{\zeta }_3^{(f)}\\ C_{12}\frac{\partial u_1^{o(f)}}{\partial X_1}+C_{23}{\xi }_2^{(f)}+C_{22}{\zeta }_3^{(f)}\\ C_{66}\left[{\xi }_3^{(f)}+{\zeta }_2^{(f)}\right]\\ C_{44}\left[\frac{\partial u_3^{o(f)}}{\partial X_1}+{\zeta }_1^{(f)}\right]\\ C_{44}\left[\frac{\partial u_2^{o(f)}}{\partial X_1}+{\xi }_1^{(f)}\right]\end{array}\right\}$$

(27)

The average stresses in the matrix (subcell “m”) are determined using the following relations:

$${\overline{S}}_{ij}^{(m)}=\frac{1}{{\left(R+\frac{h}{2}\right)}^2-\frac{\pi R^2}{4}}\left({\displaystyle {\int }_0^{R+h/2}{\displaystyle {\int }_0^{R+h/2}{\sigma }_{ij}^{(m)}d{x}_{2}d{x}_3-{\displaystyle {\int }_0^{\pi /2}{\displaystyle {\int }_0^R{\sigma }_{ij}^{(m)}rdrd\theta }}}}\right)$$

(28)

Equation (15) together with Equations (5)–(10) yields

$$\frac{\partial u_1^o}{\partial X_1}=D_n\left[1+\nu (t)\right]{\overline{S}}_{11}^{(m)}-D_n\left[\nu (t)\right]{\overline{S}}_{kk}^{(m)}$$

(29)

$${\xi }_2^{(m)}=D_n\left[1+\nu (t)\right]{\overline{S}}_{22}^{(m)}-D_n\left[\nu (t)\right]{\overline{S}}_{kk}^{(m)}$$

(30)

$${\zeta }_3^{(m)}=D_n\left[1+\nu (t)\right]{\overline{S}}_{33}^{(m)}-D_n\left[\nu (t)\right]{\overline{S}}_{kk}^{(m)}$$

(31)

$${\xi }_3^{(m)}+{\zeta }_2^{(m)}=2D_n\left[1+\nu (t)\right]{\overline{S}}_{23}^{(m)}$$

(32)

$$\frac{\partial u_3^{o(m)}}{\partial X_1}+{\zeta }_1^{(m)}=2D_n\left[1+\nu (t)\right]{\overline{S}}_{31}^{(m)}$$

(33)

$$\frac{\partial u_2^{o(m)}}{\partial X_1}+{\xi }_1^{(m)}=2D_n\left[1+\nu (t)\right]{\overline{S}}_{12}^{(m)}$$

(34)

2.3. Continuity Conditions

In a RUC, the conditions of continuity of the movements at the interface between sub-cells must be ensured. These conditions must be assured in both the X₂ and X₃ directions. From Figure 4, the following relations hold:

$$x_2^{(f)}=X_2^i\mp X_2^{(f)}=\mp R\cos \theta$$

(35)

$$x_2^{(m)}=X_2^i\mp X_2^{(m)}=\pm \left(\frac{h}{2}+R-R\cos \theta \right)$$

(36)

where there are $\theta$ located points on the interface.

Introducing Equations (35) and (36) into Equation (1) for the cases when $\lambda =f$ and $\lambda =m$ results in

$$u_i^{(f)}=u_i^{o(f)}\mp R\cos \theta {\xi }_i^{(f)}+x_3^{(f)}{\zeta }_i^{(f)}$$

(37)

$$u_i^{(m)}=u_i^{o(m)}\pm \left(\frac{h}{2}+R-R\cos \theta \right){\xi }_i^{(m)}+x_3^{(m)}{\zeta }_i^{(m)}$$

(38)

where

$$u_i^{o(f)}=u_i^{o(f)}(X_2^{(f)})$$

(39)

$$u_i^{o(m)}=u_i^{o(m)}(X_2^{(m)})$$

(40)

The continuity of the displacements at the interface is considered in the average sense. This is expressed by the following relations:

$${\displaystyle {\int }_{-\pi /2}^{\pi /2}\left[u_i^{o(f)}\mp R\cos \theta {\xi }_i^{(f)}+x_3^{(f)}{\zeta }_i^{(f)}\right]R\cos \theta d\theta =}$$

$${\displaystyle {\int }_{-\pi /2}^{\pi /2}\left[u_i^{o(m)}\pm \left(\frac{h}{2}+R-R\cos \theta \right){\xi }_i^{(m)}+x_3^{(m)}{\zeta }_i^{(m)}\right]R\cos \theta d\theta }$$

(41)

Equation (41) produces

$$u_i^{o(f)}\left(X_2^{(f)}\right)\pm \frac{\pi R}{2}{\xi }_i^{(f)}=u_i^{o(m)}\left(X_2^{(m)}\right)\mp \left(h+2R-\frac{\pi R}{2}\right){\xi }_i^{(m)}$$

(42)

The addition of the two Equations (42) offers us

$$u_i^{o(f)}=u_i^{o(m)}=u_i^o$$

(43)

which represent the continuity conditions for displacement.

2.4. Average Strain

Considering the composite specimen presented earlier, and the continuity conditions, the average of the strains over volume is

$${\overline{\epsilon }}_{ij}=\frac{1}{V}{\displaystyle {\int }_V{\epsilon }_{ij}dV}$$

(44)

For the representative cell studied, the following formula is obtained:

$${\overline{\epsilon }}_{ij}=\frac{1}{A}{\displaystyle \sum _{\lambda =f,m}{\overline{\epsilon }}_{ij}^{(\lambda )}{A}_{\lambda }}$$

(45)

Here $A=A_m+A_f$ and ${\overline{\epsilon }}_{ij}^{(\lambda )}$ are the strains obtained using Equations (5)–(10) ($\lambda =f,m$). Considering Equations (1) and (3) (considering $i=j=1$), we obtain

$${\overline{\epsilon }}_{ij}^{(\lambda )}=\frac{\partial u_1^o}{\partial X_1}$$

(46)

Equation (46) together with Equation (45) yields

$${\overline{\epsilon }}_{11}=\frac{1}{{\left(R+\frac{h}{2}\right)}^2}\left\{\frac{\pi R^2}{4}\frac{\partial u_1^o}{\partial X_1}+\left[{\left(R+\frac{h}{2}\right)}^2-\frac{\pi R^2}{4}\right]\frac{\partial u_1^o}{\partial X_1}\right\}$$

(47)

or

$${\overline{\epsilon }}_{11}=\frac{\partial u_1^o}{\partial X_1}$$

(48)

The octahedral shear stress in the matrix is

$$S_{oct}^{(m)}={\left\{\frac{1}{2}\left[{\left({\overline{S}}_{11}^{(m)}-{\overline{S}}_{22}^{(m)}\right)}^2+{\left({\overline{S}}_{22}^{(m)}-{\overline{S}}_{33}^{(m)}\right)}^2+{\left({\overline{S}}_{33}^{(m)}-{\overline{S}}_{11}^{(m)}\right)}^2\right]+3\left[{\left({\overline{S}}_{12}^{(m)}\right)}^2+{\left({\overline{S}}_{23}^{(m)}\right)}^2+{\left({\overline{S}}_{31}^{(m)}\right)}^2\right]\right\}}^{\frac{1}{2}}$$

(49)

The unknowns are calculated using an incremental procedure. The unknowns are

• the six values of stresses in the two subcells: ${\overline{S}}_{11}^{(\lambda )}$, ${\overline{S}}_{22}^{(\lambda )}$, and ${\overline{S}}_{33}^{(\lambda )}$;
• the four micro-variables in the two subcells: ${\xi }_2^{(\lambda )}$ and ${\zeta }_3^{(\lambda )}$;
• the three strains in the composite: ${\overline{\epsilon }}_{11}$, ${\overline{\epsilon }}_{22}$, and ${\overline{\epsilon }}_{33}$.

Now, the connection between the stress in the matrix and in the fiber of a RUC must be determined. Using the assumptions proposed in [36,37], the shear stress in the X₂ direction is

$${\overline{S}}_{22}^{(f)}={\alpha }_f{\overline{\sigma }}_{22}$$

(50)

for the fiber and

$${\overline{S}}_{22}^{(m)}={\alpha }_m{\overline{\sigma }}_{22}$$

(51)

for the matrix from where it results:

$${\overline{S}}_{22}^{(f)}=\frac{{\alpha }_f}{{\alpha }_m}{\overline{S}}_{22}^{(m)}$$

(52)

In the X₃ direction, a similar equation is obtained:

$${\overline{S}}_{33}^{(f)}=\frac{{\beta }_f}{{\beta }_m}{\overline{S}}_{33}^{(m)}$$

(53)

The concentration factors ${\alpha }_{\lambda }$ and ${\beta }_{\lambda }$ are weighting coefficients and should satisfy the following relations:

$${\alpha }_f{v}_f+{\alpha }_m{v}_m=1$$

(54)

and

$${\beta }_f{v}_f+{\beta }_m{v}_m=1$$

(55)

In a particular case, considering that the composite loaded is in only one of the directions X₂ or X₃, the relation (38) in direction X₃ (unloaded) becomes

$${\overline{S}}_{33}^{(f)}=-\frac{v_m}{v_f}{\overline{S}}_{33}^{(m)}$$

(56)

and

$${\overline{S}}_{22}^{(f)}=-\frac{v_m}{v_f}{\overline{S}}_{22}^{(m)}$$

(57)

Consider now the case of a uniaxial load. Therefore, we obtain a linear system with 13 equations and 13 unknowns:

$${\overline{S}}_{11}^{(f)}=C_{11}{\overline{\epsilon }}_{11}+C_{12}\left({\xi }_2^{(f)}+{\zeta }_3^{(f)}\right)$$

(58)

$${\overline{S}}_{22}^{(f)}=C_{12}{\overline{\epsilon }}_{11}+C_{22}{\xi }_2^{(f)}+C_{23}{\zeta }_3^{(f)}$$

(59)

$${\overline{S}}_{33}^{(f)}=C_{12}{\overline{\epsilon }}_{11}+C_{23}{\xi }_2^{(f)}+C_{22}{\zeta }_3^{(f)}$$

(60)

$${\overline{\epsilon }}_{11}=D_n\left[1+\nu (t)\right]{\overline{S}}_{11}^{(m)}-D_n\left[\nu (t)\right]{\overline{S}}_{kk}^{(m)}$$

(61)

$${\xi }_2^{(m)}=D_n\left[1+\nu (t)\right]{\overline{S}}_{22}^{(m)}-D_n\left[\nu (t)\right]{\overline{S}}_{kk}^{(m)}$$

(62)

$${\zeta }_3^{(m)}=D_n\left[1+\nu (t)\right]{\overline{S}}_{33}^{(m)}-D_n\left[\nu (t)\right]{\overline{S}}_{kk}^{(m)}$$

(63)

$${\overline{\epsilon }}_{22}=\frac{1}{A}\left[A_f{\xi }_2^{(f)}+A_m{\xi }_2^{(m)}\right]$$

(64)

$${\overline{\epsilon }}_{33}=\frac{1}{A}\left[A_f{\zeta }_3^{(f)}+A_m{\zeta }_3^{(m)}\right]$$

(65)

$${\overline{S}}_{22}^{(f)}=\frac{{\alpha }_f}{{\alpha }_m}{\overline{S}}_{22}^{(m)}$$

(66)

$${\overline{S}}_{33}^{(f)}=\frac{{\beta }_f}{{\beta }_m}{\overline{S}}_{33}^{(m)}$$

(67)

$${\overline{\sigma }}_{11}=\frac{1}{A}\left[A_f{\overline{S}}_{11}^{(f)}+A_m{\overline{S}}_{11}^{(m)}\right]$$

(68)

$${\overline{\sigma }}_{22}=\frac{1}{A}\left[A_f{\overline{S}}_{22}^{(f)}+A_m{\overline{S}}_{22}^{(m)}\right]$$

(69)

$${\overline{\sigma }}_{33}=\frac{1}{\mathit{A}}\left[{\mathit{A}}_{\mathit{f}}{\overline{\mathit{S}}}_{33}^{(\mathit{f})}+{\mathit{A}}_{\mathit{m}}{\overline{\mathit{S}}}_{33}^{(\mathit{m})}\right]$$

(70)

The analysis presented in this section shows that a micromechanical model for the study of a unidirectional composite can provide good results. Thus, analytical relations are obtained, which then allow for the calculation of the mechanical constants of such a composite and for the study of its behavior in a range of applications. Schapery’s nonlinear constitutive equation for isothermal uniaxial loading conditions is used in the analysis, thus allowing us to consider the nonlinear viscoelastic response of the material. Papers that present many experimental results [55,56,57,58,59] demonstrate the potential of the method.

3. Description of Homogenization Theory

3.1. Overview

The theory of homogenization is a mathematical method used to average the physical properties of inhomogeneous materials. This method has been developed over the last eight decades and is used to analyze and solve differential equations with periodic coefficients. As essentially inhomogeneous materials that have a periodicity or certain symmetries in their structure, composite materials lend themselves very well to the application of these methods that determine the mechanical characteristics of a material. The experimental results validate the methods used by the theory of homogenization [55,56,57,58,59]. For this reason, the homogenization method has been used in numerous cases and engineering applications [33,34,36,37,39,40,41,59] to determine the mechanical properties of multiphase composites. In this method, a transition is made, through homogenization, from a periodic structure to a homogeneous and isotropic or transversely isotropic material throughout its structure [44].

In the research, several analytical and numerical methods have been proposed to solve the problems generated by the application of this method. Experimental results have always shown the predominant acceptance of such methods [45,60]. The interaction between the phases of the composite is modeled by unifying the homogenization problems for heterogeneous elasto-plastic and elasto-viscoplastic materials [61,62]. Other works address the improvement of the method, using the experience gained by different engineering applications [61,62,63,64,65,66,67,68,69]. Other related methods are considered in [70,71]. The object of this research is the development of reliable procedures that can be easily applied by designers. The following section presents the homogenization theory used to determine the mechanical quantities that characterize the viscoelastic material in question. An application is suggested for a composite reinforced with carbon fibers.

3.2. Homogenized Model

One of the advantages offered by the homogenization theory is the possibility of studying differential equations in which the coefficients have rapid variations or periodic variations. Engineering constants, which are useful in engineering practice, are obtained following averaging processes. Thus, a material with a periodic structure can be treated as a homogeneous material. A differential equation with periodic coefficients with large variations is thus replaced in the modeling with an equation with constant coefficients. This is how the continuum concept is extended to micro-structured materials (composite materials also belong to this class). The bases of this mathematical theory are presented in [72,73,74,75,76,77]. In this application, the calculation method is used to analyze the creep response of a unidirectional composite reinforced with carbon fibers.

The stress field ${\sigma }^{\delta }$ for repeating unit cells of size $\delta$ must obey the following equations:

$$\begin{gathered} \frac{\partial {\sigma }_{11}^{\delta }}{\partial x_1}+\frac{\partial {\tau }_{12}^{\delta }}{\partial x_2}+\frac{\partial {\tau }_{13}^{\delta }}{\partial x_3}=f_1(x) \\ \frac{\partial {\tau }_{21}^{\delta }}{\partial x_1}+\frac{\partial {\sigma }_{22}^{\delta }}{\partial x_2}+\frac{\partial {\tau }_{23}^{\delta }}{\partial x_3}=f_2(x) \\ \frac{\partial {\tau }_{31}^{\delta }}{\partial x_1}+\frac{\partial {\tau }_{32}^{\delta }}{\partial x_2}+\frac{\partial {\sigma }_{33}^{\delta }}{\partial x_3}=f_3(x) \end{gathered}$$

(71)

where ${\sigma }_{ij}^{\delta }={\sigma }_{ji}^{\delta }$, for $i,j=1,2,3$.

The contour conditions that must be respected by the displacements are

$${u^{\delta }|}_{{\partial }_1\Omega }=\tilde{u}$$

(72)

The boundary conditions are

$$\begin{gathered} {\sigma }_{11}^{\delta }{n}_1+{\tau }_{12}^{\delta }{n}_2+{\tau }_{13}^{\delta }{n}_3=T_1(x) \\ {\tau }_{21}^{\delta }{n}_1+{\sigma }_{22}^{\delta }{n}_2+{\tau }_{23}^{\delta }{n}_3=T_2(x) \\ {\tau }_{31}^{\delta }{n}_1+{\tau }_{32}^{\delta }{n}_2+{\sigma }_{33}^{\delta }{n}_3=T_3(x) \end{gathered}$$

(73)

on the contour ${\partial }_2\Omega$, (${\partial }_1\Omega \cup {\partial }_2\Omega =\partial \Omega$). Hook’s Law is

$${\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\tau }_{23}\\ {\tau }_{31}\\ {\tau }_{12}\end{array}\right\}}^{\delta }=\left[\begin{array}{cccccc}{C}_{1111}& C_{1122}& C_{1133}& & & \\ C_{2211}& C_{2222}& C_{2233}& & 0& \\ C_{3311}& C_{3322}& C_{3333}& & & \\ & & & C_{2323}& & \\ & 0& & & C_{3131}& \\ & & & & & C_{1212}\end{array}\right]{\left\{\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ {\gamma }_{23}\\ {\gamma }_{31}\\ {\gamma }_{13}\end{array}\right\}}^{\delta }$$

(74)

Or, using a compact notation,

$${\sigma }^{\delta }=C{\epsilon }^{\delta }$$

(75)

The elasticity matrix C is semi-positive definite:

$$C_{ijkh}{x}_{ij}{x}_{kh}\ge \alpha \hspace{0.17em}{x}_{ij}{x}_{kh}$$

(76)

for $\alpha >0$ and $\forall x_{ij},x_{kh}\in \mathit{R}$, $C_{ijkh}(x)$ is a periodical function of x with the period equal with the dimension $\delta$ of the unit cell. Considering a new function y, $y=x/\delta$:

$$C_{ijkh}(x)=C_{ijkh}(y\delta )=C_{ijkh}(y)$$

(77)

and the stress is

$${\sigma }_{ij}^{\delta }={\sigma }_{ij}^o(x,y)+{\sigma }_{ij}^1(x,y)\delta +.....$$

(78)

The dependence of stress on y is “quasi-periodical”. Introducing Equation (78) in Equation (71) produces

$${\delta }^{-1}\frac{\partial {\sigma }_{ij}^o}{\partial y_j}+\left(\frac{\partial {\sigma }_{ij}^o}{\partial x_j}+\frac{\partial {\sigma }_{ij}^1}{\partial y_j}\right)\hspace{0.17em}{\delta }^o+\left(\frac{\partial {\sigma }_{ij}^1}{\partial x_j}+\frac{\partial {\sigma }_{ij}^2}{\partial y_j}\right)\hspace{0.17em}{\delta }^1+.....=f_i(x)$$

(79)

The following relation is used:

$$\frac{d}{dx}\left(f\right)=\frac{\partial f}{\partial x}dx+\frac{\partial f}{\partial x}dy$$

(80)

but $y=x/\delta$, and, thus, $dy=dx/\delta$; so,

$$\frac{d}{dx}\left(f\right)=\frac{\partial f}{\partial x}\left(f\right)+\frac{1}{\delta }\frac{\partial }{\partial y}\left(f\right)$$

(81)

The coefficients of ${\delta }^{-1}$ in Equation (79) must be 0; therefore,

$$\frac{\partial {\sigma }_{ij}^o}{\partial y_j}=0$$

(82)

Equation (82) is called the “local equation”. Identifying the terms of $\delta {\hspace{0.17em}}^0$ produces

$$\frac{\partial {\sigma }_{ij}^o}{\partial x_j}+\frac{\partial {\sigma }_{ij}^1}{\partial y_j}=f_i(x)\hspace{1em}i=1,2,3$$

(83)

Applying the average operator to Equation (83) results in the following equation:

$$\frac{\partial \langle {\sigma }_{ij}^o\rangle }{\partial x_j}+\langle \frac{\partial {\sigma }_{ij}^1}{\partial y_j}\rangle =f_i(x)\hspace{1em}i=1,2,3$$

(84)

but

$$\langle \frac{\partial {\sigma }_{ij}^1}{\partial y_j}\rangle =\frac{1}{V}{\displaystyle {\int }_V\partial {\sigma }_{ij}^{1}dV}=\frac{1}{V}{\displaystyle {\int }_{\partial V}{\sigma }_{ij}^1{n}_{j}dS}=0$$

(85)

The stresses ${\sigma }_{ij}^1$ take equal values on the corresponding points of the boundary of the cell $\Gamma$ (due to the property of periodicity), so

$$\frac{\partial \langle {\sigma }_{ij}^o\rangle }{\partial x_j}=f_i(x)\hspace{1em}i=1,2,3$$

(86)

We state that

$${\epsilon }_{ij,x}(w)=\frac{1}{2}\left(\frac{\partial w_i}{\partial x_j}+\frac{\partial w_j}{\partial x_i}\right)\hspace{1em};\hspace{1em}i,j=1,2,3$$

(87)

$${\epsilon }_{ij,y}(w)=\frac{1}{2}\left(\frac{\partial w_i}{\partial y_j}+\frac{\partial w_j}{\partial y_i}\right)\hspace{1em};\hspace{1em}i,j=1,2,3$$

(88)

The displacement field can be expressed by the series:

$$u(x,y)=u^o(x)+u^1(x,y)\hspace{0.17em}\delta +u^2(x,y)\hspace{0.17em}{\delta }^2+....$$

(89)

$u^o(x)$ is a function on x (only). The terms $u^1(x,y)\hspace{0.17em}$, $u^2(x,y)$ are considered to be quasi-periodical. Using (87)–(89), it can be written as

$$\begin{gathered} {\epsilon }_{\mathit{k}\mathit{h},\mathit{x}}(\mathit{u})=\frac{1}{2}\left(\frac{\partial {\mathit{u}}_{\mathit{k}}}{\partial {\mathit{x}}_{\mathit{h}}}+\frac{\partial {\mathit{u}}_{\mathit{h}}}{\partial {\mathit{x}}_{\mathit{k}}}\right)\hspace{1em} \\ =\frac{1}{2}\left(\frac{\partial u_k^o}{\partial x_h}+\frac{\partial u_h^o}{\partial x_k}\right)+\frac{\delta }{2}\left(\frac{\partial u_k^1}{\partial x_h}+\frac{\partial u_h^1}{\partial x_k}\right)+\frac{\delta }{2}\left(\frac{\partial u_k^2}{\partial y_h}+\frac{\partial u_h^2}{\partial y_k}\right)+....\hspace{1em} \\ ={\epsilon }_{kh,x}(u^o)+{\epsilon }_{kh,y}(u^1)+\delta \left[{\epsilon }_{kh,x}(u^1)+{\epsilon }_{kh,y}(u^2)\right]+{\delta }^2[\dots ].+\dots \hspace{1em};\hspace{1em}k,h=1,2,3 \end{gathered}$$

(90)

or

$${\epsilon }_{kh,x}(u)={\epsilon }_{kh}^o+\delta \hspace{0.17em}{\epsilon }_{kh}^1+\dots \hspace{1em};\hspace{1em}k,h=1,2,3$$

(91)

where

$${\epsilon }_{kh}^o={\epsilon }_{kh,x}(u^o)+{\epsilon }_{kh,y}(u^1)\hspace{1em}k,h=1,2,3$$

(92)

$${\epsilon }_{kh}^1=\left[{\epsilon }_{kh,x}(u^1)+{\epsilon }_{kh,y}(u^2)\right]\hspace{1em};\hspace{1em}k,h=1,2,3$$

(93)

Applying the linear Hooke’s law results in

$${\sigma }_{ij}^o=C_{ijkh}{\epsilon }_{kh}^o\hspace{1em},\hspace{1em}i,j,k,h=1,2,3$$

(94)

From Equation (94), it follows that

$$\frac{\partial \left(C_{ijkh}{\epsilon }_{kh}^o\right)}{\partial y_j}=0\hspace{1em},\hspace{1em}i,j,k,h=1,2,3$$

(95)

or

$$\frac{\partial \left[C_{ijkh}\left({\epsilon }_{kh,x}^{}(u^o)+{\epsilon }_{kh,y}^{}(u^1)\right)\right]}{\partial y_j}=0\hspace{1em},\hspace{1em}i,j,k,h=1,2,3$$

(96)

The terms ${\epsilon }_{kh,x}^{}(u^o)$ depend only on x. Equation (96) can be written as

$$-\frac{\partial \left[C_{ijkh}\hspace{0.17em}{\epsilon }_{kh,y}^{}(u^1)\right]}{\partial y_j}={\epsilon }_{kh,x}(u^o)\frac{\partial C_{ijkh}}{\partial y_j}\hspace{1em},\hspace{1em}i,j,k,h=1,2,3$$

(97)

introducing

$$u^1=w^{kh}{\epsilon }_{kh,x}^{}(u^o)+k(x)$$

(98)

Using Equations (97) and (98) with k(x) an arbitrary function on x, we can obtain

$${\epsilon }_{lm,y}^{}(u^1)={\epsilon }_{kh,x}^{}(u^o)\frac{1}{2}\left(\frac{\partial w_l^{kh}}{\partial y_m}+\frac{\partial w_m^{kh}}{\partial y_l}\right)={\epsilon }_{kh,x}^{}(u^o){\epsilon }_{lm,y}^{}(w^{kh})$$

(99)

Equation (99) becomes

$$-{\epsilon }_{kh,x}^{}(u^o)\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}(w^{kh})\right]}{\partial y_j}={\epsilon }_{kh,x}(u^o)\frac{\partial C_{ijkh}}{\partial y_j}\hspace{1em},\hspace{1em}i,j,k,h=1,2,3\hspace{1em}=1$$

(100)

Equation (100) is valid for any strain field ${\epsilon }_{kh,x}(u^o)$, so Equation (100) becomes

$$-\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}(w^{kh})\right]}{\partial y_j}=\frac{\partial C_{ijkh}}{\partial y_j}\hspace{1em},\hspace{1em}i,j,k,h=1,2,3$$

(101)

Using Green’s theorem, we obtain

$$\begin{array}{l}{\displaystyle {\int }_{\Gamma }\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}(w^{kh})\right]}{\partial y_j}{v}_{i}dV}+{\displaystyle {\int }_{\Gamma }{C}_{ijlm}{\epsilon }_{lm,y}(w^{kh})\frac{\partial v_i}{\partial y_j}dV}\\ ={\displaystyle {\int }_{\partial \Gamma }{u}_i{C}_{ijlm}{\epsilon }_{lm,y}(w^{kh})v_i^{}dS}=0\hspace{1em},\hspace{1em}i,j,k,h=1,2,3\end{array}$$

(102)

$$\begin{array}{l}{\displaystyle {\int }_{\Gamma }\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}(w^{kh})\right]}{\partial y_i}{v}_{j}dV}+{\displaystyle {\int }_{\Gamma }{C}_{ijlm}{\epsilon }_{lm,y}(w^{kh})\frac{\partial v_j}{\partial y_i}dV}=\\ ={\displaystyle {\int }_{\partial \Gamma }{u}_j{C}_{ijlm}{\epsilon }_{lm,y}(w^{kh})v_{j}dS}=0\hspace{1em},\hspace{1em}i,j,k,h=1,2,3\end{array}$$

(103)

In Equation (103), the indices i and j have been interchanged and the property $C_{ijlm}=C_{jilm}$ has been considered. From (102) and (103), it follows that

$$\begin{gathered} {\displaystyle {\int }_{\Gamma }\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}(w^{kh})\right]}{\partial y_j}{v}_{i}dV}+{\displaystyle {\int }_{\Gamma }\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}(w^{kh})\right]}{\partial y_i}{v}_{j}dV} \\ =2{\displaystyle {\int }_{\Gamma }{C}_{ijlm}{\epsilon }_{lm,y}(w^{kh})\frac{1}{2}\left(\frac{\partial v_i}{\partial y_j}+\frac{\partial v_j}{\partial y_i}\right)dV}=2{\displaystyle {\int }_{\Gamma }{C}_{ijlm}{\epsilon }_{ij,y}(v){\epsilon }_{lm,y}(w^{kh})dV},\hspace{1em}i,j,k,h=1,2,3 \end{gathered}$$

(104)

Because $C_{ijlm}=C_{jilm}$, multiplying Equation (104) by $v$ results in

$$\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}(w^{kh})\right]}{\partial y_j}{v}_i=\frac{\partial C_{ijkh}}{\partial y_j}{v}_i$$

(105)

Interchanging the indices i and j results in

$$\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}(w^{kh})\right]}{\partial y_i}{v}_j=\frac{\partial C_{ijkh}}{\partial y_i}{v}_j$$

(106)

The integration and addition of Equations (105) and (106) using Equation (104) offer

$${\displaystyle {\int }_{\Gamma }{C}_{ijlm}{\epsilon }_{ij,y}(v){\epsilon }_{lm,y}(w^{kh})dV}={\displaystyle {\int }_{\Gamma }\frac{\partial C_{ijkh}}{\partial y_j}{v}_{i}dV}$$

(107)

We must then find $w^{kh}$ in $V_y$ such that $\forall v\in V_y$, which verifies Equation (107). If $w^{kh}$ is obtained, then

$$\begin{gathered} {\sigma }_{ij}^o=C_{ijkh}\left[{\epsilon }_{kh,x}^{}(u^o)+{\epsilon }_{kh,y}^{}(u^1)\right] \\ =C_{ijkh}\left[{\epsilon }_{kh,x}^{}(u^o)+{\epsilon }_{kh,x}^{}(u^o)\hspace{0.17em}{\epsilon }_{lm,y}\left(w^{kh}\right)\right] \end{gathered}$$

(108)

By applying the average operator, it results in

$$\begin{gathered} \langle {\sigma }_{ij}^o\rangle =\langle C_{ijkh}{\epsilon }_{kh,x}(u^o)\rangle +\langle C_{ijkh}{\epsilon }_{kh,x}(u^o){\epsilon }_{lm,y}(u^{kh})\rangle \\ =\langle C_{ijkh}\rangle {\epsilon }_{kh,x}(u^o)+\langle C_{ijkh}{\epsilon }_{lm,y}(u^{kh})\rangle {\epsilon }_{kh,x}(u^o) \end{gathered}$$

(109)

This produces

$$\langle {\sigma }_{ij}^o\rangle =\left[\langle C_{ijkh}\rangle +\langle C_{ijkh}{\epsilon }_{lm,y}(u^{kh})\rangle \right]{\epsilon }_{kh,x}(u^o)$$

(110)

A comparison of Equation (110) can be made with

$$\langle {\sigma }_{ij}^o\rangle =C_{ijkh}^o{\tilde{\epsilon }}_{kh}(u^o)$$

(111)

and, if we denote ${\epsilon }_{kh,x}(u^o)\equiv {\tilde{\epsilon }}_{kh}(u^o)$, the homogenized coefficients can be obtained:

$$C_{ijkh}^o=\langle C_{ijkh}\rangle +\langle C_{ijkh}{\epsilon }_{lm,y}(u^{kh})\rangle$$

(112)

Therefore, there are two ways to obtain the homogenized coefficients:

• Using the local equations, the strain and stress field and the averages are determined, obtaining the homogenized coefficients;
• Using the variational formulation and determining the function $w^{kh}$ can also help us to determine the homogenized coefficients.

For the fiber-reinforced composite, there is a class of solutions ${\mathrm{w}}^{kh}$, with k,h = 1,2,3 satisfying

$$-\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}\left(w\right)\right]}{\partial y_j}=\frac{\partial C_{ijkh}}{\partial y_j}\hspace{1em},\hspace{1em}i=1,2,3$$

(113)

with the boundary conditions

$${w^{kh}|}_{\partial \Gamma }=0$$

(114)

and

$$\langle w^{kh}\rangle =0$$

(115)

If $(x_1,x_2,x_3)$ are the principal material axes, we state

$$\begin{gathered} C_{1111}=C_{11}; C_{2222}=C_{22}; C_{1122}=C_{1133}=C_{12}; C_{2211}=C_{3311}=C_{21} \\ {C}_{3322}=C_{2233}=C_{23}; C_{3333}=C_{33}; C_{4444}=\left(C_{22}-C_{23}\right)/2; C_{5555}=C_{44}; \\ {C}_{6666}=C_{44} \end{gathered}$$

(116)

The other components of $C_{ijkl}$ are zero (we work with a transversely isotropic material). The stress–strain relation becomes

$$\left\{\begin{array}{c}{\sigma }_{22}\\ {\sigma }_{33}\\ {\tau }_{23}\end{array}\right\}=\left[\begin{array}{ccc}{C}_{22}& C_{23}& 0\\ C_{23}& C_{33}& 0\\ 0& 0& \frac{C_{22}-C_{23}}{2}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{22}\\ {\epsilon }_{33}\\ {\gamma }_{23}\end{array}\right\}$$

(117)

or

$$\left\{\sigma \right\}=\left[C\right]{\left\{\epsilon \right\}}^o$$

(118)

The equilibrium conditions in Equation (71) are

$$\left[\begin{array}{ccc}\frac{\partial }{\partial y_2}& 0& \frac{\partial }{\partial y_3}\\ 0& \frac{\partial }{\partial y_3}& \frac{\partial }{\partial y_2}\end{array}\right]\left\{\begin{array}{c}{\sigma }_{22}\\ {\sigma }_{33}\\ {\tau }_{23}\end{array}\right\}=0$$

(119)

Equation (119) can be expressed in compact form:

$$[\partial ]\left\{\sigma \right\}=[\partial ]\left[C\right]{\left\{\epsilon \right\}}^o=0$$

(120)

Equation (92) becomes

$${\left\{\epsilon \right\}}^o={\left\{\epsilon \right\}}_{,x}^{u^o}+{\left\{\epsilon \right\}}_{,y}^{u^1}$$

(121)

and (120) becomes

$$\begin{gathered} -\left[\begin{array}{ccc}\frac{\partial }{\partial y_2}& 0& \frac{\partial }{\partial y_3}\\ 0& \frac{\partial }{\partial y_3}& \frac{\partial }{\partial y_2}\end{array}\right]\left[\begin{array}{ccc}{C}_{22}^{(\lambda )}& C_{23}^{(\lambda )}& 0\\ C_{23}^{(\lambda )}& C_{33}^{(\lambda )}& 0\\ 0& 0& \frac{C_{22}^{(\lambda )}-C_{23}^{(\lambda )}}{2}\end{array}\right]{\left\{\epsilon \right\}}_{,y}^{u^1} \\ =\left[\begin{array}{ccc}\frac{\partial }{\partial y_2}& 0& \frac{\partial }{\partial y_3}\\ 0& \frac{\partial }{\partial y_3}& \frac{\partial }{\partial y_2}\end{array}\right]\left[\begin{array}{ccc}{C}_{22}^{(\lambda )}& C_{23}^{(\lambda )}& 0\\ C_{23}^{(\lambda )}& C_{33}^{(\lambda )}& 0\\ 0& 0& \frac{C_{22}^{(\lambda )}-C_{23}^{(\lambda )}}{2}\end{array}\right]{\left\{\epsilon \right\}}_{,x}^{u^o} \end{gathered}$$

(122)

Considering the plane strain loading conditions, we obtain

$$\left[\begin{array}{ccc}\frac{\partial }{\partial y_2}& 0& \frac{\partial }{\partial y_3}\\ 0& \frac{\partial }{\partial y_3}& \frac{\partial }{\partial y_2}\end{array}\right]\left[\begin{array}{ccc}{C}_{22}^{(\lambda )}& C_{23}^{(\lambda )}& 0\\ C_{23}^{(\lambda )}& C_{33}^{(\lambda )}& 0\\ 0& 0& \frac{C_{22}^{(\lambda )}-C_{23}^{(\lambda )}}{2}\end{array}\right]{\left\{\epsilon \right\}}_{,y}^{u^1}=0$$

(123)

Using the determined functions ${\mathit{w}}^{kh}$, it can be deduced that

$${\left\{\epsilon \right\}}_{,y}^{u^1}=\left[\begin{array}{ccc}{\epsilon }_{22}\left(w^{22}\right)& {\epsilon }_{22}\left(w^{33}\right)& {\epsilon }_{22}\left(w^{23}\right)\\ {\epsilon }_{33}\left(w^{22}\right)& {\epsilon }_{33}\left(w^{33}\right)& {\epsilon }_{33}\left(w^{23}\right)\\ {\epsilon }_{23}\left(w^{22}\right)& {\epsilon }_{23}\left(w^{33}\right)& {\epsilon }_{23}\left(w^{23}\right)\end{array}\right]{\left\{\epsilon \right\}}_{,x}^{u^o}$$

(124)

or, in an alternative form,

$${\left\{\epsilon \right\}}_{,y}^{u^1}=\left[\begin{array}{ccc}\left\{\epsilon \left(w^{22}\right)\right\}& \left\{\epsilon \left(w^{33}\right)\right\}& \left\{\epsilon \left(w^{23}\right)\right\}\end{array}\right]{\left\{\epsilon \right\}}_{,x}^{u^o}$$

(125)

Additionally,

$$\left[\begin{array}{ccc}\frac{\partial }{\partial {\mathit{y}}_2}& 0& \frac{\partial }{\partial {\mathit{y}}_3}\\ 0& \frac{\partial }{\partial {\mathit{y}}_3}& \frac{\partial }{\partial {\mathit{y}}_2}\end{array}\right]\left[\begin{array}{ccc}{\mathit{C}}_{22}^{(\lambda )}& {\mathit{C}}_{23}^{(\lambda )}& 0\\ {\mathit{C}}_{23}^{(\lambda )}& {\mathit{C}}_{33}^{(\lambda )}& 0\\ 0& 0& \frac{{\mathit{C}}_{22}^{(\lambda )}-{\mathit{C}}_{23}^{(\lambda )}}{2}\end{array}\right]\left[\begin{array}{ccc}\left\{\epsilon \left(w^{22}\right)\right\}& \left\{\epsilon \left(w^{33}\right)\right\}& \left\{\epsilon \left(w^{23}\right)\right\}\end{array}\right]{\left\{\epsilon \right\}}_{,x}^{u^o}=0$$

(126)

$$\left[\begin{array}{ccc}\left\{\epsilon \left(w^{22}\right)\right\}& \left\{\epsilon \left(w^{33}\right)\right\}& \left\{\epsilon \left(w^{23}\right)\right\}\end{array}\right]{\left\{\epsilon \right\}}_{,x}^{u^o}=0$$

(127)

Equation (127) should remain valid for all ${\left\{\epsilon \right\}}_{,x}^{u^o}$. Equation (113) becomes

$$-\frac{C_{ijlm}^{(\lambda )}\partial {\epsilon }_{lm,y}^{(\lambda )}(w^{kh})}{\partial y_j}=0\hspace{1em},\hspace{1em}i=1,2,3$$

(128)

In we consider the case of plane strain $\mathit{i}=2,3$ and $j=2,3$,

$${\mathit{C}}_{22}^{(\lambda )}\frac{\partial {\epsilon }_{22}\left(w^{kh}\right)}{\partial {\mathit{y}}_2}+{\mathit{C}}_{23}^{(\lambda )}\frac{\partial {\epsilon }_{33}\left(w^{kh}\right)}{\partial {\mathit{y}}_2}+\frac{1}{2}\left[{\mathit{C}}_{22}^{(\lambda )}-{\mathit{C}}_{23}^{(\lambda )}\right]\frac{\partial {\epsilon }_{23}\left(w^{kh}\right)}{\partial {\mathit{y}}_3}=0$$

(129)

and

$$C_{23}^{(\lambda )}\frac{\partial {\epsilon }_{22}\left(w^{kh}\right)}{\partial y_3}+C_{22}^{(\lambda )}\frac{\partial {\epsilon }_{33}\left(w^{kh}\right)}{\partial y_3}+\frac{1}{2}\left[C_{22}^{(\lambda )}-C_{23}^{(\lambda )}\right]\frac{\partial {\epsilon }_{23}\left(w^{kh}\right)}{\partial y_2}=0$$

(130)

The solution is

$$w^{kh}=\{\begin{array}{c}{w}^{kh,(f)}\hspace{1em}for\hspace{1em}y\in V_f\\ w^{kh,(m)}\hspace{1em}for\hspace{1em}y\in V_m\end{array}$$

(131)

satisfying the boundary conditions

$${w^{kh,(f)}|}_{\partial \Gamma }=\hspace{1em}{w^{kh,(m)}|}_{\partial \Gamma }; {\sigma }_{ij}^{(f)}{n}_j={\sigma }_{ij}^{(m)}{n}_j$$

(132)

The boundary conditions for the RUC are: ${\mathit{u}}_i^{}={\alpha }_{ij}{\mathit{y}}_j$. It is possible to show that the average strain is: $\langle {\epsilon }_{ij}\rangle ={\overline{\epsilon }}_{ij}={\alpha }_{ij}$. Let us denote the displacement field by $w^{*}$ having the property ${w^{*}|}_{\partial \mathit{\Gamma }}={\mathit{u}\hspace{0.17em}|}_{\partial \mathit{\Gamma }}$ and ${\overline{\epsilon }}_{kh}\left(w^{*}\right)={\alpha }_{ij}$. Due to the existing symmetry in the distribution of the unit cell it can be concluded that $\langle w^{*}\rangle =0$. The field $w$ is introduced as

$$w=w^{*}-u$$

(133)

with the boundary conditions

$${w\hspace{0.17em}|}_{\partial \Gamma }=0$$

(134)

$$\langle w\rangle =\langle w^{*}-u\rangle =\langle w^{*}\rangle -\langle u\rangle =0$$

(135)

This function $\left(w\right)$ verifies the condition of zero average and value zero on the contour, and it verifies Equation (130). For the “quasi-periodical fields” $u{\hspace{0.17em}}_1$, it follows that

$$u_1={\alpha }_{22}\left\{\begin{array}{c}\frac{w_{22}^{*}}{{\alpha }_{22}}-y_2\\ 0\end{array}\right\}+{\alpha }_{33}\left\{\begin{array}{c}0\\ \frac{w_{22}^{*}}{{\alpha }_{22}}-y_2\end{array}\right\}$$

(136)

The strain field is

$${\epsilon }_{22}\left(w^{22}\right)=\frac{{\epsilon }_{22}\left(w^{*}\right)}{{\alpha }_{22}}-1\hspace{1em};\hspace{1em}{\epsilon }_{33}\left(w^{22}\right)=0$$

(137)

$${\epsilon }_{22}\left(w^{33}\right)=\frac{{\epsilon }_{22}\left(w^{*}\right)}{{\alpha }_{33}}-1\hspace{1em};\hspace{1em}{\epsilon }_{33}\left(w^{33}\right)=0$$

(138)

and

$${\epsilon }_{22}\left(u_1\right)={\epsilon }_{22}\left(w^{*}\right)-{\alpha }_{22}\hspace{1em};{\epsilon }_{33}\left(u_1\right)={\epsilon }_{33}\left(w^{*}\right)-{\alpha }_{33}\hspace{1em}$$

(139)

or

$${\overline{\epsilon }}_{22}\left(w^{22}\right)=\frac{{\overline{\epsilon }}_{22}\left(w^{*}\right)}{{\alpha }_{22}}-1=0\hspace{1em};\hspace{1em}{\overline{\epsilon }}_{33}\left(w^{33}\right)=\frac{{\overline{\epsilon }}_{33}\left(w^{*}\right)}{{\alpha }_{33}}-1=0\hspace{1em}$$

(140)

For the fiber-reinforced unidirectional composite, the homogenized coefficients can be obtained with the following relations:

$$\begin{gathered} C_{ijkh}^o=\langle C_{ijkh}\rangle +\langle C_{ijkh}{\epsilon }_{lm,y}\left(w^{kh}\right)\rangle \\ =\frac{1}{V}{\displaystyle {\int }_{\Gamma }{C}_{ijkh}}dV+\frac{1}{V}{\displaystyle {\int }_{\Gamma }{C}_{ijkh}{\epsilon }_{lm}\left(w^{*}\right)}dV \\ =\frac{1}{V}\left(C_{ijkh}^{(f)}{V}_f+C_{ijkh}^{(m)}{V}_m\right)+\frac{1}{V}\left(C_{ijlm}^{(f)}{\overline{\epsilon }}_{lm}^{(f)}\left(w\right)V_f+C_{ijlm}^{(m)}{\overline{\epsilon }}_{lm}^{(m)}\left(w\right)V_m\right) \end{gathered}$$

(141)

where

$${\alpha }_{lm}{\overline{\epsilon }}_{lm}^{(f)}\left(w\right)={\overline{\epsilon }}_{lm}^{(f)}\left(w^{*}\right)-{\alpha }_{lm}\hspace{1em};\hspace{1em}{\alpha }_{lm}{\overline{\epsilon }}_{lm}^{(m)}\left(w\right)={\overline{\epsilon }}_{lm}^{(m)}\left(w^{*}\right)-{\alpha }_{lm}$$

(142)

Thus, we have

$$C_{ijkh}^o=v_f{C}_{ijkh}^{(f)}+v_m{C}_{ijkh}^{(m)}+v_f{C}_{ijkh}^{(f)}\left[\frac{{\overline{\epsilon }}_{lm}^{(f)}\left(w*\right)}{{\alpha }_{lm}}-1\right]+v_m{C}_{ijkh}^{(m)}\left[\frac{{\overline{\epsilon }}_{lm}^{(m)}\left(w*\right)}{{\alpha }_{lm}}-1\right]$$

(143)

As a result (considering the plane strain loading conditions),

$$\begin{gathered} C_{22}^o=v_f{C}_{22}^{(f)}+v_m{C}_{22}^{(m)}+v_f{C}_{22}^{(f)}\left[\frac{{\overline{\epsilon }}_{22}^{(f)}\left(w*\right)}{{\alpha }_{22}}-1\right]+v_m{C}_{22}^{(m)}\left[\frac{{\overline{\epsilon }}_{22}^{(m)}\left(w*\right)}{{\alpha }_{22}}-1\right] \\ =v_f{C}_{22}^{(f)}\frac{{\overline{\epsilon }}_{22}^{(f)}\left(w*\right)}{{\alpha }_{22}}+v_m{C}_{22}^{(m)}\frac{{\overline{\epsilon }}_{22}^{(m)}\left(w*\right)}{{\alpha }_{22}} \end{gathered}$$

(144)

and

$$\begin{gathered} C_{23}^o=v_f{C}_{23}^{(f)}+v_m{C}_{23}^{(m)}+v_f{C}_{23}^{(f)}\left[\frac{{\overline{\epsilon }}_{33}^{(f)}\left(w*\right)}{{\alpha }_{33}}-1\right]+v_m{C}_{23}^{(m)}\left[\frac{{\overline{\epsilon }}_{33}^{(m)}\left(w*\right)}{{\alpha }_{33}}-1\right] \\ =v_f{C}_{23}^{(f)}\frac{{\overline{\epsilon }}_{33}^{(f)}\left(w*\right)}{{\alpha }_{33}}+v_m{C}_{23}^{(m)}\frac{{\overline{\epsilon }}_{33}^{(m)}\left(w*\right)}{{\alpha }_{33}} \end{gathered}$$

(145)

4. The Mori–Tanaka Model

4.1. Mathematical Model

In the following section, the mathematical model proposed by Mori and Tanaka is applied to obtain the engineering parameters that define Hooke’s law for a one-dimensional fiber-reinforced composite [34]. We consider an epoxy matrix with a visco-elastic response, reinforced with monotonous and parallel aligned carbon fibers that are uniformly distributed inside the resin (Figure 5). The resulting material has an orthotropic behavior. However, there are applications where the fibers are elliptical cylinders. These cylinders are randomly distributed, and the behavior of the material is a transverse isotropic.

The theory developed in [78] is applied in [28] for a reinforced material with continuous cylindrical fibers with an elliptic section. To solve this problem, Mori–Tanaka’s [34] mean-field theory is used. In [79,80], the two phases of the composite are two isotropic materials.

We consider a comparison material (CM). In the CM, there is a linear relation between the mean strain field $\epsilon$° and the mean stress field $\overline{\sigma }$:

$$\overline{\sigma }=C_m{\epsilon }^0$$

(146)

The average strain field in the RUC is ${\epsilon }^m={\epsilon }^0+\overline{\epsilon }$ and the mean stress field is ${\sigma }^m=\overline{\sigma }+\tilde{\sigma }$. This results in

$${\sigma }^m=\overline{\sigma }+\tilde{\sigma }=C_m({\epsilon }^0+\overline{\epsilon })$$

(147)

The mean strain fields in the fiber and in the matrix are differentiated through an additional term ${\epsilon }^{pt}$, hence ${\epsilon }^f={\epsilon }^m+{\epsilon }^{pt}={\epsilon }^0+\tilde{\epsilon }+{\epsilon }^{pt}$. In a similar way, the average stress field differs by the term ${\sigma }^{pt}$ and, therefore, ${\sigma }^f=\overline{\sigma }+\tilde{\sigma }+{\sigma }^{pt}$. The generalized Hooke law becomes

$${\sigma }^f=\overline{\sigma }+\tilde{\sigma }+{\sigma }^{pt}=C_f({\epsilon }^0+\tilde{\epsilon }+{\epsilon }^{pt})$$

(148)

or

$${\sigma }^f=\overline{\sigma }+\tilde{\sigma }+{\sigma }^{pt}=C_f({\epsilon }^0+\tilde{\epsilon }+{\epsilon }^{pt})=C_m({\epsilon }^0+\tilde{\epsilon }+{\epsilon }^{pt}-{\epsilon }^{\ast })$$

(149)

We introduce ${\epsilon }^{pt}$ in Equation (149).

$${\epsilon }^{pt}=P{\epsilon }^{\ast }$$

(150)

The Eshelby transformation tensor P from Equation (150) is presented in Appendix A (where P_ikjl = P_jikl = P_ijtk). The average stress in the whole RUC is

$$\begin{array}{l}\overline{\sigma }=v_f{\sigma }_f+v_m{\sigma }_m=v_f(\overline{\sigma }+\tilde{\sigma }+{\sigma }^{pt})+v_m(\overline{\sigma }+\tilde{\sigma })\\ =(v_f+v_m)\overline{\sigma }+(v_f+v_m)\tilde{\sigma }+v_f{\sigma }^{pt}=\overline{\sigma }+\tilde{\sigma }+v_f{\sigma }^{pt}\end{array}$$

(151)

which reduces to

$$\tilde{\sigma }=-v_f{\sigma }^{pt}$$

(152)

In a similar way, we can obtain

$$\overline{\epsilon }=-v_f({\epsilon }^{pt}-{\epsilon }^{\ast })=-v_f(P{\epsilon }^{\ast }-{\epsilon }^{\ast })=-v_f(P-I){\epsilon }^{\ast }$$

(153)

I denotes the unit tensor. Equations (147) and (149) yield

$$C_f\left[{\epsilon }^0-v_f\left(P-I\right){\epsilon }^{\ast }+P{\epsilon }^{\ast }\right]=C_m\left[{\epsilon }^0-v_f\left(P-I\right){\epsilon }^{\ast }+P{\epsilon }^{\ast }-{\epsilon }^{\ast }\right]$$

(154)

or

$$\left[C_f\left(-v_f\left(P-I\right)+P\right)+C_m\left(v_f\left(P-I\right)-P+I\right)\right]{\epsilon }^{\ast }+\left(C_f-C_m\right){\epsilon }^0=0$$

(155)

or

$$\left[C_f\left(v_{m}P+v_{f}I\right)-C_m{v}_m(P-I)\right]{\epsilon }^{\ast }+\left(C_f-C_m\right){\epsilon }^0=0$$

(156)

and

$$\left[v_m\left(C_f-C_m\right)P+v_f\left(C_f-C_m\right)+C_m\right]{\epsilon }^{\ast }+\left(C_f-C_m\right){\epsilon }^0=0$$

(157)

The final form is

$$\left[\left(C_f-C_m\right)\left(v_{m}P+v_{f}I\right)+C_m\right]{\epsilon }^{\ast }+\left(C_f-C_m\right){\epsilon }^0=0$$

(158)

This offers

$$\begin{array}{l}{\epsilon }_{11}^{\ast }=\frac{1}{A}\left(A_{11}{\epsilon }_{11}^o+A_{12}{\epsilon }_{22}^o+A_{13}{\epsilon }_{33}^o\right)\hspace{1em};\hspace{1em}\\ {\epsilon }_{11}^{\ast }=\frac{1}{A}\left(A_{21}{\epsilon }_{11}^o+A_{22}{\epsilon }_{22}^o+A_{23}{\epsilon }_{33}^o\right)\hspace{1em};\\ {\epsilon }_{11}^{\ast }=\frac{1}{A}\left(A_{31}{\epsilon }_{11}^o+A_{32}{\epsilon }_{22}^o+A_{33}{\epsilon }_{33}^o\right)\hspace{1em}.\end{array}$$

(159)

The coefficients $A_{ij}$ are presented in Appendix A [30]. The shear strain is [30]

$${\epsilon }_{12}^{\ast }=\frac{\left(G_{12,f}-G_m\right)}{\left(G_{12,f}-G_m\right)\left(2v_m{P}_{1212}+v_f\right)+G_m}{\epsilon }_{12}^0$$

(160)

$${\epsilon }_{23}^{\ast }=\frac{\left(G_{23,f}-G_m\right)}{\left(G_{23,f}-G_m\right)\left(2v_m{P}_{2323}+v_f\right)+G_m}{\epsilon }_{23}^0$$

(161)

$${\epsilon }_{31}^{\ast }=\frac{\left(G_{31,f}-G_m\right)}{\left(G_{31,f}-G_m\right)\left(2v_m{P}_{3131}+v_f\right)+G_m}{\epsilon }_{31}^0$$

(162)

Equations (158)–(162) can now be used to determine the elastic/viscoelastic parameters of a composite, which is considered as an orthotropic body. To compute the Young’s modulus $E_m$, the composite specimen is subjected to a pure traction ${\overline{\sigma }}_{11}$. This results in the following equation: ${\overline{\sigma }}_{11}=E_{11}{\overline{\epsilon }}_{11}$ and ${\overline{\sigma }}_{11}=E_m{\overline{\epsilon }}_{11}^0$; ${\overline{\epsilon }}_{22}^0={\overline{\epsilon }}_{33}^0=-{\nu }_m{\overline{\epsilon }}_{11}^0$.

Equation (158) produces

$$\begin{array}{ll}\hfill {\overline{\epsilon }}_{11}={\overline{\epsilon }}_{11}^0+& v_f{\overline{\epsilon }}_{11}^{\ast }={\overline{\epsilon }}_{11}^0+v_f\left(\frac{A_{11}}{A}{\overline{\epsilon }}_{11}^0+\frac{A_{12}}{A}{\overline{\epsilon }}_{22}^0+\frac{A_{13}}{A}{\overline{\epsilon }}_{33}^0\right)\hfill \\ & \hspace{0.17em}={\overline{\epsilon }}_{11}^0\left(1+v_f{a}_{11}\right)-v_f{a}_{12}{v}_m{\overline{\epsilon }}_{11}^0-v_f{a}_{13}{v}_m{\overline{\epsilon }}_{11}^0\hfill \\ & \hspace{0.17em}={\overline{\epsilon }}_{11}^0\left[1+v_f\left[a_{11}-v_m\left(a_{12}+a_{13}\right)\right]\right]\hfill \end{array}$$

(163)

Here, we show that $a_{ij}^{}=A_{ij}^{}/A$, $A_{ij}^{}$, and $A$ is presented in Appendix A; see rel. (A6). This results in

$$E_{11}=\frac{{\overline{\epsilon }}_{11}^0}{{\overline{\epsilon }}_{11}}{E}_m=\frac{E_m}{1+v_f\left[a_{11}-v_m\left(a_{12}+a_{13}\right)\right]}$$

(164)

For the other directions, in a similar way, we obtain the following equations:

$$E_{22}=\frac{{\overline{\epsilon }}_{22}^0}{{\overline{\epsilon }}_{22}}{E}_m=\frac{E_m}{1+v_f\left[a_{22}-v_m\left(a_{21}+a_{23}\right)\right]}$$

(165)

and

$$E_{33}=\frac{{\overline{\epsilon }}_{33}^0}{{\overline{\epsilon }}_{33}}{E}_m=\frac{E_m}{1+v_f\left[a_{33}-v_m\left(a_{31}+a_{32}\right)\right]}$$

(166)

Considering the shear moduli, we have

$${\overline{\sigma }}_{12}=2G_{12}{\overline{\epsilon }}_{12\hspace{0.17em}}\hspace{0.17em};\hspace{0.17em}{\overline{\sigma }}_{12}=2G_m{\overline{\epsilon }}_{12}^0\hspace{0.17em}$$

(167)

However,

$${\overline{\epsilon }}_{12}={\overline{\epsilon }}_{12}^0+v_f{\overline{\epsilon }}_{12}^{\ast }={\epsilon }_{12}^{\ast }-v_f\frac{G_{12,f}-G_m}{\left(G_{12,f}-G_m\right)\left(2v_m{P}_{1212}+v_f\right)+G_m}{\epsilon }_{12}^0$$

(168)

Using Equations (167) and (168) produces $G_{12}$:

$$G_{12}=G_m\left(1+\frac{v_f}{\frac{G_m}{G_{12,f}-G_m}+2v_m{P}_{1212}}\right)$$

(169)

In the same way, we obtain

$$G_{23}=G_m\left(1+\frac{v_f}{\frac{G_m}{G_{23,f}-G_m}+2v_m{P}_{2323}}\right)$$

(170)

and

$$G_{31}=G_m\left(1+\frac{v_f}{\frac{G_m}{G_{31,f}-G_m}+2v_m{P}_{3131}}\right)$$

(171)

The Poisson ratio is computed using the formulas

$${\overline{\epsilon }}_{22}=-v_m{\overline{\epsilon }}_{11}\hspace{0.17em};\hspace{0.17em}{\overline{\epsilon }}_{22}^0=\hspace{0.17em}{\overline{\epsilon }}_{33}^0=-v_m{\overline{\epsilon }}_{11}^0\hspace{0.17em}$$

(172)

Note that

$$\begin{array}{l}{\overline{\epsilon }}_{11}={\overline{\epsilon }}_{11}^0+v_f{\overline{\epsilon }}_{11}^{\ast }={\overline{\epsilon }}_{11}^0+v_f{a}_{11}{\overline{\epsilon }}_{11}^0+v_f{a}_{12}{\overline{\epsilon }}_{22}^0+v_f{a}_{13}{\overline{\epsilon }}_{33}^0=\\ ={\overline{\epsilon }}_{11}^0\left(1+v_f{a}_{11}\right)+v_f{a}_{12}{\overline{\epsilon }}_{22}^0+v_f{a}_{13}{\overline{\epsilon }}_{33}^0\end{array}$$

(173)

and

$${\overline{\epsilon }}_{22}={\overline{\epsilon }}_{22}^0+v_f{\overline{\epsilon }}_{22}^{\ast }=v_f{a}_{21}{\overline{\epsilon }}_{11}^0+{\overline{\epsilon }}_{22}^0\left(1+v_f{a}_{22}\right)+v_f{a}_{23}{\overline{\epsilon }}_{33}^0$$

(174)

or

$${\overline{\epsilon }}_{11}=\left[\left(1+v_f{a}_{11}\right)-v_f{a}_{12}{v}_m-v_f{a}_{13}{v}_m\right]{\overline{\epsilon }}_{11}^0=\left[v_f{a}_{21}-v_m\left(1+v_f{a}_{22}\right)-v_m{v}_f{a}_{23}\right]{\overline{\epsilon }}_{11}^0$$

(175)

Introducing Equation (174) into Equation (175) produces

$$v_{12}=-\frac{{\overline{\epsilon }}_{22}}{{\overline{\epsilon }}_{11}}=-\frac{v_f{a}_{21}-v_m\left(1+v_f{a}_{22}\right)-v_m{v}_f{a}_{23}}{1+v_f{a}_{11}-v_f{a}_{12}{v}_m-v_f{a}_{13}{v}_m}$$

(176)

which can be written as

$$v_{12}=\frac{v_m-v_f\left[a_{22}-v_m\left(a_{21}+a_{23}\right)\right]}{1+v_f\left[a_{11}-v_m\left(a_{12}+a_{13}\right)\right]}$$

(177)

In the same way, it produces

$$v_{23}=\frac{v_m-v_f\left[a_{22}-v_m\left(a_{21}+a_{23}\right)\right]}{1+v_f\left[a_{33}-v_m\left(a_7+a_8\right)\right]}$$

(178)

and

$$v_{31}=\frac{v_m-v_f\left[a_{33}-v_m\left(a_{31}+a_{32}\right)\right]}{1+v_f\left[a_{11}-v_m\left(a_{12}+a_{13}\right)\right]}$$

(179)

5. The Finite Element Method Used to Obtain the Creep Response

Recently, FEM has become the main method used for the study of elastic systems, as it is able to address a multitude of situations and types of materials, including composite materials [81]. Specialized problems are also studied, such as the influence of temperature on the stresses that appear in the analyzed structures [82]. In [83], a model is presented for the study of a composite reinforced with silicon carbide fibers. A similar model is addressed in [84]. Bodies with transverse isotropy were also studied, as in [11,85]. If we are dealing with microstructured systems, where a unit cell can be identified, the geometric symmetry allows the analysis to be conducted only on a quarter or half of the unit cell, on a unit previously defined as the “representative unit cell” (RUC). The unit cell model with finite elements is presented in Figure 6 and Figure 7; two models of a RUC that are used in various applications are also presented (Models 1 and 2).

The mechanical constants used in the application are

E_m = 4.14 GPa ; v_m = 0.22; E_f = 86.90 GPa; v_f = 0.34

(180)

The results of the analysis are shown in

Table 1. Average values of stress and strain (Case 1).

$\overline{\mathit{\sigma }}$	Fiber	Matrix	RUC	$\overline{\mathit{\epsilon }}$	Fiber	Matrix	RUC
${\overline{\sigma }}_{22}$	0.146 × 103	0.122 × 103	0.138 × 103	${\overline{\epsilon }}_{22}$	0.16 × 100	0.263 × 101	0.106 × 101
${\overline{\sigma }}_{33}$	0.71 × 100	−0.906 × 100	0.118 × 100	${\overline{\epsilon }}_{33}$	−0.445 × 10−1	−0.137 × 101	−0.529 × 100
${\overline{\sigma }}_{11}$	0.324 × 102	0.413 × 103	0.357 × 102	${\overline{\epsilon }}_{11}$	0.0 × 100	0.0 × 100	0.0 × 100
${\overline{\sigma }}_{23}$	0.194 × 10−4	0.412 × 102	0.359 × 10−4	${\overline{\epsilon }}_{23}$	0.539 × 10−7	0.447 × 10−5	0.167 × 10−5

,

Table 2. Computed values of elastic moduli (Case 1).

Modulus [MPa]	Matrix	Fiber	Average
$E_{11}$	4140.0	86,900.0	56,278.0
$E_{23}=E_{13}$	4140.0	86,899.0	12,741.0
${\nu }_1$	0.34	0.22	0.259
${\nu }_{23}$	0.34	0.22	0.475
$G_{23}$	1544.0	35,614.7	4318.2
$K_{23}$	4827.4	63,597.7	12,886.2

,

Table 3. Average values of stress and strain (Case 2).

$\overline{\mathit{\sigma }}$	Fiber	Matrix	RUC	$\overline{\mathit{\epsilon }}$	Fiber	Matrix	RUC
${\overline{\sigma }}_{22}$	0.147 × 103	0.122 × 103	0.138 × 103	${\overline{\epsilon }}_{22}$	0.134 × 100	0.183 × 101	0.757 × 100
${\overline{\sigma }}_{33}$	0.857 × 102	0.702 × 102	0.801 × 102	${\overline{\epsilon }}_{33}$	0.485 × 10−1	0.157 × 100	0.881 × 10−1
${\overline{\sigma }}_{11}$	0.512 × 102	0.653 × 102	0.564 × 102	${\overline{\epsilon }}_{11}$	0.0 × 100	0.0 × 100	0.0 × 100
${\overline{\sigma }}_{23}$	0.992 × 10−5	0.322 × 10−4	0.181 × 10−4	${\overline{\epsilon }}_{23}$	0.306 × 10−7	0.190 × 10−5	0.717 × 10−6

,

Table 4. Computed values of the elastic moduli (Case 2).

Modulus	Matrix	Fiber	Average
$E_{11}$	4140.0	86,900.0	56,278.0
$E_{23}=E_{13}$	4140.0	86,900.0	12,741.0
${\nu }_1$	0.34	0.22	0.259
${\nu }_{23}$	0.34	0.22	0.475
$G_{23}$	1544.0	35,614.8	4318.2
$K_{23}$	4827.4	63,597.8	12,886.2

,

Table 5. Average values of stress and strain (Case 3).

$\overline{\mathit{\sigma }}$	Fiber	Matrix	RUC	$\overline{\mathit{\epsilon }}$	Fiber	Matrix	RUC
${\overline{\sigma }}_{22}$	0.147 × 103	0.123 × 103	0.138 × 103	${\overline{\epsilon }}_{22}$	0.160 × 100	0.263 × 101	0.106 × 10−1
${\overline{\sigma }}_{33}$	0.644 × 100	−0.983 × 100	0.049 × 10−1	${\overline{\epsilon }}_{33}$	−0.446 × 10−1	−0.137 × 101	0.530 × 100
${\overline{\sigma }}_{11}$	0.324 × 102	0.414 × 102	0.357 × 102	${\overline{\epsilon }}_{11}$	0.0 × 100	0.0 × 100	0.0 × 100
${\overline{\sigma }}_{23}$	0.899 × 10−5	−0.979 × 10−4	−0.301 × 10−4	${\overline{\epsilon }}_{23}$	0.258 × 10−7	−0.638 × 10−5	0.232 × 10−5

and

Table 6. Computed values of the elastic moduli (Case 3).

Modulus	Matrix	Fiber	Average
$E_{11}$	4140.0	86,900.0	56,279.0
$E_{23}=E_{13}$	4140.0	86,899.0	12,754.0
${\nu }_1$	0.34	0.22	0.259
${\nu }_{23}$	0.34	0.22	0.475
$G_{23}$	1544.0	35,614.7	4322.2
$K_{23}$	4827.4	63,597.7	12,900.8

(in these tables, $\overline{\sigma }$ is the average stress and $\overline{\epsilon }$ is the average strain).

In this paper, we used a three-dimensional model to obtain the shear modulus and Poisson’s ratios in a plane perpendicular to $x_2{x}_3$.

A few of the foregoing models are listed in

Table 7. Finite element models and associated boundary conditions (BCs).

Case	Model	B.C. (x2 Direction)	B.C. (x3 Direction)
1	Model 1-a	px = 137.90 (MPa)	py = 0.00 (MPa)
2	Model 1-b	px = 137.90 (MPa)	py = 80.0 (MPa)
3	Model 2-a	px = 137.90 (MPa)	py = 0.00 (MPa)
4	Model 2-b	ux = 0.01 (mm)	uy = 0.01 (mm)
5	Model 2-c	ux = 0.01 (mm)	uy = 0.01 (mm)

, for which the results using finite element analysis are obtained.

There are some discrepancies between the present FE results and those presented in [49]. With respect to these discrepancies, the following verification should be considered. If the boundary condition for the FE model is taken as $u_i={\alpha }_{ij}{x}_j$ (where ${\alpha }_{ij}={\alpha }_{ji}$), the average strain should be equal to ${\overline{\epsilon }}_{ij}={\alpha }_{ij}$. This can be demonstrated as follows:

$${\overline{\epsilon }}_{ij}=\frac{1}{V}{\displaystyle {\int }_{\Gamma }{\epsilon }_{ij}dV=}\frac{1}{2V}{\displaystyle {\int }_{\Gamma }\left(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\right)dV}$$

(181)

By applying Green’s theorem, it follows that

$$\begin{gathered} {\overline{\epsilon }}_{ij}=\frac{1}{2V}{\displaystyle {\int }_{\partial \hspace{0.17em}\Gamma }\left(n_i{u}_j+n_j{u}_i\right)ds}= \\ =\frac{1}{2V}\left({\displaystyle {\int }_{\partial \hspace{0.17em}\Gamma }{n}_i{\alpha }_{jk}{x}_{k}ds}+{\displaystyle {\int }_{\partial \hspace{0.17em}\Gamma }{n}_j{\alpha }_{il}{x}_{l}ds}\right)= \\ =\frac{1}{2V}\left({\alpha }_{jk}{\displaystyle {\int }_{\partial \hspace{0.17em}\Gamma }{n}_i{x}_{k}ds}+{\alpha }_{il}{\displaystyle {\int }_{\partial \hspace{0.17em}\Gamma }{n}_j{x}_{l}ds}\right) \end{gathered}$$

(182)

or

$$\begin{gathered} {\overline{\epsilon }}_{ij}=\frac{1}{2V}\left({\alpha }_{jk}{\displaystyle {\int }_{\hspace{0.17em}\Gamma }\frac{\partial x_k}{\partial x_i}dV}+{\alpha }_{il}{\displaystyle {\int }_{\hspace{0.17em}\Gamma }\frac{\partial x_l}{\partial x_j}dV}\right)= \\ =\frac{1}{2V}\left({\alpha }_{jk}{\displaystyle {\int }_{\hspace{0.17em}\Gamma }{\delta }_{ki}dV}+{\alpha }_{il}{\displaystyle {\int }_{\hspace{0.17em}\Gamma }{\delta }_{lj}dV}\right)=\frac{1}{2V}\left({\alpha }_{ji}+{\alpha }_{ij}\right)={\alpha }_{ij} \end{gathered}$$

(183)

The discrepancy identified with the results of [49] can be attributed to the different type of finite elements used.

Therefore, we obtain average strains and stresses, viz., ${\overline{\sigma }}_{22},{\overline{\sigma }}_{33},{\overline{\sigma }}_{11},{\overline{\sigma }}_{23}={\overline{\tau }}_{23},{\overline{\epsilon }}_{22},$ ${\overline{\epsilon }}_{33},{\overline{\epsilon }}_{11},{\overline{\epsilon }}_{23}=1/2\hspace{0.17em}{\gamma }_{23}$. Using these values, it is now possible to obtain the mechanical constants of the studied composite [56]. To determine the longitudinal elastic modulus $E_{11}$, we use the well-known rule of mixture:

$$E_{11}=E_f{\nu }_f+E_m{\nu }_m$$

(184)

where

$${\nu }_f=\frac{A_f}{A}\hspace{1em};\hspace{1em}{\nu }_f=\frac{A_m}{A}\hspace{1em}$$

(185)

The following relations exist:

$$\begin{gathered} {\overline{\sigma }}_{22}=C_{22}{\overline{\epsilon }}_{22}+C_{23}{\overline{\epsilon }}_{33}; \\ {\overline{\sigma }}_{33}=C_{23}{\overline{\epsilon }}_{22}+C_{22}{\overline{\epsilon }}_{33}; \\ {\overline{\sigma }}_{11}=C_{12}\left({\overline{\epsilon }}_{22}+{\overline{\epsilon }}_{33}\right); \\ {\overline{\tau }}_{23}=C_{66}\hspace{0.17em}{\overline{\gamma }}_{23} \end{gathered}$$

(186)

from which results

$$\left[\begin{array}{cc}{\overline{\epsilon }}_{22}& {\overline{\epsilon }}_{33}\\ {\overline{\epsilon }}_{33}& {\overline{\epsilon }}_{22}\end{array}\right]\left\{\begin{array}{c}{C}_{22}\\ C_{23}\end{array}\right\}=\left\{\begin{array}{c}{\overline{\sigma }}_{22}\\ {\overline{\sigma }}_{33}\end{array}\right\}$$

(187)

and

$$\left\{\begin{array}{c}{C}_{22}\\ C_{23}\end{array}\right\}=\frac{1}{{\overline{\epsilon }}_{22}^2-{\overline{\epsilon }}_{33}^2}\left[\begin{array}{cc}{\overline{\epsilon }}_{22}& -{\overline{\epsilon }}_{33}\\ -{\overline{\epsilon }}_{33}& {\overline{\epsilon }}_{22}\end{array}\right]\left\{\begin{array}{c}{\overline{\sigma }}_{22}\\ {\overline{\sigma }}_{33}\end{array}\right\}$$

(188)

This results in the following:

$$C_{22}=\frac{{\overline{\sigma }}_{22}{\overline{\epsilon }}_{22}-{\overline{\sigma }}_{33}{\overline{\epsilon }}_{33}}{{\overline{\epsilon }}_{22}^2-{\overline{\epsilon }}_{33}^2}\hspace{1em};\hspace{1em}{C}_{23}=\frac{{\overline{\sigma }}_{33}{\overline{\epsilon }}_{22}-{\overline{\sigma }}_{22}{\overline{\epsilon }}_{33}}{{\overline{\epsilon }}_{22}^2-{\overline{\epsilon }}_{33}^2}\hspace{1em}$$

(189)

For $C_{12}$ and $C_{66}$,

$$C_{12}=\frac{{\overline{\sigma }}_{11}}{{\overline{\epsilon }}_{22}+{\overline{\epsilon }}_{33}}\hspace{1em};\hspace{1em}{C}_{66}=\frac{{\overline{\tau }}_{23}}{{\overline{\gamma }}_{23}}\hspace{1em}$$

(190)

To determine the bulk modulus $K_{23}$ is used in the relation:

$$K_{23}=\frac{C_{22}+C_{33}}{2}=\frac{{\overline{\sigma }}_{22}+{\overline{\sigma }}_{33}}{2\left({\overline{\epsilon }}_{22}+{\overline{\epsilon }}_{33}\right)}\hspace{1em}.$$

(191)

The longitudinal Poisson’s ratio is calculated via the following relation:

$${\nu }_1={\nu }_{21}={\nu }_{31}=\frac{1}{2}{\left(\frac{C_{11}-E_{11}}{K_{23}}\right)}^{1/2}=\frac{C_{12}}{C_{22}+C_{33}}=\frac{{\overline{\sigma }}_{11}}{\left({\overline{\sigma }}_{22}+{\overline{\sigma }}_{33}\right)}\hspace{1em}.$$

(192)

and the shear modulus

$$G_{23}=\frac{C_{22}-C_{33}}{2}=\frac{{\overline{\sigma }}_{22}-{\overline{\sigma }}_{33}}{2\left({\overline{\epsilon }}_{22}-{\overline{\epsilon }}_{33}\right)}\hspace{1em}.$$

(193)

or from

$$G_{23}=C_{66}=\frac{{\overline{\sigma }}_{23}}{2{\overline{\epsilon }}_{23}}\hspace{1em}.$$

(194)

By introducing the following parameter,

$$\psi =1+\frac{4{\nu }_1^2{K}_{23}}{E_{11}}\hspace{1em},$$

(195)

the transverse moduli and Poisson’s ratio are obtained using the following relations:

$$E_{22}=E_{33}=\frac{4G_{23}{K}_{23}}{K_{23}+\psi G_{23}}\hspace{1em}$$

(196)

and

$${\nu }_{23}=\frac{K_{23}-\psi G_{23}}{K_{23}+\psi G_{23}}\hspace{1em}$$

(197)

As such, the expressions for $E_{11},\hspace{0.17em}{E}_{22}=E_{33},\hspace{0.17em}{\nu }_{12}={\nu }_{13},\hspace{0.17em}{\nu }_{23},\hspace{0.17em}{G}_{23},K_{23}$ were determined. From

$$C_{22}+C_{33}=2K_{23}\hspace{1em}; C_{22}-C_{33}=2G_{23}\hspace{1em}$$

(198)

one may obtain

$$C_{22}=K_{23}+G_{23}\hspace{1em};\hspace{1em}{C}_{23}=K_{23}-G_{23}\hspace{1em}.$$

(199)

Recall that

$$C_{44}=G_1=G_{12}=G_{13}; C_{12}={\nu }_1\left(C_{22}+C_{23}\right)=2{\nu }_1{K}_{23}\hspace{1em}$$

(200)

and

$$C_{11}=E_{11}+\frac{2C_{12}^2}{C_{22}+C_{23}}=E_{11}+4{\nu }_1^2{K}_{23}=\psi E_{11}\hspace{1em}.$$

(201)

In a similar way, FEM was used to determine the average stresses and strains in a 3D elastic solid. These are:${\overline{\sigma }}_{11},{\overline{\sigma }}_{22},{\overline{\sigma }}_{33},{\overline{\sigma }}_{12}={\overline{\tau }}_{12},{\overline{\sigma }}_{23}={\overline{\tau }}_{23},{\overline{\sigma }}_{31}={\overline{\tau }}_{31},$ ${\overline{\epsilon }}_{11},{\overline{\epsilon }}_{22},{\overline{\epsilon }}_{33},{\overline{\epsilon }}_{12}=1/2\hspace{0.17em}{\gamma }_{12},$ ${\overline{\epsilon }}_{23}=1/2\hspace{0.17em}{\gamma }_{23},{\overline{\epsilon }}_{31}=1/2\hspace{0.17em}{\gamma }_{31}$. The general Hooke’s Law can be written as follows:

$$\begin{gathered} {\overline{\sigma }}_{11}=C_{11}{\overline{\epsilon }}_{11}+C_{12}{\overline{\epsilon }}_{22}+C_{12}{\overline{\epsilon }}_{33} \\ {\overline{\sigma }}_{22}=C_{12}{\overline{\epsilon }}_{11}+C_{22}{\overline{\epsilon }}_{22}+C_{23}{\overline{\epsilon }}_{33} \\ {\overline{\sigma }}_{33}=C_{12}{\overline{\epsilon }}_{11}+C_{23}{\overline{\epsilon }}_{22}+C_{22}{\overline{\epsilon }}_{33} \\ {\overline{\sigma }}_{23}={\overline{\tau }}_{23}=\left(C_{11}-C_{23}\right){\overline{\epsilon }}_{23}=\frac{1}{2}\left(C_{11}-C_{23}\right){\overline{\gamma }}_{23} \\ {\overline{\sigma }}_{31}={\overline{\tau }}_{31}=2C_{44}{\overline{\epsilon }}_{31}=C_{66}{\overline{\gamma }}_{31} \\ {\overline{\sigma }}_{12}={\overline{\tau }}_{12}=2C_{44}{\overline{\epsilon }}_{12}=C_{66}{\overline{\gamma }}_{12} \end{gathered}$$

(202)

From the last part of Equation (202), we can obtain

$$C_{44}=\frac{{\overline{\sigma }}_{12}}{2\hspace{0.17em}{\overline{\epsilon }}_{12}}=\frac{{\overline{\tau }}_{12}}{{\overline{\gamma }}_{12}}=G_{12}=G_{13}=G_1$$

(203)

Equation (202) yields

$${\overline{\sigma }}_{22}-{\overline{\sigma }}_{33}=\left(C_{22}-C_{33}\right)\left({\overline{\epsilon }}_{22}-{\overline{\epsilon }}_{33}\right)$$

(204)

The law of mixture offers us

$$E_{11}=E_f{v}_f+E_m{v}_m$$

(205)

Using Equation (205), one can replace the redundant fourth relation from Equation (201) with

$$E_{11}=C_{11}-\frac{2C_{12}^2}{C_{22}+C_{23}}$$

(206)

The addition of the second and third equations in Equation (202) yields

$${\overline{\sigma }}_{22}+{\overline{\sigma }}_{33}-\frac{2C_{12}{\overline{\epsilon }}_{11}}{{\overline{\epsilon }}_{22}+{\overline{\epsilon }}_{33}}=C_{22}+C_{23}$$

(207)

From the first part of Equation (202), one can show that

$$C_{11}=\frac{{\overline{\sigma }}_{11}-C_{12}\left({\overline{\epsilon }}_{22}+{\overline{\epsilon }}_{33}\right)}{{\overline{\epsilon }}_{11}}$$

(208)

The substitution of Equation (208) into $E_{11}$ yields

$$E_{11}=\frac{{\overline{\sigma }}_{11}}{{\overline{\epsilon }}_{11}}+C_{12}\frac{{\overline{\epsilon }}_{22}+{\overline{\epsilon }}_{33}}{{\overline{\epsilon }}_{11}}-\frac{2C_{12}^2\left({\overline{\epsilon }}_{22}+{\overline{\epsilon }}_{33}\right)}{\left({\overline{\sigma }}_{22}+{\overline{\sigma }}_{33}-2C_{12}{\overline{\epsilon }}_{11}\right)}$$

(209)

from which it is possible to compute $C_{12}$.

Figure 8 and Figure 9 present two creep curves for a composite carbon/epoxy at two different temperatures [56,57].

6. Conclusions

The method presented in this paper proves to be a calculus method suitable for obtaining the general mechanical constants of a multiphase composite material. The material constants required by designers are obtained using the average of the values obtained by applying FEM. The results obtained experimentally verified the models proposed by different researchers. The tests and measurements conducted here show a good concordance between the results obtained using the proposed models and the experimental verifications. Thus, FEM proves to be a powerful tool for determining the engineering constants of composite materials. Compared to the methods described in the other sections, this method proves to be a useful and relatively simple means of identifying the constitutive laws. The results were also applied to a study of the creep behavior of a composite material. This case is more complicated because, in the case of creep phenomena, the influences of temperature prove to be nonlinear. All the presented models can replace expensive methods of determining the engineering constants of a viscoelastic material by experimental measurements with calculation-based methods.

This review focuses on the behavior of unidirectional fiber composites.