Version 1.2

Version 1.2 - Transonic correction and analytical gradients

This version mainly features a transonic solver and analytical pressure derivatives.

Changes

• Add transonic correction
• Add analytical gradients
• Add feature to deform the body
• Fix calculation of motion-induced velocity and unsteady pressure
• Refactor APIs and utilities

Notes

Transonic correction
The unsteady pressure solution can now be corrected by higher-fidelity data (typically obtained from CFD) to include transonic and viscous effects.
Analytical gradients
Although CoDiPack works well, its memory footprint is too large for practical unsteady cases (see #6). Analytical pressure and aerodynamic forces gradients have therefore been implemented so that flutter constrained optimization can be performed within a reasonable computational time.
Body deformation
The body nodes can now be set to new positions for aeroelastic and optimization purposes.
Formulation fixes
The calculation of the surface velocity induced by a rigid pitching motion has been corrected to use original instead of compressible coordinates. Additionally, a mistake in the calculation of the unsteady pressure coefficient has been fixed.
Refactoring
The default values in the API have been refactored. The slicing utility now also outputs z-coordinates.

Tests

• Ubuntu 22.04 (gcc 11, python 3.10)

sdpm/src/sdpmBody.cpp

@@ -43,13 +43,13 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
     std::string wkName = _tag->getName() + "Wake";
     // get TE nodes and sort them according to their y-coordinate (ascending)
     Group te(_msh, teName);
-    std::vector<Node *> teNodes;
     for (auto e : te.getElements())
         for (auto n : e->getNodes())
-            teNodes.push_back(n);
-    std::sort(teNodes.begin(), teNodes.end(), [](Node *a, Node *b) -> bool
+            _teNodes.push_back(n);
+    std::sort(_teNodes.begin(), _teNodes.end(), [](Node *a, Node *b) -> bool
               { return a->getCoords()(1) < b->getCoords()(1); });
-    teNodes.erase(std::unique(teNodes.begin(), teNodes.end()), teNodes.end());
+    _teNodes.erase(std::unique(_teNodes.begin(), _teNodes.end()), _teNodes.end());
+    size_t nTe = _teNodes.size() - 1;
     // create wake geometry if it does not already exist
     std::map<std::string, Tag *> tags = _msh.getTags();
     auto it = tags.find(wkName);
@@ -62,38 +62,38 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
             err << "sdpm::Body: zero or negative characteristic chord length given for " << *_tag << "!\n";
             throw std::runtime_error(err.str());
         }
+        // duplicate TE nodes
+        std::map<Node *, Node *> teMap;
+        for (auto n : _teNodes)
+        {
+            Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords());
+            _msh.addNode(nN);
+            teMap[n] = nN;
+        }
         // create tag and add it to the mesh
         Tag *tagp = new Tag(tags.size() + 1, wkName, 2);
         _msh.addTag(wkName, tagp);
         // create wake nodes and elements
-        std::vector<Node *> wkNodes = teNodes;
+        std::vector<Node *> wkNodes = _teNodes;
         for (size_t i = 0; i < 10 * nShedDiv; ++i)
         {
             // translate TE nodes along x-coordinate
             sdpmVector3d delta((i + 1) * shedLgth / nShedDiv, 0., 0.);
-            for (auto n : teNodes)
+            for (auto n : _teNodes)
             {
                 Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords() + delta);
                 wkNodes.push_back(nN);
                 _msh.addNode(nN);
             }
             // create elements
-            for (size_t j = 0; j < teNodes.size() - 1; ++j)
+            for (size_t j = 0; j < nTe; ++j)
             {
-                std::vector<Node *> qnodes = {wkNodes[i * teNodes.size() + j], wkNodes[(i + 1) * teNodes.size() + j], wkNodes[(i + 1) * teNodes.size() + j + 1], wkNodes[i * teNodes.size() + j + 1]};
+                std::vector<Node *> qnodes = {wkNodes[i * (nTe + 1) + j], wkNodes[(i + 1) * (nTe + 1) + j], wkNodes[(i + 1) * (nTe + 1) + j + 1], wkNodes[i * (nTe + 1) + j + 1]};
                 _msh.addElement(new Quad4(_msh.getElements().back()->getId() + 1, tagp, qnodes));
             }
         }
-        // duplicate TE nodes
-        std::map<Node *, Node *> teMap;
-        for (auto n : teNodes)
-        {
-            Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords());
-            _msh.addNode(nN);
-            teMap[n] = nN;
-        }
         // create wake
-        _wake = new Wake(_msh, wkName, _tag, teNodes.size() - 1);
+        _wake = new Wake(_msh, wkName, _tag, nTe);
         // swap nodes
         std::unordered_set<Element *> lwEl; // get set of unique wing TE elements on the pressure side
         for (auto p : _wake->getElMap())
@@ -126,7 +126,7 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
         getUniqueNodes();
         // map the lower TE nodes to their upper nodes
         std::map<Node *, Node *> iTeMap;
-        for (auto ten : teNodes)
+        for (auto ten : _teNodes)
             for (auto wn : _nodes)
                 if (ten->getCoords() == wn->getCoords() && ten->getId() != wn->getId())
                 {
@@ -134,11 +134,11 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
                     break;
                 }
         // create wake and print
-        _wake = new Wake(_msh, wkName, _tag, teNodes.size() - 1, iTeMap);
+        _wake = new Wake(_msh, wkName, _tag, nTe, iTeMap);
         std::cout << *_wake << "already existing, nothing done." << std::endl;
         // check if parameters match existing wake geometry
-        std::vector<Element *> wEls = _wake->getElements();
-        if (std::abs(wEls[0]->getNodes()[1]->getCoords()(0) - wEls[0]->getNodes()[0]->getCoords()(0) - shedLgth / nShedDiv) > 1e-12 || wEls.size() != 10 * nShedDiv * (teNodes.size() - 1))
+        std::vector<Node *> wNods = _wake->getNodes();
+        if (std::abs(wNods[nTe + 1]->getCoords()(0) - wNods[0]->getCoords()(0) - shedLgth / nShedDiv) > 1e-12 || _wake->getElements().size() != 10 * nShedDiv * nTe)
         {
             std::stringstream err;
             err << "sdpm::Body: parameters \'shedLgth\' and \'nShedDiv\' do not match existing wake geometry!\n";
@@ -152,7 +152,8 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
             _neMap[n].push_back(e);
 
     // Size vectors
-    _nLoads.resize(_nodes.size());
+    _dGradPressure.resize(_msh.getNodes().size(), 0.);
+    _nLoads.resize(_nodes.size(), sdpmVector3d::Zero());
     _nLoads1Re.resize(nM, std::vector<std::vector<sdpmVector3d>>(nF, std::vector<sdpmVector3d>(_nodes.size(), sdpmVector3d::Zero())));
     _nLoads1Im.resize(nM, std::vector<std::vector<sdpmVector3d>>(nF, std::vector<sdpmVector3d>(_nodes.size(), sdpmVector3d::Zero())));
     _iVelocityS.resize(nM, std::vector<sdpmVector3d>(getElements().size(), sdpmVector3d::Zero()));
@@ -168,6 +169,26 @@ Body::~Body()
     std::cout << "~sdpm::Body()\n";
 }
 
+/**
+ * @brief Set new coordinates for all nodes
+ */
+void Body::setNodes(std::vector<std::vector<double>> const &coords)
+{
+    // Check size
+    if (coords.size() != _nodes.size())
+        throw std::runtime_error("sdpm::Problem: size of coordinates vector must be equal to the number of nodes!\n");
+    // Set nodes coordinates
+    for (size_t i = 0; i < _nodes.size(); ++i)
+        _nodes[i]->setCoords(sdpmVector3d(coords[i][0], coords[i][1], coords[i][2]));
+    // Update wake
+    sdpmDouble dx = _wake->getLagMap().at(_wake->getElements()[0]); // dx = dt of first cell
+    size_t nTe = _teNodes.size();
+    std::vector<Node *> wkNodes = _wake->getNodes();
+    for (size_t i = 1; i < wkNodes.size() / nTe; ++i)
+        for (size_t j = 0; j < nTe; ++j)
+            wkNodes[i * nTe + j]->setCoords(_teNodes[j]->getCoords() + sdpmVector3d(i * dx, 0., 0.));
+}
+
 /**
  * @brief Add motion
  */
@@ -179,9 +200,9 @@ void Body::addMotion()
 /**
  * @brief Set rigid motion parameters
  */
-void Body::setRigidMotion(size_t m, double aAmp, double hAmp, double xRef, double zRef)
+void Body::setRigidMotion(size_t m, double aAmp, double hAmp)
 {
-    _motions[m]->setRigid(aAmp, hAmp, xRef, zRef);
+    _motions[m]->setRigid(aAmp, hAmp);
 }
 
 /**
@@ -195,17 +216,44 @@ void Body::setFlexibleMotion(size_t m, double mAmp, std::vector<std::vector<doub
     _motions[m]->setFlexible(mAmp, xMod, _nodes);
 }
 
+/**
+ * @brief Set transonic pressure derivative wrt to angle of attack on surface nodes
+ */
+void Body::setTransonicPressureGrad(std::vector<double> const &dCpA)
+{
+    for (size_t i = 0; i < _nodes.size(); ++i)
+        _dGradPressure[_nodes[i]->getId() - 1] = dCpA[i];
+}
+
+/**
+ * @brief Return the elements on the trailing edge on the suction and lower sides
+ */
+std::pair<std::vector<Element *>, std::vector<Element *>> Body::getTrailingEdgeElements() const
+{
+    std::vector<Element *> we = _wake->getElements();
+    std::map<Element *, std::pair<Element *, Element *>> wwMap = _wake->getElMap();
+    size_t nTe = _teNodes.size() - 1;
+    std::vector<Element *> up(nTe), lw(nTe);
+    for (size_t i = 0; i < nTe; ++i)
+    {
+        std::pair<Element *, Element *> p = wwMap.at(we[i]);
+        up[i] = p.first;
+        lw[i] = p.second;
+    }
+    return std::make_pair(up, lw);
+}
+
 /**
  * @brief Compute the motion-induced velocity on the body surface
  */
-void Body::computeVelocity(size_t m, sdpmVector3d const &ui)
+void Body::computeVelocity(size_t m, sdpmVector3d const &ui, sdpmDouble xRef, sdpmDouble zRef)
 {
     std::vector<Element *> elems = getElements();
     std::vector<sdpmVector3d> ucS(elems.size()), ucH(elems.size());
     for (size_t i = 0; i < elems.size(); ++i)
     {
         ucS[i] = _motions[m]->computeSteady(*elems[i], ui);
-        ucH[i] = _motions[m]->computeHarmonic(*elems[i]);
+        ucH[i] = _motions[m]->computeHarmonic(*elems[i], xRef, zRef);
     }
     _iVelocityS[m] = ucS;
     _iVelocityH[m] = ucH;