# An Improved Galerkin Framework for Solving Unsteady High-Reynolds Navier-Stokes Equations
## 3 Numerical Methods  
### Taylor-Hood Element Discretization  
In this paper, the Taylor-Hood element space is considered such that  
\[\mathbf{v}_{h}\in V_{h}=\{\mathbf{v}_{h}\in V|\mathbf{v}_{h}|_{K}\in P_{2}(K),\forall K\in\mathcal{T}_{h}\},\mathbf{q}_{h}\in Q_{h}=\{\mathbf{q}_{h}\in Q| \mathbf{q}_{h}|_{K}\in P_{1}(K),\forall K\in\mathcal{T}_{h}\},\]and select two sets of basis functions \(\{\phi_{i}\}_{i=1}^{N_{u}}\) and \(\{\psi_{i}\}_{i=1}^{N_{p}}\), respectively. Since any element in a finite-dimensional space can be expressed as a linear combination of its basis, the velocity field and pressure field can be approximated as  
\[\mathbf{u}_{h}=\sum_{i=1}^{N_{h}}\mathbf{u}_{i}\phi_{i}(\mathbf{x}),\quad p_{h}= \sum_{j=1}^{N_{p}}p_{j}\psi_{j}(\mathbf{x}),\]  
with the piecewise quadratic basis functions \(\phi_{i}\) and linear basic function \(\psi_{i}\). Substituting these approximations into the weak form (2) and choosing the test functions from the same discrete spaces results in a large, coupled system of nonlinear algebraic equations. This system can be expressed in the following block matrix form as  
\[ \begin{bmatrix}\mathbf{A}+\mathbf{C}(\mathbf{u}_{h})&\mathbf{B}^{T}\\ \mathbf{B}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\mathbf{u}\\ p\end{bmatrix} = \begin{bmatrix}\mathbf{F}\\ \mathbf{0}\end{bmatrix} , \tag{4}\]  
Among which are  
\[\mathbf{A}_{ij}=\
u\int_{\Omega}\
abla\phi_{i}\cdot\
abla\phi_{j}\,d\Omega,\]  
\[\mathbf{C}(\mathbf{u}_{h})_{ij}=\int_{\Omega}(\mathbf{u}_{h}\cdot\
abla\phi_{ j})\cdot\phi_{i}\,d\Omega,\]  
\[\mathbf{B}_{ij}=-\int_{\Omega}\psi_{i}\cdot\
abla\phi_{j}\,d\Omega.\]  
For the error analysis of this element (4), the velocity field accuracy is \(O(h^{2})\) and the pressure field accuracy is \(O(h)\), with a detailed analysis process found in [PERSON] et al. (2006) [28]. For the temporal discretization, a generalized \(\theta\)-scheme is employed. This allows for a unified framework to control the accuracy and stability of the time integration. In this work, it has been focused on the second-order accurate and Crank-Nicolson scheme for its excellent long-term stability properties. First, consider discretizing the time term using a generalized \(\theta\)-scheme, rewriting the momentum equation in the weak form (2) as  
\[ \begin{cases}\int_{\Omega}\frac{\mathbf{u}_{h}^{n+1}-\mathbf{u}_{h}^{n}}{\Delta t }\cdot\mathbf{v}_{h}\,d\Omega+\theta[\int_{\Omega}(\mathbf{u}_{h}^{n+1}\cdot \
abla)\mathbf{u}_{h}^{n+1}\cdot\mathbf{v}_{h}+\
u\
abla\mathbf{u}_{h}^{n+1} \cdot\
abla\mathbf{v}_{h}-p_{h}^{n+1}\cdot\
abla\mathbf{v}_{h}\,d\Omega]+\\ (1-\theta)[\int_{\Omega}(\mathbf{u}_{h}^{n}\cdot\
abla)\mathbf{u}_{h}^{n} \cdot\mathbf{v}_{h}+\
u\
abla\mathbf{u}_{h}^{n}\cdot\
abla\mathbf{v}_{h}-p_{h }^{n}\cdot\
abla\mathbf{v}_{h}\,d\Omega]=\int_{\Omega}\mathbf{f}^{n+\theta} \cdot\mathbf{v}_{h}\,d\Omega,\\ \int_{\Omega}\mathbf{q}_{h}\cdot\
abla\cdot(\theta\mathbf{u}_{h}^{n+1}+(1- \theta)\mathbf{u}_{h}^{n})=0.\end{cases} \]  
Combining with (4), the final [PERSON] system of equations can be obtained as  
\[ \begin{bmatrix}\frac{1}{\Delta t}\mathbf{M}+\theta(\mathbf{A}+ \mathbf{C}(\mathbf{u}_{h}^{n+1}))&\theta\mathbf{B}^{T}\\ \theta\mathbf{B}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\mathbf{u}^{n+1}\\ p^{n+1}\end{bmatrix} \] \[= \begin{bmatrix}\frac{1}{\Delta t}\mathbf{M}\mathbf{u}^{n}-(1- \theta)(\mathbf{A}+\mathbf{C}(\mathbf{u}_{h}^{n}))\mathbf{u}^{n}-(1-\theta) \mathbf{B}^{T}p^{n}+\mathbf{F}^{n+\theta}\\ -(1-\theta)\mathbf{B}\mathbf{u}^{n}\end{bmatrix} , \tag{5}\]  
Among which we have  
\[\mathbf{M}_{ij}=\int_{\Omega}\phi_{i}\cdot\phi_{j}\,d\Omega.\]  
When \(\theta=\frac{1}{2}\), the temporal accuracy is \(O(\Delta t^{2})\), and otherwise it is \(O(\Delta t)\). The choice of the time-step size (\(\Delta t\)) is crucial for both the accuracy and stability of unsteady simulations.