Numerical approach for the one stationary nonlinear
problem governing the flow of incompressible viscous
fluid in 𝐿-shaped domain
Alexey V. Rukavishnikova , Viktor A. Rukavishnikovb
a
  Institute of Applied Mathematics of Far-Eastern Branch, Russian Academy of Sciences, Dzerzhinsky Str., 54, Khabarovsk,
680000, Russia
b
  Computing Center of the Far Eastern Branch of the Russian Academy of Sciences, Kim Yu Chen Str., 65, Khabarovsk, 680000,
Russia


                                          Abstract
                                          The steady Navier-Stokes equations governing the flow of an incompressible viscous fluid in the rotation form
                                          in 𝐿-shaped domain is considered. The weighted finite element method based on the definition of an 𝑅𝜈 -
                                          generalized solution is constructed. The advantage of the proposed approach over classical approximations
                                          is numerically established. The modern elements of computational technologies to find the optimal parameters
                                          of the proposed method are used.

                                          Keywords
                                          corner singularity, weighted finite element method, preconditioning, steady Navier-Stokes equations




1. Introduction
Most of the mathematical models representing natural processes are described using boundary value
problems for the systems of partial differential equations with a singularity. The peculiarity of the
solution is as follows systems in a bounded, connected domain of the Euclidean space 𝐑2 can be
attributed to the presence of obtuse corners on its boundary, to the degeneration of initial data, or
to internal characteristics of the solution. If the solution of the boundary value problem does not
belong to the Sobolev space 𝑊21 (Ω), then it is called strong singular. If the solution of the boundary
value problem belongs to 𝑊21 (Ω), but does not belong to 𝑊22 (Ω), then the boundary value problem is
called weakly singular. The generalized solution of such problems in Ω with a boundary containing a
reentrant corner 𝜎 ∈ (𝜋, 2𝜋] belongs to the space 𝑊21+𝛼−𝜖 (Ω), 𝛼 < 1. Moreover, the approximate finite
element or finite difference solution by classical method converges to the exact one with a (ℎ𝛼 ) rate.
   In [1], it was proposed to define the solution of elliptic boundary value problems with a singularity
as an 𝑅𝜈 -generalized one. The approach allows us to introduce a weight space or a set, depending
on the geometry of the domain and input data (right-hand sides, equation coefficients, boundary
and initial data) to which an 𝑅𝜈 -generalized solution belongs. In [2, 3, 4], the existence, uniqueness
and differential properties of the elliptic problems solution are proved. In [5], a weighted analogue
of the Ladyzhenskaya-Babuska-Brezzi condition for the Stokes problem is established. In [6, 7, 8, 9,
10], a weighted finite element method (FEM) for an approximate solution of elliptic problems with a
singularity has been developed.

Far Eastern Workshop on Computational Technologies and Intelligent Systems, March 2–3, 2021, Khabarovsk, Russia
" 78321a@mail.ru (A.V. Rukavishnikov); vark0102@mail.ru (V.A. Rukavishnikov)
 0000-0002-7585-4559 (A.V. Rukavishnikov); 0000-0002-3702-1126 (V.A. Rukavishnikov)
                                       © 2021 Copyright for this paper by its authors.
                                       Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
    CEUR
    Workshop
    Proceedings
                  http://ceur-ws.org
                  ISSN 1613-0073       CEUR Workshop Proceedings (CEUR-WS.org)
   In the paper, an 𝑅𝜈 -generalized solution of the steady Navier-Stokes equations governing the flow of
an incompressible viscous fluid in the rotation form in 𝐿-shaped domain is defined. We use Picard’s
iterative procedure [11] to find a solution of a nonlinear problem. Then, we construct a weighted
finite element scheme based on the definition of an 𝑅𝜈 -generalized solution: 1) the functions of the
finite element spaces satisfy the mass conservation law in a strong sense — Scott-Vogelius (SV) pair
of 2nd order [12]; 2) basis functions are the product of SV functions and weight functions in some
degree. This construction allows us better take into account the behavior of the solution in the 𝛿-
neighborhood of the singularity point and increase the convergence rate of the approximate solution
to the exact one to the first order with respect to the grid step ℎ, i. e. (ℎ) rate in 𝑊2,𝜈1 (Ω) norm.

Thus, we develop the numerical method overcomes the so-called pollution effect (see [13]). The same
advantage for other hydrodynamics problems was achieved in [14, 15, 16]. The optimal values 𝜈, 𝜈 ∗
and 𝛿 of the presented weighted FEM using the modern elements of computational technologies were
derived numerically.
   The paper consists of six sections. Section 2 is devoted to the definition of an 𝑅𝜈 -generalized solu-
tion. In Section 3, we present the weighted FEM. The iterative procedure for solving the systems of
linear algebraic equations is constructed in Section 4. In Section 5, we shaw and discuss the results of
computational experiments. Necessary conclusions are made in last section.


2. The problem statement
Let Ω be a bounded, connected domain in the Euclidean space 𝐑2 . Denote by Ω̄ and Γ the closure and
boundary
      √ of Ω, respectively, Ω = Ω ∪ Γ. Let 𝐱 = (𝑥1 , 𝑥2 ) be an element of 𝐑 , where 𝑑𝐱 = 𝑑𝑥1 𝑑𝑥2 and
                                 ̄                                               2

‖𝐱‖ = 𝑥1 + 𝑥2 are the measure and norm of 𝐱, respectively.
         2   2

   We write the steady Navier-Stokes equations governing the flow of an incompressible viscous fluid
in the convection form: find a velocity field 𝐮 = 𝐮(𝐱) = (𝑢1 (𝐱), 𝑢2 (𝐱)) and a kinematic pressure 𝑝 = 𝑝(𝐱)
from
                 − 𝜈̄ △ 𝐮 + (𝐮 ⋅ ∇)𝐮 + 𝛼𝐮 + ∇𝑝 = 𝐟,              div 𝐮 = 𝟎        in     Ω,             (1)
                                                                     𝐮=𝐠         on      Γ,               (2)
where 𝜈̄ > 0 is the kinematic viscosity coefficient (inversely proportional to the Reynolds number),
𝛼 > 0, 𝐟 = 𝐟(𝐱) = (𝑓1 (𝐱), 𝑓2 (𝐱)) and 𝐠 = 𝐠(𝐱) = (𝑔1 (𝐱), 𝑔2 (𝐱)) are given force field in Ω and boundary
data on Γ, respectively. Denote by △, ∇ and div the Laplace, gradient and divergence operators in 𝐑2 ,
respectively.
  Further, we introduce the necessary notation. Let 𝐰 = (𝑤1 , 𝑤2 ), 𝐯 = (𝑣1 , 𝑣2 ) and 𝑏 – scalar, then
                                                                                     𝜕𝑤1 𝜕𝑤2
            𝐰 ⋅ 𝐯 = 𝑤1 𝑣1 + 𝑤2 𝑣2 ,      𝑏 × 𝐰 = (−𝑏𝑤2 , 𝑏𝑤1 )𝑇 ,      curl 𝐰 = −       +    .
                                                                                     𝜕𝑥2 𝜕𝑥1
  We have the identity
                                                                 1
                                     (𝐮 ⋅ ∇)𝐮 = ( curl 𝐮) × 𝐮 + ∇𝐮2 .                                     (3)
                                                                 2
  It follows from the equality (𝐰 ⋅ ∇)𝐮 + (𝐮 ⋅ ∇)𝐰 = ∇(𝐰 ⋅ 𝐮) + ( curl 𝐰) × 𝐮 + ( curl 𝐮) × 𝐰 and assumption
that 𝐰 = 𝐮.
  Using (3), with 𝑃 = 𝑝 + 12 𝐮2 for the system (1), (2) we get the rotation form of the steady Navier-
Stokes equations: find a velocity field 𝐮 and a Bernoulli pressure 𝑃 such that

                  − 𝜈̄ △ 𝐮 + ( curl 𝐮) × 𝐮 + 𝛼𝐮 + ∇𝑃 = 𝐟,        div 𝐮 = 0      in       Ω,               (4)
                                                                                            𝐮=𝐠            on          Γ.        (5)
   The system (4), (5) as well as (1), (2) is nonlinear due to the presence of the rotation term ( curl 𝐮)×𝐮
in the momentum equations. The system on (4), (5) and term in particular we linearized by Picard’s
procedure (see [11]).
   Starting with an initial approximation 𝐮(0) for which

                          div 𝐮(0) = 0             in          Ω          and         𝐮(0) = 𝐠          on         Γ             (6)

Picard’s iteration constructs a sequence of solutions (𝐮(𝑘) , 𝑃 (𝑘) ) by solving the linear system:

             − 𝜈̄ △ 𝐮(𝑘) + ( curl 𝐮(𝑘−1) ) × 𝐮(𝑘) + 𝛼𝐮(𝑘) + ∇𝑃 (𝑘) = 𝐟,                   div 𝐮(𝑘) = 0            in        Ω,   (7)

                                                                                               𝐮(𝑘) = 𝐠           on        Γ.   (8)
    Note that the initial Bernoulli pressure in (7) need not be specified. If 𝜈̄ be a not too small and 𝐟 be a
not too large, the steady Navier-Stokes equations (4), (5) have a unique solution (𝐮, 𝑃) and the iterates
(𝐮(𝑘) , 𝑃 (𝑘) ), 𝑘 = 1, 2, in (7), (8) converge to it as 𝑘 → ∞ for any choice of the arbitrary 𝐮(0) satisfying
(6) (see [11]).
    Note that for a linearized system (7), (8) the conservation laws of the mass and momentum remain
in force.
    In the article, we consider the special case of a bounded polygon domain Ω. Let Ω be a 𝐿-shaped
domain with one reentrant obtuse corner equals to 3𝜋2 on the boundary and its vertex coincides with
the origin. We define an 𝑅𝜈 -generalized solution in each Picard’s iteration of the problem (7), (8) and
construct the effective weighted FEM. Thus, we solve the nonlinear problem (4), (5) governing the
flow of a incompressible viscous fluid in the rotation form and show the advantage of our approximate
method over the classical approaches in a 𝐿-shaped domain by the computational simulations.
    Let us introduce the notation and define necessary spaces of generalized functions. Denote by
                            ≤ 𝛿 < 1, 𝛿 > 0} a part of a 𝛿-neighborhood of a point (0, 0) contained in Ω̄ . Let
  ′
Ω𝛿 = {𝐱           ̄
          { ∈ Ω ∶ ‖𝐱‖    ′
             ‖𝐱‖, 𝐱 ∈ Ω𝛿 ,
𝜌(𝐱) =                       ′     be a weight function.
               𝛿 , 𝐱 ∈ Ω̄ ⧵ Ω𝛿
    Denote by 𝐿2,𝛽 (Ω) and 𝑊2,𝛽      1 (Ω) the spaces of functions 𝑣(𝐱) with a bounded norms

                                                        √
                                                   ‖𝑣‖𝐿2,𝛽 (Ω) =         2𝛽   2
                                                                      ∫ 𝜌 (𝐱)𝑣 (𝐱)𝑑𝐱
                                                                      Ω

and                                                 √
                               ‖𝑣‖𝑊2,𝛽
                                   1 (Ω) =              ‖𝜌 𝛽 (𝐱)|𝑣(𝐱)|‖2𝐿2 (Ω) + ‖𝜌 𝛽 (𝐱)|𝐷 1 𝑣(𝐱)|‖2𝐿2 (Ω) ,
                                             |𝑚|
respectively, where 𝐷 𝑚 𝑣(𝐱) = 𝜕𝑥 𝑚𝜕 1 𝜕𝑥𝑣 𝑚2 , |𝑚| = 𝑚1 + 𝑚2 , 𝑚𝑖 ≥ 0 - integer.
                                         1          2
  Let 𝑊2,𝛽
       1 (Ω, 𝛿) for 𝛽 > 0 be a set of functions from the space 𝑊 1 (Ω), meets the conditions
                                                                2,𝛽

                                                                                              𝛿 𝛽+𝑚                    ′
                   ∫       𝜌 2𝛽 (𝐱)𝑣 2 (𝐱)𝑑𝐱 ≥ 𝐶1 > 0,                 |𝐷 𝑚 𝑣(𝐱)| ≤ 𝐶2 (                        𝐱 ∈ Ω𝛿 ,         (9)
                                                                                             𝜌(𝐱) )
                      ′
                 Ω̄ ⧵Ω𝛿

where 𝑚 = 0, 1 and 𝐶2 a positive constant which is not depend on 𝑚, with the norm of a space
   1 (Ω). Denote by 𝐿 (Ω, 𝛿) a set of functions from the space 𝐿 (Ω) which subject to conditions (9)
𝑊2,𝛽                 2,𝛽                                             2,𝛽
(only for 𝑚 = 0) with a norm of a space 𝐿2,𝛽 (Ω). Let 𝐿02,𝛽 (Ω, 𝛿) = {𝑞 ∈ 𝐿2,𝛽 (Ω, 𝛿) ∶ ∫ 𝜌 𝛽 𝑞𝑑𝐱 = 0}.
                                                                                                                   Ω
         𝑜              𝑜
   Let 𝑊2,𝛽
          1 (Ω, 𝛿) (𝑊 1 (Ω, 𝛿) ⊂ 𝑊 1 (Ω, 𝛿)) be a closure by 𝑊 1 (Ω) norm of a set of the infinitely-
                          2,𝛽            2,𝛽                      2,𝛽
differentiable functions with a compact support in Ω comply with the conditions (9). We will say
𝜑(𝐱) ∈ 𝑊2,𝛽 1/2
                (Γ, 𝛿), if exists a function Φ(𝐱) ∈ 𝑊2,𝛽
                                                     1 (Ω, 𝛿) such that Φ(𝐱)| = 𝜑(𝐱) and ‖𝜑‖
                                                                             Γ                1/2
                                                                                             𝑊2,𝛽 (Γ,𝛿) =
 inf ‖Φ‖𝑊2,𝛽
           1 (Ω) .
Φ|Γ =𝜑
   For the vector field 𝐯 = (𝑣1 , 𝑣2 ) we define sets 𝐋2,𝛽 (Ω, 𝛿) and 𝐖12,𝛽 (Ω, 𝛿) such that 𝑣𝑖 ∈ 𝐿2,𝛽 (Ω, 𝛿) and
                                                                         √
        1 (Ω, 𝛿), respectively, with a bounded norms ‖𝐯‖
                                                               𝐋2,𝛽 (Ω) = ‖𝑣1 ‖𝐿2,𝛽 (Ω) + ‖𝑣2 ‖𝐿2,𝛽 (Ω) for the first set
𝑣𝑖 ∈ 𝑊2,𝛽                                                                      2               2
                   √
and ‖𝐯‖𝐖12,𝛽 (Ω) = ‖𝑣1 ‖2𝑊 1 (Ω) + ‖𝑣2 ‖2𝑊 1 (Ω) for the second one. Similarly, we define the sets of vector
                             2,𝛽            2,𝛽
                        𝑜
fields 𝐖1/2
        𝛽   (Γ, 𝛿) and 𝐖 1 (Ω, 𝛿) on Γ and in Ω, respectively.
                         2,𝛽
   We introduce the concept of an 𝑅𝜈 -generalized solution for the linearized problem (7), (8).
                           𝜈 ∈ 𝐖2,𝜈 (Ω, 𝛿) and 𝑃𝜈 ∈ 𝐿2,𝜈 (Ω, 𝛿) is called an 𝑅𝜈 -generalized solution
   Definition 1. The pair 𝐮(𝑘)      1              (𝑘)    0
                                                                                                   𝑜
                                 𝜈 satisfies a condition (8) on Γ for any pair 𝐯 ∈𝐖2,𝜈 (Ω, 𝛿) and 𝑞 ∈
of the problem (7), (8), where 𝐮(𝑘)                                                1
 0
𝐿2,𝜈 (Ω, 𝛿)
                            𝑎𝑘 (𝐮(𝑘)            (𝑘)
                                 𝜈 , 𝐯) + 𝑏(𝐯, 𝑃𝜈 ) = 𝑙(𝐯),   𝑐(𝐮(𝑘)
                                                                 𝜈 , 𝑞) = 0

hold, where bilinear and linear forms are as follows

                  𝜈 , 𝐯) = ∫ [𝛼𝜌 𝐮𝜈 ⋅ 𝐯 + 𝜈̄ ∇𝐮𝜈 ⋅ ∇(𝜌 𝐯) + 𝜌 (( curl 𝐮𝜈
             𝑎𝑘 (𝐮(𝑘)           2𝜈 (𝑘)         (𝑘)    2𝜈     2𝜈        (𝑘−1)
                                                                             ) × 𝐮(𝑘)
                                                                                  𝜈 ) ⋅ 𝐯]𝑑𝐱,
                            Ω


 𝑏(𝐯, 𝑃𝜈(𝑘) ) = − ∫ 𝑃𝜈(𝑘) div (𝜌 2𝜈 𝐯)𝑑𝐱,            𝜈 , 𝑞) = − ∫ (𝜌 𝑞) div 𝐮𝜈 𝑑𝐱,
                                                  𝑐(𝐮(𝑘)            2𝜈       (𝑘)
                                                                                                𝑙(𝐯) = ∫ 𝜌 2𝜈 𝐟 ⋅ 𝐯𝑑𝐱
                Ω                                              Ω                                       Ω

and 𝐟 ∈ 𝐋2,𝛽 (Ω, 𝛿), 𝐠 ∈ 𝐖1/2
                          2,𝛽 (Γ, 𝛿), 𝜈 ≥ 𝛽 ≥ 0.


3. The weighted finite element method
Perform triangulation Υℎ based on the barycentric partition of the elements 𝐿𝑖 of the quasi-uniform
triangulation 𝑇ℎ of the domain Ω. Then, we divide each element 𝐿𝑖 ∈ 𝑇ℎ (macroelement) into three
triangles 𝐾𝑖𝑗 (finite element), 𝐾𝑖𝑗 ∈ Υℎ (their common vertex is in the barycenter of the macroelement
𝐿𝑖 ). Let 𝑅𝑙 and 𝑆𝑚 be the vertices and midpoints of the sides 𝐾𝑠 ∈ Υℎ , respectively. Introduce the
notation of sets:
1)𝑍 𝑣𝑒𝑙 = 𝑍Ω𝑣𝑒𝑙 ∪ 𝑍Γ𝑣𝑒𝑙 = {𝑅𝑙 ∪ 𝑆𝑚 }, where 𝑍Ω𝑣𝑒𝑙 and 𝑍Γ𝑣𝑒𝑙 are triangulation nodes subsets for the velocity
field components in Ω and on Γ, respectively;
2)𝑍 𝑝𝑟𝑒𝑠 = {𝑄𝑙 } of triangulation nodes for the Bernoulli pressure, where the node 𝑄𝑙 an exact match
to the node 𝑅𝑚 at the appropriate 𝐾𝑖𝑗 .
    We denote by Ωℎ = ⋃ 𝐾𝑠 the totality of the finite elements with sides of order ℎ. Next, we
                             𝐾𝑠 ∈Υℎ
describe the Scott-Vogelius (SV) element pair (see [12]). For the components of the velocity field, we
use polynomials of the second degree (𝑋 ℎ ), and for the pressure — the first one (𝑍 ℎ ):
𝑋 ℎ = {𝑣ℎ ∈ 𝐶(Ω) ∶ 𝑣ℎ |𝐾 ∈ 𝑃2 (𝐾 ), ∀𝐾 ∈ Υℎ }, 𝐗ℎ = 𝑋 ℎ × 𝑋 ℎ ;
𝑍 ℎ = {𝑧ℎ ∈ 𝐿2 (Ω) ∶ 𝑧ℎ |𝐾 ∈ 𝑃1 (𝐾 ), ∀𝐾 ∈ Υℎ , ∫ 𝑧ℎ 𝑑𝐱 = 0}.
                                                   Ω
  The SV pair has useful feature, namely div 𝐗ℎ ⊂ 𝑍 ℎ . Next, we represent special basis functions and
construct a scheme of the weighted finite element method. To each node 𝑀𝑚 ∈ 𝑍Ω𝑣𝑒𝑙 (𝑄𝑙 ∈ 𝑍 𝑝𝑟𝑒𝑠 ) we
associate the basis function
                                   ∗                           ∗
                  𝜃𝑚 (𝐱) = 𝜌 𝜈 (𝐱) ⋅ 𝜑𝑚 (𝐱), (𝜒𝑙 (𝐱) = 𝜌 𝜇 (𝐱) ⋅ 𝜓𝑙 (𝐱)), 𝑚 = 0, 1, … , ( 𝑙 = 0, 1, …),
                                                                                                    {
                                                                                                      1, 𝑙 = 𝑗,
where 𝜑𝑚 ∈ 𝑋 ℎ , 𝜑𝑚 (𝑀𝑗 ) = 𝛿𝑚𝑗 𝑚, 𝑗 = 0, 1, … (𝜓𝑙 ∈ 𝑍 ℎ , 𝜓𝑙 (𝑄𝑗 ) = 𝛿𝑙𝑗 , 𝑙, 𝑗 = 0, 1, …); 𝛿𝑙𝑗 =                 , 𝜈∗
                                                                                                      0, 𝑙 ≠ 𝑗,
and 𝜇 ∗ are real parameters.
   The spaces 𝑉 ℎ and 𝑄 ℎ for the components of the velocity field and pressure are defined as linear
span of the basis functions {𝜃𝑚 }𝑚 and {𝜒𝑙 }𝑙 , respectively. Let 𝑉0ℎ be a subspace of 𝑉 ℎ ∶ 𝑉0ℎ = {𝑣 ℎ ∈
𝑉 ℎ ∶ 𝑣ℎ (𝑀𝑚 )|𝑀𝑚 ∈𝑍 𝑣𝑒𝑙 = 0}. The approximate components of the velocity field 𝐮(𝑘)      𝜈,ℎ = (𝑢𝜈,ℎ,1 , 𝑢𝜈,ℎ,2 ) and
                                                                                                   (𝑘)     (𝑘)
                      Γ

pressure 𝑃𝜈,ℎ
          (𝑘)
              we seek as a

                 (𝑘)            (𝑘)          (𝑘)            (𝑘)           (𝑘)
                𝑢𝜈,ℎ,1 (𝐱) = ∑ 𝑑2𝑚  𝜃𝑚 (𝐱), 𝑢𝜈,ℎ,2 (𝐱) = ∑ 𝑑2𝑚+1 𝜃𝑚 (𝐱), 𝑃𝜈,ℎ (𝐱) = ∑ 𝑒𝑙(𝑘) 𝜒𝑙 (𝐱),                      (10)
                                  𝑚                                𝑚                                  𝑙

                                       (𝑘)
                                                          𝑖 . The coefficients 𝑑𝑗
where 𝑑𝑗(𝑘) = 𝜌 −𝜈 (𝑀[𝑗/2] ) 𝑑̃ 𝑗 , 𝑒𝑖(𝑘) = 𝜌 −𝜇 (𝑄𝑖 ) 𝑒̃ (𝑘)                       and 𝑒𝑖(𝑘) in (10) are found as
                      ∗                               ∗                         (𝑘)

a result of solving a system (11), (12) (see below). Let 𝐕ℎ = 𝑉 ℎ × 𝑉 ℎ , 𝐕ℎ0 = 𝑉0ℎ × 𝑉0ℎ and 𝐕ℎ ⊂
                          𝑜
𝐖12,𝜈 (Ωℎ , 𝛿), 𝐕ℎ0 ⊂𝐖12,𝜈 (Ωℎ , 𝛿), 𝑄 ℎ ⊂ 𝐿02,𝜈 (Ωℎ , 𝛿).
  Definition 2. The pair 𝐮(𝑘)   𝜈,ℎ ∈ 𝐕 and 𝑃𝜈,ℎ ∈ 𝑄 is called an approximate 𝑅𝜈 -generalized solution of
                                        ℎ         (𝑘)      ℎ

the problem (7), (8) for any pair 𝐯ℎ ∈ 𝐕ℎ0 and 𝑞 ℎ ∈ 𝑄 ℎ the equalities

                              𝑎𝑘 (𝐮(𝑘)    ℎ       ℎ (𝑘)          ℎ
                                   𝜈,ℎ , 𝐯 ) + 𝑏(𝐯 , 𝑃𝜈,ℎ ) = 𝑙(𝐯 )          and           𝑐(𝐮(𝑘)    ℎ
                                                                                              𝜈,ℎ , 𝑞 ) = 0
                                                                                                                         (11)

hold, where 𝐟 ∈ 𝐋2,𝛽 (Ω, 𝛿), 𝐠 ∈ 𝐖1/2
                                  2,𝛽 (Γ, 𝛿), 𝜈 ≥ 𝛽 ≥ 0.
  Thus, we construct a weighted FEM to find an 𝑅𝜈 -generalized solution for the problem (7), (8). We
get a system of linear algebraic equation:

                                         𝐀𝑘 𝐝(𝑘) + 𝐁𝐞(𝑘) = 𝜔           and         𝐂𝑇 𝐝(𝑘) = 𝟎,                          (12)

where 𝐝(𝑘) = (𝑑0(𝑘) , 𝑑2(𝑘) , … , 𝑑1(𝑘) , 𝑑3(𝑘) , …)𝑇 , 𝐞(𝑘) = (𝑒0(𝑘) , 𝑒1(𝑘) , 𝑒2(𝑘) , …)𝑇 and 𝜔 be a vector of values of the
linear form 𝑙(𝜃𝑚 ).


4. Iterative procedure
Now, we present an iterative procedure for solving the sequences of systems view (12), 𝑘 = 1, 2, 3, …
and thus we will approximately solve the original nonlinear problem in rotation form (4), (5):
1. Let 𝐝(0) and 𝐞(0) be an arbitrary vectors such that 𝐮(0)
                                                        𝜈,ℎ (𝑀𝑙 )|𝑀𝑙 ∈𝑍 𝑣𝑒𝑙 = 𝐠(𝑀𝑙 ), div 𝐮𝜈,ℎ (𝑀𝑙 )|𝑀𝑙 ∈𝑍 𝑣𝑒𝑙 = 0
                                                                                           (0)
                                                                                       Γ

(for example 𝐮(0)𝜈,ℎ (𝑀𝑙 )|𝑀𝑙 ∈𝑍Ω𝑣𝑒𝑙 = 𝟎) and 𝑃𝜈,ℎ (𝑀𝑙 )|𝑀𝑙 ∈𝑍
                                               (0)
                                                               𝑝𝑟𝑒𝑠 = 0.

2. Realize the Picard’s procedure 𝑘 = 0, 1, 2, … until the stopping condition is fulfilled:
a) Let 𝜁0(𝑘) ∶= 𝐝(𝑘) and 𝜂(𝑘)0 ∶= 𝐞 ;
                                       (𝑘)

b) We construct an internal convergent iterative process (see [17]). For 𝑛 = 0, 1, … , 𝑁𝑘 − 1 ∶
                                              (𝑘)              −1
                                             𝜁𝑛+1 = 𝜁𝑛(𝑘) + 𝐀̂ 𝑘 (𝜔 − 𝐀𝑘 𝜁𝑛(𝑘) − 𝐁𝜂(𝑘)
                                                                                   𝑛 )
                                                     (𝑘)
                                                    𝜂𝑛+1 = 𝜂(𝑘) ̂ −1 𝑇 (𝑘)
                                                            𝑛 + 𝐒𝑘 𝐂 𝜁𝑛+1 ;
c) Let 𝐝(𝑘+1) ∶= 𝜁𝑁(𝑘)
                    𝑘
                       and 𝐞(𝑘+1) ∶= 𝜂(𝑘)
                                      𝑁𝑘 ,
where 𝐀𝑘 and 𝐒𝑘 are the preconditioning matrices to 𝐀𝑘 and 𝐒𝑘 = 𝐂𝑇 𝐀−1
        ̂       ̂
                                                                           𝑘 𝐁, respectively.
   At first, we build a preconditioner 𝐀̂ 𝑘 applying an incomplete LU factorization. We employ the
                                                            −1
GMRES(5)-method (see [18]). If we have error 𝐫0 = 𝐀̂ 𝑘 (𝐬 − 𝐀𝑘 𝐯) for the problem 𝐀𝑘 𝐯 = 𝐬, then the
                                                                                  −1              −1
Arnoldi procedure will build an orthonormal basis of the subspace: Span{𝐫0 , 𝐀̂ 𝑘 𝐀𝑘 𝐫0 , … , (𝐀̂ 𝑘 𝐀𝑘 )4 𝐫0 }.
   Further, we construct an auxiliary matrix 𝐒̃ 𝑘 to 𝐒̂ 𝑘 , which is the weight mass matrix 𝐌𝜈,𝜇 ,𝜈̄ , such
                                                                                                     ∗


that on each 𝐿 ∈ Υℎ ∶
                                  ∗           1           ∗
                            (𝐌𝜈,𝜇 ,𝜈̄ )𝑙𝑚 =      ∫ 𝜌 2(𝜈+𝜇 ) 𝜓𝑙 (𝐱) 𝜓𝑚 (𝐱)𝑑𝐱, 𝑙, 𝑚 = 0, 1, … .
                                              𝜈̄
                                                𝐿
                                                                ∗                     ∗
  After that, we define a diagonal matrix 𝐒̄ 𝑘 = 𝐌
                                                 ̄ 𝜈,𝜇 ,𝜈̄ , where (𝐌
                                                                    ̄ 𝜈,𝜇 ,𝜈̄ ) = ∑ (𝐌𝜈,𝜇 ∗ ,𝜈̄ ) .
                                                                               𝑖𝑖                𝑖𝑙
                                                                                               𝑙
   It is known (see [19]), that such diagonal lumping 𝐒̄ 𝑘 is a good preconditioner to matrix 𝐒̃ 𝑘 . In order
                                        −1
to determining the vector 𝜓 ⋄ ∶= 𝐒̂ 𝑘 𝜃 we need to find a solution of the internal procedure:
1) 𝜙0 = 𝟎;
2) 𝜙𝑚 = 𝜙𝑚−1 + 𝐒̄ 𝑘 (𝜃 − 𝐒̃ 𝑘 𝜙𝑚−1 ) (𝑚 = 1, … , 𝑀);
                  −1

3) 𝜓 = 𝜙𝑀 .
      ⋄

We use the GMRES(5)-method: (Span{𝐫̄ , (𝐒̄ 𝑘 𝐒̃ 𝑘 )𝐫̄ , … , (𝐒̄ 𝑘 𝐒̃ )4𝑘 𝐫̄ }, 𝐫̄ = 𝐒̄ 𝑘 (𝜃 − 𝐒̃ 𝑘 𝜙𝑚−1 )).
                                                 −1             −1                     −1



5. Results of numerical experiments
Let Ω = (−1, 1)×(−1, 1)⧵[0, 1]×[−1, 0]. Then we split Ω̄ by the horizontal and vertical lines 𝑥1(𝑗) = −1+𝑗 ℎ,
and 𝑥2(𝑖) = −1 + 𝑖 ℎ, respectively, into elementary squares 𝑆𝑙 , where 𝑗, 𝑖 = 0, … , 𝑁 , ℎ = 𝑁2 , 𝑁 – even
number. After that, we divide each 𝑆𝑙 by the diagonal (the lower left corner connects to the upper
right corner) into two triangles 𝐿𝑚 (macroelements). Further, each macroelement 𝐿𝑚 is partitioned
into three triangles 𝐾𝑠 (barycentric partition). Consider a solution (𝐮, 𝑃) of nonlinear problem (4), (5)
which has a singularity in the vicinity of the reentrant corner 𝜎 = 32𝜋 with apex at the origin (0, 0) ∶
                                      𝜆                                       𝑥2                            𝑥2
          𝑢1 (𝑥1 , 𝑥2 ) = (𝑥12 + 𝑥22 ) 2 ((1 + 𝜆) 𝐸(𝑥1 , 𝑥2 ) ⋅ sin(arctg        ) + 𝐺(𝑥1 , 𝑥2 ) ⋅ cos(arctg )),
                                                                              𝑥1                            𝑥1
                                      𝜆                               𝑥2                                    𝑥2
          𝑢2 (𝑥1 , 𝑥2 ) = (𝑥12 + 𝑥22 ) 2 (𝐺(𝑥1 , 𝑥2 ) ⋅ sin(arctg        ) − (1 + 𝜆) 𝐸(𝑥1 , 𝑥2 ) ⋅ cos(arctg )),
                                                                      𝑥1                                    𝑥1
                                                      𝜆−1   (1 + 𝜆)2 𝐺(𝑥1 , 𝑥2 ) + 𝐻 (𝑥1 , 𝑥2 )
                           𝑃(𝑥1 , 𝑥2 ) = (𝑥12 + 𝑥22 ) 2 (                                       ),
                                                                        𝜆−1
where
                                                                𝑥2                      𝑥2
                             𝐸(𝑥1 , 𝑥2 ) = cos((1 − 𝜆) arctg       ) − cos((1 + 𝜆) arctg )+
                                                                𝑥1                      𝑥1
                          sin((1 + 𝜆) arctg 𝑥𝑥21 ) ⋅ cos(𝜆 𝜎)       sin((1 − 𝜆) arctg 𝑥𝑥21 ) ⋅ cos(𝜆 𝜎)
                      +                                         +                                         ,
                                          𝜆+1                                     𝜆−1
                                                      𝑥2                              𝑥2
                   𝐺(𝑥1 , 𝑥2 ) = (1 + 𝜆) sin((1 + 𝜆) arctg
                                                         ) − (1 − 𝜆) sin((1 − 𝜆) arctg )+
                                                      𝑥1                              𝑥1
                                       𝑥2                                 𝑥2
                     +cos((1 + 𝜆) arctg ) ⋅ cos(𝜆 𝜎) − cos((1 − 𝜆) arctg ) ⋅ cos(𝜆 𝜎),
                                       𝑥1                                 𝑥1
Table 1
The error between the generalized solution 𝐮ℎ and exact one 𝐮 in the 𝐖12 (Ω) norm.

                                      N=        280             140            70
                                             1.379e-1      1.988e-1       2.848e-1


Table 2
The error between an 𝑅𝜈 -generalized solution 𝐮𝜈,ℎ and exact one 𝐮 in the 𝐖12,𝜈 (Ω) norm, for different values
𝜈, 𝛿 and 𝜈 ∗ (𝜇 ∗ = 𝜈 ∗ ).
                                 (𝜈, 𝜈 ∗ , 𝛿), 𝑁 =        280           140           70
                                (2.0, −0.5, 0.015)    1.627e-5        3.295e-5      6.629e-5
                               (2.0, 𝜆 − 1, 0.015)    1.394e-5        2.796e-5      5.626e-5
                                (2.0, −0.4, 0.015)    1.181e-5        2.355e-5      4.760e-5
                                (1.9, −0.5, 0.016)    2.785e-5        5.624e-5      1.131e-4
                               (1.9, 𝜆 − 1, 0.016)    2.331e-5        4.672e-5      9.368e-5
                                (1.9, −0.4, 0.016)    1.869e-5        3.759e-5      7.549e-5


Table 3
The number of grid nodes 𝑀𝑖 ∈ 𝑍Ω𝑣𝑒𝑙 (in percentage of their total number), where the errors 𝛽𝑗𝑖 are more than
the given limit values 𝜀𝑙 , 𝑙 = 1, 2, of the generalized solution 𝜈 = 𝜈 ∗ = 0, 𝛿 = 1.
                                                     N=     280          140
                                           𝜀1 = 10−5      64.51%       77.82%
                                           𝜀2 = 10−6      90.97%       94.81%


Table 4
The number of grid nodes 𝑀𝑖 ∈ 𝑍Ω𝑣𝑒𝑙 (in percentage of their total number), where the errors 𝛽𝜈,𝑗     𝑖 are more than
                                                                                ∗
the given limit values 𝜀𝑙 , 𝑙 = 1, 2, of an 𝑅𝜈 - generalized solution 𝜈 = 1.8, 𝜈 = −0.31, 𝛿 = 0.014.
                                                     N=     280          140
                                           𝜀1 = 10−5      28.33%       54.30%
                                           𝜀2 = 10−6      77.32%       85.76%


                                                         𝑥2                               𝑥2
                  𝐻 (𝑥1 , 𝑥2 ) = (1 − 𝜆)3 sin((1 − 𝜆) arctg ) − (1 + 𝜆)3 sin((1 + 𝜆) arctg )+
                                                         𝑥1                               𝑥1
                                           𝑥                                          𝑥
            +(𝜆 − 1)2 cos((1 − 𝜆) arctg ) ⋅ cos(𝜆 𝜎) − (𝜆 + 1)2 cos((1 + 𝜆) arctg ) ⋅ cos(𝜆 𝜎).
                                            2                                           2
                                           𝑥1                                         𝑥1
   Let 𝛼 = 𝜈̄ = 1 and 𝜆 = 0.54448. The pair of functions (𝐮, 𝑃) is analytic in Ω̄ ⧵ (0, 0), but 𝐮 ∉ 𝐖22 (Ω)
and 𝑃 ∉ 𝑊21 (Ω). It is a typical situation in non-convex polygonal domains.
   Numerical experiments were carried out on grids with different steps ℎ. The errors of the gen-
eralized (classical FEM with 𝜈 = 0, 𝛿 = 1, 𝜈 ∗ = 𝜇 ∗ = 0) and 𝑅𝜈 -generalized (presented weighted
FEM) solutions were determined using the modulus of the difference between the exact solution
and approximate one at the nodes 𝑀𝑖 , i. e. 𝛽𝑗𝑖 = |𝑢𝑗 (𝑀𝑖 ) − 𝑢ℎ,𝑗 (𝑀𝑖 )| for the generalized solution and
                  𝜈,ℎ,𝑗 (𝑀𝑖 )| for an 𝑅𝜈 -generalized one, where 𝑀𝑖 ∈ 𝑍Ω , 𝑗 = 1, 2, and also in the norms
 𝑖 = |𝑢 (𝑀 ) − 𝑢
𝛽𝜈,𝑗                                                                        𝑣𝑒𝑙
        𝑗   𝑖
of generalized functions. See Figures 1-2 and Tables 1-4. The optimal values of parameters 𝜈, 𝜈 ∗ and
𝛿 were derived numerically.
Figure 1: Distribution of the points 𝑀𝑖 with errors 𝛽 of the generalized solution (𝜈 = 0, 𝛿 = 1, 𝜈 ∗ = 𝜇 ∗ = 0):
𝑎) 𝑁 = 140, 𝑐) 𝑁 = 280 and 𝑏) 𝑁 = 140, 𝑑) 𝑁 = 280 for the 1st and 2nd components of 𝐮ℎ , respectively.


6. Conclusions
The results of computational experiments for the steady Navier-Stokes equations (4), (5) lead to the
following conclusions:
   1) An approximate 𝑅𝜈 -generalized solution by the weighted FEM converges to the exact one with
a (ℎ) rate in the 𝐖12,𝜈 (Ω) norm (see Table 2), while the approximate generalized solution by the
classical FEM converges to the exact one with a (ℎ0.54 ) rate in the 𝐖12 (Ω) norm (see Table 1). In
other words, the proposed method suppresses the so-called pollution effect [13].
   2) For all values of 𝛿, 𝜈 and 𝜈 ∗ from the range of optimal values (𝛿 ∼ ℎ, 𝜈 ∼ 2 and 𝜈 ∗ ∼ 1 − 𝜆) an
approximate 𝑅𝜈 -generalized solution converges to the exact one with a (ℎ) rate in the 𝐖12,𝜈 (Ω) norm.
   3) The number of nodes and their surroundings by using a weighted FEM, in which the values of
the absolute errors 𝛽𝜈,𝑗
                      𝑖 , 𝑗 = 1, 2, do not exceed the given values, increases with 𝑁 and is much more

then by using the classical FEM (see Tables 3-4) and Figures 1-2.


Acknowledgments
The reported study was supported by RSF according to the research project No. 21-11-00039. Compu-
tational resources for the numerical experiments were provided by the Shared Services Center "Data
Center of FEB RAS".
Figure 2: Distribution of the points 𝑀𝑖 with errors 𝛽𝜈 of the 𝑅𝜈 -generalized solution (𝜈 = 1.9, 𝛿 = 0.014, 𝜈 ∗ =
𝜇 ∗ = −0.35): 𝑎) 𝑁 = 140, 𝑐) 𝑁 = 280 and 𝑏) 𝑁 = 140, 𝑑) 𝑁 = 280 for the 1st and 2nd components of 𝐮𝜈,ℎ ,
respectively.


References
 [1] V. A. Rukavishnikov, On the differential properties of 𝑅𝜈 -generalized solution of Dirichlet prob-
     lem, Dokl. Akad. Nauk SSSR 309 (1989) 1318–1320.
 [2] V. A. Rukavishnikov, On the existence and uniqueness of an 𝑅𝜈 -generalized solution of a bound-
     ary value problem with uncoordinated degeneration of the input data, Doklady Mathematics 90
     (2014) 562–564. doi:10.1134/S1064562414060155.
 [3] V. A. Rukavishnikov, E. V. Kuznetsova, The 𝑅𝜈 -generalized solution of a boundary value problem
     with a singularity belongs to the space 𝑊2,𝜈+𝛽/2+𝑘+1
                                                𝑘+2       (Ω, 𝛿), Differential Equations 45 (2009) 913–
     917. doi:10.1134/S0012266109060147.
 [4] V. A. Rukavishnikov, S. G. Nikolaev, On the 𝑅𝜈 -generalized solution of the Lamé system with cor-
     ner singularity, Doklady Mathematics 92 (2015) 421–423. doi:10.1134/S1064562415040080.
 [5] V. A. Rukavishnikov, A. V. Rukavishnikov, Weighted finite element method for the Stokes prob-
     lem with corner singularity, Journal of Computational and Applied Mathematics 341 (2018)
     144–156. doi:10.1016/j.cam.2018.04.014.
 [6] V. A. Rukavishnikov, A. Y. Bespalov, An exponential rate of convergence of the finite element
     method for the Dirichlet problem with a singularity of the solution, Doklady Mathematics 62
     (2000) 266–270.
 [7] V. A. Rukavishnikov, E. V. Kuznetsova, A finite element method scheme for boundary value
     problems with noncoordinated degeneration of input data, Numerical Analysis and Applications
     2 (2009) 250–259. doi:10.1134/S1995423909030069.
 [8] V. A. Rukavishnikov, E. I. Rukavishnikova, Numerical method for Dirichlet problem with de-
     generation of the solution on the entire boundary, Symmetry 11 (2019) 1455. doi:10.3390/
     sym11121455.
 [9] V. A. Rukavishnikov, A. O. Mosolapov, E. I. Rukavishnikova, Weighted finite element method
     for elasticity problem with a crack, Computers and Structures 243 (2021) 106400. doi:10.1016/
     j.compstruc.2020.106400.
[10] V. A. Rukavishnikov, O. P. Tkachenko, Dynamics of a fluid-filled curvilinear pipeline, Applied
     Mathematics and Mechanics 39 (2018) 905–922. doi:10.1007/s10483-018-2338-9.
[11] M. Benzi, H. Golub, J. Liesen, Numerical solution of saddle point problems, Acta Numerica 14
     (2005) 1–137. doi:10.1017/S0962492904000212.
[12] L. R. Scott, M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator
     in spaces of piecewise polynomials, Mathematical Modelling and Numerical Analysis 19 (1985)
     111–143.
[13] H. Blum, The influence of reentrant corners in the numerical approximation of viscous flow
     problems, volume 30 of Numerical Treatment of the Navier-Stokes Equations, Springer, 1990.
     doi:110.1007/978-3-663-14004-7-4.
[14] V. A. Rukavishnikov, A. V. Rukavishnikov, New approximate method for solving the Stokes
     problem in a domain with corner singularity, Bulletin of the South Ural State University. Ser.
     Mathematical Modelling, Programming & Computer Software 11 (2018) 95–108. doi:10.14529/
     mmp180109.
[15] V. A. Rukavishnikov, A. V. Rukavishnikov, New numerical method for the rotation form of the
     Oseen problem with corner singularity, Symmetry 11 (2019) 54. doi:10.3390/sym11010054.
[16] V. A. Rukavishnikov, A. V. Rukavishnikov, The method of numerical solution of the one sta-
     tionary hydrodymics problem in convective form in 𝐿-shaped domain, Computer Research and
     Modeling 12 (2020) 1291–1306. doi:10.20537/2076-7633-2020-12-6-1291-1306.
[17] J. H. Bramble, J. E. Pasciak, A. T. Vassilev, Analysis of the inexact Uzawa algorithm for sad-
     dle point problems, SIAM Journal on Numerical Analysis 34 (1997) 1072–1092. doi:10.1137/
     S0036142994273343.
[18] Y. Saad, Iterative methods for sparse linear systems, 2nd. ed., SIAM, Philadelphia, 2003.
[19] M. Olshanskii, A. Reusken, Analysis of a Stokes interface problem, Numerische Mathematik 103
     (2006) 129–149. doi:10.1007/S00211-005-0646-X.