USCSDPM v1.0

Version 1.0 - Unsteady compressible source and doublet panel method

USCSDPM, shortened SDPM, is an implementation of the original source and doublet panel method proposed by Morino (L. Morino and P. Chen, Steady and oscillatory subsonic and supersonic aerodynamics around complex configurations, AIAA Journal, 1975) and revisited by Dimitriadis (M. Martínez and G. Dimitriadis, Subsonic source and doublet panel methods, Journal of Fluids and Structures, 2022). The main purpose of SDPM is to provide steady and unsteady pressure distributions on fixed and moving wings. The code is currently capable of computing steady or unsteady, compressible or incompressible flows on isolated or several wing surfaces.

Features

• C++ code interfaced in Python using SWIG
• Steady or unsteady flows caused by forced harmonic pitching motion of the boundary
• Incompressible or compressible flows
• Isolated or multiple wings

Limitations

• Flow are currently limited to the compressible subsonic, but not transonic regime
• No gradient computation method is currently featured
• Only rigid pitching motion can be currently imposed on the moving boundary

These limitations will be addressed in future developments.

Related issues

• #4 (closed) implemented
• #5 (closed) fixed

Tests

• Ubuntu 22.04 (gcc 11, python 3.10)
• Windows 10 (vs 2022, python 3.8)

sdpm/src/sdpmBody.cpp

+/*
+ * Copyright 2023 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 "sdpm.h"
+#include "sdpmBody.h"
+#include "sdpmMesh.h"
+#include "sdpmElement.h"
+#include "sdpmQuad4.h"
+#include "sdpmTag.h"
+#include "sdpmGmshExport.h"
+#include <unordered_set>
+#include <iostream>
+using namespace sdpm;
+
+Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
+           double shedLgth, size_t nShedDiv,
+           size_t nM, size_t nF,
+           double ux) : Group(msh, name),
+                        _cl(0), _cd(0), _cs(0), _cm(0),
+                        _cl1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.))),
+                        _cd1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.))),
+                        _cs1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.))),
+                        _cm1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.)))
+{
+    // Create wake
+    std::cout << "Creating wake... " << std::flush;
+    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
+              { return a->getCoords()(1) < b->getCoords()(1); });
+    teNodes.erase(std::unique(teNodes.begin(), teNodes.end()), teNodes.end());
+    // create wake geometry if it does not already exist
+    std::map<std::string, Tag *> tags = _msh.getTags();
+    auto it = tags.find(wkName);
+    if (it == tags.end())
+    {
+        // sanity check
+        if (shedLgth <= 0.)
+        {
+            std::stringstream err;
+            err << "sdpm::Body: zero or negative characteristic chord length given for " << *_tag << "!\n";
+            throw std::runtime_error(err.str());
+        }
+        // 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;
+        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)
+            {
+                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)
+            {
+                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]};
+                _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, ux);
+        // swap nodes
+        std::unordered_set<Element *> lwEl; // get set of unique wing TE elements on the pressure side
+        for (auto p : _wake->getElMap())
+            lwEl.insert(p.second.second);
+        for (auto e : lwEl)
+        {
+            std::vector<Node *> lwNods = e->getNodes();
+            for (size_t i = 0; i < lwNods.size(); ++i)
+            {
+                try
+                {
+                    e->setNode(i, teMap.at(lwNods[i]));
+                }
+                catch (const std::out_of_range &)
+                {
+                    // std::cout << lwe->nodes[i]->no << "not found in map!\n";
+                }
+            }
+        }
+        // get body nodes (now that they are duplicated and swaped at lower TE)
+        getUniqueNodes();
+        // print and save new mesh
+        std::cout << *_wake << "created. " << std::flush;
+        GmshExport writer(_msh);
+        writer.save();
+    }
+    else
+    {
+        // get body nodes
+        getUniqueNodes();
+        // map the lower TE nodes to their upper nodes
+        std::map<Node *, Node *> iTeMap;
+        for (auto ten : teNodes)
+            for (auto wn : _nodes)
+                if (ten->getCoords() == wn->getCoords() && ten->getId() != wn->getId())
+                {
+                    iTeMap[wn] = ten;
+                    break;
+                }
+        // create wake and print
+        _wake = new Wake(_msh, wkName, _tag, teNodes.size() - 1, ux, 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::stringstream err;
+            err << "sdpm::Body: parameters \'shedLgth\' and \'nShedDiv\' do not match existing wake geometry!\n";
+            throw std::runtime_error(err.str());
+        }
+    }
+
+    // Associate each node to its adjacent elements
+    for (auto e : _tag->getElements())
+        for (auto n : e->getNodes())
+            _neMap[n].push_back(e);
+
+    // Size vectors
+    _nLoads.resize(_nodes.size());
+    _nLoads1.resize(nM, std::vector<std::vector<sdpmVector3cd>>(nF, std::vector<sdpmVector3cd>(_nodes.size(), sdpmVector3cd::Zero())));
+}
+
+Body::~Body()
+{
+    delete _wake;
+    for (auto m : _motions)
+        delete m;
+    _motions.clear();
+    std::cout << "~sdpm::Body()\n";
+}
+
+/**
+ * @brief Add motion
+ */
+void Body::addMotion(double ux, double amp, double xr)
+{
+    _motions.push_back(new Motion(ux, amp, xr));
+}
+
+/**
+ * @brief Compute the motion-induced velocity on the body surface
+ */
+void Body::computeVelocity(size_t m)
+{
+    std::vector<Element *> elems = getElements();
+    std::vector<sdpmVector3cd> ucS(elems.size()), ucH(elems.size());
+    for (size_t i = 0; i < elems.size(); ++i)
+    {
+        ucS[i] = _motions[m]->computeSteady(*elems[i]);
+        ucH[i] = _motions[m]->computeHarmonic(*elems[i]);
+    }
+    _iVelocityS.push_back(ucS);
+    _iVelocityH.push_back(ucH);
+}
+
+/**
+ * @brief Return a vector of unique nodes belonging to the tag
+ */
+void Body::getUniqueNodes()
+{
+    for (auto e : _tag->getElements())
+        for (auto n : e->getNodes())
+            _nodes.push_back(n);
+    std::sort(_nodes.begin(), _nodes.end());
+    _nodes.erase(std::unique(_nodes.begin(), _nodes.end()), _nodes.end());
+}
+
+void Body::write(std::ostream &out) const
+{
+    out << "sdpm::Body " << *_tag << "with " << *_wake;
+}