Method of fundamental solutions for the problem of doubly-periodic potential flow

The method of fundamental solutions, or the charge simulation method, [7, 16] is a fast numerical solver for potential problems

$$\begin{cases}-\triangle u=0&\mbox{in}\ \ \mathscr{D}\\ u=f&\mbox{on}\ \ \partial\mathscr{D},\end{cases}$$

(1)

where $\mathscr{D}$ is a domain in the $n$-dimensional Euclidean space $\mathbb{R}^{n}$, and $f$ is a function given on $\partial\mathscr{D}$. In two dimensional problems $(n=2)$, equalizing the Euclid plane $\mathbb{R}^{2}$ with the complex plane $\mathbb{C}$, the method of fundamental solutions gives an approximate solution of (1) of the form ²²2The scheme shown here is the invariant scheme [17, 18] proposed by Murota.

$$u(z)\simeq u_{N}(z)=Q_{0}-\frac{1}{2\pi}\sum_{j=1}^{N}Q_{j}\log|z-\zeta_{j}|\quad(\>z=x+\mathrm{i}y\>),$$

(2)

where $Q_{0},Q_{1},\ldots,Q_{N}$ are unknown real coefficients such that

$$\sum_{j=1}^{N}Q_{j}=0,$$

(3)

and $\zeta_{1},\ldots,\zeta_{N}$ are points given in $\mathbb{C}\setminus\overline{\mathscr{D}}$. We call $Q_{j}$ the “charges” and $\zeta_{j}$ the “charge points”. We remark that the approximate solution $u_{N}(z)$ exactly satisfies the Laplace equation in $\mathscr{D}$. Regarding the boundary condition, we pose the collocation condition on $u_{N}(z)$, namely, we assume that $u_{N}(z)$ satisfies the equations

$$u_{N}(z_{i})=f(z_{i}),\quad i=1,\ldots,N$$

(4)

with $z_{1},\ldots,z_{N}$ given on the boundary $\partial\mathscr{D}$, which are called the “collocation points”. The equations (4) are rewritten as

$$Q_{0}-\frac{1}{2\pi}\sum_{j=1}^{N}Q_{j}\log|z_{i}-\zeta_{j}|=f(z_{i}),\quad i=1,\ldots,N,$$

(5)

which, together with (3), form a system of linear equations with respect to $Q_{0},Q_{1},\ldots,Q_{N}$. We determine the unknowns $Q_{j}$ by solving the linear system (3) and (5) and obtain the approximate solution $u_{N}(z)$. The method of fundamental solutions has the advantages that it is easy to program, its computational cost is low, and it shows fast convergence such as exponential convergence [13] under some condition. It was first used for studies of electric field problems [26, 27], and now it is widely used in science and engineering, for example, problem of scattering of earthquake waves [25].

The method of fundamental solutions is also used for the approximation of complex analytic functions. Let $f(z)$ be an analytic function in a domain $\mathscr{D}\subset\mathbb{C}$. The real part of $f(z)$, which is a harmonic function in $\mathscr{D}$, can be approximated using the form (2), and the imaginary part of $f(z)$, which is the conjugate harmonic function of $\operatorname{Re}f(z)$, is approximated using

$$-\frac{1}{2\pi}\sum_{j=1}^{N}Q_{j}\arg(z-\zeta_{j}).$$

Then, the analytic function $f(z)$ is approximated using a linear combination of the complex logarithmic functions

$$Q_{0}-\frac{1}{2\pi}\sum_{j=1}^{N}Q_{j}\log(z-\zeta_{j}).$$

(6)

From this point of view, Amano applied the method of fundamental solution to numerical conformal mappings [1, 2, 3].

In this paper, we examine the problem of two-dimensional potential flow past an infinite doubly-periodic array of obstacles as shown in Figure 1. A two-dimensional potential flow is characterized by a complex velocity potential $f(z)$ [15], which is an analytic function in the flow domain $\mathscr{D}$ such that it gives the velocity field $\boldsymbol{v}=(u,v)$ by $f^{\prime}(z)=u-\mathrm{i}v$, and its imaginary part satisfies the boundary condition

$$\operatorname{Im}f=\mbox{constant}\quad\mbox{on}\ \ \partial\mathscr{D}.$$

(7)

Physically, the boundary condition (7) means that the fluid flows along the boundary $\partial\mathscr{D}$ since the contour lines of $\operatorname{Im}f$ are the streamlines. Therefore, in order to obtain the potential flow of our problem, we have to find an analytic function $f(z)$ in $\mathscr{D}$ satisfying the boundary condition (7).

We have, however, one problem in applying the method of fundamental solutions to our problem. The velocity field $\boldsymbol{v}$ is obviously a doubly-periodic function due to the the double periodicity of the flow domain $\mathscr{D}$. Then, the complex velocity potential $f(z)$ involves a doubly-periodic function, and it is difficult to approximate it using the form (6) of the conventional method. To overcome this challenge, we propose to solve our problem using a doubly-periodic fundamental solution, that is, a logarithmic potential with a doubly-periodic array of sources. We construct a doubly-periodic fundamental solution using the theta functions and approximate the complex velocity potential by a linear combination of the doubly-periodic fundamental solutions.

The previous works related to this paper are as follows. As studies of problems of periodic flow, Zick and Homsy [28] proposed an integral equation method for three-dimensional Stokes flow problems with a three-dimensional periodic array of spheres, where the solution is given by an integral including the periodic fundamental solution. Greengard and Kropinski [8] proposed an integral equation method for two-dimensional Stokes flow problems in doubly-periodic domains, where an approximate solution is given as a complex variable formulation and the fast multipole method is used. Liron [14] studied Stokes flow due to infinite array of Stokeslets and applied it to the problems of fluid transport by cilia. As studies on methods of fundamental solutions applied to problems with periodicity, Ogata et al. proposed a method of numerical conformal mappings of complex domains with a single periodicity [24], where the mapping function, an analytic function involving a singly-periodic function, is approximated by the method of fundamental solution using the singly-periodic fundamental solutions. They proposed a method of fundamental solutions also to the problems of two-dimensional periodic Stokes flows [21, 20], where the solutions are approximated by the periodic fundamental solutions of the Stokes equation, that is, the Stokes flows induced by a periodic array of point forces. The author proposed a method of fundamental solutions for the two-dimensional elasticity problem with one-dimensional periodicity [19], where the solution is approximated using the periodic fundamental solutions of the elastostatic equation, that is, the displacements induced by concentrated forces in a periodic array.

The remainder of this paper is structured as follows. Section 2 proposes a method of fundamental solutions for our problem of potential flow with double periodicity. Section 3 shows some numerical examples which show the effectiveness of our method. In Section 4, we conclude this paper and present problems related to future studies.

We consider a potential flow past the doubly-periodic array of obstacle. The flow domain is mathematically given by

$$\mathscr{D}=\mathbb{C}\setminus\bigcup_{m,n\in\mathbb{Z}}\overline{D_{mn}},$$

where $D_{0}$ is one of the obstacles, a simply-connected domain in $\mathbb{C}$ and

$$D_{mn}=\left\{\>z+m\omega_{1}+n\omega_{2}\>\left|\>z\in D_{00}\>\right.\right\},\quad m,n\in\mathbb{Z}$$

with complex numbers $\omega_{1}$ and $\omega_{2}$ giving the periods of the obstacle array such that $\operatorname{Im}(\omega_{2}/\omega_{1})>0$ and $\overline{D_{mn}}\cap\overline{D_{m^{\prime}n^{\prime}}}=\emptyset$ if $(m,n)\neq(m^{\prime},n^{\prime})$. We assume that the spatial average of the velocity field $\boldsymbol{v}=(u,v)$ is a uniform flow, in whose direction the real axis is taken, that is,

$$\langle\boldsymbol{v}\rangle=\frac{1}{|\mathscr{D}_{0}|}\iint_{\mathscr{D}_{0}}\boldsymbol{v}\mathrm{d}x\mathrm{d}y=(U,0),$$

(8)

where $U$ is the magnitude of the velocity of the unit flow, $\mathscr{D}_{0}$ is the set defined by

	$\displaystyle\mathscr{D}_{0}=\>$	$\displaystyle\mbox{(a period parallelogram)}\setminus\bigcup_{m,n\in\mathbb{Z}}\overline{D_{mn}}$
	$\displaystyle=\>$	$\displaystyle\left\{\>z_{0}+a_{1}\omega_{1}+a_{2}\omega_{2}\>\left|\right.\>0\leq a_{1},a_{2}\leq 1\>\right\}\setminus\bigcup_{m,n\in\{0,1\}}\overline{D_{mn}}$		(9)

with a point $z_{0}\in D_{00},$ (see Figure 2), and $|\mathscr{D}_{0}|$ is the area of $\mathscr{D}_{0}$.

As mentioned in the previous section, a two-dimensional potential flow is characterized by a complex velocity potential, $f(z)$, which is an analytic function in the flow domain $\mathscr{D}$ such that it gives the velocity field $\boldsymbol{v}=(u,v)$ by

$$f^{\prime}(z)=u-\mathrm{i}v,$$

(10)

and its imaginary part is constant on the boundary of the domain $\mathscr{D}$. The condition (8) is rewritten as

$$\frac{1}{|\mathscr{D}_{0}|}\iint_{\mathscr{D}_{0}}f^{\prime}(z)\mathrm{d}x\mathrm{d}y=U.$$

(11)

Therefore, our potential flow problem is mathematically reduced to the problem to find an analytic function $f(z)$ in $\mathscr{D}$ satisfying

$$\operatorname{Im}f=\mbox{constant}\quad\mbox{on}\ \ \partial\mathscr{D}$$

(12)

and (11).

In solving our problem by the method of fundamental solution method, however, we have one problem. The complex velocity potential $f(z)$ involves a periodic function since the velocity field $\boldsymbol{v}=(u,v)$ given by (10) is obviously a periodic function, and it is impossible to approximate $f(z)$ by the conventional method using the form (2). To overcome this challenge, we propose to approximate the complex velocity potential $f$ by

$$f(z)\simeq f_{N}(z)=Uz-\frac{\mathrm{i}}{2\pi}\sum_{j=1}^{N}Q_{j}\left\{\log\vartheta_{1}\left(\left.\frac{z-\zeta_{j}}{\omega_{1}}\right|\tau\right)-u_{j}z\right\},$$

(13)

where $\vartheta_{1}(v|\tau)$ is the theta function [4] ³³3 In [4], the theta function $\vartheta_{1}(z|\tau)$ is defined by $\displaystyle\vartheta_{1}(z|\tau)=$ $\displaystyle 2\sum_{n=0}^{\infty}q^{(n+1/2)^{2}}\sin(2n+1)z$ $\displaystyle=$ $\displaystyle 2q^{1/4}\sin z\prod_{n=1}^{\infty}(1-q^{2n})(1-2q^{2n}\cos 2z+q^{4n}),\quad q=\mathrm{e}^{\mathrm{i}\pi\tau}.$

	$\displaystyle\vartheta_{1}(v|\tau)=$	$\displaystyle 2\sum_{n=0}^{\infty}q^{(n+1/2)^{2}}\sin(2n+1)\pi v$
	$\displaystyle=$	$\displaystyle 2q^{1/4}\sin\pi v\prod_{n=1}^{\infty}(1-q^{2n})(1-2q^{2n}\cos 2\pi v+q^{4n})$

with $\tau=\omega_{2}/\omega_{1}$ and $q=\mathrm{e}^{\mathrm{i}\pi\tau}$,

$$u_{j}=\frac{1}{|\mathscr{D}_{0}|}\iint_{\mathscr{D}_{0}}\frac{1}{\omega_{1}}\frac{\vartheta_{1}^{\prime}((z-\zeta_{j})/\omega_{1}|\tau)}{\vartheta_{1}((z-\zeta_{j})/\omega_{1}|\tau)}\mathrm{d}x\mathrm{d}y,\quad j=1,\ldots,N,$$

(14)

$\zeta_{1},\ldots,\zeta_{N}$ are points given in $D_{00}$ and $Q_{1},\ldots,Q_{N}$ are unknown real coefficients such that

$$\sum_{j=1}^{N}Q_{j}=0.$$

(15)

We call $Q_{1},\ldots,Q_{N}$ the “charge” and $\zeta_{1},\ldots,\zeta_{N}$ the “charge points”. The theta function $\vartheta_{1}(v|\tau)$ is an entire function which has simple zeros at $m+n\tau$, $m,n\in\mathbb{Z}$ and satisfies the pseudo-periodicity

$$\vartheta_{1}(v+1|\tau)=-\vartheta_{1}(v|\tau),\quad\vartheta_{1}(v+\tau|\tau)=-q^{-1}\mathrm{e}^{-2\pi\mathrm{i}v}\vartheta_{1}(v|\tau).$$

(16)

Therefore, the functions $-(1/(2\pi))\log\vartheta_{1}((z-\zeta_{j})/\omega_{1}|\tau)$ appearing on the rightmost side of (13) can be regarded as the complex logarithmic potential with sources at the points $m\omega_{1}+n\omega_{2}+\zeta_{j},\ m,n\in\mathbb{Z},$ and it is a periodic fundamental solution of the two-dimensional Poisson equation ⁴⁴4 Hasimoto [10] presented a periodic fundamental solution of the two-dimensional Poisson equation in terms of the Weierstrass elliptic functions . The term $-u_{j}z$ is added in each term of the rightmost side of (13) so that $f_{N}(z)$ satisfies the condition (11).

The approximate potential $f_{N}(z)$ satisfies the pseudo-periodicity

	$\displaystyle f_{N}(z+\omega_{1})-f_{N}(z)=\>$	$\displaystyle\omega_{1}\left(U+\frac{\mathrm{i}}{2\pi}\sum_{j=1}^{N}Q_{j}u_{j}\right),$		(17)
	$\displaystyle f_{N}(z+\omega_{2})-f_{N}(z)=\>$	$\displaystyle U\omega_{2}+\sum_{j=1}^{N}Q_{j}\left(\frac{\zeta_{j}}{\omega_{1}}+\frac{\mathrm{i}}{2\pi}\omega_{2}u_{j}\right).$		(18)

In fact,

		$\displaystyle f_{N}(z+\omega_{1})-f_{N}(z)$
	$\displaystyle=\>$	$\displaystyle U\omega_{1}-\frac{\mathrm{i}}{2\pi}\sum_{j=1}^{N}Q_{j}\left\{\log\vartheta_{1}\left(\left.\frac{z-\zeta_{j}}{\omega_{1}}+1\right|\tau\right)-\log\vartheta_{1}\left(\left.\frac{z-\zeta_{j}}{\omega_{1}}\right|\tau\right)-u_{j}\omega_{1}\right\}$
	$\displaystyle=\>$	$\displaystyle U\omega_{1}-\frac{\mathrm{i}}{2\pi}\sum_{j=1}^{N}Q_{j}\left\{\log(-1)-u_{j}\omega_{1}\right\}=\omega_{1}\left(U+\frac{\mathrm{i}}{2\pi}\sum_{j=1}^{N}Q_{j}u_{j}\right),$
where (16) is used on the second equality and (15) on the third equality, and
		$\displaystyle f_{N}(z+\omega_{2})-f_{N}(z)$
	$\displaystyle=\>$	$\displaystyle U\omega_{2}-\frac{\mathrm{i}}{2\pi}\sum_{j=1}^{N}Q_{j}\left\{\log\vartheta_{1}\left(\left.\frac{z-\zeta_{j}}{\omega_{1}}+\tau\right|\tau\right)-\log\vartheta_{1}\left(\left.\frac{z-\zeta_{j}}{\omega_{1}}\right|\tau\right)-u_{j}\omega_{2}\right\}$
	$\displaystyle=\>$	$\displaystyle U\omega_{2}-\frac{\mathrm{i}}{2\pi}\sum_{j=1}^{N}Q_{j}\left\{\log(-q^{-1})-2\pi\mathrm{i}\frac{z-\zeta_{j}}{\omega_{1}}-u_{j}\omega_{2}\right\}$
	$\displaystyle=\>$	$\displaystyle U\omega_{2}+\sum_{j=1}^{N}Q_{j}\left(\frac{\zeta_{j}}{\omega_{1}}+\frac{\mathrm{i}}{2\pi}\omega_{2}u_{j}\right),$

where (16) is used on the second equality and (15) on the third equality. Then, the complex velocity $f_{N}^{\prime}(z)=u_{N}-\mathrm{i}v_{N}$ satisfies the double periodicity

$$f_{N}^{\prime}(z+\omega_{1})=f_{N}^{\prime}(z),\quad f_{N}^{\prime}(z+\omega_{2})=f_{N}^{\prime}(z),$$

(19)

which means that $f_{N}^{\prime}(z)$ is an elliptic function of periods $\omega_{1}$ and $\omega_{2}$.

The approximate velocity potential $f_{N}(z)$ in (13) is an analytic function in $\mathscr{D}$ such that it satisfies the condition (11) and gives the complex velocity $f_{N}^{\prime}(z)$ which is an elliptic function of periods $\omega_{1}$ and $\omega_{2}$. Regarding the boundary condition (12), we pose a collocation condition on $f_{N}(z)$, namely, we assume the equation

$$\operatorname{Im}f_{N}(z_{i})=C,\quad i=1,\ldots,N,$$

(20)

where $z_{1},\ldots,z_{N}$ are given points on $\partial D_{00}$ called the “collocation points” and $C$ is an unknown real constant. The equations (20) are rewritten as

$$-\frac{1}{2\pi}\sum_{j=1}^{N}Q_{j}\left\{\log\left|\vartheta_{1}\left(\left.\frac{z_{i}-\zeta_{j}}{\omega_{1}}\right|\tau\right)\right|-\operatorname{Re}(u_{j}z_{i})\right\}-C=-U(\operatorname{Im}z_{i}),\\ i=1,\ldots,N.$$

(21)

The equations (15) and (21) form a system of linear equations with respect to the unknowns $Q_{1},\ldots,Q_{N}$ and $C$. We determine the unknown charges $Q_{j}$ by solving this linear system and obtain the approximate velocity potential $f_{N}(z)$. Owing to the pseudo-periodicity (17) and (18), the approximate potential $f_{N}(z)$ automatically satisfies the collocation condition $\operatorname{Im}f=\mbox{constant}$ on the boundary of other obstacle $D_{mn}$, $m,n\in\mathbb{Z}$, that is,

$$\operatorname{Im}f_{N}(z_{i}+m\omega_{1}+n\omega_{2})=\mbox{constant},\quad i=1,\ldots,N.$$

In this section, we show some numerical examples which show the effectiveness of our method. All the computations were performed using programs coded in C++ with double precision.

Figure 3 shows the examples of flows past a doubly-periodic array of cylinders, that is, flows in the domain

	$$\displaystyle\mathscr{D}=\mathbb{C}\setminus\bigcup_{m,n\in\mathbb{Z}}\overline{D_{mn}},$$
where
	$$\displaystyle D_{mn}=\left\{\>z\in\mathbb{C}\>\left|\>|z-m\omega_{1}-n\omega_{2}|<r\>\right.\right\},\quad m,n\in\mathbb{Z}$$

with a positive constant $r$ and periods $\omega_{1},\omega_{2}\in\mathbb{C}$ such that $\operatorname{Im}(\omega_{2}/\omega_{1})>0$ and $\overline{D_{mn}}\cup\overline{D_{m^{\prime}n^{\prime}}}=\emptyset$ if $(m,n)\neq(m^{\prime},n^{\prime})$. The charge points $\zeta_{j}$ and the collocation points $z_{j}$ are respectively taken as

$$\zeta_{j}=qr\exp\left(\mathrm{i}\frac{2\pi(j-1)}{N}\right),\quad z_{j}=r\exp\left(\mathrm{i}\frac{2\pi(j-1)}{N}\right),\quad j=1,\ldots,N,$$

(22)

where $q$ is a constant such that $0<q<1$ and was taken as $q=0.7$ in the examples, and $u_{j}$ in (14) are computed by the Monte Carlo method with a million points in $\mathscr{D}_{0}$.

To estimate the accuracy of our method, we evaluated the value

$$\epsilon_{N}=\frac{1}{Ur}\max_{z\in\partial D_{00}}|\operatorname{Im}f_{N}(z)-C|,$$

where $C$ is the constant obtained in solving the system of linear equations (15) and (20). The value $\epsilon_{N}$ shows how accurately the approximate potential $f_{N}(z)$ satisfies the boundary condition (7). Figure 4 shows $\epsilon_{N}$ evaluated for the above numerical example with $\omega_{1}=4r$ and $\omega_{2}=4r\mathrm{i}$ for $q=0.4,0.5,0.6,0.7$, where $q$ is the parameter appearing in (22). The figure shows that the approximation of our method converges exponentially as the number of unknowns $N$ increases. Table 1 shows the decay rates of the error estimates $\epsilon_{N}$, which are also shown in Figure 4 in broken lines. The table shows that the error estimate $\epsilon_{N}$ roughly obeys

$$\epsilon_{N}=\mathrm{O}(q^{N})$$

for $q=0.5,0.6,0.7$ in the examples of $(\omega_{1},\omega_{2})=(4r,4r\mathrm{i}),(4r,4r\mathrm{e}^{\mathrm{i}\pi/3}),(4r\mathrm{e}^{\mathrm{i}\pi/6},4r\mathrm{i})$ and for $q=0.6,0.7$ in the example of $(\omega_{1},\omega_{2})=(4r,2r(1+\mathrm{i}))$.

In the actual computations, the condition (8) that the average of the velocity is $\langle\boldsymbol{v}\rangle=(U,0)$ is not exactly satisfied because the integrals on $\mathscr{D}_{0}$ giving $u_{j}$ by (13) are approximately computed by the Monte-Carlo method. We computed the average velocity $\langle\boldsymbol{v}\rangle$ by evaluating the integral in (11) using the Monte-Carlo method with a million points. Table 2 shows the result. The table shows that the condition (8) is approximately satisfied with error $\sim 10^{-4}$ to $10^{-5}$.

In this paper, we proposed a method of fundamental solutions for the problems of two-dimensional potential flow past a doubly-periodic array of obstacles. In terms of mathematics, our problem is to find the complex velocity potential, an analytic function in the flow domain with double periodicity. The solution obviously involves a doubly-periodic function, and it is difficult to approximate it by the conventional method. Then, we proposed a new method of fundamental solution for this problem using the periodic fundamental solutions, which is the logarithmic potential with a doubly-periodic array of sources constructed using the theta functions. The proposed method inherits the advantages of the conventional method of fundamental solutions and approximates the solution of our problem with double periodicity well. The numerical examples showed the effectiveness of our method.

We have two issues regarding this paper for future study. The first problem is a theoretical error estimate of our method. Theoretical studies on the accuracy of the method of fundamental solution have been presented for special cases such as two-dimensional Laplace equation and Helmholtz equation in a disk [13, 5, 22] and two-dimensional Laplace equation in a domain with an analytic boundary [12, 23]. However, theoretical error estimate still remains unknown for many types of problems and methods including our method. The author believes it one of the most important works on the method of fundamental solutions to give a theoretical error estimate.

The second problem is to extend our method to other problems with periodicity than two-dimensional potential problems. The author has already given methods of fundamental solutions for Stokes flow problems [20, 21] and elasticity problems [19] with periodicity. In the work on periodic Stokes flow [21], the author et al. proposed to use the periodic fundamental solutions, which was presented by Hasimoto [9] and is given by a Fourier series. It is expected to construct a method of fundamental solutions using periodic fundamental solutions given by the theta functions or elliptic functions [11, 6] as in this paper.

The author thanks Professor Yusaku Yamamoto for his advice in numerical computations.