Separation Modeling

[antoniosgeme]

Hello Cam,

This is an awesome project. I was wondering at what flapping frequency do you expect the unsteady VLM solver to stop giving accurate results? Have you attempted any validation on the force outputs? I would imagine that at higher frequencies you'd need some leading edge vortex shedding to capture what is happening (but maybe not?!)

[camUrban]

Hi Antonios,

Thanks for the kind words about my project. I'm familiar with some of the awesome work you and Dr. Jones do at UMD, so it means a lot coming from you!

I have done some work to validate the solver, but more work in this area would be valuable! Check out /validation/validation_paper.pdf for a write-up comparing Ptera Software's simulations to wind tunnel data of a flapping-wing robot. The paper is a bit out of date (many bugs have been found and fixed since). However, you can try the simulation with the most up-to-date code by running /validation/validation_study.py. Here's an example of the output:

[camUrban]

You also make an excellent point regarding LEVs and separated flow! Many animals use stable LEVs to enhance their force production during flapping-wing flight [1]. Other forms of flow separation are probably prominent as well. I am extremely interested in incorporating these phenomena into this project's solvers!

I have read a few papers that attempt to do something similar [2][3][4]. Unfortunately, I don't think their solvers are publicly available or open-source. Additionally, I'm not sure I have the expertise required to vet their approaches. Given your past work, I'd love to hear your opinions on if/how we could incorporate LEVs and separation into Ptera Software.

[1] Chin, D. D., & Lentink, D. (2016). Flapping wing aerodynamics: from insects to vertebrates. Journal of Experimental Biology, 219(7), 920–932. https://doi.org/10.1242/JEB.042317
[2] Roccia, B., Preidikman, S., Massa, J., & Mook, D. (2013). Modified unsteady vortex-lattice method to study flapping wings in hover flight. AIAA Journal, 51(11), 2628–2642. https://doi.org/10.2514/1.J052262
[3] Nguyen, A., Kim, J., Han, J., & Han, J.-H. (2016). Extended unsteady vortex-lattice method for insect flapping wings. Journal of Aircraft, 53(6), 1709–1718. https://doi.org/10.2514/1.C033456
[4] Levin, D., & Katz, J. (1981). Vortex-lattice method for the calculation of the nonsteady separated flow over delta wings. Journal of Aircraft, 18(12), 1032–1037. https://doi.org/10.2514/3.57596

[antoniosgeme]

Great job on the validation study. You bring up some interesting points on the comparison with experiments.

Indeed including LEV effects would be awesome. Here are some thoughts:

I do not think there is a standard "accepted" way of including leading edge separation in 3D potential flows. The simplest and perhaps least rigorous way is what ref. 3 did (from the references you cited). Basically decide whether the flow is separated by comparing the effective angle of attack (AoA) to your static stall AoA. Then take the suction force (tangential to the wing) and rotate it by 90 degrees and recompute lift. This is based on Polhamus' theory of LEV lift and I'd love to see how well it works but personally I am not convinced of its theoretical validity.

The next level of fidelity is to simply to shed vorticity at the leading edge by imposing the Kutta condition on it. It sounds funky but it actually gives decent results (in 2D as far as I know) and it is simple to implement numerically. You can even have an off/on switch again based on your effective angle of attack. However, this method is agnostic of the exact geometry of the leading edge, which as you can imagine plays a pretty significant role on the separation dynamics. This leads us to the next paradigm.

The "state-of-the-art" right now for leading edge separation in vortex models is the Leading Edge Suction Parameter (LESP). The idea is that for some given operating conditions and geometry, the leading edge can hold some critical amount of suction which is motion independent. If the suction exceeds that critical value, vorticity must be shed at the leading edge to bring it back down. A critical LESP value of 0 corresponds to the kutta condition (basically the edge cannot hold any suction). I have not seen this implemented in 3D so I am not sure if there are any implementation subtleties.

I think its worth giving any of these methods a try! Let me know what you think.

[camUrban]

Thank you @antoniosgeme for your insights! Your comments led me to read some papers on modeling LEV separation using potential methods. For this conversation, let's refer to the rotated leading-edge suction technique as Method 1, the alpha-based leading-edge shedding technique as Method 2, and the LESP-based leading-edge shedding technique as Method 3.

Also, is it okay with you if I change the name of this issue to separation modeling?

[camUrban]

One thing I like about Method 1 is that reference 3 validated it with both experimental and CFD data. This validation is shown in figures 8, 10, and 11. Furthermore, figure 8 shows that Method 1's results are closer to experimental data than results from a standard UVLM. Other positives for Method 1 are its ease of implementation and its low computational overhead.

One downside is that I haven't found many other research groups besides the authors of reference 3 using it. Another con is the need for a heuristic stall angle-of-attack.

[camUrban]

Regarding Method 2, I was initially unconvinced of its accuracy due to some issues in reference 2 (which I believe implements this technique). Notably, figure 7 shows that a standard UVLM is more accurate than a UVLM with leading-edge shedding. Therefore, I reached out to the lead author, Dr. Roccia, with this question.

According to him, the downstroke, which is the first half of figure 7, is contaminated because the simulation is just starting. Once we ignore the first half, figure 7 shows that Method 2 is slightly more accurate than a standard UVLM.

In his response, he also said that if we want to model LEVs, we may benefit from using a hybrid UVLM-VPM (vortex particle method) to avoid numerical instability and speed up the solver. What are your thoughts on this idea? During a preliminary search, I found one example of someone else trying this.

He also mentioned the Polhamus technique (which I'm calling Method 1). He said that this technique is much more straightforward but that it only works for analyzing insect flight and falling Samara seeds. Do you have any ideas why Method 1 would be more applicable to these cases?

[camUrban]

Finally, I think I found a paper that implements Method 3 [5]! I haven't had time to read it yet, but I thought I'd leave the paper here in case you find it interesting.

[5] Hirato, Y., Shen, M., Gopalarathnam, A., & Edwards, J. R. (2021). Flow criticality governs leading-edge-vortex initiation on finite wings in unsteady flow. Journal of Fluid Mechanics, 910. https://doi.org/10.1017/JFM.2020.896

[antoniosgeme]

Separation modeling is a much better name, go for it.

There is definitely nothing to lose by trying out method 1. It should be straight forward to implement (famous last words). However, this is definitely something which would need validation in my opinion. Perhaps some CFD since experimental data sets are so hard to come by. I am not sure why Dr. Roccia suggested that it would only work for insects and samara seeds. Polhamus first discussed this theory with delta wings in mind. The LEV on delta wings is stabilized by the flux of vorticity along its core. Similar mechanisms exist in wings with high rotation rates so maybe this is what Dr. Roccia was referring to.

VPM is an interesting thought. I have never dealt with instabilities using panel methods. You probably have more experience on this given that you implemented the 3D UVLM solver. I can see the issue with rings being the TE-LE interaction. Particles are much further apart and can be regularized to not yield infinite velocities when the get close to one another. I read in the README that you use vortex decay? How does that work? Are you accounting for the spurious forces due to the vortex circulation change?

As for [5] I do not think they actually implemented it but just studied the validity of the method and it seems they found that the theory holds. I am sure they are currently working out the details of the implementation and a paper will be out soon.

[camUrban]

There is definitely nothing to lose by trying out method 1. It should be straight forward to implement (famous last words). However, this is definitely something which would need validation in my opinion. Perhaps some CFD since experimental data sets are so hard to come by.

Alright! Are you interested in implementing any of the methods and submitting a PR? Once we have it implemented we can find some previous CFD studies to compare against. No pressure if you're too busy, I just want to make sure we don't do double the work. 😊

I am not sure why Dr. Roccia suggested that it would only work for insects and samara seeds. Polhamus first discussed this theory with delta wings in mind. The LEV on delta wings is stabilized by the flux of vorticity along its core. Similar mechanisms exist in wings with high rotation rates so maybe this is what Dr. Roccia was referring to.

I wonder if Dr. Roccia may have meant that Method 1 only works for stabilized LEVs. There has been research demonstrating that some birds and bats also fly with stabilized LEVs, but it is quite recent [1]. In your opinion, would Method 1 also work for transient LEVs or just stabilized LEVs?

VPM is an interesting thought. I have never dealt with instabilities using panel methods. You probably have more experience on this given that you implemented the 3D UVLM solver. I can see the issue with rings being the TE-LE interaction. Particles are much further apart and can be regularized to not yield infinite velocities when the get close to one another. I read in the README that you use vortex decay? How does that work? Are you accounting for the spurious forces due to the vortex circulation change?

I think I may have used the term "vortex decay" incorrectly. I'm using a modified Biot-Savart law that eliminates the singularity issues at the wake vortex cores by giving them a finite radius. Then, I use the method described by Nguyen et al. on page 1711 to allow this finite radius to increase over time [3]. The increasingly thick radius prevents instabilities caused by wake vortices rolling up to be super close to each other and then shooting off to infinity. I'm not familiar with the spurious forces you mentioned. Could you elaborate on that?

[antoniosgeme]

Alright! Are you interested in implementing any of the methods and submitting a PR? Once we have it implemented we can find some previous CFD studies to compare against. No pressure if you're too busy, I just want to make sure we don't do double the work. 😊

I can probably take a crack at method 1 some time next week. Could you perhaps share some info on the implementation of the ring method? How are you enforcing the Kutta condition?

I wonder if Dr. Roccia may have meant that Method 1 only works for stabilized LEVs. There has been research demonstrating that some birds and bats also fly with stabilized LEVs, but it is quite recent [1]. In your opinion, would Method 1 also work for transient LEVs or just stabilized LEVs?

I do not think method 1 will work for LEVs that are in relative motion with the wing. Once the LEV detaches, stall conditions commence. However, insects and some birds flap in such a way so that this never happens, so perhaps for practical flows it could work well. It all depends on what we are trying to capture.

I think I may have used the term "vortex decay" incorrectly. I'm using a modified Biot-Savart law that eliminates the singularity issues at the wake vortex cores by giving them a finite radius. Then, I use the method described by Nguyen et al. on page 1711 to allow this finite radius to increase over time [3]. The increasingly thick radius prevents instabilities caused by wake vortices rolling up to be super close to each other and then shooting off to infinity. I'm not familiar with the spurious forces you mentioned. Could you elaborate on that?

Ah I see, so you regularized the vortices (great work!). Spurious forces arise when the circulation of a free vortex changes and the so-called the Brown-Michael Correction accounts for just that. Brown-Michael [6] argued that the change in circulation strength results in unbalanced forces in the flow, and this concept was later formalized by Michelin and Llewellyn Smith [7]. In order to maintain a force-free branch cut joining the vortex and its associated edge, a new term must be added to the convection velocity. In your case, as long as the circulation of the vortex remains fixed, even if its radius is changing, you should be on the clear.

[6] CE Brown and WH Michael. Effect of leading-edge separation on the lift of a delta wing. Journal of the Aeronautical Sciences (Institute of the Aeronautical Sciences), 21(10), 1954.

[7] Sebastien Michelin and Stefan G Llewellyn Smith. An unsteady point vortex method for coupled fluid{solid problems. Theoretical and Computational Fluid Dynamics, 23(2):127{153, 2009.

[camUrban]

Great! Thanks so much for offering to help out. 😄

In my opinion, the best description of the UVLM is in "Low-Speed Aerodynamics" by Joseph Katz and Allen Plotkin [7]. You should be able to find PDFs of it online (let me know if you're having trouble though). A broad overview is provided in chapters 12-14, and section 13.12 is especially relevant. Also, 13.11 discusses the Kutta condition. For a quicker, less detailed overview, the vortex lattice method Wikipedia page is great.

One area where I diverge from Katz and Plotkin is in the calculation of forces from the panel vortex strengths. Instead of using a modified version of Bernoulli's equation to calculate forces, I use the Kutta-Jokowski theorem to find them directly. Lambert, 2015 provides a good overview of the differences [8]. However, it may be necessary to implement Katz and Plotkins' technique for Method 1.

Feel free to keep me updated and continue to ask questions along the way!

[7] Katz, J., & Plotkin, A. (2001). Low-speed aerodynamics (2nd ed.). Cambridge University Press. https://doi.org/10.1017/cbo9780511810329
[8] Lambert T. (2015). Modeling of aerodynamic forces in flapping flight with the unsteady vortex lattice method. University of Liege.

[camUrban]

Hey @antoniosgeme! How is implementing Method 1 going? Let me know if you need any help or want to discuss any details regarding the code! 😄

[antoniosgeme]

Hey @camUrban I have not done any coding yet but I have looked through your code. As you said you are using vortex impulse to calculate the loads and not unsteady bernoulli. This is problematic as method 1 will require the load for each wing separately and vortex impulse can only calculate the total force experienced by the airplane as a whole. This will be the first step in the implementation.

Next we need to define a tangential direction. This is a bit awkward as wings can have 'washout' or twist but we'll figure something out.

Am I correct in assuming that the library only supports rigid wings?

[camUrban]

As you said you are using vortex impulse to calculate the loads and not unsteady bernoulli. This is problematic as method 1 will require the load for each wing separately and vortex impulse can only calculate the total force experienced by the airplane as a whole. This will be the first step in the implementation.

I think this first step may be easier than it seems! While it is true that I'm using the Kutta-Joukowski equation, I am doing it for each panel. Also, each wing object contains a list of all its panels. Therefore, we could iterate through the list of panels for each wing, solve the unsteady Bernoulli equation for each of its panels, and then sum the result. The Kutta-Joukowski force calculation is described in equation 2.12, and the unsteady Bernoulli calculation in equation 2.16 [8].

In essence, this will involve writing a method that's a Bernoulli alternative for calculate_near_field_forces_and_moments in unsteady_ring_vortex_lattice_method.py.

Next we need to define a tangential direction. This is a bit awkward as wings can have 'washout' or twist but we'll figure something out.

Hmm yes, this could be a bit tricky. Perhaps we could approximate it as the quarter-chord line of each panel?

Am I correct in assuming that the library only supports rigid wings?

Yes! At least for now, all surfaces are rigid and all motion is prescribed.

[antoniosgeme]

I see. You are calculating the force for each panel and then summing up. I wonder if this force would be the same locally on the wing using both methods. I think the answer is no, but the total force will be the same. We can easily test this once implemented.

I am looking at your force calculation, and for the unsteady term I see the difference in circulation but it is not divided by the timestep. Am I missing something?

unsteady_near_field_forces_geometry_axes = (
            self.current_operating_point.density
            * np.expand_dims(
                (self.current_vortex_strengths - self.last_panel_vortex_strengths),
                axis=1,
            )
            * np.expand_dims(self.panel_areas, axis=1)
            * self.panel_normal_directions
        )

[antoniosgeme]

Also I have noticed that during the first few timesteps of the simulation the airplane experiences large negative lift and drag coefficients. Below is a simulation of a wing with no camber, sweep or twist, very large aspect ratio at 20 degrees angle of attack. In general its behavior is correct! There is a delay in lift build up due to the wake. However I do no think that initial negative lift coefficient is physical. Any ideas as to why this happens?

[camUrban]

Hi @antoniosgeme! Sorry for the delay in getting back to you. These past few weeks have been crazy.

Regarding the unsteady force contribution, I think you just found a significant bug in the code. Thank you so much! Tomorrow, I'll publish a hot fix release that addresses this bug.

Out of curiosity, I corrected the equation and reran my test cases. The results still look reasonable. However, this doesn't fix the negative lift coefficients during the first few steps. I've noticed this result before and agree that it is strange. To be honest, I always brushed it off as an unimportant numerical issue because we are typically interested in results after the wake has fully developed. However, I'll reread how other UVLM implementations deal with the first timesteps and get back to you!

[antoniosgeme]

No problem, glad I could help.

I will look into the negative forces as well. I have some ideas.

[camUrban]

I just released a hotfix! The bug in the unsteady force calculation is fixed in v2.2.1.

Regarding the negative forces, I think the issue is pretty recent. I noticed that this image shows results from the unsteady static example without an initially negative lift coefficient. However, when I run the example now, I see the negative lift coefficient we discussed. I last updated this file on 4/11/2022, so the bug must be from the last three-ish months.

[antoniosgeme]

Hey @camUrban just wanted to let you know I haven't forgotten about this and will hopefully pick it up again next week.

One comment that I have regarding that initial spike is that in general we are expecting singular behavior at the first instant but in the positive lift direction. This is because of the infinite acceleration we are commanding (and therefore an infinite added-mass force). This added mass force comes in through the "unsteady term" that you have here, so we could start by looking at that calculation for a missing minus sign or something along those lines.

[camUrban]

Hey @antoniosgeme! No worries at all! I'm in the final two weeks of my current job, so I'm working overtime to wrap up several projects. Therefore, I think I'll be a bit MIA myself until Monday, 8/1. I hope that's okay, and again, thank you so much for all the effort on your part! Starting in August, I'll be excited to contribute more to tracking down this issue and implementing separation modeling.

[camUrban]

Hey @antoniosgeme, I hope you're doing well! I just finished up at my job, so I'll have more time to work on this in the next few weeks.

I'm trying to root out that negative-transient-lift bug you discovered. I think you're right that the issue comes from the unsteady_near_field_forces_geometry_axes calculation. This value is found via equation 2.11 [8]:

I'm approximating the vortex time derivative as a finite difference. I noticed that for the first time step, I assume that the bound vortices previously had zero vorticity. Therefore, the derivative is just the current vorticity divided by the time step.

Do you think that this is an appropriate assumption? Do you think it would be better to set this equation equal to zero for the first time step?

[antoniosgeme]

Hello @camUrban I believe it is okay to assume zero initial vorticity. This approximates the dirac delta response we expect in the first timestep (but for some reason in the wrong direction). A way to resolve this is to calculate the forces using bernoulli and compare. I gave it a try and was having trouble but then went on vacation and pretty much forgot everything 😅 . I will get back to you when I have something.

[antoniosgeme]

Force behavior should be an easy fix! #27