Exact solutions and superposition rules for Hamiltonian systems generalizing time-dependent SIS epidemic models with stochastic fluctuations

1.   Introduction

The time-dependent IMNSE model is a parabolic partial differential equation, which couples the linear velocity vector and pressure with the angular velocity. Eringen ^[1] incorporated an account of the microstructure fluid into the well-known Navier-Stokes (NS) fluid equations, and obtained the IMNSE model in 1965. The IMNSE model constitutes a framework to describe the dynamics of continuum media ^[2] where the material particles have both translational and rotational degrees of freedom.

Given a final time $T > 0$ and a bounded, regular polygonal/polyhedral domain $\Omega \subset {\mathbb{R}^{{d}}}$ ($d = 2$ or $3$) along with Lipschitz continuous boundary $\partial \Omega$, we consider the numerical approximation of the following IMNSE model ^{[3,4,5,6,7,8,9,10,11,12]}:

where ($u$, $p$, $w$) denote the linear velocity, pressure, and the micro-rotation field, respectively. $({f}, g) = ({f}(x, t), g(x, t))$ are given body force functions. The physical parameters $j$ is inertia density, and $\nu_0$, $v_r$, $c_1$, and $c_2$ are the kinematic viscosities, which are assumed to be constant and positive. The three nondimensional numbers appearing in the equations are the Newtonian kinematic viscosity $\nu_0$, the dynamic micro-rotation viscosity $v_r$, and the angular viscosity $c_1$.

When $d = 3$,

when $d = 2$,

and to simplify notation, we do not differentiate between the 'curl' differences when the spatial dimension is $d = 2$ or $3$.

Because of their significant applications in modern industry^[3], such as bearing lubrication^[4], the IMNSE model has received considerable attention in both theoretical and numerical algorithms in recent years. There are many numerical approaches for solving the IMNSE system, such as that by Ortega-Torres et al. in ^[5], where they presented the fully discrete penalty FE method and obtained optimal error estimates. Jiang et al. analyzed two decoupled time-stepping FE methods in ^[8], where the method decouples velocity and angular velocity. Maimaiti et al.^[9] analyzed the first-order and second-order pressure-correction (PC) projection methods for the MNSE, which decouple velocity and pressure. Zhang et al. in ^[11] combined the scalar auxiliary variable (SAV) approach for the convective terms and some subtle implicit-explicit (IMEX) treatments for the coupling terms, as well as proposed a decoupled, linear, and unconditionally energy stable scheme for MNS system. A decoupling FE method with different time steps for the model was proposed and analyzed in ^[12]. Some projection methods (or fractional-step methods) were proposed and analyzed for the MNSE in ^[13]. From the expressions of the IMNSE system, it is evident that the system encompasses characteristics of incompressibility, strong nonlinearity, and coupling between the angular equation of motion and the fluid equation of motion. This complexity renders solving the MNSE system a challenging endeavor.

In recent years, the grad-div stabilization method for incompressible flow problems has attracted much attention because it can improve the quality of the numerical solution for the velocity vector ${u}$ ^[14]. The usual technique is to add $0 = -\gamma \nabla(\nabla\cdot u)$ to the continuous momentum equation. Because the grad-div term is nonzero in discrete cases, the grad-div term has a positive impact on discrete equations. As the grad-div parameter $\gamma$ increases, it leads to reduced computational efficiency. The matrix resulting from the grad-div term becomes singular for larger stabilized parameter $\gamma$. This will result in solver breakdown. Recently, Qin Y et al. proposed a standard grad-div stabilization (SGD) in ^[15], which involves adding the grad-div stabilization terms $- \beta \nabla \nabla \cdot {u_t}$ and $-\gamma \nabla(\nabla\cdot u)$($\beta \ge 0, \gamma \ge 0$) in the governing equations directly, where $\gamma$ and $\beta$ are stabilized parameters. It is found that the SGD algorithm can also effectively improve the computational accuracy, but, with the increase of parameters, the computational efficiency decreases. Fortunately, Fiordilino, et al. in ^[16] proposed an MGD stabilization method for the NS equations to weaken this problem. The MGD algorithm constitutes two steps, which involve a post-processing step for linear velocity, in which the solution of the equation is not affected by the stabilization parameters. Therefore, the MGD method is widely used in different equations, such as in 2D/3D nonstationary incompressible magnetohydrodynamic equations ^[17], time-dependent thermally coupled MHD equations^[18], and the incompressible non-isothermal fluid flows^[19]. The BDF2-MGD stabilization method was proposed for the unsteady NS equations ^[20] and the Stokes/Darcy equations ^[21]. A series of new expansions for the MGD stabilization method have been presented, including the Crank-Nicolson LeapFrog MGD stabilization method for the time-dependent Navier-Stokes equations in ^[22], and the rotational pressure-correction method based on the MGD in ^[23]. In this paper, we devised a BDF2 modular grad-div stabilization method for the micropolar Navier-Stokes equations. This method not only avoids solver breakdown, but also improves the efficiency of the solution and preserves the advantages of the grad-div stabilization method.

Suppressing the spacial discretization for the moment, we propose a two-step time semi-discrete method for the IMNSE model. Let $\Delta t$ denote the size of the time steps, and $t_n = n\Delta t$, where $n = 0, 1, 2, ..., N$ with $T: = N\Delta t$.

Algorithm 1.1. Modular grad-div stabilization

Step 1. Given $u^{-1} = u^{0};w^{-1} = w^0$, find ${\hat u^{n + 1}}, w^{n+1}$ and $p^{n+1}$ for $n = 0, 1, 2, \cdots, N-1$ satisfying:

Step 2. Given ${\hat u^{n + 1}}$, find $u^{n+1}$ satisfying:

In the provided equations, $\beta \ge 0$ and $\gamma \ge 0$, and their specific values depend on the application at hand. The combined result of Steps 1 and 2 is the precise discretization of the model by using the consistent BDF2 method.

Algorithm 1.1 is essentially equivalent to Eq (1.6). If Eq (1.6) is discretized by using a decoupled linearly extrapolated BDF2 scheme, we refer to it as standard grad-div (SGD) stabilization. If the stabilized parameters $\beta = \gamma = 0$, the (1.6) reduces to the BDF2 scheme without stabilization.

Our main contributions to the work are four-fold:

(1) By combining MGD with a linearly extrapolated BDF2 scheme, we have proposed a decoupled, linear, and second-order scheme for the IMNSE model. Step 1 of this scheme only requires solving two smaller sub-physics problems at each time step, which reduces the size of the linear systems to be solved and allows for parallel computing of the two sub-physics problems.

(2) For the SGD stabilization method (1.6), the module grad-div stable terms are directly added to the momentum equation (the first equation of (1.6)). The grad-div stable terms of the MGD method are separated as the second step (post-processing) in this paper, making it easy to implement. Not only can the MGD stabilization method keep the advantages of the grad-div algorithm, but it also avoids the influence of large parameters $\beta$ and $\gamma$ on the solution.

(3) We have conducted rigorous unconditional stability and error estimates for the proposed MGD linearly extrapolated BDF2 scheme.

(4) We have provided several numerical examples in 2D and 3D settings to verify the theoretical findings and to demonstrate the benefits of the MGD stabilization method.

The organization of this article is outlined below. Section 2 furnishes an introduction to fundamental mathematical concepts. Section 3 offers an detailed outline of the MGD algorithm for the IMNSE model, along with an analysis of its stability. Section 4 is devoted to convergence analysis of the MGD algorithm. A series of numerical tests are used to validate the theoretical findings in Section 5. Finally, in Section 6, the paper is concluded.

2.   Preliminaries

This work employs the conventional and well-established Hilbertian and Sobolev spaces, equipped with their respective inner product and norm definitions.

For a given function $a(x, t)$ defined on $[0, T]$, we establish the norm as enumerated below:

Additionally, we introduce the following discrete norms:

The functional spaces are outlined below:

The dual space of $H_0^1(\Omega)$ is defined as ${H^{-1}}(\Omega)$ and endowed with the negative norm.

The expression governing the curl operator ^[24] is

The skew-symmetrized trilinear form is indicated as

where $b^*(u, v, w)$ satisfies the bounds ^[24]

Furthermore, we introduce the finite element spaces. Let ${\tau _h} = \cup{K}$ be a uniform triangulation/tetrahedron partition of domain $\Omega^{d}(d = 2$ or $3)$, with $h: = \mathop {\max }\limits_{\forall \:K \in {\tau _h}} \; \; diameter(K)$. Further, we introduce the finite element subspaces $\left({{X_h}, {Q_h}, {Y_h}} \right) \subset \left({X, Q, Y} \right)$ for linear velocity, pressure, and angular velocity, respectively. We assume that $\left({{X_h}, {Q_h}} \right)$ satisfies the discrete $LBB^h$ ^[24] condition, meaning there exists a constant $M_0$ such that

There are finite elements that satisfy the stability of $LBB^h$, such as the Tayler-Hood finite elements $(P_k, P_{k-1})$ and Mini finite elements $(P^{b}_1, P_1)$ on the quai-uniform triangular(2D)/tetrahedral(3D) meshes, where $P_k$ is defined by $k$-order polynomials on element $K$.

The discrete divergence-free space $V_h$ is defined as:

3.   Numerical scheme and its stability

We first present the fully discrete linearly extrapolated BDF2 MGD scheme, then conduct stability analysis.

3.1.   Numerical scheme

Algorithm 3.1 (MGD-BDF2 full discrete decoupled scheme).

Step 1. Given $u_h^{ - 1} = u_h^{0} = u^{0}\in V_{h}$ and $w_h^{ - 1} = w_h^{0} = w^{0}\in Y_{h}$, for all $({v_h}, {q_h}, {z_h}) \in ({X_h}, {Q_h}, {Y_h})$ find $(\hat u_h^{n + 1}, p_h^{n + 1}, w_h^{n + 1}) \in ({X_h}, {Q_{h, }}{Y_h})$, satisfying:

Step 2. Given $\hat u_h^{n + 1} \in {X_h}$, for all $v_h \in X_h$ find $u_h^{n + 1}\in{X_{h}}$, satisfying:

Remark 3.2. In Step 1, we decouple the fully coupled IMNSE model into two smaller sub-physics problems at each time step (one is for the linear velocity and pressure, the other is for the angular velocity), which reduces the size of the linear systems to be solved and allows for parallel computing of the two sub-physics problems.

Remark 3.3. Since Algorithm 3.1 is a two-step scheme, it requires that the initial values $(u^{0}_{h}, w^{0}_{h})$ and $(u^{1}_{h}, w^{1}_{h})$ are second-order accuracy. For simplicity, here we use ghost points $(u^{-1}_{h}, w^{-1}_{h})$ and set $(u^{-1}_{h}, w^{-1}_{h})$ = $(u^{0}_{h}, w^{0}_{h})$ = $(u^{0}, w^{0})$. This ensures that $(u^{0}_{h}, w^{0}_{h})$ and $(u^{1}_{h}, w^{1}_{h})$ are of second-order accuracy.

3.2.   Stability analysis

In this subsection, we will conduct a stability analysis.

Lemma 3.4. ^[20] For the MGD-BDF2 scheme, the following identities hold for Eq (3.4):

and

The subsequent statement certifies its unconditional stability.

Theorem 3.5. Assume that $f, g \in {L^2}(0, T;{H^{ - 1}}{(\Omega)^d})$, then the following holds for $\Delta t > 0$:

Proof. Substituting ${v_h} = \hat u_h^{n + 1}, {q_h} = p_h^{n + 1}$ into (3.1) and (3.2), and ${z_h} = w_h^{n + 1}$ into (3.3), respectively, we derive the following equations:

We first add Eqs (3.8) and (3.9), and apply the identity $2(3a-4b+c)a = a^2-b^2+(2a-b)^2-(2b-c)^2+(a-2b+c)^2$ and Lemma 3.4. Then we sum the equation from $n = 0$ to $N-1$ to obtain

By employing Young's inequality and the Cauchy-Schwartz inequality, we establish the following relation for the terms on the right-hand sides (RHS) of (3.10):

Substituting (3.11) and (3.12) into (3.10), we end the proof and obtain Theorem 3.5.

4.   Error analysis

In what follows, an account of the error estimates of the MGD-BDF2 algorithm shall be provided.

Assumption 4.1 To obtain the most precise error estimates for the MGD-BDF2 Algorithm 3.1, we need the regularity

Additionally, we define two types of projections^[25,26]. Denote by $(s_{h}, l_{h})$ the Stokes projection of $(u, p)\in X\times Q$ such that

and the Ritz projection operator ${P_h}:w(t) \in Y \to {P_h}w \in {Y_h}, \forall t \in \left[ {0, T} \right],$

Furthermore, for $0\le m\le k$, the following approximation properties hold:

Following ^[20], the following inequality holds:

As for the Ritz projection (4.2), it satisfies the approximation property ^[27]

To facilitate the convergence analysis, we introduce error decompositions. Let ${{\tilde u}_h}$ be an approximation of $u$ in $V_h$. We denote ${u^n} = u({t^n}), {{\tilde u}^n} = \tilde u({t^n})$, ${w^n} = w({t^n}), {{\tilde w}^n} = {P_h}\tilde w({t^n})$. The following error decompositions are

Lemma 4.2. Under regularity Assumption 4.1, we have the following consistency error relation: $\forall\; \sigma > \; 0$,

where

Proof. Utilizing the integral version of Taylor's theorem and the Cauchy-Schwarz inequality, we can derive the following expressions:

The proof of the formulas $\left({\frac{{3{u^{n + 1}} - 4{u^n} + {u^{n - 1}}}}{{2\Delta t}} - {u_t}({t_{n + 1}})} \right)$ and $\left({\nabla ({u^{n + 1}} - 2{u^n} + {u^{n - 1}})} \right)$ can be found in the references ^[28,29]. Similarly, we can establish an expression for ${\tau ^{n + 1}}({z_h})$.

To establish convergence, a critical lemma concerning Step 2 of Algorithm 3.1 is required.

Lemma 4.3. ^[20] The following inequality holds:

Theorem 4.4. Under the condition of Assumption 4.1 and $K\Delta t < \theta\in (0, 1)$, where $K = {\max _{0 \le t \le T}}\left\{ {1 + C{\gamma ^{ - 1}}\left\| {u(t)} \right\|_2^2, \left({C\left\| {\nabla u(t)} \right\|_{{L^\infty }}^2 + \frac{{v_r^2}}{{{c_1}}}} \right), {\nu _0}} \right\}$, the error between the solutions of Algorithm 3.1 and the exact solution of (1.1) satisfies

Proof. Subtracting (3.1)–(3.4) from (1.1)–(1.5), respectively, we get the error equations

Taking ${{\tilde u}_h} = {s_h}, {v_h} = \Lambda _{u, h}^{n + 1} \in {V_h}$ in (4.9), we obtain the equation

Setting ${v_h} = \frac{{3\phi _{u, h}^{n + 1} - 4\phi _{u, h}^n + \phi _{u, h}^{n - 1}}}{3} \in {V_h}$ in (4.12), we have the following equations;

Integrating (4.13) and (4.14) into one equation, using the error decompositions (4.4), multiplying both sides of the integrated equation by $4\Delta t$, and applying Lemma 4.3, we obtain

Applying the Cauchy-Schwarz and Young's inequalities on the RHS of (4.15), we obtain

and

Using Lemma 4.2, we obtain

For the trilinear terms, we handle them as follows:

Subsequently, we proceed to perform bound estimation on the nonlinear terms after the above processing:

and for the last term as in Eq (4.21), we use a similar treatment in ^[20] to obtain Eq (4.24),

Combining inequalities (4.16)–(4.25) into (4.15), we get

Similarly, taking ${z_h} = \phi _{w, h}^{n + 1}$ in (4.11), the angular velocity equation is similarly processed as (4.26). A series of inequalities related to angular velocity can be derived by making use of the Cauchy-Schwarz-Young inequality, and integrating all inequalities yields

The equations for linear velocity (4.26) and angular velocity (4.27) are combined and generalized into a single equation. The sum of $n$ from $0$ to $N-1$ is given by

When $\Delta t$ is small enough, such that $K\Delta t \le 1$, we use the discrete Gronwall lemma ^[16] to get

where ${K^{n+1}} = \max \left\{ {1 + C{\gamma ^{ - 1}}\left\| {{u^{n + 1}}} \right\|_2^2, \left({C\left\| {\nabla {u^{n + 1}}} \right\|_{{L^\infty }}^2 + \frac{{v_r^2}}{{{c_1}}}} \right), {\nu _0}} \right\}.$ Finally, we conclude the proof by applying the triangle inequality.

Corollary 4.5. Under the conditions of Theorem 4.4 with $k = 2$ and $m = 1$, we have

Remark 4.6. Similarly, when using the $(P^{b}_1, P_1, P_1)$ elements to approximate $(u, p, w)$, the error is $C(h+(\Delta t)^2)$.

5.   Numerical tests

We present four illustrative examples in this section. The initial pair of examples serves to compute error and convergence rates in three and two dimensional cases, respectively. The subsequent two examples involve the simulation of square cavity-driven flow and the bearing lubrication problem. For the implementation of the algorithm, we have opted to utilize the publicly available finite element software FreeFem++^[30].

5.1.   Convergence order validation

5.1.1.   Analytical solution in three dimensions

The fluid domain is $\Omega = {[0, 1]^3}$. The parameters of the MGD-BDF2 method are $\gamma = 1.0$, $\beta = 0.2$, ${\nu_0} = 2.0$, $\nu_r = 1.0$, $j = 1.0$, $c_1 = 2.0$, and $c_2 = 1.0$. We take the following analytical solutions:

Analytical solution 1:

Analytical solution 2:

The convergence rates are determined by analyzing the errors at two consecutive mesh sizes. Both the spatial convergence rates (SCR) and temporal convergence rates (TCR) will be evaluated in this study. The errors and convergence rates for $u, p$ and $w$ in the $L_2$ norm and $H_1$ norm are recorded by using the $(P^{b}_1, P_1, P_1)$ and $(P_2, P_1, P_2)$ elements, respectively.

The data results of SCR are placed in Tables 1 and 2 by using analytical solution 1 (5.1). The errors and convergence rates of TCR are placed in Tables 3 and 4 using analytical solution 2 (5.2).

From Tables 1–4, it is easy to see that these findings satisfy the convergence rates of Theorem 4.4 (i.e., Corollary 4.5 and Remark 4.6).

5.1.2.   Analytical solution in two dimensions

We will check the convergence performance of the IMNSE model in two dimensions by using analytical solution (5.3). We set $\Omega = {(0, 1)^2}$, $\nu_0 = 2.0$, $c_1 = 1.0$, $c_2 = 1.0$, $\beta = 0.2$, $\gamma = 1.0$, $j = 1.0$, $\alpha = 1.0$.

Tables 5 and 6 are the numerical results of spatial convergence rates (SCR), while Table 7 is the numerical results of the temporal convergence rates (TCR). We see that these findings of Tables 5–7 satisfy the convergence rates of Theorem 4.4 (i.e., Corollary 4.5.).

5.2.   Computational efficiency

In order to demonstrate the superiority of the MGD algorithm, we conducted two numerical experiments. The first numerical experiment is to calculate the error of $\nabla\cdot u$ by using the analytical solution (5.3) at $T = 0.01$, $\Delta t = 0.0001$, for the BDF2 scheme (without grad-div terms), the standard grad-div (SGD-BDF2) scheme, and the modular grad-div (MGD-BDF2) scheme. Numerical results are presented in Table 8. It is evident that the SGD-BDF2 and MGD-BDF2 schemes are superior to the BDF2 scheme without grad-div terms.

Next, our experiment focuses on increasing $Re$, where $Re \approx \nu _0^{ - 1}$ is the Reynolds number. We set $\Delta t$ = $h$ = $\frac{1}{32}$, $\beta$ = $0.2$, and $\gamma = 10$. In Table 9, we compare the errors for divergence velocity with increasing $Re$ of BDF2, SGD-BDF2, and MGD-BDF2. We notice that the error of the proposed algorithm barely increases. It can be seen that the addition of stabilization can enhance the conservation of mass sufficiently.

The third numerical experiment is to test the computational efficiency by comparing the computational time of the SGD and MGD stabilization methods by using analytical solution (5.3) in 2D, and (5.2) for 3D. Numerical results are given in Tables 10 and 11.

In the 2D test, we fixed $\Delta t = h = \frac{1}{32}$ while varying the parameters $\gamma$ and $\beta$ and conducted a comparison of the elapsed times for the SGD-BDF2 and MGD-BDF2 schemes. For the 3D test, $\Delta t = \frac{1}{80}$ and $h = \frac{1}{10}$ were selected. Notably, when both $\gamma$ and $\beta$ are set to $0$, the MGD-BDF2 scheme essentially transforms into the classical Galerkin BDF2 scheme. We used the standard GMRES solver to calculate the SGD-BDF2 scheme and Step 1 of the MGD-BDF2 scheme, and a CG solver to calculate Step 2 of the MGD-BDF2 scheme. When the GMRES solver fails to converge after a single iteration for certain $\gamma$ and $\beta$, we denote this as "Failed" in tables. From Table 10 for 2D and Table 11 for 3D, we can observe that the MGD-BDF2 scheme is superior to the SGD-BDF2 scheme.

Regarding the SGD-BDF2 algorithm, as the stabilization parameter increases, the overall CPU time generally exhibits a growing trend. Notably, in the cases of $\beta = 0$, $\gamma = 40$ in 2D and 3D, $\beta = 0.4$, $\gamma = 0.3$ in 2D, and $\beta = 4$, $\gamma = 0.3$ in 3D, the SGD-BDF2 algorithm failed to converge successfully, and is therefore marked as "Failed".

In contrast, the CPU time performance of the MGD-BDF2 algorithm remains relatively stable, without showing a significant increase as the stabilization parameter increases. Furthermore, despite significant variations in the stabilization parameter, the MGD-BDF2 algorithm demonstrates remarkable adaptability, enabling it to consistently maintain a relatively stable solving time.

5.3.   The Lid-Driven flow

Next we simulated the flow on $\Omega = {[0, 1]^2}$ by applying similar no-slip boundary conditions at the lid. The upper boundary, where $y = 1$ and $0 < x < 1$, is subject to the conditions: $u_1 = 1$, $u_2 = 0$, and $w = 0$. The normal component of velocity is assumed to be zero on $\partial \Omega$, while the tangential component is zero, except at $y = 1$ where it is set to 1. Figures 1–3 display the components of velocity ($u_1$, $u_2$) and angular velocity ($w$) contours for the MGD-BDF2 scheme and different viscosity coefficients $\nu_0$ by using $\frac{1}{h} = 48$ and $(P^{b}_1, P_1, P_1)$ elements.

Upon observation of Figures 1–3, it becomes evident that as the viscosity coefficient $\nu_0$ diminishes, the components of velocity ($u_1$, $u_2$) and angular velocity ($w$) gradually lose their symmetry and deviate toward one side.

5.4.   The Bearing lubrication problem

The IMNSE model is used to simulate lubrication problems ^[4]. The relevant fluid domain encompasses a ring-shaped region bounded by the outer boundary ${\Gamma _1}$ with a radius of ${r_1}$, and the inner boundary ${\Gamma_2}$ with a radius of ${r_2}$. In Figures 4–6, $(u_1, u_2)$ and $w$ adhere to the homogeneous dirichlet boundary condition on $\Gamma_1$; while on ${\Gamma_2}$, $w = 0$ and the following conditions hold: ${u_1} = - {r_2}{w_r}\sin \theta, {u_2} = - {r_2}{w_r}\cos \theta$, $\theta = arctan\frac{y}{x}$. Here, ${w_r}$ represents the rotational angular velocity. We will explore three distinct scenarios, where ${w_r}$ takes on the values of $200,600,$ and $1000$. The parameters are selected as ${r_1} = 0.1$, ${r_2} = 0.06$, $\nu_0 = 2.0$, $c_1 = c_2 = j = 1.0$, $\beta = 0.2$, $\gamma = 1.0$, $h = \frac{1}{150}$, and $\Delta t = 0.001$.

Figures 4–6 show that the velocity components ($u_1$, $u_2$) and angular velocity ($w$) increase as the rotational angular velocity $w_r$ increases. Consequently, the bearing can endure higher loads. However, upon comparing the velocity components and angular velocity, it is evident that the polarity effect of the fluid remains minimal, with the angular velocity tending to rotate around the inner circle. This implies that, as the rotational angular velocity $w_r$ rises, the micro-polarity effect of the fluid becomes more pronounced.

6.   Conclusions

We have introduced a novel method (linearly extrapolated fully discrete MGD-BDF2) for solving the IMNSE model in $\mathbb{R}^{d}, d = 2, 3$. Our scheme is second-order, fully discrete, and semi-implicit, ensuring unconditional stability and optimal convergence rates in both time and space. By employing this approach, we only need to solve linear systems at each time step, which simplifies the implementation process and improves efficiency. Additionally, the study compared the computational time differences between the MGD method and SGD methods with varying stabilization parameters. Standard implementations showed a rapid increase in costs as the parameters increased. In addition, the paper also simulated the flow of cavity fluids and bearing lubrication issues. These simulations all yielded good results. The results show that the MGD-BDF2 method has a significant advantage in computational time and keep mass conservation.

The MGD method decouples the velocity and angular velocity, yet the velocity and pressure remain coupled, which is typically a key challenge in solving the NS equations. In the future, we will continue to work on the expansion of this model to a wider range of applications. For example, in combination with pressure correction and the scalar auxiliary variable (SAV), we will explore whether it can maintain the advantages of several algorithms.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (12171141).

Conflict of interest

The authors declare no conflict of interest.