# Bearing Capacity of Karst Cave Roof under File Foundation Load Using Limit Analysis
## 3 Mathematical Model
### Flow Law Constraints  
For ideal rigid plastic materials, the strain component obtained from the deformation compatibility condition should be equal to the plastic strain rate component obtained from the associated flow rule and yield condition. Since the strain within each element is continuous, the flow law constraint can be expressed as  
\[\left\{ \begin{array}{l}\mathbf{B}^{e}\mathbf{u}^{e}-\lambda^{e}\
abla f^{e}( \mathbf{\sigma}^{e})=0;(e=1,\ldots,\text{N}_{1})\\ \lambda^{e}f^{e}(\mathbf{\sigma}^{e})=0;(e=1,\ldots,\text{N}_{1})\\ \lambda^{e}\geq 0;(e=1,\ldots,\text{N}_{1})\end{array} \right. \tag{5}\]  
where \(\text{N}_{1}\) is the number of elements, \(\mathbf{B}^{e}\) is the compatibility matrix of the element, \(\mathbf{u}^{e}\) is the velocity vector of the element, \(\lambda^{e}\) is the non-negative plastic multiplier of the element, \(f^{e}(\sigma)\) is the yield function, and the specific expressions for \(\mathbf{B}^{e}\) and \(\mathbf{u}^{e}\) are as follows  
\[\left\{ \begin{array}{l}\mathbf{B}^{e}=\left[\mathbf{B}_{1},\mathbf{B}_{2}, \mathbf{B}_{3},\mathbf{B}_{4}\right]\\ \\ \mathbf{B}_{i}=\frac{1}{6\
abla^{e}}\left[\begin{array}{cccccc}b_{i}&0&0&c_{ i}&0&d_{i}\\ 0&c_{i}&0&b_{i}&d_{i}&0\\ 0&0&d_{i}&0&c_{i}&b_{i}\\ \end{array}\right]^{\text{T}}(i=1,2,3,4)\\ \\ \mathbf{u}^{e}=\left[\mathbf{u}_{1}^{e},\mathbf{u}_{2}^{e},\mathbf{u}_{3}^{e},\mathbf{u}_{4}^{e}\right] \\ \
abla=\left\{\partial/\partial\sigma_{xx},\partial/\partial\sigma_{yy}, \partial/\partial\sigma_{zz},\partial/\partial\gamma_{xy},\partial/\partial \gamma_{yz},\partial/\partial\gamma_{zz}\right\}^{\text{T}}\end{array} \right. \tag{6}\]  
Figure 1: Upper-bound tetrahedral element: (**a**) Velocity and Stress of tetrahedral element, (**b**) Velocity discontinuity. where \(\mathbf{B}_{i}(i=1,2,3,4)\) is the compatibility matrix of the node \(i\) of the element. On the global kinematically admissible velocity field, the flow rule can be expressed as  
\[\left\{ \begin{array}{l}\mathbf{B}\boldsymbol{u}-\sum\limits_{e=1}^{E}\lambda^ {e}\
abla f^{e}(\boldsymbol{\sigma})=0;(e=1,\ldots,\mathrm{N}_{1})\\ \lambda^{e}\
abla f^{e}(\boldsymbol{\sigma})=0;(e=1,\ldots,\mathrm{N}_{1})\\ \lambda^{e}\geq 0;(e=1,\ldots,\mathrm{N}_{1})\end{array} \right. \tag{7}\]  
where \(\mathbf{B}\) is the global compatibility, \(\boldsymbol{u}\) is the global continuous velocity vector, and \(\boldsymbol{\sigma}\) is the global stress vector. In order to satisfy the flow law for a velocity discontinuity occurring on a common surface, it is imperative to apply constraints on the nodes associated with the discontinuity. The introduction of a node velocity variable in the local coordinate system between any adjacent nodes allows for the description of the velocity discontinuity, thereby guaranteeing that both tangential and normal velocities traverse the discontinuity to satisfy the flow law. Within the context of the Mohr-Coulomb yield criterion, the correlation between the jump in normal velocity and the jump in tangential velocity is established. \[\Delta u_{n}= |\Delta u_{t}|\tan\varphi \tag{8}\]  
where \(\Delta u_{n}\) and \(\Delta u_{t}\) are the normal and tangential velocity jump, respectively, \(\varphi\) is the internal friction angle of the material. The tangential velocity jump vector \(\Delta\boldsymbol{u^{\prime}}_{t}^{lm}\) of each pair of nodes \((l,m)\) on the discontinuity is expressed by the sum of two non-negative velocity vectors \(\boldsymbol{v}^{+lm}\) and \(\boldsymbol{v}^{-lm}\). \[\Delta\boldsymbol{u^{\prime}}_{t}^{lm}=\boldsymbol{v}^{+lm}-\boldsymbol{v}^{ -lm} \tag{9}\]  
where \(\boldsymbol{v}^{+lm}=\left\{v_{x^{\prime}}^{+lm},v_{y^{\prime}}^{+lm}\right\} ^{\mathrm{T}}\), \(\boldsymbol{v}^{-lm}=\left\{v_{x^{\prime}}^{-lm},v_{y^{\prime}}^{-lm}\right\} ^{\mathrm{T}}\), \(\Delta\boldsymbol{u^{\prime}}_{t}^{lm}=\left\{\Delta u^{\prime lm},\Delta u^{ \prime lm}_{y^{\prime}}\right\}^{\mathrm{T}}\).