44 #include "FGFDMExec.h" 45 #include "FGAerodynamics.h" 46 #include "input_output/FGPropertyManager.h" 47 #include "input_output/FGXMLElement.h" 53 IDENT(IdSrc,
"$Id: FGAerodynamics.cpp,v 1.60 2016/05/30 19:05:58 bcoconni Exp $");
54 IDENT(IdHdr,ID_AERODYNAMICS);
63 Name =
"FGAerodynamics";
73 AxisIdx[
"NORMAL"] = 2;
81 AeroFunctions =
new AeroFunctionArray[6];
82 AeroFunctionsAtCG =
new AeroFunctionArray[6];
84 impending_stall = stall_hyst = 0.0;
85 alphaclmin = alphaclmax = 0.0;
86 alphaclmin0 = alphaclmax0 = 0.0;
87 alphahystmin = alphahystmax = 0.0;
90 bi2vel = ci2vel = 0.0;
92 vDeltaRP.InitMatrix();
106 for (j=0; j<AeroFunctions[i].size(); j++)
107 delete AeroFunctions[i][j];
109 for (j=0; j<AeroFunctionsAtCG[i].size(); j++)
110 delete AeroFunctionsAtCG[i][j];
112 delete[] AeroFunctions;
113 delete[] AeroFunctionsAtCG;
122 bool FGAerodynamics::InitModel(
void)
124 if (!FGModel::InitModel())
return false;
126 impending_stall = stall_hyst = 0.0;
127 alphaclmin = alphaclmin0;
128 alphaclmax = alphaclmax0;
129 alphahystmin = alphahystmax = 0.0;
132 bi2vel = ci2vel = 0.0;
134 vDeltaRP.InitMatrix();
135 vForces.InitMatrix();
136 vMoments.InitMatrix();
145 if (Holding)
return false;
147 unsigned int axis_ctr;
148 const double twovel=2*in.Vt;
153 if ( in.Qbar > 1.0) {
157 clsq = (vFw(eLift) + vFwAtCG(eLift))/ (in.Wingarea*in.Qbar);
166 bi2vel = in.Wingspan / twovel;
167 ci2vel = in.Wingchord / twovel;
169 alphaw = in.Alpha + in.Wingincidence;
170 qbar_area = in.Wingarea * in.Qbar;
172 if (alphaclmax != 0) {
173 if (in.Alpha > 0.85*alphaclmax) {
174 impending_stall = 10*(in.Alpha/alphaclmax - 0.85);
180 if (alphahystmax != 0.0 && alphahystmin != 0.0) {
181 if (in.Alpha > alphahystmax) {
183 }
else if (in.Alpha < alphahystmin) {
189 vFwAtCG.InitMatrix();
190 vFnative.InitMatrix();
191 vFnativeAtCG.InitMatrix();
193 for (axis_ctr = 0; axis_ctr < 3; ++axis_ctr) {
194 AeroFunctionArray::iterator f;
196 AeroFunctionArray* array = &AeroFunctions[axis_ctr];
197 for (f=array->begin(); f != array->end(); ++f) {
202 (*f)->cacheValue(
true);
203 vFnative(axis_ctr+1) += (*f)->GetValue();
206 array = &AeroFunctionsAtCG[axis_ctr];
207 for (f=array->begin(); f != array->end(); ++f) {
208 (*f)->cacheValue(
true);
209 vFnativeAtCG(axis_ctr+1) += (*f)->GetValue();
228 vFw = in.Tb2w*vFnative;
230 vFw(eDrag)*=-1; vFw(eLift)*=-1;
232 vFwAtCG = in.Tb2w*vFnativeAtCG;
233 vForcesAtCG = vFnativeAtCG;
234 vFwAtCG(eDrag)*=-1; vFwAtCG(eLift)*=-1;
238 vFw(eDrag)*=-1; vFw(eLift)*=-1;
239 vForces = in.Tw2b*vFw;
240 vFw(eDrag)*=-1; vFw(eLift)*=-1;
242 vFwAtCG = vFnativeAtCG;
243 vFwAtCG(eDrag)*=-1; vFwAtCG(eLift)*=-1;
244 vForcesAtCG = in.Tw2b*vFwAtCG;
245 vFwAtCG(eDrag)*=-1; vFwAtCG(eLift)*=-1;
248 vFw = in.Tb2w*vFnative;
249 vFnative(eX)*=-1; vFnative(eZ)*=-1;
252 vFwAtCG = in.Tb2w*vFnativeAtCG;
253 vFnativeAtCG(eX)*=-1; vFnativeAtCG(eZ)*=-1;
254 vForcesAtCG = vFnativeAtCG;
257 cerr << endl <<
" A proper axis type has NOT been selected. Check " 258 <<
"your aerodynamics definition." << endl;
263 if ( fabs(vFw(eDrag) + vFwAtCG(eDrag)) > 0.0)
264 lod = fabs( (vFw(eLift) + vFwAtCG(eLift))/ (vFw(eDrag) + vFwAtCG(eDrag)));
271 if (AeroRPShift) vDeltaRP(eX) = AeroRPShift->
GetValue()*in.Wingchord;
273 vDXYZcg(eX) = in.RPBody(eX) - vDeltaRP(eX);
274 vDXYZcg(eY) = in.RPBody(eY) + vDeltaRP(eY);
275 vDXYZcg(eZ) = in.RPBody(eZ) - vDeltaRP(eZ);
277 vMomentsMRC.InitMatrix();
279 for (axis_ctr = 0; axis_ctr < 3; axis_ctr++) {
280 AeroFunctionArray* array = &AeroFunctions[axis_ctr+3];
281 for (AeroFunctionArray::iterator f=array->begin(); f != array->end(); ++f) {
286 (*f)->cacheValue(
true);
287 vMomentsMRC(axis_ctr+1) += (*f)->GetValue();
290 vMoments = vMomentsMRC + vDXYZcg*vForces;
292 vForces += vForcesAtCG;
293 vFnative += vFnativeAtCG;
306 string scratch_unit=
"";
307 Element *temp_element, *axis_element, *function_element;
315 DetermineAxisSystem(document);
319 if ((temp_element = document->
FindElement(
"alphalimits"))) {
321 if (scratch_unit.empty()) scratch_unit =
"RAD";
324 alphaclmin = alphaclmin0;
325 alphaclmax = alphaclmax0;
328 if ((temp_element = document->
FindElement(
"hysteresis_limits"))) {
330 if (scratch_unit.empty()) scratch_unit =
"RAD";
335 if ((temp_element = document->
FindElement(
"aero_ref_pt_shift_x"))) {
336 function_element = temp_element->
FindElement(
"function");
337 AeroRPShift =
new FGFunction(PropertyManager, function_element);
341 while (axis_element) {
342 AeroFunctionArray ca;
343 AeroFunctionArray ca_atCG;
345 function_element = axis_element->
FindElement(
"function");
346 while (function_element) {
348 bool apply_at_cg =
false;
350 if (function_element->
GetAttributeValue(
"apply_at_cg") ==
"true") apply_at_cg =
true;
354 ca.push_back(
new FGFunction(PropertyManager, function_element) );
355 }
catch (
const string& str) {
356 cerr << endl <<
fgred <<
"Error loading aerodynamic function in " 357 << current_func_name <<
":" << str <<
" Aborting." <<
reset << endl;
362 ca_atCG.push_back(
new FGFunction(PropertyManager, function_element) );
363 }
catch (
const string& str) {
364 cerr << endl <<
fgred <<
"Error loading aerodynamic function in " 365 << current_func_name <<
":" << str <<
" Aborting." <<
reset << endl;
371 AeroFunctions[AxisIdx[axis]] = ca;
372 AeroFunctionsAtCG[AxisIdx[axis]] = ca_atCG;
376 PostLoad(document, PropertyManager);
391 void FGAerodynamics::DetermineAxisSystem(
Element* document)
395 while (axis_element) {
397 if (axis ==
"LIFT" || axis ==
"DRAG") {
398 if (axisType == atNone) axisType = atLiftDrag;
399 else if (axisType != atLiftDrag) {
400 cerr << endl <<
" Mixed aerodynamic axis systems have been used in the" 401 <<
" aircraft config file. (LIFT DRAG)" << endl;
403 }
else if (axis ==
"SIDE") {
404 if (axisType != atNone && axisType != atLiftDrag && axisType != atAxialNormal) {
405 cerr << endl <<
" Mixed aerodynamic axis systems have been used in the" 406 <<
" aircraft config file. (SIDE)" << endl;
408 }
else if (axis ==
"AXIAL" || axis ==
"NORMAL") {
409 if (axisType == atNone) axisType = atAxialNormal;
410 else if (axisType != atAxialNormal) {
411 cerr << endl <<
" Mixed aerodynamic axis systems have been used in the" 412 <<
" aircraft config file. (NORMAL AXIAL)" << endl;
414 }
else if (axis ==
"X" || axis ==
"Y" || axis ==
"Z") {
415 if (axisType == atNone) axisType = atBodyXYZ;
416 else if (axisType != atBodyXYZ) {
417 cerr << endl <<
" Mixed aerodynamic axis systems have been used in the" 418 <<
" aircraft config file. (XYZ)" << endl;
420 }
else if (axis !=
"ROLL" && axis !=
"PITCH" && axis !=
"YAW") {
421 cerr << endl <<
" An unknown axis type, " << axis <<
" has been specified" 422 <<
" in the aircraft configuration file." << endl;
427 if (axisType == atNone) {
428 axisType = atLiftDrag;
429 cerr << endl <<
" The aerodynamic axis system has been set by default" 430 <<
" to the Lift/Side/Drag system." << endl;
438 string AeroFunctionStrings =
"";
439 bool firstime =
true;
440 unsigned int axis, sd;
442 for (axis = 0; axis < 6; axis++) {
443 for (sd = 0; sd < AeroFunctions[axis].size(); sd++) {
447 AeroFunctionStrings += delimeter;
449 AeroFunctionStrings += AeroFunctions[axis][sd]->GetName();
455 if (FunctionStrings.size() > 0) {
456 if (AeroFunctionStrings.size() > 0) {
457 AeroFunctionStrings += delimeter + FunctionStrings;
459 AeroFunctionStrings = FunctionStrings;
463 return AeroFunctionStrings;
472 for (
unsigned int axis = 0; axis < 6; axis++) {
473 for (
unsigned int sd = 0; sd < AeroFunctions[axis].size(); sd++) {
474 if (buf.tellp() > 0) buf << delimeter;
475 buf << AeroFunctions[axis][sd]->GetValue();
481 if (FunctionValues.size() > 0) {
482 if (buf.str().size() > 0) {
483 buf << delimeter << FunctionValues;
485 buf << FunctionValues;
494 void FGAerodynamics::bind(
void)
499 PropertyManager->
Tie(
"forces/fby-aero-lbs",
this, 2, (PMF)&FGAerodynamics::GetForces);
500 PropertyManager->
Tie(
"forces/fbz-aero-lbs",
this, 3, (PMF)&FGAerodynamics::GetForces);
502 PropertyManager->
Tie(
"moments/m-aero-lbsft",
this, 2, (PMF)&FGAerodynamics::GetMoments);
503 PropertyManager->
Tie(
"moments/n-aero-lbsft",
this, 3, (PMF)&FGAerodynamics::GetMoments);
505 PropertyManager->
Tie(
"forces/fwy-aero-lbs",
this, 2, (PMF)&FGAerodynamics::GetvFw);
506 PropertyManager->
Tie(
"forces/fwz-aero-lbs",
this, 3, (PMF)&FGAerodynamics::GetvFw);
509 PropertyManager->
Tie(
"aero/qbar-area", &qbar_area);
510 PropertyManager->
Tie(
"aero/alpha-max-rad",
this, &FGAerodynamics::GetAlphaCLMax, &FGAerodynamics::SetAlphaCLMax,
true);
511 PropertyManager->
Tie(
"aero/alpha-min-rad",
this, &FGAerodynamics::GetAlphaCLMin, &FGAerodynamics::SetAlphaCLMin,
true);
512 PropertyManager->
Tie(
"aero/bi2vel",
this, &FGAerodynamics::GetBI2Vel);
513 PropertyManager->
Tie(
"aero/ci2vel",
this, &FGAerodynamics::GetCI2Vel);
514 PropertyManager->
Tie(
"aero/alpha-wing-rad",
this, &FGAerodynamics::GetAlphaW);
515 PropertyManager->
Tie(
"systems/stall-warn-norm",
this, &FGAerodynamics::GetStallWarn);
516 PropertyManager->
Tie(
"aero/stall-hyst-norm",
this, &FGAerodynamics::GetHysteresisParm);
538 void FGAerodynamics::Debug(
int from)
540 if (debug_lvl <= 0)
return;
546 cout << endl <<
" Aerodynamics (Lift|Side|Drag axes):" << endl << endl;
548 case (atAxialNormal):
549 cout << endl <<
" Aerodynamics (Axial|Side|Normal axes):" << endl << endl;
552 cout << endl <<
" Aerodynamics (X|Y|Z axes):" << endl << endl;
555 cout << endl <<
" Aerodynamics (undefined axes):" << endl << endl;
560 if (debug_lvl & 2 ) {
561 if (from == 0) cout <<
"Instantiated: FGAerodynamics" << endl;
562 if (from == 1) cout <<
"Destroyed: FGAerodynamics" << endl;
564 if (debug_lvl & 4 ) {
566 if (debug_lvl & 8 ) {
568 if (debug_lvl & 16) {
570 if (debug_lvl & 64) {
572 cout << IdSrc << endl;
573 cout << IdHdr << endl;
double GetLoD(void) const
Retrieves the lift over drag ratio.
std::string GetAttributeValue(const std::string &key)
Retrieves an attribute.
bool HasAttribute(const std::string &key)
Determines if an element has the supplied attribute.
std::string GetAeroFunctionValues(const std::string &delimeter) const
Gets the aero function values.
static char reset[5]
resets text properties
Encapsulates the aerodynamic calculations.
Element * FindElement(const std::string &el="")
Searches for a specified element.
double FindElementValueAsNumberConvertFromTo(const std::string &el, const std::string &supplied_units, const std::string &target_units)
Searches for the named element and converts and returns the data belonging to it. ...
std::string GetFunctionStrings(const std::string &delimeter) const
Gets the strings for the current set of functions.
double GetValue(void) const
Retrieves the value of the function object.
const FGColumnVector3 & GetvFw(void) const
Retrieves the aerodynamic forces in the wind axes.
static char fgred[6]
red text
virtual bool Run(bool Holding)
Runs the model; called by the Executive.
double GetClSquared(void) const
Retrieves the square of the lift coefficient.
void Tie(const std::string &name, bool *pointer, bool useDefault=true)
Tie a property to an external bool variable.
const FGColumnVector3 & GetForces(void) const
Gets the total aerodynamic force vector.
std::string GetAeroFunctionStrings(const std::string &delimeter) const
Gets the strings for the current set of aero functions.
Base class for all scheduled JSBSim models.
Represents a mathematical function.
~FGAerodynamics()
Destructor.
Element * FindNextElement(const std::string &el="")
Searches for the next element as specified.
Encapsulates the JSBSim simulation executive.
virtual bool Load(Element *el)
Loads this model.
bool Load(Element *element)
Loads the Aerodynamics model.
std::string GetFunctionValues(const std::string &delimeter) const
Gets the function values.
bool Run(bool Holding)
Runs the Aerodynamics model; called by the Executive Can pass in a value indicating if the executive ...
const FGColumnVector3 & GetMoments(void) const
Gets the total aerodynamic moment vector about the CG.