From f9e3ff23896c2b656f57121dac3a64ff7cfcaddd Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Tue, 4 Apr 2023 13:47:27 +0200
Subject: [PATCH] Update unsteady pressure computation
---
sdpm/src/sdpmSolver.cpp | 21 ++++++++++++---------
sdpm/tests/lann.py | 8 ++++----
2 files changed, 16 insertions(+), 13 deletions(-)
diff --git a/sdpm/src/sdpmSolver.cpp b/sdpm/src/sdpmSolver.cpp
index bcb8c29..fda7e15 100644
--- a/sdpm/src/sdpmSolver.cpp
+++ b/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];
}
diff --git a/sdpm/tests/lann.py b/sdpm/tests/lann.py
index 5b0da28..9975adb 100644
--- a/sdpm/tests/lann.py
+++ b/sdpm/tests/lann.py
@@ -69,10 +69,10 @@ def main():
tests.add(CTest('CD0', sol.getDragCoeff(), 0.0020, 5e-2)) # MATLAB: 0.0031
tests.add(CTest('CS0', sol.getSideforceCoeff(), 0.008, 5e-2))
tests.add(CTest('CM0', sol.getMomentCoeff(), -0.079, 5e-2))
- tests.add(CTest('Abs(CL1)', abs(sol.getUnsteadyLiftCoeff(0, 0)) / amp / math.pi, 1.54, 5e-2)) # MATLAB: 1.62
- tests.add(CTest('Ang(CL1)', cmath.phase(sol.getUnsteadyLiftCoeff(0, 0)) * 180 / math.pi, 0.66, 5e-2)) # MATLAB: 2.6
- tests.add(CTest('Abs(CD1)', abs(sol.getUnsteadyDragCoeff(0, 0)) / amp / math.pi, 0.038, 5e-2)) # MATLAB: 0.041
- tests.add(CTest('Ang(CD1)', cmath.phase(sol.getUnsteadyDragCoeff(0, 0)) * 180 / math.pi, 14.6, 5e-2)) # MATLAB: 16.8
+ tests.add(CTest('Abs(CL1)', abs(sol.getUnsteadyLiftCoeff(0, 0)) / amp / math.pi, 1.61, 5e-2)) # MATLAB: 1.69
+ tests.add(CTest('Ang(CL1)', cmath.phase(sol.getUnsteadyLiftCoeff(0, 0)) * 180 / math.pi, 0.46, 5e-2)) # MATLAB: 0.64
+ tests.add(CTest('Abs(CD1)', abs(sol.getUnsteadyDragCoeff(0, 0)) / amp / math.pi, 0.036, 5e-2)) # MATLAB: 0.037
+ tests.add(CTest('Ang(CD1)', cmath.phase(sol.getUnsteadyDragCoeff(0, 0)) * 180 / math.pi, 16.4, 5e-2)) # MATLAB: 19.9
tests.add(CTest('Abs(CS1)', abs(sol.getUnsteadySideforceCoeff(0, 0)) / amp / math.pi, 0.037, 5e-2))
tests.add(CTest('Ang(CS1)', cmath.phase(sol.getUnsteadySideforceCoeff(0, 0)) * 180 / math.pi, -10.0, 5e-2))
tests.add(CTest('Abs(CM1)', abs(sol.getUnsteadyMomentCoeff(0, 0)) / amp / math.pi, 0.33, 5e-2))
--
GitLab