An Optimal Finite Element Method with Uzawa Iteration for Stokes Equations including Corner Singularities

.e Uzawa method is an iterative approach to find approximated solutions to the Stokes equations. .is method solves velocity variables involving augmented Lagrangian operator and then updates pressure variable by Richardson update. In this paper, we construct a new version of the Uzawa method to find optimal numerical solutions of the Stokes equations including corner singularities..e proposedmethod is based on the dual singular functionmethod which was developed for elliptic boundary value problems. We estimate the solvability of the proposed formulation and special orthogonality form for two singular functions. Numerical convergence tests are presented to verify our assertion.


Introduction
To study the solution of a partial differential equation, the equation is sometimes interpreted in a weak (variational) sense and we can define the regular problem in this manner. Consider a variational problem: where H m 0 ⊂ V ⊂ H m (Ω) and a(·, ·) is bilinear form with continuity and coercivity. We call the problem is H s -regular if, for every f ∈ H s− 2m (Ω), there is a solution u ∈ H s (Ω) with a constant c such that It is well known that the approximated solution of regular problems shows optimal convergence by using common numerical methods such as the finite difference method or the finite element method.
However, if a computational domain of an elliptic boundary value problem is a polygon including reentrant corners, then the problem is not H s -regular, s ≥ 2, and it is sometimes called a corner singular problem, and the singular solution is hard to approximate optimally. erefore, there have been many numerical approaches developed to solve the problem efficiently.
ese methods aim to improve accuracy and to resolve the convergence difficulties. e hp version of the finite element method is a typical method. It employs elements of variable sizes h and polynomial degrees p to improve the convergence rate and accuracy in [1]. e postprocessing method has been proposed in [2][3][4][5][6]. ese methods calculate the coefficients of the singular coefficients from the finite element solution. Singular function boundary integral method also has been proposed to treat singular problems in [7][8][9]. ese methods calculate the unknown coefficients of singular functions directly. e solution is approximated by the leadinig terms of the local asymptotic solution expansion, and the Dirichlet boundary conditions are weakly enforced utilizing Lagrange multipliers.
It is well known that the solution of an elliptic boundary value problem can be decomposed of a finite number of socalled singular functions, which come from the neighborhood of a corner point, and a smooth remainder function which is called a regular solution (see, for instance, [10][11][12]). e singular function method is one approach to use this fact. It takes into account the form of the singular solution by finding stress intensity factors and regular solution (see, for example, [13][14][15]). ese methods consist of augmenting the spline space by the singular functions. e dual singular function is based on the observation that the regular part of the solution and the stress intensity factor are related to each other by testing a dual singular function in [16][17][18]. It was implemented as an iterative procedure that iterates back and forth between these equations. Also, this approach was extended to multigrid versions.
In [19], Cai and Kim developed and analyzed a finite element method for the accurate computation of the solution and intensity factors. If a regular part of the solution is smoother than the solution itself, then approximated solutions of a standard finite element lose accuracy. erefore, the new method finds approximated regular solution first and then computes stress intensity factors and solution. is method decoupled variational formulations by testing the dual singular function and disjoint cut-off function. Consequently, the regular solution is uniquely determined by a well-posed variational problem, and stress intensity factors can be expressed by a regular solution. We call this method the finite element dual singular function method (FE-DSFM).
is method was extended to other problems. Poisson problems with mixed boundary conditions and interface problems were applied in [20,21]. Also, this algorithm was studied to corner singularities of the Helmholtz equation and heat equation (see [22,23]). Some studies have been designed to target the singular solution of the Stokes equation which is much more complicated. We refer the interested reader to the papers [24,25] for a formula for corner singularity of the Stokes equations. A mixed finite element method-based FE-DSFM was developed for Stokes equations in [26,27]. ese proved the accuracy and wellposedness of the algorithm including two singularities at each corner.
In this paper, we proposed a new algorithm for solving Stokes equations. e proposed algorithm is based on the Uzawa algorithm and finite element method in [28]. We estimate the solvability of the proposed formulation and special orthogonality form for two singular functions. We also give numerical convergence tests to verify our assertion.

Singular Functions of Stokes Equations
Let the computational domain Ω be an open and bounded concave polygon in R 2 especially having one reentrant corner. e steady-state Stokes equation is where f is a given external force field causing an acceleration of the flow in H − 1 (Ω), Ω is a computational domain in R 2 , and μ � Re − 1 is the reciprocal of the Reynolds number. e unknowns are the (vector) velocity field u ∈ H 1 0 (Ω) and the (scalar) pressure p ∈ L 2 0 (Ω). e pressure gradient plays a role in an additional force, which prevents a change in the density. In particular, high pressure builds up at points, where, otherwise, a source of the sink would be created. Mathematically, the pressure can be considered as a Lagrange multiplier. Besides, the weak formulation of the Stokes equations leads to a saddle point problem with the restriction ∇ · u � 0.
ere are mainly two approaches to find finite element approximation of solutions of (3). One approach is called the mixed finite element method which solves velocities and pressure simultaneously by constructing a big linear system. Another approach is an iteration method called the Uzawa method which solves velocity variables involving augmented Lagrangian operator and then updates pressure variable by Richardson update. e advantage of the Uzawa method is that it uses less memory because it solves velocity and pressure separately. However, the iteration process sometimes takes more computational time.
e Uzawa iteration can be extended to projection-type methods such as the gauge-Uzawa method to solve unsteady incompressible Navier-Stokes equations.
If the solution of (3) is smooth enough, namely, (u, p) ∈ H s+1 (Ω) × H s (Ω) with s ≥ 1, and if a suitable finite element pair is imposed for velocity and pressure, then the finite element solution (u h , p h ) using the standard mixed finite element method has optimal error bounds as shown in [14]: where h is the biggest mesh size. However, if s < 1, then the error bounds only become For the case s < 1, we call the solution (u, p) a singular solution, otherwise a regular solution. Since singularities are due to re-entrant corners of a computational domain Ω, we assume that Ω is open and bounded polygonal domain in R 2 with one re-entrant corner.
To derive singular and dual singular functions for the Stokes equations, the polar coordinate (r, θ) of homogeneous Stokes system (3) should be considered, and we can find the singular and dual singular solution via solving the homogeneous system with the separation of variables. en, we arrive at the singular function of system (3): Journal of Mathematics where λ is solution of and C 1 and C 2 are Similarly, the dual singular functions are derived by where λ is solution of (7) and D 1 and D 2 are e following lemma describes the number of singular functions depending on the value of ω and the relationship between values of λ and ω is shown in Figure 1.
It is shown that 2 singular functions of solution (3) can be involved in each reentrant corner. If there is one reentrant corner, the solution of (3) can be written by the form where α 1 and α 2 are the stress intensity factors, , are singular functions of (6), and (w, q) ∈ H 2 (Ω) × H 1 (Ω) is the regular solution.
Lemma 2 (see [26]). The singular function respectively. e boundary conditions of u s i and u d i vanish on Γ in , but the boundary value of u d i is not defined at the origin. Both of u s i and u d i are not 0 on Γ out � zΩ/Γ in .

FE-DSFM and Uzawa Iteration
Uzawa method solves velocity variables involving augmented Lagrangian operator and then updates the pressure variable by Richardson update. Algorithm 1 is the Uzawa method for solving stokes equation (3). ere are some theoretical works of literature on the Uzawa method and augmented Lagrangian method. In [29][30][31][32], the convergence of Algorithm 1 is proved. Especially, it is shown that κ � 0 and β � 1 have optimal Journal of Mathematics convergence in [31]. So, we choose κ � 0 and β � 1 for our algorithms.
Our goal is to construct FE-DSFM with the Uzawa iteration. We consider Stokes equations with polygonal domain including one reentrant corner. To represent FE-DSFM, we first introduce the cut-off function. Let ω be the internal angle of Ω satisfying π < ω < 2π. Now, we set and B(r 1 ) � B(0; r 1 ) in the polar coordinates (r, θ) and define cut-off function where δ(r) is very smooth function, ρ is a parameter in (0, 2], and R is a fixed real number which will be determined later so that η 2ρ has 0 on whole Γ out . We note that the singular function (6) Let η be a smooth cut-off function (14) with smooth function which is equal to the one identically in the neighborhood of origin, and the support of η is small enough so that the functions ηu s i , i � 1, 2, vanish identically on zΩ.
en, in general, the solution (u, p) including singular parts of (3) can be rewritten in the form (11): where α 1 and α 2 are the stress intensity factors and e strategy of FE-DSFM is to find the regular solution (w, q) ∈ H 2 (Ω) × H 1 (Ω) and stress intensity factors α 1 and α 2 by applying the standard Uzawa method with an additional equation. To make simpler formula, we assume that there is only one singular solution, that is, one of α 1 or α 2 is zero. Let u n+1 � w n+1 + α n+1 η ρ u s and p n � q n + α n η ρ p s . en, Step 1 of the Uzawa algorithm can be Even though − μΔη ρ u s and ∇η ρ p s ∉ L 2 , we know that − μ△(η ρ u s ) + ∇(η ρ p s ) ∈ L 2 . So, (16) can be rewritten as For the sake of a clear explanation, we note that the inner product of vectors a � (a 1 , a 2 ) and en, the weak formulation of (17) is (19) is solvable, and it has unique solution w n+1 ∈ H 1 0 (Ω) if we know α n+1 ∈ R. However, since α n+1 is unknown, we need one more equation, and it should be (i) A single equation because α is a real number (ii) Linearly independent to (19) (iii) α which is not disappeared (iv) Solvable near the singularity corner Since equation (19) is a Poisson-type problem to find w n+1 , we can make the additional equation by testing the dual singular function of the following Poisson equation, instead of that of Stokes equation (9): Step 1: find u n+1 ∈ H 1 0 (ω) as the solution of − μΔu n+1 − κ∇(∇ · u n+1 ) + ∇p n � f, in Ω Step 2: update p n+1 ∈ L 2 (Ω) from the Richardson update p n+1 � p n − βμ∇ · u n+1 ALGORITHM 1: (Uzawa method). Given any initial guess p 0 ∈ L 2 (Ω), repeat the following steps until ‖u n+1 − u n ‖ ≤ tolerance: 4 Journal of Mathematics Lemma 3 (see [19]). The function s d L in (20) has the following properties: (17) and apply Lemma 3 to obtain Let β s � − △η ρ u s + 〈∇η ρ p s , η 2ρ S d L 〉 and β f � 〈f, η 2ρ S d L 〉 − 〈∇q n , η 2ρ S d L 〉. en, we have the system We first check solvability of system (22).
Proof. By Lemma 4, contract mapping theorem provides the existence of the unique fixed point of F ∘ T. We now define finite element space similarly as in previous sections to construct fully discrete FE-DSFM. Let T � K { } be a shape-regular quasi-uniform partition of Ω of mesh size h into closed elements K. e vector and scalar finite element spaces are Journal of Mathematics where P(K) and Q(K) are spaces of polynomials with degree bounded uniformly with respect to K ∈ T. We stress that the space P h is composed of continuous functions to use integration by parts: for all q h ∈ P h . en, the finite element approximation for (22) is to find w h ∈ V h and α n+1 h ∈ R such that And, the matrix form of the coupled system (30) becomes where A is symmetric positive definite square matrix and a and b are column vectors. System (31) can be solved by the Sherman-Morrison formula: and α n+1 h can be computed by second equation of (30). Sherman-Morrison formula is a rigorous approach to solve system (30), but it is sometimes complicated to apply for some problems including many singular functions. Since (F ∘ T)(α n ) � α n+1 by Lemma 4, we choose α n h instead of α n+1 h in first equation of (30). Algorithm 2 is the proposed method for one singular function.
□ Remark 1. Compared with the original Uzawa iteration (Algorithm 1), the proposed method (Algorithm 2) needs only one more linear solver for the initial process.

Algorithm for Two Singular Functions
In this section, we construct an algorithm for two singularities in one corner. We consider the solution is With the implicit formation of two stress intensity factors, step 1 of Uzawa algorithm is Here, the unknowns are w n+1 , α n+1 1 , and α n+1 2 . However, our test function which is not in H 1 0 is only η 2ρ S d L . Instead of finding another test function, we use the following property.

(38)
Since Step We can conclude that equation (35) holds if However, we could not show equations (40)-(43) are valid by algebraically computation. erefore, we calculate them numerically.
For numerical substitution, we first compute the above integral forms: ω λω − π − ω λω + π C 1 2λ 2 + 3λ (1 + cos(λω)) − C 2 2λ 2 + 7λ + 2 sin(λω) We remark that the values λ 1 and λ 2 are from the equation and it is not possible to solve directly. So, we use the numerical root-finding method and bisection method for various tolerances of approximated λ. For internal re-entrant corner angle ω, we choose 400 angles from π to 2π by 0.0025π interval. Denote that (43) We choose three tolerances 1e − 12 , 1e − 13 , and 1e − 14 for solving λ 1 and λ 2 . By substituting λ 1 , λ 2 , and ω to (44)-(47). Figure 2 shows numerical substitution values of U A V A , U A V B , U B V A , and U B V B . We can see that the values strongly depend on the tolerances. If λ 1 and λ 2 are exact real values, we can claim our assertion. Now, we construct the Uzawa method for two singular solutions. We assume that the form (33) is valid. e weak formulation of Step 1 in the Uzawa method test by η 2ρ S L is in the form Let w � (w 1 , w 2 ) and f � (f 1 , f 2 ). en, x component of (49) is and y component of (49) is If we use the following notations then (50) and (51) can be rewritten by (53) erefore, we obtain Finally, combining (54)

Numerical Experiments
In this section, we provide numerical simulations using the standard Uzawa method and Algorithm 3. Our goal is to check their performance for the numerical solution u h and p h , which is composed of a regular solution w h and q h and stress intensity factors α 1h and α 2h . All numerical integration is calculated by the Gaussian quadrature 6 points. We notate Γ out � (r, θ) ∈ zΩ : θ ≠ 0 and θ ≠ ω { }.
Minielement chooses V h and P h are the piecewise P1bubble and piecewise P1 space pair satisfying discrete infsup condition. In this finite dimensional space, we desire ‖u − u h ‖ 0 ≤ Ch 2 , ‖Iu − u h ‖ ∞ ≤ Ch 2 , ‖u − u h ‖ 1 ≤ Ch, and ‖p − p h ‖ 0 ≤ Ch as an optimal convergence for some constant C which does not depend on numerical solution (u, p). In Figure 6, standard Uzawa iteration loses accuracy of u and p because the solutions contain singular functions. On the contrary, proposed Uzawa algorithm shows optimal accuracy of w, q, α 1 , and α 2 . Since u only depends w, α 1 , and α 2 and similarly p only depends q, α 1 , and α 2 . e figure guarantees convergence rate of u and p. e errors of stress intensity factors α 1 and α 2 decrease somewhat irregularly. However, we can see that the overall decay of both errors follows the optimal rate. e next finite element space is the Taylor-Hood space which chooses V h and P h as the piecewise P2 and piecewise P1 space pair.
is also satisfies the discrete inf-sup condition. In this space, the desired convergence rates are ‖u − u h ‖ 0 ≤ Ch 3 , ‖Iu − u h ‖ ∞ ≤ Ch 3 , ‖u − u h ‖ 1 ≤ Ch 2 , and ‖p − p h ‖ 0 ≤ Ch 2 , and Figure 7 shows the results calculated by each method. Even in this case, standard Uzawa iteration shows insufficient performance, which is     almost the similar convergence rate as in the minielement case, despite higher-order finite element space. In contrast, the proposed Uzawa algorithm presents optimal results.

Conclusions
We have presented a finite element using dual singular functions with Uzawa iteration for the steady Stokes equations including corner singularities. e proposed method uses the fact that the solution of an elliptic boundary value problem can be decomposed of a finite number of singular functions and a smoother remainder function which is called a regular solution. For the decoupled variational formulations with Uzawa iteration, we test the dual singular function of the Laplace problem and another cut-off function. Also, orthogonality of singular function of Stokes equation was presented, which imply that the proposed algorithm is available for two singular functions in the same corner. A numerical example indicates that the scheme is optimally accurate in both velocity and pressure. Besides, the scheme requires only one additional linear solver in the initial step, so it is much more cost-effective than the conventional dual singular function method for solving the Sherman-Morrison formulation. Since our method is based on the Uzawa iteration, it can be an extended version to solve the Navier-Stokes equation later.

Data Availability
No data were used to support the findings of the study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.