Hi!

I am writing a 2-D Euler structured solver in frequency domain using "harmonic balance" method. One of the references I am using is

"Thomas, J. P., Dowell, E. H., and Hall, K. C., "Further Investigation of Modeling Limit Cycle Oscillation Behavior of the F-16 Fighter Using a Harmonic Balance Approach," AIAA Paper 2005-1917, 46th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials (SDM) Conference, Austin, TX, April, 2005."

I am testing this code for the AGARD test case CT1 (pitching NACA0012 airfoil), but I am not getting good match with the experimental values. I know it should match from reading other publications implementing this method. Has anyone else implemented this method before or is working in a related area? I will be glad to share snippets of code if it helps me and other users resolve similar problems.

Manish

Are you actually performing the iterations in the frequency domain, or are you just evaluating the time derivatives in the frequency domain and then iterating in the time-domain? The harmonic balance method can be used both ways, but I think the latter treatment is preferred. The reason, why I am asking this, is that it would be relatively easy to build a dual-time marching method, if your iteration is already performed in the time domain. You should first validate your Euler code by a time-marching method before applying the more sophisticated HBM.

I am evaluating the time derivatives in frequency domain and then iterating in time domain. My Euler code is validated for the dual time marching method. My dual time implementation (time domain) gives very good agreement with the experimental values for AGARD CT1 test case. My results for harmonic balance method show similar trend but I was expecting a better match based on previously published results. I tested the harmonic balance code with 1 harmonic for steady flow case and for omega = 0 case and it gives good agreement. Somehow for oscillating case it doesn't match. I would like to send you my current results if you would like to take a look. It would be great if you can you share a sample code for how you implement harmonic balance using dual time marching?

If you have validated your time domain code against the experimental results then if you are implementing your harmonic balance method correctly your results should converge to these as you increase the number of harmonics you use. Do you have a converged result with respect to the number of harmonics?

For harmonic balance case any more than 5 harmonics had no effect on the solution for the case I described. I tried upto 15 harmonics but they fell on the same curve. I tried the following which might give some clue. I Fourier transformed my lift coeffcient values for both time domain code and frequency domain code and compared the mean and first cosine and sine components as higher harmonics were negligible. I noticed that the mean and first sine (1s) components were close for both methods, but there was a discrepancy in the first cosine component. The values were as follows:

Time domain: 0 -> 0.4131, 1c -> 0.0504, 1s -> 0.2725 Frequency domain: 0 -> 0.4120, 1c -> 0.0158, 1s -> 0.2657

I calculated the lift curve using frequency domain coefficients by replacing the 1c component with the time domain 1c component and the results agreed very well. So, I know the discrepancy in 1c component is what I should look at. But, I am not sure where to go from here and which part of the code is the culprit for this.

Hmmm..check to make sure your fourier transform matrices(matrix E and it's inverse I believe is the way that they are represented in Jeff's papers) are setup correctly. I can give you representative code if you need it.

I got the Matrices from Jeff and checked them just now, and my matrices matched perfectly with his. I even tried putting his matrices directly in my code, but with no change in result. I think I have hit a wall and now don't know which way to proceed. It will be a great help if you can share your code with me. My email is manishk@gmail.com Thanks a millions, Manish

Well my code for the fourier matrices is going to be the same as Jeff's since I got it from him when I was at Duke. What happens if you drop the amplitude of the pitching? Do the results start to agree better? We may want to take this offline. Feel free to e-mail me to continue the discussion.

There are a few comments to add on convergence. As Peter said, you should expect your harmonic balance result to converge to the time marching result as the number of harmonics is increased (for a case like this I would use no more than 10). But there is one condition: Your time-marching method should also be accurate in the sense that additional temporal resolution (increasing time steps per period, or using higher order time difference) will not change the result. Make sure to compare a converged time-marching result to a converged harmonic balance result. Otherwise you cannot draw any conclusions, because both methods will have different discretization errors. I am assuming you did just that, and you still see significant discrepancies, where the time marching matches the experiments but the harmonic balance is off. In that case it's pretty safe to say that there is some problem in your implementation of the HB method.

The major difference between time-marching and HB is the evaluation of the time derivatives, so that's the part in the code that you should focus on. All you need to do for HB is to (1) take the Fourier transform of all flow variables at all grid cells, (2) take their analytical derivatives in Fourier space, and then (3) take the inverse Fourier transforms of those derivatives (as opposed to forming the finite time differences for the dual-time method). So I see three likely sources of error, namely the above three operations.

To test your code you can initialize the whole flow field (in space-time) with a known solution, for example let all flow variables be sinusoidal functions of time with a specified frequency and amplitude. Without applying boundary conditions or doing anything with the flow, simply use your subroutines to get the time derivatives of this initial condition and compare them with the known exact solution for the time derivatives.

Mani,

You said "simply use your subroutines to get the time derivatives of this initial condition and compare them with the known exact solution for the time derivatives".

Isn't this a test for just the time derivative matrix which is referred to in the paper as omega*D matrix. Anyway I did that and came up with exact agreement with the analytically calculated velocities. Basically what I did was just choose one cell. Assume the flow field was given by 5+4*Sin(omega*t). So the analytical velocity is 4*omega*Cos(omega*t). I initialized the cell value at 11 time steps. Then I multiplied my omega*D matrix with that flow field vector and calculated the velocities. They turned out to be an exact match.

I am using a chracteristic based boundary condition (using Riemann Invariants) for the farfield which is the same as implemented in my steady flow solver and my unsteady dual time solver. So, I am not using the nonreflecting boundary condition formulation as described in other HB papers. Do you think that might have an effect on the solution that somehow those boundary conditions are not suitable for the HB method?

>Isn't this a test for just the time derivative matrix which >is referred to in the paper as omega*D matrix. Anyway I did >that and came up with exact agreement with the analytically >calculated velocities. Basically what I did was just choose >one cell. Assume the flow field was given by >5+4*Sin(omega*t). So the analytical velocity is >4*omega*Cos(omega*t). I initialized the cell value at 11 time >steps. Then I multiplied my omega*D matrix with that flow >field vector and calculated the velocities. They turned out >to be an exact match.

Did you do this test with your CFD solver or with a separate code? It should be done with the actual solver.

The boundary condition is irrelevant. If it works with time-marching, it should work with the harmonic balance method. I use the same boundary conditions for either case.

One more thing you can try: Instead of using the Fourier transform (or matrix), just apply a finite difference to obtain the time derivatives, solving the whole period simultaneously. Let's say your velocity is v(x,t), define dvdt(x,t) = (v(x,t) -v(x,t-1))/dt for first order. For the first time level in the period you need to apply it in a cycled way, i.e. dvdt(x,t1) = (v(x,t1) -v(x,tn))/dt, where t1 is your first time level and tn is the last level in the period. You can use the same difference that is used in your dual-time application, and in that case your solution should be absolutely identical to the (converged!) time-marching solution! If that's not the case, then your problem is not related to the evaluation of the time derivatives. (did you take care of grid velocities?)

Mani, Thanks for your response. I calculated the grid velocities with the actual solver. I had already tested the calculation of grid velocties both ways, using finite difference and using omega*D matrix. My final solution remained unaffected. So, I am sure the problem is not with the time derivatives. Do you have any results (CL vs alpha) for a simple oscillating airfoil using your harmonic balance code? I would like to get some idea on what kind of results I should expect. Since the discrepancy in my results is in out-of-phase mode, there could be a problem in some nonlinear coupling. In your opinion which portion of the code could be responsible for this?

Manish