Finite thermo‐elastic decoupled two‐scale analysis

Engineering structures are often characterized by different structural properties, depending on the length scale under consideration. Fiber reinforced composites are determined, eg, by a heterogeneous microstructure, but are sufficiently described by homogeneous characteristics at their macroscopic level. Moreover, different loading situations, eg, of thermal or mechanical nature, require the consideration of multiphysical equilibrium states. The challenging engineering task is the computation of the effective material properties of these different loading scenarios. The contribution at hand introduces a finite thermo‐elastic two‐scale analysis, where the effective macroscopic material properties are computed in a decoupled manner with respect to the different length scales.


INTRODUCTION
A numerical analysis of structures or substructures with respect to their behavior due to any external load applied is often required in engineering tasks. Utilizing numerical methods and approaches enables to tailor the geometry and the properties of the structure. Mechanical or thermal loads are often dominating. An interaction of the thermal and mechanical phenomena is not negligible for many engineering tasks. Especially, structures at large deformations and thermal conduction or dissipative phenomena require the consideration of heat generation and convection inside or at the surface of the structure. For example, synthetic materials, such as polymers, are significantly deformable. These materials can be reinforced by particles or fibers, in order to improve their strength or stiffness for certain purposes during their operating life. Reinforced material consisting of different bulk materials is heterogeneous. In specific cases, it is anisotropic. Considering a structure, eg, bearings, car bumpers, or damping devices, made of reinforced synthetic materials, the reinforcement is usually characterized by a microscopic length scale.
In general, homogenization methods are used in continuum mechanics to investigate heterogeneous materials or composites. Homogenization is motivated by large length scale differences of composite structures between their microstructure and their macrostructure. Structures made of fiber reinforced polymers are typical composites, see Figure 1 for an example. The fiber material and the geometry at the microscale as well as its orientation influence the structural properties and behavior. The macroscale can be considered as homogeneous, but both scales are linked to each other. The interaction of the material structure at different scales of a heterogeneous body can be evaluated by different approaches. Especially, computational homogenization is one of the most suitable scale bridging methods to predict structural behavior based on a heterogeneous microstructure.
Knowing the construction and geometry of a representative volume element (RVE) of the microstructure, a common way to predict the structural behavior of the macroscale is to determine effective properties of such an RVE. Having obtained these properties, the constitutive behavior of the macroscale is computed. Another way of investigating the macrostructure is to scale down to the representative part of the microscale at the position of the macroscale, where constitutive properties are needed and therewith to decompose the macroscale locally. After localization and determination of the microscopic structural response, a homogenization procedure scales the information up again to its related macroscopic position.
An important aspect of homogenization methods and their prediction accuracy is the determination of the structure of the representative microvolume as well as its spatial extent. Following the work of Hill 2 or the work of Kröner, 3 an RVE represents a part of the microstructure, which leads to the equivalent macroscopic properties independently of boundary conditions or location. Generally, uncertainty investigations (considering, eg, stochastics) lead to the determination of the representativity of an RVE with respect to its size or heterogeneity. Another important aspect of the concept of RVEs is the length scale difference of the spatial extent of the RVE with respect to the spatial extent of the macrostructure. As stated in the work of Hashin, 4 the ratio of both dimensions can be used as a measure for the scale separation and should fulfill the condition → 0. Separated scales enable the pointwise linkage of different scales.
Common homogenization methods are introduced, eg, in related works. [4][5][6] Further approaches to homogenization that make use of upscaling and downscaling of the properties and driving boundary conditions are given in, eg, other works. [7][8][9][10] In recent years, remarkable progress has been made in the field of computational homogenization in structural engineering, see, eg, related works [11][12][13] and the references therein for a review. The main challenge of computational homogenization is still the requirement that the cost of numerical homogenization solution should be less than the cost of the numerical solution on the fully discretized structure. Various solutions have been proposed in literature in order to derive homogenization methods for purely mechanical problems. Temizer and Wriggers 14 developed an adaptive method based on a material map from the microscale to the macroscale. This database methodology is used in order to determine the effective macroscopic behavior for orthotropic nonlinear elasticity.
Most approaches are based on the key assumption of separate scales, which means that the volume of the microscale tends to zero, viewed from the macroscale. Guedes and Kikuchi 15 and Ghosh et al 16 derive effective properties of composite microstructures for linear elasticity by averaging them over the microstructural domain. They also use the homogenized quantities for the macroscopic stress analysis within an finite element method (FEM) framework. Further approaches, see, eg, other works, 7,17-23 use a microscopic fine scale fluctuation field in order to determine the microstructural deformation state of the mechanical heterogeneity. This microstructural state determination is based on certain assumptions for the boundary of such a microstructure. Once the microstructure has reached a mechanical equilibrium, stresses and stiffness for the related material point within the macrostructure are computed, based on the prescribed boundary information. Higher-order approximations for the scale-bridging deformation gradient are introduced and used in related works. 8,24,25 Similar to the approaches aforementioned, the microscopic fluctuation field needs to be determined by assumptions, which allow the microscopic equilibrium to be achieved. In the work of Temizer and Wriggers, 26 the computation of the macroscopic tangent in terms of the microscopic fluctuation field is investigated and a procedure for determining the stiffness from a stress controlled testing method at the microscale is developed.
Multiphysical computational homogenization methods are under development. A few publications try to bridge the scales for thermo-mechanical problems. In the work of Özdemir et al, 9 a staggered algorithm is introduced, which determines the macroscopic stresses and heat fluxes from the average over the microscale within a nested FEM procedure. A model reduction technique combined with a homogenization approach is used in the work of Monteiro et al 27 in order to solve thermal and electric conduction problems efficiently and precisely. This solution is based on periodic microstructures as well as orthogonal decomposition methods. In the work of Temizer and Wriggers, 10 macroscopic quantities are identified by averaging the appropriate microscopic quantities separately over the microstructure in finite thermo-elasticity. In the work of Zhang et al, 28 steady-state microscopic thermo-dynamic wave propagation equations are derived, by means of an asymptotic multiscale homogenization approach from higher-order methods in homogenization. Furthermore, using the asymptotic expansion for thermo-mechanical homogenization, in the work of Temizer, 29 two uncoupled cells are obtained at the microscale, one for the mechanical field and one for the thermal field, which allow an efficient numerical solution by ensuring thermo-dynamical consistency. In the work of Fleischhauer et al, 30 another approach is introduced to multiphysical homogenization based on the assumption of continuity of the appropriate kinematic fields across both scales and the homogenized balance equations that drive the considered field quantities. The formulation is applied to the fully coupled partial differential equations of thermo-mechanics.
Not only multiphysical homogenization strategies have been developed in recent years, multiscale approaches on cracks and delamination phenomena are also under investigation. In the work of Fish and Yu, 31 a nonlocal theory for damage of brittle composites is developed by applying an asymptotic expansion in order to obtain the overall macroscopic properties. The computation of effective macroscopic properties based on certain mesoscopic structures incorporating damage and fracture processes of concrete is proposed in the work of Wriggers and Moftah. 32 Recently, in the work of Matouš et al, 33 a multiscale cohesive method is developed for coupling the microscale damage evolution to the resulting macroscopic properties. A scale bridging method for a cohesive interface between the macroscale and an RVE at the microscale is proposed in the work of Hirschberger et al, 34 where the macroscopic properties are computed via homogenization. In the work of Holl et al, 35 a homogenization technique for predicting crack propagation and coalescence of cracks at the microscale and macroscale contributes to the extended finite element framework.
The paper at hand is structured as follows, in order to consistently introduce the finite thermo-elastic decoupled two-scale analysis. Section 2 depicts the kinematic basis considered. Section 3 summarizes the stress and thermal measures used, and in Section 4, the thermo-elastic equilibrium conditions as well as the necessary constitutive descriptions are introduced. The main development is presented in Section 5, where effective quantities, energetic volume averages, as well as the related Hill-Mandel criterion and the decoupled approach are introduced. A representative numerical example in order to verify and validate the developments is shown in Section 6. Finally, conclusions, an outlook, and appendices are given.

KINEMATIC BASIS
The subsequent section introduces the kinematics of a homogeneous material body, deformation measures and their splits, as well as the absolute temperature.

Configuration, motion, and deformation gradient
Consider a body B, composed of material points P ∈ B that are identically connected to the domain  in the Euclidean space R 3 . The configuration is a bijective mapping of P of B onto a position x ∈ R 3 . The motion t is a function of configurations, where t ∈ R denotes the time. t can be expressed as Observing a material point P, Equation (2) describes the path of P in R 3 while B moves. Equation (3) introduces the current or spatial configuration  t at current time t. At a fixed time t 0 , the body represents the configuration  0 , which is generally denoted by an undistorted, stress-free state, with a homogeneous temperature distribution. Therefore,  0 is commonly named reference or material configuration and reads Equation (5) introduces the reference position or Lagrangian coordinates X of P, while x of P denotes the current or Eulerian coordinates. Considering now only the Euclidean space and eliminating the variable P from the set of Equations (1) to (6), the nonlinear deformation map (X, t) = x, is defined. maps the reference coordinates X ∈  0 of P onto their current coordinates x ∈  t . For a schematic two-dimensional depiction of the configurations and their mathematical dependencies in R 3 , see Figure 2. The displacement vector can then be expressed as u ∶= x − X.
The deformation gradient is used for pull-back and push-forward operations, in order to transform physical quantities between different configurations. The identity J ∶= det F is introduced. For later use, the spatial velocity gradient is defined. Based on the deformation gradient, it can be rephrased as

Temperature change and absolute temperature
The absolute temperature is described as where > 0 is the absolute temperature of a point P ∈  t measured in Kelvin [K]. The change of temperature with respect to the reference temperature 0 is denoted as . Equation (6) depicts the reference configuration (usually) with a homogeneous temperature distribution 0 .

Volumetric and isochoric split of the deformation gradient
As outlined in the work of Kaliske, 36 the split of the deformation gradient into an isochoric and a volumetric part, implies the decomposition of the total deformation into a volume preserving and volume changing part. This assumption is locally valid in the infinitesimal vicinity of the material point P under observation. Observing one can conclude that any change in volume of  t is only caused by F vol .

Multiplicative thermo-elastic split of the deformation gradient
The following multiplicative split of the deformation gradient F is used, in order to describe the deformation behavior at the macroscopic level. Based on the assumptions made in the work of Terada et al, 37 the multiplicative decomposition of F into a dilatoric and into a volumetric part is used, where J = det F. At the same time, a decomposition of F into a thermal expansion part F and a stress-inducing elastic part F e is carried out, similar to the decompositions in other works. [38][39][40] The thermal expansion part can be specified by for an anisotropic material. The scalar-valued entries in the anisotropic thermal coefficient tensor depict the thermal expansion coefficients in the appropriate directions. The temperature change = − 0 with respect to a reference temperature 0 is introduced in Equation (13). By Equation (17), it follows that the stress-inducing elastic part of the deformation gradient takes the form

STRESS MEASURES AND THERMAL QUANTITIES
The following section introduces stresses, heat flow, and heat flux vectors.

Stresses
The Cauchy stress theorem defines a stress vector t at the surface  t of the current configuration and, thus, postulates an infinitesimal surface force tda according to t ∶= n , where are the Cauchy stresses and n is the outward normal at the material point x ∈  t . For a schematic illustration of the current surface force in the reference and current configuration, see Figure 3. The configuration  t can be any subdomain of the real body under observation. Subsequently, is defined for any point x ∈  t and denotes the true stresses, since it relates the stresses of the current configuration to the area element of the same configuration. Connecting the current stresses to the surface of the reference configuration  0 leads to the first Piola-Kirchhoff stresses P nda = JF −T NdA =∶ PNdA. For later use, the first Piola-Kirchhoff stresses are defined according to

Thermal quantities
The heat flow q n out of the surface  t of the current configuration can be expressed as where q denotes the spatial heat flux through point x ∈  t and n the current outward normal at the observed point. Equation (24) can be used to postulate the infinitesimal transport of thermal energy per time through point x ∈  t according to Within Equation (25), Q denotes the material heat flux, which depends on q through Figure 3 depicts the current and reference heat fluxes schematically. The spatial heat flux vector, describing the ability of the material to conduct heat, is assumed to follow Fourier's law in case of isotropic material where k > 0 represents the thermal conduction coefficient of an isotropic conductive material behavior. Next, the anisotropic heat flux vectors at the reference and current configuration are introduced as compare the work of Miehe. 38 The thermal conductivity coefficients k ij in the appropriate coordinate directions are introduced for the case of anisotropic material, reading where the vectors A and B are the direction vectors for the orthotropic material characteristics. Accordingly, the thermal expansion coefficients ij in Equation (18) read

THERMO-ELASTIC EQUILIBRIUM AND CONSTITUTIVE SPECIFICATIONS AT FINITE DEFORMATIONS
In order to simulate finite thermo-elastic problems via the FEM, the following underlying coupled problem has to be solved. The local formulations read Equation (32) is the quasi-static equilibrium condition with respect to the current configuration, resulting from the local balance of linear momentum, and Equation (33) is the thermo-elastic part of the transient heat conduction equation, following from the local balance of energy. For details of deriving this set of equations, the reader is referred to the work of Miehe. 38 Equation (32) depicts the evolution equation for the displacements u, and Equation (33) is responsible for the evolution of the temperature change .

Anisotropic thermo-elastic material model
Based on the developments made in the work of Terada et al, 37 the material model proposed in the work of Kaliske 41 is used and extended to finite thermo-elasticity. The Helmholtz free energy at the macroscale is extended to thermo-elasticity and takes generally the form The temperature dependency of the stresses is included into the volumetric part of the free energy, while the isochoric part is chosen to be independent of temperature, in order to prevent nonvolumetric deformation due to excitation by the temperature field. The specific constitutive formulations read The function C( ) is given in order to preserve the heat capacity at constant deformation, defined by c ∶= − 2 , since the volumetric part of the free energy is nonlinear with respect to temperature. Determining the second derivative of the Helmholtz free energy function, the corrector function C( ) has to ensure the result c = − 2 . No explicit computation is necessary within this framework, since the coupling between temperature and displacement is assumed within W vol (J e ( )). Equation (35) shows that the specific choice of the kinematic dependency of the volumetric energy part is equivalent to other approaches to consider thermo-mechanical coupling, see, eg, other works. [42][43][44][45][46] The invariants used in W iso (C) are defined in the work of Kaliske. 41 What remains is the specification of the Cauchy stresses. Due to the split of the free energy, see Equation (34), the Cauchy stresses are derived by They are a function of the current metric tensor g. Here, g is the covariant representation of the symmetric metric of the current configuration and the right Cauchy-Green deformation tensor C = F T gF denotes its pull-back to the reference configuration.
Further expansion leads to In order to transform Cauchy stresses to first Piola-Kirchhoff stresses, the standard transformation of Equation (23) can be applied. Equation (33) requires the computation of the following formulas: Their derivatives, necessary for an FEM-implementation, have a general form of introducing the current tangent moduli ℂ = 1∕J g ( J ). Furthermore, the sensitivity m of the Cauchy stresses with respect to the temperature is introduced, computed by The different linearized external power terms with respect to g, u, and are named w g , w u , and w t and they are derived by Last, heat storage is a function of temperature only and needs to be linearized according to The detailed derivation and the explicit formulations are given in Appendix A. Concluding the macroscopic material model, a full list of material constants to be identified according the numerical material testing procedure, introduced in the work of Terada et al, 37 is given in Table 1.

Isotropic thermo-elastic material model
Basis for the homogenization is the heterogeneous unit cell, where the different materials are discretized. The proposed thermo-elastic formulation used here is based on the developments in the work of Miehe 38 and is briefly summarized by the following equations. Similar as the kinematics introduced in Section 2, the determinant of the thermal part of the deformation gradient is constitutively described by where t is the isotropic thermal expansion coefficient and the isochoric part of the deformation gradient is given bȳ The Cauchy stresses are derived from the finite thermo-elastic stress potential from the work of Miehe 38 and yield the form where and are the bulk modulus and the shear modulus, respectively.b =FF T is the isochoric part of the left Cauchy-Green deformation tensor. The temperature sensitivity of the Cauchy stresses is given by (60) Additionally, the heat capacity quantitiesĉ are given in order to complete the necessary terms for evaluating the thermo-elastic equilibrium. Their derivatives, necessary for FEM-implementation, follow schematically Equations (50) up to (55). Their specific formulations are given in Appendix B. Concluding the isotropic Neo-Hookean material model, a full list of material constants is given in Table 1.

HOMOGENIZATION AT FINITE THERMO-ELASTICITY
Generally, macroscopic systems might be constituted by mesoscopically or microscopically heterogeneous material structures. At different length scales, different structural properties are observed in such cases, eg, the macroscopic level can be considered as homogeneous while the mesoscopic or microscopic scale consists of different materials. In order to computationally model such systems at the engineering (macroscopic) length scale, homogenization is an adequate method. The following section explains the procedure to identify the parameters given in Table 1 of the microscopic structural behavior for an arbitrary unit cell and introduces the appropriate homogenization procedure with respect to finite thermo-elasticity. The proposed approach given in the following is developed first in a consistent manner for fully coupled FE 2 methods, and subsequently, the two length scales are decoupled in order to preserve efficiency for macroscopic structural analyses.

Effective thermo-elastic quantities
At a point P of a solid body B with a position X ∈  0 in the reference configuration  0 , effective quantities are introduced based on the developments made in the works of Miehe 7 and Bayreuther. 47 Thus, the decomposition of the RVE boundary  0 into parts of prescribed displacements  0 u , parts of prescribed surface traction vectors  0 T , and parts of prescribed temperatures  0 is important. Different combinations of the prescribed boundaries are possible, A schematic representation of the boundary decomposition is given in Figure 4. First, effective deformations are defined and transformed into surface data by Furthermore, an effective first Piola-Kirchhoff stress measure is given and transformed into an RVE surface integral. Basis for the dual mechanical field variables P and F within this homogenization framework is the Helmholtz free energy 0 (F, ), since Having a look at the thermal field, no analogous dual variables can be defined with respect to the dependency of the Helmholtz free energy on the absolute temperature . Following the works of Coleman et al, 48,49 the Helmholtz free energy cannot directly depend on the temperature gradient  = GRAD( ), which is a direct consequence of the Clausius-Planck inequality, regarding dissipative effects. Thus, the derivative is zero. Considering this result, no dual thermal variable to the temperature gradient  exists. The direct dependency of the Helmholtz free energy on the thermal field variable and knowing the definition of the entropy = − (F, ) of a point in a thermo-dynamic system implies the duality of and . However, there is no effective entropy to be defined for an RVE in a generalized thermo-mechanical sense, since the entropy evolution for reversible processes depends on heat storage as well as on thermal expansion effects and, for irreversible processes, additionally on dissipative mechanisms. Nevertheless, an effective entropy production is defined, in order to ensure the positiveness of the second law of thermo-dynamics at both scales. The reference representation of the local form of the thermo-dynamic restriction with respect to the specific entropy production reads where internal heat sources are neglected. A further expansion, a manipulation by multiplying with the strictly positive absolute temperature and the incorporation of Equation (69) can be further recast into by using Equation (65) and the definition c = − 2 . At this point, volume averaging is carried out, in order to ensure that the entropy production of the RVE is identical to the entropy production at the appropriate local point of the macroscale. The effective form of the second law of thermo-dynamics reads Equation (71) determines the effective thermal material behavior at the macroscopic length scale. Furthermore, an effective energy can be defined as a volume average of the energy of an RVE.

Energetic volume averaging and the generalized Hill-Mandel criterion
Starting point for the subsequent developments is the Legendre transformation of Gibbs' equation that defines the relation of the Helmholtz energy , the internal energy e, and the entropy of a thermo-dynamical system, according to Following other works, 48-50 the mass-specific internal energy density changes according to the mechanical and the thermal power, reading where Equation (68) and the specific heat capacity c = − 2 are incorporated. Assuming a pointwise linkage of an RVE to a macroscopic point P ∈  0 , one can define the volume average of the rate of the internal energy of the RVE equal to the rate of the internal energy density 0ė by Under consideration of Equation (75), Equation (76) can be recast into Therewith, the transition from microscale to macroscale is achieved.

Hill-Mandel criterion
The linkage of the two scales and the necessary boundary conditions for the field variables x and is achieved by the proposed generalized thermo-elastic Hill-Mandel criterion The developed thermo-elastic homogenization boundary conditions are schematically depicted in Figure 5.

Dirichlet-type boundary conditions or linear displacement boundary conditions
In order to derive displacement and temperature boundary conditions for the RVE, Equation (64) is substituted into the right side of Equation (78), and Equation (77) is incorporated into the left side of Equation (78), which leads to To formulate surface integrals on the left side of Equation (80), further manipulations are performed, reading An assignment of terms as a requirement to compare equivalent terms between left side and right side of Equation (80) The requirement to prescribe the whole RVE boundary by displacements and temperature in order to achieve equality in Equation (83) and Equation (84). A closer look onto Equation (86) leads to the result of linear deformations at the boundary  0 of the RVE driven by the macroscopic deformation gradient F, compare the work of Miehe 7 and Figure 5A, since no microscopic displacement-type fluctuations are allowed. Imposing the macroscopic temperature on  0 (see Equation (87)) ensures the identity together with Equation (86). What remains is to fulfill for the volume of the RVE. First fact to notice is that Equation (89) cannot be expressed as a surface integral. Due to the nature of the heat storage effect, this is reasonable, since only convective or insulated phenomena are present at the surface of an RVE. A second fact to observe is that volume averaging of the heat capacity would lead to a strict requirement of a unique temperature distribution equal to the macroscopic temperature over the whole domain of the RVE. This is due to an equality of the average of the product of two functions with the product of two averages of the appropriate functions. Having achieved steady state for the evolution of the temperature, ie, . = 0 ∀ X ∈  0 and . = 0, the effective macroscopic heat capacity cannot be computed. Therefore, the triviality of the solution 0 = 0 in Equation (89) for steady states has to be avoided. At temperature rates that are not zero, the effective macroscopic heat capacity is obtained by integrating Equation (89) in time. Integration by parts leads to for the assumption of constant material properties in time. By Equation (91), one can compute the effective macroscopic heat capacity by implying the nonzero property > 0.

Neumann-type boundary conditions or uniform traction boundary conditions (UTBC)
In order to derive traction and temperature boundary conditions of the RVE, Equation (63) is substituted into the right side of Equation (78) and Equation (77) is incorporated in the left side of Equation (78), which leads to A further manipulation with respect to the effective deformation gradient on the right side and the surface transformations from Equations (81) and (82) on the left side of Equation (93) yields The requirement to prescribe the whole RVE boundary by traction vectors and temperatures in order to achieve equality in Equations (94) and (95). A closer look to Equation (97) leads to the result of homogeneous traction vectors on the boundary  0 of the RVE driven by the macroscopic first Piola-Kirchhoff stresses P, compare the work of Miehe 7 and Figure 5C. To prescribe the macroscopic temperature on  0 ensures the identity together with Equation (97). In order to ensure Equation (96), Equation (91) has to be applied.

Perodicity-type boundary conditions or periodic displacement boundary conditions
Next, periodic displacement and temperature boundary conditions are derived. In order to excite the unknown field quantities of the RVE, displacement jumps, regarding opposite points P + ∈  + 0 and P − ∈  − 0 at opposite surfaces  + 0 and  − 0 of an RVE, are introduced. Thus, the difference between x + − x − = ⟦x⟧ ∀ x + ∈  + 0 ∧ x − ∈  − 0 is computed, compare Figure 5B, with the property N + = −N − for the normals at the appropriate reference positions. Subsequently, Equation (64) is substituted into the right side of Equation (78), which leads to A further manipulation with respect to the displacement jump and the effective first Piola-Kirchhoff stresses as well as the requirement to prescribe the whole RVE boundary by displacements jumps and temperature in order to achieve equality in Equations (102) and (103). A closer look to Equation (105) leads to the result of periodic deformations on the boundary  0 of the RVE, driven by the macroscopic deformation gradient F, compare the work of Miehe 7 and Figure 5B. To prescribe the macroscopic temperature on  0 (see Equation (106)) ensures the identity in Equation (103)

Decoupled algorithmic treatment of homogenization considering the different boundary conditions
Following the work of Terada et al 37 and the developments made therein, the following section explains the proposed approach of a decoupled two-scale analysis in case of finite thermo-elasticity. Figure 6 depicts the decoupled two-scale analysis schematically.

Preprocessing of RVE-data
In the first step of the decoupled thermo-elastic homogenization approach, preprocessing regarding the input data of the heterogeneous RVE is carried out. Data of the RVE are read and stored for later use. Furthermore, nodal coordinates as well as element information, such as nodal connections and material parameter information are part of the input.

Numerical material test
The following section describes the numerical material tests. According to the boundary conditions applied, either uniform traction vectors (UTBC), linear (LDBC), or periodic displacements (PDBC) and their loading function, the homogenized quantities of the heterogeneous RVE response are computed. The deformation F is prescribed, such that mixed loading cases are applied to the RVE. The values F i are separated according to the three spatial coordinate directions. Therefore, the first mixed loading case is defined by F 1 | = 1, … , 3 (loading in x 1 direction) and the two other cases (x 2 and x 3 ) accordingly. The three loading sets are separately applied to the boundary data of the RVE. The advantage of prescribing mixed loading situations is capturing of coupled thermo-elastic RVE responses for the anisotropic mechanical, thermal, and thermo-mechanical effects inside the RVE. Similarly, mixed loading cases are constructed for the traction vectors, if P is prescribed.
In case of prescribed boundary temperatures, the temperature changes at opposite sides of the heterogeneous RVE are prescribed for separate loading in x 1 -, x 2 -, and x 3 -directions.
Equation (108) is basis for computing effective elastic constants and properties, respectively, at the macroscale. In order to ensure Equation (71), reading each of its terms is evaluated separately, according to The volume averages are stored for determining the effective specific heat capacity by Equation (110), the effective thermal expansion coefficients by Equation (111), and the effective thermal conductivity coefficients by Equation (113). Equation (112) cannot be used as an indicator for the thermal conductivity behavior, since the heat flow per unit surface area Q n is independent of the thermal conductivity coefficients. The general form of the algorithm to apply the boundary data is based on the differentiation between inner degrees of freedom x a and outer or prescribed degrees of freedom x b . The general FE-residual is given by as difference between internal and external contributions of the unknown fields Consistent linearization of Equation (114), according to with respect to the nodal unknowns in Equation (115), leads to In order to generalize the algorithmic treatment of the three different boundary conditions, the increments of the nodal unknowns are recast into where R is a reduction matrix with respect to the nodal unknowns Δx to be computed and K is a constraint allocation matrix with respect to the constraints per unknown degree of freedom in C.

Specification for linear displacement at the RVE boundary (LDBC)
According to a full prescription of the boundary data, Δx takes the form and, consequently, R is constructed as Equations (86) and (87) require to construct C as where H = F − 1 is the displacement gradient of the prescribed deformation F and b is the prescribed boundary temperature. What follows is the construction of the constraint allocation matrix K according to

Specification for periodic displacements at the RVE boundary (PDBC)
According to a separation of the data into partner degrees of freedom at positive  + 0 and negative  − 0 boundaries, a detailed expansion of Equation (115) leads to Consequently, the increments of Equation (123) read In order to satisfy Equations (105) and (106), the unknown degrees of freedom are computed and, consequently, R is constructed as Equations (105) and (106) furthermore require to construct C as where H = F − is the displacement gradient of the prescribed deformation F and b is the prescribed boundary temperature. Finally, the constraint allocation matrix K in case of PDBC is given by

Specification for uniform traction vectors at the RVE boundary (UTBC)
In case of UTBC, nodal forces at the boundary  0 T are constructed according to Equation (97) by and assemble, under consideration of the inner RVE domain, the external forces F u ext for the RVE. F u ext read compare Equation (114). Consequently, the increments of the nodal unknowns Δx with respect to Equation (115) are given by and the reduction matrix R is constructed as Since only the temperature degrees of freedom at the boundary of the RVE are prescribed, C reads and K in case of UTBC is given as Since all the mechanical degrees of freedom are not prescribed and no fixed boundary conditions exist, artificial rigid body motions of the RVE are prevented by construction of a mass-type perturbation vector R and matrix M for every finite element having a surface at the boundary of the RVE. Both artificial contributions read and where N I = N I ( ) denotes the shape functions. The perturbation value in Equations (135) and (136) is chosen to be a small number, compared to the stiffness matrix of the appropriate element. R is added to the internal element forces and M is added to the mechanical part of the element stiffness matrix.

Finite element solution for LDBC, PDBC, and UTBC
In order to solve for the structural response of the RVE by using the FEM, Equation (118) is inserted into Equation (116), which yields After manipulation of Equation (137) by R T the transformed FE problem is given by and it is solved for Δx according to the appropriate boundary conditions, eg, by standard Newton-type solvers, in an iterative manner. After the solution of Equation (139), a back-transformation for all the real unknown nodal values of the RVE has to be applied, according to Equation (118), until equilibrium is achieved.

Macroscopic material parameter identification and macroscopic homogeneous material model application
After having computed the data from numerical material tests, the set of macroscopic material parameters in Table 1 has to be identified. Therefore, appropriate optimization techniques need to be applied, see, eg, the work of Amaran et al 51 for an overview. Common optimization methods are, eg, minimization of squared errors by gradient based approaches, see, eg, the work of Fischer, 52 or evolution strategies, see, eg, the work of Deb and Tiwari. 53 The optimization approach applied within this framework is gradient based minimization of the squared errors. The reason for this application is the efficiency of gradient based approaches compared to evolution or searching methods. Although only local minima are eventually computed, the efficiency is a major advantage of this approach. Gradient optimization methods are depending on the starting point in the parameter space. This property might lead to the computation of a local instead of a global minimum.
Here, the error functions e I ( i ) consist of the differences between the heterogeneous RVE computations (the numerical material tests for the three loading cases, compare Figure 6) and the homogeneous sample computations. The argument i represents the respective model parameter of the appropriate error function. Superscript □ I depicts the evaluation point of the error e. In order to determine the mechanical model parameters of the homogeneous anisotropic material model, the qualitative errors are the differences of the homogenized first Piola-Kirchhoff stresses (in case of LDBC and PDBC) as well as the homogenized deformation gradient (in case of UTBC) at an evaluation point I, reading Subscripts □ ho and □ he depict a quantity due to a homogeneous or heterogeneous material consideration. For computing the effective heat capacity, the error function takes the form while the anisotropic thermal expansion coefficients are determined by The thermal conductivity coefficients are computed by Exemplarily, Equations (140) and (141) can be interpreted as a function of the effective bulk modulus , the effective shear modulus equivalents a i · · ·g i , and the effective thermal expansion coefficients ij . Therefore, the derivatives of the homogeneous sample are computed. In a general nonlinear material model, and P respectively depend nonlinearly on the deformation measure. The deformation measure itself again might depend on the elastic constants vice versa. Due to the fact that Equations (140) and (141) are computed in an equilibrium state for the homogeneous samples at loading state I, only the constitutive equations for and P respectively have to be derived with respect to the effective material properties.
Generally, the error function for a set of specific effective material constants i reads within the presented approach. The term SIM ho ( ) represents the homogeneous effective sample response as a function of the appropriate effective material parameter, while EXP he denotes the heterogeneous effective RVE response. The general form of the specific derivative is depicted, in order to be utilized within the gradient based minimization approach. The optimization follows the Taylor series expansion of Equation (149) around the root Having a look at Equation (152), it has to be noticed that the number of evaluation points N needs to be at least the number of parameters p, ie, N ≥ p. This necessary condition ensures that the system of equations in Equation (152) is not under-determined, what would mean that an infinite number of solutions for the set of parameters is present. Typically, N is larger than p. Hence, an overdetermined system of equations has to be solved. Within this framework, a solution based on a singular value decomposition of the coefficient matrix K in Equation (151) is applied, see the work of Golub and Reinsch. 54 Based on the nonzero singular values of K, an optimal solution point of Equation (152) for is computed, satisfying all equations in an averaged manner, see the work of Schwarz. 55 Equation (152) is solved for until ‖ ‖ = 0. The iterative solution approach depends on the initial parameter values and might yield multiple roots. The solution of Equation (152) is found such that an admissible data range for all j is considered. This is necessary in order to prevent computation of a root yielding, eg, a negative heat capacity or any other nonadmissible parameter values. After the computation of all effective material parameters of the anisotropic effective material properties, such that the homogeneous sample has an equal structural behavior compared to the heterogeneous RVE, the macroscopic homogeneous structural applications can be computed based on the microstructural information at the RVE-level.

NUMERICAL EXAMPLE
Subsequently, the verification of the implementation as well as the validation of the formulation of the decoupled thermo-elastic homogenization approach is shown.

Structural thermo-elastic response of a fiber reinforced composite RVE at three different homogenization boundary conditions
This example firstly verifies the implementation of the constitutive equations and the finite element solver by investigation of the convergence behavior with respect to the iteration steps. Second, the overall mechanical response of the RVE according to Dirichlet-type, Neumann-type, and periodicity-type boundary conditions is compared. Figure 7 shows the RVE and its discretization and material distribution with respect to fiber orientation and matrix material. The RVE consists of 2280 linear temperature-displacement elements with 2651 nodes. The initial heterogeneous problem consists of the material distribution, shown in Figure 7A. The volume ratio between fiber and matrix material corresponds to a fiber The homogenized problem consists of the same discretization but of a homogeneous material distribution as shown in Figure 7B. The discretization of the homogeneous sample is identical to the heterogeneous RVE, in order not to introduce additional discretization errors into the computation. According to Figure 6, the RVE is initially investigated with respect to the numerical material tests at unit loading cases for the appropriate unknown fields. In the following, the stress and temperature results of the heterogeneous RVE computations at certain loading situations are used, in order to determine the material parameter of the anisotropic thermo-elastic constitutive laws, see Section 4.1. The identification is carried out by an optimization approach, as mentioned in Section 5.3.3. Table 2 depicts the material parameters used for the heterogeneous RVE. The loading function applied is given in Figure 8. In case of Neumann-type boundary conditions or uniform traction boundary condition (UTBC), the loading function for the mechanical field takes the form 1.5 · 10 4 6.0 · 10 3 2.6 · 10 −5 1.0 · 10 0 3.1 · 10 −1 Material 2 (matrix) 2.0 · 10 3 1.0 · 10 3 2.0 · 10 −4 1.5 · 10 0 1.5 · 10 −1 In case of Dirichlet-type or linear displacement boundary conditions (LDBCs) and periodic displacement boundary conditions (PDBCs), the loading function for the deformation of the boundary is

Numerical material test at the fiber reinforced RVE
Decoupling of the spatial loads in Equations (153) to (154) is carried out, in order to account for the three main deformation phenomena in the three spatial directions, ie, to capture the characteristics of the heterogeneous RVE in x 1 -direction, the first row of, eg, F is prescribed at the boundary of the RVE with l(t). This procedure ensures not to account only for single mechanical properties of the RVE with respect to, eg, F i , also combined (stretch and shear) RVE characteristics are activated. The temperature variation at the boundary follows the loading Decoupling of the temperature load at the boundary  0 is performed such that mechanical load in x i -direction is combined with thermal load = ∀ (X i = max X i ∈  0 ∨ X i = min X i ∈  0 ). The spatial decoupling of the thermal load is carried out in order to account for a clear distinction of thermo-mechanical coupling and thermal conductivity phenomena, similar to the decoupling of the spatial loads in Equations (153) to (154), compare the work of Terada et al. 37 The comparison of the homogenized quantities of Equations (107) and (108) as stress versus deformation plots is given exemplarily for the x 1 loading case from Equations (153) to (155) in Figure 9. As can be seen in Figure 9, the structural response of the RVE is in accordance with the limit cases, compare, eg, the work of Terada et al, 21 for a monotonic loading case, ie, t ≤ 300 s in Figure 8. The stiffest stress-strain response of the RVE is observed, if linear deformations (LDBC) are prescribed at the boundary of the RVE, while the softest structural response is obtained by uniform traction vectors (UTBC) as boundary conditions and periodicity constraints lead to a stress-strain behavior in between both limit cases. Due to the contrary loading effects with respect to the temperature field, the stress-strain-dependency for F ii and P ii reverses if a new loading occurs after the unloading path, ie, for t ≥ 600 s. The contrary loading effect results from the external or stress power consideration in Equation (33). Whenever negative stress is prescribed, the external power term in Equation (33) generates a positive temperature change, but the applied loading function l(t) enforces a negative temperature change at negatively prescribed stresses. Generally, all shear-stress and shear-deformation plots in Figure 9 fulfill the homogenization limit cases, since the coupling effect of volumetric thermal expansion is not present for shearing. The deformed configuration as well as the internal temperature distribution for loading case 1 (x 1 -direction) at t = 300 s is shown in Figure 10A. Figure 11 depicts the Cauchy stress and internal temperature distribution for loading case 2 (x 2 -direction) of the heterogeneous RVE and Figure 12 shows the appropriate quantities accordingly. Comparing Figures 10A, 11A, and 12A, it has to be noticed that the LDBC cases are denoted by the largest norm of the Cauchy  Figure 7A) according to the three different boundary conditions of homogenization for the depicted loading function at t = 300 s at loading scenario x 1 . LDBC, linear displacement boundary condition; PDBC, periodic displacement boundary condition; UTBC, uniform traction boundary condition stresses, while UTBC cases have the lowest Cauchy stress norm and PDBC lays between both limit cases. The high Cauchy stress norm in case of LDBC is originated in the completely prescribed boundary deformation of the RVE, which prevents the structure from reducing the stress distribution. In comparison to LDBC, the PDBC cases in Figures 10A, 11A, and 12A show less stress magnitudes, due to less strict boundary deformation condition. The RVE can reduce the Cauchy stresses by larger displacements at its surfaces compared to LDBC. The UTBC cases show the highest boundary deformation, compared to LDBC and PDBC, since the surfaces are free to deform and the high deformation yields a high stress reduction.
Additionally, a difference in the spatial direction of the reinforcement can be found for all loading scenarios and all homogenization boundary conditions. Compared to the other two equal directions, obviously the reinforcement leads to higher Cauchy stress magnitudes, due to large fiber stiffness (see Table 2), when deformed in fiber direction.
Comparing Figures 10B, 11B, and 12B, only minor differences of the internal temperature distribution are observed for three loading conditions at the three spatial directions. With respect to the influence of homogenization boundary conditions for one loading direction, it can be found that the thermo-mechanical coupling phenomena lead to significantly different RVE core temperatures, depending on prescribed displacements (LDBC, PDBC) or prescribed tractions (UTBC). Lower stress distribution yields larger core temperature and vice versa. The less restraints on the deformation are prescribed, the more homogeneous the thermo-mechanical response becomes. The numerical material test as well as the parameter identification are carried out for all three types of boundary conditions. In case of LDBC and PDBC, the objective function is composed of Equations (140) up to (144). In case of UTBC, the objective is the minimization of the differences in the volume averaged deformation gradient between the heterogeneous RVE and homogeneous sample computations, due to prescribed forces at the boundaries of the RVE. The objective is to determine the effective bulk modulus , the effective shear modulus equivalent a 1 , the effective thermal expansion coefficients A , B , and A×B (compare Equation (31)), the effective thermal conductivity coefficients k A , k B , and k A×B (compare Equation (30)), as well as the effective heat capacity 0 c, such that the homogeneous sample has an equivalent structural behavior like the heterogeneous RVE.
The orthotropic direction vectors A = [1, 0, 0] T and B = [0, 1, 0] T are given with respect to the fiber orientation and the global coordinates, see Figure 7, and are kept constant for the optimization. The effective shear modulus equivalents a 2 · · ·g 6 are set to be zero and are excluded from optimization.  Table 3 shows the identified set of effective material parameters and the starting values. The starting point of the identification of the effective thermo-elastic material parameters is based on the volume average of the heterogeneous counterparts. As can be seen, each of the different homogenization boundary conditions leads to a different set of model parameters. Recalling Figure 9, the heterogeneous RVE responses for LDBC and PDBC are similar; therefore, a similar set of model parameters is determined by the optimizer. The elastic constants and a 1 have lower values in case of UTBC compared to LDBC and PDBC, due to a much softer response of the heterogeneous RVE in this case. In contrast, the thermal expansion coefficients A , B , and A×B are identified with higher values in case of UTBC compared to the two other boundary conditions. This is caused by the higher deformation of the heterogeneous RVE, shown in Figures 10,  11, and 12. The specific heat capacity 0 c is almost identical and therewith independent of the applied boundary conditions. The thermal conductivity coefficient in fiber direction is the highest in case of UTBC, but identified to be the lowest in the two other directions. Generally, the purely thermal parameters 0 c, k A , k B , and k A×B are almost independent of the applied boundary conditions, since in all three loading cases, the prescribed temperature boundary data is identical. Figure 13 shows the results obtained by the LDBC numerical material tests (exp) at the heterogeneous RVE and the results obtained by the homogeneous sample simulations (sim) with respect to a comparison of the loading cases for the first Piola-Kirchhoff stresses, see Equation (108). Figures 13A, 13B, and 13C refer to loading case 1 and show the comparison between homogeneous simulation and numerical experiment. As can be seen, the results agree qualitatively very  well for the exemplarily depicted loading case 1. Figure 14 shows the results obtained by the PDBC numerical material tests at the heterogeneous RVE and the results obtained by the homogeneous sample simulations with respect to a comparison of the load case 1 for the first Piola-Kirchhoff stresses, see Equation (108). Similarly as the linear boundary constraints, the results agree qualitatively very well for the depicted load case. Figure 15 shows the results obtained by the UTBC numerical material tests at the heterogeneous RVE and the results obtained by the homogeneous sample simulations with respect to a comparison of the first load case for the deformation gradient, see Equation (107). Similarly as the linear and periodic boundary constraints, the identified model parameters in Table 3 agree qualitatively very well for the appropriate load case. The volume averaged external power, see Equation (109) of the heterogeneous RVE and the homogeneous sample according to the prescribed boundary temperature change and the three different constraints, are shown in Figure 16 exemplarily for the second load case. Having a look at Figures 16A and 16B, the comparison between the numerical material test data and the results computed by the set of model parameter in Table 3 is very good. Only in case of UTBC, slight differences are observable in Figure 16C. The reason for the discrepancies is the averaged solution for the overdetermined set of equations in Equation (152). The absolute values for homogenized external power are small and, therefore, underestimated in the solution. The volume averaged heat capacity, see Equation (110), of the heterogeneous RVE and the homogeneous sample according to the prescribed boundary temperature variation at the three different constraints are shown in Figure 17 and exemplarily depicted for load case 3. The effective heat capacity is similarly computed for all different mechanical load cases and for identical thermal loads. Therefore, only small differences for each load case of the homogenization boundary constraints are found. The comparison between numerical experiments and homogeneous identification results is very good for all plots in Figure 17.
The volume averaged heat conduction term, see Equation (113), of the heterogeneous RVE, and the homogeneous sample according to the prescribed boundary temperature change and the three different constraints are shown in Figure 18 for load case 1. The identified effective anisotropic heat conduction coefficients lead to almost equal results for the numerical material test and the homogeneous sample simulation for all plots of Figure 18.
As an additional information for the verification of the parameter identification, Figure 19 shows the convergence behavior of the norm of the error functions. As can be seen, the error norm converges to a steady state after some few iteration steps. Compared to evolution strategies or search algorithms, much less FEM simulations are required and the results obtained are satisfying, as shown above.  Figure 20 shows the Cauchy stress magnitude as well as the internal temperature distribution for the first mixed load case. The stress norm is much more homogeneously distributed for all types of boundary conditions, since only one material is applied to the sample. The stresses due to LDBC are nearly uniformly distributed, while the stresses due to PDBC are periodically and less uniformly arranged over the sample. In case of UTBC, the stresses are low due to larger deformations. The nonuniform distribution is caused by the anisotropy of the constitutive description inside the volume of the sample. Comparing Figure 20A to Figure 10A, the differences of the stress norm distribution is obvious, although the volume average is almost equal, see Figures 13,14,and 15. The temperature distribution in Figure 20B is almost equal to Figure 10B, which leads to almost identical results in Figures 17 and 18. Similar to loading case 1 in Figure 20, the other two mixed loading situations in Figures 21 and 22 are denoted by large differences in the stress norm distribution but show almost equally distributed internal temperature changes. The example demonstrates the applicability of the presented approach to heterogeneous two-scale analysis. The different results for three different boundary conditions need to be considered for a realistic example, since they influence the effective material behavior, as demonstrated above. Therefore, the choice of the appropriate scale-bridging constraint is an important engineering task.

CONCLUSION AND OUTLOOK
The finite thermo-elastic decoupled homogenization approach proposed is based on thermo-dynamically consistent material models. The extension of the anisotropic material model in the work of Kaliske 41 to thermo-elasticity provides the basis to describe a wide range of thermo-elastic characteristics of different RVEs, such as anisotropic thermal expansion effects or heat conduction phenomena.
A consistent theoretical homogenization approach is formulated and introduced. Effective measures for the thermo-elastic properties of the RVE are used and introduced, such as the effective entropy production. Furthermore, the scale transition from microscale to macroscale is carried out by applying an energetic volume averaging. The introduction of the generalized thermo-elastic Hill-Mandel criterion provides the means to derive three types of boundary conditions, in order to determine the thermo-elastic equilibrium of the appropriate RVE. Based on the developments in the work of Terada et al, 37 a decoupling scheme is developed, where effective mechanical and thermal properties of the RVE are computed from numerical material tests by using a gradient based optimization approach. The full set of formulations is processed and algorithmically treated to be implemented into the FEM framework.
The applicability of the presented method to numerical solutions of physical problems is demonstrated for an exemplary fiber reinforced RVE. The wide range of application options for any volume element that denotes a representative microstructure is therewith indicated. A validation using experimental data for the numerical example remains to be carried out. This implies consideration of a physical heterogeneous composite and predictive simulations with identified effective thermo-elastic composite properties.
The extension to thermo-inelasticity of the presented framework and method, in order to capture and describe dissipative phenomena, is a further step toward a much wider range of thermo-mechanical applications.
compare Equations (52) to (54). For sake of completeness, the following expressions are given: