## 2 Finite Element Formulations and Developed Procedure
### Fiber Beam Element Model with Nonlinear Interfacial Connection
For the fiber numbered \(k\) within the cross-section of the element, the specified data includes fiber centroid coordinates (\(y_{k}\), \(z_{k}\)), fiber area \(A_{k}\), fiber region flag \(\varphi_{k}^{s}\) and \(\varphi_{k}^{c}\), fiber stress \(\sigma_{k}\) and fiber tangent modulus \(E_{k}\), as shown in Figure 2b. For fibers in the concrete slab (\((y_{k},\,\mathrm{z}_{k})\in\Omega_{k}\)), \(\varphi_{k}^{s}=0\) and \(\varphi_{k}^{c}=1\); otherwise (\((y_{k},\,\mathrm{z}_{k})\in\Omega_{k}\)), \(\varphi_{k}^{s}=1\) and \(\varphi_{k}^{c}=0\). Then, the axial strain of fiber \(k\) can be yielded by Equation (7). The slip displacement at the location \(x\) can be denoted as Equation (8). \[\begin{array}{llll}\varepsilon(x)&=\big{[}\ \varphi_{s}^{k}(y,z)\ \\[\mathbf{F} =\mathbf{P}^{T}\mathbf{q} \tag{11}\] \[\mathbf{q} =\int_{0}^{L}\mathbf{B}^{T}(x)\left[\int\limits_{A}\mathbf{Q}^{T}(y,z)\sigma_{k}dA+\mathbf{h}^{T}f_{s}\right]dx=\int_{0}^{L}\mathbf{B}^{T}(x) \mathbf{r}_{s}(x)dx\]  
The element stiffness matrix \(\mathbf{K}\) can be obtained by the partial derivative of Equation (11) to \(\mathbf{u}\), as shown in Equation (12). \(\mathbf{k}^{s}\) is the geometric stiffness matrix and \(\mathbf{k}^{m}\) is the material stiffness matrix in the basic coordinate system. The derived expression of \(\mathbf{k}^{m}\) is shown in Equation (13), in which \(\mathbf{k}_{s}\) is the sectional stiffness matrix. Based on the idea of the fiber section method, the integration style of \(\mathbf{r}_{s}\) and \(\mathbf{k}_{s}\) can be calculated by the algebraic summary of the total fibers' contribution, as shown in Equations (14) and (15). The integrations of the \(\mathbf{k}^{m}\) and \(\mathbf{q}\) of the element are implemented using the Gauss-Lobatto method. \[\mathbf{K} =\frac{\partial\mathbf{F}}{\partial\mathbf{u}}=\frac{\partial}{ \partial\mathbf{u}}\big{(}\mathbf{P}^{T}\mathbf{q}\big{)}\] (12) \[=\frac{\partial\mathbf{P}}{\partial\mathbf{u}}\mathbf{q}+ \mathbf{P}^{T}\frac{\partial\mathbf{q}}{\partial\mathbf{u}}\frac{\partial  
Assuming that the initial length of the element was \(L_{0}\) which is the sum length of all segments, after deformation, the deformed element length can be expressed by the displacement vector \(\mathbf{U}_{t}\), as shown in Equation (17),  
\[\epsilon^{t}=\frac{\mathbf{U}_{t}^{T}\mathbf{P}_{t}^{T}\mathbf{\gamma}_{t}\left[ \mathbf{P}_{t}\mathbf{U}_{t}+\mathbf{R}^{T}\mathbf{L}_{0}\right]}{L_{0}} \tag{17}\]  
Figure 4: Stress and modulus update flow during the iterative calculation process. Figure 5: Configuration of the multi-node slipping cable element. in which,  
\[\mathbf{P}_{t}=\left[ \begin{array}{ccccc}-1&0&1&&&\\ &-1&0&1&&\\ &&-1&0&1&\\ &&\ddots&\ddots&\ddots&\\ &&&-1&0&1\end{array} \right]_{(2n_{p}-2)\times(2n_{p})},\] \[\mathbf{r}_{j}=\left[ \begin{array}{ccccc}\cos X&\cos Y\\ -\cos Y&\cos X\end{array} \right],\;\cos X=\frac{x_{j+1,0}-x_{j,0}}{l_{j,0}},\; \cos Y=\frac{y_{j+1,0}-y_{j,0}}{l_{j,0}},\] \[\Upsilon_{t}=diag\Big{(}\left[ \begin{array}{cc}l_{1,t}^{-1} \mathbf{I}_{2,\times 2}&l_{j,t}^{-1}\mathbf{I}_{2,\times 2}&\cdots&l_{n_{p}-1,t}^{-1} \mathbf{I}_{2,\times 2}&\end{array} \right]\Big{)},\mathbf{L}_{0}=\left[ \begin{array}{ cc}\mathbf{I}_{1,0}^{T}&\mathbf{I}_{j,0}^{T}&\cdots&\mathbf{I}_{n_{p}-1,0}^{T} \end{array} \right]^{T},\mathbf{L}_{j,0}=\left[ \begin{array}{cc}l_{j,0}&0 \end{array} \right]^{T}.\]  
By the virtual work equation and variation principle, the resisting force vector and stiffness matrix of the element can be derived and expressed in Equations (18) and (19), respectively. \(A_{p}\) is the area of the external tendon, \(\sigma_{p}\) and \(E_{p}\) denote the tendon stress and tangent modulus at time \(t\), respectively, which is calculated according to the material constitute model. \(\mathbf{k}_{t}^{S}\) and \(\mathbf{k}_{t}^{m}\) denote the geometric and material stiffness matrix, respectively.