Computing the effective response of heterogeneous materials with thermomechanically coupled constituents by an implicit fast Fourier transform‐based approach

Thermomechanical couplings are present in many materials and should therefore be considered in multiscale approaches. Specific cases of thermomechanical behavior are the isothermal and the adiabatic regime, in which the behavior of real materials differs. Based on the consistent asymptotic homogenization framework for thermomechanically coupled generalized standard materials, the present work is devoted to computing the effective thermomechanical behavior of composite materials in the context of fast Fourier transform (FFT)‐based micromechanics. Exploiting the homogeneity of the temperature on the microscale, we develop a fast implicit staggered solution scheme for the coupled problem, which is compatible to existing strain‐based micromechanics solvers. Due to its implicit formulation, the algorithm permits large time steps for computations involving strong thermomechanical coupling. We investigate the performance of modern FFT‐based algorithms combined with the proposed thermomechanical solution strategy. In this context, the Barzilai–Borwein method is identified as particularly efficient, inducing only a small overhead compared with the traditional isothermal setting. We demonstrate the effectiveness of the presented approach for short‐fiber reinforced composites with viscoelastic matrix behavior.

equation). For instance, in the vicinity of their glass transition temperature polymers are particularly sensitive to temperature variations. 1 Especially when subjected to cyclic loading, self-heating due to dissipation can critically affect the mechanical properties of materials and the life time of components, see, for example, Rittel, 2 Mortazavian-Fatemi, 3 or Katunin. 4 Thus, for the optimal use of materials, characterizing and predicting their thermomechanical behavior is of central importance. For composite materials, this proves to be a challenging task as their properties rest on their individual constituents and microstructure. In a small-strain framework, Chatzigeorgiou et al. 5 used an asymptotic homogenization approach to derive the governing thermomechanical equations on the micro-and macroscale for generalized standard materials, 6 taking into account both the microstructure and the thermomechanical material behavior. This generalized previous studies using asymptotic approaches, for example, by Terada et al. 7 for porothermoelasticity or Temizer 8 for finite thermoelasticity. A recent review on the homogenization of dissipative materials was given by Charalambakis et al. 9 A particular result of the asymptotic homogenization is that the microscopic balance of linear momentum depends only on the macroscopic temperature and is independent of temperature fluctuations on the microscale. 5 In contrast to earlier works on the homogenization of thermomechanical material properties, for example, by Willis, 10 the uniform temperature on the microscale is not an ad hoc assumption but arises as a direct consequence of first-order homogenization. As a result, the thermomechanical problem on the microscale may be solved for a homogeneous temperature and is decoupled from microscopic heat-conduction. Based on these results, Tikkarrouchine et al. 11 homogenized unidirectional short-fiber structures with temperature-independent material parameters in the context of concurrent multiscale simulations, using the finite element (FE)-software ABAQUS. Similar FE-based multiscale studies which still consider thermal conduction on the microscale were carried out by Özdemir et al. 12 for elastoplasticity and Li et al. 13 for single-crystal elastoviscoplasticity.
Motivated by the aforementioned studies, we consider solvers based on the fast Fourier transform (FFT) for the computational homogenization of thermomechanically coupled materials on the microscale. Since the pioneering work of Moulinec-Suquet, 14 these algorithms have become a well-established tool in computational micromechanics, owing to three major advantages. First, FFT-based solvers operate on a regular grid, enabling them to handle complex microstructures and making them directly compatible to modern volumetric imaging techniques such as microcomputed X-ray tomography. Second, in contrast to traditional FE-based solvers, the FFT-based methods come with a natural matrix-free implementation. This is especially relevant for finely discretized complex microstructures, where memory can become a limiting factor. Last but not least, FFT-based methods have proven to be computationally efficient, typically outperforming traditional FE-based solvers. [15][16][17] Indeed, FFT-based methods profit from optimized and freely available implementations of the eponymous FFT. 18 Furthermore, research effort was devoted to developing fast solution algorithms, which improve upon the performance of Moulinec-Suquet's basic iterative scheme. Currently, the most efficient families of solution schemes are the polarization-based methods, [19][20][21] inexact Newton methods 22,23 and quasi-Newton methods. [24][25][26] An overview of the advantages of these methods is provided in Appendix A. Comparative studies on the performance of various FFT-based solvers were carried out, for instance, by Mishra et al. 27 or Schneider. 25 FFT-based methods have found widespread application in various fields, such as fracture mechanics, 28,29 interface damage, 30 and gradient plasticity, 31 as well as coupled problems such as electromechanics 32 and recrystallization kinetics. 33 In the context of thermomechanics, FFT-based methods have been used to homogenize linear thermoelastic materials [34][35][36] and linear thermomagnetoelectroelastic 37 materials. Shantraj et al. 38 proposed a FFT-based staggered algorithm for coupled multiphysics problems, taking thermal conduction on the microscale into account.
To exploit the power of FFT-based methods for computing the effective thermomechanical behavior of nonlinear dissipative materials, we rely upon the framework of asymptotic homogenization, as pioneered by Chatzigeorgiou et al. 5 Due to the weak coupling of mechanics and thermal conduction, the cell problem on the microscale is governed only by the microscopic balance of linear momentum and the evolution of the macroscopic temperature, cf. Section 2. Based on these results, we propose a staggered solution algorithm, where strain-field and temperature are updated in an alternating fashion, cf. Section 3. The proposed solution scheme may be applied on top of any iterative strain-based solution method and can be easily integrated into existing FFT-based computational micromechanics codes. Owing to the homogeneity of the temperature on the microscale, the temperature update only involves solving a scalar equation and introduces little overhead. The usefulness of the approach is demonstrated in Section 4 for glass-fiber reinforced polypropylene composites with strong thermomechanical coupling.

FIRST-ORDER HOMOGENIZATION OF THERMOMECHANICAL COMPOSITES
Chatzigeorgiou et al. 5 introduced a framework for the asymptotic homogenization of thermomechanically coupled generalized standard materials in the quasi-static small-strain setting. As a result, they obtained governing equations for macro-and microscale. In the following, we review the equations relevant for solving the thermomechanical cell problem on the microscopic level.
Let Y ⊆ R d be a rectangular cell, with microscopic point x ∈ Y and d ∈ {1, 2, 3} spatial dimensions. We denote by Sym(d) the space of symmetric d × d matrices. For the following discussion, we consider the displacement fluctuation with a sufficiently large vector space Z and the macroscopic absolute temperature ∈ R >0 . For a heterogeneous Helmholtz free energy density which is related to the internal energy e by we express stress and entropy by the potential relations under the assumption that is differentiable in all arguments except for the first. 39 As a result from the asymptotic homogenization of Chatzigeorgiou et al., 5 only the macroscopic temperature enters . Thus, the temperature in a microstructure cell corresponding to a macroscopic point can be interpreted as homogeneous.
We assume that the free Helmholtz energy density can be additively decomposed into a component heat associated to heat storage and a component mech representing the storage of mechanical energy. This splitting does not reflect physics, but is computationally convenient, cf. Section 3. Many commonly used thermomechanical material models, such as viscoelasticity, 11 elastoplasticity, 5 and viscoplasticity 40 feature a free energy in the form of (3). The heat capacity density at constant strain is typically assumed to be independent of strain and internal state z. Under this condition, the temperature dependence of the mechanical free energy mech may at most be linear. Consequently, we also partition the entropy For generalized standard materials, the evolution of internal variables is governed by Biot's equation We assume that is convex in its third argument and (⋅, , 0) = 0 holds.
For the stress and strain field, the microscopic static balance of linear momentum without volume-force densities div = 0, (6) and the kinematic compatibility condition hold, where ⟨⋅⟩ Y = 1∕|Y |∫ Y (⋅) dV denotes the volume average over Y and ∇ s stands for the symmetrized gradient. The macroscopic temperature is determined by the macroscopic balance of internal energẏ where we neglect additional source terms and div x denotes the divergence with respect to the position x ∈ Ω in the macroscopic body Ω ⊆ R d . It is common to reformulate the balance of internal energy as a heat equation in terms of the entropẏ or the temperature Note that, in the small-strain setting, the material time derivative( ⋅) reduces to the local time derivative (⋅) t . As we are only interested in solving the cell problem, that is, we only consider a single macroscopic point, the term −div x ⟨q⟩ Y cannot be further specified and acts as a volumetric heat supply term. Hence, we denote  = −div x ⟨q⟩ Y and treat  and as boundary conditions. For a treatment in a concurrent multiscale context, cf., for example, Chatzigeorgiou et al. 5 or Tikkarrouchine et al. 11

SOLUTION SCHEME FOR THE FULLY COUPLED THERMOMECHANICAL CELL PROBLEM
Consider the Hilbert space V = L 2 (Y ; Sym(d)) of Y -periodic and square integrable stress and strain fields with inner product and the induced norm For a certain point in time, we want to find a strain field and a macroscopic temperature which solve Equations (5)- (8) for prescribed and . For the convenience of the reader, we restrict to pure strain boundary conditions, see Kabel et al. 41 for an extension to mixed boundary conditions. To solve our problem, we consider a fixed time step and apply an implicit Euler discretization in time to our system of equations. We define the operator with the mean entropy s n and internal variables z n at the last converged time step and the time increment Δt. When evaluating M( , ) or H( , ), the internal variables z are computed by solving the discretized Biot's equation for given strain-field and temperature . The thermomechanical cell problem is defined by the system of equations where (14) describes the mechanical problem for the strain-field and (15) is the thermal problem, determining the evolution of the temperature . There exist two general approaches for solving the thermomechanically coupled problem. In monolithic schemes, (14) and (15) are solved simultaneously, whereas staggered approaches treat the subproblems (14) and (15) separately. 42,43 Monolithic approaches enjoy unconditional stability, but the resulting system is usually nonsymmetric. 42 Provided each subproblem by itself is symmetric, staggered schemes circumvent this difficulty and thereby enable using more efficient solution algorithms. 44,45 Furthermore, they are convenient in terms of implementation, as existing solvers for the subproblems may be used. 38,46,47 Hence, we focus on staggered algorithms in the following.
Typically, staggered schemes are based on an isothermal split, 42,44 where the mechanical problem is solved for a fixed temperature and the thermal problem is solved for a fixed strain-field. More precisely, for given iterates k and k , where 0 = n is set to the temperature in the last converged time step, the following steps are performed: 1. Solve M( , k ) = 0, with fixed temperature k and assign the solution to k + 1 . 2. Solve H( k+1 , ) = 0, with fixed strain-field k + 1 and assign the solution to k+1 .
In this context, we distinguish between explicit and implicit staggered schemes. For explicit schemes, steps 1 and 2 are carried out only once, whereas for implicit schemes the steps are repeated until a prescribed convergence criterion is fulfilled.
Thus, explicit schemes are naturally faster. However, they suffer from lower accuracy 47,48 and are prone to instabilities 42,46,49 for problems with strong thermomechanical coupling. To address the latter difficulty, Armero-Simo 42,49 proposed an unconditionally stable adiabatic split, where (14) is solved under the conditioṅs = 0. In the present work, we do not follow this approach (see below for a discussion) and consider an implicit staggered approach with an isothermal split. Implicit staggered schemes enjoy the same accuracy as monolithic algorithms 43 and have been shown to be more stable than explicit schemes. 46 However, when repeating steps 1 and 2 until convergence, the subproblems (14) and (15) have to be solved multiple times per time step, which is computationally expensive. Therefore, we propose two simplifications to enhance the overall efficiency of the scheme.
First, suppose we have an iterative strain-based fixed point scheme which solves (14) for a fixed temperature . For better readability, we suppress the possible dependency of f on additional algorithmic parameters and the boundary conditions . Instead of solving (14) after each temperature update, we only perform a single iteration (16) of the mechanical solver. The second simplification concerns the temperature update, that is, solving (15). Evaluating H( , ) is computationally expensive, as it involves solving (13) for all points in Y to compute the mechanical entropy s mech (⋅, k , k , z k ) and the dissipation z (⋅, k , k , z k ) ⋅ (z k − z n ), cf. (12). To obtain an efficient algorithm, we wish to avoid this operation outside of (16), that is, without improving our current guess for the strain field. Thus, we propose an additive split of H( , ) into an implicit part H impl ( ) and an explicit part H expl ( , ) following our partition of the entropy. We emphasize that this splitting is not physical but computationally convenient. Instead of solving H( k+1 , ) = 0 for updating the temperature, we solve H impl ( ) + H expl ( k , k ) = 0. More precisely, we compute the effective mechanical entropy as well as the mean dissipation as part of our mechanical iteration (16), cf. Miehe. 50 Subsequently, we solve This is significantly more efficient than solving H( k+1 , ) = 0, as it only involves the effective entropy related to heat storage, which is efficiently computed by for an N-phase composite material with volume fractions c j and phase-specific entropies s heat, j . To summarize, our modified implicit algorithm involves the following steps, which are repeated until the convergence criterion of the mechanical solver is met: 1. Update the strain field k+1 = F( k , k ) with a single iteration of the mechanical solver (16). Compute s mech,k and  k as part of the iteration. 2. Solve H split ( ) = 0 and assign the solution to k+1 .
The proposed algorithm is compatible to any mechanical solver in the form of (16), including classical FE-based methods.
For our concrete implementation, we rely on FFT-based solution schemes, due to their computational efficiency. [15][16][17] In particular, we consider Moulinec-Suquet's basic scheme, 14 the Barzilai-Borwein method, 25 and the inexact Newton-CG method, 23 cf. Appendix A. Typically, the iteration scheme (16) involves applying the operator Γ = ∇ s (div ∇ s ) −1 div in Fourier space and evaluating material law = (⋅, , , z). As convergence criterion for the static equilibrium (14), we use see section 5 of Schneider et al. 51 for further details. The mean mechanical entropy s mech,k and the mean dissipation  k are computed when evaluating the material law. For the temperature update we use Newton's method, that is, we iterate is met. Thus, we set k+1 = i . If the convergence criterion (21) is met, we proceed to the next time step. Otherwise, we repeat updates (16) and (19). The algorithm is summarized in Algorithm 1. Algorithm 1. Implicit staggered solution scheme ( , maxit mech , mech , , maxit heat , heat ) 1: Set initial values for̄and 2: k ← 0 3: r mech ← 1 4: while k < maxit mech and r mech > mech do 5: k ← k + 1 6: (16) 7: while i < maxit heat and r heat > heat do ⊳ Temperature update (19) 10: 1. For the temperature update, we use the entropy-based heat equation (9) instead of the more common temperature-based formulation (10). Consider the change in s mech under the assumption that the heat capacity c , as defined in (4), depends only on the temperature. Using the implicit Euler time discretization oṅs heat in (9) yields If, alternatively, we discretized the corresponding term c ( )̇in (10), we obtain Apparently, the change in entropy is basically linearized. Hence, to obtain higher precision for large time increments we prefer using (9). 2. For the temperature update (19), we only consider the temperature dependency of the entropy related to heat storage s heat , whereas s mech and  remain fixed. Indeed, if the heat capacity c depends only on the temperature, s mech is temperature independent and  is at most a linear function of the temperature. Thus, as the strain field converges, changes in subsequent iterates k become small. As a result, the solution of (19) approaches the solution of (15). 3. Due to the homogeneity of the macroscopic temperature , computing the mean entropy (20) is comparatively inexpensive. Thus, solving the scalar equation (19) introduces no significant computational overhead. Compared with solving the isothermal mechanical problem, higher computation times may still arise for the thermomechanically coupled case, due to two factors. First, the iterative solver (16) may require more iterations to converge, depending on the thermomechanical coupling of the composite, that is, the temperature dependence of the material laws and the magnitude of s mech and . Second, evaluating s mech and  may affect the overall runtime of the algorithm, provided the associated computational effort is similar to the evaluation of the material law and the Γ operator. 4. Armero-Simo 42,49 analyzed different operator splits for thermomechanically coupled problems and found that the explicit isothermal split, that is, an isothermal mechanical step followed by a temperature update, is only conditionally stable. As an alternative, they proposed the unconditionally stable adiabatic split, where the mechanical problem is solved under the conditioṅs = 0. For the present algorithm, we rely on the isothermal split as it is more convenient from the viewpoint of implementation. Suppose we already have an existing code for a purely mechanics-based solution scheme. For the isothermal split, only an update of the temperature dependent material parameters and the computation of s mech and  have to be added to the already implemented material law. For the adiabatic split, on the other hand, simple reduced forms of the material law, which identically fulfill s = 0, can only be derived in special cases such as linear thermoelasticity. 42 For more complex material laws with arbitrary temperature dependencies, the implementation of an additional adiabatic formulation may be cumbersome or even require an iterative local solution scheme. Thus, for tackling the issue of instability, we prefer using an implicit staggered approach based on an isothermal split. 46 Indeed, we encountered no numerical instabilities in our numerical experiments in Section 4, even for a composite with strong thermomechanical coupling.

Setup
Algorithm 1 for thermomechanically coupled problems was implemented in an in-house FFT-based computational homogenization code written in Python 3.7 with FFTW 18 bindings. Applying Γ and evaluating the material law were integrated as Cython extensions 52 and parallelized using OpenMP. Throughout, we rely on the discretization by trigonometric polynomials introduced by Moulinec-Suquet. 14 As convergence criterion for the iterative FFT-based solver, we use (21) with a prescribed tolerance of mech = 10 −5 . The tolerance for the convergence criterion (22) of the temperature update For the computations on the two-dimensional microstructure in Section 4.2, a desktop computer with 32 GB RAM and a six-core Intel i7-8700K CPU was used. The computations on the three-dimensional microstructure in Section 4.3 were performed on a workstation with 512 GB RAM and two 12-core Intel Xeon(R) Gold 6146.

Continuous glass-fiber reinforced polypropylene
In the following example, we consider a composite consisting of a polypropylene matrix unidirectionally reinforced by continuous glass fibers with a volume fraction of 30%. The microstructure, cf. Figure 1, is modeled as a two-dimensional periodic cell with a resolution of 512 2 , containing 200 fibers. It was generated using the adaptive shrinking cell algorithm of Torquato-Jiao. 53 The glass fibers are modeled as an isotropic linear thermoelastic material. The free energy related to heat storage reads and corresponds to a material with a constant heat capacity c ( ) = c 0 . Typically, for solids, states of constant strain are hard to realize under fluctuating temperatures. Hence, the heat capacity c at constant strain is usually not measured experimentally. However, its value is typically close to the heat capacity c at constant stress. The mechanical part of the free energy is given by implying the stress-strain relation with a stiffness tensor C and a thermal expansion tensor ∈ Sym(d). The associated entropies read and s mech ( ) = ∶ C ∶ .
As the material is elastic, no energy is dissipated, that is,  = 0, and the thermomechanical coupling is governed solely by s mech . Changes in s mech cause self-heating under hydrostatic compression and self-cooling under hydrostatic extension. This phenomenon is commonly referred to as thermoelastic coupling effect, cf. section 13.2 in Haupt, 54 or Gough-Joule effect, cf. section 96 in Truesdell-Noll. 55 For the glass fibers, we assume that both stiffness and thermal expansion are isotropic, that is, with bulk modulus K, shear modulus G, and isotropic coefficient of thermal expansion 0 . By P 1 and P 2 we denote the projectors onto the spherical and deviatoric d × d matrices, respectively. The parameters of the model are taken from Tikkarrouchine et al. 11 and listed in Table 1.
For the polypropylene matrix, we assume a linear thermoviscoelastic model based on a generalized Maxwell model, cf. Figure 1 and section 3.5.1 in Tschoegl's book. 56 For models accounting for effects outside of the viscoelastic domain, we refer, for example, to Krairi et al. 57 and Benaarbia et al. 58 for extensions to viscoplasticity and damage, and Tscharnauter et al. 59 for a study on polypropylene. Based on the caloric data in table 18.10 in the Springer Handbook of Materials Data, 60 we assume a heat-storage related free energy of the form corresponding to a linear heat capacity The energy stored the generalized Maxwell model with N MW Maxwell elements reads Consequently, the stress computes as We assume that the viscosity tensor associated to a dashpot of the generalized Maxwell model has the form where a ∶ R >0 → R denotes a temperature-dependent shift function. The corresponding fluidity F is defined by the pseudoinverse In terms of the partial stresses the evolution equation for the viscous strains readṡv For simplicity, we assume that polypropylene is isotropic and linear elastic in dilation, cf. section 9.4 in Brinson-Brinson. 61 More precisely, the stiffness tensors and the thermal expansion tensor have the form For the present study, we restrict to linear viscoelastic behavior and focus on the effects induced by the thermomechanical coupling. In particular, we omit a possible pressure dependence of the shift factor as suggested by Fillers-Tschoegl 63 based on free-volume considerations.
For our implementation, we use the time-integration scheme of Taylor et al., 64 which is based on the partial stresses v instead of the viscous strains v . The update reads v = exp where (⋅) n denotes the value of the last converged time step and is a reduced time defined via We compute the change in reduced time Δ = − n by a five-point Gauss quadrature, assuming a constant temperature rate. The entropies and dissipation in terms of the partial stresses v read ] , The used material parameters are listed in Table 2. The caloric parameters were chosen based on tables 18.9 and 18.10 in the Springer Handbook of Materials Data 60 and the viscoelastic parameters are taken from the experimental study by Kehrer et al. 65 Note that Kehrer et al. 65 characterized the behavior of polypropylene over a wide range of frequencies and temperatures, using 27 Maxwell elements for their model. For the present study, we restrict to moderate temperature and frequency changes and only consider nine elements with time constants ∈ [10 −4 , 10 4 ] in order to reduce computation times. The shear moduli of the elements with > 10 4 are added to the elastic shear modulus G 0 , whereas the elements with < 10 −4 were omitted.

Uniaxial extension
In our first set of experiments, we take a look at the stress-strain behavior under uniaxial extension and compression. We want to assess the strength of the thermomechanical coupling for the investigated composite microstructure. In addition, we are interested in the performance of the different FFT-based solution algorithms in Appendix A in conjunction with the staggered thermomechanical solution scheme in Algorithm 1. To this end, we apply mixed boundary conditions, cf. Kabel et al., 41 corresponding to strain-controlled uniaxial extension/compression to 5% with a strain rate of 1/s at various loading angles in the xz-plane with respect to the x-direction, that is, the fiber direction. For the first set of computations, we consider adiabatic conditions, that is,  = 0, where self-heating/-cooling of the material is expected. The second set of computations is performed with a fixed temperature of ref = 293.15 K as reference.
The resulting stress-strain curves are plotted in Figure 2. In the isothermal setting, there is no distinction between tension and compression and we observe a linear relation between stresses and strains. For the loading under adiabatic conditions, however, the thermomechanical coupling induces an effectively nonlinear behavior. To be more precise, under compression, s mech decreases, which leads to a rise in temperature, resulting in the softening of the polypropylene matrix. Conversely, under tension s mech increases, leading to a lower temperature and the stiffening of polypropylene. Two factors contribute to the strength of the observed thermomechanical coupling. First, due to its high thermal expansion coefficient , the Gough-Joule effect, that is, the strain-induced change of s mech is rather pronounced for polypropylene. Second, the mechanical behavior of polypropylene is very sensitive to temperature changes in the vicinity of its glass transition temperature, as encapsulated by the WLF equation (23). Note that the computations for the 0 • load angle represent an exception to these observations. In this case, the fibers carry most of the load and we observe no difference between isothermal and adiabatic computations, due to temperature independence of their stiffness.

Performance comparison for a single load step
Next, we take a closer look at the performance of the different FFT-based solution schemes. In particular, we are interested how their convergence behavior changes in case of strong thermomechanical coupling, compared with the isothermal setting. Hence, we consider the load case of uniaxial extension at a 90 • load angle, where the coupling is most pronounced. First, the performance is evaluated for a single load step up to 5% strain. The residual is plotted as a function of iteration counts and computation time in Figure 3. Note that the convergence behavior of the Newton-CG method in the adiabatic setting is distinctly different in comparison with the isothermal computation. For the isothermal case, the decrease of the residual gradually grows in subsequent Newton iterations. Due to the adaptive forcing-term choice of Eisenstat-Walker, 66 the linear system is thus solved to higher accuracy. In contrast, the convergence rate with respect to Newton iterations is roughly constant for the adiabatic computation. This is due to the fact that we do not consider the temperature dependence of the material behavior in the computation of the Hessian. Thus, the linear approximation of the gradient is less precise than for the isothermal computation. With respect to the overall performance, this effect is somewhat alleviated by Eisenstat-Walker's forcing-term choice, 66 as the linear system is solved to lower accuracy, thereby reducing the cost of each Newton iteration. Even though, Newton-CG requires 75% more Newton iterations in the adiabatic setting, the runtime only increases by about 30%, cf. Table 3.
For the basic scheme, the convergence rate in the adiabatic and isothermal setting is nearly identical. The same is true for the Barzilai-Borwein method, which displays its characteristic nonmonotone behavior and converges in much fewer iterations than the basic scheme, cf. Table 3. Even though the iteration counts of both schemes are roughly identical for both settings, the overall computation times are slightly higher for the adiabatic computations.
A look at the computational cost of the most expensive operations, that is, the material law, the FFTs and the Γ-operator, clarifies this phenomenon. In Table 4, the average computation times per application of these operations are listed for the 0 • load case solved by the Barzilai-Borwein scheme. Notably, for the adiabatic setting, the additional computation of  and s mech in the material law increases its time per application by about 70%. Thus, the overall cost per iteration ends up 30% higher. The same is true for the basic scheme. Comparing the overall performance of the schemes, we observe that the Barzilai-Borwein method is the fastest for both the isothermal and adiabatic setting. The Newton-CG method is only slightly slower but suffers from increased iteration counts for the adiabatic case. The basic scheme is by far the slowest, taking five to six times longer than the Barzilai-Borwein method.

Performance comparison for 20 load steps
Next, we investigate the performance of the solvers, when subdividing the strain loading of 5% into 20 equally spaced load steps. An affine-linear extrapolation 14 is applied at the beginning of each load step to obtain an initial guess for the strain field. The total iteration counts and computation times of the different solvers in each load step are plotted in Figure 4. For all solvers, the iteration counts decrease up to step 5 as the affine-linear extrapolation takes effect. For the isothermal computations, the iteration counts further decrease after this point, due to the linear stress-strain behavior. In contrast, the iteration counts stagnate or, in case of the basic scheme, even increase for the adiabatic computations. This coincides with the onset of the effectively nonlinear material behavior for uniaxial strains larger than 1%, cf. Figure 2(D). Hence, the affine-linear extrapolation becomes less effective, which leads to higher iteration counts compared with the isothermal computations. Note that for a material which already behaves nonlinearly under isothermal conditions, this difference between adiabatic and isothermal computations is expected to be less pronounced ( Table 5).
As for the loading in a single step, the Barzilai-Borwein method is fastest. Its computation time for the adiabatic case increases by roughly 35%, due to higher iteration counts and the additional cost per iteration. The performance of the Newton-CG method is nearly identical to the Barzilai-Borwein method for the isothermal computation. However, it exhibits a larger decrease in performance for the adiabatic computation, with an increase in computation time by nearly 60%. For the basic scheme, the iteration counts are roughly identical for the isothermal and the adiabatic setting. In the first nine steps, the adiabatic computation converges faster, as a consequence of the stiffening due to self-cooling and the resulting reduction in material contrast. For the subsequent steps, the isothermal computation requires fewer iterations, due to the more effective affine-linear extrapolation. Fortuitously, these effects roughly cancel each other out. Overall, the basic scheme is still the slowest, taking three to four times longer than the Barzilai-Borwein method to converge.
To summarize, we observe that the convergence behavior of the basic scheme and the Barzilai-Borwein method in conjunction with Algorithm 1 is similar to their convergence behavior under isothermal conditions, even for a composite with strong thermomechanical coupling. The computation times for the thermomechanically coupled computations increase by roughly 30% for both schemes, which is mainly due to the additional cost of computing the dissipation  and mechanical entropy s mech in the material law. The Newton-CG method suffered the highest decrease  in performance for the coupled computations, as the temperature dependence is neglected in the Hessian computation. This leads to a significant increase in Newton-and CG-iterations, in addition to the higher cost per Newton iteration.
Considering the overall performance, the Barzilai-Borwein method and the Newton-CG method are the fastest solvers. Due to its lower memory requirements and its more robust convergence behavior in the thermomechanically coupled computations, we use the Barzilai-Borwein method for all following computations.

Planar short glass-fiber reinforced polypropylene
Motivated by the numerical experiments in the last section, we investigate a more complex microstructure, cf. Figure 5.
We consider a polypropylene matrix reinforced by 1130 short glass-fibers with an aspect ratio of 20. The fiber volume fraction amounts to 13.2%, corresponding to mass fraction of 30%. The microstructure was generated by the sequential addition and migration algorithm 67 and discretized by 512 × 512 × 64 voxels. The second-order fiber-orientation tensor reads A = diag(0.45, 0.45, 0.1), cf. Advani-Tucker. 68 For the following investigations, we use the same material models and parameters as in Section 4.2.

Dynamic mechanical analysis
The macroscopic behavior of viscoelastic composites is often investigated under steady-state oscillations with a fixed frequency f ∈ R ≥0 , cf. section 5.5 in  This is commonly called dynamic-mechanical analysis (DMA). Suppose a linear viscoelastic material is harmonically excited by uniaxial tension/compression where the strain component in loading direction is given by with the strain amplitude amp ∈ R ≥0 and the angular frequency = 2 f . The stress response of the material in loading direction reads with the stress amplitude amp ∈ R ≥0 and phase difference ∈ [0, ∕2]. Typical characteristics for the material are the storage modulus and loss modulus The storage modulus is related to the average elastic energy stored in a load cycle cycle = 1 4 2 amp E ′ and serves as a measure of the material's elastic stiffness. The loss modulus is proportional to the energy dissipated over a load cycle cf. section 9.1 in Tschoegl. 56 Thus, E ′′ is of particular interest in cases of harmonic loadings with high cycle counts. For instance in fatigue experiments, 69,70 the dissipated energy accumulates, leading to an increase of temperature over time.
For linear viscoelastic material models, such as the generalized Maxwell model for polypropylene, E ′ and E ′′ can be computed analytically in the isothermal setting, cf. section 11.1 in Tschoegl. 56 However, as we have seen in Section 4.2.1, F I G U R E 6 Polypropylene: relative error between analytic values and the results of the virtual DMA tests for E ′ , E ′′ , and  cycle as a function of load steps per cycle the thermomechanical coupling induces a nonlinear behavior due to self-heating and self-cooling. Thus, we characterize the viscoelastic behavior of the composite by simulating DMA tests. More precisely, we run through the following steps: 1. In analogy to tensile DMA experiments, a static uniaxial tensile load of static is applied in a single step in 1 second. 2. The loading static is held constant for 100 s. In actual experiments, the holding time is usually much shorter. However, we want to mitigate the effects of the initial stress relaxation on our numerical experiments.  (26) and (27).
Note that we do not use the affine-linear extrapolation for these computations, due to the nonmonotone loading. To validate our approach and to determine the necessary number of load steps per cycle, we apply steps 1-4 for a homogeneous polypropylene microstructure under isothermal conditions. The parameters for the sinusoidal loading are static = 0.1%, amp = 0.05%, cf. Kehrer et al., 65 and f = 10 Hz. For this frequency, the storage and loss modulus of the viscoelastic model for polypropylene are given by E ′ = 2012.22 MPa and E ′′ = 177.68 MPa. In addition to E ′ and E ′′ , we also track the effective dissipated energy (18) in our computations and compare it to the analytical formula (28). The relative errors for E ′ , E ′′ and  cycle are shown in Figure 6 as a function of the load steps per cycle.
For more than 30 load steps per cycle, the relative error for all tracked quantities falls below 1%. Indeed, E ′′ as determined by our DMA computation virtually coincides with its analytical value. Note that the error in dissipation does not tend to 0 for finer resolutions. This is a consequence of the stress relaxation under static strain loading, which still causes a small additional amount of energy dissipation. In preliminary computations, a higher number of cycles was considered as well. However, the results did not differ substantially. Hence, we choose 30 load steps per cycle for all subsequent computations.
With the established procedure, we simulate uniaxial DMA tests at various static load values for the planar short glass-fiber reinforced polypropylene microstructure, see Figure 5. In particular, the effect of the thermomechanical coupling under adiabatic conditions on E ′ and E ′′ is of interest. The loading is applied in the xz-plane at angles between 0 • and 90 • with respect to the x-direction. Static loads static between 0.1% and 1.0% are considered. The amplitude and frequency are fixed at amp = 0.05% and f = 10 Hz, respectively. In Figure 7, the results for E ′ and E ′′ are plotted alongside the mean temperature during the harmonic excitation as a function of the loading angle. First, we take a look at the storage modulus. For the 0 • load case, that is, in-plane loading, the storage modulus is at its peak value. This is due to the stiffening effect of the fibers. matrix and an increase E ′ . The effect is most pronounced for the 90 • load case, where we observe the largest temperature difference between adiabatic and isothermal conditions. For a static load of 1.0%, the relative error between the adiabatic computation and the isothermal computation is slightly below 6%. The loss modulus E ′′ displays a slightly different profile with respect to the loading angle. Its value is at its maximum between 0 • to 15 • , where the strong strain localization around the fibers leads to high amounts of dissipated energy. Subsequently, the loss modulus decreases linearly. The effect of the static loading on the loss modulus under adiabatic conditions is more pronounced than for the storage modulus. As a decrease in temperature brings the temperature of polypropylene closer to its glass transition temperature, the dissipated energy and E ′′ increase. At a load angle of 90 • , where the self-cooling is most pronounced, even the lowest static loading of static = 0.1% leads to a 5% difference in the loss moduli. The difference increases with the static loading, reaching 13% for static = 1.0%. We conclude that the thermomechanical coupling can have a significant effect when characterizing thermoplastics-based composites using DMA. Due to the Gough-Joule effect, the effective behavior of the material, in particular E ′′ , becomes load dependent, that is, nonlinear. This is particularly pronounced for high loading frequencies, when there is no time for thermal conduction or radiation to take place and the conditions are approximately adiabatic.
To obtain precise results for real-life experiments, a strict temperature control of the specimen and low static loadings are therefore necessary.

Self-heating under harmonic loading
In the previous Section 4.3.1, we considered an oscillatory loading with a small number of cycles. In this case, the observed temperature changes were mostly due to the Gough-Joule effect caused by the static loading. However, for a high number of cycles, the dissipated energy accumulates over time and becomes the main driver of the temperature evolution.
For example, such conditions frequently occur in fatigue testing, where the self-heating of the specimen poses a major challenge. 2,3 Typically, in the first hundreds of cycles, the temperature increases in a roughly linear fashion 71 and subsequently reaches an equilibrium value when dissipation and thermal conduction reach an equilibrium state. This limits, for instance, the range of viable loading frequencies for testing. 72,73 Motivated by these findings, we take a look at the effect of the thermomechanical coupling on the dissipative characteristics of the short glass-fiber reinforced composite in the initial stage of a high cycle test. More precisely, we prescribe 100 cycles of harmonic stress-controlled uniaxial tensile loading in x-direction with a frequency of f = 10 Hz. The static stress is fixed at static = 30MPa with a stress amplitude of amp = 30MPa, corresponding to a load factor of R = min ∕ max = 0. As only a short time-frame of 10 s is considered, we assume adiabatic conditions. First, we consider the evolution of the temperature and the strain amplitude. In Figure 8(A), the minimum, maximum, and average temperature are plotted for each cycle. Initially, the mean temperature is lower than the reference ref = increase and after 25 cycles the initial cool-down is compensated. Together with the temperature, the strain amplitude increases as the material softens, cf. Figure 8(B). However, the reference value of amp for the isothermal case is reached after 42 cycles when the mean temperature has already surpassed ref . Taking a look at the minimum and maximum temperature in Figure 8(A), we observe that the large stress amplitude leads to a significant fluctuation of about 1 K for each cycle. Hence, the behavior of polypropylene fluctuates within each cycle, resulting in a slight reduction of the amplitude.
Last but not least, we take a look at the dissipation and the loss modulus for each cycle. Consistent with our observations in Section 4.3.1, the magnitude of the loss modulus, cf. Figure 9(A), is initially higher than the isothermal prediction and subsequently decreases with increasing temperature. As the temperature reaches its reference value, so does E ′′ , indicating that it is mostly unaffected by the large stress amplitude and the resulting intercyclic temperature fluctuations.
The dissipation per cycle follows a similar trend. However, it barely exceeds the isothermal reference value in the first few cycles, as the higher loss modulus is partly compensated by the lower strain amplitude. Overall, we observe that the dissipative behavior of the material changes significantly in the first cycles of a long-term harmonic excitation. At the end of 100 cycles, the loss modulus and dissipation are 16% and 12% lower, respectively, than the values predicted for the isothermal setting. Thus, when predicting the temperature changes for fatigue tests based on (28), cf. Handa et al., 69 accounting for the temperature dependence of the material is mandatory.

CONCLUSIONS
The present study was devoted to enabling the efficient computational homogenization of thermomechanically coupled materials. Based on the consistent asymptotic homogenization framework for dissipative materials, 5 we presented an efficient staggered algorithm compatible to strain or displacement-based micromechanical solvers. Due to their computational power, we focused on FFT-based solution schemes and found that best performance was achieved in combination with the Barzilai-Borwein method. Even for a composite with strong thermomechanical coupling, its iteration counts and convergence behavior hardly differed from the usual isothermal setting. The powerful class of polarization-based schemes 19,20,74 was excluded from the present work, as the complexity-reduction approach by Schneider et al. 51 may prevent the evaluation of dissipation and entropy. Further studies are necessary, to make these solvers available for thermomechanically coupled problems. In our numerical experiments, we observed that the computational overhead for the temperature-update step in the proposed algorithm was negligible. The difference in runtime between thermomechanically coupled and isothermal computations was dominated by evaluating the entropy and the dissipation, as part of the material law. In particular, computing the dissipation was costly for the chosen linear viscoelastic model, as it involves applying an inverse stiffness tensor for each Maxwell element. This lead to an increase in overall computation times by 20%-30%. However, for material laws such as J 2 -plasticity, where the dissipation is readily computed, the difference is much smaller. Overall, we conclude that the proposed algorithm enables computing the effective mechanical behavior of thermomechanically coupled materials with nearly the same computational efficiency as traditional FFT-based methods in an isothermal setting.
For the investigated glass-fiber reinforced polypropylene composites we observed that the thermomechanical coupling induced an effectively nonlinear material behavior, even though the underlying material model was linear viscoelastic. In particular, the dissipative characteristics of the materials changed significantly between the isothermal and adiabatic computations. Expanding the study of similar polymer-based lightweight-materials, such as sheet-molding compounds, 75 to include thermal effects seems promising. The presented thermomechanical solver is compatible to the interpolation approach by Köbler et al., 76 enabling the development of effective (macroscopic) surrogate models for arbitrary fiber orientations. For more general structures and material models, thermomechanical FFT-based computations may enter data-driven approaches, such as deep material networks, [77][78][79] to facilitate the simulation of components on the macroscale.
With regard to the material model of the polymer, it would be interesting to apply a free-volume-based approach for the shift factor, 63 which takes into account the pressure dependence of the viscosity. Whereas a tensile loading mechanically increases the free volume, the accompanying adiabatic cooldown, observed in this study, may weaken this effect. Investigating the interaction between these phenomena seems worthwhile to enable a thorough characterization of the thermomechanical material behavior. In addition, expanding the material model to the viscoplastic domain 57 appears attractive to investigate the influence of the plastic dissipation on the self-heating behavior of the material. As self-heating effects are particularly relevant in the context of fatigue and life-time predictions, coupling the presented thermomechanical solver with FFT-based schemes for damage 30,80 or fracture 28,29 would be of interest.

APPENDIX A. FFT-BASED SOLUTION SCHEMES FOR THE ISOTHERMAL PROBLEM
Due to their computational efficiency, we chose FFT-based algorithms for the iterative solution scheme in our staggered thermomechanical algorithm 1. In this appendix, we recast the isothermal balance of linear momentum, cf. (14), as an optimization problem. Based on this variational framework, the basic scheme, 14 the Barzilai-Borwein method, 25 and the inexact Newton method 22,23 are briefly summarized in Sections A.1-A. 3.
Let w (⋅, ) = w(⋅, , ) denote the condensed incremental potential of a generalized standard material for a fixed temperature , cf. Lahellec-Suquet. 81 The material law is defined by the hyperelastic relation = w (⋅, ).
For a prescribed macroscropic strain , we seek a compatible strain field where H 1 # (Y , R d ) denotes the first-order Sobolev space of periodic and mean-free displacement fluctuations u. By the Helmholtz decomposition, cf. chapter 12.1 in Milton's book, 83 Γ = ∇ s (div ∇ s ) −1 div is an orthogonal projector onto U. Thus, we can write the differential in the form  Note that the condition ∇W = 0 for critical points of W is equivalent to the static equilibrium div = 0 in (A2).

A.1 Basic scheme
Using the gradient descent method for solving (A2) yields the iterative scheme 23 The optimal step size for the scheme is given by cf. section 1.2.3. in Nesterov's book, 84 where + and − denote the extreme eigenvalues of the material tangent 2 w 2 evaluated over all microscopic points. Rewriting (A3) as a Lippmann-Schwinger equation , with C 0 = 1 s Id, and Γ 0 = sΓ, we see that the gradient-descent method coincides with Moulinec-Suquet's basic scheme. 14 The pseudocode for the basic scheme is summarized in Algorithm 2. ← FFT −1 ( )

5: return
The basic scheme is robust and memory efficient as it operates on a single strain-like field. However, it converges exceedingly slow for composites with high material contrast = + ∕ − .

A.2 Barzilai-Borwein method
Schneider 25 exploited the quasi-Newton-based Barzilai-Borwein 85 step size for gradient descent to accelerate the basic scheme (A3). The resulting algorithm is presented in Algorithm 3 and requires storing two strain-like fields. Note that for the first iteration (k = 0), a single iteration of the basic scheme, cf. Algorithm 2, should be used to fix the mean value of .  In practice, the Barzilai-Borwein method outperforms the basic scheme significantly, both in terms of iteration count and runtime. 25,26 For problems where the computational effort of evaluating = w ( ) is comparable to the cost of