# A Semi-Explicit Multi-Step Method for Solving Incompressible Navier-Stokes Equations
## 2 Governing Equations
### Space Discretization  
The Galerkin weak form of (8) reads:  
\[\left(\boldsymbol{\omega},\rho\frac{\partial\mathbf{u}}{\partial t}\right)+( \boldsymbol{\omega},\mathcal{A}_{1})+(q,\mathcal{A}_{2})=\langle\boldsymbol{ \omega},\mathbf{f}\rangle, \tag{9}\]  
where \(\boldsymbol{\omega}\) and \(q\) are the velocity and pressure test functions. In the notation adopted, \((\cdot,\cdot)\) and \(\langle\cdot,\cdot\rangle\) stand for a bilinear and linear forms, respectively. In order to obtain stable solutions for convection-dominated flows, the equations must be stabilized. Additionally, pressure stabilization is necessary for equal order velocity-pressure interpolations, as they do not satisfy the compatibility condition [27]. The Galerkin/least squares (GLS) stabilization method allows one to circumvent this problem by summing \(\sum_{n}\int_{\Omega}\mathcal{A}(\boldsymbol{\omega},q)\cdot(\mathcal{F}- \mathcal{A}(\mathbf{u},p)-\frac{\partial\mathbf{u}}{\partial t})\) to the original weak form of the problem [28]. This means summing the product of the residual with the Navier-Stokes operator \(\mathcal{A}(\boldsymbol{\omega},q)\) (this operator may be viewed as a changed weight function). We note that for linear velocity-pressure elements, the GLS method coincides with the more recent (and less diffusive for a general case) algebraic sub-grid scale (ASGS) method [29], the latter being defined based on the adjoint operator \(\mathcal{A}^{*}\) as a weight function for the residual. Summing up the Galerkin weak form with the stabilization terms results in:  
\[ \begin{array}{l}\rho(\mathbf{\omega},\frac{\partial\mathbf{u}}{\partial t})+\rho(\mathbf{ \omega},\mathbf{u}\cdot\
abla\mathbf{u})+(\mathbf{\omega},\
abla p)-\rho\
u(\mathbf{\omega},\Delta \mathbf{u})+(q,\
abla\cdot\mathbf{u})+\\ \sum_{el}\tau\int_{\Omega_{d}}\left(\
abla q+\rho\mathbf{\omega}\cdot\
abla\mathbf{ \omega}-\rho\
u\Delta\mathbf{\omega}\right)(\mathbf{f}-\rho\frac{\partial\mathbf{u}}{ \partial t}-\rho\mathbf{u}\cdot\
abla\mathbf{u}+\rho\
u\Delta\mathbf{u}-\
abla p), \end{array}  \tag{10}\]  
where the stabilization coefficient \(\tau\) is defined as (see, e.g., [28]):  
\[\tau=\frac{2||\mathbf{u}||}{h}+\frac{4\
u}{h^{2}}, \tag{11}\]  
where \(h\) is the element size. Integrating the viscous term by parts and rearranging the terms (noting that in Equation (10), products containing the test function \(\mathbf{\omega}\) contribute to the momentum equation, while the ones containing \(q\) contribute to the mass conservation equation), the modified linear momentum and mass conservation equations are obtained:  
\[ \begin{array}{l}\rho(\mathbf{\omega},\frac{\partial\mathbf{u}}{\partial t})+\rho( \mathbf{\omega},\mathbf{u}\cdot\
abla\mathbf{u})+\rho\
u(\
abla\mathbf{\omega},\
abla\mathbf{u})- (\
abla\cdot\mathbf{\omega},p)+\\ \sum_{el}\tau\int_{\Omega_{d}}\left(\mathbf{\omega}\cdot\
abla\mathbf{\omega}\right) \left(\mathbf{f}-\rho\frac{\partial\mathbf{u}}{\partial t}-\rho\mathbf{u}\cdot \
abla\mathbf{u}-\
abla p\right)=\langle\mathbf{f},\mathbf{\omega}\rangle,\end{array}  \tag{12}\]  
\[ \begin{array}{l}(q,\
abla\cdot\mathbf{u})+\sum_{el}\tau\int_{\Omega_{d}}\
abla q \left(\mathbf{f}-\rho\frac{\partial\mathbf{u}}{\partial t}-\rho\mathbf{u}\cdot \
abla\mathbf{u}-\
abla p\right)=0.\end{array}  \tag{13}\]  
Prior to establishing the space-discretized stabilized governing system, we define the matrices and vectors corresponding to each term of Galerkin weak form:  
\[ \begin{array}{l}(\mathbf{\omega},\mathbf{u})\longrightarrow\mathbf{M}\bar{\mathbf{u}}, \end{array}  \tag{14}\]  
\[ \begin{array}{l}(\
abla\cdot\mathbf{\omega},p)\longrightarrow\mathbf{G}\bar{\mathbf{ p}},\end{array}