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/sdpmSolver.cpp

@@ -522,20 +522,23 @@ void Solver::postUnsteady(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdp
         std::vector<std::vector<sdpmComplex>> cp012(elems.size(), std::vector<sdpmComplex>(3));
         for (size_t i = 0; i < elems.size(); ++i)
         {
-            // si = phi_t + 1/2 conv(u_k,u_k)_k - 1/2 u_inf^2
-            sdpmVectorXcd si(5);
-            si << 0., std::conj(-sdpmComplex(0, 1) * omega * emu(_rows[elems[i]])), -0.5 * ui.dot(ui), -sdpmComplex(0, 1) * omega * emu(_rows[elems[i]]), 0.;
+            // phi_t, phi_x
+            sdpmVector3cd phit(std::conj(-sdpmComplex(0, 1) * omega * emu(_rows[elems[i]])), 0., -sdpmComplex(0, 1) * omega * emu(_rows[elems[i]]));
+            sdpmVector3cd phix(_u1[imd][ifq][elems[i]->getId() - 1].conjugate()(0), _u[elems[i]->getId() - 1](0) - ui(0), _u1[imd][ifq][elems[i]->getId() - 1](0));
+            // 1 - phi_t
+            sdpmVectorXcd ecp(5);
+            ecp << 0., -2. * phit(0) / ui.dot(ui), 1., -2. * phit(2) / ui.dot(ui), 0.;
+            // phi_xx, phi_tt, phi_tx
+            ecp += mi * mi / ui.dot(ui) * (convolve(phix, phix) + convolve(phit, phit) + convolve(phix, phit) / ui.norm());
+            // conv(u_k,u_k)
             for (size_t j = 0; j < 3; ++j)
             {
                 sdpmVector3cd utrm(_u1[imd][ifq][elems[i]->getId() - 1].conjugate()(j), _u[elems[i]->getId() - 1](j), _u1[imd][ifq][elems[i]->getId() - 1](j));
-                si += 0.5 * convolve(utrm, utrm);
+                ecp -= convolve(utrm, utrm) / ui.dot(ui);
             }
-            // 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) * mi * mi / (ui.dot(ui) * ui.dot(ui));
+            // extract 0th, 1st and 2nd harmonics
             for (size_t u = 0; u < 3; ++u)
-                cp012[i][u] = ecp(4 + u);
+                cp012[i][u] = ecp(2 + u);
             _cp1[imd][ifq][elems[i]->getId() - 1] = cp012[i][1];
         }