Main Body

# VI. The Panel Method: An Introduction

The panel method is an analysis method that can be used to arrive at an approximate solution for the forces acting on an object in a flow. The method, as we present it here, is based on inviscid flow analysis, so it is limited to the resultant pressure forces over the surface. The panel method is basically a numerical approximation that relies on using discrete elements on the surface of an object and then prescribing a flow element (such as a vortex or doublet or source or sink) on each element that will satisfy certain boundary conditions (like no flow crosses the surface of the object). The interaction of the elements are accounted for and must also satisfy the condition that far from the object the flow should be equal to the free stream velocity approaching the object. There are a number of books and papers written that describe the method in very general terms and even the inclusion of viscous forces to some degree. But here we are just introducing the method to get a feel for its usefulness in external flows, so we will use a simply geometry with a simply distribution of flow elements. More complicated models exist but they all are based on the simplified form presented here.

We will assume that we have potential flow such that the governing equation for the flow field is the Laplace of the velocity potential, [latex]{\mathrm{\nabla }}^2\phi =0[/latex]. The boundary condition at an impermeable surface, where the velocity normal to the surface is zero, is [latex]\mathrm{\nabla }\phi \cdot \boldsymbol{n}=0[/latex]. Also, we can put our frame of reference on the object so fluid flow approaches the object. Keep in mind that since it is inviscid there may be a nonzero velocity component tangent to the surface. Also, the goal is to determine the velocity on the surface, and once this is found the Bernoulli equation can be used to find the local pressure distribution. The pressure can then be integrated over the surface to find the force by the fluid flow.

Without deriving this it can be shown that the following defines the velocity potential at any point P in the flow field (using Green’s Identify):

[latex]\phi \left (P \right)=\frac{1}{4\pi }\int{(\frac{\mathrm{\nabla }\phi }{r}-\phi \mathrm{\nabla }\frac{1}{r})\cdot ndS}[/latex]

(6.1)

where the integral is over the surface area of the flow (assuming two dimensional flow), S. This equation indicates that to solve for the velocity potential we must evaluate the integral on the flow boundaries (both the solid surface and infinitely far away).

This can be shown to be a result of the application of the divergence theorem that says:

[latex]\int{n_ju_jdS}=\int{\frac{\partial u_i}{\partial x_i}d\forall }[/latex]

(6.2)

where [latex]{u{}_{i}}[/latex] is the velocity vector. This states that the rate of expansion of a volume of fluid, given by the divergence of the velocity vector, within a given volume equals the flux through the volume boundaries. If the fluid is incompressible then the right side is zero and the net flux through the entire boundary is zero (this is the general argument for the conservation of mass of steady flow).

All this is nice and can be a powerful tool to find [latex]\phi[/latex] and therefore the velocity in the flow. But we really don’t need to find the potential of the entire flow, what our goal is, is to set up a flow that we know satisfies the velocity boundary condition at the surface of the object (with no velocity component across the surface) and far from the surface where the velocity is known to be the freestream velocity. Once this is established the force can be found. That is to say we only need to evaluate the surface velocity and then the pressure on the surface.

The general approach is to select a “grid” which is a series of “panels” that form the surface. Here we take the panels as straight flat surfaces arranged over the real surface. In the limit of very small panels the constructed surface will simulate the actual surface. On each panel we place a distribution of flow elements (like sources, sinks, vortices, etc) that when combined together will result in a flow field that will satisfy the surface boundary condition. There are lots of ways to identify which elements to use and how they may be distributed on the panels. Here we will use vortex elements, with one placed on each panel. The net flow is the result of superposition of the flow set up by each vortex on each panel element. So at each point in the field we add together the flow caused by all of the panel elements using the superposition rule. The panels that are far away from a given point will have less and less influence on the flow because the strength of the flow caused by a flow element decreases with distance from the element origin. For instance for a “source” the velocity decreases with increasing radial position because the flow is spreading out away from the source. But the influence never really goes to zero.

In placing a series of panels over the surface we first need to specify the size of each panel. We place a vortex of some strength some where on the panel (whose location is shown below) and we must identify points on the surface where we want to make sure that the velocity is zero across the surface. To be clear, individual points on each panel are used to evaluate the element (vortex) flow field –it needs its own origin, or coordinate system, to write an equation for the flow generated by this element. We also only pick a point on the panel to check to make sure that the net sum of contributions from all elements results in zero flow across the surface. The fact that we only satisfy the condition at one point on each panel will be satisfactory if the panels are made to be reasonably small. These points are called “collocation points” on each panel. In the end each panel will have coordinates that define its location on the surface, coordinates for the element location on the surface and coordinates for the collocation point.

The method of solution for the force on the object is the determination of the magnitudes or strengths of the elements on each panel. Once we have this distribution of strengths we can calculate the total lift on the surface that results from all of these elements. Recall that the lift experienced by a surface is really the component of the pressure at the surface integrated over the surface area in a direction normal to the approach velocity vector of the flow — it is not necessarily vertical, but normal to the freestream velocity.

For illustration we are going to use flat plate panels with vortex elements, one per panel. To get an idea of how to define the vortex origins on each panel we can examine flow over a single flat plate at some angle of attack to the freestream velocity vector, [latex]\alpha[/latex].

First we define the “center of pressure”, [latex]{x{}_{cp}}[/latex]. This is the location on the surface where the resultant lift force caused by the distributed load on the surface acts such that there is no net moment on the surface (see F.M. White, Chapt. 8, for some details on this). In general, if the moment on the surface about the leading edge is [latex]{Mo}[/latex], and [latex]{L}[/latex] is the lift force then using [latex]{x}[/latex] as the coordinate measured along the flat surface from the leading edge then:

[latex]{x{}_{cp\ }= -Mo/L}[/latex].

Consider now flow over the same flat surface, the lift we have seen is [latex]L= \rho {U} \Gamma[/latex] note that we are going to define [latex]\Gamma[/latex] as __positive clockwise__ — this is opposite to what we have done previously, but it helps get rid of some negative signs, this is not necessary, but convenient. To find out the value of [latex]{x{}_{cp}}[/latex] we use the theoretical evaluation of the lift force based on “thin airfoil theory”. The interested reader can refer to F.M. White Chapt. 8 for details on this. First, assuming a two dimensional flow, we assume that there is some continuous distribution of vorticity over the surface (in contrast to discrete vortex elements), such that there is a local circulation per unit length of surface, [latex]{\gamma (x)}[/latex]. The surface integral of this distribution results in the net total circulation: [latex]\mathrm{\Gamma }=\int^c_0{\gamma dx}[/latex] where “c” is the length of the flat surface, which we call the “chord”.

The Kutta condition is imposed where we want the flow to leave the surface flowing parallel with the back (trailing) edge. If this is a surface like an airfoil we want both the top and bottom flows to leave smoothly from the trailing edge. For this to happen then we don’t expect there to be a pressure difference between the top and bottom of the surface at the trailing edge, which would cause the flow to deflect up or down. If there is no velocity or pressure difference across the flow at the trailing edge the local circulation at the trailing edge, [latex]{(x=c)}[/latex], must be zero. So now we have a boundary condition for the integration of [latex]{\gamma (x)}[/latex] to get the total circulation [latex]\Gamma[/latex], that is [latex]{\gamma = 0}[/latex] at [latex]{x=c}[/latex]. The needed distribution has been figured out for a single flat surface at an angle of attack of [latex]\alpha[/latex], the angle between the oncoming flow and the surface. The function that works is: [latex]\gamma =2 U \sin \alpha \ {(\frac{c}{x}-1)}^{1/2}[/latex]. The integral of this expression for [latex]{\gamma(x)}[/latex] yields the value of [latex]\Gamma[/latex]. Once the circulation is known then the lift, [latex]{L,}[/latex] can be found. Using this lift, and also using the above definition of [latex]{x{}_{cp}}[/latex], it turns out that the “center of pressure occurs at [latex]{x{}_{cp} = c/4}[/latex]. So by assuming the lift force occurs at [latex]{x{}_{cp}}[/latex] then there is no moment caused by the distributed pressure force.

The lift force, [latex]{L}[/latex], as calculated as stated above, can be made nondimensional by divided by [latex]{1/2 \rho U{}^{2}b}[/latex], where b is the span of the surface (in and out of the page). The result is an expression for the lift coefficient:

[latex]C{}_{L} = 2 \pi \sin \alpha[/latex]

(6.3)

This result is the thin airfoil theory result for the lift on a surface at angle of attack, [latex]\alpha[/latex]. So at this point, we take this result for a flat surface and place vortex elements at the center of pressure for each flat element, which will be at a point ΒΌ of the distance along the panel from the beginning edge of the panel. Although this is not absolutely required for the panel method it is a convenient choice. As panel size is reduced smaller and smaller the impact of this choice becomes small.

The next step is to find the location where we want to impose the condition of zero velocity crossing the surface for each panel. For this we are using a coordinate system for the panel to be ([latex]{x,z}[/latex]) where [latex]{x}[/latex] is, again, along the flat surface from the panel leading edge and [latex]{z}[/latex] is normal (upward) from the surface. To find these collocation points, we note that the the velocity set up by a vortex element at location [latex]{r}[/latex] from the origin is [latex]v_{\theta }=-\mathrm{\Gamma}/2\pi r[/latex]. Writing this in Cartesian coordinates where at the panel surface [latex]v_{\theta }=\ w[/latex] (where [latex]{w}[/latex] is the [latex]{z}[/latex] component of velocity) is:

[latex]w=-\frac{\Gamma }{2\pi }\frac{(x-x_o)}{{(z-z_o)}^2+{(x-x_o)}^2}[/latex]

(6.4)

where [latex](x_o,z_o)[/latex] is taken as the center of pressure (where the vortex element is located). Also, since the collocation point is on the surface then both [latex]{z}[/latex] and [latex]z{}_{o}[/latex] are zero.

To get the total normal component velocity at the panel surface we add together [latex]{w}[/latex] being the velocity normal to the surface induced by the vortex, with the component of the freestream flow normal to the surface which is at angle [latex]{\alpha}[/latex] to the surface, [latex]U \sin \alpha[/latex], and set this to zero:

[latex]w+U\ sin\ \alpha =0[/latex]

(6.5)

We next insert Eqn. (6.4) for [latex]{w}[/latex], with [latex]{z{}_{o} = 0}[/latex] and [latex]{x{}_{0} = c/4}[/latex]. And we define the location at which we evaluate the zero normal velocity component as the collocation point as: [latex]{x= kc}[/latex]. Finally we solve for “[latex]{k}[/latex]” to be [latex]{k=3/4}[/latex]. This says that the [latex]{x}[/latex] location along any given panel where the boundary condition is to be satisfied is at [latex]{x=3 \ c/4}[/latex]. Recall that the vortex element is located at [latex]{x_{o} = c/4 }[/latex].

In summary, once the panels are set up on the surface, each has its own angle of attack, [latex]\alpha_{\rho}[/latex], depending on the local orientation of the surface. Each panel has a vortex element located at its own center of pressure and the normal velocity must be equal to zero at [latex]{x=3 \ c/4}[/latex] measured from the front of each panel. An equation can be written for the zero velocity condition at the collocation point of a panel by taking the sum of all of the contributions from all of the panels, each with a different vortex strength [latex]\Gamma_{\rho}[/latex], and setting this sum to zero. This will result in N equations (one for each panel) with N unknowns for [latex]\Gamma_{\rho}[/latex] for each panel. We can then solve for all [latex]\Gamma_{\rho}[/latex] using this set of equations.

Now we can set this up for a simple example to help pull all of this together, see Fig. 6.3. We use a thin flat surface representing a flat airfoil at an angle of attack of [latex]\alpha[/latex], with a uniform approaching flow of U. We use two panels for illustration. Note that the entire length of the two panels together is “c”, NOT the length of each panel. Also, [latex]{x}[/latex] and [latex]{z,}[/latex] used previously for a single panel, are along and normal to a line from the start to the end of the surface, respectively. This provides a coordinate system for the “system” of panels. For this simple case the angle of attack for both panels are the same, but we allow individual panel vortex strengths, [latex]{\Gamma}_{1}[/latex] and [latex]{\Gamma}_{2}[/latex].

In determining the vortex location for each panel, using the overall coordinate system we get (c/8,0) and (5c/8,0) for panel 1 and panel 2, respectively. The collocation points are (3c/8,0) and (7c/8,0), respectively.

Note that the each outward normal, [latex]n_{i}[/latex], points in the same direction for this case since both panels have the same orientation, but in general each panel could be at some angle [latex]\beta_{i}[/latex]. Fig. 6.4 shows the general definition of the panel angle to the chord, the latter being along the [latex]{x}[/latex] axis. The outward normal to the panel is given as [latex]n_{i}[/latex]. For the general case shown in Fig. 6.4 the outward normal is at an angle [latex]\beta_{i}[/latex] to the [latex]{z}[/latex] axis.

The general equation for the boundary condition is Eqn. (6.5), which we write slightly differently as:

[latex]v_i\cdot n_i=-U\cdot n_i[/latex]

(6.6)

Where [latex]v_i[/latex] is the induced velocity by the set of vorticies evaluated at each panel [latex]j[/latex] and [latex]n_i[/latex] is the panel outward normal. The left hand side is the dot product of the induced velocity with the outward normal. The right hand side is the component of the freestream velocity aligned with the outward normal. These two balance each other to yield zero velocity crossing the panel. The induced flow from all of the panel elements needs to be added together to get the value of [latex]v_{i}[/latex] in the above equation.

The general form for the velocity, which has components identified as [latex]{u,w}[/latex],at any point in the flow with a vortex element located at ([latex]x_o[/latex], [latex]z_o[/latex]) is obtained from above and written here as:

[latex]u=\frac{\mathrm{\Gamma }}{2\pi }\frac{(z-z_o)}{{(x-x_o)}^2+{(z-z_o)}^2}[/latex]

(6.7)

[latex]w=\frac{\mathrm{-}\mathrm{\Gamma }}{2\pi }\frac{(x-x_o)}{{(x-x_o)}^2+{(z-z_o)}^2}[/latex]

(6.8)

where we have transformed [latex]v_{\theta }[/latex] into the Cartesian components [latex]{u,w}[/latex].

Writing this as a matrix we obtain a general expression for the induced velocity anywhere in the flow [latex]\left ( x, z \right )[/latex]:

[latex]\binom{u}{w}= \frac{\Gamma}{2 \pi r^2} \Big ( \frac{01}{-10} \Big ) \binom{x-x_o}{z-z_o}[/latex]

(6.9)

Note here that [latex]x-x_o[/latex] is the distance along the x axis __from a vortex element__ to the point of interest where the velocity is being evaluated. Similarly for [latex]{{z-}z_o}[/latex]. Since the values for [latex]\Gamma[/latex] are not known, it is convenient to write this set of two equations for a “unit value” of [latex]\Gamma[/latex], that is [latex]\Gamma=1[/latex]. Keep in mind that the flow caused by each vortex when evaluated on a surface results only in a component normal to the surface. We also are going to use values of [latex]x_{o}[/latex] and [latex]{x}[/latex] based on the location of the vortices and collocation points, respectively.

Consequently, for our two panel example, where each panel is along the chord line, or the [latex]{x}[/latex] axis, we can write the velocity components given below, however here these velocities assume [latex]{\Gamma}=1[/latex] in the above equations. The notation used below is that the first subscript, [latex]{i}[/latex], represents the collocation point of interest at some panel, and the second subscript, [latex]{j}[/latex], identifies the vortex element at some panel that induces the velocity at the [latex]{i}[/latex] collocation point. That is to say these subscripts do NOT represent components of vectors in space in the conventional manner, and here [latex]u[/latex] and [latex]w[/latex] are the [latex]x[/latex] and [latex]z[/latex] component of the velocity at the surface.

[latex](u_{11}, w_{11}) = \left (0, -\frac{1}{\frac{2 \pi c}4} \right ) = (0, -\frac{2}{\pi c})[/latex]

(6.10)

[latex](u_{21}, w_{21}) = \left (0, -\frac{1}{\frac{2 \pi 3 c}4} \right ) = (0, -\frac{2}{3 \pi c})[/latex]

(6.11)

[latex](u_{12}, w_{12}) = \left ( 0, \frac{2}{\pi{c}} \right )[/latex]

(6.12)

[latex](u_{22}, w_{22}) = \left ( 0, - \frac{2}{\pi{c}} \right )[/latex]

(6.13)

For instance, [latex]u_{12}[/latex] is the [latex]{x }[/latex] component of velocity at the collocation point of panel 1 caused by the vortex element on panel 2. Also, [latex]w_{12}[/latex] is the z component of velocity at collacation point 1 caused by the vortex element on panel 2.

We can define, again for [latex]\Gamma_{\rho} = {1}[/latex], for each panel with vortex [latex]{j}[/latex], the following matrix [latex]a_{ij}[/latex]:

[latex]a_{ij}={\boldsymbol{v}}_{ij}\cdot \boldsymbol{n}[/latex]

(6.14)

This represents the component of the induced velocity in the outward normal direction for a given panel [latex]\Gamma=1[/latex]. We need to be careful here with this notation, [latex]{v}_{ij}[/latex] in the equation above is the velocity vector caused by the vorticies with unit value circulation with the subscripts: [latex]{i}[/latex] represents the collocation point and [latex]{j}[/latex] the vortex location inducing that flow. Also, the dot product with [latex]{n}[/latex] (the outward normal for each panel), gives the projection of [latex]{v}[/latex] normal to the surface. For example for the above example [latex]n[/latex] is only in the [latex]z[/latex] direction.

[latex]a_{11} =\left ( 0, - \frac{2}{\pi c} \right ) \cdot (0.1)=-\frac{2}{\pi c}[/latex]

(6.15)

so the elements [latex]a _{ij}[/latex] are quantities representing a measure of the velocity contribution normal to each panel from each vortex [latex]{j}[/latex] at each collocation point [latex]{i}[/latex] for a [latex]\Gamma_{j}=1[/latex] at each panel.

Once we have identified the components [latex]a_{ij}[/latex], which, as shown above for the two panel example, are determined only by geometry, it is possible to use the impervious boundary condition to find each of the values for [latex]\Gamma _{j}[/latex]. We have a general equation for our two element surface for panels [latex]{i}=1,2[/latex] that says the sum of the outward normal velocity at each panel with the freestream component normal to that panel must be zero. The outward normal velocity due to summation of the induced velocity over both panels is for some unknown distribution of circulations, [latex]\Gamma_{j}[/latex]:

[latex]\sum^2_{j=1}{a_{ij\ }{\mathrm{\Gamma }}_j}[/latex]

(6.16)

So this expression is balanced by the contribution from the free stream velocity: [latex]-U\cdot n[/latex] (this is the component of U normal to the surface).

The final system of equations become:

[latex]\sum^2_{j=1}{a_{ij}{\mathrm{\Gamma }}_j=-(U_x,U_z)\cdot (n_x,\ n_z)}[/latex]

(6.17)

Here the right hand side value of [latex]{U}[/latex] is represented as the vector ([latex]{U}_{x},{U}_{z}[/latex]) indicating that [latex]{U}[/latex] is a vector with [latex]{x}[/latex] and [latex]{z}[/latex] components that depends on the angle of attack given by the chord line. The [latex]{x,z}[/latex] components of the outward normal will depend on the angle [latex]{\beta}[/latex] for each panel. The geometry of all of the panels fully determines [latex]a_{ij}[/latex], as shown above in the simple two panel example. Consequently, once the geometry is known the values of [latex]a_{ij}[/latex] and [latex](n_x,\ n_z)[/latex] are all determined. Then knowing [latex]{U}[/latex] the set of Eqns. (6.17) provides the means to find all of the values of [latex]{\Gamma_j}[/latex], where the index represents each of the panel vortices. In the above two panel example there are two values of [latex]{\Gamma}[/latex].

From the solution of [latex]\Gamma_j[/latex] we can find the lift on each panel, [latex]\Delta L_j=\rho U \Gamma_j[/latex] and then the total lift as the sum over all panels:

[latex]\sum^2_{j=1}{\Delta L_j}=C_L\left(1/2\rho U^2c\right)[/latex]

(6.18)

Recall that the lift is normal to the direction of the freestream velocity, so the lift of each panel, even though the panels may all be at different angles of attack, are all in the same direction.

It is also possible at this point to calculate the velocity tangent to the surface at each panel. That is to say, instead of finding the component normal to the surface to satisfy the impermeable boundary condition find the component tangent to the panel, for each panel and combine this with the component of U tangent to the panel. The induced velocity tangent to the panel can be found using the dot product of the vortex induced velocity with the unit vector tangent to the panel. The tangent unit vector is given by its [latex]x[/latex] and [latex]z[/latex] components as [latex]cos\beta_i[/latex], [latex]sin\beta_i[/latex], respectively with [latex]\beta_i[/latex] being the angle of the panel relative to the chord line which is the [latex]x[/latex] axis, show in Fig 6.4. The freestream contribution is found from the value of U. This combined velocity yields the velocity vector of the flow at each panel and can be used in the Bernoulli equation to find the pressure at each panel. The Bernoulli equation is applied from a known upstream condition (known pressure and freestream velocity) to each panel in order to calculate the panel pressure.

Fluid Mechanics by F.M. White, Chapter 8 (introduction to thin airfoil theory)

Low Speed Aerodynamics by Katz \& Plotkin (details on the panel method)

Fundamentals of Aerodynamics by J.D. Anderson (more advanced airfoil theory)