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

Authors: Jae-Hong Pyo; Deok-Kyu Jang
Journal: Journal of Mathematics
 (2021)

---

## Abstract

The Uzawa method is an iterative approach to find approximated solutions to the Stokes equations. This 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. The proposed method is based on the dual singular function method 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.


---

## Body

## 1. 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:(1)au,v=f,vL2,for allv∈V,where H0m⊂V⊂HmΩ and a⋅,⋅ is bilinear form with continuity and coercivity. We call the problem is Hs-regular if, for every f∈Hs−2mΩ, there is a solution u∈HsΩ with a constant c such that(2)us≤cfs−2m.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 notHs-regular, s≥2, and it is sometimes called a corner singular problem, and the singular solution is hard to approximate optimally. Therefore, there have been many numerical approaches developed to solve the problem efficiently. These methods aim to improve accuracy and to resolve the convergence difficulties. The 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].The postprocessing method has been proposed in [2–6]. These 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–9]. These methods calculate the unknown coefficients of singular functions directly. The 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 so-called 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–12]). The 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–15]). These methods consist of augmenting the spline space by the singular functions. The 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–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. Therefore, the new method finds approximated regular solution first and then computes stress intensity factors and solution. This 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).This 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]. These proved the accuracy and well-posedness of the algorithm including two singularities at each corner.In this paper, we proposed a new algorithm for solving Stokes equations. The 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.


## 2. Singular Functions of Stokes Equations

Let the computational domainΩ be an open and bounded concave polygon in ℝ2 especially having one reentrant corner. The steady-state Stokes equation is(3)−μ△u+∇p=f,inΩ,∇⋅u=0,inΩ,u=0,on∂Ω,where f is a given external force field causing an acceleration of the flow in H−1Ω, Ω is a computational domain in ℝ2, and μ=Re−1 is the reciprocal of the Reynolds number. The unknowns are the (vector) velocity field u∈H01Ω and the (scalar) pressure p∈L02Ω. The 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.There 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. The 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. The 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∈Hs+1Ω×HsΩ with s≥1, and if a suitable finite element pair is imposed for velocity and pressure, then the finite element solution uh,ph using the standard mixed finite element method has optimal error bounds as shown in [14]:(4)u−uh0+hu−uh1+hp−ph0≤Chs+1us+1+ps,where h is the biggest mesh size. However, if s<1, then the error bounds only become(5)u−uh0+hsu−uh1+hsp−ph0≤Ch2sus+1+ps.For the cases<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 ℝ2 with one re-entrant corner.To derive singular and dual singular functions for the Stokes equations, the polar coordinater,θ 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. Then, we arrive at the singular function of system (3):(6)usvsps=C1rλμλsinθsin1−λθrλμsinλθ−λsinθcos1−λθ−2rλ−1λcos1−λθ−C2rλμsinλθ+λsinθcos1−λθrλμλsinθsin1−λθ2rλ−1λsin1−λθ,where λ is solution of(7)λ2sin2ω=sin2λω,and C1 and C2 are(8)C1=sinλω+λsinωcos1−λω,C2=λsinωsin1−λω.Similarly, the dual singular functions are derived by(9)udvdpd=D1−r−λμλsinθsin1+λθ−r−λμsinλθ−λsinθcos1+λθ2r−λ−1λcos1+λθ+D2r−λμsinλθ+λsinθcos1+λθr−λμλsinθsin1+λθ2r−λ−1λsin1+λθ,where λ is solution of (7) and D1 and D2 are(10)D1=sinλω+λsinωcos1+λω,D2=λsinωsin1+λω.The 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.Figure 1
λ1 and λ2 values by numerical calculation by the bisection root-finding method.Lemma 1.
(see [26]). For χ:≈1.430296653124203, we have the following properties:(1)
λ=0 and λ=1 are solutions for any ω. We call these trivial solutions.(2)
There are only trivial solutions for0≤ω≤π. It means that there is no u,p∉H2Ω×H1Ω, solution of Stokes equation with homogeneous boundary condition on Γin, where Γin is the boundary near re-entrant corner.(3)
Ifπ≤ω≤χπ, then there is a unique solution λ except trivial solutions. Therefore, there is one singular solution u,p∈H1Ω×L2Ω.(4)
Ifχπ<ω<2π, then there are two solutions λ except trivial solutions. Thus, there is two singular solutions u,p∈H1Ω×L2Ω.(5)
Ifω=2π, then solution is λ=0.5 except trivial solutions.It is shown that 2 singular functions of solution (3) can be involved in each reentrant corner. If there is one re-entrant corner, the solution of (3) can be written by the form(11)up=wq+α1u1sp1s+α2u2sp2s,where α1 and α2 are the stress intensity factors, uis,pis∉H2Ω×H1Ω, i=1,2, are singular functions of (6), and w,q∈H2Ω×H1Ω is the regular solution.Lemma 2 (see [26]).
The singular functionuis,pis∈H1+λΩ×HλΩ and the dual singular function uid,pid∉H1Ω×L2Ω, i=1,2, satisfy(12)−μΔuis+∇pis=0,inΩ,∇⋅uis=0,onΓin,−μΔuid+∇pid=0,inΩ,∇⋅uid=0,onΓin,respectively. The boundary conditions of uis and uid vanish on Γin, but the boundary value of uid is not defined at the origin. Both of uis and uid are not 0 on Γout=∂Ω/Γin.


## 3. FE-DSFM and Uzawa Iteration

Uzawa method solves velocity variables involving augmented Lagrangian operator and then updates the pressure variable by Richardson update. Algorithm1 is the Uzawa method for solving stokes equation (3).Algorithm 1: (Uzawa method). Given any initial guess p0∈L2Ω, repeat the following steps until un+1−un≤tolerance:
Step 1: find un+1∈H01ω as the solution of −μΔun+1−κ∇∇⋅un+1+∇pn=f,inΩStep 2: update pn+1∈L2Ω from the Richardson update pn+1=pn−βμ∇⋅un+1There are some theoretical works of literature on the Uzawa method and augmented Lagrangian method. In [29–32], the convergence of Algorithm 1 is proved. Especially, it is shown that κ=0 and β=1 have optimal 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(13)Br1;r2=r,θ:r1<r<r2and0<θ<ω∩Ω,and Br1=B0;r1 in the polar coordinates r,θ and define cut-off function(14)ηρr=1,inB12ρR,δr,inB12ρR;ρR,0,inΩ\BρR,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) (i=1,2) are the solutions of homogeneous Stokes equation (3) with vanishing Dirichlet boundary condition on Γin. Also, λi has to be a positive real number and uis,pis∈H1+λiΩ×HλiΩ, i=1,2. 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 ηuis, i=1,2, vanish identically on ∂Ω. Then, in general, the solution u,p including singular parts of (3) can be rewritten in the form (11):(15)up=wq+α1η1u1sη1p1s+α2η2u2sη2p2s,where α1 and α2 are the stress intensity factors and w,q∈H2Ω×H1Ω.The strategy of FE-DSFM is to find the regular solutionw,q∈H2Ω×H1Ω 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 un+1=wn+1+αn+1ηρus and pn=qn+αnηρps. Then, Step 1 of the Uzawa algorithm can be(16)−μΔwn+1+αn+1ηρus+∇qn+αnηρps=f.Even though−μΔηρus and ∇ηρps∉L2, we know that −μ△ηρus+∇ηρps∈L2. So, (16) can be rewritten as(17)−μΔwn+1+∇qn+αn+1−μΔηρus+∇ηρps=f.For the sake of a clear explanation, we note that the inner product of vectorsa=a1,a2 and b=b1,b2 is(18)a,b=a1,b1+a2,b2,∇a,∇b=∂xa1,∂xb1+∂xa2,∂xb2+∂ya1,∂yb1+∂ya2,∂yb2.Then, the weak formulation of (17) is(19)μ∇wn+1,∇v+∇qn,v+αn+1−μΔηρus+∇ηρps,v=f,v,for all v∈H01Ω. Since qn∈H1Ω and −μ△ηρus+∇ηρps∈L2Ω, equation (19) is solvable, and it has unique solution wn+1∈H01Ω if we know αn+1∈ℝ. 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 cornerSince equation (19) is a Poisson-type problem to find wn+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):(20)sLdr,θ=r−π/ωsinπωθ.Lemma 3.
(see [19]). The function sLd in (20) has the following properties:(1)
sLd is not defined at origin(2)
sLd∈L2 and sLd∉H1(3)
−△sLd=0(4)
sLd=0on=Γin (except origin)(5)
For allu∈H01∩H2, −Δu,η2ρsLd=u,−Δη2ρsLd
We denote the vector form of (20), η2ρSLd=η2ρsLη2ρsL. We now test η2ρSLd to (17) and apply Lemma 3 to obtain(21)μwn+1,−Δη2ρSLd+∇qn,η2ρSLd+αn+1−μΔηρus+∇ηρps,η2ρSLd=f,η2ρSLd.
Letβs=−△ηρus+∇ηρps,η2ρSLd and βf=f,η2ρSLd−∇qn,η2ρSLd. Then, we have the system(22)μ∇wn+1,∇v+αn+1−μΔηρus+∇ηρps,v=f,v−∇qn,v,αn+1=βfβs−1βswn+1,−△η2ρSLd.We first check solvability of system (22).Lemma 4.
Define a mappingT:ℝ⟶H1 and F:H1⟶ℝ such that Tα=w, where(23)μ∇w,∇v+α−μΔηρus+∇ηρps,v=f,v−∇q,v,Fw=βfβs−1βsw,−Δη2ρSLd.
Given a functionf,∇qn∈L2,(24)F∘T=0.Therefore, there is a unique valueγ such that F∘Tγ∗=γ for all γ∗∈ℝ.Proof.
Letα1,α2∈ℝ with α1≠α2, Tα1=w1, and Tα2=w2. Define two equations:(25)μ∇a,∇v+∇qn,v=f,v,μ∇z,∇v+−μΔηρus+∇ηρps,v=0.
Then,Tα1=w1=a+α1z and Tα2=w1=a+α2z. Then,(26)F∘Tα2−F∘Tα1=1βsw1−w2,−Δη2ρSLd=1βsα1−α2z,−Δη2ρSLd=α1−α2βsz,−Δη2ρSLd.
We note that−μΔηρus+∇ηρps=0 in Ω\BρR and so is z with z=0 on ∂Ω. Since −Δη2ρSLd=0 in BρR, we have(27)z,−Δη2ρSLd=0.Theorem 1.
Formulation (22) has a unique solution w in H2Ω∩H01Ω.Proof.
By Lemma4, 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. LetT=K be a shape-regular quasi-uniform partition of Ω of mesh size h into closed elements K. The vector and scalar finite element spaces are(28)Wh≔wh∈L2Ω:wh|K∈PK,∀K∈T,Vh≔Wh∩H01Ω,ℙh≔qh∈L2Ω∩C0Ω:qh|K∈QK,∀K∈T,where PK and QK are spaces of polynomials with degree bounded uniformly with respect to K∈T. We stress that the space ℙh is composed of continuous functions to use integration by parts:(29)∇⋅vh,qh=−vh,∇qh,∀vh∈Vh,for all qh∈ℙh. Then, the finite element approximation for (22) is to find wh∈Vh and αhn+1∈ℝ such that(30)μ∇whn+1,∇vh+αhn+1−μΔηρus+∇ηρps,vh=f,vh−∇qhn,vh,αhn+1=βfhβs−1βswhn+1,−Δη2ρSLd,where βfh=f,η2ρSLd−∇qhn,η2ρSLd. And, the matrix form of the coupled system (30) becomes(31)A+a⋅bTwh=f,where A is symmetric positive definite square matrix and aandb are column vectors. System (31) can be solved by the Sherman–Morrison formula:(32)A+a⋅bT−1=A−1−A−1abTA−11+bTA−1a,and αhn+1 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 αhn instead of αhn+1 in first equation of (30). Algorithm 2 is the proposed method for one singular function.Algorithm 2: (FE-DSFM Uzawa iteration for one singular function). Start with any initial values ph0,αh0, and zh as the solution of μ∇zh,∇vh+−μΔηρus+∇ηρps,vh=0, for all vh∈Wh. Repeat the following steps until whn+1−whn≤tolerance:
Step 1: for all vh∈Wh, find w^hn+1∈Wh as the solution of μ∇w^hn+1,∇vh+∇qhn,∇vh+αhn−Δηρus+∇ηρps,vh=f,vhStep 2: find αhn+1∈ℝ from the equality αhn+1=βfh/βs−1/βsw^hn+1,−Δη2ρSLdStep 3: update whn+1∈Wh by whn+1=w^hn+1+αhn+1−αhnzStep 4: for all ph∈ℙh, update qhn+1 by qhn+1,ph=qhn,ph−μ∇⋅whn+1,ph−αhn+1∇⋅ηρus,phRemark 1.
Compared with the original Uzawa iteration (Algorithm1), the proposed method (Algorithm 2) needs only one more linear solver for the initial process.


## 4. Algorithm for Two Singular Functions

In this section, we construct an algorithm for two singularities in one corner. We consider the solution is(33)up=wq+α1ηρu1sηρp1s+α2ηρu2sηρp2s.With the implicit formation of two stress intensity factors, step 1 of Uzawa algorithm is(34)−μΔwn+1+∇qn+α1n+1−Δηρu1s+∇ηρp1s+α2n+1−Δηρu2s+∇ηρp2s=f.Here, the unknowns arewn+1, α1n+1, and α2n+1. However, our test function which is not in H01 is only η2ρSLd. Instead of finding another test function, we use the following property.Theorem 2.
(properties of singular functions). Letu1s=u1sv1s and u2s=u2sv2s be the singular functions of Stokes equation on same corner. Then,(35)−Δηρu1s+ηρp1sx,η2ρsLd−Δηρu2s+ηρp2sx,η2ρsLd+−Δηρv1s+ηρp1sy,η2ρsLd−Δηρv2s+ηρp2sy,η2ρsLd=0.Proof.
To tell the conclusion, it is too difficult to show equation (35) by algebraic calculation because the integral functions depend on some constants which are constructed by solutions of transcendental function.
First of all, from singular function (6) and cut-off function (14) with smooth function, we obtain(36)−Δηρus+∇ηρps=−C1∂2ηρ∂r2rλ+2λ+1∂ηρ∂rrλ−1λsinθsin1−λθsinλθ−λsinθcos1−λθ−C12λ∂ηρ∂rrλ−1cosθcos1−λθsinθcos1−λθ+C2∂2ηρ∂r2rλ+2λ+1∂ηρ∂rrλ−1sinλθ+λsinθcos1−λθλsinθsin1−λθ−C22λ∂ηρ∂rrλ−1cosθsin1−λθsinθsin1−λθ.
Let(37)UAλθ=−C1λsinθsin1−λθ+C2sinλθ+λsinθcos1−λθ,UBλθ=−C12λ+1λsinθsin1−λθ+2λcosθcos1−λθ+C22λ+1sinλθ+λsinθcos1−λθ−2λcosθsin1−λθ,VAλθ=−C1sinλθ−λsinθcos1−λθ+C2λsinθsin1−λθ,VBλθ=−C12λ+1sinλθ−λsinθcos1−λθ+2λsinθcos1−λθ+C22λ+1λsinθsin1−λθ−2λsinθsin1−λθ.
Then,(38)−Δηρus+∇ηρpsx=rληρrrUAλθ+rλ−1ηρrUBλθ,−Δηρvs+∇ηρpsy=rληρrrVAλθ+rλ−1ηρrVBλθ.
Since(39)−Δηρus+∇ηρpsx,η2ρsLd=∫∫0ω−Δηρus+∇ηρpsxη2ρr−π/ωsinπωθdθdr=∫∫0ωrληρrrUAλθ+rλ−1ηρrUBλθη2ρr−π/ωsinπωθdθdr=∫rληρrrη2ρr−π/ωdr∫0ωUAλθsinπωθdθ+∫rλ−1ηρrη2ρr−π/ωdr∫0ωUBλθsinπωθdθ,−Δηρvs+∇ηρpsy,η2ρsLd=∫rληρrrη2ρr−π/ωdr∫0ωVAλθsinπωθdθ+∫rλ−1ηρrη2ρr−π/ωdr∫0ωVBλθsinπωθdθ.
We can conclude that equation (35) holds if(40)∫0ωUAλ1θsinπωθdθ∫0ωUAλ2θsinπωθdθ+∫0ωVAλ1θsinπωθdθ∫0ωVAλ2θsinπωθ=0,(41)∫0ωUAλ1θsinπωθdθ∫0ωUBλ2θsinπωθdθ+∫0ωVAλ1θsinπωθdθ∫0ωVBλ2θsinπωθ=0,(42)∫0ωUBλ1θsinπωθdθ∫0ωUAλ2θsinπωθdθ+∫0ωVBλ1θsinπωθdθ∫0ωVAλ2θsinπωθ=0,(43)∫0ωUBλ1θsinπωθdθ∫0ωUBλ2θsinπωθdθ+∫0ωVBλ1θsinπωθdθ∫0ωVBλ2θsinπωθ=0.
However, we could not show equations (40)–(43) are valid by algebraically computation. Therefore, we calculate them numerically.
For numerical substitution, we first compute the above integral forms:(44)∫0ωUAλθ×sinπωθdθ=14ωλω−π−ωλω+πC1λ1+cosλω−C22+λsinλω+14ω2−λω+π−ω2−λω−π×C1λ1+cos2−λω−C2λsin2−λω,(45)∫0ωUBλθ×sinπωθdθ=14ωλω−π−ωλω+πC12λ2+3λ1+cosλω−C22λ2+7λ+2sinλω+14ω2−λω+π−ω2−λω−π×C12λ2−λ1+cos2−λω−C22λ2−λsin2−λω,(46)∫0ωVAλθ×sinπωθdθ=14ωλω−π−ωλω+πC12−λsinλω−λC21+cosλω+14ω2−λω+π−ω2−λω−π×C1λsin2−λω+C2λ1+cos2−λω,(47)∫0ωVBλθ×sinπωθdθ=14ωλω−π−ωλω+πC2λ−2λ21+cosλω−C12λ2−5λ−2sinλω+14ω2−λω+π−ω2−λω−π×C2λ−2λ21+cos2−λω−C1λ−2λ2sin2−λω.
We remark that the valuesλ1 and λ2 are from the equation(48)λ2sin2ω=sin2λω,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(i)
UAVA≔ left side of equation (40)(ii)
UAVB≔ left side of equation (41)(iii)
UBVA≔ left side of equation (42)(iv)
UBVB≔ left side of equation (43)
We choose three tolerances1e−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 UAVA, UAVB, UBVA, and UBVB. 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. The weak formulation of Step 1 in the Uzawa method test by η2ρSL is in the form(49)μwn+1,−△η2ρSLd+∇qn,η2ρSLd+α1n+1−Δηρu1s+∇ηρp1s,η2ρSLd+α2n+1−Δηρu2s+∇ηρp2s,η2ρSLd=f,η2ρSLd.
Letw=w1,w2 and f=f1,f2. Then, x component of (49) is(50)μw1n+1,−Δη2ρsLd+qnx,η2ρsLd+α1n+1−Δηρu1s+ηρp1sx,η2ρsLd+α2n+1−Δηρu2s+ηρp2sx,η2ρsLd=f1,η2ρsLd,and y component of (49) is(51)μw2n+1,−Δη2ρsLd+qny,η2ρsLd+α1n+1−Δηρv1s+ηρp1sy,η2ρsLd+α2n+1−Δηρv2s+ηρp2sy,η2ρsLd=f2,η2ρsLd.
If we use the following notations(52)Us1=−Δηρu1s+ηρp1sx,η2ρsLd,Us2=−Δηρu2s+ηρp2sx,η2ρsLd,Vs1=−Δηρv1s+ηρp1sy,η2ρsLd,Vs2=−Δηρv2s+ηρp2sy,η2ρsLd,βW1n+1=μw1n+1,−Δη2ρsLd+qnx,η2ρsLd,βW2n+1=μw2n+1,−Δη2ρsLd+qny,η2ρsLd,βf1=f1,η2ρsLd,βf2=f2,η2ρsLd,then (50) and (51) can be rewritten by(53)βW1n+1+α1n+1Us1+α2n+1Us2=βf1,βW2n+1+α1n+1Vs1+α2n+1Vs2=βf2.
Therefore, we obtain(54)α1n+1Us1⋅Us1=βf1−βW1n+1−α2n+1Us2Us1,(55)α1n+1Vs1⋅Vs1=βf2−βW2n+1−α2n+1Vs2Vs1.
Finally, combining (54) and (55) and equation (35) in Theorem 2 yields(56)α1n+1Us1⋅Us1+Vs1⋅Vs1=βf1−βW1n+1Us1+βf2−βW2n+1Vs1−α2n+1Us2Us1+Vs2Vs1=βf1−βW1n+1Us1+βf2−βW2n+1Vs1.
As a result, we can calculate one stress intensity factorα1n+1 by(57)α1n+1=βf1−βW1n+1Us1+βf2−βW2n+1Vs1Us1⋅Us1+Vs1⋅Vs1.
Similarly, another stress intensity factorα2n+1 can be computed by(58)α2n+1=βf1−βW1n+1Us2+βf2−βW2n+1Vs2Us2⋅Us2+Vs2⋅Vs2.Figure 2
UAVA, UAVB,UBVA, and UBVB values by numerical substitution with various tolerance of λ.
(a)(b)(c)(d)Remark 2.
Sinceηρus,−Δη2ρsLd=0, Lemma 4 is valid for two singular functions case.
Algorithm3 is the proposed Uzawa iteration for two singularities.Algorithm 3: (FE-DSFM Uzawa iteration for two singular functions). Start with any initial values ph0,α1h0, and α2h0. z1 and z2, for all vh∈Wh, is the solution of
μ∇z1,∇vh+−μΔηρu1s+∇ηρp1s,vh=0μ∇z2,∇vh+−μΔηρu2s+∇ηρp2s,vh=0Repeat the following steps untilwhn+1−whn≤tolerance:Step 1: for all vh∈Wh, find w^hn+1∈Whas the solution of μ∇w^hn+1,∇vh+∇qhn,∇vh+α1hn−Δηρu1s+∇ηρp1s,vh+α2hn−Δηρu2s+∇ηρp2s,vh=f,vhStep 2: by w^hn+1=W1n+1,W2n+1, find α1hn+1 and α2hn+1 from the equationsα1hn+1=βf1−βW1n+1Us1+βf2−βW2n+1Vs1/Us1⋅Us1+Vs1⋅Vs1α2n+1=βf1−βW1n+1Us2+βf2−βW2n+1Vs2/Us2⋅Us2+Vs2⋅Vs2Step 3: update whn+1∈Wh, by whn+1=w^hn+1+α1hn+1−α1hnz1+α2hn+1−α2hnz2Step 4: for all ph∈ℙh, update qhn+1 by qhn+1,ph=qhn,ph−μ∇⋅whn+1,ph−α1hn+1∇⋅ηρu1s,ph−α2hn+1∇⋅ηρu2s,ph


## 5. Numerical Experiments

In this section, we provide numerical simulations using the standard Uzawa method and Algorithm3. Our goal is to check their performance for the numerical solution uh and ph, which is composed of a regular solution wh and qh and stress intensity factors α1h and α2h. All numerical integration is calculated by the Gaussian quadrature 6 points. We notate Γout=r,θ∈∂Ω:θ≠0andθ≠ω.Example 1.
Simulation environments are(1)
PDE model: Stokes equation (3)(2)
Computational domain:−1,1×−1,1\0,1×−1,0∈ℝ×ℝ (Γ shape domain, ω=3π/2) (see Figure 3).(3)
Parameters:μ=1.0λ1=0.5444837367824639λ2=0.9085291898460987ρ=1.0 and radius R=3.0∈ℝ for cut-off functionun+1−un≤tolerance=1.0−8(4)
Exact solution (see Figures4 and 5):u=−sin2πxsin2πy+α1u1s+α2u2sv=sin2πxsin2πy+α1v1s+α2v2sp=2+cosπx2+cosπy−4+α1p1s+α2p2sα1=2.0,α2=−3.0(5)
Boundary condition: Dirichlet boundary condition.u=0 on Γin|θ=0 and Γin|θ=ω, and u≠0 on Γout.(6)
Finite elements: Taylor-Hood and minielements.
Minielement choosesVh and ℙh are the piecewise P1-bubble and piecewise P1 space pair satisfying discrete inf-sup condition. In this finite dimensional space, we desire u−uh0≤Ch2, Iu−uh∞≤Ch2, u−uh1≤Ch, and p−ph0≤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. The figure guarantees convergence rate of u and p. The 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.
The next finite element space is the Taylor-Hood space which choosesVh and ℙh as the piecewise P2 and piecewise P1 space pair. This also satisfies the discrete inf-sup condition. In this space, the desired convergence rates are u−uh0≤Ch3, Iu−uh∞≤Ch3, u−uh1≤Ch2, and p−ph0≤Ch2, 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.Figure 3
Computational domainh=1/8 of Example 1.Figure 4
Exact solutions,u (a) and v (b), of Example 1.
(a)(b)Figure 5
Exact solution,p (a) and the quiver of u (b), of Example 1.
(a)(b)Figure 6
Error decays between standard Uzawa method and Algorithm3 for Example 1. All algorithms use minielements.
(a)(b)Figure 7
Error decays between standard Uzawa method and Algorithm3 for Example 1. All algorithms use Taylor-Hood elements.
(a)(b)


## 6. Conclusions

We have presented a finite element using dual singular functions with Uzawa iteration for the steady Stokes equations including corner singularities. The 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.



---

*Source: 5592982-2021-06-09.xml*