Determining parameters in generalized thermomechanics for metamaterials by means of asymptotic homogenization

Advancement in manufacturing methods enable designing so called metamaterials with a tailor-made microstructure. Microstructure affects materials response within a length-scale, where we model this behavior by using the generalized thermomechanics. Strain gradient theory is employed as a higher-order theory with thermodynamics modeled as a first-order theory. Developing multiphysics models for heterogeneous materials is indeed a challenge and even this ``simplest'' model in generalized thermomechanics causes dozens of parameters to be determined. We develop a computational model by using a given microstructure, modeled as a periodic domain, and numerically calculate all parameters by means of asymptotic homogenization. Finite element method (FEM) is employed with the aid of open-source codes (FEniCS). Some example with symmetric and random distribution of voids in a model problem verifies the method and provides an example at which length-scale we need to consider generalized thermoeleasticity in composite materials.


Introduction
The majority of natural (e.g., polycrystals, wood, and bone) and man-made (e.g., fiber-reinforced composites, concrete, ceramics, and metallic foams) materials are heterogeneous at the micrometer (µm) length-scale. Heterogeneous materials have differing physical properties within the structure at the so-called microstructure; microstructural, crystal structural, or compositional heterogeneity exists. Heterogeneous materials own their widespread use in engineering and scientific applications (e.g., spaceflight technology, energy conversion, or energy storage) to the combination of inherent or tunable mutually beneficial properties such as low relative density, heat insulation, high heat resistance, and chemical resistance, or extreme hardness, [1,2,3,4,5]. Heterogeneous materials, despite being composed of domains possessing distinct physical properties at the microscale, may be modeled accurately at macroscale as homogeneous materials by an effective (homogenized) material-like properties, [6,7,8]. As expected, physical response on the continuum level is strongly coupled to the microscale heterogeneity for a length-scale near the microscale's length-scale.
In the case of thermo-mechanical processes, microscale heterogeneity's role may be significant within a length-scale. For example, due to a mismatch of microscale thermal expansion coefficients, materials subject to high stress or temperature environments (e.g., concrete and bedrock used in nuclear waste storage) exhibit sharp stresses at the macroscopic level [9]. Hence, damage may occur due to induced thermo-mechanical stresses and minimize the operational lifetime of the component. For this reason, efforts have been made to develop more accurate theoretical models to predict material's physical behavior [10,11,12]. The evolution of theoretical models necessitates the development of numerical homogenization methods based on averaging different physical fields to obtain the effective physical properties [13,14]. The average and the calculation of local field quantities are carried out by solving the underlying physics problem within a so-called representative volume element (RVE) to model the microstructure.
In classical Cauchy continuum mechanics, linear elastic models are based on Hooke's law that implies a linear relation between Cauchy stress and strain. In order to encapsulate thermal effects, linear elastic models have been extended to include temperature by adding an extra linear dependency of Cauchy stress, as formulated in Duhamel-Neumann extension in thermomechanics. At a macroscale with several orders of magnitude larger length-scale than the microscale, which is indeed the case in many engineering applications, the aforementioned models are perfectly admissible. However, the same models fail to account for the complex heterogeneous microstructure at the macroscale with a similar length-scale to microscale, accurately [15]. Hence, generalized continuum theories have been developed to counteract the inability of classical continuum mechanics to account for the microstructural effects, [16,17].
The generalized continuum additionally incorporates higher-order gradients of essential kinematic variables and associated length-scale parameters. The most common application of such models is the strain gradient model, where alongside strain, we have a gradient of strain as an additional state variable [16]. The addition of strain gradient introduces higher-order "hyper" stress as a work conjugate of the strain gradient [18]. However, the extension of the strain gradient models to account for temperature poses new challenges [19]. One approach follows Coleman-Noll rational thermodynamics [20] where only temperature is included in Helmholtz free energy in addition to strain and gradient strain [21]. Another approach includes temperature and gradient of temperature into the Helmholtz free energy, [22,10,23]. Both approaches lack experimental data that would provide additional material parameters arising from the inclusion of temperature and gradient of temperature into the Helmholtz free energy. This problem becomes more challenging for models with both temperature and its gradient as we are adding not one but two additional variables in the free energy formulation. Furthermore, another problem that arises from the addition of temperature gradient is an extra time derivative that appears in the flux term. This choice leads to an extension of Fourier's law into Cattaneo's equation, where alongside conductivity, we have an additional parameter coupled with the time derivative of the temperature gradient [24].
Going from the classical continuum to the generalized continuum, homogenization models go from first-order approaches dealing with strain (displacement first derivative) to second-order approaches dealing with strain and gradient of strain (displacement second derivatives) [25,26,27]. Particularly, the first-order (to be precise, first-gradient) approach requires strict separation of scales, and adherence to the concept of local action negates the ability to capture microscale geometry and deal with localization problems [28]. The second-order approach, by virtue of the generalized continuum, enables us to capture the microscale geometry by introducing length-scale into the material constitutive law [29,30,31,13,32]. Although, applications of the first-order homogenization approaches in thermoelastic problems are abundant in the literature, [33,34,35], the second-order homogenization approaches are relatively rare due to previously mentioned problems. Instead of homogenization, multiscale approaches exist, where the finite element method at both scales (FE 2 ) may be used to do a thermo-mechanical analysis of heterogeneous solids [36]. Such an approach requires high-order continuity of macroscale equations, which relies on finite element formulation that should have (at least) C 1 regularity in displacement and temperature fields. A significant characteristic of the multiscale asymptotic approach is the ability to avoid continuity requirements owing to the reestablishment of the high-order macroscale derivatives by post-processing. Thus, other researchers [37,38] used a second-order asymptotic expansion approach to analyze the coupled thermo-mechanical problems. High-order asymptotic models effectively investigate coupled problems by solving periodic functions at the microscale and generating the macroscale displacement and temperature fields [39]. Additional parameters emerge and they need to be explicitly calculated. This work aims to explicitly calculate all the higher-order material terms associated with the generalized continuum model, such as higher-order elastic constants, coupling constants, and parameters associated with temperature.
In the present study, we only include temperature in our model to avoid extending Fourier's law. In this manner, we analyze the simplest thermo-mechanical model in strain gradient elasticity. We follow the asymptotic homogenization in strain gradient elasticity as introduced in [40], verified in [41], and applied in [42]. In order to incorporate temperature in the asymptotic homogenization model, we follow existing methods [43,44,45]. In doing so, we develop a higher-order asymptotic homogenization model for thermoelastic strain gradient materials that accounts for all of the accompanying higher-order material parameters.
The rest of the paper is organized as follows. The higher-order asymptotic homogenization method and computational implementation are explained in detail in the second section. Numerical results and a discussion of higher-order parameters are presented in the third section, followed by the conclusion.

Methodology
We follow the asymptotic homogenization method [46,40] and extend it to thermomechanics. The microstructure is denoted by y and in the rest of the paper we will call it microscale; and its corresponding homogenized continuum is denoted by X, called macroscale. Their transformation is handled by a so-called homothetic ratio, . Thus, we circumvent a scale separation which enables us to use the same coordinate system for both length-scales. The approach is based on the "known" microscale leading to the "sought after" parameters at the macroscale.
We begin with thermomechanics at microscale and use balance of momentum for calculating the displacement, u, by a defined stress, σ, under a given (specific) body force, g, as follows: where we use a comma notation denoting the space derivative and ρ m is the (known) microscale mass density. Herein and henceforth, we use standard continuum mechanics formulation with summation convention over repeated indices. Similarly, the balance of internal energy reads where the specific (per mass) internal energy, u, and heat flux, q, need to be defined. Supply term, r, is the specified internal thermal source. It should be noted that strain used in Eq. (2) is defined in linear form shown below, Geometric nonlinearities are ignored such that the reference frame is equal as the current frame. Therefore, the rate is simply the partial time derivative in the reference frame that we choose as the known initial (undeformed) configuration. Generalization to higher order is adequate by using an energy formulation. By choosing the specific Helmholtz free energy: where T m is the microscale temperature, and η m is the microscale specific entropy. Here we introduce the first simplification, f m = f m (T m , ε m ), indicating that the free energy depends only on temperature and strain. This approach is valid in thermoelasticiy and we circumvent ourselves from justifications like objectivity (usually done in rational thermodynamics) and use a more direct approach of defining the free energy in an axiomatic manner (as in continuum thermodynamics or in non-equilibrium thermodynamics) where the internal energy is simply defined. By insertingu m = (f m + T m η m ) • into the Eq. (2), dividing by T m , using and since there is no dissipative stress in the system (elasticity), we obtain After rewriting the heat flux in a straight-forward manner, we obtain the balance of entropy: The right-hand side is the entropy production in thermoelasticity which is positive according to the second law of thermodynamics. This assertion results in a restriction for the heat flux, herein, we use a linear relation called Fourier's law: where κ m ij is the thermal conductivity. Furthermore, since f m = f m (T m , ε m ), we have η m = η m (T m , ε m ) as a simple mathematical fact based on Eq. (5)-often it is introduced as a principle of equipresence, but there is no need for such a principle, since there is a mathematical justification for this, we refer the readers to [47] for further details. By summing up the equations for thermoelasticity, we obtain In this manner, the whole formulation is reduced to one scalar function, Helmholtz free energy and its definition. Corresponding to the linear material model (Fourier's law) used in heat flux, we continue to model the microscale as a linear thermoelastic material. Thus, we use linear elastic model with the known stiffness tensor, C m , thermoelastic interaction, β m ij = C m ijkl α m kl , with the well-established coefficient of thermal expansion, α m . In this setting, the Helmholtz free energy is modeled as a quadratic one, To compare microscale and macroscale Helmholtz free energies, we simplify Eq. (10) by expanding logarithmic temperature function through Taylor expansion as: where ξ = T /T ref . Thus, first expression on the right hand side of Eq. (10) is expanded, where specific heat capacity is relatied to parameter a m as: If we assume that room temperature, T ref , is at 300 K, we observe that the expansion is accurate for a temperature range from 180 K to 540 K, see Appendix A. We may now introduce above Eq. (12) into Eq. (10) and obtain Furthermore, the symmetry of the strain tensor leads to minor symmetries of the stiffness matrix C m ijkl = C m jikl = C m ijlk , and without loss of generality the symmetry of thermoelastic interaction β M ij = β M ji , we obtain The equations are closed such that thermoelastic material is modeled at microscale by means of Eq. (5), as follows: Furthermore, we consider steady state condition for temperature and displacement by setting their rate terms equal to zero, For the homogenized continuum, we employ one axiom that the free energy within the RVE, Ω, is identical at micro-and macroscale as: Furthermore, we simplify the analysis by assuming a mass ratio for the following terms:

Macroscale
A higher-order macroscale model may be defined by strain ε M , gradient of strain ∇ε M , and temperature T M . In other words, we begin with a specific free energy, where we use the comma notation as the partial space derivative. We emphasize that the microstructure causes higher order in displacement because of homogenization of the structure [48]; however, we exclude temperature gradient from the free energy. We stress that temperature gradient is used in heat flux as a consequence of the second law of thermodynamics. Free energy is obtained from the internal energy, in our formulation, internal energy incorporates reversible quantities.
, is then specified as a quadratic form, where we have used the symmetry of strain, ε M ij = (u M i,j + u M j,i )/2, allowing us to consider additional minor symmetries with the usual restrictions for positive definiteness [49,50], and without loss of generality the symmetry of γ M ijk = γ M kji . In analogy with Eq. (9), the governing equations at the macroscale read by following a variational formulation [51], as follows: Now, by using the model in Eq. (20) for the free energy, in the case of steady state as in Eq. (17), we obtain The main aim is to find a relation between microscale and macroscale parameters. In other words, we start with given parameters in Eq. (17) and obtain the parameters in Eq. (22).
We introduce a geometric center of the RVE, c X, as follows: Assuming displacement and temperature field u M and T M are continuous over the microscale, we approximate the macroscale displacement and temperature by a Taylor expansion around the value at the geometric center and truncate terms with orders higher than quadratic for displacement, since we assume that the energy depends up to the second gradient [52], and truncate terms with orders higher than linear for temperature, since the energy depends on temperature, but not on its gradient. Macroscopic displacement field and its displacement gradients read since (·)| c X is evaluated at the geometric center and thus a constant vanishing by taking its derivative. We stress that there is no scale separation such that the gradient at macroscale is used in this expansion by means of the comma notatopn. Macroscopic temperature field is assumed constant over the RVE (at microscale) leading to In Eq. (24), the first and second derivatives of macroscopic deformation field are the unknowns. They are obtained by spatial averaging and exploiting the fact that terms evaluated at c X are constant within Ω, thus, Going back to the Eq. (24), we replace the displacement gradients with spatially averaged values from Eq. (26), as follows: We use the axiom in Eq. (18) and insert the above averages into the macroscale energy definition on the right-hand side of Eq. (20). All the spatial averaged terms are constant within the RVE such that they are taken out of the integral. For the sake of clarity, we write each term of the free energy at macroscale by denoting the corresponding parameter to be determined where Separating above equations and combining parameters with identical combinations of spatial averages, we obtain macroscale deformation energy,

Microscale free energy
The basic idea behind this method is to transfer between microscale, y, and macroscale, X, by means of the so-called homothetic ratio . We visualize the meaning of this rather abstract constant in Fig. 1 and emphasize that it is a constant within the RVE. In the computation, there is one single coordinate system and we may set = 1 and model the RVE in real geometric dimensions. In general, the homothetic ratio is the connection between = l L = microscale length macroscale length Homothetic ratio such that we obtain y i,j = δ ij / . Microscale displacement field for the RVE is then expanded with regard to , u m (X) = where n u(X, y) is y-periodic as an assertion from the unit cell. Furthermore, we assume that temperature is constant over the RVE, in order to ensure the y-periodicity of a scalar field. By using the chain rule, we obtain the first derivative of (microscale) displacement field, We utilize Eq. (17), use chain rule, and insert Eq. (34), Comparing coefficients in Eq. (35) of the same order of leads to Only possible solution for Eq. (36) is to define 0 u as a function of X because C m and β m depends on local variable y. This argumentation leads to the conclusion, We use a separation of variables or also called Bernoulli ansatz and rewrite, where we introduce unknown tensors ϕ, ψ, P of one rank higher so that the formulation is general. We emphasize that the temperature is only expanded upto one order less than the displacement. By inserting these into Eq. (37), we acquire The only possible general solution is to fulfill (solve) independently the following governing equations in order to obtain ϕ and for acquiring P . Repeating the same procedure for Eq. (38), we have which is rewritten Furthermore, from Eq. (22), by inserting and using the same cut-off procedure for displacement third derivative and temperature second derivative, we obtain u M a,bc δ cj C m ijkl δ ak δ bl + We separate the independent parts and enforce to fulfill in the case of T M ,a = 0 and in the case of u M a,bc = 0. At the moment that we have assumed temperature without higher order contributions, the latter is identically fulfilled. Therefore, we solve Eqs. (44), (45), (50) in order to determine ϕ, P , ψ, respectively, under the condition T M ,a = 0. By introducing these shorthand notations L abij =δ ia δ jb + ∂ϕ abi ∂y j , N abcij =ϕ abi δ jc + ∂ψ abci ∂y j , We need to fulfill By using Eq. (42) in Eq. (32), we obtain We aim for making free energies equivalent, where the macroscale energy is given in Eq. (30), now we want to find an expression for the microscale energy in order to determine the parameters, where T M ,a = 0. Moreover, third derivative in displacement vanishes as before and we obtain after inserting Eq. (27), we acquire By introducing Eq. (57) is rewritten by using Eq. (52), as follows: Using the above equation microscale energy becomes Comparing microscale energy in Eq. 60 to macroscale energy in Eq. 30, we obtain homogenized values

Heat conduction
For the sake of completeness, we additionally calculate homogenized value of thermal conductivity, κ.
We begin with Eq. (22) and use Fourier's law, q i = κ ij T ,j leading to We follow the same procedure and expand the microscale temperature field for the RVE with the same accuracy up to first order in as follows: Expanded microscale temperature field (64) As we are only interested in expansion up to 1 T , thermal conductivity is in relation up to the first-order. In other words, we start with Fourier's equation at the microscale and result in Fourier's equation at the macroscale, Hence, we obtain Substituting local coordinate y into Eq. (64), and using the chain rule, we obtain the first derivative of microscale temperature field, Inserting Eq. (67) into Eq. (63) and again using chain rule, we obtain an asymptotically expanded governing equation: Comparing coefficients in Eq. (68) of the same order of leads to Only possible solution for Eq. (69) is to define 0 T as a function of X because κ m ij depends on local variable y. This observation leads to a straightforward conclusion, By using Eq. (72) in Eq. (70), we obtain where we have utilized Eq. (66). Since the second gradient in T M has been vanishing, we obtain

Numerical implementation in FEniCS
Calculation of macroscale parameters in Eq. (62) and Eq. (75) requires the solution of P i , ϕ abi , ψ abci , and R j tensors from Eq. (45), Eq. (50), Eq. (44), and Eq. (73). The computational work has two steps. First, an RVE, Ω, is created in Salome, which has periodic boundary conditions. The solutions of y-periodic fields, P i , ϕ abi , ψ abci , and R j , have to be periodic. Corresponding surfaces must have the identical mesh such that the solution is restricted by this periodicity. In order to attain this condition adequately, projection method is used in Salome for ensuring that the node positions on corresponding boundaries are matching. Second, the weak form is implemented in a Python code to be solved by the open-source packages developed by the FEniCS project. The weak form is obtained by the standard variational formulation by multiplying the governing equation by an arbitrary test function and integrating by parts in order to reduce the regularity condition of the discrete functions. Discretization for the finite element method (FEM) is established by polynomial form functions. We use the same form functions for the fields and their test functions as known as the Galerkin procedure. For representing a vector, for example, P i in 2D i = 1, 2, we use the Hilbertian Sobolev space, H n , is of polynomial order, n, hence, we use standard (continuous) Lagrange elements in the FEM [53]. As known as the Galerkin approach, we use the same type of a functional space for test functions, where we skip testing the solution at Dirichlet boundaries, Ω D , with the known solution. The computational domain, Ω, is the image of the RVE with the Dirichlet type boundary conditions, Ω D , being basically periodic boundaries for all fields.
This calculation is done by solving the corresponding weak forms for governing equations in Eq. (45), Eq. (50), Eq. (44), and Eq. (73), respectively, Solutions of these fields are then used to construct the macroscale parameters as shown in Fig. 2.

Problem Description
Two types of problems have been designed to study the behavior of homogenized material parameters: • The first problem is intended as a homogeneous material in order to verify that higher-order parameters vanish, G M = 0, D M = 0, and γ M = 0, see Figure 3a.
• The second problem is defined as a porous material with three types of circular pore distributions, Figure 3. All the structures have the same porosity of 20%. The idea behind the second problem has been to use different pore distributions to capture previously examined behavior of the higher-order parameters, [42].
The material used in all cases is aluminum, for which we assume a linear elastic material model at the microscale, with material properties compiled in Table 1. For parameter a M which is related to heat capacity as in Eq. (13), we define for all of the cases T M = 400 K. For a better representation of parameters, we use Voigt's notation. In the case of stiffness tensor, the matrix notation reads as well as

Porous material case
We emphasize that the higher-gradient is often neglected in composite materials. During the standard homogenization by volume averaging, thermoelastic parameters are obtained in the similar sense as the results herein. Yet the difference, herein, relies on obtaining higher-order parameters as well, where their significance depends on the chosen length-scale, see for a numerical study [54].

Single central pore case
In the case of a single centrally located pore, we observe cubic material behavior by inspecting the stiffness tensors components, C M 1111 = C M 2222 . In determining mechanical parameters, we assume constant temperature, T M , over the whole RVE. Thus, thermoelastic interaction, β M , is solely a function of geometry. In other words, it affects only the volumetric component of thermal strain. Higher order parameters G M and γ M are essentially zero due to the centro-symmetry of the RVE with one centrally located pore. Yet higher-order parameters, D M , arise as follows:

Four uniformly distributed pores case
In case of four uniformly distributed pores, we observe the same stiffness tensor, since the porosity is kept the same. We emphasize that the stiffness tensor depends on the porosity; however, the higherorder parameters do depend on the homothetic ratio . The results for all of the parameters are as follows:  In case of four randomly distributed pores, the microscale structure creates an anisotropic material behavior at the macroscale. Thermoelastic interaction, β M , is again a function of geometry so that in the random case, it affects not only the volumetric component of thermal strain but also the shear component of thermal strain. Higher-order parameters G M and γ M are different from zero as the random distribution circumvents the centro-symmetry of the RVE, while D M is also affected by the random distribution of the pores and shows indeed an anisotropic behavior as well. The results for all of the parameters are as follows:

Conclusion
A higher-order asymptotic homogenization model for generalized thermomechanics has been proposed and implemented by means of strain gradient elasticity and a first-order thermodynamics modeling approach. This model incorporates the effects of microscale morphology through additional (higherorder) thermal and mechanical material parameters at the macroscale. On the mechanical side, we account for stiffness matrix C M and higher-order parameters D M and G M , while on the thermal side, we account for thermoelastic interaction β M and higher-order parameter γ M . All of the thermal and mechanical macroscale parameters are explicitly computed in this work by assuming a linear thermoelastic material model at the microscale.
In this framework, the higher-order asymptotic homogenization is implemented in the FEniCS platform and used to solve the partial differential equations generated from the homogenization procedure. The methodology is verified by using a problem without the microstructure leading to homogeneous parameters. Also a porous material is calculated providing insight how the parameters alter for three types of distributions (single, uniform, random distribution). These cases has shown that thermoelastic interaction β M has a similar reaction to pore morphology as the stiffness matrix C M while higher-order thermal parameters γ M mirror the response of the higher-order mechanical parameter G M as both of them are linearly dependent to the size of the RVE and have the same behavior with respect to the centro-symmetry of the RVE. Even though these numerical results for thermal parameters still need indepth numerical analysis and experimental validation, we provide herein a complete methodology and its implementation to encourage further research for a better understanding of the interplay between microscale morphology and thermo-mechanical material parameters.

Acknowledgment
This work was supported by a project entitled "Time-dependent THMC properties and microstructural evolution of damaged rocks in excavation damage zone" funded by the U.S. Department of Energy (DOE), Office of Nuclear Energy under award #DE-NE0008771.

Appendix A Taylor Expansion of the Logarithm Function
First parameter on the left-hand side of Eq. (10) needs to be simplified in order to compare microscale and macroscale Helmholtz free energies. This simplification is done by expanding the logarithmic term through Taylor expansion, which may have several different forms depending on the value of ξ, see Eq. (11). What follows is fully developed Eq. 12 for ξ ≥ 1 2 , as follows: As we are specifying ξ ≥ 1 2 , we need to determine the temperature range where Eq. 12 is valid. In Figure 4 we show comparisons between ln T T ref and its Taylor expansion. If we assume that room temperature, T ref , is at 300 K, from Figure 4, we see that the expansion is accurate for a temperature range from 180 K to 540 K.