[BUG] Negative Spike in Lift Forces

[MOLOCH-dev]

Bug Description

Thank you for this great software. @RaphZufferey and I am using it currently to analyse flapping flight with UVLM to incorporate Aeroelasticity effects. I observed an initial Negative spike in Lift Force after plotting results from Ptera. This spike makes no physical sense for a flapping wing, so I was wondering if this is a numerical issue that could potentially be solved.

Reproduction

Running the provided example
with the addition of

ps.output.plot_results_versus_time(
    show=True,
    save=True,
)

will reproduce the above results

Expected Behavior

No spike in the lift coefficient is observed in the README.

Desktop

Please complete the following information.

• OS: 'Windows 10 Home'
• Version: '22H2'
• PteraSoftware Version : I have forked the "develop" branch, the closest release would be 3.0.1

Additional Context

A similar issue was being discussed about a year ago here #17 . I am wondering if any resolution was achieved to the same or any work is in progress. I am open to picking up some work towards this issue as I have been making myself familiar with Ptera's UVLM solver.

[camUrban]

@MOLOCH-dev Thanks for the bug report. As you pointed out, the negative lift spike is most likely a numerical issue. The plot from the README was from an unsteady simulation of a static (non-flapping) wing. This may be the reason we don't see such a negative spike. The last time I looked into this, I hit a dead-end and, unfortunately, became sidetracked with other bugs.

However, I have a hunch about a potential culprit. During the first time step, there are a few things different from subsequent time steps:

1. There are no wake vortices. Normally, the rear line vortices of the ring vortices on a wing's trailing edge (red) are mostly canceled out by the front line vortices of the first row of wake ring vortices (blue).
   Katz, J., & Plotkin, A. (2001). Low-speed aerodynamics (2nd ed.). Cambridge University Press. https://doi.org/10.1017/cbo9780511810329. Perhaps we should exclude the red line vortices on the initial time step?
2. There is no motion-velocity. During the first time step, a function calculates the velocity at each collocation point. This velocity is the sum of the freestream velocity and the velocity induced by all the vortices. However, during future time steps, the velocity calculations also include an approximate velocity due to each panel's movement, calculated by subtracting the current position from the last and dividing by the change in time. This method won't work for the first time step, but we could calculate it directly if this is causing the problem.

[camUrban]

Please let me know if you have any luck finding the bug or need help. While I'm a bit overwhelmed with work this week, I will take a look and report back what I find!

[camUrban]

@MOLOCH-dev Okay, I have some updates on this issue. I've ruled out the motion-induced velocity calculation as being the sole reason for this bug. This is because the negative lift spike still occurs in simulations with no motion.

[camUrban]

Also, I found a few reference plots to compare our results against:

Katz, J., & Plotkin, A. (2001). Low-speed aerodynamics (2nd ed.). Cambridge University Press. https://doi.org/10.1017/cbo9780511810329

This case is something easily simulated in Ptera Software. The one catch is that Katz and Plotkin use a different method for calculating forces than I'm using. They use a modified version of Bernoulli's equation, whereas I use the unsteady Kutta-Joukowski equation. I've read a paper comparing the two approaches (referred to as the Katz method and the Joukowski method), and they generate extremely similar results for most converged simulations, but the transient lift might be different.

However, I also found the following plot from a paper that uses the Joukowski method:

Gao, W., Liu, Y., Li, Q., & Lu, B. (2023). Aerodynamic Modeling and Simulation of Multi-Lifting Surfaces Based on the Unsteady Vortex Lattice Method. Aerospace 2023, Vol. 10, Page 203, 10(2), 203. https://doi.org/10.3390/AEROSPACE10020203

As you can see, both examples show lift coefficients that begin at roughly positive infinity and then dip slightly lower than expected before finally asymptotically converging. This is the opposite of what we're seeing in our simulations.

[camUrban]

@MOLOCH-dev I just created a branch ("hotfix-3.0.1") for us to work on this issue.

[camUrban]

@MOLOCH-dev As a baseline, I created suddenly_accelerated_plate.py in a new directory, "issue testing." It runs a simulation identical to one of the simulations shown in the Katz and Plotkin figure. Here's what it looks like now:

[antoniosgeme]

I believe the "unsteady_near_field_forces_geometry_axes" term in the UnsteadyRingVortexLatticeMethodSolver.calculate_near_field_forces_and_moments() method has the wrong sign. It should be negative (or the order of subtraction of the circulations needs to be changed). After correcting this, the force behavior is correct as shown in the plot below that attempts to replicate the above plots. The large positive peak at the initial instant is actually the correct behavior due to the large acceleration. As the aspect ratio increases the solution converges to Wagner's analytical solution. The fact that early on the agreement with Wagner is not perfect is also expected. If anyone is interested, here is recent paper that explains why this is the case.

[antoniosgeme]

If the fix isn't obvious, checkout my hotfix branch.