Strojniški vestnik - Journal of Mechanical Engineering 57(2011)6, 485-494 D01:10.5545/sv-jme.2010.212 Paper received: 07.10.2010 Paper accepted: 17.02.2011 Numerical Methods for TMF Cycle Modeling Henrik Zaletelj* - Gorazd Fajdiga - Marko Nagode University of Ljubljana, Faculty of Mechanical Engineering, Slovenia Experimental tests for the endurance evaluation of the machine parts that are exposed to thermo-mechanical fatigue (TMF) require advanced and expensive testing machines. Numerical methods for the determination of stress-strain material behavior have become very frequent and known due to lower costs. There are several different approaches for the determination of stress-strain behavior. In the article three different numerical methods and their results are presented. The numerical results for different load conditions are compared with the experimental results and the accuracy of the methods can be compared. The Chaboche, Skelton and Prandtl operator approaches are presented, presuming a stabilized elastoplastic response and not including creep. The properties of the model, their weaknesses and possible improvements are also studied in the paper. ©2011 Journal of Mechanical Engineering. All rights reserved. Keywords: cyclic loading, elastoplasticity, kinematic hardening, stress-strain trajectory, thermo-mechanical fatigue 0 INTRODUCTION Many components such as internal combustion engines, turbines, nuclear reactors, etc., are subjected to thermo-mechanical fatigue (TMF) [1] and [2]. Fatigue life depends primarily on loads, material, geometry and environmental effects. Its evolution is generally based on tests of three forms [3] and [4]: • isothermal strain-controlled low cycle fatigue (LCF) tests, • TMF tests on specimens and components, and • thermal shock tests. In view of their relative simplicity [5], LCF tests are often favored. The data from LCF tests conducted on servo-controlled uniaxial testing machines have been collected and tabulated for many years [3]. The key idea is, therefore, to predict fatigue life by avoiding expensive TMF and thermal shock tests. The present work is concerned with the cyclic stress-strain response for variable temperatures. Creep and transient effects, such as cyclic hardening and cyclic softening are not considered in this paper. Thus, hysteresis loops are supposed to be stabilized and constitutive equations for elastoplasticity are applicable for the stress-strain behavior modeling. To analyze the response of different calculating models, the material 9Cr2Mo alloy is used where the material parameters are presented in [5]. The paper presents the predicted TMF cycles based on the three different approaches. The Skelton, Chaboche and spring-slider model with the temperature dependant Prandtl densities [7] are verified with the experimental results. The strain and temperature are controlled. The paper is structured as follows. After the explanation of the material and TMF tests, the constitutive equations of the Chaboche model and the definition of the parameters are presented. The following section introduces the spring-slider model with the Prandtl densities. Then, the review of different model results and verification with several TMF tests follow. Finally, the final section lists the conclusions where the characteristics of the Prandtl and Chaboche models are introduced. 1 MATERIAL AND TYPES OF TMF TESTS 1.1 Material The material used to analyze the efficiency of the predicted stress-strain curves is advanced ferritic-martensitic steel EM 12 (9% Cr, 2% Mo). It is used as material in conventional thermal plants operating at temperatures up to 600 °C [5]. The material belongs to the cyclic softening class of alloys. Further details can be found in [5]. 1.2 TMF Cycles The cycle time for each loop was five minutes, i.e. at the strain rate of 4 x 10-5 s-1 for 0.6% strain range. Fig. 1 presents TMF cycles, where it is convenient to plot strain on the vertical axis and temperature on the horizontal axis. The strain range is 0.6% and temperature varies between 270 and 570 °C. Using the paths in Fig. 1, nine TMF types of cycles were performed to compare the stress-strain hysteresis curves of numerical models. Path PXRM represents a 45° kite cycle formed in the anticlockwise direction, while the corresponding clockwise cycle is taken in the PMRX order. A similar scheme applies in the 135° kite cycle PKRZ. The parallelogram-shaped PXRZ and PKRM mark the 45° zero strain and the 135° zero strain, respectively. Finally, complex cycle (dashed line) was considered, which shows a cycle of industrial gas turbine blades. A detailed explanation and the reasons for choosing these types of TMF loops are given in [5]. 0.3- 0.2- 0.1 - 0.0- ti- to -0.1 - -0.2- -0.3- -0 4- F M z H / / / ' \ \ \ \ / ' \ \ \ \ P \ K \ \ \ R \ \ \\ \ \ X / \ J \\ \ 7 V\ \ / v\ \/> v\ ^ / G V K E 250 300 350 400 450 7TC] 500 550 600 Fig. 1. TMF cycles For the observed (experimental) values of stress-strain curves, the Ramberg-Osgood parameters were used. The parameters of cyclic hardening coefficient A and cyclic hardening exponent ft are examined for isothermal loops for the temperature span of 270 to 570 °C [5]. Experimental curves were drawn for Eq. (1) where the corresponding total strain range, Ae, at any temperature is given by: while Act stands for the stress range. 2 THE CHABOCHE MODEL Most metals approach a cyclically stable state after a certain number of cycles. Cyclically stable or half-life material properties are usually used in fatigue analysis [7]. To describe the stable state of material, the Chaboche model of kinematic hardening can be used [8] and [9]. To the Chaboche model considerable attention has been paid due to its capacity of modeling a wide range of inelastic material behavior such as cycling hardening/softening, the Bauschinger effect, stress relaxation and creep for a range of materials [10]. Cyclic stress-strain curves are modeled with the uniaxial form of the Chaboche model. Stabilized cyclic curves are defined with kinematic hardening, which corresponds to the movement of the loading surface, where ct is stress at each moment and k is the yield stress. The hardening variable x (back-stress) indicates the present position of the loading surface. Back stress x also indicates the directional dependent effects, such as the Bauschinger effect. The criterion and the equation of flow and hardening can be expressed in the form [ll]: f(P, X, k) = \a- /j - k = 0 . (2) The evolution equation of the back-stress for non-linear kinematic hardening used in the Chaboche material model was originally introduced by Armstrong and Frederick [12]: Xt = Cfip -Y,X,P- (3) The Eq. presents the nonlinear kinematic hardening where p is the accumulated plastic strain, i = l, 2, x = Xi+X2. C and yt are temperature dependant material parameters. If yt = 0, Eq. (3) presents the model of the linear kinematic hardening (Prager's kinematic hardening law) [ll]. The essence of the model is in the velocity of plastic deformation sp and in the velocity of accumulated plastic deformation p : P = \Sp (4) The model in Eq. (3) describes kinematic hardening within one load cycle as well as kinematic hardening of the accumulated plastic deformation developed over more cycles to the saturation condition. The integration of Eq. (3) gives: (5) — Xi = -7 ( - exP {-YiSp)), / i and for i = 1, 2, 2 C X = X~I1 - exP(-Yisp)) andX = X1+X2. (6) i=i Yi The first part of kinematic hardening X1 describes the transition area of inelastic deformation, while the second part x2 describes the behavior at higher inelastic deformations after x1 reaches saturation value C1/y1. In tensioncompression, and more generally, in proportional loading, the evolution equation of hardening can be integrated analytically to give [11]: Ci = v — + Yi Xo-v — Iexp (( (ep -£po )), (7) Yî where v = ±1 according to the direction of the flow, ep0 and Xo denote the initial values, for example at the beginning of each plastic flow. It is not necessary to update variables ep0 and x0 from the previous flow. At each moment the stress is given by: o = X + vk. (8) 2.1 Parameter Estimation To model the cyclic curves, kinematic hardening variable of transition area X1 is used. With regard to the strain magnitude, the saturated value of x1 is not reached, so parameters C1 and y1 have to be defined. The presented model, Eqs. (7) and (8), contains four material parameters. These are Young's modulus E, kinematic hardening parameters C1 and y1 and the initial size of yield surface k. Young's modulus is presented in [5] as a function of temperature: E = a - bT, (9) where a = 2.08 x 105 MPa and b = 97.5 for T in °C. Parameters Q and y1 can be estimated from the tension part of the cyclically stable rising hysteresis branch. They can also be estimated on the first cyclically stable stress-strain curve. The latter option is used in this paper. The estimation is made at a low level of plastic strain, where the transition kinematic hardening X\ is more obvious. As a > x it follows: £p = a - -R - k Z (10) where a presents total stress, X kinematic hardening, R isotropic hardening, Z and n viscous parameters and k yield surface. Viscosity stress is defined as: av = a - x - R - k, (11) As the stabilized cycle is observed where R = 0 by considering the constant k parameter and viscosity stress av = 0, the derivate of Eq. (11) is presented as: du ds„ dX ds„ (12) Considering Eq. (3), the logarithm of Eq. (12) is expressed as: ln ds = ln(C1) -Yisp (13) p J Parameters C1 and yx are determined from the linear regression. Parameter -yj presents the line slope and parameter ln(C1) the intersection with the ordinate axis. a and sp are evaluated with the Ramberg-Osgood equation with the parameters that are presented in the article [5]. The parameters are evaluated for the temperatures ranging from 270 to 570 °C with the interspaces of 30 °C. To obtain the material parameters that were not measured, the linear parameter interpolation is used. The estimated values are used as initial values in the optimization process, where the parameters were finally fitted on the stabilized cycle loop. The aim of the optimization is to find the minimum difference between the back-stress values of the first cycle and the back-stress of the stabilized cycle. The parameters are presented in Table 1. Table 1. parameters Chaboche and Ramberg-Osgood Temperature C1 Y1 E A [°C] [MPa] [MPa] [MPa] ß 270 104541 181675 1611 0.135 300 101844 178750 1554 0.127 330 101344 175825 1534 0.123 360 113963 172900 1763 0.147 390 92139 169975 1423 0.115 420 97456 192 167050 1505 0.126 450 89232 164125 1395 0.121 480 71252 161200 1133 0.095 510 74475 158275 1165 0.108 540 59331 155350 964 0.091 570 63820 152425 990 0.108 2.2 Temperature Dependence in the Back-Stress Evolution Equation The response and evolution of kinematic hardening is dependent on temperature. Temperature influence is considered in parameter dependence; besides, the changing of temperature vs. time is taken into account. If the temperature changing vs. time is fast, the influence is high. The influence is negligible for long periods. Fig. 2 presents the graph where the stress vs. strain is dependent on temperature and time. The temperature at s = 0 was 270 °C and at maximum load s = 0.006 it was 570 °C. Eq. (14) represents the back-stress dependence of temperature (only one back-stress is considered here). As compared to Eq. (3), temperature influence is introduced directly by the variation of parameter C. 1 dC X = C(T)sp-yxP + ——-L%T. (14) Ci(T) dT The main advantage of the model is that the equation does not contain new material parameters. It uses the parameters in dependence on the temperature where the partial derivatives upon temperature can be estimated from. 2.3 Relaxation of Mean Stress If the load is not purely alternating, additional effects can occur. In a strain-controlled test, when the mean strain is not zero, the phenomena of the relaxation of the mean stress appear. The initial asymmetry of the stress Fig. 2. Stress-strain dependence on temperature and time 0.000 0.002 0.004 0.006 [/] Fig. 3. Relaxation of mean stress disappears progressively in the first few cycles, Fig. 3. The temperature changing inside the load cycle defines the level of mean stress. Temperature load cycle defines the shape of the hysteresis curve, maximum and minimum stress. In Fig. 4, the position of a stabilize curve is shown for the changing the temperature between 270 and 570 °C. The mean stress of the stabilized curve where the temperature is constant equals zero, Fig. 4. Fig. 4. Stabilized curves and temperature variation 3 PRANDTL OPERATOR APPROACH It is standard practice in the isothermal strain-life approach to use cyclic material properties together with the Masing and Memory models to define the cyclic uniaxial response of the material, which determines the stress and strain range of the closed hysteresis loop and the mean stress associated with each loop. It has been shown in [13] that the Masing and Memory models are not valid if the temperature varies during the cycle. Masing based his finding on the rheological spring-slider model and assumed that the model parameters are time independent. As the spring-slider model supports the modeling of elastoplastic hardening, it has been adapted for variable temperatures [7]. The model is capable of modeling elastoplastic hardening solids and nonlinear kinematic hardening under strain control (see Fig. 5). The stress controlled model is given in [14]. From the equilibrium in a single springslider segment, total strain s is obtained s = sqa- + saj where slider strain | sqj- | can never exceed fictive yield strain q^ that is also known as the half-width of the play operator. Spring strain saj can now be expressed as the play operator with general initial values [7]: \e(ti )- qj, saj (t ) = max j . f. (15) j Imm {) + qj ^ ^JM Fig. 5. Rheological spring-slider model [7] for 0 < t1 < t2 < ... < ti<... . Thus, current strain state ea](ti) depends on the previous saj(ti_l) called the memory point. Presumably, there is no residual strain initially, so ea;(0) = 0 and cj(0) = 0. Determination of the parameter nq is thoroughly described in [26]. The stress in the spring-slider segment is then: a, (ti) = Ej (Ti )eaJ (ti) = a, (Ti )eaJ (ti), (16) where Tt = T(ti) and aj(Tl) is the Prandtl density. Adding the spring-slider stresses results in total stress in the form known as the operator of the Prandtl type [7]: v(tt) = (Ti )eaj (ti), (17) j=i with temperature-dependant Prandtl densities. The play operator given in Eq. (16) is independent of time and temperature. Therefore, it is modified [7] to assure equilibrium in the spring-slider: eai (ti ) = maX s(ti) - qj, 's(ti) + qj, a. (TJ mm ^ !-1/ eaJ(ti-!) a j (T ) (18) When Eq. (18) is inserted into both Eqs. (17) and (19), temperature-dependant stress-strain behavior can be modeled because the temperature-modified play operator guarantees equilibrium in the spring-slider segments at any time and temperature: ay (ti_x) = a. (Ti _l)s(ti_l). (19) The Prandtl densities can be calculated from the temperature-dependent Ramberg-Osgood curves. For a further explanation, interested readers should see papers [13] to [21]. In newer publications index j in Eq. 17 starts with 1 instead of 0 to simplify the notation. The Skelton approach is thoroughly discussed in [3], [5] and [23]. 4 VERIFICATION OF MODELS The Chaboche, Skelton and Prandtl operator approaches have been compared to several TMF tests conducted by Skelton [5] at the total strain range of 0.6 % on the 9Cr2Mo alloy. The paper is concerned with the paths given in Fig. 1. The observed hysteresis loops from the tests are plotted as crosses in Figs. 6 to 15. The circles and the dot line denote the stress-strain trajectories modeled by the Skelton and the Prandtl operator approach, respectively. The thin solid line denote the stabilized cycle of the Chaboche non-linear kinematic hardening model. The thick solid line presents the shifted Chaboche hysteresis curve 400 Fig. 6. Shift of the Chaboche hysteresis curve (load case PZRKP) 300- 200- . 100- -100- -200- -300 -400 .o'-j,—y yo-y^ /Á ¿O /.' p / fft /7+ t/+ y? p / ^^ + observed ......Prandtl Chaboche o Skelton -0.3 -0.2 -0.1 0.0 0.1 0.2 £ [%] Fig. 8. TMF cycle PMRXP 0.3 400 300 200- 100- 0- b-100- -200- -300- -400 o ,o--o«Z jji-^ri + Q^v + S^ '' + A + A P ''' + / '' + P '' .'/ + / '' /P + / P J + / ty + P x>y + / ^r^ observed - Prandtl - Chaboche K Skelton -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 £ [%] Fig. 7. TMF cycle PZRKP 400 300 200- ' 100 -100- -200 -300 -400 o ¿9 M o* sSy ' + A R /'+ A p 4- / f .6/ + observed Prandtl O Chaboche Skelton -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 £ [%] Fig. 9. TMF cycle PXRMP (thin line) that is fitted to the observed results, Fig. 6. The shift is required due to unwanted ratcheting effect built into the Chaboche model. The ratcheting effect is noticed at different load conditions. It happens due to a small amount of plastic strain in each cycle, which leads to unacceptable accumulated strain. This is true even for the material that does not intrinsically present a risk of ratcheting [24]. The ratcheting effect is dependent on the TMF load as well as the changing of the temperature vs. time. The presented Chaboche model does not consider the elimination of the ratcheting effect, so only the shape of the curve can be observed. To take into consideration the effect of ratcheting, several kinematic hardening parameters have to be defined requesting additional work and calculations. Figs. 7 to 15 present the result agreement of different numerical models to the observed values. It can be seen from the figures that the results of the Prandtl operator approach fit the Chaboche model well. A better agreement of the Prandtl operator results is noticed with the Skelton model. The stress-strain behaviors of all models 400 300 200 ■ 100- -100 -200- -300 -400 - - ° $$ /: ' * ^ ■7 R.' a i -f/ o'Y o./ 7 P ' y o'y °7 7 9'/ 7 ■7 rr y ' „ + + observed ......Prandtl Chaboche o Skelton —1—i—1—i—1—i—1—i—1—r -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 £ [%] Fig. 10. TMF cycle PKRZP 400 ~i—1—i—1—i—1—i—1—i—1—r -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 £ [%] Fig. 12. TMF cycle PKRMP 400 300- 200 100 - & 0- -100 -200- -300 .-oj R ■ jP^ 97 7 ? 7 7 • 7/ 7 ■ $/ 7 p-cy / ' + observed / + .....Prandtl rf Chaboche o Skelton ~i—1—i—1—i—1—i—1—i—1—r -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 £ [%] Fig. 11. TMF cycle PXRZP 400 300- 200- 100 0- -100- -200 -300- -400 ^77 P /' 7'' 7 P J y J-S / ^^ R +7 /7 ■ + observed ......Prandtl Chaboche xJ^9 o Skelton H—1—i—1—i—1—r -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 £ [%] Fig. 13. TMF cycle PZRXP are in good agreement with the observed values. The shapes of the stress-strain curves for the TMF loops are predicted well, as well as the minimum and maximum value of stress. The advantage of the Prandtl operator approach is the fact that it does not take into account the effect of ratcheting, which is a benefit as compared to the Chaboche model. For more precise results of the Chaboche model, the elimination of ratcheting should be included. A comparison of the results is also applied for the complex cycle, Fig. 15. The shapes of the curves deviate from the observed values. It should be noted that experimental testing cannot perform the changing of strain vs. temperature as linearly as it can in numerical calculation. 5 CONCLUSIONS The results of the three models are introduced and compared to the experimental TMF cycles. The Prandtl operator approach is compared to the Skelton and Chaboche models. The non-linear kinematic hardening model is used in the framework of time-independent plasticity to model the stabilized curves. The influence of temperature is taken into account in all models. Temperature influence is introduced as parameter dependence as well as changing temperature vs. time. The results of the Chaboche model indicate the ratcheting effect, which is difficult to eliminate. Its effect depends on the TMF cycle. The classical Chaboche constitutive equations do not describe the ratcheting effect correctly, especially the ones observed when the mean-stress is significantly lower than the stress amplitude [24]. The non-linear kinematic model greatly over predicts ratcheting when its identification is performed for normal monotonic and reversed cyclic conditions. In [24] a set of modified kinematic rules for the elimination of ratcheting is introduced. The non-linear kinematic model with a threshold presents the best choice to describe both the normal cyclic condition and the ratcheting condition [24]. The comparison of the results of different numerical models shows good agreement with the observed values. The ratcheting effect can cause higher deviation from the predicted results and for this reason the Prandtl operator approach is preferred to the Chaboche model. For precise results of the Chaboche model the correction of the ratcheting effect should be considered. 6 REFERENCES [1] Constantinescu, A., Charkaluk, E., Lederer, G., Verger, L. (2004). A computational approach to thermomechanical fatigue. International Journal of Fatigue, vol. 26, p. 805-818. Fig. 14. TMF cycle PMRKP Fig. 15. TMF cycle — dashed line [2] Muhič, M., Tušek, J., Kosel, F., Klobčar, D. (2010). Analysis of Die Casting Tool Material. Strojniški vestnik - Journal of Mechanical Engineering, vol. 56, no. 6, p. 351-356. [3] Skelton, R.P., Webster, G.A. (1996). History effects in the cyclic stress-strain response of a polycrystalline and single crystal nickel-base superalloy. Mater. Sci. Engen., vol. A 216, p. 139-154. [4] Charkaluk, E., Bignonnet, A., Constantinescu, A., Dang, V.K. (2002). Fatigue design of structures under thermomechanical loadings. Fatigue Fract. Engng. Mater. Struct., vol. 25, p. 1199-1206. [5] Skelton, R.P. (2004). Hysteresis, yield, and energy dissipation during thermo-mechanical fatigue of ferritic steel. International Journal of Fatigue, vol. 26, p. 253-264. [6] Urevc, J., Koc, P., Štok, B. (2009). Numerical simulation of stress relieving of an austenite stainless steel. Strojniški vestnik - Journal of Mechanical Engineering, vol. 55, no. 10, p. 590-598. [7] Nagode, M., Fajdiga, M. (2005). Temperature-stress-strain trajectory modeling during thermo-mechanical fatigue. Fatigue Fract. Engng Mater. Struct., vol. 29, p. 175182. [8] Chaboche, J.L. (1989). Constitutive equations for cyclic plasticity and cyclic viscoplasticity. International Journal of Plasticity, vol. 5, p. 247-302. [9] Chaboche, J.L. (2008). A review of some plasticity and viscoplasticity constitutive theories. International Journal of Plasticity, vol. 24, p. 1642-1693. [10] Tong, J., Vermeulen, B. (2003). The description of cyclic plasticity and viscoplasticity of waspaloy using unified constitutive equations. International Journal of Fatigue, vol. 25, p. 413-420. [11] Lemaitre, J., Chaboche, J.L. (1990). Mechanics of solid materials. Cambridge University Press, Cambridge. [12] Armstrong, P.J., Frederick, C.O. (1966). A mathematical representation of the multiaxial Bauschinger effect. CEGB Report RD/B/ N731. [13] Nagode, M., Zingsheim, F. (2004). An online algorithm for temperature influenced fatigue- life estimation. Strain-life approach. Int. J. Fatigue, vol. 26, p. 155-161. [14] Nagode, M., Hack, M. (2004). An online algorithm for temperature influenced fatigue-life estimation. Stress-life approach. Int. J. Fatigue, vol. 26, p. 163-171. [15] Conle, A., Oxland, T.R., Topper, T.H. (1988). Computer-based prediction of cyclic deformation and fatigue behaviour. Low Cycle Fatigue, Solomon, H.D., Halford, G.R., Kaisand, L.R., Leis, B.N. (Eds.). ASTM STP 942, p. 1218-1236. [16] Brokate, M., Sprekels, J. (1996). Hysteresis and Phase Transition. Applied Material Science 121. Springer Verlag, New York. [17] Nagode, M., Fajdiga, M. (2006). Thermo-Mechanical Modeling of Stochastic StressStrain States. Strojniški vestnik - Journal of Mechanical Engineering, vol. 52, no. 2, p. 74-84. [18] Brokate, M., Rachinskii, D. (2005). On global stability of the scalar Chaboche models. Nonlinear Anal.: Real World Appl., vol. 6, p. 67-82. [19] Hack, M. (1998). Schadigungsbasierte Hysteresefilter. Shaker Verlag, Aachen. [20] Nagode, M., Hack, M., Fajdiga, M. (2009). Low cycle thermo-mechanical fatigue. damage operator approach. Fatigue Fract Engng. Mater. Struct., vol. 33, p. 149-160. [21] Nagode, M., Langler, F., Hack, M. (2011) A time-dependant damage operator approach to thermo-mechanical fatigue of Ni-resist D-5S. International Journal of Fatigue, vol. 33, p. 692-699. [22] Jovičic, G., Živkovic, M., Jovičic, N. (2009). Numerical simulation of crack modeling using extended finite element method. Strojniški vestnik - Journal of mechanical engineering, vol. 55, no. 9, p. 549-554. [23] Skelton, R.P., Webster, G.A., de Mestral, B., Wang, C.Y. (2000). Modeling thermo-mechanical fatigue hysteresis loops from isothermal cyclic data. Thermo-Mechanical Fatigue Behaviour of Materials. Sehitoglu, H. Maier, H.J. (Eds.). ASTM STP 1371, p. 69-84. [24] Chaboche, J.L. (1991). On some modifications of kinematic hardening to improve the description of ratcheting effects. International Journal of Plasticity, vol. 7, p. 661-678. [25] Stamenkovic, D., Maksimovic, K., Nikolic-Stanojevic, V., Maksimovic, S., Stupar, S., Vasovic, I. (2010). Fatigue life estimation of notched structural components. Strojniški vestnik - Journal of mechanical engineering, vol. 56, no. 12, p. 846-852. [26] Laug, P., Borouchaki, H. (2004). Curve linearization and discretization for meshing composite parametric surfaces. Commun. Numer. Meth. Engng., vol. 20, p. 869-876.