
Orbiter SDK Orbiter software developers post your questions and answers about the SDK, the API interface, LUA, meshing, etc. 

Thread Tools 
03292019, 09:46 AM  #1 
Orbinaut

Trouble with the produced drag of airfoils
Hello everyone,
I'm currently (still) working on a lifting capsule Soyuz model for reentry simulation. As described in the API Guide I have defined two airfoils, one for horizontal and one for vertical lift with all the given parameters. The Code is as follows: Code:
double D1 = 2.2; VECTOR3 crossSecVec; crossSecVec.x = 4.052293854; // projection along xaxis into yzplane crossSecVec.y = 4.052293854; // projection along yaxis into xzplane crossSecVec.z = PI * ((D1 * D1) / 4); SetCrossSections(crossSecVec); ClearAirfoilDefinitions(); hVertAirfoil = CreateAirfoil2( LIFT_VERTICAL, dyn>CoP, SoyuzCMAerodynamicCoeffVert, 2.142, 0, D1 ); hHoriAirfoil = CreateAirfoil2( LIFT_HORIZONTAL, dyn>CoP, SoyuzCMAerodynamicCoeffHoriz, 2.142, 0, D1 ); With this the capsule is aerdoynamicly stable and flies just as it should. However the LifttoDrag ratio is at all points of the simulation half the size it should be (compared to the values given by literature) . Calcluating the drag manually for some points of the trajectory shows, that the drag is twice as high as supposed. Thus the Lift is just fine as it is. Since the difference in Lift and Drag calculation is simply the coefficient I tried working with them to solve the problem. Dividing the drag coefficient in half for each airfoil as well as setting one of them to zero seems to effect the whole aerodynmaic model causing the capsule to fly at a high slip angle of ~45° rather than the supposed one of ~0°. The angle of attack is also way higher than the old value of 17,5° (for the working model with twice the drag) Does anyone have an idea of why this is happening and how to solve the problem? I personally have run out of ideas. Changes in the chord length and wing span have, as expected, no influence on the LifttoDrag ratio. Thanks in advance for your time and effort 
03292019, 10:01 AM  #2 
Not funny anymore

What is your lift function?

03292019, 11:26 AM  #3 
Orbinaut

What do you mean by lift function?
The Coefficients are calculated this way: Code:
/* Mach Number Steps of Reference Points */ const double aero_MAref[6] = { 0.60, 0.95, 1.10, 1.78, 2.52, 5.96 }; /* Angle of Attack (AoA) Steps of Reference Points */ const double aero_AOAref[7] = { 0.0*RAD, 5.0*RAD, 10.0*RAD, 15.0*RAD, 20.0*RAD, 25.0*RAD, 30.0*RAD }; /* Axial Force Coefficient (This is NOT equal to Drag) for aoa = { 0, 5, 10, 15, 20, 25, 30 } deg and for all reference mach numbers. */ const double aero_cx[6][7] = { /* aoa 0 5 10 15 20 25 30 deg*/ { 1.01, 1.00, 1.00, 0.975, 0.958, 0.921, 0.879 }, // Ma = 0.60 { 1.27, 1.27, 1.25, 1.22, 1.21, 1.18, 1.15 }, // Ma = 0.95 { 1.40, 1.40, 1.39, 1.37, 1.34, 1.29, 1.25 }, // Ma = 1.10 { 1.52, 1.52, 1.50, 1.46, 1.43, 1.38, 1.31 }, // Ma = 1.78 { 1.52, 1.52, 1.48, 1.46, 1.43, 1.38, 1.31 }, // Ma = 2.52 { 1.57, 1.54, 1.51, 1.44, 1.38, 1.28, 1.18 }, // Ma = 5.96 }; /* Normal Force Coefficient (This is NOT equal to Lift) for aoa = { 0, 5, 10, 15, 20, 25, 30 } deg and for all reference mach numbers. */ const double aero_cz[6][7] = { /* aoa 0 5 10 15 20 25 30 deg*/ { 0.00, 0.04, 0.06, 0.07, 0.06, 0.04, 0.02 }, // Ma = 0.60 { 0.00, 0.03, 0.04, 0.01, 0.01, 0.07, 0.14 }, // Ma = 0.95 { 0.00, 0.03, 0.04, 0.00, 0.04, 0.09, 0.13 }, // Ma = 1.10 { 0.00, 0.03, 0.01, 0.08, 0.16, 0.26, 0.36 }, // Ma = 1.78 { 0.00, 0.04, 0.09, 0.14, 0.18, 0.24, 0.31 }, // Ma = 2.52 { 0.00, 0.03, 0.06, 0.09, 0.13, 0.19, 0.25 }, // Ma = 5.96 }; /* Lift to Drag Ratio for aoa = { 0, 5, 10, 15, 20, 25, 30 } deg and for all reference mach numbers. */ const double aero_ld[6][7] = { /* aoa 0 5 10 15 20 25 30 deg*/ { 0.00, 0.123, 0.249, 0.343, 0.443, 0.514, 0.614 }, // Ma = 0.60 { 0.00, 0.100, 0.209, 0.286, 0.363, 0.389, 0.420 }, // Ma = 0.95 { 0.00, 0.100, 0.200, 0.271, 0.323, 0.389, 0.446 }, // Ma = 1.10 { 0.00, 0.066, 0.171, 0.200, 0.229, 0.243, 0.257 }, // Ma = 1.78 { 0.00, 0.057, 0.114, 0.171, 0.229, 0.271, 0.300 }, // Ma = 2.52 { 0.00, 0.063, 0.129, 0.191, 0.257, 0.300, 0.329 }, // Ma = 5.96 }; /* Pitch Moment Coefficient for aoa = { 0, 5, 10, 15, 20, 25, 30 } deg and for all reference mach numbers. */ const double aero_cm[6][7] = { /* aoa 0 5 10 15 20 25 30 deg*/ { 0.0300, 0.0200, 0.0100, +0.0035, +0.0153, +0.0247, +0.0318 }, // Ma = 0.60 { 0.0359, 0.0241, 0.0141, 0.0059, +0.0012, +0.00711, +0.0094 }, // Ma = 0.95 { 0.0388, 0.0282, 0.0188, 0.0094, 0.0012, +0.0053, +0.0082 }, // Ma = 1.10 { 0.0447, 0.0400, 0.0300, 0.0176, 0.0041, +0.0082, +0.0112 }, // Ma = 1.78 { 0.0500, 0.0412, 0.0300, 0.0176, 0.0041, +0.0082, +0.0112 }, // Ma = 2.52 { 0.0535, 0.0435, 0.0329, 0.0224, 0.0118, +0.0000, +0.0100 } // Ma = 5.96 }; void SoyuzCMAerodynamicCoeffVert(double aoa, double Ma, double Re, double *cl, double *cm, double *cd) { /* We enter backwards, so rotate the AoA by PI */ double aoa2 = 0; if ( aoa < 0) aoa2 = (aoa + PI); if ( aoa > 0) aoa2 = (aoa  PI); double aoa3 = aoa2; //the actual AoA is needed for the correct lift and drag Coeffiecents aoa2 = abs (aoa2 ); // use only absolute value here /* Get lower AoA index */ int i_AoA1=1, i_AoA2=0; double d_AoA; // 0 .. 1 point between i_AoA1 and i_AoA2 for (int i=0; i<7; i++) if (aero_AOAref[i] < aoa2){ i_AoA1 = i;} else break; i_AoA2 = i_AoA1 + 1; d_AoA = (aoa2  aero_AOAref[i_AoA1]) / (aero_AOAref[i_AoA2]  aero_AOAref[i_AoA1]); /* If AoA is higher than 30 deg, limit.*/ if (aoa2 > aero_AOAref[6]) { i_AoA1 = 5; i_AoA2 = 6; d_AoA = 1; } if (aoa3 > 30 * (PI / 180)) aoa3 = 30 * (PI / 180); else if (aoa3 < 30 * (PI / 180)) aoa3 = 30 * (PI / 180); /* Get next lower mach number index */ int i_Ma1=1, i_Ma2=0; double d_Ma = 0; // 0 .. 1 point between i_Ma1 and i_Ma2 for (int i=0; i<6; i++) if (aero_MAref[i]<Ma) i_Ma1 = i; else break; i_Ma2 = i_Ma1 + 1; d_Ma = (Ma  aero_MAref[i_Ma1]) / (aero_MAref[i_Ma2]  aero_MAref[i_Ma1]); /* If mach number is higher, set to highest mach index instead of extrapolation. "[...] apparent changes of the arodynamic coefficients will occur for hypersonic Mach numbers up to 30, [...] the differences are probably low." [1] */ if (Ma > aero_MAref[5]) { i_Ma1 = 4; i_Ma2 = 5; d_Ma = 1; } /* If mach number is lower, set to the lowest mach index instead of extrapolation.*/ if (i_Ma1 == 1) { i_Ma1 = 0; i_Ma2 = 1; d_Ma = 0; } /*Calculation of the actual pure axial/normal Coefficient*/ /* Axial Component */ double cx11 = aero_cx[i_Ma1][i_AoA1]; // lower Ma, lower AoA double cx21 = aero_cx[i_Ma2][i_AoA1]; // higher Ma, lower AoA double cx12 = aero_cx[i_Ma1][i_AoA2]; // lower Ma, higher AoA double cx22 = aero_cx[i_Ma2][i_AoA2]; // higher Ma, higher AoA /* Linear interpolation over Ma */ double cx1 = cx11 + d_Ma * (cx21  cx11); // for lower AoA double cx2 = cx12 + d_Ma * (cx22  cx12); // for higher AoA /* Linear interpolation over AoA */ double cx = cx1 + d_AoA * (cx2  cx1); /* Normal Component */ double cz11 = aero_cz[i_Ma1][i_AoA1]; // lower Ma, lower AoA double cz21 = aero_cz[i_Ma2][i_AoA1]; // higher Ma, lower AoA double cz12 = aero_cz[i_Ma1][i_AoA2]; // lower Ma, higher AoA double cz22 = aero_cz[i_Ma2][i_AoA2]; // higher Ma, higher AoA /* Linear interpolation over Ma */ double cz1 = cz11 + d_Ma * (cz21  cz11); // for lower AoA double cz2 = cz12 + d_Ma * (cz22  cz12); // for higher AoA /* Linear interpolation over AoA */ double cz = cz1 + d_AoA * (cz2  cz1); /*Drag and Lift Coefficient*/ *cd = cx * cos(aoa3) + cz * sin(aoa3); //*cd = *cd/2; *cl = cx * sin(aoa3) + cz * cos(aoa3); *cl *= 1; // we fly backward in our coordinate system, so reverse sign /* Pitch Moment Coefficient */ double cm11 = aero_cm[i_Ma1][i_AoA1]; // lower Ma, lower AoA double cm21 = aero_cm[i_Ma2][i_AoA1]; // higher Ma, lower AoA double cm12 = aero_cm[i_Ma1][i_AoA2]; // lower Ma, higher AoA double cm22 = aero_cm[i_Ma2][i_AoA2]; // higher Ma, higher AoA /* Linear interpolation over Ma */ double cm1 = cm11 + d_Ma * (cm21  cm11); // for lower AoA double cm2 = cm12 + d_Ma * (cm22  cm12); // for higher AoA /* Linear interpolation over AoA */ *cm = cm1 + d_AoA * (cm2  cm1); *cm *= 1; // we fly backward in our coordinate system, so reverse sign if (aoa < 0) *cm *= 1; // for negative angles of attack, use different sign } void SoyuzCMAerodynamicCoeffHoriz(double aoa, double Ma, double Re, double *cl, double *cm, double *cd) { /* We enter backwards, so rotate the AoA by PI */ double aoa2 = 0; if (aoa < 0) aoa2 = (aoa + PI); if (aoa > 0) aoa2 = (aoa  PI); double aoa3 = aoa2; //the actual AoA is needed for the correct lift and drag Coeffiecents aoa2 = abs(aoa2); // use only absolute value here /* Get lower AoA index */ int i_AoA1 = 1, i_AoA2 = 0; double d_AoA; // 0 .. 1 point between i_AoA1 and i_AoA2 for (int i = 0; i<7; i++) if (aero_AOAref[i] < aoa2){ i_AoA1 = i; } else break; i_AoA2 = i_AoA1 + 1; d_AoA = (aoa2  aero_AOAref[i_AoA1]) / (aero_AOAref[i_AoA2]  aero_AOAref[i_AoA1]); /* If AoA is higher than 30 deg, limit.*/ if (aoa2 > aero_AOAref[6]) { i_AoA1 = 5; i_AoA2 = 6; d_AoA = 1; } if (aoa3 > 30 * (PI / 180)) aoa3 = 30 * (PI / 180); else if (aoa3 < 30 * (PI / 180)) aoa3 = 30 * (PI / 180); /* Get next lower mach number index i_Ma==1 oder i_Ma2>5 bedeutet extrapolation */ int i_Ma1 = 1, i_Ma2 = 0; double d_Ma = 0; // 0 .. 1 point between i_Ma1 and i_Ma2 for (int i = 0; i<6; i++) if (aero_MAref[i]<Ma) i_Ma1 = i; else break; i_Ma2 = i_Ma1 + 1; d_Ma = (Ma  aero_MAref[i_Ma1]) / (aero_MAref[i_Ma2]  aero_MAref[i_Ma1]); /* If mach number is higher, set to highest mach index instead of extrapolation. "[...] apparent changes of the arodynamic coefficients will occur for hypersonic Mach numbers up to 30, [...] the differences are probably low." [1] */ if (Ma > aero_MAref[5]) { i_Ma1 = 4; i_Ma2 = 5; d_Ma = 1; } /* If mach number is lower, set to the lowest mach index instead of extrapolation. */ if (i_Ma1 == 1) { i_Ma1 = 0; i_Ma2 = 1; d_Ma = 0; } /* Axial Component (Drag) */ double cx11 = aero_cx[i_Ma1][i_AoA1]; // lower Ma, lower AoA double cx21 = aero_cx[i_Ma2][i_AoA1]; // higher Ma, lower AoA double cx12 = aero_cx[i_Ma1][i_AoA2]; // lower Ma, higher AoA double cx22 = aero_cx[i_Ma2][i_AoA2]; // higher Ma, higher AoA /* Linear interpolation over Ma */ double cx1 = cx11 + d_Ma * (cx21  cx11); // for lower AoA double cx2 = cx12 + d_Ma * (cx22  cx12); // for higher AoA /* Linear interpolation over AoA */ double cx = cx1 + d_AoA * (cx2  cx1); /* Normal Component (Lift) */ double cz11 = aero_cz[i_Ma1][i_AoA1]; // lower Ma, lower AoA double cz21 = aero_cz[i_Ma2][i_AoA1]; // higher Ma, lower AoA double cz12 = aero_cz[i_Ma1][i_AoA2]; // lower Ma, higher AoA double cz22 = aero_cz[i_Ma2][i_AoA2]; // higher Ma, higher AoA /* Linear interpolation over Ma */ double cz1 = cz11 + d_Ma * (cz21  cz11); // for lower AoA double cz2 = cz12 + d_Ma * (cz22  cz12); // for higher AoA /* Linear interpolation over AoA */ double cz = cz1 + d_AoA * (cz2  cz1); /*Drag and Lift Coefficient calculation*/ *cd = cx * cos(aoa3) + cz * sin(aoa3); //*cd = *cd/2; *cl = cx * sin(aoa3) + cz * cos(aoa3); *cl *= 1; // we fly backward in our coordinate system, so reverse sign /* Pitch Moment Coefficient */ double cm11 = aero_cm[i_Ma1][i_AoA1]; // lower Ma, lower AoA double cm21 = aero_cm[i_Ma2][i_AoA1]; // higher Ma, lower AoA double cm12 = aero_cm[i_Ma1][i_AoA2]; // lower Ma, higher AoA double cm22 = aero_cm[i_Ma2][i_AoA2]; // higher Ma, higher AoA /* Linear interpolation over Ma */ double cm1 = cm11 + d_Ma * (cm21  cm11); // for lower AoA double cm2 = cm12 + d_Ma * (cm22  cm12); // for higher AoA /* Linear interpolation over AoA */ *cm = cm1 + d_AoA * (cm2  cm1); *cm *= 1; // we fly backward in our coordinate system, so reverse sign if (aoa < 0) *cm *= 1; // for negative angles of attack, use different sign } 
03292019, 12:03 PM  #4 
Orbiter Founder

This is a stab in the dark, but is it possible that the problem is that both your airfoils are generating drag, but (effectively) only one is generating nonzero lift?

03292019, 12:41 PM  #5 
Orbinaut

Quote:
Or do you mean, that the two airfoils create the amount of lift that only one should (compared to the drag they produce)? Is there a way to get the value of the Wing Area that's calculated from the cross sections that the airfoil works with? Like a GetWingArea() type of function that returns the wing area for that moment with the current attitute of the sapcecraft? Then I could manually check, if at least the amount of lift is correct. Changing the Lift would only be a temporary solution for me, since the drag value is obviously not what it should be and, as this is for my Bachlor Thesis, it's fairly important that I work with the real figures. Thanks for the Idea though 
03292019, 12:50 PM  #6 
Not funny anymore

Quote:

03292019, 01:51 PM  #7 
Orbinaut

Quote:
The (real) Soyuz is not only aerodynamicly stable but also stabilized by a control system. The Slip angle during reentry is very small and doen't really exceed + 2°. That induces very little Lift. If I was willing to simply "ignore" that small horitzontal lift, couldn't I remove the Airfoil with Horizontal lift? Then the drag coming from just one arifoil should have the value I need and the Lift almost stays the same, as most of it is vertical anyways. 
03292019, 02:14 PM  #8 
Orbiter Founder

The point I was making is that if you have a body that is rotationally symmetric w.r.t. its longitudinal axis (as I guess your capsule approximately is) then your two airfoils probably will also have identical and symmetric lift characteristics.
Therefore the body's rotation along its longitudinal axis shouldn't change the lift (except that is rotates in the local vessel frame). So without loss of generality, in the case of negligible slip angle during entry, you can assume that the capsule is rotated such that one of the airfoils is generating all the lift, while the other is generating zero lift. Rotate the capsule by 90 degrees, and it's the other airfoil that is now generating the lift. Anywhere in between, both airfoils generate part of the lift, but the magnitude of the resulting lift vector remains the same. This is what I meant when saying that effectively only one is generating lift. Whereas the drag force from both airfoils is added up in all orientations. So I would argue that the desired drag characteristics must be split evenly over the two airfoils (except induced drag, since that is directly linked to the lift), while each of the airfoils must be defined with the full lift characteristics. Does this make sense? Not sure if my argument is very clear. 
03292019, 02:38 PM  #9 
Not funny anymore

I am not sure if it would work with just one foil  the horizontal lift function also has a stabilizing function. You could describe it as zerolift for simplification, but I would not ignore the torque function.

03292019, 05:06 PM  #10 
Orbinaut

Quote:
Quote:
Thank you very much for your help, I'll let you know if the changes helped  Post added at 05:06 PM  Previous post was at 03:28 PM  Quote:
Once again thank you very much, this helped me greatly And now, onwords to the next thing to fix 
03292019, 05:23 PM  #11 
Orbiter Founder

Good to hear that it works now!
Although the fact that lift and drag for the horizontal foil are now zero makes me slightly ... uneasy What happens if you remove that airfoil altogether? Does it behave the same? Just to confirm: Your lift, moment and drag coefficients _are_ symmetric around AOA=0, correct? That is, cl(aoa=0) = 0 cl(aoa) = cl(aoa) same for moment coefficients, and cd(aoa) = cd(aoa) In that case, if your slip angle is 0, lift will be 0 as well, so it should behave the same as if you force cl=0 outright. What happens if you flip the capsule on its side for the reentry, so that the roles of the horizontal and vertical airfoils are swapped? Does all hell break loose? It might be worth putting a debug output into your airfoil coefficient callback function that echos the input AOA and resulting output cl values on the screen and to monitor that during entry to make sure that the values make sense. 

Tags 
airfoil, drag, sdk 
Thread Tools  


Quick Links  Need Help? 