**1. Introduction**

Many mathematical models of natural processes are described by the boundary value problems for systems of partial differential equations with a singularity. The singularity of the solution to such systems in the two-dimensional closed domain Ω may be due to the degeneration of initial data, to the presence of reentrant corners on a boundary, or to internal features of the solution. The boundary value problem has a strong singularity if its solution does not belong to the Sobolev space *W*<sup>1</sup> 2 (Ω). In short, the Dirichlet integral from the solution diverges. In the case when the solution belongs to the space *W*<sup>1</sup> 2 (Ω), but it does not belong to the *W*<sup>2</sup> 2 (Ω), a boundary value problem is called weakly singular. The generalized solution of a boundary value problem in the two-dimension domain with a boundary containing an initial angle *ω* belongs to the space *W*1+*α*− 2 (Ω), where 0.25 ≤ *α* < 1 for *π* < *ω* ≤ 2*π* and  is an arbitrary positive real number. Therefore, the approximate solution produced by the classical finite difference or finite element methods converges to an exact one no faster than at the O(*hα*) rate [1].

For the boundary value problem with singularity, there are various numerical approaches founded on the separation of singular and regular components of the generalized solution, on mesh refinement toward singularity points, and on the multiplicative identification of singularities. These methods slow down the convergence rate of the approximate solution to an exact one or to the significant complication of the finite element scheme, which in total influences the computational process speed and accuracy of the result.

In reference [2], we suggested to define the solution of the boundary value problem with weak or strong singularity as an *<sup>R</sup>ν*-generalized one in the weighted Sobolev space or set. Relying on this approach, numerical methods were created with a convergence rate independent of the value (size) of a singularity. In the papers [3–5] for the boundary value problems with a strong singularity, the weighted finite element method (FEM) and the weighted edge-based FEM were built. The approximate solution converges to an exact one with the second and first order rates (under the mesh step *h*) in the norms of the weighted Lebesgue and Sobolev spaces, respectively. In references [6,7], a weighted FEM for the

Lame system in a domain with the reentrant corner on the boundary was built. The rate of convergence is equal to O(*h*) and independent of the size of a reentrant corner.

We study the incompressible Navier–Stokes equations in the two-dimensional polygonal domain Ω with one internal corner greater than 180◦ on its boundary. The nonlinearity in this system can be written in several equivalent forms. For one case, if we regard these equations in the velocity field and kinematic pressure variables, then this leads to the convection form of nonlinear terms. For another case, if we consider these equations in the velocity field and total pressure variables, then it gives nonlinear terms in the rotation form. In order to meet the non-stationary incompressible system, we must be able to find the solution of a steady linearized one. The stationary Navier–Stokes system we can linearize in different manners. We use a scheme that is based on Picard's iterative procedure (see [8] and the references therein). Starting with an arbitrary vector as a velocity field, which satisfies the law of conservation of mass, Picard's iterative procedure forms the sequence of solutions of the corresponding linear Oseen system. We note that linearizations of convection and rotation forms of nonlinear terms tend to the systems of linear algebraic equations with various features. In the paper, we study the Oseen system in the rotation form. The fact is that the rotation form allows us (using a skew-symmetric of the resulting matrix) to construct a Schur complement preconditioner, which is acceptable to all parameters of the Oseen problem and becomes more effective for large Reynolds numbers (see [9] and the references therein). For the convection form of the Oseen problem, this is not so.

As usual, to solve a fluid problem, the explorer has freedom and can construct a method in different manners by selecting various discretization algorithms for the system of linear algebraic equations. There are many opportunities to solve the considered system. The researcher can select various finite difference, finite volume, or finite element methods. However, the chosen method is effective if it gives the best result in terms of the convergence rate under certain restrictions on the input data and geometric singularities of the domain Ω.

In the paper, we consider a special case, where Ω is a polygon with one internal corner greater than 180◦ on its boundary. The flow of the viscous fluid in a *δ*-neighborhood of a reentrant angle was studied in [10]. It is not a secret that the velocity field and pressure, as a weak solution of a problem for the domain with corner singularity, do not belong to Sobolev spaces **W**<sup>2</sup> 2(Ω) and *W*<sup>1</sup> 2 (Ω), respectively [11]. Therefore, the rate of convergence of the approximate solution to an exact one is equal to O(*hα*), *α* < 1, in the norm of standard and weighted Sobolev spaces (see [12] and the references therein) for different classical finite difference and finite element methods. Earlier, for the Stokes problem, we defined the *<sup>R</sup>ν*-generalized solution; in [13], we formulated and proved the weighted LBBinequality (inf-sup condition [14]); and in [15], we showed the advantage of our method over classical approaches.

The aim of the paper is to present a new numerical approach for the rotation form of the Oseen problem using (see [16]) a mass conservation space pair; to show that the rate of convergence of the approximate solution to an exact one (the velocity field) is equal to O(*h*) for all considered sizes of the internal corner greater than 180◦ on the boundary in the norm of the space **<sup>W</sup>**12,*ν*(<sup>Ω</sup>*k*); so that this rate is much better than if using the classical finite difference or finite element methods.

The article consists of six sections. Section 2 is devoted to the definition of the *<sup>R</sup>ν*-generalized solution for the rotation form of the Oseen system in a domain Ω with one internal corner greater than 180◦ on its boundary. In Section 3, we construct the presented FEM. The iterative algorithm for the resulting system of linear algebraic equations is built in Section 4. In Section 5, we discuss the numerical results of computational experiments. Necessary conclusions are made in Section 6.

### **2.** *Rν***-Generalized Solution of the Oseen Problem**

Let **x** = (*<sup>x</sup>*1, *<sup>x</sup>*2) be an element of the Euclidean space **R**2, where **x** = *x*2 1 + *x*2 2 1/2 and *d***x** = *dx*1 *dx*2 are the norm and measure of **x**, respectively. Denote by Ω a bounded domain in **R**2. Let Γ and Ω¯ be the boundary and closure of Ω, respectively, where Ω¯ = Ω ∪ Γ.

At first, we write incompressible Navier–Stokes equations in such a form: find a velocity field **<sup>u</sup>**(**<sup>x</sup>**, *<sup>t</sup>*)=(*<sup>u</sup>*1(**<sup>x</sup>**, *t*), *<sup>u</sup>*2(**<sup>x</sup>**, *t*)) and a kinematic pressure *p*(**<sup>x</sup>**, *t*) from:

$$\frac{\partial \mathbf{u}}{\partial t} - \overline{v} \triangle \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p = \mathbf{f} \quad \text{and} \qquad \text{div } \mathbf{u} = \mathbf{0} \quad \text{in} \qquad \Omega \times (0, T], \tag{1}$$

with given force field **f** = (*f*1, *f*2) and viscosity *ν*¯ = 1*Re* > 0. Let , div , and ∇ be the Laplace, divergence, and gradient operators in **R**2, respectively. The equations in (1) are the convection form of Navier–Stokes equations.

We supplement the system (1) with a boundary and initial conditions:

$$\mathbf{u} = \mathbf{g} \quad \text{on} \quad \Gamma \times (0, T], \qquad \mathbf{u}(\mathbf{x}, 0) = \mathbf{u}^0(\mathbf{x}) \quad \text{in} \quad \Omega,\tag{2}$$

where **g** = (*g*1, *g*2) is given vector function on Γ and **<sup>u</sup>**<sup>0</sup>(**x**)=(*u*01(**x**), *u*02(**x**)) — in Ω. We introduce the following notation:

$$\mathbf{v} \cdot \mathbf{u} = \sum\_{i=1}^{2} u\_i v\_{i\prime} \quad \text{curl } \mathbf{u} = -\frac{\partial u\_1}{\partial \mathbf{x}\_2} + \frac{\partial u\_2}{\partial \mathbf{x}\_1}, \quad a \times \mathbf{u} = \left( \begin{array}{c} -au\_2 \\ au\_1 \end{array} \right).$$

We have a formal equality:

$$\nabla(\mathbf{u}\cdot\mathbf{v}) + (\operatorname{curl}\,\mathbf{u})\times\mathbf{v} + (\operatorname{curl}\,\mathbf{v})\times\mathbf{u} = (\mathbf{u}\cdot\nabla)\mathbf{v} + (\mathbf{v}\cdot\nabla)\mathbf{u}.\tag{3}$$

If **u** = **v** in (3), then we have a relation:

$$(\operatorname{curl} \mathbf{v}) \times \mathbf{v} + \frac{1}{2} \nabla \mathbf{v}^2 = (\mathbf{v} \cdot \nabla) \mathbf{v}. \tag{4}$$

.

Let *P* = *p* + 12**u**2, using (4), for vector function **u**; we ge<sup>t</sup> the rotation form of the Navier–Stokes system for an incompressible flow:

$$\frac{\partial \mathbf{u}}{\partial t} - \overline{\nu} \triangle \mathbf{u} + (\text{curl } \mathbf{u}) \times \mathbf{u} + \nabla P = \mathbf{f} \text{ and } \quad \text{div } \mathbf{u} = \mathbf{0} \text{ in } \Omega \times (0, T]. \tag{5}$$

We supplement the system (5) with the boundary and initial conditions (2). Using implicit time integration of (5) compared to explicit methods reduces accuracy, stability, and flexibility in selecting the step size for a time variable.

In our research, on each time level, we solve the following system of equations:

$$-\overline{\nu}\triangle\mathbf{u} + \text{curl } \mathbf{u} \times \mathbf{u} + \mathfrak{a}\ \mathbf{u} + \nabla P = \mathbf{f} \text{ and } \quad \text{div } \mathbf{u} = \mathbf{0} \quad \text{in} \quad \Omega,\tag{6}$$

$$\mathbf{u} = \mathbf{g} \qquad \text{on} \qquad \Gamma\_{\prime} \tag{7}$$

and parameter *α* is a known positive constant.

The system (6) and (7) is nonlinear due to the fact that there is a rotation term curl **u** × **u** in the first Equation (6). This term and the system as a whole we linearized by Picard's iterative procedure (see [8] and the references therein).

At each iteration, we need to solve the following problem:

$$-\nabla \triangle \mathbf{u} + \mathbf{u} \times \mathbf{u} + \mathbf{a} \cdot \mathbf{u} + \nabla P = \mathbf{f}, \quad \text{and} \quad \text{div } \mathbf{u} = \mathbf{0} \quad \text{in} \quad \Omega,\tag{8}$$

$$
\mathbf{u} = \mathbf{g} \qquad \text{on} \qquad \Gamma\_{\prime} \tag{9}
$$

which is called the Oseen system in a rotation form, where *w* = curl **U** and **U** is some approximation to **u**.

The linearization of convection and rotation forms of nonlinear terms tends to the systems of linear algebraic equations with various features. In the paper, we study the Oseen system in the rotation form. The fact is that the rotation form allows us (using a skew-symmetric of the resulting matrix) to construct the Schur complement preconditioner, which is acceptable to all parameters of the Oseen problem and becomes more effective when *ν*¯ → 0 (see [9] and the references therein). For the convection form of the Oseen problem, this is not so.

We note that for the linearized system (8) and (9), the laws of the conservation of momentum and mass remain valid.

In the paper, we consider a special case, where Ω is a bounded non-convex polygonal domain with one internal corner greater than 180◦ on Γ. Let its vertex be located at the origin. We define an *<sup>R</sup>ν*-generalized solution of the Oseen problem (8) and (9) with a corner singularity and construct the weighted FEM. We demonstrate the advantage of the proposed approach over the classical finite element methods for all sizes of the reentrant corner.

Let Ω*δ* = {**x** ∈ Ω¯ : **x** ≤ *δ*, *δ* ∈ (0, 1)} be a part of a *δ*-neighborhood, with a vertex located at the origin, which is in Ω ¯ . Denote by *ρ*(**x**) a weight function: *ρ*(**x**) = - **x** , **x** ∈ <sup>Ω</sup>*δ*, *δ* , **x** ∈ Ω¯ \ <sup>Ω</sup>*δ*.

Let *Dmv*(**x**) = *∂*|*m*| *v*(**x**) *∂xm*<sup>1</sup> 1 *∂xm*<sup>2</sup> 2 be the *m*th order generalized derivatives of a function *v*(**x**) in Ω, where |*m*| = *m*1 + *m*2, *mi*, nonnegative integers. For the function *<sup>v</sup>*(**x**), we define the following inequalities:

$$\int\_{\Omega \cap \Omega\_{\delta}'} \rho^{2a} v^2 d\mathbf{x} \ge C\_1 > 0,\tag{10}$$

$$|D^m v(\mathbf{x})| \le C\_2 \left(\frac{\delta}{\rho(\mathbf{x})}\right)^{a+m} \text{ for } \mathbf{x} \in \Omega\_\delta^{'} \text{ and } m = 0, 1,\tag{11}$$

where *α* > 0 and constant *C*2 > 0 do not depend on *m* and *α*.

Denote by *<sup>L</sup>*2,*α*(Ω) a space of functions *<sup>v</sup>*(**x**), such that:

$$\|v\|\_{L\_{2,a}(\Omega)} = \left(\int\_{\Omega} \rho^{2\alpha} v^2 d\mathbf{x}\right)^{1/2} < \infty.$$

If **w** = (*<sup>w</sup>*1, *<sup>w</sup>*2) is a vector function, then we define the weighted vector function space **<sup>L</sup>**2,*α*(Ω) with a norm **w <sup>L</sup>**2,*α*(Ω) = *<sup>w</sup>*1 2*<sup>L</sup>*2,*α*(Ω) + *<sup>w</sup>*2 2*<sup>L</sup>*2,*α*(Ω)1/2.

Further, denote by *<sup>L</sup>*2,*α*(<sup>Ω</sup>, *<sup>δ</sup>*), *α* > 0, a set of elements *v*(**x**) from the *<sup>L</sup>*2,*α*(Ω) space for which Inequalities (10) and (11) (the case *m* = 0) are valid with a bounded *<sup>L</sup>*2,*α*(Ω) norm. Let *<sup>L</sup>*02,*α*(<sup>Ω</sup>, *δ*) be a subset of functions *<sup>v</sup>*(**x**), such that *<sup>L</sup>*02,*α*(<sup>Ω</sup>, *δ*) = {*v* ∈ *<sup>L</sup>*2,*α*(<sup>Ω</sup>, *δ*) : Ω *ρα vd***x** = <sup>0</sup>}. If **w** = (*<sup>w</sup>*1, *<sup>w</sup>*2) is a

vector function, then we define a set **<sup>L</sup>**2,*α*(<sup>Ω</sup>, *δ*) = {**w** : *wi* ∈ *<sup>L</sup>*2,*α*(<sup>Ω</sup>, *δ*)} with a bounded **<sup>L</sup>**2,*α*(Ω) norm. Let *<sup>W</sup>*12,*α*(Ω) be a weighted space of functions *<sup>v</sup>*(**x**), such that:

$$\|\|\boldsymbol{\upsilon}\|\|\_{W^1\_{2,a}(\Omega)} = \left(\sum\_{|\boldsymbol{m}| \le 1} \|\|\boldsymbol{\rho}^a |D^{\boldsymbol{m}}\boldsymbol{\upsilon}|\|\_{L\_2(\Omega)}^2\right)^{1/2} < \infty.$$

If **w** = (*<sup>w</sup>*1, *<sup>w</sup>*2) is a vector function, then we denote by **<sup>W</sup>**12,*α*(Ω) the weighted vector function space with a norm **w <sup>W</sup>**12,*α*(Ω) = *<sup>w</sup>*1 2*<sup>W</sup>*12,*α*(Ω) + *<sup>w</sup>*2 2*<sup>W</sup>*12,*α*(Ω)1/2.

Let *<sup>W</sup>*12,*α*(<sup>Ω</sup>, *<sup>δ</sup>*), *α* > 0, be a set of functions *v*(**x**) from the space *<sup>W</sup>*12,*α*(Ω), that meet the conditions (10) and (11) with a bounded *<sup>W</sup>*12,*α*(Ω) norm. We denote by *o <sup>W</sup>*12,*α* (<sup>Ω</sup>, *δ*)( *o <sup>W</sup>*12,*α* (<sup>Ω</sup>, *δ*) ⊂ *<sup>W</sup>*12,*α*(<sup>Ω</sup>, *δ*)) a closure, with respect to the *<sup>W</sup>*12,*α*(Ω) norm, of the set of infinitely-differentiable functions with compact support in Ω that meet the inequalities (10) and (11). Then, we denote by *W*1/2 2,*α* (<sup>Γ</sup>, *δ*) the set of functions *<sup>θ</sup>*(**x**) on Γ : *<sup>θ</sup>*(**x**) ∈ *W*1/2 2,*α* (<sup>Γ</sup>, *<sup>δ</sup>*), if there exists a function <sup>Θ</sup>(**x**) from the set *<sup>W</sup>*12,*α*(<sup>Ω</sup>, *<sup>δ</sup>*), such that <sup>Θ</sup>(**x**)|Γ = *<sup>θ</sup>*(**x**) and *θ W*1/2 2,*α* (<sup>Γ</sup>,*<sup>δ</sup>*) = inf <sup>Θ</sup>|Γ=*<sup>θ</sup>* Θ *<sup>W</sup>*12,*α*(<sup>Ω</sup>,*<sup>δ</sup>*).

If **w** = (*<sup>w</sup>*1, *<sup>w</sup>*2) is a vector function, then we define a set **<sup>W</sup>**12,*α*(<sup>Ω</sup>, *δ*) = {**w** : *wi* ∈ *<sup>W</sup>*12,*α*(<sup>Ω</sup>, *δ*)} with a norm of space **<sup>W</sup>**12,*α*(Ω). Similarly, we define the set *o* **<sup>W</sup>**12,*α* (<sup>Ω</sup>, *δ*) of vector functions in Ω and **W**1/2 2,*α*(<sup>Γ</sup>,*<sup>δ</sup>*),onΓ.

Let known functions *w*, **f** = (*f*1, *f*2) and **g** = (*g*1, *g*2) in (8) and (9) meet the following conditions:

$$w \in L\_{2,\gamma}(\Omega, \delta), \quad \mathbf{f} \in \mathbf{L}\_{2,\gamma}(\Omega, \delta), \quad \mathbf{g} \in \mathbf{W}^{1/2}\_{2,\gamma}(\Gamma, \delta), \ \gamma \ge 0. \tag{12}$$

Bilinear and linear forms are as follows:

$$a(\mathbf{u}\_{\nu}, \mathbf{v}) = \int\_{\Omega} \left[ \mathbb{P} \nabla \mathbf{u}\_{\nu} \cdot \nabla (\rho^{2\nu} \mathbf{v}) + \rho^{2\nu} (w \times \mathbf{u}\_{\nu}) \cdot \mathbf{v} + a \rho^{2\nu} \mathbf{u}\_{\nu} \cdot \mathbf{v} \right] d\mathbf{x},$$

$$b(\mathbf{v}, P\_{\nu}) = -\int\_{\Omega} P\_{\nu} \text{ div } (\rho^{2\nu} \mathbf{v}) d\mathbf{x}, \ c(\mathbf{u}\_{\nu}, q) = -\int\_{\Omega} \rho^{2\nu} q \text{ div } \mathbf{u}\_{\nu} \, d\mathbf{x}, \ l(\mathbf{v}) = \int\_{\Omega} \rho^{2\nu} \mathbf{f} \cdot \mathbf{v} d\mathbf{x}.$$

**Definition 1.** *The pair* (**<sup>u</sup>***ν*(**x**), *<sup>P</sup>ν*(**x**)) ∈ **<sup>W</sup>**12,*ν*(<sup>Ω</sup>, *δ*) × *<sup>L</sup>*02,*ν*(<sup>Ω</sup>, *δ*) *is called the <sup>R</sup>ν-generalized solution for an Oseen system in the rotation form (8) and (9) such that for all pairs* (**v**(**x**), *q*(**x**)) ∈ *o* **<sup>W</sup>**12,*ν* (<sup>Ω</sup>, *δ*) × *<sup>L</sup>*02,*ν*(<sup>Ω</sup>, *<sup>δ</sup>*)*, theequalities:*

$$\begin{aligned} a(\mathbf{u}\_{\nu}, \mathbf{v}) + b(\mathbf{v}\_{\nu} P\_{\nu}) &= l(\mathbf{v})\_{\nu}, \\ c(\mathbf{u}\_{\nu}, q) &= 0 \end{aligned}$$

*hold, where functions w*, **f** *and* **g** *satisfy the conditions (12) and ν* ≥ *γ*.

Note that the bilinear and linear forms in the definition of an *<sup>R</sup>ν*-generalized solution include a weight function *ρ*(**x**). The introduction of the weight function into integral identities suppresses the influence of the singularity in the solution and ensures that **u***ν* and *Pν* belong to the weighted sets **<sup>W</sup>**22,*ν*(<sup>Ω</sup>, *δ*) and *<sup>W</sup>*12,*ν*(<sup>Ω</sup>, *<sup>δ</sup>*), respectively. This property of the *<sup>R</sup>ν*-generalized solution allows one to construct a finite element scheme with a O(*h*) rate. This rate is significantly higher than in the classical finite element method for the Oseen problem in a polygonal domain with the internal corner greater than 180◦ on the boundary.

### **3. The Weighted Finite Element Scheme**

Now, we construct a finite element scheme for an Oseen problem in the rotation form (8) and (9) based on the definition of an *<sup>R</sup>ν*-generalized solution.

We would like to use the finite element space pair, which satisfies the law of mass conservation not in the weak (like the well-known Taylor–Hood (TH) element pair [14]), but in the strong sense. The fact is that the implementation of the mass conservation law in a weak sense combines pressure and velocity field errors and does not eliminate possible instabilities [17]. In the paper, we apply the Scott–Vogelius (SV) element pair [16] that will help us to obtain strong mass conservation of the approximate solution.

First, we divide Ω into a finite quantity of triangles *Li*, which we call macro-elements. The set of elements *Li* represents a quasi-uniform (see [1]) triangulation *Th* of Ω. Then, we divide each macro-element *Li* ∈ *Th* into three triangles *Kij* using the barycenter of *Li*. Thus, we construct a triangulation Υ*h*, which is based on a barycenter refinement of a triangulation *Th*. Denote by Ω*h* the set of resulting triangles (which are called finite elements) with sides of order *h*, i.e., Ω*h* = *Kij*<sup>∈</sup>Υ*h Kij* =

$$\bigcup\_{L\_i \in T\_h} \left( \bigcup\_{j=1}^3 \mathcal{K}\_{i\_j} \right) = \bigcup\_{\substack{L\_i \in T\_h\\1 \le D < 1}} L\_i.$$

Let *Am* and *Bl* be vertices and midpoints of the finite element sides *Kij* ∈ Υ*h*, respectively. Then, for the components of a velocity field and pressure, we define sets of nodes *G* and *H*, respectively, such that *G* = *G* Ω *G*Γ = {*Am* ∪ *Bl*}, where *G* Ω is a totality of Υ*h* nodes in Ω, *G*Γ, on Γ, and *H* = {*Ck*}, where *Ck* coincide with a node *Am* on the appropriate element *Kij*∈ Υ*h* (see Figure 1).

Now, we define spaces of the SV element pair. The space *Xh*, for the components of the velocity field, coincides with the corresponding space of degree two of the THelement pair, i.e., *X<sup>h</sup>* = {*w<sup>h</sup>* ∈ *C*(Ω) : *<sup>w</sup><sup>h</sup>*|*Ki j* ∈ *<sup>P</sup>*2(*Kij* ), ∀ *Kij* ∈ <sup>Υ</sup>*h*} and for a velocity field **X***<sup>h</sup>* = *X<sup>h</sup>* × *Xh*. The space *<sup>Y</sup>h*, for the pressure, differs from the corresponding space degree one of the TH element pair by the fact that it is discontinuous in Ω, i.e., *Y<sup>h</sup>* = {*y<sup>h</sup>* ∈ *<sup>L</sup>*2(Ω) : *<sup>y</sup><sup>h</sup>*|*Ki j* ∈ *<sup>P</sup>*1(*Kij* ), ∀ *Kij* ∈ Υ*h*, Ω*yhd***x** = <sup>0</sup>}.

**Figure 1.** The macro-element *Li*: squares and dots are the velocity and pressure nodes on *Kij* , *j* = 1, 2, 3, respectively.

The SV element pair has an important property, namely div **X***<sup>h</sup>* ⊂ *Yh*. This means that there exists a function *yh* ∈ *Y<sup>h</sup>* equal to div **w***<sup>h</sup>* such that: from the condition for performing mass conservation in a weak sense, i.e., Ω div **w***<sup>h</sup>ψhd***x** = 0 ∀ *ψh* ∈ *<sup>Y</sup>h*, we ge<sup>t</sup> a pointwise mass conservation, i.e., div **w***<sup>h</sup> <sup>L</sup>*2(Ω) = 0. Moreover, in [18], it was established that spaces of the SV element pair before us thecondition.Note,thatobtainedthe˘

satisfy Ladyzhenskaya–Babuska–Brezzi approximations using TH element, pair unlike the SV element pair, in general, do not achieve pointwise mass conservation. Then, we define the weighted basis functions and describe a special finite element method for

the Oseen system in the rotation form (8) and (9). For components of the velocity field, for each node *Mk* ∈ *G* Ω, we will match a function:

$$\Phi\_k(\mathbf{x}) = \rho^{\nu^\diamond}(\mathbf{x}) \cdot q\_k(\mathbf{x}), \ k = 0, 1, \dots, n$$

where *ϕk*(**x**) ∈ *<sup>X</sup>h*, *<sup>ϕ</sup>k*(*Mj*) = -1, *k* = *j*, 0, *k* = *j*, *k*, *j* = 0, 1, . . . ; *ν* is a parameter.

We define a set *<sup>V</sup>h*, for components of the velocity field, such that for any velocity field **v***<sup>h</sup>* = (*v<sup>h</sup>* 1, *vh* 2), *vh i*∈ *<sup>V</sup>h*, we have:

$$v\_1^h(\mathbf{x}) = \sum\_k d\_{2k} \, \Phi\_k(\mathbf{x}), \ v\_2^h(\mathbf{x}) = \sum\_k d\_{2k+1} \, \Phi\_k(\mathbf{x}), \tag{13}$$

where *dl* = *ρ*<sup>−</sup>*ν* (*<sup>M</sup>*[*l*/2]) ˜*dl*.

Let *Vh*0 be a subset in *V<sup>h</sup>* such that *Vh*0 = {*w<sup>h</sup>* ∈ *V<sup>h</sup>* : *<sup>w</sup><sup>h</sup>*(*Mk*)|*Mk*∈*G*Γ = <sup>0</sup>}. Moreover, we define velocity field sets **V***<sup>h</sup>* = *V<sup>h</sup>* × *V<sup>h</sup>* and **V***h*0 = *Vh*0 × *Vh*0 .

For the pressure, for each node *Nl* ∈ *H*, we will match a function:

$$\Theta\_m(\mathbf{x}) = \rho^{\mu^\diamond}(\mathbf{x}) \cdot \theta\_m(\mathbf{x}), \ m = 0, 1, \dots, \ell$$

where *<sup>θ</sup>m*(**x**) ∈ *<sup>Y</sup>h*, *<sup>θ</sup>m*(*Nj*) = -1, *m* = *j*, 0, *m* = *j*, *m*, *j* = 0, 1, . . . ; *μ* is a parameter.

Then, we define a set *Q<sup>h</sup>*, for the pressure, such that for any *qh* ∈ *Q<sup>h</sup>*, we have:

$$q^h(\mathbf{x}) = \sum\_{\mathfrak{m}} \varepsilon\_{\mathfrak{m}} \Theta\_{\mathfrak{m}}(\mathbf{x}),\tag{14}$$

where *em* = *ρ*<sup>−</sup>*μ* (*Nm*)*e*˜*m*.

**Remark 1.** *The coefficients dj and ei in (13) and (14) are defined as a solution of a system (17) (see below).*

**Remark 2.** *The following embedding of sets is valid:*

$$\mathbf{V}^{\mathrm{h}} \subset \mathbf{W}^{1}\_{2,\nu}(\Omega\_{\mathrm{h}},\delta), \mathbf{V}^{\mathrm{h}}\_{0} \subset \mathbf{W}^{1}\_{2,\nu}(\Omega\_{\mathrm{h}},\delta), \ Q^{\mathrm{h}} \subset L^{0}\_{2,\nu}(\Omega\_{\mathrm{h}},\delta).$$

**Definition 2.** *The pair* (**u***hν*(**x**), *Phν* (**x**)) ∈ **V***<sup>h</sup>* × *Qh is called an approximate <sup>R</sup>ν-generalized solution for an Oseen system in the rotation form (8) and (9) obtained by the weighted FEM if the equalities:*

$$a(\mathbf{u}\_{\nu}^{h}, \mathbf{v}^{h}) + b(\mathbf{v}^{h}, P\_{\nu}^{h}) = l(\mathbf{v}^{h}),\tag{15}$$

$$\varepsilon(\mathbf{u}\_{\nu}^h \boldsymbol{q}^h) = 0 \tag{16}$$

*hold for any pair* (**v***<sup>h</sup>*(**x**), *q<sup>h</sup>*(**x**)) ∈ **V***h*0 × *Q<sup>h</sup>*, *where* **u***hν* = (*uh<sup>ν</sup>*,1, *<sup>u</sup>h<sup>ν</sup>*,<sup>2</sup>) *and ω* ∈ *<sup>L</sup>*2,*γ*(<sup>Ω</sup>, *<sup>δ</sup>*),**<sup>f</sup>** ∈ **<sup>L</sup>**2,*γ*(<sup>Ω</sup>, *<sup>δ</sup>*), **g** ∈ **W**1/2 2,*γ* (<sup>Γ</sup>, *<sup>δ</sup>*), *ν* ≥ *γ*.

Thus, we construct a weighted FEM to find an *<sup>R</sup>ν*-generalized solution for the rotation form of the Oseen problem (8) and (9).

Then, using (15) and (16), we ge<sup>t</sup> a system of linear algebraic equations:

$$\mathbf{A}\mathbf{d} + \mathbf{B}\mathbf{e} = \omega, \qquad \mathbf{C}^{\mathrm{T}}\mathbf{d} = \mathbf{z}, \tag{17}$$

where **d** = (*d*0, *d*2, *d*4,..., *d*1, *d*3, *<sup>d</sup>*5,...)*<sup>T</sup>*, **e** = (*<sup>e</sup>*0,*e*1,*e*2,...)*<sup>T</sup>*, *ω* = **<sup>F</sup>***h*, **z** = **0**.

### **4. Iterative Algorithm**

Now, we present an iterative procedure for solving the system of equations (17). Note that the system (17), which needs to be solved, has a large dimension, and moreover, its matrix is sparse. Finding the solution of the system by the direct method is not possible, so that we will construct a convergen<sup>t</sup> iterative process of the following type [19]:


where **A** ˆ is a preconditioning matrix to **A** and **S** ˆ is a preconditioning matrix to **S** = **<sup>C</sup>***T***A**−1**B**, which is called the Schur complement matrix. Next, we describe the process of constructing preconditioning matrices **A** ˆ and **S** ˆ .

*Symmetry* **2019**, *11*, 54

At first, we build a preconditioner **A**ˆ applying an incomplete **LU** factorization, where **L** and **U** are low unitriangular and upper triangular matrices respectively. At each iteration in Item 2, we employ the GMRES(*l*) method (see [20]) as the solution of a problem **Av** = **s** with the left preconditioner **A** ˆ . The method is designed so that it approximates the solution in an *l*th order Krylov subspace. In our research, the dimension of a Krylov subspace is equal to 10; so that if **r**0 = **A** ˆ <sup>−</sup><sup>1</sup>(**s** − **Av**), then the Arnoldi procedure will build an orthonormal basis of the subspace: Span{**<sup>r</sup>**0,(**A**<sup>ˆ</sup> <sup>−</sup>1**A**)<sup>1</sup>**r**0,...,(**A**<sup>ˆ</sup> <sup>−</sup>1**A**)<sup>9</sup>**r**0}.

Secondly, we build an intermediate matrix **S** ˜ to **S** ˆ . The matrix **S** ˜ represents a mass matrix **<sup>M</sup>***<sup>ν</sup>*,*μ*,*<sup>ν</sup> P* of a special view, such that on all elements *K* ∈ Υ*h* :

$$(\mathbf{M}\_P^{\overline{\eta}, \mu^{\diamond} \nu})\_{ij} = \frac{1}{\overline{\nu}} \int\_K \rho^{2(\nu + \mu^{\diamond})} \theta\_j(\mathbf{x}) \, \theta\_i(\mathbf{x}) d\mathbf{x}, \ \theta\_j(\mathbf{x}), \ \theta\_i(\mathbf{x}) \in Y^h, j, i = 0, 1, \dots, n$$

After that, we determine a matrix **S** ¯ , which is equal to a diagonal matrix **M** ¯ *<sup>ν</sup>*,*μ*,*<sup>ν</sup> P* with elements **M**¯ *<sup>ν</sup>*,*μ*,*<sup>ν</sup> P ii* = ∑*k* **M***<sup>ν</sup>*,*μ*,*<sup>ν</sup> P ik*. In other words, **S**¯ *ii* = ∑*k* **S**˜ *ik*. It is known (see [9] and the references therein)thatsuchdiagonallumping **S** ¯ isagoodpreconditionertotheinitialmatrix **S** ˜ .

Therefore, in order to determine the vector Ψ- := **S**ˆ <sup>−</sup>1*χ*, at each iteration of Item 3, we must find a solution to the following internal procedure: (1) *φ*0 = **0**; (2) *φm* = *φ<sup>m</sup>*−<sup>1</sup> + **S**¯ <sup>−</sup><sup>1</sup>(*χ* − **S**˜ *φ<sup>m</sup>*−<sup>1</sup>) (*m* = 1, . . . , *<sup>M</sup>*); (3) Ψ- = *φ<sup>M</sup>*.

We apply the GMRES(5) method, where Span{**r**¯,(**S**¯ <sup>−</sup>1**S**˜)<sup>1</sup>**r**¯, ... ,(**S**¯ −1**S**˜)<sup>4</sup> **<sup>r</sup>**¯}, and **r** ¯ = **S** ¯ <sup>−</sup><sup>1</sup>(*χ* − **S**˜ *φ<sup>m</sup>*−<sup>1</sup>).

### **5. Numerical Experiments**

Now, we present numerical results for the Oseen system in the rotation form (8) and (9) and show the advantage of the proposed method.

Let Ω*k* = (−*l*; *l*) × (−*l*; *l*) \ *G*¯ *k* be a polygon with one internal corner greater than 180◦ on Γ*k* whose vertex is at the origin. We will consider the following sizes of the reentrant corner: *ωk* = <sup>2</sup>*k*+<sup>1</sup> 2*k π*, *k* = 1, 2, 3. The triangulation Υ*h* (see Section 3) of each Ω ¯ *k*, *k* = 1, 2, 3 and *l* = 1 we present in Figure 2.

**Figure 2.** The triangulation Υ*h* of a domain Ω ¯ *k*.

In a test problem, we consider the solution of the problem (8), (9), which has a singularity in a neighborhood of a point located at the origin. Let *α* = *ν*¯ = 1, *w* = *b* · curl **u**, *b* = 0.95, and for each corner *ωk* in polar coordinates (*r*, *ϕ*), we have an auxiliary function:

$$\Psi\_k(\varphi) = \frac{\sin((1+\lambda\_k)\varphi)\cos(\lambda\_k\omega\_k)}{1+\lambda\_k} - \frac{\sin((1-\lambda\_k)\varphi)\cos(\lambda\_k\omega\_k)}{1-\lambda\_k} + \cos((1-\lambda\_k)\varphi) - \cos((1+\lambda\_k)\varphi).$$

Then, the exact solution **u** = (*<sup>u</sup>*1, *<sup>u</sup>*2) and *P* of the problem (8) and (9) for each corner *ωk*, *k* = 1, 2, 3, in polar coordinates has the following form:

$$\mu\_1(r,\,\,\boldsymbol{\varrho}) = r^{\lambda\_k} \cdot ( (\lambda\_k + 1) \sin \boldsymbol{\varrho} \cdot \boldsymbol{\Psi}\_k(\boldsymbol{\varrho}) + \cos \boldsymbol{\varrho} \cdot \boldsymbol{\Psi}\_k'(\boldsymbol{\varrho}) ),$$

$$\mu\_2(r,\,\,\boldsymbol{\varrho}) = r^{\lambda\_k} \cdot (\sin \boldsymbol{\varrho} \cdot \boldsymbol{\Psi}\_k'(\boldsymbol{\varrho}) - (\lambda\_k + 1) \cos \boldsymbol{\varrho} \cdot \boldsymbol{\Psi}\_k(\boldsymbol{\varrho})),$$

$$P(r,\,\,\,\boldsymbol{\varrho}) = r^{\lambda\_k - 1} \cdot \frac{(\lambda\_k + 1)^2 \boldsymbol{\Psi}\_k'(\boldsymbol{\varrho}) + \boldsymbol{\Psi}\_k^{\prime\prime\prime}(\boldsymbol{\varrho})}{\lambda\_k - 1},$$

where *λk* = min{*λ* : sin(*λ <sup>ω</sup>k*) + *λ* sin *ωk* = 0 and *λ* > <sup>0</sup>}.

Thus, for the corner *ω*1 = 3*π* 2 , we have *λ*1 ≈ 0.544483, for *ω*2 = 5*π* 4 , *λ*2 ≈ 0.673583, and for *ω*3 = 9*π* 8 , *λ*3 ≈ 0.800766. The proposed solution is analytical in Ω¯ *k* \ (0, <sup>0</sup>), but unfortunately, *P* ∈ *W*<sup>1</sup> 2 (<sup>Ω</sup>*k*), **u** ∈ **W**<sup>2</sup> <sup>2</sup>(<sup>Ω</sup>*k*).

In numerical experiments, we use meshes with a various step size *h* and number *N*, where *N* · *h* equals two. The approximate generalized solution (velocity field) by classical FEM converges to the exact one in the **W**<sup>1</sup> 2(<sup>Ω</sup>*k*) norm with a rate depending on the size of reentrant corner *ω*, the so-called pollution effect (see [12] and the references therein): for a corner *ω*1 = 3*π* 2 , we have the rate of convergence, which is equal to <sup>O</sup>(*h*0.55), for a corner *ω*2 = 5*π* 4 , <sup>O</sup>(*h*0.67), and for a corner *ω*3 = 9*π* 8 , O(*h*0.8) (see Table 1); whereas, the approximate *<sup>R</sup>ν*-generalized solution by the presented weighted FEM converges to the exact one in the **<sup>W</sup>**12,*ν*(<sup>Ω</sup>*k*) norm with a rate that is independent of the value of the internal angle *ω* and has the first order by *h* (see Table 2), where we derive computationally the optimal parameters *δ*, *ν* = *<sup>ν</sup>op<sup>t</sup>* and *ν* . Both errors for the *<sup>R</sup>ν*-generalized and generalized solutions visually are represented in Figure 3 for different values of a number *N*.

**Table 1.** The generalized solution error (**u***<sup>h</sup>* − **u**) in the norm of a space **W**<sup>1</sup> <sup>2</sup>(<sup>Ω</sup>*k*).


**Table 2.** The *<sup>R</sup>ν*-generalized solution error (**u***<sup>h</sup> ν* − **<sup>u</sup>***ν*) in the norm of a space **<sup>W</sup>**12,*ν*(<sup>Ω</sup>*k*), where *ν* = *μ* = *λk* − 1 and *ν* = *μ* = *<sup>ν</sup>opt*.


Let *δji* = |*uj*(*Mi*) − *uh j* (*Mi*)| and *δji* = |*uj*(*Mi*) − *uh <sup>ν</sup>*,*j* (*Mi*)|, *j* = 1, 2, *Mi* ∈ *G* Ω be errors for the generalized and *<sup>R</sup>ν*-generalized solutions, respectively. Then, we show the percentage of nodes, where *δ*1*i* and *δ*1*i* are less than a given value ¯ *l*. The quantity of points *Mi* ∈ *G*Ω, where *δ*1*j* < ¯ *l* (for the classical FEM), is significantly less in relation to the quantity of points *Mi* ∈ *G*Ω, where *δ*1*j* < ¯ *l* (for the proposed weighted FEM) for all sizes of the reentrant corner *ω* (see Table 3). Moreover, in numerical experiments, the number of nodes *Mi*, where *δ*2*i* < ¯ *l* and *δ*2*i* < ¯ *l* are approximately equal to the number of nodes *Mi*, where *δ*1*i*< ¯ *l* and *δ*1*i* < ¯ *l*, *l* = 1, 2, respectively.


**Table 3.** The percentage of points *Mi* ∈ *G*Ω, where the values *δ*1*i* and *δ*1*i*are less than ¯ *l*, *l* = 1, 2.

Then, we present the distribution of errors *δji* and *δji* in the points *Mk* for components *uhν*,*j* and *uhj* for all sizes *ωl*, *l* = 1, 2, 3, *j* = 1, 2 , and *h*, such that *N* = 148 and *N* = 296. The weighted finite element method allows us to perform computations with high accuracy both inside of the domain and near the point of singularity. Moreover, the error of the proposed FEM is localized near the point of singularity and does not extend into the interior of the domain, in contrast to the error of the classical FEM for all values of the internal corner *ω* (see Figures 4–15).

In Figures 16–18, we show the dependence of error in the **<sup>W</sup>**12,*ν*(<sup>Ω</sup>*k*) norm on the parameter *ν* (*μ* = *ν*), where each minimum is compatible with the best value *<sup>ν</sup>opt*. Any value from the interval (*<sup>λ</sup>k* − 1, 0) can be taken as an exponent *ν* for the presented FEM in the domain Ω*k* with a reentrant corner *ωk*. Moreover, if the exponent *μ* does not coincide with *<sup>ν</sup>*, then we ge<sup>t</sup> substantially worse results. This research was supported in through computational research provided by the Shared Facility Center "Data Center of FEB RAS".

**Figure 3.** The errors of (**left**) a classical FEM in the **W**12 norm and (**right**) a weighted FEM in the **<sup>W</sup>**12,*ν* norm, where *ν* = 1.6, *δ* = 0.01375 : *ω*1 = 3*π*2 , *ν* = *<sup>ν</sup>op<sup>t</sup>* = −0.35; *ω*2 = 5*π*4 , *ν* = *<sup>ν</sup>op<sup>t</sup>* = −0.25; *ω*3 = 9*π*8 , *ν* = *<sup>ν</sup>op<sup>t</sup>* = −0.125, for different values of a number *N*.

**Figure 4.** The errors *δ*1*i* of the approximate generalized solution (*uh*1) : *ω*1 = 3*π*2 , (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 5.** The errors *δ*1*i* of the approximate *<sup>R</sup>ν*-generalized solution (*uh<sup>ν</sup>*,<sup>1</sup>) : *ω*1 = 3*π*2 , *ν* = 1.6, *δ* = 0.01375, *ν* = *μ* = −0.35, (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 6.** The distribution of the errors *δ*2*i* of the approximate generalized solution (*uh*2) : *ω*1 = 3*π*2 , (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 7.** The errors *δ*2*i* of the approximate *<sup>R</sup>ν*-generalized solution (*uh<sup>ν</sup>*,<sup>2</sup>) : *ω*1 = 3*π*2 , *ν* = 1.6, *δ* = 0.01375, *ν* = *μ* = −0.35, (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 8.** The errors *δ*1*i* of the approximate generalized solution (*uh*1) : *ω*2 = 5*π*4 , (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 9.** The errors *δ*1*i* of the approximate *<sup>R</sup>ν*-generalized solution (*uh<sup>ν</sup>*,<sup>1</sup>) : *ω*2 = 5*π*4 , *ν* = 1.6, *δ* = 0.01375, *ν* = *μ* = −0.25, (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 10.** The errors *δ*2*i* of the approximate generalized solution (*uh*2) : *ω*2 = 5*π*4 , (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 11.** The errors *δ*2*i* of the approximate *<sup>R</sup>ν*-generalized solution (*uh<sup>ν</sup>*,<sup>2</sup>) : *ω*2 = 5*π*4 , *ν* = 1.6, *δ* = 0.01375, *ν* = *μ* = −0.25, (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 12.** The errors *δ*1*i* of the approximate generalized solution (*uh*1) : *ω*3 = 9*π*8 , (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 13.** The errors *δ*1*i* of the approximate *<sup>R</sup>ν*-generalized solution (*uh<sup>ν</sup>*,<sup>1</sup>) : *ω*3 = 9*π*8 , *ν* = 1.6, *δ* = 0.01375, *ν* = *μ* = −0.125, (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 14.** The errors *δ*2*i* of the approximate generalized solution (*uh*2) : *ω*3 = 9*π*8 , (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 15.** The errors *δ*2*i* of the approximate *<sup>R</sup>ν*-generalized solution (*uh<sup>ν</sup>*,<sup>2</sup>) : *ω*3 = 9*π*8 , *ν* = 1.6, *δ* = 0.01375, *ν* = *μ* = −0.125, (**left**) *N* = 148, (**right**) *N* = 296.

**Figure 16.** The dependence of error (**u***hν* − **u**) in the **<sup>W</sup>**12,*ν*(<sup>Ω</sup>1) norm on the degree *<sup>ν</sup>*, *ω*1 = 3*π*2.

**Figure 17.** The dependence of error (**u***hν* − **u**) in the **<sup>W</sup>**12,*ν*(<sup>Ω</sup>2) norm on the degree *<sup>ν</sup>*, *ω*2 = 5*π*4 .

**Figure 18.** The dependence of error (**u***hν* − **u**) in the **<sup>W</sup>**12,*ν*(<sup>Ω</sup>3) norm on the degree *<sup>ν</sup>*, *ω*3 = 9*π*8.
