Weak Galerkin method for the Navier-Stokes equation with nonlinear damping term

1.   Introduction

Humanity has always been attracted to celestial bodies. However, our ability to observe the universe is significant biased. The light released by the atmosphere; as by a star is refracted as it enters the human eyes or optical telescopes, it results in the formation of the our perceived star image. Solar radiation affects the erratic motion of molecules in the atmosphere, resulting in unpredictable fluctuations in physical characteristics such as temperature, pressure, and humidity. This unpredictable thermal motion causes a disruption in the refraction of light, resulting in varying light routes and random fluctuations in the amplitude and phase of the light reaching the surface. As a consequence, stars appear to sparkle and flicker. Atmospheric turbulence, the described atmospheric motion phenomena, disrupts the refraction of light. Due to atmospheric turbulence, imaging findings appear blurred and deteriorated even if weather conditions are optimal and the technology is enhanced, not meeting the expected theoretical outcomes. Consequently, the research on turbulence has consistently been a prominent topic in the realm of fluid mechanics.

Atmospheric turbulence refers to the occurrence of turbulence in gaseous fluids; however, turbulence can also be observed in other types of fluids. Once the Reynolds number exceeds a critical threshold, the fluid streamlines become disorganized, and numerous small eddies emerge within the flow field. As a consequence, there is sliding and blending between neighboring layers of flow, which induces the fluid to generate a velocity component that is perpendicular to the direction of the streamlined flow, ultimately leading to an erratic motion. The Navier-Stokes (NS) equation is widely acknowledged as the prevailing model for describing turbulence processes. We consider the steady incompressibility Navier-Stokes equation with a nonlinear damping term. Find $(u,p)$ such that

where $\Omega\in\mathbb{R}^d\; (d = 2,\; 3)$ is an open bounded domain with Lipschitz continuous boundary. ${\boldsymbol{u}}=(u_1,...,u_d)^T$ represents the fluid velocity, $p$ represents the pressure, and ${\boldsymbol{f}}$ represents the external force. $\mu>0$ is the viscosity coefficient. $2< r<\infty$ and $\alpha>0$ are two damping parameters. Nonlinear damping describes the damping force of a flowing fluid that is not directly proportional to its velocity. This phenomenon is frequently observed in porous medium flow, friction effects, and dissipation mechanisms ^[1,9].

The Navier-Stokes equation with damping is a complex nonlinear second-order partial differential equation with multiple variables and unknown parameters, making it challenging to determine its solution. At present, the mainstream method used to tackle this problem is the mixed finite element method (MFEM) ^[6]. The Navier-Stokes-type variational inequality with nonlinear damping term and friction boundary conditions has been discussed in ^[16,17]. Moreover, Li et al. developed a two-level mixed finite element approach and demonstrated its ability to significantly reduce computational time without compromising accuracy^[11]. A stabilized mixed finite element method was proposed for the NS equations with damping in ^[12]. The existence and uniqueness of the weak solutions were proven by the Brouwer fixed-point theorem, and numerical examples were implemented to confirm the theoretical analysis. Furthermore, ^[15] studied the multi-level stabilized algorithms, which combine the stabilized finite element technique with the multi-level method. Using a regularized mortar-type finite element discretization, ^[5] proposed a partitioned Dirichlet-Neumann algorithm for the fluid fow, which is described by the incompressible NS equations. ^[4] adopted an incremental pressure fractional step method and proposed a hybrid vertex-centered finite volume/finite element method for the NS equations on unstructured grids.

In the context of these Eqs (1.1)–(1.3), namely the Stokes equation with a damping term, there have been more studies on its numerical solution. Its existence and uniqueness have been proved, and the conforming MFEM has been developed to discretize the model. Additionally, superclose and superconvergence results for the MFEM applied to the Stokes equations with damping were obtained in ^[18]. In ^[27] the interior penalty discontinuous Galerkin (IPDG) schemes on general mesh were established, and the existence, boundedness, and uniqueness of the discrete solutions were analyzed. With the help of a local pressure projection stabilization, low-order finite element pairs were applied to approximate the velocity and pressure in ^[10]. Furthermore, Burman and Hansbo presented a continuous interior penalty finite element method with edge stabilization for the incompressible Navier-Stokes equations and proved energy-type a priori error estimates independent of the local Reynolds number^[2].

In this paper, we are concerned with the weak Galerkin (WG) finite element method. The WG finite element method is a novel numerical technique that has emerged in recent years for solving partial differential equations ^[3]. The main idea behind this method is the employment of totally discontinuous weak functions on the partitions and the replacement of classical derivatives by weak derivatives. In ^[20], the WG scheme was extended to general polygon mesh by introducing a stabilizer. When comparing with mixed finite element methods, the WG finite element methodis more flexible due to the availability of low-order weak Galerkin finite element pairs. In the past few years, this method has been utilized to address the Stokes equation ^[13,21], Brinkman equation ^[23,25], Helmholtz equation ^[22], and linear elastic equations ^[19,24]. In contrast to the references that focused on the WG scheme applied to the Navier-Stokes equations^[8,26], we introduce a nonlinear damping term to the NS equation, making it more complex. Dealing with this damping term presents a challenge due to the need for more mathematical techniques to prove error estimates, including the topological degree lemma, which is detailed in the following section.

We propose the weak Galerkin finite element method for solving the Navier-Stokes equation with nonlinear damping. The paper is structured as follows. Section 2 outlines the variational formulation of the problem and introduces the notion of weak finite element spaces as well as weak differential operators such as weak gradient and weak divergence. In addition, we establish the WG finite element numerical scheme and prove the existence and uniqueness of the numerical solutions. The error equation and estimates for velocity and pressure are derived and examined in Section 3. In Section 4, we present numerical examples to verify the theoretical analysis conducted before.

2.   Weak Galerkin finite element method

2.1.   Variational form

Let $\Omega \in\mathbb{R}^d\; (d = 2,\; 3)$ be an open bounded domain with Lipschitz continuous boundary. We denote by $H^m(\Omega)$ the standard Sobolev space, where the associated inner product, boundary product, norm, and semi-norm are given by $\left(\cdot,\cdot\right)_{m}$, $\langle\cdot,\cdot\rangle_{m,\partial \Omega}$, $\|\cdot\|_{m}$, and $|\cdot|_{m}$. When $m=0,$ the space $H^m(\Omega)$ coincides with $L^2(\Omega)$, where the inner product and norm are denoted by $(\cdot,\cdot)_{\Omega}$ and $\|\cdot\|_{\Omega}$, respectively. We denote the subspace of $L^2\left(\Omega\right)$ by $L^2_0\left(\Omega\right)$, where the mean value is zero.

We introduce some bilinear forms and trilinear forms. For any ${\boldsymbol{v}},\; {{\boldsymbol{w}}}\in V$, $p\in W$, we define

If $\nabla\cdot{\boldsymbol{u}}=0$, by the Green's formula,

Then, the variational form of Eqs (1.1)–(1.3) is to find $({\boldsymbol{u}},p)\in V\times W$ such that

for any ${\boldsymbol{v}}\in V$, $q\in W$.

2.2.   Numerical scheme of WG

The main objective of this part is to develop the weak Galerkin finite element numerical scheme for the problems (1.1)–(1.3). The essential step in establishing a numerical scheme for the WG method is to replace the differential operators in the variational form with weak differential operators and incorporate a stabilizer to improve the weak continuity of the approximation functions. This part aims to provide an introduction to fundamental topics in WG techniques, including the local weak function space, discrete weak function space, weak differential operators, and the stabilizer.

Let ${\mathcal{T}}_h$ be the regularity partition of $\Omega$ satisfying conditions^[20]. $T$ denotes the partition unit and $\partial T$ denotes the boundary of $T$. For each $T\in {\mathcal{T}}_h$, $h_T$ denotes its diameter and $h=max_{T\in\mathcal{T}_h}$ denotes the mesh size of ${\mathcal{T}}_{h}$. For any $k\geq 1$, $P_{k}(T)$ denotes the set of polynomials with degree no more than $s$ on $T$. The function ${\boldsymbol{v}}=\{{\boldsymbol{v}}_0,{\boldsymbol{v}}_b\}$ is called the weak vector valued function, where ${\boldsymbol{v}}_0$ and ${\boldsymbol{v}}_b$ are the values of the weak function ${\boldsymbol{v}}$ in the interior $T$ and the boundary $\partial T$, respectively.

Define the local weak function space $S(T)$, discrete velocity space, and pressure space $V_h$, $V_h^0$, and $W_h$, as follows.

Define two types of weak differential operators, namely the discrete weak gradient and discrete weak divergence.

Definition 1. For any $\tau\in [P_{k-1}(T)]^{d\times d}$, discrete weak gradient of vector valued functions $\nabla_w \boldsymbol{v}\in{[P_{k-1}(T)]^{d\times d}}$ is defined as the unique polynomial satisfying the following equation

Definition 2. For any $q\in P_{k-1}(T)$, discrete weak divergence of vector valued functions $\nabla_w\cdot \boldsymbol{v}\in{P_{k-1}(T)}$ is defined as the unique polynomial satisfying the following equation

We denote the stabilizer by bilinear form $s({\boldsymbol{u}}_h,{\boldsymbol{v}}_h)$, which is defined as

Here, we provide the weak Galerkin finite element numerical approach for the problems (1.1)–(1.3) in the following manner:

Algorithm 1. Find $\boldsymbol{u}_h\in V_h^0+H^{1},\ \boldsymbol{p}_h\in W_h+H^{1}\ $such that

for any $(\boldsymbol{v}_h,q_h)\in V_h^0\times W_h$, where

2.3.   Well-posedness of the WG scheme

This subsection provides a proof of the well-posedness for the WG numerical scheme. Next, we present the subsequent norm definition and lemmas from the documents ^{[11,14,21,26]}.

Definition 3. For any $\boldsymbol{v}_h\in{V_h^0}$, define

Definition 4. For any $\boldsymbol{v}_h\in{V_h^0}$, define a mesh-dependent norm

Lemma 1. For any $\boldsymbol{v}_h\in{V_h^0}\ $we have

Lemma 2. $\left[\mathit{\text{Discrete Sobolev Inequality}}\right]$ For any $\boldsymbol{v}_h\in{V_h^0}$, there exists a positive constant $C_0$ independent of $h$ such that

Lemma 3. For any $\rho_{h}\in W_{h}$, there exists a positive constant $\beta$ independent of $h$ such that

Lemma 4. For all $a,b\in{R^n},\ r>2$, there holds that

Lemma 5. Let $Z$ be a finite dimensional functional space equipped with a norm $\|\cdot\|_Z$, let $\theta>0$, and let $\Psi:\; Z\times[0,1]\rightarrow Z$, satisfying the following assumptions

1. $\Psi$ is continuous;

2. For any $(z,\rho)\in Z\times[0,1],\; \Psi(z,\rho)=0$ implies $\|z\|_Z\neq \theta$;

3. $\Psi(\cdot,0)$ is an affine function and the equation $\Psi(Z,0)=0$ has a solution $z\in Z$ such that $\|z\|_Z<\theta$.

There exists $z\in Z$ such that $\Psi(z,1)=0$ and $\|z\|_Z<\theta$.

The next propositions can be deduced from Definition 1, 3 and Cauchy-Schwarz inequality.

Proposition 1. For any $\boldsymbol{u}_h, \boldsymbol{v}_h\in{V_h^0}$, we have

where $C$ and $\ \nu\ $are constant.

The upper bound of the trilinear form $c_h(\cdot;\cdot,\cdot)$ and $d_h(\cdot; \cdot,\cdot)$ can be proved.

Proposition 2. For any $\boldsymbol{u}_h,\boldsymbol{v}_h,\boldsymbol{w}_h\in{V_h^0}\ $, we have

where $C_1,\ C_2\ $are constants independent of $h$.

Proof. From the Lemmas 1 and 2, we have

In a similar manner, we may obtain

□

The theorem presented aims to establish the well-posedness of the WG scheme.

Theorem 1. $\left[\mathit{\text{Existence and Boundedness}}\right]$ There exists a solution $(\boldsymbol{u}_h,p_h)\in V^0_h\times W_h$ to the WG schemes (2.1) and (2.2) that satisfies the following estimates

Proof. Let $Z_h:=V^0_h\times W_h$. It is obvious that $Z_h$ is a finite dimensional functional space. For any $({{\boldsymbol{w}}}_h, r_h)\in Z_h$, the norm $\|({{\boldsymbol{w}}}_h, r_h)\|_{Z_h}$ is given by

Define the continuous map $\Psi: Z_h\times [0,1]\rightarrow Z_h$ such that, for $(({\boldsymbol{w}}_h,r_h),\rho)\in Z_h\times [0,1]$, $\; (\mathit{\boldsymbol{\xi}},\zeta_h)=\Psi(({\boldsymbol{w}}_h,r_h),\rho)$ is defined as the unique polynomial pairs satisfying the following equations

for any ${\boldsymbol{v}}_h\in V^0_h,\; q_h\in W_h$.

We are now in a position to check the conditions of Lemma 5 one by one.

It is easy to know that $\Psi$ is continuous. Let $(\mathit{\boldsymbol{\xi}}_h,\zeta_h)=\Psi(({\boldsymbol{w}}_h,r_h),\rho)=(\boldsymbol{0},0)$. Taking ${\boldsymbol{v}}_h={\boldsymbol{w}}_h$ in Eq (2.12) and $q_h=r_h$ in Eq (2.13), and using Eqs (2.6) and (2.7), we have

The following can be obtaine} from Eq (2.5)

From Eq (2.3), we can further have ${|||} {\boldsymbol{w}}_h{|||}\leq \mu^{-1}\|{\boldsymbol{f}}\|_{*,h}.$ According to Lemma 4,

which implies

Taking

where $\delta$ is a positive constant, $\|({\boldsymbol{w}}_h,r_h)\|_{Z_h}\neq\theta$ .

Note that $\Psi({\boldsymbol{z}}_h,0)$ is an affine map when $\rho=0$. The equation $\Psi({{\boldsymbol{w}}}_h,0)=\boldsymbol{0}$ is well-posed because in the result of WG finite element method for the Stokes problem ^[21], we denote it by $\boldsymbol{z}_h=({\boldsymbol{w}}_h,r_h)\in Z_h$ and the solution satisfies the estimate $\|({\boldsymbol{w}}_h,r_h)\|_{Z_h}< \theta$. There exists $\boldsymbol{z}_h^*=({\boldsymbol{u}}_h,\; p_h)\in Z_h$ such that $\Psi(\boldsymbol{z}_h,1)=0$ and $\|\boldsymbol{z}_h^*\|_{Z_h}< \theta$ from Lemma 5. Thus the WG scheme solution is $({{\boldsymbol{u}}_h,\; p_h})$, which satisfies the estimates (2.10) and (2.11). □

Theorem 2. The solution of the WG schemes (2.1) and (2.2) is unique under the condition $C_2\|\boldsymbol{f}\|_{*,h}< \mu^2$.

Proof. Let $({\boldsymbol{u}}_{h1},p_{h1})\in V^0_h\times W_h, ({\boldsymbol{u}}_{h2},p_{h2})\in V^0_h\times W_h$ solve Eqs (2.1) and (2.2) and satisfy the estimates (2.10) and (2.11). Then, we have

Observe that

So we arrive at

Taking ${\boldsymbol{v}}_h={\boldsymbol{u}}_{h1}-{\boldsymbol{u}}_{h2}$, we have the next inequalities based on Lemma 4, Proposition 1, and Theorem 1.

Thus,

The condition $C_2\|{\boldsymbol{f}}\|_{*,h}< \mu^2$ yields that ${\boldsymbol{u}}_{h1}={\boldsymbol{u}}_{h2}=0$. It follows from Eq (2.11)

Combining with Eq (2.5), we know that $p_{h1}=p_{h2}$. □

3.   Error analysis

In this section, we construct the error equation of the WG numerical scheme and obtain the optimal error estimates for the velocity in an energy norm and the pressure in the $L^2$ norm, which shows the efficiency of the WG method in theory.

We first introduce several $L^2$ projection operators. Let $Q_{h}=\{Q_{0},Q_{b}\}$, where $Q_{0}$ is the $L^2$ space projection onto $[P_{k}(T)]^d$, $Q_{b}$ is the $L^2$ space projection onto $[P_{k}(e)]^d$. Moreover, ${{\boldsymbol{Q}}}_{h}, R_{h}$, and ${{\mathbb{Q}}}_{h}$ represent the projection operators from $L^2$ space onto $[P_{k-1}(T)]^{d\times d}$, $[P_{k-1}(T)]^{d}$, and $P_{k-1}(T)$, respectively.

Proposition 3. Let $1\leq r\leq k$, $\boldsymbol{w}|_{\Omega}\in [H^{r+1}(\Omega)]^d$, $\rho|_{\Omega}\in H^{r}(\Omega)$, $\boldsymbol{v}_h\in V_h$. Thus, under the condition of regular subdivision, we have

where

Proof. Based on the definition of $Q_b$, trace inequality, Cauchy-Schwarz inequality, and projection inequality (A.5), we can get

□ Based on the lemma mentioned above, we provide the subsequent error equations.

Theorem 3. Let $(\boldsymbol{u},p)\in{([H^{k+1}(\Omega)]^d\cap[H^1_0(\Omega)]^d)\times(H^k(\Omega)\cap L^2_0(\Omega))}$ be the exact solution of Eqs (1.1)–(1.3), $(\boldsymbol{u}_h,p_h)\in V_h^0\times W_h$ be the numerical solution of the WG schemes (2.1) and (2.2), and we define $\boldsymbol{e}_{h}=Q_h\boldsymbol{u}-\boldsymbol{u}_{h}$ $\varepsilon_h={{\mathbb{Q}}}_h p-p_{h},$ the following equations

hold for any $\boldsymbol{v}_h \in V_h^0$ and $\ q_h \in W_h$, where

Proof. Multiplying by ${\boldsymbol{v}}_0$ in ${\boldsymbol{v}}=\{{\boldsymbol{v}}_0,{\boldsymbol{v}}_b\}$ to both sides of Eq (1.1), we have

Using integration by parts to get

It follows from Eqs (A.1) and (A.2)

Noting that

Further simplify as

Subtracts Eq (3.3) from Eq (2.1) as

Replacing $({\boldsymbol{u}}\cdot\nabla{\boldsymbol{u}},{\boldsymbol{v}}_0)$ by Eq (A.6), we have

Next, we deal with the nonlinear terms $\delta_d({\boldsymbol{u}};{\boldsymbol{u}},{\boldsymbol{v}}_h)-d_{w}({\boldsymbol{u}}_h;{\boldsymbol{u}}_h,{\boldsymbol{v}}_h)$. First, we get

Furthermore,

Combining Eqs (3.4) and (3.5), we obtain

Multiply both sides of Eq (1.2) by $q_h\in W_h$ and use Lemma A1 as follows

The difference of Eqs (2.2) and (3.6) yields the following equation

□

Next, we are ready to derive the optimal error estimates for velocity and pressure. Taking the velocity function as an example, we can in fact separately estimate the errors $Q_h{\boldsymbol{u}} -{\boldsymbol{u}}_h$ and $Q_h{\boldsymbol{u}} -{\boldsymbol{u}}$ in order to estimate ${\boldsymbol{u}} -{\boldsymbol{u}}_h$. Because the interpolation error $Q_h{\boldsymbol{u}} -{\boldsymbol{u}}$ is proved by the finite element theory in ^[20], we need to demonstrate the error estimates theorem as follows.

Theorem 4. Let $(\boldsymbol{u},p)\in ([H^{k+1}(\Omega)]^2\cap [H^1_0(\Omega)]^2)\times (H^k(\Omega)\cap L^2_0(\Omega))$ be the solution of Eqs (1.1)–(1.3), and $({${\textbf{u}}$}_h,{${\textbf{v}}$}_h)$ be the solution of the WG schemes (2.1) and (2.2). Then, the following error estimate holds true

Proof. Taking ${\boldsymbol{v}}_h={\boldsymbol{e}}_h$ in Eq (3.1), $q_h=\varepsilon_h$ in Eq (3.2), and adding $\alpha(|Q_0{\boldsymbol{u}}|^{r-2}Q_0{\boldsymbol{u}},{\boldsymbol{e}}_0)_h$ to both sides of Eq (3.1), we observe

Based on Eqs (2.8) and (2.9), and Lemma 4, we have

The following is an estimate of each item of the right side of the above formula. First, we deal with the linear terms. According to Eqs (2.3) and (2.9), we have

From Proposition 3, we have

Using $H\ddot{o}lder$ inequality, Discrete Sobolev inequality (2.4), and projection inequalities (A.3) and (A.4),

Furthermore, using trace inequality, we have

At last, we deal with $\alpha(|{\boldsymbol{u}}|^{r-2}{\boldsymbol{u}},{\boldsymbol{e}}_0)_h+\alpha(|Q_0{\boldsymbol{u}}|^{r-2}Q_0{\boldsymbol{u}},{\boldsymbol{e}}_0)_h$. Using $H\ddot{o}lder$ inequality, Lemma (2.4), and projection inequality (2.10), we get

Thus, we get the following error estimate of velocity.

In order to get the estimates of $p$, we can start with Eq (3.2)

On the right-hand side of the above equation, we only need to estimate the terms $\alpha(|{\boldsymbol{u}}|^{r-1}{\boldsymbol{u}},{\boldsymbol{v}}_0)-c_h({\boldsymbol{u}}_h;{\boldsymbol{u}}_h,{\boldsymbol{v}}_0)$ and $d_{h}({\boldsymbol{e}}_h;Q_h{\boldsymbol{u}};{\boldsymbol{v}}_h)+d_{hw}({\boldsymbol{u}}_h;{\boldsymbol{e}}_h,{\boldsymbol{v}}_h)$ because estimates of other terms have been estimated in the previous proof. Using Lemma 1, Proposition 1, and $H\ddot{o}lder$ inequality,

From Sobolev Inequality (2.4), we have

Thus,

Combing with inf-sup condition (2.5), the estimate of $p$ is given as

□

4.   Numerical experiments

Three numerical examples are used to assess the efficacy of the WG finite element numerical solution in this section. Taking into account the nonlinear term of the WG numerical scheme of Navier-Stokes equation with the nonlinear terms (2.1) and (2.2), we use the Oseen iterative method ^[11] to linearize it first: assuming that we have $({\boldsymbol{u}}^n_h,p^n_h)$, find $({\boldsymbol{u}}^{n+1}_h,p^{n+1}_h)\in V^0_h\times W_h$ such that

where the termination condition is set to $1.0E-6$.

The first two simulations are provided with precise, smooth solutions, while the third simulation involves a stable nonlinear fluid and is considered a quasi-basis problem. Each of these simulations is executed using MATLAB.

Example 1. Let $\Omega=(0,1)\times (0,1)$, viscosity coefficient $\mu=1$, damping coefficient $\; r=3,\; \alpha=1$. The source term $f$ is selected to ensure that the exact solution is

We employ a uniform rectangular grid partitioning of the region $\Omega$ in Example 1 and subsequently establish connections between their diagonals to produce triangulated meshes. We proceed to mesh in steps of $h=1/4, 1/8, 1/16, 1/32, 1/64$. The WG numerical technique approximates the velocity function using a polynomial function of degree one within the element and a polynomial function of degree zero on the boundary of the element. The pressure function is approximated by a first-degree polynomial function within the element. Namely, the WG scheme is in case of $P_{1}P_{1}-P_{0}$.

The results of Example 1 are illustrated in Figures 1–4 for a mesh size of $h=1/64$. Figure 1(a) displays the numerical solutions of the velocity function ${\boldsymbol{u}}_{h}$ in the first direction, whereas Figure 1(b) displays the exact solutions of the velocity function ${\boldsymbol{u}}$ in the first direction. Figure 2(a),(b) shows the numerical solutions and the exact solutions of the velocity function in the second direction. The numerical solutions $p_h$ and exact solutions $p$ of the pressure function are displayed in Figure 3. The left and right subgraphs of Figure 4 illustrate the two contour maps of pressure. Observing the images above, it becomes apparent that the numerical solutions for the velocity and pressure functions are able to accurately approximate the respective analytical solutions.

Subsequently, we show the error and convergence of the velocity and pressure functions. The convergence rates with respect to $h$ are calculated using the formula $log(e_{i}/e_{i+1})/log(h_{i+1}/h_{i})$, where $e_{i}$ and $e_{i+1}$ are the relative errors corresponding to the mesh sizes $h_{i}$ and $h_{i+1}$. Table 1 denotes the obtained error and convergence order. The third and fifth columns of the table clearly demonstrate that the rate of convergence for the $H^1$-norm and $L^2$-norm of velocity is of first and second order, respectively. The outcomes coincide precisely with the predictions made by Eq (3.7) in Theorem 4. Furthermore, it is evident that the seventh column of Table 1 demonstrates a convergence rate over order 1 for the $L^2$-norm of the pressure function. Also, we show the image of the convergence order in Figure 5, which is compared with lines of first and second order visually. This indicates that the convergence is superlative and surpasses the expectations set by the theoretical study.

Example 2. Let $\Omega=(0,1)\times (0,1)$, viscosity coefficient $\mu=1$, damping coefficient $\; r=5,\; \alpha=2$. The source term $f$ is selected to ensure that the exact solution is

In Example 2, we ensure that the mesh partitions, mesh sizes, and the weak Galerkin finite element function space are identical to those in Example 1. In this example, we calculate the numerical solutions for the velocity ${\boldsymbol{u}}_{h}$ and pressure $p_{h}$ when the value of $h$ is set to $1/64$. The numerical solutions of the velocity functions in the first and second direction are shown in Figures 6(a) and 7(a). The exact solutions of the velocity functions in the first and second direction are shown in Figures 6(b) and 7(b). The numerical solutions $p_h$ are presented in Figure 8(a), and the exact solutions $p$ are displayed in Figure 8(b). Figure 9(a),(b) displays the pressure contours. By analyzing these figures, we may obtain numerical solutions for both velocity and pressure that closely approximate the corresponding exact results. In addition, we compute the errors and convergence and present them in Table 2 and Figure 5(b). The table displays the rate of convergence for the $H^1$-norm error of velocity in the third, fifth, and seventh columns, as well as the $L^2$-norm error of velocity. The $H^1$-norm error of the pressure exhibits orders of 1, 2, and 1, which aligns with the theoretical analysis presented in Theorem 3.7.

We are also interested in determining whether variations in the damping coefficient impact the numerical results of the WG method proposed in this paper. By a number of experiments in which the coefficients change by orders of magnitude, it shows that when the damping coefficient changes are extremely small, even up to $10E-06$, the effect on the solution of the equation still remains small, which is up to $4.07\%$. However, when $\alpha$ increases, the numerical results are effected significantly. In practice, the damping parameters describe the damping force that is not directly proportional to the velocity of a flowing fluid. Large damping coefficients lead to an increase in the proportion of nonlinear terms in the equation, and this causes change in calculation results of ${\boldsymbol{u}}_h$.

Example 3. Consider the lid-driven cavity flow problem ^[7]. Let $\Omega=(0,1)\times (0,1)$, viscosity coefficient $\mu=0.1$, damping coefficient $\; r=3,\; \alpha=1$. The source term $f$ is zero. The velocity function satisfies the Dirichlet boundary condition, namely, we impose the normal component of the velocity to be zero on $\partial\Omega$ and the tangential component to be zero except along the top boundary where it is set to one.

In Example 3, we standardize the triangulation partitioning to match that of the initial two examples. The WG scheme is in case of $P_{1}P_{1}-P_{0}$. When mesh size $h=1/64$, we represent the numerical velocity vector diagram as Figure 10(a) and the streamline diagram as Figure 10(b). By analyzing Figure 10(a), it is evident that the initial velocity of the top border in the tangential direction induces a clockwise flow of the fluid in the cavity. Upon analyzing Figure 10(b), it becomes evident that, as one moves closer to the top cover, the streamline becomes thicker, indicating a higher velocity. This has the same calculation results as ^[11]. The density of streamlines decreases and the velocity of the flow decreases as the flow approaches the three boundaries.

5.   Results

In this paper, we propose the weak Galerkin finite element method for the Navier-Stokes equation with a nonlinear damping term. The theoretical analysis demonstrates the existence and uniqueness of numerical solutions obtained by the WG finite element method. We analyse the convergence of the velocity error in the energy norm and the pressure error in the $L^2$-norm. The numerical experiments show the validity of the theoretical study and the efficacy of the WG method. For improved computational efficiency, this two-level technique can be explored in the future, and we will make an effort for the time-dependent models.

Use of AI tools declaration

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

Conflict of interest

The authors declare there is no conflict of interest.

Authors contribution

Yue Tai: Methodology and writing original draft; Xiuli Wang: Formal analysis; Weishi Yin: Review and editing; Pinchao Meng: Supervision and validation.

Appendix A

In order to give the error estimates, we need to estimate the remainder $s(Q_h\textbf{w},\textbf{v}_h),\varphi_{\textbf{w}}(\textbf{v}_h)$ and $\theta_{\rho}(\textbf{v}_h)$ with the help of the lemmas in ^[21] and ^[26]. In this section, we shall introduce these techniques.

Lemma A1. The projection operators $Q_h, {{\textbf{Q}}}_h,$ and ${{\mathbb{Q}}}_h$ satisfy the following commutative properties

Lemma A2. For any $\boldsymbol{w}\in [H^1(\Omega)]^d,$ $\rho\in H^1(\Omega), \boldsymbol{v}_h\in V^0_h$, it follows that

Lemma A3. Let ${\mathcal{T}}_h$ be a finite element partition of domain $K$ satisfying the shape regularity assumptions and $\boldsymbol{w}\in[H^{r+1}(K)]^d$ and $\rho\in H^r(K)$ with $1\le r\le k$. Then, for $0\le s\le 1$, we have

here, $C$ denotes a generic constant independent of the mesh size $h$ and the functions in the estimates.

Lemma A4. Let $\boldsymbol{u}\in [H^1_0(\Omega)]^d$ and $\nabla\cdot\boldsymbol{u}=0$, we have

where