Many engineering structures operate under severe conditions, subjected to both mechanical and thermal loads. High-performance systems demand advanced design optimization methods to meet stringent requirements for performance metrics such as weight, strength, reliability, and material cost. For instance, components in turbojet engines must endure extreme temperature gradients, as depicted in Fig. 1, ranging from ambient air temperatures to over 1000 °C. These temperature gradients often occur over very short distances, particularly in the high-pressure turbine receiving hot gas from the combustor (Lagow 2016). In such structures, temperature variation can greatly affect the thermal and mechanical properties that ultimately determine the response and performance during operation. Additionally, thermal stress fracture can lead to failure and must be mitigated.

FGMs offer a unique solution to these challenges with material properties that vary throughout a single component. FGMs can be achieved by varying either microscale geometry or material composition through space (Hofmann et al. 2014). Herein, we focus on varying material compositions to exploit variations in the material properties of multiple alloys. Distinct material compositions may exhibit differences in their temperature-dependent thermal and mechanical behaviors.

The intelligent gradation of material properties in a thermal-elastic structure permits the simultaneous management of temperature, thermal stress, and mechanical loading paths. The accommodation of continuously varying material composition during optimization can improve performance, enhance sustainability, and reduce the utilization of critical or expensive materials compared to designs composed of a single material or even discrete phases of multiple materials.

By considering the phenomenological dependence of constitutive properties on temperature during a thermal-elastic structure optimization, we can better exploit the massive design freedom offered by FGMs. Analyses and design optimizations that ignore this temperature dependency may suffer in accuracy, leading to sub-optimal designs and ultimately impacting the performance and safety of realized structures. The optimization of the geometry and material distribution in a thermal-elastic structure is an intricate problem that motivates our development of a holistic TO approach.

TO has emerged as a powerful tool for structure design, optimizing the layout of material in a design domain to best achieve performance goals under constraints. The vast literature on the topic spans methods including solid isotropic material with penalization (SIMP) (Sigmund 2001), evolutionary structural optimization (Xie and Steven 1993), level-set representations (Wang et al. 2003), and moving morphable component (Guo et al. 2014) approaches. There have been many developments in TO for both mechanical and thermal systems, which optimize mechanical response and heat transfer performance, respectively. More recently, structure optimizations for multi-physics systems, particularly thermal-mechanical coupling, have been explored and developed (Shishir and Tabarraei 2024; Qing Li and Xie 2001; Gao and Zhang 2010; Deaton and Grandhi 2016; Sigmund and Torquato 1997).

Thermal-mechanical coupling in TO typically accounts for stresses induced by thermal expansion and the influence of temperature on material properties. A weakly coupled model, where temperature impacts displacement, is considered sufficient for static equilibrium (Zhu et al. 2019). Early studies on thermal-elastic TO focused on minimizing compliance under uniform temperature fields, with or without considering strains from thermal expansion (Rodrigues and Fernandes 1995; Yang and Li 2013; Zhang et al. 2014). However, realistic scenarios involve temperature fields governed by heat conduction, influenced by the structural configuration. Additionally, methods incorporating design-dependent temperature and stress constraints, using sensitivity-based level-set methods (Deng and Suresh 2017a) or stabilizing control schemes (Meng et al. 2020), have been reported. Some studies further impose constraints on temperature fields alongside thermally induced stress constraints (Kambampati et al. 2020; Meng et al. 2022). Despite these advancements, many thermal-elastic TO approaches assume constant material properties, limiting their applicability to small temperature ranges. Research has shown that neglecting temperature-dependent properties can lead to inaccuracies even in purely thermal TO (Tang et al. 2019). Recent efforts have addressed this by incorporating temperature-dependent elastic modulus and yield strength under stress-constrained formulations, albeit with constant temperature fields (Wang et al. 2022). More comprehensive approaches account for temperature-dependent material properties and design-dependent temperature fields but often treat thermal conductivity as constant in linear thermal analyses (Chen et al. 2022).

Alongside thermal-elastic formulations, multi-material TO is also evolving and is driven by advancements in manufacturing. Multi-material TO can improve design performance with larger design spaces that accommodate multiple materials within a single structure. Combinations of different materials can be tailored to enhance desirable performance metrics such as strength, flexibility, conductivity, or material cost (Hofmann et al. 2014). Classic methods like SIMP have been extended to handle multiple phases associated with multiple compositions (Bendsøe and Sigmund 1999), but challenges remain, primarily due to computational cost and severe nonlinearities for more than two material phases. Innovative schemes like discrete material optimization (Stegmann and Lund 2005), shape function with penalization (Bruyneel 2011), and ordered SIMP (Zuo and Saitou 2017) have been developed to address these issues. Most multi-material TO works are concerned with structures consisting of discrete material phases (Huang and Li 2021; Gibiansky and Sigmund 2000), even for thermal-elastic cases (Gao et al. 2016), where continuous material gradation could mitigate thermal stress concentrations.

In contrast, advancements in the field of advanced manufacturing have led to the development of FGMs, which are capable of being fabricated into ceramics (Besisa and Ewais 2016), concretes (Torelli et al. 2020), metals (Ghanavati and Naffakh-Moosavy 2021), and polymers (Lipkowitz et al. 2022) with gradations in properties. Continuously varying materials, achieved through either composition gradients or composite structures, can exhibit smooth property transitions across a structure, allowing for the management of stress, thermal expansion, or other factors leading to structural weakness or failure. There has been considerable interest in the experimental characterization of continuously graded alloys. Notable experiments include the printing of graded stainless steel 316 L to Inconel 718 (Li et al. 2022; Yang et al. 2022), and the printing of Fe-Co and Fe-Ni systems (Chaudhary et al. 2020). The computational design of material gradient paths has also been explored within Fe–Ni–Cr systems (Kirk et al. 2018).

The ability to deposit graded material compositions during additive manufacturing introduces the possibility of treating the spatial distribution of material composition as an explicit and continuous design variable, alongside geometry. In this context, FGM TO approaches have been proposed for elastic systems, parameterized to include up to three solid phases (two solids and void) (Xia and Wang 2008; Conde et al. 2022; Silva et al. 2024). One approach employs a hybrid design representation, where geometry is implicitly defined by a level-set function and material type explicitly defined as a scalar field (Xia and Wang 2008). Another method focuses on minimizing mechanical compliance and maximum stress using a density-based FGM representation (Conde et al. 2022). More recently, TO research has extended to FGMs under thermal-elastic loading, addressing energy and stress objectives while considering thermal expansion due to imposed temperature fields (Silva et al. 2025). However, a critical gap remains in simultaneously addressing these challenges alongside heat conduction and temperature-dependent material properties.

Combining thermal-mechanical and graded material TO into one paradigm is a promising avenue for designing high-performance FGM structures. We present a density-based TO framework with automatic differentiation for improved computational efficiency and extensibility. The proposed TO approach integrates coupled thermal-elasticity, temperature-dependent material properties, and continuously graded material compositions. We summarize and contrast the capabilities of the proposed work against previous studies in Table 1. The main contributions of this paper are as follows:

• A thermal-elastic TO framework with nonlinear thermal analysis
• Concurrent design of geometry and graded material composition with consideration of design-dependent temperature fields and loading conditions
• An RBF-based interpolation scheme for thermal-elastic material properties as functions of both temperature and material composition
• Computational efficiency and adaptability through extensibility to GPUs and the adjoint method integrated with automatic differentiation

An overview of the procedure is visualized in Fig. 2. This method enables the co-design of material distributions and geometric configurations. Thermal conductivity, coefficient of thermal expansion, elastic stiffness, and yield strength are all formulated as explicit functions of temperature and material composition. Mapping alloy composition to constitutive properties has been explored with data-driven models based on experimental data (Shaheen et al. 2023; Xiong et al. 2020; Yu et al. 2021). We utilize radial basis function interpolation to evaluate material properties, facilitating continuously varying material combinations with minimal computational overhead. The accommodation of temperature-dependent thermal conductivity results in nonlinear thermal analyses, ensuring accurate temperature solutions. This is crucial since the elastic properties are sensitive to temperature as well.

Section 2 outlines the utilized and proposed methods, giving context for the numerical examples presented in Sect. 3. The numerical studies are chosen to demonstrate the effects of temperature-dependent material properties on optimized structures and the advantages of continuously grading structures under thermal-elastic loads. A case study focusing on turbine blade design for aerospace applications demonstrates the multiple advantages of our approach. Continuous material design efficiently balances the interplay between temperature variation and mechanical stresses, achieving significant improvements over single-material designs. Benefits include limiting the utilization of critical, expensive, or unsustainable materials in high-performance engineering structures.

We begin by describing the assumptions and methods developed to build our CGA design framework in preparation for the numerical results presented in Sect. 3.

2.1 Governing equations for thermal-elasticity

In the case of coupled thermal-elasticity, the governing partial differential equations (PDEs) describe the mechanical behavior of a solid body under thermal and mechanical loads. We consider a design domain  with a fixed temperature boundary  and fixed displacement boundary , such as that in Fig. 3. Heat flux  and traction loads  are applied on the boundaries if specified, as well as a body force  throughout the domain. The remaining boundaries are thermally insulated and traction-free. We consider only the equilibrium state, and all materials are assumed to be isotropic. The state equations are

where T is the temperature,  is the prescribed temperature on , k is the thermal conductivity,  is the displacement vector,  is the prescribed displacement on ,  is the stress tensor, and  is the outward normal unit vector on boundaries  and . The following constitutive and kinematic equations relate the stress tensor to the displacement and temperature fields:

where  and  are the Lamé parameters,  is the strain tensor,  is the strain due to thermal expansion,  is the coefficient of thermal expansion,  is a reference temperature which is taken to be 0 °C in this work, and  is the identity tensor. Note that the material parameters, thermal conductivity (k), thermal expansion coefficient (), and the Lamé parameters ( and ), depend on both the material distribution and temperature (T) throughout the design domain.

2.2 Forward and backward problems

These governing relations can be efficiently solved using the finite element method (FEM). This work is built upon an open-source FEM library called JAX-FEM (Xue et al. 2023), which leverages JAX (Bradbury et al. 2018), a Python library for efficient scientific computing and automatic gradient computation. The functional mapping from design configuration to system performance metrics in Fig. 2, with the intermediate FEA solution schemes, constitutes the forward problem.

TO is a case of PDE-constrained optimization, and to leverage gradient-based optimization techniques, we require sensitivities or gradients of the objective function and constraints with respect to the design variables. The calculation of these gradients constitutes the backward problem. The performance metrics typically depend on implicitly defined solution fields, in our case temperature  and displacement , but naive automatic differentiation through each step of a solution procedure for these quantities is often prohibitively expensive.

Instead, we apply the implicit function theorem (IFT) to compute sensitivities analytically. For a converged PDE solution  (either  or ) that satisfies a residual equation , where  denotes design variables (e.g., ), the total derivative of  with respect to  is given by

The end-to-end total derivative of a scalar-valued objective or constraint function  can be expressed as

For a given problem, the manual derivation of concrete expressions for these partial derivatives is tedious, so in JAX-FEM we leverage automatic differentiation to compute these quantities. To efficiently retrieve the total derivative, we evaluate the last term from left to right, utilizing a vector-Jacobian product (VJP) corresponding to the adjoint method:

Here,  is the adjoint variable, which solves a linear system of equations with the transposed Jacobian of the residual. This formulation enables efficient gradient computation without backpropagating through the entire PDE solve. We solve these adjoint linear systems for each objective and constraint and implement the corresponding VJPs for both the thermal and structural solvers directly within the automatic differentiation infrastructure of JAX. A more thorough discussion of the discrete adjoint method is provided by Cao et al. (2003), and the practical implementation of VJPs in JAX is illustrated by Xue et al. (2023).

The gradients of objective values and constraints are computed and passed to an optimizer, which proposes a new design configuration to be evaluated in the next iteration of the optimization loop. Our current choice of optimizer is the Method of Moving Asymptotes (Svanberg 1987). TO generally tends to be computationally demanding because it almost always involves solving the forward and backward problems many times to evaluate each proposed design.

2.3 Design parameterization

Consider that a design configuration is defined simultaneously by its geometry and its material composition. Geometry will be represented as a vector, , with a continuous scalar corresponding to each cell in a mesh of the entire design domain. While geometry is only meaningful as a binary variable,  for solid or  for void, relaxing the discrete problem to a continuous one permits gradient calculation through the simulation scheme, which is desirable because gradients can guide the search for an optimal geometry. To prevent singularity of the FEM stiffness matrices, the values of  are never allowed to be exactly zero, rather all values will be between 0.001 and 1.0. For a mesh with N cells,

Material composition is represented similarly to geometry as a matrix, , where each row corresponds to a cell in a mesh and the columns correspond to different basis material compositions. The continuity of material composition, however, is physically meaningful in the sense that we can conceive of a smooth chemical transition between, for example, two basis alloys. The stored image of material composition has a dimension corresponding to the number of basis materials, M, being considered, analogous to the multiple channels of a digital RGB-color image:

In a discretized domain with N cells and M basis materials, the design is parameterized by  variables: a vector  encoding geometry and a matrix  encoding material distributions. Each entry  corresponds to the unnormalized contribution of material j in a given cell i. We normalize these values within each cell to obtain physically meaningful material fractions  as

This enforces that the physical material fractions lie on the unit simplex, i.e., they represent valid barycentric coordinates, summing to one. The normalized concentrations  are then used in the constitutive interpolations to follow.

2.4 Constitutive modeling

2.4.1 Material interpolation through composition and temperature

Material behavior is inherently sensitive to temperature, exhibiting significant variation in properties such as thermal conductivity, thermal expansion, elastic stiffness, and strength. These changes can critically impact structural performance under thermal loads, potentially leading to sub-optimal or even unsafe configurations if overlooked. A comprehensive optimization framework, therefore, necessitates the integration of temperature-dependent characteristics to ensure structural integrity and functionality under anticipated thermal loads. For instance, Fig. 4 illustrates the properties of iron, nickel, and cobalt-containing alloys as functions of temperature. The functional grading of these alloys has been investigated experimentally (Li et al. 2022; Yang et al. 2022; Chaudhary et al. 2020) because they exhibit chemical compatibility, potentially making them suitable for FGM parts. It is apparent that the constitutive parameters vary just as much, if not more, across temperature as they do across chemical composition. These parameters govern the physics of the system and the performance of the structure; so to have confidence in our design performance evaluation, we must account for material behavior as a function of both composition and temperature.

The TO literature largely ignores the temperature dependence of constitutive properties. However, it is necessary to accurately evaluate design performance for systems subjected to large temperature ranges. We consider thermal conductivity k, thermal expansion , elastic stiffness E, and yield strength  as functions of composition and temperature in a unified way using radial basis function interpolation.

Radial basis function (RBF) interpolation is an approximation theory for building smooth and accurate interpolants of scattered data. It is particularly well suited for this application as it is computationally inexpensive and easily extendable to arbitrarily many dimensions, allowing for the accommodation of as many basis materials for which data are available. For the constitutive properties acknowledged as functions of temperature, we interpolate through each material basis dimension and temperature simultaneously, with the temperature dimension normalized to  such that the mean distance between adjacent data points corresponding to pure component materials is 1. Mass density and material cost are treated as constants through temperature, so the interpolants for these properties are only constructed through material composition space from the data in Table 2. For the construction of all RBF interpolants, we utilize a multiquadric kernel,

with shape parameter . For K data points, the RBF interpolant takes the form,

with  chosen such that the values of the properties f at observed points  are exactly recovered by the interpolant s,

Note that each interpolant s is formulated as an explicit function, and gradients are computed analytically and propagated through the automatic differentiation chain. RBF interpolation generally provides a smoother surface than other interpolation methods, and interpolation, as opposed to regression, offers greater accuracy, recovering true data labels when evaluated at observed data points. Material property values evaluated at intermediate compositions, such as those shown in Fig. 5, are obtained via radial basis function interpolation using a multiquadric kernel, constructed from data for the pure basis materials. This purely data-driven approach is favored over mixture rules because some relevant properties, such as yield strength, lack rigorous mixture models. Additionally, classical mixture rules such as the Hashin–Shtrikman (HS) bounds (Hashin and Shtrikman 1963) are generally inapplicable to CGAs, where material behavior is governed not only by mechanical interaction but also by underlying phase evolution and chemical interaction. A well-known counterexample is Nitinol, whose Young’s modulus is significantly less than that of both constituent element, Nickel and Titanium (Lee et al. 2018), thereby violating the HS bounds. The current interpolated trends may reflect smooth transitions, but these are artifacts of the interpolation kernel and not necessarily grounded in mixture theory.

In real applications, this design procedure could be coupled with high-throughput experimentation (Vecchio et al. 2021; Sommer et al. 2023; Feng et al. 2021) that samples material properties throughout the composition and temperature space to construct globally accurate interpolants or regression models for each material property. While the current interpolation is based on sparse data from pure material compositions and may not capture true trends at intermediate compositions, these interpolants serve as essential surrogates that enable gradient-based optimization under coupled thermal-mechanical physics. They provide smooth, differentiable approximations where classical mixture rules are unavailable or inapplicable. Ongoing work focuses on the experimental characterization of intermediate compositions, with future datasets to be integrated into material property models using regression techniques or neural networks to improve physical fidelity and generalization across composition-temperature space.

2.4.2 Material interpolation through geometry

To promote convergence to design configurations with binary, 0 and 1, density values, there are interpolation schemes that penalize intermediate density values. The well-known SIMP (Solid Isotropic Material with Penalization) interpolation scheme presents numerical challenges in the presence of design-dependent loads such as thermal expansion, inertia, or gravity, largely due to the vanishing gradient for low-density values. This makes the re-materialization of void cells highly unlikely under SIMP, which can be a drawback. Gao et al. (2016) suggest that RAMP (Rational Approximation of Material Properties) is a more appropriate material interpolation scheme in the presence of design-dependent boundary conditions. The non-zero derivative of the RAMP scheme with  at zero density, shown in Fig. 6, provides a mechanism for the optimizer to re-materialize void elements compared to the SIMP form  where .

Thus, the stiffness and thermal conductivity are interpolated according to

where  is the design density at cell i of the mesh.  and  are the RAMP penalization parameter for Young’s modulus, E, and thermal conductivity, k, respectively.  denotes the RBF interpolant for E evaluated at the physical concentration and normalized temperature of cell i, and this notation for RBF interpolants persists for all material properties. The rest of the material properties are interpolated according to

The coefficient of thermal expansion, , is modeled as constant with respect to  to prevent artificial thermal stresses from forming at the interfaces between solid and void cells (Silva et al. 2025). Yield strength is also modeled as constant across , although the von Mises stress response is computed with a penalization of Young’s modulus. The material cost and mass are linearly interpolated with respect to the design density .

2.5 Optimization formulation

Generally, there are two approaches to computational PDE-constrained optimization: optimize-then-discretize or discretize-then-optimize. We utilize the latter, discretizing the system via the finite element method before solving the PDEs. A general formulation for the CGA TO is written below.

The thermal and elastic governing PDE constraints are satisfied implicitly by the solution schemes for  and , with  and  corresponding to the global finite element stiffness matrices of the thermal and elastic problems, respectively. We seek a design geometry  and material composition distribution  that minimize an objective function J while satisfying N additional constraints . The objective function and constraints may be expressed explicitly as functions of the displacement and temperature fields, as well as the underlying design variables. The definitions of the performance metrics to be minimized or constrained are outlined in the next section.

2.6 Performance metrics — objectives and constraints

Mechanical compliance () is a measure of the elastic strain energy contained in a structure subject to a load. It quantifies the global stiffness of a structure and can be expressed as

where  is the global stiffness matrix of the elastic FEA and  is the global displacement vector evaluated without thermal expansion, achieved by letting .

Thermal compliance () is a measure of the heat transfer capacity of a structure subject to a thermal flux. It can be expressed as

where  is the final global stiffness matrix of the thermal FEA.

Target solution difference () measures the difference between the response of the structure and prescribed target solution values at certain degrees of freedom. It can be expressed as

where  is a vector encoding the target displacements at degrees of freedom of interest and with zeros at all other entries, and  is a vector of ones at the target degrees of freedom and with zeros elsewhere.

The stress design factor (DF) quantifies the load-carrying capacity of the structure before yielding anywhere, and it is the inverse of the factor of safety (FOS). A lesser design factor indicates that a design is less likely to fail under a load. Stress is an intrinsically local phenomenon, but we seek to prevent yielding everywhere in the design domain, so the metric of interest is the maximum taken over all elements. The stress design factor can be written as

where  is the von Mises stress response in cell i, and  is the yield limit in cell i. Note that the yield limit varies across cells in the mesh due to variations in both material composition and temperature.

Volume fractions (vf) and total cost are simply volume-weighted averages over all cells. They can be expressed as

where i is an index over cells in the mesh, j is an index over basis material compositions, and  corresponds to the determinant of the Jacobian matrix for the finite element corresponding to cell i in the mesh.

2.7 Stress constraint

For a structure to be feasible, the von Mises stress must not exceed the local yield strength at any element. The enforcement of these local conditions for each mesh cell can be computationally challenging (Da Silva et al. 2021; Senhora et al. 2020). Prior works have addressed this by aggregating many local constraints into a single criterion (Wang et al. 2022; Deaton and Grandhi 2016). Equation 18 represents the aggregation of local constraints imposed on each individual cell into a single global constraint via the maximum operation. To integrate this into a gradient-based optimization framework, we utilize LogSumExp (LSE) as a smooth and differentiable approximation to the maximum function:

A greater value of the parameter P corresponds to a more accurate approximation of the maximum function but reduces numerical smoothness, leading to gradient arrays that are mostly zeros even if many local constraints are violated. This can introduce instability during gradient-based optimization. Conversely, small values of P (such as  used in this work) provide greater smoothness and more stable optimization trajectories, albeit with a looser approximation of the true maximum. To correct the approximation error so that the approximate constraint is only active when the exact constraint is active, we use a correction factor c,

To stabilize the movement of the constraint boundary due to the correction factor across iterations, we follow (Yang et al. 2018) and use a stability transformation method to dampen the changes in the correction factor. This can be written as

where  corresponds to no correction of the stress constraint boundary,  corresponds to exact correction, and  corresponds to the case of damped stress constraint correction. The damped correction strategy has the advantage of enforcing the exact constraint boundary value once converged while reducing the movement of the boundary value between adjacent design iterations. The latter can cause instability and oscillation during an optimization procedure if not addressed.

In this section, we conduct several numerical studies to illustrate the proposed framework, using the Python libraries JAX (Bradbury et al. 2018) and JAX-FEM (Xue et al. 2023). The default parameters are as follows:

• The reference temperature is  .
• The RAMP penalization parameters are  for conductivity and  for Young’s modulus during the FEA to determine temperature and displacement responses, and  is used to evaluate the von Mises stress response afterward. These choices align with existing recommendations (Deaton and Grandhi 2016).
•  in the stress constraint aggregation of Equation 20.
• The damping parameter is  for the stress correction scheme of Equation 22.
• The 2D examples use fully integrated 4-node quadrilateral elements under the assumption of plane stress. The 3D examples use fully integrated 8-node hexahedral elements.
• Optimizations terminate when the relative change in objective value is less than  over five consecutive iterations, or when 200 iterations have been completed.

To optimize the structures shown throughout this work, we use an open-source Python implementation (Chandrasekhar et al. 2022) of the Method of Moving Asymptotes (Svanberg 1987). We also implement sensitivity filtering (Sigmund and Petersson 1998) with a radius of 1.5 times the mean element length. Optimizations are initialized with spatially uniform distributions of the design variables  and  such that any volume fraction constraints are initially exactly satisfied, and the composition variables  correspond to an equal mixture of each constituent alloy, .

3.1 Single material design with temperature-dependent properties

Figure 7 demonstrates the significance of considering temperature-dependent material properties, namely thermal conductivity, during the structure design loop. We consider a design domain of 20 cm by 10 cm and restrict ourselves to a single material, Incoloy 800. The design problem is formulated as the constrained optimization problem in Eq. 23, where the objective is to minimize mechanical compliance (Eq. 23a) under thermal equilibrium (Eq. 23b) and mechanical equilibrium (Eq. 23c) while satisfying a thermal compliance constraint of 2 MJ (Eq. 23d) and a volume fraction constraint that limits material usage to 60% of the design domain (Eq. 23e).

The boundary conditions are visualized in Fig. 7. For the thermal simulation, the temperature is fixed at 0°C on the left edge, and a heat flux of 100 kW/m is applied to the central 25% of the right edge. For the elastic simulation, the displacement is fixed on the left edge, and a downward traction of 100 MPa is applied to the same location as the heat source. We perform three optimizations that are identical except in their consideration of the temperature dependence of material properties. Case (a) is the simplest of the group, treating all material properties as constants through temperature, using the property values evaluated at 500C. Case (b) considers the temperature dependence of all material properties except for thermal conductivity, which is taken as a constant value corresponding to the value evaluated at 500C. Case (c) considers all material properties as functions of temperature, including thermal conductivity; therefore, the associated thermal problem is nonlinear. Note that the thermal FEA is linear for the first two examples because thermal conductivity is invariant with respect to the temperature solution. After these optimizations converge, the final designs are all evaluated with full material property temperature dependencies to recover their true performance metrics. The final structure designs are shown in Fig. 7, and their fairly evaluated objective and constraint values are found in Table 3.

The similarity between designs in Fig. 7(a) and (b) arises from the fact that both treat thermal conductivity as constant, which yields nearly identical temperature fields. As a result, the optimizer produces comparable material distributions since thermal gradients, rather than temperature-dependent elastic variation, dominate the structural performance in this case.

The similarity of  values across cases (a), (b), and (c) suggests that the temperature dependence of elastic material properties alone has a negligible effect on optimized designs in this example. However, the designs from (a) and (b) violate the thermal compliance constraint when re-evaluated with fully temperature-dependent properties. This highlights the need to include temperature-dependent thermal conductivity during the optimization loop, as it significantly affects the temperature field and overall design outcome.

The consideration of thermal nonlinearity leads to a more accurate evaluation of the temperature field and state of heat flux. Even in this simple example without stress constraints, the resulting morphology is significantly different in case (c), which considers the decrease in thermal conductivity of Incoloy 800 at low temperatures. For this reason and to satisfy the thermal compliance constraint, material was allocated to effectively sink heat into the left boundary condition. Cases (a) and (b), which used constant thermal conductivity during optimization, easily satisfied the thermal compliance constraint due to the inaccuracy of the thermal model. When re-evaluated using temperature-dependent material properties, designs (a) and (b) violate the thermal compliance constraint and exhibit a 2% increase in mechanical compliance relative to design (c), due to elevated temperatures captured in the higher-fidelity simulation. These structures are therefore less effective at dissipating heat than the one resulting from case (c).

3.2 Co-design of CGA with temperature-invariant properties

The consideration of coupled thermal and mechanical physics results in a holistic structure optimization. Even with material properties held constant with respect to temperature, different properties can be exploited in different ways, and the optimal utilization of material is highly dependent on the differences between the basis materials and the magnitudes of the loading conditions. Through these parametric studies, the multifaceted nature of multi-material systems is demonstrated. The influences of distinct material properties, mass density, conductivity, stiffness, and thermal expansion are isolated in Fig. 8. Materials A and B are defined to be identical except for one material property at a time, for each of mass density, conductivity, stiffness, and thermal expansion. The loading conditions are 0°C and 0 displacement on the bottom, while 800°C and an upward traction are applied on the top. Additionally, there is an inertial body force proportional to the mass density throughout the design domain. The corresponding optimization problem is given in Eq. 24, where the objective is to minimize the sum of squared displacements at the traction boundary (eq. 24a), subject to thermal and mechanical equilibrium (eq. 24b, eq. 24c), and per-material volume fraction constraints limiting each of the two materials to at most 30% of the design domain (eq. 24d).

Near discrete material transitions may arise even with continuous design freedom, particularly when differences in mechanical properties are examined in isolation. This is because the optimizer selects one material over the other in regions where its dominant property provides maximum benefit, for instance, low-density regions to reduce inertia or high-stiffness zones to resist deformation. For the case of thermal conductivity difference, the optimizer insulates the majority of the structure from the hot end with concentrated thermally resistive material and sinks heat into the cool end by utilizing thermally conductive material throughout the rest of the part. This reduces the average temperature of the structure and therefore the magnitude of thermal expansion. For the case of thermal expansion mismatch, the more thermally expansive material is allocated to the widest section of the geometry to reduce the net vertical component of dilation.

We extend our parametric study to examine the effects of different loading terms on the optimal structure configuration. To generate the results shown in Fig. 9, the temperature, traction, and inertial load magnitudes are manipulated such that one dominates the others in each example.

Obviously, the optimal geometry and material composition distribution throughout a structure depend on the optimization formulation. This includes each of the underlying material properties and the boundary conditions on the design domain. Figures 8 and 9 show how each of these factors in isolation can produce starkly different results, which motivates the consideration of them all to ensure a balanced design that accounts for competing effects in a holistic manner. These studies provide insight into the optimized distribution of material properties across a structure, ensuring that each region is tailored to withstand its specific loading conditions effectively.

3.3 Co-design of CGA with temperature-dependent properties

CGAs offer a pathway to achieve superior performance by continuously varying material properties, creating a seamless transition between different materials within a single component. They are especially promising for thermal-mechanical systems, where tailored thermal and mechanical behaviors can lead to structures that are inherently more resilient to thermal stresses, minimizing the risk of thermal shock and degradation over time. Moreover, CGAs can enhance the effectiveness of heat dissipation mechanisms, directly contributing to the longevity and reliability of the optimized structure.

Consider Stainless Steel 316 and Inconel 718; the former is relatively affordable, while the latter is relatively strong. There is broad interest in fabricating functionally graded parts that exploit the advantages of both (Yang et al. 2022; Zhang et al. 2021; The Minerals 2022). However, their thermal expansion coefficients differ substantially, and an abrupt interface between the two alloys can result in thermal stress concentrations at elevated temperatures. To mitigate such effects, we perform structure optimizations that maximize the von Mises factor of safety under thermal-mechanical loading. The problem is formulated in Equation 25, where the objective function (eq. 25a) to minimize is the smooth maximum of the design factor (inverse FOS). This formulation is subject to thermal and mechanical equilibrium constraints (eqs. 25b, 25c) and alloy-specific volume fraction limits of 15% each (eq. 25d).

The boundary conditions impose 0°C and zero displacement on the left edge, a leftward traction on the right edge, and a rightward body force acting throughout the domain in opposition to the traction, as visualized in Fig. 10.

The three optimizations demonstrated in Fig. 10 differ in their handling of material composition. Case (a) is a multi-material design restricted to discrete material compositions via a penalty on intermediate compositions. In other words, only pure SS316 or IN718 is permitted. Case (b) composes the design of a homogenous material throughout the domain; in particular, a 50–50 mixture of SS316 and IN718. Case (c) offers the most flexibility, allowing the material composition to continuously vary through space.

The optimized design for case (a) exhibits a safety factor to 0.98, with large thermal stress concentrations at the alloy interfaces. Furthermore, case (b) exhibits a safety factor of 2.5. However, the design does not fully exploit the superior high-temperature properties of Inconel 800. In contrast, we observe the best performance in case (c) with a safety factor of 4.8. The design provides an optimal compromise between the two extremes of material uniformity and variation.

3.4 Stress-constrained cost minimization of a turbine blade

Here, we consider the cost minimization of a turbine blade while ensuring structural integrity under thermal-mechanical loading. The optimization problem is defined in Equation 26, where the objective is to minimize the total cost of utilized material (eq. 26a) while satisfying a constraint that the von Mises factor of safety be at least 1 (eq. 26d) under thermal and mechanical equilibrium (eqs. 26b, 26c).

The design utilizes Inconel 718 and Stainless Steel 316 as basis materials; the former is stronger, while the latter is more economical. The interpolated material properties for these alloys and their mixtures are shown in Fig. 5.

The boundary conditions impose 800°C and 200 MPa at the top of the blade, and 0°C with fixed displacement at the bottom. Two types of material composition optimizations are conducted: one allowing spatially varying CGA compositions across the domain, and another enforcing a homogeneous mixture where all cells share identical composition values. Both cases are initialized with a uniform 50/50 mixture of IN718 and SS316. The blade geometry remains fixed; only the composition is optimized. The boundary conditions and resulting temperature and displacement fields for the optimized CGA case are visualized in Fig. 11.

The progression of the material cost objective and the factor of safety constraint for both optimizations are visualized in Fig. 12, and the resulting material distributions are shown in Fig. 13. Note that the objective values increase in the first several steps of the optimization procedure because the initialized configurations violated the stress constraint. This corresponds to the allocation of additional IN718 to strengthen the structures, increasing the material expense.

Figure 14 compares the results of both optimizations with a few benchmark homogeneous material designs. Both the homogeneous and graded material optimized designs outperform the pure IN718 and pure SS316 turbine blades. The IN718 case is the most expensive of the five, while the SS316 case violates the stress constraint, failing under the imposed loading conditions. The optimized homogeneous material design improves upon both, utilizing merely enough IN718 to maintain feasibility with respect to the stress constraint. Furthermore, the continuously graded material design case offers more than a 30% reduction in cost compared to the optimal homogeneous material case. With continuous material design freedom, IN718 is selectively distributed only where needed to strengthen the part near high temperatures and geometric stress concentrations. With the single homogeneous material design case, IN718 must be added everywhere, even though it may only be a local stress concentration driving this change. We also observe that a homogeneous material turbine blade with the same cost as the functionally graded design violates the stress constraint with a factor of safety of 0.42, demonstrating the advantage of FGM design.

These findings highlight the significant potential of graded alloy structures to outperform their single-material counterparts, offering critical benefits in terms of both cost efficiency and sustainability. If the cost function is adapted to represent material criticality instead of financial cost, this framework could be instrumental in reducing the use of environmentally harmful or geopolitically sensitive materials. For instance, while cobalt is known to enhance the high-temperature performance of alloys, its extraction and processing are linked to significant environmental damage, water pollution, and human rights abuses (Bamana et al. 2021). By leveraging compositional grading, it would be possible to strategically minimize cobalt usage in high-temperature components, advancing both sustainability and ethical material sourcing.

In this work, we introduced a thermal-elastic TO framework that simultaneously optimizes geometric configurations and material distributions in structures composed of CGAs. By integrating temperature-dependent material properties into the optimization process, our approach addresses key limitations of previous thermal-elastic TO methods and provides a more realistic structural performance evaluation under varying thermal conditions, and leading to better performance.

The numerical studies performed demonstrate the following critical insights:

• Significance of Temperature-Dependent Material Properties: Our results highlight the substantial impact of considering temperature-dependent properties such as thermal conductivity, thermal expansion, elastic stiffness, and yield strength. The optimized structures that account for these dependencies exhibit significantly improved thermal and mechanical performance compared to those assuming constant properties.
• Advantages of Functionally Graded Materials: The study underscores the benefits of continuous material grading. Structures optimized with FGMs show enhanced performance, namely with higher safety factors. This flexibility in material composition allows for efficient utilization of materials, balancing performance with cost and sustainability considerations.

Overall, the proposed thermal-elastic TO approach with temperature-dependent properties and FGM optimization represents a significant advancement in the design of high-performance engineering structures. The accuracy of thermal-elastic structure optimizations is enhanced by considering temperature-dependent material properties, including thermal conductivity in the nonlinear thermal analysis. The acquisition and utilization of more complete experimental material property datasets is the focus of ongoing work. Our approach is efficient and extensible, capable of utilizing GPU hardware and the automatic differentiation infrastructure provided by JAX (Bradbury et al. 2018) and JAX-FEM (Xue et al. 2023). Additionally, the development of continuous material design freedom leads to designs where compositionally graded material is strategically utilized to reduce the amount of critical material and enhance sustainability. A limitation of this method is the susceptibility of convergence to local minima when the material property interpolants are not monotonic. Future work could address this with design space exploration strategies that the current framework lacks, such as optimizing from multiple initializations or applying continuation schemes on material property models. Furthermore, advancing manufacturing constraints, mitigating deleterious phases in graded alloys, and incorporating creep and oxidation effects remain critical areas of exploration to enhance the design of high-performance thermal-elastic structures.