Hello Folks,

I've been working on the numerical solution of the transonic small disturbance equation and have recently discovered that depending upon the way in which I non-dimensionalise my equations I get a different solution.As far as I am aware this should not be the case - if someone could please advise me on this matter I'd be very grateful.

I have attempted to clarify the problem below. I must apologise for the length of my post but as you will all understand I want to make my problem as clear as possible - after all I am hoping someone will be kind enough to help me! )

I am trying to solve the problem of transonic 2d flow around a thin airfoil which is undergoing prescribed harmonic pitching motion. The airfoil angle of attack is prescribed as: a=a0*sin(W*T) where W is the cyclic frequency and T is dimensional time. In all that follows I have taken upper case letters to be dimensional variables and lower case to be their corresponding non-dimensional form.

The TSD equations are quasi-linear (in x). The non-dimensionalised, unsteady TSD equation can take one of the following forms (depending upon how we non-dimensionalise the time 'T' term when deriving the equation).

------------------------------------------- OPTION #1: (NON-DIMENSIONALISE USING FREQUENCY)

TSD BECOMES: k*p_xt + (p_x -constant)*p_xx - p_yy = 0

Here '_xt' denotes a second derivative with respect to space and time once each; '_xx' corresponds to second derivative with respect to x and so on. The non-dimensional time is scaled according to the cyclic frequency 'W': t=W*T. The term 'k' is the reduced frequency equal to WC/U (C=chord length and U=upstream velocity, both are taken equal to 1 for simplicity).

I choose an increment of, say, dt=0.01 as the non-dimensional timestep between time levels n and n+1. I evaluate the solution of the equation for 2e4 timesteps.At each timestep the angle of attack is prescribed using the non-dimensional time values as: a=a0*sin(WT)=a0*sin(t) {due to the way in which time has been non-dimensionalised}. The code produces values of aerodynamic moment and lift acting upon the airfoil for every timestep. Prior to plotting graphs I re-scale back to physical time by dividing the time axis values by the cyclic frequency W => T=t/W.

------------------------------------------- OPTION #2: (NON-DIMENSIONALISE WITH VELOCITY & CHORD)

TSD BECOMES: p_xt + (p_x -constant)*p_xx - p_yy = 0

This time the temporal term is scaled using t=(U/C)*T which removes the reduced frequency term from the equation. The timestep is, again, taken as dt=0.01. The characteristic velocity & length scales are both taken equal to 1 in order to simplify things and to make comparison with option #1. The airfoil motion is prescribed, as before according to a=a0*sin(WT) = a0*sin(W*(C/U)*(U/C)*T) = a0*sin(wt)

The solutions for moment and lift have once again been plotted vs time. In this case there is no-need to rescale as the time-scaling factor is 1 exactly.

My problem is that there are large differences observed in the two solutions The amplitudes and phase values of the plots differs between those obtained using option #1 and #2. I was expecting that the solutions would, once suitably re-scaled back to dimensional co-ordinates, be virtually identical. There are only two temporal terms used in the method: when I calculate the derivative 'p_xt' and when I prescribe the airfoil motion. I expect the problem is in the way I have dealt with these terms although I am not sure where. I also have suspicions that the problem may be caused by the nonlinearity of my equation (in x) and scaling factor of the p_xt term.

If someone could offer some advice then I'd greatly appreciate it. If not then thank you for reading this far in my rather lengthy post!

Bren

Are both solutions converged with respect to the timestep used? Also what happens if you use a cyclic frequency W=1?? Do the solutions compare favorably for this?

Hi Peter,

I have tried the W=1 case and both sets of solutions are identical. When W is not equal to one I find that the solutions are dependent upon the way in which I non-dimensionalise the time variable. The two different methods provide solutions with different amplitudes and phase-offsets.

The code is converging in the sense that the output values for moment and lift is periodic (which is what I expect). If I run the unsteady algorithm to solve a steady flow (by, say, fixing the angle of attack) then it does not matter how I non-dimensionalise the time variable - I always get the same solution.

Many thanks,

Bren

The fact that when W=1 the solutions are the same is the answer to your problem...What I meant by converged was does your solution converge to one answer when you vary dt from 0.01 to 0.005 to 0.001,etc. The problem has some characteristic time scale which you must resolved. In the case of the oscillating airfoil with prescribed motion this will be approximately 2*pi/W. In your first non-dimensionalization by taking dt=0.01 your are basically taking 100 steps per period of oscillation. This is different than the second non-dimensionalization where you are saying your characteristic time scale is C/U. When you are using dt=0.01 for this you are taking 100 steps based upon this. If W=2*pi then these are the same and you will get the same answers. Otherwise you won't. For the nonlinear fluid-structure interaction problem you are working on the non-dimensionalization usually used is the second one since you don't know the reduced frequency apriori.

Hi Peter,

The solutions I have calculated are both converged. If I reduce the timestep to 0.001 or even 0.0001 then I get the same 2 solutions; they just takes a lot longer to compute.

I don't understand your reasoning behind your timestep explanation. If I have prescribed a frequency of, say, 0.2 and define dt=0.01 and use the same (U/C)=1 as before, then my 2 options of non-dimensionalising are:

Frequency: t=0.2*T => a=a0*sin(T) and the number of timesteps in one oscillation is 2*pi/dt.

Velocity/Chord: t=(U/C)*T => a=a0*sin(wT) and the number of timesteps in one oscillation is, 2*pi/(w*dt).

My understanding if the timestep is such that both solutions are converged then by rescaling back to physical co-ordinates should provide the same solutions. The systems are undegoing the same oscillation at the same rate and so I don't understand why the solutions are different.

I'm possibly missing something very obvious here! Many thanks, Bren

Hmm..ok. I thought maybe in one case you were resolving well enough. But if BOTH methods show convergence with respect to dt then that is not the case. So the only difference in the solutions is the phase and amplitude? Do both give the correct frequency behavior? ie the period of CL should be nearly 2*pi/W. Also do you get better agreement when the airfoil motion is small which you should expect to give a linear response.

Hi Peter,

Both methods are giving the correct frequency behaviour. I have used a spacing of dt=0.01 using each choice of non-dimensionalisation. There is a phase and amplitude difference between the two sets of results. The angle of attack is also very small (only 0.5 degrees). If I lower it further then I still get differences between the two methods.

The 2 codes are virtually identical - the only differences are in the definition of the equations (i.e whether or not the 'k' term appears) and the prescription of the angle of attack (i.e. sin(T) vs sin(wT) depending upon the scaling.

Thanks Peter, I really appreciate your help in this matter

Bren

When you have eliminated the impossible, whatever remains, however improbable, must be the truth.

Just kidding! Anyway, here is what my intuition says. The dimensional freestream velocity you have taken to be unity, without loss of generality. However, then the nondimensional velocity must differ in the two nondimensionalizations. In the second option, the nondimensional velocity will of course also be unity, since you select the freestream velocity as the reference speed for that option. In the first option, however, you choose length and time scale reference values. You are then NOT free to pick the (dimensional) reference velocity, which instead must be calculated as V = C*W.

I believe that this reference velocity must enter your formulation in some way that you have not recognized yet. Perhaps as a freestream boundary condition. Even more likely is that it enters into the nondimensionalization of pressure, which is normally p = P/(RHO*V*V), where RHO is the (dimensional) reference density. I would check your derivation of the TSD to trace the influence of the reference velocity. You will notice that your TSD as stated is not dimensionally homogeneous. The second term is dimensionally quadratic in P (while the other terms are linear) and dimensionally inversely cubic in C (while the other terms are inversely quadratic). It has been many (many!) years since I looked at the full-potential transonic flow equations, but I imagine that the factor (p_x - constant) arises from the convective term of the Euler equations and has not been correctly nondimensionalized. Most textbooks will not bother to give you nondimensional equations; their listed equations are usually dimensional. This would explain the difference you see in amplitudes. The phase error I do not have any clue about. You have correctly accounted for the time scale differences between the two options, and therefore your adjusted frequencies match.

Thanks very much for your advice Ananda, I really appreciate your guidance in this matter.

The TSD equation I am solving is derived from the non-dimensional Euler equations using the small disturbance approximations. An asymptotic approach is used to derive the equations of motion based on a small parameter delta (which here is taken as the airfoil thicknessratio). Eventually the analysis provides (again using lower case for dimensionless terms): du/dy=dv/dx as one of the equations (where u and v are dimensionless velocity terms in the x & y directions, and the x & y terms are themselves dimensionless). This equation suggests that we can introduce a velocity potential p such that p_x=u and p_y=v. The potential is dimensionless as are it's derivatives with respect to x,y and t. The constant term that appears in the equations is actually a combination of dimensionless coefficients. The TSD equation should therefore be completely dimensionless: p_xt + (p_x -constant)*p_xx - p_yy = 0

I may have missed something though! I have been through my derivation several times and am pretty sure that I haven't dropped any dimensional terms anywhere. The boundary conditions in the problem are also stated in non-dimensional form. My equation & BC's also agree with those found in the literature.

I think I understand what you mean when you talk about the velocity terms. As far as I am aware though I am free to choose the freestream velocity as I choose. The freestream velocity doesn't even appear in my code as a variable and so I don't think that the problem can lie there. The values of lift and moment coefficient don't rely upon the freestream velocity - I have calculated them directly in terms of the potential and so I don't think that the freestream velocity scaling will affect the amplitude of the plotted solutions..... it really is a weird problem!

It is driving me crazy...er )

Many thanks, Bren

Sorry Ananda I've just realised that what I said earlier wasn't technically correct.

While the velocity scaling doesn't appear directly in any of my equations it does appear indirectly in the reduced frequency term

k=W*C/U

The problem is that if I take U=C*W to be my velocity scale then I always get a value of k=1. I'm therefore unable to account for the frequency in my model.

Yes, my mistake(s). I took "p" to be a nondimensional pressure, and further, I misidentified your "Option 1" with a third option which would have been the way I nondimensionalized things had I based my reference time scale on the pitching frequency. This third option is just as valid as the other two, though, and it is NOT true that you would be unable to account for frequency in that model. But let us ignore this third option, and focus instead on the two you listed.

First of all, let me say that the flaw in your stated TSDE lies not in the second term, but in the third term (-p_yy). My guess is that there is a coefficient to that term that is unity when nondimensionalizing one way but would be frequency (or velocity) dependent when nondimensionalizing the other way. However, I cannot prove this unless I derive the TSDE from the dimensional Euler equations, and I do not have the time to do this, sorry. As I said in the previous post, your stated TSDE is not dimensionally homogeneous, even in light of the fact that p is a velocity potential.

Let me illustrate what I expect out of a consistent formulation. Let us follow your convention of uppercase letters representing dimensional quantities and lowercase representing nondimensional ones. In each case, let L denote the reference length scale, S the reference time scale, and Q the reference speed scale. L, S, Q are all dimensional quantities. Further, let F denote the freestream speed, and C the airfoil chord. Then

X=xL, Y=yL, T=tS, U=uQ, V=vQ ______(E1)

Let me pick a fictitious dimensional TSDE, which looks like a hybrid of your TSDE and the Euler x-momentum equation:

U_T + U*U_X - U*V_Y = 0 ______(E2)

with an angle-of-attack given by

a = a0*sin(WT) ______(E3)

where W (actually W/(2*Pi)) is the pitching frequency. The virtue of (E2) is that it is dimensionally homogeneous, as is (E3). Being dimensionally homogeneous is necessary for the equation to be valid when measurements are made in arbitrary consistent unit systems.

YOUR OPTION 1 FOR NONDIMENSIONALIZATION:

We pick L=C, Q=F, S=1/W. Then [using (E1)], (E2) and (E3) become

k*u_t + u*u_x - u*v_y = 0 ______(E4)

a = a0*sin(t) ______(E5)

where k = C*W/F is the reduced frequency of pitching (which shows that what essentially matters is the ratio of the pitch period to the C/F time period, and not the absolute value of the pitch period).

YOUR OPTION 2 FOR NONDIMENSIONALIZATION:

We pick L=C, Q=F, S=L/Q=C/F. Then [using (E1)], (E2) and (E3) become

u_t + u*u_x - u*v_y = 0 ______(E6)

a = a0*sin(k*t) ______(E7)

If you pick the same C and F for the two codes (you picked C=1 and F=1), you should get the same answer with both codes after accounting for the differences in S. Keep in mind that the velocity potential P=p*Q*L, which would make a difference if you went with the third option for nondimensionalization. That is it. Should be as simple as the above.

Now, just for fun, you can try to work out what your "nondimensional" equation looks like if you started out with a dimensionally inhomogeneous (and therefore incorrect) equation such as

U_T + U*U_X - V_Y = 0 _______(E8)

Also for fun, you could try to transform (E2) and (E3) with a third option for nondimensionalization (you are right, k disappears from the time term, but do you know where else it pops up?):

L=C, S=1/W, Q=L/S=W/C.

Hope this clears things up. It is as simple as I can make it without actually deriving the TSDE with a small parameter expansion.

Hey guys,

I finally fixed the problem. I had made a simple mistake in my code that essentially acted as an additonal re-scaling of the problem. It was such a small thing I never even noticed it when repeatedly checking through the code. The two different scaling methods now give solutions that are identical. I still have no idea why the rescaling caused a phase change in the results but it's working now so I'm happy )

Thanks Peter and Ananda - your advice proved invaluable.

Bren

When it comes to writing code..the problem is always something simple...especially after you have figure it out!!! : ) I am glad you fixed it.

If you don't mind me asking are you working at a company or is this graduate work and if so what school do you attend?

Glad you caught it. It is a good thing you tried the two different nondimensionalizations, or else you might never have stumbled across the error. I am still puzzled though as to how that third term in your TSDE does not throw the answer off under different nondimensionalizations. As long as you are getting consistent results, I guess what you've put in the code must be okay now.

Hi Bren, The solutions you get out of the code are dependent on two factors only. The frequency, and the amplitude of the prescribed pitch motion.

For the first case, t = WT. Hence the frequency of the prescribed pitching motion is 1. Thus, for dt = 0.01, after 200*pi timesteps, the prescribed motion should repeat itself. As you rightly pointed out, when the code has been computed the aero load as a function of t, the t-axis values are multiplied by (1/W), to get a T-axis. k should still be entered as W.

For the second case, w = W, bcos U/C = 1. Clearly, the frequency of the prescribed pitching motion is W. Thus, for dt = 0.01, after 100*(2*pi/W) timesteps, the prescribed motion should repeat itself. When the code has computed the aero load as a function of t, the t-axis values are multiplied by (C/U), to get a T-axis.

The most likely source of the problem is the definition of the reduced frequency, k. k = W*C/(2*U) and nu = W*C/U. Hence, if a value of k is assigned in the first nondimensional case, Try W = 2*k in the second nondimensional case.

Good Luck

Hi Ananda, I've a VERY limited amount of experience with the Transonic Small Disturbance Theory (and only in steayd flows) but I'm pretty sure that the equation stated is the one I've seen in the literature. I'm unsure as to what you mean by "dimensionally homogeneous". Would you mind giving a brief explanation as to what this term means and why it is important for an equation to be such? I'm afraid I'm still fairly new to this CFD business and some of the theory is new to me.

If the term means what I think it does then surely many equations are of this form, for example the linear wave equation: du/dt + (du/dx) = 0 would surely be such an equation?

S

Hi Bren, Good to hear you've cracked that problem. One thing though if you don't me asking. How did you non-dimensionalise the TSD equation to get

p_xt + (p_x - const)*p_xx - p_yy = 0

or k*p_xt + (p_x - const)*p_xx - p_yy = 0

The non-dimensional forms of the TSD i keep stumbling upon are of the form

p_xt + (((1+Y)/2)*p_x - (1-(M^2))/(2*(M^2)))*p_xx - (1/2)*(M^2)*p_yy = 0

This makes sense, bcos the TSDE must reflect the dependence of the unsteady flow on the Mach number. Also, I realised that you solve for the harmonic components of the lift and the moment coefficients as a function of discrete time from the two formulations you presented. As a suggestion, you could try solving for p in the frequency domain for prescribed periodic motions, by setting p = ps + pu*(e^(i*k*t)) in the unsteady TSD equation and solving two TSD equations for the steady component of p, ps and the first harmonic of the unsteady component of p, pu.

All the best Nick

Hi Nick, I avoided including a few terms in my TSD equation because they weren't really relevant to the question I was asking (and they make the equation look scary!). The actual equation I've derived is written:

[2*k/tratio^(2/3)]*phi_xt = [(1-Minf^2)/tratio^(2/3) - (Y+1)*phi_x]*phi_xx + phi_y'y'

which is pretty nasty-looking ) The terms I avoided in my original post are dimensionless constants (involving either the heat ratio Y or the mach number Minf) and so weren't important when it came to the non-dimensionalisation problems I was having. The above equation is in similarity form but it's validity can be improved by changing it slightly to include additional mach number terms:

[2*k*(Minf^2)/tratio^(2/3)]*phi_xt = [(1-Minf^2)/tratio^(2/3) - (Y+1)*(Minf^2)*phi_x]*phi_xx + phi_y'y'

These additional mach number terms provide better agreement between the numerical solutions and experimental results although the similarity property is lost.

I'm not familiar with the particular form of the TSD equation you are using although as you can see it's very similar to the one i'm solving. I suggest that the main source of the differences is due to the way that the equations have been derived. I have seen a few different forms of the TSD and the differences are usually due to the way in which the equation has been derived and the underlying assumptions (axial symmetry, small parameter expansion etc). The traditional way to derive the equations is from the full potential equations & flow assumptions (such as isentropy etc) and then to eliminate small terms to arrive at the TSD equation. The way I have derived the equation is by using a small parameter expansion, with the thickness ratio as the small parameter. This involves a bit more work and lots of simple, but annoying, algebra but seems a lot less hand-wavy than some other methods.

The constant that appears before your phi_yy term is probably due to the way that the equation has been scaled. The perturbation field caused by the airfoil grows rapidly in the y-direction as the flow approaches mach 1 (due to the angle of the sonic cone) and so the y-coordinate must be rescaled in order to correctly account for the flow behaviour. In my equation this is represented by using the phi_y'y' term, where y' is rescaled using the thickness ratio. In your equation the y-derivative is multiplied by the mach number, essentially the two equations say the same thing. The TSD equation is not uniquely defined and there are several different forms of it. This can make the literature very confusing although in general the differences are small.

I had previously thought about using the steady flow solution and introducing an oscillating perturbation to the flow as you suggest. It will be difficult to include within my current solver, however. I will possibly give it a try in the near future though, if I do I'll let you know how it works out.

Many thanks, Bren

Yes, you are right, the linear wave equation with unit advective speed du/dt + du/dx = 0 in a euclidean spacetime is such a dimensionally inhomogeneous equation. For that reason, if you consider the equation to be in its nondimensional form, and solve it, the nondimensional solution u = f(x,t) needs to be mapped with care into the (X,T,U) dimensional space. There is of course no difficulty with mapping u to U, since the equation is dimensionally homogeneous with respect to the occurrence of u. However, the relationship between the (x,t) plane and the (X,T) plane is where the ambiguities arise. Of course, in this particularly simple case, the researcher is able to easily see what must be scaled and how it must be scaled, but in more complex equations this may mean that you are not solving the equation you intended to solve. In the particularly simple nondimensional form above, the researcher has effectively specialized the physics to consider advection along x = t+c lines in the (x,t) plane. Equivalently, the researcher has picked the velocity scale to give unit coefficient to the second term, and thus picked the relationship between the space and time scales, leaving only one of them to be picked freely.

One may view the wave equation above as du/dt + 1*du/dx = 0, which is a dimensionally homogeneous equation if we treat the "1" as having the numerical value of unity, but having dimensions of velocity (of dx/dt). Effectively, this is choosing the velocity scale, and as mentioned in the previous paragraph, leaves us the free choice of either the spatial scale or the temporal scale, but not both, if we wish the solution to apply to all linear advection situations in the dimensional (X,T) space. Alternatively, if we freely pick both the space and time scales, then one and the same nondimensional solution u(x,t) represents solutions of differing physics in the (X,T,U) space, namely advection along parallel lines of differing slope which depends on the choice of length and time scales. Provided that the researcher keeps all this clearly in mind, he/she can get away with this approach to nondimensionalization, as must be the case with Bren's TSDE and other equations.

A better way out of all these ambiguities is to solve instead the dimensionally homogeneous (linear) advection equation du/dt + a*du/dx = 0, where "a" has the dimensions of velocity. Then one can freely pick both space and time scales(which then determine the velocity scale) so as to easily set initial/boundary conditions and frequencies, and also freely pick the physics, namely, the angle at which the initial conditions are advected (which amounts to setting the value a = A*S/L, where S is the time scale, and L the length scale).

You can see from Nick and Bren's posts of 19 Aug that the TSDE is more complex than initially presented, and that Bren's y direction has been scaled with a different reference length than the x direction, which is how the airfoil thickness ratio enters into the picture. I am unwilling to hazard a guess as to exactly how Bren's nondimensionalization leads to the final TSDE used, but starting from the dimensionally homogeneous Euler equations, there must be some approximations such as treating the Mach number as being approximately unity (so that 1 - M^2 is approximately zero) and rescaling with the thickness ratio. These choices may correspond to either picking a time scale or a length scale in a particular direction, or to picking a particular subset of the physics as applicable. But that is how the TSDE as stated ends up being dimensionally inhomogeneous.

All of this post is just my thoughts on nondimensionalization, and I cannot attribute it to any published paper or book, though some author may have published a similar explanation for students, if my views are not mistaken. The main reason that I find dimensionally homogeneous equations more comfortable to work with than inhomogeneous ones is that the scaling is more transparent in the former (making it easier to go back and forth between dimensional and nondimensional quantities), and I am less likely to be puzzled by what particularization or approximation has been made in the physics.

Dear Ananda, I followed the posts exchange between you and Bren with interest. I understand you are a full potential equation expert. I then would like a little advice from you. I programmed a Finite Element code to solve the full potential equation (steady) on unstructured grid. I am having lots of troubles when the shocks get strong, ie the solution diverges. I am using density upwinding but still handling the nonlinearity updating the density with the solution at the previous step. I understood that this cannot work since I violate the flow physics when the equation gets hyperbolic. What is the best way out of this problem, in your opinion? I was thinking about solving the unsteady equation to steady steate implicitly in time, but how can I linearize it? Any of your suggestions would be of great help! Thanks a lot Alberto

PS Anyone can answer, of course