fly/flightgear
1
0
Fork 0
Code Activity

Added the core of the AI wake computations with its tests.

Browse source 
At the moment, this is dead code: only the tests are compiled. FG is still compiled without this code.

A new directory is created that contains all the numerical computations made to estimate the wake induced by AI aircrafts. This is based on the venerable Vortex Lattice Method (VLM) which was all the rage in the 60's Computational Fluid Dynamics (CFD).

Even though quite old, the method is relevant to compute aircrafts wake in real time since 3D Navier-Stokes (NS3D) is out of reach for real time computations even with modern multicore personal computers and their GPUs.
This commit is contained in:
Bertrand Coconnier 2017-06-10 18:27:19 +02:00
parent 8667300483
commit c1313f2ecb
12 changed files with 994 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

121
src/FDM/AIWake/AIWakeGroup.cxx Normal file
View file

@ -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 <math.h>
#include <vector>
#include <simgear/structure/SGSharedPtr.hxx>
#include <simgear/math/SGVec3.hxx>
#include <simgear/math/SGGeod.hxx>
#include <simgear/math/SGGeoc.hxx>
#include <simgear/math/SGQuat.hxx>
#ifdef FG_TESTLIB
#include <map>
#include <simgear/props/props.hxx>
#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;
}

55
src/FDM/AIWake/AIWakeGroup.hxx Normal file
View file

@ -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<int, AIWakeData> _aiWakeData;
SGPropertyNode_ptr _density_slugft;
public:
AIWakeGroup(void);
void AddAI(FGAIAircraft* ai);
SGVec3d getInducedVelocityAt(const SGVec3d& pt) const;
// Garbage collection
void gc(void);
};
#endif

82
src/FDM/AIWake/AeroElement.cxx Normal file
View file

@ -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 <cmath>
#include <simgear/structure/SGSharedPtr.hxx>
#include <simgear/math/SGVec3.hxx>
#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;
}

45
src/FDM/AIWake/AeroElement.hxx Normal file
View file

@ -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> AeroElement_ptr;
#endif

100
src/FDM/AIWake/AircraftMesh.cxx Normal file
View file

@ -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 <vector>
#include <cmath>
#include <simgear/structure/SGSharedPtr.hxx>
#include <simgear/math/SGVec3.hxx>
#include <simgear/math/SGGeoc.hxx>
#include <simgear/math/SGGeod.hxx>
#include <simgear/math/SGQuat.hxx>
#include "AircraftMesh.hxx"
#include <FDM/flight.hxx>
#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; i<nelm; ++i) {
SGVec3d pt = elements[i]->getCollocationPoint();
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<double> rhs;
rhs.resize(nelm, 0.0);
for (int i=0; i<nelm; ++i)
rhs[i] = dot(elements[i]->getNormal(),
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; i<nelm; ++i) {
SGVec3d mp = elements[i]->getBoundVortexMidPoint();
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;
}

49
src/FDM/AIWake/AircraftMesh.hxx Normal file
View file

@ -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<SGVec3d> collPt, midPt;
SGQuatd Te2b;
SGVec3d moment;
};
typedef SGSharedPtr<AircraftMesh> AircraftMesh_ptr;
#endif

106
src/FDM/AIWake/WakeMesh.cxx Normal file
View file

@ -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 <vector>
#include <cmath>
#include <simgear/structure/SGSharedPtr.hxx>
#include <simgear/math/SGVec3.hxx>
#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; i<nelm; i++) {
double y2 = y1 + ds;
elements.push_back(new AeroElement(SGVec3d(-chord, y1, 0.),
SGVec3d(0., y1, 0.),
SGVec3d(0., y2, 0.),
SGVec3d(-chord, y2, 0.)));
y1 = y2;
}
influenceMtx = nr_matrix(1, nelm, 1, nelm);
Gamma = nr_matrix(1, nelm, 1, 1);
for (int i=0; i < nelm; ++i) {
SGVec3d normal = elements[i]->getNormal();
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; i<nelm; ++i)
f += rho*Gamma[i+1][1]*cross(v, elements[i]->getBoundVortex());
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; i<nelm; ++i)
v += Gamma[i+1][1] * elements[i]->getInducedVelocity(at);
return v;
}

46
src/FDM/AIWake/WakeMesh.hxx Normal file
View file

@ -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<AeroElement_ptr> elements;
double **influenceMtx, **Gamma;
};
typedef SGSharedPtr<WakeMesh> WakeMesh_ptr;
#endif

17
tests/CMakeLists.txt
View file

@ -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)

58
tests/fakeAIAircraft.hxx Normal file
View file

@ -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

145
tests/testAeroElement.cxx Normal file
View file

@ -0,0 +1,145 @@
#include <simgear/constants.h>
#include <simgear/misc/test_macros.hxx>
#include <simgear/structure/SGSharedPtr.hxx>
#include <simgear/math/SGVec3.hxx>
#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();
}

170
tests/testAeroMesh.cxx Normal file
View file

@ -0,0 +1,170 @@
#include <vector>
#include <map>
#include <iostream>
#include <simgear/constants.h>
#include <simgear/misc/test_macros.hxx>
#include <simgear/math/SGVec3.hxx>
#include <simgear/math/SGGeod.hxx>
#include <simgear/math/SGQuat.hxx>
#include <simgear/math/SGGeoc.hxx>
#include <simgear/props/props.hxx>
#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();
}