From 50670d87e45db674d32ff95a1abaa25e299eae91 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Wed, 8 Mar 2023 11:39:51 +0100
Subject: [PATCH] Remove speed of sound. Add output. Refactor names.
---
sdpm/src/sdpmProblem.cpp | 3 +-
sdpm/src/sdpmProblem.h | 12 ++---
sdpm/src/sdpmSolver.cpp | 110 +++++++++++++++++++++++++++++----------
3 files changed, 88 insertions(+), 37 deletions(-)
diff --git a/sdpm/src/sdpmProblem.cpp b/sdpm/src/sdpmProblem.cpp
index 50fcb8f..fe7db83 100644
--- a/sdpm/src/sdpmProblem.cpp
+++ b/sdpm/src/sdpmProblem.cpp
@@ -22,7 +22,7 @@
using namespace sdpm;
Problem::Problem(Mesh &msh, size_t nmods, std::vector<double> const &fs) : _msh(msh),
- _aoa(0.), _aos(0.), _mach(0.), _beta(1.0), _sound(0.), _uinf(sdpmVector3d::Zero()),
+ _aoa(0.), _aos(0.), _mach(0.), _beta(1.0), _uinf(sdpmVector3d::Zero()),
_dDir(sdpmVector3d::Zero()), _sDir(sdpmVector3d::Zero()), _lDir(sdpmVector3d::Zero()),
_sRef(0.), _cRef(0.), _xRef(sdpmVector3d::Zero()), _ySym(false),
_nModes(nmods), _omega(fs),
@@ -54,7 +54,6 @@ void Problem::updateFreestream(double aoa, double aos, double mch, double ui)
// Compressiblity
_mach = mch;
_beta = std::sqrt(1 - _mach * _mach);
- _sound = ui / _mach;
// Velocity and directions
_uinf = ui * sdpmVector3d(cos(_aoa) * cos(_aos), sin(_aos), sin(_aoa) * cos(_aos));
_dDir = sdpmVector3d(cos(_aoa) * cos(_aos), sin(_aos), sin(_aoa) * cos(_aos));
diff --git a/sdpm/src/sdpmProblem.h b/sdpm/src/sdpmProblem.h
index bee2d6f..a6bcbaa 100644
--- a/sdpm/src/sdpmProblem.h
+++ b/sdpm/src/sdpmProblem.h
@@ -39,7 +39,6 @@ class SDPM_API Problem : public sdpmObject
sdpmDouble _aos; ///< Angle of sideslip
sdpmDouble _mach; ///< Mach number
sdpmDouble _beta; ///< Compressibility factor
- sdpmDouble _sound; ///< Speed of sound
sdpmVector3d _uinf; ///< Freestream velocity
sdpmVector3d _dDir; ///< Direction of drag force
sdpmVector3d _sDir; ///< Direction of side force
@@ -71,17 +70,16 @@ public:
sdpmDouble getAos() const { return _aos; }
sdpmDouble getMach() const { return _mach; }
sdpmDouble getCompressibility() const { return _beta; }
- sdpmDouble getSpeedOfSound() const { return _sound; }
+ sdpmVector3d const &getVelocity() const { return _uinf; }
+ sdpmVector3d const &getDragDir() const { return _dDir; }
+ sdpmVector3d const &getSideDir() const { return _sDir; }
+ sdpmVector3d const &getLiftDir() const { return _lDir; }
sdpmDouble getRefArea() const { return _sRef; }
sdpmDouble getRefLgth() const { return _cRef; }
sdpmVector3d const &getRefCtr() const { return _xRef; }
sdpmDouble getSym() const { return _ySym; }
- std::vector<double> const &getFreqs() const { return _omega; }
size_t getNumModes() const { return _nModes; }
- sdpmVector3d const &getFrstrmVel() const { return _uinf; }
- sdpmVector3d const &getDragDir() const { return _dDir; }
- sdpmVector3d const &getSideDir() const { return _sDir; }
- sdpmVector3d const &getLiftDir() const { return _lDir; }
+ std::vector<double> const &getFrequencies() const { return _omega; }
#endif
void addBody(Body &b);
void updateFreestream(double aoa, double aos, double mch, double ui);
diff --git a/sdpm/src/sdpmSolver.cpp b/sdpm/src/sdpmSolver.cpp
index e830ff1..5e0e1eb 100644
--- a/sdpm/src/sdpmSolver.cpp
+++ b/sdpm/src/sdpmSolver.cpp
@@ -52,7 +52,7 @@ Solver::Solver(Problem &pbl) : _pbl(pbl), _cl(0), _cd(0), _cs(0), _cm(0)
// Get sizes
size_t nm = _pbl.getNumModes();
- size_t nf = _pbl.getFreqs().size();
+ size_t nf = _pbl.getFrequencies().size();
size_t ne = _pbl.getElements().size();
// Resize AIC matrics
@@ -84,13 +84,13 @@ Solver::Solver(Problem &pbl) : _pbl(pbl), _cl(0), _cd(0), _cs(0), _cm(0)
<< "Angle of attack: " << _pbl.getAoa() * 180 / 3.14159 << " deg\n"
<< "Angle of sideslip: " << _pbl.getAos() * 180 / 3.14159 << " deg\n"
<< "Mach number: " << _pbl.getMach() << "\n"
- << "Velocity: " << _pbl.getFrstrmVel().norm() << "\n";
+ << "Velocity: " << _pbl.getVelocity().norm() << "\n";
std::cout << "--- Reference data ---\n"
<< "Surface area: " << _pbl.getRefArea() << "\n"
<< "Chord length: " << _pbl.getRefLgth() << "\n"
<< "Origin: (" << _pbl.getRefCtr().transpose() << ")\n";
std::cout << "--- Unsteady parameters ---\n"
- << "Number of frequencies: " << _pbl.getFreqs().size() << "\n"
+ << "Number of frequencies: " << _pbl.getFrequencies().size() << "\n"
<< "Number of motions: " << _pbl.getNumModes() << "\n"
<< std::endl;
}
@@ -107,9 +107,9 @@ void Solver::update()
{
// Update elements and get values
_pbl.updateElements();
- std::vector<double> omega = _pbl.getFreqs();
- sdpmDouble mach = _pbl.getMach();
- sdpmDouble sound = _pbl.getSpeedOfSound();
+ std::vector<double> omega = _pbl.getFrequencies();
+ sdpmDouble ui = _pbl.getVelocity().norm();
+ sdpmDouble mi = _pbl.getMach();
sdpmDouble beta = _pbl.getCompressibility();
sdpmComplex img(0, 1);
@@ -137,10 +137,10 @@ void Solver::update()
_B0(_rows.at(ei), _rows.at(ej)) = mutau[1];
for (size_t fk = 0; fk < omega.size(); ++fk)
{
- sdpmDouble om = omega[fk] / sound / beta;
+ sdpmDouble om = omega[fk] * mi / ui / beta;
sdpmVector3d r = Builder::computeDistVec(false, ei, ej);
- sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
- _A[fk](_rows.at(ei), _rows.at(ej)) = (mutau[0] * (1. + img * om * r.norm()) - mutau[1] * img * om * mach * nx) * exp;
+ sdpmComplex exp = std::exp(-img * om * (r.norm() - mi * r(0)));
+ _A[fk](_rows.at(ei), _rows.at(ej)) = (mutau[0] * (1. + img * om * r.norm()) - mutau[1] * img * om * mi * nx) * exp;
_B[fk](_rows.at(ei), _rows.at(ej)) = mutau[1] * exp;
}
}
@@ -163,9 +163,9 @@ void Solver::update()
_A0(_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= mutau[0];
for (size_t fk = 0; fk < omega.size(); ++fk)
{
- sdpmDouble om = omega[fk] / sound / beta;
+ sdpmDouble om = omega[fk] * mi / ui / beta;
sdpmVector3d r = Builder::computeDistVec(false, ei, ew);
- sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
+ sdpmComplex exp = std::exp(-img * om * (r.norm() - mi * r(0)));
sdpmComplex waktrm = mutau[0] * (1. + img * om * r.norm()) * exp * std::exp(-img * omega[fk] * lagMap.at(ew));
_A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).first)) += waktrm;
_A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= waktrm;
@@ -194,10 +194,10 @@ void Solver::update()
_B0(_rows.at(ei), _rows.at(ej)) -= mutau[1];
for (size_t fk = 0; fk < omega.size(); ++fk)
{
- sdpmDouble om = omega[fk] / sound / beta;
+ sdpmDouble om = omega[fk] * mi / ui / beta;
sdpmVector3d r = Builder::computeDistVec(true, ei, ej);
- sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
- _A[fk](_rows.at(ei), _rows.at(ej)) -= (mutau[0] * (1. + img * om * r.norm()) - mutau[1] * img * om * mach * nx) * exp;
+ sdpmComplex exp = std::exp(-img * om * (r.norm() - mi * r(0)));
+ _A[fk](_rows.at(ei), _rows.at(ej)) -= (mutau[0] * (1. + img * om * r.norm()) - mutau[1] * img * om * mi * nx) * exp;
_B[fk](_rows.at(ei), _rows.at(ej)) -= mutau[1] * exp;
}
}
@@ -220,9 +220,9 @@ void Solver::update()
_A0(_rows.at(ei), _rows.at(wwMap.at(ew).second)) += mutau[0];
for (size_t fk = 0; fk < omega.size(); ++fk)
{
- sdpmDouble om = omega[fk] / sound / beta;
+ sdpmDouble om = omega[fk] * mi / ui / beta;
sdpmVector3d r = Builder::computeDistVec(true, ei, ew);
- sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
+ sdpmComplex exp = std::exp(-img * om * (r.norm() - mi * r(0)));
sdpmComplex waktrm = mutau[0] * (1. + img * om * r.norm()) * exp * std::exp(-img * omega[fk] * lagMap.at(ew));
_A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).first)) -= waktrm;
_A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).second)) += waktrm;
@@ -242,7 +242,7 @@ void Solver::run()
{
std::cout << "Running steady SDPM solver... " << std::endl;
std::cout << "computing sources... " << std::flush;
- sdpmVector3d ui = _pbl.getFrstrmVel(); // freestream velocity
+ sdpmVector3d ui = _pbl.getVelocity(); // freestream velocity
sdpmVectorXd etau(_nP);
for (auto body : _pbl.getBodies())
for (auto e : body->getElements())
@@ -273,7 +273,7 @@ void Solver::runUnsteady()
for (auto body : _pbl.getBodies())
body->computeVelocity(imd);
// Solve
- for (size_t ifq = 0; ifq < _pbl.getFreqs().size(); ++ifq)
+ for (size_t ifq = 0; ifq < _pbl.getFrequencies().size(); ++ifq)
{
std::cout << "Mode #" << imd << ", frequency #" << ifq << std::endl;
std::cout << "computing sources... " << std::flush;
@@ -284,7 +284,7 @@ void Solver::runUnsteady()
std::vector<sdpmVector3cd> uc1 = body->getHarmonicVelocity(imd);
std::vector<Element *> elems = body->getElements();
for (size_t i = 0; i < elems.size(); ++i)
- etau(_rows[elems[i]]) = (uc0[i] + uc1[i] * sdpmComplex(0., 1.) * _pbl.getFreqs()[ifq]).conjugate().dot(elems[i]->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
+ etau(_rows[elems[i]]) = (uc0[i] + uc1[i] * sdpmComplex(0., 1.) * _pbl.getFrequencies()[ifq]).conjugate().dot(elems[i]->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
}
std::cout << "done." << std::endl;
@@ -306,15 +306,52 @@ void Solver::runUnsteady()
*/
void Solver::save(GmshExport const &writer, std::string const &suffix)
{
+ // Get data
+ size_t nE = _pbl.getElements().size();
+ size_t nM = _pbl.getNumModes();
+ std::vector<sdpmDouble> omega = _pbl.getFrequencies();
+ size_t nF = omega.size();
+ sdpmDouble toDeg = 180. / 3.14159;
+
// Write files
std::cout << "Saving files... " << std::endl;
// set up results
Results results;
+ // steady
results.scalarsAtElems["phi"] = &_phi;
results.scalarsAtElems["mu"] = &_mu;
results.scalarsAtElems["tau"] = &_tau;
results.vectorsAtElems["u"] = &_u;
results.scalarsAtElems["cp"] = &_cp;
+ // unsteady
+ std::vector<std::vector<std::vector<sdpmVector3d>>> u1_abs(nM, std::vector<std::vector<sdpmVector3d>>(nF, std::vector<sdpmVector3d>(nE)));
+ std::vector<std::vector<std::vector<sdpmVector3d>>> u1_arg(nM, std::vector<std::vector<sdpmVector3d>>(nF, std::vector<sdpmVector3d>(nE)));
+ std::vector<std::vector<std::vector<sdpmDouble>>> cp1_abs(nM, std::vector<std::vector<sdpmDouble>>(nF, std::vector<sdpmDouble>(nE)));
+ std::vector<std::vector<std::vector<sdpmDouble>>> cp1_arg(nM, std::vector<std::vector<sdpmDouble>>(nF, std::vector<sdpmDouble>(nE)));
+ for (size_t imd = 0; imd < nM; ++imd)
+ {
+ for (size_t ifq = 0; ifq < nF; ++ifq)
+ {
+ // set names
+ std::stringstream u1_abs_name, u1_arg_name, cp1_abs_name, cp1_arg_name;
+ u1_abs_name << "|u1|_" << imd << "_" << ifq;
+ u1_arg_name << "< u1_" << imd << "_" << ifq;
+ cp1_abs_name << "|cp1|_" << imd << "_" << ifq;
+ cp1_arg_name << "< cp1_" << imd << "_" << ifq;
+ // set values
+ for (size_t i = 0; i < nE; ++i)
+ {
+ u1_abs[imd][ifq][i] = _u1[imd][ifq][i].array().abs();
+ u1_arg[imd][ifq][i] = _u1[imd][ifq][i].array().arg() * toDeg;
+ cp1_abs[imd][ifq][i] = std::abs(_cp1[imd][ifq][i]);
+ cp1_arg[imd][ifq][i] = std::arg(_cp1[imd][ifq][i]) * toDeg;
+ }
+ results.vectorsAtElems[u1_abs_name.str()] = &u1_abs[imd][ifq];
+ results.vectorsAtElems[u1_arg_name.str()] = &u1_arg[imd][ifq];
+ results.scalarsAtElems[cp1_abs_name.str()] = &cp1_abs[imd][ifq];
+ results.scalarsAtElems[cp1_arg_name.str()] = &cp1_arg[imd][ifq];
+ }
+ }
// save results
writer.save(results, suffix);
// save aerodynamic loads breakdown
@@ -323,13 +360,29 @@ void Solver::save(GmshExport const &writer, std::string const &suffix)
outfile.open("aeroloads_breakdown.dat");
for (auto b : _pbl.getBodies())
{
+ // steady
outfile << "# Body - " << b->getName() << " (" << b->getElements().size() << " elements)" << std::endl;
outfile << std::fixed << std::setprecision(5)
<< "Cl = " << b->getLiftCoeff() << "\n"
<< "Cd = " << b->getDragCoeff() << "\n"
<< "Cs = " << b->getSideforceCoeff() << "\n"
- << "Cm = " << b->getMomentCoeff() << std::endl
- << std::endl;
+ << "Cm = " << b->getMomentCoeff() << "\n";
+ // unsteady
+ for (size_t imd = 0; imd < nM; ++imd)
+ {
+ for (size_t ifq = 0; ifq < nF; ++ifq)
+ {
+ outfile << "|Cl1| (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::abs(b->getUnsteadyLiftCoeff(imd, ifq)) << "\n"
+ << "< Cl1 (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::arg(b->getUnsteadyLiftCoeff(imd, ifq)) * toDeg << "\n"
+ << "|Cd1| (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::abs(b->getUnsteadyDragCoeff(imd, ifq)) << "\n"
+ << "< Cd1 (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::arg(b->getUnsteadyDragCoeff(imd, ifq)) * toDeg << "\n"
+ << "|Cs1| (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::abs(b->getUnsteadySideforceCoeff(imd, ifq)) << "\n"
+ << "< Cs1 (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::arg(b->getUnsteadySideforceCoeff(imd, ifq)) * toDeg << "\n"
+ << "|Cm1| (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::abs(b->getUnsteadyMomentCoeff(imd, ifq)) << "\n"
+ << "< Cm1 (mode #" << imd << ", omega=" << omega[ifq] << ") = " << std::arg(b->getUnsteadyMomentCoeff(imd, ifq)) * toDeg << "\n"
+ << std::endl;
+ }
+ }
}
outfile.close();
std::cout << "done." << std::endl;
@@ -366,8 +419,8 @@ void Solver::post(sdpmVectorXd const &etau, sdpmVectorXd const &emu)
}
// Compute surface flow functionals
- sdpmVector3d ui = _pbl.getFrstrmVel(); // freestream velocity
- sdpmDouble mi = _pbl.getMach(); // freestream speed of sound
+ sdpmVector3d ui = _pbl.getVelocity(); // freestream velocity
+ sdpmDouble mi = _pbl.getMach(); // freestream speed of sound
for (auto body : bodies)
{
for (auto e : body->getElements())
@@ -446,13 +499,14 @@ void Solver::postUnsteady(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdp
}
// Compute flow and loads
- sdpmVector3d ui = _pbl.getFrstrmVel();
- double omega = _pbl.getFreqs()[ifq];
+ sdpmVector3d ui = _pbl.getVelocity();
+ sdpmDouble mi = _pbl.getMach();
+ sdpmDouble omega = _pbl.getFrequencies()[ifq];
sdpmDouble sRef = _pbl.getRefArea();
sdpmDouble cRef = _pbl.getRefLgth();
sdpmVector3d xRef = _pbl.getRefCtr();
- double ca = cos(_pbl.getAos());
- double sa = sin(_pbl.getAos());
+ sdpmDouble ca = cos(_pbl.getAos());
+ sdpmDouble sa = sin(_pbl.getAos());
_cl1[imd][ifq] = sdpmComplex(0., 0.);
_cd1[imd][ifq] = sdpmComplex(0., 0.);
_cs1[imd][ifq] = sdpmComplex(0., 0.);
@@ -480,7 +534,7 @@ void Solver::postUnsteady(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdp
// ecp = -2 / u_inf^2 * si + 1/u_inf^2/a_inf^2 * conv(si, si)
sdpmVectorXcd ecp(9);
ecp << 0., 0., -2 * si / ui.dot(ui), 0., 0.;
- ecp += convolve(si, si) / (_pbl.getSpeedOfSound() * _pbl.getSpeedOfSound() * ui.dot(ui));
+ ecp += convolve(si, si) * mi * mi / (ui.dot(ui) * ui.dot(ui));
for (size_t u = 0; u < 3; ++u)
cp012[i][u] = ecp(4 + u);
_cp1[imd][ifq][elems[i]->getId() - 1] = cp012[i][1];
--
GitLab