# Linear-Strength Vortex Unsteady Panel Method (2D)

 Register Blogs Members List Search Today's Posts Mark Forums Read

 April 3, 2012, 13:44 #2 New Member   james nathman Join Date: May 2011 Posts: 22 Rep Power: 8 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.

 April 3, 2012, 15:41 #3 New Member   Anthony Join Date: Jul 2011 Posts: 9 Rep Power: 8 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.

 April 4, 2012, 13:52 #4 New Member   james nathman Join Date: May 2011 Posts: 22 Rep Power: 8 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.

 April 4, 2012, 14:15 #5 New Member   james nathman Join Date: May 2011 Posts: 22 Rep Power: 8 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.

 April 4, 2012, 16:17 #6 New Member   Anthony Join Date: Jul 2011 Posts: 9 Rep Power: 8 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).

 April 5, 2012, 14:04 #7 New Member   james nathman Join Date: May 2011 Posts: 22 Rep Power: 8 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.

 April 6, 2012, 12:34 #8 New Member   Anthony Join Date: Jul 2011 Posts: 9 Rep Power: 8 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?)

 April 6, 2012, 14:45 #9 New Member   james nathman Join Date: May 2011 Posts: 22 Rep Power: 8 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.

April 6, 2012, 14:49
#10
New Member

james nathman
Join Date: May 2011
Posts: 22
Rep Power: 8
Attachment here
Attached Files
 phi10.txt (4.3 KB, 28 views)

 April 6, 2012, 17:26 #11 New Member   Anthony Join Date: Jul 2011 Posts: 9 Rep Power: 8 Pitch axis is c/4, yes.

 April 7, 2012, 05:41 #12 New Member   Anthony Join Date: Jul 2011 Posts: 9 Rep Power: 8 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.

 April 11, 2012, 11:00 #13 New Member   Anthony Join Date: Jul 2011 Posts: 9 Rep Power: 8 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.

 April 13, 2012, 05:46 #14 New Member   Anthony Join Date: Jul 2011 Posts: 9 Rep Power: 8 I could send my code to anyone of you and you could give it a check... just for some suggestions... Thanks Last edited by Sonik_87; April 13, 2012 at 06:13.

 June 30, 2012, 17:56 Code Request #15 New Member   Mido Join Date: Mar 2011 Posts: 22 Rep Power: 8 I running through some problems with my code if Possible could you send me your code my email : night_wing7@yahoo.com

 February 4, 2013, 14:48 potential calculation #16 New Member   Pappu Join Date: Feb 2013 Posts: 1 Rep Power: 0 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

November 4, 2015, 13:11
Code request
#17
New Member

Simone Simeone
Join Date: Nov 2015
Posts: 1
Rep Power: 0
Quote:

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.

 November 5, 2015, 14:49 #18 New Member   james nathman Join Date: May 2011 Posts: 22 Rep Power: 8 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.

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post subash OpenFOAM 0 May 29, 2010 01:23 Gerrit Groot Main CFD Forum 17 September 17, 2007 09:13 gocarts Main CFD Forum 2 April 22, 2007 22:50 Daniel Main CFD Forum 0 September 15, 2006 07:51 Daniel Main CFD Forum 5 January 6, 2003 09:09

All times are GMT -4. The time now is 16:31.