New Numerical Method for the Rotation form of the Oseen Problem with Corner Singularity

: In the paper, a new numerical approach for the rotation form of the Oseen system in a polygon Ω with an internal corner ω greater than 180 ◦ on its boundary is presented. The results of computational simulations have shown that the convergence rate of the approximate solution (velocity ﬁeld) by weighted FEM to the exact solution does not depend on the value of the internal corner ω and equals O ( h ) in the norm of a space W 12, ν ( Ω ) .


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 1 2 (Ω). In short, the Dirichlet integral from the solution diverges. In the case when the solution belongs to the space W 1 2 (Ω), but it does not belong to the W 2 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 R ν -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][4][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 2 2 (Ω) and W 1 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 R ν -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 W 1 2,ν (Ω 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 R ν -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.

R ν -Generalized Solution of the Oseen Problem
Let x = (x 1 , x 2 ) be an element of the Euclidean space R 2 , where x = x 2 1 + x 2 2 1/2 and dx = 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 u(x, t) = (u 1 (x, t), u 2 (x, t)) and a kinematic pressure p(x, t) from: 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: where g = (g 1 , g 2 ) is given vector function on Γ and u 0 (x) = (u 0 1 (x), u 0 2 (x)) -in Ω. We introduce the following notation: We have a formal equality: If u = v in (3), then we have a relation: Let P = p + 1 2 u 2 , using (4), for vector function u; we get the rotation form of the Navier-Stokes system for an incompressible flow: 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: 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: 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 R ν -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: be the m th order generalized derivatives of a function v(x) in Ω, For the function v(x), we define the following inequalities: where α > 0 and constant C 2 > 0 do not depend on m and α. Denote by L 2,α (Ω) a space of functions v(x), such that: If w = (w 1 , w 2 ) is a vector function, then we define the weighted vector function space (10) and (11) : If w = (w 1 , w 2 ) is a vector function, then we denote by W 1 2,α (Ω) the weighted vector function space with a norm w W 1 2,α (Ω) = w 1 , that meet the conditions (10) and (11) with a bounded W 1 2,α (Ω) norm. We denote by ) a closure, with respect to the W 1 2,α (Ω) norm, of the set of infinitely-differentiable functions with compact support in Ω that meet the inequalities (10) and (11). Then, we denote by (8) and (9) meet the following conditions: Bilinear and linear forms are as follows: Oseen system in the rotation form (8) and (9) such that for all where functions w, f and g satisfy the conditions (12) and ν ≥ γ.
Note that the bilinear and linear forms in the definition of an R ν -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 W 2 2,ν (Ω, δ) and W 1 2,ν (Ω, δ), respectively. This property of the R ν -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.

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 R ν -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 L i , which we call macro-elements. The set of elements L i represents a quasi-uniform (see [1]) triangulation T h of Ω. Then, we divide each macro-element L i ∈ T h into three triangles K i j using the barycenter of L i . Thus, we construct a triangulation Υ h , which is based on a barycenter refinement of a triangulation T h . Denote by Ω h the set of resulting triangles (which are called finite elements) with sides of order h, i.e., Ω h = Let A m and B l be vertices and midpoints of the finite element sides K i j ∈ Υ h , respectively. Then, for the components of a velocity field and pressure, we define sets of nodes G and H, respectively, where C k coincide with a node A m on the appropriate element K i j ∈ Υ h (see Figure 1). Now, we define spaces of the SV element pair. The space X h , for the components of the velocity field, coincides with the corresponding space of degree two of the THelement pair, i.e.,  The SV element pair has an important property, namely div X h ⊂ Y h . This means that there exists a function y h ∈ Y h equal to div w h such that: from the condition for performing mass conservation in a weak sense, i.e., Ω div w h ψ h dx = 0 ∀ψ h ∈ Y h , we get a pointwise mass conservation, i.e., div w h L 2 (Ω) = 0. Moreover, in [18], it was established that spaces of the SV element pair before us satisfy the Ladyzhenskaya-Babuska-Brezzi condition. Note, that approximations obtained using the 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 M k ∈ G Ω , we will match a function: We define a set V h , for components of the velocity field, such that for any velocity field where For the pressure, for each node N l ∈ H, we will match a function: Then, we define a set Q h , for the pressure, such that for any q h ∈ Q h , we have: where e m = ρ −µ (N m )ẽ m . (13) and (14) are defined as a solution of a system (17) (see below).
Thus, we construct a weighted FEM to find an R ν -generalized solution for the rotation form of the Oseen problem (8) and (9).
Then, using (15) and (16), we get a system of linear algebraic equations:

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 convergent iterative process of the following type [19]: (1) Let (d 0 , e 0 ) be an initial guess for the system (17). We iterate (n = 0, 1, 2, . . .) until the stopping condition is fulfilled; whereÂ is a preconditioning matrix to A andŜ is a preconditioning matrix to S = C T A −1 B, which is called the Schur complement matrix. Next, we describe the process of constructing preconditioning matricesÂ andŜ.
At first, we build a preconditionerÂ 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Â. 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 =Â −1 (s − Av), then the Arnoldi procedure will build an orthonormal basis of the subspace: Span{r 0 , (Â −1 A) 1 r 0 , . . . , (Â −1 A) 9 r 0 }.
Secondly, we build an intermediate matrixS toŜ. The matrixS represents a mass matrix M ν,µ ,ν P of a special view, such that on all elements K ∈ Υ h : After that, we determine a matrixS, which is equal to a diagonal matrixM It is known (see [9] and the references therein) that such diagonal lumpingS is a good preconditioner to the initial matrixS.

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) \Ḡ 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 = 2 k +1 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. 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: Then, the exact solution u = (u 1 , u 2 ) and P of the problem (8) and (9) for each corner ω k , k = 1, 2, 3, in polar coordinates has the following form: where λ k = min{λ : sin(λ ω k ) + λ sin ω k = 0 and λ > 0}. 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, 0), but unfortunately, P ∈ W 1 2 (Ω k ), u ∈ W 2 2 (Ω 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 1 2 (Ω 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 O(h 0.55 ), for a corner ω 2 = 5π 4 , O(h 0.67 ), and for a corner ω 3 = 9π 8 , O(h 0.8 ) (see Table 1); whereas, the approximate R ν -generalized solution by the presented weighted FEM converges to the exact one in the W 1 2,ν (Ω 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 δ, ν = ν opt and ν . Both errors for the R ν -generalized and generalized solutions visually are represented in Figure 3 for different values of a number N.  Table 2. The R ν -generalized solution error (u h ν − u ν ) in the norm of a space W 1 2,ν (Ω k ), where ν = µ = λ k − 1 and ν = µ = ν opt . Let δ ji = |u j (M i ) − u h j (M i )| and δ ji = |u j (M i ) − u h ν,j (M i )|, j = 1, 2, M i ∈ G Ω be errors for the generalized and R ν -generalized solutions, respectively. Then, we show the percentage of nodes, where δ 1i and δ 1i are less than a given value¯ l . The quantity of points M i ∈ G Ω , where δ 1j <¯ l (for the classical FEM), is significantly less in relation to the quantity of points M i ∈ G Ω , where δ 1j <¯ l (for the proposed weighted FEM) for all sizes of the reentrant corner ω (see Table 3). Moreover, in numerical experiments, the number of nodes M i , where δ 2i <¯ l and δ 2i <¯ l are approximately equal to the number of nodes M i , where δ 1i <¯ l and δ 1i <¯ l , l = 1, 2, respectively. Then, we present the distribution of errors δ ji and δ ji in the points M k for components u h ν,j and u h j 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 . In Figures 16-18, we show the dependence of error in the W 1 2,ν (Ω k ) norm on the parameter ν (µ = ν ), where each minimum is compatible with the best value ν opt . Any value from the interval (λ 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 ν , then we get substantially worse results. This research was supported in through computational research provided by the Shared Facility Center "Data Center of FEB RAS".

Conclusions
The main results of the numerical experiments for the Oseen problem (8) and (9) lead to the following conclusions:

•
The approximate generalized solution (velocity field) by classical FEM converges to the exact one in the W 1 2 (Ω k ) norm with a rate O(h λ ), λ < 1, where the exponent λ depends on the size of reentrant corner ω, the so-called pollution effect (see [12] and the references therein), while the approximate R ν -generalized solution by the presented weighted FEM converges to the exact one in the W 1 2,ν norm with a rate that is independent of the value of the internal angle ω and has the first order by h for various values of ν, δ (see Tables 1 and 2 and Figure 3). • Thanks to Theorem 3.1 in [13], there exists a limitation on the radius δ * of the neighborhood of a reentrant corner ω and ρ(x) exponent ν * in Definition 1, that for all δ < δ * and ν > ν * , a weighted inf-sup condition holds. After a series of computational experiments, we conclude that ν * ∼ 1 and δ * ∼ h.

•
The proposed approach allows us to compute the approximate solution by the weighted FEM with a given accuracy 10 −3 , for example in a case when the internal corner ω is equal to 3π 2 , about 10 6 -times faster than using classical FEM. Note that in implementing the weighted FEM, one can spend about 10 6 -times less computing resources and energy consumption. • The weighted finite element method enables us to perform computations with high accuracy, both inside of the domain and near the point of singularity.
Author Contributions: V.A.R. and A.V.R. contributed equally in each stage of the work. All authors read and approved the final version of the paper.

Funding:
The reported study was supported by RFBR and RSF according to the research project Nos. 19-01-00007-a and 19-71-20006, respectively.