A Numerical Procedure for Shakedown Analysis of Thick Cylindrical Vessels with Crossholes under Dual Cyclic Loadings

A modified numerical procedure for the shakedown analysis of structures under dual cyclic loadings, based on the Abdalla method, is proposed in this paper. Based on the proposed numerical procedure, the shakedown analysis of the thick cylindrical vessels with crossholes (TCVCs) under cyclic internal pressure and cyclic thermal loading was carried out. The effects of material parameters (elastic modulus and thermal expansion coefficient) and crosshole radius on the elastic shakedown limit of TCVCs are discussed and, finally, normalized and formularized. Furthermore, the obtained shakedown limit boundary formulation is compared with FEA results and is verified to evaluate the shakedown behavior of TCVCs under cyclic internal pressure and cyclic thermal loading.


Introduction
Pressure-bearing components in petrochemical and nuclear power plants are usually operated under cyclic loadings. Throughout their lifetime, pressure-bearing components may face the problem of reversed plasticity or ratcheting failure. Therefore, it is very important to establish an effective shakedown analysis method to prevent early failure due to reversed plasticity or ratcheting; the corresponding Bree loading zoning diagram is shown in Figure 1, where σ 0 denotes the mechanical stress, σ t is the thermal stress, and σ s is the yield stress [1].
With the increasing emphasis on shakedown analysis in engineering, the theoretical study of the classical upper and lower bound shakedown theorems has gradually matured, and the corresponding numerical algorithms for shakedown analysis have also been rapidly developed [1][2][3][4][5][6][7][8][9][10][11][12][13][14]. The researchers have proposed some shakedown evaluation methods, such as the cycle-by-cycle method (CBC) [15][16][17], elastic compensation method (ECM) [2][3][4], and nonlinear superposition methods (Muscat and Mackenzie method, min{P L , 2P e } method, and Abdalla method) [5][6][7][8][9][10][11]. The CBC method has higher accuracy in shakedown assessment, but it is very time-consuming, so it is mainly used for verification of other shakedown analysis methods. The ECM method calculates the shakedown limits through the elastic finite element iteration, and each iteration needs to adjust the elastic modulus of the elements to obtain the redistribution of stress. The ECM method can achieve high accuracy for simple structures, but it often has a large calculation error for complex structures [18,19]. The Muscat and Mackenzie method can only be used for shakedown analysis of structures under mechanical loadings, while the min{P L , 2P e } method just supports the proportional loadings. It is necessary to find a simple and accurate method applicable to structures subjected to thermomechanical loadings as shown in Table 1. Furthermore, in the existing studies  considerable research has been undertaken on the shakedown of pressure-bearing structures under constant pressure and cyclic thermal loading, and those considering the shakedown behavior of TCVC and other structures under the cyclic pressure and cyclic thermal loading simultaneously (dual cyclic loadings) have not been reported. The establishment of a numerical flow for calculating the shakedown analysis of structures under dual cyclic loadings is crucial.
Materials 2023, 16, x FOR PEER REVIEW 2 of 12 can achieve high accuracy for simple structures, but it often has a large calculation error for complex structures [18,19]. The Muscat and Mackenzie method can only be used for shakedown analysis of structures under mechanical loadings, while the min{PL, 2Pe} method just supports the proportional loadings. It is necessary to find a simple and accurate method applicable to structures subjected to thermomechanical loadings as shown in Table 1. Furthermore, in the existing studies  considerable research has been undertaken on the shakedown of pressure-bearing structures under constant pressure and cyclic thermal loading, and those considering the shakedown behavior of TCVC and other structures under the cyclic pressure and cyclic thermal loading simultaneously (dual cyclic loadings) have not been reported. The establishment of a numerical flow for calculating the shakedown analysis of structures under dual cyclic loadings is crucial.  Thick cylindrical vessels are one kind of important pressure-bearing component in the process equipment. Meanwhile, thick cylindrical vessels often need crossholes to ensure equipment maintenance, material transfer, and connection between devices [25,26]. During service, thick cylindrical vessels with crossholes are often subjected to simultaneous cyclic pressure and cyclic temperature changes due to equipment shutdowns, startups, and peak regulation. In addition, the introduction of crossholes changes the original stress distribution of the structures, which causes significant stress concentration [27]. The structures are then prone to losing their initial shakedown state and may undergo reversed plasticity or ratcheting failure. Therefore, it is vital to carry out research on safety analysis of thick cylindrical vessels with crossholes (TCVCs) under dual cyclic loadings.
In summary, there is an urgent need for the development of a shakedown analysis procedure for the TCVC under dual cyclic loadings in which the influencing factors of material properties and structural dimensions are considered. In this paper, the obtained shakedown limits of four classical methods, namely, the min{PL, 2Pe} method, the Abdalla method, the stress analysis design method (SAD), and the cycle-by-cycle method (CBC), are first introduced and a precision comparison is carried out. According to the comparison results, a modified numerical procedure for the shakedown analysis of structures  Thick cylindrical vessels are one kind of important pressure-bearing component in the process equipment. Meanwhile, thick cylindrical vessels often need crossholes to ensure equipment maintenance, material transfer, and connection between devices [25,26]. During service, thick cylindrical vessels with crossholes are often subjected to simultaneous cyclic pressure and cyclic temperature changes due to equipment shutdowns, startups, and peak regulation. In addition, the introduction of crossholes changes the original stress distribution of the structures, which causes significant stress concentration [27]. The structures are then prone to losing their initial shakedown state and may undergo reversed plasticity or ratcheting failure. Therefore, it is vital to carry out research on safety analysis of thick cylindrical vessels with crossholes (TCVCs) under dual cyclic loadings.
In summary, there is an urgent need for the development of a shakedown analysis procedure for the TCVC under dual cyclic loadings in which the influencing factors of material properties and structural dimensions are considered. In this paper, the obtained shakedown limits of four classical methods, namely, the min{P L , 2P e } method, the Abdalla method, the stress analysis design method (SAD), and the cycle-by-cycle method (CBC), are first introduced and a precision comparison is carried out. According to the comparison results, a modified numerical procedure for the shakedown analysis of structures under dual cyclic loadings is established based on the Abdalla method by making a distinction between constant loadings and cyclic loadings and introducing the loading ratios. Then, considering the effects of material parameters (elastic modulus and thermal expansion coefficient) and crosshole radius, the shakedown limits of the TCVC under dual cyclic loadings are studied and formularized based on the proposed numerical procedure, which was verified to be effective. The research results have important reference and guiding significance for studying the shakedown behavior of structures under dual cyclic loadings.

Cycle-by-Cycle Method
The cycle-by-cycle method is the most fundamental technique for determining the shakedown limit of structures, and it is frequently employed to examine the effectiveness of other shakedown analysis methods [15]. Using the cycle-by-cycle method, the finite element method is used to simulate the stress-strain behavior of structures under cyclic loadings. The shakedown behavior of structures is judged according to the convergence of accumulated plastic strain, and the approximate solution of the shakedown limit load is obtained by the loading approximation method. Using the cycle-by-cycle method, Zheng et al. investigated the shakedown limit of thick cylinders with radial openings subjected to thermomechanical loadings [16]. Camilleri et al. investigated the shakedown and ratcheting behavior of a thin cylinder, a thick cylinder, and a thick cylinder with a radial crosshole [17].
The cycle-by-cycle method is applicable to structures under non-proportional loadings. It requires a considerable number of cycles and a series of loading combinations to find the shakedown limit load, which greatly affects the efficiency of the calculation. Therefore, this method is often used to check the accuracy of other methods.

min{P L , 2P e } Method
For the min {P L , 2P e } method, the elastic shakedown limit of a structure under proportional loadings can be further simplified to the lesser of the limit load and twice the elastic limit load. The limit load P L of the structure is determined according to the elastic-perfectly plastic analysis, while the elastic limit P e of the structure is determined according to the yield strength σ s of the material and the maximum elastic equivalent stress |σ ei | max of the structure under any proportional loading P i : The elastic shakedown limit P s is: The min {P L , 2P e } method is convenient in calculation, but it can only be used to determine the shakedown limit of structures under proportional loadings.

Abdalla Method
Abdalla split the operating loading under non-proportional loading into two parts: constant loading and cyclic loading, and proposed a simplified method to determining the shakedown limit of the structures [9]. The simplified method utilizes small displacement formulation and determines the shakedown limit load by performing elastic analysis and elastic-plastic analysis. In the elastic analysis, only the cyclic loading type is applied individually to obtain the elastic stress field σ E of the structures. In the elastic-plastic analysis, both the constant and cyclic loading are applied in two consecutive analysis steps to obtain the elastic-plastic stress field σ EPi of the structures. Finally, the two stress fields are superimposed to obtain the maximum residual stress field satisfying the Melan shakedown theorem: Materials 2023, 16, 3364 4 of 11 The shakedown limit load is the loading that corresponds to the residual stress field. Abdalla et al. applied this method to solve two benchmark shakedown problems, namely, the two-bar structure subjected to constant axial force and cyclic thermal loading, and the Bree cylinder subjected to constant internal pressure and cyclic high-temperature variation across its wall [10]. The advantages of the Abdalla method are high accuracy and efficiency, while the method is also applicable to structures under non-proportional loadings.

Stress Analysis Design Method
In the ASME Boiler and Pressure Vessel Code Sec. VIII Division 2, the stress analysis design method is supplied [28]. Using this method, the maximum shear stress theory is adopted, and the stress is divided into primary stress, secondary stress, and peak stress based on the location of stresses, loading conditions, stress properties, and other factors. Also in the stress analysis design method, the stress intensity checks are classified into five categories according to the stress conditions in different paths of the pressure vessel, which are listed as follows: where P m is the primary membrane stress, P L means the local primary membrane stress, P b is the primary bending stress, Q is the secondary stress, F is the peak stress, and S a is the fatigue strength. The stress analysis design method can be used to quickly judge whether the common component shakes down or not; however, in complex engineering practice, there are still disputes in how to classify the stresses. In summary, the Abdalla method has higher accuracy and generality. In this study, the Abdalla method is further developed for dual cyclic loadings in the proposed shakedown analysis procedure of the TCVC, as introduced in the following section.

The Modified Shakedown Analysis Procedure for Dual Cyclic Loadings
Since most of the existing shakedown studies just consider the pressure-bearing structures under constant pressure and cyclic thermal loading, in this study a modified numerical procedure for the shakedown analysis of engineering structures subjected to dual cyclic loadings was developed. In the proposed shakedown analysis procedure, the Abdalla method was further developed by making a distinction between constant loadings and cyclic loadings, while a concept of loading ratio was also introduced. Furthermore, two cyclic loading types were applied along each loading path to calculate the residual stress field. The proposed procedure can be divided into the following five steps:

1.
Establish an initial polygonal loading domain, in which the domain corners are corresponding to different loading ratios, as shown in Figure 2.

2.
Apply two cyclic loadings on the structure monotonically along each loading path determined by a certain loading ratio, and calculate the elastic stress field of the structure. 3.
Increase the maximum loading (T 0 + ∆T and P 0 + ∆P) incrementally along the loading path in N steps, and calculate the elastic-plastic stress field of the structure at each step.

4.
Then assume that the total stress field consists of the elastic stress field and residual stress field. Calculate the residual stress field at each step by removing the elastic stress field, and output the maximum loading value of the residual stress field not exceeding the yield strength. Then the corresponding elastic shakedown limit load under the loading ratio is determined.

5.
Change the loading path for different loading ratios and repeat steps 2 to 4 in order to determine the elastic shakedown loading region of the structure under dual cyclic loadings as the orange region shown in Figure 2.
determined by a certain loading ratio, and calculate the elastic stress field of the structure. 3. Increase the maximum loading (T0 + ∆T and P0 + ∆P) incrementally along the loading path in N steps, and calculate the elastic-plastic stress field of the structure at each step. 4. Then assume that the total stress field consists of the elastic stress field and residual stress field. Calculate the residual stress field at each step by removing the elastic stress field, and output the maximum loading value of the residual stress field not exceeding the yield strength. Then the corresponding elastic shakedown limit load under the loading ratio is determined. 5. Change the loading path for different loading ratios and repeat steps 2 to 4 in order to determine the elastic shakedown loading region of the structure under dual cyclic loadings as the orange region shown in Figure 2. In this study, a square loading domain was established according to the shakedown limit loads under single cyclic internal pressure and single cyclic thermal loading. The maximum combined load is taken as the value of the coordinates on the square boundary. The accuracy of the obtained results can be only affected by the increment of ∆P and ∆T at each step. In addition, it needs to be noted that the maximum combined load should exceed the yield strength of the structure.

Finite Element Model and the Cyclic Loadings
In this study, 316 stainless steel thick cylindrical vessels with crossholes (TCVCs) were considered for the shakedown analysis. The elastic-perfectly plastic constitutive model was utilized to simulate the shakedown behavior of the TCVC subjected to dual cyclic loadings. The material property parameters of 316 stainless steel at different temperatures are given in Table 2. The initial geometric dimensions of the TCVC are shown in Figure 3, for which the inner radius Ri = 300 mm, outer radius Ro = 450 mm, half length of the whole model L = 800 mm, and crosshole radius r = 60 mm. Considering the structure geometry and loading symmetry, a quarter of the TCVC is adopted for modeling, as shown in Figure 4. For the quarter finite element model, the symmetry constraints are In this study, a square loading domain was established according to the shakedown limit loads under single cyclic internal pressure and single cyclic thermal loading. The maximum combined load is taken as the value of the coordinates on the square boundary. The accuracy of the obtained results can be only affected by the increment of ∆P and ∆T at each step. In addition, it needs to be noted that the maximum combined load should exceed the yield strength of the structure.

Finite Element Model and the Cyclic Loadings
In this study, 316 stainless steel thick cylindrical vessels with crossholes (TCVCs) were considered for the shakedown analysis. The elastic-perfectly plastic constitutive model was utilized to simulate the shakedown behavior of the TCVC subjected to dual cyclic loadings. The material property parameters of 316 stainless steel at different temperatures are given in Table 2. The initial geometric dimensions of the TCVC are shown in Figure 3, for which the inner radius R i = 300 mm, outer radius R o = 450 mm, half length of the whole model L = 800 mm, and crosshole radius r = 60 mm. Considering the structure geometry and loading symmetry, a quarter of the TCVC is adopted for modeling, as shown in Figure 4. For the quarter finite element model, the symmetry constraints are applied on the planes of X = 0 and Z = 0, while the coupling displacement is applied on the end plane of the TCVC with Z = L in the Z direction.
The loading conditions of the finite element model primarily include a cyclic internal pressure loading and a cyclic thermal loading, as shown in Figure 5. The cyclic pressure varies from P 0 to P 0 + ∆P, and is applied as the inner pressure of the TCVC with an initial internal pressure of P 0 = 0 MPa. When the environmental pressure remains P 0 , ∆P means the maximum pressure difference between the inner surface and the outer surface of the TCVC. Then the cyclic thermal loading is distributed linearly along the thickness with a constant outer surface temperature of T 0 and a cyclic inner surface temperature of T 0 + ∆T (T 0 = 300 • C), in which ∆T denotes the maximum temperature difference between the inner surface and the outer surface of the TCVC. Hence, the elastic shakedown limit load can be characterized jointly by the maximum pressure difference ∆P and temperature difference ∆T. Furthermore, it should be mentioned that for the condition of TCVC subjected to cyclic thermomechanical loadings, Formula (8) must be used to limit the sum of the primary plus secondary stress to within two times the yield strength in order to ensure the shakedown of the structure. applied on the planes of X = 0 and Z = 0, while the coupling displacement is applied o the end plane of the TCVC with Z = L in the Z direction.   The loading conditions of the finite element model primarily include a cyclic interna pressure loading and a cyclic thermal loading, as shown in Figure 5. The cyclic pressur varies from P0 to P0 + ∆P, and is applied as the inner pressure of the TCVC with an initia internal pressure of P0 = 0 MPa. When the environmental pressure remains P0, ∆P mean the maximum pressure difference between the inner surface and the outer surface of th TCVC. Then the cyclic thermal loading is distributed linearly along the thickness with constant outer surface temperature of T0 and a cyclic inner surface temperature of T0 + ∆T (T0 = 300 °C), in which ∆T denotes the maximum temperature difference between the inne surface and the outer surface of the TCVC. Hence, the elastic shakedown limit load can b characterized jointly by the maximum pressure difference ∆P and temperature differenc ∆T. Furthermore, it should be mentioned that for the condition of TCVC subjected to cycli thermomechanical loadings, Formula (8) must be used to limit the sum of the primar plus secondary stress to within two times the yield strength in order to ensure the shake down of the structure. applied on the planes of X = 0 and Z = 0, while the coupling displacement is applied on the end plane of the TCVC with Z = L in the Z direction.   The loading conditions of the finite element model primarily include a cyclic internal pressure loading and a cyclic thermal loading, as shown in Figure 5. The cyclic pressure varies from P0 to P0 + ∆P, and is applied as the inner pressure of the TCVC with an initial internal pressure of P0 = 0 MPa. When the environmental pressure remains P0, ∆P means the maximum pressure difference between the inner surface and the outer surface of the TCVC. Then the cyclic thermal loading is distributed linearly along the thickness with a constant outer surface temperature of T0 and a cyclic inner surface temperature of T0 + ∆T (T0 = 300 °C), in which ∆T denotes the maximum temperature difference between the inner surface and the outer surface of the TCVC. Hence, the elastic shakedown limit load can be characterized jointly by the maximum pressure difference ∆P and temperature difference ∆T. Furthermore, it should be mentioned that for the condition of TCVC subjected to cyclic thermomechanical loadings, Formula (8) must be used to limit the sum of the primary plus secondary stress to within two times the yield strength in order to ensure the shakedown of the structure.

Method Comparison
In order to compare the effectiveness of different shakedown analysis methods, the shakedown analysis of the TCVC subjected to only the cyclic internal pressure was first carried out, as shown in this section. Figure 6 shows the elastic shakedown limits for the

Method Comparison
In order to compare the effectiveness of different shakedown analysis methods, the shakedown analysis of the TCVC subjected to only the cyclic internal pressure was first carried out, as shown in this section. Figure 6 shows the elastic shakedown limits for the structure with different crosshole radii r obtained by the min{P L , 2P e } method, Abdalla method, stress analysis design method (SAD) and cycle-by-cycle method (CBC). With the comparison with the results from the cycle-by-cycle method, the most fundamental and accurate method, the most accurate to least accurate are shown to be: min{P L , 2P e } method, Abdalla method, and stress analysis design method. In comparison with the stress analysis design method, the results obtained by the min{P L , 2P e } method and Abdalla method are more accurate. Meanwhile, compared with the min{P L , 2P e } method, the Abdalla method can be used to further analyze the shakedown of structures under nonproportional loadings, which has contributed to its widespread adoption.

Method Comparison
In order to compare the effectiveness of different shakedown analysis methods, the shakedown analysis of the TCVC subjected to only the cyclic internal pressure was first carried out, as shown in this section. Figure 6 shows the elastic shakedown limits for the structure with different crosshole radii r obtained by the min{PL, 2Pe} method, Abdalla method, stress analysis design method (SAD) and cycle-by-cycle method (CBC). With the comparison with the results from the cycle-by-cycle method, the most fundamental and accurate method, the most accurate to least accurate are shown to be: min{PL, 2Pe} method, Abdalla method, and stress analysis design method. In comparison with the stress analysis design method, the results obtained by the min{PL, 2Pe} method and Abdalla method are more accurate. Meanwhile, compared with the min{PL, 2Pe} method, the Abdalla method can be used to further analyze the shakedown of structures under non-proportional loadings, which has contributed to its widespread adoption.

The Shakedown Analysis Considering Dual Cyclic Loadings
Based on the proposed shakedown analysis flow, the elastic shakedown limit of the TCVC with different material parameters and crosshole radius under dual cyclic loadings are investigated by taking the cyclic thermal loading into account. Firstly, the elastic modulus and thermal expansion coefficient were multiplied by different coefficients Ec and αc on the original basis to study the effects of material parameters on the elastic shakedown limit of the TCVC. In the simulations, 11 loading ratios were adopted. The simulated results are shown in Figure 7. It can be found that when the yield strength is constant, the elastic shakedown limit load decreases with the increase in elastic modulus and thermal expansion coefficient. In order to conduct an in-depth study on the shakedown of the TCVC, the dimensionless treatment was carried out for

The Shakedown Analysis Considering Dual Cyclic Loadings
Based on the proposed shakedown analysis flow, the elastic shakedown limit of the TCVC with different material parameters and crosshole radius under dual cyclic loadings are investigated by taking the cyclic thermal loading into account. Firstly, the elastic modulus and thermal expansion coefficient were multiplied by different coefficients E c and α c on the original basis to study the effects of material parameters on the elastic shakedown limit of the TCVC. In the simulations, 11 loading ratios were adopted. The simulated results are shown in Figure 7. It can be found that when the yield strength is constant, the elastic shakedown limit load decreases with the increase in elastic modulus and thermal expansion coefficient. In order to conduct an in-depth study on the shakedown of the TCVC, the dimensionless treatment was carried out for Figure 7. In the dimensionless treatment, the maximum pressure difference ∆P is divided by the limit pressure P L of the TCVC as the vertical coordinate, and the elastic thermal loading σ t = Eα∆T/2(1 − v) divided by the yield stress σ s replaces ∆T as the horizontal coordinate. Figure 8 shows the results of dimensionless treatment. It can be seen from the figures that the elastic shakedown limits of the TCVC with different material parameters almost coincide, and a general rule of shakedown for the TCVC at a certain size is obtained.
Based on the research above, the effect of different crosshole radii (r/R i = 0.1, 0.2, 0.3, 0.4) on the elastic shakedown limit of the TCVC with a constant inner radius is analyzed. As an illustration, the model for r/R i = 0.2 in Figure 4 is used and the crosshole radius changes to 30 mm, 90 mm, and 120 mm. In order to establish the shakedown region in Bree's diagram, the elastic shakedown limit of the TCVC with 13 loading ratios is calculated based on the proposed numerical flow. The simulated results are shown in Figure 9. The elastic shakedown limit can be approximately expressed by piecewise linear functions. It can be seen from Figure 9 that when σ t /σ s is small, the elastic shakedown limit decreases as the crosshole radius increases; however, when σ t /σ s reaches a specific value, the crosshole radius has little effect on the elastic shakedown limit. treatment, the maximum pressure difference ∆P is divided by the limit pressure PL of the TCVC as the vertical coordinate, and the elastic thermal loading divided by the yield stress σs replaces ∆T as the horizontal coordinate. Figure 8 shows the results of dimensionless treatment. It can be seen from the figures that the elastic shakedown limits of the TCVC with different material parameters almost coincide, and a general rule of shakedown for the TCVC at a certain size is obtained.  Based on the research above, the effect of different crosshole radii (r/Ri = 0.1, 0.2, 0.3, 0.4) on the elastic shakedown limit of the TCVC with a constant inner radius is analyzed. As an illustration, the model for r/Ri = 0.2 in Figure 4 is used and the crosshole radius changes to 30 mm, 90 mm, and 120 mm. In order to establish the shakedown region in Bree's diagram, the elastic shakedown limit of the TCVC with 13 loading ratios is calculated based on the proposed numerical flow. The simulated results are shown in Figure 9. The elastic shakedown limit can be approximately expressed by piecewise linear functions. It can be seen from Figure 9 that when σt/σs is small, the elastic shakedown limit decreases as the crosshole radius increases; however, when σt/σs reaches a specific value, the crosshole radius has little effect on the elastic shakedown limit. treatment, the maximum pressure difference ∆P is divided by the limit pressure PL of the TCVC as the vertical coordinate, and the elastic thermal loading divided by the yield stress σs replaces ∆T as the horizontal coordinate. Figure 8 shows the results of dimensionless treatment. It can be seen from the figures that the elastic shakedown limits of the TCVC with different material parameters almost coincide, and a general rule of shakedown for the TCVC at a certain size is obtained.  Based on the research above, the effect of different crosshole radii (r/Ri = 0.1, 0.2, 0.3, 0.4) on the elastic shakedown limit of the TCVC with a constant inner radius is analyzed. As an illustration, the model for r/Ri = 0.2 in Figure 4 is used and the crosshole radius changes to 30 mm, 90 mm, and 120 mm. In order to establish the shakedown region in Bree's diagram, the elastic shakedown limit of the TCVC with 13 loading ratios is calculated based on the proposed numerical flow. The simulated results are shown in Figure 9. The elastic shakedown limit can be approximately expressed by piecewise linear functions. It can be seen from Figure 9 that when σt/σs is small, the elastic shakedown limit decreases as the crosshole radius increases; however, when σt/σs reaches a specific value, the crosshole radius has little effect on the elastic shakedown limit. In summary, after the normalization of Bree's diagram considering the material parameters, the elastic shakedown limit of TCVC in Figure 9 still receives a dispersion of the crosshole radius influence. In order to obtain a more flexible shakedown limit description, we performed a segmented linear fitting to the shakedown boundary of the TCVC for different crosshole radii in Figure 9, and the resulting equation is listed as following: In summary, after the normalization of Bree's diagram considering the material parameters, the elastic shakedown limit of TCVC in Figure 9 still receives a dispersion of the crosshole radius influence. In order to obtain a more flexible shakedown limit description, we performed a segmented linear fitting to the shakedown boundary of the TCVC for different crosshole radii in Figure 9, and the resulting equation is listed as following: Further, it can also be seen from Figure 9 that the shakedown limit value corresponding to σ t /σ s = 0 and the slope of the linear function before the turning point are related to the crosshole radius ratio r/R i , which can be fitted as a formula as: Therefore, the uniform description of the elastic shakedown limit of the TCVC considering different crosshole radii can be expressed as: In order to validate the reliability of the above fitting formula, the shakedown limit boundary obtained from Formula (12) is compared with the numerical solution result at r/R i = 0.35 simulated with the finite element method, as shown in Figure 10. It can be seen from Figure 10 that the calculated results of the formulation match well with the finite element analysis results obtained from the proposed procedure. Meanwhile, the shakedown limit boundary determined by the cycle-by-cycle (CBC) method is also shown in Figure 10, which agrees well with that obtained from the proposed procedure. Therefore, the proposed numerical procedure for the shakedown analysis of the TCVC under dual cyclic loadings established in this study is verified to be effective and reliable.  (12) In order to validate the reliability of the above fitting formula, the shakedown limit boundary obtained from Formula (12) is compared with the numerical solution result at r/Ri = 0.35 simulated with the finite element method, as shown in Figure 10. It can be seen from Figure 10 that the calculated results of the formulation match well with the finite element analysis results obtained from the proposed procedure. Meanwhile, the shakedown limit boundary determined by the cycle-by-cycle (CBC) method is also shown in Figure 10, which agrees well with that obtained from the proposed procedure. Therefore, the proposed numerical procedure for the shakedown analysis of the TCVC under dual cyclic loadings established in this study is verified to be effective and reliable.

Conclusions
In this work, a developed numerical procedure for the shakedown analysis of TCVCs under dual cyclic loadings is proposed, and the shakedown behavior of TCVCs under dual cyclic loadings is studied to verify the effectiveness and reliability of the proposed procedure. It is shown that:

Conclusions
In this work, a developed numerical procedure for the shakedown analysis of TCVCs under dual cyclic loadings is proposed, and the shakedown behavior of TCVCs under dual cyclic loadings is studied to verify the effectiveness and reliability of the proposed procedure. It is shown that: (1) Comparing different shakedown analysis methods for only the application of cyclic internal pressure, it is shown that the Abdalla method is more accurate and is also effective for non-proportional loading. (2) The elastic shakedown limit load of the TCVC decreases with the increase in the elastic modulus and the thermal expansion coefficient of the material. After dimensionless treatment, the influence of material parameters can be eliminated. (3) The modified shakedown evaluation method and procedure proposed in this work can be used to accurately evaluate the shakedown of TCVCs under cyclic internal pressure and cyclic thermal loading.