Reliability based topology optimization of thermoelastic structures using bi-directional evolutionary structural optimization method

The aim of this paper is to propose a novel computational technique of applying reliability-based design to thermoelastic structural topology optimization. Therefore, the optimization of thermoelastic structures' topology based on reliability-based design is considered by utilizing geometrical nonlinearity analysis. For purposes of introducing reliability-based optimization, the volume fraction parameter is viewed as a random variable with a normal distribution having a mean value and standard deviation. The Monte Carlo simulation approach for probabilistic designs is used to calculate the reliability index, which is used as a constraint related to the volume fraction constraint of the deterministic problem. A new bi-directional evolutionary structural optimization scheme is developed, in which a geometrically nonlinear thermoelastic model is applied in the sensitivity analysis. The impact of changing the constraint of a defined volume of the required design in deterministic problems is examined. Additionally, the impact of altering the reliability index in probabilistic problems is investigated. The effectiveness of the suggested approach is shown using a benchmark problem. Additionally, this research takes into account probabilistic thermoelastic topology optimization for a 2D L-shaped beam problem.


Introduction
Topology optimization (TO) has become a dynamic discipline, thus attracting researchers from different fields: mechanical, civil, and aerospace fields. It is supposed to be one of the most rewarding methods (Li et al. 2020;Zhou and Rozvany 1991;Zhu et al. 2016). Depending on the definition of the design variables, TO may generally be divided into two groups. The density approach falls under the first group. It connects each design variable to a finite element and uses that connection to determine if a solid material is present (1) or not (0). For such category, various approaches have been developed such as evolutionary structural optimization (ESO) method (Xie and Steven 1993) and solid isotropic material with penalization (SIMP) (Bendsøe 1989). While the second category involves design variables related to the structural boundaries such as the phase field method (Wang and Zhou 2004) and level set method (Sethian 1999).
When structures undergo significant deformations, the tangential stiffness matrix may no longer be positive definite. As a result, various studies have arisen in the past two decades that take geometrically nonlinear topology into consideration (Abdi et al. 2018;Buhl et al. 2000).
Examining contemporary researches in the field of topology optimization that considers geometrically nonlinear, a topology optimization code was developed by (Chen et al. 2019) where an additive hyperelasticity approach is utilized to solve topology optimization of structures going through large deformations. (Zhu et al. 2021) presented topology optimization algorithm by adopting SIMP topology optimization method. By converting the strain energy between small and large deformation theories, an energy interpolation technique was proposed by . Using a nonlinear FE formulation, (Buhl et al. 2000) addressed the topology optimization of geometrically nonlinear problems by coupling SIMP with it. In many circumstances, the solutions from the nonlinear modeling are just slightly different from the linear ones, as shown by the examples supplied by Buhl and his colleagues. Nonetheless, a big disparity may emerge, if snap-through effects are at play in the issues at hand.
Considering reliability designs might be regarded as an important strategy since the engineer should deal with the existence of uncertainties in various cases including differences in dimensions and material properties. Thus, structural topology optimization is more applicable (Chun et al. 2016;Logo 2007;Zhao and Wang 2014). (Jung and Cho 2004) proposed reliability based design optimization technique by considering probabilistic constraints related to displacement using the performance measure approach. Based on probabilistic methods, a robust topology optimization technique was proposed to solve problems considering uncertainties by (Meng et al. 2021). In the study of (Lógó 2012), probabilistic topology optimization problem of continuum structure was proposed by considering random location of the applied loads. (Dunning et al. 2011) considered uncertainties in performing topology optimization by randomly assuming the value and the direction of the applied forces. A reliability based optimization technique was formed to improve the material distribution in the presence of microstructure uncertainties in the study of (Gao and Liu 2021). Reliability based design optimization is also studied by (Ghasemi et al. 2015b), who present a probabilistic computational optimization technique for internal cooling channels in Ceramic Matrix Composite (CMC) subject to thermal and mechanical loadings. Moreover, (Ghasemi et al. 2015a) proposed an optimization procedure that uses two stages and a sequential approach to discover the best fiber quantity and distribution in solid composites while accounting for uncertainty in the design parameters.
Most prior structural optimization research focused on externally loaded systems. Nonetheless, there is still opportunity to get better concerning instabilities associated with thermoelastic problems of structures. Considering thermoelastic models in topology optimization can sufficiently be useful practice to face the challenges of complex design. (Rodrigues and Fernandes 1995) pioneered thermoelastic structural topology optimization by using a homogenization technique to meet the aim of structure compliance reduction under combined temperature and mechanical loads. An algorithm for optimizing the topology of structures exposed to thermoelastic and mechanical loads at the same time was described by (Deaton and Grandhi 2016). Considering multiple materials structure, (Gao et al. 2016) suggested a topology optimization technique taking steady-state temperature and mechanical loads into consideration. An evolutionary technique considering thermoelastic problems was developed by (Li et al. 1999) for minimization of displacement. The buckling performance of infill structure was considered to propose topological optimization algorithm by (Gan and Wang 2022), in which the authors made use of the SIMP technique to study the thermoelastic coupling effect.
Reliability-based design according to the assumption of volume fraction as a random variable considering linear and geometrical nonlinearities as well as elastic and elastoplastic models were adopted by BESO method recently by (Movahedi Rad et al. 2021). Consequently, the current paper is a continuation of prior research to propose a novel computational algorithm by applying geometrically nonlinear reliability-based design to thermoelastic structural topology optimization in which BESO method is developed to fulfill the aim of this research.
The remaining portions of the arrangement of this publication are coordinated as: The theoretical context of the problem is discussed in Sect. 2. The technique of the enhanced BESO approach is discussed in Sect. 3. This study's numerical examples are introduced in Sect. 4. In Sect. 5 the work is finally summarized.

Finite element analysis of geometric nonlinearities
Nonlinear Lagrangian FE model is utilized for performing the analysis of large displacements: where u stands for point-wise displacement, and i, j and k indicate the coordinate axes.
where the displacement vector is represented by U, and B is the matrix that transforms the change in displacement dU into the change in strain. Equation 3 is used to write the Hooke's law for materials of intermediate densities: In which the second stress defined by Piola-Kirchhoff is denoted by s ij , the power of penalization is given as p , while the Péclet number is represented by p e and strains proposed by Green-Lagrange are given as kl , also, solid isotropic material produces a constitutive tensor denoted by C 0 ijkl . Consequently, the residual may be conceptualized as the deviation from the achieved equilibrium by following: where the externally applied load is represented by P , and Piola-Kirchhoff stress in this equation is denoted by the vector s . Taking into account that when the residual equal 0, the state of equilibrium is found. Generally, Equation (4) can be iteratively solved by adopting Newton-Raphson iterative technique.
where K T indicates the tangential stiffness matrix.

Implementing thermoelastic analysis in topology optimization
According to earlier research (Huang and Xie 2009;Querin et al. 1998), the topology optimization adopting BESO approach may be described as: where C represents the compliance of the considered structure and K represents the tensor of global stiffness. N is the total number of elements allowed in the design space, V * is the volume of the whole structure and V i stands for element's volume. Taking into consideration that the binary design variable x i is (1) or (0), referring to material and void, respectively. And V f is the volume fraction. Besides, f and u are loads and displacements vector, respectively. It should be mentioned that both mechanical loading f m and thermomechanical loading f heat are applied and combined in the loading vector as f = f m + f heat . Consequently, the displacements measured and shown by the objective function already account for thermal expansion. The semi-coupled thermoelasticity hypothesis is utilized. The temperature field is obtained by first solving the heat regulating equations. Next, the forces on the body caused by the temperature field are added to the other applied forces to determine the whole response of the elastic body (Ghasemi et al. 2015b). According to the Bubnov-Galerkin weak form: where Ω denotes the domain, Γ represents the boundary condition, stands for the temperature field, b is the body force, t Γ is the traction on the borders, represents the arising strain, the variational operator, C is the structural compliance, and u is displacements vector.
The w -function is a test chosen from a set of functions for B-splines. They are also used to approximatively model the temperature and displacement fields.
where u represents the vector of nodal displacements and represents temperatures, and the shape functions are represented by N , , are the knots, and p, q are the degrees of B-spline curve.
Both the strain-displacement and heat flux-temperature gradient expressions may be constructed as: and The elastic and thermal problems are represented by the matrices B e and B heat , which store the derivatives of the shape functions N.
The matrix form of the discretized system of equations is obtained by substituting the B-spline approximation function into Bubnov-Galerkin weak form is expressed as: The heat force vector, f heat , and the local conduction matrix, Kc , are calculated in the following ways: The heat conduction matrix, H , is denoted by the superscript T , which denotes transpose. Equation (16)'s first and second integrals, represent heat conduction (in volume Ω) and convection (on surface Γ3), respectively.
Finally, the matrix form of the discretized linear system of equations for the thermoelasticity issue is obtained by plugging in the test function and its derivatives into Eq. (11): where K is the global stiffness matrix of the elastic problem.
The BESO approach is based on concurrently adding and eliminating parts according to its sensitivities, resulting in the optimal design of the structure. To consider heat conduction (Xia et al. 2018) in the optimization problem, this is how the sensitivity number is calculated: The e − th element's design variable X e might be absent (X e = min ) or present X e = 1 depending on the design variable. The sensitivity number, denoted by the symbol k e , of the e − th element in the k − th optimization loop is determined by:  (19): From Eqs 21 and 22, we obtain: For purpose of avoiding mesh dependency and checkerboard patterns, a spatial filtering technique is applied Xie 2007, 2010) and defined as: where w ej is the weight which is defined as: where d(e, j) represents the Euclidean distance between centers of the e -th and j th elements and the filter radius is denoted by r min .
In order to figure out how the final results will be obtained in BESO method, the target volume for the next iteration (Vk + 1) must be provided before any changes are made to the existing design (such as the removal or addition of elements). The target volume in each iteration may drop or rise until the constraint volume is attained, depending on whether the volume constraint ( V * ) is larger or lower than the volume of the original estimate design. When expressing the change in volume, we may use the expression: where k is the current iteration number and ER represents the evolutionary ratio.
Once the volume limit has been met, the following formula guarantees that the volume of the structure will not change in subsequent iterations: Then, the elements are arranged from most sensitive to least sensitive.
In the case of solid components, deletion will occur if: In the case of void elements, it is added if: where th del and th add are the threshold sensitivity numbers of erasing and adding elements, respectively.
In the (BESO) approach, optimization is carried out in an iterative process by repeatedly adding and removing components until the convergence condition is satisfied. According to this research, the following criteria for convergence was used: The objective function is denoted by F , the tolerance for convergence is denoted by , the current iteration is denoted by k , and the integer number N that yields constant compliance for at least ten iterations.

Introducing reliability based design into deterministic optimization problem
In this study, for purpose of calculating reliability index ( ) which is based on probability of failure ( P f ) values, Monte-Carlo sampling method is considered. The idea of Monte-Carlo technique generally depends on the concept of producing realizations x according to X which is random vector of PDF-f X (x) . Therefore, to calculate P f , we compare the number of points generated to the total number of points within the failure domain (Stanton et al. 2000). By adding a D f , indicator function, this idea may be expressed as: To get the mean and standard deviation of Df (X) , compute as follows: When calculating P f , an estimate of the mean value is adopted as follows: where X (z) are set of independent random vectors having f X (x) taking into account that ( z = 1, … , Z).
It is important to note that in probabilistic models, V f is treated as a random variable. As a result, we can calculate its mean value as well as its standard deviation. In addition, its mean is and its standard deviation is ar . It should be noted that Gaussian distribution (Normal distribution) model is utilized in this research due to its simplicity, so the entire distribution can be indicated by only specifying two parameters: mean and variance. To get the estimator's mean and standard deviation, we do the following: The reliability constraint then can be presented by expressing as: The method will end when the calculated reliability index for each iteration, denoted by calc , achieves the desired value, denoted by target .
For determining target and calc , the following terms are used: for where (1) stands as a sign for the inverse of the normal distribution Φ.
Accordingly, the problem of reliability based optimization is constructed in the presence of reliability constraint as: Here, Eqs. (42, 45 and 46) serve the same purpose as Eqs. 6, 7, 8 and 9. While Equation 47 displays the new circumstances associated with the volume fraction reliability limit.

The expanded BESO method
After a short mathematical explanation of the problem, the approach for the reliability based thermoelastic optimization problem with the given constraints adopting BESO method is provided in Fig. 1, the technique is summed up as follows: 1. Identifying the design domain. 2. FEA execution followed by sensitivity calculation. 3. Averaging the sensitivity values. 4. Setting the value of targeted volume for the next iteration of the process. 5. The subtraction and addition of constituent elements. 6. repeat the procedures from the steps (2 to 5) until all of the conditions have been met.
and convergence criteria.

Numerical examples
In this part, we will look at two different numerical examples. The optimization of a linear thermoelastic topology using a probabilistic approach is used as a starting point while the second example is considered for probabilistic geometrically nonlinear thermoelastic topology optimization. Thus, the second example is considered to show the differences between the results of the linear and geometrically nonlinear thermoelastic optimum designs. A rectangular plate clamped from both sides is counted as a first example, and the results are evaluated against a (45) target − calc ≤ 0 Reliability based topology optimization of thermoelastic structures using bi-directional… set of benchmark problems which were done by (Li et al. 1999;Rodrigues and Fernandes 1995) in order to establish its viability. This research also includes a second numerical example, this time based on the L-shaped beam and the geometrically nonlinear thermoelastic reliability it entails. For reliability evaluation, Monte-Carlo simulation is used with the assumption that V f is random to simulate the probabilistic nature.

Linear problem: rectangular plate clamped from both sides
Initially, an optimization problem for a rectangular plate that is clamped on both sides is investigated. Figure 2 demonstrates this optimization problem where the considered applied loads are (mechanical load) which is F = 10 kN and three different (temperature loads) which are ΔT = 0 o C, 5 o C and 10 o C . The design is 10mm in thickness. The chosen material also has a thermal expansion coefficient of 12 × 10 −6 , Poisson's ratio of 0.3 , and Young's modulus is assumed 210000MPa . The considered parameters of BESO are ER = 2% , AR max = 1% , r min = 30mm and = 1% . Considering that V f has probabilistic nature with mean value = 4 0% and variance = 5% . The assumed sample points for probabilistic Monte-Carlo simulations are ( Z = 1.0 × 10 8 ).
As previously mentioned, to approve the proficiency of the proposed work, the findings of this example are assessed with a benchmark problems (Li et al. 1999;Rodrigues and Fernandes 1995). Table 1 represents the obtained optimal solution according to two different values of V f . Furthermore, for probabilistic design, Table 2 shows the results of two different values of target . Various topologies are achieved for deterministic and probabilistic designs when different thermal loads ΔT are considered, as seen by the layouts produced. In addition, by considering each case of the applied thermal loads ΔT , it can be noticed that when V f is changed, the resulted optimized shape will change too for deterministic case. Also, a significant difference between the resulted layouts for each thermal load ΔT considering different target by adopting probabilistic technique. On the other hand, inserting a reliability constraint into a deterministic design clearly affects the results, given that the results of deterministic designs vary from those of probabilistic designs.
According to the obtained optimal topologies, we can say that the filter scheme which discussed earlier in Sect. 2.2 functions as a low-pass filter, removing non-essential elements of the structure below a certain length scale. The optimal topology will no longer be sensitive to the mesh sizes, which is a major advantage of using this filter method.
All models' displacements are examined based on complementary work. Thus, another comparison according to the resulted complementary work for deterministic and probabilistic designs are considered. Figure 3 helps to illustrate that as the applied thermal load increases the complementary work increases too. For instance, by considering V f = 45% , the complementary work is increased by 59.05% from 135.74kJ in case of ΔT = 0 • C to 331.52kJ in case of ΔT = 10 • C . Also, by considering V f = 40% , the complementary work is raised by 52.16% from 143.27kJ in case of ΔT = 0 • C to 299.51kJ in case of ΔT = 10 • C.
Complementary work values for probabilistic design as a result are represented in Fig. 4 in which we can say that the complementary work increases as the applied thermal load increases for each value of target . By considering target = 5.36 , the complementary work is increased by 62.03% from 128.88kJ in case of ΔT = 0 • to 339.46kJ in case of ΔT = 10 o C . Besides, by considering target = 3.13 , the complementary work is increased by 59.88% from 133.18kJ in case of ΔT = 0 • C to 331.87kJ in case of SSSS-CXSΔT = 10 • C.

Geometrically nonlinear problem: L-shaped beam
The second topology optimization problem is geometrically nonlinear optimization of 2D L-shaped beam to learn more about how the suggested formulation might be used to the design of nonlinear thermoelastic structures. The design space that has to be optimized is shown in Fig. 5. The considered applied loads are (mechanical load) which is F = 12 kN and four different (temperature loads) which are ΔT = 0 • C,5 • C,10 • C and15 • C . The design is 10mm in thickness. Both 45% and 40% of the design domain are accounted for by the values of V f . Furthermore, the chosen material has a thermal expansion coefficient of 12 × 10 −6 , a Poisson's ratio of 0.3 , and a Young's modulus of E = 70000MPa . The following criteria are taken into account by BESO, ER = 1% , AR max = 1% , r min = 18mm and = 0.1% . Considering that V f has probabilistic nature with mean value = 4 0% and variance = 5% .  Similar to what we have done in the previous problem, a comparison between the obtained layouts of deterministic and probabilistic linear designs is considered in this example too. Besides, a new comparison according to the obtained layouts Tables 3 and 4 indicate the deterministic results of the obtained layouts according to two different values of V f for linear and geometrically nonlinear designs, respectively. When applying the nonlinear topology design technique, it is clear that completely different solutions are achieved for each thermal load case than with the linear design. It has also been discovered that the optimal designs are able to regulate the thermal effect to the point that the provided temperature change can counterbalance the mechanical load, even to the point of preventing buckling and nonlinear snap-through behavior. This result stands in sharp contrast to linear elastic optimization, where the compliance is temperature-dependent. It can be concluded from these optimized shapes that the optimized structures for greater thermoelastic loads contain material layouts with slits and struts, which enhance the buckling capacity in areas where the design is vulnerable to buckling. Moreover, we can say here also that the filter scheme approach which discussed earlier in Sect. 2.2 functions perfectly to make the optimum topologies no longer sensitive to the mesh sizes.
While Tables 5 and 6 indicate the probabilistic results of the obtained layouts according to two different values of target . Again, here we can say that introducing reliability constraint into deterministic design will change the optimized layouts of the structures. In deterministic design, for each specified applied thermal load, the resulted optimized shapes considering geometric nonlinearities are not the same as which are obtained in linear design based on V f value. Also, the resulted optimized shapes indicate that considering different target will affect the final optimized layout for each case of the applied thermal load according to probabilistic design.
Optimized structures for greater thermoelastic loads contain material layouts with slits and struts, which enhance the buckling capacity in areas where the design is vulnerable to buckling.
All models' displacements are examined based on complementary work here also. Thus, another comparison based on the resulted complementary work for deterministic and probabilistic designs are considered. Figure 6 represents the resulted values of complementary work considering deterministic linear and geometric nonlinearity designs. The nonlinear designs are clearly superior to the linear designs in terms of performance under the load conditions under which they were developed, as seen by the smaller magnitudes of the complementary work required by each. This highlights the need for a nonlinear topology optimization technique for large-displacement situations like those incorporating snap-through effects, where the complementary work of the linear design combining buckling and snap-through effects is much greater than for the nonlinear one.
Considering linear and geometric nonlinear designs, for each value of ΔT , as V f decreases the complementary work increases. Also, the obtained complementary work values considering linear design are greater than considering geometrically nonlinear design if we compare it according to the same value of V f and ΔT . Furthermore, by considering linear and geometric nonlinear designs, for each value of V f , as the applied thermal load increases the complementary work increases too. For instance, by considering V f = 45% and linear design, the complementary work is increased by 1.37% from 25.92kJ considering ΔT = 0 • C to 26.28kJ considering ΔT = 15 • C . Also, considering V f = 45% and a geometrically nonlinear design increases the complementary work by 3.2% from 25.41kJ when ΔT = 15 • C is considered to 26.25kJ when ΔT = 15 • C is considered. Figure 7 represents the resulted values of complementary work considering probabilistic linear and geometric nonlinearity designs. Considering linear and geometric nonlinear designs, for each value of ΔT , as target, decreases the complementary work increases. In addition, the obtained complementary work values considering linear design are greater than considering geometrically nonlinear design if we compare it according to the same value of target and ΔT.
Furthermore, by considering linear and geometric nonlinear designs, for each value of target , as the applied thermal load increases the complementary work increases too. For instance, by considering target = 5.61 and linear design, the complementary work is increased by 1.62% from 23.1kJ considering ΔT = 15 • C to 23.49kJ considering ΔT = 15 • C . In addition, by considering target = 5.61 and geometrically nonlinear design, the complementary work increases by 1.84% from 23kJ when ΔT = 15 • C is considered to 23.43kJ when ΔT = 15 • C is considered.

Conclusions
This study aimed to optimize thermoelastic models with geometric nonlinearity and randomness using extended BESO to achieve reliability based design. 3.23 Table 6 The optimal layouts for a probabilistic geometrically nonlinear design Geometrically nonlinear designprobabilistic target Optimized shape 5.61

3.23
Also, to examine how thermoelastic loads affect a structure's material distribution.
Throughout the optimization process, for incorporating reliability based design, the V f was treated random. In other words, by applying the reliability theory, the proposed method was developed in order to get the optimal solution by satisfying the reliability constraint. And P f were found by adopting Monte-Carlo approach, and therefore the values. In optimization problems, the design-dependent thermoelastic stress loads are examined. The suggested approach for optimizing the topology of thermoelastic structures may be regarded as an effective tool for designers who encounter considerable difficulties while minimizing the mean compliance of thermoelastic structures. The findings of the examined benchmark illustrate the worth and usability of the offered methodology. In addition, the findings of the examined problems reveal the link between reliabilitybound and objective function.
The following main points can conclude the proposed work: 1. When introducing reliability constraint into deterministic design, the optimal layouts of the Geometrically nonlinear probabilistic designs 0℃ 5℃ 10℃ 15℃ 5. As it is shown in the obtained results, by considering a specified value of ΔT , there is negative relation between the complementary work and V f for deterministic designs of both models (linear models and models with geometric nonlinearity). Besides, by considering a specified value of ΔT , there is negative relation between the complementary work and target for probabilistic designs of both models (linear models and models with geometric nonlinearity).
The proposed work settles the major concerns that emerge by dealing with thermal loads. And therefore, it is obvious to see how this is a major step forward towards more reasonable framework for thermoelastic topological designs of geometrically nonlinear problems under reliability constraint. Future work is anticipated to address other topology optimization problems by defining extra constraints within problem formulation.
Funding Open access funding provided by Széchenyi István University (SZE). The authors did not receive support from any organization for the submitted work.

Data availability
The datasets generated during and/or analyzed during the current study are available in the main manuscript.

Conflicts of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.