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

@@ -34,7 +34,7 @@ Solver::Solver(Problem &pbl, int vrb) : _pbl(pbl), _vrb(vrb), _cl(0), _cd(0), _c
     // Say hi
     std::cout << std::endl;
     std::cout << "*******************************************************************************" << std::endl;
-    std::cout << "** Hi! My name is SDPM v1.1-2401                                             **" << std::endl;
+    std::cout << "** Hi! My name is SDPM v1.2-2412                                             **" << std::endl;
     std::cout << "** ULiege 2023-2024                                                          **" << std::endl;
     std::cout << "*******************************************************************************" << std::endl;
 
@@ -558,7 +558,7 @@ void Solver::post(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdpmVectorX
             sdpmComplex phit = -sdpmComplex(0., 1.) * omega * emu(_rows[elems[i]]);
             sdpmComplex phix = u1[i](0) - (uc0[i](0) + uc1[i](0) * sdpmComplex(0., 1.) * omega);
             sdpmComplex phix0(_u[elems[i]->getId() - 1](0) - ui(0), 0.);
-            sdpmComplex cp1 = -sdpmComplex(2., 0.) * phit + sdpmComplex(mi * mi, 0.) * (sdpmComplex(2., 0.) * phix0 * phix + phix0 * phit);
+            sdpmComplex cp1 = -sdpmComplex(2., 0.) * phit + sdpmComplex(2 * mi * mi, 0.) * (phix0 * phix + phix0 * phit);
             for (size_t j = 0; j < 3; ++j)
                 cp1 -= sdpmComplex(2 * _u[elems[i]->getId() - 1](j), 0.) * u1[i](j);
             // cp
@@ -616,7 +616,9 @@ void Solver::computeGAF(size_t ifq)
     // Get sizes
     size_t n = _pbl.getNodes().size();
     size_t nm = _pbl.getNumModes();
-    // Get z-component of modal displacements and nodal loads
+    // Get x-component of reference center
+    sdpmDouble xRef = _pbl.getRefCtr()(0);
+    // Get z-component of displacements and nodal loads
     sdpmMatrixXd mz1 = sdpmMatrixXd::Zero(nm, n);
     sdpmMatrixXd lz1Re = sdpmMatrixXd::Zero(n, nm);
     sdpmMatrixXd lz1Im = sdpmMatrixXd::Zero(n, nm);
@@ -624,12 +626,15 @@ void Solver::computeGAF(size_t ifq)
     {
         for (auto b : _pbl.getBodies())
         {
+            sdpmDouble aAmp = b->getPitchAmplitude(imd);
+            sdpmDouble hAmp = b->getPlungeAmplitude(imd);
             mz1.row(imd) = Eigen::Map<const sdpmVectorXd>(b->getModeZ(imd).data(), n).transpose();
             std::vector<Node *> nods = b->getNodes();
             std::vector<sdpmVector3d> l1Re = b->getUnsteadyLoadsReal(imd, ifq);
             std::vector<sdpmVector3d> l1Im = b->getUnsteadyLoadsImag(imd, ifq);
             for (size_t i = 0; i < nods.size(); ++i)
             {
+                mz1(imd, nods[i]->getId() - 1) += aAmp * (nods[i]->getCoords()(0) - xRef) + hAmp;
                 lz1Re(nods[i]->getId() - 1, imd) = -l1Re[i](2);
                 lz1Im(nods[i]->getId() - 1, imd) = -l1Im[i](2);
             }