## 4 Integration Algorithm  
The numerical simulation of cyclic elastoplastic problems requires the constitutive equations in the rate form to be integrated so as to reproduce the effective material response. In this research, the well-known return mapping algorithm (RMA), originally proposed by [PERSON] and [PERSON] [31], is implemented to provide an efficient numerical solution of the cyclic elasto-plastic constitutive model. In such a problem, the variables expressed in the rate form are required to be integrated within the time interval [\(t_{n}\), \(t_{n+1}\)], where \(n\) and \(n\)+1 subscripts are the previous and current load steps, respectively. Using the implicit backward Euler method, the set of constitutive equations are discretized and all the rate quantities are replaced with corresponding incremental unknowns. Assuming that the loading histories illustrated in Figures 1b and 2b are prescribed in the total time interval length [0, \(t_{t}\)] that, in turn, is partitioned into \(N\) sub-intervals by a sequence of discrete time steps \(t_{0}<t_{1}<t_{  
With a little algebra and rearranging the terms in the Chaboche nonlinear kinematic hardening model of Equations (7-8), the following expression is yielded:  
\[\mathbf{X}_{n+1}\ =\ \left(\mathbf{X}_{n}/\left(1+\sum\limits_{n^{\prime}\ =\ 1}^{3}\gamma^{n \prime}\Delta\lambda_{n+1}\right)\right)+\left(\left(\frac{2}{3}\sum\limits_{n ^{\prime}\ =\ 1}^{3}C^{n\prime}\Delta\lambda_{n+1}\right)/\left(1+\Delta\lambda_{n+1}\sum \limits_{n^{\prime}\ =\ 1}^{3}\gamma^{n\prime}\right)\right)\left(\beta_{n+1}/\left\|\beta_{n+1}\right\|\right) \tag{34}\]  
Subtracting Equation (33) from Equation (34), gives:  
\[ \begin{array}{l}\mathbf{\sigma}^{{}^{\prime}}_{n+1}-\mathbf{X}_{n+1}\ =\ \left(\mathbf{\sigma}^{{}^{\prime}}_{n+1}-\mathbf{X}_{n}/\left(1+\Delta \lambda_{n+1}\sum\limits_{n^{\prime}\ =\ 1}^{3}\gamma^{n\prime}\right)\right)-\Delta\lambda_{n+1}\sqrt{3/2}(\beta_{n+1} /\left\|\beta_{n+1}\right\|)(2G+\\ \left(\left(\frac{2}{3}\sum\limits_{n^{\prime}\ =\ 1}^{3}C^{n\prime}\right)/ \left(1+\Delta\lambda_{n+1}\sum\limits_{n^{\prime}\ =\ 1}^{3}\gamma^{n\prime}\right)\right)\end{array}  \tag{35}\]  
It is evident from Equation (35) that the following identity is established:  
\[\frac{\mathbf{\sigma}^{{}^{\prime}}_{n+1}-\mathbf{X}_{n+1}}{\|\mathbf{\sigma} ^{{}^{\prime}}_{n+1}-\mathbf{X}_{n+1}\|}\ =\ \frac{\mathbf{j}_{n+1}^{trial}}{\|\mathbf{j}_{n+1}^{trial}\|} \tag{36}\]  
where  
\[\mathbf{j}_{n+1}^{trial}\ =\ \sqrt{3/2}\|\mathbf{\sigma}^{{}^{\prime}}_{n+1}- \mathbf{X}_{n}/\left(1+\Delta\lambda_{n+1}\sum\limits_{n^{\prime}\ =\ 1}^{3}\gamma^{n\prime}\right)\| \tag{37}\]  
Using the relation of Equation (36), the following expression for Equation (35) is obtained:  
\[ \begin{array}{l}\mathbf{\sigma}^{{}^{\prime}}_{n+1}-\mathbf{X}_{n+1}\ =\ \left(\mathbf{\sigma}^{{}^{\prime}}_{n+1}-\mathbf{X}_{n}/\left(1+\sum \limits_{n^{\prime}\ =\ 1}^{3}\gamma^{n\prime}\Delta\lambda_{n+1}\right)\right)\left(1-\Delta \lambda_{n+1}\sqrt{3/2}\left(2G+\left(\frac{2}{3}\sum\limits_{n^{\prime}\ =\ 1}^{3}C^{n\prime}\right.\right.\\ \left.\left./1+\sum\limits_{n^{\prime}\ =\ 1}^{3}\gamma^{n\prime}\ \Delta\lambda_{n+1} \right)\right)/\left\|\mathbf{\sigma}^{{}^{\prime}}_{n+1}-\mathbf{X}_{n}/ \left(1+\Delta\lambda_{n+1}\sum\limits_{n^{\prime}\ =\ 1}^{3}\gamma^{n\prime}\right)\right\|\right)\end{array}  \tag{38}\]  
Finally, using Equation (38), the yield criterion, i.e., \(F_{n+1}\) = 0, is written as:  
\[F_{n+1}\ =\ \sqrt{3/2}\|\mathbf{\sigma}^{{}^{\prime}}_{n+1}-\mathbf{X}_{n}/ \left(1+\sum\limits_{n^{\prime}\ =\ 1}^{3}\gamma^{n\prime}\Delta\lambda_{n+1}\right)\|-\Delta \lambda_{n+1}(3G\ +\left(\sum\limits_{n^{\prime}\ =\ 1}^{3}C^{n\prime}/1+\sum\limits_{n^{\prime}\ =\ 1}^{3}\gamma^{n\prime}\ \Delta\lambda_{n+1} \right))-\sigma_{y}\ =\ 0 \tag{39}\]  
Equation (39) is effectively solved by the iterative method of Newton-Raphson procedure with \(p\)th iterations and designated tolerance (tol = \(10^{-6}\)) [25] to obtain the incremental plastic multiplier \(\Delta\lambda_{n+1}\).