Topology optimization of turbine disk considering maximum stress prediction and constraints

Cheng YAN, Ce LIU, Hn DU, Cunfu WANG,*, Zeyong YIN,c

a School of Aerospace Engineering ,Xiamen University, Xiamen 361102, China

b National Key Laboratory of Science and Technology on Helicopter Transmission, AECC Hunan Aviation Powerplant Research Institute, Zhuzhou 412002, China

c Aero Engine Academy of China, Beijing 101304, China

KEYWORDS Centrifugal loads;Maximum stress prediction;Stress constraints;Topology optimization;Turbine disk

Abstract For the stress-constrained topology optimization of a turbine disk under centrifugal loads, the jagged boundaries of the mesh and the gray densities on the solid/void interfaces could make the calculated stress field inconsistent with the actual value.It may result in overestimating the maximum stress and thus affect the effectiveness of stress constraints.This paper proposes a new method for predicting the maximum stress to overcome the difficulty.In the process, a predicted density is newly defined to obtain stable boundaries with thin layers of gray elements, a transition factor is innovatively proposed to evaluate the effects of intermediate-density elements, two different stiffness penalty schemes are flexibly used to calculate the elastic modulus of elements, and a linear stress penalty is further adopted to relax the stress field of the structure.The proposed approach for predicting the maximum stress value is verified by the analysis of a structure with smooth boundaries and the topology optimization of a turbine disk.An updating scheme of the stress constraint in the topology optimization is also developed using the predicted maximum stress.Some key ingredients affecting the optimization results are discussed in detail.The results prove the effectiveness and efficacy of the proposed maximum stress prediction and developed stress constraint methods.

1.Introduction

As a powerful conceptual design method, topology optimization has been widely used in industrial designs.1Topology optimization can find the optimal material distribution in a given design domain to maximize or minimize the given objective function under assigned constraints.In the structural design process, stress is an essential indicator since it determines whether the loads exceed the bearing capacity of the material.Therefore, topology optimization design considering stress constraints has received continuous attention and extensive research, and it is still an essential research direction.Topology optimization problems considering stress constraints under fixed loads have been well studied.However,the designdependent loads, such as the centrifugal loads, are essential factors affecting the structural stress field for rotating machinery.Topology optimization for such problems considering stress constraints is still challenging.As a typically rotating part, the turbine disk is primarily subjected to centrifugal loads.Topology optimization of the turbine disk is significant for lightweight design.2This paper mainly studies the stressconstrained topology optimization of a turbine disk under centrifugal loads.

Structural optimization has been widely used to improve the performance of the turbine disk.Shape optimization of the turbine disk is relatively mature.Bhavikatti and Ramakrishnan3used the nonlinear programming method to optimize the shape of the rotating disk.Hattori et al.4used the sequential linear programming method to design a turbine rotor under stress limit constraints,considering centrifugal and thermal loads.Cheu5considered geometric and stress constraints and used the feasible direction method and sequential linear programming to optimize the shape of turbine disks.O¨ zakc?a et al.6introduced a robust algorithm for optimizing turbine disks based on the automated Structural Shape Optimization(SSO).Liu et al.7developed a new method named Metamorphic Development to optimize the profile of a turbine disk and used the axisymmetric element to display the structure.Meng et al.8considered the disk’s centrifugal loads and developed a finite cell method to optimize the shape of an axisymmetric disk.However, there are few publications on topology optimization of turbine disks considering design-dependent loads.Lopes et al.9took the optimality criteria to optimize rotors under self-weight and inertial forces, and the original domain is a three-dimensional structure.Shen et al.10optimized the turbine disk using the Evolutionary Structural Optimization(ESO)method.Zheng and Gea11selected compliance and kinetic energy as objective functions and optimized the turbine disk.Wang and Yu12optimized the 1/4 model of the turbine disc and obtained a hollow turbine disk.Liu13utilized the Bi-directional Evolutionary Structural Optimization(BESO) method to optimize the turbine disk considering centrifugal and thermal loads.Wang et al.14considered the influence of the temperature field on the structure and made a lightweight design of the turbine disk.Until now, the density method is less studied for the topology optimization of the turbine disk.Considering that the density method theory and algorithm are mature and widely used in designing complex industrial products, the topology optimization of the turbine disk based on the density method still has essential research significance.

However, solving topology optimization problems with stress constraints faces many challenges in the density method.15–17The first challenge is the singularity phenomenon in the optimization process.Sved and Ginos18first observed the singularity phenomenon when bars could not be removed in the truss optimization.In the density-based topology optimization, the local stress is discontinuous when the density approaches zero.Thus, the optimization space is discontinuous in the stress-constrained topology optimization problem.The degenerate space would hinder the gradient-based optimization algorithm from converging to the optimal value.In order to make the optimization space continuous, the stress relaxation method has been widely used.The first relaxation method is the ε-relaxation method.19Based on the εrelaxation method, the qp approach and stress interpolation schemes were proposed.20,21The second challenge is a large number of stress constraints.The nature of stress constraints is local, and each element corresponds to a stress constraint.Therefore,it is difficult and time-consuming to solve the sensitivity by either the direct or adjoint methods.In order to reduce the number of stress constraints, aggregation methods,such as the p-norm, Kreiselmeier-Steinhauser (KS) functions and the Heaviside projection integral22–24,have been proposed to establish a global stress constraint.It could significantly improve the calculation efficiency.However, the stress aggregation function reduces the ability to control local stress.In order to make a trade-off between global and local stressconstrained methods, block and cluster aggregation methods are proposed to improve the control performance of local stress.25–26The third challenge is the high nonlinearity of stress.The stress of the element is susceptible to the change of density, and the severe stress concentration often leads to a large stress gradient.27To improve the convergence stability,Le et al.21adopted the density filtering method.

In addition to the three challenges mentioned above, the fourth challenge in stress-constrained topology optimization is the calculation accuracy of the stress field.In densitybased topology optimization, the accurate finite element analysis of the structure is critical, especially for calculating the stress field.28However, due to the existence of jagged boundaries and blurred regions (namely gray densities), it is difficult to calculate the stress field accurately through the fixed mesh.Especially when the design-dependent loads exist, the intermediate-density elements significantly influence the accurate calculation of the stress field.Usually, the structure needs to be post-processed to eliminate the jagged boundaries and blurred regions after topology optimization.Then the finite element analysis with body-fitted mesh would be conducted for the optimized structure with smooth boundaries.Unfortunately,the maximum stress calculated by the body-fitted mesh is usually inconsistent with that calculated by the fixed mesh in topology optimization and may even exceed the bearing capacity of the material.It will affect the effectiveness of stress constraints in topology optimization and even lead to the failure of optimization and the unavailability of results.To avoid these difficulties, the stress calculation on the jagged boundaries and the maximum stress calculation of the structure have been studied.28–31Sva¨rd28used the Internal Value Extrapolation(IVE)method to calculate the boundary stress more accurately.De Troya and Tortorelli29refined the mesh in the highstress concentration area yet increased the calculation cost.Silva et al.30improved the calculation accuracy of stress on the jagging boundaries by selecting appropriate Heaviside parameters and establishing an intermediate material layer between solid and void.Amir31used a sharper projection method to calculate the maximum stress of the structure and reduced the error with the actual maximum stress.However,all these studies are stress-constrained topology optimization problems under fixed loads.

For the stress-constrained topology optimization under the design-dependent loads, the jagging boundaries of the mesh and the gray densities on the solid/void interface would significantly affect the stress field distribution of the structure,which is different from that under the fixed loads.To our knowledge,the accuracy of using the density method to calculate the stress field under the design-dependent loads has not been investigated.In this paper,the p-norm based aggregation is first used to approximate the maximum stress, and then normalized by the actual maximum stress to construct the global stress constraint.Hence, the calculation of the actual maximum stress of the structure is directly related to the topology optimization results.We develop a method to predict the actual maximum stress in this paper.The Heaviside function is used to map the filtered density to obtain the newly defined predicted density, enabling stable boundaries with thin layers of gray elements and a more stable and accurate maximum stress.A transition factor is innovatively proposed to help the flexible selection of the linear or nonlinear stiffness penalty schemes to calculate the elastic modulus of elements.The linear stress penalty is further employed to relax the structure’s stress field and calculate the predicted maximum stress.The method is verified both in a structure with smooth boundaries and in the topology optimization of a turbine disk.The Rational Approximation of Material Properties (RAMP) interpolation is used in both stiffness and stress penalization for the structural stress field calculation and the sensitivity analysis.A practical updating scheme of the stress constraint in topology optimization is developed using the predicted maximum stress.Parameters that might affect the optimization results are thoroughly discussed.An optimized turbine disk is reconstructed and analyzed using the ANSYS software to further verify the effectiveness of the proposed maximum stress prediction method and the developed stress constraint method.Finally,a comprehensive analysis of the different topological configurations of the turbine disk is conducted in order to provide more support for engineering design.

The structure of the remainder of this paper is given as follows.The studied turbine disk model and the axisymmetric finite element analysis method are given in Section 2.The proposed maximum stress prediction method is introduced in Section 3.The formulations of two kinds of optimization problems and the corresponding stress measure, stress constraint relaxation method, and sensitivity analysis are presented in Section 4.The example of turbine disk topology optimization is discussed in Section 5.The summary and conclusion of this paper are given in Section 6.

2.Turbine disk model and axisymmetric structure analysis

2.1.Turbine disk model

The turbine disk is a crucial load-bearing and connecting component in an aircraft engine.The actual geometry of the turbine disk is complex, which is challenging to be directly used for topology optimization.For the conceptual design of the turbine disk, the actual geometric structure should be reasonably simplified.The mounting edges and turbine blades could be removed,while the disk body should be retained in the simplification process.In this paper, the topology optimization design of the disk body is carried out.The simplified threedimensional model of the turbine disk is shown in Fig.1.

Fig.1 Simplified three-dimensional model of a turbine disk.

When the aero-engine is working,the turbine disk rotates at high speeds in the high-temperature environment.Accordingly, the loads on the turbine disk mainly come from three aspects.On the one hand, the high-speed rotation of the turbine disk itself can lead to the design-dependent centrifugal loads.On the other hand, the fixed tensile loads on the outer radius of the turbine disk are caused by the high-speed rotation of turbine blades.Additionally, the thermal loads are caused by the non-uniform temperature field on the turbine disk.The tensile loads are design-independent loads, which remain unchanged during the topology optimization of the turbine disk.The thermal and centrifugal loads belong to designdependent loads.The change in disk structure will cause the change of loads, showing the characteristics of ‘‘no load without material”.It brings great difficulties to the topology optimization of the turbine disk.The introduction of thermal loads mainly faces the challenges of the real-time updating of convective heat transfer boundary conditions, the accurate calculation of the temperature field, and the stable convergence of the stress-constrained topology optimization.The introduction of centrifugal loads is mainly faced with the following challenges: the non-monotonic behavior of the compliance, the parasitic effects for low densities, the accurate calculation of maximum stress, and the stable convergence of the stressconstrained topology optimization.As far as we know, no research has solved the above problems simultaneously.Considering that the centrifugal and tensile loads are more dominant than the thermal loads.Especially for the low-pressure turbine disk,the thermal loads are obviously less than the centrifugal and tensile loads.In the conceptual design stage of the turbine disk, topology optimization is usually utilized to obtain the initial disk configuration without considering the thermal loads.32,33In the detailed design stages, the shape and size of the turbine disk will be re-optimized by considering more comprehensive thermal loads and even fluid-thermalstructure coupling analysis.34,35The primary purpose of this paper is to solve the difficulties caused by centrifugal loads.Some novel topological configurations will be obtained,which is a critical consideration for the turbine disk’s design and is meaningful for engineering.Therefore, only the centrifugal and fixed tensile loads are considered in this paper, and the temperature field of the turbine disk is set to be uniform.The material properties at a constant temperature are used.

The centrifugal loads are symmetrical to the central rotation axis.The tensile loads are approximately evenly distributed on the disk rim surface.Therefore, the tensile loads are also considered to be symmetrical to the central rotation axis.The turbine disk is constrained only by axial displacement.Thus,the studied turbine disk is an axisymmetric model.The original three-dimensional model can be further simplified to a two-dimensional axisymmetric model, which can significantly save computing resources and improve the efficiency of optimization design.Considering that topology optimization is a conceptual design, the existing design domain can be expanded to some extent.Fig.2 shows a schematic diagram of the axisymmetric model of the expanded turbine disk.The cylindrical coordinate system is used to describe the structure,where r represents the radial direction and z represents the axial direction.The shaded areas are the non-design domains,and the blank area is the design domain.The geometry of the turbine disk is defined using five parameters,which are R1,R2,W1, W2, and T.In detail, R1is the inner radius of the turbine disk;R2is the outer radius of the turbine disk;W1is the inner width of the turbine disk; W2is the outer width of the turbine disk;T is the thickness of the non-design domains.F is the tensile load.

2.2.Finite element analysis of axisymmetric structure

This paper uses the rectangular axisymmetric element with four nodes to discretize the expanded axisymmetric turbine disk.The axisymmetric element is similar to the plane element but has particular properties.36Each axisymmetric element essentially represents an annular structure, as shown in Fig.3.Accordingly, the circumferential strain and stress must be considered.

Because the model is axisymmetric, the circumferential displacement of each node of the axisymmetric element is zero.So, each node has two degrees of freedom, and the displacement vector unof each node could be expressed as

Fig.2 Schematic diagram of the axisymmetric model of expanded turbine disk.

By combining the node displacement vector ueof the element and the shape function matrix N(r,z), the element displacement field uf（r,z ) could be uniquely expressed by the polynomial.

where B is the strain matrix,εris the radial normal strain,εθis the circumferential normal strain,εzis the axial normal strain,γrzis the shear strain between the r direction and the z direction.

The selected material in this paper is isotropic.The elastic matrix D of the axisymmetric element could be written as

where σris the radial normal stress, σθis the circumferential normal stress, σzis the axial normal stress, σrzis the shear stress between the r direction and the z direction.

Then, the element stiffness matrix can be calculated by employing the volume integral on the annular structure of the element.The expression is given as

where Ωerepresents the annular structure of the axisymmetric element, and Aeis the cross-section of the axisymmetric element.Note that Keis not a constant matrix and will change with the element’s position.

As for the body forces(such as the centrifugal force and the gravity)and the surface forces(such as the pressure)of the element,the equivalent nodal load vector of the element could be calculated by

Fig.3 Rectangular axisymmetric element with four nodes.

where b is the unit body force,p is the unit surface force,and leis the edge line where the surface force acts.

By respectively assembling the element stiffness matrices Keand the equivalent nodal load vectors fe, the global stiffness matrix K of the structure and the equivalent load vector f of the structure could be obtained.

where G is the transformation matrix of element nodal degrees of freedom and structural nodal degrees of freedom, Neis the amount of elements.

By employing the principle of virtual work,the equilibrium equation of the axisymmetric structure could be written as

where u is the displacement vector of the structure.

By introducing the displacement constraints and solving the equilibrium equation, the displacement vector ueof each element could be obtained.Accordingly, the element strain field ε and the element stress field σ could be calculated.

The widely used von Mises stress criterion is selected as the failure criterion.The expression of von Mises stress σVMis

3.Maximum stress prediction method

In density-based topology optimization, the fixed finite element mesh is used instead of the body-fitted mesh.After optimization, the boundaries of the structure are usually jagged rather than smooth, and there are apparent blurred regions on the solid/void interfaces.The blurred regions often contain high-density elements with density values of 0.7–0.9, mediumdensity elements with density values of 0.4–0.6, and lowdensity elements with density values of 0.1–0.3, as shown in Fig.4.When the intermediate-density elements and jagged boundaries exist simultaneously,it may be inaccurate to calculate the maximum stress of the structure under designdependent loads.For the topology optimization with stress constraints, the accuracy of the calculated maximum stress of the structure would directly affect the optimization results of the structure.The optimized structure will be conservative if the calculated maximum stress is too large compared with the actual maximum stress.On the contrary, if the calculated maximum stress is too small, the optimized structure’s maximum stress may exceed the material’s allowable stress.Therefore, it is significant to use appropriate methods to compute the maximum stress of the structure accurately.

3.1.Calculation method of predicted maximum stress

The popular IVE method extrapolates the boundary stress of the structure under fixed loads through the internal elements near the boundaries.28However,in the presence of centrifugal loads, the maximum stress of the structure would be changed by the intermediate-density elements at the boundaries.In the finite element analysis, the element properties are consistent with the actual materials when the element density is 1 or 0.However, the material properties of the intermediatedensity elements are obtained by the interpolation functions.They are inconsistent with the actual materials, which is the main reason for the inaccurate stress calculation.

Fig.4 Jagged boundaries and blurred regions in density-based topology optimization.

This paper proposes a novel maximum stress prediction method for accurately predicting the actual maximum stress of the structure under centrifugal loads.The core idea is to:(A)newly define the predicted density to obtain stable boundaries with thin layers of gray elements (the predicted density is obtained by using the Heaviside function to map the filtered density); (B) innovatively propose a transition factor φ(namely the ratio of elements with predicted density values greater than 0.9 or less than 0.1 in all elements) to evaluate the impact of intermediate-density elements;(C)flexibly introduce two different stiffness penalty schemes to calculate the elastic modulus of elements according to the feedback of the transition factor φ; (D) further use the linear stress penalty to relax the stress field of the structure and obtain the predicted maximum stress.

Fig.5 shows the flowchart for calculating the predicted maximum stress.It mainly involves the following steps.

(1) Use the density filtering method to filter the initial design variables x and obtain the filtered density ρ～.A detailed description of the density filtering method is shown in Section 4.2.

(2) Use the Heaviside project function to map the filtered density ρ～and obtain the predicted density xH.The expression is given as

where βHcontrols the steepness of the Heaviside function.To be noted,compared with the filtered density ρ～,the value of the predicted density xHis much closer to 0 or 1.

(3) Introduce a transition factor φ to control the calculation method of the elastic modulus.The transition factor φ could be written as

where Ngis the number of elements with a predicted density greater than 0.9,Nlis the number of elements with a predicted density less than 0.1, and N is the total number of elements.

a.When the transition factor φ is smaller than a predefined threshold value φ*, the nonlinear stiffness penalty scheme is adopted to calculate the elastic modulus EHof each element.It could be expressed as

where Emaxis the elastic modulus of the actual material, and Emin(=10-6Emax) is a small value added to avoid the matrix singularity.

b.When the transition factor is greater than the threshold value φ*, the linear stiffness penalty scheme will be adopted to calculate the elastic modulus EHof each element.The expression is given as

Fig.5 Flowchart for calculating the predicted maximum stress.

To be noted,the predicted maximum stress will be used for the stress constraint relaxation method shown in Section 4.4.

3.2.Verification of axisymmetric finite element analysis program

The adopted finite element analysis program for the axisymmetric structure is developed in MATLAB.To verify the accuracy of the developed program,the axisymmetric model of the expanded turbine disk shown in Fig.2 is selected as the benchmark test case.The geometric parameters of the expanded turbine disk are set as R1=83 mm,R2=237 mm,W1=92 mm,and W2=40 mm.GH4169 is selected as the calculation material in this paper.The material properties at a constant temperature are used.The density is 8240 kg/m3,the Poisson’s ratio is 0.3, and the elastic modulus is 192 GPa.Note that the yield stress of the material is 1000 MPa.The selected element type is the 4-node axisymmetric rectangular element.The element size is set to 1 mm.The axisymmetric structure is finally discretized into 14168 elements.The rotation angular velocity is set to 1047.2 rad/s (namely 1 × 104r/min).The tensile loads caused by the high-speed rotation of turbine blades are set to 100 MPa and applied at the edge with a length of 40 mm on the right side of the turbine disk.The axial displacement constraint is applied at point 1 in the lower-left corner of the turbine disk.

Fig.6 shows the von Mises stress distribution of the expanded axisymmetric turbine disk, respectively calculated using the ANSYS software and the self-developed program.It can be seen that the stress distribution is entirely consistent.Table 1 shows the maximum von Mises stress and the maximum radial displacement of the expanded axisymmetric turbine disk, respectively calculated by the ANSYS software and the self-developed program.We can see that the maximum radial displacement of the axisymmetric structure calculated by the self-developed program is the same as that calculated by the ANSYS software;the maximum von Mises stress calculated by the self-developed program is only 0.246%lower than that calculated by the ANSYS software.

Therefore, the self-developed finite element analysis program for the axisymmetric structure is accurate and reliable.

3.3.Verification of maximum stress prediction method

Fig.6 Von Mises stress distribution of expanded axisymmetric turbine disk.

Table 1 Maximum von Mises stress and maximum radial displacement of expanded axisymmetric turbine disk respectively calculated by ANSYS software and self-developed program.

Fig.7 Schematic diagram of a Y-type axisymmetric structure.

To verify the proposed maximum stress prediction method, a Y-type axisymmetric structure with jagged boundaries and blurred regions is selected as the benchmark test case.Fig.7 shows the detailed shape and size information of the selected Y-type axisymmetric structure.The physical density of elements in the void and solid areas are set to 0 and 1, respectively.The elements denoted by 0.3, 0.5, and 0.7 mean that the physical density of these elements are set to 0.3, 0.5, and 0.7, respectively.GH4169 is selected as the calculation material.The rotation angular velocity is set to 1047.2 rad/s(namely 1 × 104r/min).The tensile loads are set to 100 MPa and applied at the edge with a length of 6 mm on the right side of the Y-type structure.

The structure could be reconstructed by using lines to connect diagonals of elements with a density of 0.5 and boundaries of elements with a density of 1.As shown in Fig.7, the structure surrounded by the red lines is the reconstructed geometry model.The ANSYS software is used to generate the bodyfitted mesh shown in Fig.8(a) and then conduct the finite element analysis.The corresponding von Mises stress distribution of the Y-type axisymmetric structure is shown in Fig.8(b).It can be seen that the actual maximum stress of the Y-type axisymmetric structure is 367.0 MPa.

Note that the transition factor φ of this Y-type axisymmetric structure is 0.9025.It is not necessary to verify the accuracy of the maximum stress prediction method when the transition factor is less than 0.9, because the transition factor φ must be greater than 0.9 at the end of optimization.Table 2 and Fig.9 show the predicted maximum stress of the Y-type axisymmetric structure using different values of the Heaviside parameter βH,respectively calculated using the linear and nonlinear stiffness penalty schemes.From them, we can see that: (A) compared with the actual maximum stress calculated by the ANSYS software, the difference of the predicted maximum stress calculated using the nonlinear stiffness penalty scheme varies from 2.081% to 5.364%, while that calculated using the linear stiffness penalty scheme varies from - 0.273%to - 0.109%; (B) when the Heaviside parameter βHincreases from 4 to 16, the predicted maximum stress calculated using the nonlinear stiffness penalty scheme decreases gradually;when βHincreases from 16 to 25, the change of the predicted maximum stress is not apparent; (C) when the Heaviside parameter βHincreases from 4 to 25, the predicted maximum stress calculated using the linear stiffness penalty scheme is almost unchanged and approximately equal to the actual maximum stress.

Therefore, the proposed maximum stress prediction method for the axisymmetric structure is effective and accurate.

4.Topology optimization method of turbine disk

4.1.Formulation of two kinds of optimization problems

Two kinds of topology optimization problems with stress constraints are considered in this paper.One is to minimize the structure’s mass under stress constraints, and the other is to minimize the structure’s compliance under stress and mass constraints.The structure is discretized into N rectangular elements to generate a vector x containing N design variables.Each design variable xerepresents the state of the corresponding structural element.A value of 1 indicates that the element has material, while a value of 0 indicates that the element is void.In order to facilitate the solution and fully use the popular mathematical programming algorithm, the design variables are defined as continuous parameters ranging from 0 to 1.

For the design of the turbine disk,the structure’s stress is a critical evaluation factor to be considered.If the structure’s maximum stress exceeds the material’s yield stress,the material will fail and cannot satisfy the service requirements.Therefore,the optimization problem of minimizing structural mass with stress constraints is always a hot issue.

The optimization formulation could be written as

where Mfis the mass fraction of the turbine disk, ρeis the physical density of the e-th element, veis the volume of the e-th element, σeis the stress of the e-th element, and σlis the yield stress of the material.

It is also significant to seek the maximum stiffness design of the structure with stress and mass constraints.The corresponding optimization formulation could be written as

where C is the structural compliance,and Mf*is the allowable maximum mass fraction.

Fig.8 Body-fitted mesh and von Mises stress distribution of Y-type axisymmetric structure generated or calculated by ANSYS software.

Table 2 Different settings of βH and the corresponding predicted maximum stress of Y-type axisymmetric structure.

Fig.9 Predicted maximum stress of Y-type axisymmetric structure using different values of Heaviside parameter βH.

To be noted, the different formulations for the minimum mass problem and the minimum compliance problem considering stress constraints have been introduced by Holmberg et al.26The relevant sensitivity analysis has also been fully discussed.However,the centrifugal loads are not considered.This paper mainly studies the topology optimization of the turbine disk, considering the prediction and constraints of maximum stress under centrifugal loads.The updating scheme of the relax coefficient c is similar to that in the work of Le et al.21By comparison, we use the prediction maximum stress instead of the maximum stress in the calculation of c.

4.2.Density filtering and material property interpolation function

where rminis the given filter radius, ciis the center of the i-th element, and ceis the center of the e-th element.

To obtain a black-and-white layout without gray elements,the Heaviside projection function is employed to map the filtered density ρ～to obtain the physical density ρe:39–41

where β controls the steepness of the Heaviside function; η is the threshold value, which controls the inflection point of the threshold projection.Referring to previous studies17,41, η is set to 0.5 in this paper.

The interpolation function is used to establish the relationship between the properties of the intermediate-density elements and the actual material.The Solid Isotropic Material with Penalization (SIMP) and RAMP methods are the two most popular interpolation functions.Due to the centrifugal loads on the turbine disk, the SIMP method would lead to low sensitivity in the low-density area and hence numerical difficulties in optimization.42Therefore, the RAMP method is selected as the stiffness penalty function in this paper and written as43

where qcis the stiffness penalty factor.The elastic modulus of the intermediate-density element could be finally written as

For the intermediate-density elements,the von Mises stress σVMcalculated by Eq.(13)still needs to be relaxed.The SIMPlike relaxation method has been widely used in the stressconstrained topology optimization of structures under fixed loads.However, stress-constrained topology optimization under the centrifugal load has rarely been studied.Considering the numerical difficulties of the SIMP method,42this paper adopts the RAMP method to penalize the stress of the intermediate-density elements of the turbine disk under centrifugal loads.It could be expressed as44

4.3.Stress measure

When the structure is discretized into N elements,the stress of each element needs to be calculated and less than the material yield stress,which results in N stress constraints.It is challenging to calculate the sensitivity values of the N stress constraints using either the direct or adjoint methods.Moreover,so many constraints often lead to expensive computation costs and even make the optimization difficult to converge.To reduce the number of constraints and improve the calculation efficiency,the maximum von Mises stress is used to convert the N local stress constraints into a global stress constraint.

However, the maximum von Mises stress is nondifferentiable and challenging to use directly in the densitybased topology optimization that relies on gradient information.Hence, there is an urgent need for a differentiable stress index to replace the maximum von Mises stress.The global aggregation functions are usually utilized to overcome this difficulty.23The most commonly used aggregation methods are the p-norm and Kreisselmeier-Steinhauser (KS) functions.17Verbart et al.45studied and compared the p-norm and KS functions in detail.The main difference between the p-norm function and the KS function is that the p-norm function only deals with non-negative value problems,while the KS function does not have the limitation.45–46In this paper,the stress value is non-negative, so both the p-norm and KS functions can be used.For simplicity and generality, the p-norm function is adopted in this paper to aggregate the N stress and approximate the maximum von Mises stress.The stress constraint is then transformed as

where σ～is the approximate maximum stress, and P is the norm factor.When P approaches 1, σ～is close to the sum of the stress values of the N elements.In this situation, the aggregation function is relatively smooth, and the nonlinearity degree is relatively low, yet it deviates greatly from the maximum von Mises stress.When P is large enough, σ～is approximately equal to the maximum stress of the N elements, yet it is highly nonlinear and difficult to optimize.In order to balance the accuracy of the approximation and the nonlinearity of the p-norm function, the value of P is set to 8 in this paper.

4.4.Stress constraint relaxation method

The p-norm function aggregates the local stress constraints into a global stress constraint.The approximate maximum stress σ～is used to replace the structure’s actual maximum von Mises stress.Although P is set to 8,there is still an apparent error between the approximate and actual maximum stress values, which may affect the effectiveness of stress constraints.21Inspired by Le et al.21, a new normalized global stress method is developed to deal with the stress constraint by employing the predicted maximum stress.The stress constraint is finally transformed into the following form:

In order to avoid the rapid oscillation of the stress constraint and make the iterative process more stable, the relaxation coefficient c is updated every five steps, and the iterative maximum stress is gradually updated according to the predicted maximum stress in the current step and the iterative maximum stress in the last step.The updating scheme could be expressed as follows:

(1) In the first step,

To be noted, the stress constraint essentially is that the predicted maximum stress should be less than or equal to the material’s yield stress.

4.5.Sensitivity analysis

The popular gradient-based mathematical programming algorithm, the Method of Moving Asymptotes (MMA), is adopted to perform the topology optimization of the axisymmetric turbine disk under centrifugal loads.Thus,the first-order sensitivity information of the objective function and constraints should be derived to update the design variables.

The derivative of the mass fraction Mfof the turbine disk with respect to the physical density ρeof the e-th element could be written as

where ω denotes the angular velocity, ρ0is the material density.

It is challenging to directly calculate the derivative of the approximate maximum stress σ～with respect to the physical density.Therefore, the adjoint method is adopted, and the augmented form of the approximate maximum stress σ～is written as

Finally, in terms of the chain rule, the derivatives of the mass fraction Mf,the structural compliance C,and the approximate maximum stress σ～with respect to the design variable xecould be easily obtained.

5.Results and discussion

After verifying the effectiveness of the axisymmetric finite element analysis program and the maximum stress prediction method, the topology optimization of the turbine disk is conducted using the self-developed program.Fig.2 shows the axisymmetric model of the expanded turbine disk.The involved geometry, material, element, and loading parameters have been introduced in Section 3.2.The yield stress of the material is set to 1000 MPa.MMA is adopted as the optimization solver.47To improve the convergence performance, the initial values of design variables in the design domain are set to 0.5,and the move limit of design variables in each iteration is set to 0.1.Additionally, the filter radius rminin the density filter function is set to 3.5.The initial value of the Heaviside parameter βHin Eq.(15)is set to 4 and is increased by 4 every 20 iterations.The maximum value of βHis set to 20.The initial value of the Heaviside parameter β in Eq.(26) is set to 4 and multiplied by a factor of 1.2 every 20 iterations.The maximum value of β is set to 8.The stiffness penalty factor qcin Eq.(27)is set to 4.Moreover,the maximum number of iterations of the optimization is set to 400.The optimization is terminated when the maximum change of the design variables is less than 0.05 and the maximum predicted stress is smaller than the yield stress σlof the material, or when the number of iterations reaches the predefined maximum value of 400.

Section 5.1 discusses the selection of the penalty factor in the calculation of stress of intermediate densities.Then, the effects of the threshold on the transition factor φ,which is used in the proposed maximum stress prediction method, are discussed in Section 5.2.In Section 5.3,the effects of the fixed tensile and centrifugal loads on the optimization results of the turbine disk are studied.Subsequently, two popular optimization formulations are compared,which are the mass minimization with the stress constraints and the compliance minimization with the stress and mass constraints.Finally,the optimized turbine disk is reconstructed to obtain smooth boundaries and analyzed by the ANSYS software to illustrate the effectiveness of the proposed maximum stress prediction method and the developed stress constraint method.

5.1.Effects of stress penalty factor

In Ref.21, only the fixed loads are considered, the SIMP method is adopted for stiffness penalty, and the SIMP-like method is adopted for stress penalty.However, the structure will have a parasitic effect if the design-dependence loads are considered.42This section investigates the RAMP method for the stress penalty of the topology optimization of the turbine disk under centrifugal loads.The RAMP method has unique properties.Different stress penalty factors will affect the topology optimization of the structure under the stress constraint.Fig.10 (a) shows the graph of the RAMP function in Eq.(29) using different stress penalty factors.At the same time,in order to show the punishment intensity of the nonlinear penalty compared with the linear penalty, we establish an intensity function I(ρe) to describe the punishment intensity,and the expression is

Fig.10(b)shows the graph of the intensity functions using different stress penalty factors.It can be seen that the higher the peak value of the curve, the greater the nonlinear punishment intensity.For the listed stress penalty factors -0.99,-0.97, -0.95, -0.93, -0.9, -0.8, -0.7, -0.6, and -0.5, the corresponding density values at the peak points of the curves are 0.09,0.15,0.18,0.21,0.24,0.31,0.35,0.39,and 0.41 respectively.The smaller the penalty factor is, the farther the peak point is from the x-axis and the closer it is from the y-axis.Additionally, the peak point is a transition point of the stress penalty.On the left side of the peak point, the punishment intensity increases rapidly as the density value increases; on the right side of the peak point, the punishment intensity decreases gradually.Obviously, the penalty factors directly affect the calculation of the stress field of the structure.

Fig.11 shows the configurations of the optimized turbine disk when the stress penalty factors are set to -0.5, -0.6,-0.7, -0.8, -0.9, -0.93, -0.95, -0.97, and -0.99, respectively.Note that the threshold value φ*of the transition factor of the proposed maximum stress prediction method is set to 0.9 in this section.From Fig.11, we can see that: (A) when the stress penalty factor is greater than -0.9, the boundaries of the optimized turbine disk are not clear;(B)when the stress penalty factor is less than or equal to -0.9, the configurations of the optimized turbine disk are almost similar.Table 3 lists the mass fraction of the optimized turbine disk when the stress penalty factors are set to-0.5,-0.6,-0.7,-0.8,-0.9,-0.93,-0.95,-0.97,and-0.99,respectively.It can be seen that when the stress penalty factor is-0.95,the mass fraction of the optimized turbine disk is the smallest.It is 0.1006.

Therefore, the selection of the stress penalty factor has a significant impact on the optimization results of the turbine disk.The inappropriate penalty factor will lead to the failure of the topology optimization.The optimization performs best when the stress penalty factor is between-0.95 and-0.93.For ease of comparison, the stress penalty factor is set to -0.95 in the following sections.

5.2.Effects of threshold value of transition factor

The threshold value φ*of the transition factor in the proposed maximum stress prediction method determines the predicted maximum stress.It is of great significance to explore the effects of the threshold value of the transition factor on optimization results.Fig.12 shows the convergence curves of the mass fraction and the configurations of the optimized turbine disk when the threshold values are set to 0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,and 1.0,respectively.Table 4 lists the mass fraction of the optimized turbine disk using different threshold values of the transition factor.From them, we can see that: (A) when the threshold value is less than or equal to 0.4, the configurations of the optimized turbine disk are unreasonable,and there are large fluctuations in the convergence curves; (B) when the threshold value is greater than or equal to 0.5, the configurations of the optimized turbine disk are reasonable and innovative, the convergence curves are smooth, and the mass reductions of the optimized turbine disk are significant; (C)when the threshold value is increased from 0.5 to 0.6, there is a noticeable difference in the optimized configuration,changing from a single inner-ring-cavity structure to a double inner-ring-cavity structure; (D) when the threshold value is between 0.6 and 0.9,the mass reductions of the optimized turbine disk are approximately the same; (E) when the threshold value is 0.9,the convergence curves are pretty smooth,and the convergence speed is the fastest.

It is interesting to explore why the configuration transitions from a single inner-ring-cavity structure to a double inner-ringcavity structure.Fig.13 shows the convergence curves of the transition factor and the predicted maximum stress at a threshold of 0.5.It can be seen that: (A) in the 20th iteration, the transition factor exceeds 0.5, the stiffness penalty in the maximum stress prediction is changed from nonlinear penalty to linear penalty, and the predicted maximum stress is reduced from 1046 MPa to 660 MPa; in the current iteration, the predicted maximum stress value is less than the material’s yield stress; (B) in the subsequent iterations, the transition factor gradually increases, and the mass fraction gradually decreases until the predicted maximum stress value meets the stress constraint.Fig.14 shows the density fields and stress contour of the turbine disk at typical iteration steps with a threshold value of 0.5.It can be seen that: (A) the maximum stress is concentrated in the middle area on the left side;(B)in subsequent iterations, the turbine disk gradually evolves into a single innerring-cavity structure.Fig.15 shows the convergence curves of the transition factor and the predicted maximum stress at a threshold of 0.6.It can be seen that:(A)in the 21st iteration,the transition factor exceeds 0.6, and the predicted maximum stress is reduced from 1063 MPa to 654 MPa;(B)in the subsequent iterations, the transition factor gradually increases to 1,and the predicted maximum stress value is gradually approaching the material’s yield stress.Fig.16 shows the density fields and stress contour of the turbine disk at typical iteration steps with a threshold value of 0.6.It can be seen that:(A)the maximum stress is distributed in the upper and lower areas on the left side;(B)it changes the optimization direction,and the turbine disk gradually evolves into a double inner-ring-cavity structure in the subsequent iterations.That is to say,the initial distribution of the maximum stress will affect the final configuration of the optimized turbine disk.

Fig.10 Graphs of RAMP and intensity functions using different stress penalty factors.

Fig.11 Configurations of optimized turbine disk using different stress penalty factors.

Table 3 Mass fraction of optimized turbine disk using different stress penalty factors.

Fig.12 Convergence curves of mass fraction and configurations of optimized turbine disk using different threshold values of transition factor.(Left: Configurations of optimized turbine disk.Right: Convergence curves of mass fraction).

It is also of great significance to explore the interaction between the threshold value of the transition factor and the stress penalty factor.Table 5 lists the mass fraction of the optimized turbine disk using different threshold values of the transition factor and the stress penalty factors.It can be seen that:(A) when the stress penalty factor is -0.8, the minimum mass fraction of 0.1508 is obtained at the threshold value of 0.9;(B)when the stress penalty factor is-0.9,the minimum mass fraction of 0.1024 is obtained at the threshold value of 0.9; (C)when the stress penalty factor is -0.95 and the threshold values are between 0.6 and 0.9, the mass reductions of the optimized turbine disk are approximately the same; (D) when the stress penalty factor is -0.99 and the threshold values are between 0.7 and 0.9, the mass reductions are approximately the same; (E) when the stress penalty factor is greater than or equal to -0.9, the threshold value has a more significant impact on the results; (F) when the stress penalty factor is less than or equal to -0.95, the threshold value has less impact on the results; (G) when the threshold values are respectively set to 0.5, 0.6, 0.7, 0.8, 0.9, and 1.0, the minimum mass fraction is obtained throughout at the stress penalty factor of -0.95.

Table 4 Mass fraction of optimized turbine disk using different threshold values of transition factor.

Fig.13 Convergence curves of transition factor and predicted maximum stress at a threshold value of 0.5.

Fig.14 Density fields and stress distributions of turbine disk at typical iteration steps at a threshold value of 0.5.(Left: Density fields.Right: Stress distribution contour).

Therefore,the selection of the threshold value of the transition factor of the proposed maximum stress prediction method significantly affects the topology optimization results of the turbine disk.The threshold value can control the evolution of the optimized turbine disk from the single inner-ring-cavity structure to the double inner-ring-cavity structure.The inappropriate threshold value will lead to the failure of the topology optimization.Considering the optimization effect,convergence speed, and convergence stability, the threshold value of the transition factor is set to 0.9 in the following sections.

5.3.Effects of loads

The turbine disk mainly bears the centrifugal loads and the fixed tensile loads.The loads directly determine the distribution of the stress field of the turbine disk.Exploring the influence of the centrifugal and fixed tensile loads on the optimization results is necessary.

Fig.15 Convergence curves of transition factor and predicted maximum stress at a threshold value of 0.6.

Fig.16 Density fields and stress distributions of turbine disk at typical iteration steps at a threshold value of 0.6.(Left: Density fields.Right: Stress distribution contour).

Table 5 Mass fraction of optimized turbine disk using different threshold values of transition factor and stress penalty factors.

5.3.1.Effects of fixed tensile loads

In this section,the fixed tensile loads are changed from 0 MPa to 600 MPa, while the rotating speed is set to 1 × 104r/min from beginning to end.Fig.17 shows the configurations of the optimized turbine disk under different fixed tensile loads F of 0 MPa, 20 MPa, 40 MPa, 60 MPa, 100 MPa, 150 MPa,200 MPa, 250 MPa, 300 MPa, 400 MPa, 500 MPa, and 600 MPa.Table 6 shows the mass fraction and maximum stress of the optimized turbine disk under different fixed tensile loads.From them, we can see that: (A) when the fixed tensile loads are less than 40 MPa,the configurations of the optimized turbine disk are unreasonable, the material distribution is discontinuous, and the final design tends to be material free; (B)when the fixed tensile loads are between 60 MPa and 500 MPa,reasonable and innovative configurations are obtained, the mass fraction increases gradually with the increase of the fixed tensile loads, and the turbine disk gradually evolves from the double inner-ring-cavity structure into the single inner-ringcavity structure and further evolves into a solid structure;(C) when the fixed tensile loads are greater than 600 MPa,the optimization could not converge because the stress fails to meet the constraint.

Therefore,the optimized configuration of the turbine disk is closely related to the fixed tensile loads.When the fixed tensile loads are small enough,the stress distribution mainly depends on the centrifugal loads,and the final design tends to be material free.When the fixed tensile loads increase, the mass fraction increases so as to meet the stress constraint, and the optimized configuration evolves from the double inner-ringcavity structure to the single inner-ring-cavity structure and a solid structure.That is to say,the solid structure performs best in the bearing capacity for the tensile loads, followed by the single inner-ring-cavity and the double inner-ring-cavity structure.

5.3.2.Effects of centrifugal loads

Compared with the fixed tensile loads, the centrifugal loads affect the calculation of each element’s equivalent loads and stresses.In this section, the rotating speed is changed from 0 r/min to 2.5 × 104r/min, while the fixed tensile loads are set to 100 MPa.Fig.18 shows the configurations of the optimized turbine disk at different rotating speeds ω of 0 r/min,3×103r/min,6×103r/min, 8×103r/min, 1×104r/min, 1.2×104r/min, 1.4 × 104r/min, 2 × 104r/min, and 2.5 × 104r/min.Table 7 shows the mass fraction and maximum stress of the optimized turbine disk at different rotating speeds.From them,we can see that:(A)when the rotating speed is less than or equal to 1.4×104r/min,reasonable and innovative configurations are obtained, and the mass fraction increases gradually with the increase of the rotating speed;(B)the turbine disk gradually evolves from the double inner-ring-cavity structure into the triple inner-ring-cavity structure and further evolves into a single inner-ring-cavity structure; (C) when the rotating speed is greater than or equal to 2 × 104r/min, the optimizations cannot converge since the maximum stress exceeds the material’s yield stress.

Table 6 Mass fraction and maximum stress of optimized turbine disk under different fixed tensile loads.

Fig.17 Configurations of optimized turbine disk under different fixed tensile loads.

Therefore, the magnitude of the centrifugal loads significantly affects the optimized configuration of the turbine disk.The single inner-ring-cavity structure performs best in the bearing capacity for the centrifugal loads, followed by the triple inner-ring-cavity structure and the double inner-ring-cavity structure.

Table 7 Mass fraction and maximum stress of optimized turbine disk at different rotating speeds.

5.4.Comparison of two optimization formulations

In order to obtain the structure with the maximum stiffness,the structural compliance minimization problem with mass constraints has been widely researched and is still a hotspot.15,26The optimization results of the turbine disk would be more valuable and practical if the stress constraint is introduced.Thus,it is of great significance to compare the two optimization formulations,namely the mass minimization with the stress constraints and the compliance minimization with the stress and mass constraints.For ease of comparison,the stress penalty factors of the two optimization formulations are set to -0.95, the threshold values of the transition factors of the two optimization formulations are set to 0.9, the fixed tensile loads of the two optimization formulations are set to 100 MPa,and the rotating speeds of the two optimization formulations are set to 1×104r/min.Additionally,the allowable maximum mass fraction of the compliance minimization problem is set to 0.1006,which is the mass fraction of the optimized turbine disk of the mass minimization problem.Fig.19 shows the convergence curves of the two optimization formulations and the corresponding configurations of the optimized turbine disk.Table 8 compares the optimization results of the two optimization formulations.Combining Fig.19 and Table 8,we can see that:(A)the optimized mass fraction of the compliance minimization problem is 0.1033, which violates the given mass constraint slightly; (B) the configurations of the optimized turbine disk of the two optimization formulations are similar; (C) compared with the mass minimization problem,the convergence curve of the compliance minimization problem is not smooth enough; (D) compared with the mass minimization problem, the mass fraction of the compliance minimization problem is increased by 2.684%, yet the structural compliance is reduced by 1.728%; (E) although, strictly speaking, the optimization of the compliance minimization problem does not converge, the optimized scheme is still valuable.

Fig.18 Configurations of optimized turbine disk at different rotating speeds.

Fig.19 Convergence curves of two optimization formulations and corresponding configurations of optimized turbine disk.

Table 8 Comparison of optimization results of two optimization formulations.

Therefore,similar results are obtained by the two optimization formulations.The mass minimization problem performs better in achieving the minimum mass.The compliance minimization problem performs better in order to obtain the minimum compliance.However,selecting the allowable maximum mass fraction for the compliance minimization problem is challenging.Inappropriate selection may lead to convergence issues.

5.5.Model reconstruction and performance analysis

The optimized turbine disk contains jagged boundaries and blurred regions and is difficult to use directly for processing and manufacturing.Therefore, it is necessary to reconstruct the structural boundaries to obtain smooth boundaries and accurately analyze the performance of the optimized turbine disk.The optimal scheme of the mass minimization problem introduced in Section 5.4 is selected as the test case for reconstruction.Fig.20(a) shows the density distribution before reconstruction.The fourth-order polynomial curve is used to fit the boundaries of the optimized turbine disk by employing open-source software for 2D topology optimization postprocessing.48Fig.20(b) shows the stress contour of the reconstructed turbine disk.The ANSYS software is used to discretize the reconstructed structure with the same element size in topology optimization.Fig.20(c) shows the generated body-fitted mesh of the reconstructed turbine disk.The total number of elements is 3012.The finite element analysis for the reconstructed structure is conducted in ANSYS software.Fig.20(d) shows the turbine disk’s von Mises stress distribution after reconstruction.The maximum stress of the reconstructed turbine disk is 1004 MPa, which is only 0.4% larger than the predicted maximum stress of the turbine disk before reconstruction.

Therefore,the maximum stress of the reconstructed turbine disk is consistent with the predicted maximum stress in topology optimization.It further verifies that the proposed maximum stress prediction method for the axisymmetric structure is accurate and that the developed stress constraint method for the topology optimization of the turbine disk is effective.

5.6.Comprehensive analysis of topological configurations

It has been proved that the optimized configuration of the turbine disk is closely related to the threshold value of the transition factor, the fixed tensile loads, and the centrifugal loads.Four topological configurations of the turbine disk, including the solid structure, the single inner-ring-cavity structure, the double inner-ring-cavity structure, and the triple inner-ringcavity structure, are generated during the topology optimization.The topological configuration is a critical consideration for the turbine disk’s design and is meaningful for engineering.Therefore, this section attempts to make a comprehensive analysis of this phenomenon.

Fig.21 shows the stress distributions of the optimized turbine disk at different threshold values of 0.5 and 0.6.Here,the rotating speed and the fixed tensile load are set to 1 × 104r/min and 100 MPa,respectively.According to Fig.21 and Section 5.2,we can see that:(A)the threshold value affects the optimization direction by changing the predicted maximum stress,which leads to different local optimal solutions for the two schemes;(B)the stress of the two schemes are mainly distributed between 700 MPa and 1000 MPa,which are relatively uniform;(C) for the single inner-ring-cavity structure, the stress on the left middle side is very large, while the stress on the left upper side and the left lower side are very small; (D) for the double inner-ring-cavity structure, the appearance of the left cavity causes the redistribution of the loads, the high stress migrates from the left middle side to the left upper and left lower sides,and the mass of the turbine disk decreases accordingly.

Fig.20 Model reconstruction and performance analysis of optimized turbine disk.

Fig.21 Stress distributions of optimized turbine disk at different threshold values of 0.5 and 0.6.

According to Section 5.3.1, four types of turbine disks are obtained under different optimization conditions, namely different fixed tensile loads of 100 MPa, 250 MPa, 400 MPa,and 500 MPa.Fig.22 shows the stress distributions of the four optimized turbine disks under different working conditions,namely fixed tensile loads.Here, the threshold value of the transition factor and the rotating speed are set to 0.9 and 1 × 104r/min, respectively.According to the second column of Fig.22,we can see that:(A)with the increase of the tension loads of the optimization condition,the mass of the optimized turbine disk increases gradually, and the structure changes obviously, evolving from a large double inner-ring-cavity structure to a small double inner-ring-cavity structure,a single inner-ring-cavity structure, and a solid structure; (B) from Subgraph (b) to Subgraph (e), the left-ring-cavity becomes smaller, and the high stress area (greater than 900 MPa) on the left middle side becomes larger,extending to the left upper side and the left lower side;(C)from Subgraph(e)to Subgraph(h), the left-ring-cavity disappears and evolves to a low stress area (less than 700 MPa), and accordingly the high stress area(greater than 900 MPa) on the left middle side becomes smaller;the right ring-cavity becomes smaller and the stress around it decreases; (D) from Subgraph (h) to Subgraph (k), the high stress area (greater than 900 MPa) on the left middle side becomes larger, extending to the left upper side and left lower side; the right-ring-cavity disappears and evolves to a low stress area (less than 700 MPa).By comparing the second and third columns of Fig.22, we can see that: (A) when the tensile loads of the working condition are greater than those of the optimization condition, the stress of the whole turbine disk increases synchronously, and the stress of the most areas of the turbine disk exceeds the stress limit of the material; (B)most positions of the turbine disk need to be thickened to bear the larger tensile loads of the working condition; this is why the turbine disk evolves from a large double inner-ring-cavity structure to a small double inner-ring-cavity structure,a single inner-ring-cavity structure, and a solid structure with the increase of the tension loads of the optimization condition.By comparing the first and second columns of Fig.22, we can see that: (A) when the tensile loads of the working condition are less than those of the optimization condition, the stress of the whole turbine disk synchronously decreases, and some extremely-low-stress areas appear; (B) most positions of the turbine disk need to be thinned to obtain a lighter structure;this is why the turbine disk evolves from a solid structure to a single inner-ring-cavity structure, a small double innerring-cavity structure, and a large double inner-ring-cavity structure with the decrease of the tension loads of the optimization condition.These phenomena show that:(A)the existence of inner-ring-cavity helps to reduce the mass of the turbine disk,but changes the stress distributions and increases the stress around it; (B) with the increase of the tensile loads,the inner ring cavity gradually reduces or even disappears due to the limitation of the design area, and the stress around it decreases accordingly, which makes the structure have greater tensile bearing capacity.

Fig.22 Stress distributions of four optimized turbine disks under different working conditions, namely different fixed tensile loads.(Note: the four types of turbine disks are obtained under different optimization conditions, namely different fixed tensile loads).

Fig.23 Stress distributions of three optimized turbine disks under different working conditions,namely different rotating speeds.(Note:the three turbine disks are obtained under different optimization conditions, namely different rotating speeds).

According to Section 5.3.2, three types of turbine disks are obtained under different optimization conditions, namely different rotating speeds of 1 × 104r/min, 1.2 × 104r/min, and 1.4 × 104r/min.Fig.23 shows the stress distributions of the three optimized turbine disks under different working conditions, namely different rotating speeds.Here, the threshold value of the transition factor and the fixed tensile loads are set to 0.9 and 100 MPa, respectively.According to the second column of Fig.23,we can see that:(A)with the increase of the rotating speed of the optimization condition, the mass of the optimized turbine disk increases gradually, and the structure changes obviously, evolving from a double inner-ring-cavity structure to a triple inner-ring-cavity structure and a large single inner-ring-cavity structure; (B) from Subgraph (b) to Subgraph (e), the right-ring-cavity becomes larger, the left-ringcavity also becomes larger and evolves into two ring-cavities,and the material on the left side of the structure increases;(C) from Subgraph (e) to Subgraph (h), the three small ringcavities evolve into a large single ring-cavity, and the material on the left side of the structure increases obviously, extending to the left upper side and left lower side.By comparing the second and third columns of Fig.23, we can see that: (A) when the rotating speeds of the working condition are greater than those of the optimization condition, the closer to the left side of the structure(namely the rotational center axis),the greater the stress increases;only small parts of stress on the left side of the turbine disk exceed the stress limit of the material; (B) the material on the right side needs to be further removed so that the material on the left side could bear larger rotating speeds of the working condition; this is why the right-ring-cavity becomes larger and finally the turbine disk evolves into a large single ring-cavity structure with the increase of the rotating speeds of the optimization condition.By comparing the first and second columns of Fig.23, we can see that: (A) when the rotating speeds of the working condition are less than that of the optimization condition, the stress on the left side of the turbine disk decreases significantly,while the stress on the right side decreases slightly;(B)the material on the left side needs to be removed and the new ring-cavity needs to be generated to obtain a lighter structure; this is why the turbine disk evolves from a large single ring-cavity structure to a triple innerring-cavity structure and a double inner-ring-cavity structure with the decrease of the rotating speeds of the optimization condition.

Therefore, the three factors, namely the threshold value of the transition factor,the fixed tensile loads,and the centrifugal loads, together determine the final configuration of the optimized turbine disk.The threshold value could affect the optimization direction by changing the predicted maximum stress,which may lead to different local optimal solutions.Both the fixed tensile loads and the centrifugal loads could affect the optimized configuration of the turbine disk by changing the stress level and the stress distributions.The difference is that the effect of fixed tensile loads on the stress level of the structure is global, while the effect of the centrifugal loads tends to be regional.The stress of the whole turbine disk increases synchronously with the increase of fixed tensile loads, while with the increase of the rotating speeds, the closer to the rotational center axis of the structure, the greater the stress increases.

6.Conclusions

To reduce the influence of jagged boundaries of the mesh and gray densities near the solid/void interfaces on the calculation of structural stress field under centrifugal loads,a new method to predict the maximum stress is proposed in this paper.In the method, the Heaviside function is used to project the filtered density to obtain the newly defined predicted density,enabling stable boundaries with thin layers of gray elements and a more stable and accurate maximum stress; the transition factor is proposed to flexibly select the linear or nonlinear stiffness penalty schemes for calculating the elastic modulus of elements;the linear stress penalty is further employed to relax the structure’s stress field, and the predicted maximum stress is finally obtained.The effectiveness of the proposed method is verified by the analysis of the Y-type structure with smooth boundaries and the topology optimization of the turbine disk.The RAMP interpolation is used in both stiffness and stress penalization for the structural stress field calculation and the sensitivity analysis.The p-norm function is adopted to aggregate the local stresses and approximate the maximum von Mises stress.A practical updating scheme of the stress constraint in the topology optimization is developed using the predicted maximum stress.The influences of the stress penalization factor, the threshold value of the transition factor in the maximum stress prediction,and the loads on the optimization results are investigated.Two kinds of optimization formulations, namely the mass minimization with the stress constraints and the compliance minimization with the stress and mass constraints, are also compared.An optimized turbine disk is reconstructed to obtain smooth boundaries and analyzed using the ANSYS software to further illustrate the effectiveness of the proposed maximum stress prediction method and the developed stress constraint method.A comprehensive analysis of the different topological configurations of the turbine disk is conducted to provide more support for engineering design.From the results of this study,some findings of this work could be summarized as follows:

(1) The maximum radial displacement of the axisymmetric structure calculated by the self-developed axisymmetric finite element analysis program is the same as that calculated by the ANSYS software; the maximum von Mises stress calculated by the self-developed program is only 0.246% lower than that calculated by the ANSYS software.It illustrates that the self-developed finite element analysis program for the axisymmetric structure is accurate and reliable.

(2) For the Y-type axisymmetric structure with jagged boundaries and blurred regions,the predicted maximum stress of the proposed method is approximately equal to the actual maximum stress calculated by the ANSYS software.It illustrates that the proposed maximum stress prediction method for the axisymmetric structure is effective and accurate.

(3) The selection of the stress penalty factor has a significant impact on the optimization results of the turbine disk.The inappropriate penalty factor will lead to the failure of the topology optimization.The optimization performs best when the stress penalty factor is between-0.95 and -0.93.For ease of comparison, the stress penalty factor is set to -0.95 in this paper.

(4) The selection of the threshold value of the transition factor in the proposed maximum stress prediction method significantly affects the topology optimization results of the turbine disk.The threshold value can control the evolution of the optimized turbine disk from the single inner-ring-cavity structure to the double inner-ringcavity structure.The inappropriate threshold value would lead to the failure of the topology optimization.Considering its effects on the optimization results, convergence speed and convergence stability, the threshold value of the transition factor is set to 0.9 in this paper.

(5) The optimized configuration of the turbine disk is closely related to the fixed tensile loads.When the fixed tensile loads are small enough, the stress distribution mainly depends on the centrifugal loads, and the final design tends to be material free.When the fixed tensile loads increase, the mass fraction increases so as to meet the stress constraint, and the optimized configuration evolves from the double inner-ring-cavity structure to the single inner-ring-cavity structure and further evolves into a solid structure.That is to say, the solid structure performs best in the bearing capacity for the tensile loads, followed by the single inner-ring-cavity and the double inner-ring-cavity structure.

(6) Compared with the fixed tensile loads, the centrifugal loads affect the calculation of each element’s equivalent loads and stresses.The magnitude of the centrifugal loads significantly affects the optimized configuration of the turbine disk.When the centrifugal loads are large enough, the optimizations cannot converge since the maximum stress exceeds the material’s yield stress.The large single inner-ring-cavity structure performs best in the bearing capacity for the centrifugal loads, followed by the triple inner-ring-cavity structure and the double inner-ring-cavity structure.

(7) The two optimization formulations obtain similar results.The mass minimization problem performs better in designing the turbine disk for the minimum mass.The compliance minimization problem performs better in designing the turbine disk for the minimum compliance.However, selecting the allowable maximum mass fraction for the compliance minimization problem is challenging.Inappropriate selection may lead to convergence problems.

(8) The maximum stress of the reconstructed turbine disk is only 0.4% larger than the predicted maximum stress of the turbine disk before reconstruction.It further verifies that the proposed maximum stress prediction method for the axisymmetric structure is accurate and that the developed stress constraint method for the topology optimization of the turbine disk is effective.

(9) The threshold value could affect the optimization direction by changing the predicted maximum stress, which may lead to different local optimal solutions.Both the fixed tensile loads and the centrifugal loads could affect the optimized configuration of the turbine disk by changing the stress level and the stress distributions.The difference is that the effect of fixed tensile loads on the stress level of the structure is global, while the effect of the centrifugal loads tends to be regional.

It is worth noting that this paper carries out two finite element analyses in one iteration.The first finite element analysis is used for structural stress field calculation and sensitivity analysis, and the second finite element analysis is used to predict the maximum stress of the structure.It leads to a slight increase in computational cost, which is acceptable for 2D problems and would be further improved for 3D problems in our future work.

Declaration of Competing 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.

Acknowledgement

This work was co-supported by the National Natural Science Foundation of China (Nos.52005421 and 12102375), the Natural Science Foundation of Fujian Province of China(No.2020J05020),the National Science and Technology Major Project,China (No.J2019-I-0013-0013), the Fundamental Research Funds for the Central Universities,China (No.20720210090), and the project funded by the China Postdoctoral Science Foundation (Nos.2020M682584 and 2021T140634).