## 2 Numerical Model
### Solution of Fluid Sub-Problem  
The following projection scheme is used to further decouple (7) and (8) as follows:  
1. Velocity Helmholtz equation \[\frac{\mathbf{u}^{\#}-\mathbf{u}^{n}}{\Delta t}-\frac{1}{2 Re}\
abla^{2}\mathbf{u}^{\#}= \frac{1}{2 Re}\
abla^{2}\mathbf{u}^{n}+\mathbf{H}^{n},\] (13) where \(\mathbf{H}^{n}=-\
abla p^{n}-\frac{1}{2}[3(\mathbf{u}\cdot\
abla\mathbf{u})^{n}-(\mathbf{u} \cdot\
abla\mathbf{u})^{n-1}]+\mathbf{\lambda}^{n}\). 2. Pressure Poison equation \[\
abla^{2}\varphi=\frac{\
abla\cdot\mathbf{u}^{\#}}{\Delta t}.\] (14)
3. Velocity and pressure correction \[\frac{\mathbf{u}^{*}-\mathbf{u}^{\#}}{\Delta t}=-\
abla\varphi.\] (15)  
\[p^{n+1}=p^{n}+\varphi. \tag{16}\]  
The velocity Helmholtz Equation (13) is solved with the Douglas Gunn alternating-direction-implicit (ADI) scheme [36]:  
\[\frac{\mathbf{u}^{n*}-\mathbf{u}^{n}}{\Delta t}-\frac{1}{2 Re}\
abla_{x}^{2}\mathbf{u}^{n*} =\frac{1}{2 Re}\
abla_{x}^{2}\mathbf{u}^{n}+\frac{1}{Re}(\
abla_{y}^{2}\mathbf{u}^{n}+ \
abla_{z}^{2}\mathbf{u}^{n})+\mathbf{H}^{n}, \tag{17}\]  
\[\frac{\mathbf{u}^{n**}-\mathbf{u}^{n*}}{\Delta t}-\frac{1}{2 Re}\
abla_{y}^{2}\mathbf{u}^{n** }=-\frac{1}{2 Re}\
abla_{y}^{2}\mathbf{u}^{n}, \tag{18}\]  
\[\frac{\mathbf{u}^{\#}-\mathbf{u}^{n**}}{\Delta t}-\frac{1}{2 Re}\
abla_{z}^{2}\mathbf{u}^{ \#}=-\frac{1}{2 Re}\
abla_{z}^{2}\mathbf{u}^{n}. \tag{19}\]  
In the above equations, \(\
abla_{x}^{2}\), \(\
abla_{y}^{2}\), and \(\
abla_{z}^{2}\) represent the second-order derivative in the \(x\), \(y\), and \(z\) directions, respectively. In the absence of \(\mathbf{H}^{n}\), the [PERSON] scheme above is second-order accurate and unconditionally stable. A homogeneous grid is used. For our previous FD method, a half-staggered grid was employed, which means that the velocity components are collocated but the pressure nodes are staggered with the velocity nodes. Somehow this scheme was found unstable for the present problem, and thus we adopt the fully staggered grid. The velocity derivatives in the convection and diffusion terms are discretized with a second-order central difference scheme. The most time-consuming step for the incompressible flow is the solution of the pressure Poison Equation (14). Here, we develop a combined FFT and tri-diagonal-system fast solver for the present problem. The pressure is defined at the cell center and the discretized pressure equation at the pressure node \(i,j,k\) can be written as follows:  
\[p_{i+1,j,k}+p_{i-1,j,k}+p_{i,j+1,k}+p_{i,j-1,k}+p_{i,j,k+1}+p_{i,j,k-1}-6p_{i,j,k}=f_{i,j,k}, \tag{20}\]  
for \(i=0,M-1;j=0,N-1;k=0,L-1\). Here \(M\), \(N\), and \(L\) represent the cell numbers of the computational domain in the \(x\), \(y\), and \(z\) directions, respectively. Since the pressure is periodic in the \(y\) direction and its normal gradient on the side wall is zero, we perform an FFT in the \(y\) direction and a fast cosine transformation in the spanwise (z) direction, and the pressure in the physical space \(p_{i,j,k}\) can be expanded with its spectrum \(P_{i,n,l}\):  
\[p_{i,j,k}=\frac{2}{NL}\sum_{n=0}^{N-1}\sum_{l=0}^{L-1}P_{i,n,l}e^{-\frac{2\pi i \pi}{N}}\cos\frac{\pi(k+\frac{1}{2})l}{L}, \tag{21}\]  
where \(i^{\prime}\) denotes the imaginary number \(i^{\prime}=\sqrt{-1}\). Expanding all terms in (20) with (21), one can obtain the following equation for the pressure spectrum:  
\[P_{i+1,n,l}+P_{i-1,n,l}+P_{i,n,l}(2\cos\frac{2\pi n}{N}+2\cos\frac{\pi l}{L}-6) =F_{i,n,l}, \tag{22}\]  
which is a tri-diagonal system and holds for \(i=1,M-2\). For the shifted periodic boundary condition, the equation for \(i=0\) and \(M-1\) needs to be treated specifically. We assume that the row shift satisfies  
\[\
abla^{2}\varphi=\frac{\
abla\cdot\mathbf{u}^{\#}}{\Delta t}.\]0  
where \(h\) denotes the mesh size.