Semi-Analytical Method for Computing Effective Thermoelastic Properties in Fiber-Reinforced Composite Materials

: Effective elastic and thermal properties for isotropic or transversely isotropic thermoelastic ﬁbrous composite materials are obtained. Fibers are distributed with the same periodicity along the two perpendicular directions to the ﬁber orientation. The periodic cell of the composite has a square or hexagonal distribution. Perfect contact between the ﬁber and the matrix is presented. The effective properties are calculated using a semi-analytical method. The semi-analytical method consists of obtaining the differential equations that describe the local problems using the Asymptotic Homogenization Method. Then, these equations are solved using the Finite Element Method. Effective elastic coefﬁcient ( C ), effective thermal expansion coefﬁcient ( α ) and the effective thermal conductivity ( κ ) are obtained. The numerical results are compared with the semi-analytical solution and with results reported by other authors. Additionally, the effective properties for a ﬁber with an elliptical cross section are calculated. Distributions of the ﬁber’s cross section with different orientations are also studied. A MATLAB program for computing the effective coefﬁcients is presented.


Introduction
A composite material is a combination of two or more materials with different physical and (or) chemical properties. A composite material usually has characteristics that do not correspond to those of its isolated components [1]. Thermoelastic composites constitute an important class of materials with a wide variety of applications ranging from aerospace structures and electronic printed circuit boards to recreational and commercial equipment. Some of the most important and useful properties of these composites are lightweight, high strength and stiffness, excellent frictional properties, good resistance to fatigue and retention of these properties at high temperatures [2]. The effective thermomechanical properties of the composite depend on properties of the constituents and the fiber volume fraction.
Previous knowledge of physical properties of composite materials is one of the problems that science faces. Success in practical applications of composites largely depends on the possibility to predict their properties and behavior through the development of appropriate models. Numerous studies and investigations have been carried out to obtain the properties of composite materials [3][4][5][6]. The prediction of global homogenized properties based on knowledge of the composite constituents provides a useful tool for the design of composite materials. Homogenization methods have proven to be powerful techniques for the study of heterogeneous media. Among the homogenization methods, one of the most used methods is the Asymptotic Homogenization Method (AHM). The AHM can provide overall material constants which describe the global behavior of the composite structure. First, the periodic cell is identified; then, by means of the AHM [7][8][9][10][11][12], the original constitutive relations with rapidly oscillating material coefficients are transformed in two sets of mathematical problems. The first set has new physical relations over the periodic cell with constant coefficients which represent the properties of an equivalent homogeneous medium and they are called the effective coefficients of the composite under study. The second set allows the calculation of some functions which are solutions of the so-called local problems. The local problems are solved over the periodic cell. A state-of-the-art review where the AHM is applied to composites is provided in ref. [13]. A variety of existing homogenization methods are presented with their advantages and disadvantages.
Previous research has been performed for the estimation of the effective properties in composite materials. The elastic properties of carbon nanotube/polymer composites are computed in ref. [14]. A combination of the Finite Element Method (FEM) and the Strain Energy Method was used to obtain the effective elastic constants. The nanotubes are assumed to follow a helical path along the length of the array. The finite element software ANSYS 7.1 was used for the analysis. A review on the application of the FEM in the analysis of composite materials is performed in ref. [15]. A numerical homogenization approach to calculate the effective properties of fiber and particle reinforced materials including smart and multifunctional materials with a focus on piezoelectric fiber composites applied to control vibration and noise radiation of structures was presented in ref. [16]. The authors employ the FEM based on homogenization to optimize the material distribution at the micro-scale by applying an evolutionary approach to receive a desired global behavior of a structure at the macro-scale level. The FEM based on the homogenization of a fiber-reinforced composite has also been applied in ref. [17] by using MAPLE symbolic calculations. Improved versions of the Stochastic perturbation-based Boundary Element Method (SBEM) have been formulated and implemented in ref. [18]. The new approach is compared with a previous version of the SBEM in both isotropic and composite cantilever panels.
Additionally, analytical expressions for the effective properties have been found. Apart from their theoretical importance, exact analytical expressions prove to be very useful for the validation of numerical results. Simple closed analytical expressions for the effective elastic coefficients were found in refs. [10,19] for square and hexagonal cell distributions, respectively. The AHM was applied to a fiber-reinforced composite with a cylindrical shape periodically distributed in the matrix. Both, matrix and fiber, had transversely isotropic properties and they were in perfect contact. The formulas obtained are also valid for isotropic constituents.
The effective moduli of a periodic elastic composite material reinforced by straight parallel circular fibers in a square array were determined in ref. [3]. The composite is made of a transversely isotropic material with axial, tangential and normal imperfect contacts at the interface. The investigation was carried out by adopting the Semi-Analytical Method (SAFEM) which is based on the FEM and the AHM. Quadrilateral quadratic Lagrange functions were used as shape functions in the FEM formulation. In ref. [6], the effective properties of composites materials with fibers of different geometries (circular, square and rhombic) distributed in a square array were studied and three different FEM approaches were implemented: 4-node quadrilaterals, 8-node quadrilaterals and 12-node quadrilaterals. An extension of this SAFEM for coupled field problems was reported in ref. [20]. In addition, a three-dimensional SAFEM was implemented in ref. [21] to calculate the effective properties of periodic elastic-reinforced nanocomposites. Different inclusions were considered, such as discs, ellipsoids, spheres, carbon nanotubes and carbon nanowires.
The effective thermoelastic properties in composite materials have also been extensively studied. A simple closed form of the effective thermoelastic properties was obtained in ref. [8]. Three-phase unidirectional transversely isotropic composites were studied. Composites with hexagonal and square arrangements of cylindrical fibers were considered.
The effective elastic properties and thermal expansion coefficients of unidirectional cylindrical fiber composites with imperfect interface conditions on the basis of the generalized Self Consistent Scheme model were evaluated in ref. [5]. Imperfect interface conditions were defined in terms of spring type coupling constants relating interface displacement discontinuities with tractions. In ref. [22], the thermal analysis of ceramic matrix composites (CMC) were performed. A Representative Volume Element model with periodic boundary conditions and a full-size model with the actual thickness were developed to investigate the temperature field, the heat flux field, and the effective thermal conductivity of the CMC. Ref. [23] considers the thermoelastic response of a composite body reinforced by coated cylindrical fibers and oriented in different directions. Lower bounds of the effective moduli were obtained, which provide a more accurate estimation of the overall composite properties. In ref. [2], a unit cell model was used to predict the effective thermo-mechanical properties of three-phase coated unidirectional cylindrical fibers through homogenization techniques for different fiber volume fractions. The numerical approach was based on the FEM. Longitudinal and transversal effective thermo-mechanical coefficients were calculated with the FEM and compared with analytical solutions based on the AHM.
In the present study, thermoelastic effective properties are obtained for a thermoelastic isotropic or transversely isotropic fiber-reinforced composite. The composite is formed by fibers periodically distributed within a matrix and exhibits perfect contact at the interface between the fiber and matrix, which means that tractions and displacements are continuous across the interface [24]. The results are obtained by means of the SAFEM. The local problems obtained by the AHM are solved through the FEM. The effective coefficients are computed for two different geometrical configurations of cells: square and hexagonal. The FEM implementation is proposed by using 4-node and 8-node quadrilaterals.
The main contribution of this research consists of the generalization of the SAFEM proposed in refs. [3,6] to the field of thermoelastic composites with perfect contact. Additionally, a computational tool for calculating the effective thermoelastic properties of composite materials from the properties of their constituents was created. Finally, the effective properties of composite materials based on elliptical cross-section fibers with different orientations, were estimated. Following this introduction, the reminder of the paper is organized as follows: the statement of the basic equations of thermoelasticity and the fundamental steps of the AHM are presented in Section 2. In Section 3, linear and quadratic finite element formulations are developed. Meshes are constructed to discretize the regions corresponding to the quarter of the periodic cell, and FEM formulation is applied to two reference local problems. Numerical results are shown in Section 4. Effective thermoelastic properties are obtained for different fiber volume fractions. In order to validate the numerical results, comparisons with the semi-analytical expressions reported in ref. [8] are performed. The results are compared with those reported in ref. [23] and they are also shown in the form of figures. Section 5 shows the effective properties for an elliptical cross-section fiber embedded into a square cell distribution. Distributions of the fiber's cross section with different orientations are also studied. In Section 6, some details considered during the implementation of the program are explained. Finally, the conclusions are presented in Section 7.

Problem Statement: Local Problems
We consider a two-phase uniaxial reinforced material where the fibers and matrix have isotropic or transversely isotropic thermoelastic properties. Let the position of a typical point of the body be denoted by three coordinates (x 1 , x 2 , x 3 ) in a Cartesian system of coordinates. Axis x 3 corresponds to the axis of transverse symmetry which points along the direction of the fiber. Additionally, the fibers are periodically distributed along the directions of the axes x 1 and x 2 . Figure 1a illustrates the fiber distribution for square (up) and hexagonal (down) cells.
The characteristic size of the inhomogeneities is assumed to be much smaller than the global size L of the whole structure ( L). The AHM is used to solve the problem. The AHM allows to determine the effective properties of composite materials from a doublescale asymptotic development. The study of the initial problem in a heterogeneous medium ( Figure 1a) becomes an equivalent problem in a homogeneous medium. The problem in the homogeneous medium requires the solution of local problems on a representative cell of the material: the periodic cell Y. Figure 1b represents the periodic cell for the square (up) and the hexagonal (down) distributions. A dimensionless small parameter = /L, 1 characterizing the rate of heterogeneity of the composite structure is introduced. In order to separate the micro and macroscale components of the solution, the so-called slow (x i ) and fast (y i ) coordinates are also introduced The fast variable y i represents the rapid change of properties of the composite or the degree of the heterogeneity. The linear equation of equilibrium in the quasi-static and stationary problem is given by where the subscripts assume the values 1, 2 and 3; the comma denotes partial differentiation and the summation convention is applied; the body forces are represented by X i ; and σ is the stress tensor. The stress (σ) and strain (ε) tensors are related to the temperature change (θ) by the constitutive relations also known as the Duhamel-Neumann law [8] where x = (x 1 , x 2 , x 3 ) and y = (y 1 , y 2 , y 3 ). The constitutive laws contain the components of the elasticity fourth-order tensor C ijkl (y), the components of the thermal stress second-order tensor β ij (y) and α kl (y), which represent the thermal expansion coefficients. The strain tensor is defined by the Cauchy relations where u k are the components of the displacement. The heat flow vector is expressed in terms of the temperature as stated by Fourier's law where κ ij (y) denotes the thermal conductivity components. Using Equation (6), the heat balance equation is written as follows [25] where f is the density of the internal heat sources. Equations (2)-(7) form a coupled linear system of thermoelasticity equations with rapidly oscillating coefficients for an anisotropic inhomogeneous solid. Appropriate boundary and initial conditions should be specified for this system. All coefficients C ijkl , β ij , α ij and κ ij are considered to be picewise-smooth Y-periodic functions of the fast coordinates y i . The idea of the homogenization approach is to obtain a closed system of equations with constant coefficients, equivalent to the system given by Equation (2)-(7). These new coefficients are considered to represent the global properties of the composite. Through the AHM [8], the original constitutive relations with rapidly oscillating material coefficients are transformed in two sets of mathematical problems. The first set has new physical relations over Y with constant coefficients C, β and κ which represent the properties of an equivalent homogeneous medium, and they are called the effective coefficients of the composite. They are calculated by using the following relations [8] where * represents the average operator such that f = Y f dY |Y| , and the subscripts i, j, k, l, p and q assume the values 1, 2 and 3. Equation (8) shows that the effective elastic coefficients can be obtained in the same manner for elastic [10] and thermoelastic materials.
The second set of resulting problems from the AHM allows the calculation of the functions, pq U k (y), pp R k (y) and q M(y) which are solutions of the so-called local problems. The local problems for the elastic ( pq L) and thermal ( pp T, q I) problems are solved over the periodic cell Y and they were derived in ref. [8] as with perfect contact conditions at the interface Γ shown as The vector n represents the outward unit normal vector to the interface Γ and n j are the vector's components. The discontinuity of the function g through the interface is denoted by ||g||, and pq U k (y), pp R k (y) and q M(y) constitute the Y−periodic functions related to the mechanical displacements, the thermal stress tensor and thermal conductivity, respectively. The material coefficients C ijpq (y), β ij (y) and κ ij (y) are known functions of the constituent phases in the local variable y on the periodic cell Y. More details explaining the AHM procedure and the mathematical background can be found in the classical works [13,[26][27][28]. Considering only 1 /4 of the periodic cell's surface to perform the calculations, the periodic boundary conditions must be transformed into mixed boundary conditions as illustrated in Appendix A. Table 1 shows the relationship between the local problems and the effective properties obtained through their solution. Table 1. Relationship between the local problems and the effective properties.
Given the existing symmetry between p and q, there are just six non-homogeneous local problems pq L that contribute to the estimation of the effective elastic properties: the plane problems 11 L, 12 L, 22 L, 33 L and the anti-plane problems 13 L, 23 L. The expressions to compute the elastic coefficients C ijkl are obtained from Equation (8) where Ω represents the periodic cell and A corresponds to the cell's area. Thermal problems pp T are mathematically analogous to purely elastic plane problems. The solutions obtained from problem 11 T are sufficient to determine the effective coefficients β 11 , β 22 and β 33 . Therefore, it is only necessary to solve the differential equations that result from the local problem 11 T. The effective values of the thermal stress tensor β ij are derived from Equation (9) and they can be obtained as Once the C and the β coefficients have been calculated, it is possible to obtain the effective thermal expansion coefficients α through Equation (4).
For thermal capacity problems q I, it can be deduced that problems 1 I and 2 I are analogous to purely elastic anti-plane problems. From Equation (10) it follows that κ 33 is consistent with the Voigt average [8]. The effective thermal conductivity is derived from Equation (10) and can be estimated as follows

Finite Element Method
The thermoelastic effective properties are given by Equation (8)- (10). Additionally, to calculate the effective coefficients it is necessary to solve the local problems given by Equation (11)- (13) to obtain the values: pq U k,l (y), pp R k,l (y) and q M ,k (y). The FEM is one of the most popular numerical methods used to solve differential equations [29]. The region of interest must be discretized into elements in order to apply the FEM. The strong formulation of the differential problem is transformed into a weak formulation or integral problem. The matrix formulation is obtained on each element and then it is assembled into a global system of equations. The shape functions are chosen and incorporated into the formulation. Finally, the resulting system of equations is solved.
The composite material is formed by a matrix with cylindrical fibers periodically distributed in the x 3 direction as mentioned previously. This makes it possible to reduce the problem into a two-dimensional problem. The region of interest will be the periodic cell, since the fibers are periodically distributed. In this work, square and hexagonal cells are studied, as shown in Figure 1a.
Since the representative cell has a symmetric geometry, calculations are performed using only 1 /4 of the cell, which reduces significantly the number of operations and CPU time. In this section, a finite element formulation is performed to solve differential Equations (12) and (13) for the 1 I and 11 T local problems. The remaining local problems are implemented in a similar fashion.

Discretization of the Region
The first step in the FEM implementation is to discretize the region where the differential equation is to be solved, which corresponds to the quarter of the cell. For the square distribution, let us assume that the periodic cell is 1 unit long; therefore, the width and height of the quarter square cell are 0.5 units long. The symmetric quarter of the cell corresponding to a hexagonal distribution is made of one rectangle and two quarters of fibers. The hexagon's side length is assumed to be 1 unit, so the width of the periodical hexagonal cell is 3 /4 units and its height is equal to √ 3/4 units. The (0, 0) coordinate is located in the lower left corner of the cell in both cases.
In this work, the meshes that represent the quarter of the cell are automatically generated. The method described in ref. [30] was used to construct these meshes. First, a description of the geometry to be discretized is provided. In order to generate the initial geometry of the cell, the user must specify the radius of the fiber and the number of elements that the final mesh will have in the y 2 direction. Once these two parameters are defined, the initial blocks are created. Blocks are 8-node quadrilaterals that discretize the region. Each block is made of a single material. The number of subdivisions to be performed in the y 1 and y 2 directions is defined for each block. Additionally, if there are elements with curved sides, then for each of those elements, the coordinates of a point located close to the middle section and over the curved side must be specified. The connectivity of the blocks must also be defined, which is the way in which the nodes that constitute a specific block are related.
Once the initial geometry has been determined, the corresponding mesh is generated. The correspondence between a block and the parent element must be established in order to generate the corresponding mesh. This correspondence can be achieved through the same isoparametric representation that will be used later during the FEM formulation. Figure 2a,b show the meshes obtained using the method described in ref. [30] for square and hexagonal cells, respectively. Linear and quadratic FEM formulations are carried out throughout this investigation. Thus, one more step must be performed to obtain the 4-node quadrilateral mesh when solving the problem by using the linear FEM formulation. The 8-node mesh should be transformed into a 4-node mesh, through division of each 8-node element into four 4-node elements. Figure 3a,b show the meshes obtained by transforming the 8-node quadrilateral meshes presented in Figure 2a,b, respectively. Let Ω be the region where the differential problem must be solved. The region Ω is represented by the mesh previously determined.

Shape Functions
As mentioned in Section 3.1, the region Ω is divided into a finite number of quadrilateral elements. Two types of quadrilateral elements are used: 4-node quadrilateral elements (linear case) and 8-node quadrilateral elements (quadratic case).
The main idea of the FEM is to approximate the unknown function at each discretized element using some approximation functions (also known as shape functions). Once the shape functions have been chosen, they are incorporated into the problem's formulation. There are many ways to choose a set of functions that form a vector basis on which to approximate the exact solution of the problem. Lagrange functions are the most common type of shape functions in the FEM. These functions are piecewise polynomials defined at each element e. In this work, linear and quadratic Lagrange shape functions are used. Therefore, in the linear case, each element will have 4 nodes, while in the quadratic case each element will have 8 nodes. Expressions for the linear and quadratic Lagrange shape functions can be found in ref. [29]. Let N = {N 1 N 2 · · · N n } with n = 4, 8 be the basis of the Lagrange functions in the space where we seek to approximate the solution.
The elements of the mesh that are used to discretize the region Ω are defined by the local coordinates (y 1 , y 2 ). However, the shape functions are generally defined in parent coordinates (ξ, η). Therefore, it is necessary to perform a transformation from local coordinates to parent coordinates and vice versa. The transformation from local coordinates (y 1 , y 2 ) to parent coordinates (ξ, η) obeys the following mapping where vectors y e 1 = [y 11 , y 12 , ... , y 1n ] t and y e 2 = [y 21 , y 22 , ... , y 2n ] t are the known nodal coordinates of the element, and n = 4 (n = 8) has been used for the linear (quadratic) case. Using expressions (18) and the chain rule, the ξ and η derivatives of the ith shape function could be expressed in terms of the y 1 and y 2 derivatives as follows Equation (19) can be written in a more compact form as where is the Jacobian matrix. The Jacobian matrix is shown below for an n-node element (n = 4 or n = 8) In general, each element of the discretization has a different Jacobian matrix, and entries of the Jacobian matrices are functions of (ξ, η).
When the Jacobian determinant |J| is positive at any point in the region, it is possible to invert the transformation in order to determine ξ(y 1 , y 2 ) and η(y 1 , y 2 ). This possibility means that for a given point (y 1 , y 2 ) within the quadrilateral element, there is a single point (ξ, η) in the parent element. Provided that the transformation is invertible, the following can be written where J −1 is the inverse of the Jacobian matrix. Additionally, the surface differentials at each coordinate system are related as follows 3.3. The 1 I Local Problem: Weak Formulation The differential equation used to calculate the thermal conductivity for the problem 1 I can be obtained from Equation (13) and it is given by Perfect contact conditions (14) for problem qI are given by ||(κ 11 + κ 11 1 M ,1 )n 1 + κ 22 1 M ,2 n 2 || Γ = 0, ||κ 11 n 1 + τ 11 n 1 + τ 22 n 2 || Γ = 0, ||τ 11 n 1 + τ 22 n 2 || Γ = −||κ 11 n 1 || Γ , where τ 11 = κ 11 1 M ,1 , τ 22 = κ 22 1 M ,2 and ||g|| denotes the discontinuity or jump of the function g across the interface. The jump at the interface Γ for the conductivity problem is defined by where κ m 11 and κ f 11 are the matrix and fiber thermal conductivity κ 11 coefficient, respectively. By using symmetry and periodicity conditions, the boundary conditions for problem 1 I are obtained where ∂Ω L , ∂Ω R , ∂Ω U and ∂Ω D represent the left, right, upper and lower boundaries of the cell quarter.
Assuming that we look for the approximate solutionM ≈ M of the problem (25) within a certain space of functions H with dimension n, where N = [N 1 N 2 · · · N n ] is a basis of H and q = [q 1 q 2 · · · q n ] t are the nodal values of the unknown field. The weak formulation of problem 1 I is obtained by applying the weighted residual method and Galerkin's method [30,31]. The differential problem (25) Expressions (29)-(31) are the element stiffness matrix, the element force vector and the integral vector that corresponds to the boundary of each element. The vector B e must be evaluated only at the edges that belong to the interface Γ. Equation (26) is used to determine the vector B e . Once the matrices and vectors of each element have been calculated, they are assembled to form the global system of equations which results in the following expression
The discontinuity at the interface Γ is defined by where the superscripts m and f are used for the matrix and fiber coefficients, respectively. Additionally, using symmetry and periodicity conditions, the obtained boundary conditions for problem 11 T are given by where ∂Ω L , ∂Ω R , ∂Ω U and ∂Ω D represent the left, right, upper and lower boundaries of the cell quarter. Suppose that we look for the approximate solutionR = R 1R2 t ≈ R of the problem (33) within certain space of functions H with dimension n such that where {N 1 N 2 · · · N n } is a basis of the space H and N = N 1 0 N 2 0 · · · N n 0 0 N 1 0 N 2 · · · 0 N n .
Vector q = [q 11 q 21 q 12 q 22 · · · q 1n q 2n ] t represents the nodal values of the unknown field. The positions 2i − 1 and 2i (i = 1, ..., n) in vector q are in agreement with the y 1 and y 2 directions of the ith node in the discretization. In order to develop the finite element formulation, the weak formulation of the problem 11 T is obtained by applying the weighted residual method and Galerkin's method [30,31] to Equation (33), resulting in the following integral equation where The vector B e must be evaluated only at the edges that belong to the interface Γ. Equation (36) will be used to calculate the vector B e . Once the element matrices and vectors have been calculated, they are assembled into the following global system of equations that must be solved

Numerical Results
Once the local problems have been solved, it is possible to compute the effective thermoelastic properties of the composite material. The effective elastic coefficients C ijpq , the effective thermal stress β ij and the effective thermal conductivity κ iq are obtained by substituting the calculated values of pq U k , pp R k and q M into Equations (15), (16) and (17). The effective thermal expansion coefficient α pq can be calculated by substituting β ij and C ijpq into Equation (4).
In this section, C ijpq , α pq and κ iq are computed for two isotropic and one transversely isotropic composite materials. The first isotropic material is formed by Nicalon fibers embedded into a matrix of Barium Magnesium Aluminosilicate (Nicalon/BMAS). The second isotropic composite material is formed by Aluminum, Boron and Silicon fibers embedded within an Epoxy matrix (Al-B-Si/Epoxy). Table 2 shows the values of the elastic (Young's modulus E, Poisson's ratio ν, shear modulus G) and thermal (thermal expansion α, thermal conductivity κ) constants for the fiber and matrix of the two isotropic composite materials.
The transversely isotropic material is constituted by circular graphite fibers embedded into an epoxy matrix (Graphite/Epoxy). Table 3 shows the elastic and thermal constants for the fiber and the matrix of the Graphite/Epoxy composite material. The three experiments were solved for square and hexagonal cell distributions of circular fibers composite materials. Different fiber volume fraction values γ were used to perform the calculations. Linear (4-nodes) and quadratic (8-nodes) quadrilateral approximations were used for the FEM implementation. The results obtained were compared with the values provided by the semi-analytical expressions proposed in ref. [8].
The calculated thermoelastic effective properties are illustrated in this section. In each figure, the effective coefficients obtained by using the linear and quadratic implementations of the FEM for the square and hexagonal fiber distributions are compared. Results calculated for the square and hexagonal cells are shown in blue and red with different intensities, respectively. The values of the effective coefficients determined through the semi-analytical formulas, those calculated by using the linear FEM implementation and those obtained through the quadratic FEM formulation are shown with a square, a triangle and a point, respectively. The superscripts A, 4 and 8 in the legend correspond to the values of the semi-analytical expressions, and those obtained through the linear and quadratic FEM formulations, respectively.
In general, Figures 4-6 show that the effective coefficients obtained using both FEM implementations are close to the semi-analytical values. Furthermore, it can be seen that the values calculated by means of the quadratic FEM formulation are closer to the semianalytical values than those obtained by using the linear FEM formulation. Additionally, some figures show an inset with a smaller scale to show relative differences between the semi-analytic and numerical solutions. Finally, the results obtained when performing the calculations using the square and hexagonal distributions are very similar. Due to the symmetry of the tensors and the periodic cell geometry, there are effective coefficients that yield the same results. Only effective coefficients that have different values are illustrated. Additionally, supplementary data have been included in ref. [32].
The first numerical experiment corresponds to a Nicalon/BMAS composite material. Calculated effective properties for different fiber volume fraction values γ are shown. As illustrated in Figure 4, semi-analytical expressions and both FEM implementations provide similar results. The results obtained for the square and hexagonal cell distributions are also similar, differing slightly for large volume fraction values.   From the computed elastic coefficients, it is possible to calculate the effective Young's modulus and the effective Poisson's ratio. Table 4 shows the effective Young's modulus (E 4 ) and the effective Poisson's ratio (ν 4 ) values obtained by using the linear FEM approximation and the results reported in ref. [23] (E and ν). As illustrated in Table 4, the estimated values of E 4 and ν 4 are very similar to those reported in ref. [23].
The effective thermal expansion for the Nicalon/BMAS composite material was also calculated in ref. [23] for a square cell distribution. Dots were used in Figure 7 to illustrate the effective thermal expansion coefficients reported in ref. [23]. The results calculated through the linear SAFEM model presented in this work are illustrated with square symbols. The relative difference between the value of α 22 and α 11 can not be appreciated within the scale shown in Figure 7. Note that for a volume fraction of γ = 0, the effective thermal expansion values correspond to the thermal expansion coefficient of the matrix. Table 4. Young's modulus (E 4 ) and Poisson's ratio (ν 4 ) for the Nicalon/BMAS composite material calculated by using the linear FEM approximation and those reported in [23].  The second example corresponds to the Al-B-Si/Epoxy composite material. This is also an isotropic composite material. Effective elastic properties C, effective thermal expansion α and effective thermal conductivity κ are illustrated in Figure 5 for square and hexagonal distributions. Calculations for the linear and quadratic FEM approximations, and the semi-analytical expressions are preformed for each coefficient. Results obtained from both FEM implementations are very similar to those estimated when applying the expressions proposed in [8].
The last example belongs to a transversely isotropic composite material formed by Graphite fibers embedded into an Epoxy matrix. Effective thermal properties for Graphite/Epoxy are shown in Figure 6. As in the previous experiments, the values obtained from the quadratic FEM formulation are closer to the semi-analytical values than those obtained by using the linear FEM formulation. In addition, the figure illustrates that there is little variation between the results obtained when performing the calculations using the square and hexagonal distributions.
From the calculated effective coefficients, it is possible to find the effective Young's modulus and Poisson's ratio. The linear FEM approximation was used to compute the effective E 4 a and ν 4 31 values. In Table 5, the E 4 a and ν 4 31 values are compared with those given in ref. [23] (E a and ν 31 ).
The effective thermal expansion coefficients values for the Graphite/Epoxy composite material have been reported in ref. [23]. Figure 8 shows a comparison between the calculated values of α by using the 4-node FEM implementation (square marks) and those reported in ref. [23] (dot marks). Note that if the material is constituted by a pure matrix material, the effective thermal expansion values correspond to the thermal expansion coefficient of the matrix.

Elliptical Fiber Cross Section
In this section, the thermoelastic effective coefficients are calculated for a composite material formed by fibers with an elliptical cross section embedded into a matrix with a square distribution as illustrated in Figure 9a. Figure 9b shows the periodic cell Y for the material described and the quarter of the periodic cell is indicated with dashed lines. The ellipse is defined by the length of its semi-axes, denoted by a and b in Figure 9b. It is possible to calculate the volume fraction occupied by the fiber using the following relation γ = abπ. The isotropic compound formed by Nicalon fibers in a BMAS matrix was used to perform the calculations. The values of the elastic and thermal properties of the fiber and the matrix are listed in Table 2. The values summarized in the following tables show the results calculated through the model presented in this work. The first column shown in each table corresponds to the values of the fiber volume fraction (γ) used during the numerical simulations. The results obtained by applying the linear and quadratic FEM are shown for each effective coefficient. Table 6 shows the elastic effective coefficients calculated by using linear and quadratic FEM approximations for different values of γ. The calculations were performed using an aspect ratio for the ellipse of a /b = 2. Finally, the differences between the different approximations used in this work are shown in each table.  Table 7 shows the calculated values of the thermal effective coefficients (thermal expansion and thermal conductivity) for the Nicalon/BMAS composite material formed by elliptical fibers. Additionally, the difference between FEM approximations are also illustrated. On the one hand, composites with circular cross-section fibers have the same values for the effective coefficients C 2323 and C 1313 . On the other hand, when the composite material is made with elliptical cross-section fibers, C 2323 and C 1313 do not have the same value. Similarly, the thermal expansion coefficients α 11 and α 22 and thermal conductivity coefficients κ 11 and κ 22 have the same values for composites with circular cross-section fibers and different values for composites with elliptical cross-section fibers, as observed in Tables 6 and 7. The aforementioned difference between effective coefficients is expected because the fiber is not evenly distributed in the y 1 and y 2 directions. To study the effects of the fiber geometry on the values of the effective properties, several experiments were performed. The fiber volume fraction value γ = 0.1 was kept constant, and the aspect ratio of the fiber was changed. The shape of the fiber varied by considering zero eccentricity values where a /b = 1 to large eccentricity values where a /b = 7. Figure 10 shows the effective thermoelastic properties calculated through both FEM implementations: linear and quadratic. Linear values are shown using dots and quadratic values are illustrated using square symbols. Each subfigure represents two different effective coefficients.

Distributions of the Fiber's Cross Section with Different Orientations
In this section, the effect of the elliptical cross-section orientation on the effective thermoelastic properties is analyzed. The calculated results of the effective coefficients are closely related to the fiber geometry, as mentioned in the previous section. For the elliptical cross section, it has been shown that C 2323 and C 1313 , and κ 11 and κ 22 , do not have the same values. Square periodic cells composed of thirty-six ellipses were considered to estimate the relation between the effective properties of the composite material and the orientation of the elliptical fibers. Consequently, the cell quarter must have a square with nine ellipses, as shown in Figure 11.
Effective thermoelastic properties were calculated using an aspect ratio of a /b = 2 and γ = 0.2. Figure 12a shows the effective coefficient values C 1313 and C 2323 for different orientations θ. The effective coefficient C 1313 (C 2323 ) increases (decreases) monotonically as a function of θ until C 1313 has the exact same value as C 2323 for θ = 45 • , where the composite material shows a transversely isotropic behavior as illustrated in Figure 12a. The thermal conductivity coefficients κ 11 and κ 22 show a similar monotonic behavior as a function of θ, reaching the expected transversely isotropic behavior for θ = 45 • , as observed in Figure 12b.

Computational Implementation
A MATLAB program was developed to calculate the thermoelastic effective coefficients by using square or hexagonal cell distributions. Since the representative cell has symmetric geometry, calculations are performed by using only 1 /4 of the cell. This allows to reduce the number of operations and CPU time.
The program receives as input the type of cell that constitutes the composite material as well as the elastic (Young's modulus, Poisson's ratio) and thermal properties (conductivity and thermal expansion) of the materials that correspond to the fiber and the matrix. The fiber volume fraction values γ, for which the calculations of the effective properties will be performed, are selected. For each value of γ a mesh that discretizes the quarter of the cell is generated. The method described in ref. [30] was used to construct an 8-node quadrilateral mesh. The 4-node quadrilateral meshes are created by using the obtained 8-node quadrilateral meshes. The effective thermoelastic properties, such as the elastic coefficients C, the thermal expansion α and thermal conductivity κ coefficients, are obtained for each value of γ. Linear and quadratic FEM formulations are used to solve the local problems obtained from the AHM. The generated meshes are used to perform calculations through the linear and quadractic FEM formulations, respectively. Graphics are created and tables are estimated for each value of γ, where the numerical values of the thermoelastic effective coefficients are shown.
As a result of applying the AHM and in order to obtain the thermoelastic effective coefficients, the local problems must be solved for each value of γ. The values of γ vary from 0 to 1 with a step size of 0.05 as a default value. For each local problem, every element within the mesh is analyzed. The local matrices of each element are calculated and the global matrices are assembled. The dimension of the matrices is twice the number of nodes in the mesh, for the elastic plane local problems and the pp T local problem. The dimension of the matrices is the same as the number of nodes in the mesh, for the elastic anti-plane local problems and the q I local problems. The matrices have several zero values, and due to the large volume of calculations that must be performed with large matrices where a large number of zeros is present, the MATLAB's sparse data structure was used. The sparse data structure stores only the positions and the values of the elements that are not zero, occupying a significantly less amount of memory. Additionally, MATLAB redefines the sparse array operations so that only nonzero elements are considered. Finally, approximately 30% of the execution time can be saved through this strategy. Figure 13 shows a flowchart that describes the computational implementation steps.

Conclusions
The generalization of the SAFEM for thermoelastic composite materials formed by elliptical cross-section fibers with perfect contact between the interface and the matrix was proposed. Square and hexagonal cell distributions were considered. Effective properties were obtained for three different composite materials. The calculated results by using the linear and quadratic FEM formulations are very similar to the semi-analytical values reported by other authors. Additionally, the calculated effective properties for the square and hexagonal distributions are very similar. The effective coefficients for a composite material with several orientations for an elliptic cross section fiber show a non-isotropic behavior according to the results of this work. Additionally, as the angle of orientation increases, the behavior of the composite material approaches smoothly to a transversely isotropic behavior.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.
Appendix A Table A1. Boundary conditions for local problems.