Calculation of Homogenized Mechanical Coefficients of Fiber-Reinforced Composite Using Finite Element Method

Abstract

Determining the mechanical properties of a composite material represents an important stage in its design and is generally a complicated operation. These values are influenced by the topology and geometry of the resulting composite and the values of the elastic constants of the components. Due to the importance of this subject and the increasing use of composite materials, different calculation methods have been developed over the last fifty years. Some of the methods are theoretical, with results that are difficult to apply in practice due to difficulties related to numerical calculation. In the current paper, using theoretical results offered by the homogenization theory, values of engineering elastic constants are obtained. The finite element method (FEM) is used to determine the stress and strain field required in these calculations; this is an extremely powerful and verified calculation tool for the case of a material with any type of structure and geometry. In order to minimize errors, the paper proposes the method of least squares, a mathematical method that provides the best estimate for the set of values obtained by calculating FEM. It is useful to consider as many load cases as possible to obtain the best estimates. The elastic constants for a transversely isotropic material (composite reinforced with cylindrical fibers) are thus determined for a real case.

1. Introduction

Almost all industrial branches have had an explosive growth in the last few decades in the use of new and composite materials. It has been practically demonstrated that they represent in many cases a better alternative to the use of classic materials. This has determined, as a consequence, the development of manufacturing technologies for obtaining composite materials and structures made from these materials for industrial, medical, transport, military, aerospace, etc., applications. One of the major applications of composites has turned out to be structural applications. Structures made of composites have greater resistance and less weight. The existence of these materials has determined numerous studies that have dealt with the properties of these materials, their behavior in different conditions and the possibilities of improving properties needed for practical applications. A composite is composed of two or more phases of material, one playing the role of matrix and the other of reinforcing elements. Fiber-reinforced polymers stand out among matrix materials. The properties of these materials are defined by the matrix (usually polyester resin) and the reinforcing materials (fibers) [1].

In designing and manufacturing a structure that has composite materials in its composition, knowing the engineering constants of the materials is a major requirement. This necessity has led to numerous studies in this direction with the aim of providing the designer with efficient and fast calculation methods. This objective (knowing the mechanical properties of a composite) remains important in the current context of the continuous development of the use of composite materials. Up to now, a systematization of this research has not yet been carried out. Experimental measurements represent the safest method to obtain useful and credible results. Unfortunately, experimental methods are time-consuming and involve large expenses. As a result, the development of methods for the rapid estimation of elastic constants remains a current and continuous desire. There are a number of methods that will briefly be presented in the following. Most are focused on the study of linear behavior. But alongside these, improved methods have been developed that can analyze the non-linear behavior [2] of composites. Different mathematical principles are used. For example, a variational principle that also takes into account a time factor allows a relatively simple mathematical description of such cases [3,4,5]. The stress field induced by low temperatures in the body of a polymer composite at low temperatures (a spacecraft application is considered) is presented in [6], where the geometry of the composite microstructure is studied.

One of the methods used, in order to avoid the calculation of stresses and deformations in complex cases, is to consider particular loading cases [7]. Thus, in the mentioned paper, the material constants were calculated for an orthotropic composite and for an isotropic transverse composite. The disadvantage of this approach is that it provides lower and upper bounds for the values of the physical constants. Sometimes these methods can be somewhat or even totally useless because of the large range in which these values can be found [8,9,10,11]. Methods based on the micromechanical analysis of composites were among the first to be used [12,13]. Experimental results show good agreement with the results obtained using this method [4,5]. For the study of the most used composites, those consisting of a matrix and reinforced with fibers, numerous studies have been carried out, and the results have been presented in numerous papers published in recent years [14,15,16,17,18,19,20,21,22,23,24,25]. Other results related to the research of these types of problems are described in [26,27,28,29,30].

To analyze the visco-elastic response of a composite with short fibers, an analytical model is used in [31]. Thus, FEM proves to be a very useful tool in the case of determining the global mechanical properties of composite materials [32].

Special situations that may appear, such as the influence of temperature and humidity on the stress field for a unidirectional graphite/epoxy material, were also studied. These factors were found to lead to a decrease in the tensile strength of a silicon carbide unidirectional fiber composite. Brinson et al. used a similar model in [33,34,35] to determine material constants and how these values depend on the values of the component phase constants of a composite. In the present work, FEM analysis is presented to determine the stress and strain field in a transversely isotropic composite, in order to finally determine the values of the engineering constants for the studied material. Other studies concerning the behavior of these composites are presented in [36,37,38].

Micromechanical models used together with FEM to determine the stress and strain field were successfully used in [39] to determine the Young’s modulus and thermal expansion coefficient of a nanocomposite. For this type of approach, a three-dimensional volume-representative element with randomly distributed SiC reinforcement particles inside was used. Based on the results obtained from the calculation, conclusions that are useful in practical applications were formulated. The performed tests confirm the correctness of the performed calculation. A study of carbon fiber-reinforced polymer composites containing unwanted voids using FEM is presented in [40]. Voids will cause a strong decrease in the mechanical properties of a composite. The results of the study are in agreement with experimental measurements. A micromechanical model of intelligent composites reinforced with piezoelectric constituents is developed in [41]. A hybrid composite reinforced with carbon fibers based on carbon nanotubes where FEM is used for analysis is presented in [42]. The damping characteristics [43] of a unidirectional composite material reinforced with polymer fibers are determined in [44]. Again, micromechanical models show their efficiency in performing such calculations.

In the paper, the representative volume element (RVE) is used for microstructure analysis. This method lends itself well to the study of a large class of materials with very different topologies. Obviously, at the microstructure level, any composite material studied is a heterogeneous material. Its properties are determined by the structure of the composite and the properties of the component materials, matrix and reinforcement material. In the case of a finite element analysis model, these inhomogeneities are neglected. The properties of a single finite element are assumed to be constant.

The homogenized parameters of a composite that define the engineering constants of the material are determined based on a representative volume element using the finite element method. The obtained values agree very well with experimentally determined values. FEM also allows the analysis of some material variants, regarding the geometry and the percentage of the reinforcement material used. Other results concerning the application of FEM in the problems of determining the elastic constant of the composite material multiphasic are presented in [45,46,47,48,49,50,51]. In any of the proposed methods, FEM was always used to determine the field of stresses and deformations, which is used in different ways. Fiber-reinforced composites, used in a wide class of applications, have attracted the attention of researchers, and there are numerous studies on the subject. Some more interesting results obtained in the field can be mentioned. Thus, research has dealt with a wide class of materials, including cement [52]. FEM was used to obtain a field of stresses and strains, which was then used with homogenization theory to determine two material constants for cement reinforced with metal fibers. The obtained values proved to be in accordance with the values obtained using other calculation methods. Polymeric composites reinforced with natural fibers are used in many applications at the moment. A study that aims to determine the visco/elastic constants of these materials is presented in [53]. Natural fibers have viscoelastic behavior and superior specific properties. They also have good ecological characteristics and a low cost. The mentioned work proposes an analytical model of homogenization. Experimental checks demonstrate the accuracy of the model used. A method for determining the mechanical properties of composites reinforced with fibers (carbon and glass) is presented in [54]. This type of material has viscoplastic behavior. The creep behavior of such a material subjected to high temperatures is determined using the presented theory and is verified experimentally. Other research in the field for different situations that can be encountered in practice is presented in [55,56].

It is expected that the calculation accuracy is that normally obtained by applying FEM. In general, FEM has proven to be a reliable method that provides credible and sufficiently accurate results for engineering applications.

In this research, based on micromechanical models, the authors used FEM to calculate the mechanical constants of a composite material reinforced with unidirectional cylindrical fibers. Results are obtained for a composite in which the epoxy matrix is reinforced with carbon fibers. The method can be applied to any type of composite reinforced with fibers, if as a whole it can be considered as a homogeneous transversely isotropic body. In the paper, the authors considered several loading cases for the studied specimen, thus obtaining more information regarding the elastic behavior of the material. In order to obtain the best possible result, the least squares method was applied. Thus, an objective criterion for determining the elastic constants is provided. It can be specified that these elastic constants are obtained and no upper or lower bounds are obtained for them, as happens in a series of energetic methods, previously presented.

The presented method is extremely suitable for the design phase, where it is desirable to know sufficiently precise estimates of the mechanical properties of a material in order to create a mechanical structure. Obviously, after the adoption of a solution and a project and the realization of a prototype, experimental measurements can be made to provide numerical values for the estimated properties that allow possible corrections of the project (if it proves necessary).

The estimation method proposed in the paper is based on models previously established by researchers, which establish formula relationships that can be used to help to determine the mechanical properties of a material if the topology, geometry and properties of the component phases are taken into account. Most of the methods are definitively based on the determination of the field of stress and strain for certain load states of the material. This can be achieved, in most cases, only using numerical methods, which are approximate methods. During the work, the least squares method was proposed to find the best possible approximation, based on the calculation of the stress–strain field in a composite body in several loading situations. Thus, a theoretical basis is provided for this type of estimation.

Models used up to now deal with the simulation of particular loading cases on the basis of which one of the elastic constants can be calculated. The disadvantage of such analysis is that the models generally provide upper and lower limits for elastic constants, which are often very imprecise (the given interval is very large). Other methods, such as the homogenization method, imply as a preliminary stage knowledge of the stress and strain field, a difficult thing to calculate in general, which leads to the use of numerical calculation methods, implicitly FEM. After determining these values, quite difficult calculations follow to obtain these coefficients. The least squares method offers a more direct approach; the method is well known, and the consideration of a larger number of analyzed cases leads to more accurate values of the elastic constants with the help of the least squares criterion.

2. Materials and Methods

Let us consider a composite material made of an epoxy matrix reinforced with aligned cylindrical fibers (Figure 1). The orientation of the fibers defines the axis of symmetry. If a reference frame is considered, as in Figure 1, the properties along the x direction are different if the properties along the y and z directions are considered. In the y and z directions, the properties are the same.

Plane yz is an isotropic plane. The properties of such a material are defined by five parameters (the longitudinal (axial) Young’s modulus E_A, the transverse Young’s modulus E_T, the transverse and longitudinal Poisson’s ratio (ν_A and ν_T) and the shear modulus G_A).

The composite material as a whole is treated as a transversely isotropic material. The properties that are used in strength calculations represent the global properties of the composite, which is composed of two phases (fiber and matrix) that have different properties. Obviously, practical applications require knowledge of the global properties of the material. Starting from the individual properties of the two phases, it is necessary to determine the properties of the material. Various calculation methods have been used to determine them. The first methods applied were methods that used micromechanical models [10,11,12,13]. However, these led to the need to know the field of stress and strain in the resulting material and to obtain the average of the values for stresses and strains. In general, this method, which involves numerical calculations, was used to determine only some engineering constants, simply obtained through special loading states [1]. Variational methods also require the knowledge of functions that define the stress and deformation field. Other calculation methods [1] operate with simplifications of the model and approximations for the performed calculations. In general, in most of the works that deal with these problems, lower and higher bounds for engineering constants (buck modulus, Young’s moduli and Poisson’s ratios) are determined. Although useful, for certain areas of fiber concentration (generally towards the middle of the range), the errors can be huge. This requires the use of methods that bring us closer to the real values.

In this work, the authors proposed to use a micromechanical model for an elementary cell (which takes into account a single fiber of a given length embedded in a matrix). This elementary cell is an RVE, and the material in the assembly can be conceived as an overlap of such RVEs. Using the results from micromechanical analysis to determine the equivalent engineering constants of the RVE, it is necessary to determine the stress and strain field for an arbitrary load and then to determine the average of these quantities [1]. Considering fundamental relationships in the mechanics of transversely isotropic continuous media, the values of the equivalent engineering constants can then be determined.

The study is carried out for an RVE and for a package of elementary cells (Figure 2).

Equations that provide the relationship between stresses and strains, assuming isothermal conditions, considering a transversely isotropic material, are as follows [13]:

$$\left\{\begin{array}{c}{o}_{xx}\\ o_{yy}\\ o_{zz}\\ o_{yz}\\ o_{zx}\\ o_{xy}\end{array}\right\}=\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_{44}& 0& 0\\ 0& 0& 0& 0& C_{66}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{zz}\\ 2{\epsilon }_{yz}\\ 2{\epsilon }_{zx}\\ 2{\epsilon }_{xy}\end{array}\right\}$$

(1)

where:

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

(2)

or:

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

(3)

Considering microstructural models [14,15], it can be considered that any part of the composite behaves as a homogeneous transversely isotropic body [57,58].

As a result, it is considered that if an RVE is taken, its behavior to any demands is, on average, identical to the behavior of the entire isotropic transverse body. In essence, the stress and strain field for an RVE is calculated, averages of the stress and strain values for the entire RVE are calculated, and then it is considered that these values must comply with Hooke’s generalized law for transversely isotropic materials. In this way, a set of equations is obtained that connects the average stress and strain to the material constants. The obtained set of equations is used to determine these constants.

Obviously, FEM, representing an approximation method, introduces errors into calculations made. For this reason, a single calculation of the requests that appear in the RVE could bring errors in the determination of the material constants. For this reason, a large number of calculations are made with different types of load, after which the least square method is used to determine the best approximation for the material constants.

In the following, the average of a size will be denoted with a bar over the symbol of the size.

With FEM, the values of stresses and strains are determined numerically for an RVE consisting of a fiber incorporated in a matrix. A case of loading the RVE is considered. If $n_f$ is the number of finite elements through which the fiber that is included in the RVE was discretized, $n_m$ is the number of finite elements in which the matrix was discretized and $n_t=n_f+n_m$ is the total number of finite elements; after FEM modeling, sets of equations with stress values and strains ${\left\{\sigma \right\}}_k,{\left\{\epsilon \right\}}_k;k=\overline{1,n_f}$, ${\left\{\sigma \right\}}_l,{\left\{\epsilon \right\}}_l;l=\overline{1,n_m}$ are used for each considered finite element. Vectors $\left\{\sigma \right\}$ and $\left\{\epsilon \right\}$ have the following components:

$$\left\{\sigma \right\}=\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{zz}\\ {\sigma }_{yz}\\ {\sigma }_{zx}\\ {\sigma }_{xy}\end{array}\right\};\left\{\epsilon \right\}=\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{zz}\\ {\gamma }_{yz}\\ {\gamma }_{zx}\\ {\gamma }_{xy}\end{array}\right\}=\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{zz}\\ 2{\epsilon }_{yz}\\ 2{\epsilon }_{zx}\\ 2{\epsilon }_{xy}\end{array}\right\}.$$

(4)

Averages of these sizes for the entire RVE are given by:

$$\left\{\overline{\sigma }\right\}={\nu }_f{\left\{\overline{\sigma }\right\}}^{\left(f\right)}+{\nu }_m{\left\{\overline{\sigma }\right\}}^{\left(m\right)};$$

(5)

$$\left\{\overline{\epsilon }\right\}={\nu }_f{\left\{\overline{\epsilon }\right\}}^{\left(f\right)}+{\nu }_m{\left\{\overline{\epsilon }\right\}}^{\left(m\right)},$$

(6)

where index (f) is used to denote a property for a fiber and index (m) is used for a matrix. The average stresses in the fiber and matrix, respectively, are:

$${\left\{\overline{\sigma }\right\}}^{\left(f\right)}=\frac{{\displaystyle \sum _{i=1}^{n_f}{\sigma }_i^{(f)}}}{n_f};{\left\{\overline{\sigma }\right\}}^{\left(m\right)}=\frac{{\displaystyle \sum _{i=1}^{n_m}{\sigma }_i^{(m)}}}{n_m},$$

(7)

and the average strains in the fiber and matrix, respectively, are:

$${\left\{\overline{\epsilon }\right\}}^{\left(f\right)}=\frac{{\displaystyle \sum _{i=1}^{n_f}{\epsilon }_i^{(f)}}}{n_f};{\left\{\overline{\epsilon }\right\}}^{\left(m\right)}=\frac{{\displaystyle \sum _{i=1}^{n_m}{\epsilon }_i^{(m)}}}{n_m}.$$

(8)

Here, ${\sigma }_i^{(f)}$, ${\sigma }_i^{(m)}$, ${\epsilon }_i^{(f)}$, ${\epsilon }_i^{(m)}$ represent the vectors of stress and strain in the RVE, labelled with index i for fiber and matrix, respectively. Compared to the exact theoretical value that exists when considering the RVE as an isotropic transverse body, there is a deviation that can be written as follows:

$$\left\{\delta \right\}=\left\{\overline{\sigma }\right\}-\left[C\right]\left\{\overline{\epsilon }\right\},$$

(9)

with components:

$${\delta }_1={\overline{\sigma }}_{xx}-C_{11}{\overline{\epsilon }}_{xx}-C_{12}\left({\overline{\epsilon }}_{yy}+{\widehat{\epsilon }}_{zz}\right),$$

(10)

$${\delta }_2={\overline{\sigma }}_{yy}-C_{12}{\overline{\epsilon }}_{xx}-C_{22}{\overline{\epsilon }}_{yy}-C_{23}{\overline{\epsilon }}_{zz},$$

(11)

$${\delta }_3={\overline{\sigma }}_{zz}-C_{12}{\overline{\epsilon }}_{xx}-C_{23}{\overline{\epsilon }}_{yy}-C_{22}{\overline{\epsilon }}_{zz},$$

(12)

$${\delta }_4={\overline{\sigma }}_{yz}-\left(C_{22}-C_{23}\right){\overline{\epsilon }}_{yz},$$

(13)

$${\delta }_5={\overline{\sigma }}_{zx}-2C_{66}{\overline{\epsilon }}_{zx},$$

(14)

$${\delta }_6={\dot{\sigma }}_{xy}-2C_{66}{\overline{\epsilon }}_{xy},$$

(15)

The total squared deviation for the considered load case is

$$\left\{{\delta }^2\right\}={\left\{\delta \right\}}^T\left\{\delta \right\}={\left\{\overline{\sigma }\right\}}^T\left\{\overline{\sigma }\right\}-{\left\{\overline{\sigma }\right\}}^T\left[C\right]\left\{\overline{\epsilon }\right\}-{\left\{\overline{\epsilon }\right\}}^T{\left[C\right]}^T\left\{\overline{\sigma }\right\}+{\left\{\overline{\epsilon }\right\}}^T{\left[C\right]}^T\left[C\right]\left\{\overline{\epsilon }\right\},$$

(16)

Let us consider that the FEM calculation is now conducted for a particular case of loading. Let us index with i the considered case. In the case under study, different loading cases will be considered, so $i=\overline{1,n}$. In this case, the components of the deviation vector for the considered state of load will be noted as f:

$${\delta }_1(i)={\overline{\sigma }}_{xx,i}-C_{11}{\overline{\epsilon }}_{xx,i}-C_{12}\left({\overline{\epsilon }}_{yy,i}+{\overline{\epsilon }}_{zz,i}\right),$$

(17)

$${\delta }_2(i)={\overline{\sigma }}_{yy,i}-C_{12}{\overline{\epsilon }}_{xx,i}-C_{22}{\overline{\epsilon }}_{yy,i}-C_{23}{\overline{\epsilon }}_{zz,i},$$

(18)

$${\delta }_3(i)={\overline{\sigma }}_{zz,i}-C_{12}{\overline{\epsilon }}_{xx,i}-C_{23}{\overline{\epsilon }}_{yy,i}-C_{22}{\overline{\epsilon }}_{zz,i},$$

(19)

$${\delta }_4(i)={\overline{\sigma }}_{yz,i}-\left(C_{22}-C_{23}\right){\overline{\epsilon }}_{yz,i},$$

(20)

$${\delta }_5(i)={\overline{\sigma }}_{zx,i}-2C_{66}{\overline{\epsilon }}_{zx,i},$$

(21)

$${\delta }_6(i)={\overline{\sigma }}_{xy,i}-2C_{66}{\overline{\epsilon }}_{xy,i},$$

(22)

The squared deviation for the loading state labelled i is:

$$\begin{array}{c}{\delta }^2(i)={\delta }_1^2(i)+{\delta }_2^2(i)+{\delta }_3^2(i)+{\delta }_4^2(i)+{\delta }_5^2(i)+{\delta }_6^2(i)\\ ={\left[{\overline{\sigma }}_{xx,i}-C_{11}{\overline{\epsilon }}_{xx,i}-C_{12}\left({\overline{\epsilon }}_{yy,i}+{\overline{\epsilon }}_{zz,i}\right)\right]}^2+{\left[{\overline{\sigma }}_{yy,i}-C_{12}{\overline{\epsilon }}_{xx,i}-C_{22}{\overline{\epsilon }}_{yy,i}-C_{23}{\overline{\epsilon }}_{zz,i}\right]}^2\\ +{\left[{\overline{\sigma }}_{zz,i}-C_{12}{\overline{\epsilon }}_{xx,i}-C_{23}{\overline{\epsilon }}_{yy,i}-C_{22}{\overline{\epsilon }}_{zz,i}\right]}^2+{\left[{\overline{\sigma }}_{yz,i}-\left(C_{22}+C_{23}\right){\overline{\epsilon }}_{yz,i}\right]}^2\\ +{\left[{\overline{\sigma }}_{zx,i}-2C_{66}{\overline{\epsilon }}_{zx,i}\right]}^2+{\left[{\overline{\sigma }}_{xy,i}-2C_{66}{\overline{\epsilon }}_{xy,i}\right]}^2\end{array}$$

(23)

The total squared deviation is:

$$\begin{array}{c}{\delta }^2={\displaystyle \sum _{i=1}^n\delta {(i)}^2=}{\displaystyle \sum _{i=1}^n{\left[{\overline{\sigma }}_{xx,i}-C_{11}{\overline{\epsilon }}_{xx,i}-C_{12}\left({\overline{\epsilon }}_{yy,i}+{\overline{\epsilon }}_{zz,i}\right)\right]}^2}\\ +{\displaystyle \sum _{i=1}^n{\left[{\overline{\sigma }}_{yy,i}-C_{12}{\overline{\epsilon }}_{xx,i}-C_{22}{\overline{\epsilon }}_{yy,i}-C_{23}{\overline{\epsilon }}_{zz,i}\right]}^2}\\ +{\displaystyle \sum _1^n{\left[{\overline{\sigma }}_{zz,i}-C_{12}{\overline{\epsilon }}_{xx,i}-C_{23}{\overline{\epsilon }}_{yy,i}-C_{22}{\overline{\epsilon }}_{zz,i}\right]}^2}+{\displaystyle \sum _1^n{\left[{\overline{\sigma }}_{yz,i}-\left(C_{22}+C_{23}\right){\overline{\epsilon }}_{yz,i}\right]}^2}\\ +{\displaystyle \sum _1^n{\left[{\overline{\sigma }}_{zx,i}-2C_{66}{\overline{\epsilon }}_{zx,i}\right]}^2}+{\displaystyle \sum _1^n{\left[{\overline{\sigma }}_{xy,i}-2C_{66}{\overline{\epsilon }}_{xy,i}\right]}^2}.\end{array}$$

(24)

The condition for this squared deviation to be minimal is given by relationships offered by mathematical analysis:

$$\frac{\partial {\delta }^2}{\partial C_{11}}=-2\left[{\displaystyle \sum _{i=1}^n{\overline{\sigma }}_{xx,i}}{\overline{\epsilon }}_{xx,i}-C_{11}{\displaystyle \sum _{i=1}^n{\overline{\epsilon }}_{xx,i}^2}-C_{12}\left({\displaystyle \sum _{i=1}^n{\overline{\epsilon }}_{yy,i}}{\overline{\epsilon }}_{xx,i}+{\displaystyle \sum _{i=1}^n{\overline{\epsilon }}_{zz,i}{\overline{\epsilon }}_{xx,i}}\right)\right]=0,$$

(25)

$$\begin{array}{c}\frac{\partial {\delta }^2}{\partial C_{12}}=-2\left[{\displaystyle \sum _{i=1}^n{\overline{\sigma }}_{xx,i}}\left({\overline{\epsilon }}_{yy,i}+{\overline{\epsilon }}_{zz,i}\right)-C_{11}{\displaystyle \sum _{i=1}^n{\overline{\epsilon }}_{xx,i}}\left({\overline{\epsilon }}_{yy,i}+{\overline{\epsilon }}_{zz,i}\right)-C_{12}{\displaystyle \sum _{i=1}^n{\left({\overline{\epsilon }}_{yy,i}+{\overline{\epsilon }}_{zz,i}\right)}^2}\right]\\ -2\left[{\displaystyle \sum _1^n{\overline{\sigma }}_{yy,i}{\overline{\epsilon }}_{xx,i}}-C_{12}{\displaystyle \sum _1^n{\overline{\epsilon }}_{xx,i}^2}-C_{22}{\displaystyle \sum _1^n{\overline{\epsilon }}_{xx,i}{\overline{\epsilon }}_{yy,i}}-C_{23}{\displaystyle \sum _1^n{\overline{\epsilon }}_{xx,i}{\overline{\epsilon }}_{zz,i}}\right]\\ -2\left[{\displaystyle \sum _1^n{\overline{\sigma }}_{zz,i}{\overline{\epsilon }}_{xx,i}}-C_{12}{\displaystyle \sum _1^n{\overline{\epsilon }}_{xx,i}^2}-C_{23}{\displaystyle \sum _1^n{\overline{\epsilon }}_{xx,i}{\overline{\epsilon }}_{yy,i}}-C_{22}{\displaystyle \sum _1^n{\overline{\epsilon }}_{xx,i}{\overline{\epsilon }}_{zz,i}}\right]=0,\end{array}$$

(26)

$$\begin{array}{c}\frac{\partial {\delta }^2}{\partial C_{22}}=-2{\left[{\displaystyle \sum _i^n{\overline{\sigma }}_{yy,i}{\overline{\epsilon }}_{yy,i}}-C_{12}{\displaystyle \sum _1^2{\overline{\epsilon }}_{xx,i}{\overline{\epsilon }}_{yy,i}}-C_{22}{\displaystyle \sum _1^n{\overline{\epsilon }}_{yy,1}^2}-C_{23}{\displaystyle \sum _1^n{\overline{\epsilon }}_{yy,i}{\overline{\epsilon }}_{zz,i}}\right]}^2\\ -2\left[{\displaystyle \sum _1^n{\overline{\sigma }}_{zz,i}{\overline{\epsilon }}_{zz,i}}-C_{12}{\displaystyle \sum _1^n{\overline{\epsilon }}_{xx,i}{\overline{\epsilon }}_{zz,i}}-C_{23}{\displaystyle \sum _1^n{\overline{\epsilon }}_{yy,i}{\overline{\epsilon }}_{zz,i}}-C_{22}{\displaystyle \sum _1^n{\widehat{\epsilon }}_{zz,i}^2}\right]\\ -2\left[{\displaystyle \sum _1^n{\overline{\sigma }}_{yz,i}{\overline{\epsilon }}_{yz,i}}-\left(C_{22}+C_{23}\right){\displaystyle \sum _1^n{\overline{\epsilon }}_{yz,i}^2}\right]=0,\end{array}$$

(27)

$$\begin{array}{c}\frac{\partial {\delta }^2}{\partial C_{23}}=-2\left[{\displaystyle \sum _1^n{\overline{\sigma }}_{yy,i}{\overline{\epsilon }}_{zz,i}}-C_{12}{\displaystyle \sum _i^n{\overline{\epsilon }}_{xx,i}{\overline{\epsilon }}_{zz,i}}-C_{22}{\displaystyle \sum _1^n{\overline{\epsilon }}_{yy,i}{\overline{\epsilon }}_{zz,i}}-C_{23}{\displaystyle \sum _1^n{\overline{\epsilon }}_{zz,i}^2}\right]\\ -2\left[{\displaystyle \sum _1^n{\overline{\sigma }}_{zz,i}{\overline{\epsilon }}_{yy,i}}-C_{12}{\displaystyle \sum _1^n{\overline{\epsilon }}_{xx,i}{\overline{\epsilon }}_{yy,i}}-C_{23}{\displaystyle \sum _1^n{\overline{\epsilon }}_{yy,i}^2}-C_{22}{\displaystyle \sum _1^n{\overline{\epsilon }}_{zz,i}{\overline{\epsilon }}_{yy,i}}\right]\\ -2\left[{\displaystyle \sum _1^n{\overline{\sigma }}_{yz,i}{\overline{\epsilon }}_{yz,i}}-\left(C_{22}+C_{23}\right){\displaystyle \sum _1^n{\overline{\epsilon }}_{yz,i}^2}\right]=0,\end{array}$$

(28)

$$\frac{\partial {\delta }^2}{\partial C_{66}}=-2\left[{\displaystyle \sum _1^n{\overline{\sigma }}_{zx,i}{\overline{\epsilon }}_{zx,i}}-2C_{66}{\displaystyle \sum _i^n{\overline{\epsilon }}_{zx,i}^2}\right]-2\left[{\displaystyle \sum _1^n{\overline{\sigma }}_{xy,i}{\overline{\epsilon }}_{xy,i}}-2C_{66}{\displaystyle \sum _i^n{\overline{\epsilon }}_{xy,i}^2}\right]=0,$$

(29)

The system obtained by applying these conditions represents a linear system of five equations with five unknowns. By solving it, the five unknowns can be obtained, which represent material constants with the help of which the usual elastic engineering constants can then be determined: Young’s moduli, shear moduli, bulk modulus, Poisson’s ratios, etc.

Equation (29) immediately results in the following:

$$C_{66}=\frac{{\displaystyle \sum _1^n{\overline{\sigma }}_{zx,i}{\overline{\epsilon }}_{zx,i}}+{\displaystyle \sum _1^n{\overline{\sigma }}_{xy,i}{\overline{\epsilon }}_{xy,i}}}{2\left({\displaystyle \sum _i^n{\overline{\epsilon }}_{zx,i}^2}+{\displaystyle \sum _i^n{\overline{\epsilon }}_{xy,i}^2}\right)},$$

(30)

The system remaining to be solved is presented in Appendix A and offers other constants of materials, $C_{11},C_{12},C_{22},C_{23}$. Using these values now, it is possible to obtain the desired engineering elastic coefficients.

Bulk modulus can be computed by the following relationship:

$$K_{23}=C_{22},$$

(31)

The shear modulus is:

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

(32)

Other relationships between the elastic coefficients give us all the values that are necessary for the resistance calculation of a structure made of the proposed material.

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

(33)

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

(34)

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

(35)

$${\nu }_{23}=\frac{K_{23}-\psi G_{23}}{K_{23}+\psi G_{23}};$$

(36)

$$C_{66}=\frac{{\overline{\sigma }}_{12}}{2{\overline{\epsilon }}_{12}}=\frac{{\overline{\tau }}_{12}}{{\overline{\gamma }}_{12}}=G_{12}=G_{13}=G_1.$$

(37)

3. Results

A composite bar made up of an epoxy matrix reinforced with a circular fiber is considered. In this way, a rectangular parallelepiped is obtained, which represents an RVE. The results obtained for this RVE will be used in the calculations. In a perfected version, a rectangular parallelepiped made up of 16 such RVEs is considered. In this case, a calculation can be made more precisely. To see the potential application of the proposed method, a micromechanical FEM model is developed for a composite bar made of a glass/epoxy matrix reinforced with carbon fibers [59].

Values of the main material properties for the two phases of the composite material are:

E_m = 4.14 GPa; ν_m = 0.22; E_f = 86.90 GPa; ν_f = 0.34.

The two materials are considered homogeneous and isotropic. They are defined by two material constants. For the situation studied in the paper, the fiber volume ration of the total volume is $v_f=54.54\%$ and represents an important parameter in all methods of homogenization.

Five distinct load cases of the bar are considered (Figure 3). The load is considered to be a force that acts on one of the ends of the bar, the other end being embedded. In the five cases, the force acts under angles equal to 0, 30°, 45°, 60° and 90° from the x axis. For each of these cases, the stress and strain field is calculated and then the average of these stresses and strain is calculated for the entire RVE. Based on these values, elastic constants are then determined.

In the framework of analysis carried out in the work, the stresses and deformations are determined with the help of Altair software 2020 and the Hypermesh preprocessor and Hyperview postprocessor, respectively. These values are then used to determine the average values of stresses and strains for the fiber, the matrix and then for the RVE.

Table 1. Average strain in fiber for four loading cases.

Loading Case	${\overline{\mathit{\epsilon }}}_{\mathit{x}\mathit{x}}$	${\overline{\mathit{\epsilon }}}_{\mathit{y}\mathit{y}}$	${\overline{\mathit{\epsilon }}}_{\mathit{z}\mathit{z}}$	${\overline{\mathit{\epsilon }}}_{\mathit{y}\mathit{z}}$	${\overline{\mathit{\epsilon }}}_{\mathit{z}\mathit{x}}$	${\overline{\mathit{\epsilon }}}_{\mathit{x}\mathit{y}}$
1	3.60 × 10−4	−1.22 × 10−4	−1.22 × 10−4	−1.60 × 10−13	−1.01 × 10−17	−3.86 × 10−17
2	3.60 × 10−4	−1.22 × 10−4	−1.22 × 10−4	−8.54 × 10−8	−3.27 × 10−9	−5.43 × 10−8
3	3.11 × 10−4	−1.05 × 10−4	−1.05 × 10−4	−2.45 × 10−4	−7.56 × 10−9	2.41 × 10−10
4	1.79 × 10−4	−6.03 × 10−5	−6.02 × 10−5	−4.25 × 10−4	−1.46 × 10−8	2.33 × 10−9
5	−1.46 × 10−6	4.46 × 10−7	4.56 × 10−7	−4.91 × 10−4	−1.67 × 10−8	4.06 × 10−11

,

Table 2. Average strain in matrix for four loading cases.

Loading Case	${\overline{\mathit{\epsilon }}}_{\mathit{x}\mathit{x}}$	${\overline{\mathit{\epsilon }}}_{\mathit{y}\mathit{y}}$	${\overline{\mathit{\epsilon }}}_{\mathit{z}\mathit{z}}$	${\overline{\mathit{\epsilon }}}_{\mathit{y}\mathit{z}}$	${\overline{\mathit{\epsilon }}}_{\mathit{z}\mathit{x}}$	${\overline{\mathit{\epsilon }}}_{\mathit{x}\mathit{y}}$
1	3.60 × 10−4	−8.83 × 10−5	−8.83 × 10−5	−4.25 × 10−17	4.38 × 10−11	−1.11 × 10−17
2	3.60 × 10−4	−8.70 × 10−5	−9.15 × 10−5	−1.08 × 10−7	−2.74 × 10−6	−1.42 × 10−7
3	3.11 × 10−4	−7.95 × 10−5	−7.23 × 10−5	−2.48 × 10−4	4.86 × 10−9	−4.83 × 10−10
4	1.79 × 10−4	−4.56 × 10−5	−4.16 × 10−5	−4.29 × 10−4	8.71 × 10−9	−1.05 × 10−9
5	−1.33 × 10−6	3.69 × 10−7	2.32 × 10−7	−4.95 × 10−4	1.43 × 10−8	−1.53 × 10−9

,

Table 3. Average strain in RVE for four loading cases.

Loading Case	${\overline{\mathit{\epsilon }}}_{\mathit{x}\mathit{x}}$	${\overline{\mathit{\epsilon }}}_{\mathit{y}\mathit{y}}$	${\overline{\mathit{\epsilon }}}_{\mathit{z}\mathit{z}}$	${\overline{\mathit{\epsilon }}}_{\mathit{y}\mathit{z}}$	${\overline{\mathit{\epsilon }}}_{\mathit{z}\mathit{x}}$	${\overline{\mathit{\epsilon }}}_{\mathit{x}\mathit{y}}$
1	3.60 × 10−4	−1.09 × 10−4	−1.09 × 10−4	−2.84 × 10−7	−9.09 × 10−7	−2.84 × 10−7
2	3.60 × 10−4	−1.06 × 10−4	−1.08 × 10−4	−9.58 × 10−8	−1.25 × 10−6	−9.40 × 10−8
3	3.11 × 10−4	−9.33 × 10−5	−9.01 × 10−5	−2.46 × 10−4	−1.92 × 10−9	−8.84 × 10−11
4	1.79 × 10−4	−5.36 × 10−5	−5.18 × 10−5	−4.27 × 10−4	−4.01 × 10−9	7.93 × 10−10
5	−1.40 × 10−6	4.11 × 10−7	3.54 × 10−7	−4.93 × 10−4	−2.58 × 10−9	−6.73 × 10−10

,

Table 4. Average stress in fiber for four loading cases (in GPa).

Loading Case	${\overline{\mathit{\sigma }}}_{\mathit{x}\mathit{x}}$	${\overline{\mathit{\sigma }}}_{\mathit{y}\mathit{y}}$	${\overline{\mathit{\sigma }}}_{\mathit{z}\mathit{z}}$	${\overline{\mathit{\sigma }}}_{\mathit{y}\mathit{z}}$	${\overline{\mathit{\sigma }}}_{\mathit{z}\mathit{x}}$	${\overline{\mathit{\sigma }}}_{\mathit{x}\mathit{y}}$
1	3.14 × 10	4.71 × 10−2	4.71 × 10−2	−4.38 × 10−8	−5.65 × 10−13	−2.24 × 10−12
2	3.60 × 10−4	−1.22 × 10−4	−1.22 × 10−4	−8.54 × 10−8	−3.27 × 10−9	−5.43 × 10−8
3	3.11 × 10−4	−1.05 × 10−4	−1.05 × 10−4	−2.45 × 10−4	−7.56 × 10−9	2.41 × 10−10
4	1.79 × 10−4	−6.03 × 10−5	−6.02 × 10−5	−4.25 × 10−4	−1.46 × 10−8	2.33 × 10−9
5	−1.46 × 10−6	4.46 × 10−7	4.56 × 10−7	−4.91 × 10−4	−1.67 × 10−8	4.06 × 10−11

,

Table 5. Average stress in matrix for four loading cases (in GPa).

Loading Case	${\overline{\mathit{\sigma }}}_{\mathit{x}\mathit{x}}$	${\overline{\mathit{\sigma }}}_{\mathit{y}\mathit{y}}$	${\overline{\mathit{\sigma }}}_{\mathit{z}\mathit{z}}$	${\overline{\mathit{\sigma }}}_{\mathit{y}\mathit{z}}$	${\overline{\mathit{\sigma }}}_{\mathit{z}\mathit{x}}$	${\overline{\mathit{\sigma }}}_{\mathit{x}\mathit{y}}$
1	1.46 × 100	−5.72 × 10−2	−5.72 × 10−2	1.08 × 10−18	−3.93 × 10−14	−1.60 × 10−13
2	3.60 × 10−4	−8.70 × 10−5	−9.15 × 10−5	−1.08 × 10−7	−2.74 × 10−6	−1.42 × 10−7
3	3.11 × 10−4	−7.95 × 10−5	−7.23 × 10−5	−2.48 × 10−4	4.86 × 10−9	−4.83 × 10−10
4	1.79 × 10−4	−4.56 × 10−5	−4.16 × 10−5	−4.29 × 10−4	8.71 × 10−9	−1.05 × 10−9
5	−1.33 × 10−6	3.69 × 10−7	2.32 × 10−7	−4.95 × 10−4	1.43 × 10−8	−1.53 × 10−9

,

Table 6. Average stress in RVE for four loading cases (in GPa).

Loading Case	${\overline{\mathit{\sigma }}}_{\mathit{x}\mathit{x}}$	${\overline{\mathit{\sigma }}}_{\mathit{y}\mathit{y}}$	${\overline{\mathit{\sigma }}}_{\mathit{z}\mathit{z}}$	${\overline{\mathit{\sigma }}}_{\mathit{y}\mathit{z}}$	${\overline{\mathit{\sigma }}}_{\mathit{z}\mathit{x}}$	${\overline{\mathit{\sigma }}}_{\mathit{x}\mathit{y}}$
1	1.78 × 10	−7.15 × 10−3	−7.14 × 10−3	−9.64 × 10−4	−3.08 × 10−3	−9.64 × 10−4
2	3.60 × 10−4	−1.06 × 10−4	−1.08 × 10−4	−9.58 × 10−8	−1.25 × 10−6	−9.40 × 10−8
3	3.11 × 10−4	−9.33 × 10−5	−9.01 × 10−5	−2.46 × 10−4	−1.92 × 10−9	−8.84 × 10−11
4	1.79 × 10−4	−5.36 × 10−5	−5.18 × 10−5	−4.27 × 10−4	−4.01 × 10−9	7.93 × 10−10
5	−1.40 × 10−6	4.11 × 10−7	3.54 × 10−7	−4.93 × 10−4	−2.58 × 10−9	−6.73 × 10−10

and

Table 7. Computed values of elastic moduli.

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.28
${\nu }_{23}$	0.34	0.22	0.31
$G_{23}$	1544.0	35,614.7	12,318.2
$K_{23}$	4827.4	63,597.7	27,953.2

present the results of the analysis.

The same calculations were also carried out considering a material cell made up of 16 fibers (Figure 4).

For a load case, the stresses in the sections made in the first, second, third and fourth rows are represented (Figure 5, Figure 6, Figure 7 and Figure 8). It is observed that these tensions vary quite strongly at the ends, unrelated to the overall stress of the bar. This occurs at the ends of the bar: both the loaded end of the bar and the recessed end (Figure 9 and Figure 10). Due to these variations, when the averages were calculated, the values obtained in this part of the bar were removed. In the calculations made in the study, 200 finite elements were removed from one end and 200 finite elements were removed from the other end.

4. Discussion

The theoretical models developed by the researchers presented in the introduction section allow the determination of the elastic constants of a material considering a single RVE. A linear system of five equations with five unknowns is obtained that provides these values. However, these values are obtained for a single loading case of the material, based on the stress and deformation field obtained for that loading, with errors arising from the determination of this field. The method used is FEM, which, being an approximation method, inherently introduces some errors. The obtained constants will be affected by these errors. If different loading situations are considered, the errors will be different and, obviously, will depend on these loading cases. To minimize the importance of these errors in the calculation of the elastic constants, several load cases will be considered and the condition is set to look for those values that minimize the quadratic deviation. It is obvious that a large number of considered load cases will minimize the errors that may occur in some situations where the errors may be too large. It is therefore recommended to consider very different loading situations, so as to cover as wide a range of situations as possible. In the case of the present study, these were considered as five distinct load cases.

Most of the methods used to determine the properties of a composite material obtain results in the form of inequalities, which provide the sought constants in a certain interval. For some concentrations of the phases, these ranges may be wide enough so that the information provided is ineffective. The presented method provides values that are affected only by the errors involved in the FEM. This leads to much more useful results for the designer. Experimental checks for simpler situations, where only a single loading state is considered, are presented in [40] and show good agreement between the computational and the experimental results.

The studies that dealt with the determination of the elastic constants of a material generally considered particular cases of RVE loading and, based on the particular strain and stress fields, determined some of these constants based on direct calculation [7,8,9,10,11,12,13]. Some of the methods obtained only lower and upper bounds of these constants, providing a relatively large range for these constants. In the case of this paper, several loading cases were considered and, based on them, elastic constants were considered using the least squares method. In this way, considering more loading cases (as many as possible) can minimize the error with which the elastic constants are determined. The accuracy with which these values are obtained should be higher than the calculation accuracy when considering a single loading state [60,61,62,63,64].

Other effects, such as thermal, humidity, etc., can be taken into account. But in this case, the variant corresponding to the considered case must be used in modeling. In this case, too, homogenization of the elastic coefficients can be carried out by averaging, but other physical parameters necessary for the studied model can also be obtained. So the method can be applied in any complex situation, but each time a corresponding model must be considered, which usually involves additional parameters.

5. Conclusions

The presented method gives us values of the engineering constants of a composite if these values are known for the components of the composite and if its topology and geometry are known. It is the case often encountered in practice that it is necessary to estimate the properties of a newly composed material. These values obtained are approximate and are affected by errors occurring in the use of FEM. In order to minimize these errors, the paper proposes the method of least squares, a mathematical method that provides the best estimate for a set of values obtained by calculating FEM. It is useful to consider as many load cases as possible to obtain the best estimate. The method has the advantage of using FEM, which has become a common method, verified by countless applications and frequently used in engineering. In this way, the sought values can be obtained with relative ease.

The main novelty is the application of the least squares method, which is a well-known method. With its help, several elastic constants of a material can be determined (in the case of isotropic materials, two; in the case of transversely isotropic materials, five; etc.). The precision in the case of this approximation method increases with the number of cases considered. For a model of a representative structure of the material, that model is made with finite elements for which the field of stresses and deformations is determined for several different loading cases. Based on these values, it is now possible to determine the engineering constants of the material using the least squares method, a simple and easy-to-apply method.

It is mentioned that the study of as many loading situations as possible should give us, according to the least squares theory, the best approximation of the elastic constants respecting the conditions presented in the paper. So the goal when applying the method is to perform as many calculations of the stress and strain field as possible, so as to obtain the most accurate results of the engineering constants sought.