Implementing the gust encounter for 3D meshes and for free-flying, elastic aircraft

[ArneVoss]

Proposed Changes

I would like to be able to calculate a gust encounter of a free-flying, elastic aircraft. The rigid body reaction (in terms of velocities u, v, w and rotation rates p, q, r) as well as the elastic structural dynamic reaction (displacement of all surface nodes of the CFD mesh) are calculated by our own aeroelastic solver, coupled via the python interface, while SU2 shall provide the aerodynamic side, as a higher fidelity alternative to panel methods such as VLM & DLM.

Update, 13.04.2023: The combination of grid velocities from the moving + deforming grid turned out to be more complicated as expected (at least for me) and the moving grid is not fully implemented for unsteady simulations (according to Pedro, some source terms are missing). The new approach uses a farfield onflow condition at alpha=0.0 and translates and rotates the whole aircraft in the elastic mesh for the rigid body + structural dynamic reaction. This is more laborious on my side (transformation of the surface, the deformation vector and the forces about the Euler angles) but simplifies the work on the SU2 side:

• extend the gust to 3D meshes (current implementation only for 2D meshes)
• update the grid movement during the time domain simulation via the python interface
  Involved parameters:
  GRID_MOVEMENT= ROTATING_FRAME
  TRANSLATION_RATE= u v w
  ROTATION_RATE= p q r
• update the surface mesh deformation during the time domain simulation via the python interface
  Involved parameters:
  DEFORM_MESH= YES
  MARKER_DEFORM_MESH= ( list of markers )

Summing up, I will need to combine the grid velocities due to aircraft movement, elasticity and gust. Currently, this is not possible in SU2.
I split the work into the following steps. Here is what is working / not working:

• gust with 3D meshes --> plausibility check with rectangular, 3 meter wing at 2° + long gust equivalent to 2° = approx doubling of the lift coefficient
• gust + moving grid --> plausibility check by comparison with farfield onflow condition from above
• new python interface for grid movement update --> plausibility check with a) sudden jump in heave velocity w, b) sinusoidal velocity w at different frequencies to observe influence of unsteady aerodynamics (phase shift and reduction of resulting lift)
• mesh deformation via python interface --> plausibility check with sinusoidal displacement z, scaled and shifted by 90° to correspond to velocity w from above, leads to same results
• mesh deformation + moving grid --> differences with respect to farfield onflow condition / does NOT lead to same results as above. The grid velocities look right though, both on the temporary screen output and when plotting them by slicing through the volume results, but the amplitude of the lift coefficient is too small. I guess that the velocities are not set properly / mixed up / require a call to some update function or similar. I'm stuck at this point. Does someone have any idea what I'm doing wrong?? I feel that I am missing an important step here...
• final test: mesh deformation + moving grid + gust
• final test: mesh deformation + gust

Apart from the steps above, the following questions / tasks remain:

• The the Split Velocity Method (classes SetWind_GustField and CSourceWindGust) is not used in SU2 as far a I understand and for all gust types, the inputs are zero. Instead, the Field Velocity Method is used, but when trying to remove the functions SetWindGust and SetWindGustDer, I broke the code elsewhere. In order to enable 3D meshes, I simply set the outputs to zero for now. I believe that we should clean up the code, but my C++ skills are the limiting factor here, so I need some help with that. --> In a new, dedicated pull request

Related Work

The gust encounter for 3D meshes was requested in #1319 and in the Forum https://www.cfd-online.com/Forums/su2/234021-gusts-3d.html and https://www.cfd-online.com/Forums/su2/225914-no-slip-condition-3d-gust.html

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

• I am submitting my contribution to the develop branch.
• My contribution generates no new compiler warnings (try with --warnlevel=3 when using meson).
• My contribution is commented and consistent with SU2 style (https://su2code.github.io/docs_v7/Style-Guide/).
• I have added a test case that demonstrates my contribution, if necessary.
• I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp), if necessary.

[pcarruscag]

How are you combining moving grid and deformation?

[pcarruscag]

Is it possible the lift coefficient is different because the reference velocity is not the same in the 2 configurarions? Can you output the forces breakdown and see what is the force factor?

[ArneVoss]

You mean that only the coefficients are calculated wrong while the solution itself is correct?

Actually, I checked that yesterday by calculating the lift coefficient manually from the nodal forces, which I obtained through the python interface with GetFlowLoad. I found that the CL from SU2 is correct, the results were numerically equivalent. It might be that the Vertex Tractions are calculated with a wrong reference value?? That would be something that happens deep in the solver, which is beyond my abilities...

For the farfield, I set MACH_NUMBER= 0.0 and used MACH_MOTION= 0.2, which should be used for the scaling of coefficients, forces, etc. as far I understood. The velocity are u=68.0588 m/s, v=0.0, w=+/- 2.37 m/s, which corresponds exactly to Ma=0.2 at sea level with pressure=101325.0, density=1.225 and temperature=288.15. I checked the values by slicing through the volume mesh, and I can see exactly these numbers in far in front of the wing / close to the farfield.

[pcarruscag]

I think that is what is happening, the logic to use either MACH_MOTION or the farfield velocity as the reference is in CEulerSolver::SetReferenceValues, I think if you have fluid load markers we will assume FSI and not use MACH_MOTION.

Give this a try and if it fixes the issue we need more robust logic.

[ArneVoss]

I tried changing line 1387 in CEulerSolver::SetReferenceValues to simply if (dynamic_grid) {...}, but that had no effect.

[ArneVoss]

In the meantime, I did a few additional comparisons. In the following plot, you can see a gust encounter of my 3m wing, in green color farfield onflow + gust and in violet color the same but with moving grid + gust. As we can see, there is a significant difference, which shouldn't be there. What is happening here??
My current understanding is that the green curve (farfield onflow + gust) is correct, because it looks similar to other plots I found in literature and from asking my colleagues. I have the feeling that in my implementation, I am violating some physics / broke the numerics which describe the physics, but I have not idea what I'm doing wrong...

[pcarruscag]

I don't have an answer here unfortunately, have you ruled out issues or differences due to the GCL correction we make when computing the dual time source term?

[ArneVoss]

Do you have any idea where / in which file I should look?

[pcarruscag]

In CFVMFlowSolverBase<V, FlowRegime>::SetResidual_DualTime (CFVMFLowSolverBase.inl::~1600) there is different treatment for dynamic grid.

And for rotating frame (which you have in some cases) we call CSolver::SetRotatingFrame_GCL from CEulerSolver::Source_Residual. Since you are combining things I don't know if this makes some weird interaction.

[ArneVoss]

Great, thank you! I will look into this within the next days.

[ArneVoss]

Hello Pedro,
thank you for your last comment, I think that some unintended double accounting was happening.
I removed two sections:

• the first with if (rotating_frame) is captured by if (dynamic_grid) in CFVMFLowSolverBase
• the second if (windgust) should be zero anyway

The results look plausible to me: the green line shows farfield onflow + gust while the violet line shows moving grid + gust. If you zoom in (second plot) you see some small differences, for which I have no explanation so far.

What do you think?
(Side note: I will be in the Alps for the next few day and be back in office on Thursday. )

[ArneVoss]

Hi Pedro,
The following two configs are for a sinus motion using an elastic mesh in combination with farfield or rotating frame.
Previous results: Looking at the CL, there used to be a difference of about 10%.
Current Results: Good agreement with only a small difference.

Side note: Although the Naca0009 is a symmetric airfoil, my mesh appears to be not perfectly symmetric, resulting in non-zero lift at alpha=0.0 deg. I didn't look into this any further as the mesh is only a quick test case and symmetry probably not important for our problem.

NACA0009_elastic.zip
NACA0009.zip

[ArneVoss]

... and here are the two gust cases.

NACA0009_gust.zip

[pcarruscag]

This is a somewhat tricky problem.
Moving reference frames are usually used in steady state, to do unsteady in a moving frame with variable velocity there are source terms due to the acceleration of the frame (we do not have those terms).
Then, the rotating frame source term is due to steady rotation. So based on your configs it doesn't have influence. As a side note, to combine rotation and translation the source term needs to be modified (I did not know this at the time you added the feature). See slide 11 https://courses.ansys.com/wp-content/uploads/2020/09/Flows-with-Moving-and-Rotating-Objects-Lesson-4-Handouts-1.pdf
On the GCL, the GCL term we have in the source term function, and the one in the dual time function are slightly different. The second one is normal for deforming meshes (see Tom's thesis). The first one (in sources) is to make steady-state rotating frame preserve free-stream better, and so it should probably not be called for unsteady problems.

Then the windgust derivatives do not look very healthy to be honest, commented out code is never a good sign. Have you looked at the paper mentioned in the code?

[ArneVoss]

Thank you Pedro for the information, I will definitely have a look at the documents, although I admit that it sounds pretty complicated to me...

Then the windgust derivatives do not look very healthy to be honest, commented out code is never a good sign. Have you looked at the paper mentioned in the code?

In the paper it says that the split velocity method (which uses the derivatives) didn't work out. So it is kind of a left-over in the code which I believe should clean up.

[ArneVoss]

I forgot to mention: I found that the small offset (delta Cl = 0.0001) between farfield onflow and grid velocities observed in the previous plots already exist in a normal time simulation with alpha=0.0 deg. So the offset has nothing to do with a) the mesh deformation or b) the gust.

[pcarruscag]

Maybe we can approach it from the minimal set of features you need.
Mesh deformation seems to work fine, then you are applying the gust via the farfield?
Then for the free-flying part maybe the best would be a translating moving frame (you'd need to implement the acceleration term). (Translating the mesh over large distances would lead to floating point issues)
For rotation maybe rigid rotation would be the easiest route? This would not require the rotating frame source term.

[ArneVoss]

Hi Pedro,
I'm slightly confused right now and not sure if I'm able to do the necessary modifications. I understand that it is unlikely someone else might implement this but I think that this is beyond my abilities right now. If you don't mind, do you have some time to talk through my options in a video call? I hope that will help me to get a clearer picture. Next week, I'm generally available except for Monday, which is a public holiday here. Germany is 9 hours in advance of California, so your morning / my evening might work?

[pcarruscag]

We have our developers meeting in 1h 30m, happy to chat about this https://meet.jit.si/SU2_DevMeeting

[ArneVoss]

Hi Pedro,
As discussed last week, I now translate and rotate the whole aircraft in the elastic mesh in combination with a farfield onflow. Implementing and doing the coordinate transformations right took me a few hours, but now everything seems to work properly and fast :)

1. Currently, activating the gust resets/overwrites the grid velocities due to the deformed mesh, but I haven't found the place yet. Any ideas?

2. Should I clean up / remove the split velocity approach as described in the first post or would you like to keep it?

3. How to handle the new approach, should I close this pull request and open a new one? There are a few commits which I needed to undo.

4. Generally, I still need the rotating frame approach for steady maneuver load cases, e.g. to calculate the pitching, rolling or yawing aircraft in a steady simulation. The acceleration terms are zero in this case, but I understood that the Coriolis-Term with omega x velocity is missing. Is that correct? I guess they are probably important for objects like a propeller which spins at a couple of thousand RPMs but maybe it is justified to neglect them for slow objects?

[pcarruscag]

Hi Arne,
I'm sorry for the complicated meeting schedules.
SetWind_GustField seems to be resetting the grid velocity based on some conditions.
We should split the cleanup of obsolete code into another PR if possible.
I think you can continue on this PR.
The rotating frame term we have assumes 0 translation velocity, we have w x u and based on that Ansys presentation it should be w x (u - u_translation) if the translation is parallel to the rotation axis there is no issue, if not the steady state assumption is questionable, so it may be fine to use it as-is.

[ArneVoss]

SetWind_GustField seems to be resetting the grid velocity based on some conditions.

OK, found it :)

[ArneVoss]

Where do the meshes come from or where should I place the new mesh and the restart files? For example, there is no ./TestCases/gust/mesh_NACA0012_inv.su2 in the repository...

[pcarruscag]

Is it possible to use the Onera M6 mesh that we already have in the repo?
Meshes and restart files are kept in this repo https://github.com/su2code/TestCases/pulls

[ArneVoss]

I see that the Onera M6 mesh is a half-wing with a symmetry plane and I'm not sure how the mesh deformation behaves in this case (for example, there should be zero out-of-plane movement).

Personally, I like really simple examples where I know what should happen (symmetric airfoil, subsonic flow, no sweep) and I did all my testing during the last weeks with the NACA0012 3m wing. At the same time, it is closer to how I would model a "real" aircraft (aircraft in the center of a spherical farfield) than the Onera M6, which is more wind-tunnel-like. The mesh and everything is there, so it's no additional work for me and I would be happy to contribute this as a new example.

What's the main argument in favor of the Onear M6 wing / objection against a new example?

[pcarruscag]

Just keeping the size of the repo on the small side. But go for it if it is better for the test, please use a binary restart instead of csv.

[ArneVoss]

please use a binary restart instead of csv

How do I do this? To my understanding, SU2 currently needs both files.

[ArneVoss]

I haven't fully figured out how the two repositories interact, but if I copy the data from su2code/TestCases#124 into ./TestCases/gust and launch the hybrid_regression.py, the new testcase succeeds.

[pcarruscag]

Ah yes, that problem... I thought the csv was only needed for something else

[pcarruscag]

You'll have to take out the tecplot output from the config https://github.com/su2code/SU2/actions/runs/5006027201/jobs/8971286794

[pcarruscag]

Lgtm, now that there are some tests maybe you can remove the seemingly dead code in another PR

[ArneVoss]

Thank you Pedro, the code removal is on my To-Do list 👍

[pcarruscag]

Sounds good, go ahead and merge this when you can.

[ArneVoss]

I think I don't have the authority to merge, it says: "This pull request can be automatically merged by project collaborators"
and "Only those with write access to this repository can merge pull requests."

[pcarruscag]

Ah right, I thought you were a member already, I've sent an invite.