From 1f8b1200432c56bc89047dbb6ed1ceac6cc8acff Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Fri, 21 Apr 2023 14:33:55 +0200
Subject: [PATCH] Added forced transition capability in vii.
---
vii/api/core.py | 9 +++++++-
vii/src/DBoundaryLayer.cpp | 1 +
vii/src/DBoundaryLayer.h | 1 +
vii/src/DDiscretization.cpp | 22 +++++++++++++++++++
vii/src/DDiscretization.h | 5 +++++
vii/src/DDriver.cpp | 44 +++++++++++++++++++++++++++++--------
vii/src/DDriver.h | 4 +++-
vii/tests/bli2D.py | 10 +++++----
vii/utils.py | 10 +++++++--
9 files changed, 89 insertions(+), 17 deletions(-)
diff --git a/vii/api/core.py b/vii/api/core.py
index af26fbf..31bda72 100644
--- a/vii/api/core.py
+++ b/vii/api/core.py
@@ -28,6 +28,7 @@ def initVII(cfg, icfg, iSolverName='DART'):
- Minf (float): Freestream Mach number (only used for time step computation).
- CFL0 (float): Initial CFL number.
- Verb (int): Verbosity level of the viscous solver.
+ - xtrF ([float, float]): Forced transition location [upper, lower].
- couplIter (int): Maximum number of coupling iterations.
- couplTol (float): Tolerance to end computation (drag count between 2 iterations).
@@ -76,6 +77,12 @@ def initVII(cfg, icfg, iSolverName='DART'):
_verbose = cfg['Verb']
else:
_verbose = 1
+ _xtrF = [-1, -1]
+ if 'xtrF' in cfg:
+ for i, xtr in enumerate(cfg['xtrF']):
+ if not(0 <= xtr <= 1) and xtr != -1:
+ raise RuntimeError("Incorrect forced transition location.")
+ _xtrF[i] = xtr
_span = 0
_nSections = 1
@@ -98,7 +105,7 @@ def initVII(cfg, icfg, iSolverName='DART'):
__resetInv = False
# Viscous solver object.
- vSolver = vii.Driver(_Re, _Minf, _CFL0, _nSections, _span, verbose =_verbose)
+ vSolver = vii.Driver(_Re, _Minf, _CFL0, _nSections, _xtrF[0], _xtrF[1], _span, verbose =_verbose)
# Solvers interface.
solversAPI = interface(iSolverObjects['sol'], vSolver, [iSolverObjects['blwb'], iSolverObjects['blww']])
# Coupler
diff --git a/vii/src/DBoundaryLayer.cpp b/vii/src/DBoundaryLayer.cpp
index 6301a13..a690d71 100644
--- a/vii/src/DBoundaryLayer.cpp
+++ b/vii/src/DBoundaryLayer.cpp
@@ -19,6 +19,7 @@ BoundaryLayer::BoundaryLayer(double _xtrF)
xtrF = _xtrF;
xtr = 1.0;
+ xoctr = 1.0;
closSolver = new Closures();
mesh = new Discretization();
blEdge = new Edge();
diff --git a/vii/src/DBoundaryLayer.h b/vii/src/DBoundaryLayer.h
index 7c8e223..5d3853c 100644
--- a/vii/src/DBoundaryLayer.h
+++ b/vii/src/DBoundaryLayer.h
@@ -45,6 +45,7 @@ public:
size_t transMarker; ///< Marker id of the transition location.
double xtrF; ///< Forced transition location in the global frame of reference.
double xtr; ///< Transition location in the global frame of reference.
+ double xoctr; ///< Transition location in % of the chord.
BoundaryLayer(double _xtrF = -1.0);
~BoundaryLayer();
diff --git a/vii/src/DDiscretization.cpp b/vii/src/DDiscretization.cpp
index 29b7033..8f4e561 100644
--- a/vii/src/DDiscretization.cpp
+++ b/vii/src/DDiscretization.cpp
@@ -1,5 +1,6 @@
#include "DDiscretization.h"
#include <cmath>
+#include <algorithm>
using namespace vii;
@@ -19,6 +20,7 @@ void Discretization::setGlob(std::vector<double> &_x, std::vector<double> &_y, s
nMarkers = _x.size();
x.resize(nMarkers, 0.);
+ xoc.resize(nMarkers, 0.);
y.resize(nMarkers, 0.);
z.resize(nMarkers, 0.);
nElm = nMarkers - 1;
@@ -27,6 +29,26 @@ void Discretization::setGlob(std::vector<double> &_x, std::vector<double> &_y, s
z = _z;
computeBLcoord();
+ computeOCcoord();
+}
+
+/**
+ * @brief Compute x/chord coordinate.
+*/
+void Discretization::computeOCcoord()
+{
+ // Find min and max of array x
+ auto min_it = std::min_element(x.begin(), x.end());
+ double xMin = *min_it;
+ auto max_it = std::max_element(x.begin(), x.end());
+ double xMax = *max_it;
+
+ // Compute xoc
+ chord = xMax - xMin;
+ if (chord <= 0)
+ throw;
+ for (size_t iPoint = 0; iPoint < x.size(); ++iPoint)
+ xoc[iPoint] = (x[iPoint] - xMin) / chord;
}
/**
diff --git a/vii/src/DDiscretization.h b/vii/src/DDiscretization.h
index 1f9855b..1f40cb8 100644
--- a/vii/src/DDiscretization.h
+++ b/vii/src/DDiscretization.h
@@ -15,6 +15,7 @@ class VII_API Discretization
{
private:
std::vector<double> x; ///< x coordinate in the global frame of reference.
+ std::vector<double> xoc; ///< x/c coordinate in the global frame of reference.
std::vector<double> y; ///< y coordinate in the global frame of reference.
std::vector<double> z; ///< z coordinate in the global frame of reference.
std::vector<double> loc; ///< xi coordinate in the local frame of reference.
@@ -23,6 +24,7 @@ private:
std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference.
size_t nMarkers; ///< Number of points in the domain.
size_t nElm; ///< Number of cells in the domain.
+ double chord; ///< Chord of the section.
public:
size_t prevnMarkers; ///< Number of points in the domain at the previous iteration.
@@ -34,16 +36,19 @@ public:
size_t getnMarkers() const { return nMarkers; };
double getLoc(size_t iPoint) const { return loc[iPoint]; };
double getx(size_t iPoint) const { return x[iPoint]; };
+ double getxoc(size_t iPoint) const { return xoc[iPoint]; };
double gety(size_t iPoint) const { return y[iPoint]; };
double getz(size_t iPoint) const { return z[iPoint]; };
double getExt(size_t iPoint) const { return ext[iPoint]; };
double getAlpha(size_t iElm) const { return alpha[iElm]; };
size_t getDiffnElm() const { return nMarkers - prevnMarkers; };
+ double getChord() const { return chord; };
// Setters.
void setGlob(std::vector<double> &_x, std::vector<double> &_y, std::vector<double> &_z);
void setExt(std::vector<double> extM);
void computeBLcoord();
+ void computeOCcoord();
};
} // namespace vii
#endif // WDISCRETIZATION_H
diff --git a/vii/src/DDriver.cpp b/vii/src/DDriver.cpp
index f226c03..8000696 100644
--- a/vii/src/DDriver.cpp
+++ b/vii/src/DDriver.cpp
@@ -17,16 +17,22 @@ using namespace vii;
* @param _Minf Freestream Mach number.
* @param _CFL0 Initial CFL number for pseudo-time integration.
* @param nSections Number of sections in the computation (1 for 2D cases).
+ * @param xtrF_top Forced transition location on the upper side.
+ * @param xtrF_bot Forced transition location on the lower side.
* @param _Span Span of the considered wing (only used for drag computation).
* @param _verbose Verbosity level of the solver 0<=verbose<=1.
*/
-Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double _Span, unsigned int _verbose)
+Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
+ double xtrF_top, double xtrF_bot, double _span, unsigned int _verbose)
{
// Freestream parameters.
Re = _Re;
CFL0 = _CFL0;
verbose = _verbose;
+ // Vector of forced transition {upper, lower, wake}.
+ std::vector<double> xtrF = {xtrF_top, xtrF_bot, 0.};
+
// Initialzing boundary layer
sections.resize(nSections);
convergenceStatus.resize(nSections);
@@ -35,7 +41,7 @@ Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double
convergenceStatus[iSec].resize(3);
for (size_t i = 0; i < 3; ++i)
{
- BoundaryLayer *region = new BoundaryLayer();
+ BoundaryLayer *region = new BoundaryLayer(xtrF[i]);
sections[iSec].push_back(region);
}
sections[iSec][0]->name = "upper";
@@ -49,7 +55,7 @@ Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double
Cdt = 0.;
Cdf = 0.;
Cdp = 0.;
- span = _Span;
+ span = _span;
// Pre/post processing flags.
stagPtMvmt.resize(sections.size(), false);
@@ -199,7 +205,9 @@ int Driver::run(unsigned int const couplIter)
// Solve equations
tms["1-Solver"].start();
pointExitCode = tSolver->integration(iPoint, sections[iSec][iRegion]);
- // sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar()+2] = 9.;
+ // Force transition if required by the user
+ if (sections[iSec][iRegion]->xtrF != -1 && sections[iSec][iRegion]->mesh->getxoc(iPoint) >= sections[iSec][iRegion]->xtrF)
+ sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar()+2] = 9.;
tms["1-Solver"].stop();
if (pointExitCode != 0)
@@ -210,7 +218,9 @@ int Driver::run(unsigned int const couplIter)
if (!lockTrans)
{
// Check if perturbation waves are growing or not.
- checkWaves(iPoint, sections[iSec][iRegion]);
+ // Avoided in the case of a forced transition.
+ if (sections[iSec][iRegion]->xtrF != -1 && sections[iSec][iRegion]->mesh->getxoc(iPoint) <= sections[iSec][iRegion]->xtrF)
+ checkWaves(iPoint, sections[iSec][iRegion]);
// Amplification factor is compared to critical amplification factor.
if (sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar() + 2] >= sections[iSec][iRegion]->getnCrit())
@@ -220,7 +230,7 @@ int Driver::run(unsigned int const couplIter)
{
std::cout << std::fixed;
std::cout << std::setprecision(2);
- std::cout << sections[iSec][iRegion]->xtr
+ std::cout << sections[iSec][iRegion]->xoctr
<< " (" << sections[iSec][iRegion]->transMarker << ") ";
}
lockTrans = true;
@@ -292,6 +302,7 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
// Compute transition location.
bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getx(iPoint - 1);
+ bl->xoctr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getxoc(iPoint) - bl->mesh->getxoc(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getxoc(iPoint - 1);
// Percentage of the subsection corresponding to a turbulent flow.
double avTurb = (bl->mesh->getx(iPoint) - bl->xtr) / (bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1));
@@ -471,7 +482,7 @@ int Driver::outputStatus(double maxMach) const
std::vector<double> meanTransitions(2, 0.);
for (size_t iSec = 0; iSec < sections.size(); ++iSec)
for (size_t iReg = 0; iReg < 2; ++iReg)
- meanTransitions[iReg] += sections[iSec][iReg]->xtr;
+ meanTransitions[iReg] += sections[iSec][iReg]->xoctr;
for (size_t iReg = 0; iReg < meanTransitions.size(); ++iReg)
meanTransitions[iReg] /= sections.size();
@@ -694,10 +705,10 @@ void Driver::setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue)
std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
{
size_t nMarkersUp = sections[iSec][0]->mesh->getnMarkers();
- size_t nMarkersLw = sections[iSec][1]->mesh->getnMarkers() - 1; // No stagnation point doublon.
+ size_t nMarkersLw = sections[iSec][1]->mesh->getnMarkers(); // No stagnation point doublon.
size_t nVar = sections[iSec][0]->getnVar();
- std::vector<std::vector<double>> Sln(16);
+ std::vector<std::vector<double>> Sln(20);
for (size_t iVector = 0; iVector < Sln.size(); ++iVector)
Sln[iVector].resize(nMarkersUp + nMarkersLw + sections[iSec][2]->mesh->getnMarkers(), 0.);
// Upper side.
@@ -722,6 +733,11 @@ std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
Sln[14][iPoint] = sections[iSec][0]->delta[revPoint];
Sln[15][iPoint] = sections[iSec][0]->mesh->getx(revPoint);
+ Sln[16][iPoint] = sections[iSec][0]->mesh->getxoc(revPoint);
+ Sln[17][iPoint] = sections[iSec][0]->mesh->gety(revPoint);
+ Sln[18][iPoint] = sections[iSec][0]->mesh->getz(revPoint);
+ Sln[19][iPoint] = sections[iSec][0]->blEdge->getVt(revPoint);
+
--revPoint;
}
// Lower side.
@@ -745,6 +761,12 @@ std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
Sln[14][nMarkersUp + iPoint] = sections[iSec][1]->delta[iPoint];
Sln[15][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getx(iPoint);
+ Sln[16][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getxoc(iPoint);
+ Sln[17][nMarkersUp + iPoint] = sections[iSec][1]->mesh->gety(iPoint);
+ Sln[18][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getz(iPoint);
+ Sln[19][nMarkersUp + iPoint] = sections[iSec][1]->blEdge->getVt(iPoint);
+
+
}
// Wake.
for (size_t iPoint = 0; iPoint < sections[iSec][2]->mesh->getnMarkers(); ++iPoint)
@@ -767,6 +789,10 @@ std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
Sln[14][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->delta[iPoint];
Sln[15][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getx(iPoint);
+ Sln[16][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getxoc(iPoint);
+ Sln[17][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->gety(iPoint);
+ Sln[18][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getz(iPoint);
+ Sln[19][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->blEdge->getVt(iPoint);
}
if (verbose > 2)
diff --git a/vii/src/DDriver.h b/vii/src/DDriver.h
index 3651560..3690a68 100644
--- a/vii/src/DDriver.h
+++ b/vii/src/DDriver.h
@@ -37,13 +37,15 @@ protected:
fwk::Timers tms; ///< Internal timers
public:
- Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double _span = 0., unsigned int verbose = 1);
+ Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
+ double xtrF_top=-1., double xtrF_bot=-1., double _span = 0., unsigned int verbose = 1);
~Driver();
int run(unsigned int const couplIter);
// Getters.
double getxtr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xtr; };
double getRe() const { return Re; };
+ double getxoctr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xoctr; };
std::vector<std::vector<double>> getSolution(size_t iSec);
// Setters.
diff --git a/vii/tests/bli2D.py b/vii/tests/bli2D.py
index e661a8b..f6b29a2 100644
--- a/vii/tests/bli2D.py
+++ b/vii/tests/bli2D.py
@@ -25,6 +25,7 @@
# - Transition routines.
# - Compressible flow routines.
# - Laminar and turbulent flow.
+# - Forced transition.
#
# Untested functionalities.
# - Separation routines.
@@ -54,6 +55,7 @@ def main():
Re = 1e7
alpha = 2.*math.pi/180
M_inf = 0.2
+ xtrF = [-1, 0.45]
# Numerical parameters.
CFL0 = 1
@@ -74,7 +76,7 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
- vSolver = viscUtils.initBL(Re, M_inf, CFL0, 1)
+ vSolver = viscUtils.initBL(Re, M_inf, CFL0, 1, xtrF=xtrF)
iSolverAPI = viscUtils.initDart(pbl, blw, vSolver)
Coupler = viiCoupler.Coupler(iSolverAPI, vSolver, _maxCouplIter=20, _couplTol=1e-4, _iterPrint=1, _resetInv=False)
Coupler.run()
@@ -110,11 +112,11 @@ def main():
# Test for n0012 (le=0.01, te=0.01), alpha=2, M=.2, Re=1e7.
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- tests.add(CTest('Cl', iSolverAPI.getCl(), 0.234, 5e-3)) # XFOIL 0.2325
+ tests.add(CTest('Cl', iSolverAPI.getCl(), 0.2344, 5e-3)) # XFOIL 0.2325
tests.add(CTest('Cd', vSolution['Cdt'], 0.0055, 1e-3, forceabs=True)) # XFOIL 0.00531
- tests.add(CTest('Cdp', vSolution['Cdp'], 0.0012, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
+ tests.add(CTest('Cdp', vSolution['Cdp'], 0.0010, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
tests.add(CTest('xtr_top', vSolution['xtrT'], 0.21, 0.03, forceabs=True)) # XFOIL 0.1877
- tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.50, 0.03, forceabs=True)) # XFOIL 0.4998
+ tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.45, 0.01, forceabs=True)) # XFOIL 0.4998
tests.run()
# eof
diff --git a/vii/utils.py b/vii/utils.py
index 04a1178..c090dc2 100644
--- a/vii/utils.py
+++ b/vii/utils.py
@@ -3,7 +3,7 @@ from fwk.wutils import parseargs
from fwk.coloring import ccolors
import numpy as np
-def initBL(Re, Minf, CFL0, nSections, span=0, verb=None):
+def initBL(Re, Minf, CFL0, nSections, xtrF=[-1, -1], span=0, verb=None):
""" Initialize boundary layer solver.
Params
@@ -12,6 +12,7 @@ def initBL(Re, Minf, CFL0, nSections, span=0, verb=None):
- Minf: Freestream Mach number.
- CFL0: Initial CFL number for time integration.
- nSections: Number of sections in the domain.
+ - xtrF: Forced transition location [upper, lower].
- span: Wing span (not used for 2D computations.
- verb: Verbosity level of the solver.
@@ -39,7 +40,12 @@ def initBL(Re, Minf, CFL0, nSections, span=0, verb=None):
raise RuntimeError("Invalid parameter")
else:
_verbose = verb
- solver = vii.Driver(Re, Minf, CFL0, nSections, span, verbose=_verbose)
+ _xtrF = [-1, -1]
+ for i, xtr in enumerate(xtrF):
+ if not(0 <= xtr <= 1) and xtr != -1:
+ raise RuntimeError("Incorrect forced transition location.")
+ _xtrF[i] = xtr
+ solver = vii.Driver(Re, Minf, CFL0, nSections, xtrF[0], xtrF[1], span, verbose=_verbose)
return solver
def initDart(pbl, blw, vSolver, iters = 25):
--
GitLab