From 3f0ef99c71058595c2257d87c911d60f22ad349b Mon Sep 17 00:00:00 2001 From: Henning Stahlke Date: Mon, 3 Apr 2017 21:42:41 +0200 Subject: [PATCH] YASim: add MAC (mean aerodynamic chord) and c.g. calculations. --- src/FDM/YASim/Airplane.cpp | 81 ++++++++++++++++++++++++------------ src/FDM/YASim/Airplane.hpp | 58 ++++++++++++++++---------- src/FDM/YASim/FGFDM.cpp | 10 ++++- src/FDM/YASim/FGFDM.hpp | 1 + src/FDM/YASim/Wing.cpp | 7 ++++ src/FDM/YASim/Wing.hpp | 40 +++++++++++------- src/FDM/YASim/yasim-test.cpp | 22 +++++++--- 7 files changed, 146 insertions(+), 73 deletions(-) diff --git a/src/FDM/YASim/Airplane.cpp b/src/FDM/YASim/Airplane.cpp index cdc9ef8cc..4098ee4b6 100644 --- a/src/FDM/YASim/Airplane.cpp +++ b/src/FDM/YASim/Airplane.cpp @@ -57,8 +57,10 @@ Airplane::Airplane() _failureMsg = 0; _wingsN = 0; - _cgMaxX = -1e6; - _cgMinX = 1e6; + _cgMax = -1e6; + _cgMin = 1e6; + _cgDesiredMax = 0.33f; // FIXME find reasonable default value + _cgDesiredMin = 0.1f; // FIXME find reasonable default value } Airplane::~Airplane() @@ -115,7 +117,7 @@ void Airplane::calcFuelWeights() } } -void Airplane::getPilotAccel(float* out) +const void Airplane::getPilotAccel(float* out) { State* s = _model.getState(); @@ -125,14 +127,13 @@ void Airplane::getPilotAccel(float* out) Math::vmul33(s->orient, out, out); out[0] = -out[0]; - // The regular acceleration - float tmp[3]; + float acceleration[3]; // Convert to aircraft coordinates - Math::vmul33(s->orient, s->acc, tmp); - tmp[1] = -tmp[1]; - tmp[2] = -tmp[2]; + Math::vmul33(s->orient, s->acc, acceleration); + acceleration[1] = -acceleration[1]; + acceleration[2] = -acceleration[2]; - Math::add3(tmp, out, out); + Math::add3(acceleration, out, out); // FIXME: rotational & centripetal acceleration needed } @@ -240,10 +241,6 @@ void Airplane::addGear(Gear* gear) g->gear = gear; g->surf = 0; _gears.add(g); - float pos[3]; - g->gear->getPosition(pos); - if (pos[0] > _cgMaxX) _cgMaxX = pos[0]; - if (pos[0] < _cgMinX) _cgMinX = pos[0]; } void Airplane::addThruster(Thruster* thruster, float mass, float* cg) @@ -371,6 +368,9 @@ float Airplane::compileWing(Wing* w) _wingsN->getNode("wing-area", true)->setFloatValue(w->getArea()); _wingsN->getNode("aspect-ratio", true)->setFloatValue(w->getAspectRatio()); _wingsN->getNode("standard-mean-chord", true)->setFloatValue(w->getSMC()); + _wingsN->getNode("mac", true)->setFloatValue(w->getMAC()); + _wingsN->getNode("mac-x", true)->setFloatValue(w->getMACx()); + _wingsN->getNode("mac-y", true)->setFloatValue(w->getMACy()); } float wgt = 0; @@ -582,7 +582,6 @@ void Airplane::compile() RigidBody* body = _model.getBody(); int firstMass = body->numMasses(); SGPropertyNode_ptr baseN = fgGetNode("/fdm/yasim/model/wings", true); - SGPropertyNode_ptr n; // Generate the point masses for the plane. Just use unitless // numbers for a first pass, then go back through and rescale to @@ -594,6 +593,15 @@ void Airplane::compile() { if (baseN != 0) _wingsN = baseN->getChild("wing", 0, true); aeroWgt += compileWing(_wing); + + // convert % to absolute x coordinates + _cgDesiredFront = _wing->getMACx() - _wing->getMAC()*_cgDesiredMin; + _cgDesiredAft = _wing->getMACx() - _wing->getMAC()*_cgDesiredMax; + if (baseN != 0) { + SGPropertyNode_ptr n = fgGetNode("/fdm/yasim/model", true); + n->getNode("cg-range-front", true)->setFloatValue(_cgDesiredFront); + n->getNode("cg-range-aft", true)->setFloatValue(_cgDesiredAft); + } } if (_tail) { @@ -648,33 +656,35 @@ void Airplane::compile() ThrustRec* tr = (ThrustRec*)_thrusters.get(i); tr->handle = _model.addThruster(tr->thruster); } - - // Ground effect - // If a double tapered wing is modelled with wing and mstab, wing must - // be outboard to get correct wingspan. + if(_wing) { - float gepos[3]; + // Ground effect + // If a double tapered wing is modelled with wing and mstab, wing must + // be outboard to get correct wingspan. + float pos[3]; float gespan = 0; gespan = _wing->getSpan(); - _wing->getBase(gepos); + _wing->getBase(pos); if(!isVersionOrNewer( Version::YASIM_VERSION_2017_2 )) { //old code //float span = _length * Math::cos(_sweep) * Math::cos(_dihedral); //span = 2*(span + Math::abs(_base[2])); - gespan -= 2*gepos[1]; // cut away base (y-distance) - gespan += 2*Math::abs(gepos[2]); // add (wrong) z-distance + gespan -= 2*pos[1]; // cut away base (y-distance) + gespan += 2*Math::abs(pos[2]); // add (wrong) z-distance } if (baseN != 0) baseN->getChild("wing", 0)->getNode("gnd-eff-span", true)->setFloatValue(gespan); // where does the hard coded factor 0.15 come from? - _model.setGroundEffect(gepos, gespan, 0.15f); + _model.setGroundEffect(pos, gespan, 0.15f); } - + // solve function below resets failure message // so check if we have any problems and abort here if (_failureMsg) return; solveGear(); + calculateCGHardLimits(); + if(_wing && _tail) solve(); else { @@ -742,11 +752,24 @@ void Airplane::solveGear() } } +void Airplane::calculateCGHardLimits() +{ + _cgMax = -1e6; + _cgMin = 1e6; + for (int i = 0; i < _gears.size(); i++) { + GearRec* gr = (GearRec*)_gears.get(i); + float pos[3]; + gr->gear->getPosition(pos); + if (pos[0] > _cgMax) _cgMax = pos[0]; + if (pos[0] < _cgMin) _cgMin = pos[0]; + } +} + void Airplane::initEngines() { for(int i=0; i<_thrusters.size(); i++) { ThrustRec* tr = (ThrustRec*)_thrusters.get(i); - tr->thruster->init(); + tr->thruster->init(); } } @@ -835,7 +858,6 @@ void Airplane::runApproach() Math::vmul33(_approachState.orient, wind, wind); setFuelFraction(_approachFuel); - setupWeights(true); // Run the thrusters until they get to a stable setting. FIXME: @@ -1094,4 +1116,11 @@ void Airplane::solveHelicopter() } +const float Airplane::getCGMAC() +{ + float cg[3]; + _model.getBody()->getCG(cg); + return (_wing->getMACx() - cg[0]) / _wing->getMAC(); +} + }; // namespace yasim diff --git a/src/FDM/YASim/Airplane.hpp b/src/FDM/YASim/Airplane.hpp index 13dcdc836..690517bd6 100644 --- a/src/FDM/YASim/Airplane.hpp +++ b/src/FDM/YASim/Airplane.hpp @@ -33,11 +33,12 @@ public: void setPilotPos(float* pos) { Math::set3(pos, _pilotPos); } void getPilotPos(float* out) { Math::set3(_pilotPos, out); } - void getPilotAccel(float* out); + const void getPilotAccel(float* out); void setEmptyWeight(float weight) { _emptyWeight = weight; } void setWing(Wing* wing) { _wing = wing; } + Wing* getWing() { return _wing; } void setTail(Wing* tail) { _tail = tail; } void addVStab(Wing* vstab) { _vstabs.add(vstab); } @@ -64,56 +65,62 @@ public: void addSolutionWeight(bool approach, int idx, float wgt); - int numGear() { return _gears.size(); } + const int numGear() { return _gears.size(); } Gear* getGear(int g) { return ((GearRec*)_gears.get(g))->gear; } Hook* getHook() { return _model.getHook(); } - int numHitches() { return _hitches.size(); } + const int numHitches() { return _hitches.size(); } Hitch* getHitch(int h); Rotorgear* getRotorgear() { return _model.getRotorgear(); } Launchbar* getLaunchbar() { return _model.getLaunchbar(); } - int numThrusters() { return _thrusters.size(); } + const int numThrusters() { return _thrusters.size(); } Thruster* getThruster(int n) { return ((ThrustRec*)_thrusters.get(n))->thruster; } - int numTanks() { return _tanks.size(); } + const int numTanks() { return _tanks.size(); } void setFuelFraction(float frac); // 0-1, total amount of fuel // get fuel in kg - float getFuel(int tank) { return ((Tank*)_tanks.get(tank))->fill; } + const float getFuel(int tank) { return ((Tank*)_tanks.get(tank))->fill; } // set fuel in kg float setFuel(int tank, float fuel) { return ((Tank*)_tanks.get(tank))->fill = fuel; } // get fuel density in kg/m^3 - float getFuelDensity(int tank) { return ((Tank*)_tanks.get(tank))->density; } - float getTankCapacity(int tank) { return ((Tank*)_tanks.get(tank))->cap; } + const float getFuelDensity(int tank) { return ((Tank*)_tanks.get(tank))->density; } + const float getTankCapacity(int tank) { return ((Tank*)_tanks.get(tank))->cap; } void compile(); // generate point masses & such, then solve void initEngines(); void stabilizeThrust(); // Solution output values - int getSolutionIterations() { return _solutionIterations; } - float getDragCoefficient() { return _dragFactor; } - float getLiftRatio() { return _liftRatio; } - float getCruiseAoA() { return _cruiseAoA; } - float getTailIncidence() { return _tailIncidence; } - float getApproachElevator() { return _approachElevator.val; } + const int getSolutionIterations() { return _solutionIterations; } + const float getDragCoefficient() { return _dragFactor; } + const float getLiftRatio() { return _liftRatio; } + const float getCruiseAoA() { return _cruiseAoA; } + const float getTailIncidence() { return _tailIncidence; } + const float getApproachElevator() { return _approachElevator.val; } const char* getFailureMsg() { return _failureMsg; } static void setupState(const float aoa, const float speed, const float gla, yasim::State* s); // utility void loadApproachControls() { loadControls(_approachControls); } void loadCruiseControls() { loadControls(_cruiseControls); } - float getCGMinX() { return _cgMinX; } - float getCGMaxX() { return _cgMaxX; } + const float getCGHardLimitXMin() { return _cgMin; } // get min x-coordinate for c.g (from main gear) + const float getCGHardLimitXMax() { return _cgMax; } // get max x-coordinate for c.g (from nose gear) + const float getCGMAC(); // return c.g. x as fraction of MAC + // set desired range for C.G. in % of MAC, 0% = leading edge, 100% trailing edge + void setDesiredCGRangeInPercentOfMAC(float MACPercentMin, float MACPercentMax) { _cgDesiredMin = MACPercentMin; _cgDesiredMax = MACPercentMax; } + const float getCGSoftLimitXMin() { return _cgDesiredAft; } // get x-coordinate limit calculated from MAC and setCGRange values + const float getCGSoftLimitXMax() { return _cgDesiredFront; } // get x-coordinate limit calculated from MAC and setCGRange values + void setAutoBallast(bool allowed) { _autoBallast = allowed; } private: struct Tank { float pos[3]; float cap; float fill; - float density; int handle; }; + float density; int handle; }; struct Fuselage { float front[3], back[3], width, taper, mid, _cx, _cy, _cz, _idrag; - Vector surfs; }; + Vector surfs; }; struct GearRec { Gear* gear; Surface* surf; float wgt; }; struct ThrustRec { Thruster* thruster; - int handle; float cg[3]; float mass; }; + int handle; float cg[3]; float mass; }; struct Control { int control; float val; }; struct WeightRec { int handle; Surface* surf; }; struct SolveWeight { bool approach; int idx; float wgt; }; @@ -136,7 +143,8 @@ private: float normFactor(float f); void updateGearState(); void setupWeights(bool isApproach); - + void calculateCGHardLimits(); + Model _model; ControlMap _controls; @@ -187,9 +195,13 @@ private: Control _approachElevator; const char* _failureMsg; - // hard limits for cg from gear positions - float _cgMaxX; - float _cgMinX; + float _cgMax; // hard limits for cg from gear position + float _cgMin; // hard limits for cg from gear position + float _cgDesiredMax; // desired cg max in %MAC from config + float _cgDesiredMin; // desired cg min in %MAC from config + float _cgDesiredFront; // calculated desired cg x max + float _cgDesiredAft; // calculated desired cg x min + bool _autoBallast = false; }; }; // namespace yasim diff --git a/src/FDM/YASim/FGFDM.cpp b/src/FDM/YASim/FGFDM.cpp index d981eb516..3a604ef28 100644 --- a/src/FDM/YASim/FGFDM.cpp +++ b/src/FDM/YASim/FGFDM.cpp @@ -133,11 +133,12 @@ void FGFDM::init() // write some compile time information to property tree _yasimN->getNode("config-version",true)->setIntValue(_airplane.getVersion()); - _yasimN->getNode("model/cg-x-min",true)->setFloatValue(_airplane.getCGMinX()); - _yasimN->getNode("model/cg-x-max",true)->setFloatValue(_airplane.getCGMaxX()); + _yasimN->getNode("model/cg-x-min",true)->setFloatValue(_airplane.getCGHardLimitXMin()); + _yasimN->getNode("model/cg-x-max",true)->setFloatValue(_airplane.getCGHardLimitXMax()); // prepare nodes for write at runtime _cg_x = _yasimN->getNode("cg-x-m", true); + _cg_xmacN = _yasimN->getNode("cg-x-mac", true); _cg_y = _yasimN->getNode("cg-y-m", true); _cg_z = _yasimN->getNode("cg-z-m", true); _vxN = _yasimN->getNode("velocities/v-x", true); @@ -258,6 +259,10 @@ void FGFDM::startElement(const char* name, const XMLAttributes &atts) if( !_airplane.isVersionOrNewer( Version::YASIM_VERSION_CURRENT ) ) { SG_LOG(SG_FLIGHT, SG_DEV_ALERT, "This aircraft does not use the latest yasim configuration version."); } + _airplane.setDesiredCGRangeInPercentOfMAC(attrf(a, "cg-min", 0.1f), attrf(a, "cg-max", 0.3f)); //FIXME find reasonable defaults + if (attrb(a, "auto-ballast")) { + _airplane.setAutoBallast(true); + } } else if(eq(name, "approach")) { float spd, alt = 0; if (a->hasAttribute("speed")) { spd = attrf(a, "speed") * KTS2MPS; } @@ -691,6 +696,7 @@ void FGFDM::setOutputProperties(float dt) _cg_x->setFloatValue(cg[0]); _cg_y->setFloatValue(cg[1]); _cg_z->setFloatValue(cg[2]); + _cg_xmacN->setFloatValue(_airplane.getCGMAC()); State* s = _airplane.getModel()->getState(); float v[3], acc[3], rot[3], racc[3]; diff --git a/src/FDM/YASim/FGFDM.hpp b/src/FDM/YASim/FGFDM.hpp index 066862cfd..c888202d2 100644 --- a/src/FDM/YASim/FGFDM.hpp +++ b/src/FDM/YASim/FGFDM.hpp @@ -118,6 +118,7 @@ private: SGPropertyNode_ptr _arxN; SGPropertyNode_ptr _aryN; SGPropertyNode_ptr _arzN; + SGPropertyNode_ptr _cg_xmacN; }; }; // namespace yasim diff --git a/src/FDM/YASim/Wing.cpp b/src/FDM/YASim/Wing.cpp index 12a20728c..68489231a 100644 --- a/src/FDM/YASim/Wing.cpp +++ b/src/FDM/YASim/Wing.cpp @@ -212,6 +212,13 @@ void Wing::compile() // wingspan in y-direction (not for vstab) _wingspan = Math::abs(2*_tip[1]); _aspectRatio = _wingspan / _meanChord; + + _netSpan = Math::abs(2*(_tip[1]-_base[1])); + // http://www.nasascale.org/p2/wp-content/uploads/mac-calculator.htm + const float commonFactor = _chord*(0.5+_taper)/(3*_chord*(1+_taper)); + _mac = _chord-(2*_chord*(1-_taper)*commonFactor); + _macRootDistance = _netSpan*commonFactor; + _macX = _base[0]-Math::tan(_sweep) * _macRootDistance + _mac/2; // The wing's Y axis will be the "left" vector. The Z axis will // be perpendicular to this and the local (!) X axis, because we diff --git a/src/FDM/YASim/Wing.hpp b/src/FDM/YASim/Wing.hpp index d757f8236..9560b9a98 100644 --- a/src/FDM/YASim/Wing.hpp +++ b/src/FDM/YASim/Wing.hpp @@ -17,28 +17,28 @@ public: // Do we mirror ourselves about the XZ plane? void setMirror(bool mirror) { _mirror = mirror; } - bool isMirrored() { return _mirror; }; + const bool isMirrored() { return _mirror; }; // Wing geometry in local coordinates: // base point of wing void setBase(const float* base) { Math::set3(base, _base); } - void getBase(float* base) { Math::set3(_base, base); }; + const void getBase(float* base) { Math::set3(_base, base); }; // dist. ALONG wing (not span!) void setLength(float length) { _length = length; } - float getLength() { return _length; }; + const float getLength() { return _length; }; // at base, measured along X axis void setChord(float chord) { _chord = chord; } - float getChord() { return _chord; }; + const float getChord() { return _chord; }; // fraction of chord at wing tip, 0..1 void setTaper(float taper) { _taper = taper; } - float getTaper() { return _taper; }; + const float getTaper() { return _taper; }; // radians void setSweep(float sweep) { _sweep = sweep; } - float getSweep() { return _sweep; }; + const float getSweep() { return _sweep; }; // radians, positive is "up" void setDihedral(float dihedral) { _dihedral = dihedral; } - float getDihedral() { return _dihedral; }; + const float getDihedral() { return _dihedral; }; void setIncidence(float incidence); void setTwist(float angle) { _twist = angle; } @@ -67,27 +67,31 @@ public: // Compile the thing into a bunch of Surface objects void compile(); - void getTip(float* tip) { Math::set3(_tip, tip);}; + const void getTip(float* tip) { Math::set3(_tip, tip);}; // valid only after Wing::compile() was called - float getSpan() { return _wingspan; }; - float getArea() { return _wingspan*_meanChord; }; - float getAspectRatio() { return _aspectRatio; }; - float getSMC() { return _meanChord; }; + const float getSpan() { return _wingspan; }; + const float getArea() { return _wingspan*_meanChord; }; + const float getAspectRatio() { return _aspectRatio; }; + const float getSMC() { return _meanChord; }; + const float getMAC() { return _mac; }; // get length of MAC + const float getMACx() { return _macX; }; // get x-coord of MAC leading edge + const float getMACy() { return _base[1]+_macRootDistance; }; // get y-coord of MAC leading edge - int numSurfaces() { return _surfs.size(); } + + const int numSurfaces() { return _surfs.size(); } Surface* getSurface(int n) { return ((SurfRec*)_surfs.get(n))->surface; } - float getSurfaceWeight(int n) { return ((SurfRec*)_surfs.get(n))->weight; } + const float getSurfaceWeight(int n) { return ((SurfRec*)_surfs.get(n))->weight; } // The overall drag coefficient for the wing as a whole. Units are // arbitrary. void setDragScale(float scale); - float getDragScale() { return _dragScale; } + const float getDragScale() { return _dragScale; } // The ratio of force along the Z (lift) direction of each wing // segment to that along the X (drag) direction. void setLiftRatio(float ratio); - float getLiftRatio() { return _liftRatio; } + const float getLiftRatio() { return _liftRatio; } private: void interp(const float* v1, const float* v2, const float frac, float* out); @@ -114,6 +118,10 @@ private: // calculated from above float _tip[3]; float _meanChord; // std. mean chord + float _mac; // mean aerodynamic chord length + float _macRootDistance; // y-distance of mac from root + float _macX; // x-coordinate of mac (leading edge) + float _netSpan; float _wingspan; float _aspectRatio; diff --git a/src/FDM/YASim/yasim-test.cpp b/src/FDM/YASim/yasim-test.cpp index 39cac4eba..9e40dc233 100644 --- a/src/FDM/YASim/yasim-test.cpp +++ b/src/FDM/YASim/yasim-test.cpp @@ -237,7 +237,9 @@ int main(int argc, char** argv) } } else { - printf("Solution results:"); + printf("==========================\n"); + printf("= YASim solution results =\n"); + printf("==========================\n"); float aoa = a->getCruiseAoA() * RAD2DEG; float tail = -1 * a->getTailIncidence() * RAD2DEG; float drag = 1000 * a->getDragCoefficient(); @@ -247,18 +249,26 @@ int main(int argc, char** argv) float SI_inertia[9]; a->getModel()->getBody()->getInertiaMatrix(SI_inertia); + float MAC = a->getWing()->getMAC(); + float MACx = a->getWing()->getMACx(); + float MACy = a->getWing()->getMACy(); printf(" Iterations: %d\n", a->getSolutionIterations()); printf(" Drag Coefficient: %f\n", drag); printf(" Lift Ratio: %f\n", a->getLiftRatio()); - printf(" Cruise AoA: %f\n", aoa); - printf(" Tail Incidence: %f\n", tail); - printf("Approach Elevator: %f\n", a->getApproachElevator()); - printf(" CG: x:%.3f, y:%.3f, z:%.3f\n\n", cg[0], cg[1], cg[2]); - printf("Inertia tensor [kg*m^2], origo at CG:\n"); + printf(" Cruise AoA: %f deg\n", aoa); + printf(" Tail Incidence: %f deg\n", tail); + printf("Approach Elevator: %f\n\n", a->getApproachElevator()); + printf(" CG: x:%.3f, y:%.3f, z:%.3f\n", cg[0], cg[1], cg[2]); + printf(" Wing MAC (*1): x:%.2f, y:%.2f, length:%.1f \n", MACx, MACy, MAC); + printf(" CG-x rel. MAC: %.3f\n", a->getCGMAC()); + printf(" CG-x desired: %.3f < %.3f < %.3f \n", a->getCGSoftLimitXMin(), cg[0], a->getCGSoftLimitXMax()); + + printf("\nInertia tensor [kg*m^2], origo at CG:\n\n"); printf(" %7.3f, %7.3f, %7.3f\n", SI_inertia[0], SI_inertia[1], SI_inertia[2]); printf(" %7.3f, %7.3f, %7.3f\n", SI_inertia[3], SI_inertia[4], SI_inertia[5]); printf(" %7.3f, %7.3f, %7.3f\n", SI_inertia[6], SI_inertia[7], SI_inertia[8]); + printf("\n(*1) MAC calculation works on only! Numbers will be wrong for segmented wings, e.g. +.\n"); } delete fdm; return 0;