A general high order two-dimensional panel method

TL;DR

This paper introduces a high-order, quadrature-free two-dimensional panel method optimized for airfoil design, balancing accuracy and computational efficiency, and provides an open-source implementation.

Contribution

It presents a novel high-order panel method that is quadrature-free and integrates with a fast multipole method for improved efficiency.

Findings

Method achieves high accuracy with reduced computational cost

Trade-offs between approximation order and complexity are analyzed

Open-source implementation available for practical use

Abstract

We develop an efficient and high order panel method with applications in airfoil design. Through the use of analytic work and careful considerations near singularities our approach is quadrature-free. The resulting method is examined with respect to accuracy and efficiency and we discuss the different trade-offs in approximation order and computational complexity. A reference implementation within a package for a two-dimensional fast multipole method is distributed freely.

Figures

Equations

$$\vec{V}=\nabla \varphi +{\vec{V}}_{\omega },$$

$$\vec{V}\cdot \hat{n}=0.$$

$$\overline{V\left(z\right)}=\frac{1}{2\pi }\frac{Q+i\Gamma }{z-z_v},$$

$$z^{\prime }=z^{\prime }(z)=\left(z-z_1\right)e^{-i\theta },$$

$$\lambda =\text{mbox}Re\left\{\left(z_2-z_1\right)e^{-i\theta }\right\}.$$

$$z_p=z_p(\zeta )=\zeta +i\eta \left(\zeta \right)=\zeta +i\sum _{k=0}^M{a}_k{\zeta }^k,$$

$$\gamma =\gamma (\zeta )=\sum _{j=0}^N{b}_j{\zeta }^j.$$

$$\overline{V\left(z\right)}$$

$$=\frac{e^{-i\theta }}{2\pi }\frac{-i}{a_M}\underset{0}{\overset{\lambda }{\int }}\frac{{\sum }_{j=0}^N{b}_j{\zeta }^j}{{\sum }_{k=0}^M{c}_k{\zeta }^k}d\zeta ,$$

$$c_0$$

$$\frac{{\sum }_{j=0}^N{b}_j{\zeta }^j}{{\sum }_{k=0}^M{c}_k{\zeta }^k}=\sum _{m=0}^{N-M}{d}_m{\zeta }^m+\frac{{\sum }_{j=0}^{M-1}{h}_j{\zeta }^j}{{\sum }_{k=0}^M{c}_k{\zeta }^k}$$

$$\sum _{k=0}^M{c}_k{\zeta }^k=\prod _{k=1}^M\left(\zeta -x_k\right).$$

$$\frac{{\sum }_{j=0}^N{b}_j{\zeta }^j}{{\sum }_{k=0}^M{c}_k{\zeta }^k}$$

$$=\sum _{m=0}^{N-M}{d}_m{\zeta }^m+\sum _{k=1}^M\frac{A_k}{\left(\zeta -x_k\right)},$$

$$A_k=\frac{{\sum }_{m=0}^{M-1}{h}_m{x}_k^m}{{\prod }_{n\ne k}\left(x_k-x_n\right)}.$$

$$\overline{V\left(z\right)}=\frac{e^{-i\theta }}{2\pi }\frac{-i}{a_M}\left(\sum _{m=0}^{N-M}\frac{d_m}{m+1}{\lambda }^{m+1}+\sum _{k=1}^M{A}_k\log \left(\frac{x_k-\lambda }{x_k}\right)\right).$$

$$\overline{V\left(z\right)}=\frac{1}{2\pi }\sum _{n=1}^{\infty }\frac{f_n}{{\left(z-z_0\right)}^n},$$

$$f_n$$

$$\underset{Panel}{\int }\vec{V}\cdot \hat{n}dS=0.$$

$$\underset{S}{\oint }\vec{V}\cdot \hat{n}dS=\underset{A}{\iint }\nabla \cdot \vec{V}dA.$$

$$\underset{S}{\oint }\vec{V}\cdot \hat{n}dS=\underset{Panel}{\int }\vec{V}\cdot \hat{n}dS+\underset{Line}{\int }\vec{V}\cdot \hat{n}dS,$$

$$\underset{Panel}{\int }\vec{V}\cdot \hat{n}dS=-\underset{Line}{\int }\vec{V}\cdot \hat{n}dS+\sum _{q\in A}q,$$

$$\underset{Line}{\int }\vec{V}\cdot \hat{n}dS=\text{mbox}Re\left(\frac{q}{2\pi }\log \left(\frac{z^{\prime }-\lambda }{z^{\prime }}\right)\right)=\text{mbox}Re\left(\frac{q}{2\pi }\log \left(\frac{z-z_2}{z-z_1}\right)\right).$$

$$\overline{V\left(z\right)}=\sum _{n=0}^P{p}_n{\left(z-z_0\right)}^n,$$

$$\underset{Line}{\int }\vec{V}\cdot \hat{n}dS=\sum _{n=0}^P{p}_n\frac{{\left(z_2-z_0\right)}^{n+1}-{\left(z_1-z_0\right)}^{n+1}}{n+1}.$$

$$\displaystyle\int\limits_{Line}\vec{V}\cdot\hat{n}\,dS=\Biggl{[}f_{1}\log\left|z-z_{0}\right|-\sum_{n=2}^{P}\frac{f_{n}}{\left(n-1\right)\left(z-z_{0}\right)^{n-1}}\Biggr{]}_{z=z_{1}}^{z=z_{2}}.$$

$$\underset{0}{\overset{\lambda }{\int }}\frac{1}{2\pi }\gamma \log \left(\frac{z-z_2}{z-z_1}\right)$$

$$-\frac{1}{2\pi }$$

$$\frac{1}{2\pi }\underset{0}{\overset{\lambda }{\int }}\gamma \log \left(\frac{z-z_2}{z-z_1}\right)=$$

$$+\sum _{k=1}^M\left(\sum _{m=0}^N\left(s_m\left(x_{2,k}\right)-s_m\left(x_{1,k}\right)\right){\lambda }^m\right)\}$$

$$s_0\left(x_k\right)$$

$$s_m\left(x_k\right)$$

$$C=\pm i\underset{0}{\overset{\lambda }{\int }}q$$

Full text

Keywords: Boundary element method; Vortex method; Airfoil design; Potential flow; Fast multipole method.

AMS subject classification: Primary: 76M15, 76M23; Secondary: 65M38, 65M80.

1 Introduction

Most methods to numerically solve partial differential equations (PDEs) fall into one of two categories. The first is volume discretization methods, including, for example, finite element and finite volume methods. Here the resulting set of equations is large but sparse since the discretization nodes are connected only locally. If the PDE has a known fundamental solution, the full solution may instead be obtained in a discretization procedure involving only the boundary. This approach is commonly known as boundary element methods (BEMs) and generally involves fewer unknowns which, however, are connected globally.

In the current work, focus will be on fluid mechanical applications where the Laplace equation is used to calculate potential flow solutions. For many aerodynamics simulations, the target is a small object in the form of, e.g., an airfoil in a large domain which makes the BEM particularly attractive [1, 2, 3]. For time-dependent calculations, the method is commonly combined with the release and subsequent advection of vortices, effectively discretization points which approximates the flow. This is the method known as the vortex method [4]. Here the flow velocity $\vec{V}$ is calculated as a combination of a potential flow and of vortex contributions,

[TABLE]

where $\phi$ is the solution to a Laplace equation and where $\vec{V}_{\omega}$ is the contribution from the vortices in the flow [3, 5, 6]. Note that the contribution to the flow velocity from each vortex has to be calculated at each vortex position and at every time step. This is an $N$-body problem for which, as discussed below, fast algorithms should be employed.

One common method for solving the potential flow problem is to discretize the boundaries using panels. For fluid mechanical applications, these panels are constructed such that the flow satisfy the no-penetration boundary condition,

[TABLE]

The standard approach in 2D is to discretize the boundary using linear panels consisting of straight line segments, and to use a panel strength which is either constant or linearly varying along this line. These panels can be constructed from source and vortex sheets, but other possibilities also exist [7].

An obvious issue with panels that have a linear shape is that there will be sharp corners at the transitions between panels such that the solution to Laplace’s equation approaches infinity at the corner. For a pure potential flow solution, this is often not a problem, since the no penetration boundary condition is only satisfied at the centers of the panels. For vortex methods, however, this yields large numerical errors when evaluating the velocity in the vicinity of such transition points. A remedy is to use panels with higher order shapes to make the boundary smooth. For three-dimensional implementations, such panels have been discussed in [8, 9, 10]. Here, numerical integration is required to solve the flow contribution from the panel. This can be time consuming, especially for velocity evaluations close to the singularity of the fundamental solution, and it is therefore desirable to avoid numerical integration whenever possible. Indeed, for two-dimensional calculations, Ramachandran and co-authors have derived a solution involving panels with cubic shape and a linear distribution of the panel strength [11].

To improve the computational speed, vortex methods commonly rely on the fast multipole method (FMM) [12], and the same method can also be used to accelerate the solution of the dense BEM matrix with the influence coefficients of the panels [13, 11].

In the present paper we extend the work in [11] and design a general framework for two-dimensional panels with high order in both shape and strength and which does not require numerical integration. We develop the necessary analytic relations in §2, where we also discuss practical implementation issues allowing the method to be evaluated via the FMM. The performance of our framework is evaluated in §3 and conclusions are summarized in §4.

Our method has been implemented within the 2D FMM-software described in [14, 15], and is distributed as open source. See §4.1 for details.

2 Theory and implementation

We construct the panels from point sources/vortices in §2.1. The procedure for calculating the contribution from a panel is developed in §2.2, including both near and far-field evaluations. How the boundary conditions can be solved is described in §2.3, which includes how to integrate the contribution of a source point/panel over the panel surface. Finally, corrections to ensure continuity of the source strength are discussed in §2.4

2.1 Representation of panels

For a two dimensional flow using the complex number representation, the velocity $V$ at position $z$ from a source/vortex at position $z_{v}$ is given by

[TABLE]

where $\overline{V\left(z\right)}$ denotes the complex conjugate of the velocity, $Q$ is the source strength, and $\Gamma$ is the vortex strength.

Assume that a part of the boundary of an object extends between the points $z_{1}$ and $z_{2}$. To generate a high order panel that models this boundary, we need to select a baseline along which the panel is parameterized. Since only the direction of this reference line will be of importance, we let it start at $z_{1}$ with some angle $\theta$ pointing in the general direction of $z_{2}$, see Figure 2.1. The natural choice is to choose the baseline between $z_{1}$ and $z_{2}$, but to allow for the panel to be split (as relied upon in the FMM), it is necessary to allow that $z_{2}$ is not on the baseline.

Following the notation of [11], if we are interested in the flow velocity at position $z$, we apply the transformation

[TABLE]

which will make the reference line parallel to the real axis and the panel will extend between 0 and $\lambda$, where

[TABLE]

Using this reference system we can write a position $z_{p}$ on the panel as

[TABLE]

for a panel shape defined by a polynomial of degree $M$. To ensure that the shape of the panel remains numerically reasonable, it will be assumed that the coefficients $a_{k}$ are real.

The panel strength distribution is similarly defined by an $N$th degree polynomial,

[TABLE]

The coefficients $b_{j}$ are complex numbers, where the real part represents the source strength and the imaginary part the vortex strength according to (2.1).

2.2 Evaluation of velocities

The velocity contribution from the panel can be found by integration,

[TABLE]

where $z^{\prime}$ is defined in (2.2) and

[TABLE]

Hence the task can been reduced to an integration of the quotient of two complex polynomials. First we ensure that the order of the nominator is lower than the denominator (polynomial division),

[TABLE]

Now, since $c_{M}=1$, we can factorize the denominator,

[TABLE]

This allows a partial fraction decomposition, reducing the order of the denominator,

[TABLE]

where the constants $A_{k}$ are given by

[TABLE]

The final expression for the flow velocity from one panel thus becomes

[TABLE]

This is the general expression for evaluating the flow velocity from panels with arbitrary order polynomial shapes and strength distributions. So far, the only numerical step is to find the roots of the polynomial in the denominator, a well-studied problem for which efficient algorithms exist [16, 17].

For numerical issues to consider when using this expression, see A.

2.2.1 Outgoing expansion

Although (2.12) provides a solution for the velocity from an arbitrary polynomial panel, the expression is time consuming to calculate and can also give rise to cancellation errors when the distance $z^{\prime}$ in (2.2) is much larger than the panel length $\lambda$. For a distance sufficiently far away from the panel, a faster and more accurate way is to evaluate the velocity through a series expansion. Let $z_{0}$ be a suitably chosen expansion point for the panel. If we denote the radius of the smallest circle that encloses the panel by $R$, then if $\left|z-z_{0}\right|>R$, the evaluation can be carried out through the expansion

[TABLE]

Hence the coefficients $f_{n}$ can be evaluated through straightforward integration of a polynomial. The number of coefficients that has to be included in the series will depend on the distance $\left|z-z_{0}\right|/R$ and on the desired precision. Note that this expansion is on the same form as the expansions used in the fast multipole method. Hence, the panels can easily be included in the multipole calculations by shifting these expansions to the center of the FMM box containing the panel. Since the expansion is the same as used in the FMM, the same error estimates can be used to determine the required number of coefficients, see [12] for details. The current implementation uses these series expansions if $R<1.8\left|z-z_{0}\right|$, thus ensuring a relatively fast convergence of the series.

2.3 Solving the no-penetration boundary condition

The most common way of implementing the no-penetration boundary condition is to define one control point at the center of each panel and enforce that the normal flow velocity at this point is zero. For some applications it is desirable that the flow does not leak out of the domain, and it can be more suitable to satisfy the boundary condition that the net flow through the entire panel should be zero,

[TABLE]

This condition also makes the solution less sensitive to vortices located close to a control point.

2.3.1 Contribution from a point source

The contribution of one single vortex on a panel can be determined as follows. Start by applying Gauss’ law,

[TABLE]

As a closed contour $S$ is required, we close the contour by simply inserting a straight line between the end points. That is,

[TABLE]

which will give us the flow through the panel as

[TABLE]

where $q$ represents the strength of the vortex and where $A$ is the area enclosed by $S$. This shows that the integral only has to be calculated along the line, and then a correction is applied whenever the vortex is situated within $A$. In turn, the line integral is given by

[TABLE]

Note that the imaginary part of the expression is the flow parallel to the panel.

If so-called smoothing kernels of the vortices are used [4], and if the kernel overlaps with the panel, numerical integration may be the best option.

2.3.2 Evaluation from local- and far series expansions

In the fast multipole method, each multipole box will be associated with a local field expansion of the form

[TABLE]

where $z_{0}$ is the center of the expansion. This expansion contains the contribution from all vortices far away, and it is hence much faster to evaluate the net panel flow through this expansion rather than from (2.19). If we carry out the integral over the line, we find

[TABLE]

When evaluating the contribution from another panel far away, yet close enough to be in the near-field of the fast multipole method, it is preferable to carry out the evaluation from the far-field expansion around some point $z_{0}$. The general form of these expansions is given in (2.13), here truncated to $P$ coefficients. We find

[TABLE]

This expression can be used when (2.13) is valid over the line of integration.

2.3.3 Panel to panel interaction

To use the flow through the panel when generating the dense BEM interaction matrix, we need to integrate the flow from one panel over another panel. This expression can be obtained by integrating (2.19) over the source panel. By rewriting the panel the same way as in (2.6), we find

[TABLE]

The constant $N_{1}$ has been included to compensate for the possible changes in the branch of the logarithm. Equation (2.23) shows that the solution is reduced to the solution of two logarithmic potentials. A brief derivation of the solution to the logarithmic potential is included in B.

Using the expression for the logarithmic potential and by assuming that the panels do not intersect each other, the panel to panel interaction can be written as

[TABLE]

where

[TABLE]

Here, $x_{1,k}$ and $x_{2,k}$ are the roots to the polynomials defined by coefficients $c_{1,k}$ and $c_{2,k}$ respectively.

The above expression experiences a problem when the end point $z_{2}$ of the source panel coincide with any of the end points of the evaluation panel. If this is the case, the easiest solution is to generate a new panel which is parameterized in the opposite direction, hence moving the root from $\lambda$ to 0.

A correction has to be added whenever the panel end up within the area that is enclosed by the line between the evaluation panel end points and the evaluation panel itself. In this case, a constant with the magnitude

[TABLE]

should be added (see the last term in (2.23)), where the sign is chosen depending on the direction of the integral over the line and in relation to the direction in (2.18).

2.4 Continuity of the source strength

As stated in §1, to generate a smooth velocity field, the shape of the panels should have a continuous derivative and the strength should be continuous. In the current definition of the panels, however, the source strength is distributed along the reference line rather than along the tangential direction of the panel. It follows that the strength of the panel (in the tangential direction of the panel) will depend on the slope of the panel. To ensure continuity when determining the panel strength, it is convenient to define the panel by its strength (and derivatives of the strength for higher order continuity) in the end points. Hence, we want to give the strength in the end points along the tangential direction of the panel when defining the panel. One possible approach is to modify (2.5) to read

[TABLE]

Assuming that the slope of the panel is small, it is possible to expand (2.26) into a series, effectively a polynomial, where (2.12) applies anew. Any panels with large slope can be split into smaller panels until the series expansion is convergent. The disadvantages of this approach is that it generates additional panels and that the resulting polynomial in the nominator of (2.6) is of high degree. As an alternative approach, one can choose to only apply a correction factor to the given input values for the strength at the panel end points. Assume that the desired strength in the tangential direction of the panel is given by $\gamma_{input}$, then one can apply the correction at this point

[TABLE]

and then use (2.5) without any modification. This will give a different solution, than using (2.26), but continuity will still be fulfilled. The corrections for higher order continuity are obtained by differentiating (2.27) with respect to $\zeta$.

3 Numerical experiments

This section will illustrate the performance of the higher order panels for two different cases. In all cases, the panels are specified by giving the position, the derivative, the second order derivative and so on at each end point, where the number of derivatives are chosen depending on the order of the panel. We have tested panel shapes of all odd orders up to the 7th order. The panel strength is also chosen by giving the value, derivative, etc. at each end point, hence only odd orders are included here as well. This also means that panels of order 3 will have continuous first derivatives, panels of order 5 will additionally have continuous second derivatives etc. Consequently, 1st order panels imply 1 unknown per panel, 3rd order panels 2 unknowns, 5th order 3 unknowns and so on. All simulations will apply the correction to the panel strength at the end points according to (2.27) and its derivatives.

If more unknowns are added, there will also be a need for satisfying the zero flow condition on more panels. To generate more panels, the source panels are split into several evaluation panels while maintaining the same shape. Also, due to issues with BEM matrices that are close to singular, in the general case, it was concluded that better solutions are obtained if the source panels are split into one evaluation panel too much, and then the zero flow conditions is solved in a least square sense. All solutions in this section will apply this technique.

3.1 Flow around a deformed circle

For the first test case, a smooth shape in the form of a deformed circle is chosen. The deformation is accomplished by using the conformal mapping

[TABLE]

of a unit circle, and the resulting shape is illustrated in Figure 3.1. The deformation is added to the circle to avoid testing the panels on the symmetric circle, which artificially can give higher convergence than the general case.

The flow velocity when a conformal mapping is applied to a unit circle (without circulation) has the well-known analytic solution

[TABLE]

The measurement of error that will be used here is

[TABLE]

where $V$ is the velocity calculated from the panels, $S$ is a curve around the deformed circle at a fixed distance from the object and $L$ is the length of the curve $S$. This expression is evaluated through numerical integration with Gauss-Kronrod quadrature in Matlab [18].

For this case, panels with pure source strength are chosen (i.e. no vorticity) and the additional criterion that the total source strength of all panels should be zero is also enforced. All panels are distributed evenly around the unit circle in the $s$-plane.

As a first evaluation, the error is given in Figure 3.2 as a function of the distance from the deformed circle when using 50 panels. Some distance away from the circle (in the far-field of the panels), the error decreases with approximately the same factor for all panels, but the higher order panels give significantly better results for the same amount of panels. Very close to the circle (starting when the distance from the circle is about the same as the length of the panel), there is a significant change in the slope of the error curves and the error grows as the distance approaches the object. One general trend is that the order of the strength of the panels appears to be more important close to the surface of the object.

For linear strength ($N=1$), compared to the standard approach of using linear panels ($M=1$), there are significant improvements for increasing to cubic shape, but very modest improvements for the quintic shape. For cubic strength, a quintic shape is required to obtain significant improvements compared to the linear strength (while septic shape gives modest improvements) and a quintic shape shows some improvements for both the quintic and the septic shape. These trends are also visible in Figure 3.3, where the convergence in the far-field with respect to the number of panels is shown. Here, it can be seen that panels with linear shape have 2nd order convergence (regardless of the order of the strength), panels with cubic shape have 4th order convergence if strength is linear or higher, quintic panels with have 6th order convergence if strength is cubic or higher and septic panels initially have 8th order convergence if strength is quintic or higher. The trend change for the septic panels are likely related to numerical errors in the solution of the BEM equations. The maximum possible convergence rates in the far-field is according to the derivation in C equal to $\min(M+1,2N+2)$ under the condition that the optimal strength coefficients are obtained when solving the BEM equations. However, this assumes that the number of unknowns for the BEM equations is equal to $N+1$. Here, higher order continuity is enforced, which reduces the number of unknowns by a factor 2. Second, the BEM equations are applied in the near-field, and solves the no-penetration boundary condition, while the error is measured as the least square error in the far-field in Figure 3.3. The trend seen here is that the theoretical estimate applies to constant and linear strength, while when increasing the strength to cubic or higher, the reduction by a factor 2 applies (panels with constant strength are not shown in Figure 3.3 but have quadratic convergence with similar performance as $M=1$, $N=1$ for all tested shapes).

For the deformed circle case, the conclusion is that significant improvements can be gained from increasing the order of the panels. It is also clear that having a good approximation of the shape is more important in this case, than a good approximation of the strength.

3.2 Flow around an airfoil

In the second reference case, the flow around an airfoil is studied, which is a common application for both panel methods and vortex methods. Here, two aspects are important. The first is how well the flow close to the boundary is approximated. This is important for solutions of the boundary layer and therefore have many applications [6, 19]. The second aspect is how well the circulation around the airfoil is calculated, which will determine the calculated lift force. For this case, a Karman Trefftz airfoil is chosen, which has an analytic solution through the use of the conformal mapping

[TABLE]

which maps a cylinder in the complex $s$-plane with center in $\mu$ and radius $R_{ref}=\left|1-\mu\right|$ to an airfoil. The velocity for zero degrees pitch angle is obtained from

[TABLE]

where

[TABLE]

is the circulation. Here, the error in circulation is defined as

[TABLE]

and the error in flow velocity is calculated according to (3.3) by integrating along a line with with a fixed distance from the surface (an arc is added around the trailing edge to connect the upper and lower sides). The circulation around the airfoil is calculated through the Kutta condition [20] by setting the vortex strength to zero at the trailing edge [7].

Here, two different profiles (A and B) are tested, see Figure 3.4. The panels are generated by evenly separating the points in the $s$-plane with start and end panels on the trailing edge. This will generate smaller panels close to the trailing edge, which is the most difficult part to solve. The BEM matrix is then normalized to give all evaluation panels the same weight before the least square fit (which makes the flow close to the trailing edge more important). As the circulation is of interest, the airfoils are modeled with pure vortex panels.

For numerical reasons, defining higher order panels through their higher order derivatives does not work properly for quintic panels, as the second order derivative grows close to the trailing edge. Hence the order of shape for all points above 1.9 in the $x$-coordinate has to be reduced to cubic. Note that this only affects a very small part of the panels close to the trailing edge.

The calculated circulation for the two airfoils are shown in Figure 3.5. This shows that the convergence is similar for all panels, hence the high order convergence obtained with the circle is not present here. This can probably be related to the large gradients of the solution close to the trailing edge, which causes numerical difficulties. The oscillations seen for the linear panels are related to the modeling of the leading edge. The high accuracy values are obtained when the leading edge panel is perpendicular to the incoming flow, while the low accuracy values are obtained when the transition between two panels are at the leading edge, giving a sharp edge there. It should be noted that it is not the same panel orders that give the best results in both cases. Pure linear panels do for example give among the best results for airfoil A, but not for airfoil B. Apparently, determining the best panel order a priori is not so easy.

The convergence with respect to distance (Figure 3.6) shows similar trends as for the circle. Some distance away, most panel orders tend to have the same convergence. The exception is panels with linear shape and strength, which show significant improvements with respect to distance. This probably comes from the fact that the linear panels obtained a good value for the circulation, and having the correct total circulation becomes more important with distance. It is also clearly seen that close to the panel surface, the higher order panel shapes show significantly improved results. To further illustrate the origin of this, the error as a function of distance along the integration curve is shown in Figure 3.7 at a distance of 0.01 from the airfoil surface. Here, panels with linear shape have large errors over the entire airfoil surface, while the higher order panels show good values over most of the surface, but with large errors close to the trailing edge (and somewhat larger errors close to the leading edge).

Finally, the error in velocity at distance 0.01 as a function of the number of panels is illustrated in Figure 3.8. Similar to the error in circulation, all different panel orders tend to have relatively similar convergence behavior (compared to the large differences for the circle). Again, the quintic panels with quintic strength show the best results, but only quadratic convergence is obtained. Panels of septic order were tested as well, but without any improvements to the results.

3.2.1 Comparison of different panel implementations

One question is how big the difference is between solving for zero flow through the panels, or solving for zero flow velocity at control points in the center of the panels. To investigate this, a comparison between the two have been added in Figure 3.9 for airfoil A (error in flow velocity at distance 0.01). A small increase in performance is seen for quintic panels when using the flow through the panel while for cubic panels, no difference was seen. This shows that a control point approach is a decent approximation for this particular case. In the same figure, a comparison between using the correction to the panel strength at the end points according to (2.27), or by using (2.26) with series expansions of the correction to make the source distribution really follow the equation (red line). According to this case, the accuracy from using the series expansions is actually decreased, showing that corrections to the end points are a valid approximation, and is hence the recommended implementation as it is numerically more stable and faster, as the order of the strength polynomial is kept low and panels do not have to be split.

3.3 Computational complexity

We have seen that under certain conditions, the accuracy can be increased by using higher order panels. However, the evaluation speed also needs to be considered before choosing the order of the panel. Here we will assume that the panels is to be used in a vortex method, which means that it will be the panel-vortex interaction that will be the most time-consuming part, since this has to be performed every time-step, while the BEM matrices often only have to be calculated once at the start of the simulation.

3.3.1 Near field evaluation

In the first case, the evaluation within the near field of the panel will be considered. Here, 100 panels extending over the diagonal of the unit square are used, and 1 000 evaluation points are distributed randomly within the unit square. The time values measured are the wall clock times for all steps of the calculations; see Figure 3.10 for the results. Clearly, it is the order of the shape that dominates the time of the evaluation. For orders above 3, where numerical methods are required to calculate the roots in (2.9), there is a relatively rapid increase in the evaluation time. From a theoretical point of view, one can expect the growth both to have a linear component, see (2.12), and also a quadratic component from the root finding algorithm.

3.3.2 Fast multipole evaluation

In the second case, which is intended to illustrate applications in a vortex method, the unit circle from §3.1 will be used. Here, 100 000 vortices are randomly distributed using a homogeneous distribution in the polar plane between $0<\theta\leq 2\pi,\,1<r\leq 10$. The fast multipole method is configured to have about 35 particles in the lowest level boxes, which according to [15] is a suitable configuration for the current implementation for vortex dominated systems. The tolerance of the FMM is set to $10^{-9}$. In a real vortex method application, the velocity from the mutual interaction between the vortices also has to be calculated. This should be carried out in the same FMM evaluation as the panel interaction to reduce computational overhead. Hence, evaluation times with the panels included are compared with evaluation times with only the vortices present. This evaluation method was chosen to illustrate the added time of using the different panels and to exclude the time to build the FMM tree of the vortices (which otherwise will dominate over the vortex evaluation time). As Figure 3.10 showed that it is mainly the panel shape that determines the evaluation speed, the order of the panel strength is fixed to 3 in this evaluation. The results are presented in Figure 3.11 and shows that the time indeed increases with increasing panel order. However, by comparing with the particle to particle interaction time of 0.29s, the added time of the panels is still a small part of the total time.

Note that when comparing the evaluation times with Figure 3.10, 100 times more vortices are included here, implying about three orders of magnitude faster evaluation. Without the use of the FMM, the time would be about 0.81s (100 panels) for all panel orders, as the time for evaluating the far-field is independent of the panel order.

One trend seen for $M=6$ is that the evaluation time stops to decrease for low numbers of panels. This is a consequence of both that large panels give larger regions, where direct evaluation is necessary, and that large panels will cover several FMM boxes, and hence have to be split (this is likely the reason for the peak at about 50 panels).

3.4 Discussion

One of the main questions is which order of the panels should be used and as seen from the results this depends on the application. For the circle case, higher order panels generally give better results, but for the airfoil, the most suitable panels depend on what the code has to be used for. If only the circulation is of interest, it can be suitable to choose a low order panel method that is fast, but higher order panels have benefits for evaluations close to the airfoil surface. It is recommended to perform studies of the particular application before choosing panel method.

One other aspect of the choice between many panels of low order, or few of high order, is what happens for a real vortex application, where there are lots of vortices close to the boundary. Here, it can be beneficial to have many panels, as they more easily can adapt the strength to fulfill the no-penetration boundary condition. It is seen in Figure 3.11 that for high order panels, using too few panels will not decrease the evaluation time further. Hence, it may be a good idea not to use too high panel orders with very few panels, especially since the high order panels may experience difficulties when the solution has large derivatives, as is seen at the trailing edge of the airfoil case.

One additional question, not addressed in this work, is if one should use source or vortex panels, or a combination of both of them. From an evaluation point of view, as the whole algorithm is calculated with complex numbers, using both source and vortex strength should not give any noticeable increase in computational speed for the velocity evaluation, however the amount of unknowns will increase. This is an area where further studies can be performed for best performance.

4 Conclusions

We have presented a general high order panel method where, except for finding the roots of a polynomial, all parts have analytic solutions. The method hence allows for the evaluation of the panel contribution with high accuracy for all points in space. We have also presented how the flow through panels of high order can be calculated, both from point sources and from other high order panels. The necessary equations for implementing the panels in a fast multipole method has been presented as well. Together, these should provide the necessary basis for including higher order panels in vortex method applications.

We have demonstrated that for certain applications, a high order convergence can be obtained with higher order panels. However, in certain cases, especially close to points with large derivatives of the solution, numerical difficulties may reduce this high order convergence, showing that care has to be taken in these cases.

We have also demonstrated that the evaluation time mainly depend on the order of the shape of the panels, and high order panels are significantly more expensive to evaluate in the near field. However, the cost of evaluating the contribution from the higher order panels quickly decreases for evaluation further away when the fast multipole method can be applied. If the number of vortices is large, it is likely that the vortex to vortex interaction will dominate over the panel to vortex interaction in terms of computational cost.

4.1 Reproducibility

The full source code of the fast multipole high order panel solver is available as open source, and comes with a Matlab interface. All test codes used to generate the results in this work is also included, thus allowing for easy reproduction of the results. The code is available at www.stenglib.org.

Acknowledgment

This work was financially supported by the Swedish Energy Agency and was conducted within StandUp for Wind, a part of the StandUp for Energy strategic research framework (A. Goude), and by the Swedish Research Council within the UPMARC Linnaeus center of Excellence (S. Engblom).

Appendix A Numerical considerations

There are many numerical considerations that has to be taken into account for accurate use of (2.12). This section will briefly go through the most important numerical special cases that have to be handled.

A.1 Large roots

One case that can cause numerical issues is if one or several of the roots of the polynomial in (2.9) is large ($\left|x_{k}\right|\gg\lambda$) and this becomes particularly significant in the case when $N\geq M$. The problem can be illustrated by looking at the expression for polynomial division

[TABLE]

which has the recursive solution

[TABLE]

Hence the coefficients grow as $R^{j}$ which for large $R$ will cause large numerical errors (and increasing with the order of the strength $N$).

To reduce this problem one uses Taylor expansions for the large roots. Assume that $x_{i}$ is the large root and write

[TABLE]

This expression can be approached in the same way as (2.6), but now without the issue with the large root.

A.2 Double roots

Another case which may produce numerical errors is when two roots are close to each other, causing division by a small inaccurate number in (2.11). This can be handled by studying these two roots together. Assume that roots $j$ and $k$ are nearly the same. Then we have

[TABLE]

By expanding the logarithm of one root around the other we get

[TABLE]

Also, the sum can be expanded as

[TABLE]

This allows us to write the whole expansion as

[TABLE]

For convergence of $B_{j}$ in (A.7), a sufficient condition is that the distance between the roots is smaller than the distance for any of the roots to the panel end points.

A.3 Evaluation close to the panel end

For connected panels with continuous derivative in the shape and continuous strength, the velocity solution should be smooth on the panel surface when shifting between panels. However, (2.12) contains a weak logarithmic singularity at this point (when $x_{k}=0$ or $\lambda$). Although it is possible to derive similar techniques as for the double root in this case, as long as the overall method avoids evaluations in the nearest vicinity of the end points, this does not appear to be critical within the target tolerances considered here.

Appendix B Logarithmic potential

The logarithmic potential is generally given by

[TABLE]

Using the same procedure as for the potential in (2.1) we get

[TABLE]

In terms of

[TABLE]

where

[TABLE]

Combining (B.2) with (B.3) we get

[TABLE]

which is the final expression for the logarithmic potential. As in §A.1, in case of large roots, a series expansion might be employed for numerical reasons.

Appendix C Convergence with respect to number of panels

This section will cover the convergence estimate for the far-field region. The flow velocity of the panel is an expansion on the form

[TABLE]

In the far-field, where $z^{\prime}$ is large enough, one can readily expand the denominator in a geometric series containing negative powers of $z^{\prime}$. Hence by the truncation of polynomials in both the nominator and the denominator, one finds that the far-field convergence rate is $\min\left(M+1,N+1\right)$.

However, in terms of the strength, a higher convergence rate can be obtained than the suggested value $N+1$. In contrast to the shape, the source strength is not obtained from an explicit value. Instead, it is calculated implicitly by inverting the BEM matrix. For the far-field, we can look at this problem in terms of series expansions. Write (C.1) as an expansion around some point $z_{0}$,

[TABLE]

Similarly, write the expansion of the true source as

[TABLE]

It follows that the best possible coefficients, in the sense of the highest possible truncation order, are given by (cf. (2.14))

[TABLE]

Note that the constant terms $z_{0}$, $z_{1}$ and $a_{0}$ in (C.4) represent the distance between the starting point and the center of the panel, and will be of the order $\lambda/2$.

Combining (C.3) and (C.4), we can write the values for the coefficients as

[TABLE]

where $h_{n,k}(\lambda)=O(1)+O(\lambda)$. Note that for panels with linear shape, $h_{n,k}$ would be independent of $\lambda$. Now, assume that we want to approximate the source with a polynomial

[TABLE]

Instead of truncating $G$ in (C.3) to $N$ coefficients, we will instead choose the coefficients of $B$ to exactly match the $N+1$ lowest terms in the series expansion (C.5). This gives the system of equations

[TABLE]

which we can write in matrix notation as

[TABLE]

We can now split the vector $c$ into two parts

[TABLE]

where $m$ is the element index, starting from 0. Noting that $c_{1}=Ag$ we get

[TABLE]

where $b_{2}$ can be thought of as the error coefficients. As the lower order terms up to $\lambda^{N+1}$ belongs to $c_{1}$, while $b_{2}$ contains the remaining terms from $c_{2}$, the lowest order term in $b_{2}$ will be of order $b_{2,0}=O(\lambda^{N+1})$. Now, for any coefficient $f_{n}$, we have the approximation

[TABLE]

The error is just the difference between the approximation and the correct solution, thus,

[TABLE]

The lower order terms up to $n=N+1$ in (C.12) are 0 thanks to the matching in (C.7). Also, since the lowest order term in $b_{2}$ is $O(\lambda^{N+1})$, the lowest order term in (C.12) is $O(\lambda^{n+N+1})$. It follows that the lowest order term we can find in any of the coefficients is of order $\lambda^{2N+3}$. Finally, one notes that the number of panels increases with $1/\lambda$, and hence the total convergence rate in the far-field for the source strength is bounded by $\lambda^{2N+2}$.