From 67b15a9ce084924b088b92d4ca51163e5edc22bc Mon Sep 17 00:00:00 2001 From: acrovato <a.crovato@uliege.be> Date: Tue, 14 Jul 2020 13:54:53 +0200 Subject: [PATCH] Add Builder and Solver --- fpm/_src/fpmw.h | 2 + fpm/_src/fpmw.i | 10 +- fpm/src/CMakeLists.txt | 6 -- fpm/src/fBody.cpp | 17 ++- fpm/src/fBody.h | 21 ++-- fpm/src/fBuilder.cpp | 162 +++++++++++++++++++++++++++++ fpm/src/fBuilder.h | 71 +++++++++++++ fpm/src/fProblem.cpp | 62 ++++++++++- fpm/src/fProblem.h | 16 ++- fpm/src/fProblem.inl | 31 ++++++ fpm/src/fSolver.cpp | 228 +++++++++++++++++++++++++++++++++++++++++ fpm/src/fSolver.h | 69 +++++++++++++ fpm/src/fpm.h | 4 + fpm/tests/basic.py | 11 +- 14 files changed, 684 insertions(+), 26 deletions(-) create mode 100644 fpm/src/fBuilder.cpp create mode 100644 fpm/src/fBuilder.h create mode 100644 fpm/src/fProblem.inl create mode 100644 fpm/src/fSolver.cpp create mode 100644 fpm/src/fSolver.h diff --git a/fpm/_src/fpmw.h b/fpm/_src/fpmw.h index 59e1f6a..0045215 100644 --- a/fpm/_src/fpmw.h +++ b/fpm/_src/fpmw.h @@ -15,6 +15,8 @@ */ #include "fBody.h" +#include "fBuilder.h" #include "fProblem.h" +#include "fSolver.h" #include "fTimers.h" #include "fWake.h" diff --git a/fpm/_src/fpmw.i b/fpm/_src/fpmw.i index b68248d..c09a59a 100644 --- a/fpm/_src/fpmw.i +++ b/fpm/_src/fpmw.i @@ -74,5 +74,13 @@ def __getitem__(self, name): %include "fWake.h" %shared_ptr(fpm::Problem); -%immutable fpm::Problem::msh; // avoids the creation of the setter method +%immutable fpm::Problem::msh; %include "fProblem.h" + +%shared_ptr(fpm::Builder); +%immutable fpm::Builder::pbl; +%include "fBuilder.h" + +%shared_ptr(fpm::Solver); +%immutable fpm::Solver::aic; +%include "fSolver.h" diff --git a/fpm/src/CMakeLists.txt b/fpm/src/CMakeLists.txt index 54ef33c..c583cfd 100644 --- a/fpm/src/CMakeLists.txt +++ b/fpm/src/CMakeLists.txt @@ -21,12 +21,6 @@ MACRO_DebugPostfix(fpm) TARGET_INCLUDE_DIRECTORIES(fpm PUBLIC ${PROJECT_SOURCE_DIR}/fpm/src) -# -- TBB -FIND_PACKAGE(TBB REQUIRED) -SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${TBB_CXX_FLAGS_DEBUG}") -TARGET_INCLUDE_DIRECTORIES(fpm PUBLIC ${TBB_INCLUDE_DIRS}) -TARGET_LINK_LIBRARIES(fpm ${TBB_LIBRARIES}) - # -- Eigen FIND_PACKAGE(EIGEN 3.3.4 REQUIRED) TARGET_INCLUDE_DIRECTORIES(fpm PUBLIC ${EIGEN_INCLUDE_DIRS}) diff --git a/fpm/src/fBody.cpp b/fpm/src/fBody.cpp index 10e9200..571fb58 100644 --- a/fpm/src/fBody.cpp +++ b/fpm/src/fBody.cpp @@ -30,6 +30,17 @@ using namespace fpm; Body::Body(std::shared_ptr<MshData> _msh, std::string const &id, std::vector<std::string> const &teIds, double xF) : Group(_msh, id), Cl(0), Cd(0), Cs(0), Cm(0) { + // Get nodes + for (auto e : tag->elems) + for (auto n : e->nodes) + nodes.push_back(n); + std::sort(nodes.begin(), nodes.end()); + nodes.erase(std::unique(nodes.begin(), nodes.end()), nodes.end()); + // Associate each node to its adjacent elements + for (auto e : tag->elems) + for (auto n : e->nodes) + neMap[n].push_back(e); + // If wakes are requested, check if they are already available, otherwise create them bool hasChanged = false; size_t j = 0; @@ -94,9 +105,9 @@ Body::Body(std::shared_ptr<MshData> _msh, std::string const &id, } // Size load vectors - cLoadX.resize(tag->elems.size()); - cLoadY.resize(tag->elems.size()); - cLoadZ.resize(tag->elems.size()); + cLoadX.resize(nodes.size()); + cLoadY.resize(nodes.size()); + cLoadZ.resize(nodes.size()); } /** diff --git a/fpm/src/fBody.h b/fpm/src/fBody.h index 7cdea5a..181af38 100644 --- a/fpm/src/fBody.h +++ b/fpm/src/fBody.h @@ -20,25 +20,28 @@ #include "fpm.h" #include "wGroup.h" #include <vector> +#include <map> namespace fpm { /** - * @brief Manage a body immersed in the fluid + * @brief Manage a body immersed in the fluid (lifting surface) * @authors Adrien Crovato */ class FPM_API Body : public tbox::Group { public: - double Cl; ///< lift coefficient - double Cd; ///< drag coefficient - double Cs; ///< sideforce coefficient - double Cm; ///< pitch moment coefficient (positive nose-up) - std::vector<double> cLoadX; ///< x-component of aerodynamic load coefficient on boundary - std::vector<double> cLoadY; ///< y-component of aerodynamic load coefficient on boundary - std::vector<double> cLoadZ; ///< z-component of aerodynamic load coefficient on boundary - std::vector<Wake *> wakes; ///< wake(s) attached to the lifting surface + double Cl; ///< lift coefficient + double Cd; ///< drag coefficient + double Cs; ///< sideforce coefficient + double Cm; ///< pitch moment coefficient (positive nose-up) + std::vector<double> cLoadX; ///< x-component of aerodynamic load normalized by dynamic pressure + std::vector<double> cLoadY; ///< y-component of aerodynamic load normalized by dynamic pressure + std::vector<double> cLoadZ; ///< z-component of aerodynamic load normalized by dynamic pressure + std::vector<tbox::Node *> nodes; ///< nodes of the surface + std::map<tbox::Node *, std::vector<tbox::Element *>> neMap; ///< map between nodes and adjacent elements + std::vector<Wake *> wakes; ///< wake(s) attached to the lifting surface Body(std::shared_ptr<tbox::MshData> _msh, std::string const &name, std::vector<std::string> const &teNames, double xF); ~Body() { std::cout << "~Body()\n"; } diff --git a/fpm/src/fBuilder.cpp b/fpm/src/fBuilder.cpp new file mode 100644 index 0000000..fc7bfdd --- /dev/null +++ b/fpm/src/fBuilder.cpp @@ -0,0 +1,162 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "fBuilder.h" +#include "fProblem.h" +//#include "fField.h" +#include "fBody.h" +#include "fWake.h" + +#include "wTag.h" +#include "wElement.h" +#include "wMem.h" + +#define PI 3.14159 + +using namespace tbox; +using namespace fpm; + +/** + * @brief Initialize the solver and perform sanity checks + * @authors Adrien Crovato + */ +Builder::Builder(std::shared_ptr<Problem> _pbl) : pbl(_pbl), valid(false) +{ + // Check the problem and update element memory + pbl->check(); + pbl->updateMem(); + + // Setup AIC matrices + nF = 0; //pbl->field->tag->elems.size(); + nP = 0; + for (auto body : pbl->bodies) + nP += body->tag->elems.size(); + A.resize(nP, nP); + B.resize(nP, nP); + C.resize(nP, nF); + Cx.resize(nP, nF); + Cy.resize(nP, nF); + Cz.resize(nP, nF); +} + +/** + * @brief Fill AIC matrices + */ +void Builder::run() +{ + // global->local id + std::map<Element *, int> rows; + int i = 0; + for (auto body : pbl->bodies) + for (auto e : body->tag->elems) + { + rows[e] = i; + i++; + } + + // body to body + for (auto ibody : pbl->bodies) + for (auto ei : ibody->tag->elems) + { + // body panels + for (auto jbody : pbl->bodies) + for (auto ej : jbody->tag->elems) + { + A(rows.at(ei), rows.at(ej)) = mu(ei, ej); + B(rows.at(ei), rows.at(ej)) = tau(ei, ej); + } + // wake panels + for (auto jbody : pbl->bodies) + for (auto wake : jbody->wakes) + for (auto ew : wake->tag->elems) + { + double wmu = mu(ei, ew); + A(rows.at(ei), rows.at(wake->wkMap.at(ew).first)) += wmu; + A(rows.at(ei), rows.at(wake->wkMap.at(ew).second)) -= wmu; + } + } + + // field to body + //for (auto body : pbl->bodies) + // for (auto e : body->tag->elems) + // for (size_t j = 0; j < field->tag->elems.size(); ++j) + // { + // C(rows.at(ei), j) = sigma(ei, field->tag->elems[j]); + // Cx(rows.at(ei), j) = sigmaX(ei, field->tag->elems[j]); + // Cy(rows.at(ei), j) = sigmaY(ei, field->tag->elems[j]); + // Cz(rows.at(ei), j) = sigmaY(ei, field->tag->elems[j]); + // } + + valid = true; +} + +/** + * @brief Compute doublet AIC from body panel ej to body panel ei + */ +double Builder::mu(Element *ei, Element *ej) +{ + if (ei == ej) + return 0.5; + else + return 0; +} +/** + * @brief Compute source AIC from body panel ej to body panel ei + */ +double Builder::tau(Element *ei, Element *ej) +{ + if (ei == ej) + return 0; + else + return 0; +} +/** + * @brief Compute source AIC from field panel ej to body panel ei + */ +double Builder::sigma(Element *ei, Element *ej) +{ + return 0; +} +/** + * @brief Compute source x-velocity AIC from field panel ej to body panel ei + */ +double Builder::sigmaX(Element *ei, Element *ej) +{ + return 0; +} +/** + * @brief Compute source y-velocity AIC from field panel ej to body panel ei + */ +double Builder::sigmaY(Element *ei, Element *ej) +{ + return 0; +} +/** + * @brief Compute source z-velocity AIC from field panel ej to body panel ei + */ +double Builder::sigmaZ(Element *ei, Element *ej) +{ + return 0; +} + +void Builder::write(std::ostream &out) const +{ + out << "fpm::Builder with AIC matrices:" + << "\n\tA(" << A.rows() << "," << A.cols() << "), B(" << B.rows() << "," << B.cols() << ")" + << "\n\tC(" << C.rows() << "," << C.cols() << ")" + << "\n\tCx(" << Cx.rows() << "," << Cx.cols() << "), Cy(" << Cy.rows() << "," << Cy.cols() << "), Cz(" << Cz.rows() << "," << Cz.cols() << ")" + << "\nwith\t" << *pbl; +} diff --git a/fpm/src/fBuilder.h b/fpm/src/fBuilder.h new file mode 100644 index 0000000..bb290cc --- /dev/null +++ b/fpm/src/fBuilder.h @@ -0,0 +1,71 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#ifndef FBUILDER_H +#define FBUILDER_H + +#include "fpm.h" +#include "wObject.h" + +#include <iostream> +#include <vector> +#include <memory> +#include <Eigen/Dense> + +namespace fpm +{ + +/** + * @brief Aerodynamic Influence Coefficients builder + * @authors Adrien Crovato + * @todo add symmetry support + * @todo implement AIC computation + */ +class FPM_API Builder : public fwk::wSharedObject +{ +public: + std::shared_ptr<Problem> pbl; ///< problem definition + + bool valid; ///< matrices are up-to-date + int nP; ///< number of body panels + int nF; ///< number of field panels + Eigen::MatrixXd A; ///< body-body doublet matrix + Eigen::MatrixXd B; ///< body-body source matrix + Eigen::MatrixXd C; ///< field-body source matrix + Eigen::MatrixXd Cx; ///< field-body x-velocity source matrix + Eigen::MatrixXd Cy; ///< field-body y-velocity source matrix + Eigen::MatrixXd Cz; ///< field-body z-velocity source matrix + + Builder(std::shared_ptr<Problem> _pbl); + ~Builder() { std::cout << "~Builder()\n"; } + + void run(); + +#ifndef SWIG + virtual void write(std::ostream &out) const override; +#endif + +private: + double mu(tbox::Element *ei, tbox::Element *ej); + double tau(tbox::Element *ei, tbox::Element *ej); + double sigma(tbox::Element *ei, tbox::Element *ej); + double sigmaX(tbox::Element *ei, tbox::Element *ej); + double sigmaY(tbox::Element *ei, tbox::Element *ej); + double sigmaZ(tbox::Element *ei, tbox::Element *ej); +}; + +} // namespace fpm +#endif //FBUILDER_H diff --git a/fpm/src/fProblem.cpp b/fpm/src/fProblem.cpp index 5ac5429..c2b666e 100644 --- a/fpm/src/fProblem.cpp +++ b/fpm/src/fProblem.cpp @@ -19,12 +19,12 @@ #include "wTag.h" #include "wElement.h" #include "wMem.h" +#include "wCache.h" using namespace tbox; using namespace fpm; -Problem::Problem(std::shared_ptr<MshData> _smsh, double aoa, double aos, double sref, double cref, - double xref, double yref, double zref) : smsh(_smsh), alpha(aoa), beta(aos), - sR(sref), cR(cref) +Problem::Problem(std::shared_ptr<MshData> _msh, double aoa, double aos, double mch, + double sref, double cref, double xref, double yref, double zref) : msh(_msh), alpha(aoa), beta(aos), mach(mch), sR(sref), cR(cref) { xR(0) = xref; xR(1) = yref; @@ -39,6 +39,62 @@ void Problem::add(std::shared_ptr<Body> b) bodies.push_back(b); } +/** + * @brief Compute velocity at element + */ +Eigen::Vector3d Problem::U(tbox::Element &e, std::vector<double> const &phi) +{ + Eigen::Vector3d v(0, 0, 0); + // @todo to be implemented if usefull + return v; + std::cout << std::endl; +} + +/** + * @brief Compute density at element + */ +double Problem::Rho(tbox::Element &e, std::vector<double> const &phi) +{ + // particularized: pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1)); + Eigen::Vector3d u = U(e, phi); + double a = 1 + 0.2 * mach * mach * (1 - u.dot(u)); + return sqrt(a * a * a * a * a); +} + +/** + * @brief Compute Mach number at element + */ +double Problem::Mach(tbox::Element &e, std::vector<double> const &phi) +{ + if (mach == 0) + return 0; + else + { + // particularized: pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1)); + Eigen::Vector3d u = U(e, phi); + double a2 = 1 / (mach * mach) + 0.2 * (1 - u.dot(u)); // speed of sound squared + return u.norm() / sqrt(a2); + } +} + +/** + * @brief Compute pressure coefficient at element + */ +double Problem::Cp(tbox::Element &e, std::vector<double> const &phi) +{ + Eigen::Vector3d u = U(e, phi); + if (mach == 0) + { + return 1 - u.dot(u); + } + else + { + //particularized: 2 / (gamma * mInf * mInf) * (pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), gamma / (gamma - 1)) - 1); + double a = 1 + 0.2 * mach * mach * (1 - u.dot(u)); + return 2 / (1.4 * mach * mach) * (sqrt(a * a * a * a * a * a * a) - 1); + } +} + /** * @brief Update the elements memory (Jacobian) */ diff --git a/fpm/src/fProblem.h b/fpm/src/fProblem.h index 54cf015..7c87145 100644 --- a/fpm/src/fProblem.h +++ b/fpm/src/fProblem.h @@ -28,26 +28,36 @@ namespace fpm /** * @brief Manage the problem * @authors Adrien Crovato + * @todo check and implement velocity computation */ class FPM_API Problem : public fwk::wSharedObject { public: - std::shared_ptr<tbox::MshData> smsh; ///< surface mesh data structure + std::shared_ptr<tbox::MshData> msh; ///< mesh data structure #ifndef SWIG std::vector<std::shared_ptr<Body>> bodies; ///< bodies in fluid #endif // Reference values double alpha; ///< Angle of attack double beta; ///< Angle of sideslip + double mach; ///< Mach number double sR; ///< Reference surface double cR; ///< Reference chord Eigen::Vector3d xR; ///< Reference center point (for moment computation) - Problem(std::shared_ptr<tbox::MshData> _smsh, double aoa, double aos, double sref, double cref, double xref, double yref, double zref); + Problem(std::shared_ptr<tbox::MshData> _msh, double aoa, double aos, double mch, double sref, double cref, double xref, double yref, double zref); ~Problem() { std::cout << "~Problem()\n"; } void add(std::shared_ptr<Body> b); + // Functions + inline double PhiI(Eigen::Vector3d const &pos); + inline Eigen::Vector3d Ui(); + Eigen::Vector3d U(tbox::Element &e, std::vector<double> const &phi); + double Rho(tbox::Element &e, std::vector<double> const &phi); + double Mach(tbox::Element &e, std::vector<double> const &phi); + double Cp(tbox::Element &e, std::vector<double> const &phi); + #ifndef SWIG void check() const; void updateMem(); @@ -55,6 +65,8 @@ public: #endif }; +#include "fProblem.inl" + } // namespace fpm #endif //FPROBLEM_H diff --git a/fpm/src/fProblem.inl b/fpm/src/fProblem.inl new file mode 100644 index 0000000..175e1a3 --- /dev/null +++ b/fpm/src/fProblem.inl @@ -0,0 +1,31 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +/** + * @brief Compute freestream potential + */ +inline double Problem::PhiI(Eigen::Vector3d const &pos) +{ + return cos(alpha) * cos(beta) * pos(0) + sin(beta) * pos(1) + sin(alpha) * cos(beta) * pos(2); +} + +/** + * @brief Compute freestream velocity + */ +inline Eigen::Vector3d Problem::Ui() +{ + return Eigen::Vector3d(cos(alpha) * cos(beta), sin(beta), sin(alpha) * cos(beta)); +} diff --git a/fpm/src/fSolver.cpp b/fpm/src/fSolver.cpp new file mode 100644 index 0000000..1769580 --- /dev/null +++ b/fpm/src/fSolver.cpp @@ -0,0 +1,228 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "fSolver.h" +#include "fBuilder.h" +#include "fProblem.h" +#include "fBody.h" + +#include "wMshData.h" +#include "wNode.h" +#include "wElement.h" +#include "wMem.h" +#include "wTag.h" +#include "wResults.h" +#include "wMshExport.h" +#include "wCache.h" + +using namespace tbox; +using namespace fpm; + +#define ANSI_COLOR_YELLOW "\x1b[1;33m" +#define ANSI_COLOR_RESET "\x1b[0m" + +/** + * @brief Initialize the solver and perform sanity checks + * @authors Adrien Crovato + */ +Solver::Solver(std::shared_ptr<Builder> _aic) : aic(_aic), Cl(0), Cd(0), Cs(0), Cm(0) +{ + // Say hi + std::cout << std::endl; + std::cout << "*******************************************************************************" << std::endl; + std::cout << "** \\_/ **" << std::endl; + std::cout << "** \\_______O(_)O_______/ **" << std::endl; + std::cout << "** **" << std::endl; + std::cout << "*******************************************************************************" << std::endl; + std::cout << "** Hi! My name is FPM v0.1-2007 **" << std::endl; + std::cout << "** Adrien Crovato & Romain Boman **" << std::endl; + std::cout << "** ULiege 2020-2021 **" << std::endl; + std::cout << "*******************************************************************************" << std::endl + << std::endl; + + // Setup variables (unkwown vector) + phi.resize(aic->pbl->msh->nodes.size(), 0.); + vPhi.resize(aic->pbl->msh->nodes.size(), 0.); + U.resize(aic->pbl->msh->elems.size(), Eigen::Vector3d::Zero()); + rho.resize(aic->pbl->msh->elems.size(), 0.); + M.resize(aic->pbl->msh->elems.size(), 0.); + Cp.resize(aic->pbl->msh->elems.size(), 0.); +} + +/** + * @brief Solve the (full) potential equation + */ +void Solver::run() +{ + if (!aic->valid) + throw std::runtime_error("fpm::Solver:run() AIC matrices are not up-to-date!\n"); + + // @todo compute sigma from an interpolator and compute volume flow + std::cout << "Computing field sources... " << std::flush; + Eigen::VectorXd sigma = Eigen::VectorXd::Zero(aic->nF); + Eigen::VectorXd uX = aic->Cx * sigma; + Eigen::VectorXd uY = aic->Cy * sigma; + Eigen::VectorXd uZ = aic->Cz * sigma; + std::cout << "done!" << std::endl; + + std::cout << "Computing body sources... " << std::flush; + Eigen::VectorXd mu(aic->nP), tau(aic->nP); + int ig = 0; + for (auto body : aic->pbl->bodies) + for (int i = 0; i < body->tag->elems.size(); ++i, ++ig) + tau(ig) = (aic->pbl->Ui() + Eigen::Vector3d(uX(ig), uY(ig), uZ(ig))).dot(body->tag->elems[i]->normal()); + std::cout << "done!" << std::endl; + + std::cout << "Solving for body doublets... " << std::flush; + Eigen::PartialPivLU<Eigen::MatrixXd>(aic->A).solve(aic->B * tau); + std::cout << "done!" << std::endl; + + std::cout << "Computing flow and loads on bodies... " << std::flush; + computeFlow(mu, tau, sigma); + computeLoad(); + std::cout << "done!" << std::endl << std::endl; +} + +/** + * @brief Write the results + */ +void Solver::save(std::shared_ptr<MshExport> mshWriter, int n) +{ + // Write files + std::cout << "Saving files... " << std::endl; + // setup results + Results results; + results.scalars_at_nodes["phi"] = φ + results.scalars_at_nodes["varPhi"] = &vPhi; + results.vectors_at_elems["U"] = &U; + results.scalars_at_elems["rho"] = ρ + results.scalars_at_elems["Mach"] = &M; + results.scalars_at_elems["Cp"] = &Cp; + // save (mesh and boundary surface) + if (n > 0) + { + mshWriter->save(aic->pbl->msh->name + "_" + std::to_string(n), results); + for (auto body : aic->pbl->bodies) + body->save(body->tag->name + "_" + std::to_string(n)); + } + else + { + mshWriter->save(aic->pbl->msh->name, results); + for (auto body : aic->pbl->bodies) + body->save(body->tag->name); + } +} + +/** + * @brief Compute flow solution on bodies + */ +void Solver::computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau, const Eigen::VectorXd &sigma) +{ + // @todo should we + // - compute velocity from phi at nodes or mu at nodes? if mu, then store phi at elements only + + // Compute phi at elements and transfer at nodes + Eigen::VectorXd phie = aic->A * mu + aic->B * tau + aic->C * sigma; + auto pbl = aic->pbl; + // reset all values + for (auto n : pbl->msh->nodes) + vPhi[n->no - 1] = 0; + int ig = 0; + for (auto body : pbl->bodies) + { + // associate the potential value to its element + std::map<Element *, double> mapPhi; + for (auto e : body->tag->elems) + { + mapPhi[e] = phie[ig]; + ig++; + } + // average at nodes + for (auto n : body->nodes) + for (auto e : body->neMap.at(n)) + vPhi[n->no - 1] += mapPhi.at(e) / body->neMap.at(n).size(); // maybe use area average? + } + + // Compute surface flow + for (auto body : pbl->bodies) + { + for (auto n : body->nodes) + phi[n->no - 1] = pbl->PhiI(n->pos) + vPhi[n->row]; // total potential + + for (auto e : body->tag->elems) + { + U[e->no - 1] = pbl->U(*e, phi); // velocity + rho[e->no - 1] = pbl->Rho(*e, phi); // density + M[e->no - 1] = pbl->Mach(*e, phi); // Mach number + Cp[e->no - 1] = pbl->Cp(*e, phi); // pressure coefficient + } + } +} + +/** + * @brief Compute aerodynamic load coefficients on bodies + */ +void Solver::computeLoad() +{ + Cl = 0; + Cd = 0; + Cs = 0; + Cm = 0; + auto pbl = aic->pbl; + for (auto body : pbl->bodies) + { + // Compute normalized load (i.e. load / pressure) + for (size_t i = 0; i < body->nodes.size(); ++i) + { + body->cLoadX[i] = 0; + body->cLoadY[i] = 0; + body->cLoadZ[i] = 0; + for (auto e : body->neMap.at(body->nodes[i])) + { + Eigen::Vector3d cLoad = Cp[e->no - 1] * e->getVMem().getVol() * e->normal() / e->nodes.size(); + body->cLoadX[i] += cLoad(0); + body->cLoadY[i] += cLoad(1); + body->cLoadZ[i] += cLoad(2); + } + } + + // Compute aerodynamic load coefficients (i.e. Load / pressure / surface) + // compute coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D) + body->Cm = 0; + Eigen::Vector3d Cf(0, 0, 0); + for (auto e : body->tag->elems) + { + Eigen::Vector3d l = e->getVMem().getCg() - pbl->xR; // lever arm + Eigen::Vector3d cf = Cp[e->no - 1] * e->getVMem().getVol() * e->normal(); // force coefficient + Cf -= cf; + body->Cm -= (l(2) * cf(0) - l(0) * cf(2)) / (pbl->sR * pbl->cR); // moment is positive along y (3D) and negative along z (2D) => positive nose-up + } + // rotate to flow direction + body->Cd = (Cf(0) * cos(pbl->alpha) * cos(pbl->beta) + Cf(1) * sin(pbl->beta) + Cf(2) * sin(pbl->alpha) * cos(pbl->beta)) / pbl->sR; + body->Cs = (-Cf(0) * cos(pbl->alpha) * sin(pbl->beta) + Cf(1) * cos(pbl->beta) - Cf(2) * sin(pbl->alpha) * sin(pbl->beta)) / pbl->sR; + body->Cl = (-Cf(0) * sin(pbl->alpha) + Cf(2) * cos(pbl->alpha)) / pbl->sR; + // compute total + Cl += body->Cl; + Cd += body->Cd; + Cs += body->Cs; + Cm += body->Cm; + } +} + +void Solver::write(std::ostream &out) const +{ + out << "fpm::Solver\nwith\t" << *aic; +} diff --git a/fpm/src/fSolver.h b/fpm/src/fSolver.h new file mode 100644 index 0000000..a7c9084 --- /dev/null +++ b/fpm/src/fSolver.h @@ -0,0 +1,69 @@ +/* + * Copyright 2020 University of Liège + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#ifndef FSOLVER_H +#define FSOLVER_H + +#include "fpm.h" +#include "wObject.h" + +#include <iostream> +#include <vector> +#include <memory> +#include <Eigen/Dense> + +namespace fpm +{ + +/** + * @brief Field panel solver + * @authors Adrien Crovato + * @todo check and implement velocity computation + */ +class FPM_API Solver : public fwk::wSharedObject +{ +public: + std::shared_ptr<Builder> aic; ///< AIC builder + + std::vector<double> phi; ///< full potential + std::vector<double> vPhi; ///< perturbation potential + std::vector<Eigen::Vector3d> U; ///< velocity + std::vector<double> rho; ///< density + std::vector<double> M; ///< mach number + std::vector<double> Cp; ///< pressure coefficient + double Cl; ///< lift coefficient + double Cd; ///< drag coefficient + double Cs; ///< sideforce coefficient + double Cm; ///< pitch moment coefficient (positive nose-up) + +public: + Solver(std::shared_ptr<Builder> _aic); + ~Solver() { std::cout << "~Solver()\n"; } + + void run(); + void save(std::shared_ptr<tbox::MshExport> mshWriter, int n = 0); + +#ifndef SWIG + virtual void write(std::ostream &out) const override; +#endif + +private: + void computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau, const Eigen::VectorXd &sigma); + void computeLoad(); +}; + +} // namespace fpm +#endif //FSOLVER_H diff --git a/fpm/src/fpm.h b/fpm/src/fpm.h index a56c566..b9b586e 100644 --- a/fpm/src/fpm.h +++ b/fpm/src/fpm.h @@ -43,6 +43,10 @@ class Problem; class Body; class Wake; +// solver +class Builder; +class Solver; + // misc class Edge; class EquEdge; diff --git a/fpm/tests/basic.py b/fpm/tests/basic.py index af8a994..8228946 100644 --- a/fpm/tests/basic.py +++ b/fpm/tests/basic.py @@ -33,10 +33,17 @@ def main(): pars = {'spn' : 5, 'nC' : 80, 'nS' : 10} msh = gmsh.MeshLoader('../models/n0012.geo', __file__).execute(**pars) # problem - pbl = fpm.Problem(msh, 0, 0, 5., 1., 0., 0., 0.) + pbl = fpm.Problem(msh, 0, 0, 0, 5., 1., 0., 0., 0.) wing = fpm.Body(msh, 'wing', ['te'], 5) pbl.add(wing) - print(pbl) + # AIC builder + aic = fpm.Builder(pbl) + aic.run() + # solver + slv = fpm.Solver(aic) + slv.run() + slv.save(tbox.GmshExport(msh)) + print('\n', slv) # end timer tms['total'].stop() print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) -- GitLab