From b90464fd10e5015670cbe817106eedc6c41f0ce4 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Sat, 26 Oct 2024 13:24:42 +0200
Subject: [PATCH] (fix) Fixed bug in transition
deltaStar was not recomputed for the transition point after averaging solutions, leading to oscillations in the transition region. Test cases have been updated accordingly and deltaStar computation is now tested in the 2D subsonic test case
---
blast/src/DDriver.cpp | 21 +++++++++------------
blast/tests/dart/naca0012_2D.py | 1 +
blast/validation/raeValidation.py | 23 +++++++++++++++--------
3 files changed, 25 insertions(+), 20 deletions(-)
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index 8a4fb03..050c6dc 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -257,13 +257,9 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced)
for (size_t k = iPoint; k < bl->nNodes; ++k)
bl->regime[k] = 1; // Set Turbulent regime.
- // Compute transition location.
- // if (forced)
- // bl->xtr = bl->xtrF;
- // else
- bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->x[iPoint] - bl->x[iPoint-1]) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->x[iPoint-1];
+ bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->xoc[iPoint] - bl->xoc[iPoint-1]) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->xoc[iPoint-1];
- double avTurb = (bl->x[iPoint] - bl->xtr) / (bl->x[iPoint] - bl->x[iPoint-1]);
+ double avTurb = (bl->xoc[iPoint] - bl->xtr) / (bl->xoc[iPoint] - bl->xoc[iPoint-1]);
double avLam = 1. - avTurb;
// std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl;
@@ -293,16 +289,19 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced)
// Solve point in turbulent configuration.
int exitCode = tSolver->integration(iPoint, bl);
- if (exitCode != 0 && verbose > 1)
+ if (exitCode != 0)
std::cout << "Warning: Transition point turbulent computation terminated with exit code " << exitCode << std::endl;
// Average both solutions (Now solution @ iPoint is the turbulent solution).
bl->u[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->u[(iPoint)*nVar + 0]; // Theta.
bl->u[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->u[(iPoint)*nVar + 1]; // H.
- bl->u[iPoint * nVar + 2] = 9.; // N.
+ bl->u[iPoint * nVar + 2] = bl->getnCrit(); // N.
bl->u[iPoint * nVar + 3] = avLam * lamSol[3] + avTurb * bl->u[(iPoint)*nVar + 3]; // ue.
bl->u[iPoint * nVar + 4] = avTurb * bl->u[(iPoint)*nVar + 4]; // Ct.
+ bl->deltaStar[iPoint-1] = bl->u[(iPoint-1)*nVar + 0] * bl->u[(iPoint-1)*nVar + 1];
+ bl->deltaStar[iPoint] = bl->u[iPoint * nVar + 0] * bl->u[iPoint * nVar + 1];
+
// Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations).
std::vector<double> lamParam(8, 0.);
bl->closSolver->laminarClosures(lamParam, bl->u[(iPoint - 1) * nVar + 0], bl->u[(iPoint - 1) * nVar + 1], bl->Ret[iPoint - 1], bl->Me[iPoint-1], bl->name);
@@ -328,7 +327,7 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced)
}
/**
- * @brief Friction, pressure and total drag computation qt each station of the configuration.
+ * @brief Friction, pressure and total drag computation at each station of the configuration.
*
* @tparam Cdt = theta * Ue^((H+5)/2).
* @tparam Cdf = integral(Cf dx)_upper + integral(Cf dx)_lower.
@@ -357,11 +356,9 @@ void Driver::computeSectionalDrag(std::vector<BoundaryLayer *> const &sec, size_
double Cdf_lower = 0.0;
for (size_t k = 1; k < sec[1]->nNodes; ++k)
Cdf_lower += (frictioncoeff_lower[k] + frictioncoeff_lower[k - 1]) * (sec[1]->x[k] - sec[1]->x[k-1]) / 2;
-
- // Friction drag coefficient.
Cdf_sec[i] = Cdf_upper + Cdf_lower; // No friction in the wake
- // Total drag coefficient. Compute from upper and lower TE if there is no wake.
+ // Total drag coefficient. Computed from upper and lower TE if there is no wake.
if (sec.size() > 2 && sec[2]->getName() == "wake")
{
size_t lastWkPt = (sec[2]->nNodes - 1) * nVar;
diff --git a/blast/tests/dart/naca0012_2D.py b/blast/tests/dart/naca0012_2D.py
index 04ab790..a5cd8f1 100644
--- a/blast/tests/dart/naca0012_2D.py
+++ b/blast/tests/dart/naca0012_2D.py
@@ -143,6 +143,7 @@ def main():
tests.add(CTest('Cdf', vsol.Cdf, 0.0047, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
tests.add(CTest('xtr_top', vsol.getAverageTransition(0), 0.196, 0.03, forceabs=True)) # XFOIL 0.1877
tests.add(CTest('xtr_bot', vsol.getAverageTransition(1), 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
+ tests.add(CTest('error dStar', np.linalg.norm(vSolution[0]['deltaStar'] - vSolution[0]['theta']*vSolution[0]['H']), 0., 0., forceabs=True))
tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 19, 0, forceabs=True))
tests.run()
diff --git a/blast/validation/raeValidation.py b/blast/validation/raeValidation.py
index ef462b8..19a1e1e 100644
--- a/blast/validation/raeValidation.py
+++ b/blast/validation/raeValidation.py
@@ -29,6 +29,7 @@ import fwk
from fwk.testing import *
from fwk.coloring import ccolors
import os.path
+import tbox
import math
@@ -42,7 +43,7 @@ def cfgInviscid(nthrds, verb):
'Verb' : verb, # verbosity
# Model (geometry or mesh)
'File' : os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/rae_2.geo', # Input file containing the model
- 'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF': 10, 'msTe' : 0.0075, 'msLe' : 0.001}, # parameters for input file model
+ 'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF': 10, 'msTe' : 0.01, 'msLe' : 0.001}, # parameters for input file model
'Dim' : 2, # problem dimension
'Format' : 'gmsh', # save format (vtk or gmsh)
# Markers
@@ -64,13 +65,13 @@ def cfgInviscid(nthrds, verb):
'y_ref' : 0.0, # reference point for moment computation (y)
'z_ref' : 0.0, # reference point for moment computation (z)
# Numerical
- 'LSolver' : 'SparseLu', # inner solver (Pardiso, MUMPS or GMRES)
+ 'LSolver' : 'PARDISO', # inner solver (Pardiso, MUMPS or GMRES)
'G_fill' : 2, # fill-in factor for GMRES preconditioner
'G_tol' : 1e-5, # tolerance for GMRES
'G_restart' : 50, # restart for GMRES
'Rel_tol' : 1e-6, # relative tolerance on solver residual
'Abs_tol' : 1e-8, # absolute tolerance on solver residual
- 'Max_it' : 50 # solver maximum number of iterations
+ 'Max_it' : 100 # solver maximum number of iterations
}
def cfgBlast(verb):
@@ -84,7 +85,7 @@ def cfgBlast(verb):
'iterPrint': 5, # int, number of iterations between outputs
'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
'sections' : [0],
- 'xtrF' : [0., 0.],# Forced transition location
+ 'xtrF' : [0.03, 0.03],# Forced transition location
'interpolator' : 'Matching',
'nDim' : 2
}
@@ -98,6 +99,12 @@ def main():
icfg = cfgInviscid(args.k, args.verb)
vcfg = cfgBlast(args.verb)
+ try:
+ _ = tbox.Pardiso()
+ except:
+ print('PARDISO not found, using SparseLu')
+ icfg['LSolver'] = 'SparseLu'
+
gr = 1.1
if gr == 1.0:
icfg['Pars']['msF'] = icfg['Pars']['msLe']
@@ -132,12 +139,12 @@ def main():
# Test solution
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- tests.add(CTest('Cl', isol.getCl(), 0.730, 5e-2))
- tests.add(CTest('Cd wake', vsol.Cdt, 0.0095, 1e-3, forceabs=True))
+ tests.add(CTest('Cl', isol.getCl(), 0.740, 5e-2))
+ tests.add(CTest('Cd wake', vsol.Cdt, 0.0096, 1e-3, forceabs=True))
tests.add(CTest('Cd integral', isol.getCd() + vsol.Cdf, 0.0132, 1e-3, forceabs=True))
- tests.add(CTest('Cdf', vsol.Cdf, 0.0068, 1e-3, forceabs=True))
+ tests.add(CTest('Cdf', vsol.Cdf, 0.0067, 1e-3, forceabs=True))
if icfg['LSolver'] == 'SparseLu':
- tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 21, 0, forceabs=True))
+ tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 34, 0, forceabs=True))
tests.run()
expResults = np.loadtxt(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/references/rae2822_AR138_case6.dat')
--
GitLab