From cd548558319dced7b13020e3038a9246ac1a18f3 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Wed, 16 Jun 2021 15:15:06 +0200
Subject: [PATCH] Consistent signs for forces due to pressure
---
flow/src/wAdjoint.cpp | 14 +++++++-------
flow/src/wLoadFunctional.cpp | 8 ++++----
flow/src/wSolver.cpp | 4 ++--
3 files changed, 13 insertions(+), 13 deletions(-)
diff --git a/flow/src/wAdjoint.cpp b/flow/src/wAdjoint.cpp
index a5d9be24..11eebf98 100644
--- a/flow/src/wAdjoint.cpp
+++ b/flow/src/wAdjoint.cpp
@@ -267,8 +267,8 @@ void Adjoint::buildGradientLoadsFlow(Eigen::VectorXd &dL, Eigen::VectorXd &dD)
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t i = 0; i < e->nodes.size(); ++i)
{
- Eigen::RowVectorXd AeL = -dirL * Ae.block(i * 3, 0, 3, eV->nodes.size());
- Eigen::RowVectorXd AeD = -dirD * Ae.block(i * 3, 0, 3, eV->nodes.size());
+ Eigen::RowVectorXd AeL = dirL * Ae.block(i * 3, 0, 3, eV->nodes.size());
+ Eigen::RowVectorXd AeD = dirD * Ae.block(i * 3, 0, 3, eV->nodes.size());
for (size_t j = 0; j < eV->nodes.size(); ++j)
{
Node *nodj = eV->nodes[j];
@@ -342,7 +342,7 @@ void Adjoint::buildGradientLoadsAoa(double &dL, double &dD)
Eigen::Vector3d Cf(0, 0, 0);
for (auto bnd : sol->pbl->bnds)
for (size_t i = 0; i < bnd->nodes.size(); ++i)
- Cf -= bnd->nLoads[i];
+ Cf += bnd->nLoads[i];
// Rotate to gradient of flow direction
dD = Cf.dot(dirD);
dL = Cf.dot(dirL);
@@ -533,10 +533,10 @@ void Adjoint::buildGradientLoadsMesh(Eigen::RowVectorXd &dL, Eigen::RowVectorXd
for (size_t i = 0; i < e->nodes.size(); ++i)
{
// rotate to flow direction
- Eigen::RowVectorXd AesL = -dirL * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
- Eigen::RowVectorXd AesD = -dirD * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
- Eigen::RowVectorXd AevL = -dirL * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
- Eigen::RowVectorXd AevD = -dirD * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
+ Eigen::RowVectorXd AesL = dirL * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
+ Eigen::RowVectorXd AesD = dirD * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
+ Eigen::RowVectorXd AevL = dirL * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
+ Eigen::RowVectorXd AevD = dirD * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
for (size_t m = 0; m < sol->pbl->nDim; ++m)
{
// surface
diff --git a/flow/src/wLoadFunctional.cpp b/flow/src/wLoadFunctional.cpp
index 6c0bbe17..1113e3d8 100644
--- a/flow/src/wLoadFunctional.cpp
+++ b/flow/src/wLoadFunctional.cpp
@@ -36,7 +36,7 @@ Eigen::VectorXd LoadFunctional::build(Element const &e, Element const &eV, std::
Gauss &gauss = cache.getVGauss();
// Pressure force coefficient in associated volume
- Eigen::Vector3d fcp = cp.eval(eV, phi, 0) * e.getNormal();
+ Eigen::Vector3d fcp = -cp.eval(eV, phi, 0) * e.getNormal(); // "-" because force on body not on fluid
// Integrate cp on the surface
Eigen::VectorXd b = Eigen::VectorXd::Zero(3 * e.nodes.size());
@@ -62,7 +62,7 @@ Eigen::MatrixXd LoadFunctional::buildGradientFlow(Element const &e, Element cons
Gauss &gauss = cache.getVGauss();
// Gradient of pressure force coefficient in associated volume
- Eigen::MatrixXd dfcp = e.getNormal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getJinv(0) * eV.getVCache().getDsf(0);
+ Eigen::MatrixXd dfcp = -e.getNormal() * cp.evalD(eV, phi, 0) * eV.computeGradient(phi, 0).transpose() * eV.getJinv(0) * eV.getVCache().getDsf(0); // "-" because force on body not on fluid
// Integrate cp on the surface
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(3 * e.nodes.size(), eV.nodes.size());
@@ -91,10 +91,10 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> LoadFunctional::buildGradientMesh(E
Gauss &gauss = cache.getVGauss();
// Variables
- double cP = cp.eval(eV, phi, 0); // pressure coefficient
+ double cP = -cp.eval(eV, phi, 0); // pressure coefficient ("-" because force on body not on fluid)
Eigen::Vector3d nrm = e.getNormal(); // unit normal
// Gradients
- double dCp = cp.evalD(eV, phi, 0); // pressure coefficient
+ double dCp = -cp.evalD(eV, phi, 0); // pressure coefficient ("-" because force on body not on fluid)
Eigen::VectorXd dPhi = eV.computeGradient(phi, 0); // potential
Eigen::MatrixXd const &iJ = eV.getJinv(0); // inverse Jacobian
std::vector<Eigen::MatrixXd> const &dJ = eV.getGradJ(0); // Jacobian
diff --git a/flow/src/wSolver.cpp b/flow/src/wSolver.cpp
index 1b15b049..7fae0493 100644
--- a/flow/src/wSolver.cpp
+++ b/flow/src/wSolver.cpp
@@ -197,8 +197,8 @@ void Solver::computeLoad()
{
Eigen::Vector3d l = sur->nodes[i]->pos - pbl->x_ref; // lever arm
Eigen::Vector3d cf = sur->nLoads[i]; // normalized force
- Cf -= cf;
- sur->Cm -= (l(pbl->nDim - 1) * cf(0) - l(0) * cf(pbl->nDim - 1)) / (pbl->S_ref * pbl->c_ref); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
+ Cf += cf;
+ sur->Cm += (l(pbl->nDim - 1) * cf(0) - l(0) * cf(pbl->nDim - 1)) / (pbl->S_ref * pbl->c_ref); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
}
// rotate to flow direction
sur->Cd = Cf.dot(dirD);
--
GitLab