Multi-Material Topology Optimization of Thermo-Elastic Structures with Stress Constraint

Abstract

This paper proposes a multi-material topology optimization formulation for thermo-elastic structures considering coupled mechanical and uniform thermal loads. The ordered-SIMP multiple materials interpolation model is introduced, combined with examples considering structural volume minimization under stress constraints. The p-norm function with the adjusted coefficient was adopted to measure the global maximum stress. The adjoint variable method is presented to discuss the sensitivity of stress constraints, and the method of moving asymptotes (MMA) is utilized to update the design variables. The results demonstrate that clear topologies are obtained for complicated multiple material combinations with various temperature values. Meanwhile, the optimized configuration with stress constraints has clear sensitivity to uniform temperature variations. Therefore, the proposed model demonstrates the necessity of a thermo-elastic model influenced by temperature in optimization.

1. Introduction

Various mechanical structures have complex working environments, and often consider the effects of coupled multi-physical fields. Thermo-elastic structures coupled with mechanical and thermal loads are more general research objects, such as battery dissipation, engine exhaust mechanisms, high-speed aircraft thermal protection systems, etc. To meet their structural strength requirements and achieve lightweight demand, the topology optimization method of thermo-elastic structures has emerged and gradually been applied as an advanced design tool.

Topology optimization aims to find the optimal distribution of structural materials under specified boundary conditions. In recent years, based on the development of heat conduction structure research [1,2], topology optimization has been expanded rapidly to generate various innovative thermo-elastic configurations. Rodrigues et al. [3] first proposed thermo-elastic topology optimization by the homogenization method for structural compliance minimization. Li et al. [4] further developed the evolutionary optimization method based on element thickness to minimize displacement in thermo-elastic structures. Xia et al. [5] used the level set method to study the topological optimization design with volume constraints and compliance objectives. Zhu et al. [6] proposed the topology optimization of a thermodynamic formulation considering temperature constraints. It can be observed that thermo-elastic topology optimization is usually performed with the object of minimizing structural compliance. However, the topology optimization models may be inefficient to achieve the structural strength criterion. Pedersen et al. [7] questioned the topological optimization model of thermo-elastic structures with compliance minimization, but the structural strength was not necessarily optimized. The work conducted by Deaton et al. [8] also concluded that the compliance topology optimization model is not suitable for thermo-elastic topology optimization, especially when the thermal load is comparable to the mechanical load.

Topology optimization with stress constraints currently suffers from three main problems: the singularity phenomenon, the local nature of stresses, and the highly nonlinear behavior of stress constraints. Cheng and Jiang [9] presented a mathematical explanation for the singular problem in topology optimization, which is mainly due to stress discontinuity. Duysinx and Bendsøe [10] pointed out that a singular phenomenon exists in stress-constrained topology optimization. Recommended approaches to deal with a singular phenomenon include the ε relaxation technique [11,12], qp relaxation technique [13,14], and the stress penalty function [15]. The existing solutions for the computational surge of sensitivity analysis due to local stresses are generally through aggregation functions, such as the p-norm and K-S functions, which transform the local stress constraint into a global stress constraint measurement to reduce the computational effort [16,17,18]. The highly nonlinear behavior of stress refers to the stress concentration areas that can lead to large stress gradients, such as corners, depressions, and other special areas of the structure [14,19,20], for which suitable optimization algorithms are needed. The abovementioned research illustrated the challenge of stress-constrained topology optimization for thermo-elastic structures.

However, the stress-constrained topology optimization of thermo-elastic structures predicated on a single material is often challenging to achieve integrated optimal performance. The emergence of additive manufacturing technology has made it possible to manufacture the mechanical structures of multiple materials obtained by topology optimization [21,22,23,24]. In this context, Scholars have carried out extensive research in multi-material topology optimization [25,26]. Doan et al. [27] proposed a new multi-material topology optimization algorithm, which can economically and effectively select various material combinations to meet structural stiffness requirements. Alfouneh [28] studied the topology optimization problem of multilayer composite structures under static load. The proposed methods have been applied to three objective functions: minimizing compliance, maximizing mutual strain energy, and minimizing full-stress designs. Wang et al. [29] developed a topology optimization model for the multi-material level set method to solve various complex multiple material optimization problems. Gao et al. [30] combined the homogeneous multiple material interpolation to achieve the optimal topology optimization. Long et al. [31] used a sequential quadratic programming algorithm to solve multi-material topology optimization for the transient heat conduction problem. Li et al. [32] integrated functional gradient constraints into multi-material topology optimization for a transient heat conduction structure. Therefore, multi-material topology optimization is considered in the stress-constrained thermos-elastic design.

This paper proposes a mathematical model of multi-material topology optimization based on the stress constraints of thermo-elastic structures. Compared with the traditional topology optimization model of dimensionless calculations, the actual parameter values of materials are used for analysis. The ordered-SIMP multiple materials interpolation model [33] is combined with structural volume minimization. The stress singularity is solved by applying a penalty format to the stress values and constructing a p-norm aggregation function to establish a global stress constraint. The adjoint variable method is used to derive the sensitivity of stress constraint, and the MMA is selected to iteratively update the design variables. Finally, the effectiveness of the proposed method is demonstrated through several numerical examples.

2. Multi-Material Interpolation Model

For the multi-material topology optimization of thermo-elastic designs, it is necessary to discretize the design domain and use a normalization scheme to deal with the material properties of each element. The normalization process is expressed as follows:

$$\{\begin{array}{l}{\rho }_j^N={\rho }_j/{\rho }_0\\ E_j^N=E_j/E_0\\ {\xi }_j^N={\xi }_j/{\xi }_0\end{array} (j=1,2,\cdots ,n)$$

(1)

where ${\rho }_j^N$, $E_j^N$ and ${\xi }_j^N$ are the normalized density, elastic modulus and thermal stress coefficient of material j, respectively. ρ_j, E_j and ξ_j are the actual density, elastic modulus and thermal stress coefficient of material j, respectively, n is the total number of materials, and ρ₀, E₀ and ξ₀ are the maximum values of density, elastic modulus and thermal stress coefficient of all materials, respectively, where ξ₀ is defined as:

$${\xi }_0=c_0\times E_0$$

(2)

where c₀ is the maximum value of the thermal expansion coefficient.

Based on the ordered-SIMP model, the multiple material interpolation scheme for elastic modulus is defined as:

$$\{\begin{array}{l}{E}_e^i(x_e)=A_E^i{(x_e)}^{{\mu }_i}+B_E^i\hspace{1em}{x}_e\in \left[{\rho }_j^N,{\rho }_{j+1}^N\right]\\ A_E^i=\frac{E_j^N-E_{j+1}^N}{{({\rho }_j^N)}^{{\mu }_i}-{({\rho }_{j+1}^N)}^{{\mu }_i}} (i=\mathrm{I},\mathrm{II})\\ B_E^i=E_j^N-A_E^i{({\rho }_j^N)}^{{\mu }_i}\end{array}$$

(3)

where μ_i denotes the ordered-SIMP penalty parameter, μ_I and μ_II indicate the penalty factors for elastic modulus and element stress interpolation, respectively, $E_e^i$ represent the interpolated multi-material elastic modulus, $A_E^i$ and $B_E^i$ are the scale coefficient and translation coefficient, respectively, and ${\rho }_{j+1}^N$ represent the normalized density of material j and j + 1.

Similarly, for the thermal stress coefficient, the multi-material interpolation model can be expressed as:

$$\{\begin{array}{l}{\xi }_e(x_e)=A_{\xi }{(x_e)}^{{\mu }_{\xi }}+B_{\xi }\hspace{1em}{x}_e\in \left[{\rho }_j^N,{\rho }_{j+1}^N\right]\\ A_{\xi }=\frac{{\xi }_j^N-{\xi }_{j+1}^N}{{({\rho }_j^N)}^{{\mu }_{\xi }}-{({\rho }_{j+1}^N)}^{{\mu }_{\xi }}}\\ B_{\xi }={\xi }_j^N - A_{\xi }{({\rho }_j^N)}^{{\mu }_{\xi }}\end{array}$$

(4)

where μ_ξ represents the penalty factor for thermal stress coefficient, and ξ_e represents the interpolated multi-material thermal stress coefficient.

3. Multi-Material Topology Optimization Method for Thermo-Elastic Structures

3.1. Description of Thermo-Elastic Structures

Figure 1 depicts a generalized two-dimensional thermo-elastic structure consisting of non-design and design domains. A displacement boundary condition, Γ_d, and an external mechanical load, ${\mathit{F}}^m$, are presented. A uniform temperature variation field, ΔT, is specified in the design domain that generates thermal load. The design domain needs to be discretized into a finite number of elements.

For the thermo-elastic structure, the overall equilibrium equation can be expressed by:

$$\mathit{K}(\mathit{x})\mathit{U}(\mathit{x})={\mathit{F}}^m+{\mathit{F}}^{th}(\mathit{x})$$

(5)

where K(x), U(x) and ${\mathit{F}}^m$ are the respective global stiffness matrix, node displacement vector and node-independent external mechanical load vector, respectively, and ${\mathit{F}}^{th}(\mathit{x})$ is the overall thermal load vector.

The global stiffness matrix is assembled from the element stiffness matrix that can be expressed as:

$$\mathit{K}(\mathit{x})={\displaystyle \sum _{e=1}^{N_e}{\mathit{K}}_e}$$

(6)

where N_e denotes the total number of elements.

The element stiffness matrix is given by:

$${\mathit{K}}_e={\displaystyle {\int }_{{\Omega }_e}{\mathit{B}}_e^{\mathrm{T}}}{\mathit{D}}_e{\mathit{B}}_{e}hd{\Omega }_e$$

(7)

where B_e is the strain displacement matrix of elements obtained from the element shape function matrix, h is the thickness of the planar element, and D_e is the material elasticity matrix of the element, expressed as:

$${\mathit{D}}_e=E_e^{\mathrm{I}}(x_e)E_0{\mathit{D}}_0=(A_E^{\mathrm{I}}{(x_e)}^{{\mu }_{\mathrm{I}}}+B_E^{\mathrm{I}})E_0{\mathit{D}}_0$$

(8)

where D₀ is the coefficient matrix for an element with unit elastic modulus.

Substituting Equations (7) and (8) into Equation (6), the expression of K(x) can be obtained as:

$$\mathit{K}(\mathit{x})={\displaystyle \sum _{e=1}^{N_e}{E}_e^{\mathrm{I}}(x_e)E_0{\displaystyle {\int }_{{\Omega }_e}{\mathit{B}}_e^{\mathrm{T}}}{\mathit{D}}_0{\mathit{B}}_{e}hd{\Omega }_e}$$

(9)

The overall thermal load vector is assembled from the element heat load, defined as:

$${\mathit{F}}^{th}(\mathit{x})={\displaystyle \sum _{e=1}^{N_e}{\mathit{F}}_e^{th}(x_e)}$$

(10)

Considering uniform temperature variation, the thermal load can be expressed as:

$${\mathit{F}}_e^{th}(x_e)={\xi }_e(x_e){\xi }_0\Delta T{\displaystyle {\int }_{{\Omega }_e}{\mathit{B}}_e^{\mathrm{T}}}{\mathit{D}}_0\mathit{\omega }d{\Omega }_e=(A_{\xi }{(x_e)}^{{\mu }_{\xi }}+B_{\xi }){\xi }_0\Delta T{\displaystyle {\int }_{{\Omega }_e}{\mathit{B}}_e^{\mathrm{T}}}{\mathit{D}}_0\mathit{\omega }d{\Omega }_e$$

(11)

where ΔT is the element temperature value, ω is the thermal expansion coefficient vector that is defined as [1, 1, 0]^T for 2D problems and [1, 1, 1, 0, 0, 0]^T for 3D problems.

The stress singularity phenomenon exists in the stress-constrained topology optimization. Therefore, a different power exponent is used to penalize the stress values. The penalized element stress is established as:

$${\mathit{\sigma }}_e=E_e^{\mathrm{II}}(x_e){\mathit{\sigma }}_e{}_0=(A_E^{\mathrm{II}}{(x_e)}^{{\mu }_{\mathrm{II}}}+B_E^{\mathrm{II}}){\mathit{\sigma }}_e{}_0$$

(12)

where σ_e is the interpolated element stress, σ_e0 is the stress vector of the eth element using the stress value at the center of each element, and can be calculated as:

$${\mathit{\sigma }}_{e0}=E_0[{\mathit{D}}_0{\mathit{B}}_e{\mathit{U}}_e-{\mathit{D}}_0{c}_0\mathit{\omega }\Delta T]$$

(13)

In the 2D problems, the element stress vector σ_e0 consists of the x-direction stress σ_ex, the y-direction stress σ_ey and the shear stress of the element τ_exy.

$${\mathit{\sigma }}_{e0}={\left[{\sigma }_{ex} {\sigma }_{ey} {\tau }_{exy}\right]}^{\mathrm{T}}$$

(14)

Similarly, the 3D element stress vector can be expressed as:

$${\mathit{\sigma }}_{e0}={\left[{\sigma }_{ex} {\sigma }_{ey} {\sigma }_{ez} {\tau }_{exy} {\tau }_{eyz} {\tau }_{ezx}\right]}^{\mathrm{T}}$$

(15)

Adopting the fourth strength theory as the material failure criterion, it is necessary to find the von Mises equivalent stress, which can be obtained from the three components of the element stress vector. To facilitate subsequent calculations in the stress constraint sensitivity, it can be rewritten into matrix form as:

$${\sigma }_e^{vm}=\sqrt{{({\mathit{\sigma }}_e)}^{\mathrm{T}}\mathit{M}{\mathit{\sigma }}_e}$$

(16)

where ${\sigma }_e^{vm}$ denotes the von Mises stress at the centroid of the eth element, M is the matrix of stress coefficients; in the 2D plane stress case, this takes the following form:

$$\mathit{M}=\left[\begin{array}{ccc}1& -1/2& 0\\ -1/2& 1& 0\\ 0& 0& 3\end{array}\right]$$

(17)

Additionally, in the 3D spatial stress case, M can be expressed as:

$$\mathit{M}=\left[\begin{array}{cccccc}1& -1/2& -1/2& 0& 0& 0\\ -1/2& 1& -1/2& 0& 0& 0\\ -1/2& -1/2& 1& 0& 0& 0\\ 0& 0& 0& 3& 0& 0\\ 0& 0& 0& 0& 3& 0\\ 0& 0& 0& 0& 0& 3\end{array}\right]$$

(18)

3.2. Mathematical Model of Topology Optimization

For the topology optimization of thermo-elastic structures, the structural volume depreciation is taken as the objective function, and the maximum stress value is less than the yield stress limit as the constraint to satisfy the static strength failure criterion. Therefore, the stress-constrained topology optimization model for thermo-elastic structures is described as follows:

$$\{\begin{array}{l}\mathrm{find}:\mathit{x}={\left\{x_1,x_2,\cdots ,x_{N_e}\right\}}^{\mathrm{T}}\hfill \\ \min :V(\mathit{x})={\displaystyle \sum _{e=1}^{N_e}{x}_e{v}_e}\hfill \\ \mathrm{s}.\mathrm{t}. \mathit{K}(\mathit{x})\mathit{U}(\mathit{x})={\mathit{F}}^m+{\mathit{F}}^{th}(\mathit{x})\hfill \\ \hspace{1em}\hspace{1em}{\left[{\sigma }_e^{vm}\right]}_{\max }\le {\sigma }_s\hfill \\ \hspace{1em}\hspace{1em}0<x_{min}\le x_e\le 1(e=1,2,\dots ,N_e)\hfill \end{array}$$

(19)

where x is the design variable vector, V is the overall structural volume, v_e is the element volume, N_e is the total number of elements, and σ_s is the material yield stress limit.

3.3. Global p-Norm Stress Measure

If stress constraints are set on each element, the number of constraints is enormous, which leads to the inefficient computation of stress-constrained topology optimization. Meanwhile, the maximum stress is not fixed during the optimization process; thus, the optimization problem has difficulty converging quickly. Therefore, the aggregation p-norm function is used to integrate all von Mises stresses into a global stress model which is approximately equal to the maximum von Mises stress, expressed as:

$${\sigma }^{PN}={\left({\displaystyle \sum _{e=1}^{N_e}{\left(\frac{{\sigma }_e^{vm}}{{\sigma }_s}\right)}^p}\right)}^{\frac{1}{p}}$$

(20)

where σ^PN and p denote the integrated global stress and the stress norm parameter. When p tends to infinity, σ^PN is approximately equal to the maximum von Mises stress, expressed as:

$$\max ({\sigma }_e^{vm})/{\sigma }_s\le {\sigma }^{PN}\le 1$$

(21)

However, the larger value of p increases the degree of nonlinearity of the aggregation function, leading to oscillation problems in the optimization process. Otherwise, the value of p is small, and the aggregation function cannot capture the maximum value of the stress constraint. An improved constraint equation with a revised factor is introduced to reduce the deviation of the global stress constraint, expressed as:

$${\overline{\sigma }}^{PN}=cp{\sigma }^{PN}\le 1$$

(22)

where cp is the revised factor. Before each iteration of optimization, the factor can be calculated as:

$$cp=\frac{\max \left({\sigma }_e^{vm}\right)}{{\sigma }_s\cdot {\sigma }^{PN}}$$

(23)

4. Sensitivity Analysis

It is worth noting that the revised factor updates the sensitivity of the global stress relative to the element design variable. According to the chain rule, the derivative information can be obtained as follows:

$$\frac{\partial {\overline{\sigma }}^{PN}}{\partial x_e}=cp\frac{\partial {\sigma }^{PN}}{\partial x_e}={\displaystyle \sum _{e=1}^{N_e}cp\frac{\partial {\sigma }^{PN}}{\partial {\sigma }_e^{vm}}{\left(\frac{\partial {\sigma }_e^{vm}}{\partial {\mathit{\sigma }}_e}\right)}^{\mathrm{T}}\left(\frac{\partial {\mathit{\sigma }}_e}{\partial x_e}+\frac{\partial {\mathit{\sigma }}_e}{\partial {\mathit{U}}_e}\frac{\partial {\mathit{U}}_e}{\partial x_e}\right)}$$

(24)

The sensitivity information for solving the global stress can be obtained by combining the derivative of the p-norm function relative to the von Mises stress, the derivative of the von Mises stress relative to the stress component and the derivative of the stress component relative to the design variable.

The derivative of the p-norm function concerning the von Mises stress can be obtained from Equation (20), written as:

$$\frac{\partial {\sigma }^{PN}}{\partial {\sigma }_e^{vm}}={\left({\displaystyle \sum _{e=1}^{N_e}{\left(\frac{{\sigma }_e^{vm}}{{\sigma }_s}\right)}^p}\right)}^{\frac{1}{p}-1}{\left(\frac{{\sigma }_e^{vm}}{{\sigma }_s}\right)}^{p-1}\frac{1}{{\sigma }_s}$$

(25)

For 2D and 3D problems, the derivative of the element stress for the stress component can be defined as:

$$\{\begin{array}{c}\frac{\partial {\sigma }_e^{vm}}{\partial {\sigma }_{ex}}=\frac{1}{2{\sigma }_e^{vm}}\left(2{\sigma }_{ex}-{\sigma }_{ey}\right)\\ \frac{\partial {\sigma }_e^{vm}}{\partial {\sigma }_{ey}}=\frac{1}{2{\sigma }_e^{vm}}\left(2{\sigma }_{ey}-{\sigma }_{ex}\right)\\ \frac{\partial {\sigma }_e^{vm}}{\partial {\tau }_{exy}}=\frac{3{\tau }_{exy}}{{\sigma }_e^{vm}}\end{array}$$

(26)

$$\{\begin{array}{c}\frac{\partial {\sigma }_e^{vm}}{\partial {\sigma }_{ex}}=\frac{1}{2{\sigma }_e^{vm}}\left(2{\sigma }_{ex}-{\sigma }_{ey}-{\sigma }_{ez}\right)\\ \frac{\partial {\sigma }_e^{vm}}{\partial {\sigma }_{ey}}=\frac{1}{2{\sigma }_e^{vm}}\left(2{\sigma }_{ey}-{\sigma }_{ex}-{\sigma }_{ez}\right)\\ \frac{\partial {\sigma }_e^{vm}}{\partial {\sigma }_{ez}}=\frac{1}{2{\sigma }_e^{vm}}\left(2{\sigma }_{ez}-{\sigma }_{ex}-{\sigma }_{ey}\right)\\ \frac{\partial {\sigma }_e^{vm}}{\partial {\tau }_{exy}}=\frac{3{\tau }_{exy}}{{\sigma }_e^{vm}}\\ \frac{\partial {\sigma }_e^{vm}}{\partial {\tau }_{exz}}=\frac{3{\tau }_{exz}}{{\sigma }_e^{vm}}\\ \frac{\partial {\sigma }_e^{vm}}{\partial {\tau }_{ezx}}=\frac{3{\tau }_{ezx}}{{\sigma }_e^{vm}}\end{array}$$

(27)

The derivative of the stress component to the design variable can be written as:

$$\frac{\partial {\mathit{\sigma }}_e}{\partial x_e}=\frac{d{E}_e^{\mathrm{II}}(x_e)}{d{x}_e}{E}_0({\mathit{D}}_0{\mathit{B}}_e{\mathit{U}}_e-{\mathit{D}}_{0}cp\mathit{\omega }\Delta T)+E_e^{\mathrm{II}}(x_e)E_0{\mathit{D}}_0{\mathit{B}}_e\frac{\partial {\mathit{U}}_e}{\partial \mathit{U}}\frac{\partial \mathit{U}}{\partial x_e}$$

(28)

By introducing the Lagrange product factor, the Lagrange function of the stress constraint is constructed as:

$$L={\overline{\sigma }}^{PN}-{\mathit{\lambda }}^{\mathrm{T}}\left(\mathit{K}(\mathit{x})\mathit{U}-{\mathit{F}}^m-{\mathit{F}}^{th}(\mathit{x})\right)$$

(29)

The sensitivity to the design variable x_e is defined as:

$$\frac{\partial L}{\partial x_e}=\frac{\partial {\overline{\sigma }}^{PN}}{\partial x_e}-{\mathit{\lambda }}^{\mathrm{T}}\left(\frac{\partial \mathit{K}(\mathit{x})}{\partial x_e}\mathit{U}+\mathit{K}(\mathit{x})\frac{\partial \mathit{U}}{\partial x_e}-\frac{\partial {\mathit{F}}^m}{\partial x_e}-\frac{\partial {\mathit{F}}^{th}(\mathit{x})}{\partial x_e}\right)$$

(30)

The total sensitivity can be obtained by combining the previous sensitivity information, defined as:

$$\begin{array}{ll}\frac{\partial L}{\partial x_e}=& {\displaystyle \sum _{e=1}^{N_e}cp\frac{\partial {\sigma }^{PN}}{\partial {\sigma }_e^{vm}}{\left(\frac{\partial {\sigma }_e^{vm}}{\partial {\mathit{\sigma }}_e}\right)}^{\mathrm{T}}}\frac{d{E}_e^{\mathrm{II}}(x_e)}{d{x}_e}{E}_0({\mathit{D}}_0{\mathit{B}}_e{\mathit{U}}_e-{\mathit{D}}_{0}cp\mathit{\omega }\Delta T)-{\mathit{\lambda }}^{\mathrm{T}}\left(\frac{\partial \mathit{K}(\mathit{x})}{\partial x_e}\mathit{U}-\frac{\partial {\mathit{F}}^m}{\partial x_e}-\frac{\partial {\mathit{F}}^{th}(\mathit{x})}{\partial x_e}\right)\\ & +\left[{\displaystyle \sum _{e=1}^{N_e}cp\frac{\partial {\sigma }^{PN}}{\partial {\sigma }_e^{vm}}{\left(\frac{\partial {\sigma }_e^{vm}}{\partial {\mathit{\sigma }}_e}\right)}^{\mathrm{T}}}{E}_e^{\mathrm{II}}(x_e)E_0{\mathit{D}}_0{\mathit{B}}_e\frac{\partial {\mathit{U}}_e}{\partial \mathit{U}}-{\mathit{\lambda }}^{\mathrm{T}}\mathit{K}(\mathit{x})\right]\frac{\partial \mathit{U}}{\partial x_e}\end{array}$$

(31)

The adjoint variables can be obtained through the adjoint variable equation, defined by:

$$\mathit{K}(\mathit{x})\mathit{\lambda }={\displaystyle \sum _{e=1}^{N_e}cp\frac{\partial {\sigma }^{PN}}{\partial {\sigma }_e^{vm}}{E}_e^{\mathrm{II}}(x_e)E_0{\left(\frac{\partial {\mathit{U}}_e}{\partial \mathit{U}}\right)}^{\mathrm{T}}{\mathit{B}}_e^{\mathrm{T}}{\mathit{D}}_0^{\mathrm{T}}\left(\frac{\partial {\sigma }_e^{vm}}{\partial {\mathit{\sigma }}_e}\right) }$$

(32)

Considering the design independence of the external mechanical load, the sensitivity is expressed as:

$$\frac{\partial L}{\partial x_e}={\displaystyle \sum _{e=1}^{N_e}cp\frac{\partial {\sigma }^{PN}}{\partial {\sigma }_e^{vm}}{\left(\frac{\partial {\sigma }_e^{vm}}{\partial {\mathit{\sigma }}_e}\right)}^{\mathrm{T}}}\frac{d{E}_e^{\mathrm{II}}(x_e)}{d{x}_e}{E}_0({\mathit{D}}_0{\mathit{B}}_e{\mathit{U}}_e-{\mathit{D}}_{0}cp\mathit{\omega }\Delta T)-{\mathit{\lambda }}^{\mathrm{T}}\left(\frac{\partial \mathit{K}(\mathit{x})}{\partial x_e}\mathit{U}-\frac{\partial {\mathit{F}}^{th}(\mathit{x})}{\partial x_e}\right)$$

(33)

According to information from Equations (9) and (11), the sensitivity of the stiffness matrix and the thermal load vector for the design variable can be written as:

$$\frac{\partial \mathit{K}(\mathit{x})}{\partial x_e}={\displaystyle \sum _{e=1}^{N_e}\frac{d{E}_e^{\mathrm{I}}(x_e)}{d{x}_e}{E}_0{\displaystyle {\int }_{{\Omega }_e}{\mathit{B}}_e^{\mathrm{T}}}{\mathit{D}}_0{\mathit{B}}_{e}hd{\Omega }_e}$$

(34)

$$\frac{\partial {\mathit{F}}^{th}(\mathit{x})}{\partial x_e}={\displaystyle \sum _{e=1}^{N_e}\frac{d{\xi }_e(x_e)}{d{x}_e}{\xi }_0\Delta T}{\displaystyle {\int }_{{\Omega }_e}{\mathit{B}}_e^{\mathrm{T}}}{\mathit{D}}_0\mathit{\omega }d{\Omega }_e$$

(35)

Based on Equations (3) and (4), the sensitivity of the interpolated multi-material elastic modulus and thermal stress coefficient concerning the design variable can be obtained, respectively.

$$\frac{d{E}_e^i(x_e)}{d{x}_e}=\frac{E_j^N-E_{j+1}^N}{{({\rho }_j^N)}^{{\mu }_i}-{({\rho }_{j+1}^N)}^{{\mu }_i}} {\mu }_i{x}_e{}^{{\mu }_i-1} (i=\mathrm{I},\mathrm{II})$$

(36)

$$\frac{d{\xi }_e(x_e)}{d{x}_e}=\frac{{\xi }_j^N-{\xi }_{j+1}^N}{{({\rho }_j^N)}^{{\mu }_{\xi }}-{({\rho }_{j+1}^N)}^{{\mu }_{\xi }}}{\mu }_{\xi }{x}_e{}^{{\mu }_{\xi }-1}$$

(37)

The sensitivity of the total volume with respect to the design variable can be derived by:

$$\frac{\partial V}{\partial x_e}={\displaystyle \sum _{e=1}^{N_e}\frac{\partial x_e{v}_e}{\partial x_e}{v}_e}$$

(38)

The method of moving asymptotes (MMA) is used to solve the stress-constrained topology optimization for thermo-elastic structure problems. Due to the highly nonlinear behavior of the topology optimization problem [34], an external design variable shift limit of 0.05 is imposed on the MMA algorithm. The flowchart of the stress-constrained topology optimization method for thermo-elastic structures is illustrated in Figure 2.

5. Numerical Examples

This section demonstrates the feasibility of the stress-constrained topology optimization model for multi-material thermo-elastic structures through several 2D and 3D numerical examples. The initial values of design variables were 1, and the optimization was terminated if the relative variation of the objective function was less than 0.1%. The values of the penalty factors μ_I, μ_II and μ_ξ in the SIMP interpolation model were set as 3, 0.8 and 2, respectively. The value of parameter p for the p-norm function was set to 8. For multiple-materials design, the material C was considered as steel, with an elastic modulus of 2.1 GPa, Poisson’s ratio of 0.3, yield strength of 235 MPa (to consider the influence of different stress constraint values, it was set to 280 MPa in the Example 1), and thermal expansion coefficient of 1.21 × 10⁻⁵/°C. The three selected materials were normalized by the ordered SIMP method, and the normalization scheme with material density, elastic modulus, yield strength, and thermal stress coefficient was shown in

Table 1. Material parameter values for multi-material interpolation.

Material	ρN	EN	ξN	Color Setting
Void	0	0	0
A	0.5	0.3	0.5
B	0.7	0.5	0.7
C	1.0	1.0	1.0

.

5.1. 2D L-Shaped Beam Design

The first example considers the classical 2D L-shaped beam optimization problem, the most frequently considered test example in topology optimization. The geometric dimensions and boundary conditions are shown in Figure 3. It can be seen that the most prominent feature is the corner in the structure, which may lead to a certain stress concentration around the corner, which may be the problem that needs to be solved for stress-constrained topology optimization for thermo-elastic structures. The thickness of the design domain was set to 1 mm. The top left side of the L-shaped beam was restrained, and the upper right corner was subjected to a total load of F^m = 400 N. To avoid stress concentration near the load position, the concentrated load was equally distributed to the five adjacent nodes in the upper right corner. The structure was divided into 6400 four-node quadrilateral elements, and the maximum number of optimization iterations was set to 200.

Figure 4, Figure 5 and Figure 6 show the optimized topological configurations of the two materials (B and C), as well as the stress distribution for temperature variations of ΔT = 10 °C, ΔT = 50 °C, and ΔT = 90 °C, respectively. When ΔT = 10 °C, it can be seen from Figure 4 that the proportion of material B was higher, and there was a small distribution of material C near the corner of the L-shaped beam. With the condition of ΔT = 50 °C, the thermal load increases, and it can be seen that the proportion of material B increases compared with ΔT = 10 °C. There is a uniform distribution of relatively thin beams in the middle area of the structure, and there is also a distribution of material C in the middle part of these beams to accommodate the local stress increase. When ΔT = 90 °C, the value of the thermal load rises significantly, which is in the same order of magnitude as the applied external mechanical load, and thermal load accounts for the main effect. A high percentage of material C can be seen in Figure 6, to accommodate the higher global stress levels, and its stress distribution is more uniform. However, the stress values are slightly higher near the corners of the L-shaped beam than in other structural domains due to the external mechanical loads.

Considering the effect of different stress constraints, the topology optimization results for a stress constraint of 280 MPa are given in Figure 7. As the stress constraint value increased, the overall stiffness of the structure decreased, did the proportion of material C. The design volume also decreased from

Table 2. Topology optimization results for L-shaped beam with different values of ΔT.

Case	σs/MPa	ΔT/°C	Volume/%	Max Von Mises/MPa
1	235	10	26.53	234.86
2	235	50	30.67	234.90
	280	50	24.32	299.32
3	235	90	44.94	234.87

, relative to the stress constraint value of 235 MPa. The single-material thermo-elastic structure topology optimization results obtained by Deaton et al. [34] are shown in Figure 8. We can see that the multi-material optimization results can avoid stress concentrations more effectively and have a more uniform stress distribution by comparing Figure 5 and Figure 8.

The detailed evolution of the topology, together with the material and von Mises stress distribution for the case of ΔT = 50 °C, are shown in Figure 9 and Figure 10. It can be observed that the stress concentration appears initially at the corner of the L-shaped structure in the initial phase of the iteration. As the iterative process continues, by removing material at the smooth boundary of the L-shaped structure, the stress distribution of the topology gradually becomes uniform and free of stress concentration.

Figure 11 gives an iterative history of the volume ratio and maximum stress value for the L-shaped beam under three different temperature variations. The final volume fractions and the maximum von Mises stress values are listed in Table 2. As can be seen from Table 2, the target volume fraction of the optimized structure increases as the temperature change ΔT is enlarged, which clearly illustrates that the target structure requires more material to guarantee the effect of the higher thermal load. It can also be seen that the maximum stress value of the final topologies obtained from the minimized volume optimization model based on the stress constraint is always within the prescribed limit, and the final iterative process tends to be converged.

5.2. 2D T-Shaped Beam Design

The second example considers the T-shaped beam problem with the design domain and boundary conditions shown in Figure 12. The left and right sides of the structure are fully constrained, and the top-left point is subjected to a concentrated load $F_x^m$ = 500 N in the x-direction and $F_y^m$ = 400 N in the y-direction. The selected steel material, with a yield strength of 235 MPa, and the structure, were discretized into 8800 four-node quadrilateral elements.

Figure 13 and Figure 14 show the optimized topologies associated with three materials (A, B and C) at ΔT = 10 °C and ΔT = 60 °C, respectively, and it can be seen that the T-beam has an overall roughly triangular structure under the influence of external mechanical and thermal loads. The proportion of material A is higher at ΔT = 10 °C, and the proportion of material B is the second highest, mainly concentrated around material C. The proportion of material C is mainly concentrated in the right corner of the T-beam, because the horizontal force creates a bending moment at the corner, causing the stress value there to be higher. Due to the low thermal load at ΔT = 10 °C, the left side of the structure is only supported by a relatively thin beam, but the stress distribution shows that the maximum von Mises stress is always within the predefined limit and there is no stress concentration. When the value of ΔT increases to 60 °C, the thermal load is significantly enlarged and the shape of the optimized topology changes significantly; as can be seen, the number of beams increases, resulting in more triangular beams to guarantee the higher thermal load. The maximum stresses are still located at the right corner of the T-beam, and the percentage of material C is significantly higher. The part of the structure that is relatively less affected by mechanical loads is composed of a mixture of material A and material B, which is in line with our expected concept.

Figure 15 shows the iterative process with a fast convergence rate for the volume fraction and maximum von Mises stress at ΔT = 10 °C and ΔT = 60 °C, respectively.

Table 3. Topology optimization results for T-shaped beams with different values of ΔT.

Case	ΔT/°C	Volume/%	Max Von Mises/MPa
1	10	11.27	234.75
2	60	13.23	234.91

reveals the final volume fraction and maximum von Mises stress, and it can be seen that the volume fraction also increases significantly for increasing values of ΔT.

5.3. 3D T-Shaped Beam Design

In this section, a 3D L-shaped beam model with the mechanical and thermal loads and boundary conditions is shown in Figure 16. The initial geometrical dimensions are set as 60 mm × 60 mm × 8 mm. The top plane of the 3D L-shaped beam is fixed, and an external load with a magnitude of 140 N is applied vertically downward on its right side. The stress constraint is set to 235 MPa, and the maximum number of optimization iterations is set to 150.

The optimized topologies of the thermo-elastic structure with three materials (A, B and C) at ΔT = 5 °C and ΔT = 30 °C are shown in Figure 17 and Figure 18, respectively. For case ΔT = 5 °C, when the effect of mechanical load dominates, a higher proportion of material A and smaller proportions of material B and material C are presented. The stress distribution shows that the maximum von Mises stress is always within the constraint limit and there are no stress concentration phenomena. When the ΔT increases to 30 °C, the thermal load increases significantly. Meanwhile, the configuration of the optimized structure changes, and the proportions of the higher-strength materials B and C increase significantly. This follows the variation patterns of the two-dimensional L-shaped structure.

Figure 19 shows the iterative process for volume fraction and maximum von Mises stress at ΔT = 5 °C and ΔT = 30 °C, respectively.

Table 4. Topology optimization results for 3D L-shaped beams with different values of ΔT.

Case	ΔT/°C	Volume/%	Max Von Mises/MPa
1	5	19.24	234.32
2	30	31.43	234.51

demonstrates the final volume fraction and maximum von Mises stress. It can be seen that the minimum volume fraction increases significantly for increasing values of ΔT to make the structural strength meet the increased thermal load, whereas the maximum von Mises stress is maintained at the stress constraint limit.

6. Conclusions

This paper integrates the ordered-SIMP material interpolation model into the stress-constrained topology optimization for thermo-elastic structures. The coupled effect of mechanical and temperature loads on structural stresses is considered. A global stress constraint is adopted with the p-norm function, and the sensitivity of stresses is discussed using the Lagrange multiplier function and adjoint variable method. The proposed approach is discussed by three classical numerical examples, and the following conclusions can be drawn:

(1)

The maximum structural stress can be controlled by applying the global stress constraint to meet the material strength requirement. The stress concentration can be effectively avoided to achieve a uniform stress distribution. From the iterative curves of the objective function and the stress constraint, it can be seen that the proposed method can achieve stable convergence;

(2)

The optimized topologies under different temperature loads are significantly different. When the temperature load is small, the results are similar to the topology obtained without a temperature load and the proportion of high-performance material is low, whereas when the temperature load is high, the topologies change, and the proportion of high-performance material increases significantly to guarantee the higher thermal loads;

(3)

As the temperature loads increase, the volume ratio of the final topology increases, i.e., an increase in temperature load leads to the need for more material to satisfy the stress constraint. Future work will extend this approach to other optimization problems involving other multiple physical fields.