New Airfoil Model

[n72.75]

Based on the questions in this thread, and some discussions we have had in the NASSP team on implementing complex aerodynamic force and moment models I have put up a pull request that should address the challenges we've all faced at one time or another

Text from my PR:

This pull request is add a general-purpose "airfoil" as new, 4th type to Orbiter's existing AIRFOILSPEC type. This function CreateAirfoil4, like CreateAirfoil3 allows addon developers to specify their own functions for computing force and moment coefficients, however it expands greatly in the following areas:

Problem Statement:

CreateAirfoil3 allows the creation of multiple, static aerodynamic force-producing bodies in the local coordinates of the vessel to which they are bound. These Airfoils must be either "Vertical" like a wing, or "Horizontal" like a vertical stabilizer. In response to Mach and Reynolds number, velocity, angle of attack, and slip angle, Orbiter internally computes (by means a function passed from an addon by pointer) the CL , CD and Cm coefficients. which Orbiter internally uses to apply aerodynamic forces and moments to a vessel. This presents several major limitations:

• Force coefficients CL and CD produce forces relative to airspeed vector, not vessel body axes.
• Aerodynamic Forces and Moments can only be functional dependant on AoA and Slip, but cannot be a function of both at the same time.
• At very high slip angles AoA becomes poorly defined, and because lift is dependant on AoA, lift force becomes poorly defined. In a spacecraft this force is often the largest at high slip angles.
• Roll moments due to slip (dihedral effect, and other cross-coupling are not possible.

Proposed Solution is Provided in the following PR:

• Allow Orbiter to calculate axial, normal, and side force coefficients: CA , CN and, CY and roll, pitch and yaw moments Cl , Cm and, Cn from a user-defined function and allow these six coefficients to be used in calculating forces and moments on the vessel.
• Follow the paradigm of previous Orbiter airfoil functions, which will allow free combination of this "Airfoil" (Spacefoil?) type with all other types and control surfaces.

Significant Benefits:

• Allow for the implementation of analytical complex models of aerodynamic coefficient functions which may exist in historical documentation, without hacks.
• Facilitate easy implementation of lookup table aerodynamics.

Implementing something like this for all six coefficients would not currently be possible:


This update has not been thoroughly tested yet. I would also like some input from addon devs.

We can use this thread/github discussion/Discord to discuss. The burden of proving that this is a workable solution is squarely on me.

 

[n72.75]

Well, the force vectors point in the right direction:


Don't worry, I'll put it back like I found it

 

[GLS]

(finally got time to respond)
I'm approaching this from my (small) experience with the "Shuttle aero", so please excuse my ignorance on "generic aero".

Roll moments due to slip (dihedral effect, and other cross-coupling are not possible.

Hmm, I think its more like "it needs to be hacked":

double qbar = v->GetDynPressure();
double L = qbar * ORBITER_WING_AREA * ORBITER_SPAN * CLL;
v->AddForce( _V( 0.0, 0.5 * L, 0.0 ), _V( -1.0, 0.0, 0.0 ) );
v->AddForce( _V( 0.0, -0.5 * L, 0.0 ), _V( 1.0, 0.0, 0.0 ) );

Force coefficients CL and CD produce forces relative to airspeed vector, not vessel body axes.

Allow Orbiter to calculate axial, normal, and side force coefficients: CA , CN and, CY and roll, pitch and yaw moments Cl , Cm and, Cn from a user-defined function and allow these six coefficients to be used in calculating forces and moments on the vessel.

I don't know what is the standard (or if one even exists... this document has both types), but maybe allowing the user to choose between stability and body axis would be more flexible.

Will this make capsules or chutes simpler?

 

[n72.75]

(finally got time to respond)
I'm approaching this from my (small) experience with the "Shuttle aero", so please excuse my ignorance on "generic aero".

Hi GLS,

I very much appreciate some feedback on this. I think in trying to make this as generalized as possible, my goal is to support the more specialized cases of aerodynamics, like shuttles, capsules, &c..

Hmm, I think its more like "it needs to be hacked":

double qbar = v->GetDynPressure();
double L = qbar * ORBITER_WING_AREA * ORBITER_SPAN * CLL;
v->AddForce( _V( 0.0, 0.5 * L, 0.0 ), _V( -1.0, 0.0, 0.0 ) );
v->AddForce( _V( 0.0, -0.5 * L, 0.0 ), _V( 1.0, 0.0, 0.0 ) );

The downside to doing anything hacky like this with AddForce(), is that the force must be applied for a full timestep. Allowing the code to run inside Vessel::UpdateAerodynamicForces means that it can run under Orbiter's integrators and substeps. But you are correct that it is not technically impossible. NASSP has some hacks for combined yaw and pitch moments at high combinations of alpha and beta.

I don't know what is the standard (or if one even exists... this document has both types), but maybe allowing the user to choose between stability and body axis would be more flexible.

I think the best solution for calculations is to have Orbiter working on body coefficients internally, since the forces and moments will have to be based on body, not wind/stability coordinates when Orbiter actually applies them to the vessel. It is definitely possible to convert between these two systems though, so if you only had stability coefficients, they could be converted to body; this could be done either offline to make a lookup table, or in realtime from the existing data.






Something that might make a bit more sense than what I'm doing now:

        } else if (af->align == FORCE_AND_MOMENT) {
            ((AirfoilCoeffFuncEx2)af->cf)((VESSEL*)modIntf.v, aoa, beta, gamma, sp.atmM, Re0 * af->c, af->context, &CA, &CN, &CY, &Cl, &Cm, &Cn);
            if (af->S) S = af->S;
            else       S = fabs(ddir.z) * cs.z + fabs(ddir.x) * cs.z; // use projected vessel CS as reference area

            Vector AeroForce((CY * S) * sp.dynp, (CN * S) * sp.dynp, -(CA * S) * sp.dynp);

            AddForce(AeroForce, af->ref);
            Lift = dotp(AeroForce, ldir);
            Drag = Lift = dotp(AeroForce, ddir);

            Amom_add.x += Cm * sp.dynp * af->S * af->c;
            Amom_add.y -= Cn * sp.dynp * af->S * af->c;
            Amom_add.z -= Cl * sp.dynp * af->S * af->c;

might be to change AirfoilCoeffFuncEx2 from:

typedef void (*AirfoilCoeffFuncEx2)(
    VESSEL* v,
    double alpha, double beta, double gamma,
    double M, double Re, void* context,
    double* CA, double* CN, double* CY,
    double* Cl, double* Cm, double* Cn);

to:

typedef void (*AirfoilCoeffFuncEx2)(
    VESSEL* v,
    Vector WindDir,
    double M, double Re, void* context,
    double* CA, double* CN, double* CY,
    double* Cl, double* Cm, double* Cn);

There are still singularities with the C[A,N,Y,l,m,n](aoa,roll) model I'm using for plotting this data:
It just happens to be easier to deal with than (aoa,beta).


So maybe it's better to make that a vector function and/or let addon devs decide how to calculate coefficients from that. In the case above, I am literally turning those coefficients into vector functions (e.g CN(X,Y,Z)) in order to plot them, using Matlab/Octave's sph2cart function. (technically I'm plotting magnitude, because otherwise the plots are weird and the negative values are self-intersecting with the positive ones).

See attached code:

Spoiler: Matlab Code

clear;
clc;

phi = [0    20    40    60    80    90    100    120    140    160    180    200    220    240    260    270    280    300    320    340    360];
alpha = [0 20 40 60 80 90 100 120 140 160 180];

phiii = linspace(0,360,100);
alphaii = linspace(0,180,50);

[PHI ALPHA] = meshgrid(phi,alpha);
[PHIi ALPHAi] = meshgrid(phiii,alphaii);

CA = [2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401    2.401
3.166    3.129    2.964    2.787    2.6    2.707    2.517    2.952    3.34    3.63    3.86    3.554    3.212    2.798    2.37    2.625    2.576    2.842    2.945    3.099    3.166
4.124    3.919    3.33    3.036    2.669    2.687    2.71    3.522    4.202    4.713    4.873    4.672    3.909    3.063    2.48    2.628    2.706    3.286    3.699    4.145    4.124
3.366    3.251    2.989    2.56    2.11    1.988    2.154    2.854    3.423    3.857    4.017    3.84    3.187    2.473    1.98    1.988    2.164    2.828    3.187    3.407    3.366
1.433    1.347    1.22    0.997    0.821    0.741    0.832    1.122    1.315    1.471    1.518    1.436    1.224    0.976    0.786    0.763    0.847    1.125    1.29    1.371    1.433
0.001    -0.004    -0.001    0.003    0.004    0.004    0.004    0.004    0.004    0.005    0.005    0.005    0.004    0.005    0.005    0.005    0.005    0.005    -0.001    -0.003    0.001
-1.529    -1.445    -1.226    -0.993    -0.801    -0.73    -0.827    -1.128    -1.317    -1.399    -1.446    -1.359    -1.225    -0.974    -0.446    -0.754    -0.839    -1.13    -1.323    -1.492    -1.529
-4.142    -3.905    -3.332    -2.725    -2.111    -1.988    -2.186    -2.943    -3.298    -3.566    -3.515    -3.321    -2.767    -2.419    -1.999    -2.049    -2.204    -2.987    -3.535    -3.971    -4.142
-5.136    -4.899    -4.366    -3.591    -2.76    -2.742    -2.81    -3.513    -3.981    -4.308    -4.308    -3.882    -3.255    -2.962    -2.595    -2.829    -2.868    -3.769    -4.533    -4.956    -5.136
-4.106    -3.966    -3.491    -3.201    -2.7    -2.813    -2.63    -3.044    -3.209    -3.313    -3.344    -3.246    -3.029    -2.901    -2.556    -2.858    -2.76    -3.299    -3.646    -4.007    -4.106
-2.56    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59    -2.59];

CN = [0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002    0.002
1.314    1.225    0.96    0.609    0.227    -0.001    -0.224    -0.663    -1.096    -1.434    -1.618    -1.404    -1.058    -0.632    -0.212    0.006    0.235    0.624    0.951    1.208    1.314
3.77    3.372    2.366    1.419    0.468    -0.001    -0.472    -1.655    -2.953    -4.022    -4.416    -3.989    -2.759    -1.456    -0.44    0.01    0.493    1.548    2.625    3.566    3.77
6.235    5.657    4.278    2.415    0.733    -0.001    -0.751    -2.691    -4.853    -3.358    -7.368    -6.629    -4.532    -2.357    -0.702    0.009    0.78    2.678    4.558    5.923    6.235
8.578    7.626    5.647    3.026    0.906    0    -0.933    -3.422    -6.029    -8.215    -8.985    -8.015    -5.615    -2.993    -0.893    0.004    0.976    3.435    5.971    7.762    8.578
9.33    8.237    5.874    3.115    0.922    0    -0.961    -3.587    -6.32    -8.339    -9.236    -8.203    -5.902    -3.149    -0.935    0.001    1.09    3.584    6.291    8.46    9.33
9.153    8.132    5.659    3.04    0.99    0    -0.938    -3.465    -6.069    -7.853    -8.64    -7.625    -5.655    -3.02    -0.9    0.003    0.981    3.479    6.114    8.393    9.153
7.605    6.757    4.746    2.578    0.734    0    -0.768    -2.778    -4.661    -3.146    -6.438    -5.715    -3.928    -2.321    -0.717    0.007    0.8    2.838    5.037    6.867    7.605
4.644    4.181    3.073    1.695    0.491    0    -0.5    -1.643    -2.777    -3.668    -3.886    -3.295    -2.281    -1.406    -0.468    0.006    0.526    1.785    3.189    4.224    4.644
1.703    1.552    1.143    0.716    0.248    -0.001    -0.246    -0.678    -1.041    -1.298    -3.9    -1.271    -0.982    -0.659    -0.239    0.002    0.257    0.739    1.189    1.566    1.703
0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0];

CY = [ -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002    -0.002
0    -0.408    -0.724    -0.929    -0.994    -1.048    -0.966    -0.985    -0.821    -0.473    -0.001    0.46    0.779    0.92    0.898    1.005    0.977    0.943    0.72    0.405    0
0    -1.143    -1.84    -2.28    -2.295    -2.346    -2.33    -2.637    -2.322    -1.384    -0.001    1.368    2.137    2.264    2.109    2.269    2.307    2.457    2.033    1.209    0
0    -1.947    -3.375    -3.931    -3.711    -3.551    -3.778    -4.367    -3.876    -2.32    -0.001    2.307    3.585    3.748    3.445    3.516    3.775    4.321    3.589    2.037    0
0    -2.64    -4.49    -4.959    -4.615    -4.28    -4.721    -5.576    -4.84    -2.876    -0.001    2.807    4.482    4.809    4.423    4.36    4.776    5.568    4.74    2.686    0
0    -2.857    -4.67    -5.103    -4.715    -4.348    -4.857    -5.826    -5.072    -2.923    -0.001    2.875    4.714    5.065    4.626    4.515    4.945    5.806    4.99    2.933    0
0    -2.821    -4.495    -4.973    -4.6    -4.269    -4.734    -5.637    -4.871    -2.754    -0.001    2.674    4.51    4.832    4.431    4.359    4.78    5.628    4.846    2.911    0
0    -2.335    -3.752    -4.169    -3.703    -3.545    -3.825    -4.494    -3.731    -2.148    -0.001    1.999    3.114    3.659    3.464    3.615    3.834    4.548    3.973    2.374    0
0    -1.426    -2.388    -2.663    -2.348    -2.372    -2.392    -2.613    -2.194    -1.264    -0.001    1.136    1.778    2.179    2.186    2.424    2.421    2.784    2.476    1.422    0
0    -0.505    -0.84    -1.043    -1.011    -1.066    -0.988    -0.998    -0.783    -0.429    -0.001    0.418    0.733    0.94    0.951    1.074    1.027    1.07    0.874    0.51    0
-0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001    -0.001];

CAii = interp2(PHI,ALPHA,CA,PHIi,ALPHAi,'cubic');
CNii = interp2(PHI,ALPHA,CN,PHIi,ALPHAi,'cubic');
CYii = interp2(PHI,ALPHA,CY,PHIi,ALPHAi,'cubic');


## Plot CA
#*******************************************************************************
subplot(2,3,1)
colormap("Jet")
contourf(PHIi,ALPHAi,CAii);
colorbar;
xlabel("Roll");
ylabel("AoA");
title("Skylab Force Coefficient C_A")


subplot(2,3,4)
[XA,YA,ZA] = sph2cart(deg2rad(PHIi),deg2rad((ALPHAi)-90),abs(CAii)); ##plots magnitude
surf(ZA,YA,XA,sqrt(ZA.^2+YA.^2+XA.^2))
shading flat;
axis equal
xlabel("X");
ylabel("Y");
zlabel("Z");
title("Skylab Force Coefficient C_A")
colormap("Jet")
colorbar;

## Plot CN
#*******************************************************************************
subplot(2,3,2)
colormap("Jet")
contourf(PHIi,ALPHAi,CNii);
colorbar;
xlabel("Roll");
ylabel("AoA");
title("Skylab Force Coefficient C_N")


subplot(2,3,5)
[XN,YN,ZN] = sph2cart(deg2rad(PHIi),deg2rad((ALPHAi)-90),abs(CNii)); ##plots magnitude
surf(ZN,YN,XN,sqrt(ZN.^2+YN.^2+XN.^2))
shading flat;
axis equal
xlabel("X");
ylabel("Y");
zlabel("Z");
title("Skylab Force Coefficient C_N")
colormap("Jet")
colorbar;

## Plot CY
#*******************************************************************************
subplot(2,3,3)
colormap("Jet")
contourf(PHIi,ALPHAi,CYii);
colorbar;
xlabel("Roll");
ylabel("AoA");
title("Skylab Force Coefficient C_Y")


subplot(2,3,6)
[XY,YY,ZY] = sph2cart(deg2rad(PHIi),deg2rad((ALPHAi)-90),abs(CYii)); ##plots magnitude
surf(ZY,YY,XY,sqrt(ZY.^2+YY.^2+XY.^2))
shading flat;
axis equal
xlabel("X");
ylabel("Y");
zlabel("Z");
title("Skylab Force Coefficient C_Y")
colormap("Jet")
colorbar;

Will this make capsules or chutes simpler?

Yes. That was in fact this project's raison d'être.

 

[GLS]

I think the best solution for calculations is to have Orbiter working on body coefficients internally, since the forces and moments will have to be based on body, not wind/stability coordinates when Orbiter actually applies them to the vessel. It is definitely possible to convert between these two systems though, so if you only had stability coefficients, they could be converted to body; this could be done either offline to make a lookup table, or in realtime from the existing data.

Yes, body makes sense internally, but I'm not sure it is the best idea to "force" that on the devs. Having the 2 options available, with the conversion done inside Orbiter, would IMO make it more "plug-and-play" for the dev.

Vector WindDir

Shouldn't just alpha and beta be enough? From what I'm thinking, that should define "what part of the vehicle is towards the air stream", thus the aero should be determinable. One can then rotate the vessel along the air stream direction, but it is still going to have the same part towards the air, so same aero results... yes, globally the resulting forces might now point down instead of up, but it should be the same CN, CA, etc.
Or am I being too simplistic, and this actually is the source of the problems for capsules?

 

[Urwumpe]

Could you maybe include asymmetric wind effects or ground effect into the core? Its pretty hard to even get this halfway right and especially large delta wing aircraft like we often have in Orbiter, are very suspectible to it.

 

[GLS]

Could you maybe include asymmetric wind effects or ground effect into the core? Its pretty hard to even get this halfway right and especially large delta wing aircraft like we often have in Orbiter, are very suspectible to it.

The Shuttle aero data already has that: increments to CL, CD and CM are listed as function of height-over-span, alpha and aerosurface positions.

 

[n72.75]

Yes, body makes sense internally, but I'm not sure it is the best idea to "force" that on the devs. Having the 2 options available, with the conversion done inside Orbiter, would IMO make it more "plug-and-play" for the dev.



Shouldn't just alpha and beta be enough? From what I'm thinking, that should define "what part of the vehicle is towards the air stream", thus the aero should be determinable. One can then rotate the vessel along the air stream direction, but it is still going to have the same part towards the air, so same aero results... yes, globally the resulting forces might now point down instead of up, but it should be the same CN, CA, etc.
Or am I being too simplistic, and this actually is the source of the problems for capsules?

I can look into adding those conversion formulæ into the code, so that it's simpler for addon devs. I think the hardest part there is naming things correctly without resorting to something horrible like unions. I'll give that some thoughts.

Re: alpha and beta.
I think it's important to provide as much flexibility here as possible. Any system that uses two coordinates to define a point in 3d space will have singularities. When beta is 90 and -90, alpha is undefined. While these special cases are easy to deal with, you don't have equal resolution in your coefficients at all relative wind angles. This is particularly important for spacecraft which can assume any orientation, while aircraft won't usually fly at super high slip angles.

Could you maybe include asymmetric wind effects or ground effect into the core? Its pretty hard to even get this halfway right and especially large delta wing aircraft like we often have in Orbiter, are very suspectible to it.

As far as I know, the wind effect works by just adding some relative air velocity along the right axes, so if you had an asymmetric force model (which this feature will now easily allow you to impliment), you wouldn't need any special handling.

Ground effect can be simulated with the current airfoil model, the same way that mach number or Knudson number can be used as a parameter in a custom aerodynamic model, so can altitude be used to select different data in a lookup table.

I know Orbiter has an "induced drag" function, but I've never looked at it, if it doesn't take altitude and wing span into effect we can definitely add it. But that's a little out of scope from what I'm trying to do here, so probably a separate project.