A multi-material topology optimization with temperature-dependent thermoelastic properties

This is a study on the development of a novel multi-material topology optimization scheme considering temperature-dependent thermoelastic properties for engineering structure design. Two cases, a three-point-bending beam under a uniform temperature field and a cantilever beam under a non-uniform temperature field, are investigated for the effects of thermoelastic properties on topology optimization. The proposed optimization scheme is compared with two existing topology optimization approaches: conventional topology optimization and thermoelastic topology optimization. The results show that the temperature-dependent elastic modulus dominantly influences the design optimization outcomes, in terms of material distribution, structural shape and compliance, while the temperature-dependent thermal expansion coefficient has a much more crucial impact on determining the material distribution and structural geometry than on compliance. Taken together, the findings demonstrate that the developed topology optimization scheme can be used to design thermally sensitive multi-materials in industrial applications, e.g. aerospace structures under high temperature and polymers in additive manufacturing.


Introduction
Structural optimization is attracting increasing attention in the rapidly developing field of computational mechanics. In the past several decades, significant progress has been made in structural optimization to acquire the best structural performance with appropriate material distribution (Chen, Cheng, and Fu 2020;Li, Steven, and Xie 1999) and specific shape under various mechanical loading conditions (Bendsøe and Sigmund 2004;Chen et al. 2018Li and Kim 2018;Sigmund 2001a;Zhang et al. 2021). While thermal effects often play a significant role in structural responses and topological design in practice, as indicated in the literature (Li, Steven, and Xie 2001), where engineering structures in a thermal environment usually experience transient heat conduction, operational temperatures and stresses may change with time, heat sources, and thermal or kinematic boundary conditions. Thus, in structural design it is essential to take into account the materials' thermoelastic behaviours.
Conventional topology optimization (CTO) for engineering structures and architecture without considering the thermal impact can be widely seen, such as the design of an upper 'beam' to span  (Beghini et al. 2014), (b) the Messerschmitt-Bölkow-Blohm (MBB) beam considering thermal loading (Madsen et al. 2016), and (c) an engine exhaust-washed structure considering thermoelastic loading (Deaton 2014).
several towers (Beghini et al. 2014), as shown in Figure 1(a). However, as mentioned above, the thermal effect often significantly influences the design results when considering practical scenarios, e.g. fire safety (Madsen et al. 2016) and high-temperature environments (Deaton 2014). Therefore, recent studies have paid considerable attention to topologically optimizing structures under thermoelastic considerations (Chung, Amir, and Kim 2020;Gao, Xu, and Zhang 2016;Li, Steven, and Xie 2001;Xia and Wang 2008;Yang et al. 2019). Li, Steven, and Xie (2001) developed an evolutionary optimization approach using an iterative calculation integrating finite element heat analysis, finite element thermoelastic analysis and topology design based on thermoelasticity. In contrast, Xia and Wang (2008) adopted the level-set method in the topology optimization of thermoelastic structures. In their study, minimization of the mean compliance of a structure was defined as the objective and the material volume was set as the constraint. A series of methods was adopted and implemented, e.g. sensitivity analysis based on a continuum model and an augmented Lagrangian multiplier method, while processing the thermoelastic design problem. Furthermore, Chung, Amir, and Kim (2020) performed a topology optimization of structures undergoing thermal and mechanical loads by incorporating thermoelastic nonlinearity into the level-set topology optimization scheme. Design examples demonstrate that temperature changes significantly influence the optimal design results, indicating that the large deformations due to thermoelastic loads provide new possibilities for controlling and manipulating the thermoelastic response via topology optimization. Yang et al. (2019) performed topology optimization of thermal protection systems for low thermal conductivity and high mechanical performance. After optimization, their designs were reconstructed from the optimized material layout and corresponding thermoelastic analyses were implemented.
Most of these studies were performed for single-material thermoelastic topology optimization (TTO), whereas nowadays increasing attention is being paid to multi-material thermoelastic design (Gao, Xu, and Zhang 2016;Vantyghem et al. 2019;Wu, Fang, and Li 2019). For instance, based on the rational approximation of material properties (RAMP) scheme, topology optimization with multiple materials under the conditions of steady-state temperature and mechanical loading was performed (Gao, Xu, and Zhang 2016). In the topology optimization of a practical case, a brickwork support bracket, it was concluded that multi-material topology design combining structural and thermal computations can provide innovative design solutions (Vantyghem et al. 2019). To solve a thermoelastic buckling problem, a density-based model integrating a sequence of topological techniques was adopted to obtain the optimal multi-material structure (Wu, Fang, and Li 2019).
Recent decades have witnessed the rapid development of new methods for multi-material topology optimization, e.g. solid isotropic material with penalization (SIMP) (Hvejsel and Lund 2011), RAMP (Gao, Xu, and Zhang 2016), level-set (Luo et al. 2009;Wang et al. 2015Wang et al. , 2017 and the alternating active-phase algorithm (Tavakoli and Mohseni 2014). However, there are still few relevant investigations on multi-material thermoelastic optimization. More importantly, engineering material properties are usually simultaneously sensitive to the environmental temperatures (e.g. engine exhaust-washed structure under high temperature, most polymers for additive manufacturing), which may lead to different structural responses (Deaton 2014;Giraldo-Londoño et al. 2020;Joshua and Ramana 2011;Li et al. 2020). For instance, the elastic modulus of one three-dimensional (3D) printed polymer (VeroWhitePlus RGD835) can be reduced by over six times when the temperature increases from 20°C to 40°C (Takezawa and Kobashi 2017), so temperature-dependent thermoelastic properties were considered to achieve accurate results (Giraldo-Londoño et al. 2020;Joshua and Ramana 2011;Li et al. 2020). Accordingly, with the temperature variation, multi-material topological optimization may produce different design results. However, to the best knowledge of the authors, the study of topological design optimization of multi-materials considering temperature-dependent thermoelastic properties is still relatively unexplored.
This research develops a TTO for multi-materials with temperature-dependent properties. A multi-material optimization interpolation scheme with an iterative updating process of temperaturedependent thermoelastic properties is built using coupled mechanical and thermal governing equations. To study the significant development of temperature-dependent thermoelastic properties in the design stage, three industrial 3D printing polymer materials with thermally distinct behaviours are adopted as sample materials. Two study cases for the typical temperature fields are accounted for: the Messerschmitt-Bölkow-Blohm (MBB) beam under the steady-state uniform temperature field and the cantilever beam subjected to the steady-state non-uniform temperature field. The effects of temperature-dependent thermoelastic properties on multi-material topological design are investigated and the temperature-dependent thermoelastic model is compared with two widely used conventional schemes: conventional topology optimization (CTO) and thermoelastic topology optimization (TTO).

Multi-material temperature-dependent thermoelastic topology optimization (TTTO)
In this study, the thermoelastic topology design with multiple materials, i.e. two design materials and one void for a solid structure, is developed using the density-based method, aiming for minimum compliance. The thermoelastic computational models with detailed discretized form for finite element analysis are given in Appendix A1 in the supplementary material. Then, the optimization problem of this study can be calculated as: where U is the global displacement vector, T is the global temperature vector, K M and K T are the global stiffness and conductivity matrices, respectively, and F T , F M and F E are the global thermal load, mechanical load and thermal expansion load (caused by thermal expansion behaviour) vectors, respectively. More details about these vectors with their calculations can be found in Appendix A1 in the supplementary material. ρ is the vector of design variables/relative densities ρ 1 and ρ 2 . During optimization calculation, ρ min and ρ max are the minimum and maximum relative density, respectively. The ρ min is set to avoid calculation singularity caused by zero solution [in this study, ρ min = 0.001 (Bendsøe and Sigmund 2004) and ρ max = 1]. v i (ρ), i = 1,2,3, represents the constraint equations for material volume fraction, which are specifically shown in Equations (2)-(4). In this study, to investigate the thermal effects on the design result and optimal material proportions, the constraint equations for the material volumes are not seriously restricted using the method proposed in Thurier et al. (2019). Assuming that the proportion of solid material (including design materials 1 and 2) is f 0 , then the void percentage is 1 − f 0 . Furthermore, a trade-off ratio f 1 is defined to determine the presence of design materials 1 and 2. Then, the volume constraints can be mathematically written as follows. The volume constraint for the void: The volume constraint for design material 1: The volume constraint for design material 2: In the above equations, V e and V are the elemental volume and total volume, respectively. This study adopts f 0 = 0.4 and f 1 = 0.625 as the volume constraints, i.e. the volume constraints for void (V f 3 ), and design materials 1 (V f 1 ) and 2 (V f 2 ) are 0.6, 0.25 and 0.15, respectively. As mentioned, these constraints are not linearly independent owing to one constraint being a function of the other two. Then, in this study, the constraints on the material volume fractions are relaxed to choose between one phase and another (Thurier et al. 2019) so as to particularly investigate the effects of temperaturedependent thermoelastic properties on the topological design results.
In the topology optimization, the design of a structure is considered to be the spatial distribution of elemental candidates throughout the design domain (Swan and Kosaka 1997). The power-law method of SIMP (Bendsøe 1989) to interpolate the material properties for element with intermediate densities is widely recommended and utilized for single-material optimization design (Bendsøe and Sigmund 2004;Kunakote and Bureerat 2011;Li, Gao, and Li 2014;Strömberg 2020). In contrast, for a multi-material problem, modified interpolation formulations based on the power-law theory for multi-material interpolation with two materials and one void can be developed (Sigmund 2001b;Sigmund and Torquato 1997).
For interpolation of the elastic modulus of element e: For interpolation of the thermal conductivity of element e: For interpolation of the thermal expansion coefficient of element e: In Equations (5)- (7), E(ρ 1 e , ρ 2 e ) is the effective (interpolated) elastic modulus as a function of design densities or variables of design materials 1 (ρ 1 e ) and 2 (ρ 2 e ), by which the dependence of the elastic modulus can be denoted to design densities. Based on the power-law theory, both ρ 1 e and ρ 2 e range from 0 to 1 in relation to the existence of void and solid material for design materials 1 and 2, respectively. E 1 and E 2 are the elastic moduli of design materials 1 and 2, respectively. Likewise, k(ρ 1 e , ρ 2 e ) and α(ρ 1 e , ρ 2 e ) are the interpolated thermal conductivity and thermal expansion coefficient, respectively. Corresponding to design materials 1 and 2, k 1 , k 2 are the thermal conductivities and α 1 , α 2 are thermal expansion coefficients. p 1 and p 2 are the penalization power coefficients, which should be greater than 1. In this study, they were defined as 3 and 1, respectively, for all cases (Sigmund 2001;Thurier et al. 2019).
Furthermore, the temperature-dependent multi-material properties, i.e. elastic moduli, can be expressed as where the temperature of the eth element T e can be calculated by Equation (A21) in Appendix A1. Then, the effective modulus of Equation (5) is modified for the temperature-dependent elastic modulus as: In addition, for interpolation of the temperature-dependent multi-material thermal conductivity of element e: Likewise, for interpolation of the temperature-dependent multi-material thermal expansion coefficient of element e: In Equations (9) and (10), k 1 T , k 2 T and α 1 T , α 2 T are the temperature-dependent thermal conductivities and thermal expansion coefficients (of design materials 1 and 2), respectively.
When incorporating the temperature-dependent thermoelastic properties into multi-material topology optimization, an updating of T e is calculated and called iteratively, leading to a simultaneous iterative updating of thermoelastic properties when implementing Equations (8)-(10). This updating is iteratively carried out until the convergence has been reached by a maximum fractional change in all elemental design variables between two consecutive iterations. If the change is less than 1%, then the design optimization is converged. In optimization, the method of moving asymptotes (MMA), developed by Svanberg (1987), is utilized for solving the optimization problem, i.e. Equation (1), owing to its stability and computational efficiency (Jung et al. 2020;Lin, Luo, and Tong 2010). The sensitivity analysis and filtering methods for multi-material temperature-dependent thermoelastic topology optimization (TTTO) are accordingly introduced in Appendix A2.

Temperature-dependent thermoelastic properties
Polymer and polymer-based multi-material composites are widely utilized in a broad range of industrial fields. Additive manufacturing techniques for 3D printing of multi-material polymers have matured and undergone rapid development (Chen and He 2020;Takezawa and Kobashi 2017). The commercial printer Stratasys J750 TM is a successful example of a multi-material printer.
In this study, three polymer-based sample materials-VeroWhitePlus RGD835 (material 1) (Takezawa and Kobashi 2017), SMPU45 (material 2) (Cheng et al. 2019) and FLX9895-DM (material 3) (Takezawa and Kobashi 2017)-were selected for optimal design and analysis. The specific temperature-dependent thermoelastic material properties (i.e. the elastic modulus and thermal expansion coefficient) of these three materials can be found in other studies (Cheng et al. 2019;Takezawa and Kobashi 2017). The general elastic modulus of polymers can be expressed as a function of temperature exponentially (Cheng et al. 2019;Tobushi et al. 1997) and, in line with the descriptions in Section 2, the temperature-dependent elastic modulus within the element e can be mathematically formulated as: where E j T (j = 1, 2, 3) stands for the temperature-dependent elastic moduli of materials 1, 2 and 3, respectively; E j represents the elastic moduli of materials 1, 2 and 3, respectively, at 20°C. a j and b j are the material fitting parameters based on the experimental results (Cheng et al. 2019;Takezawa and Kobashi 2017). The final elastic modulus versus temperature curves for sample materials are presented in Figure 2, with the specific parameters and corresponding R-squared values listed in Table 1.
The temperature-dependent thermal expansion coefficient is formulated as: In Equation (12), α j T (j = 1, 2, 3) represents the temperature-dependent thermal expansion coefficients of materials 1, 2 and 3, respectively; α j symbolizes the thermal expansion coefficients of materials 1, 2 and 3, respectively, at 20°C. m j and n j are the material fitting parameters, shown in Table 2 with R-squared values. The thermal expansion coefficient-temperature curves for the three materials are shown in Figure 3. Based on Figures 2 and 3, the elastic modulus of material 1 is similar to that of material 2 but obviously different from material 3 (while the thermal expansion coefficient of material 1 is similar to material 3 but very different from material 2), such that they can be appropriately applied to interactively compare and analyse the effects of the temperature-dependent elastic modulus and thermal expansion coefficient on the optimization results.  The thermal conductivity k of these three sample polymers is insensitive to thermal perturbations and has a similar value in all three polymers, of about 0.2376 W/m°C (Mikkelson 2014). Therefore, k = 0.2376 W/m°C was used for all three materials, interpolated in the topology optimization without considering temperature effects. The Poisson's ratios for these sample materials were assumed to be 0.3 in the topology optimization.

MBB beam under a uniform temperature field
The first case is the MBB beam (Li and Kim 2018) subjected to a steady-state uniform temperature field (Rodrigues and Fernandes 1995). The initial geometric parameters for the beam are 200 mm × 60 mm (w × h), with a uniform thickness of 1 mm, as shown in Figure 4. As for the MBB beam under three-point bending compression (Chen and Ye 2021), a concentrated load F m of 300 N was imposed at the top middle of the design space and the two bottom ends were simply supported. A uniform temperature of T was prescribed throughout the design domain. The penalization power coefficients (p 1 = 3, p 2 = 1) and the filtering radius (r min = 1. 2) were applied as constants for all numerical cases. In the subsequent figures showing the results of topology optimization, the material with the highest modulus of elasticity (E 1 ) is shown in black in reference to material 1; materials 2 and 3, with the lower modulus of elasticity (E 2 and E 3 ) are shown in blue and red, respectively. The void is represented in white.

CTO versus TTO
To begin with, a comparative study of the CTO (without considering thermal effects) and the TTO was performed, and the results are presented in Figure 5. It can be found that the final compliance between the CTO and TTO is marginal, irrespective of the materials used, while the final optimized configurations using materials 1 and 2 change rapidly when taking into account the thermal effects, but those using materials 1 and 3 vary only slightly. This is because materials 1 and 3 have similar thermal expansion coefficients, contributing to the same-level priority in optimization, even under thermal influences, while materials 1 and 2 have significantly different thermal expansion coefficients, resulting in apparent different topological configurations and material distributions. In addition, owing to larger elastic moduli overall, the final compliances using materials 1 and 2 are smaller than those using materials 1 and 3. It can be concluded that in multi-material topological design, the TTO must be considered for these sample materials when they have distinctive thermal properties.

Effects of temperature-dependent thermoelastic properties on multi-material topology optimization
The effects of temperature-dependent thermoelastic properties on multi-material topology optimization were investigated using the TTO (E j and TTTO with only elastic moduli (E T j = a j E j e (b j T e ) , α j T = α j ). As shown in Table 3, either the temperature-dependent thermal expansion coefficients or the temperature-dependent elastic moduli influence the topology optimization design results significantly, even though all calculations were performed at T = 40°C. In particular, for the temperature-dependent elastic moduli, their properties decrease rapidly as the environmental temperature is considered, leading to increased compliance with much higher volumes (using materials 1 and 2) or stiffer materials (materials 1 and 3) being prioritized in optimization. Specifically, the material volumes of TTTO using materials 1 and 2 are both increased to over 0.32 and 0.23, respectively, while the stiffer material 1 becomes predominant (with a volume of approximately 0.4) and the volume of material 3 is greatly decreased (to approximately 0.06) for TTTO using materials 1 and 3. The temperature-dependent thermoelastic properties thus exert a significant impact on the topology optimization results and should not be ignored, especially for thermally sensitive materials with alternative distinguishable characteristics. Figure 6 shows the deformation results at T = 40°C of TTO and TTTO, respectively. By considering the temperature-dependent thermoelastic properties using TTTO, the deflection along the y-axis is increased from below 8 mm to over 22 mm, accounting for the physical thermal material behaviour in the design stage, which further indicates the effectiveness and significance of developing temperature-dependent thermoelastic properties in the topology optimization. Table 3. Effects of temperature-dependent thermoelastic properties on multi-material topology optimization of the Messerschmitt-Bölkow-Blohm (MBB) beam at T = 40°C.    Figure 7 illustrates the effects of temperature (T = 20°C, 30°C and 40°C) on temperature-dependent thermoelastic multi-material topology optimization. Based on Figures 2 and 3, as the temperature increases, the decreased elastic moduli with increased thermal expansion coefficients of the materials lead to maximized compliance of the MBB beam (the compliance increased more than three times as the temperature increased from 20°C to 40°C). Besides, it can be observed that material 1 is always positioned as the main upper frame to resist the load, regardless of the temperature variation or material combination, owing to its having the largest elastic modulus with the stiffest performance.

Cantilever beam in non-uniform temperature field
The second case is the cantilever beam in a steady-state non-uniform temperature field. As shown in Figure 8, the design scope has a rectangular geometry of 200 mm × 120 mm (w × h) with a uniform thickness of 1 mm. A concentrated downward force of 300 N is imposed at the bottom right corner, while the left-hand side edges are totally constrained, and a steady-state uniform temperature T h is imposed, while the other edges are defined with a constant temperature of T l = 20°C. (Note that Figure 9. Optimization process and results of thermoelastic topology optimization (TTO) and temperature-dependent thermoelastic topology optimization (TTTO) for the cantilever beam using (a) materials 1 and 2, and (b) materials 1 and 3. the conductivity in the left-hand-side edge is enhanced by 1000 times to investigate the temperature effects.) Figure 9 presents a comparison of the results of TTO and TTTO for a cantilever beam at T h = 40°C. It is evident that the optimized configurations are different when using these two models and different material combinations. As such, material 2 is mainly selected for the upper structural components, while material 1 is chosen primarily for the bottom components when using a combination of materials 1 and 2, whereas material 1 forms the main structural frame when using materials 1 and 3. In addition, the structural compliance using TTTO is slightly lower than that using TTO. Figure 10 presents the deformation of a cantilever beam using TTO and TTTO. The deflection along the y-axis differs insignificantly, with a relative variation ratio of 1-5%, but the optimal material distribution changes noticeably.

Effects of temperature-dependent thermoelastic properties on multi-material topology optimization
The effects of temperature-dependent thermoelastic properties on multi-material topology optimization of a cantilever beam subjected to a non-uniform temperature field at T h = 40°C were investigated, as shown in Table 4. When considering the effects of temperature-dependent multimaterial properties, particularly the temperature-dependent elastic modulus, a significant difference can be observed in the optimization results in terms of material distribution, geometric profile and structural compliance. Specifically, with TTTO, the structural compliances dropped by 8.7% and 11.2% using materials 1 and 2 and materials 1 and 3, respectively, compared to those with TTO. With regard to the temperature-dependent thermal expansion coefficients, they still have relatively little influence on the compliance, but they obviously affect the material distribution and structural configuration. Figure 11 shows the non-uniform temperature fields at T h = 20°C, 30°C and 40°C, and the corresponding optimized results of the cantilever beam using TTTO. By calculation, with the increment Figure 10. Deformation results of thermoelastic topology optimization (TTO) and temperature-dependent thermoelastic topology optimization (TTTO) for the cantilever beam at T = 40°C. Table 4. Effects of temperature-dependent thermoelastic properties on multi-material topology optimization of the cantilever beam at T h = 40°C. Figure 11. Non-uniform temperature fields and effects of temperature (T = 20, 30 and 40°C) on temperature-dependent thermoelastic multi-material topology optimization of the cantilever beam.

Effects of temperature on TTTO
in temperature T h from 20°C to 40°C, the compliance decreases by 9.5% and 7.1% for materials 1 and 2 and materials 1 and 3, respectively.

Conclusions
In this study, a TTTO scheme has been developed successfully for multi-material structural parts. Some main conclusions can be drawn: • The temperature-dependent elastic modulus affects the multi-material optimization results significantly, in terms of the material distribution, geometric profile and structural compliance. In contrast, the thermal expansion coefficient has a greater influence on the material distribution and structural shape than on the compliance. For a reliable and accurate design, the proposed thermoelastic multi-material topology optimization scheme, TTTO, should be employed, especially for structures made of thermally sensitive multi-materials where component materials have different properties from each other. • The comparative study proved that TTO must be considered for materials with distinguishable thermal properties in multi-material topological design compared to CTO. Furthermore, the yaxis deflection clearly increased with TTTO compared to that using TTO. • The effects of temperature on TTTO show that increasing T from 20°C to 40°C results in decreased elastic moduli, although with increased thermal expansion coefficients of the materials. This leads to increased compliance of the MBB beam (over three times greater as the temperature increases from 20°C to 40°C). All of the results consolidate the importance of developing temperature-dependent multi-material properties in multi-material structural optimization. • A steady-state non-uniform temperature field was developed for designing a cantilever beam. The results show that the structural compliance dropped by 8.7% and 11.2% for materials 1 and 2 and materials 1 and 3, respectively, when comparing the TTTO method with the TTO method. The effect of temperature was studied for TTTO, demonstrating that the compliance decreased by over 7% as the temperature along the predefined edge T h increased from 20°C to 40°C.
With regard to future work and direct applications of the proposed thermoelastic multi-material topology optimization scheme, extensions of the current model with more materials (i.e. more than three) for particular applications would be appropriate, and experimental studies on 3D printed sample materials are ongoing and will be reported by the authors soon.