Hi, i wrote my personal 3D panel code for external aerodynamics and i'm having some trouble.

Because this is my first code i made use of vortex rings instead of source-doublet panels (no local frame of ref., no coord. transform, etc.); this is based on the fact that (i hope this is not a mistake) vortex rings are equivalent to costant doublet panels. The boundary conditions are of neumann type so, for each panel, i can write:

sum ( aij * Gj, j=1..Npanels)=-Vinf * ni

where aij is the inf. coeff. of the j-th vortex ring on the i-th panel (the control point is in the middle of the panel, in example xc=(x1+x2+x3+x4)/4 yc=...), Gj is the circulation of the j-th vortex ring and ni is the outward normal of the i-th panel. Following Katz & Plotkin (and i think they're right), when the j-th v.ring is on the upper surface of a wing and is at the trailing edge i have to add to its inf.coeff. that of the adiacent wake panel, when the ring is at the t.e. and on the lower surface to subtract to its inf.coeff. that of the same adiacent wake panel.

The wake is build by time stepping but having the wing fixed so, each time step, i move the v. rings of the wake by the sum of all induced velocity + the free stream velocity Vinf, obviously preserving their intensity.

The geometry is build by surfaces so, for example, the sphere is only one surface, a complete wing is build by 2 tip surfaces and 2 main surfaces (is build by 2 semi-wing and so i have two wakes per complete wing).

Now that i explained something about my code (thanks for reading all this) i can talk about my trouble.

First i tryed to simulate the potential flow about a sphere of unitary radius (To do this i modified the last eqtn. of my system and i imposed sum(Gj, j=1..Npanels)=0). Integrating the pressure over the surface i take CL=0 and CD=0 but if i plot the surface distribution of V i have that the maximum value of it is about 0.7*Vinf instead of 1.5*Vinf, as should be. I used 1600 panels so...

The second i tryed to simulate the flow about a rectangular wing. The first trouble is that the matrix of the system is rank deficient so, if i have 300 panels and 1 wing the rank is 299. Moreover if i have m panels and n wings the rank of the matrix is m-n. Even if this was completely compromising my whole code i chosed to calculate the M.P. inverse of the matrix (it produces the least squares solution of the system) just to check if at least the wake shedding subroutine was well written. Now start the fun because this approach produce every time a solution such that sum(Gj,j=1...Npanels)=0 but, if the angle of attack of my wing is not 0 (i used symmetric profiles over all the wing) it produces the tip vortices and the Gj in the wake goes to zero approximately in the right way. Moreover i have that CL and CD are not 0 but none of theese values match those of a right (for sure) calculation made by the lifting line method.

I also tryed to enforce the Kelvin condition on the last row of the matrix but it was still singular.

Thanks to anyone read all this. Any kind of help is appreciated. P.S.(Sorry for my bad english).

Maybe this is not the right place to post about panel method, i know. But maybe someone of you could help me in find a similar forum where i can ask about my problem because i'm going crazy and i really don't know what to do. Thanks

The best thing to do, first, is to check whether you can get your potential flow over a sphere working. If this doesn't work what makes you think a much more complicated case will? Two questions come to my mind:

(1) are you evaluating the a_ij terms analytically or numerically? If analytically make sure you account for special cases correctly. If numerically, then you need quite high order accuracy specially for panels that are adjacent to each other. In particular, how are you accounting for self-interaction; i.e., a_jj? Using integration or assigning a fixed value for the "solid angle"? These minor issues makes a significant difference in results...

(2) The method you're using is fairly low-order (in accuracy) ... You should do a convergence study to examine the trend. If the trend is correct, then perhaps there's nothing wrong with your implementation and it's only a matter of accuracy. If not, then obviously you have an implementation bug.

Rank deficiency should have been taken care of by removing (an arbitrary) row and replacing it with conservation of circulation.

Adrin Gharakhani

Thanks for your time Adrin.

1) In determining the a_ij terms i followed Katz & Plotkin. More precisely i'm using vortex ring as panels and velocity - neumann formulation so i determine the velocity induced by a vortex ring summing up the 4 velocity coming from the four sides of it, each determined as the velocity induced by a straight vortex segment. Just to be sure i derived by myself the routine "vortxl" in K. & P. and it's right. About the self-interaction...in effect i didn't pay too much attention to it because, as stated in K. & P. (and it seems to be right), this formulation is valid everywhere.

So i think the approach can be considered analytical. Just one parameter i had to choose, it was the cut-off distance. I choosed it accounting for the smallest panel i have in the geometry so that it's self-induced velocity can be computed and the matrix preserve it's diagonally-dominant character.

2) I didn't do any convergence study. In effect i dont'know what is the approximation of such a discretization and how the error should goes down when i refine the grid over the body (maybe first order). Moreover, if i choose the sphere analysis to check the convergence of the method, what do you think should be the control parameter ? Do you think that the maximum velocity achieved over the surface of the sphere could be the right one ? (Considering that this is the only parameter on which i have difference from the exact solution, all the other being with good accuracy right).

About the Kelvin condition (cons. of circ.) i tryed the substitution in randomly choosed rows of the grid but whatever row i choosed the matrix was still singular. Because the singularity comes out only when i add the kutta condition, my worry is that there's something wrong with it. Do you know if there's something wrong in applying the same kutta condition that comes out in source-doublet panel methods in a vortex ring panel method?

Thanks again for your time and sorry for my poor english.

Paolo

> So i think the approach can be considered analytical. Just one parameter i had to choose, it was the cut-off distance. I choosed it accounting for the smallest panel i have in the geometry so that it's self-induced velocity can be computed and the matrix preserve it's diagonally-dominant character.

I assume by cutoff you mean distance beyond which you take the panel effect on a collocation point to be negligible. This is strictly speaking not true - the effect of each vortex loop persists quite a bit. Try the matrix setup without any cutoff and see what happens to your results. What do the sum of rows add up to?

As for convergence the rate of convergence is immaterial for you - what matters is whether your solution converges to a known/exact result as you increase the grid resolution. That's all. For potential flow over a sphere you have an analytical solution for the surface velocity, so you can compare your velocities at the collocation points to this exact velocity profile. Make sure you compare _only_ at the collocation point - as you go away from this point the velocity values will be quite inaccurate.

The question for me is whether you see a 50% error for the maximum velocity at the top and fairly good comparison elsewhere - or do you see large errors everywhere. Again, check against the velocity profile from upstream to downstream for all collocation points.

As for the rank deficiency. If you replace one equation with conservation of circulation your matrix cannot/must not be rank deficient. Even if you don't replace one row with the conservation equation your matrix should still "invert"; that is, provide the correct result for _singly connected_ geometries, such as the sphere. For multiply connected geometries, such as a donut, you absolutely must include conservation of circulation; otherwise, as you increase grid resolution you'll end up diverging from the correct solution.

Adrin

What kind of information do you need , exactly?

Thanks again for your time Adrin. I've been investigating about my code.

About the cut-off distance, i changed it between 10^-8 and 10^(-16) and nothing changed for my regular paneled cylinder. So i think it is ininfluent regarding the errors i'm experiencing.

About the error for the velocity around the cylinder and the sphere i discovered that just half panel width away from the body surface i reach approximately the right velocities. So, could this be due to the fact that a vortex ring badly (or not at all) influence the tangential velocity in its control point? However the velocity surface distribution achieved is such that it is the right one but scaled by the factor 0.5 everywhere on the surface. This should suggest that there's a wrong 0.5 factor somewhere in the code but the fact that all over there is the right velocity distribution makes this statement wrong (i think).

About the single-multiple connected regions...i'm having opposite results. In example, if i simulate the flow over a cave cylinder (that is, without tip surfaces) i don't have rank deficiency. But whatever kind of geometry that is simple connected (a sphere, a wing with tip surfaces, closed cylinder) give me a singular matrix with a rank deficiency that is the number of simply connected bodies in my geometry. When for the cylinder and the sphere i sobstitute the condition

sum(Gj,j=1..Npanels)=0

then i remove the rank deficiency and the matrix is not still singular. But this doesn't work if i use a Kelvin condition for a wing, so..do you know what should be a good statement for a kelvin condition over a wing on which i already have imposed the kutta condition?

Thanks

I don't know exactly what's wrong with my code. The questions are:

1) Is the vortex approach correct, even if i'm simulating non-lifting bodies?

2)If it is correct, is a neumann condition all over the body(ies) the right boundary condition for every kind of problem?

3) It is correct the kutta condition that i'm using, that is:

Gw=Gu-Gl

?

4) What should be a right kelvin condition (cons. of circ.) for a lifting 3-D body?

Reading some books (i.e. Katz & Plotkin) i've seen that the answer to the first question is "yes" but is no used all over the world (and i don't undersytand why); about the second question...i found that a lot of codes uses this condition but what i think is that it should be wrong, because fixing all over the boundary of the domain the neumann (d/dn) condition gives inf solutions, each differing from the others by a costant. about the third question, this kind of kutta condition is the most widely used, even if it is used with source-doublet codes, and i used it too.

The only thing i know is that i would not convert my whole code in a source-doublet one, so i'd like to know if theoretically i used the right approach and, if someone of you already had this kind of problem so you could help me.

Thanks

I write just to thank you again for your time Adrin. Finally i find what i was looking for on the net. If you are interested, you can check at:

http://www.aoe.vt.edu/~devenpor/aoe5104/

lessons 21 and 22.

This is a site from a Professor at Virginia Tech. Prof. Davenpor used the same method i'm using for my code, with the only difference in the way he handles the wake. There are also two different codes written in MATLAB, one for non-lifting bodies and the other for 4-digit NACA wings.

The velocity distribution on the surfaces of lifting and non lifting bodies, as he writes, must be computed adding a principal value term computed as -0.5*grad(G).

About the singularity, he solves the problem in two ways for lifting and non-lifting bodies. I'm tryng to understand what he does but i've substantially resolved all my problems so, thanks again.

Best regards. Paolo

Paolo,

The factor 0.5 (difference) that you see on the boundary suggests to me that you're not computing the self-interactions correctly. You're either missing a term in your integrals or your integrals don't solve the correct integrals accurately. Note that if you're right _on_ the surface you'll get a 0.5 factor (related to the solid angle), but if you're epsilon away that term should go to zero (any point in the fluid would have a zero degree solid angle - this is a crude explanation but it hopefully will make sense when you look at the equations). The best way to handle this is to put the factor 0.5 on the diagonal (bypassing the integrals for the self interaction terms - also you have to be careful with the choice of 4*Pi, etc.)

As for rank deficiency, I'd forgot that you're enforcing the Neuman BC. Yes, for the latter the matrix is rank deficient. If you were to enforce the tangential BC, the matrix would not be rank deficient!

Good luck

Adrin

Paolo,

I posted my previous response essentially the same time as you posted yours. The 0.5*grad(G) term is the self-interaction term I was talking about...

Good luck

Adrin

[pranavladkat]

Hi,

I'm also need to code the 3D panel method for my course project for a wing. Can you share the link to the notes that you provided link in one of your response? Also can you redirect me to the references which will be helpful for the understand and coding 3D panel method?

Thank,
Pranav

[sbaffini]

At the moment i don't have access to those files. In the meanwhile consider also the following resources:

http://www.cfd-online.com/Forums/blo...el-method.html

Note that the codes of Prof. Davenport are actually available, only the notes are not. A classical reference for panel method is the book of Katz and Plotkin, Low Speed Aerodynamics.

For the vortex ring method applied to closed lifting surfaces the main reference is from Ashook Srivastava:

http://nal-ir.nal.res.in/1536/
http://nal-ir.nal.res.in/2933/
http://nal-ir.nal.res.in/946/

I also have these files but not at the moment.