diff --git a/src/FDM/AIWake/AIWakeGroup.cxx b/src/FDM/AIWake/AIWakeGroup.cxx new file mode 100644 index 000000000..98889012c --- /dev/null +++ b/src/FDM/AIWake/AIWakeGroup.cxx @@ -0,0 +1,121 @@ +// AIWakeGroup.hxx -- Group of AI wake meshes for the computation of the induced +// wake. +// +// Written by Bertrand Coconnier, started April 2017. +// +// Copyright (C) 2017 Bertrand Coconnier - bcoconni@users.sf.net +// +// This program is free software; you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation; either version 2 of the License, or (at your option) any later +// version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// this program; if not, write to the Free Software Foundation, Inc., 51 +// Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id$ + +#include +#include +#include +#include +#include +#include +#include +#ifdef FG_TESTLIB +#include +#include +#include "tests/fakeAIAircraft.hxx" +Globals g; +Globals *globals = &g; +#else +#include "Main/globals.hxx" +#include "AIModel/performancedata.hxx" +#include "AIModel/AIAircraft.hxx" +#endif +#include "FDM/AIWake/AIWakeGroup.hxx" + +AIWakeGroup::AIWakeGroup(void) +{ + SGPropertyNode* _props = globals->get_props(); + _density_slugft = _props->getNode("environment/density-slugft3", true); +} + +void AIWakeGroup::AddAI(FGAIAircraft* ai) +{ + int id = ai->getID(); + PerformanceData* perfData = ai->getPerformance(); + + if (_aiWakeData.find(id) == _aiWakeData.end()) { + double span = perfData->wingSpan(); + double chord = perfData->wingChord(); + _aiWakeData[id] = AIWakeData(new WakeMesh(span, chord)); + + SG_LOG(SG_FLIGHT, SG_DEV_ALERT, + "Created mesh for " << ai->_getName() << " ID: #" << id + // << "Type:" << ai->getAcType() << endl + // << "Name:" << ai->_getName() << endl + // << "Speed:" << ai->getSpeed() << endl + // << "Vertical speed:" << ai->getVerticalSpeed() << endl + // << "Altitude:" << ai->getAltitude() << endl + // << "Heading:" << ai->_getHeading() << endl + // << "Roll:" << ai->getRoll() << endl + // << "Pitch:" << ai->getPitch() << endl + // << "Wing span (ft):" << span << endl + // << "Wing chord (ft):" << chord << endl + // << "Weight (lbs):" << perfData->weight() << endl + ); + } + + AIWakeData& data = _aiWakeData[id]; + data.visited = true; + + data.position = ai->getCartPos() * SG_METER_TO_FEET; + SGGeoc geoc = SGGeoc::fromCart(data.position); + SGQuatd Te2l = SGQuatd::fromLonLatRad(geoc.getLongitudeRad(), + geoc.getLatitudeRad()); + data.Te2b = Te2l * SGQuatd::fromYawPitchRollDeg(ai->_getHeading(), + ai->getPitch(), 0.0); + + double hVel = ai->getSpeed()*SG_KT_TO_FPS; + double vVel = ai->getVerticalSpeed()/60; + double gamma = atan2(vVel, hVel); + double vel = sqrt(hVel*hVel + vVel*vVel); + double weight = perfData->weight(); + _aiWakeData[id].mesh->computeAoA(vel, _density_slugft->getDoubleValue(), + weight*cos(gamma)); +} + +SGVec3d AIWakeGroup::getInducedVelocityAt(const SGVec3d& pt) const +{ + SGVec3d vi(0.,0.,0.); + for (auto item : _aiWakeData) { + AIWakeData& data = item.second; + if (!data.visited) continue; + + SGVec3d at = data.Te2b.transform(pt - data.position); + vi += data.Te2b.backTransform(data.mesh->getInducedVelocityAt(at)); + } + return vi; +} + +void AIWakeGroup::gc(void) +{ + for (auto it=_aiWakeData.begin(); it != _aiWakeData.end(); ++it) { + if (!(*it).second.visited) { + SG_LOG(SG_FLIGHT, SG_DEV_ALERT, "Deleted mesh for aircraft #" + << (*it).first); + _aiWakeData.erase(it); + break; // erase has invalidated the iterator, other dead meshes (if + // any) will be erased at the next time steps. + } + } + + for (auto &item : _aiWakeData) item.second.visited = false; +} diff --git a/src/FDM/AIWake/AIWakeGroup.hxx b/src/FDM/AIWake/AIWakeGroup.hxx new file mode 100644 index 000000000..e4fe934e2 --- /dev/null +++ b/src/FDM/AIWake/AIWakeGroup.hxx @@ -0,0 +1,55 @@ +// AIWakeGroup.hxx -- Group of AI wake meshes for the computation of the induced +// wake. +// +// Written by Bertrand Coconnier, started April 2017. +// +// Copyright (C) 2017 Bertrand Coconnier - bcoconni@users.sf.net +// +// This program is free software; you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation; either version 2 of the License, or (at your option) any later +// version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// this program; if not, write to the Free Software Foundation, Inc., 51 +// Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id$ + +#ifndef _FG_AIWAKEGROUP_HXX +#define _FG_AIWAKEGROUP_HXX + +#include "FDM/AIWake/WakeMesh.hxx" + +class FGAIAircraft; + +class AIWakeGroup { +#ifdef FG_TESTLIB +public: +#endif + struct AIWakeData { + explicit AIWakeData(WakeMesh* m = nullptr) : mesh(m) {} + + SGVec3d position {SGVec3d::zeros()}; + SGQuatd Te2b {SGQuatd::unit()}; + bool visited {false}; + WakeMesh_ptr mesh; + }; + + std::map _aiWakeData; + SGPropertyNode_ptr _density_slugft; + +public: + AIWakeGroup(void); + void AddAI(FGAIAircraft* ai); + SGVec3d getInducedVelocityAt(const SGVec3d& pt) const; + // Garbage collection + void gc(void); +}; + +#endif diff --git a/src/FDM/AIWake/AeroElement.cxx b/src/FDM/AIWake/AeroElement.cxx new file mode 100644 index 000000000..8fe4d83cc --- /dev/null +++ b/src/FDM/AIWake/AeroElement.cxx @@ -0,0 +1,82 @@ +// AeroElement.cxx -- Mesh element for the computation of a wing wake. +// +// Written by Bertrand Coconnier, started March 2017. +// +// Copyright (C) 2017 Bertrand Coconnier - bcoconni@users.sf.net +// +// This program is free software; you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation; either version 2 of the License, or (at your option) any later +// version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// this program; if not, write to the Free Software Foundation, Inc., 51 +// Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id$ + +#include + +#include +#include +#include "AeroElement.hxx" + +AeroElement::AeroElement(const SGVec3d& n1, const SGVec3d& n2, + const SGVec3d& n3, const SGVec3d& n4) +{ + p1 = 0.75*n2 + 0.25*n1; + p2 = 0.75*n3 + 0.25*n4; + normal = normalize(cross(n3 - n1, n2 - n4)); + collocationPt = 0.375*(n1 + n4) + 0.125*(n2 + n3); +} + +SGVec3d AeroElement::vortexInducedVel(const SGVec3d& p, const SGVec3d& n1, + const SGVec3d& n2) const +{ + SGVec3d r0 = n2-n1, r1 = p-n1, r2 = p-n2; + SGVec3d v = cross(r1, r2); + double vSqrNorm = dot(v, v); + double r1SqrNorm = dot(r1, r1); + double r2SqrNorm = dot(r2, r2); + + if ((vSqrNorm < 1E-6) || (r1SqrNorm < 1E-6) || (r2SqrNorm < 1E-6)) + return SGVec3d::zeros(); + + r1 /= sqrt(r1SqrNorm); + r2 /= sqrt(r2SqrNorm); + v *= dot(r0, r1-r2) / (4.0*M_PI*vSqrNorm); + + return v; +} + +SGVec3d AeroElement::semiInfVortexInducedVel(const SGVec3d& point, + const SGVec3d& vEnd, + const SGVec3d& vDir) const +{ + SGVec3d r = point-vEnd; + double rSqrNorm = dot(r, r); + double denom = rSqrNorm - dot(r, vDir)*sqrt(rSqrNorm); + + if (fabs(denom) < 1E-6) + return SGVec3d::zeros(); + + SGVec3d v = cross(r, vDir); + v /= 4*M_PI*denom; + + return v; +} + +SGVec3d AeroElement::getInducedVelocity(const SGVec3d& p) const +{ + const SGVec3d w(-1.,0.,0.); + SGVec3d v = semiInfVortexInducedVel(p, p1, w); + v -= semiInfVortexInducedVel(p, p2, w); + v += vortexInducedVel(p, p1, p2); + + return v; +} diff --git a/src/FDM/AIWake/AeroElement.hxx b/src/FDM/AIWake/AeroElement.hxx new file mode 100644 index 000000000..b424639b3 --- /dev/null +++ b/src/FDM/AIWake/AeroElement.hxx @@ -0,0 +1,45 @@ +// AeroElement.hxx -- Mesh element for the computation of a wing wake. +// +// Written by Bertrand Coconnier, started March 2017. +// +// Copyright (C) 2017 Bertrand Coconnier - bcoconni@users.sf.net +// +// This program is free software; you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation; either version 2 of the License, or (at your option) any later +// version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// this program; if not, write to the Free Software Foundation, Inc., 51 +// Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id$ + +#ifndef _FG_AEROELEMENT_HXX +#define _FG_AEROELEMENT_HXX + +class AeroElement : public SGReferenced { +public: + AeroElement(const SGVec3d& n1, const SGVec3d& n2, const SGVec3d& n3, + const SGVec3d& n4); + const SGVec3d& getNormal(void) const { return normal; } + const SGVec3d& getCollocationPoint(void) const { return collocationPt; } + SGVec3d getBoundVortex(void) const { return p2 - p1; } + SGVec3d getBoundVortexMidPoint(void) const { return 0.5*(p1+p2); } + SGVec3d getInducedVelocity(const SGVec3d& p) const; +private: + SGVec3d vortexInducedVel(const SGVec3d& p, const SGVec3d& n1, + const SGVec3d& n2) const; + SGVec3d semiInfVortexInducedVel(const SGVec3d& point, const SGVec3d& vEnd, + const SGVec3d& vDir) const; + SGVec3d p1, p2, normal, collocationPt; +}; + +typedef SGSharedPtr AeroElement_ptr; + +#endif diff --git a/src/FDM/AIWake/AircraftMesh.cxx b/src/FDM/AIWake/AircraftMesh.cxx new file mode 100644 index 000000000..cfcae8013 --- /dev/null +++ b/src/FDM/AIWake/AircraftMesh.cxx @@ -0,0 +1,100 @@ +// AircraftMesh.cxx -- Mesh for the computation of the wake induced force and +// moment on an aircraft. + +// Written by Bertrand Coconnier, started March 2017. +// +// Copyright (C) 2017 Bertrand Coconnier - bcoconni@users.sf.net +// +// This program is free software; you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation; either version 2 of the License, or (at your option) any later +// version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// this program; if not, write to the Free Software Foundation, Inc., 51 +// Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id$ + +#include +#include + +#include +#include +#include +#include +#include +#include "AircraftMesh.hxx" +#include +#ifndef FG_TESTLIB +#include "AIModel/AIAircraft.hxx" +#else +#include "AIWakeGroup.hxx" +#include "fakeAIAircraft.hxx" +#endif +extern "C" { +#include "../LaRCsim/ls_matrix.h" +} + +AircraftMesh::AircraftMesh(double _span, double _chord) + : WakeMesh(_span, _chord) +{ + collPt.resize(nelm, SGVec3d::zeros()); + midPt.resize(nelm, SGVec3d::zeros()); +} + +void AircraftMesh::setPosition(const SGVec3d& _pos, const SGQuatd& orient) +{ + SGVec3d pos = _pos * SG_METER_TO_FEET; + SGGeoc geoc = SGGeoc::fromCart(pos); + SGQuatd Te2l = SGQuatd::fromLonLatRad(geoc.getLongitudeRad(), + geoc.getLatitudeRad()); + Te2b = Te2l * orient; + + for (int i=0; igetCollocationPoint(); + collPt[i] = pos + Te2b.backTransform(pt); + pt = elements[i]->getBoundVortexMidPoint(); + midPt[i] = pos + Te2b.backTransform(pt); + } +} + +SGVec3d AircraftMesh::GetForce(const AIWakeGroup& wg, const SGVec3d& vel, + double rho) +{ + std::vector rhs; + rhs.resize(nelm, 0.0); + + for (int i=0; igetNormal(), + Te2b.transform(wg.getInducedVelocityAt(collPt[i]))); + + for (int i=1; i<=nelm; ++i) { + Gamma[i][1] = 0.0; + for (int k=1; k<=nelm; ++k) + Gamma[i][1] += influenceMtx[i][k]*rhs[k-1]; + } + + SGVec3d f(0.,0.,0.); + moment = SGVec3d::zeros(); + + for (int i=0; igetBoundVortexMidPoint(); + SGVec3d v = Te2b.transform(wg.getInducedVelocityAt(midPt[i])); + v += getInducedVelocityAt(mp); + + // The minus sign before vel to transform the aircraft velocity from the + // body frame to wind frame. + SGVec3d Fi = rho*Gamma[i+1][1]*cross(v-vel, + elements[i]->getBoundVortex()); + f += Fi; + moment += cross(mp, Fi); + } + + return f; +} diff --git a/src/FDM/AIWake/AircraftMesh.hxx b/src/FDM/AIWake/AircraftMesh.hxx new file mode 100644 index 000000000..884b7cb32 --- /dev/null +++ b/src/FDM/AIWake/AircraftMesh.hxx @@ -0,0 +1,49 @@ +// AircraftMesh.hxx -- Mesh for the computation of the wake induced force and +// moment on an aircraft. +// +// Written by Bertrand Coconnier, started March 2017. +// +// Copyright (C) 2017 Bertrand Coconnier - bcoconni@users.sf.net +// +// This program is free software; you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation; either version 2 of the License, or (at your option) any later +// version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// this program; if not, write to the Free Software Foundation, Inc., 51 +// Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id$ + +#ifndef _FG_AEROMESH_HXX +#define _FG_AEROMESH_HXX + +#include "WakeMesh.hxx" + +class FGAIAircraft; +class AIWakeGroup; + +class AircraftMesh : public WakeMesh { +public: + AircraftMesh(double _span, double _chord); + void setPosition(const SGVec3d& pos, const SGQuatd& orient); + SGVec3d GetForce(const AIWakeGroup& wg, const SGVec3d& vel, double rho); + const SGVec3d& GetMoment(void) const { return moment; }; + +#ifndef FG_TESTLIB +private: +#endif + std::vector collPt, midPt; + SGQuatd Te2b; + SGVec3d moment; +}; + +typedef SGSharedPtr AircraftMesh_ptr; + +#endif diff --git a/src/FDM/AIWake/WakeMesh.cxx b/src/FDM/AIWake/WakeMesh.cxx new file mode 100644 index 000000000..57dd7c356 --- /dev/null +++ b/src/FDM/AIWake/WakeMesh.cxx @@ -0,0 +1,106 @@ +// WakeMesh.cxx -- Mesh for the computation of a wing wake. +// +// Written by Bertrand Coconnier, started March 2017. +// +// Copyright (C) 2017 Bertrand Coconnier - bcoconni@users.sf.net +// +// This program is free software; you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation; either version 2 of the License, or (at your option) any later +// version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// this program; if not, write to the Free Software Foundation, Inc., 51 +// Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id$ + +#include +#include + +#include +#include + +#include "WakeMesh.hxx" +extern "C" { +#include "../LaRCsim/ls_matrix.h" +} + +WakeMesh::WakeMesh(double _span, double _chord) + : nelm(10), span(_span), chord(_chord) +{ + double y1 = -0.5*span; + double ds = span / nelm; + + for(int i=0; igetNormal(); + SGVec3d collPt = elements[i]->getCollocationPoint(); + + for (int j=0; j < nelm; ++j) + influenceMtx[i+1][j+1] = dot(elements[j]->getInducedVelocity(collPt), + normal); + } + + // Compute the inverse matrix with the Gauss-Jordan algorithm + nr_gaussj(influenceMtx, nelm, 0, 0); +} + +WakeMesh::~WakeMesh() +{ + nr_free_matrix(influenceMtx, 1, nelm, 1, nelm); + nr_free_matrix(Gamma, 1, nelm, 1, 1); +} + +double WakeMesh::computeAoA(double vel, double rho, double weight) +{ + for (int i=1; i<=nelm; ++i) { + Gamma[i][1] = 0.0; + for (int k=1; k<=nelm; ++k) + Gamma[i][1] -= influenceMtx[i][k]; + Gamma[i][1] *= vel; + } + + // Compute the lift only. Velocities in the z direction are discarded + // because they only produce drag. This include the vertical component + // vel*sin(alpha) and the induced velocities on the bound vortex. + // This assumption is only valid for small angles. + SGVec3d f(0.,0.,0.); + SGVec3d v(-vel, 0.0, 0.0); + + for (int i=0; igetBoundVortex()); + + double sinAlpha = -weight/f[2]; + + for (int i=1; i<=nelm; ++i) + Gamma[i][1] *= sinAlpha; + + return asin(sinAlpha); +} + +SGVec3d WakeMesh::getInducedVelocityAt(const SGVec3d& at) const +{ + SGVec3d v(0., 0., 0.); + + for (int i=0; igetInducedVelocity(at); + + return v; +} diff --git a/src/FDM/AIWake/WakeMesh.hxx b/src/FDM/AIWake/WakeMesh.hxx new file mode 100644 index 000000000..6605d6ae0 --- /dev/null +++ b/src/FDM/AIWake/WakeMesh.hxx @@ -0,0 +1,46 @@ +// WakeMesh.hxx -- Mesh for the computation of a wing wake. +// +// Written by Bertrand Coconnier, started March 2017. +// +// Copyright (C) 2017 Bertrand Coconnier - bcoconni@users.sf.net +// +// This program is free software; you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free Software +// Foundation; either version 2 of the License, or (at your option) any later +// version. +// +// This program is distributed in the hope that it will be useful, but WITHOUT +// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +// details. +// +// You should have received a copy of the GNU General Public License along with +// this program; if not, write to the Free Software Foundation, Inc., 51 +// Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +// +// $Id$ + +#ifndef _FG_WAKEMESH_HXX +#define _FG_WAKEMESH_HXX + +#include "AeroElement.hxx" + +class WakeMesh : public SGReferenced { +public: + WakeMesh(double _span, double _chord); + virtual ~WakeMesh(); + double computeAoA(double vel, double rho, double weight); + SGVec3d getInducedVelocityAt(const SGVec3d& at) const; + +#ifndef FG_TESTLIB +protected: +#endif + int nelm; + double span, chord; + std::vector elements; + double **influenceMtx, **Gamma; +}; + +typedef SGSharedPtr WakeMesh_ptr; + +#endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 79c90055b..07e964c1f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -113,3 +113,20 @@ flightgear_test(test_flightplan test_flightplan.cxx) add_executable(test_ls_matrix test_ls_matrix.cxx ${CMAKE_SOURCE_DIR}/src/FDM/LaRCsim/ls_matrix.c) target_link_libraries(test_ls_matrix SimGearCore) add_test(test_ls_matrix ${EXECUTABLE_OUTPUT_PATH}/test_ls_matrix) + +add_executable(testAeroElement testAeroElement.cxx ${CMAKE_SOURCE_DIR}/src/FDM/AIWake/AeroElement.cxx) +target_link_libraries(testAeroElement SimGearCore) +add_test(testAeroElement ${EXECUTABLE_OUTPUT_PATH}/testAeroElement) + +add_executable(testAeroMesh testAeroMesh.cxx + ${CMAKE_SOURCE_DIR}/src/FDM/AIWake/AeroElement.cxx + ${CMAKE_SOURCE_DIR}/src/FDM/AIWake/AircraftMesh.cxx + ${CMAKE_SOURCE_DIR}/src/FDM/AIWake/WakeMesh.cxx + ${CMAKE_SOURCE_DIR}/src/FDM/AIWake/AIWakeGroup.cxx + ${CMAKE_SOURCE_DIR}/src/FDM/LaRCsim/ls_matrix.c + ) +set_target_properties (testAeroMesh PROPERTIES COMPILE_DEFINITIONS "FG_TESTLIB") +target_include_directories(testAeroMesh PRIVATE ${CMAKE_SOURCE_DIR}/tests + ${CMAKE_SOURCE_DIR}/src/FDM/JSBSim) +target_link_libraries(testAeroMesh SimGearCore JSBSim) +add_test(testAeroMesh ${EXECUTABLE_OUTPUT_PATH}/testAeroMesh) diff --git a/tests/fakeAIAircraft.hxx b/tests/fakeAIAircraft.hxx new file mode 100644 index 000000000..d4d53b5d7 --- /dev/null +++ b/tests/fakeAIAircraft.hxx @@ -0,0 +1,58 @@ +#ifndef FAKEFGAIAIRCRAFT +#define FAKEFGAIAIRCRAFT + +struct Globals { + Globals(void) { root.getNode("environment/density-slugft3", true); } + SGPropertyNode* get_props(void) { return &root; } + SGPropertyNode root; +}; + +struct PerformanceData { + double _span, _chord, _weight; + PerformanceData(double s, double c, double w) + : _span(s), _chord(c), _weight(w) {} + double wingSpan(void) const { return _span; } + double wingChord(void) const { return _chord; } + double weight(void) const { return _weight; } +}; + +class FGAIAircraft { +public: + explicit FGAIAircraft(int id) + : _id(id), _pos(SGVec3d::zeros()), _heading(0.0), _pitch(0.0), + _vel(0.0), _perf(0.0, 0.0, 0.0) {} + void setPosition(const SGVec3d& pos) { _pos = pos; } + void setOrientation(double heading, double pitch) + { + _heading = heading; + _pitch = pitch; + } + void setGeom(double s, double c, double w) + { + _perf._span = s; + _perf._chord = c; + _perf._weight = w; + } + void setVelocity(double v) { _vel = v; } + const SGVec3d& getCartPos(void) const { return _pos; } + double _getHeading(void) const { return _heading; } + double getPitch(void) const { return _pitch; } + int getID(void) const { return _id; } + PerformanceData* getPerformance(void) { return &_perf; } + std::string getAcType(void) const { return "The AC type"; } + std::string _getName(void) const { return "The AC name"; } + double getSpeed(void) const { return _vel * SG_FPS_TO_KT; } + double getVerticalSpeed(void) const { return 0.0; } + double getAltitude(void) const { return 3000.; } + double getRoll(void) const { return 0.0; } + +private: + int _id; + SGVec3d _pos; + double _heading, _pitch, _vel; + PerformanceData _perf; +}; + +extern Globals* globals; + +#endif diff --git a/tests/testAeroElement.cxx b/tests/testAeroElement.cxx new file mode 100644 index 000000000..96597a3cc --- /dev/null +++ b/tests/testAeroElement.cxx @@ -0,0 +1,145 @@ +#include +#include +#include +#include + +#include "FDM/AIWake/AeroElement.hxx" + +void testNormal() +{ + AeroElement_ptr el = new AeroElement(SGVec3d(-1., -0.5, 0.), + SGVec3d(0., -0.5, 0.), + SGVec3d(0., 0.5, 0.), + SGVec3d(-1., 0.5, 0.)); + SGVec3d n = el->getNormal(); + SG_CHECK_EQUAL_EP(n[0], 0.0); + SG_CHECK_EQUAL_EP(n[1], 0.0); + SG_CHECK_EQUAL_EP(n[2], -1.0); +} + +void testCollocationPoint() +{ + AeroElement_ptr el = new AeroElement(SGVec3d(-1., -0.5, 0.), + SGVec3d(0., -0.5, 0.), + SGVec3d(0., 0.5, 0.), + SGVec3d(-1., 0.5, 0.)); + SGVec3d cp = el->getCollocationPoint(); + SG_CHECK_EQUAL_EP(cp[0], -0.75); + SG_CHECK_EQUAL_EP(cp[1], 0.0); + SG_CHECK_EQUAL_EP(cp[2], 0.0); +} + +void testBoundVortexMidPoint() +{ + AeroElement_ptr el = new AeroElement(SGVec3d(-1., -0.5, 0.), + SGVec3d(0., -0.5, 0.), + SGVec3d(0., 0.5, 0.), + SGVec3d(-1., 0.5, 0.)); + SGVec3d mp = el->getBoundVortexMidPoint(); + SG_CHECK_EQUAL_EP(mp[0], -0.25); + SG_CHECK_EQUAL_EP(mp[1], 0.0); + SG_CHECK_EQUAL_EP(mp[2], 0.0); +} + +void testBoundVortex() +{ + AeroElement_ptr el = new AeroElement(SGVec3d(-1., -0.5, 0.), + SGVec3d(0., -0.5, 0.), + SGVec3d(0., 0.5, 0.), + SGVec3d(-1., 0.5, 0.)); + SGVec3d v = el->getBoundVortex(); + SG_CHECK_EQUAL_EP(v[0], 0.0); + SG_CHECK_EQUAL_EP(v[1], 1.0); + SG_CHECK_EQUAL_EP(v[2], 0.0); +} + +void testInducedVelocityOnBoundVortex() +{ + AeroElement_ptr el = new AeroElement(SGVec3d(-1., -0.5, 0.), + SGVec3d(0., -0.5, 0.), + SGVec3d(0., 0.5, 0.), + SGVec3d(-1., 0.5, 0.)); + SGVec3d mp = el->getBoundVortexMidPoint(); + SGVec3d v = el->getInducedVelocity(mp); + SG_CHECK_EQUAL_EP(v[0], 0.0); + SG_CHECK_EQUAL_EP(v[1], 0.0); + SG_CHECK_EQUAL_EP(v[2], 1.0/M_PI); +} + +void testInducedVelocityOnCollocationPoint() +{ + AeroElement_ptr el = new AeroElement(SGVec3d(-1., -0.5, 0.), + SGVec3d(0., -0.5, 0.), + SGVec3d(0., 0.5, 0.), + SGVec3d(-1., 0.5, 0.)); + SGVec3d cp = el->getCollocationPoint(); + SGVec3d v = el->getInducedVelocity(cp); + SG_CHECK_EQUAL_EP(v[0], 0.0); + SG_CHECK_EQUAL_EP(v[1], 0.0); + SG_CHECK_EQUAL_EP(v[2], (1.0+sqrt(2.0)/M_PI)); +} + +void testInducedVelocityAtFarField() +{ + AeroElement_ptr el = new AeroElement(SGVec3d(-1., -0.5, 0.), + SGVec3d(0., -0.5, 0.), + SGVec3d(0., 0.5, 0.), + SGVec3d(-1., 0.5, 0.)); + SGVec3d mp = el->getBoundVortexMidPoint(); + SGVec3d v = el->getInducedVelocity(mp+SGVec3d(-1000.,0.,0.)); + SG_CHECK_EQUAL_EP(v[0], 0.0); + SG_CHECK_EQUAL_EP(v[1], 0.0); + SG_CHECK_EQUAL_EP(v[2], 2.0/M_PI); +} + +void testInducedVelocityAbove() +{ + AeroElement_ptr el = new AeroElement(SGVec3d(-1., -0.5, 0.), + SGVec3d(0., -0.5, 0.), + SGVec3d(0., 0.5, 0.), + SGVec3d(-1., 0.5, 0.)); + SGVec3d mp = el->getBoundVortexMidPoint(); + SGVec3d v = el->getInducedVelocity(mp+SGVec3d(0.,0.,-0.5)); + SG_CHECK_EQUAL_EP(v[0], -1.0/(sqrt(2.0)*M_PI)); + SG_CHECK_EQUAL_EP(v[1], 0.0); + SG_CHECK_EQUAL_EP(v[2], 0.5/M_PI); +} + +void testInducedVelocityAboveWithOffset() +{ + AeroElement_ptr el = new AeroElement(SGVec3d(-1., -0.5, 0.), + SGVec3d(0., -0.5, 0.), + SGVec3d(0., 0.5, 0.), + SGVec3d(-1., 0.5, 0.)); + SGVec3d mp = el->getBoundVortexMidPoint(); + SGVec3d v = el->getInducedVelocity(mp+SGVec3d(0.0, 0.5, -1.0)); + SG_CHECK_EQUAL_EP(v[0], -1.0/(4.0*M_PI*sqrt(2.0))); + SG_CHECK_EQUAL_EP(v[1], -0.125/M_PI); + SG_CHECK_EQUAL_EP(v[2], 0.125/M_PI); +} + +void testInducedVelocityUpstream() +{ + AeroElement_ptr el = new AeroElement(SGVec3d(-1., -0.5, 0.), + SGVec3d(0., -0.5, 0.), + SGVec3d(0., 0.5, 0.), + SGVec3d(-1., 0.5, 0.)); + SGVec3d mp = el->getBoundVortexMidPoint(); + SGVec3d v = el->getInducedVelocity(mp+SGVec3d(0.5, 0.0, 0.0)); + SG_CHECK_EQUAL_EP(v[0], 0.0); + SG_CHECK_EQUAL_EP(v[1], 0.0); + SG_CHECK_EQUAL_EP(v[2], (1.0-sqrt(2.0))/M_PI); +} + +int main(int argc, char* argv[]) +{ + testNormal(); + testCollocationPoint(); + testBoundVortexMidPoint(); + testBoundVortex(); + testInducedVelocityOnBoundVortex(); + testInducedVelocityAtFarField(); + testInducedVelocityAbove(); + testInducedVelocityAboveWithOffset(); + testInducedVelocityUpstream(); +} diff --git a/tests/testAeroMesh.cxx b/tests/testAeroMesh.cxx new file mode 100644 index 000000000..9a445ae93 --- /dev/null +++ b/tests/testAeroMesh.cxx @@ -0,0 +1,170 @@ +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "fakeAIAircraft.hxx" +#include "FDM/AIWake/AircraftMesh.hxx" +#include "FDM/AIWake/AIWakeGroup.hxx" +extern "C" { +#include "src/FDM/LaRCsim/ls_matrix.h" +} + +#include "FDM/JSBSim/math/FGLocation.h" +#include "FDM/JSBSim/math/FGQuaternion.h" + +using namespace std; +using namespace JSBSim; + +double rho = 2.0E-3; + +void testLiftComputation() +{ + double b = 10.0; + double c = 2.0; + AircraftMesh_ptr mesh = new AircraftMesh(b, c); + SGGeod geodPos = SGGeod::fromDeg(0.0, 0.0); + SGVec3d pos; + double vel = 100.; + double weight = 50.; + + SGGeodesy::SGGeodToCart(geodPos, pos); + mesh->setPosition(pos, SGQuatd::unit()); + + FGAIAircraft ai(1); + ai.setPosition(pos); + ai.setGeom(b, c, weight); + ai.setVelocity(vel); + + AIWakeGroup wg; + wg.AddAI(&ai); + + SGVec3d force = mesh->GetForce(wg, SGVec3d(vel, 0., 0.), rho); + + SG_CHECK_EQUAL_EP(force[1], 0.0); + SG_CHECK_EQUAL_EP(force[2], -weight); + + SGVec3d moment = mesh->GetMoment(); + SG_CHECK_EQUAL_EP(moment[0], 0.0); + SG_CHECK_EQUAL_EP(moment[1], -0.5*weight); + SG_CHECK_EQUAL_EP(moment[2], 0.0); + + for (int i=1; i<= mesh->nelm; ++i) + SG_CHECK_EQUAL_EP(wg._aiWakeData[1].mesh->Gamma[i][1], + mesh->Gamma[i][1]); +} + +void testFourierLiftingLine() +{ + double b = 10.0; + double c = 2.0; + double vel = 100.; + double weight = 50.; + + WakeMesh_ptr mesh = new WakeMesh(b, c); + int N = mesh->nelm; + double **mtx = nr_matrix(1, N, 1, N); + double **coef = nr_matrix(1, N, 1, 1); + + mesh->computeAoA(vel, rho, weight); + + for (int m=1; m<=N; ++m) { + double vm = M_PI*m/(N+1); + double sm = sin(vm); + coef[m][1] = c*M_PI/(2*b); + for (int n=1; n<=N; ++n) + mtx[m][n] = sin(n*vm)*(1.0+coef[m][1]*n/sm); + } + + nr_gaussj(mtx, N, coef, 1); + + double S = b*c; + double AR = b*b/S; + double lift = 0.5*rho*S*vel*vel*coef[1][1]*M_PI*AR; + double sinAlpha = weight / lift; + lift *= sinAlpha; + + cout << "y, Lift (Fourier), Lift (VLM), Corrected lift (VLM)" << endl; + + for (int i=1; i<=N; ++i) { + double y = mesh->elements[i-1]->getBoundVortexMidPoint()[1]; + double theta = acos(2.0*y/b); + double gamma = 0.0; + for (int n=1; n<=N; ++n) + gamma += coef[n][1]*sin(n*theta); + + gamma *= 2.0*b*vel*sinAlpha; + + cout << y << ", " << gamma << ", " << mesh->Gamma[i][1] << ", " + << mesh->Gamma[i][1] / gamma - 1.0 << endl; + } + + nr_free_matrix(mtx, 1, N, 1, N); + nr_free_matrix(coef, 1, N, 1, 1); +} + +void testFrameTransformations() +{ + double b = 10.0; + double c = 2.0; + double yaw = 80. * SGD_DEGREES_TO_RADIANS; + double pitch = 5. * SGD_DEGREES_TO_RADIANS; + double roll = -10. * SGD_DEGREES_TO_RADIANS; + SGQuatd orient = SGQuatd::fromYawPitchRoll(yaw, pitch, roll); + SGGeod geodPos = SGGeod::fromDeg(45.0, 10.0); + SGVec3d pos; + + SGGeodesy::SGGeodToCart(geodPos, pos); + SGGeoc geoc = SGGeoc::fromCart(pos); + + FGLocation loc(geoc.getLongitudeRad(), + geoc.getLatitudeRad(), + geoc.getRadiusFt()); + + SG_CHECK_EQUAL_EP(pos[0] * SG_METER_TO_FEET, loc(1)); + SG_CHECK_EQUAL_EP(pos[1] * SG_METER_TO_FEET, loc(2)); + SG_CHECK_EQUAL_EP(pos[2] * SG_METER_TO_FEET, loc(3)); + + AircraftMesh_ptr mesh = new AircraftMesh(b, c); + mesh->setPosition(pos, orient); + + FGQuaternion qJ(roll, pitch, yaw); + FGMatrix33 Tb2l = qJ.GetTInv(); + FGColumnVector3 refPos = loc.GetTec2l() * loc; + + for (int i=0; i < 4; ++i) + SG_CHECK_EQUAL_EP(orient(i), qJ((i+1) % 4 + 1)); + + for (int i=0; i < mesh->nelm; ++i) { + SGVec3d pt = mesh->elements[i]->getBoundVortexMidPoint(); + FGColumnVector3 ptJ(pt[0], pt[1], pt[2]); + FGColumnVector3 p = loc.GetTl2ec() * (refPos + Tb2l * ptJ); + SG_CHECK_EQUAL_EP(mesh->midPt[i][0], p(1)); + SG_CHECK_EQUAL_EP(mesh->midPt[i][1], p(2)); + SG_CHECK_EQUAL_EP(mesh->midPt[i][2], p(3)); + + pt = mesh->elements[i]->getCollocationPoint(); + ptJ.InitMatrix(pt[0], pt[1], pt[2]); + p = loc.GetTl2ec() * (refPos + Tb2l * ptJ); + SG_CHECK_EQUAL_EP(mesh->collPt[i][0], p(1)); + SG_CHECK_EQUAL_EP(mesh->collPt[i][1], p(2)); + SG_CHECK_EQUAL_EP(mesh->collPt[i][2], p(3)); + } +} + +int main(int argc, char* argv[]) +{ + globals->get_props()->getNode("environment/density-slugft3", false) + ->setDoubleValue(rho); + + testFourierLiftingLine(); + testLiftComputation(); + testFrameTransformations(); +}