## 2 Development of Elastoplastic SW-DVIB Method
### Elastoplastic SW-DVIB Model  
The SW potential is built from two components: two- and three-body potentials. In developing the elastoplastic SW-DVIB model, plasticity is incorporated into the two-body interaction, whereas the three-body interaction remains elastobrittle [63]. Consequently, the common constitutive model of the SW-DVIB retains the following structure, with the nodal force vector and stiffness matrix defined accordingly:  
\[F_{i}=\frac{\partial W}{\partial u_{i}}=\frac{1}{2}\sum\limits_{I}^{N(N-1)} \frac{\partial\Phi_{2}}{\partial l_{I}}\cdot\frac{\partial l_{I}}{\partial u_{ i}}+\frac{1}{2}\sum\limits_{I}^{N(N-1)}\sum\limits_{I}^{N-2}\frac{\partial \Phi_{3}}{\partial\theta_{IJ}}\cdot\frac{\partial\theta_{IJ}}{\partial u_{i}} \tag{6}\]  
\[ \begin{split} K_{ij}=\frac{\partial^{2}W}{\partial u_{i}\partial u _{j}}=\frac{1}{2}\sum\limits_{I}^{N(N-1)}\Bigg{(}\frac{\partial^{2}\Phi_{2}}{ \partial l_{I}^{2}}\cdot\frac{\partial l_{I}}{\partial u_{j}}\cdot\frac{ \partial l_{I}}{\partial u_{i}}+\frac{\partial\Phi_{2}}{\partial l_{I}}\cdot \frac{\partial^{2}l_{I}}{\partial u_{i}\partial u_{j}}\Bigg{)}\\ +\frac{1}{2}\sum\limits_{I}^{N(N-1)}\sum\limits_{I}^{N-2}\Bigg{(} \frac{\partial^{2}\Phi_{3}}{\partial\theta_{IJ}^{2}}\cdot\frac{\partial \theta_{IJ}}{\partial u_{j}}\cdot\frac{\partial\theta_{IJ}}{\partial u_{i}}+ \frac{\partial\Phi_{3}}{\partial\theta_{IJ}}\cdot\frac{\partial^{2}\theta_{ IJ}}{\partial u_{i}\partial u_{j}}\Bigg{)}\end{split}  \tag{7}\]  
[PERSON] et al. [50] incorporated plasticity into the two-body potential using the following equation presented in [65]:  
\[\Phi^{\prime}(l)=f=\left\{ \begin{matrix}\Gamma^{\prime}(l)&\text{if} \quad l\leq l_{b}\\ \frac{(l-l_{c})f_{l}}{l_{b}-l_{c}}&\text{if}\quad l_{b}<l\leq l_{p}\\ \frac{(l-l_{c})f_{l}}{l_{p}-l_{c}}&\text{if}\quad l_{p}<l\leq l_{c}\\ 0&\text{if}\quad l_{c}<l\end{matrix} \right. \tag{8}\]  
Figure 2: Description of bond crack criterion: (**a**) normal bond force vs bond length; (**b**) bond moment \(m_{IJ}\) vs. bond angle. Here, we take Equation (8) as the two-body interaction of the elastoplastic SW-DVIB. The first derivative of the bond potential function relative to the bond length can be expressed as follows:  
\[\Phi^{\prime}{}_{2}(l)=\frac{\partial\Phi_{2}(l)}{\partial l}= \begin{cases}k_{e}( l-l_{0})&l\leq l_{y}&\text{Elastic}\\ k_{e}\big{(}l_{y}-l_{0}\big{)}+k_{p}\big{(}l-l_{y}\big{)}&l_{y}\leq l\leq l_{b} &\text{Plastic}\\ k_{e}(l-l_{0})-\big{(}k_{e}-k_{p}\big{)}\big{(}l_{uP}-l_{y}\big{)}&l_{y}\leq l \leq l_{uP}&\text{Unloading-P}\\ k_{e}\big{(}l_{b}-l_{p}^{*}\big{)}\cdot\Big{(}\frac{l_{f}-l}{l-l_{b}}\big{)} &l_{b}\leq l\leq l_{f}&\text{ fracturinguring}\\ k_{e}\Big{(}l_{b}-l_{p}^{*}\big{)}\cdot\Big{(}\frac{l_{f}-l_{eF}}{l-l_{b}} \Big{)}\cdot\Big{(}\frac{l-l_{p}^{*}}{l_{AF}-l_{p}^{*}}\Big{)}&l_{p}^{*}\leq l \leq l_{uF}&\text{Unloading-F}\\ 0&l\geq l_{f}&\text{Failure}\end{cases}  \tag{9}\]  
\[\Phi^{\prime\prime}{}_{2}(l)= \begin{cases}k_{e}&l\leq l_{y}&\text{Elastic}\\ k_{p}&l_{y}\leq l\leq l_{b}&\text{Plastic}\\ k_{e}&l_{y}\leq l\leq l_{uP}&\text{Unloading}-\text{P}\\ -\frac{k_{e}\big{(}l_{b}-l_{p}^{*}\big{)}}{l_{f}-l_{b}}&l_{b}\leq l\leq l_{f}& \text{ fracturinguring}\\ \cdot\frac{k_{e}\big{(}l_{b}-l_{p}^{*}\big{)}}{l_{AF}-l_{p}^{*}}\cdot\frac{l_{ f}-l_{AF}}{l_{f}-l_{b}}&l_{p}^{*}\leq l\leq l_{uF}&\text{Unloading}-\text{F}\\ 0&l\geq l_{f}&\text{Failure}\end{cases}  \tag{10}\]  
A three-body potential representing linear elastic-brittle behavior was formulated in [65], and is presented in Equation (3). For the current elastoplastic SW-DVIB model, the three-body interaction used is as follows:  
\[\Phi_{3}= \begin{cases}\frac{1}{2}\lambda\big{(}\theta_{IJ}-\theta_{IJ0}\big{)} ^{2}&\text{if}\quad\max\Big{(}\frac{l_{I}}{l_{0}},\frac{l_{J}}{l_{J0}}\Big{)} \leq 1+\varepsilon_{r}\text{ and }\big{|}\theta_{IJ}-\theta_{IJ0}\big{|}\leq \varepsilon_{\theta}\theta_{IJ0}\\ \text{constant}&\text{else}\end{cases}