diff --git a/dart/tests/bliHighLift.py b/dart/tests/bliHighLift.py
index e5783674b232189a839e262b2d1b4dd011e07ee1..022ee4b8991cabf951e915eff74b94b73be394d5 100755
--- a/dart/tests/bliHighLift.py
+++ b/dart/tests/bliHighLift.py
@@ -141,7 +141,7 @@ def main():
floD.initViewer(pbl)
tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
- val=validation.Validation(case)
+ val=validation.Validation()
val.main(isolver.Cl,vsolver.wData)
# eof
diff --git a/dart/tests/bliLowRe.py b/dart/tests/bliLowRe.py
index 98f0ce3c72238c5ec4935c6a1f9f324b40f9f980..5e67e4e77c461ae34053f05b6c9f2deeb7aaeff2 100755
--- a/dart/tests/bliLowRe.py
+++ b/dart/tests/bliLowRe.py
@@ -86,7 +86,7 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.0001}
#pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.002, 'msLe' : 0.0001}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
@@ -131,13 +131,13 @@ def main():
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
if Re == 5e5 and M_inf == 0.2 and alpha == 2*math.pi/180:
- tests.add(CTest('Cl', isolver.Cl, 0.2303, 5e-2)) # XFOIL 0.2135
- tests.add(CTest('Cd', vsolver.Cd, 0.0089, 1e-3, forceabs=True)) # XFOIL 0.00706
- tests.add(CTest('Cdp', vsolver.Cdp, 0.0041, 1e-3, forceabs=True)) # 0.00238
- tests.add(CTest('xtr_top', vsolver.xtr[0], 0.574, 0.03, forceabs=True)) # XFOIL 0.5642
- tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.89, 0.03, forceabs=True)) # XFOIL 0.9290
+ tests.add(CTest('Cl', isolver.Cl, 0.2206, 5e-2)) # XFOIL 0.2135
+ tests.add(CTest('Cd', vsolver.Cd, 0.007233, 1e-3, forceabs=True)) # XFOIL 0.00706
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 1e-3, forceabs=True)) # 0.00238
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.58, 0.03, forceabs=True)) # XFOIL 0.5642
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.90, 0.03, forceabs=True)) # XFOIL 0.9290
else:
- raise Exception('Test not defined for this flow')
+ print(ccolors.ANSI_RED + 'Test not defined for this flow', ccolors.ANSI_RESET)
tests.run()
diff --git a/dart/tests/bliNonLift.py b/dart/tests/bliNonLift.py
index 9fdd94cf0929dc6ed73a31cfebc55911d1570b4e..674660638eb822e759637b6cd0f6c4ebe166d8ad 100755
--- a/dart/tests/bliNonLift.py
+++ b/dart/tests/bliNonLift.py
@@ -60,20 +60,20 @@ def main():
# define flow variables
Re = 1e7
- alpha = 5.*math.pi/180
- M_inf = 0.2
+ alpha = 0.*math.pi/180
+ M_inf = 0.
#user_xtr=[0.41,0.41]
#user_xtr=[0,None]
user_xtr=[None,None]
user_Ti=None
- mapMesh = 1
+ mapMesh = 0
nFactor = 2
# Time solver parameters
Time_Domain = True # True for unsteady solver, False for steady solver
- CFL0 = 2
+ CFL0 = 5
# ========================================================================================
diff --git a/dart/tests/bliPolar.py b/dart/tests/bliPolar.py
new file mode 100755
index 0000000000000000000000000000000000000000..c4d6fd560b10d6f31071377c703826d7c9925471
--- /dev/null
+++ b/dart/tests/bliPolar.py
@@ -0,0 +1,167 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+# -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+# -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+# -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+# -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+from numpy.core.function_base import linspace
+import dart.utils as floU
+import dart.default as floD
+
+import dart.viscous.Master.DCoupler as floVC
+import dart.viscous.Drivers.DDriver as floVSMaster
+
+import dart.viscous.Solvers.Steady.StdCoupler as floVC_Std
+import dart.viscous.Solvers.Steady.StdSolver as floVS_Std
+
+import dart.viscous.Physics.DValidation as validation
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+import numpy as np
+from matplotlib import pyplot as plt
+
+import math
+
+def main():
+ # timer
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ # define flow variables
+ Re = 6e6
+ #alpha = 0.*math.pi/180
+ M_inf = 0.15
+
+ #user_xtr=[0.41,0.41]
+ #user_xtr=[0,None]
+ user_xtr=[0, 0]
+ user_Ti=None
+
+ mapMesh = 0
+ nFactor = 2
+
+ # Time solver parameters
+ Time_Domain = True # True for unsteady solver, False for steady solver
+ CFL0 = 1
+
+ # ========================================================================================
+
+ #U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+ c_ref = 1
+ dim = 2
+ tol = 1e-4 # tolerance for the VII (usually one drag count)
+ case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
+
+ # mesh the geometry
+ print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+ tms['msh'].start()
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.0001}
+ #pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.002, 'msLe' : 0.0001}
+ msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+ tms['msh'].stop()
+
+ alphaVect = [
+ 15.2600*math.pi/180,
+ 16.3000*math.pi/180,
+ 16.5*math.pi/180,
+ 17.1300*math.pi/180,
+ 17.2*math.pi/180,
+ 17.5*math.pi/180,
+ 17.7*math.pi/180,
+ 18.0200*math.pi/180]
+
+ alphaSuc = []
+ clSuc = []
+ cdSuc = []
+ alphaFailed = []
+ clFailed = []
+ cdFailed = []
+
+ for alpha in alphaVect:
+ # set the problem and add medium, initial/boundary conditions
+ tms['pre'].start()
+ pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+ tms['pre'].stop()
+
+ # solve viscous problem
+
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+ isolver = floD.newton(pbl)
+
+ # Choose between steady and unsteady solver
+ if Time_Domain is True: # Run unsteady solver
+ vsolver=floVSMaster.Driver(Re, user_xtr, CFL0, M_inf, case, mapMesh, nFactor)
+ coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter, bnd)
+
+ elif Time_Domain is False: # Run steady solver
+ vsolver=floVS_Std.Solver(Re)
+ coupler = floVC_Std.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
+ try:
+ coupler.run()
+ tms['solver'].stop()
+
+ # extract Cp
+ Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+ tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+ # display results
+ print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+ print(' Re M alpha Cl Cd Cdp Cdf Cm')
+ print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+
+ # display timers
+ tms['total'].stop()
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+
+ alphaSuc.append(alpha*180/math.pi)
+ clSuc.append(isolver.Cl)
+ cdSuc.append(vsolver.Cd)
+ except KeyboardInterrupt:
+ quit()
+ except:
+ alphaFailed.append(alpha*180/math.pi)
+ clFailed.append(isolver.Cl)
+ cdFailed.append(vsolver.Cd)
+
+
+ print(alphaSuc, clSuc, cdSuc)
+ print(alphaFailed, clFailed, cdFailed)
+
+if __name__ == "__main__":
+ main()
diff --git a/dart/tests/bliSeparated.py b/dart/tests/bliSeparated.py
index 7b3ef4441ce5ae68f19fb4d944e70a14c601d6cd..49822dec4b639b2d9676cbc7bcccc5075654de51 100755
--- a/dart/tests/bliSeparated.py
+++ b/dart/tests/bliSeparated.py
@@ -60,9 +60,9 @@ def main():
# define flow variables
Re = 1e7
- alpha = 15.*math.pi/180
+ alpha = 15*math.pi/180
M_inf = 0.2
- user_xtr=[None,None]
+ user_xtr=[None, None]
user_Ti=None
mapMesh = 1
@@ -127,11 +127,17 @@ def main():
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
if Re == 1e7 and M_inf == 0.2 and alpha == 15*math.pi/180:
- tests.add(CTest('Cl', isolver.Cl, 1.60, 5e-2)) # Xfoil 1.6654
- tests.add(CTest('Cd', vsolver.Cd, 0.0154, 1e-3, forceabs=True)) # XFOIL 0.01640
- tests.add(CTest('Cdp', vsolver.Cdp, 0.0115, 1e-3, forceabs=True)) # XFOIL 0.01255
- tests.add(CTest('xtr_top', vsolver.xtr[0], 0.00676, 0.003, forceabs=True)) # XFOIL 0.0067
- tests.add(CTest('xtr_bot', vsolver.xtr[1], 1, 0.03, forceabs=True)) # XFOIL 1.0000
+ tests.add(CTest('Cl', isolver.Cl, 1.615, 5e-2)) # Xfoil 1.6654
+ tests.add(CTest('Cd', vsolver.Cd, 0.0143, 1e-3, forceabs=True)) # XFOIL 0.01460
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0103, 1e-3, forceabs=True)) # XFOIL 0.01255
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.00660, 0.003, forceabs=True)) # XFOIL 0.0067
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 1, 0.03, forceabs=True)) # XFOIL 1.0000
+ elif Re == 5e6 and M_inf == 0.2 and alpha == 15*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 1.58, 5e-2)) # Xfoil 1.6109
+ tests.add(CTest('Cd', vsolver.Cd, 0.0201, 1e-3, forceabs=True)) # XFOIL 0.01944
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0130, 1e-3, forceabs=True)) # XFOIL 0.01544
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.008, 0.003, forceabs=True)) # XFOIL 0.0081
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 1, 0.03, forceabs=True)) # XFOIL 1.0000
else:
raise Exception('Test not defined for this flow')
@@ -141,7 +147,7 @@ def main():
floD.initViewer(pbl)
tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
- val=validation.Validation(case)
+ val=validation.Validation()
val.main(isolver.Cl,vsolver.wData)
# eof
diff --git a/dart/tests/bliTransonic.py b/dart/tests/bliTransonic.py
index f67b686932c44b948c254f0c09d6122f8c6532dd..c29ee906ffb6d1829ea00f41e0ff26c0e5e64264 100755
--- a/dart/tests/bliTransonic.py
+++ b/dart/tests/bliTransonic.py
@@ -50,6 +50,7 @@ import tbox.utils as tboxU
import fwk
from fwk.testing import *
from fwk.coloring import ccolors
+import numpy as np
import math
@@ -61,24 +62,24 @@ def main():
# define flow variables
Re = 1e7
alpha = 2*math.pi/180
- M_inf = 0.715
+ M_inf = 0.73
- user_xtr=[None, None]
+ user_xtr=[0, 0]
user_Ti=None
mapMesh = 1
- nFactor = 5
+ nFactor = 2
# Time solver parameters
Time_Domain = True # True for unsteady solver, False for steady solver
- CFL0 = 1
+ CFL0 = 2
# ========================================================================================
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
c_ref = 1
dim = 2
- tol = 1e-3 # tolerance for the VII (usually one drag count)
+ tol = 1e-4 # tolerance for the VII (usually one drag count)
case='Case1FreeTrans.dat' # File containing results from XFOIL of the considered case
# mesh the geometry
@@ -130,17 +131,31 @@ def main():
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- tests = CTests()
- if Re == 1e7 and M_inf == 0.715 and alpha == 2*math.pi/180:
- tests.add(CTest('Cl', isolver.Cl, 0.686, 5e-2))
- tests.add(CTest('Cd', vsolver.Cd, 0.0057, 1e-3, forceabs=True))
- tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 1e-3, forceabs=True))
- tests.add(CTest('xtr_top', vsolver.xtr[0], 0.34, 0.03, forceabs=True))
- tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.5, 0.03, forceabs=True))
- else:
- raise Exception('Test not defined for this flow')
-
- tests.run()
+ try:
+ tests = CTests()
+ if Re == 1e7 and M_inf == 0.715 and alpha == 2*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.72, 5e-2))
+ tests.add(CTest('Cd', vsolver.Cd, 0.0061, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.34, 0.08, forceabs=True))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.5, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0.72 and alpha == 2*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.674, 5e-2))
+ tests.add(CTest('Cd', vsolver.Cd, 0.0061, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.36, 0.08, forceabs=True))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.5, 0.03, forceabs=True))
+ elif Re == 1e7 and M_inf == 0.73 and alpha == 2*math.pi/180 and user_xtr==[0, 0]:
+ tests.add(CTest('Cl', isolver.Cl, 0.671, 5e-2))
+ tests.add(CTest('Cd', vsolver.Cd, 0.0089, 1e-3, forceabs=True))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0024, 1e-3, forceabs=True))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0, 0.08, forceabs=True))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0, 0.03, forceabs=True))
+ else:
+ raise Exception('Test not defined for this flow')
+ tests.run()
+ except Exception as e:
+ print(e)
# visualize solution and plot results
floD.initViewer(pbl)
diff --git a/dart/viscous/Drivers/DDriver.py b/dart/viscous/Drivers/DDriver.py
index 2eea39b95aee5137292aa60d1afb7407d43f7e82..eb4e4b2a88aa3ca4caf2f9f4f4154aa5b608dc8b 100644
--- a/dart/viscous/Drivers/DDriver.py
+++ b/dart/viscous/Drivers/DDriver.py
@@ -33,227 +33,238 @@ import numpy as np
import matplotlib.pyplot as plt
class Driver:
- """Driver of the time-dependent boundary layer equation solver (viscous solver)
- """
-
- def __init__(self,_Re, user_xtr, _CFL0, _Minfty, _case, _mapMesh, _nFactor):
- self.Re=_Re
- self.Minfty = _Minfty
- self.xtrF=user_xtr
- self.Ti=None
- self.Cd=0
- self.it=0
- self.uPrevious=None
- self.mapMesh = _mapMesh
- self.nFactor = _nFactor
- self.CFL0 = _CFL0
- self.res=None
- self.upper=0
- self.case=_case
- self.xtr=[1,1]
- self.nbrElmUpper = -1
- pass
-
- def dictionary(self):
- """Create dictionary with the coordinates and the boundary layer parameters
- """
-
- wData = {}
- wData['x'] = self.data[:,0]
- wData['y'] = self.data[:,1]
- wData['z'] = self.data[:,2]
- wData['x/c'] = self.data[:,0]/max(self.data[:,0])
- wData['Cd'] = self.blScalars[0]
- wData['Cdf'] = self.blScalars[1]
- wData['Cdp'] = self.blScalars[2]
- wData['delta*'] = self.blVectors[0]
- wData['theta'] = self.blVectors[1]
- wData['H'] = self.blVectors[2]
- wData['Hk'] = self.blVectors[3]
- wData['H*'] = self.blVectors[4]
- wData['H**'] = self.blVectors[5]
- wData['Cf'] = self.blVectors[6]
- wData['Cdissip'] = self.blVectors[7]
- wData['Ctau'] = self.blVectors[8]
- wData['delta'] = self.blVectors[9]
- if self.xtrF[0] is not None:
- wData['xtrTop'] = self.xtrF[0]
- else:
- wData['xtrTop'] = self.xtr[0]
- if self.xtrF[1] is not None:
- wData['xtrBot'] = self.xtrF[1]
- else:
- wData['xtrBot'] = self.xtr[1]
- self.wData=wData
- return wData
-
- def writeFile(self):
- """Write the results in airfoil_viscous.dat
- """
-
- wData = self.dictionary()
- toW = ['delta*', 'H', 'Cf', 'Ctau'] # list of variable to be written (should be a user input)
- # Write
- print('Writing file: airfoil_viscous.dat...', end = '')
- f = open('airfoil_viscous.dat', 'w+')
-
- f.write('$Aerodynamic coefficients\n')
- f.write(' Re Cd Cdp Cdf xtr_top xtr_bot\n')
- f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(self.Re, wData['Cd'], wData['Cdp'],
- wData['Cdf'], wData['xtrTop'],
- wData['xtrBot']))
- f.write('$Boundary layer variables\n')
- f.write(' x y z x/c')
-
- for s in toW:
- f.write(' {0:>15s}'.format(s))
- f.write('\n')
-
- for i in range(len(self.data[:,0])):
- f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f}'.format(wData['x'][i], wData['y'][i], wData['z'][i], wData['x/c'][i]))
- for s in toW:
- f.write(' {0:15.6f}'.format(wData[s][i]))
- f.write('\n')
-
- f.close()
- print('done.')
-
-
- def run(self, groups):
- """Runs viscous solver driver.
- - Data preprocessing and initialization of the different components.
- - Runs the pseudo-time solver
- - Data postprocessing to ensure proper communication between the solvers.
-
- Args
- ----
-
- - groups (data structure): Structures data containing the necessary information
- related to each group (Airfoil upper surface, airfoil lower surface and wake).
- """
-
- nFactor = self.nFactor
-
- # Initialize computation structures.
-
- dataIn = [None,None] # Merge upper, lower sides and wake
- for n in range(0, len(groups)):
- groups[n].connectListNodes, groups[n].connectListElems, dataIn[n] = groups[n].connectList()
- fMeshDict, cMeshDict, data = GroupConstructor().mergeVectors(dataIn, self.mapMesh, nFactor)
-
- var = Variables(fMeshDict, self.xtrF, self.Ti, self.Re) # Initialize data container for flow config.
-
- if self.it!=0 and self.nbrElmUpper == var.bNodes[1]:
- var.extPar[4::var.nExtPar] = fMeshDict['xxExt'] # Extenal flow local markers.
-
- DirichletBC=Conditions.Boundaries(self.Re) # Boundary condition (Airfoil/Wake).
-
- gradComp=Discretization(var.xx, var.extPar[4::var.nExtPar], var.bNodes) # Gradients computation.
-
- sysSolver = sysMatrix(var.nN, var.nVar, gradComp.fou,gradComp.fou) # Boundary layer equations solver that can provide
- # Residuals and Jacobian computation.
-
- transSolver = Transition(var.xtrF) # Transition solver.
-
- # Initialize time Integrator
-
- tIntegrator = Integration(sysSolver, self.Minfty, self.CFL0)
- tSolver=TSolver(tIntegrator, transSolver, DirichletBC.wakeBC)
-
- # Output setup at the first iteration.
-
- if self.it == 0:
- print('+------------------------------- Solver setup ---------------------------------+')
- print('- Reynolds number: {:>1.2e}'.format(self.Re))
- if self.mapMesh:
- print('- Mesh mapped from {:>5.0f} to {:>5.0f} points.'.format(len(cMeshDict['locCoord']), var.nN))
- print('- Number of elements: {:>5.0f}.'.format(var.nN))
- if var.xtrF[0] is None:
- print('- Free transition on the upper side.')
- else:
- print('- Forced transition on the upper side; x/c = {:<.4f}.'.format(var.xtrF[0]))
- if var.xtrF[1] is None:
- print('- Free transition on the lower side.')
- else:
- print('- Forced transition on the lower side; x/c = {:<.4f}.'.format(var.xtrF[1]))
- print('- Critical amplification ratio: {:<.2f}.'.format(var.Ncrit))
-
- print('')
- print('Numerical parameters')
- print('- CFL0: {:<3.2f}'.format(tSolver.integrator.GetCFL0()))
- print('- Tolerance: {:<3.0f}'.format(np.log10(tSolver.integrator.Gettol())))
- print('- Solver maximum iterations: {:<5.0f}'.format(tSolver.integrator.GetmaxIt()))
- print('+------------------------------------------------------------------------------+')
-
- # Output preprocessing satut.
-
- print('- Mach max: {:<.2f}. IC:'.format(np.max((var.extPar[0::var.nExtPar]))), end = ' ')
-
- # Solver setup.
-
- if self.uPrevious is None or self.nbrElmUpper != var.bNodes[1]: # Typically for the first iteration.
-
- tSolver.initSln = 1 # Flag to use time solver initial condition.
- tSolver.integrator.SetCFL0(1) # Low starting CFL.
- print('Restart.')
- if self.it != 0 and self.nbrElmUpper != var.bNodes[1]:
- print(ccolors.ANSI_BLUE, '- Stagnation point moved.',ccolors.ANSI_RESET)
+ """Driver of the time-dependent boundary layer equation solver (viscous solver)
+ """
+
+ def __init__(self,_Re, user_xtr, _CFL0, _Minfty, _case, _mapMesh, _nFactor):
+ self.Re = _Re
+ self.Minfty = _Minfty
+
+ self.xtrF = user_xtr
+ self.Ti = None
+ self.mapMesh = _mapMesh
+ self.nFactor = _nFactor
+ self.xtr = [1,1]
+ self.CFL0 = _CFL0
+ self.case = _case
+
+ self.Cd = 0
+ self.it = 0
+ self.uPrevious = None
+ self.nbrElmUpper = -1
+ self.stagPtMvmt = 0
+ pass
+
+ def dictionary(self):
+ """Create dictionary with the coordinates and the boundary layer parameters
+ """
+
+ wData = {}
+ wData['x'] = self.data[:,0]
+ wData['y'] = self.data[:,1]
+ wData['z'] = self.data[:,2]
+ wData['x/c'] = self.data[:,0]/max(self.data[:,0])
+ wData['Cd'] = self.blScalars[0]
+ wData['Cdf'] = self.blScalars[1]
+ wData['Cdp'] = self.blScalars[2]
+ wData['delta*'] = self.blVectors[0]
+ wData['theta'] = self.blVectors[1]
+ wData['H'] = self.blVectors[2]
+ wData['Hk'] = self.blVectors[3]
+ wData['H*'] = self.blVectors[4]
+ wData['H**'] = self.blVectors[5]
+ wData['Cf'] = self.blVectors[6]
+ wData['Cdissip'] = self.blVectors[7]
+ wData['Ctau'] = self.blVectors[8]
+ wData['delta'] = self.blVectors[9]
+ if self.xtrF[0] is not None:
+ wData['xtrTop'] = self.xtrF[0]
+ else:
+ wData['xtrTop'] = self.xtr[0]
+ if self.xtrF[1] is not None:
+ wData['xtrBot'] = self.xtrF[1]
+ else:
+ wData['xtrBot'] = self.xtr[1]
+ self.wData = wData
+ return wData
+
+ def writeFile(self):
+ """Write the results in airfoil_viscous.dat
+ """
+
+ wData = self.dictionary()
+ toW = ['delta*', 'H', 'Cf', 'Ctau'] # list of variable to be written (should be a user input)
+ # Write
+ print('Writing file: airfoil_viscous.dat...', end = '')
+ f = open('airfoil_viscous.dat', 'w+')
+
+ f.write('$Aerodynamic coefficients\n')
+ f.write(' Re Cd Cdp Cdf xtr_top xtr_bot\n')
+ f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(self.Re, wData['Cd'], wData['Cdp'],
+ wData['Cdf'], wData['xtrTop'],
+ wData['xtrBot']))
+ f.write('$Boundary layer variables\n')
+ f.write(' x y z x/c')
+
+ for s in toW:
+ f.write(' {0:>15s}'.format(s))
+ f.write('\n')
+
+ for i in range(len(self.data[:,0])):
+ f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f}'.format(wData['x'][i], wData['y'][i], wData['z'][i], wData['x/c'][i]))
+ for s in toW:
+ f.write(' {0:15.6f}'.format(wData[s][i]))
+ f.write('\n')
+
+ f.close()
+ print('done.')
+
+
+ def run(self, groups):
+ """Runs viscous solver driver.
+ - Data preprocessing and initialization of the different components.
+ - Runs the pseudo-time solver
+ - Data postprocessing to ensure proper communication between the solvers.
+
+ Args
+ ----
+
+ - groups (data structure): Structures data containing the necessary information
+ related to each group (Airfoil upper surface, airfoil lower surface and wake).
+ """
+
+ nFactor = self.nFactor
+
+ # Initialize computation structures.
+
+ dataIn = [None,None] # Merge upper, lower sides and wake
+ for n in range(0, len(groups)):
+ groups[n].connectListNodes, groups[n].connectListElems, dataIn[n] = groups[n].connectList()
+ fMeshDict, cMeshDict, data = GroupConstructor().mergeVectors(dataIn, self.mapMesh, nFactor)
+
+ var = Variables(fMeshDict, self.xtrF, self.Ti, self.Re) # Initialize data container for flow config.
+
+ if self.it != 0:
+ var.extPar[4::var.nExtPar] = fMeshDict['xxExt'] # Extenal flow local markers.
+
+ DirichletBC=Conditions.Boundaries(self.Re) # Boundary condition (Airfoil/Wake).
+
+ gradComp=Discretization(var.xx, var.extPar[4::var.nExtPar], var.bNodes) # Gradients computation.
+
+ sysSolver = sysMatrix(var.nN, var.nVar, gradComp.fou,gradComp.fou) # Boundary layer equations solver that can provide
+ # Residuals and Jacobian computation.
+
+ transSolver = Transition(var.xtrF) # Transition solver.
+
+ # Initialize time Integrator
+
+ tIntegrator = Integration(sysSolver, self.Minfty, self.CFL0)
+ tSolver=TSolver(tIntegrator, transSolver, DirichletBC.wakeBC)
+
+ # Output setup at the first iteration.
+
+ if self.it == 0:
+ print('+------------------------------- Solver setup ---------------------------------+')
+ print('- Reynolds number: {:>1.2e}'.format(self.Re))
+ if self.mapMesh:
+ print('- Mesh mapped from {:>5.0f} to {:>5.0f} points.'.format(len(cMeshDict['locCoord']), var.nN))
+ print('- Number of points: {:>5.0f}.'.format(var.nN))
+ if var.xtrF[0] is None:
+ print('- Free transition on the upper side.')
+ else:
+ print('- Forced transition on the upper side; x/c = {:<.4f}.'.format(var.xtrF[0]))
+ if var.xtrF[1] is None:
+ print('- Free transition on the lower side.')
+ else:
+ print('- Forced transition on the lower side; x/c = {:<.4f}.'.format(var.xtrF[1]))
+ print('- Critical amplification ratio: {:<.2f}.'.format(var.Ncrit))
+ if var.xtrF == [0, 0]:
+ print('- Turbulent stagnation point. Ct = 0.03')
+
+ print('')
+ print('Numerical parameters')
+ print('- CFL0: {:<3.2f}'.format(tSolver.integrator.GetCFL0()))
+ print('- Tolerance: {:<3.0f}'.format(np.log10(tSolver.integrator.Gettol())))
+ print('- Solver maximum iterations: {:<5.0f}'.format(tSolver.integrator.GetmaxIt()))
+ print('+------------------------------------------------------------------------------+')
- else: # If we use a solution previously obtained.
+ # Output preprocessing satut.
- var.u = self.uPrevious.copy() # Use previous solution.
+ print('- Mach max: {:<.2f}.'.format(np.max((var.extPar[0::var.nExtPar]))), end = ' ')
- print('Updating.')
+ # Solver setup.
- self.nbrElmUpper = var.bNodes[1]
- DirichletBC.airfoilBC(var) # Reset boundary condition.
-
- # Run the time integration
- #print('+------------------------------- Solver begin ---------------------------------+')
- convergedPoints = tSolver.Run(var)
- #print('+-------------------------------- Solver Exit ---------------------------------+')
+ if self.uPrevious is None or self.nbrElmUpper != var.bNodes[1] or self.stagPtMvmt:
- # Output time solver status.
- if np.any(convergedPoints != 0):
- print(ccolors.ANSI_YELLOW + 'Some point(s) did not converge.', ccolors.ANSI_RESET)
+ tSolver.integrator.itFrozenJac0 = 1
+ tSolver.initSln = 1 # Flag to use time solver initial condition.
+ tSolver.integrator.SetCFL0(1) # Low starting CFL.
+ print('Restart sln.')
- # Save solution to use it as initial condition the next iteration.
+ if self.it != 0 and self.nbrElmUpper != var.bNodes[1]:
- self.uPrevious=var.u.copy() # Copy solution.
+ self.stagPtMvmt = 1
+ print(ccolors.ANSI_BLUE, '- Stagnation point moved.',ccolors.ANSI_RESET)
- # Compute blowing velocities and sort the boundary layer parameters.
+ else:
- self.blScalars, blwVelUp, blwVelLw, blwVelWk, AFblVectors, WKblVectors, AFdeltaStarOut, WKdeltaStarOut = PostProcess().sortParam(var, self.mapMesh, cMeshDict, nFactor)
-
- self.xtr = var.xtr.copy()
- self.blVectors = AFblVectors
+ self.stagPtMvmt = 0
- AFblwVel = np.concatenate((blwVelUp, blwVelLw))
+ else: # If we use a solution previously obtained.
- # Group modification to ensure communication between the solvers.
+ var.u = self.uPrevious.copy() # Use previous solution.
- groups[0].TEnd = [AFblVectors[1][var.bNodes[1]-1], AFblVectors[1][var.bNodes[2]-1]]
- groups[0].HEnd = [AFblVectors[2][var.bNodes[1]-1], AFblVectors[2][var.bNodes[2]-1]]
- groups[0].CtEnd = [AFblVectors[8][var.bNodes[1]-1], AFblVectors[8][var.bNodes[2]-1]]
- groups[0].deltaStar = AFdeltaStarOut
- groups[0].xx = cMeshDict['locCoord'][:cMeshDict['bNodes'][2]]
+ print('Updating sln.')
- # Sort the following results in reference frame.
- groups[0].deltaStar = groups[0].deltaStar[groups[0].connectListNodes.argsort()]
- groups[0].xx = groups[0].xx[groups[0].connectListNodes.argsort()]
- groups[0].u = AFblwVel[groups[0].connectListElems.argsort()]
- self.data = data[:var.bNodes[2],:]
+ self.nbrElmUpper = var.bNodes[1]
+ DirichletBC.airfoilBC(var) # Reset boundary condition.
- # Drag is measured at the end of the wake.
- self.Cd = self.blScalars[0]
- self.Cdf = self.blScalars[1]
- self.Cdp = self.blScalars[2] # Cdf comes from the airfoil as there is no friction in the wake.
- groups[1].deltaStar = WKdeltaStarOut
- groups[1].xx = cMeshDict['locCoord'][cMeshDict['bNodes'][2]:]
+ # Run the time integration
+ convergedPoints = tSolver.Run(var)
+ # Output time solver status.
+ if np.any(convergedPoints != 0):
+ print(ccolors.ANSI_YELLOW + 'Some point(s) did not converge.', ccolors.ANSI_RESET)
- # Sort the following results in reference frame.
- groups[1].deltaStar = groups[1].deltaStar[groups[1].connectListNodes.argsort()]
- groups[1].xx = groups[1].xx[groups[1].connectListNodes.argsort()]
- groups[1].u = blwVelWk[groups[1].connectListElems.argsort()]
+ self.uPrevious=var.u.copy() # Save solution to use it as initial condition the next iteration.
- pass
\ No newline at end of file
+ # Compute blowing velocities and sort the boundary layer parameters.
+
+ self.blScalars, blwVelUp, blwVelLw, blwVelWk, AFblVectors, WKblVectors, AFdeltaStarOut, WKdeltaStarOut = PostProcess().sortParam(var, self.mapMesh, cMeshDict, nFactor)
+
+ self.xtr = var.xtr.copy()
+ self.blVectors = AFblVectors
+
+ AFblwVel = np.concatenate((blwVelUp, blwVelLw))
+
+ # Group modification to ensure communication between the solvers.
+
+ groups[0].TEnd = [AFblVectors[1][var.bNodes[1]-1], AFblVectors[1][var.bNodes[2]-1]]
+ groups[0].HEnd = [AFblVectors[2][var.bNodes[1]-1], AFblVectors[2][var.bNodes[2]-1]]
+ groups[0].CtEnd = [AFblVectors[8][var.bNodes[1]-1], AFblVectors[8][var.bNodes[2]-1]]
+ groups[0].deltaStar = AFdeltaStarOut
+ groups[0].xx = cMeshDict['locCoord'][:cMeshDict['bNodes'][2]]
+
+ # Sort the following results in reference frame.
+ groups[0].deltaStar = groups[0].deltaStar[groups[0].connectListNodes.argsort()]
+ groups[0].xx = groups[0].xx[groups[0].connectListNodes.argsort()]
+ groups[0].u = AFblwVel[groups[0].connectListElems.argsort()]
+ self.data = data[:var.bNodes[2],:]
+
+ # Drag is measured at the end of the wake.
+ self.Cd = self.blScalars[0]
+ self.Cdf = self.blScalars[1]
+ self.Cdp = self.blScalars[2] # Cdf comes from the airfoil as there is no friction in the wake.
+ groups[1].deltaStar = WKdeltaStarOut
+ groups[1].xx = cMeshDict['locCoord'][cMeshDict['bNodes'][2]:]
+
+ # Sort the following results in reference frame.
+ groups[1].deltaStar = groups[1].deltaStar[groups[1].connectListNodes.argsort()]
+ groups[1].xx = groups[1].xx[groups[1].connectListNodes.argsort()]
+ groups[1].u = blwVelWk[groups[1].connectListElems.argsort()]
+
+
+ """if self.it >= 24:
+ self.writeFile()
+ Validation().main(1,self.wData)"""
+
+ pass
\ No newline at end of file
diff --git a/dart/viscous/Drivers/DGroupConstructor.py b/dart/viscous/Drivers/DGroupConstructor.py
index 281a60aecb07056ea56fd9639fe0aef96de6f7df..44877c7d62e24180c383d8319156baa755c36ed9 100755
--- a/dart/viscous/Drivers/DGroupConstructor.py
+++ b/dart/viscous/Drivers/DGroupConstructor.py
@@ -78,11 +78,12 @@ class GroupConstructor():
fdStarExt_up = self.interpolateToFineMesh(dataIn[0][0][:,7], fMeshUp, nFactor)
fdStarExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,7], fMeshLw, nFactor)
fdStarExt_wk = self.interpolateToFineMesh(dataIn[1][:,7] , fMeshWk, nFactor)
- fdStarext = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk))
fxxExt_up = self.interpolateToFineMesh(dataIn[0][0][:,8], fMeshUp, nFactor)
fxxExt_lw = self.interpolateToFineMesh(dataIn[0][1][:,8], fMeshLw, nFactor)
fxxExt_wk = self.interpolateToFineMesh(dataIn[1][:,8] , fMeshWk, nFactor)
+
+ fdStarext = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk))
fxxext = np.concatenate((fxxExt_up, fxxExt_lw[1:], fxxExt_wk))
fdv = np.concatenate((dv_up, dv_lw[1:], dv_wk))
diff --git a/dart/viscous/Drivers/DPostProcess.py b/dart/viscous/Drivers/DPostProcess.py
index 4649568cd2dca6f514863f6770cf443e74327ec8..bc49d67c2f01c86ea17fd309e90cfc26a3c8d8fc 100644
--- a/dart/viscous/Drivers/DPostProcess.py
+++ b/dart/viscous/Drivers/DPostProcess.py
@@ -1,208 +1,208 @@
import numpy as np
class PostProcess:
- def __init__(self):
- pass
+ def __init__(self):
+ pass
- def sortParam(self,var, mapMesh, cMeshDict, nFactor):
+ def sortParam(self,var, mapMesh, cMeshDict, nFactor):
- # Post process delta star
+ # Post process delta star
- # Recompute deltaStar.
- var.blPar[9::var.nBlPar] = var.u[0::var.nVar] * var.u[1::var.nVar]
+ # Recompute deltaStar.
+ var.blPar[9::var.nBlPar] = var.u[0::var.nVar] * var.u[1::var.nVar]
- # Map delta star to the coarse mesh and divide between the airfoil
- # and wake related points
- if mapMesh:
+ # Map delta star to the coarse mesh and divide between the airfoil
+ # and wake related points
+ if mapMesh:
- deltaStar = self.inverseMeshMapping(var.blPar[9::var.nBlPar], var.bNodes, cMeshDict['globCoord'], cMeshDict['bNodes'], nFactor)
+ deltaStar = self.inverseMeshMapping(var.blPar[9::var.nBlPar], var.bNodes, cMeshDict['globCoord'], cMeshDict['bNodes'], nFactor)
- AFdeltaStarOut = deltaStar[:cMeshDict['bNodes'][2]]
- WKdeltaStarOut = deltaStar[cMeshDict['bNodes'][2]:]
+ AFdeltaStarOut = deltaStar[:cMeshDict['bNodes'][2]]
+ WKdeltaStarOut = deltaStar[cMeshDict['bNodes'][2]:]
- else:
-
- AFdeltaStarOut = var.blPar[9:var.bNodes[2]*var.nBlPar:var.nBlPar]
- WKdeltaStarOut = var.blPar[var.bNodes[2]*var.nBlPar + 9::var.nBlPar]
+ else:
+
+ AFdeltaStarOut = var.blPar[9:var.bNodes[2]*var.nBlPar:var.nBlPar]
+ WKdeltaStarOut = var.blPar[var.bNodes[2]*var.nBlPar + 9::var.nBlPar]
- # Compute blowing velocity on each cell.
+ # Compute blowing velocity on each cell.
- blwVelUp = np.zeros(var.bNodes[1] - 1)
- blwVelLw = np.zeros(var.bNodes[2] - var.bNodes[1])
- blwVelWk = np.zeros(var.nN - var.bNodes[2] - 1)
+ blwVelUp = np.zeros(var.bNodes[1] - 1)
+ blwVelLw = np.zeros(var.bNodes[2] - var.bNodes[1])
+ blwVelWk = np.zeros(var.nN - var.bNodes[2] - 1)
- # Upper
- i = 0
- for iPoint in range(1, var.bNodes[1]):
+ # Upper
+ i = 0
+ for iPoint in range(1, var.bNodes[1]):
- iPrev = iPoint - 1
+ iPrev = iPoint - 1
- blwVelUp[i] = (var.extPar[iPoint*var.nExtPar+1] * var.u[iPoint*var.nVar+3] * var.blPar[iPoint*var.nBlPar+9]
- - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[iPoint*var.nExtPar+1] * (var.xx[iPoint] - var.xx[iPrev]))
+ blwVelUp[i] = (var.extPar[iPoint*var.nExtPar+1] * var.u[iPoint*var.nVar+3] * var.blPar[iPoint*var.nBlPar+9]
+ - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[iPoint*var.nExtPar+1] * (var.xx[iPoint] - var.xx[iPrev]))
- i += 1
+ i += 1
- # Lower
- i = 0
- for iPoint in range(var.bNodes[1], var.bNodes[2]):
+ # Lower
+ i = 0
+ for iPoint in range(var.bNodes[1], var.bNodes[2]):
- iPrev = 0 if iPoint == var.bNodes[1] else iPoint - 1
+ iPrev = 0 if iPoint == var.bNodes[1] else iPoint - 1
- blwVelLw[i] = (var.extPar[iPoint*var.nExtPar+1] * var.u[iPoint*var.nVar+3] * var.blPar[iPoint*var.nBlPar+9]
- - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[iPoint*var.nExtPar+1] * (var.xx[iPoint] - var.xx[iPrev]))
+ blwVelLw[i] = (var.extPar[iPoint*var.nExtPar+1] * var.u[iPoint*var.nVar+3] * var.blPar[iPoint*var.nBlPar+9]
+ - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[iPoint*var.nExtPar+1] * (var.xx[iPoint] - var.xx[iPrev]))
- i += 1
+ i += 1
- # Wake
- i = 0
- for iPoint in range(var.bNodes[2] + 1, var.nN):
+ # Wake
+ i = 0
+ for iPoint in range(var.bNodes[2] + 1, var.nN):
- iPrev = iPoint - 1
+ iPrev = iPoint - 1
- blwVelWk[i] = (var.extPar[iPoint*var.nExtPar+1] * var.u[iPoint*var.nVar+3] * var.blPar[iPoint*var.nBlPar+9]
- - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[iPoint*var.nExtPar+1] * (var.xx[iPoint] - var.xx[iPrev]))
+ blwVelWk[i] = (var.extPar[iPoint*var.nExtPar+1] * var.u[iPoint*var.nVar+3] * var.blPar[iPoint*var.nBlPar+9]
+ - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[iPoint*var.nExtPar+1] * (var.xx[iPoint] - var.xx[iPrev]))
- i += 1
+ i += 1
- if mapMesh:
- blwVelUp, blwVelLw, blwVelWk = self.mapBlowingVelocities(blwVelUp, blwVelLw, blwVelWk, var.bNodes, cMeshDict['globCoord'], cMeshDict['bNodes'], var.xGlobal, nFactor)
+ if mapMesh:
+ blwVelUp, blwVelLw, blwVelWk = self.mapBlowingVelocities(blwVelUp, blwVelLw, blwVelWk, var.bNodes, cMeshDict['globCoord'], cMeshDict['bNodes'], var.xGlobal, nFactor)
- # Postprocess the boundary layer parameters.
+ # Postprocess the boundary layer parameters.
- # Normalize Friction and dissipation coefficients.
- frictionCoeff=var.u[3::var.nVar]**2 * var.blPar[2::var.nBlPar]
- dissipCoeff=var.u[3::var.nVar]**2 * var.blPar[1::var.nBlPar]
+ # Normalize friction and dissipation coefficients.
+ frictionCoeff=var.u[3::var.nVar]**2 * var.blPar[2::var.nBlPar]
+ dissipCoeff=var.u[3::var.nVar]**2 * var.blPar[1::var.nBlPar]
- # Compute total drag coefficient.
- CdTot=2*var.u[-5]*(var.u[-2]**((var.u[-4]+5)/2))
+ # Compute total drag coefficient.
+ CdTot=2*var.u[-5]*(var.u[-2]**((var.u[-4]+5)/2))
- # Compute friction and pressure drag coefficients.
- Cdf_upper=np.trapz(frictionCoeff[:var.bNodes[1]],var.xGlobal[:var.bNodes[1]])
- # Insert stagnation point in the lower side.
- Cdf_lower=np.trapz(np.insert(frictionCoeff[var.bNodes[1]:var.bNodes[2]],0,frictionCoeff[0]),
- np.insert(var.xGlobal[var.bNodes[1]:var.bNodes[2]],0,var.xGlobal[0]))
- Cdf=Cdf_upper+Cdf_lower # No friction in the wake
- Cdp=CdTot-Cdf
+ # Compute friction and pressure drag coefficients.
+ Cdf_upper=np.trapz(frictionCoeff[:var.bNodes[1]],var.xGlobal[:var.bNodes[1]])
+ # Insert stagnation point in the lower side.
+ Cdf_lower=np.trapz(np.insert(frictionCoeff[var.bNodes[1]:var.bNodes[2]],0,frictionCoeff[0]),
+ np.insert(var.xGlobal[var.bNodes[1]:var.bNodes[2]],0,var.xGlobal[0]))
+ Cdf=Cdf_upper+Cdf_lower # No friction in the wake
+ Cdp=CdTot-Cdf
- # [CdTot, Cdf, Cdp]
- blScalars=[CdTot, Cdf, Cdp]
+ # [CdTot, Cdf, Cdp]
+ blScalars=[CdTot, Cdf, Cdp]
- # Store boundary layer final parameters in lists.
- # [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta]
- blVectors_airfoil=[var.blPar[9:var.bNodes[2]*var.nBlPar:var.nBlPar],
- var.u[0:var.bNodes[2]*var.nVar:var.nVar],
- var.u[1:var.bNodes[2]*var.nVar:var.nVar],
- var.blPar[6:var.bNodes[2]*var.nBlPar:var.nBlPar],
- var.blPar[4:var.bNodes[2]*var.nBlPar:var.nBlPar],
- var.blPar[5:var.bNodes[2]*var.nBlPar:var.nBlPar],
- frictionCoeff[:var.bNodes[2]], dissipCoeff[:var.bNodes[2]],
- var.u[4:var.bNodes[2]*var.nVar:var.nVar],
- var.blPar[3:var.bNodes[2]*var.nBlPar:var.nBlPar]]
+ # Store boundary layer final parameters in lists.
+ # [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta]
+ blVectors_airfoil=[var.blPar[9:var.bNodes[2]*var.nBlPar:var.nBlPar],
+ var.u[0:var.bNodes[2]*var.nVar:var.nVar],
+ var.u[1:var.bNodes[2]*var.nVar:var.nVar],
+ var.blPar[6:var.bNodes[2]*var.nBlPar:var.nBlPar],
+ var.blPar[4:var.bNodes[2]*var.nBlPar:var.nBlPar],
+ var.blPar[5:var.bNodes[2]*var.nBlPar:var.nBlPar],
+ frictionCoeff[:var.bNodes[2]], dissipCoeff[:var.bNodes[2]],
+ var.u[4:var.bNodes[2]*var.nVar:var.nVar],
+ var.blPar[3:var.bNodes[2]*var.nBlPar:var.nBlPar]]
- blVectors_wake=[var.blPar[var.bNodes[2]*var.nBlPar+9::var.nBlPar],
- var.u[var.bNodes[2]*var.nVar+0::var.nVar],
- var.u[var.bNodes[2]*var.nVar+1::var.nVar],
- var.blPar[var.bNodes[2]*var.nBlPar+6::var.nBlPar],
- var.blPar[var.bNodes[2]*var.nBlPar+4::var.nBlPar],
- var.blPar[var.bNodes[2]*var.nBlPar+5::var.nBlPar],
- frictionCoeff[var.bNodes[2]:], dissipCoeff[var.bNodes[2]:],
- var.u[var.bNodes[2]*var.nVar+4::var.nVar],
- var.blPar[var.bNodes[2]*var.nBlPar+3::var.nBlPar]]
+ blVectors_wake=[var.blPar[var.bNodes[2]*var.nBlPar+9::var.nBlPar],
+ var.u[var.bNodes[2]*var.nVar+0::var.nVar],
+ var.u[var.bNodes[2]*var.nVar+1::var.nVar],
+ var.blPar[var.bNodes[2]*var.nBlPar+6::var.nBlPar],
+ var.blPar[var.bNodes[2]*var.nBlPar+4::var.nBlPar],
+ var.blPar[var.bNodes[2]*var.nBlPar+5::var.nBlPar],
+ frictionCoeff[var.bNodes[2]:], dissipCoeff[var.bNodes[2]:],
+ var.u[var.bNodes[2]*var.nVar+4::var.nVar],
+ var.blPar[var.bNodes[2]*var.nBlPar+3::var.nBlPar]]
- return blScalars, blwVelUp, blwVelLw, blwVelWk, blVectors_airfoil, blVectors_wake, AFdeltaStarOut, WKdeltaStarOut
+ return blScalars, blwVelUp, blwVelLw, blwVelWk, blVectors_airfoil, blVectors_wake, AFdeltaStarOut, WKdeltaStarOut
- def mapBlowingVelocities(self, blwVelUp, blwVelLw, blwVelWk, fbNodes, cMesh, cMeshbNodes, fMesh, nFactor):
- """Maps blowing velocties from cell to cell using weighted average interpolation.
+ def mapBlowingVelocities(self, blwVelUp, blwVelLw, blwVelWk, fbNodes, cMesh, cMeshbNodes, fMesh, nFactor):
+ """Maps blowing velocties from cell to cell using weighted average interpolation.
- Args
- ----
+ Args
+ ----
- - blwVelUp (numpy.ndarray): Blowing velocities on the fine mesh of the upper side.
+ - blwVelUp (numpy.ndarray): Blowing velocities on the fine mesh of the upper side.
- - blwVelLw (numpy.ndarray): Blowing velocities on the fine mesh of the lower side.
-
- - blwVelWk (numpy.ndarray): Blowing velocities on the fine mesh of the wake.
+ - blwVelLw (numpy.ndarray): Blowing velocities on the fine mesh of the lower side.
+
+ - blwVelWk (numpy.ndarray): Blowing velocities on the fine mesh of the wake.
- - fbNodes (numpy.ndarray): Boundary nodes of the fine mesh.
+ - fbNodes (numpy.ndarray): Boundary nodes of the fine mesh.
- - cMesh (numpy.ndarray): Coarse mesh coordinates in the global frame of reference.
+ - cMesh (numpy.ndarray): Coarse mesh coordinates in the global frame of reference.
- - cMeshbNodes (numpy.ndarray): Boundary nodes of the coarse mesh.
+ - cMeshbNodes (numpy.ndarray): Boundary nodes of the coarse mesh.
- - fMesh (numpy.ndarray): Fine mesh coordinates in the global frame of reference.
+ - fMesh (numpy.ndarray): Fine mesh coordinates in the global frame of reference.
- - nFactor (int): Multiplicator underlying mesh adaptation.
+ - nFactor (int): Multiplicator underlying mesh adaptation.
- Example
- -------
+ Example
+ -------
- Fine mesh: |-----------|-----------|-----------|----------|
- \ /
- (1/2) \ / (1/2)
- \ /
- \ /
- Coarse mesh: |-----------------------|----------------------|
- """
+ Fine mesh: |-----------|-----------|-----------|----------|
+ \ /
+ (1/2) \ / (1/2)
+ \ /
+ \ /
+ Coarse mesh: |-----------------------|----------------------|
+ """
- cblwVelUp = np.empty(len(cMesh[:cMeshbNodes[1]]) - 1)
- cblwVelLw = np.empty(len(cMesh[cMeshbNodes[1]:cMeshbNodes[2]]))
- cWKblwVel = np.empty(len(cMesh[cMeshbNodes[2]:]) - 1)
+ cblwVelUp = np.empty(len(cMesh[:cMeshbNodes[1]]) - 1)
+ cblwVelLw = np.empty(len(cMesh[cMeshbNodes[1]:cMeshbNodes[2]]))
+ cWKblwVel = np.empty(len(cMesh[cMeshbNodes[2]:]) - 1)
- for cCell in range(len(cblwVelUp)):
+ for cCell in range(len(cblwVelUp)):
- cblwVelUp[cCell] = np.mean(blwVelUp[cCell*nFactor:(cCell + 1)*nFactor])
-
- for cCell in range(len(cblwVelLw)):
+ cblwVelUp[cCell] = np.mean(blwVelUp[cCell*nFactor:(cCell + 1)*nFactor])
+
+ for cCell in range(len(cblwVelLw)):
- cblwVelLw[cCell] = np.mean(blwVelLw[cCell*nFactor:(cCell + 1)*nFactor])
+ cblwVelLw[cCell] = np.mean(blwVelLw[cCell*nFactor:(cCell + 1)*nFactor])
- for cCell in range(len(cWKblwVel)):
+ for cCell in range(len(cWKblwVel)):
- cWKblwVel[cCell] = np.mean(blwVelWk[cCell*nFactor:(cCell + 1)*nFactor])
+ cWKblwVel[cCell] = np.mean(blwVelWk[cCell*nFactor:(cCell + 1)*nFactor])
- return cblwVelUp, cblwVelLw, cWKblwVel
+ return cblwVelUp, cblwVelLw, cWKblwVel
- def inverseMeshMapping(self, fVector, fbNodes, cMesh, cbNodes, nFactor):
- """Maps the variables required for blowing velocity computation from
- the fine mesh to the coarse mesh. Group separation is performed
- to avoid confusion between the points at the same place on the mesh
- but with different physical meaning (e.g First wake point).
+ def inverseMeshMapping(self, fVector, fbNodes, cMesh, cbNodes, nFactor):
+ """Maps the variables required for blowing velocity computation from
+ the fine mesh to the coarse mesh. Group separation is performed
+ to avoid confusion between the points at the same place on the mesh
+ but with different physical meaning (e.g First wake point).
- The mapping is based on a weigthed average with the nearest neighbors.
+ The mapping is based on a weigthed average with the nearest neighbors.
- Args
- ----
+ Args
+ ----
- - fVector (numpy.ndarray): Vector containing the function on the fine mesh.
+ - fVector (numpy.ndarray): Vector containing the function on the fine mesh.
- - fbNodes (numpy.ndarray): Id of the boundary nodes on the fine mesh.
+ - fbNodes (numpy.ndarray): Id of the boundary nodes on the fine mesh.
- - cMesh (numpy.ndarray): Coarse mesh.
+ - cMesh (numpy.ndarray): Coarse mesh.
- - cbNodes (numpy.ndarray): Id of the boundary nodes on the coarse mesh.
+ - cbNodes (numpy.ndarray): Id of the boundary nodes on the coarse mesh.
- - nFactor (int): Multiplicator underlying mesh adaptation.
+ - nFactor (int): Multiplicator underlying mesh adaptation.
- Return
- ------
+ Return
+ ------
- - cVector (numpy.ndarray): Function contained in 'fVector' mapped to the coarse mesh.
+ - cVector (numpy.ndarray): Function contained in 'fVector' mapped to the coarse mesh.
- Example
- -------
+ Example
+ -------
- Fine mesh: |-----------|-----------|-----------|----------|
- \ | /
- (1/4) \ | (1/2) / (1/4)
- \ V /
- -----> <-----
- Coarse mesh: |-----------------------|----------------------|
- """
+ Fine mesh: |-----------|-----------|-----------|----------|
+ \ | /
+ (1/4) \ | (1/2) / (1/4)
+ \ V /
+ -----> <-----
+ Coarse mesh: |-----------------------|----------------------|
+ """
- cVector = np.zeros(len(cMesh))
+ cVector = np.zeros(len(cMesh))
- for cPoint in range(len(cMesh)-1):
- cVector[cPoint] = fVector[cPoint * nFactor]
- cVector[-1] = fVector[-1]
- return cVector
+ for cPoint in range(len(cMesh)-1):
+ cVector[cPoint] = fVector[cPoint * nFactor]
+ cVector[-1] = fVector[-1]
+ return cVector
diff --git a/dart/viscous/Master/DAirfoil.py b/dart/viscous/Master/DAirfoil.py
index bf1431f50f62accb19598b66a6f82c4c0820dc69..6773c30a7e2a1b5ff1fde5da2da28db1bfec85e0 100755
--- a/dart/viscous/Master/DAirfoil.py
+++ b/dart/viscous/Master/DAirfoil.py
@@ -74,6 +74,11 @@ class Airfoil(Boundary):
eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
# Find the stagnation point
+ """if it == 0:
+ self.idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
+ idxStag = self.idxStag
+ else:
+ idxStag = self.idxStag"""
idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
# Find TE nodes (assuming suction side element is above pressure side element in the y direction)
diff --git a/dart/viscous/Master/DCoupler.py b/dart/viscous/Master/DCoupler.py
index 91c383198ce2fe862cb8454d1b51de7b930c7900..4d83869c771a6eec9dbaf8dca1f3a9adbbe7eefa 100755
--- a/dart/viscous/Master/DCoupler.py
+++ b/dart/viscous/Master/DCoupler.py
@@ -43,6 +43,9 @@ class Coupler:
"""Viscous inviscid coupling.
"""
+ # Save history.
+ cvgMonitoring = []
+
# Initialize loop
it = 0
converged = 0 # temp
@@ -55,6 +58,7 @@ class Coupler:
# Run inviscid solver
print('+----------------------------- Inviscid solver --------------------------------+')
+ self.isolver.reset()
self.isolver.run()
for n in range(0, len(self.group)):
@@ -72,6 +76,7 @@ class Coupler:
# Write inviscid pressure distribution on the first iteration
Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
tboxU.write(Cp, 'Cpinv_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+
# Run viscous solver
print('+------------------------------ Viscous solver --------------------------------+')
self.vsolver.run(self.group)
@@ -80,16 +85,21 @@ class Coupler:
for i in range(0, self.group[n].nE):
self.group[n].boundary.setU(i, self.group[n].u[i])
+ error = abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd
+ cvgMonitoring.append(error)
# Check convergence and output coupling iteration.
- error = abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd
-
if error <= self.tol:
print('')
print(ccolors.ANSI_GREEN, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} < {:<2.3f}'
.format(it, self.vsolver.Cd, np.log10(error), np.log10(self.tol)), ccolors.ANSI_RESET)
print('')
+ """plt.plot(cvgMonitoring[:it+1])
+ plt.yscale('log')
+ plt.xlabel('Number of iterations (-)')
+ plt.ylabel('log[Res] (-)')
+ plt.show()"""
converged = 1
else:
@@ -106,13 +116,8 @@ class Coupler:
# Save results
- # Inviscid solver files.
- self.isolver.save(self.writer)
-
- # Viscous solver files.
- self.vsolver.writeFile()
-
- # Pressure coefficient.
- Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
+ self.isolver.save(self.writer) # Inviscid solver files.
+ self.vsolver.writeFile() # Viscous solver files.
+ Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp) # Pressure coefficient.
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
print('\n')
diff --git a/dart/viscous/Physics/DBIConditions.py b/dart/viscous/Physics/DBIConditions.py
index c1a145b9127aa156d64104d976e0f187ef75d598..1f2d5c5777b02fa6a439253b6a5d373b829c0203 100755
--- a/dart/viscous/Physics/DBIConditions.py
+++ b/dart/viscous/Physics/DBIConditions.py
@@ -123,77 +123,83 @@ class Boundaries:
self.Re = _Re
- # Airfoil boundary conditions
- def airfoilBC(self, var):
+ def airfoilBC(self, cfgFlow):
"""Boundary conditions at the stagnation point.
Twaithes solution method values used for the
stagnation point.
"""
- if var.dv[0] > 0:
+ if cfgFlow.dv[0] > 0:
+ T0 = np.sqrt(0.075 / (self.Re * cfgFlow.dv[0]))
+ else:
+ T0 = 1e-6
- T0 = np.sqrt(0.075 / (self.Re * var.dv[0]))
+ H0 = 2.23 # One parameter family Falkner Skan.
+ n0 = 0 # Initial log of amplification factor.
+ ue0 = cfgFlow.extPar[2] # Inviscid flow velocity.
+ Ct0 = 0 # Stagnation point is laminar.
- else:
+ cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar + 0] = ue0*T0*cfgFlow.Re # Reynolds number based on the momentum thickness
- T0 = 1e-6
+ if cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar+0] < 1:
+ cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar+0] = 1
+ ue0 = cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar + 0] / (T0*cfgFlow.Re)
- H0 = 2.23 # One parameter family Falkner Skan.
- n0 = 0 # Initial log of amplification factor.
- ue0=var.extPar[2] # Inviscid flow velocity.
- Ct0 = 0 # Stagnation point is laminar.
+ cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar + 9] = T0 * H0 # deltaStar
- # Compute Reynolds number based on the momentum thickness
- var.blPar[var.bNodes[0]*var.nBlPar + 0] = ue0*T0*var.Re
+ if cfgFlow.xtrF == [0, 0]:
+ # Turbulent closures.
- if var.blPar[var.bNodes[0]*var.nBlPar+0] < 1:
- var.blPar[var.bNodes[0]*var.nBlPar+0] = 1
- ue0 = var.blPar[var.bNodes[0]*var.nBlPar + 0] / (T0*var.Re)
+ Ct0 = 0.03
+ cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar+1:cfgFlow.bNodes[0]*cfgFlow.nBlPar+9] = Closures.turbulentClosure(T0, H0,
+ Ct0, cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar+0],
+ cfgFlow.extPar[cfgFlow.bNodes[0]*cfgFlow.nExtPar+0], 'wake')
+ else:
+ # Laminar closures.
- var.blPar[var.bNodes[0]*var.nBlPar + 9] = T0 * H0 # deltaStar
- var.blPar[(var.bNodes[0]*var.nBlPar + 1):(var.bNodes[0]*var.nBlPar+7)]=Closures.laminarClosure(T0, H0,
- var.blPar[var.bNodes[0]*var.nBlPar+0],
- var.extPar[var.bNodes[0]*var.nExtPar+0],
- 'airfoil')
- #var.blPar[var.bNodes[0]*var.nBlPar+8]=(var.blPar[var.bNodes[0]*var.nBlPar+4]/2)*(1-4*(var.blPar[var.bNodes[0]*var.nBlPar+6]-1)/(3*H0)) # Equivalent normalized wall slip velocity
+ cfgFlow.blPar[(cfgFlow.bNodes[0]*cfgFlow.nBlPar + 1):(cfgFlow.bNodes[0]*cfgFlow.nBlPar+7)]=Closures.laminarClosure(T0, H0,
+ cfgFlow.blPar[cfgFlow.bNodes[0]*cfgFlow.nBlPar+0],
+ cfgFlow.extPar[cfgFlow.bNodes[0]*cfgFlow.nExtPar+0],
+ 'airfoil')
# Set boundary condition.
- var.u[var.bNodes[0]*var.nVar:(var.bNodes[0]+1)*var.nVar] = [T0, H0, n0, ue0, Ct0]
+
+ cfgFlow.u[cfgFlow.bNodes[0]*cfgFlow.nVar:(cfgFlow.bNodes[0]+1)*cfgFlow.nVar] = [T0, H0, n0, ue0, Ct0]
pass
# Wake boundary conditions
- def wakeBC(self, var):
+ def wakeBC(self, cfgFlow):
"""Boundary conditions at the first point of the wake
based on the parameters at the last point on the upper and lower sides.
Drela (1989).
"""
- xEnd_up=var.u[(var.bNodes[1]-1)*var.nVar:var.bNodes[1]*var.nVar]
- xEnd_lw=var.u[(var.bNodes[2]-1)*var.nVar:var.bNodes[2]*var.nVar]
+ # Extract last points on the upper and lower sides
- T0=xEnd_up[0]+xEnd_lw[0] # (Theta_up+Theta_low).
- H0=(xEnd_up[0]*xEnd_up[1]+xEnd_lw[0]*xEnd_lw[1])/T0 # ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low).
- n0=9
- ue0=var.extPar[var.bNodes[2]*var.nExtPar+2] # Inviscid velocity.
- Ct0=(xEnd_up[0]*xEnd_up[4]+xEnd_lw[0]*xEnd_lw[4])/T0 # ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low).
+ xEnd_up = cfgFlow.u[(cfgFlow.bNodes[1]-1)*cfgFlow.nVar:cfgFlow.bNodes[1]*cfgFlow.nVar]
+ xEnd_lw = cfgFlow.u[(cfgFlow.bNodes[2]-1)*cfgFlow.nVar:cfgFlow.bNodes[2]*cfgFlow.nVar]
- # Reynolds number based on the momentum thickness.
- var.blPar[var.bNodes[2]*var.nBlPar+0] = ue0*T0*var.Re
+ T0 = xEnd_up[0]+xEnd_lw[0] # (Theta_up+Theta_low).
+ H0 = (xEnd_up[0]*xEnd_up[1]+xEnd_lw[0]*xEnd_lw[1])/T0 # ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low).
+ n0 = 9
+ ue0 = cfgFlow.extPar[cfgFlow.bNodes[2]*cfgFlow.nExtPar+2] # Inviscid velocity.
+ Ct0 = (xEnd_up[0]*xEnd_up[4]+xEnd_lw[0]*xEnd_lw[4])/T0 # ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low).
+
+ cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+0] = ue0*T0*cfgFlow.Re # Reynolds number based on the momentum thickness.
- if var.blPar[var.bNodes[2]*var.nBlPar+0] < 1:
- var.blPar[var.bNodes[2]*var.nBlPar+0] = 1
- ue0 = var.blPar[var.bNodes[2]*var.nBlPar+0] / (T0*var.Re)
-
- var.blPar[var.bNodes[2]*var.nBlPar+9]=T0*H0 # deltaStar
-
- # Laminar reset
- var.blPar[(var.bNodes[2]*var.nBlPar+1):(var.bNodes[2]*var.nBlPar+7)]=Closures.laminarClosure(T0, H0,
- var.blPar[var.bNodes[2]*var.nBlPar+0],
- var.extPar[var.bNodes[2]*var.nExtPar+0],
- 'wake')
- #var.blPar[var.bNodes[2]*var.nBlPar+7]=0 # Cteq stays = 0
- #var.blPar[var.bNodes[2]*var.nBlPar+8]=(var.blPar[var.bNodes[2]*var.nBlPar+4]/2)*(1-4*(var.blPar[var.bNodes[2]*var.nBlPar+6]-1)/(3*H0)) if H0!=0 else 0 # Equivalent normalized wall slip velocity
+ if cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+0] < 1:
+ cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+0] = 1
+ ue0 = cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+0] / (T0*cfgFlow.Re)
+
+ cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+9] = T0*H0 # deltaStar
+ # Turbulent closures.
+
+ cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+1:cfgFlow.bNodes[2]*cfgFlow.nBlPar+9] = Closures.turbulentClosure(T0, H0,
+ Ct0, cfgFlow.blPar[cfgFlow.bNodes[2]*cfgFlow.nBlPar+0],
+ cfgFlow.extPar[cfgFlow.bNodes[2]*cfgFlow.nExtPar+0], 'wake')
+
# Set boundary condition.
- var.u[var.bNodes[2]*var.nVar:(var.bNodes[2]+1)*var.nVar] = [T0, H0, n0, ue0, Ct0]
- pass
+
+ cfgFlow.u[cfgFlow.bNodes[2]*cfgFlow.nVar:(cfgFlow.bNodes[2]+1)*cfgFlow.nVar] = [T0, H0, n0, ue0, Ct0]
+ pass
\ No newline at end of file
diff --git a/dart/viscous/Physics/DClosures.py b/dart/viscous/Physics/DClosures.py
index 1acad305c0fc9e3053338ac73b4dc5a49e9336a1..4ffcd75dbfde6d3e1d8ce12e2846b4fd500e43ce 100755
--- a/dart/viscous/Physics/DClosures.py
+++ b/dart/viscous/Physics/DClosures.py
@@ -15,7 +15,6 @@
# Imports
# ------------------------------------------------------------------------------
import numpy as np
-from math import exp
class Closures:
def __init__(self) -> None:
@@ -61,13 +60,16 @@ class Closures:
def turbulentClosure(theta, H, Ct, Ret, Me, name):
"""Turbulent closure derived from Nishida and Drela (1996)
"""
+
+ H = max(H, 1.0005)
# eliminate absurd transients
Ct = max(min(Ct, 0.30), 0.0000001)
Hk = (H - 0.29*Me**2)/(1+0.113*Me**2)
Hk = max(min(Hk, 10), 1.00005)
Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me**2
- gam = 1.4 - 1
- Fc = np.sqrt(1+0.5*gam*Me**2)
+ #gam = 1.4 - 1
+ #Fc = np.sqrt(1+0.5*gam*Me**2)
+ Fc = np.sqrt(1+0.2*Me**2)
if Ret <= 400:
H0 = 4
@@ -82,9 +84,9 @@ class Closures:
Hstar = (Hk-H0)**2*(0.007*np.log(Ret)/(Hk-H0+4/np.log(Ret))**2 + 0.015/Hk) + 1.5 + 4/Ret
Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction
- logRt = np.log(Ret/Fc)
- logRt = max(logRt,3)
- arg = max(-1.33*Hk, -20)
+ #logRt = np.log(Ret/Fc)
+ #logRt = max(logRt,3)
+ #arg = max(-1.33*Hk, -20)
Us = (Hstar/2)*(1-4*(Hk-1)/(3*H)) # Equivalent normalized wall slip velocity
delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta)
Ctcon = 0.5/(6.7**2*0.75)
@@ -97,11 +99,13 @@ class Closures:
Cdw = 0 # wall contribution
Cdd = (0.995-Us)*Ct**2 # turbulent outer layer contribution
Cdl = 0.15*(0.995-Us)**2/Ret # laminar stress contribution to outer layer
+ Cteq = np.sqrt(4*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)) #Here it is Cteq^0.5
else:
if Us > 0.95:
Us = 0.98
n = 1
- Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt/2.3026)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-1))
+ #Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt/2.3026)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-1))
+ Cf = 1/Fc*(0.3*np.exp(-1.33*Hk)*(np.log10(Ret/Fc))**(-1.74-0.31*Hk) + 0.00011*(np.tanh(4-Hk/0.875) - 1))
Hkc = max(Hk-1-18/Ret, 0.01)
# Dissipation coefficient
Hmin = 1 + 2.1/np.log(Ret)
@@ -110,82 +114,12 @@ class Closures:
Cdw = 0.5*(Cf*Us) * Dfac
Cdd = (0.995-Us)*Ct**2
Cdl = 0.15*(0.995-Us)**2/Ret
- CteqSquared = n*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)
- if CteqSquared >= 0:
- Cteq = np.sqrt(CteqSquared) #Here it is Cteq^0.5
- else:
- Cteq = 0.03
+ Cteq = np.sqrt(Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)) #Here it is Cteq^0.5
Cd = n*(Cdw+ Cdd + Cdl)
- Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2)
+ #Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2)
return Cd, Cf, delta, Hstar, Hstar2, Hk, Cteq, Us
-
- def amplRoutine(Ha,Ret,Ret_previous,dt,theta,name, Ncrit=9):
- '''Compute the amplitude of the TS waves envelop for the transition
- '''
- Ret=max(Ret,1)
- Tu=None
- nu=1e-5
- if Ncrit is not None and Tu is None:
- Tu_prime=100*exp((Ncrit+8.43)/(-2.4))
- Tu=2.7*np.arctanh(Tu_prime/2.7)
- # No empirical relations for Ret/dt
- dRet_dt=(Ret-Ret_previous)/dt
-
- #lambdat=theta**2/nu*due_dx # See also Thwaites H − λθ relation (Ye Thesis)
- lambdat=0.058*(Ha-4)**2/(Ha-1)-0.068
- lambdat=min(lambdat,0.1)
- lambdat=max(-0.1,lambdat)
- if lambdat<=0:
- F_lambdat=1-(-12.986*lambdat-123.66*lambdat**2-405.689*lambdat**3)*exp(-(Tu/1.5)**1.5)
- else: # if lambdat>0
- F_lambdat=1+0.275*(1-exp(-35*lambdat))*exp(-Tu/0.5)
- A_coeff=0.1
- B_coeff=0.3
- if Tu<=1.3:
- Ret_onset=(1173.51-589.428*Tu+0.2196/(Tu**2))*F_lambdat
- else: # if Tu>1.3
- Ret_onset=(331.5*(Tu - 0.5658)**(-0.671))*F_lambdat
- # Numerical robustness
- Ret_onset=max(Ret_onset,20)
-
- # Comment one of the expression for Recrit
- # Recrit Drela
- #Recrit=10**(2.492*(1/(Ha-1))**0.43+0.7*(np.tanh(14*1/(Ha-1)-9.24)+1))
- # Recrit Arnal's (More digits more fun)
- Hmi = (1/(Ha-1))
- #logRet_crit=(0.267659/(Ha-1)+0.394429)*np.tanh(12.7886/(Ha-1)-8.57463)+3.04212/(Ha-1)+0.6660931
- logRet_crit2 = 2.492*(1/(Ha-1))**0.43 + 0.7*(np.tanh(14*Hmi-9.24)+1.0) # value at which the amplification starts to grow
-
- r=1/B_coeff*(Ret/Ret_onset-1)+1/2
- if r<=0:
- g=0
- elif 0<r<1:
- g=A_coeff*(3*r**2-2*r**3)
- else: #if r>=1
- g=A_coeff
-
- rnorm=(np.log10(Ret)-(logRet_crit2-0.08))/(2*0.08)
- if rnorm<=0:
- rfac=0
- elif 0<rnorm<1:
- rfac=3*rnorm**2-2*rnorm**3
- else: # if rnorm>=1
- rfac=1
-
- dRet_dx=-0.05+2.7*(1/(Ha-1))-5.5*(1/(Ha-1))**2+3*(1/(Ha-1))**3
-
- dN_dRet=0.028*(Ha-1)-0.0345*exp(-(3.87*1/(Ha-1)-2.52)**2)
-
- dTu_dx=0 #???
-
- dN_dTu=43/((-8.43-2.4*np.log(Tu/100))**2*Tu)
- Ax=(dN_dRet*(dRet_dx+dRet_dt)+dN_dTu*dTu_dx)*rfac+g/theta
- if name=='wake':
- Ax=0
- return Ax
-
- def amplRoutine2(Hk, Ret, theta):
+ def amplRoutine(Hk, Ret, theta):
'''Compute the amplitude of the TS waves envelop for the transition
'''
Dgr = 0.08
diff --git a/dart/viscous/Physics/DDiscretization.py b/dart/viscous/Physics/DDiscretization.py
index 85272a1924dada038365ce681f7c11bca1c734f5..5d7d5b954b09ac91795ff108558d2e45f1c9cf97 100755
--- a/dart/viscous/Physics/DDiscretization.py
+++ b/dart/viscous/Physics/DDiscretization.py
@@ -53,13 +53,14 @@ class Discretization:
dxExt[i]=xxExt[i]-xxExt[0] if i==bNodes[1] else xxExt[i]-xxExt[i-1]
self.dxExt=dxExt
- def fou(self, var, x, idxPrev2, idxPrev,idx, idxNext):
+ def fou(self, var, x, idxPrev, idx):
"""First order backward difference.
"""
du=x-var.u[var.nVar*idxPrev:var.nVar*(idxPrev+1)]
dx = var.xx[idx] - var.xx[idxPrev]
dxExt = var.extPar[idx*var.nExtPar + 4] - var.extPar[idxPrev*var.nExtPar + 4]
+
du_dx=du/dx
dHstar_dx=(var.blPar[idx*var.nBlPar+4]-var.blPar[idxPrev*var.nBlPar+4])/dx
dueExt_dx=(var.extPar[idx*var.nExtPar+2]-var.extPar[idxPrev*var.nExtPar+2])/dxExt
diff --git a/dart/viscous/Physics/DValidation.py b/dart/viscous/Physics/DValidation.py
index 5d2777316a70e8e570c2c2d0ba155ca4e4171dc5..b367785cba7522b6f79f6c2a0fba5ff85453facc 100755
--- a/dart/viscous/Physics/DValidation.py
+++ b/dart/viscous/Physics/DValidation.py
@@ -53,8 +53,8 @@ class Validation:
# Aerodynamic coefficients
print('+---------------------------- Compare with XFOIL -----------------------------+')
- print('| CL; Dart: {:<.5f}, Xfoil: {:<.5f}. {:>43s}'.format(dartcl, float(cl_xfoil), '|'))
- print('| CD; Dart: {:<.5f}, Xfoil: {:<.5f}. {:>43s}'.format(wData['Cd'], float(cd_xfoil), '|'))
+ print('- CL; Dart: {:<.5f}, Xfoil: {:<.5f}.'.format(dartcl, float(cl_xfoil)))
+ print('- CD; Dart: {:<.5f}, Xfoil: {:<.5f}.'.format(wData['Cd'], float(cd_xfoil)))
# Plot variables
plot_dStar=plt.figure(1)
diff --git a/dart/viscous/Solvers/DIntegration.py b/dart/viscous/Solvers/DIntegration.py
index 2c7f63a5313fb1703e04f4159314575e41566bae..423b169411de9d3a6fdbeb0fcdd95a82e8737d03 100644
--- a/dart/viscous/Solvers/DIntegration.py
+++ b/dart/viscous/Solvers/DIntegration.py
@@ -14,222 +14,227 @@
# ------------------------------------------------------------------------------
# Imports
# ------------------------------------------------------------------------------
+from matplotlib.colors import Normalize
from dart.viscous.Physics.DClosures import Closures
import math
import numpy as np
+from matplotlib import pyplot as plt
from fwk.coloring import ccolors
class Integration:
- """Pseudo-time integration of the integral boundary layer equations.
- The non-linear equations are solved for one point using Newton's method
-
- Attributes
- ----------
- - CFL0 (float): Starting CFL.
- - maxIt (int): Maximum number of iterations.
- - tol (float): Convergence criterion.
- - itFrozenJac0 (int): Number of iterations during which the Jacobian matrix is frozen.
- - Minf (float): Free-stream Mach number.
- - gamma (float): Fluid heat capacity ratio.
- """
-
- def __init__(self, _sysMatrix, _Minf, _Cfl0) -> None:
- self.sysMatrix = _sysMatrix
-
- self.__CFL0=_Cfl0
- self.__maxIt = 1500
- self.__tol = 1e-6
- self.itFrozenJac0 = 5
-
- self.__Minf = max(_Minf, 0.1)
- self.gamma = 1.4
- pass
-
- def GetCFL0(self): return self.__CFL0
- def SetCFL0(self, _cfl0): self.__CFL0 = _cfl0
-
- def Gettol(self): return self.__tol
- def Settol(self, _tol): self.__tol = _tol
-
- def GetmaxIt(self): return self.__maxIt
- def SetmaxIt(self, _maxIt): self.__maxIt = _maxIt
+ """Pseudo-time integration of the integral boundary layer equations.
+ The non-linear equations are solved for one point using Newton's method
+
+ Attributes
+ ----------
+ - CFL0 (float): Starting CFL.
+ - maxIt (int): Maximum number of iterations.
+ - tol (float): Convergence criterion.
+ - itFrozenJac0 (int): Number of iterations during which the Jacobian matrix is frozen.
+ - Minf (float): Free-stream Mach number.
+ - gamma (float): Fluid heat capacity ratio.
+ """
+
+ def __init__(self, _sysMatrix, _Minf, _Cfl0) -> None:
+ self.sysMatrix = _sysMatrix
+
+ self.__CFL0=_Cfl0
+ self.__maxIt = 1000
+ self.__tol = 1e-8
+ self.itFrozenJac0 = 5
+
+ self.__Minf = max(_Minf, 0.1)
+ self.gamma = 1.4
+ pass
+
+ def GetCFL0(self): return self.__CFL0
+ def SetCFL0(self, _cfl0): self.__CFL0 = _cfl0
+
+ def Gettol(self): return self.__tol
+ def Settol(self, _tol): self.__tol = _tol
+
+ def GetmaxIt(self): return self.__maxIt
+ def SetmaxIt(self, _maxIt): self.__maxIt = _maxIt
+
+ def TimeMarching (self, cfgFlow, iMarker, screenOutput = 0) -> int:
+ """Pseudo-time marching solver.
+
+ Args
+ ----
+ - DFlow (class type): Data structure containing the flow configuration.
+ - iMarker (int): Id of the point to be treated.
+
+ Returns
+ -------
+ - Exit code info (int):
+ - 0 -> sucessfull exit.
+ - >0 -> convergence to tolerance not achieved, number of iterations.
+ - <0 -> illegal input or breakdown.
+
+ Opts
+ ----
+ - screenOutput (bool, optional): Output pseudo-time marching convergence.
+ """
- def TimeMarching (self, DFlow, iMarker, screenOutput = 0) -> int:
- """Pseudo-time marching solver.
-
- Args
- ----
- - DFlow (class type): Data structure containing the flow configuration.
- - iMarker (int): Id of the point to be treated.
-
- Returns
- -------
- - Exit code info (int):
- - 0 -> sucessfull exit.
- - >0 -> convergence to tolerance not achieved, number of iterations.
- - <0 -> illegal input or breakdown.
-
- Opts
- ----
- - screenOutput (bool, optional): Output pseudo-time marching convergence.
- """
+ # Save initial condition.
- # Save initial condition.
+ sln0 = cfgFlow.u.copy()
- sln0 = DFlow.u.copy()
+ # Retreive parameters.
- # Retreive parameters.
+ nVar = cfgFlow.nVar # Number of variables of the differential problem.
+ nExtPar = cfgFlow.nExtPar # Number of boundary layer parameters used for vector indexing.
+ bNodes = cfgFlow.bNodes # Boundary nodes of the computational domain.
+ iInvVel = cfgFlow.extPar[iMarker*nExtPar + 2] # Inviscid velocity on the considered cell.
+ iMarkerPrev = 0 if iMarker == bNodes[1] else iMarker - 1 # Point upstream of the considered point.
+ dx = cfgFlow.xx[iMarker] - cfgFlow.xx[iMarkerPrev] # Cell size.
- nVar = DFlow.nVar # Number of variables of the differential problem.
- nExtPar = DFlow.nExtPar # Number of boundary layer parameters used for vector indexing.
- bNodes = DFlow.bNodes # Boundary nodes of the computational domain.
- iInvVel = DFlow.extPar[iMarker*nExtPar + 2] # Inviscid velocity on the considered cell.
- iMarkerPrev = 0 if iMarker == bNodes[1] else iMarker - 1 # Point upstream of the considered point.
- dx = DFlow.xx[iMarker] - DFlow.xx[iMarkerPrev] # Cell size.
+ # Initial residual.
- # Initial residual.
+ sysRes = self.sysMatrix.buildResiduals(cfgFlow, iMarker)
+ normSysRes0 = np.linalg.norm(sysRes)
- sysRes = self.sysMatrix.buildResiduals(DFlow, iMarker)
- normSysRes0 = np.linalg.norm(sysRes)
+ # Initialize pseudo-time integration parameters.
- # Initialize pseudo-time integration parameters.
+ CFL = self.__CFL0 # Initialized CFL number.
+ itFrozenJac = self.itFrozenJac0 # Number of iterations for which Jacobian is frozen.
+ timeStep = self.__SetTimeStep(CFL, dx, iInvVel) # Initial time step.
+ IMatrix = np.identity(nVar) / timeStep # Identity matrix used to damp the system.
+ convPlot = []
+ # Main loop.
- CFL = self.__CFL0 # Initialized CFL number.
- itFrozenJac = self.itFrozenJac0 # Number of iterations for which Jacobian is frozen.
- timeStep = self.__SetTimeStep(CFL, dx, iInvVel) # Initial time step.
- IMatrix = np.identity(nVar) / timeStep # Identity matrix used to damp the system.
+ nErrors = 0 # Count for the number of errors identified.
+ innerIters = 0 # Number of inner iterations to solver the non-linear system.
+ while innerIters <= self.__maxIt:
+
+ try:
- # Main loop.
+ # Jacobian computation.
- nErrors = 0 # Counter for the number of errors identified.
- innerIters = 0 # Number of inner iterations to solver the non-linear system.
- while innerIters <= self.__maxIt:
-
- try:
+ if innerIters % itFrozenJac == 0:
+ Jacobian=self.sysMatrix.buildJacobian(cfgFlow, iMarker,
+ sysRes, order = 1, eps = 1e-8)
- # Jacobian computation.
+ # Linear solver.
- if innerIters % itFrozenJac == 0:
- Jacobian=self.sysMatrix.buildJacobian(DFlow, iMarker, sysRes, order = 1, eps = 1e-8)
+ slnIncr = np.linalg.solve((Jacobian + IMatrix), - sysRes)
- # Linear solver.
+ # Increment solution.
- slnIncr = np.linalg.solve((Jacobian + IMatrix), - sysRes)
+ cfgFlow.u[iMarker*nVar : (iMarker+1)*nVar] += slnIncr
- # Increment solution.
+ # Residual computation.
- DFlow.u[iMarker*nVar : (iMarker+1)*nVar] += slnIncr
+ sysRes = self.sysMatrix.buildResiduals(cfgFlow, iMarker)
+ normSysRes = np.linalg.norm(sysRes + slnIncr/timeStep)
+ convPlot.append(normSysRes/normSysRes0)
- # Residual computation.
+ # Check convergence.
+
+ if normSysRes <= self.__tol * normSysRes0:
+ # Convergence criterion satisfied.
+ return 0
- sysRes = self.sysMatrix.buildResiduals(DFlow, iMarker)
- normSysRes = np.linalg.norm(slnIncr/timeStep - sysRes)
+ except KeyboardInterrupt:
+ print(ccolors.ANSI_RED + 'Program terminated by user.', ccolors.ANSI_RESET)
+ quit()
+
+ except Exception:
- # Check convergence.
-
- if normSysRes <= self.__tol * normSysRes0:
- # Convergence criterion satisfied.
- return 0
+ nErrors += 1
- except KeyboardInterrupt:
- print(ccolors.ANSI_RED + 'Program terminated by user.', ccolors.ANSI_RESET)
- quit()
- except Exception:
+ if nErrors >= 5:
+ self.__ResetSolution(cfgFlow, iMarker, sln0)
+ return -1
- nErrors += 1
+ cfgFlow.u[iMarker*nVar : (iMarker+1)*nVar] = cfgFlow.u[(iMarker-1)*nVar : (iMarker)*nVar] # Reset solution.
- if nErrors >= 5:
- self.__ResetSolution(DFlow, iMarker, sln0)
- return -1
+ # Lower CFL and recompute time step
- DFlow.u[iMarker*nVar : (iMarker+1)*nVar] = DFlow.u[(iMarker-1)*nVar : (iMarker)*nVar] # Reset solution.
+ if math.isnan(CFL):
+ CFL = self.__CFL0
+ CFL = min(max(CFL / 2,0.3),1)
+ timeStep = self.__SetTimeStep(CFL, dx, iInvVel)
+ IMatrix = np.identity(nVar) / timeStep
- # Lower CFL and recompute time step
+ sysRes = self.sysMatrix.buildResiduals(cfgFlow, iMarker)
- if math.isnan(CFL):
- CFL = self.__CFL0
- CFL = min(max(CFL / 2,0.3),1)
- timeStep = self.__SetTimeStep(CFL, dx, iInvVel)
- IMatrix = np.identity(nVar) / timeStep
+ itFrozenJac = 1
+ continue
- sysRes = self.sysMatrix.buildResiduals(DFlow, iMarker)
+ # CFL adaptation.
- itFrozenJac = 1
- continue
+ CFL = self.__CFL0* (normSysRes0/normSysRes)**0.7
- # CFL adaptation.
+ CFL = max(CFL, 1)
+ timeStep = self.__SetTimeStep(CFL, dx, iInvVel)
+ IMatrix = np.identity(nVar) / timeStep
- CFL = self.__CFL0* (normSysRes0/normSysRes)**0.7
- CFL = max(CFL, 0.5)
- timeStep = self.__SetTimeStep(CFL, dx, iInvVel)
- IMatrix = np.identity(nVar) / timeStep
+ innerIters += 1
+ else: # Maximum number of iterations reached.
- innerIters += 1
+ if normSysRes >= 1e-3 * normSysRes0:
+ self.__ResetSolution(cfgFlow, iMarker, sln0)
+ return innerIters
- else: # Maximum number of iterations reached.
- if normSysRes >= 1e-3 * normSysRes0:
- self.__ResetSolution(DFlow, iMarker, sln0)
- return innerIters
+ def __GetSoundSpeed(self, invVelSquared) -> float:
+ """Compute the speed of sound in a cell.
+ Args
+ ----
+ - invVelSquared (float): Square of the inviscid velocity.
+ """
+
+ return np.sqrt(1 / (self.__Minf * self.__Minf) + 0.5 * (self.gamma - 1) * (1 - invVelSquared))
- def __GetSoundSpeed(self, invVelSquared) -> float:
- """Compute the speed of sound in a cell.
+ def __SetTimeStep(self, CFL, dx, invVel) -> float:
+ """Compute the pseudo-time step in a cell for a given
+ CFL number.
- Args
- ----
- - invVelSquared (float): Square of the inviscid velocity.
- """
-
- return np.sqrt(1 / (self.__Minf * self.__Minf) + 0.5 * (self.gamma - 1) * (1 - invVelSquared))
+ Args
+ ----
+ - CFL (float): Current CFL number.
+ - dx (float): Cell size.
+ - invVel (float): Inviscid velocity in the cell.
- def __SetTimeStep(self, CFL, dx, invVel) -> float:
- """Compute the pseudo-time step in a cell for a given
- CFL number.
+ Returns
+ -------
+ - timeStep (float): Pseudo-time step.
+ """
- Args
- ----
- - CFL (float): Current CFL number.
- - dx (float): Cell size.
- - invVel (float): Inviscid velocity in the cell.
+ # Speed of sound.
+ gradU2 = (invVel)**2
+ # Velocity of the fastest travelling wave
+ soundSpeed = self.__GetSoundSpeed(gradU2)
+ advVel = soundSpeed + invVel
- Returns
- -------
- - timeStep (float): Pseudo-time step.
- """
+ # Time step computation
+ return CFL*dx/advVel
- # Speed of sound.
- gradU2 = (invVel)**2
- # Velocity of the fastest travelling wave
- soundSpeed = self.__GetSoundSpeed(gradU2)
- advVel = soundSpeed + invVel
+ def __ResetSolution(self, DFlow, iMarker, sln0) -> None:
- # Time step computation
- return CFL*dx/advVel
+ nVar = DFlow.nVar
+ nBlPar = DFlow.nBlPar
+ nExtPar = DFlow.nExtPar
+
+ DFlow.u[iMarker*nVar:(iMarker+1)*nVar] = sln0[(iMarker)*nVar: (iMarker+1)*nVar].copy()
+ name = 'airfoil' if iMarker <= DFlow.bNodes[2] else 'wake'
- def __ResetSolution(self, DFlow, iMarker, sln0) -> None:
+ if DFlow.flowRegime[iMarker]==0:
- nVar = DFlow.nVar
- nBlPar = DFlow.nBlPar
- nExtPar = DFlow.nExtPar
-
- DFlow.u[iMarker*nVar:(iMarker+1)*nVar] = sln0[(iMarker)*nVar: (iMarker+1)*nVar].copy()
- name = 'airfoil' if iMarker <= DFlow.bNodes[2] else 'wake'
+ # Laminar closure relations
+ DFlow.blPar[iMarker*nBlPar+1:iMarker*nBlPar+7] = Closures.laminarClosure(DFlow.u[iMarker*nVar + 0], DFlow.u[iMarker*nVar + 1],
+ DFlow.blPar[iMarker*nBlPar+0], DFlow.extPar[iMarker*nExtPar+0],
+ name)
- if DFlow.flowRegime[iMarker]==0:
+ elif DFlow.flowRegime[iMarker]==1:
- # Laminar closure relations
- DFlow.blPar[iMarker*nBlPar+1:iMarker*nBlPar+7] = Closures.laminarClosure(DFlow.u[iMarker*nVar + 0], DFlow.u[iMarker*nVar + 1],
- DFlow.blPar[iMarker*nBlPar+0], DFlow.extPar[iMarker*nExtPar+0],
- name)
-
- elif DFlow.flowRegime[iMarker]==1:
-
- # Turbulent closures relations
- DFlow.blPar[iMarker*nBlPar+1:iMarker*nBlPar+9] = Closures.turbulentClosure(DFlow.u[iMarker*nVar + 0], DFlow.u[iMarker*nVar + 1],
- DFlow.u[iMarker*nVar + 4], DFlow.blPar[iMarker*nBlPar+0],
- DFlow.extPar[iMarker*nExtPar+0], name)
+ # Turbulent closures relations
+ DFlow.blPar[iMarker*nBlPar+1:iMarker*nBlPar+9] = Closures.turbulentClosure(DFlow.u[iMarker*nVar + 0], DFlow.u[iMarker*nVar + 1],
+ DFlow.u[iMarker*nVar + 4], DFlow.blPar[iMarker*nBlPar+0],
+ DFlow.extPar[iMarker*nExtPar+0], name)
diff --git a/dart/viscous/Solvers/DSolver.py b/dart/viscous/Solvers/DSolver.py
index bd7493c25311d2983398a533ea6c3fca05223cbd..fcefa76f14ca96d735a6e8aea1a872f3d87610dd 100644
--- a/dart/viscous/Solvers/DSolver.py
+++ b/dart/viscous/Solvers/DSolver.py
@@ -20,159 +20,162 @@ from fwk.coloring import ccolors
import numpy as np
class Solver:
- """Performs the downstream marching process. Calls the pseudo-time integration based
- on Newton's method on demand.
-
- Attributes
- ----------
-
- - integrator (class type): Module able to perform pseudo-time integration of one point.
- - transSolver (class type): Module able to treat the transition related computation.
- - wakeBC (class type): Wake boundary condition.
- - initSln (bool): Flag to initialize the initial condition or to use the value in input.
- """
-
- def __init__(self, _integrator, transition, wakeBC) -> None:
-
- # Initialize time parameters
-
- self.integrator = _integrator
- self.transSolver = transition
- self.wakeBC = wakeBC
- self.initSln = 0
-
- def Run(self, DFlow):
- """Runs the downstream marching solver.
- """
-
- # Retreive parameters.
-
- nMarkers = DFlow.nN # Number of points in the computational domain.
- bNodes = DFlow.bNodes # Boundary nodes of the computational domain.
- convergedPoints = -1*np.ones(DFlow.nN) # Convergence control of each point.
- convergedPoints[0] = 0 # Boundary condition.
-
- self.transSolver.initTransition(DFlow) # Initialize the flow regime on each cell.
-
- lockTrans = 0 # Flag to remember if the transition is already found or not.
-
- for iMarker in range(1, nMarkers):
- displayState = 0 # Flag to output simulation state.
-
- if self.initSln:
- iMarkerPrev = 0 if iMarker == DFlow.bNodes[1] else iMarker - 1
- DFlow.u[iMarker*DFlow.nVar: (iMarker + 1)*DFlow.nVar] = DFlow.u[iMarkerPrev*DFlow.nVar: (iMarkerPrev + 1)* DFlow.nVar]
- if DFlow.flowRegime[iMarker] == 1 and DFlow.u[iMarker * DFlow.nVar + 4] == 0:
- DFlow.u[iMarker * DFlow.nVar + 4] = 0.03
-
- if iMarker == bNodes[0] + 1: # Upper side start.
- print('- Upper side,', end = ' ')
-
- if iMarker == bNodes[1]: # Lower side start.
- if lockTrans == 0:
- print('no transition.')
- lockTrans = 0
- print('- Lower side,', end = ' ')
-
- if iMarker == bNodes[2]: # First wake point
-
- self.wakeBC(DFlow) # Wake boundary condition.
- convergedPoints[iMarker] = 0
- if lockTrans == 0:
- print('no transition.')
- lockTrans = 1
-
- continue # Ignore remaining instructions for the first wake point.
-
- # Call Newton solver for the point.
-
- convergedPoints[iMarker] = self.integrator.TimeMarching(DFlow, iMarker, displayState)
-
- if convergedPoints[iMarker] != 0:
- print(ccolors.ANSI_RED + 'Marker {:<4.0f} (x/c={:<1.4f}): PTI terminated with exit code {:<4.0f}.'
- .format(iMarker, DFlow.xGlobal[iMarker], convergedPoints[iMarker]), ccolors.ANSI_RESET)
-
- # Check for transition.
-
- if not lockTrans:
-
- self.__CheckWaves(DFlow, iMarker)
-
- # Free transition.
- if DFlow.u[iMarker * DFlow.nVar + 2] >= DFlow.Ncrit:
- self.__AverageTransition(DFlow, iMarker, displayState)
- if iMarker < DFlow.bNodes[1]:
- print('free transition x/c = {:<1.5f} {:>4.0f}.'.format(DFlow.xtr[0],iMarker))
- elif iMarker < DFlow.bNodes[2]:
- print('free transition x/c = {:<1.5f} {:>4.0f}.'.format(DFlow.xtr[1],iMarker))
- else:
- print('error in wake')
- lockTrans = 1
-
- # Forced transition.
- elif (DFlow.xtrF[0] is not None and DFlow.xtrF[0] !=0 and DFlow.xtrF[0] != 1) or (DFlow.xtrF[1] is not None and DFlow.xtrF[1] !=0 and DFlow.xtrF[1] != 1):
- if DFlow.flowRegime[iMarker] == 1 and DFlow.flowRegime[iMarker - 1] == 0:
- print('forced transition near x/c = {:<2.3f} {:>4.0f}'.format(DFlow.xGlobal[iMarker],iMarker))
- self.__AverageTransition(DFlow, iMarker, displayState)
- lockTrans = 1
-
- return convergedPoints
-
- def __CheckWaves(self, var, iMarker) -> None:
- """Determine if the amplification waves start growing or not.
- """
-
- logRet_crit = 2.492*(1/(var.blPar[(iMarker)*var.nBlPar + 6]-1))**0.43
- + 0.7*(np.tanh(14*(1/(var.blPar[(iMarker)*var.nBlPar + 6]-1))-9.24)+1.0)
-
- if var.blPar[(iMarker)*var.nBlPar + 0] > 0: # Avoid log of negative number.
- if np.log10(var.blPar[(iMarker)*var.nBlPar + 0]) <= logRet_crit - 0.08:
- var.u[iMarker*var.nVar + 2] = 0
-
- pass
-
-
- def __AverageTransition(self, var, iMarker, displayState) -> None:
- """Averages laminar and turbulent solution at the transition location.
- The boundary condition of the turbulent flow portion is also imposed.
- """
-
- # Save laminar solution.
- lamSol = var.u[(iMarker)*var.nVar : (iMarker+1)* var.nVar].copy()
-
- # Set up turbulent portion boundary condition.
- avTurb = self.transSolver.solveTransition(var, iMarker)
-
- # Compute turbulent solution.
- var.flowRegime[iMarker] = 1
- iMarkerPrev = 0 if iMarker == var.bNodes[1] else iMarker - 1
-
- # Avoid starting with 0 for the shear stress and high H.
- var.u[iMarker*var.nVar + 4] = var.u[iMarkerPrev*var.nVar + 4]
- var.u[iMarker*var.nVar + 1] = 1.515
-
- self.integrator.TimeMarching(var, iMarker, displayState)
-
- # Average solutions.
- avLam = 1 - avTurb
- thetaTrans = avLam * lamSol[0] + avTurb * var.u[(iMarker)*var.nVar + 0]
- HTrans = avLam * lamSol[1] + avTurb * var.u[(iMarker)*var.nVar + 1]
- nTrans = 0
- ueTrans = avLam * lamSol[3] + avTurb * var.u[(iMarker)*var.nVar + 3]
- CtTrans = avTurb * var.u[(iMarker)*var.nVar + 4]
-
- var.u[(iMarker)*var.nVar : (iMarker+1)*var.nVar] = [thetaTrans, HTrans, nTrans, ueTrans, CtTrans]
-
- # Recompute closures.
- var.blPar[(iMarker)*var.nBlPar+1:(iMarker)*var.nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker)*var.nVar+0],
- var.u[(iMarker)*var.nVar+1],
- var.u[(iMarker)*var.nVar+4],
- var.blPar[(iMarker)*var.nBlPar+0],
- var.extPar[(iMarker)*var.nExtPar+0],
- 'airfoil')
- var.blPar[(iMarker-1)*var.nBlPar+1:(iMarker-1)*var.nBlPar+7]=Closures.laminarClosure(var.u[(iMarker-1)*var.nVar+0],
- var.u[(iMarker-1)*var.nVar+1],
- var.blPar[(iMarker-1)*var.nBlPar+0],
- var.extPar[(iMarker-1)*var.nExtPar+0],
- 'airfoil')
- pass
\ No newline at end of file
+ """Performs the downstream marching process. Calls the pseudo-time integration based
+ on Newton's method on demand.
+
+ Attributes
+ ----------
+
+ - integrator (class type): Module able to perform pseudo-time integration of one point.
+ - transSolver (class type): Module able to treat the transition related computation.
+ - wakeBC (class type): Wake boundary condition.
+ - initSln (bool): Flag to initialize the initial condition or to use the value in input.
+ """
+
+ def __init__(self, _integrator, transition, wakeBC) -> None:
+
+ # Initialize time parameters
+
+ self.integrator = _integrator
+ self.transSolver = transition
+ self.wakeBC = wakeBC
+ self.initSln = 0
+
+ def Run(self, DFlow):
+ """Runs the downstream marching solver.
+ """
+
+ # Retreive parameters.
+
+ nMarkers = DFlow.nN # Number of points in the computational domain.
+ bNodes = DFlow.bNodes # Boundary nodes of the computational domain.
+ convergedPoints = -1*np.ones(DFlow.nN) # Convergence control of each point.
+ convergedPoints[0] = 0 # Boundary condition.
+
+ self.transSolver.initTransition(DFlow) # Initialize the flow regime on each cell.
+
+ lockTrans = 0 # Flag to remember if the transition is already found or not.
+
+ for iMarker in range(1, nMarkers):
+ displayState = 0 # Flag to output simulation state.
+
+ if self.initSln:
+ iMarkerPrev = 0 if iMarker == DFlow.bNodes[1] else iMarker - 1
+ DFlow.u[iMarker*DFlow.nVar: (iMarker + 1)*DFlow.nVar] = DFlow.u[iMarkerPrev*DFlow.nVar: (iMarkerPrev + 1)* DFlow.nVar]
+ if DFlow.flowRegime[iMarker] == 1 and DFlow.u[iMarker * DFlow.nVar + 4] == 0:
+ DFlow.u[iMarker * DFlow.nVar + 4] = 0.03
+
+ if iMarker == bNodes[0] + 1: # Upper side start.
+ print('- Upper side,', end = ' ')
+
+ if iMarker == bNodes[1]: # Lower side start.
+ if lockTrans == 0:
+ print('no transition.')
+ lockTrans = 0
+ print('- Lower side,', end = ' ')
+
+ if iMarker == bNodes[2]: # First wake point
+
+ self.wakeBC(DFlow) # Wake boundary condition.
+ convergedPoints[iMarker] = 0
+ if lockTrans == 0:
+ print('no transition.')
+ lockTrans = 1
+
+ continue # Ignore remaining instructions for the first wake point.
+
+ # Call Newton solver for the point.
+
+ convergedPoints[iMarker] = self.integrator.TimeMarching(DFlow, iMarker, displayState)
+
+ if convergedPoints[iMarker] != 0:
+ print(ccolors.ANSI_RED + 'Marker {:<4.0f} (x/c={:<1.4f}): PTI terminated with exit code {:<4.0f}.'
+ .format(iMarker, DFlow.xGlobal[iMarker], convergedPoints[iMarker]), ccolors.ANSI_RESET)
+
+ # Check for transition.
+
+ if not lockTrans:
+
+ self.__CheckWaves(DFlow, iMarker)
+
+ # Free transition.
+ if DFlow.u[iMarker * DFlow.nVar + 2] >= DFlow.Ncrit:
+ self.__AverageTransition(DFlow, iMarker, displayState)
+ if iMarker < DFlow.bNodes[1]:
+ print('free transition x/c = {:<1.5f} {:>4.0f}.'.format(DFlow.xtr[0],iMarker))
+ elif iMarker < DFlow.bNodes[2]:
+ print('free transition x/c = {:<1.5f} {:>4.0f}.'.format(DFlow.xtr[1],iMarker))
+ else:
+ print('error in wake')
+ lockTrans = 1
+
+ # Forced transition.
+ elif DFlow.xtrF[0] is not None or DFlow.xtrF[1] is not None:
+ if iMarker != 1 and DFlow.flowRegime[iMarker] == 1 and DFlow.flowRegime[iMarker - 1] == 0:
+ print('forced transition near x/c = {:<2.3f} {:>4.0f}'.format(DFlow.xGlobal[iMarker],iMarker))
+ self.__AverageTransition(DFlow, iMarker, displayState, forced = 1)
+ lockTrans = 1
+
+ return convergedPoints
+
+ def __CheckWaves(self, var, iMarker) -> None:
+ """Determine if the amplification waves start growing or not.
+ """
+
+ logRet_crit = 2.492*(1/(var.blPar[(iMarker)*var.nBlPar + 6]-1))**0.43
+ + 0.7*(np.tanh(14*(1/(var.blPar[(iMarker)*var.nBlPar + 6]-1))-9.24)+1.0)
+
+ if var.blPar[(iMarker)*var.nBlPar + 0] > 0: # Avoid log of negative number.
+ if np.log10(var.blPar[(iMarker)*var.nBlPar + 0]) <= logRet_crit - 0.08:
+ var.u[iMarker*var.nVar + 2] = 0
+
+ pass
+
+
+ def __AverageTransition(self, var, iMarker, displayState, forced = 0) -> None:
+ """Averages laminar and turbulent solution at the transition location.
+ The boundary condition of the turbulent flow portion is also imposed.
+ """
+
+ # Save laminar solution.
+ lamSol = var.u[(iMarker)*var.nVar : (iMarker+1)* var.nVar].copy()
+
+ # Set up turbulent portion boundary condition.
+ if not forced:
+ avTurb = self.transSolver.solveTransition(var, iMarker)
+ else:
+ avTurb = self.transSolver.forcedTransition(var, iMarker)
+
+ # Compute turbulent solution.
+ var.flowRegime[iMarker] = 1
+ iMarkerPrev = 0 if iMarker == var.bNodes[1] else iMarker - 1
+
+ # Avoid starting with 0 for the shear stress and high H.
+ var.u[iMarker*var.nVar + 4] = var.u[iMarkerPrev*var.nVar + 4]
+ var.u[iMarker*var.nVar + 1] = 1.515
+
+ self.integrator.TimeMarching(var, iMarker, displayState)
+
+ # Average solutions.
+ avLam = 1 - avTurb
+ thetaTrans = avLam * lamSol[0] + avTurb * var.u[(iMarker)*var.nVar + 0]
+ HTrans = avLam * lamSol[1] + avTurb * var.u[(iMarker)*var.nVar + 1]
+ nTrans = 0
+ ueTrans = avLam * lamSol[3] + avTurb * var.u[(iMarker)*var.nVar + 3]
+ CtTrans = avTurb * var.u[(iMarker)*var.nVar + 4]
+
+ var.u[(iMarker)*var.nVar : (iMarker+1)*var.nVar] = [thetaTrans, HTrans, nTrans, ueTrans, CtTrans]
+
+ # Recompute closures.
+ var.blPar[(iMarker)*var.nBlPar+1:(iMarker)*var.nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker)*var.nVar+0],
+ var.u[(iMarker)*var.nVar+1],
+ var.u[(iMarker)*var.nVar+4],
+ var.blPar[(iMarker)*var.nBlPar+0],
+ var.extPar[(iMarker)*var.nExtPar+0],
+ 'airfoil')
+ var.blPar[(iMarker-1)*var.nBlPar+1:(iMarker-1)*var.nBlPar+7]=Closures.laminarClosure(var.u[(iMarker-1)*var.nVar+0],
+ var.u[(iMarker-1)*var.nVar+1],
+ var.blPar[(iMarker-1)*var.nBlPar+0],
+ var.extPar[(iMarker-1)*var.nExtPar+0],
+ 'airfoil')
+ pass
\ No newline at end of file
diff --git a/dart/viscous/Solvers/DSysMatrix.py b/dart/viscous/Solvers/DSysMatrix.py
index 47ccfd13b23cfec1351a5d8b817686dd74581821..1da7e742c2434f4bf2a4c3a851a04e37f78b26e9 100644
--- a/dart/viscous/Solvers/DSysMatrix.py
+++ b/dart/viscous/Solvers/DSysMatrix.py
@@ -20,276 +20,309 @@ import numpy as np
class flowEquations:
- def __init__(self, _setGradients, _fouScheme) -> None:
- self.setGradients = _setGradients
- self.fouScheme = _fouScheme
- pass
-
- def _blLaws(self, dataCont, iMarker, x = None) -> np.ndarray:
- """Evaluates the boundary layer laws at one point of the computational domain.
-
- Args:
- ----
- - dataCont (classType): Data container. Contains the geometry, the solution, BL parameters...
-
- - iMarker (integer): Id of the point.
-
- - x (np.ndarray, optional): Prescribed solution, usually a perturbed one for Jacobian computation purposes. Defaults to None.
-
- Returns:
- -------
- - linSysRes (np.ndarray): Residuals of the linear system at the considered point.
- """
-
- if iMarker == dataCont.bNodes[2]:
- return np.zeros(dataCont.nVar)
-
- timeMatrix = np.empty((dataCont.nVar, dataCont.nVar))
- spaceMatrix = np.empty(dataCont.nVar)
-
- if x is None: x = dataCont.u[iMarker*dataCont.nVar : (iMarker+1)*dataCont.nVar]
-
- iMarkerPrev2, iMarkerPrev, iMarkerNext, name, xtr, xtrF = self.__evalPoint(iMarker, dataCont.bNodes, dataCont.xtr, dataCont.xtrF, dataCont.nN)
-
- dissipRatio = 0.9 if name == 'wake' else 1 # Wall/wake dissipation length ratio
-
- # Reynolds number based on momentum thickness theta
- if dataCont.flowRegime[iMarker] == 0: dataCont.blPar[iMarker*dataCont.nBlPar+0] = max(x[3] * x[0] * dataCont.Re, 1)
- else: dataCont.blPar[iMarker*dataCont.nBlPar+0] = x[3] * x[0] * dataCont.Re
-
- dataCont.blPar[iMarker * dataCont.nBlPar + 9] = x[0] * x[1] # deltaStar
-
- if dataCont.flowRegime[iMarker]==0:
-
- # Laminar closure relations.
- dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+7] = Closures.laminarClosure(x[0], x[1],
- dataCont.blPar[iMarker*dataCont.nBlPar+0],
- dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
-
- elif dataCont.flowRegime[iMarker]==1:
-
- # Turbulent closures relations.
- dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+9] = Closures.turbulentClosure(x[0], x[1],
- x[4], dataCont.blPar[iMarker*dataCont.nBlPar+0],
- dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
-
- # Rename variables for clarity.
- Ta, Ha, _, uea, Cta = x
- Cda, Cfa, deltaa, Hstara, Hstar2a, Hka, Cteqa, Usa = dataCont.blPar[iMarker * dataCont.nBlPar + 1 : iMarker * dataCont.nBlPar + 9]
-
- # Amplification factor for e^n method.
- if name=='airfoil' and xtrF is None and dataCont.flowRegime[iMarker]==0:
-
- Axa = Closures.amplRoutine2(dataCont.blPar[iMarker * dataCont.nBlPar + 6], dataCont.blPar[iMarker * dataCont.nBlPar + 0], Ta)
-
- else:
-
- Axa=0
- dN_dx=0
-
- if iMarker == 1 or iMarker == dataCont.bNodes[1] or iMarker == dataCont.bNodes[2]+1:
- [dTheta_dx, dH_dx, dN_dx, due_dx, dCt_dx], dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt = self.fouScheme(dataCont, x,
- iMarkerPrev2, iMarkerPrev,
- iMarker, iMarkerNext)
- else:
- [dTheta_dx, dH_dx, dN_dx, due_dx, dCt_dx], dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt = self.setGradients(dataCont, x,
- iMarkerPrev2, iMarkerPrev,
- iMarker, iMarkerNext)
-
- # +------------------------------------------- Equations -----------------------------------------------------------+
- # | |
- # | |
- # | Von Karman momentum integral equation 0-th (momentum integral) |
- # | -------------------------------------------------------------- |
- # | |
- # | H/ue*dTheta/dt + Theta/ue*dH/dt + Theta*H/(ue^2)*due/dt + dTheta/dx + (2+H)*(Theta/ue)*due/dx - Cf/2 = 0 |
- # | |
- # | 1-st (kinetic-energy integral) moments of momentum equation |
- # | ----------------------------------------------------------- |
- # | |
- # | (1+H*(1-Hstar))/ue*dTheta/dt + (1-Hstar)*Theta/ue*dH/dt + (2-Hstar*H)*Theta/(ue^2)*due/dt |
- # | + Theta*(dHstar/dx) + (2*Hstar2 + Hstar*(1-H))*Theta/ue*due_dx - 2*Cd + Hstar*Cf/2 = 0 |
- # | |
- # | Unsteady e^n method |
- # | ------------------- |
- # | |
- # | dN/dt + dN/dx - Ax = 0 |
- # | |
- # | Interaction law |
- # | --------------- |
- # | |
- # | due/dt - c*H*dTheta/dt - c*Theta*dH/dt + due/dx-c*(H*dTheta/dx+Theta*dH/dx) + (-due/dx + c*deltaStar)_Old = 0 |
- # | |
- # | Unsteady shear lag equation (ignored in the laminar case) |
- # | --------------------------------------------------------- |
- # | |
- # | 2*delta/(Us*ue^2)*due/dt + 2*delta/(ue*Ct*Us)*dCt/dt + (2*delta/Ct)*(dCt/dx) |
- # | - 5.6*(Cteq-Ct*dissipRatio)-2*delta*(4/(3*H*Theta)*(Cf/2-((Hk-1)/(6.7*Hk*dissipRatio))^2)-1/ue*due/dx) = 0 |
- # | |
- # +-----------------------------------------------------------------------------------------------------------------+
-
-
- # Spacial part
- spaceMatrix[0] = dTheta_dx + (2+Ha) * (Ta/uea) * due_dx - Cfa/2
- spaceMatrix[1] = Ta * dHstar_dx + ( 2*Hstar2a + Hstara * (1-Ha)) * Ta/uea * due_dx - 2*Cda + Hstara * Cfa/2
- if xtrF is not None or dataCont.flowRegime[iMarker] == 1:
- spaceMatrix[2] = 0
- else:
- spaceMatrix[2] = dN_dx - Axa
- spaceMatrix[3] = due_dx -c*(Ha * dTheta_dx + Ta * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx
- if dataCont.flowRegime[iMarker] == 1:
-
- spaceMatrix[4] = (2*deltaa/Cta) * (dCt_dx) - 5.6*(Cteqa - Cta * dissipRatio)
- - 2*deltaa*(4/(3*Ha * Ta) * (Cfa/2 - ((Hka-1)/(6.7*Hka * dissipRatio))**2) - 1/uea * due_dx)
-
- elif dataCont.flowRegime[iMarker] == 0: spaceMatrix[4] = 0
-
- # Temporal part
- timeMatrix[0] = [Ha/uea, Ta/uea, 0, Ta * Ha/(uea**2), 0]
- timeMatrix[1] = [(1+Ha * (1-Hstara))/uea, (1-Hstara) * Ta/uea, 0, (2-Hstara * Ha)*Ta/(uea**2), 0]
- timeMatrix[2] = [0, 0, 1, 0, 0]
- timeMatrix[3] = [-c*Ha, -c*Ta, 0, 1, 0]
- if dataCont.flowRegime[iMarker]==1:
- timeMatrix[4] = [0, 0, 0, 2*deltaa/(Usa * uea**2), 2*deltaa/(uea * Cta * Usa)]
- else:
- # Shear lag equation ignored in the laminar portion of the flow
- timeMatrix[4] = [0, 0, 0, 0, 1]
+ def __init__(self, _setGradients, _fouScheme) -> None:
+ self.setGradients = _setGradients
+ self.fouScheme = _fouScheme
+ pass
+
+ def _blLaws(self, dataCont, iMarker, x) -> np.ndarray:
+ """Evaluates the boundary layer laws at one point of the computational domain.
+
+ Args:
+ ----
+ - dataCont (classType): Data container. Contains the geometry, the solution, BL parameters...
+
+ - iMarker (integer): Id of the point.
+
+ - x (np.ndarray, optional): Prescribed solution, usually a perturbed one for Jacobian computation purposes. Defaults to None.
+
+ Returns:
+ -------
+ - linSysRes (np.ndarray): Residuals of the linear system at the considered point.
+ """
+
+ timeMatrix = np.empty((dataCont.nVar, dataCont.nVar))
+ spaceMatrix = np.empty(dataCont.nVar)
+
+ iMarkerPrev, name, xtrF = self.__evalPoint(dataCont, iMarker)
+
+ dissipRatio = 0.9 if name == 'wake' else 1 # Wall/wake dissipation length ratio
+
+ dataCont.blPar[iMarker*dataCont.nBlPar+0] = max(x[3] * x[0] * dataCont.Re, 1) # Reynolds number based on momentum thickness theta
+
+ dataCont.blPar[iMarker * dataCont.nBlPar + 9] = x[0] * x[1] # deltaStar
+
+ if dataCont.flowRegime[iMarker]==0:
+
+ # Laminar closure relations.
+ dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+7] = Closures.laminarClosure(x[0], x[1],
+ dataCont.blPar[iMarker*dataCont.nBlPar+0],
+ dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
+
+ elif dataCont.flowRegime[iMarker]==1:
+
+ # Turbulent closures relations.
+ dataCont.blPar[iMarker*dataCont.nBlPar+1:iMarker*dataCont.nBlPar+9] = Closures.turbulentClosure(x[0], x[1],
+ x[4], dataCont.blPar[iMarker*dataCont.nBlPar+0],
+ dataCont.extPar[iMarker*dataCont.nExtPar+0], name)
+
+ # Rename variables for clarity.
+ Ta, Ha, _, uea, Cta = x
+ Mea = dataCont.extPar[iMarker*dataCont.nExtPar+0]
+ Cda, Cfa, deltaa, Hstara, Hstar2a, Hka, Cteqa, Usa = dataCont.blPar[iMarker * dataCont.nBlPar + 1 : iMarker * dataCont.nBlPar + 9]
+
+ # Amplification factor for e^n method.
+ if name=='airfoil' and xtrF is None and dataCont.flowRegime[iMarker]==0:
+
+ Axa = Closures.amplRoutine(dataCont.blPar[iMarker * dataCont.nBlPar + 6], dataCont.blPar[iMarker * dataCont.nBlPar + 0], Ta)
+
+ else:
+
+ Axa=0
+ dN_dx=0
+
+ [dTheta_dx, dH_dx, dN_dx, due_dx, dCt_dx], dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt = self.fouScheme(dataCont, x,
+ iMarkerPrev,
+ iMarker)
+
+ # +------------------------------------------- Equations -----------------------------------------------------------+
+ # | |
+ # | |
+ # | Von Karman momentum integral equation 0-th (momentum integral) |
+ # | -------------------------------------------------------------- |
+ # | |
+ # | H/ue*dTheta/dt + Theta/ue*dH/dt + Theta*H/(ue^2)*due/dt + dTheta/dx + (2+H)*(Theta/ue)*due/dx - Cf/2 = 0 |
+ # | |
+ # | 1-st (kinetic-energy integral) moments of momentum equation |
+ # | ----------------------------------------------------------- |
+ # | |
+ # | (1+H*(1-Hstar))/ue*dTheta/dt + (1-Hstar)*Theta/ue*dH/dt + (2-Hstar*H)*Theta/(ue^2)*due/dt |
+ # | + Theta*(dHstar/dx) + (2*Hstar2 + Hstar*(1-H))*Theta/ue*due_dx - 2*Cd + Hstar*Cf/2 = 0 |
+ # | |
+ # | Unsteady e^n method |
+ # | ------------------- |
+ # | |
+ # | dN/dt + dN/dx - Ax = 0 |
+ # | |
+ # | Interaction law |
+ # | --------------- |
+ # | |
+ # | due/dt - c*H*dTheta/dt - c*Theta*dH/dt + due/dx-c*(H*dTheta/dx+Theta*dH/dx) + (-due/dx + c*deltaStar)_Old = 0 |
+ # | |
+ # | Unsteady shear lag equation (ignored in the laminar case) |
+ # | --------------------------------------------------------- |
+ # | |
+ # | 2*delta/(Us*ue^2)*due/dt + 2*delta/(ue*Ct*Us)*dCt/dt + (2*delta/Ct)*(dCt/dx) |
+ # | - 5.6*(Cteq-Ct*dissipRatio)-2*delta*(4/(3*H*Theta)*(Cf/2-((Hk-1)/(6.7*Hk*dissipRatio))^2)-1/ue*due/dx) = 0 |
+ # | |
+ # +-----------------------------------------------------------------------------------------------------------------+
+
+
+ # Spacial part
+ spaceMatrix[0] = dTheta_dx + (2+Ha-Mea**2) * (Ta/uea) * due_dx - Cfa/2
+ spaceMatrix[1] = Ta * dHstar_dx + ( 2*Hstar2a + Hstara * (1-Ha)) * Ta/uea * due_dx - 2*Cda + Hstara * Cfa/2
+ if xtrF is not None or dataCont.flowRegime[iMarker] == 1:
+ spaceMatrix[2] = 0
+ else:
+ spaceMatrix[2] = dN_dx - Axa
+ spaceMatrix[3] = due_dx -c*(Ha * dTheta_dx + Ta * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx
+ if dataCont.flowRegime[iMarker] == 1:
+
+ spaceMatrix[4] = (2*deltaa/Cta) * (dCt_dx) - 5.6*(Cteqa - Cta * dissipRatio)
+ - 2*deltaa*(4/(3*Ha * Ta) * (Cfa/2 - ((Hka-1)/(6.7*Hka * dissipRatio))**2) - 1/uea * due_dx)
+
+ elif dataCont.flowRegime[iMarker] == 0: spaceMatrix[4] = 0
+
+ # Temporal part
+ timeMatrix[0] = [Ha/uea, Ta/uea, 0, Ta * Ha/(uea**2), 0]
+ timeMatrix[1] = [(1+Ha * (1-Hstara))/uea, (1-Hstara) * Ta/uea, 0, (2-Hstara * Ha)*Ta/(uea**2), 0]
+ timeMatrix[2] = [0, 0, 1, 0, 0]
+ timeMatrix[3] = [-c*Ha, -c*Ta, 0, 1, 0]
+ if dataCont.flowRegime[iMarker]==1:
+ timeMatrix[4] = [0, 0, 0, 2*deltaa/(Usa * uea**2), 2*deltaa/(uea * Cta * Usa)]
+ else:
+ # Shear lag equation ignored in the laminar portion of the flow
+ timeMatrix[4] = [0, 0, 0, 0, 1]
- sysRes = np.linalg.solve(timeMatrix, spaceMatrix).flatten()
- return sysRes
-
-
- def __evalPoint(self,iMarker, bNodes, xtrIn, xtrFIn, nMarker) -> tuple:
- """Evaluates where the point is, the adjactent points
- and the state of transition on the corresponding side
-
- Args:
- ----
- - iMarker (integer): Id of the cell to evaluate.
-
- - bNodes (np.ndarray): Boundary nodes of the computational domain.
-
- - xtrIn (np.ndarray): Current transiton location on [upper, lower] sides.
-
- - xtrFIn (np.ndarray): Forced transition state on [upper, lower] sides.
-
- - nMarker ([type]): Number of cells in the computational domain.
-
- Returns:
- -------
- - tuple: Id of the adjacent points, name of the corresponding group ('airfoil', 'wake'),
- transition state on the corresponding side.
- """
- if iMarker == bNodes[0] + 1:
- iMarkerPrev = 0
- iMarkerPrev2 = None
-
- if iMarker == bNodes[1]: # First point of the lower side (next to the stagnation point).
- iMarkerPrev = 0
- iMarkerPrev2 = None
-
- elif iMarker == bNodes[1]+1:
- iMarkerPrev = iMarker-1
- iMarkerPrev2 = 0
- else:
- iMarkerPrev = iMarker-1
- iMarkerPrev2 = iMarker-2
-
- if iMarker == bNodes[1]-1 or iMarker == bNodes[2]-1 or iMarker == nMarker-1:
- iMarkerNext = None
- else: iMarkerNext = iMarker+1
-
- if iMarker >= bNodes[2]:
- name='wake'
- xtr=0
- xtrF=None
-
- elif iMarker < bNodes[1]:
- name = 'airfoil'
- xtr = xtrIn[0]
- xtrF = xtrFIn[0]
- elif bNodes[1] <= iMarker < bNodes[2]:
- name = 'airfoil'
- xtr = xtrIn[1]
- xtrF = xtrFIn[1]
- return iMarkerPrev2, iMarkerPrev, iMarkerNext, name, xtr, xtrF
+ sysRes = np.linalg.solve(timeMatrix, spaceMatrix).flatten()
+
+ """sysRes = np.empty(5)
+ if dataCont.flowRegime[iMarker] == 0: # Laminar case
+ sysRes[0] = (((-((c*Ha*Ta**2)/uea**2) - Ta/uea)*(-2*Cda + (Cfa*Hstara)/2. + Ta*dHstar_dx +
+ ((2*Hstar2a + (1 - Ha)*Hstara)*Ta*due_dx)/uea))/(-((c*Ha*Ta**2)/uea**3) - Ta/uea**2) +
+ (((2*c*Ta**2)/uea**2 - (c*Ha*Hstara*Ta**2)/uea**2 + Ta/uea - (Hstara*Ta)/uea)*
+ (-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/(-((c*Ha*Ta**2)/uea**3) - Ta/uea**2) +
+ (((2*Ta**2)/uea**3 - (Ha*Ta**2)/uea**3)*(cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) +
+ due_dx - dueExt_dx))/(-((c*Ha*Ta**2)/uea**3) - Ta/uea**2))
+
+ sysRes[1] = ((((c*Ha**2*Ta)/uea**2 + Ha/uea)*(-2*Cda + (Cfa*Hstara)/2. + Ta*dHstar_dx + ((2*Hstar2a + (1 - Ha)*Hstara)*Ta*due_dx)/uea))/
+ (-((c*Ha*Ta**2)/uea**3) - Ta/uea**2) + (((-2*c*Ha*Ta)/uea**2 + (c*Ha**2*Hstara*Ta)/uea**2 - 1/uea - Ha/uea + (Ha*Hstara)/uea)*
+ (-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/(-((c*Ha*Ta**2)/uea**3) - Ta/uea**2) +
+ ((-((Ha*Ta)/uea**3) + (Ha**2*Ta)/uea**3)*(cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) +
+ due_dx - dueExt_dx))/(-((c*Ha*Ta**2)/uea**3) - Ta/uea**2))
+
+ sysRes[2] = (-Axa + dN_dx)
+
+ sysRes[3] = (-((c*Ta*(-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/((-((c*Ha*Ta**2)/uea**3) - Ta/uea**2)*uea)) -
+ (Ta*(cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + due_dx -
+ dueExt_dx))/((-((c*Ha*Ta**2)/uea**3) - Ta/uea**2)*uea**2))
+
+ sysRes[4] = 0
+
+ elif dataCont.flowRegime[iMarker] == 1: # Turbulent case
+
+ sysRes[0] = ((((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**3*Usa) - (2*deltaa*Ta)/(Cta*uea**2*Usa))*
+ (-2*Cda + (Cfa*Hstara)/2. + Ta*dHstar_dx + ((2*Hstar2a + (1 - Ha)*Hstara)*Ta*due_dx)/uea))/
+ ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)) +
+ (((4*c*deltaa*Ta**2)/(Cta*uea**3*Usa) - (2*c*deltaa*Ha*Hstara*Ta**2)/(Cta*uea**3*Usa) + (2*deltaa*Ta)/(Cta*uea**2*Usa) -
+ (2*deltaa*Hstara*Ta)/(Cta*uea**2*Usa))*(-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/
+ ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)) +
+ (((4*deltaa*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa))*
+ (cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + due_dx - dueExt_dx))/
+ ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)))
+
+ sysRes[1] = ((((2*c*deltaa*Ha**2*Ta)/(Cta*uea**3*Usa) + (2*deltaa*Ha)/(Cta*uea**2*Usa))*
+ (-2*Cda + (Cfa*Hstara)/2. + Ta*dHstar_dx + ((2*Hstar2a + (1 - Ha)*Hstara)*Ta*due_dx)/uea))/
+ ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)) -
+ (2*deltaa*((2*c*Ha*Ta)/uea**2 - (c*Ha**2*Hstara*Ta)/uea**2 + 1/uea + Ha/uea - (Ha*Hstara)/uea)*
+ (-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/
+ (Cta*uea*((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa))*Usa) +
+ (((-2*deltaa*Ha*Ta)/(Cta*uea**4*Usa) + (2*deltaa*Ha**2*Ta)/(Cta*uea**4*Usa))*
+ (cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + due_dx - dueExt_dx))/
+ ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)))
+
+ sysRes[2] = 0
+
+ sysRes[3] = ((-2*c*deltaa*Ta*(-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/
+ (Cta*uea**2*((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa))*Usa) -
+ (2*deltaa*Ta*(cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + due_dx -
+ dueExt_dx))/(Cta*uea**3*((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa))*Usa))
+
+
+ sysRes[4] = ((2*c*deltaa*Ta*(-0.5*Cfa + dTheta_dx + ((2 + Ha - Mea**2)*Ta*due_dx)/uea))/
+ (uea**3*((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa))*Usa) +
+ ((-((c*Ha*Ta**2)/uea**3) - Ta/uea**2)*(-5.6*(Cteqa - Cta*dissipRatio) + (2*deltaa*dCt_dx)/Cta -
+ 2*deltaa*((4*(Cfa/2. - (0.02227667631989307*(-1 + Hka)**2)/(dissipRatio**2*Hka**2)))/(3.*Ha*Ta) - due_dx/uea)))/
+ ((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa)) +
+ (2*deltaa*Ta*(cExt*ddeltaStarExt_dx - c*(Ta*dH_dx + Ha*dTheta_dx) + due_dx -
+ dueExt_dx))/(uea**4*((-2*c*deltaa*Ha*Ta**2)/(Cta*uea**4*Usa) - (2*deltaa*Ta)/(Cta*uea**3*Usa))*Usa))"""
+
+
+ return sysRes
+
+
+ def __evalPoint(self, dataCont, iMarker) -> tuple:
+ """Evaluates where the point is, the adjactent points
+ and the state of transition on the corresponding side
+
+ Args:
+ ----
+ - iMarker (integer): Id of the cell to evaluate.
+
+ - bNodes (np.ndarray): Boundary nodes of the computational domain.
+
+ - xtrIn (np.ndarray): Current transiton location on [upper, lower] sides.
+
+ - xtrFIn (np.ndarray): Forced transition state on [upper, lower] sides.
+
+ - nMarker ([type]): Number of cells in the computational domain.
+
+ Returns:
+ -------
+ - tuple: Id of the adjacent points, name of the corresponding group ('airfoil', 'wake'),
+ transition state on the corresponding side.
+ """
+ iMarkerPrev = 0 if iMarker == dataCont.bNodes[1] else iMarker - 1
+
+ if iMarker < dataCont.bNodes[1]:
+ name='airfoil'
+ xtrF= dataCont.xtrF[0]
+
+ elif iMarker < dataCont.bNodes[2]:
+ name = 'airfoil'
+ xtrF = dataCont.xtrF[1]
+
+ else:
+ name = 'wake'
+ xtrF = None
+ return iMarkerPrev, name, xtrF
class sysMatrix(flowEquations):
- def __init__(self, _nMarker, _nMarkers, setGradients, fouScheme) -> None:
- flowEquations.__init__(self, setGradients, fouScheme)
- self._nMarker = _nMarker
- self._nMarkers = _nMarkers
- pass
+ def __init__(self, _nMarker, _nMarkers, setGradients, fouScheme) -> None:
+ flowEquations.__init__(self, setGradients, fouScheme)
+ self._nMarker = _nMarker
+ self._nMarkers = _nMarkers
+ pass
- def buildJacobian(self, dataCont, marker, sysRes, order = 1, eps = 1e-10) -> np.ndarray:
- """Contruct the Jacobian used for solving the liner system.
-
- Args:
- ----
- - dataCont (classType): Data container (solution, BL parameters, invicid flow parameters, geometry...)
+ def buildJacobian(self, dataCont, marker, sysRes, order = 1, eps = 1e-10) -> np.ndarray:
+ """Contruct the Jacobian used for solving the liner system.
+
+ Args:
+ ----
+ - dataCont (classType): Data container (solution, BL parameters, invicid flow parameters, geometry...)
- - marker (integer, optional): Cell id on the computational domain. Default : -1 for whole domain.
+ - marker (integer, optional): Cell id on the computational domain. Default : -1 for whole domain.
- - order (unsigned int, optional): Order of the discretization used to obtain the Jacobian (1st or 2nd order implemented).
- If order == 1, linSysRes has to be specifed
+ - order (unsigned int, optional): Order of the discretization used to obtain the Jacobian (1st or 2nd order implemented).
+ If order == 1, linSysRes has to be specifed
- - linSysRes (np.ndarray, optional): Residuals of the linear system at the current iteration. Used (and mandatory) if order == 1.
+ - linSysRes (np.ndarray, optional): Residuals of the linear system at the current iteration. Used (and mandatory) if order == 1.
- - eps (double, optional): Perturbation on the variables. Default : 10^(-8).
+ - eps (double, optional): Perturbation on the variables. Default : 10^(-8).
- Returns
- -------
+ Returns
+ -------
- - JacMatrix (np.ndarray): Jacobian of the system.
+ - JacMatrix (np.ndarray): Jacobian of the system.
- """
-
- # Build Jacobian for one point
- if order == 1 and sysRes is not None:
-
- JacMatrix = np.zeros((dataCont.nVar, dataCont.nVar))
- xUp = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
+ """
+
+ # Build Jacobian for one point
+ if order == 1 and sysRes is not None:
+
+ JacMatrix = np.zeros((dataCont.nVar, dataCont.nVar))
+ xUp = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
- for iVar in range(dataCont.nVar):
- varSave = xUp[iVar]
- xUp[iVar] += eps
+ for iVar in range(dataCont.nVar):
+ varSave = xUp[iVar]
+ xUp[iVar] += eps
- JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x=xUp) - sysRes) / eps
+ JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x=xUp) - sysRes) / eps
- xUp[iVar] = varSave
- return JacMatrix
+ xUp[iVar] = varSave
+ return JacMatrix
- elif order == 2:
- JacMatrix = np.zeros((dataCont.nVar, dataCont.nVar))
- xUp = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
- xDwn = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
+ elif order == 2:
+ JacMatrix = np.zeros((dataCont.nVar, dataCont.nVar))
+ xUp = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
+ xDwn = dataCont.u[marker*dataCont.nVar : (marker+1)*dataCont.nVar].copy()
- for iVar in range(dataCont.nVar):
- varSave = xUp[iVar]
- xUp[iVar] += eps
- xDwn[iVar] -= eps
+ for iVar in range(dataCont.nVar):
+ varSave = xUp[iVar]
+ xUp[iVar] += eps
+ xDwn[iVar] -= eps
- JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x = xUp) - self._blLaws(dataCont, marker, x = xDwn)) / (2*eps)
+ JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x = xUp) - self._blLaws(dataCont, marker, x = xDwn)) / (2*eps)
- xUp[iVar] = varSave
- xDwn[iVar] = varSave
- return JacMatrix
+ xUp[iVar] = varSave
+ xDwn[iVar] = varSave
+ return JacMatrix
- else:
- raise RuntimeError('Only first and second orders Jacobian can be computed.')
+ else:
+ raise RuntimeError('Only first and second orders Jacobian can be computed.')
- def buildResiduals(self, dataCont, marker) -> np.ndarray:
- """Construct residuals of the linear system.
-
- Args:
- ----
-
- - dataCont (classType): Data container (solution, BL parameters, invicid flow parameters, geometry...)
+ def buildResiduals(self, dataCont, marker) -> np.ndarray:
+ """Construct residuals of the linear system.
+
+ Args:
+ ----
+
+ - dataCont (classType): Data container (solution, BL parameters, invicid flow parameters, geometry...)
- - marker (int, optional): Cell id on the computational domain. default: -1 for whole domain.
- """
+ - marker (int, optional): Cell id on the computational domain. default: -1 for whole domain.
+ """
- # Build Residual for one point
- return self._blLaws(dataCont, marker, x=None)
\ No newline at end of file
+ # Build Residual for one point
+ return self._blLaws(dataCont, marker, x=dataCont.u[marker*dataCont.nVar:(marker+1)*dataCont.nVar])
\ No newline at end of file
diff --git a/dart/viscous/Solvers/DTransition.py b/dart/viscous/Solvers/DTransition.py
index e59b9a9598e9e3817a06340c83ae26d202224987..36f741fcfa21f81556299f74c6941bfd171c1002 100644
--- a/dart/viscous/Solvers/DTransition.py
+++ b/dart/viscous/Solvers/DTransition.py
@@ -48,56 +48,71 @@ class Transition:
def solveTransition(self,dataCont,transMarker):
- # Wake
- dataCont.flowRegime[dataCont.bNodes[2]:] = 1
-
- # Upper side
- if dataCont.bNodes[0]<transMarker<dataCont.bNodes[1]:
-
- # Set regime to turbulent on remaining upper side points
- dataCont.transMarkers[0] = transMarker
- dataCont.flowRegime[transMarker : dataCont.bNodes[1]] = 1
-
- # Compute upper side transition location
- dataCont.xtr[0] = (dataCont.Ncrit-dataCont.u[(transMarker-1)*dataCont.nVar+2])*((dataCont.xGlobal[transMarker]-dataCont.xGlobal[transMarker-1])/(dataCont.u[transMarker*dataCont.nVar+2]-dataCont.u[(transMarker-1)*dataCont.nVar+2]))+dataCont.xGlobal[transMarker-1]
- avTurb = (dataCont.xGlobal[transMarker]-dataCont.xtr[0])/(dataCont.xGlobal[dataCont.transMarkers[0]]-dataCont.xGlobal[dataCont.transMarkers[0]-1]) # percentage of the subsection corresponding to a turbulent flow
-
- # Fully turbulent flow on lower side
- elif transMarker == dataCont.bNodes[1]:
- dataCont.transMarkers[1] = dataCont.bNodes[1]
- dataCont.xtr[1] = 0
- dataCont.flowRegime[dataCont.bNodes[1] : dataCont.bNodes[2]] = 1
-
- # Lower side
- elif dataCont.bNodes[1] <= transMarker < dataCont.bNodes[2]:
-
- # Set regime to turbulent on remaining lower side points
- dataCont.transMarkers[1] = transMarker
- dataCont.flowRegime[transMarker : dataCont.bNodes[2]] = 1
-
- # Compute lower side transition location
- dataCont.xtr[1] = (dataCont.Ncrit - dataCont.u[(transMarker-1)*dataCont.nVar+2])*((dataCont.xGlobal[transMarker]-dataCont.xGlobal[transMarker-1])/(dataCont.u[transMarker*dataCont.nVar+2]-dataCont.u[(transMarker-1)*dataCont.nVar+2])) + dataCont.xGlobal[transMarker-1]
- avTurb = (dataCont.xGlobal[transMarker] - dataCont.xtr[1])/(dataCont.xGlobal[dataCont.transMarkers[1]]-dataCont.xGlobal[dataCont.transMarkers[1]-1]) # percentage of the subsection corresponding to a turbulent flow
-
- # Boundary condition
- avLam = 1- avTurb # percentage of the subsection corresponding to a laminar flow
- dataCont.blPar[(transMarker-1)*dataCont.nBlPar + 7] = Closures.turbulentClosure(dataCont.u[(transMarker-1)*dataCont.nVar+0], dataCont.u[(transMarker-1)*dataCont.nVar+1],0.03, dataCont.blPar[(transMarker-1)*dataCont.nBlPar+0], dataCont.extPar[(transMarker-1)*dataCont.nExtPar+0],'airfoil')[6]
- dataCont.u[(transMarker-1)*dataCont.nVar + 4] = 0.7*avLam*dataCont.blPar[(transMarker-1)*dataCont.nBlPar +7]
-
- return avTurb
-
+
+ dataCont.flowRegime[dataCont.bNodes[2]:] = 1 # Wake
+
+ # Upper side
+ if dataCont.bNodes[0]<transMarker<dataCont.bNodes[1]:
+
+ # Set regime to turbulent on remaining upper side points
+ dataCont.transMarkers[0] = transMarker
+ dataCont.flowRegime[transMarker : dataCont.bNodes[1]] = 1
+
+ # Compute upper side transition location
+ dataCont.xtr[0] = (dataCont.Ncrit-dataCont.u[(transMarker-1)*dataCont.nVar+2])*((dataCont.xGlobal[transMarker]-dataCont.xGlobal[transMarker-1])/(dataCont.u[transMarker*dataCont.nVar+2]-dataCont.u[(transMarker-1)*dataCont.nVar+2]))+dataCont.xGlobal[transMarker-1]
+ avTurb = (dataCont.xGlobal[transMarker]-dataCont.xtr[0])/(dataCont.xGlobal[dataCont.transMarkers[0]]-dataCont.xGlobal[dataCont.transMarkers[0]-1]) # percentage of the subsection corresponding to a turbulent flow
+
+ # Fully turbulent flow on lower side
+ elif transMarker == dataCont.bNodes[1]:
+ dataCont.transMarkers[1] = dataCont.bNodes[1]
+ dataCont.xtr[1] = 0
+ dataCont.flowRegime[dataCont.bNodes[1] : dataCont.bNodes[2]] = 1
+
+ # Lower side
+ elif dataCont.bNodes[1] <= transMarker < dataCont.bNodes[2]:
+
+ # Set regime to turbulent on remaining lower side points
+ dataCont.transMarkers[1] = transMarker
+ dataCont.flowRegime[transMarker : dataCont.bNodes[2]] = 1
+
+ # Compute lower side transition location
+ dataCont.xtr[1] = (dataCont.Ncrit - dataCont.u[(transMarker-1)*dataCont.nVar+2])*((dataCont.xGlobal[transMarker]-dataCont.xGlobal[transMarker-1])/(dataCont.u[transMarker*dataCont.nVar+2]-dataCont.u[(transMarker-1)*dataCont.nVar+2])) + dataCont.xGlobal[transMarker-1]
+ avTurb = (dataCont.xGlobal[transMarker] - dataCont.xtr[1])/(dataCont.xGlobal[dataCont.transMarkers[1]]-dataCont.xGlobal[dataCont.transMarkers[1]-1]) # percentage of the subsection corresponding to a turbulent flow
+
+ # Boundary condition
+ avLam = 1- avTurb # percentage of the subsection corresponding to a laminar flow
+ dataCont.blPar[(transMarker-1)*dataCont.nBlPar + 7] = Closures.turbulentClosure(dataCont.u[(transMarker-1)*dataCont.nVar+0], dataCont.u[(transMarker-1)*dataCont.nVar+1],0.03, dataCont.blPar[(transMarker-1)*dataCont.nBlPar+0], dataCont.extPar[(transMarker-1)*dataCont.nExtPar+0],'airfoil')[6]
+ dataCont.u[(transMarker-1)*dataCont.nVar + 4] = 0.7*avLam*dataCont.blPar[(transMarker-1)*dataCont.nBlPar +7]
+
+ return avTurb
+
+ def forcedTransition(self, dataCont, marker):
+
+ if marker < dataCont.bNodes[1]:
+ avTurb = (dataCont.xGlobal[marker]-dataCont.xtrF[0])/(dataCont.xGlobal[marker]-dataCont.xGlobal[marker-1]) # Percentage of the subsection corresponding to a turbulent flow
+ elif marker < dataCont.bNodes[2]:
+ avTurb = (dataCont.xGlobal[marker]-dataCont.xtrF[1])/(dataCont.xGlobal[marker]-dataCont.xGlobal[marker-1]) # Percentage of the subsection corresponding to a turbulent flow
+ else:
+ print('Trying to impose transition in the wake')
+ avLam = 1- avTurb # Percentage of the subsection corresponding to a laminar flow
+
+ # Boundary condition
+ dataCont.blPar[(marker-1)*dataCont.nBlPar + 7] = Closures.turbulentClosure(dataCont.u[(marker-1)*dataCont.nVar+0], dataCont.u[(marker-1)*dataCont.nVar+1],0.03, dataCont.blPar[(marker-1)*dataCont.nBlPar+0], dataCont.extPar[(marker-1)*dataCont.nExtPar+0],'airfoil')[6]
+ dataCont.u[(marker-1)*dataCont.nVar + 4] = 0.7*avLam*dataCont.blPar[(marker-1)*dataCont.nBlPar +7]
+
+ return avTurb
def transitionBC(self, dataCont, marker):
- if marker < dataCont.bNodes[1]:
- xtr = dataCont.xtr[0]
- elif dataCont.bNodes[1] <= marker < dataCont.bNodes[2]:
- xtr = dataCont.xtr[1]
-
- avTurb = (dataCont.xGlobal[marker]-dataCont.xtr[0])/(dataCont.xGlobal[marker]-dataCont.xGlobal[marker-1]) # Percentage of the subsection corresponding to a turbulent flow
- avLam = 1- avTurb # Percentage of the subsection corresponding to a laminar flow
-
- # Boundary condition
- dataCont.blPar[(marker-1)*dataCont.nBlPar + 7] = Closures.turbulentClosure(dataCont.u[(marker-1)*dataCont.nVar+0], dataCont.u[(marker-1)*dataCont.nVar+1],0.03, dataCont.blPar[(marker-1)*dataCont.nBlPar+0], dataCont.extPar[(marker-1)*dataCont.nExtPar+0],'airfoil')[6]
- dataCont.u[(marker-1)*dataCont.nVar + 4] = 0.7*avLam*dataCont.blPar[(marker-1)*dataCont.nBlPar +7]
-
- return avTurb
\ No newline at end of file
+ if marker < dataCont.bNodes[1]:
+ xtr = dataCont.xtr[0]
+ elif dataCont.bNodes[1] <= marker < dataCont.bNodes[2]:
+ xtr = dataCont.xtr[1]
+
+ avTurb = (dataCont.xGlobal[marker]-dataCont.xtr[0])/(dataCont.xGlobal[marker]-dataCont.xGlobal[marker-1]) # Percentage of the subsection corresponding to a turbulent flow
+ avLam = 1- avTurb # Percentage of the subsection corresponding to a laminar flow
+
+ # Boundary condition
+ dataCont.blPar[(marker-1)*dataCont.nBlPar + 7] = Closures.turbulentClosure(dataCont.u[(marker-1)*dataCont.nVar+0], dataCont.u[(marker-1)*dataCont.nVar+1],0.03, dataCont.blPar[(marker-1)*dataCont.nBlPar+0], dataCont.extPar[(marker-1)*dataCont.nExtPar+0],'airfoil')[6]
+ dataCont.u[(marker-1)*dataCont.nVar + 4] = 0.7*avLam*dataCont.blPar[(marker-1)*dataCont.nBlPar +7]
+
+ return avTurb
\ No newline at end of file
diff --git a/dart/viscous/Solvers/Steady/StdSolver.py b/dart/viscous/Solvers/Steady/StdSolver.py
index 5ae44a320e37c9caa6f0e4adbb4666dff968bcf9..e7cdd1902ddc11f0851757936783ee3f883e9b7d 100755
--- a/dart/viscous/Solvers/Steady/StdSolver.py
+++ b/dart/viscous/Solvers/Steady/StdSolver.py
@@ -482,6 +482,8 @@ class Solver:
blwVel, _, group.xx, blScalars, blVectors = self.__boundaryLayer(data, group)
print('WKblwVel',len(blwVel))
print('WKdeltaStarOut', len(blVectors[0]))
+ plt.plot(blwVel)
+ plt.show()
group.deltaStar = blVectors[0]
# The drag is measured at the end of the wake
self.Cd = blScalars[0]