Divergence-free meshless local Petrov–Galerkin method for Stokes flow

Abstract

The purpose of the present paper is development of an efficient meshless solution of steady incompressible Stokes flow problems with constant viscosity in two dimensions, with algebraic order of accuracy. This is achieved by employing a weak formulation with divergence-free matrix-valued quadratic Matérn (QM) radial basis function (RBF) for the shape function and divergence-free matrix-valued compactly supported Gaussian (CSG) RBF for the weight function on the computational domain and its boundary. The continuity equation is inherently built-in in the formulation and the pressure is eliminated from the formulation with the aid of divergence theorem and the choice of divergence-free weight function. The developed method is thus iteration free, and results in a banded system of equations to be solved jointly for both velocity components. Gauss–Legendre cell integration is performed in the current investigation. The characteristics of the method are assessed by changing its free parameters, i.e., weight functions’ sub-domain radius and shape functions’ support domain radius and the shape parameter. A sensitivity test for several choices of shape functions with regular centers arrangement is done to identify the appropriate support size for the shape and weight functions and stagnation errors are reported accordingly. To the best of our knowledge, this article is initiative in introducing the application of divergence-free MLPG method to incompressible flows, aiming at elimination of pressure from the governing equations in primitive variables, with the aid of divergence-free RBFs through weak formulation. Only the momentum equation needs to be solved. Hence, the formulation of the problem is much simpler than the building of divergence-free elements in the related mesh-based methods.

Appendices

Appendix A: Tensor–vector identities

Assume vectors $\vec {V}, \vec {W}: {\mathbb {R}}^2 \rightarrow {\mathbb {R}}^2$ defined as

$$\begin{aligned} {\vec {V}} = (V_x,V_y)^T, \, {\vec {W}} = (W_x,W_y)^T, \end{aligned}$$

(5.1)

with $V_{x,x} = \displaystyle {\frac{\partial V_{x}}{\partial x}}$, $V_{x,y} = \displaystyle {\frac{\partial V_{x}}{\partial y}}$, $V_{y,x} = \displaystyle {\frac{\partial V_{y}}{\partial x}}$, $V_{y,y} = \displaystyle {\frac{\partial V_{y}}{\partial y}}$, $V_{x,xx} = \displaystyle {\frac{\partial ^2 V_{x}}{\partial x^2}}$, $V_{x,xx} = \displaystyle {\frac{\partial ^2 V_{x}}{\partial x^2}}$, $V_{x,yy} = \displaystyle {\frac{\partial ^2 V_{x}}{\partial y^2}}$, $V_{y,xx} = \displaystyle {\frac{\partial ^2 V_{y}}{\partial x^2}}$, $V_{y,yy} = \displaystyle {\frac{\partial ^2 V_{y}}{\partial y^2}}$, and also

$$\begin{aligned} (\cdot )_{,x} = \displaystyle {\frac{\partial (\cdot )}{\partial x}},\qquad (\cdot )_{,y} = \displaystyle {\frac{\partial (\cdot )}{\partial y}}, \end{aligned}$$

then the following relations hold

$$\begin{aligned} \nabla {\vec {V}}= & {} \begin{bmatrix} V_{x,x} &{} V_{y,x} \\ V_{x,y} &{} V_{y,y} \\ \end{bmatrix}, \end{aligned}$$

(5.2)

$$\begin{aligned} \Delta {\vec {V}}= & {} \nabla \cdot \left( {\nabla {\vec {V}}}\right) = \begin{bmatrix} V_{x,xx} + V_{x,yy}\\ V_{y,xx} + V_{y,yy}\\ \end{bmatrix}, \end{aligned}$$

(5.3)

$$\begin{aligned} \nabla \vec {V} \cdot \vec {W}= & {} \begin{bmatrix} V_{x,x} \, W_{x} + V_{y,x} \, W_{y}\\ V_{x,y} \, W_{x} + V_{y,y} \, W_{y}\\ \end{bmatrix}, \end{aligned}$$

(5.4)

$$\begin{aligned} \nabla \cdot \left( {\nabla \vec {V} \cdot \vec {W}}\right)= & {} \left( {V_{x,x} \, W_x + V_{y,x} \, W_y }\right) _{,x} \nonumber \\&+ \left( {V_{x,y} \, W_x + V_{y,y} \, W_y }\right) _{,y} \nonumber \\= & {} V_{x,xx} \, W_x + V_{x,x} \, W_{x,x} \nonumber \\&+ V_{y,xx} \, W_y + V_{y,x} \, W_{y,x} \nonumber \\&+ V_{x,yy} \, W_x + V_{x,y} \, W_{x,y} + V_{y,yy} \, W_y\nonumber \\&+ V_{y,y} \, W_{y,y}, \end{aligned}$$

(5.5)

$$\begin{aligned} \Delta \vec {V} \cdot \vec {W}= & {} \left( {V_{x,xx} + V_{x,yy}}\right) W_x \nonumber \\&+ \left( {V_{y,xx} + V_{y,yy}}\right) W_y \nonumber \\= & {} V_{x,xx} \, W_x + V_{x,yy} \, W_{x} \nonumber \\&+ V_{y,xx} \, W_y + V_{y,yy} \, W_{y}, \end{aligned}$$

(5.6)

$$\begin{aligned} \nabla \vec {V} : \nabla \vec {W}= & {} V_{x,x} \, W_{x,x} + V_{y,x} \, W_{y,x} \nonumber \\&+ V_{x,y} \, W_{x,y} + V_{y,y} \, W_{y,y}. \end{aligned}$$

(5.7)

According to (5.5), (5.6) and (5.7), Eq. (3.13) is trivial.

Appendix B: Analytical solution at the outlet for the laminar flow between parallel plates

Assume $\Omega = [0,L]\times [0,h]$. The Stokes Eq. (2.4) holds at the outlet, hence,

$$\begin{aligned} \begin{array}{c} V_{x,x} + V_{y,y} = 0,\\ \quad \; V_{x,xx} + V_{x,yy} = P_{,x},\\ \quad \; V_{y,xx} + V_{y,yy} = P_{,y}. \end{array} \end{aligned}$$

(5.8)

On the other hand, due to the flow property at the outlet,

$$\begin{aligned} V_y = 0. \end{aligned}$$

(5.9)

The first result is obtained from the continuity equation

$$\begin{aligned} V_{x,x} = 0, \end{aligned}$$

(5.10)

next, is due to the momentum equation that the pressure is not a function of the independent variable y, i.e., $P_{,y} = 0$. Then, $V_{x,x} = 0$ requires that

$$\begin{aligned} V_{x,yy} = P_{,x}. \end{aligned}$$

(5.11)

Two times integration of (5.11) with respect to y gives

$$\begin{aligned} V_x(L,y) = \displaystyle \frac{1}{2} P_{,x}y^2 + c_1 y + c_2, \end{aligned}$$

(5.12)

with $c_1$ and $c_2$ as the integration constants that are identified with applying the no-slip boundary conditions on the walls:

• at $y = 0$: $V_x = 0\quad \rightarrow \quad c_2 = 0$.

• at $y = h$: $V_x = 0\quad \rightarrow \quad \displaystyle \frac{h^2}{2}P_{,x} + c_1 h = 0 \quad \rightarrow \quad c_1 = -\displaystyle \frac{h}{2}P_{,x}.$

Substitution of $c_1$ and $c_2$ in (5.12), results in

$$\begin{aligned} V_{,x}(L,y) = \displaystyle \frac{1}{2}P_{,x}\left( {y^2 - hy}\right) . \end{aligned}$$

(5.13)

The volume flow rate at the outlet cross section is equal to the integration of velocity profile over the width [0, h]

$$\begin{aligned} Q= & {} \int _A{\vec {V}.\vec {n}\,dA} = \int _0^h{V_x\,dy} = \int _0^h{\displaystyle \frac{1}{2}P_{,x}\left( {y^2 - hy}\right) \,dy}\nonumber \\= & {} -\displaystyle \frac{1}{12}\,P_{,x}\,h^3. \end{aligned}$$

(5.14)

Considering the volume flow rate for the average velocity, i.e., the total flow rate per unit that is identical to $u_{\mathrm{in}}$

$$\begin{aligned} Q = h\, u_{\mathrm{in}}, \end{aligned}$$

(5.15)

and finally, comparing the right hand side of Eqs. (5.14) and (5.15), gives

$$\begin{aligned} u_{\mathrm{in}} = -\displaystyle \frac{h^2}{12}\,P_{,x}. \end{aligned}$$

(5.16)

Consequently, $V_x(L,y) = 6\, u_{\mathrm{in}}\left( {\displaystyle \frac{hy - y^2}{h^2}}\right)$.

Appendix C: A note on the entrance length

In our case of study, i.e., the Stokes flow, a simple approach to determine the entrance length is considering the parabolic shape of the fully developed velocity vector at the outlet. According to the analytical solution, the x-component of the velocity is of the form

$$\begin{aligned} V_x(L, y) = 6\, u_{\mathrm{in}}\left( {\displaystyle \frac{hy - y^2}{h^2}}\right) , \end{aligned}$$

with L denoting the channel length. $V_x(L,y)$ reaches it’s maximum value $\displaystyle \frac{3}{2}u_{\mathrm{in}}$ at $y = \displaystyle \frac{h}{2}$. As a result, the value of $|V_x(L, \frac{h}{2}) - \frac{3}{2}u_{\mathrm{in}}|$ is a good indicator whether we have considered the sufficient channel width. Please note that the relative error in maximum norm at the outlet confirms that the chosen length of the computational domain is sufficient.