Creep Response of Carbon-Fiber-Reinforced Composite Using Homogenization Method

Abstract

The homogenization theory, used for the study of differential equations with periodic coefficients, with a rapid variation, is used in the paper for the analysis of the creep phenomenon of composite materials, reinforced with fibers. Generally, a polymer composite having a matrix with a viscoelastic response manifests a creep behavior. A good knowledge of mechanical constants allows us to predict the time response under the action of a load, which is important in engineering. The homogenization method is used to determine the engineering constants for a composite reinforced with carbon fibers. The method is applied for the particular case of fiber-reinforced unidirectional composites to obtain the equations that finally offer the required values. The epoxy matrix Fibredux 6376C is reinforced with carbon fibers T800 and the thermoplastic specimens made by APC2 material is reinforced with carbon fibers of the type IM6. The experimental results give a good concordance with the theoretical predictions.

1. Introduction

Creep phenomenon in viscoelastic materials represents a permanent deformation under a mechanical stress. This phenomenon can occur after a long-term exposure of the material to a high stress and is generally time-dependent. The rate of deformation depends on time exposure, temperature, level of stress and properties of the material [1,2]. Temperature is an important factor that increases significantly the rate of deformation. The material can increase in length and in some engineering application this can be undesirable. The creep behavior of the material can occur usually near the melting point, but for some materials this phenomenon can be manifested at the room temperature. Polymers, including plastics and reinforced plastics can reveal significant creep behavior even at room temperatures. The creep phenomenon in composite materials has been studied in detail by various authors, mainly due to the fact that composite materials are currently very widely used and because, in general, the materials used as matrices are visco-elastic.

The creep is characterized by creep strain, defined as the slope of the creep strain-time curve. In the design phase of a project, if the creep behavior occurs, it is necessary to know the rate of deformation. This is obtained by creep tests, when the creep strain–time curve is performed and can be implemented directly in the project [2,3]. A brief history of the domain and the basis of the creep studies are presented in [4,5,6].

Some results in the study of the creep phenomenon are presented below, to present the general framework of research in the field. In [7], an accelerated characterization scheme for the determination of the viscoelastic behavior of composites is proposed. This kind of approach leads to a reduced number of experimental measurements. Short-term observation done on a fiber reinforced composite at different temperatures allow the prediction of the long-term response. The studies use the time–temperature superposition principle (TTSP) [8,9,10].

Other studies concerning the nonlinear viscoelastic behavior of fiber reinforced composites has been made by Schaffer and Adams [11]. Their model uses a finite element micromechanics procedure, simplifying the analysis while using the symmetry pattern of such composite. The method can be applied not only for a bi-phasic material but also for a complex topology of a multiphasic composite.

Schapery [12,13] made a study of nonlinear viscoelastic materials, where thermal- and moisture-dilatational effects were considered through initial stresses. Mohan and Adams [14] applied this approach to study the graphite- and glass-reinforced composites and obtained a correlation between the method proposed by [11] and those proposed by [12,13].

Findley et al [15,16] proposed an empirical model to study the nonlinear creep compliance of a different composite, based on the creep power law. The method is well suited to numerical procedures. A nonlinear viscoelastic model using the Findley approach is applied to study of creep behavior of graphite/epoxy composites, by Dillard et al [17,18]. These theories were verified by measurements by Walrath [19].

An important stage in these analyses is to determine the constitutive law. These represents difficulties, as the material is inhomogeneous due to the reinforcing elements. The problem is to benefit from methods for calculating engineering constants which, for practical applications, must be simple to use. To achieve this goal, different averaging methods were utilized, using theorems from the mechanics of continuous media. Some results are summarized in the following.

Hashin et al, in series of investigations [20,21,22], used a method to pass from elastic to time-dependent behavior using a variational principle. It is assumed that the fiber volume ratio remain constant in different part of the composite despite the variation of diameters of fibers. A composite is considered to be a collection of cylindrical fibers with differing sizes. The assumption is that the strain energy of a representative element (a fiber surrounded by matrix) is equal to that of a homogeneous element.

Zhao et al [23] obtained results concerning the engineering constants which characterize an orthotropic and a transversely isotropic composite. Hill [24,25,26] determined the upper and lower bounds to estimate the elastic moduli of a fiber reinforced composite. Other researcher extended the Mori–Tanaka method [27,28] to study the viscoelastic response of an isotropic composite.

Significant results concerning the overall viscoelastic behavior of a composite reinforced with parallel fibers were obtained by Aboudi [29,30] in his micromechanical approach of one-dimensional composites. Another way to study this type of problem is that of the finite element analysis, with which a good agreement with the micromechanical model was obtained by Schaffer and Adams [11]. In other work, a micromechanical model to study the creep response of a composite carbon fiber reinforced composite is proposed [31]. To verify the model, experimental tests at different load and temperatures offer the creep curves.

The study of the mechanical properties of biphasic materials is made in numerous papers in the last few years [32,33,34,35,36,37,38,39], works that lead to the subject of the paper; namely, the use of the homogenization method for the calculation of elastic constants for such materials.

Over the last few decades, different mathematical estimation methods are used to determine the viscoelastic/elastic/plastic properties of a composite material. The mean field homogenization procedure is a precise method used in determining the mechanical properties of different types of composites, being able to analyze a wide class of materials, from unidirectional composites to textile composites, short fiber reinforced composites or nanotube-reinforced composites. The application of this optimization method for some practical cases is presented in [40], making a micro and meso-mechanical analysis of these types of materials. The estimation of the mechanical properties of the obtained material must be done using specific calculation methods. The homogenization method proves to be a very suitable method for solving this problem. In [41], this method is applied and the obtained results are compared with experimental results.

Other works refer to the theoretical aspect of the homogenization method, which involves the transition from a periodic structure of a composite material to a homogeneous structure [42]. Proposed analytical and numerical models have been used to determine the mechanical constants of a composite material. The obtained results were verified by the traction/compression and shear test. The experiments indicated that both obtained theoretical and numerical prediction values are in agreement with the results of experimental verifications confirming the validity of the methodology in providing a reliable reference for the structural design of pultruded fiber-reinforced polymer composite structures (FRP) [43].

A new method of homogenization of elasto-viscoplastic composites is developed in [44]. If Euler’s integration algorithm is used, the nonlinear ordinary differential equations in the mechanics of the elastic solid that describe the phase deformation of a composite can be discretized for a numerical approach. A total of three classical methods of approach can be used in this case: direct, secant and tangent. The interaction between the different phases of the composite material and the matrix in which these phases are incorporated is achieved by defining a new second-order tensor. In this way, the problems of homogenization of heterogeneous elasto-plastic and elasto-viscoplastic materials can be unified [45]. An example in which the engineering constants can be obtained analytically is presented and it is found that the results coincide with the values known in the literature [46]. During the research, experimental measurements are performed to verify the computed values. It is found that, for most of the tests performed for the specimens of carbon fiber reinforced composites, the results are very close to the experimentally determined values.

The research carried out at present deals mainly with the refinement of the proposed methods and with applications to different situations that may be encountered in practice [47,48,49,50,51,52,53,54]. Related methods of studying the problem are presented in [55,56]. Interest in the application of the theory of homogenization for this type of problem exists, but specific models for the application of homogenization are rarely found in the literature. The development of models with experimentally verified results is useful for researchers and manufacturers of such composites.

In the paper, the homogenization theory previously presented is applied to obtain the mechanical constants of a composite material. Thus, the general theory is applied for the special case of fiber-reinforced unidirectional composites. The equations are obtained that finally offer the required values. In general, obtaining these values involves numerical procedures, which was applied for the materials studied in the paper.

2. Homogenized Model of Carbon Fiber Composite

Generally, in mathematics, homogenization represents a method for the study of partial differential equations having rapidly oscillating coefficients. This becomes important for some types of continuous materials such composite materials with a periodical structure. In this case the coefficients of differential equations describing the response of material a periodical. Such material is convenient to treat as homogeneous, despite the fact that the material possess a particular microstructure. The homogenization procedure replaces an equation having high oscillatory coefficients with an equation with constant coefficients and represents an extension of the continuum concept to a new class of materials: materials with microstructure. The foundations of homogenization method are presented in a series of papers, i.e., [57,58,59,60,61,62]. In the following, this theory is applied for a fiber reinforced composite.

2.1. Macroscopic Equations

The stress field ${\sigma }^{\delta }$ for the case when the entire material is subdivided into repeating unit cells of dimension $\delta$ is:

$$\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}$$

(1)

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

The field of displacements $u^{\delta }$ should satisfy the contour conditions:

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

(2)

on the other hand, the stresses should satisfy the boundary conditions:

$$\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}$$

(3)

on the contour ${\partial }_2\Omega$,(${\partial }_1\Omega \cup {\partial }_2\Omega =\partial \Omega$). The stress–strain relation in the local domain can be written according to the Hook’s Law as:

$${\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 }$$

(4)

or in the compact form as:

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

(5)

The above elasticity matrix is semi-positive definite with coefficient $C_{ijkh}$ being functions of x. One can write:

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

(6)

for $\alpha >0$ and $\forall x_{ij},x_{kh}\in R$. In the analysis of composite materials, $C_{ijkh}(x)$ are periodical function of x, with the period having the dimension $\delta$ of the unit cell. Having introduced the notation $y=x/\delta$, these coefficients become a function of y in domain $\mathsf{\Gamma }$ of the unit periodic cell:

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

(7)

Let us illustrate the applicability of the above system of equations by introducing a model which can be used to obtain the solution of Equation (1) written earlier. Let further the stress field in the unit cell be expressed in the following form:

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

(8)

Note that the dependence of stress on y is “quasi-periodical”.

By applying Equation (8) to the equilibrium Equation (1) it follows that:

$${\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+\mathrm{.....}=f_i(x)$$

(9)

Note that in the derivation of the above equation, use is made of the following property. For a function f which depends on x and y where y itself is a function of x, one can write:

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

(10)

but since $y=x/\delta$ and thus $dy=dx/\delta$; therefore:

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

(11)

Identification of the terms with the coefficients of ${\delta }^{-1}$ in Equation (9) yields:

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

(12)

which is called the “local equation”. In this equation, ${\sigma }_{ij}^o$ is a function of y and x. However, since the dependence of ${\sigma }_{ij}^o$ on x is very weak, x can be taken as a parameter, constant in the domain $\mathsf{\Gamma }$. On the other hand, the function ${\sigma }_{ij}^o$ depends on y in a periodical way. This in turn means that one can consider this equation only within the domain of a cell $\mathsf{\Gamma }$. By solving this differential system “the homogenized elastic matrix” can be established.

Similarly, if one identifies the corresponding terms of ${\delta }^o$, following equations are obtained:

$$\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$$

(13)

These are referred to as “microscopic equations” and they contain both variations on x as well as on y. By applying the average operator to this equation, it follows that:

$$\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$$

(14)

The average operator represents the average of the considered size over the domain V. The main idea here is to “average” the variation on y which is considered to be very small in comparison with the variation on x in the domain $\mathsf{\Gamma }$. On the other hand:

$$〈\frac{\partial {\sigma }_{ij}^1}{\partial y_j}〉=\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$$

(15)

Note that ${\sigma }_{ij}^1$ take equal values on the corresponding points of the boundary of the cell $\mathsf{\Gamma }$ (following from the property of periodicity). At the same points, $n_j$ takes opposite values as a results of which:

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

(16)

By integration of this equation, the homogenized displacement field $u^o$ is obtained for the entire domain Ω.

2.2. Evaluation of the Homogenized Coefficients

In order to evaluate the homogenized coefficients that are functions of the microstructure in a heterogeneous medium, following notations are introduced:

$${\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$$

(17)

$${\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$$

(18)

Let the displacement field be approximated by:

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

(19)

In this equation $u^o(x)$ depends only on x and terms $u^1(x,y)\hspace{0.17em}\delta$ and $u^2(x,y)$ are taken as quasi-periodical. Now using the above notations, one may write:

$$\begin{gathered} {\epsilon }_{kh,x}(u)=\frac{1}{2}\left(\frac{\partial u_k}{\partial x_h}+\frac{\partial u_h}{\partial x_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)+\mathrm{....}\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\left[\dots \right].+\dots \hspace{1em};\hspace{1em}k,h=1,2,3 \end{gathered}$$

(20)

which can be simplified to:

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

(21)

where:

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

(22)

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

(23)

In other words, the infinitesimal term $u^1(x,y)\delta$ represents a finite component of ${\epsilon }_{kh}$, the term that should be taken into account when applying the Hooke’s law:

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

(24)

From the local equation in (24), 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$$

(25)

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$$

(26)

The terms ${\epsilon }_{kh,x}^{}(u^o)$ depend only on x and thus may be considered constant for the current problem. The last of the above equations is therefore written in the alternative form 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$$

(27)

By substituting:

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

(28)

In Equation (27) with k(x) being an arbitrary function on x, it can be shown that:

$${\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})$$

(29)

Therefore, the differential Equation (29) 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$$

(30)

This relation should remain valid for any strain field ${\epsilon }_{kh,x}(u^o)$. Equation (30) can be simplified to:

$$-\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$$

(31)

Let us now make use of the relations:

$${\displaystyle {\int }_{\mathsf{\Gamma }}\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}(w^{kh})\right]}{\partial y_j}{v}_{i}dV}+{\displaystyle {\int }_{\mathsf{\Gamma }}{C}_{ijlm}{\epsilon }_{lm,y}(w^{kh})\frac{\partial v_i}{\partial y_j}dV}={\displaystyle {\int }_{\partial \mathsf{\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$$

(32)

$${\displaystyle {\int }_{\mathsf{\Gamma }}\frac{\partial \left[C_{ijlm}{\epsilon }_{lm,y}(w^{kh})\right]}{\partial y_i}{v}_{j}dV}+{\displaystyle {\int }_{\mathsf{\Gamma }}{C}_{ijlm}{\epsilon }_{lm,y}(w^{kh})\frac{\partial v_j}{\partial y_i}dV}={\displaystyle {\int }_{\partial \mathsf{\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$$

(33)

Note that in the second equation, the indices i and j have been interchanged and the property $C_{ijlm}=C_{jilm}$ has been used. It should be mentioned that here the Green’s theorem was applied while making use of the property that the functions $u_i,C_{ijlm},{\epsilon }_{lm}$ and $v_i$ are Y-periodical. From these two relations, one can write:

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

(34)

Now, multiplying both sides of Equation (34) by $v$ while taking into consideration the property $C_{ijlm}=C_{jilm}$, it follows that:

$$\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.$$

(35)

Interchanging the indices i and j yields:

$$\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$$

(36)

Integration and addition of the above two relations, while making use of the Equation (34) leads to:

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

(37)

In variational formulation, the problem is now to find $w^{kh}$ in $V_y$ such that $\forall v\in V_y$, the previous relation holds. When $w^{kh}$ is determined, 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}$$

(38)

By applying the average operator, it follows that:

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

(39)

Furthermore:

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

(40)

Comparison of this equation with:

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

(41)

while making use of the notation ${\epsilon }_{kh,x}(u^o)\equiv {\tilde{\epsilon }}_{kh}(u^o)$, yields the homogenized coefficients:

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

(42)

This means that in order to determine the homogenized coefficient, one has to calculate the Y-periodic function $w^{kh}$ with $〈w^{kh}〉=0$ which verifies the differential Equation (30).

In summary, there are two ways to evaluate the homogenized coefficients:

• If one starts from local equations, it is possible to determine the strain and stress field. By using the averages, the homogenized coefficients can be evaluated;
• One can use the variational formulation and find the function ${\mathrm{w}}^{kh}$ which allows computation of the homogenized coefficients.

For illustrative purpose, the above theory is applied to calculate the effective Young’s modulus of a composite in a one dimensional case (Appendix A).

2.3. Evaluation of Homogenized Coefficients for FRP

In the case of a fiber reinforced composite, there are two distinct phases present. These phases are the fiber and matrix materials denoted by f and m, respectively. Let $v_f=V_f/V$ and $v_m=V_m/V$ represent the volume fractions of the fiber and matrix in the composite where $V=V_f+V_m$ is the total volume of the representative unit cell. Due to the existing regular and periodic pattern of the fiber packing, it seems appropriate to limit the analysis to a repeating unit cell. Conforming to the previously described theory, there exists a class of solution ${\mathrm{w}}^{kh}$, with k,h = 1,2,3 which satisfies the differential equations:

$$-\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$$

(43)

with the boundary conditions:

$${\left.w^{kh}\right|}_{\partial \mathsf{\Gamma }}=0$$

(44)

and the supplementary conditions:

$$〈{\mathrm{w}}^{kh}〉=0$$

(45)

Let us take $(x_1,x_2,x_3)$ as the principal material axis, then, for a transversely isotropic material:

$$\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}$$

(46)

and all other components of $C_{ijkl}$ are identically zero. Under the plain strain loading condition, following stress-strain relations can be written:

$$\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\}$$

(47)

or alternatively as:

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

(48)

The equilibrium conditions written previously become:

$$\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$$

(49)

Or in a compact form:

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

(50)

The Equation (50) can then be written:

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

(51)

In this case, the equilibrium equations become:

$$\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}$$

(52)

The term ${\left\{\epsilon \right\}}_{,x}^{u^o}$ does not depend on y. The coefficients $C_{ij}^{(\lambda )}$ are constants in the two phases. It should be mentioned again that in the above equations and those which follow, $\left(\lambda \right)$ represents both fiber $\left(\lambda =f\right)$ and matrix $\left(\lambda =m\right)$ constituents in the unit cell. In the case of plane strain loading conditions and using the matrix notation, it follows that:

$$\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$$

(53)

Now, when the functions $w^{kh}$ are determined for the above conditions, one may write:

$${\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}$$

(54)

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}$$

(55)

Additionally:

$$\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[\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$$

(56)

$$\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$$

(57)

This relation should remain valid for all ${\left\{\epsilon \right\}}_{,x}^{u^o}$. With constant coefficients in both phases namely in the fiber and matrix, Equation (43) becomes:

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

(58)

For plane strain $i=2,3$ and $j=2,3$, therefore:

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

(59)

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$$

(60)

It should be noted that the above equations do not depend on k and j. This means that if one determines for example $w^{22}$, the remaining functions $w^{33}=w^{23}=w^{32}=w$ are immediately obtained. The solution of the above differential equation is composed of:

$$w^{kh}=\left\{\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}\right.$$

(61)

which should satisfy the boundary conditions:

$${\left.w^{kh,(f)}\right|}_{\partial \mathsf{\Gamma }}=\hspace{1em}{\left.w^{kh,(m)}\right|}_{\partial \mathsf{\Gamma }}$$

(62)

and:

$${\sigma }_{ij}^{(f)}{n}_j={\sigma }_{ij}^{(m)}{n}_j$$

(63)

These correspond to the continuity conditions of the displacements and of the stresses at the boundary $\partial \mathsf{\Gamma }$ of the two phases.

The following problem will now be considered: Let the repeating periodic cell be subjected to the boundary condition: $u_i^{}={\alpha }_{ij}{y}_j$. If the material is homogeneous, the average strain in the materials can be shown to be: $〈{\epsilon }_{ij}〉={\overline{\epsilon }}_{ij}={\alpha }_{ij}$. This will be demonstrated later in the chapter. Now, let us solve the problem of the elasticity where a periodical cell is subjected to the condition indicated above. Furthermore, let us denote the displacement field by $w^{*}$ with the property ${\left.w^{*}\right|}_{\partial \mathsf{\Gamma }}={\left.u\hspace{0.17em}\right|}_{\partial \mathsf{\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 $〈w^{*}〉=0$. Let the field $\mathrm{w}$ also be introduced as:

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

(64)

Corresponding to the conditions:

$${\left.w\hspace{0.17em}\right|}_{\partial \mathsf{\Gamma }}=0$$

(65)

and:

$$〈w〉=〈w^{*}-u〉=〈w^{*}〉-〈u〉=0$$

(66)

This function $\left(w\right)$ verifies the condition of zero average and has the value zero on the contour and it is also a verification of Equation (60). 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\}$$

(67)

The strain fields for this displacement field are:

$${\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$$

(68)

$${\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$$

(69)

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}$$

(70)

which can be written in a simplified form as:

$${\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}$$

(71)

Recall that for fiber reinforced unidirectional composite $C_{ijkh}^{(f)}$ and $C_{ijkh}^{(m)}$ represent the material coefficients for the fiber and matrix, respectively. The homogenized coefficients can be obtained using the following formula:

$$\begin{gathered} C_{ijkh}^o=〈C_{ijkh}〉+〈C_{ijkh}{\epsilon }_{lm,y}\left(w^{kh}\right)〉 \\ =\frac{1}{V}{\displaystyle {\int }_{\mathsf{\Gamma }}{C}_{ijkh}}dV+\frac{1}{V}{\displaystyle {\int }_{\mathsf{\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}$$

(72)

with:

$${\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}$$

(73)

Thus, one can show that:

$$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]$$

(74)

and for the plane strain loading conditions, it yields:

$$\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}$$

(75)

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}$$

(76)

Now, in order to determine the coefficients $C_{ij}^o$ in the above equations, it is necessary to obtain the field of strain and finally compute the average values of strain in both phases of the composite. This is analytically feasible for very simple problems which are usually presented in literature. However, in the case of fiber reinforced composites, analytical computation of the strain field is a complex and tedious task. Nevertheless, by using the method of finite elements, evaluation of the field if strain and/or stress becomes a relatively easy process.

3. Experimental Creep Response of Fiber Reinforced Composite

If the engineering constants for a material have been determined, the creep behavior of the materials can be studied. The equipment belongs to the mechanical testing laboratory of the Technical University of Munich.

In order to verify the model proposed in the previous section, experiments on the behavior of carbon fiber reinforced composite have been made (Figure 1). The variable parameters are the stress acting on the test specimen and the temperature. The results are presented in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9.

To perform the experiment commercially available composites have been used. The epoxy material used is Fibredux 6376C, (manufactured by Novartis AG, Basel, Switzerland). This material is reinforced with carbon fibers T800 (made Toray International, Inc, Tokyo, Japan). The thermoplastic specimens was made by APC2 material. In this case the PEEK matrix (in a semicrystalline state) has a glass transition to a temperature of 145 °C. This material is reinforced with carbon fibers (Ashland Inc, Lexington, KY, USA) of the type IM6.

The dimensions of the specimens are: length of 150 mm, width of 10 mm and thickness of 1 mm. All the specimens, made by different materials have the same dimensions. The specimens have been stored in desiccants filled chamber in the period before testing (to protect the specimens from humidity—the relative humidity in this chamber was 35%).

4. Conclusions and Discussion

It has been shown in the paper that a very good approximation of the creep response of a material can be obtained if the homogenization theory is used to determine the mechanical constants of a fiber-reinforced composite material. In this way, the homogenization theory becomes a useful method for the study and determination of the mechanical properties of such a material. Inside the paper, two materials were studied; namely, Carbon/Epoxy and Carbon/PEEK at different temperatures and different loads. Each time very good values were obtained, close to the values obtained experimentally. Therefore, homogenization theory is proving to be a useful tool for carbon fiber composites.

The influence of temperature on the behavior of the material proves to be nonlinear. As a result, in practice, special attention must be paid to a good modeling of the mechanical behavior of the material. The present study solves the problem which involves unidirectional composites with elastic but transversal isotropic fibers. This is important when we use carbon fiber, that exhibit this property (in contrast to fiberglass).

As a conclusion, the homogenization method proves to be a very powerful and accurate tool in estimating the creep response of a composite material. The measurements proved that the theoretical results obtained using the homogenization method fit very well with the results obtained from some experimental measurements.

The use of the method in the mechanical identification of the material represents an effective tool for determining the time dependent overall response of any unidirectional composite under various loading conditions with reasonable accuracy.