## 3 Results
### Flow over a Square  
The first test case considered herein was the flow over a square with periodic boundary conditions in both horizontal and vertical directions. The geometry of the problem is shown in Figure 3, which consists of a single square cylinder of side length \(D=0.04\) m and its associated volume within the square lattice of side length \(L=0.1\) m.  
The flow is driven by a body force \(\mathbf{F}\) along the horizontal direction. On the wall of the square cylinder, the no-slip boundary condition is imposed. The velocity field is initialised with zero values everywhere. The reference pressure was set to zero at the left-bottom corner. The kinematic viscosity was set to \(v=10^{-6}\) m\({}^{2}\) s\({}^{-1}\), and the body force was set to \(F=1.5\times 10^{-5}\) m s\({}^{-2}\). The maximum horizontal velocity \(U_{max}\) can reach about \  
comparison, the relative difference for a flow variable between the two numerical solutions was computed as the root-mean-squared difference normalised by the root-mean-squared value of the variable from the reference solution. For the steady-state results, the relative differences for \(U\), \(V\), and \(P\) were \(1.2\times 10^{-4}\), \(2.8\times 10^{-3}\), and \(5.4\times 10^{-3}\), respectively. Figure 5 shows the time histories of \(U\) and \(V\) at Monitoring Point 1. Figure 6 shows the time histories of \(P\) at the two monitoring points 1 and 2. It turned out that the agreement between the numerical solutions from the two methods was also excellent with regard to the time evolution of these flow variables. The relative differences for \(U_{1}\), \(V_{1}\), \(P_{1}\), and \(P_{2}\) were \(1.5\times 10^{-4}\), \(5.9\times 10^{-4}\), \(3.3\times 10^{-3}\), and \(5.4\times 10^{-3}\), respectively. To further verify the correctness of the two methods above, the flow problem here was also simulated by using the commercial CFD software ANSYS Fluent. The spatial discretisation was performed by the standard finite volume method, and the time integration was performed by a second-order implicit scheme. The results obtained by ANSYS Fluent are not shown because they visually overlapped with those obtained by the hybrid and  
Figure 4: Comparison of numerical solutions obtained by the reference method (**left**) and the hybrid method (**right**) for the first test case at \(t=4000\) s. Shown from top to bottom are the contours of the horizontal velocity component \(U\), the vertical velocity component \(V\), and the pressure \(P\), respectively. reference methods. The relative differences between the hybrid and Fluent solutions of \(U\), \(V\), and \(P\) were \(1.6\times 10^{-3}\), \(1.1\times 10^{-2}\), and \(5.3\times 10^{-3}\), respectively. The relative differences between the reference and Fluent solutions of \(U\), \(V\), and \(P\) were \(1.5\times 10^{-3}\), \(9.7\times 10^{-3}\), and \(4.5\times 10^{-3}\), respectively. To assess the spatial accuracy of the hybrid and reference methods, simulations with larger grid sizes of \(L/200\), \(L/100\), and \(L/50\) were performed. Using the steady-state solutions obtained with the smallest grid size \(L/400\) (i.e., those shown in Figure 4) as references, the root-mean-squared errors for the solutions of coarser grids were computed. Figure 7 shows the log-log plot of the errors as a function of the grid size for the horizontal velocity component \(U\), the vertical velocity component \(V\), and the pressure \(P\). As indicated by the fit slope, the convergence rates for those flow variables were around 1.5. Therefore, for both the hybrid and reference methods, the spatial order of accuracy was lower than second order, although the spatial discretisation was based on the second-order Taylor's expansion. It is reasonable to infer that the loss of accuracy was due to the WLS approximation. To assess the temporal accuracy of the hybrid and reference methods, simulations with the grid size fixed to \(L/200\) and the time step varying from 0.00625 s to 0.05 s were performed.