Acceleration techniques for the numerical simulation of the cyclic plasticity behaviour of mechanical components under thermal loads

Numerical simulations of components subjected to cyclic thermo-mechanical loads require an accurate modelling of their cyclic plasticity behaviour. Combined models permit to capture monotonic hardening as well as cyclic hardening/softening phenomena, that occur in reality. In principle the durability assessment of a component under thermal loads can be performed only if the cyclic behaviour is simulated until complete material stabilization. As materials stabilize approximately at half the number of cycles to failure, it follows that in case of small plastic strains a huge number of cycles must be considered and an unfeasible simulation time would be required. Accelerated models have thus been proposed in literature. The aim of this work is that of comparing the different acceleration techniques in the case a round mould for continuous casting loaded thermo-mechanically. It can be observed that the usual approach of using the stabilized stress-strain curve already from the first cycle could lead to relevant errors. An alternative method is that of increasing the value of the parameter that controls the speed of stabilization in the combined model. This approach permits the number of cycles to reach stabilization to be drastically reduced, without affecting the overall mechanical behaviour. Based on this approach, a simple design rule, that can be adopted, particularly when relatively small plastic strains occur, is finally proposed. 1 Numerical simulation of mechanical components under thermal loads The choice of a suitable cyclic plasticity model to be used in numerical simulation of components subjected to cyclic thermo-mechanical loads is a crucial aspect in design. Some of typical examples are discussed in [1, 2, 3]. In fact, to perform a reliable thermal fatigue analysis a model able to capture accurately the behaviour of the material observed during experimental testing is required. On the other hand, durability assessment of a component under thermal loads can be done only if complete material stabilization is reached [4]. As materials stabilize approximately at half number of cycles to failure, it follows that in case of small plastic strain a huge number of cycles must be considered and an unfeasible computational time would be required. In a recent work [5] dealing with a thermomechanical analysis of a squared mould showed that a huge number of cycles (around 60567) must be performed to reach the stabilized condition. Since the geometry of the squared mould requires a 3D finite element (FE) simulation, even taking into account symmetries and optimizing the mesh, ≈15 min is needed to compute 1 cycle or ≈630 day to obtain final results, according to the value (b≈4) of the speed of stabilization of the mould alloy, see Fig. 1. Therefore, to overcome this problem, some alternative models (accelerated, stabilized), proposed in [6], have been adopted and compared. It has been shown that accelerated models gave the most conservative results, as the lowest value of number of cycle to failure was obtained. On the other hand, in that work it was not possible to make a comparison with the combined (complete) model because a too long computational time was necessary to reach stabilization. Fig. 1. Possible strategy to reduce the computational time. Fig. 1 presents schematically the desirable diagram required to perform a correct design of the component. In principle the basic idea of the accelerated techniques ≈ 60567 (≈630 days) N (log) Sp ee d of st ab ili za tio n, b Combined model ≈ 4 Accelerated models N=100 (≈25 hours) UNREALISTIC Δεpl UNFEASIBLE 1500 A B MATEC Web of Conferences 165, 19010 (2018) https://doi.org/10.1051/matecconf/201816519010


Numerical simulation of mechanical components under thermal loads
The choice of a suitable cyclic plasticity model to be used in numerical simulation of components subjected to cyclic thermo-mechanical loads is a crucial aspect in design.Some of typical examples are discussed in [1,2,3].In fact, to perform a reliable thermal fatigue analysis a model able to capture accurately the behaviour of the material observed during experimental testing is required.On the other hand, durability assessment of a component under thermal loads can be done only if complete material stabilization is reached [4].As materials stabilize approximately at half number of cycles to failure, it follows that in case of small plastic strain a huge number of cycles must be considered and an unfeasible computational time would be required.
In a recent work [5] dealing with a thermomechanical analysis of a squared mould showed that a huge number of cycles (around 60567) must be performed to reach the stabilized condition.Since the geometry of the squared mould requires a 3D finite element (FE) simulation, even taking into account symmetries and optimizing the mesh, ≈15 min is needed to compute 1 cycle or ≈630 day to obtain final results, according to the value (b≈4) of the speed of stabilization of the mould alloy, see Fig. 1.Therefore, to overcome this problem, some alternative models (accelerated, stabilized), proposed in [6], have been adopted and compared.It has been shown that accelerated models gave the most conservative results, as the lowest value of number of cycle to failure was obtained.On the other hand, in that work it was not possible to make a comparison with the combined (complete) model because a too long computational time was necessary to reach stabilization.  is that of moving from point A to point B (i.e. from the purple "unfeasible" zone, characterized by too long computation time, to the white area, where this time is acceptable).This goal can be obtained by increasing a parameter of the material model, which governs the speed of stabilization.On the other hand, also too high speed of stabilization must be avoided, as it could lead to convergence problems and "unrealistic" results.
In order to perform a complete validation, a similar component with a simpler geometry (i.e.round mould) is considered in this work.Due to axi-symmetry, a plane model can be used thus permitting a faster simulation to be performed.The component behaviour can be followed till stabilization adopting the combined model.This case can be then taken as a reference to compare different accelerated models proposed in literature.The aim of this work is to validate the "design diagram" shown in Fig. 1, which permits a suitable speed of stabilization to be set-up, in order to reach stabilized conditions in a feasible computational time without affecting the consequent durability assessment.

Accelerated cyclic plasticity models
A combined (kinematic and isotropic) model captures simultaneously monotonic hardening (elasto-plastic) and cyclic hardening/softening behaviour of a material.The yield surface can be represented as [7]: where S is the deviatoric stress tensor, X is the back stress tensor, R is the drag stress and σ 0 is the initial yield stress.In this case, the yield surface both translates (controlled by X) and expands (controlled by R) with the plastic strain ε p .
Different kinematic models have been proposed in literature to relate X to ε pl [7].Chaboche model (nonlinear kinematic) assumes that the increment of the back stress dX is expressed as a function of the plastic strain increment dε pl and the accumulated plastic strain dε pl,acc [7]: where C is the hardening modulus and γ is the nonlinear recovery parameter that controls the decrease rate of the hardening modulus as plastic strain accumulates.The model with one pair of (C 1 , γ 1 ) is known as the Armstrong and Frederick model.In addition, considering γ=0 the Prager model (i.e.linear kinematic) is obtained and relation (2) can be expressed as [7]: The nonlinear isotropic model controls the homothetic expansion of the yield surface [7]: where R ∞ is the saturated drag stress, b is the parameter that controls the speed of hardening (R ∞ >0) or softening (R ∞ <0).Integration of (4) gives expression that links R ∞ to ε pl,acc [7]: Stabilized condition is obtained when R reaches R ∞ .Cyclic hardening/softening behaviour of a material is mainly governed by the speed of stabilization b and the accumulated plastic strain ε pl,acc .When the accumulated plastic strain is relatively small, a huge number of cycles is required to reach the stabilized condition as already shown in Fig. 1.Some accelerated methods have thus been proposed in literature.In particular, in presence of creep and thermal fatigue, authors such in [8,9,10] suggest only a limited number of cycles to be simulated; even if this procedure seems not well defined, it has to be considered that the presence of visco-elasticity generally strongly reduces the stabilization time.On the other hand, some authors [11,12] use the kinematic model with stabilized material condition from the beginning of simulation, at the same time neglecting the initial transient cyclic deformation.If creep constitutes the damage criteria, a more rigorous approach was proposed in [13], where an extrapolation technique is developed to speed up the simulation.An alternative approach was suggested in [6].Since the accumulated plastic strain (see Eq. 5) depends on loading condition and cannot be changed, it was proposed to accelerate material stabilization by increasing (10-200 times) the speed of stabilization b.

Case study: round mould under variable thermal flux
A mould is a part of the continuous casting process where solidification of steel starts.During solidification, a huge thermal flux, q, is transferred from the molten steel, that is in contact to the inner surface, to the water-cooled outer surface of the mould, see Fig. 2. The thermal flux distribution varies between two conditions: q=0 when the plant is turned-off due to maintenance and q=q max when the plant operates [5].

Fig. 2. Description of a mould and of its working condition.
mal Flux -q a r ϴ Generally, moulds are constituted by a hollow tube of copper alloys while according to the geometry of the final product, different types of cross sections (square, rectangular or rounded) can be adopted.Copper alloys are used because they provide the best combination of high thermal conductivity and adequate mechanical properties.
Temp.In this work, a round mould 1000 mm long, 16 mm thick and with 200 mm inner diameter was studied, similarly to the component described in [14].As material, it was assumed a CuAg0.1 alloy, whose properties, listed in Tab. 1, were experimentally obtained in [15,16,17].

Numerical simulation -FEM model
A two dimensional (2D) model was adopted due to axisymmetry, thus strongly decreasing the computational time.The finite element model (FEM) had 760 elements and 2487 nodes.As can be seen in Fig. 3, the mesh was refined in the meniscus area close to the liquid free surface where the maximum thermal gradients occur

Thermal analysis and results
A thermal analysis was performed considering the thermal flux acting on the inner surface and imposing convection on the outer surface of the mould to simulate a water cooling.In the FEM model 8-node elements were used.The thermal flux proposed in [14] was 50% increased to reach a maximum temperature close to 300 o C for which material parameters are available in [15].The convection coefficient is 48000 W/m 2 K and temperature of the cooling water is 40 °C.Variation of the thermal flux schematically presented in Fig. 2 was simulated by a sequence of steady state analyses.A nonlinear solution was carried out to simulate the temperature dependence of thermal properties.The highest temperatures and thermal gradients occur in the region close to meniscus.The maximum temperature 298 o C occurs in the position labelled with "A" (see Fig. 3b) that is also a representative point for the calculation presented in the following section.

Mechanical analysis and results
Temperature distribution calculated previously in the thermal analysis was the input data for the mechanical analysis.Since the component is free to expand, no mechanical constraints were imposed on the numerical model.The nonlinear mechanical analysis was performed assuming temperature dependence of material parameters and considering combined (nonlinear kinematic + nonlinear isotropic) model with parameters taken from [15].
The component is cyclically loaded until stabilization by using different material models.The comparison is performed considering the equivalent strain range ∆ε eq .Before proceeding with simulations, a criterion that defines when material stabilization actually occurs needs to be established.
Normalized maxima of the stress evolution as a function of the accumulated plastic strain are considered [7,18]: where σ max,i , σ max,1 σ max,s are the maximum stress for the N th number of cycle, first and stabilized, respectively.Assuming that the plastic strain range ∆ε pl is constant over cycles (see [7]): the following relation proposed in [18] could be introduced: Eq. ( 8) assumes that a value close to 5 is good saturation criterion of the exponent in Eq. ( 6).The criterion adopted to determine the stabilized condition of a material is:  Results obtained with the combined model were compared with those achieved considering accelerated models with 8 different values of parameter b a : 10b, 20b, 30b, 100b, 200b, 300b, 421b and 0.01b, covering a wider range with respect to that proposed in [6,18].Prager and nonlinear kinematic models with initial (E, σ 0 ) and stabilized (E s , σ 0* ) condition (where yield stress and Young's modulus were determined from 1 st and stabilized cycle, respectively) were also implemented for comparison.
Firstly, the combined model was analysed.Adopting the stabilization criterion given by Eq. ( 9), the material stabilizes within 537 cycles.This results can be argued also from Fig. 4 and 6.Fig. 4a presents hoop, axial, radial and von Mises maximum stresses (computed at the critical point A) versus the number of cycles.It can be observed that a biaxial state of stress occurs.In Fig. 4b, the three components of strain range and the equivalent strain range, determined according to [4], versus the number of cycles are presented.It can be noticed that the hoop strain range is almost constant, whereas the other two components increase until stabilization.This behaviour can be clearly observed also in Fig. 5, where hoop and axial stress-strain evolution is presented.For the sake of clarity, only the first five cycles, the 100 th , the 200 th , the 300 th , the 400 th and the final stabilized cycles are presented.
Results obtained with several accelerated models are under investigation hereinafter.According to [6], stabilization of material can be reached in a shorter simulation time by increasing the original value of speed of stabilization b to b a .On the other hand, in the same work [6] it was observed that in some cases, particularly when softening occurs, the accelerated solution can also be completely incorrect.
In this work a wide span of cases with different values of b a were thus considered (b a =10b, 20b, 30b, 100b, 200b, 300b).The lower bound case b a =0.01b was also investigated.In this case, stabilization is reached in 8 cycles.This is due to the criterion defined by Eq. ( 9).Theoretically, models with speed of stabilization that tends towards zero requires a number of cycles to stabilization that tends towards infinite.Nevertheless, as in this case the contribution of the nonlinear isotropic model is negligible, cyclic decrease of stress from cycle to cycle is smaller than tolerance and therefore a fast stabilization seems to occur.The upper bound case (b a =421b) was also investigated.
Fig. 6 shows a comparison between these cases in term of the maximum von Mises stress in each cycle.The number of cycles of each curve are normalized with respect to the corresponding number of cycles to stabilization N stab in order to emphasize differences among the models occurring in first few cycles.Despite the different behaviour at the beginning, all the accelerated models reach almost the same value of von Mises stress.Obviously in the case b a =0.01b a higher value of Von Mises stress is reached.In fact, as already explained, softening does not occur.The evolution of the equivalent strain range from first to stabilized cycle is presented in Fig. 7.The equivalent strain range is calculated for all accelerated models and compared with results calculated considering the combined model (see Tab. 2).Results obtained with the combined model show a well know "s-shaped" curve which corresponds to the adopted nonlinear isotropic model.Quite similar "s-shaped" curves are also obtained with the accelerated models with a moderately increased speed of stabilization (b a =10b÷30b), while this shape is lost for higher values of b a (b a =100b÷300b).Similar behaviour is observed also for the bound case b a =0.01b, for which, as already shown, the model reaches stabilization in few cycles.In Fig. 7 also the case b a =421b is reported, representing the upper bound.In fact, for higher values of b a numerical problems do not permit the convergence to be reached.Due to the high speed of stabilization, the model wants to reach stabilization almost immediately somehow neglecting the nonlinear kinematic model that controls the monotonic hardening.
Fig. 8 shows the correlation between speed of stabilization and number of cycles to stabilization.As can be observed, plotting the data on a log-log scale gives almost a linear relation between b and N stab .Accelerated models with b a =0.01b and b a =421b obviously deviate from this trend.8) which correlates b, N and ∆ε pl , which is computed in the first cycles.All the three components of plastic strain range (∆ε pl,a =0.0554% axial direction, ∆ε pl,θ =0.0734% hoop direction and ∆ε pl,r =0.1289% radial direction) calculated by FE after the 5 th cycle are considered and therefore three curves are presented.As can be noticed results obtained with all accelerated models except the two bound cases (b a =0.01b and b a =421b) show a good agreement with the adopted relation.As pointed out in Section 2, very often in literature it is suggested to perform the simulation only for few cycles.As shown in Fig. 7, this approach in this case would introduce a significant error.For instance, the value of the equivalent strain range obtained after 10 cycles is about 15% lower than the value obtained in stabilized condition.It follows that, entering in the Manson-Coffin curve [16], the number of cycles to failure would be overestimated.Also the usual approach of adopting the kinematic model in the stabilized condition could rise significant errors.Fig. 10 shows the equivalent strain range for all the adopted models.The case of the Prager and nonlinear kinematic (initial and stabilized) models also reported.The Prager model is often proposed in literature because parameter C lin can be estimated by simply using monotonic uniaxial test.In the initial and stabilized nonlinear kinematic models the yield stress and the Young's modulus are estimated at the beginning of loading and at stabilization, respectively.Therefore, model parameters E and σ 0 are replaced by E s and σ 0* estimated from stabilized cycles, while kinematic variables (C i , γi) remain unaffected.All these approaches are based on a kinematic model and therefore capture only the monotonic hardening behaviour.Fig. 10.Comparison between all adopted cyclic plasticity models in terms of equivalent strain range.
Tab. 2 summarizes the results shown in Fig. 10.As can be noticed, neglecting the bound cases (b a =0.01b and b a =421b), Δε eq does not show a dependence with the speed of stabilization and its value substantially coincides with that obtained with the combined model.A relative error Δe=(Δε eq,a -Δε eq,c )/Δε eq,c is calculated where Δε eq,a and Δε eq,c are the equivalent strain ranges for the considered accelerated and combined model, respectively.In the case of b a =10b÷300b the error remains always in the range of -0.48÷0.5%.Result obtained with the nonlinear kinematic model that considers initial condition of a material (Non.Kin.initial) and neglects stabilized condition gives a high relative error (-15.26%) with respect to all adopted models.A quite similar error (-15.3%) is calculated with the accelerated model when b a =0.01b, since the model stabilized after 8 cycles.Due to the strong simplification, Prager and nonlinear kinematic-stabilized (Non.Kin.stabilized) models provide, instead, a smaller but still significant relative error (-8.6% and -10.71%).This result is in good agreement with the conclusions in [6], where it was observed that the direct use of the stabilized model could lead to heavy mistakes.Based on the obtained results, it is possible to conclude that neglecting the initial as well as stabilized condition of material could give misleading results that can heavily affect durability assessment, in particular over-estimating the life of the component.
As the computation time is almost proportional to the number of cycles to stabilization, it is possible to conclude that the accelerated model permits a strong reduction of the computational effort, keeping the same accuracy as the combined model.

Design rule for practical use in FEM modelling
As Eq. ( 8) shows, for decreasing values of plastic strain range ∆ε pl , the number of cycles to reach stabilization N stab increases proportionally.This case obviously falls in the "unfeasible" purple area depicted in Fig. 1.It is of practical interest to assess a design rule to deal with these cases.For this purpose, the thermal flux proposed in [14] was only 30% increased to reach a maximum temperature close to 250 °C for which material properties are available.
Firstly, the numerical simulation with the combined model was performed up to the complete stabilization (1033 cycles) in order to have the reference case, see Tab. 3. The ratio between the equivalent plastic strain range ∆ε eq,pl and the equivalent strain range ∆ε eq was 0.57, while in the previous case this ratio was 0.64 (as expected, in the latter case the component is less plastically strained).The three components (axial, hoop and radial) of the plastic strain range computed after 5 cycles were used to evaluate the corresponding curves, according to Eq. ( 8), see Fig. 11.Entering in the design diagram with the number of cycles to stabilization considered feasible to perform the numerical simulation, the speed of stabilization to be used in the accelerated model could be obtained.For example, assuming N stab =40 gives b a =200 (the curve relative to the axial strain component was used), see Fig 11 .Running the FE analysis with the accelerated model with b a =200 and performing a simulation of 40 cycles, a value of equivalent strain range ∆ε eq =0.3144% was obtained.With respect to the combined model, which requires to simulate 1033 cycles, in this case an error of only 2.27% was obtained with a huge reduction of the computational effort.

Conclusions
The choice of the adopted cyclic plasticity model in numerical simulations is an important design aspect, particularly when dealing with components under cyclic thermo-mechanical loading.In principle, the material model should follow quite accurately the cyclic evolution of the material behaviour and the numerical simulation has to be performed up to stabilization.Quite often this goal is unfeasible and it is thus necessary to refer to simplified models or to acceleration techniques.
In this work, several models (combined, accelerated, Prager, nonlinear kinematic initial and stabilized models) have been considered and compared in terms of the equivalent strain range.Based on the obtained results, it is possible to conclude that the use of too simplified models (Prager, nonlinear kinematic initial and stabilized) neglecting initial or stabilized conditions may be inadequate, as they could lead to overestimate the fatigue life of the component.On the other hand, accelerated models give results that are always close to the fully-stabilized combined model, assumed as reference.
A design rule is finally proposed.According to this approach, particularly useful when the plastic strain range induced by the thermo-mechanical loads is relatively small, few cycles of simulation have to be firstly performed.According to the obtained value of plastic strain range and to the computational time, a curve can be calculated, permitting a suitable value of increased speed of stabilization to be introduced in the model to perform the final FEM analysis.

Fig. 1 .
Fig. 1.Possible strategy to reduce the computational time.

Fig. 1
Fig. 1 presents schematically the desirable diagram required to perform a correct design of the component.In principle the basic idea of the accelerated techniques

Fig. 4 .
Fig. 4. Combined model: maximum stress a) and strain range b) versus number of cycles.

Fig. 7 .
Fig. 7. Equivalent strain range versus number of cycles to stabilization.

Fig. 8 .
Fig. 8. Correlation between the speed of stabilization and number of cycles to stabilization.

Fig. 9
Fig. 9 presents a "design diagram" calculated with Eq. (8) which correlates b, N and ∆ε pl , which is computed in the first cycles.All the three components of plastic strain range (∆ε pl,a =0.0554% axial direction, ∆ε pl,θ =0.0734% hoop direction and ∆ε pl,r =0.1289% radial direction) calculated by FE after the 5 th cycle are considered and therefore three curves are presented.

Table 2 .
Number of cycles to stabilization and ∆ε eq for point A (T=300 o C).

Table 3 .
Number of cycles to stabilization and ∆ε eq for point A (T=250 o C).