A nontraditional method for reducing thermoelastic stresses of variable thickness rotating discs

Stresses reductions and/or raising the load-carrying capacity for a mechanical structure are always great dilemmas for researchers. In this article, a novel method is proposed, and its efficiency is examined for achieving these goals on functionally graded rotating nonuniform thickness discs. The originality of this method relies on comprising a geometrically well-defined area, into the whole structure, with certain homogeneous properties including density, thermal expansion coefficient, and elasticity matrix. This area acts as a reducer of the maximum values of various stress components. The solution of the magnetoelastic/magneto-thermoelastic problem is accomplished using the finite element method. The disc is subjected to partial uniform outer pressure, whereas, upon applying thermal loads; the thermal boundary conditions are considered symmetric. The proposed method is found to be beneficial as the obtained results demonstrated the ability to reduce the maximum stresses with different percentages depending on the location, angular width, and properties of the predefined area. This is reflected by an attainable decrease in the maximum compressive tangential stress and the von Mises stress by approximately 20.7% and 12.5%, respectively, under certain conditions.

Functionally graded materials (FGMs) have gained significant attention in the engineering community for their exceptional properties compared to traditional laminated composites and are widely used in severe working conditions 1,2 .Rotating discs, which are crucial components in many engineering systems, are often made of FGMs.Numerous studies have explored their elastic 3,4 , thermoelastic 5,6 , and elastoplastic static behaviors [7][8][9] .The dynamic performance of rotating discs has also been investigated in several studies [10][11][12] .Additionally, the effects of asymmetric loading conditions on their performance has been examined in some research 13,14 .
On the other hand, other scholars considered the magnetoelastic (ME) and magneto-thermoelastic (MTE) behaviors of discs (and their counterparts: cylinders and spheres), which are core items of many applications, such as aerospace industry, supersonic airplanes, rockets/missiles, submarine structures, nuclear energy, electronics, biomedical sector, geophysics, turbines/pumps/compressors rotors, optics, pipes, brake discs, internal combustion engines, ME sensors and actuators [15][16][17][18] .For example, Rad and Shariyat 19 developed a 3D model for FGM discs prone to nonuniform loadings.Additionally, different solution techniques were proposed to obtain the MTE performance of symmetric structures (disc/cylinder/sphere) 20,21 .These methods include for instance the finite element method (FEM) that was used by Zenkour and Abbas 16 to study the effects of the FGM's heterogeneity index on the disc's stresses.Moreover, Zenkour 22 and Dini et al. 18 obtained the steady MTE behaviors of multi-layer discs.
In the same vein, different means are sought to alleviate the stresses and/or enhance the sustainability (e.g., loading limits and failure likelihood) of any structure.For example, performing optimization for some parameters is found to be advantageous [23][24][25] .Also, the concept of grading the material properties along two simultaneous directions has its positive gain [26][27][28] .In addition, controlling a material's microstructure, to have a new substance with certain aspired properties 29 (e.g., yield strength 30 ), is also beneficial.The latter seems to be complicated; however, with the grand advancement in technology and fabrication techniques, it became attainable and easily accomplished.For instance, two methods are being used for this sake: modified slip-casting method 31 and fused deposition modeling method 32 .
In light of the extensive literature, efforts are being exerted on exploring methods to mitigate the stress and increase the loading capacity of rotating discs.Thus, this article proposes a novel and atypical approach to material customization under asymmetric loading conditions.This unconventional method involves incorporating a predefined area within the disc's FGM domain that has constant values for selected properties.The primary objective of this novel area is to reduce induced stress levels.The article explores through the FEM the effects of this area, including its properties, angular width, and position, on the disc's MT and MTE behaviors.The judgment on the success or failure of the proposed idea will be evaluated based on the reduction of the maximum stress components including the von Mises stress, to avoid plasticity occurrence.

Problem formulation
In polar coordinates ( rθz ), Fig. 1 depicts a disc with inner and outer radii r i and r o , respectively.The disc has a nonuniform thickness h = h(r) , and rotates with angular speed ω .It is comprised of two portions with distinct properties.One portion is formed of FGM (i.e., ψ > θ > ψ + φ ), while the other, named as the prescribed area ( ψ ≤ θ ≤ ψ + φ ), is made of a material with selected constant properties.The angles ψ and φ will be defined in the following section, and θ P is the angle of the applied outer pressure.
The thermomechanical equilibrium equations are written as 33 : where the comma denotes partial differentiation, T is the temperature, and k = k r 0 0 k θ is the thermal conductivity.F r = ρhω 2 r 2 is the centrifugal force, and F θ = ρhω ,t r 2 is the tangential force that exists when ω ,t = 0 , where t is the time variable, and ρ represents the material's density.In terms of h , it is proposed to change with the radial position r through the following relation 34 : with H is an imaginary thickness at r = 0 , and n j ( j = 1, 2 ) is a geometrical parameter.
Furthermore, since the disc has a large diameter-to-thickness ratio; the state of plane stress becomes a valid hypothesis; hence, the constitutive relation is 35 : where σ is the stress field of three components: σ rr , σ θθ , and σ rθ that are the radial, tangential, and shear stresses, respectively.Also, C is the elasticity matrix with the following non-zero components:

2
, where E and v resemble the material's elastic modulus and Poisson's ratio, respectively.Besides, ε refers to the strain components given through the kinematics relation (Eq.( 6)), with u and ϑ are the radial and (1) Schematic drawing of the proposed rotating annular variable thickness disc.The gradient-filled brown area represents the FGM, while the gray area, with an angular width of φ , is made of homogeneous material with constant selected properties.The outer pressure angle is denoted by θ P .

Material gradation
In the current analyses, the disc is made of two parts with dissimilar material compositions.The first part is made of FGMs (the gradient-filled brown area shown in Fig. 1).The material properties are radially graded as follows 38 : where β denotes a generic material property, and j refers to the heterogeneity index.Also, the two subscripts m and c stand for metal and ceramic, respectively.Conversely, the other part, or area, of the disc, which has an angular width of φ and named as the prescribed area (shaded in gray in Fig. 1), is located at angle ψ .This area is assumed to be perfectly bonded to the FGM area.It is introduced to be examined as being effective as an unconventional method in reducing the disc's induced stresses.In addition, it has special material characteristics compared to the FGM area.It has specific and constant selected properties ( β * ).The rest of the properties, other than the selected ones, would follow Eq. ( 9).In other words, this area has: where x is named as the property fraction, which is a positive number ( 0.25 , 0.5 or 1 ), and ξ β * is a symbol stand- ing for either m or c .For simplicity, β * would include either of ρ , α and C .In practical applications, this idea can be achieved through controlling the microstructure of the materials that can be conducted by the material's processing approach with the aid of the great advancement and technology available globally.

Finite element formulation
A finite element (FE) algorithm is built by the authors through the MATLAB software to solve for X = U T , where U = u ϑ .In FEM, X is connected to its nodal values through the shape functions ( N) 39,40 : where n n is the number of nodes per element, and N e i is read as: N at the i th node of the element e.Thereafter, the FE symbolic equation ( KX = R ) is derived using the standard Galerkin's procedures 39 .In that equation, K represents the global stiffness matrix, and R denotes the external force vector.That symbolic equation is expanded as below 39 : with n e is the total number of elements.Each term in Eq. ( 12) is calculated as follows: where  The integrations of Eq. ( 13) are accomplished through the gauss quadrature method with nine gauss points within each element that has an area of e , such that d� e = drdθ .In addition, σ n identifies the traction on a certain part Ŵ of the boundary 39,41 .Furthermore, B U represents the strain-displacement matrix, B T is the gradient matrix (see for instance Refs 39,42 .for their definitions.),and B µ is related to the magnetic field and is determined as follows: Finally, X can be easily obtained using matrix calculus.Then, ε and σ are computed directly after applying the proper boundary conditions 41 .

Algorithm verification
In this section, the validity of the proposed FEM formulation is examined through regenerating the results of three examples that exist in previous literature.
The first example is a nonuniform thickness clamped-free disc studied using the Runge-Kutta's method by Hassani et al. 43 .It rotates with ω = 300rad/s , and is prone to a temperature field defined by Eq. ( 15).Geometri- cally, the disc's thickness is defined via the power-law given in Eq. ( 16) with r i = 0.1m and r o = 0.6m .Table 1 in Ref 43 .lists the numerical values of the material properties that were graded according to Eq. ( 9) with j = 2 .Since this problem is completely symmetric ( σ rθ = 0 ); there is no need to investigate the 360 • model of the disc.Only its quarter is modeled with proper symmetric boundary conditions (BCs) at θ = 0 • and 90 •41 .Figure 2 shows the distribution of stresses obtained by the current FEM solution and the method used by Hassani et al. 43 .
The second example is presented to verify the solution ability in obtaining the ME behaviors of a nonrotating cylinder with r i = 0.2r o .It was subjected to uniform magnetic field ( H z = 2.23 × 10 9 A/m ), with the magnetic permeability varying as µ = µ 0 × r −221 .The outer surface was free-of-stress, while at r = r i a uniform and ( 14) Table 1.Mechanical properties and boundary conditions of FGM cylinder 38 .† Eq. ( 9) is used to describe the variation of E and v with j = 1.

Boundary conditions Elastic modulus ( GPa)
Poisson's ratio www.nature.com/scientificreports/symmetric pressure with a value of P 0 existed; thus, symmetric BCs prevails.Also, the following quantities were used:

and
Ev (1+v)(1−2v) = 28.8GPa 21; so that, after some mathematical manipulations, it yields: However, in this example the plane strain conditions were used; therefore, this conversion is conducted to Eq. ( 5) 44 : Results are presented in a dimensionless form according to Eq. ( 19), which is used henceforth, and are depicted in Fig. 3.The third and final example is a nonrotating cylinder (plane strain) with r o = 3r i studied by Li and Liu 38 .Table 1 lists the numeric values of the mechanical properties, and the associated BCs with P 0 = 100MPa .This (17) E = 25.5 × r −2 GPa, v = 0.378 Comparison between the radial variation of σrr , σθθ , and σrθ calculated through the current FE algorithm and in Ref. 38 .σrr and σθθ are drawn at θ = 0 • , and σrθ is calculated at θ = 45 • .
Vol:.( 1234567890) www.nature.com/scientificreports/example is re-examined to benchmark the capability of solving models with asymmetric BCs. Figure 4 portrays the resulting dimensionless stresses.
It is observed from the previous three verification examples that the obtained outcomes using the present model are in coherence with the results published in the literature.Therefore, the current solution scheme is robust and can be used for predicting the behaviors of discs with different loads and various boundary conditions.

Results and discussion
As intended in this work, mitigation of the induced ME and MTE stresses of a nonuniform thickness FGM disc is sought.This is accomplished through fixing one property or more within the prescribed area of the disc's domain.On the other hand, the disc has inner and outer radii of 0.2m and 0.5m , respectively.The parameters of Eq. ( 4) are: H = 0.1m , n 1 = 0.415196 , and n 2 = 3 8 .It rotates with ω = 700rad/s ( ω ,t = 0 ; thus F θ = 0 ), and is subjected to an axial magnetic field of H z = 10 6 A/m.
Table 2 lists the material properties assigned in this work.Equation ( 9) labels the FGM gradation model ( j = 1 ), and Eq. ( 10) describes the selected properties variation within the predefined area.Also, φ is set to 45 • , and ψ has a value of 0 • , 90 • , or 180 • .
Concerning the FE scheme, the eight-node ( n n = 8 ) isoparametric 2D element is used to discretize the disc's domain.It is found that convergence occurs at n e = 20000 elements after examining different numbers of elements.
It is observed from Fig. 5a that σrr max and σrr min stood at ∼ 1.5 and −1 , respectively.The critical tensile σθθ is found to be almost 0.92 (Fig. 5b), while the maximum compressive value is negligible (henceforth max and min denote the greatest values of tensile and compressive stresses in turns).Also, due to the lack of shear traction on any surface; the upper and lower limits of σrθ (Fig. 5c) have similar value of approximately ± 0.27 .Therefore, it can be concluded that the role of σrr min , σθθ min and σrθ is negligible compared to σrr max and σθθ max .
In order to decrease the values of the extreme stresses, Eq. ( 10) is used along with Eq. ( 9) yielding an area with certain homogeneous material properties.Table 3 lists the ME stress values at different values of ψ , and x β * while β * included ρ only.
It is seen that the variation of ψ and x ρ led to an increase in the values of σrr max due to change in F r .For exam- ple, this growth reached 2.5% at ψ = 0 • , x ρ = 1 and ξ ρ = c .Alternatively, at ψ = 90 • , x ρ = 0.5 and ξ ρ = c , a marginal increase is detected ( < 0.2% ).In terms of σθθ max , a swing in behavior can be witnessed.At x ρ = 0.5 or 1, and ξ ρ = m , σθθ max dropped by around 5% at ψ = 180 • , while it increased to 1.074 at ψ = 0 • with x ρ = 0.5 and ξ ρ = m .Also, using ψ = 0 • would yield larger stresses compared to ψ = 90 • and 180 • under this state of loading.Additionally, it is noticed that using x ρ = 0.5 with ξ ρ = c would yield nearly analogous results to the case of using ξ ρ = m at x ρ = 1.Table 2. Material properties of the FGM constituents 16,18   Afterwards, seeking further reductions in the stresses and/or raising the load-carrying capacity; C is added to ρ to follow Eq.(10).Table 4 illustrates the stresses values at different values of ψ and ξ β * .
Overall, results show that there would be a possibility to reduce the stresses by higher percentages compared to the case presented in Table 3.It can be emphasized that there is a fluctuation in the stresses' readings.However, at ψ = 90 • with ξ β * = m , σθθ max can be reduced by ∼ 7.3% while σrr max witnesses an insignificant increase ( < 1% ).This percentage of decline is greater by ∼ 2% than the case where β * only included ρ .Thus, including both ρ and C in Eq. ( 10) is more beneficial for the disc's behaviors if compared to changing ρ only, and the idea itself can produce encouraging outcomes (reducing failure/crack probability and raising the loading capacity).This showcases the originality of the proposed unconventional method.
On the other hand, it is not logical to investigate finite values of ψ and x β * , with different material proper- ties, loads, and geometry.Therefore, it can be stated that performing an adequate optimization for the current problem's parameters is to produce the utmost stresses lessening levels, which is beyond the scope of the study.

Case (2): Magneto-thermoelastic (MTE) loading.
In this case, more complications are seen compared to the previous problem.This is achieved by considering the existence of thermal loading which resembles a more practical scenario.For that sake, the following BCs are used: T(r i , θ) = 100 • C , and T(r o , θ) = 500 • C with the previous ones listed in "Case (1): Magnetoelastic (ME) loading" Section.As seen in Fig. 6, σrr max stood at 2.6 , σθθ min = −4.289, and σrθ = ±0.27 .The last one would be neglected owing to its minor value 28 .The concentra- tion would be directed towards the first two ones even though the absolute value of σθθ min is much larger than σrr max .Nonetheless, it can be said that both are important for ceramics that have tensile strength far less than the corresponding compressive strength.Thus, it can be stated that even if results yielded moderate growth for σθθ max and some dilatation for σrr max , this can be pointed out to as a positive point.Also, in this case, σVM (Fig. 6d) is included and discussed.It is calculated according to Eq. ( 20) 8 , and its dimensionless form is σVM = σ VM /P 0 .Here, σVM max peaked at 3.87.It should be noted that σVM max is excluded from the previous case study since the objective is to prove the adequacy of the idea.
Similar to the previous case, in order to apply the novel idea for stress reduction; Eq. ( 10) is used along with Eq. ( 9), where β * would include: (1) ρ , (2) ρ and α , and (3) ρ , α and C .This is done in order to obtain the utmost applicable alleviations for the three stress values: σrr max , σθθ min and σVM max .
Starting with ρ only in Eq. ( 10), a diverse variety of outcomes generated by altering ψ and x ρ , and for brev- ity, selected cases are debated ( ψ = 180 • and 0 • ).At ψ = 180 • with x ρ = 1 and ξ ρ = c , a tiny reduction of around 3% occurred for σθθ min (Fig. 7b), whereas σrr max grew by about 17% according to Fig. 7a.Despite the fact that the increase is nearly six times the value of the decrease, σVM max is impacted moderately, as it declined by (20) Table 3. Stress values obtaibed at different ψ by varying ρ only according to Eq. (10).

Stress Reference value
Equation (10) ψ nearly 3.5% as seen in Fig. 7c.This value experienced nearly a doubling hitting 6.2% at ψ = 0 • with x ρ = 1 and ξ ρ = c (Fig. 8c), where σθθ min (Fig. 8b) dropped by nearly 5.1% in contrast to σrr max (Fig. 8a) that rose slightly by approximately 1.4% .On the contrary, results revealed that both σrr max and σθθ min increased by < 1% and 8.5% , respectively, at similar ψ while x ρ = 0.5 and ξ ρ = m .Accordingly, σVM max jumped by ∼ 9% .Moreover, another interesting finding should be stated which is, while β * only included ρ , results at ξ ρ = c are better than that at ξ ρ = m .This is traced back to ρ c > ρ m ; hence F r becomes larger, and this works in opposite direction to the outer pressure especially at ψ = 0 • .Based on these selected results, a maximum feasible decline of σVM max is nearly 6.2% while β * included ρ only.Thus, the novel proposed method within the manuscript play an essential role in alleviating the stresses that can in turns lead to a reduction in the failure likelihood and an increase in the loading limits.Moreover,  www.nature.com/scientificreports/seeking larger mitigation for the stresses; another trial is conducted through letting α join ρ in Eq. (10).Table 5 lists the variation percentage detected for both σrr max and σθθ min according to the authors' records.
As seen, all the percentages are positive with different levels at x β * = 1 .Conversely, at x β * = 0.5 , σrr max and σθθ min experienced some minor declines.Nevertheless, at the two values of x β * , σVM max experienced gigantic surges that are traced back to the interaction between the loads, θ P , material properties, and the geometry.
Another finding can be comprehended from Table 5 that is, at any ψ and x β * = 0.5 , using ξ β * = c yields smaller growths for σθθ min compared to using ξ β * = m .This concludes that using ξ α = c would yield enhanced results compared to ξ α = m .In other words, the reduction of the prescribed area α is beneficial for the disc's performance as it reduces the thermal strains, and accordingly, the thermal stresses.However, generally, it can be drawn that under the current conditions, permitting β * to include α along with ρ has deleterious consequences on the disc's performance.
Afterwards, as done in Case (1), Eq. ( 10) is extended to C as well as ρ and α so as to try achieving higher stress alleviation levels.Unfortunately, records disclose that both σrr max and σθθ min increase substantially by using any value of ψ at x β * = 1 .Hence, σVM max goes up as well.For instance, at ψ = 90 • with x β * = 1 and ξ β * = c (Fig. 9), σrr max , σθθ min and σVM max significantly soared by 87% , 23% and 36% , respectively.
In contrast, using x β = 0.5 produced enhanced reductions for the three stresses as outlined in Table 6.It is seen that there are advantageous falls in the values of σrr max and σθθ min leading to noteworthy plunges of σVM regardless of ψ .For example, at ψ = 0 • , σVM max dropped by approximately 8.9% and 9.9% at ξ β * = m and c , respectively, and these two values surpass the applicable reductions levels at ψ = 90 • and 180 • .These results again prove that there would be stresses' reductions and an increase for the loading capacity through using the method proposed within this study.
Finally, and opposing to the previous results where x β * has an equal value for the three properties; hav- ing different configurations of it are examined at only ψ = 0 • .For that sake, the following values were used: As portrayed in Fig. 10a, σrr max witnessed substantial increase ( ∼ 8% ), while in Figs.11a and 12a, it experi- enced slight declines ( < 2% ).Conversely, σθθ min fell by approximately 19.8% at x C = 1 (Fig. 10b), and continued to decrease reaching more or less than 20.7% , which can be considered as a great achievement, at x C = 0.5 and 0.25 as seen in Figs.11b and 12b, respectively.Thus, according to the readings of σrr max and σθθ min , the smaller the value of x C , the higher the stresses' reductions are attainable.However, this is not dominant for the values  of σ VM max .At x C = 1 (Fig. 10c), it dropped by a large value ( ∼ 10.8% ), and continued to declined hitting nearly 12.5% at x C = 0.5 (Fig. 11c).As x C shrinks to 0.25 , a lessening of 10.5% is experienced as seen in Fig. 12c.These results prove the valuable impacts of the idea of keeping some the properties constant through a selected area within the disc to reduce the failure likelihood and/or raise the working load capacity and safety.They also confirm that the problem in hand is nonlinear; thus, the behaviors are hardly predicted.Accordingly, performing an adequate optimization is central to obtain the most suitable values of x β * and ψ .Additionally, it opens new horizons through considering more properties within the prescribed area that can have different values of ψ other the one used within the current study.

Conclusion
The main contribution of this article is trying to open new horizons for lowering the magnetoelastic/magnetothermoelastic stresses and/or increasing the working limits of rotating FGM discs.This was sought through an atypical idea that presumes the existence of an area within the disc with definite dimensions, location, and particular homogeneous properties.The assumption of perfect bonding between this area and the FGM area was utilized in the analyses.The solution of the governing equations was performed through FEM.
Overall, results showed that some constructive stress reductions were achievable depending on the location and the chosen constant properties of the prescribed area.A decline of ∼ 7.3% in the maximum tensile tangential stress occurred for the ME loading case.This percentage jumped to nearly 20.7% at the MTE case, which led to mitigation of the von Mises stress by about 12.5% at certain composition for the homogenous area.
These results were indisputably beneficial for decreasing the crack formation/propagation likelihood and/ or raising the payload limits, which proves the originality of the proposed method.However, optimizing many parameters is recommended to obtain the utmost stresses suppression especially for such nonlinear problems.These parameters include, for instance, thickness, property and its fraction, location, and the dimension of the homogenous area. https://doi.org/10.1038/s41598-023-39878-wwww.nature.com/scientificreports/

Figure 3 .Figure 4 .
Figure 3. Dimensionless stresses distributions for a cylinder subjected to magnetic field.

Table 4 .
Stress values obtained at different ψ by varying ρ and C according to Eq.