A multiple-scale Pascal polynomial for 2D Stokes and inverse Cauchy–Stokes problems

Accepted Manuscript A multiple-scale Pascal polynomial for 2D Stokes and inverse Cauchy–Stokes problems

Chein-Shan Liu, Der-Liang Young

PII: DOI: Reference:

S0021-9991(16)00071-1 http://dx.doi.org/10.1016/j.jcp.2016.02.017 YJCPH 6392

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

18 May 2015 1 February 2016 4 February 2016

Please cite this article in press as: C.-S. Liu, D.-L. Young, A multiple-scale Pascal polynomial for 2D Stokes and inverse Cauchy–Stokes problems, J. Comput. Phys. (2016), http://dx.doi.org/10.1016/j.jcp.2016.02.017

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Highlights • • • •

We propose a third-first order system and a third-third order system to solve direct and the inverse Cauchy problems in Stokes flows. A multiple-scale Pascal polynomial method is developed. The scales are determined by the collocation points. The multiple-scale Pascal polynomial expansion method is accurate and stable against large noise.

A multiple-scale Pascal polynomial for 2D Stokes and inverse Cauchy-Stokes problems Chein-Shan Liua,b , Der-Liang Younga,∗ a

Department of Civil Engineering, National Taiwan University, Taipei 106-17, Taiwan b

Center for Numerical Simulation Software in Engineering and Sciences,

College of Mechanics and Materials, Hohai University, Nanjing, Jiangsu 210098, China ∗

Corresponding author, E-mail: [email protected]

Abstract The polynomial expansion method is a useful tool for solving both the direct and inverse Stokes problems, which together with the pointwise collocation technique is easy to derive the algebraic equations for satisfying the Stokes differential equations and the specified boundary conditions. In this paper we propose two novel numerical algorithms, based on a third-first order system and a third-third order system, to solve the direct and the inverse Cauchy problems in Stokes flows by developing a multiple-scale Pascal polynomial method, of which the scales are determined a priori by the collocation points. To assess the performance through numerical experiments, we find that the multiple-scale Pascal polynomial expansion method (MSPEM) is accurate and stable against large noise. Keywords: Stokes problem, Inverse Cauchy-Stokes problem, Meshless method, Pascal triangle, Multiple-scale Pascal polynomial expansion method 1. Introduction The flows with very low Reynolds number form a class of incompressible flows, which have the common physical events described as follows. For a domain characterized by a length L and a motion characterized by the velocity U , the Reynolds number Re = LU ρ/μ is small either because the fluid is very viscous, or because the inertia or density is very small. The viscous incompressible flow is governed by the nonlinear Navier-Stokes partial differential equations (PDEs), of which the closed-form solutions are usually not available if without making any simplification. Under the assumption of negligible inertia force, the NavierStokes equations can be reduced to the Stokes equations. Although the governing equations can be simplified to linear ones, it is still very difficult to obtain the analytical solutions. In the past several decades, the numerical simulations are implemented for understanding the Stokes fluid behavior, which can be classified into four types of formulations: the primitivevariable formulation [22, 25], the velocity-vorticity formulation [9, 26], the stream functionvorticity formulation and the stream function formulation [24]. Among them, the primitivevariable formulation can directly acquire the pressure distribution. 1

What causes the situation of one flow different from another is the boundary conditions, which include the location and motion of walls, imposed pressure difference, and prescribed velocities. The formulation of boundary conditions follows a few rules, for example, the nonslip condition, but in general the physical intuition is used to make reasonable assumptions. When the velocity-vorticity formulation and the stream function-vorticity formulation are considered, the required boundary conditions are quite complex and not easy to obtain. Therefore, it is very important to develop a better numerical algorithm for directly solving the primitive-variable formulation of the Stokes equations. For this purpose, many numerical methods were employed to solve the Stokes flow, such as the finite difference method [22], the finite element method, and the boundary element method [9, 27]. They belong to the mesh-dependent numerical methods. In their numerical implementations a time-consuming mesh generation and a troublesome numerical quadrature are needed. In contrast to those mesh-dependent methods, there are numerical algorithms which are classified as the meshindependent methods. The most well-known meshless methods are the radial basis function (RBF) collocation method [23, 5], the method of fundamental solutions (MFS) [25, 26, 24, 3], and the Trefftz method [10]. The polynomial expansion method is easily implemented; however, it is seldom used in the numerical solution of Stokes problems. When the MFS and the Trefftz method are boundary-type meshless methods, both the RBF collocation method and the polynomial expansion method are of domain-type meshless methods. In any way, when the boundary conditions are specified the information has to come from measurements. For realistic applications of Stokes equations, some boundary conditions cannot be directly measured, for example the pressure, or sometimes it is estimated before numerical solution. Hence, it will result in one tpye of the inverse Stokes problems. This sort problem is known to be severely ill-posed since the solution does not exist for any pair of data, and even if such a solution exists, it does not always depend continuously on the prescribed data. Regarding the solution of inverse Cauchy-Stokes problem, the related papers are a few, of which we may mention the works [28, 29, 6, 4, 1, 2, 11, 12]. The velocity-vorticity formulation of the Stokes flow is well-known to be [23]: Δω = 0, Δv = −∇ × ω,

(1) (2)

where v and ω := ∇ × v are, respectively, the velocity and vorticity vectors, and Δ is the Laplacian operator. To solve the above equations pair one needs to impose the boundary conditions of v and ω. However, it does not guarantee that the solutions satisfy the continuity equation: ∇ · v = 0. (3) Inserting ω = ∇ × v into Eq. (1) we can obtain Δ(∇ × v) = 0,

(4)

which are third-order PDEs. In this paper we will employ the polynomial expansion method to solve Eqs. (3) and (4) directly by imposing the boundary conditions of v only, which constitute a third-first order system. Eqs. (3) and (4) are somewhat non-symmetric that when the first equation is a first-order PDE, the latter are third-order PDEs. Taking a divergence of Eq. (2) it follows that Δ(∇ · v) = 0, 2

(5)

which is again a third-order PDE. According to the minimum and maximum theorem of Laplace equation, when we impose the boundary condition ∇· v = 0 on the whole boundary, Eq. (5) guarantees that the continuity equation (3) is automatically satisfied in the whole domain. Now Eqs. (5) and (4) are more symmetric than the pair (3) and (4). Alternatively for the second option we can solve Eqs. (5) and (4) under the specified boundary conditions of v and an extra condition ∇ · v = 0 on the whole boundary, which is a third-third order system. The polynomial interpolation is an ill-posed problem and it makes the interpolation by higher-order polynomials not being easy to be numerically implemented. In order to overcome these difficulties, Liu and Atluri [20] have introduced a characteristic length into the high-order polynomials expansion, which improved the numerical accuracy for the applications to solve the ill-conditioned linear interpolation problems. At the same time, Liu et al. [21] have developed a multi-scale Trefftz-collocation Laplacian conditioner to deal with the ill-conditioned linear systems. The concept of multi-scale Trefftz-collocation method has been employed by Chen et al. [7] to solve the sloshing wave problem. Liu [15] has proposed a multi-scale half-order polynomial interpolation method, and Kuo et al. [14] have used the modified two characteristic lengths Pascal triangle to solve the inverse heat source problem. In this paper we propose a new multiple-scale expansion technique by permitting the expansion through the use of higher-order polynomials, which can overcome the above-mentioned ill-conditioned behavior. The purpose of this paper is to develop a powerful polynomial expansion method, endowing with the advantage of easy numerical implementation and having a great flexibility applied to the Stokes problem and the inverse Cauchy-Stokes problem defined in an arbitrary plane domain. This paper is arranged as follows. In Section 2 we introduce a simple modification of the Pascal polynomial expansion method by considering multiple scales, which according to the concept of equilibrated matrix can be fully determined by the collocation points. The main results are shown in Section 3, where two novel systems of third-first order system and third-third order system are proposed to solve the two-dimensional Stokes equations. The numerical examples for direct Stokes problems are solved in Section 4, while the inverse Cauchy-Stokes problems are solved in Section 5. Finally, we draw some conclusions in Section 6. 2. A multiple-scale Pascal polynomial The use of polynomial expansion as a trial solution of linear PDEs is simple, which can be used straightforward to derive the required algebraic equations to determine the expansion coefficients after a suitable collocation of spatial points in the problem domain. However, it is seldom used as a major numerical tool to solve the linear PDEs. The main difficulty is that the resultant linear algebraic equations (LAEs) are often highly ill-conditioned. How to reduce the condition number of the linear system becomes an important issue in the application of polynomials expansion method to solve the linear PDEs.

3

The elements in following polynomial matrix: ⎡ 1 y y 2 · · · y m−1 ym 2 m−1 ⎢ x xy xy · · · xy xy m ⎢ 2 2 2 2 2 m−1 ⎢ x x y x y ··· x y x2 y m ⎢ ⎢ .. .. .. .. .. ⎣ . . . ··· . . xm xm y xm y 2 · · · xm y m−1 xm y m

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

(6)

are often used to expand the solution of u(x, y). If the elements are restricted in the left-upper triangle, then such an expansion is known as the Pascal polynomial expansion: 1 x y x2 xy y 2 x3 x2 y xy 2 y 3 4 x x3 y x2 y 2 xy 3 y 4 x5 x4 y x3 y 2 x2 y 3 xy 4 y 5 ... ... ... ... ... ... ... ...

(7)

with u(x, y) being expanded by u(x, y) =

i m  

cij xi−j y j−1 ,

(8)

i=1 j=1

where the coefficients cij are to be determined. The total number of all elements is n = m(m + 1)/2, and the highest order of the above polynomial is m − 1. Because x and y in the problem domain Ω may be arbitrarily large quantity, the above expansion would lead to a divergence of the powers xi−j y j−1 . In order to obtain accurate solution of the inverse Cauchy-Stokes problem by using the modified Pascal polynomial expansion method, we have to develop a more effective and accurate solution method to solve the resulting LAEs by significantly reducing the condition number. For the purpose of demonstration, we begin with the following linear Poisson equation: Δu(x, y) = F (x, y), (x, y) ∈ Ω, u(x, y) = H(x, y), (x, y) ∈ Γ,

(9) (10)

where Δ is the Laplacian operator, Γ is the boundary of the problem domain Ω, and F and H are given functions. We consider a novel multiple-scale Pascal polynomial expansion of u(x, y) by i m   u(x, y) = cij sij xi−j y j−1 , (11) i=1 j=1

where the scales sij will be determined below. From Eq. (11) it is straightforward to write Δu(x, y) =

i m  

cij sij (i − j)(i − j − 1)xi−j−2 y j−1

i=1 j=1

+(j − 1)(j − 2)xi−j y j−3 . 4

(12)

Inserting it into Eqs. (9) and (10), and selecting n1 and n2 collocation points on the boundary and in the domain, to satisfy the boundary condition and the field equation, respectively, we can obtain a system of linear algebraic equations (LAEs) to solve the n coefficients cij . Now we use the above Poisson problem to demonstrate how to determine sij . First set sij = 1 and we can obtain a set of linear equations. It is convenient to express the resulting LAEs in terms of a matrix-vector product form: Ac = b.

(13)

Usually, Eq. (13) is an over-determined system for that we may collocate more points to produce more equations with number nc , which are used to find n coefficients in c with n  nc . Here the dimension of A is nc × n. First the coefficients cij used in the expansion (8) can be expressed as an n-dimensional vector c with the components ck , k = 1, . . . , n generated by k=0 Do i = 1, m Do j = 1, i k =k+1 ck = cij Enddo.

(14)

¯ the term u(x, y) can be expressed as an inner product of Then for a generic point (x, y) ∈ Ω a vector a with c, i.e., ⎤ ⎡ c1 ⎢ c2 ⎥ ⎥ ⎢

⎥ ⎢ u(x, y) = 1 x y x2 xy y 2 x3 x2 y xy 2 y 3 . . . ⎢ c3 ⎥ ⎢ .. ⎥ ⎣ . ⎦ cn = aT c.

(15)

Similarly, for a generic point (x, y) ∈ Ω the term Δu(x, y) can be expressed as an inner product of a vector d with c, where the components dk are generated by k=0 Do i = 1, m Do j = 1, i k =k+1 dk = (i − j)(i − j − 1)xi−j−2 y j−1 + (j − 1)(j − 2)xi−j y j−3 Enddo.

(16)

Then, when we select n1 points (xi , yi ), i = 1, . . . , n1 on the boundary Γ to satisfy the boundary condition, and n2 points (xi , yi ), i = 1, . . . , n2 on the domain Ω to satisfy the field

5

equation, we have

⎡ ⎤ H(x1 , y1 ) aT 1 ⎢ .. ⎢ .. ⎥ ⎢ . ⎢ . ⎥ ⎢ ⎢ T ⎥ ⎢ H(x n1 , y n1 ) ⎢ a n1 ⎥ ⎥, b = ⎢ A=⎢ ⎢ ⎢ dT ⎥ ⎢ F (xn1 +1 , yn1 +1 ) ⎢ 1 ⎥ ⎢ ⎢ . ⎥ .. ⎢ ⎣ .. ⎦ . ⎣ T F (x n +n dn2 1 2 , yn1 +n2 ) ⎡

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎥ ⎦

(17)

where n1 + n2 = nc . Instead of Eq. (13), we can solve a normal linear system: AT Ac = AT b,

(18)

where we will use the algorithm of conjugate gradient method (CGM) to solve Eq. (18). Indeed the scaling of linear algebraic equations is an important topic that has a long history of development. A matrix is called equilibrated if all its rows or columns have the same norm, and under this condition the matrix is better conditioned than that without using the scaling technique. The problem is the search of some suitable diagonal matrices D1 and D2 , such that the condition number of D1 AD2 is reduced as much as possible. Liu [18] has proposed a simple procedure to find D1 and D2 only through a few operations, which are derived explicitly. Previously, Liu [17] has used the concept of equilibrated matrix to choose the best source points for the method of fundamental solutions, and Liu and Atluri [19] have employed the concept of equilibrated matrix to find the best multiplescale of the Trefftz method used in the solution of inverse Cauchy problems for the Laplace equation, whose resulting linear system is less ill-conditioned. Liu [16] has developed a general purpose optimally scaled vector regularization method to treat ill-conditioned linear problems, according to the idea of equilibrated matrix. Before the application of scaling technique the norms of all columns in the coefficient matrix are in general not equal. If we demand that the norm of each column of the coefficient matrix of A is equal, the multiple-scale sij is determined by k=0 Do i = 1, m Do j = 1, i k =k+1 amax sij = ak  Enddo,

(19)

where ak denotes the kth column of A in Eq. (17), and amax is the maximum column norm of all columns of matrix A. Sometimes we can replace amax by a1  for saving time; hence, we have s11 = 1. Through this simple arrangement, in the new system Bc = b,

(20)

the n column norms of the new coefficient matrix B are equal. It can be seen that sij are fully determined by the collocation points. 6

Let k=0 Do i = 1, m Do j = 1, i k =k+1 Dk = sij Enddo,

(21)

and we can introduce a diagonal post-conditioning matrix by D := diag(D1 , . . . , Dn ),

(22)

such that the above equilibrated multiple-scale technique is equivalent to derive the new coefficient matrix B by B = AD, (23) where the dimension of B is still nc × n. 3. The 2D Stokes problem The governing law of the viscous incompressible flow is the Navier-Stokes equations. When the viscous force is very large in comparing with the inertial force, the equations can be reduced to the 2D Stokes equations: ∂u ∂v + = 0, (x, y) ∈ Ω, ∂x ∂y ∂p , (x, y) ∈ Ω, Δu = ∂x ∂p Δv = , (x, y) ∈ Ω, ∂y

(24) (25) (26)

where u(x, y) and v(x, y) are x-directional and y-directional velocity components, p(x, y) is the pressure and Δ is the Laplacian operator. Ω denotes the computational domain, whose enclosed boundary is denoted by Γ := ∂Ω. The above system of equations can be solved under suitable boundary conditions. While Eq. (24) is the continuity equation, Eqs. (25) and (26) are the x-directional and ydirectional momentum equations, respectively. By observing the system of above equations, the velocity components and pressure are coupled in these equations. So, it is not easy to analyze this system of equations directly. Since the unknowns in Eqs. (24)-(26) are velocities and pressure of the flow field, these three equations form the well-known primitive-variable formulation of the Stokes flow equations. From Eqs. (25) and (26) we can consider a new formulation: ∂ 2 p(x, y) ∂ 3u ∂ 3u = + , ∂x∂y ∂x2 ∂y ∂y 3 ∂ 2 p(x, y) ∂ 3v ∂ 3v = . + ∂y∂x ∂x3 ∂y 2 ∂x 7

(27) (28)

In view of pxy = pyx and Eq. (24), we come to a two-coupled equations for u and v: ∂ 3u ∂ 3v ∂ 3v ∂ 3u + = 0, (x, y) ∈ Ω, − − ∂x2 ∂y ∂y 3 ∂x3 ∂y 2 ∂x ∂u ∂v + = 0, (x, y) ∈ Ω. ∂x ∂y

(29) (30)

The above two equations under prescribed boundary conditions for u and v on Γ are our objective to be solved, which are third-first order system. On the other hand, by reducing Eq. (5) to the 2D Stokes flow and using Eq. (29) we can obtain another pair of third-order PDEs: ∂ 3u ∂ 3v ∂ 3u ∂ 3v − − + = 0, (x, y) ∈ Ω, ∂x2 ∂y ∂y 3 ∂x3 ∂y 2 ∂x ∂ 3u ∂ 3v ∂ 3v ∂u ∂v ∂ 3u + = 0, (x, y) ∈ Γ. + + + = 0, (x, y) ∈ Ω, 3 2 2 3 ∂x ∂x∂y ∂y∂x ∂y ∂x ∂y

(31) (32)

The above two equations by endowing with prescribed boundary conditions for u and v on Γ will be solved, which are third-third order system. In the multiple-scale Pascal polynomial expansion of u(x, y) and v(x, y) we can take u(x, y) =

i m  

c1ij s1ij xi−j y j−1 ,

(33)

c2ij s2ij xi−j y j−1 ,

(34)

i=1 j=1

v(x, y) =

i m   i=1 j=1

where the scales s1ij and s2ij are determined as that in Section 2 from the resulting linear system. The expansion coefficients c1ij and c2ij are to be determined. Inserting Eqs. (33) and (34) into Eqs. (29) and (30) we can obtain i m  

c1ij s1ij [(i − j)(i − j − 1)(j − 1)xi−j−2 y j−2 + (j − 1)(j − 2)(j − 3)xi−j y j−4 ]

(35)

i=1 j=1

−

i m  

c2ij s2ij [(i − j)(j − 1)(j − 2)xi−j−1 y j−3 + (i − j)(i − j − 1)(i − j − 2)xi−j−3 y j−1 ] = 0,

i=1 j=1 i m  

c1ij s1ij (i

− j)x

i−j−1 j−1

y

+

i=1 j=1

i m  

c2ij s2ij (j − 1)xi−j y j−2 = 0.

(36)

i=1 j=1

As that done in Eq. (14) we can let ⎡ ⎢ ⎢ ⎢ c1 = ⎢ ⎢ ⎣

c11 c12 c13 .. . c1n0

⎤

⎡

⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥=⎢ ⎥ ⎢ ⎦ ⎣

8

c111 c121 c122 .. . c1mm

⎤ ⎥ ⎥ ⎥ ⎥, ⎥ ⎦

(37)

⎡ ⎢ ⎢ ⎢ c2 = ⎢ ⎢ ⎣

c21 c22 c23 .. . c2n0

⎤

⎡

⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥=⎢ ⎥ ⎢ ⎦ ⎣

c211 c221 c222 .. . c2mm

⎤ ⎥ ⎥ ⎥ ⎥, ⎥ ⎦

(38)

where n0 = m(m+1)/2. When we collocate some points in the problem domain to satisfy the governing equations (35) and (36), and some points on the boundary to satisfy the boundary conditions, we can obtain a linear system with first setting s1ij = 1 and s2ij = 1: ⎤ ⎤ ⎡ ⎡ 0 G1 0  ⎢ 0 ⎥ ⎢ 0 G2 ⎥ c1 ⎥ ⎥ ⎢ ⎢ (39) ⎣ B1 0 ⎦ c2 = ⎣ b1 ⎦ , b2 0 B2 where G1 and G2 are linear operator due to the governing equations, and B1 and B2 are boundary linear operator due to the boundary conditions whose right-hand side data are represented by b1 and b2 . Then we apply the scaling technique as developed in Section 2 to find the scales s1ij and s2ij . When u and v are solved we can compute p by  i m   (i − j)(i − j − 1) i−j−2 j 2 2 i−j j−2 x + C, (40) cij sij y + (j − 1)x y p= j i=1 j=1 where C is a constant, which can be determined by a reference pressure at one point. The constant term does not appear in the field equations. Eq. (40) can be derived as follows. From Eqs. (25) and (33) it follows that px =

i m  

c1ij s1ij [(i − j)(i − j − 1)xi−j−2 y j−1 + (j − 1)(j − 2)xi−j y j−3 ],

(41)

i=1 j=1

where the subscript x in px denotes the partial differential with respect to x. Integrating it we have  i m   (j − 1)(j − 2) i−j+1 j−3 1 1 i−j−1 j−1 p= cij sij (i − j)x y + y + F (y), (42) x i−j+1 i=1 j=1 where F (y) is a function of y to be determined. On the other hand, from Eqs. (26) and (34) it follows that i m   c2ij s2ij [(i − j)(i − j − 1)xi−j−2 y j−1 + (j − 1)(j − 2)xi−j y j−3 ], (43) py = i=1 j=1

by equating which to the y-differential of Eq. (42) yields 

F (y) =

i m  

c2ij s2ij [(i − j)(i − j − 1)xi−j−2 y j−1 + (j − 1)(j − 2)xi−j y j−3 ]

i=1 j=1

−

i m  



c1ij s1ij

(i − j)(j − 1)x

i−j−1 j−2

y

i=1 j=1

9

(j − 1)(j − 2)(j − 3) i−j+1 j−4 + y x . (44) i−j+1

Integrating it with respect to y we can obtain (i − j)(i − j − 1) i−j−2 j i−j j−2 x y + (j − 1)x y F (y) = j i=1 j=1  i m   (j − 1)(j − 2) i−j+1 j−3 1 1 i−j−1 j−1 − + C. x cij sij (i − j)x y + y i−j+1 i=1 j=1 i m  



c2ij s2ij

(45)

Finally, inserting it into Eq. (42) we can derive Eq. (40). Alternatively, we can find p by solving the following equations when u and v are available: ∂p(x, y) ∂ 2 u(x, y) ∂ 2 u(x, y) = + , ∂x ∂x2 ∂y 2 ∂p(x, y) ∂ 2 v(x, y) ∂ 2 v(x, y) + . = ∂y ∂x2 ∂y 2 Upon taking p(x, y) =

m0  i 

c3ij s3ij xi−j y j−1 ,

(46) (47)

(48)

i=1 j=1

the scaling technique is used similarly to solve Eqs. (46) and (47) after collocating some points in Ω. Inserting the obtained expansion coefficients c3ij into the above equation we can find p. Here, we must emphasize that we do not need to specify the boundary condition of p on Γ. 4. Numerical examples of Stokes flows In order to demonstrate the performance of the multiple-scale Pascal polynomial expansion method (MSPEM), we solve some direct Stokes problems in this section and then the inverse Cauchy-Stokes problems in the next section. When Eqs. (29) and (30) are used we will call it the first method with an abbreviation MSPEM1 for solving the third-first order system, and when Eqs. (31) and (32) are used we have the second method with an abbreviation MSPEM2 for solving the third-third order system. On the other hand, for some cases we impose a disturbance on the specified boundary data for observing the sensitivity of the proposed new algorithms MSPEM1 and MSPEM2 to noise. 4.1. Example 1 We first consider the following exact solutions [8]: u(x, y) = 2xy +

y3 , 6

v(x, y) = x2 − y 2 + p(x, y) = xy,

x3 , 6 (49)

where the square computational domain is used and the velocity components, which are derived directly from the exact solutions, are specified along the boundary. 10

In this test, under the Dirichlet boundary conditions we solve this problem by using the MSPEM1: nc = 384 and n = 20 with m = 4. It means that the coefficient matrix in Eq. (39) has dimension 384 × 20. Under the convergence criterion ε = 10−10 , the CGM is convergence with 72 steps as shown in Fig. 1(a). The solutions u, v and p obtained by the MSPEM1 are very accurate with the maximum error being 8.64 × 10−12 , 2.40 × 10−12 and 2.94 × 10−11 , respectively. The logarithmic errors are defined by Log10 |u − uexact |, Log10 |v − vexact |, and Log10 |p − pexact |, which are plotted in Figs. 1(b)-1(d), respectively. In order to investigate the sensitivity of the MSPEM1 we impose a disturbance on the boundary conditions for u and v by adding a random noise with absolute value 0.01. Under the same parameter conditions and the same convergence criterion the MSPEM1 is convergence with 90 steps and with the maximum errors of the solutions u, v and p being 4.81 × 10−3 , 5.46 × 10−3 and 3.44 × 10−2 , respectively. The algorithm MSPEM1 is stable against large noise with the maximum errors being smaller than the magnitude of imposed noise. 4.2. Example 2 The second example is the well-known lid-driven cavity. The lid is moving with unit velocity in the x-direction and the boundary conditions in the other sides are assumed as non-slip [5, 10]: u(x, 0) = v(x, 0) = 0, u(1, y) = v(1, y) = 0 u(x, 1) = 1, v(x, 1) = 0, u(0, y) = v(0, y) = 0.

(50)

Under these boundary conditions we solve this problem by using the MSPEM1, where nc = 384 and n = 42 with m = 6. Under the convergence criterion ε = 10−8 , the CGM is convergence with 1512 steps for MSPEM1, and 1349 steps for MSPEM2. The MSPEM2 is slightly faster than MSPEM1. The solutions u and v obtained by the MSPEM1 and MSPEM2 are very accurate, which as usual are plotted in Fig. 2(a) of the profile of u with respect to y on the centerline x = 0.5, while that in Fig. 2(b) is the profile of v with respect to x on the centerline y = 0.5. The numerical results of u and v at the center point (x, y) = (0.5, 0.5) are (u, v) = (−0.17436, −1.4 × 10−10 ) for MSPEM1, and (u, v) = (−0.17438, −4.39 × 10−13 ) for MSPEM2. Due to the symmetry of the flow we expect that v = 0 at the center point. It can be seen that the MSPEM2 is better than the MSPEM1 from this viewpoint. Through the scaling technique the column norms of the coefficient matrix B are equal, and hence with a zero standard deviation of column norms. Furthermore, we can detect the standard deviation σ of the nc row norms of the coefficient matrix A for both MSPEM1 and MSPEM2, in order to compare the conditioned behavior for rows. Let ri , i = 1, . . . , nc denote the ith row of A, and then we have c 1  ri , nc i=1

 nc 1   (ri  − ξ)2 , σ= nc i=1

n

ξ=

11

(51)

where ξ is the mean value and σ is the standard deviation, which measures the variation of the row norms. For this problem the value of σ for MSPEM1 is 23.508, while that for MSPEM2 is 21.591. It demonstrates that the third-third order system MSPEM2 is slightly better conditioned than the third-first order system MSPEM1. 4.3. Example 3 This example is a recirculating flow in a 2D circular cavity [13, 23], whose radius is one. In the upper half of the boundary the velocities u = − sin θ and v = cos θ are prescribed, while in the lower half boundary u = v = 0 are imposed. Under these boundary conditions we solve this problem by using the MSPEM1 with nc = 744 and n = 30 with m = 10, and MSPEM2 with nc = 944 and n = 30 with m = 5. Under the convergence criterion ε = 10−10 , the CGM is convergence with 42 steps for MSPEM1, and 39 steps for MSPEM2. The solutions u and v obtained by the MSPEM1 and MSPEM2 are plotted in Fig. 3(a) of the profile of u with respect to y on the centerline x = 0, while that in Fig. 3(b) is the profile of v with respect to x on the centerline y = 0, which are compared with the analytic ones [13]. We can see that the MSPEM1 and MSPEM2 are accurate. For this problem the value of σ for MSPEM1 is 7.64, while that for MSPEM2 is 8.5. This case demonstrates that the third-third order system MSPEM2 and the third-first order system MSPEM1 almost have the same conditioned behavior. 4.4. Example 4 Then, we consider the following exact solutions [23]: u(x, y) = 2x3 + 6xy 2 + 2xy − 16x, v(x, y) = −6x2 y − 2y 3 − 3x2 − y 2 + 16y + 12, p(x, y) = 12(x2 − y 2 ) − 8y,

(52)

where the computational domain is shown in Fig. 4(a) with the contour given by ρ(θ) = exp(sin θ) sin2 (2θ) + exp(cos θ) cos2 (2θ).

(53)

Under the Dirichlet boundary conditions we solve this problem by using the MSPEM1: nc = 240 and n = 30 with m = 5, and m0 = 4. Under the convergence criterion ε = 10−10 , the CGM is convergence with 117 steps as shown in Fig. 4(b) by solid line. On the other hand the second method MSPEM2 with nc = 102 and n = 20 with m = 4 is convergence with 36 steps as shown in Fig. 4(b) by dashed line. The solutions u, v and p obtained by the MSPEM1 are very accurate with the maximum error being 8.56 × 10−13 , 7.64 × 10−13 and 9.98 × 10−12 , respectively. Furthermore, the MSPEM2 is more accurate with the maximum error being 2.06 × 10−13 , 2.90 × 10−13 and 9.01 × 10−13 , respectively. For this problem the value of σ for MSPEM1 is 17.969, while that for MSPEM2 is 4.491. It is interesting that the third-third order system as realized by the MSPEM2 yields row norms with similar magnitudes and hence better conditioned than the third-first order system as realized by the MSPEM1. In order to carry out the sensitivity analysis of the proposed algorithms we impose disturbance on the boundary conditions for u and v by adding random noise with absolute 12

value 0.01. Under the same parameter conditions and the same convergence criterion the MSPEM1 is convergence with 116 steps and with the maximum errors of the solutions u, v and p being 5.38 × 10−3 , 5.89 × 10−3 and 4.59 × 10−2 , respectively. Under the same parameter conditions and the same convergence criterion the MSPEM2 is still convergence with 36 steps and with the maximum errors of the solutions u, v and p being 9.66×10−3 , 9.32× 10−3 and 3.3×10−2 , respectively. It can be seen that the algorithms MSPEM1 and MSPEM2 are efficient and stable against large noise with the maximum errors being smaller than the magnitude of imposed noise. 5. Inverse Cauchy-Stokes problems In this section we use the MSPEM to solve two typical inverse Cauchy-Stokes problems. Consider a domain Ω, velocities ui are prescribed at Γu , and tractions ti are prescribed at Γt . If Γu and Γt satisfy Γu ∪ Γt = Γ := ∂Ω and Γu ∩ Γt = ø then a direct problem is to be solved. If both ui and ti are only specified at part of the boundary, then an inverse Cauchy-Stokes problem is to be solved [1]. In the below we use the MSPEM to solve a more ill-posed inverse Cauchy-Stokes problem without needing of overspecified traction boundary condition. We also impose a noise on the specified velocities boundary data by  s  uˆ(x, y) = u(x, y) 1 + Random × , (x, y) ∈ Γ1 , 100  s , (x, y) ∈ Γ1 , (54) vˆ(x, y) = v(x, y) 1 + Random × 100 where Random is the random numbers in [−1, 1], and s is the intensity of noise. Let Γ2 = Γ/Γ1 , and we attempt to recover the boundary data on Γ2 . 5.1. Example 5 First, we consider the following exact solutions [6]: u(x, y) = 2xy +

y3 , 6

x3 v(x, y) = x2 − y 2 + , 6 x p(x, y) = xy − , 24

(55)

where the computational domain is shown in Fig. 5(a). The contour is given by  ρ(θ) =

 12  2 cos(2θ) + 1.1 − sin (2θ) , 0 ≤ θ ≤ 2π,

(56)

and the discrete points are the collocation points. Let Γ1 be the upper half boundary of Γ and we specify u and v on Γ1 under a large noise s = 3. We attempt to recover the boundary data of u and v on Γ2 , which is the lower half boundary of Γ. We take nc = 592 and n = 35 with m = 4, and m0 = 5 for MSPEM1 and nc = 792 and n = 35 with m = 4, and m0 = 5 for MSPEM2. Under the convergence criterion ε = 10−3 , the CGM is convergence with 51 steps for MSPEM1 as shown in Fig. 5(b) by solid 13

line, and 59 steps for MSPEM2 as shown in Fig. 5(b) by dashed line. The solutions u and v on the lower half boundary obtained by the MSPEM1 and MSPEM2 are very accurate with the maximum errors being 4.04 × 10−2 and 8.56 × 10−3 for MSPEM1, and being 3.71 × 10−2 and 8.53 × 10−3 for MSPEM2. The comparisons are made in Figs. 5(c) and 5(d). 5.2. Example 6 Then we consider the exact solutions as that given in example 3, where the computational domain is doubly-connected as shown in Fig. 6(a). The outer boundary is a circle with radius 4, and the inner contour is given by Eq. (53), which is an amoeba-like irregular shape, and the discrete points are the collocation points as shown in Fig. 6(a). Let Γ1 be the outer boundary and we specify u and v on Γ1 under a large noise s = 5. We attempt to recover the boundary data of u and v on the inner boundary Γ2 . We take nc = 824 and n = 56 with m = 7 for MSPEM1, and nc = 884 and n = 56 with m = 7 for MSPEM2. Under the convergence criterion ε = 10−5 , the MSPEM1 is convergence with 165 steps as shown in Fig. 6(b) by solid line, and the MSPEM2 is convergence with 134 steps as shown in Fig. 6(b) by dashed line. The solutions u and v on the inner boundary obtained by the MSPEM1 and MSPEM2 are very accurate, of which the comparisons are made in Figs. 6(c) and 6(d). Under m0 = 3 we can recover the pressure p in the whole domain, of which the maximum error is about 1.9. By knowing the value of p is about in the order of 160, the accuracy with error 1.9 is quite good. It can be seen that the MSPEM1 and MSPEM2 are robust to recover all data although we do not give overspecified data. For this inverse problem the value of σ for MSPEM1 is 1664.85, while that for MSPEM2 is 1753.05. It demonstrates that the third-first order system MSPEM1 and the third-third order system MSPEM2 almost have the same conditioned behavior, such that they are performed similarly. 6. Conclusions In this paper we have introduced the multiple-scale in the Pascal polynomial triangle as a trial solution of two-dimensional Stokes equations. According to the concept of equilibrated matrix we can determine all introduced scales a priori when the collocation points are given. The new algorithms are called the multiple-scale Pascal polynomial expansion methods (MSPEMs). The third-third order system as presented by the MSPEM2 is found better conditioned than the third-first order system with MSPEM1. It is amazing that the accuracy of the MSPEMs for the direct Stokes problems is higher to the orders of 10−11 to 10−13 . Although under a random disturbance for detecting the sensitivity to noise, the proposed methods are quite stable against large noise. We also use the MSPEMs to recover the unknown pressure-driven Stokes flow by solving the primitive equation, which are ones of the inverse Cauchy-Stokes problems without specifying extra data. Through the tests of two typical inverse Cauchy-Stokes problems we found that the MSPEM is accurate and stable against large noise to 5%, which showed that the new numerical methods are applicable to both the direct and inverse Cauchy problems in Stokes flows. Acknowledgements:

14

The authors thank the constructive comments from anonymous referees, which improved the quality of the paper. Taiwan’s Ministry of Science and Technology project NSC-102-2221-E002-125-MY3, granted to CS, is highly appreciated. CS has been a Chair Professor of Hohai University since 2016.

References [1] R. Aboulaich, I. Ben Saad, M. Hassine, Recovering boundary data: the Cauchy Stokes system, Applied Mathematical Modelling 37 (2013) 1-12. [2] Aboulaich R, Ben Abda A, Kallel M. A control type method for solving the CauchyStokes problem. Applied Mathematical Modelling 37 (2013) 4295-4304. [3] A. Barrero-Gil, The method of fundamental solutions without fictitious boundary for solving Stokes problems, Computers & Fluids 62 (2012) 86-90. [4] G. Bastay, T. Johansson, V.A. Kozlov, D. Lesnic, An alternating method for the stationary Stokes system, ZAMM 86 (2006) 268-280. [5] S. Chantasiriwan, Performance of multiquadric collocation method in solving lid-driven cavity flow problem with low Reynolds number, Computer Modeling in Engineering and Sciences 15 (2006) 137-146. [6] C.W. Chen, D.L. Young, C.C. Tsai, K. Murugesan, The method of fundamental solutions for inverse 2D Stokes problems, Computational Mechanics 37 (2005) 2-14. [7] Y.W. Chen, W. Yeih, C.S. Liu, J.R. Chang, Numerical simulation of the twodimensional sloshing problem using a multi-scaling Trefftz method, Engineering Analysis with Boundary Elements 36 (2012) 9-29. [8] A.E. Curteanu, L. Elliott, D.B. Ingham, D. Lesnic, Laplacian decomposition and the boundary element method for solving Stokes problems, Engineering Analysis with Boundary Elements 31 (2007) 501-513. [9] C.M. Fan, D.L. Young, Analysis of the 2D Stokes flows by the non-singular boundary integral equation method, International Mathematical Journal 2 (2002) 1199-1215. [10] C.M. Fan, H.H. Li, CL Kuo, The modified collocation Trefftz method and Laplacian decomposition for solving two-dimensional Stokes problems, Journal of Marine Science and Technology 19 (2011) 522-530. [11] C.M. Fan, H.H. Li, Solving the inverse Stokes problems by the modified collocation Trefftz method and Laplacian decomposition, Applied Mathematics and Computation 219 (2013) 6520-6535. [12] C.M. Fan, H.H. Li, Solving inverse Stokes problems by modified collocation Trefftz method, Journal of Computational and Applied Mathematics 268 (2014) 68-81.

15

[13] T.Y. Hwu, D.L. Young, Y.Y. Chen, Chaotic advection for Stokes flows in circular cavity, Journal of Engineering Mechanics 123 (1997) 774-782. [14] C.L. Kuo, J.R. Chang, C.S. Liu, The modified polynomial expansion method for solving the inverse heat source problems, Numerical Heat Transfer B: Fundamentals 63 (2013) 357-370. [15] C.S. Liu, A highly accurate multi-scale full/half-order polynomial interpolation, Compututers Materials and Continua 25 (2011) 239-263. [16] C.S. Liu, Optimally scaled vector regularization method to solve ill-posed linear problems, Applied Mathematics and Computation 218 (2012) 10602-10616. [17] C.S. Liu, An equilibrated method of fundamental solutions to choose the best source points for the Laplace equation, Engineering Analysis with Boundary Elements 36 (2012) 1235-1245. [18] C.S. Liu, A two-side equilibration method to reduce the condition number of an ill-posed linear system, Computer Modeling in Engineering and Sciences 91 (2013) 17-42. [19] C.S. Liu, S.N. Atluri, Numerical solution of the Laplacian Cauchy problem by using a better postconditioning collocation Trefftz method, Engineering Analysis with Boundary Elements 37 (2013) 74-83. [20] C.S. Liu, S.N. Atluri, A highly accurate technique for interpolations using very highorder polynomials, and its applications to some ill-posed linear problems, Computer Modeling in Engineering and Sciences 43 (2009) 253-276. [21] C.S. Liu, W. Yeih, S.N. Atluri, On solving the ill-conditioned system Ax = b: generalpurpose conditioners obtained from the boundary-collocation solution of the Laplace equation, using Trefftz expansions with multiple length scales, Computer Modeling in Engineering and Sciences 44 (2009) 281-311. [22] C. Pozrikidis, Boundary Integral and Singularity Methods for Linearized Viscous Flow. Cambridge University Press, New York, 1992. [23] D.L. Young, S.C. Jane, C.Y. Lin, C.L. Chiu, K.C. Chen, Solutions of 2D and 3D Stokes laws using multiquadrics method, Engineering Analysis with Boundary Elements 28 (2004) 1233-1243. [24] D.L. Young, C.W. Chen, C.M. Fan, K. Murugesan, C.C. Tsai, The method of fundamental solutions for Stokes flow in a rectangular cavity with cylinders, European Journal of Mechanics B/Fluids 24 (2005) 703-716. [25] D.L. Young, S.J. Jane, C.M. Fan, K. Murugesan, C.C. Tsai, The method of fundamental solutions for 2D and 3D Stokes problems, Journal of Computational Physics 211 (2006) 1-8. [26] D.L. Young, C.L. Chiu, C.M. Fan, C.C. Tsai, Y.C. Lin, Method of fundamental solutions for multidimensional Stokes equations by the dual-potential formulation, European Journal of Mechanics B/Fluids 25 (2006) 877-893. 16

[27] A. Zeb, L. Elliott, D.B. Ingham, D. Lesnic, The boundary element method for the solution of Stokes equations in two-dimensional domains, Engineering Analysis with Boundary Elements 22 (1998) 317-326. [28] A. Zeb, L. Elliott, D.B. Ingham, D. Lesnic, Boundary element two-dimensional solution of an inverse Stokes problem, Engineering Analysis with Boundary Elements 24 (2000) 75-88. [29] A. Zeb, L. Elliott, D.B. Ingham, D. Lessnic, An inverse Stokes problem using interior pressure data, Engineering Analysis with Boundary Elements 26 (2002) 739-745.

17

Captions: Fig. 1.

Fig. 2.

Fig. 3.

Fig. 4.

For example 1 of the Stokes flow, showing (a) the convergence rate of MSPEM1, (b) the logarithmic error of u, (c) the logarithmic error of v, and (d) the logarithmic error of p. For example 2 of a lid-driven Stokes flow, showing (a) the velocity profile of u on the centerline x = 0.5, and (b) the velocity profile of v on the centerline y = 0.5 of a square cavity, computed by MSPEM1 and MSPEM2. For example 3 of circular cavity Stokes flow, showing (a) the velocity profile of u on the centerline x = 0, and (b) the velocity profile of v on the centerline y = 0, computed by MSPEM1 and MSPEM2. For example 4 of the Stokes flow, showing (a) the domain and collocation points, and (b) the convergence rates of MSPEM1 and MSPEM2.

Fig. 5.

For example 5 of an inverse Cauchy-Stokes problem, showing (a) the domain and collocation points, (b) the convergence rates of MSPEM1 and MSPEM2, (c) the errors of u, and (d) the errors of v on lower half boundary.

Fig. 6.

For example 6 of an inverse Cauchy-Stokes problem, showing (a) the domain and collocation points, (b) the convergence rates of MSPEM1 and MSPEM2, (c) the errors of u, and (d) the errors of v on inner boundary.

Residual

(a) 1E+3 1E+2 1E+1 1E+0 1E-1 1E-2 1E-3 1E-4 1E-5 1E-6 1E-7 1E-8 1E-9 1E-10 1E-11 0

8

16

24

32

40

48

Number of steps (b)

56

64

72

(c)

(d)

Fig. 1. For example 1 of the Stokes flow, showing (a) the convergence rate of MSPEM1, (b) the logarithmic error of u, (c) the logarithmic error of v, and (d) the logarithmic error of p.

1.0

(a)

0.8

y

0.6

0.4

MSPEM1 0.2

MSPEM2

0.0 -0.2

0.0

0.2

0.4

0.6

0.8

1.0

u 0.2

(b) MSPEM1

v

MSPEM2

0.0

-0.2 0.0

0.2

0.4

0.6

0.8

1.0

x

Fig. 2. For example 2 of a lid-driven Stokes flow, showing (a) the velocity profile of u on the centerline x=0.5, and (b) the velocity profile of v on the centerline y=0.5 of a square cavity, computed by MSPEM1 and MSPEM2.

1.0

(a)

0.8 0.6 0.4

y

0.2 0.0 -0.2

MSPEM2

-0.4

MSPEM1

-0.6

Analyt ic [13]

-0.8 -1.0 -1.0

-0.8

-0.6

-0.4

-0.2

0.0

0.2

0.4

u 0.6

(b)

0.4

v

0.2 0.0 -0.2

MSPEM2 MSPEM1

-0.4

Analyt ic [13] -0.6 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0

x

Fig. 3. For example 3 of a circular cavity Stokes flow, showing (a) the velocity profile of u on the centerline x=0, and (b) the velocity profile of v on the centerline y=0, computed by MSPEM1 and MSPEM2.

2.0

(a)

1.5 1.0

y

0.5 0.0 -0.5 -1.0 -1.5

Residual

-2

-1

0

2

3

x

(b)

1E+5 1E+4 1E+3 1E+2 1E+1 1E+0 1E-1 1E-2 1E-3 1E-4 1E-5 1E-6 1E-7 1E-8 1E-9 1E-10 1E-11

1

MSPEM1 MSPEM2

0

40

80

120

Number of steps Fig. 4. For example 4 of the Stokes flow, showing (a) the domain and collocation points, and (b) the convergence rate of MSPEM1 and MSPEM2.

1.0

(a)

y

0.5

0.0

-0.5

-1.0 -2

-1

0

1

2

x

1E+4

(b)

1E+3

Residual

1E+2 1E+1 1E+0 1E-1

MSPEM1

1E-2

MSPEM2

1E-3 1E-4 0

10

20

30

40

50

Number of steps

60

70

1.2

(c)

Exact MSPEM1 under s= 3

u(T)

MSPEM2 under s= 3

0.0

-1.2 3.0 3

3.7

4.4

5.1

5.8

6.5

5.8

6.5

T (d)

Exact MSPEM1 under s= 3

v(T)

2

MSPEM2 under s= 3

1

0

-1 3.0

3.7

4.4

5.1

T Fig. 5. For example 5 of an inverse Cauchy-Stokes problem, showing (a) the domain and collocation points, (b) the convergence rates of MSPEM1 and MSPEM2, (c) the errors of u, and (d) the errors of v on lower half boundary.

y

4

(a)

0

-4

Residual

-4

0

x

(b)

1E+9 1E+8 1E+7 1E+6 1E+5 1E+4 1E+3 1E+2 1E+1 1E+0 1E-1 1E-2 1E-3 1E-4 1E-5 1E-6

4

MSPEM1 MSPEM2

0

20

40

60

80

100 120 140 160 180

Number of st eps

20

(c)

Exact MSPEM1 under s= 5

u(T)

MSPEM2 under s= 5

0

-20 0.0 30

0.8

1.6

2.4

3.2

4.0

4.8

5.6

6.4

T (d)

Exact MSPEM1 under s= 5

20

v(T)

MSPEM2 under s= 5 10

0

-10

-20 0.0

0.8

1.6

2.4

3.2

4.0

4.8

5.6

6.4

T Fig. 6. For example 6 of an inverse Cauchy-Stokes problem, showing (a) the domain and collocation points, (b) the convergence rates of MSPEM1 and MSPEM2, (c) the errors of u, and (d) the errors of v on inner boundary.