Linear-Strength Vortex Unsteady Panel Method (2D)
 

Hi everyone,

I am working on a project for the development of a 2D Unsteady Panel Method for Airfoils. As often already suggested in this forum I have been using the book "Low Speed Aerodynamics" which has helped me a lot to produce the Steady State solver (I am writing the code in Matlab at the moment), which works very well.

However I am having great problems in validating my unsteady code.
I have followed the book, mainly the steps are:

1) compute the steady solution at first time step
2) move the airfoil via a prescribed motion path (in my case pitching A*sin(wt))
3) position a concentrated vortex at a certain distance from the T.E., which has unknown intensity GAMMA_W and which is therefore part of the solution
4) use the Kelvin theorem for completing the coefficients matrix A which will be N+2 times N+2 (N+1 gammas on the nodes of the panels) + 1 (latter) unknown shed concentrated vortex
5) compute the RHS (right-hand-side) term considering free stream velocity and wake induced velocity
6) solve the system

The problem seems that the Cp variations are not big enough. I am using a NACA report on wind tunnel experimental data (pitching motion 5 + 5sin(Wt) motion) to check the code. But also the comparison with results contained in the Katz and Plotkin (plot of CL-ALPHA of a 0012) isn't satisfactory.

I compute the Cp at each collocation point via the formula:

Cp_i = 1 - Vel_i^2

where Vel_i is the velocity at the i-th collocation point obtained from the different contributes: free stream (similar to steady state) + kinematic movement (sinusoidal motion) + induced velocities (by panels and by wake).
I consider the free stream V to be equal to 1.Honestly I am not quite sure that all of this is right (talking about Cp calculation).

Moreover I have great doubt about the Kutta condition: in the steady state case the last line of matrix A is set equal to: [1 0 0 0 . . . 0 1] so that vorticity at T.E. is 0:

gamma_1 + gamma_(N+1) = 0

I did not change this in the unsteady code only because I couldn't find any indication on how to modify it (the book says that for small oscillations the steady case condition should still work fine).

I would be really thankful if any of you could help me out, as it is more than a month that I am trying to solve the problem. I have also read many articles like Hancock and Mook's one, or similar, but they don't really say explicitly how to do things.

Thank you very much.

PS: do ask questions if you don't understand some of the procedures I have used, I will be happy to explain more thoroughly.

The unsteady Bernoulli eqn applies:

http://web.mit.edu/2.016/www/handout...ion_050921.pdf

(I'm sure there are better descriptions, but this is the first acceptable one I found).

Cp=1-V^2/2+d(phi)/dt

I strongly recommend that you try a wakeless problem such as a cylinder or ellipse. Oscillating in heave is an easier problem than oscillating in pitch.
Removing the wake eliminates the Kutta condition and convection as sources of error.

You should be able to reproduce the apparent mass calculations of Max Munk (See naca reports server for added/apparent mass).

Then use the unsteady Bernoulli equation to determine the time rate of vorticity shedding into the wake. Run a very low frequency - you should get steady-state results.

Describe the airfoil problem in detail: airfoil, frequency, amplitude, pitch axis.

The problem is I need it for airfoils. At the moment I am using a NACA 0012 and the experimental data with which I need to validate the code is this sinusoidal pitch motion:

5° + 5*sin(wt)

With k=0.2.

I'm not quite sure how to get d(phi)/dt out of my code, because I never handle the potential directly.

Sorry for the typo: Cp=1-V^2+2d(phi)/dt

I use the method of Luigi Morino, it applies to Oscillatory or unsteady or steady subsonic conditions. The time derivative is more significant as the frequency (which you haven't specified) increases.

If you have the vorticity on the airfoil surface, check me but isn't

phi=phi_at_leading_edge + integral[ gamma * ds]

The time derivative of phi_at_leading_edge will not contribute to the forces because it's integrated over a closed contour. So you can calculate forces and moments (but you haven't specified what output you are interested in, nor have you stated what cases, if any, you program does accurately calculate)

If you're not getting good results for an airfoil with wake, doesn't it make sense to simplify the problem to one which makes the error obvious? A high-frequency (k=1) wakeless case can determine if the boundary condition and Cp are correct. If you can't calculate a low-frequency (k=0.001) pitch case correctly, it's unlikely higher frequencies will work.

Your Kutta condition of equal upper and lower vorticity will not work for an oscillating airfoil which must be allowed to shed vorticity into the wake.

d(wake_circulation)/dt=-Vte*d(wake_circulation)/ds

where Vte is the mean velocity at the trailing edge.

In an article I have found a quite precise description of the Unsteady Kutta Condition for unsteady cases. The article is:

"Interaction of two-dimensional impulsively started airfoils"
(by WU Fu-bing, ZENG Nian-dong, ZHANG Liang ,and WU De-ming)

In this paper they state: (gamma_1)_k + (gamma_N+1)_k = (gamma_w)_k

(where gamma_w is the value on the extra T.E. panel of constant vorticity, and k indicates the time step)

They also use an iterative method to get delta_k (panel length) and theta_k (inclination of the panel) of the shed vorticity panel. This is then transformed into a concentrated vortex and builds up the wake.

I haven't implemented this yet but I have changed my code setting A matrix in such a way that at the line which expresses the Kutta condition I have:

[1 0 0 . . . 1 -1] which should be exactly what they suggest in the article.

Also I have found an error (which I have now corrected): i.e. I hadn't multiplied the gammas on the panel by the relative panel length to make it consisten with the wake vortices which are obviously concentrated vortices, not linear vortex distributions (this is especially important for the Kelvin condition)!

Things now are quite differente and Cp variations seem to be more realistic: for the mentioned case they go from 0 to approx -6.

Answering to a few questions posed by blackjack:

- I am interested in the Cp, above all;
- I still calculate it via the tangential component of the velocity at each collocation point, without considering the potential's derivative, also because I still don't understand how to (from the gamma or from other quantities) obtain it... should I integrate the gammas on the airfoil to obtain PHI, and then differentiate (= dPHI/dt == (PHI_k-PHI_k-1)/(t_k-t_k-1)) to get the second term of the RHS in "Cp=1-V^2+2d(phi)/dt"... or am I wrong?
- I have this file http://http://www.ar87.com/uni/19810...1981008464.pdf

and at page 7 you can see the data I am trying to get out of my code. So reduced frequency k is 0.2.

What I now get as variations of my Cp is this:
http://www.ar87.com/uni/dCp.jpg
which corresponds in the time line to these points:
http://www.ar87.com/uni/angle.jpg

(I just plotted those 4 points).

The link to the .pdf file does not work.

The scheme to find dphi/dt should work.

I have run your case. The time derivative at the leading edge is quite small - it contributes less than 0.003 to Cp but at the trailing edge contributes up to 0.09. This is the same magnitude as the boundary layer you neglect.

The most negative pressure on the airfoil I get is -4.85 without boundary layer. Less with the boundary layer.

Ok, the test case (experimental data) is this one (the link which wasn't working before):

http://www.ar87.com/uni/testCase.jpg


I have modified my code by adding a function which moves the vortices in the wake under the influence of the induced effect (panel induced velocity and wake self-induced velocity), and the same case produces this:

http://www.ar87.com/uni/dCp1WRUit.jpg

this is again the plot of 4 points on the sinusoidal motion: 0, pi/2, pi, 3/2pi.

I didn't manage to implement the dPHI/dt... how exactly should I calculate it?
What I did was something like this:

if j==1
PHI1(1) = PLEN(1)*(GAMMA(1,i)+GAMMA(2,i))/2;
else
PHI1(j) = PLEN(j)*(GAMMA(j,i)+GAMMA(j+1,i))/2 + sum(PHI1(1:j-1));
end

%dPHI/dt
dPHI(j) = PHI1(j) - PHI0(j);

PHI0(j) = PHI1(j);


where PLEN(j) is the panel length of the j-th panel, GAMMA contains the gamma distribution value at nodes of the panel, PHI should be the value of the potential I want to calculate. Could you help me out in getting it right?

Moreover, do you thing that it makes a great difference if I just use tangential component of the total velocity at each collocation point and calculate the Cp, rather than total velocity + dPHI/dt? (Because in the formula Cp=1-V^2+2d(phi)/dt, V is the module of the velocity, not the tangential velocity. Right?)

I don't know the pitch axis for your case: c/4, c/2? Moment is about c/4, c/2?

Your approach to dphi/dt sounds fine. Attached is the potential calculated for a 0012 at steady 10 degrees alpha. This should be useful for verifying that you have calculated the correct potential. The time derivative is then obtained from finite difference of two solutions.

Attachment here

Pitch axis is c/4, yes.

I'm afraid I'm not getting round to calculating the potential.
Thank you for the file, very kind of you. But still I don't get the correct values. I have calculated the 10° steady case for NACA 0012 (which is surely correct as I have correctly validate my steady state code months ago), however I don't know how to get the PHI out of my gammas.

Consider I get the gammas at nodes of the 100 panels on the surface. I have tried different calculations, but none give your results, which I presume are correct.

In Katz&Plotkin they suggest this formula for the potential (in the case of linear vortex distribution):

PHI = +/- gamma1/16*(x1^2 + 2*x1*x2 -3*x2^2)

where gamma1 should be the angular coefficient of the linear formula for the gamma:

gamma(x) = gamma0 + gamma1*(x-x1)


I also tried the formula I have written before, but that produces enormous values (10^24).
So at the moment I am a little bit confused.

My results seem to be quite correct, but not exactly like in the paper, and the minimum Cp I get is not -4.85 but a bit less.

I could send my code to anyone of you and you could give it a check... just for some suggestions... Thanks

Code Request
 

I running through some problems with my code
if Possible could you send me your code
my email : night_wing7@yahoo.com

potential calculation
 

I an dealing with linear strength vortex panel method.Can anyone of you can help in finding the potential for dphi/dt? If I am doing my programing in local coordinate system(frame fixed in wing as described in KATZ and PLOTKIN) can any body could suggest what should be the V for each panel for the formula Cp=1-V^2/2+d(phi)/dt that is I want to know V is the addition of which terms or simply the V as in steady state case as for i=1:npanel vel=0; for j=1:npanel+1 vel=vel+at(i,j)*gamma(j); end vel=vel+U*tang(i,1)+W*tang(i,2); Cp(i)=1-(vel/Q)^2; end where U=freestream X component velocity and W=freestream Y component velocity

Code request
 

Hi Anthony,

I've come across your post as I'm developing a 2D unsteady panel method for aerofoils too by following the Katz and Plotkin book. I'm doing this in a framework for gust reconstruction; do you still have your code available and if so could you please share it with me?

Thank you!

Best regards,
Simone.

Simone,

You haven't described what your program does correctly, so I propose the following course of action:

1. Add an input flag to skip the wake, so you can calculate flow without separation.

Skip the wake:
2. Calculate flow past circular cylinder in steady flow.
If the result is correct, proceed to next step.

3. Calculate flow past circular cylinder in unsteady flow (pitching or translation)
If result is correct, proceed to next step.

Activate wake:

4. Calculate flow past airfoil in steady flow.
What is lift coefficient vs. alpha?
If result is correct ...

5. Calculate flow past pitching airfoil in low-frequency oscillatory flow.
Reduced frequency should be 0.001 < k < 0.01
Answer should depend on angle of attack the same as steady flow.
If result is correct ...

6. Calculate flow past pitching airfoil for 0.01 < k < 0.5
Calculate amplitude and phase of lift - compare to Theodorsen function.

Which is the first step that the result is not correct?

The pressure coefficient in unsteady flow is Vs^2 - V^2 + d(phi)/dt
where Vs is the surface velocity, V is the fluid velocity w/r to the surface, and d(phi)/dt is the time derivative of the potential. It's trivial to evaluate if you use Luigi Morino's method. The last term provides the added mass term.
Added mass is significant to hydrodynamic analysis and aerodynamic flutter.

Hi, all
 

Can we calculate d(phi)/dt as follows?

1. Calculate the equivalent lumped vortex at the cell centre Gam = 0.5*(gamma1+gamma2) * dl

2. Calculate the potential of the point vortex phi = -Gam/(2pi)*atan(dz/dz)

3. d(phi)/dt = (phi_(ti)-phi_(ti-1))/dt