A FAST MULTIPOLE METHOD FOR HIGHER ORDER VORTEX
PANELS IN TWO-DIMENSIONS
PRABHU RAMACHANDRAN
†
, S. C. RAJAN
‡
, AND M. RAMAKRISHNA
‡
Abstract. Higher order panel methods are used to solve the Laplace equation in the presence of
complex geometries. These methods are useful when globally accurate velocity or potential fields are
desired as in the case of vortex based fluid flow solvers. This paper develops a fast multipole algorithm
to compute velocity fields due to higher order, two-dimensional vortex panels. The technique is
applied to panels having a cubic geometry and a linear distribution of vorticity. The results of the
present method are compared with other available techniques.
Key words. Potential Flow, Fast Multipole Method, Higher Order Panels, Panel Method.
AMS subject classifications. 31A99, 34B60, 35J05, 65E05, 65Y99, 76M15
1. Introduction. Panel metho ds provide a means to solve incompressible and
inviscid fluid flows in two and three dimensions. These methods reduce the dimen-
sionality of the problem by one and hence in two dimensions require a one dimensional
“grid” on the boundary. The flow of an inviscid and incompressible fluid in the pres-
ence of a solid body can be simulated by discretizing the b ody into elements called
panels, distributing some singularity having an unknown strength on them and solv-
ing for the unknowns based on the boundary condition. In two dimensions the panel
elements can be linear, parabolic or higher order. The singularity distributed can
be a source, doublet or vorticity distribution. It can be lumped at a point or dis-
tributed as a constant, linear or higher order function. Katz and Plotkin [5] provide
comprehensive details on panel methods in general. Hess and Smith [4] laid the foun-
dation for the source panel method. The idea of the vortex panel method is due to
Martensen [8] and is extended by Lewis [6]. Vortex panels unlike source panels can
be used to simulate a lifting body. It is also well known that a polynomial distribu-
tion of doublets of order q can be replaced by a vorticity distribution of order q − 1.
Therefore, the present work uses a vorticity distribution on the surface of the panels.
Only two-dimensional flows are considered. It is to be noted that the ideas developed
in this work can also be used for source and doublet distributions.
In panel methods, the b oundary conditions can be represented in terms of the ve-
locity resulting in the Neumann condition, or in terms of the potential or stream func-
tion on the boundary, called the Dirichlet condition. If N panels with one unknown
strength per panel are used to discretize the boundary, a system of N linear equations
is obtained based on the chosen boundary condition. The system of equations can
be written as Ax = b, where x represents the vector of N unknown strengths. The
boundary value problem is solved when x is obtained from these system of equations.
In panel methods, A is usually dense. There are a variety of standard methods for the
solution of dense matrices. These algorithms are usually computationally inefficient.
For example, an LU decomposition requires O(N
3
) work. Rokhlin [12] uses a Gener-
alized Conjugate Residual Algorithm (GCRA), ordinarily requiring O(N
2
) work, to
solve the equations iteratively. In this algorithm one computes increasingly accurate
approximations to x, ˜x, by computing the matrix product A˜x. For well conditioned
†
Graduate Student, Department of Aerospace Engineering, IIT-Madras, Chennai 600 036
(prabhu@aero.iitm.ernet.in).
‡
Department of Aerospace Engineering, IIT-Madras, Chennai 600 036.
1
2 RAMACHANDRAN, RAJAN AND RAMAKRISHNA
matrices this sequence of approximations converges rapidly to x. Normally, the matrix
product A˜x involves O(N
2
) work. Rokhlin reduces this to an O(N) computation by
developing a fast multipole algorithm. This considerably speeds up the solution of the
boundary value problem. Rokhlin developed the method for point sources/vortices
(single layer potential) and doublets (double layer potential). These correspond to
lumping the singularity distribution of the panel at a point. The resulting panels
generate a singular velocity/potential field in the vicinity of these lumped points.
To eliminate the singular velocity or potential field, higher order distributions
of singularity can be used. When linear panel elements are used with a vorticity
distribution of any polynomial order q ≥ 1, it can be seen [10, 11] that there is a
logarithmic singularity at the edges of the panels. This is called the “edge effect”. This
edge effect can be eliminated by using a higher order panel method where the geometry
of the panel is such that there is no discontinuity in the slope b etween two adjacent
panels. For example, this is possible with cubic panels. The authors [10] developed
a method that eliminates the edge effect by using cubic panel elements with a linear
distribution of vorticity. This method was found to be about four times slower than the
linear panel method when evaluating velo cities or potentials for a given distribution
of singularity. In [11] the authors use a fast multipole technique[3, 2] to accelerate
the computation to acceptable speeds. The method uses a linear representation of
panels for the far field fast multipole computations and uses the cubic method for the
direct computations. While this hybrid approach eliminates the edge effect, it is not
as accurate as the cubic panel method.
The present work details a procedure to perform a fast multipole summation
using panels of any order and thus provides an accurate and fast technique for two-
dimensional panel methods. The method allows one to compute the matrix product
A˜x in O(N) time as done by Rokhlin [12] and thereby allows one to solve the b oundary
value problem efficiently. The technique is also of use when the velocity or potential
field due to the panels having a known distribution of singularity is required at a very
large number of points. Such a requirement arises in the context of vortex methods
where the velocity field due to a collection of panels on a large number of vortex
blobs is to be computed. The number of vortex blobs usually exceeds the number
of panels by two orders of magnitude or more. In these cases it is possible to solve
the boundary value problem using an LU decomposition since the number of panels
involved is small. However, as mentioned in [11, p. 6], the computation of the velocity
field due to the panels on the particles needs to be performed efficiently since the
particle-particle interaction is computed using an O(N) fast multipole method. The
developed algorithm makes this possible.
The fast multipole method for higher order vortex panels developed in this work is
demonstrated using cubic panels. It is compared with the hybrid algorithm described
in [11] and also compared with Anderson’s technique [1] extended to the current
problem.
2. Mathematical details. For clarity, the velocity field due to a cubic panel
element is first derived. Subsequently, the details of the fast multipole method are
elaborated.
2.1. Cubic panels. Consider an arbitrarily oriented cubic panel as shown in
figure 2.1. The start and end points of the panel are denoted as z
1
and z
2
respectively.
The panel can be translated and rotated such that the start point is at the origin and
the chord is along the x-axis in the z
0
plane. If the complex velocity due to the panel
FMM FOR 2D HIGHER ORDER PANELS 3
θ
z
2
0
y
x
z
1
z=x+iy
z plane
cubic panel
λ
Fig. 2.1. Sketch of a single cubic panel having a chord length λ in the z plane.
0
γ
γ
1
2
η
ς
z’ plane
y’
x’
z’=x’ + iy’
ς=
λ
Fig. 2.2. Sketch of a single cubic panel having a chord length λ in the z
0
plane.
in the z
0
plane is V (z
0
), then it is easy to see that the velocity, V (z), in the z plane
can be obtained as, V (z) = V (z
0
)e
−iθ
and that z = z
1
+ z
0
e
iθ
.
The equation of a cubic panel starting at the origin, with chord oriented along
the x-axis is given by η = a
1
ζ + a
2
ζ
2
+ a
3
ζ
3
, where ζ and η respectively are the x
and y coordinates of the panel surface, as shown in figure 2.2. Note that 0 ≤ ζ ≤ λ,
where λ is the chord length of the panel. The vorticity is distributed on the surface
of the panel and is linear with respect to ζ. The equation for the velocity at a point
z
0
due to such a panel is given as
V (z
0
) =
k
2πa
3
Z
λ
0
(ζ + γ
1
/k)
³
ζ
3
+
a
2
a
3
ζ
2
+
(a
1
−i)
a
3
ζ +
iz
0
a
3
´
dζ,(2.1)
where k = (γ
2
− γ
1
)/λ. As done in [11], the cubic in the denominator is factored as
follows
ζ
3
+
a
2
a
3
ζ
2
+
(a
1
− i)
a
3
ζ +
iz
0
a
3
= (ζ − a)(ζ − b)(ζ −c),(2.2)
4 RAMACHANDRAN, RAJAN AND RAMAKRISHNA
where a, b and c are the complex cube roots of the cubic [9]. Given the roots, closed
form integrals are obtained for the panel velocity and potential. After integration and
simplification the velocity due to the cubic panel in the z plane is,
V (z) =
−γ
2
2πa
3
λ
"
a log
¡
a−λ
a
¢
(a − c)(a − b)
+
b log
¡
b−λ
b
¢
(b − c)(b − a)
+
c log
¡
c−λ
c
¢
(c − a)(c − b)
#
e
−iθ
(2.3)
−
γ
1
2πa
3
λ
"
(λ − a) log
¡
a−λ
a
¢
(a − c)(a − b)
+
(λ − b) log
¡
b−λ
b
¢
(b − c)(b − a)
+
(λ − c) log
¡
c−λ
c
¢
(c − a)(c − b)
#
e
−iθ
,
where θ is the angle between the chord of the cubic panel and the x-axis. Similarly,
the complex potential due to the panel can be obtained. The expressions for these are
not reproduced here. Using these derived expressions one can compute the velocity or
potential due to a cubic panel with a linear vorticity distribution. It it to be noted that
if the panel geometry is almost parabolic (a
3
→ 0) or linear (a
3
, a
2
→ 0), the above
expression for the velocity field would be inaccurate. In such cases, the panels should
be represented as parabolic or linear elements and the velocity integral specialized
appropriately. These are straightforward to derive. The specialization does not affect
the developed method in any way.
2.2. FMM for higher order panels. Consider a higher order panel with its
chord oriented along the x-axis as illustrated in figure 2.2. Given a vorticity distri-
bution γ(ζ) it can be easily seen that the complex velocity at a point z
0
due to the
panel is given as,
V (z
0
) = u
0
− iv
0
=
−i
2π
Z
λ
0
γ(ζ)
z
0
− (ζ + iη)
dζ.(2.4)
Substituting ξ = ζ + iη, and p erforming a binomial expansion results in,
u
0
− iv
0
=
−i
2π
Z
λ
0
γ(ζ)
z
0
− ξ
dζ
=
−i
2π
∞
X
j=1
Z
λ
0
γ(ζ)ξ
j−1
z
0j
dζ(2.5)
Without loss of generality, if z
1
of the panel is assumed to be at the origin, then
z
0
= ze
−iθ
and u − iv = e
−iθ
(u
0
− iv
0
), and the above equation reduces to,
u − iv =
−i
2π
∞
X
j=1
e
i(j−1)θ
z
j
Z
λ
0
γ(ζ)ξ(ζ)
j−1
dζ.(2.6)
This can be written as,
u − iv =
−i
2π
∞
X
j=1
A
j
z
j
,(2.7)
where,
A
j
= e
i(j−1)θ
Z
λ
0
γ(ζ)ξ(ζ)
j−1
dζ.(2.8)
FMM FOR 2D HIGHER ORDER PANELS 5
The complex potential of the higher order panel can also be obtained as follows,
Φ = φ + iψ =
−i
2π
Z
λ
0
γ ln(z
0
− ξ)dζ
=
−i
2π
Ã
ln(z
0
)
Z
λ
0
γdζ −
∞
X
k=1
1
z
0k
Z
λ
0
ξ
k
k
dζ
!
=
−i
2π
"
P
0
ln(z
0
) −
∞
X
k=1
P
k
z
0k
#
(2.9)
where,
P
0
=
Z
λ
0
γdζ ; P
k
=
1
k
Z
λ
0
ξ
k
γ dζ.(2.10)
In the present work we are interested in evaluating the velocity fields. Hence, the
multipole method is developed with that in mind. The analysis of the truncation
errors is performed only for the velocity field. The truncation error for the complex
potential can also b e easily computed in a similar fashion as done for the velocity
field.
Given ξ(ζ), γ(ζ) and θ, A
j
can be readily computed. The series (2.7) converges
if |ξ| < |z|. It is reasonable to assume that the panel is completely contained inside
a circle centered at the origin having radius λ, i.e. |ξ(ζ)| < λ. Given this, it is easy
to see from (2.8) that the error involved in truncating the series to a finite number of
terms p is,
¯
¯
¯
¯
¯
¯
V (z) +
i
2π
p
X
j=1
A
j
z
j
¯
¯
¯
¯
¯
¯
≤
Γ
2πλ(% − 1)
µ
1
%
¶
p
(2.11)
where,
Γ =
Z
λ
0
|γ(ζ)|dζ
and % = |z|/λ.
Hence, equations (2.7) and (2.8) can be used to obtain a fast multipole expansion
for higher order panels. The coefficients A
j
are to be computed by numerical integra-
tion. It is to be noted that equation (2.7) is a multipole expansion about the point z
1
(the first point) of the panel as shown in figure 2.1. For a cubic panel as used in the
earlier work[11] we have, ξ = ζ + i(a
1
ζ + a
2
ζ
2
+ a
3
ζ
3
). It is also to be noted that if
the panels do not deform or change orientation, the coefficients, A
j
, are constant and
need be computed only once. If the panels only rotate, the entire integral need not be
evaluated and the coefficients need to be multiplied by a different value of e
i(j−1)θ
. If
the panels deform, the coefficients must be recomputed. Given equations (2.7) and
(2.8), the various expressions for the fast multipole method can be derived as follows.
2.3. Multipole expansion for a collection of panels. Given n panels placed
at points z
k
inside a circle of radius R, the multipole expansion for the velocity field
of the panels is,
V (z) = u − iv =
−i
2π
n
X
k=1
∞
X
j=1
A
kj
(z − z
k
)
j
(2.12)
![Fig. 5.6. Time taken versus number of particles on a square region. 400 panels are used for the cylinder. β = 7.0. Note that the curves for the present method and the hybrid method of Ramachandran et al. [11] almost coincide.](/figures/fig-5-6-time-taken-versus-number-of-particles-on-a-square-28taq2jq.webp)
