## 3 Calculation Methods  
Computational fluid mechanics (CFD) is a method used to mathematically solve the governing equation in flow fields and to calculate series of nodes numerically in finite element meshes. For incompressible fluids, it is based on three basic laws: the law of conservation of mass, the law of conservation of momentum, and the law of conservation of energy, as shown in Equations (1) and (2). The energy conservation equation was not considered here because energy exchange was not involved in this scenario. \[\frac{\partial\rho}{\partial t}+\
abla\cdot(\rho\overset{\rightarrow}{\mathbf{u}})=0 \tag{1}\]  
\[\frac{\partial(\rho\overset{\rightarrow}{\mathbf{u}})}{\partial t}+\
abla\cdot( \rho\overset{\rightarrow}{\mathbf{u}}\overset{\rightarrow}{\mathbf{u}})=-\
abla p+ \
abla\cdot(\overset{\rightarrow}{\mathbf{\tau}})+\rho\overset{\rightarrow}{ \mathbf{g}}-\overset{\rightarrow}{\mathbf{F}} \tag{2}\]  
where \(\rho\), \(t\), and \(\overset{\rightarrow}{\mathbf{u}}\) represent the fluid density, time, and velocity vector of fluid media, respectively; subscripts x, y, and z denote the projection of the velocity vector in the x-, y-, and z-axes, respectively; \(\
abla\), \(p\), and \(\overset{\rightarrow}{\mathbf{g}}\) represent the Hamiltonian operator, the pressure of fluid media, and the acceleration vector, respectively; \(\overset{\rightarrow}{\mathbf{F}}\) refers to the interaction force between particles and fluids; and \(\overset{\rightarrow}{\mathbf{\tau}}\) is the stress tensor, which indicates the viscous force exerted by the fluids. Figure 10: Stress diagram for the rock plug. In the DEM, the particle motion follows [PERSON]'s second law of motion and involves both rotary and translational motion. The motion equations of particles (momentum conservation equation) are shown in Equations (3) and (4). \[m_{i}\frac{\partial\overrightarrow{v}_{i}}{\partial t}=\sum_{j}\overrightarrow {F}_{ij}+\overrightarrow{F}_{g}+\overrightarrow{F}_{f} \tag{3}\]  
\[I_{i}\frac{\partial\overrightarrow{\omega}_{i}}{\partial t}=\sum_{j} \overrightarrow{T}_{ij}+T_{f} \tag{4}\]  
where \(v\), \(\omega\), and \(m\) represent the velocity, angular velocity, and mass of particles, respectively; subscripts \(i\) and \(j\) indicate two particles in contact with each other; \(F_{ij}\) is the contact force between particles \(i\) and \(j\); \(F_{g}\) is the gravitational force between particles; \(F_{f}\) represents the interaction force between fluids and particles; \(I\) denotes the rotational inertia of particles; \(T\) is the moment between particles; and \(T_{f}\) indicates the torque induced by the fluids. For the interaction force between particles in the DEM, the Hertz-Mindlin contact model, considering rolling friction, was used here (Figure 11). The forces between two particles in contact with each other can be expressed using Equations (5)-(8) [34]. \[F_{n}=S_{n}\delta_{n} \tag{5}\]  
\[F_{t}=-S_{t}\delta_{t} \tag{6}\]  
\[F_{n}^{d}=-2\sqrt{\frac{5}{6}}\beta\sqrt{S_{n}m^{*}}v_{n}^{\overrightarrow{rel}} \tag{7}\]  
\[F_{t}^{d}=-2\sqrt{\frac{5}{6}}\beta\sqrt{S_{t}m^{*}}v_{t}^{\overrightarrow{rel}} \tag{8}\]  
where \(F_{n}\) and \(F_{t}\) denote the normal and tangential contact forces, respectively; \(\delta n\) and \(\delta t\) represent the normal and tangential overlaps, as shown in Figure 11; \(F_{n}^{d}\) and \(F_{t}^{d}\) represent the normal and tangential damping forces; \(v_{n}\) and \(v_{t}\) indicate the projection of the relative velocity in the normal and tangential directions; and \(m^{*}\) represents the equivalent mass. The damping coefficient \(\beta\), normal stiffness \(Sn\), and tangential stiffness \(St\) can be calculated using Equations (9)-(11).