Quadratic Midpoint Integration Method for J2 Metal Plasticity

Abstract: The quadratic variants of the generalized midpoint rule and return map algorithm for the J2 von Mises metal plasticity model are examined for the accuracy of deviatoric stress integration of the constitutive equation. The accuracy of stress integration using a strain rate vector for arbitrary direction is presented in terms of an iso-error map for comparison with the exact solution. Accuracy and stability issues of the quadratic integration method are discussed using a two-dimensional metal panel problem with a single slit-like defect in the center. The scale factor and shape factor were introduced to a quadratic integration rule for assuming a returning directional tensor from a trial stress onto the final yield surface. Luckily enough, the perfectly plastic model is the only case where the analytical solution is possible. Thus, solution accuracies were compared with those of the exact solutions. Since the standard scale factor ranges from 0 to 1, which is similar to the linear α-method, the penalty scale factors that are greater than 1 were mainly explored to examine the solution accuracies and computational efficiency. A higher value of scale factor above five shows a better computational efficiency but a decreased solution accuracy, especially in the higher plastification zone. A well-balanced scale factor for both computational efficiency and solution accuracy was found to be between one and five. The trade-off scale factor was proposed to be five. The proper shape factor was also proposed to be {1,1,4}/6 among the different combinations of weight distribution over a time interval. This proposed scale factor and shape factor is also valid for relatively long time periods.


Introduction
Calculation of crack extension and the behavior of the fracture process zone ahead of the current real void (or crack) in ductile materials (such as metals) strongly depends on the computational method used for stress analysis. Rice and Tracey [1] introduced an iterative procedure to integrate deviatoric stress under loading conditions, especially in J 2 metal plasticity. In addition, Krieg et al. [2] explored the accuracies of several solution methods, including the tangent stiffness, secant stiffness, and radial return method compared to the exact solutions obtained from the elastic perfectly plastic model. This perfectly plastic model is the only case where an analytical solution is possible. The exact solution was proposed by Krieg et al. [2] to the perfect plasticity equations for the general deviatoric stress increment. There were no approximations involved in deriving the starting differential equation of incremental plasticity. Based on of these pioneering papers, Ortiz and Popov [3] and Simo et al. [4,5] elaborated the return mapping algorithms for the plane stress problem. The linear α method was used to integrate the new stress state over the time interval. By assuming the new stress state in terms of a linear combination between current known stress and stress to be determined, the plastic corrector term has changed accordingly. This means that the direction of mapping onto the final yield surface has changed with varying α value. Thus, this directional change onto the yield surface may subsequently affect to computation time to calculate a numerical subroutine. To elaborate this issue, Willam [6] and Shim [7] proposed the higher order α method by expanding to additional time integration points and weights, while the time interval remained as a single step. Artioli et al. [8] and Jahanshahi [9] introduced multi-step methods for von Mises materials, specifically single-step midpoint method (SMPT1 and SMPT2) and double-step midpoint method (DMPT1 and DMPT2). The first two schemes used a generalized midpoint integration rule with return mapping procedure enforced by yield conditions. The latter two schemes were two-stage algorithms developed by dividing each time step into two subintervals, in which equations were solved in turn. These approaches essentially rely on sub-steps in a time interval. They develop an extensive comparison based on pointwise mixed stress-strain loading histories, iso-error maps and initial boundary value problem. They investigated the existence of solution, accuracy, stability and the algorithm behavior for long time steps. Since Willam [6] and Shim [7] proposed the higher order α method to evaluate the new stress state, they did not implement this scheme in a finite element problem. Hence, this paper examines the accuracies of the quadratic α tmethod with the return mapping algorithm for solution of a uniaxial tension problem with a slit in the center. It is instructional to examine the solution using the perfectly plastic model, which only has the analytical solution available.

J 2 Metal Plasticity
Tresca [10] concluded that yielding in metals would occur when maximum shear stress reaches a yield value. This maximum shear stress is the difference between the major and minor principal stresses. Later, von Mises [11] proposed that the yielding of metals is governed by the second invariants of the deviatoric stress, denoted as J 2 . This von Mises yield criterion is known to be in better agreement with experiments for ductile metals such as copper, nickel, aluminum, and alloy steels [12]. This paper focuses on the analytical formulation and numerical implementation of the von Mises criterion, which we shall simply refer to as J 2 plasticity. In metals, irreversible deformation is largely due to crystal dislocation, generally evident as irreversible deformation behavior when an applied load is removed.
where p = tr(σ)/3 is the mean normal stress which can be expressed as hydrostatic stress, 1 is the second-rank identity tensor, s is the deviatoric stress tensor, ε v = tr(ε) is volumetric strain, and e is the deviatoric strain tensor. For isotropic linear elastic material, the elastic constitutive equations are where, K and G are the elastic bulk and shear modulus. This gives where, the symbol ⊗ denotes dyadic product for any second-order symmetric tensors, : denotes an inner product (double contraction), and I denotes the rank-four symmetric identity tensor. J 2 is the second invariants of the deviatoric stress tensor, s is defined as where, • is the norm of second-order tensor. σ 1 , σ 2 and σ 3 are the principal stresses. The yield function of von Mises metals is defined as This yield function can be rephrased in another form such as, f = √ s : s − σ y and f = s − σ y . Here, for simplicity, the inelastic behavior of von Mises materials with perfectly plastic yielding is discussed as follows: (6)- (10). The associated flow rule for plastic strain is illustrated in Figure 1. The flow rule defines symmetric tensor of plastic strain rate, with unit direction, n and norm . λ. One can rewrite the term, tr( . e p ) by . λtr(n)m which is the trace of deviatoric strain change. Since this term is zero, there is no plastic volume change and normally referred as isochoric deformation. Here, for simplicity, the inelastic behavior of von Mises materials with perfectly plastic yielding is discussed as follows: (6)- (10). The associated flow rule for plastic strain is illustrated in Figure 1. The flow rule defines symmetric tensor of plastic strain rate, with unit direction, n and norm λ  . One can rewrite the term, ( ) p tr e  by ( ) tr λ n  m which is the trace of deviatoric strain change. Since this term is zero, there is no plastic volume change and normally referred as isochoric deformation.

( )
: : This fourth-order tangent operator ep C represents the additive relations for both isotropic bulk Elastoplastic tangent modulus: .
This fourth-order tangent operator C ep represents the additive relations for both isotropic bulk modulus, K = K1 ⊗ 1 and plastic shear operator G ep = 2G(I − (1/3) · 1 ⊗ 1 − n ⊗ n). The plastic shear operator, G ep can be rank-one updated in the form of n ⊗ n when plastic flow evolves. The rank-one update means that one of eigenvalues in G ep will be updated according to the level of plastification. Equation (10) can be written as the incremental format as in Equation (11). s n+1 − s n = C ep : e n+1 − e n or s n+1 = s n + C ep : e n+1 − e n = s n + C ep : ∆e (11) where n is the current time step which has known stress state, n + 1 is the unknown new stress state which has to be evaluated according to deviatoric strain rate, ∆e. This new deviatoric stress, s n+1 must be satisfied with the yield function. Radial return algorithm allows one to integrate the rate constitutive equations for J 2 plasticity. The typical procedure of this algorithm is to write the stress tensor in the elastic predictor and plastic corrector format as in Equation (12).
where, the trial stress, s tr n+1 = s n + ∆s = s n + 2G∆e. Using the backward implicit scheme for integrating the incremental plastic flow, λn dt ≈ ∆λn n+1 , n n+1 = s n+1 / s n+1 (13) where, ∆λ is the discrete plastic multiplier. In the meantime, . f = 0 is the consistency condition for the rate equation, the discrete consistency condition is the following Equation (14). More detailed stress integration procedure using the radial return algorithm is discussed in the following section.

J 2 Metal Plasticity with Linear α Method
The time increment over the loading interval ∆t = t n+1 − t n is constant for each loading step for [t o , · · · t n , t n+1 , · · · t f ]. Intermediate time integration (midpoint time integration) can be bounded by t n ≤ t n+α ≤ t n+1 , where n is the current known state, n + 1 is the state to be determined and α is defined as (15) and Figure 2 shows the linearization of this time integration graphically. Equation (16) defines the midpoint deviatoric stress, s n+α , in terms of α, which can vary between 0 and 1: one update means that one of eigenvalues in ep G will be updated according to the level of plastification. Equation (10) can be written as the incremental format as in Equation (11).
( )  (11) where n is the current time step which has known stress state, n + 1 is the unknown new stress state which has to be evaluated according to deviatoric strain rate, Δe . This new deviatoric stress, 1 n+ s must be satisfied with the yield function. Radial return algorithm allows one to integrate the rate constitutive equations for J2 plasticity. The typical procedure of this algorithm is to write the stress tensor in the elastic predictor and plastic corrector format as in Equation (12 (12) where, the trial stress,  (13) where, λ Δ is the discrete plastic multiplier. In the meantime, is the consistency condition for the rate equation, the discrete consistency condition is the following Equation (14). More detailed stress integration procedure using the radial return algorithm is discussed in the following section.

J2 Metal Plasticity with Linear α Method
The time increment over the loading interval  , where n is the current known state, 1 n + is the state to be determined and α is defined as (15) and Figure 2 shows the linearization of this time integration graphically. Equation (16) defines the midpoint deviatoric stress, n α + s , in terms of α , which can vary between 0 and 1:   By applying this time linearization into the midpoint strain, ε n+α , s n+α can be determined. The elastic deviatoric stress increment, ∆s = 2G∆e e = 2G(∆e − ∆e p ) is expressed by the elastic deviatoric strain increment, where G is the shear modulus and ∆e and ∆e p are the total deviatoric strain increment and plastic deviatoric strain increment respectively over the time increment, ∆t. Let s n is the known stress state in the Figure 3a. The trial stress, s tr n+1 can be assumed by adding the elastic predictor as ∆s. It must be corrected using the plastic flow rule as defined in (17). can be assumed by adding the elastic predictor as Δs . It must be corrected using the plastic flow rule as defined in (17).  Figure 3a is the specific case which the current stress state resides on the yield surface. Figure 3b is the general case which the current stress state resides inside the yield surface. This means that one needs to calculate the contact stress, c s when Δs cross through the yield surface. The unknown stress determination using the radial return method in the three-dimensional deviatoric space is shown in Figure 4. As a result, ∆e = ∆e − ∆e p can be evaluated and must satisfy the consistency condition at the new stress state as depicted in Equation (18) at t n+1 . Finally, one can describe the new stress state as (18) where s tr n+1 = s n + 2G∆e, which is the trial stress (elastic predictor) at t n+1 . By plastically correcting the trial stress state iteratively, one can finally obtain s n+1 : Figure 3a is the specific case which the current stress state resides on the yield surface. Figure 3b is the general case which the current stress state resides inside the yield surface. This means that one needs to calculate the contact stress, s c when ∆s cross through the yield surface. The unknown stress determination using the radial return method in the three-dimensional deviatoric space is shown in Figure 4.
Krieg et al. [2] evaluated the exact stress state at t n+1 in the perfectly plastic model. In their pioneering work, they compared the three techniques to correct the trial stress: (a) the tangent stiffness-radial return method, (b) the secant stiffness method, and (c) the radial return method. The errors between the exact stress and the stress determined by the radial return method with different α values were plotted in the form of an iso-error map, as shown in Figure 5. The errors were plotted in terms of the angle, θ between s Krieg E and s n+1 as illustrated in Figure 3b. The error, θ is defined by Krieg et al. [2] as cos −1 s n+1 : s Kreig E /R 2 . R is the radius of deviatoric principal stress plane as shown in Figure 3. In this case, R is the yield strength, σ y or s . If s n+1 lies between s Krieg E and s c , then the error, θ becomes positive. s c is the contact stress when elastic prediction process over a time span such that the yield surface is reached. For the monotonic loading with an arbitrary deviatoric stress increment, ∆s between first quadrant of the radial and tangential coordinates in the Figure 3a, the strain rate, ∆e = e n+1 − e n was controlled up to five times R of ∆s for both the radial and tangential direction. The stress increment, ∆s for marching toward the trial stress state, s tr n+1 is simply scaled in terms of the reference value, radius R. This is for the purpose of evaluating the solution accuracies due to a different location of elastic predictor in the first quadrant of radial and tangential coordinates. Figure 5 shows various iso-error maps for the linear α method. The accuracy improves when α = 1, known as the radial return method. The material properties used were E = 210 GPa, ν = 0.3 and σ y = 900 MPa used. Figure 6 shows the two parallel solution algorithms for finite element analysis: the rate format and the incremental format. The two formulations are equivalent when the material is in the plastic state at t n and the strain increment, ∆e is parallel to the current known deviatoric stress state, s n . Of course, when the time increment is very small, these two requirements are satisfied.   s is simply scaled in terms of the reference value, radius R. This is for the purpose of evaluating the solution accuracies due to a different location of elastic predictor in the first quadrant of radial and tangential coordinates. Figure 5 shows various iso-error maps for the linear α method. The accuracy improves when 1 α = , known as the radial return method. The material properties used were E = 210 GPa, ν = 0.3 and y σ = 900 MPa used. Figure 6 shows the two parallel solution algorithms for finite element analysis: the rate format and the incremental format. The two formulations are equivalent when the material is in the plastic state at n t and the strain increment, Δe is parallel to the current known deviatoric stress state, n s . Of course, when the time increment

Accuracy of Quadratic Midpoint Integration Method
A simple expansion of the linear α -method as discussed in the previous section was proposed by Shim [7]. A midpoint time, n t α + , bounded by the starting time, n t , and the arrival time, k is a number of interpolation point. In our case, 2 k = .

Accuracy of Quadratic Midpoint Integration Method
A simple expansion of the linear α-method as discussed in the previous section was proposed by Shim [7]. A midpoint time, t n+α , bounded by the starting time, t n , and the arrival time, t n+1 were added. These three interpolation points allow manipulation of the quadratic midpoint time integration by incorporating the Lagrangian interpolation function, L(x) = k ∑ j=0 y i · l j (x) as depicted in Equations (19) and (20). This typical sum of time interval-deviatoric stress products consists of x o , x 1 , x 2 = t n , t n+α , t n+1 , y o , y 1 , y 2 = s n , s n+α , s n+1 . k is a number of interpolation point. In our case, k = 2.
Using this polynomial expansion, one can express the plastic strain increment, ∆e p = L(t), where s n+α = (1 − α)s n + αs n+1 . For instance, α = 1/2, β 1 (t) = 2t 2 − 3t + 1, β 2 (t) = −4t 2 + 4t, and β 3 (t) = 2t 2 − t, when t n , t n+α , t n+1 = 0, 1/2, 1 . The summation of β i is equal to unity. Thus, α is the scale factor and β is the shape factor as illustrated in Figure 7. Since the backward implicit scheme, α = 1 in the linear method shows the best performance in a solution accuracy and a computation time. In this quadratic α method, we tried to explore a penalty scale factor that is higher than 1. Although this approach can be seen in the non-logical sense at a glance, it is worthwhile to examine its performance for both solution accuracy and computation time. The accuracies of the quadratic integration methods are plotted in Figure 8. When α = 1 and β 1 , β 2 , β 3 = 0, 1, 0 , the method is identical to the backward implicit scheme in the linear α-method. In addition, when α = 1/2 and β 1 , β 2 , β 3 = 0, 1, 0 , the central difference method is yielded. Thus, as long as β = 0, 1, 0 , the quadratic integration method coincides with the linear α-method. Figure 8b-f shows the results with a range of α values and indicates a violation of the interval, 0 ≤ α ≤ 1. Interestingly, the solution errors compared to Krieg's exact solution decrease as α increases. Figure 8a shows the entire negative solution error in the radial-tangential stress increment domain. However, as α increases up to 3, the solution error becomes split into negative and positive components, and its deviation from the zero error iso-line is reduced, as shown in Figure 8b,c, compared to the backward implicit scheme in Figure 8a. In addition, the level of absolute error is reduced in the penalty α method. However, when α > 5, the negative error diminishes gradually and the positive error prevails, especially in the narrow range of tangential stress increments within ∆s/R < 2. This implies that if α is greater than 1, then the positive error within ∆s/R < 2 is increased significantly. This will lose solution accuracy significantly. On the other hand, the relatively large stress increment of ∆s/R > 2 gives a more accurate estimation compared with that obtained using the standard backward Euler method. This might be helpful for relatively long time steps. In Figure 9, the effects of a range of values of β on the solution error were investigated, with a fixed value of α = 5. Better solution accuracies were found in Figure 9c,e. The similarity between these two cases arises from the higher weight for s n+1 as in the backward Euler approach. Hence, the influence of the quadratic time integration method spreads the solution error into the positive and negative ranges, while the relative error deviation is similar from the perspective of solution accuracy. Nevertheless, for the case of a large stress increment, the solution error would decay much faster than in the linear α-method. This differs from the results reported by Shim [7]. An iterative solution algorithm for implementing this quadratic time integration the finite element programming is presented in Figure 10. By altering the trial stress in terms of n n+1 , s n+1 · n n+1 = s tr n+1 · n tr n+1 − 2G∆γn n+1 , one can find ∆γ under the assumption that n n+1 n tr n+1 so that s n+1 can be defined as s tr n+1 − 2G∆γ, and the constrained equation in the n + 1 step at f tr n+1 = 0 is finally determined by ∆γ as f tr n+1 /(2G + K). K is the bulk modulus.

Implementation for Quadratic α Method
The Matlab-based finite element program developed by Ballani [13] was used as a base program for implementing the quadratic α method in the stress integration module and associated dependencies. The single-step midpoint methods (SMPTs) applied a generalized midpoint rule for integration of the differential equations of the system. These stress integration schemes were performed once per time interval. The difference between SMPT1 and SMPT2 is the choice of consistency condition at time interval, t n+1 or t n+α . As we introduced the scale factor, α and the shape factor, β earlier, one can implement the quadratic SMPT1 and SMPT2 as shown in Figure 10. A Newton-Rapson iterative solution technique was used to compute the plastic consistency parameter, λ. To calculate the consistent tangent moduli, C n+1 in both the quadratic SMPT1 and SMPT2, the formulations of Equations (21)-(26) were used.
For the quadratic SMPT 1: where, 1 ⊗ 1 is a rank-four tensor, I 4 is a fourth-order identity tensor and the values of χ 1 , χ 2 , χ 3 are defined as follows.
Finally, we get For the quadratic SMPT 2: Two scale factors were applied to examine the speed of convergence of the calculation of the plastic consistency parameter: the standard α method and the so-called penalty α method, as shown in Figure 10. The penalty means that if α is greater than 1, then the interpolation function layout will be skewed from standard quadratic method as much as the α value. Different shape factors were introduced as a weighting method, but the sum of the three components does not necessarily unite as shown in Figure 11.

2D Metal Panel Problem with Slit-Like Defect
A 1 × 1 mm steel plate with a 0.25 mm slit-like defect in the center was discretized using the Delaunay triangulation scheme with a two-dimensional mesh generator (Persson and Strang, [14]), as shown in Figure 12a,b. A total of 1423 nodes and 2606 CST triangular elements were used. A Matlab based finite element program [13] was used and modified to implement the quadratic αmethod as discussed in the previous section. An Intel Xeon 2.0 GHz processor was used. The material properties were E = 210 GPa, ν = 0.3, y σ = 900 MPa. Displacement load control was applied vertically in the upper face of the squared domain up to 0.01 mm for uniform stretching. A single element, #837, in the vicinity of the slit was monitored for von Mises stress during plastification, as shown in Figure 12a. The von Mises stress contour at the final stage of loading is plotted in Figure   12c. The quadratic time integration combinations of , α β discussed in Section 4 were implemented into the finite element program. The convergence performance and accuracy of the solution were compared with the exact solution of the uniform stretching test. The exact solution was assumed from the result of a case with 100 loading steps using the radial return mapping scheme without iteration (rate format by Simo et al. [4,5]), 30, 12, 8, and 5 loading steps were applied for the investigation. Above 30 loading steps, the accuracy of solution should not differ significantly; however, backward Euler method provided a faster convergence rate from all combination of , α β in the quadratic time integration schemes. As discussed earlier, a smaller stress increment produced more severe positive solution errors mainly localized in the narrow region of the small tangential stress increment. For this reason, the backward Euler method in the linear α -method still showed better performance.
However, for large increments of stress such as the loading step (30, 12, 8 and 5 steps), the backward Euler method lost its stability under the limit of 300 iterations for seeking global Newton-Rapson loops. In contrast, the quadratic time integration scheme with 30, 12, 8 and 5 loading steps showed a stable solution. Table 1 shows that the elapsed computation time and the iteration numbers for each loading step. Figure 13 shows the cross-plots for computational efficiency of the monitored element stress with different time integration schemes. The penalty α value reduces the computation time significantly, as shown in Table 1 and Figure 13. However, this leads to a loss in solution accuracy, Figure 11. Quadratic single-step midpoint methods (SMPT1 and SMPT2): standard α-method and penalty α-method.

2D Metal Panel Problem with Slit-Like Defect
A 1 × 1 mm steel plate with a 0.25 mm slit-like defect in the center was discretized using the Delaunay triangulation scheme with a two-dimensional mesh generator (Persson and Strang, [14]), as shown in Figure 12a,b. A total of 1423 nodes and 2606 CST triangular elements were used. A Matlab based finite element program [13] was used and modified to implement the quadratic α-method as discussed in the previous section. An Intel Xeon 2.0 GHz processor was used. The material properties were E = 210 GPa, ν = 0.3, σ y = 900 MPa. Displacement load control was applied vertically in the upper face of the squared domain up to 0.01 mm for uniform stretching. A single element, #837, in the vicinity of the slit was monitored for von Mises stress during plastification, as shown in Figure 12a. The von Mises stress contour at the final stage of loading is plotted in Figure 12c. The quadratic time integration combinations of α, β discussed in Section 4 were implemented into the finite element program. The convergence performance and accuracy of the solution were compared with the exact solution of the uniform stretching test. The exact solution was assumed from the result of a case with 100 loading steps using the radial return mapping scheme without iteration (rate format by Simo et al. [4,5]), 30, 12, 8, and 5 loading steps were applied for the investigation. Above 30 loading steps, the accuracy of solution should not differ significantly; however, backward Euler method provided a faster convergence rate from all combination of α, β in the quadratic time integration schemes. As discussed earlier, a smaller stress increment produced more severe positive solution errors mainly localized in the narrow region of the small tangential stress increment. For this reason, the backward Euler method in the linear α-method still showed better performance. However, for large increments of stress such as the loading step (30, 12, 8 and 5 steps), the backward Euler method lost its stability under the limit of 300 iterations for seeking global Newton-Rapson loops. In contrast, the quadratic time integration scheme with 30, 12, 8 and 5 loading steps showed a stable solution. Table 1 shows that the elapsed computation time and the iteration numbers for each loading step. Figure 13 shows the cross-plots for computational efficiency of the monitored element stress with different time integration schemes. The penalty α value reduces the computation time significantly, as shown in Table 1 and Figure 13. However, this leads to a loss in solution accuracy, as shown in Figure 14a. When the shape factor was set to {0,1,0} while the α value was varied from 1 to 100, the accuracy of the solution drifted in comparison to the solution from the backward Euler method at 100 time steps (exact solution). Figure 14b shows that light penalty α values such as α = 5, with a shape factor of {1,1,4}/6 yield more time-efficient results without losing solution accuracy. A comparison of computation time efficiencies between quadratic SMPT1 and SMPT2 is also given in Table 1; the quadratic SMPT2 scheme had a slightly higher efficiency.

Conclusions
The quadratic variants of the generalized midpoint rule and return map algorithm for the J 2 von Mises metal plasticity model are examined for the accuracy of deviatoric stress integration of the constitutive equation. Based on the Lagrangian interpolation expansion for the linear α method, the scale factor and shape factor were introduced to a quadratic integration rule for assuming a returning directional tensor from a trial stress onto the final yield surface. As the perfectly plastic model is the only case where the analytical solution is possible, solution accuracies can be compared with those of the exact solutions with the aid of an iso-error map. Since the standard scale factor ranges from 0 to 1, which is similar to the linear α-method, the penalty scale factors greater than 1 were mainly explored to examine the solution accuracies and computational efficiency. A higher value of scale factor above five shows a better computational efficiency but a decreased solution accuracy, especially in the higher plastification zone. A well-balanced scale factor for both computational efficiency and solution accuracy was found to be between one and five. The trade-off scale factor was proposed to be five. The proper shape factor was also proposed to {1,1,4}/6 among the different combinations of weight distribution over a time interval. Since the backward implicit scheme, α = 1 in the linear method shows the best performance in terms of solution accuracy and a computation time. In this quadratic α-method, we tried to explore a penalty scale factor that is higher than 1. Although this approach can be seen in the non-logical sense at a glance, it shows better results in either solution accuracy or computation time. A smaller stress increment produced more severe positive solution errors mainly localized in the narrow region of the small tangential stress increment. For this reason, the backward Euler method in the linear α-method still showed better performance. However, for large increments of stress such as the loading step (30, 12, 8 and 5 steps), the backward Euler method lost its stability under the limit of 300 iterations for seeking global Newton-Rapson loops. In contrast, the quadratic time integration scheme with 30, 12, 8 and 5 loading steps showed a stable solution. In addition, the quadratic SMPT2 scheme had a slightly higher efficiency compared to SMPT1.