Analysis of the Element-Free Galerkin Method with Penalty for Stokes Problems

The element-free Galerkin (EFG) method with penalty for Stokes problems is proposed and analyzed in this work. A priori error estimates of the penalty method, which is used to deal with Dirichlet boundary conditions, are derived to illustrate its validity in a continuous sense. Based on a feasible assumption, it is proved that there is a unique weak solution in the modified weak form of penalized Stokes problems. Then, the error bounds with the penalty factor for the EFG discretization are derived, which provide a rationale for choosing an efficient penalty factor. Numerical examples are given to confirm the theoretical results.


Introduction
It is well-known that Stokes equations describe low Reynolds number flow motion and play a fundamental role in the numerical modeling of incompressible viscous flows. Recently, there has been an increasing interest in solving Stokes problems by various meshfree (or meshless) methods [1][2][3][4][5] to alleviate mesh-related dilemmas, including some collocation meshless methods, such as virtual interpolation point method [6], generalized finite difference method [7], divergence-free kernel approximation method [8], as well as some Galerkin meshless methods, such as the moving least square reproducing kernel method [9], weighted extended B-spline method [10], Galerkin boundary node method [11], and the divergence-free meshless local Petrov-Galerkin method [12].
The element-free Galerkin (EFG) method [13] is a Galerkin-based meshfree discretization technique for solving partial differential equations. The trial and test functions for the EFG method are generated by the moving least squares (MLS) approximation [14]. During the past several decades, many research works have been devoted to improving and extending the MLS approximation, see [4,[15][16][17][18] for various details. To offset the lack of interpolating properties of the MLS shape functions, several interpolation-type MLS methods have been developed. We refer to [14,[18][19][20][21] and the references therein for details.
In addition to choosing the interpolation-type MLS methods, some mandatory methods, such as the Lagrange multiplier method [3,4,13], Nitsche method [22,23] and penalty method [3,4,[24][25][26][27][28] are desirable in practical applications. The important feature of these methods is that they can straightforwardly use the non-interpolating trial and test functions by modifying the traditional weak form. The penalty method seems to be more appealing because of its ease of implementation, its variable-preservation and its general framework, and these significant advantages enable numerical analysis.
In the context of the EFG method, many papers are devoted to penalty-based error analysis for elliptic problems [24,25], parabolic problems [26,28] and contact problems [27]. To the authors' knowledge, for Stokes problems, a priori errors of the meshless penalty method have not been explained, and numerical analysis with a penalty factor has not been presented either. The main difficulty may be that the modified weak form based on the penalty method is different from the standard weak form, thus the standard meshless Galerkin procedure cannot be used directly.
In order to better clarify the principle of the penalty method and determine an efficient penalty factor, the presented paper is an extension of the previous works [25,28] on the EFG method for Stokes problems. The modified weak form of penalized Stokes problems is analyzed. Based on a discrete inf-sup condition, error bounds with a penalty factor of the EFG discretization are given in H 1 norm for velocity approximation and in L 2 norm for pressure approximation, respectively. Furthermore, an error estimate with a penalty factor for velocity approximation in L 2 norm is also derived. Numerical examples are given to confirm the theoretical results.
The remainder of the paper is organized as follows. In Section 2, we introduce the Stokes problem and its standard weak formulation. Then, a priori estimates of the penalty method are determined on the Dirichlet boundary and in the problem domain in Section 3, respectively. Sections 4 and 5 present the modified weak form of the penalized Stokes problem and the EFG numerical discretization, respectively. Section 6 is devoted to the error analysis for velocity and pressure approximations. Numerical examples are presented in Section 7 and conclusions are drawn in the final Section.
Throughout this paper, the letter C, with a superscript or subscript, is used to represent a generic positive constant, independent of the characteristic distance h and could take different values at different appearances.

Penalized Stokes Problem
In the subsequent numerical discrete process, since the MLS shape functions with noninterpolating properties will be adopted, the penalty method is used to impose the Dirichlet boundary condition. In order to better illustrate the principle of the penalty method, by using a penalty factor α, (1) is approximated as the following penalized problem: in Ω, where n is the unit outward normal to Γ. By combining (1) and (5), an a priori estimate of the penalty method on the boundary Γ is first derived.
It is shown by Lemma 1 that when the penalty factor α approaches infinity, the solution u α of (5) tends to 0 with L 2 norm on the Dirichlet boundary Γ. Clearly, the boundary term pn is almost independent of the penalty factor on this point. In addition, −ν∆u α + ∇p = f is equivalent to −ν∆u α = f − ∇p, which can be regarded as the approximation of −ν∆u = f − ∇p. Then, the a priori estimate within the problem domain is exported.
It can be found that when the penalty factor α approaches infinity, the solution u α tends to u with H 1 semi-norm in the problem domain by Lemma 2. Lemmas 1 and 2 fully demonstrate the validity of the penalty method in a continuous sense.
A discrete assumption similar to (11) is used in the Nitsche method to ensure the continuity and coercivity of the bilinear form for the elliptic equation [1]. In this paper, the assumption (11) is proposed to certify the continuity of a α (·, ·) and inf-sup condition, thus deriving that A α ((·, ·); (·, ·)) is continuous and coercive.

EFG for Penalized Stokes Problem
To approximate the solution of the modified weak form (15), the approximate space of the velocity is defined as: (18) and the approximate space of the pressure is: is a set of N 2 pressure nodes. Φ i (x) and Ψ j (x) represent the MLS shape functions based on velocity nodes and pressure nodes, respectively. Now, we briefly state the MLS shape function and its approximation error by taking the velocity nodes as an example, which is similar to the pressure nodes. The MLS shape functions Φ i (x) are: in which p j (x) denotes the shifted and scaled monomial basis function [22,24,25,28], ∧(x) ={λ 1 , λ 2 , · · · , λ n x } means the global sequence numbers of nodes whose support domains cover the point x. The support domain of Assume that the configuration of velocity nodes {x i } N 1 i=1 satisfies the following conditions: (B1) Define the characteristic distance h as wherem represents the largest degree of the used monomial basis functions.
Moreover, assume that derivatives of the weight function w i (x) up to order γ are bounded and continuous such that

Lemma 4 ([24]
). Assume that w ∈ Hm +1 (Ω), conditions (B1) and (B2) are satisfied, Mw denotes the MLS approximation of w. Then The following lemma is regarded as a sufficient condition for (X h , M h ) to satisfy the discrete inf-sup condition in the meshless method.
The EFG solutions u h and p h have an estimate similar to Lemma 1.
According to Lemmas 1 and 2, theoretically, the penalty method requires that the penalty factor α tends to infinity to impose the Dirichlet boundary condition. Nevertheless, in numerical calculations, the coefficient matrix of the resulting system will become ill-conditioned when the penalty factor increases uncontrollably. By deploying α = Ch (−2m−1)/3 in (25) and (26), the optimal convergence rates are derived as: When the linear basis function is chosen in the MLS approximation, i.e.,m = 1 and the penalty factor α = Ch −1 , the corresponding convergence rates are optimal as: Clearly, in this case, the EFG solution u h converges to the exact solution u with optimal convergence rate h in H 1 (Ω), but the pressure numerical solution p h only maintains first order convergence in L 2 (Ω).
To obtain the numerical error of velocity u in L 2 norm, the following definition and lemma are required.
The following approximate error follows from the above definition. Lemma 7 ([25,31]). Let w ∈ H l (Ω) and w = 0 on Γ. If Γ is sufficiently smooth, k 0 ≥ l ≥ 2 and k 1 ≥ 1, then there exists a function ξ ∈ γ k 0 ,k 1 h 0 (Ω) such that: Clearly, the MLS shape functions satisfy the requirements of ξ in Definition 1. There- (Ω) with k 0 ≥m + 1 and k 1 ≥ 1. Now, with the aid of the duality argument, an error bound of u in the L 2 norm can be derived.

Numerical Examples
This section presents four numerical examples to illustrate the theoretical error results proposed in the previous section. The problem domain is the unit square Ω = [0, 1] × [0, 1]. An efficient discrete node configuration for velocity and pressure has been proposed to satisfy the condition of Lemma 5 [9,10], see Figure 1.

Example 1
The first example is the Stokes problem (1) with the viscosity ν = 1. The exact solution is: Figure 2 depicts the log-log plots of the errors u − u h L 2 (Ω) and u − u h H 1 (Ω) with respect to increasing penalty factors α = 10, 10 2 , 10 3 , · · · , 10 9 , 10 10 for linear basis function (m = 1). The radius of support domain of the velocity node is 1.5h and four types of equidistant nodes 11 × 11, 21 × 21, 41 × 41 and 81 × 81 are used. Clearly, a too-small or too-big penalty factor increases the numerical errors and leads to the invalidation of numerical calculations. Theorems 1 and 2 imply that the theoretical errors of velocity are dominated by α −1 for a small penalty factor. It can be observed that the numerical errors of velocity keep almost the same convergence order α −1 from the left sides of Figure 2a,b. Obviously, the numerical errors of velocity agree well with the theoretical error estimates. The condition numbers of the discrete coefficient matrix for the increasing penalty factors are shown in Figure 3. It is clear that the condition numbers increase with the increase of the penalty factor and the condition number is approximately α 2 . Therefore, a too-big penalty factor predestinates invalidate the penalty method, which in turn leads to the failure of the numerical methods using the penalty method.     From the point of view of the numerical results above, a suitably large constant penalty factor can obtain stable numerical solutions. On the other hand, the latest theoretical analysis [24,25,28] implies that the penalty factor affects the convergence order of the numerical solutions. According to Theorem 1, an option is α = C s h −1 for linear basis function, where C s is a constant related to the problem to be solved. Figure 5 shows the log-log plots of the errors u − u h H 1 (Ω) and p − p h L 2 (Ω) with respect to the nodal spacing h for different C s . Clearly, C s affects the accuracy of the numerical solutions, but hardly impacts the convergence order. By comparison, a great choice is α = 1000h −1 from a precision point of view, and the error bounds have been tabulated in Table 1. It is clear that the numerical convergence orders are consistent with the theoretical analysis.
Moreover, it can be known from Theorem 2 that a valid choice is α = C s h −5/3 for the L 2 norm of velocity errors in terms of linear basis function. Figure 6 displays the log-log plots of the errors u − u h L 2 (Ω) with respect to the nodal spacing h for α = C s h −1 and α = C s h −5/3 . Similarly, α = 1000h −5/3 is an excellent option. Meanwhile, the numerical errors of α = 1000h −1 and α = 1000h −5/3 have been shown in Table 2. Visibly, the numerical convergence order of velocity is 1/3 order higher than the theoretical result in terms of L 2 norm, but the numerical errors of α = 1000h −1 still accord with the theoretical analysis.

Example 3
The third example considers Kovasznay flow [6]. The exact solutions are: where λ = Re 2 − 4π 2 + Re 2 4 1/2 and Re = 1 ν . The EFG numerical solutions u 1h , u 2h and p h are shown in Figure 8 for Re = 40 and Figure 9 for Re = 200 using uniform 41 × 41 velocity nodes and α = 1000h −1 . Again, the EFG method gains very accurate numerical solutions. Tables 4 and 5 give the errors for Re = 40 and Re = 200, respectively. Obviously, except that the numerical convergence order of velocity of α = 1000h −5/3 is second order in L 2 norm, the numerical convergence orders are still in good agreement with the theoretical analysis for α = 1000h −1 .

Example 4
The last example considers the lid-driven cavity flow problem, which is often regarded as a typical benchmark for incompressible flows. The body force f = 0. Figure 10 shows boundary conditions. On the top side u = (1, 0) T is given, and other sides are no-slip. Figure 11 shows the EFG solutions of velocities u 1 and u 2 along the vertical centerline x 1 = 0.5 and horizontal centerline x 2 = 0.5, respectively. The numerical results are derived by using uniform 81 × 81 velocity nodes and α = 1000h −1 . Meanwhile, the results of the Galerkin boundary node method (GBNM) [11] are also plotted in this figure for comparison. Clearly, the EFG solutions are in good agreement with the GBNM results. Moreover, Figure 12 displays the computed plots of streamline, vorticity contour and pressure contour. It can be found that stable numerical results of velocity and pressure are achieved.

Conclusions
In this paper, we presented and analyzed a penalty-based EFG method for Stokes problems. The penalty method allows the use of the MLS approximation to generate trial and test functions in the modified weak form. A priori errors of the penalty method are determined on the Dirichlet boundary and in the problem domain respectively, which state the feasibility of the penalty method in a continuous sense. For the penalized Stokes problems, the existence and uniqueness of the weak solution are proved under a rational assumption, which provides a valid foundation for the EFG numerical discretization. Under the condition of discrete inf-sup, error estimates with a penalty factor are provided in H 1 and L 2 norms for velocity approximation and in L 2 norm for pressure approximation.
Numerical results reveal that the proposed method exhibits good numerical accuracy and agrees well with the theoretical prediction. Note that for linear basis functions, the numerical convergence order of velocity can reach the second order in L 2 norm, but we have only theoretically obtained a suboptimal convergence order of velocity. Therefore, more research is required on the theoretical analysis of the present method for deriving the optimal convergence order of velocity in L 2 norm. In addition, how to reduce the condition numbers is an important research topic.