## 3 Numerical Methods
### Immersed Boundary Method  
For incompressible viscous flow of a Newtonian fluid with constant density and kinematic viscosity \(\upsilon\), using tensor notation the governing equations in Cartesian coordinates can be given as:  
\[\frac{\partial u_{i}}{\partial t}+\frac{\partial\big{(}u_{i}u_{j}\big{)}}{ \partial x_{j}}=-\frac{1}{\rho}\frac{\partial p}{\partial x_{i}}+\upsilon\frac {\partial^{2}u_{i}}{\partial x_{j}\partial x_{j}}+f_{i} \tag{15}\]  
Figure 3: The value of the level set (LS) function \(\phi\): (**a**) before the re-initialization process; (**b**) after the re-initialization process. \[\frac{\partial u_{i}}{\partial x_{i}}=0 \tag{16}\]  
where \(x_{i}(i=1,\ 2,\ 3)\) represents the Cartesian coordinates \((x,y,z)\); \(u_{i}\) (or equivalently \(u\), \(v\) and \(w\)) represents the velocity components in the corresponding directions; \(t\) is the time; \(p\) is the pressure; and \(f_{i}\) is the external force, representing the effects of solid body on the fluid. For the laminar flow, Equations (15) and (16) are directly integrated in space and time without introducing any model. Large eddy simulation (LES) approach is used for turbulent flow, and the dynamic Samgorinsky model [23] is adopted to calculate the subgrid-scale (SGS) stress. The staggered grid is used in our in-house developed code. The velocity components \(u\), \(v\), and \(w\) are defined at centers of cell surfaces, and all other variables, i.e., pressure \(p\) and LS function \(\phi\) are defined at the center of the grid cell. The finite difference method is used for spatial discretization of the Navier-Stokes equations, and the second order Runge-Kutta (RK2) method is used for time integration. In each sub-step of the RK2 method, a fractional-step method of [PERSON] and [PERSON] [24] is used, and the projection method is utilized to ensure the divergence-free requirement for incompressible flow. The algorithm can be referred to [PERSON]'s work [25]. The Krylov-based multigrid solver from the PETSc library is used to solve the pressure Poisson equation. In Cartesian IB method, the Cartesian grid is not required to be precisely conformed to the fluid-solid interface. A forcing term \(\overrightarrow{f_{b}}\) is introduced into the momentum equation to represent the effect of the solid body on the flow field, and it is applied on the forcing points. At first, we advance the simulation without considering the boundary effects. Then the velocity \(\overline{u}_{b}\) on the forcing point is interpolated by the velocity of the surrounding fluid points (neighbouring to the forcing points on the \(x\)-grid line, \(y\)-grid line, or along the diagonal) and the solid point (the projection of forcing point). In present study, the solid points at the interface are found through the gradient operators of the LS function, and the velocities are decomposed into the tangential and normal components. As shown in Figure 4, the bilinear interpolation scheme is used to reconstruct the velocity \(u_{b}\) at the forcing point. For the variable \(\alpha\) in the 2D space, the interpolation coefficients for the linear interpolation can be written as:  
\[\alpha=b_{1}+b_{2}x+b_{3}y \tag{17}\]  
The coefficients \(b_{1}\), \(b_{2}\), and \(b_{3}\) in Equation (17) can be found by solving:  
\[\left[ \begin{array}{c}b_{1}\\ b_{2}\\ b_{3}\end{array} \right]=A^{-1}\left[ \begin{array}{c}\alpha_{1}\\ \alpha_{2}\\ \alpha_{3}\end{array} \right]=\left[ \begin{array}{ccc}1&x_{1}&y_{1}\\ 1&x_{2}&y_{2}\\ 1&x_{3}&y_{3}\end{array} \right]^{-1}\left[ \begin{array}{c}\alpha_{1}\\ \alpha_{2}\\ \alpha_{3}\end{array} \right] \tag{18}\]  
where \((x_{i},y_{i},z_{i})\) with \(i=1,\ 2\) and \(3\) in matrix \(A\) are the point coordinates of the interpolation stencil. For the case of swimming fish, the inversion of matrix \(A\) is calculated at each time step, because of the deformed and moving fish body.