# An Improved Galerkin Framework for Solving Unsteady High-Reynolds Navier-Stokes Equations
## 3 Numerical Methods
### Variational Multiscale Stabilization  
The fundamental idea behind the Variational Multiscale Stabilization (VMS) method is to conceptually partition the solution space into resolved (large-scale) and unresolved (small-scale) components. In the context of FEM, the large scales correspond to the discrete information that can be captured by the computational mesh, while the small scales represent the subgrid phenomena that cannot be resolved. The VMS framework provides a systematic way to model the effect of these unresolved small scales back onto the resolved large scales. Unlike traditional stabilization methods that add artificial diffusion somewhat arbitrarily, VMS derives the stabilization term directly from the underlying equations, leading to a more physically consistent and less intrusive form of stabilization. Following the VMS derivation, first, the solution has been decomposed into large scale \((h)\) and small scale \((^{\prime})\)[15; 16] as  
\[\mathbf{u}=\mathbf{u}_{h}+\mathbf{u}^{\prime},\quad p=p_{h}+p^{\prime}, \tag{6}\]  
where \(\mathbf{u}_{h}\in V_{h}\subset H^{1}(\Omega)\), \(p_{h}\in Q_{h}\subset L^{2}(\Omega)\) are the large-scale solutions, which are the numerical solutions. And \(\mathbf{u}^{\prime}\in V_{h}^{\perp}\) and \(p^{\prime}\in Q_{h}^{\perp}\) are small-scale solutions, regarded as unresolvable parts, and approximated by subgrid-scale models. Small scales are updated according to the following method  
\[\mathbf{u}^{\prime n+\theta}=\tau_{1,t}\frac{1}{\theta\Delta t}\mathbf{u}^{ \prime n}-\tau_{1,t}\Pi^{\perp}\Big{(}a\cdot\
abla\mathbf{u}_{h}^{n+\theta}+ \
abla p_{h}^{n+1}\Big{)},\quad p^{\prime n+\theta}=-\tau_{2}\Pi^{\perp}( \
abla\cdot\mathbf{u}_{h}^{n+\theta}),\]  
where  
\[\tau_{1}=\left[\left(c_{1}\frac{v}{h^{2}}\right)^{2}+\left(c_{2}\frac{|a|}{h} \right)^{2}\right]^{-1/2},\quad\tau_{2}=\frac{h^{2}}{c_{1}\tau_{1}},\quad\tau _{1,t}=\left(\frac{1}{\theta\Delta t}+\frac{1}{\tau_{1}}\right)^{-1}\]  
are stabilization parameters. \(\Pi^{\perp}=I-\Pi\) is the orthogonal projection operator and \(\Pi\) is the standard \(L^{2}\) projection. Their formulations are based on the work of [PERSON] (2002) [29], and it has been incorporated by local, cell-wise information of the velocity field magnitude \(|a|\), which makes the stabilization more sensitive to the local flow dynamics [24]. Furthermore,  
\[a=\mathbf{u}_{h}^{n+\theta}+\mathbf{u}^{\prime n+\theta},\]  
\[\mathbf{u}^{\prime 0}=\mathbf{0},\quad p^{\prime 0}=0,\]  
where \(p^{n+\theta}=\theta p^{n+1}+(1-\theta)p^{n},\mathbf{u}_{h}^{n+\theta}=\theta \mathbf{u}_{h}^{n+1}+(1-\theta)\mathbf{u}_{h}^{n}\). Substituting the decomposed solutions into the weak form (2), and separating the contributions of large and small scales (taking large-scale test functions), one has  
\[\left\{ \begin{aligned} &\int_{\Omega}\frac{\partial(\mathbf{u}_{h}+ \mathbf{u}^{\prime})}{\partial t}\cdot\mathbf{v}_{h}\,d\Omega+\int_{\Omega}(( \mathbf{u}_{h}+\mathbf{u}^{\prime})\cdot\
abla)(\mathbf{u}_{h}+\mathbf{u}^{ \prime})\cdot\mathbf{v}_{h}\,d\Omega+\
u\int_{\Omega}\
abla(\mathbf{u}_{h}+ \mathbf{u}^{\prime})\cdot\\ &\
abla\mathbf{v}_{h}\,d\Omega-\int_{\Omega}(p_{h}+p^{\prime}) \cdot\
abla\mathbf{v}_{h}\,d\Omega=\int_{\Omega}\mathbf{f}\cdot\mathbf{v}_{h} \,d\Omega,\\ &\int_{\Omega}\mathbf{q}_{h}\cdot(\
abla\cdot(\mathbf{u}_{h}+ \mathbf{u}^{\prime}))\,d\Omega=0.\end{aligned} \right.