## 3 Explicit Finite Element Method for the _u-p_ Formulation
### Spatial Discretization  
The dynamic responses of node \(k\) only are associated with the responses of the neighboring nodes. On that basis, taking the 2D plane strain problem as an example, the whole computational domain is discretized into local node systems for all nodes [44], as shown in Figure 1. For node \(k_{1}\) at the corner of the computational domain, the four nodes of element 1 adjacent to node \(k\) are renumbered as nodes 1 to 4; similarly, for the node \(k_{2}\) on the boundary and the node \(k_{3}\) located inside the computational domain, the local node systems are renumbered by the method shown in Figure 1. In the computational domain, for the local node system composed of any node, an individual finite element equation of the local node system is constructed and is computed in time domain. The dynamic responses for the whole computational domain are obtained by cycling all the nodes.  
According to Figure 1, it is assumed that there are \(n\) adjacent elements for node \(k\), and each element has four nodes. Based on Equations (6) and (7), the finite element equations of the local node system of node \(k\) are expressed as follows:  
\[\sum_{L=1}^{n}\sum_{j=1}^{4}M_{kj}^{L}\bar{u}_{j}^{L}+\sum_{L=1}^{n}\sum_{j=1} ^{4}K_{kj}^{L}u_{j}^{L}-\sum_{L=1}^{n}\sum_{j=1}^{4}Q_{kj}^{L}p_{j}^{L}=\sum_{ L=1}^{n}\mathbf{f}_{k,\mu}^{L} \tag{8}\]  
\[\sum_{L=1}^{n}\sum_{j=1}^{4}S_{kj}^{L}p_{j}^{L}+\sum_{L=1}^{n}\sum_{j=1}^{4}H_ {ij}^{L}p_{j}^{L}+\sum_{L=1}^{n}\sum_{j=1}^{4}\left(Q_{kj}^{L}\right)^{T}\bar{ u}_{j}^{L}=\sum_{L=1}^{n}f_{k,\mu}^{L} \tag{9}\]  
where \(L\) represents the influence of the adjacent element \(L\) on node \(k\), and \(j=1-4\) represent the influence of the four nodes of the element \(L\) on node \(k\). The matrices in the above equations are defined as follows:  
\[\mathbf{M}_{ij}^{L} =\int_{-1}^{1}\int_{-1}^{1}\rho\mathbf{N}_{ui}\mathbf{N}_{uj}|J|d\xi d\eta\] \[\mathbf{K}_{ij}^{L} =\int_{-1}^{1}\int_{-1}^{1}\mathbf{N}_{ui}^{L}\mathbf{T}\mathbf{D}\mathbf{L}\mathbf{N} _{uj}|J|d\xi d\eta\] \[\mathbf{Q}_{ij}^{L} =\int_{-1}^{1}\int_{-1}^{1}\mathbf{N}_{ui}^{T}\
abla p_{ij}|J|d\xi d\eta \tag{10}\] \[H_{ij}^{L} =\int_{-1}^{1}\int_{-1}^{1}\left(\
abla p_{ji}\right)^{T}\bar{K} \left(\
abla p_{ij}\right)|J|d\xi d\eta\]  
Figure 1: Local node systems.  
where \(M\) is the mass matrix of saturated porous media, \(K\) is the stiffness matrix of the solid phase, \(Q\) is the coupling matrix, \(S\) is the compressibility matrix of the fluid phase, \(H\) is the permeability matrix of the fluid phase, \(f_{u}\) and \(f_{q}\) are the external load of and the flow rate, respectively, \(T_{s}\) is the distributed force on the boundary of element \(L\) in the \(x\) and \(y\) directions, and \(q_{s}\) is the flow rate perpendicular to the boundary of element \(L\).  
In fact, Equations (8) and (9) are an alternative expression of Equations (6) and (7). Equations (6) and (7) are the global equations established by assembling the coefficient matrices of the global computational domain. In comparison, Equations (8) and (9) are the equations established for the local node system. The dynamic responses of the whole domain are obtained by cycling all nodes in the computational domain and solving the equation of local node system of each node. Therefore, for the dynamic problem with large degree of freedom, the dimension of equation of the local node system is very small and computational efficiency is higher.