## 2 Materials and Methods  
### Governing Equations and Numerical Setup  
The analysis of the vortices generated by MVGs was completed by first solving the steady-state 3D Navier-Stokes equations for incompressible Newtonian fluids, which are based on the conservation of mass and momentum. These equations are represented as  
\[\
abla\cdot\mathbf{u} =0\,, \tag{1}\] \[(\mathbf{u}\cdot\
abla)\mathbf{u} =-\frac{1}{\rho}\
abla p+\
u\
abla^{2}\mathbf{u},\]  
where \(\mathbf{u}\) [m/s] denotes the velocity vector, \(p\) [Pa] represents pressure, \(\rho\) [kg/m\({}^{3}\)] is density, and \(\
u\) [m\({}^{2}\)/s] stands for kinematic viscosity. Ansys Fluent 2020 R1 was used to solve these equations. The Reynolds number based on boundary layer thickness (\(Re_{\delta}\)) is defined as  
\[Re_{\delta}=\frac{u\delta}{\
u}, \tag{2}\]  
where \(u\) is the velocity in [m/s], \(\delta\) is the boundary layer thickness in [m], and \(\
u\) [m\({}^{2}\)/s] is the kinematic viscosity. This Reynolds number is around 300 at the end of the domain. The Reynolds number based on distance from leading edge \(Re_{x}=\frac{ux}{\
u}\) is around \(1.5\times 10^{5}\)[22].  
The computational meshes used for the various configurations are consistent with those reported by [PERSON] et al. [10], maintaining a uniform height of 13 mm in all cases. Research by [PERSON] et al. [23] indicated that, provided the computational domain is adequately large, altering its size does not significantly impact the results. Figure 1 shows a schematic of the computational domain. The domain starts at \(X_{0}\), with \(X_{MVG}\) indicating the center of the MVG relative to the domain's beginning. To optimize computational efficiency, simulations were conducted for only half of the domain, corresponding to a single MVG from each pair. A symmetry condition, as described by [PERSON] et al. [23], was applied in both sides of the domain for these simulations. A fully developed Blasius profile is implemented at the inlet in the wall-normal coordinate (y) while ensuring uniformity in the spanwise direction (z). The boundary layer thickness at the MVG location is \(\delta_{99,MVG}=3.610\,\mathrm{mm}\). At the outlet, a pressure boundary condition is set with a gauge pressure of 0 Pa, which was also applied at the top boundary, with no velocity gradient perpendicular to the boundary. To maintain consistency across different cases, careful attention was paid to avoid significant variations in mesh size or quality. The mesh is finest near the MVG and becomes progressively coarser further away. Typically, the mesh consists of approximately 3 to 5 million elements. It was verified that changes in mesh resolution have a negligible effect on the results, with a discretization error around 0.5%. Mesh independence was analogous to the study of [PERSON] et al. (2019), where the Grid Convergence Index method of [PERSON] (2018) is employed. Grids of three different sizes were tested, with the resulting investigated parameters (such as the location of the vortex center or circulation value) all falling between the 95% confidence intervals of the fitting parameters discussed in Section 2.3. The QUICK (Quadratic Upstream Interpolation for Convective Kinematics) scheme is used for discretizing momentum, while the PRESTOI scheme is employed for pressure discretization. The pressure-velocity coupling uses the segregated SIMPLEC (Semi-Implicit Method for Pressure Linked Equations-Consistent) method. Each simulation continues until the residuals achieve double-precision machine accuracy before concluding. For analyzing the downstream drag reduction abilities, a BiGlobal stability analysis is performed, detailed in [PERSON] et al. (2019).