From 2e275560dc35724502c227ec746da7ff244101ae Mon Sep 17 00:00:00 2001
From: Dechamps Paul <paul.dechamps@student.uliege.be>
Date: Thu, 18 Nov 2021 20:47:50 +0100
Subject: [PATCH] Time marching VII.
---
dart/viscous/Driver.py | 323 ++++++++++++++
dart/viscous/GroupConstructor.py | 129 ++++++
dart/viscous/PostProcess.py | 92 ++++
dart/viscous/README.md | 21 +
dart/viscous/{newton.py => Steady/Newton.py} | 0
dart/viscous/Steady/StdAirfoil.py | 131 ++++++
dart/viscous/Steady/StdBoundary.py | 36 ++
dart/viscous/Steady/StdCoupler.py | 82 ++++
.../{solver.py => Steady/StdSolver.py} | 11 +-
dart/viscous/Steady/StdWake.py | 78 ++++
dart/viscous/Validation.py | 94 +++++
dart/viscous/XfoilFiles/BlParM03A0Re6e6.dat | 186 +++++++++
dart/viscous/XfoilFiles/Case1FreeTrans.dat | 186 +++++++++
dart/viscous/XfoilFiles/case12Xfoil.dat | 186 +++++++++
dart/viscous/XfoilFiles/case15Xfoil.dat | 186 +++++++++
dart/viscous/__init__.py | 1 -
dart/viscous/airfoil.py | 5 +-
dart/viscous/boundary.py | 0
dart/viscous/coupler.py | 70 ++--
dart/viscous/model/BIConditions.py | 154 +++++++
dart/viscous/model/CSolver.py | 322 ++++++++++++++
dart/viscous/model/Closures.py | 234 +++++++++++
dart/viscous/model/DataStructures.py | 217 ++++++++++
dart/viscous/model/Discretization.py | 105 +++++
dart/viscous/model/TSolver.py | 395 ++++++++++++++++++
dart/viscous/model/Transition.py | 87 ++++
dart/viscous/model/Variables.py | 136 ++++++
dart/viscous/model/vplo.py | 87 ++++
dart/viscous/wake.py | 4 +-
29 files changed, 3526 insertions(+), 32 deletions(-)
create mode 100644 dart/viscous/Driver.py
create mode 100755 dart/viscous/GroupConstructor.py
create mode 100644 dart/viscous/PostProcess.py
create mode 100755 dart/viscous/README.md
rename dart/viscous/{newton.py => Steady/Newton.py} (100%)
mode change 100644 => 100755
create mode 100755 dart/viscous/Steady/StdAirfoil.py
create mode 100755 dart/viscous/Steady/StdBoundary.py
create mode 100755 dart/viscous/Steady/StdCoupler.py
rename dart/viscous/{solver.py => Steady/StdSolver.py} (98%)
mode change 100644 => 100755
create mode 100755 dart/viscous/Steady/StdWake.py
create mode 100755 dart/viscous/Validation.py
create mode 100644 dart/viscous/XfoilFiles/BlParM03A0Re6e6.dat
create mode 100644 dart/viscous/XfoilFiles/Case1FreeTrans.dat
create mode 100644 dart/viscous/XfoilFiles/case12Xfoil.dat
create mode 100644 dart/viscous/XfoilFiles/case15Xfoil.dat
delete mode 100644 dart/viscous/__init__.py
mode change 100644 => 100755 dart/viscous/airfoil.py
mode change 100644 => 100755 dart/viscous/boundary.py
mode change 100644 => 100755 dart/viscous/coupler.py
create mode 100755 dart/viscous/model/BIConditions.py
create mode 100644 dart/viscous/model/CSolver.py
create mode 100755 dart/viscous/model/Closures.py
create mode 100755 dart/viscous/model/DataStructures.py
create mode 100755 dart/viscous/model/Discretization.py
create mode 100644 dart/viscous/model/TSolver.py
create mode 100644 dart/viscous/model/Transition.py
create mode 100755 dart/viscous/model/Variables.py
create mode 100755 dart/viscous/model/vplo.py
mode change 100644 => 100755 dart/viscous/wake.py
diff --git a/dart/viscous/Driver.py b/dart/viscous/Driver.py
new file mode 100644
index 0000000..fc7429f
--- /dev/null
+++ b/dart/viscous/Driver.py
@@ -0,0 +1,323 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# Driver : Communication with coupler and initialization #
+# of some components of the viscous computation part #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+from dart.viscous.model.Discretization import Discretization
+from dart.viscous.model.CSolver import sysMatrix
+from dart.viscous.GroupConstructor import GroupConstructor
+from dart.viscous.model.Variables import Variables
+import dart.viscous.model.TSolver as TSolver
+import dart.viscous.model.Transition as Transition
+from dart.viscous.model.vplo import vplo
+from dart.viscous.PostProcess import PostProcess
+import dart.viscous.Validation as validation
+import dart.viscous.model.BIConditions as Conditions
+from fwk.coloring import ccolors
+import numpy as np
+import matplotlib.pyplot as plt
+import time
+import cProfile, pstats
+
+# ------------------------------------------------------------------------------
+# Driver : Input: - Data with the inviscid solution, geometry and user input for
+# time integration
+#
+# Output: - Corrected data ready to be used for viscous inviscd coupling.
+#
+# Sorts the data, allows the whole computational domain to be treated at the same
+# time (equations are solved on the three groups simultaneously) and computes
+# required data for the inviscid solver (blowing velocity aerodynamic drag,...)
+# ------------------------------------------------------------------------------
+class Driver:
+ def __init__(self,_Re, user_xtr, _timeParams, _case, _mapMesh):
+ self.Re=_Re
+ self.Cd=0
+ self.it=0
+ self.uPrevious=None
+ self.blParPrevious=None
+ self.xtrF=user_xtr
+ self.mapMesh = _mapMesh
+ self.Ti=None
+ self.timeParams=_timeParams
+ self.res=None
+ self.upper=0
+ self.case=_case
+ self.xtr=[1,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, group):
+ '''Primary component of the Diver. Receives 2 groups : 'airfoil' and 'wake'.
+ Sorts data, initialize code components and calls the time tintegrator. Then sorts the data again,
+ communicating it to the coupler.'''
+
+ # ###################################### Driver.run() ######################################
+ # ------------------------------------------------------------------------------------------
+ # Input: -'group': - List of 2 elements (airfoil/wake) containing parameters
+ # from the inviscid solver.
+ #
+ # Tasks: - Initializing: - var: Data container used throughout the code.
+ # Contains the mesh, the solution, boundary
+ # layer parameters etc.
+ # - DiricheletBC: Used to compute the boudary conditions
+ # of the airfoil (Twaithes) and the wake.
+ # - gradComp: Gradient computation using various spacial
+ # schemes.
+ # - cSolver: Evaluates the boundary layer laws at one point
+ # of the computational domain given a prescibed
+ # solution or not. Uses the closures models.
+ # - initializer: Provides the initial condition at the first
+ # iteration
+ # - tSolver: Time integration given a model, CFL and tolerences.
+ # Used the same way for explicit/implicit time
+ # marching.
+ #
+ # - Running: Runs the time integration.
+ # Tries with lower CFL if the solver crashes
+ # - Sorting: (Post Process): Calls the PostProcess class to sort the value and,
+ # from the solution, compute/correct parameters that must
+ # sent to the inviscid solver
+ # Output: - Void.
+ # Directly modifes the group that were initialy sent.
+ # ------------------------------------------------------------------------------------------
+ nFactor = 3
+ # Merge upper, lower sides and wake
+ dataIn=[None,None]
+ for n in range(0, len(group)):
+ group[n].connectListNodes, group[n].connectListElems, dataIn[n] = group[n].connectList()
+ param, xx, xGlobal, bNodes, data, dx, initMesh, bNodesInitMesh = GroupConstructor().mergeVectors(dataIn, self.mapMesh, nFactor)
+
+ # Initialize variable container
+ plotter=vplo(self.case)
+ var=Variables(data, param, xx, xGlobal, bNodes, self.xtrF, self.Ti, self.Re)
+ if self.it!=0:
+ var.extPar[4::var.nExtPar]=data[:,8]
+
+ DirichletBC=Conditions.Boundaries(self.Re) # Boundary condition (Airfoil/Wake).
+ gradComp=Discretization(var.xx, var.extPar[4::var.nExtPar], var.bNodes) # Gradients computation.
+
+ cSolver=sysMatrix(var.nN, var.nVar, gradComp.fou, gradComp.fou) # Spacial operator and Jacobian computation using the boundary layer laws.
+ transSolver = Transition.Transition(var.xtrF)
+ # Initialize time Integrator
+ tSolver=TSolver.implicitEuler(cSolver, transSolver, DirichletBC.wakeBC, self.timeParams['CFL0'])
+
+ print('+------------------------------ Preprocessing ---------------------------------+')
+ if self.mapMesh:
+ print('| - Mesh mapped from {:<.0f} to {:<.0f} elements.{:>38s}'.format(len(initMesh), var.nN, '|'))
+ print('| - Number of elements: {:<.0f}. {:>51s}'.format(var.nN,'|'))
+ print('| - Maximum Mach number: {:<.2f}. {:>49s}'.format(np.max((var.extPar[0::var.nExtPar])),'|'))
+ if var.xtrF[0] is None:
+ print('| - Free transition on the upper side. {:>41}'.format('|'))
+ else:
+ print('| - Forced transition on the upper side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[0], '|'))
+ if var.xtrF[1] is None:
+ print('| - Free transition on the lower side. {:>41}'.format('|'))
+ else:
+ print('| - Forced transition on the lower side; x/c = {:<.4f}. {:>25s}'.format(var.xtrF[1], '|'))
+ print('| - Critical amplification ratio: {:<.2f}. {:>40s}'.format(var.Ncrit,'|'))
+
+ # Initial condition: 'Driver' stores the solution form the previous coupling iteration.
+ if self.uPrevious is None or len(self.uPrevious) != len(var.u): # Typically for the first iteration
+ initializer=Conditions.Initial(DirichletBC) # Initalize boundary condition
+ initializer.initialConditions(var) # Set boundary condition
+
+ tSolver.setCFL0(0.5)
+ tSolver.setmaxIt(20000)
+
+ print('| - Using Twaithes solution as initial guess. {:>34}'.format('|'))
+ #tSolver.tol=1e-2
+ else: # If we use a solution previously obtained.
+ var.u=self.uPrevious.copy()
+ var.u[2::var.nVar] = 0
+ print('| - Using previous solution as initial guess. {:>34}'.format('|'))
+
+ print('|{:>79}'.format('|'))
+ print('| Numerical parameters {:>57}'.format('|'))
+ print('| - CFL0: {:<.2f} {:>65}'.format(tSolver.getCFL0(),'|'))
+ print('| - Tolerance: {:<.0f} {:>62}'.format(np.log10(tSolver.gettol()),'|'))
+ print('| - Solver maximum iterations: {:<.0f} {:>45}'.format(tSolver.getmaxIt(),'|'))
+ print('|{:>79}'.format('|'))
+
+ """plt.figure(1)
+ plt.plot(var.xGlobal[:var.bNodes[1]],var.extPar[2:var.bNodes[1]*var.nExtPar:var.nExtPar])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.extPar[var.bNodes[1]*var.nExtPar + 2:var.bNodes[2]*var.nExtPar:var.nExtPar])
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.extPar[var.bNodes[2]*var.nExtPar + 2::var.nExtPar])
+ plt.title('Inv Vel IN')
+ plt.xlim([0,1])
+
+ plt.figure(2)
+ plt.plot(var.xGlobal[:var.bNodes[1]],var.extPar[0:var.bNodes[1]*var.nExtPar:var.nExtPar])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.extPar[var.bNodes[1]*var.nExtPar + 0:var.bNodes[2]*var.nExtPar:var.nExtPar])
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.extPar[var.bNodes[2]*var.nExtPar + 0::var.nExtPar])
+ plt.title('Inv Mach')
+ plt.xlim([0,1])
+ plt.show()"""
+
+ print('+------------------------------- Solver begin ---------------------------------+')
+ # Run the time integration
+ nTries=1
+ convergedPoints = tSolver.timeSolve(var)
+ print('+-------------------------------- Solver Exit ---------------------------------+')
+
+ if np.any(convergedPoints!=1):
+ print('|', ccolors.ANSI_YELLOW + 'Some points did not converge.', ccolors.ANSI_RESET, '{:>55s}'.format('|'))
+
+ elif np.all(convergedPoints==1):
+ print('|', ccolors.ANSI_GREEN + 'All points converged.', ccolors.ANSI_RESET, '{:>55s}'.format('|'))
+ print('|{:>79}'.format('|'))
+
+
+ # Save solution to use it as initial condition the next iteration
+ var.u[4::var.nVar][var.flowRegime == 0] = 0.001
+ self.uPrevious=var.u.copy()
+ self.blParPrevious=var.blPar.copy()
+
+ """plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[0:var.bNodes[1]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 0:var.bNodes[2]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 0::var.nVar])
+ plt.show()
+ plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[1:var.bNodes[1]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 1:var.bNodes[2]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 1::var.nVar])
+ plt.show()
+ plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[2:var.bNodes[1]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 2:var.bNodes[2]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 2::var.nVar])
+ plt.show()
+ plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[3:var.bNodes[1]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 3:var.bNodes[2]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 3::var.nVar])
+ plt.show()
+ plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[4:var.bNodes[1]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 4:var.bNodes[2]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 4::var.nVar])
+ plt.show()"""
+
+ self.blScalars, AFblwVel, WKblwVel, AFblVectors, WKblVectors = PostProcess().sortParam(var, self.mapMesh, initMesh, bNodesInitMesh, nFactor)
+ self.xtr=var.xtr.copy()
+ self.blVectors=AFblVectors
+
+ group[0].TEnd = [AFblVectors[1][var.bNodes[1]-1], AFblVectors[1][var.bNodes[2]-1]]
+ group[0].HEnd = [AFblVectors[2][var.bNodes[1]-1], AFblVectors[2][var.bNodes[2]-1]]
+ group[0].CtEnd = [AFblVectors[8][var.bNodes[1]-1], AFblVectors[8][var.bNodes[2]-1]]
+ # Delete the stagnation point doublon for variables in the VII loop
+ group[0].deltaStar = AFblVectors[0]
+ group[0].xx = np.concatenate((var.xx[:var.bNodes[1]],var.xx[var.bNodes[1]:var.bNodes[2]]))
+
+ # Sort the following results in reference frame
+ group[0].deltaStar = group[0].deltaStar[group[0].connectListNodes.argsort()]
+ group[0].xx = group[0].xx[group[0].connectListNodes.argsort()]
+ group[0].u = AFblwVel[group[0].connectListElems.argsort()]
+ self.data=data[:var.bNodes[2],:]
+
+ group[1].deltaStar = WKblVectors[0]
+ # The 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
+
+ group[1].xx=var.xx[var.bNodes[2]:]
+ # Sort the following results in reference frame
+ group[1].deltaStar = group[1].deltaStar[group[1].connectListNodes.argsort()]
+ group[1].xx = group[1].xx[group[1].connectListNodes.argsort()]
+ group[1].u = WKblwVel[group[1].connectListElems.argsort()]
+
+ """self.writeFile()
+ val=validation.Validation('Case1FreeTrans.dat')
+ val.main(1,self.wData)"""
+
+
+ """plt.plot(var.xGlobal[:var.bNodes[1]], var.u[3:var.bNodes[1]*var.nVar:var.nVar])
+ plt.show()"""
+
+ """"plt.plot(var.xGlobal[0:var.bNodes[1]],var.u[3:var.bNodes[1]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.u[var.bNodes[1]*var.nVar + 3:var.bNodes[2]*var.nVar:var.nVar])
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar + 3::var.nVar])
+ plt.title('Ue')
+ plt.show()
+
+ plt.plot(var.xGlobal[0:var.bNodes[1]],var.blPar[9:var.bNodes[1]*var.nBlPar:var.nBlPar])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.blPar[var.bNodes[1]*var.nBlPar + 9:var.bNodes[2]*var.nBlPar:var.nBlPar])
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.blPar[var.bNodes[2]*var.nBlPar + 9::var.nBlPar])
+ plt.title('deltaStar')
+ plt.show()
+
+ plt.plot(var.xGlobal[:var.bNodes[1]],var.extPar[1:var.bNodes[1]*var.nExtPar:var.nExtPar])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],var.extPar[var.bNodes[1]*var.nExtPar + 1:var.bNodes[2]*var.nExtPar:var.nExtPar])
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.extPar[var.bNodes[2]*var.nExtPar + 1::var.nExtPar])
+ plt.title('rhoe')
+ plt.show()
+
+ plt.plot(var.xGlobal[:var.bNodes[1]],AFblwVel[:var.bNodes[1]])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],AFblwVel[var.bNodes[1]:var.bNodes[2]])
+ plt.plot(var.xGlobal[var.bNodes[1]:var.bNodes[2]],AFblwVel[var.bNodes[1]:var.bNodes[2]])
+ plt.title('BlwVel')
+ plt.show()"""
\ No newline at end of file
diff --git a/dart/viscous/GroupConstructor.py b/dart/viscous/GroupConstructor.py
new file mode 100755
index 0000000..b6fd945
--- /dev/null
+++ b/dart/viscous/GroupConstructor.py
@@ -0,0 +1,129 @@
+################################################################################
+# #
+# Flow viscous module file #
+# Constructor : Merges and sorts data #
+# of the 3 groups (upper, lower sides and wake) #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+import math
+from matplotlib import pyplot as plt
+
+class GroupConstructor():
+ """ ---------- Group Constructor ----------
+
+ Builds single vectors elements out of the 3 groups given in input.
+ """
+ def __init__(self) -> None:
+ pass
+
+ def mergeVectors(self,dataIn, mapMesh, nFactor):
+ '''Merges 3 groups upper/lower sides and wake.'''
+
+ initMesh = np.concatenate((dataIn[0][0][:,0], dataIn[0][1][1:,0], dataIn[1][:,0]))
+ bNodesInitMesh = [0, len(dataIn[0][0][:,0]), len(dataIn[0][0][:,0])+len(dataIn[0][1][1:,0])]
+ for k in range(len(dataIn)):
+
+ if len(dataIn[k])==2: # Airfoil case
+ if mapMesh:
+ dataUpper = self.refineMesh(dataIn[k][0], nFactor)
+ dataLower = self.refineMesh(dataIn[k][1], nFactor)
+ else:
+ dataUpper = dataIn[k][0]
+ dataLower = dataIn[k][1]
+
+
+ xx_up, dv_up, vtExt_up, _, alpha_up= self.__getBLcoordinates(dataUpper)
+ xx_lw, dv_lw, vtExt_lw, _, alpha_lw= self.__getBLcoordinates(dataLower)
+ else: # Wake case
+
+ if mapMesh:
+ dataWake = self.refineMesh(dataIn[k], nFactor)
+ else:
+ dataWake = dataIn[k]
+
+ xx_wk, dv_wk, vtExt_wk, _, alpha_wk= self.__getBLcoordinates(dataWake)
+
+ xx_wk += xx_up[-1]
+ # Delete stagnation point doublon
+ xx_lw=np.delete(xx_lw,0)
+ dv_lw=np.delete(dv_lw,0)
+ vtExt_lw=np.delete(vtExt_lw,0)
+ alpha_lw=np.delete(alpha_lw,0)
+
+ # Nodes at the boundary of the computational domain: [Stagnation point, First lower side point, First wake point]
+ boundaryNodes=np.empty(3, dtype=int)
+ boundaryNodes[0] = 0
+ boundaryNodes[1] = len(xx_up)
+ boundaryNodes[2] = len(xx_up)+len(xx_lw)
+
+ param={'dv':None,'vtExt':None,'alpha':None}
+ param['dv']=np.concatenate((dv_up,dv_lw,dv_wk))
+ param['vtExt']=np.concatenate((vtExt_up,vtExt_lw,vtExt_wk))
+ param['alpha']=np.concatenate((alpha_up,alpha_lw,alpha_wk))
+
+ xx=np.concatenate((xx_up,xx_lw,xx_wk))
+
+ data=np.zeros((len(xx),9))
+ data=np.concatenate((dataUpper,dataLower,dataWake),axis=0)
+ data=np.delete(data,boundaryNodes[1],0) # Delete stagnation point doublon
+
+ xGlobal = data[:,0]
+
+ dx=np.zeros(len(xx))
+ for i in range(1,len(xx)):
+ if i==boundaryNodes[1]:
+ dx[i]=xx[i]-xx[0]
+ elif i==boundaryNodes[2]:
+ dx[i]=0
+ else:
+ dx[i]=xx[i]-xx[i-1]
+ return param, xx, xGlobal, boundaryNodes, data, dx, initMesh, bNodesInitMesh
+
+
+ def __getBLcoordinates(self,data):
+ """Transform the reference coordinates into airfoil coordinates
+ """
+ nN = len(data[:,0])
+ nE = nN-1 # Nbr of element along the surface
+ vt = np.zeros(nN)
+ xx = np.zeros(nN) # Position along the surface of the airfoil
+ dx = np.zeros(nE) # Distance along two nodes
+ dv = np.zeros(nE) # Speed difference according to element
+ alpha = np.zeros(nE) # Angle of the element for the rotation matrix
+ for i in range(0,nE):
+ alpha[i] = np.arctan2((data[i+1,1]-data[i,1]),(data[i+1,0]-data[i,0]))
+ dx[i] = np.sqrt((data[i+1,0]-data[i,0])**2+(data[i+1,1]-data[i,1])**2)
+ xx[i+1] = dx[i]+xx[i]
+ vt[i] = (data[i,3]*np.cos(alpha[i]) + data[i,4]*np.sin(alpha[i])) # Tangential speed at node 1 element i
+ vt[i+1] = (data[i+1,3]*np.cos(alpha[i]) + data[i+1,4]*np.sin(alpha[i])) # Tangential speed at node 2 element i
+ dv[i] = (vt[i+1] - vt[i])/dx[i] # Velocity gradient according to element i. central scheme with half point
+ return xx, dv, vt, nN, alpha
+
+
+ def refineMesh(self, inMesh, nFactor):
+
+ outMesh = np.empty(len(inMesh[0,:]))
+
+ for marker in range(len(inMesh[:,0]) - 1):
+ dCell = inMesh[marker + 1, :] - inMesh[marker, :]
+
+ for n in range(nFactor):
+
+ outMesh = np.ma.row_stack((outMesh, inMesh[marker, :] + n*dCell/nFactor))
+
+ outMesh = np.delete(outMesh, 0, 0)
+
+ # Add last cell
+ outMesh = np.ma.row_stack((outMesh, inMesh[-1, :]))
+
+ return outMesh
\ No newline at end of file
diff --git a/dart/viscous/PostProcess.py b/dart/viscous/PostProcess.py
new file mode 100644
index 0000000..b3c11a7
--- /dev/null
+++ b/dart/viscous/PostProcess.py
@@ -0,0 +1,92 @@
+from matplotlib import pyplot as plt
+import numpy as np
+class PostProcess:
+ def __init__(self):
+ pass
+
+ def sortParam(self,var, mapMesh, initMesh, bNodesInitMesh, nFactor):
+
+ # Recompute deltaStar.
+
+ var.blPar[9::var.nBlPar] = var.u[0::var.nVar] * var.u[1::var.nVar]
+
+ # Compute blowing velocity on each cell.
+ blwVel = np.empty(var.nN - 1)
+
+ for i in range(1, var.nN): # -1 and we will remove the first wake point after.
+ iPrev = 0 if i == var.bNodes[1] else i-1
+
+ # blwVel[i-1] = (rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1]))
+
+ if i == var.bNodes[2]:
+ sep = i-1 # Separation between airfoil and wake. (-1 because the first wake point will be removed)
+ continue
+ else:
+ blwVel[i-1] = (var.extPar[i*var.nExtPar+1]*var.u[i*var.nVar+3] * var.blPar[i*var.nBlPar+9]
+ - var.extPar[iPrev*var.nExtPar+1] * var.u[iPrev*var.nVar+3] * var.blPar[iPrev*var.nBlPar+9]) / (var.extPar[i*var.nExtPar+1] * (var.xx[i] - var.xx[iPrev]))
+
+ blwVel = np.delete(blwVel, var.bNodes[2] - 1)
+
+ if mapMesh:
+ # Map blowing velocities to inviscid solver mesh.
+
+ invBlwVel = np.zeros(len(initMesh)-2)
+ for invMarker in range(len(invBlwVel)):
+ if invMarker == bNodesInitMesh[2]:
+ sep = invMarker - 1
+ invBlwVel[invMarker] = 1/4*(blwVel[nFactor*invMarker - 1] + 2* blwVel[nFactor*invMarker] + blwVel[nFactor*invMarker + 1])
+
+ AFblwVel = invBlwVel[:sep]
+ WKblwVel = invBlwVel[sep:]
+
+ else:
+ # Separate airfoil and wake blowing velocities.
+ AFblwVel=blwVel[:sep] # -1 because the first wake point was removed and indices have moved.
+ WKblwVel=blwVel[sep:]
+ """plt.figure(1)
+ plt.plot(AFblwVel)
+ plt.figure(2)
+ plt.plot(WKblwVel)
+ plt.show()"""
+
+ # 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 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]
+
+ # 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]]
+
+ return blScalars, AFblwVel, WKblwVel, blVectors_airfoil, blVectors_wake
diff --git a/dart/viscous/README.md b/dart/viscous/README.md
new file mode 100755
index 0000000..c2eed4b
--- /dev/null
+++ b/dart/viscous/README.md
@@ -0,0 +1,21 @@
+# viscous
+A dart module solving the boundary layer equations to perform viscous-inviscid coupling with the main full potential solver.
+
+
+
+## Main features
+* Solvers
+ - [x] Steady: Solves the steady boundary layer equations
+ - [x] Unsteady: Solves the unsteady boundary layer equations using (pseudo-)time marching procedures
+* Time-Marching
+ - [x] Explicit/Implicit methods
+ - [x] Newton method with direct linear solver for implicit time integration
+* Various
+ - [x] Different spacial discretizations (FOU, SOU, CENTERED)
+ - [x] Forced/Free transition
+ - [x] Free-stream turbulence intensity (affects transition location)
+ - [x] Adaptative CFL number
+
+## Sub-modules
+* [Steady](flow/viscous/Steady): Steady equations solver and components
+* [Model](flow/viscous/model): Mathematical model of the time dependent boundary layer equations
\ No newline at end of file
diff --git a/dart/viscous/newton.py b/dart/viscous/Steady/Newton.py
old mode 100644
new mode 100755
similarity index 100%
rename from dart/viscous/newton.py
rename to dart/viscous/Steady/Newton.py
diff --git a/dart/viscous/Steady/StdAirfoil.py b/dart/viscous/Steady/StdAirfoil.py
new file mode 100755
index 0000000..18db9d5
--- /dev/null
+++ b/dart/viscous/Steady/StdAirfoil.py
@@ -0,0 +1,131 @@
+#!/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.
+
+
+## Airfoil around which the boundary layer is computed
+# Amaury Bilocq
+
+from dart.viscous.Steady.StdBoundary import Boundary
+import numpy as np
+
+class Airfoil(Boundary):
+ def __init__(self, _boundary):
+ Boundary.__init__(self, _boundary)
+ self.name = 'airfoil' # type of boundary
+ self.T0 = 0 # initial condition for the momentum thickness
+ self.H0 = 0
+ self.n0 = 0
+ self.Ct0 = 0
+ self.TEnd = [0,0]
+ self.HEnd = [0,0]
+ self.CtEnd = [0,0]
+ self.nEnd = [0,0]
+
+ def initialConditions(self, Re, dv):
+ if dv > 0:
+ self.T0 = np.sqrt(0.075/(Re*dv))
+ else:
+ self.T0 = 1e-6
+ self.H0 = 2.23 # One parameter family Falkner Skan
+ self.n0 = 0
+ self.Ct0 = 0
+ return self.T0, self.H0, self.n0, self.Ct0
+
+ def connectList(self):
+ ''' Sort the value read by the viscous solver/ Create list of connectivity
+ '''
+ N1 = np.zeros(self.nN, dtype=int) # Node number
+ connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
+ connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
+ data = np.zeros((self.boundary.nodes.size(), 10))
+ i = 0
+ for n in self.boundary.nodes:
+ data[i,0] = n.no
+ data[i,1] = n.pos[0]
+ data[i,2] = n.pos[1]
+ data[i,3] = n.pos[2]
+ data[i,4] = self.v[i,0]
+ data[i,5] = self.v[i,1]
+ data[i,6] = self.Me[i]
+ data[i,7] = self.rhoe[i]
+ data[i,8] = self.deltaStar[i]
+ data[i,9] = self.xx[i]
+ i += 1
+ # Table containing the element and its nodes
+ eData = np.zeros((self.nE,3), dtype=int)
+ for i in range(0, self.nE):
+ eData[i,0] = self.boundary.tag.elems[i].no
+ 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
+ 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)
+ idxTE = np.where(data[:,1] == np.max(data[:,1]))[0] # TE nodes
+ idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], eData[:,2] == data[idxTE[0],0]))[0] # TE element 1
+ idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], eData[:,2] == data[idxTE[1],0]))[0] # TE element 2
+ ye0 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe0[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe0[0],2])[0],2]) # y coordinates of TE element 1 CG
+ ye1 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe1[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe1[0],2])[0],2]) # y coordinates of TE element 2 CG
+ if ye0 - ye1 > 0: # element 1 containing TE node 1 is above element 2
+ upperGlobTE = int(data[idxTE[0],0])
+ lowerGlobTE = int(data[idxTE[1],0])
+ else:
+ upperGlobTE = int(data[idxTE[1],0])
+ lowerGlobTE = int(data[idxTE[0],0])
+ # Connectivity
+ connectListElems[0] = np.where(eData[:,1] == globStag)[0]
+ N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes -> globStag
+ connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
+ i = 1
+ upperTE = 0
+ lowerTE = 0
+ # Sort the suction part
+ while upperTE == 0:
+ N1[i] = eData[connectListElems[i-1],2] # Second node of the element
+ connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes
+ connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
+ if eData[connectListElems[i],2] == upperGlobTE:
+ upperTE = 1
+ i += 1
+ # Sort the pressure side
+ connectListElems[i] = np.where(eData[:,2] == globStag)[0]
+ connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
+ N1[i] = eData[connectListElems[i],2]
+ while lowerTE == 0:
+ N1[i+1] = eData[connectListElems[i],1] # First node of the element
+ connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes
+ connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
+ if eData[connectListElems[i+1],1] == lowerGlobTE:
+ lowerTE = 1
+ i += 1
+ connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
+ data[:,0] = data[connectListNodes,0]
+ data[:,1] = data[connectListNodes,1]
+ data[:,2] = data[connectListNodes,2]
+ data[:,3] = data[connectListNodes,3]
+ data[:,4] = data[connectListNodes,4]
+ data[:,5] = data[connectListNodes,5]
+ data[:,6] = data[connectListNodes,6]
+ data[:,7] = data[connectListNodes,7]
+ data[:,8] = data[connectListNodes,8]
+ data[:,9] = data[connectListNodes,9]
+ # Separated upper and lower part
+ data = np.delete(data,0,1)
+ uData = data[0:np.argmax(data[:,0])+1]
+ lData = data[np.argmax(data[:,0])+1:None]
+ lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
+ return connectListNodes, connectListElems,[uData, lData]
diff --git a/dart/viscous/Steady/StdBoundary.py b/dart/viscous/Steady/StdBoundary.py
new file mode 100755
index 0000000..f3ff33a
--- /dev/null
+++ b/dart/viscous/Steady/StdBoundary.py
@@ -0,0 +1,36 @@
+#!/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.
+
+
+## Base class representing a physical boundary
+# Amaury Bilocq
+
+import numpy as np
+
+class Boundary:
+ def __init__(self, _boundary):
+ self.boundary = _boundary # boundary of interest
+ self.nN = len(self.boundary.nodes) # number of nodes
+ self.nE = len(self.boundary.tag.elems) # number of elements
+ self.v = np.zeros((self.nN, 3)) # velocity at edge of BL
+ self.u = np.zeros(self.nE) # blowing velocity
+ self.Me = np.zeros(self.nN) # Mach number at the edge of the boundary layer
+ self.rhoe = np.zeros(self.nN) # density at the edge of the boundary layer
+ self.connectListnodes = np.zeros(self.nN) # connectivity for nodes
+ self.connectListElems = np.zeros(self.nE) # connectivity for elements
+ self.deltaStar = np.zeros(self.nN) # displacement thickness
+ self.xx = np.zeros(self.nN) # coordinates in the reference of the physical surface. Used to store the values for the interaction law
diff --git a/dart/viscous/Steady/StdCoupler.py b/dart/viscous/Steady/StdCoupler.py
new file mode 100755
index 0000000..9427e71
--- /dev/null
+++ b/dart/viscous/Steady/StdCoupler.py
@@ -0,0 +1,82 @@
+#!/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.
+
+
+## Viscous-inviscid coupler (quasi-simultaneous coupling)
+# Amaury Bilocq
+
+import numpy as np
+from fwk.coloring import ccolors
+from dart.viscous.Steady.StdAirfoil import Airfoil
+from dart.viscous.Steady.StdWake import Wake
+
+class Coupler:
+ def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer):
+ self.isolver =_isolver # inviscid solver
+ self.vsolver = _vsolver # viscous solver
+ self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects
+ self.tol = _tol # tolerance of the VII
+ self.writer = _writer
+
+ def run(self):
+ ''' Perform coupling
+ '''
+ # initialize loop
+ it = 0
+ converged = 0 # temp
+ CdOld = self.vsolver.Cd
+ while converged == 0:
+ print(ccolors.ANSI_BLUE + 'Iteration: ', it, ccolors.ANSI_RESET)
+ # run inviscid solver
+ self.isolver.run()
+ print('--- Viscous solver parameters ---')
+ print('Tolerance (drag count):', self.tol * 1000)
+ print('--- Viscous problem definition ---')
+ print('Reynolds number:', self.vsolver.Re)
+ print('')
+ for n in range(0, len(self.group)):
+ print('Computing for', self.group[n].name, '...', end = ' ')
+ for i in range(0, len(self.group[n].boundary.nodes)):
+ self.group[n].v[i,0] = self.isolver.U[self.group[n].boundary.nodes[i].row][0]
+ self.group[n].v[i,1] = self.isolver.U[self.group[n].boundary.nodes[i].row][1]
+ self.group[n].v[i,2] = self.isolver.U[self.group[n].boundary.nodes[i].row][2]
+ self.group[n].Me[i] = self.isolver.M[self.group[n].boundary.nodes[i].row]
+ self.group[n].rhoe[i] = self.isolver.rho[self.group[n].boundary.nodes[i].row]
+ # run viscous solver
+ self.vsolver.run(self.group[n])
+ if self.group[n].name == 'airfoil':
+ self.group[n+1].T0 = self.group[n].TEnd[0]+self.group[n].TEnd[1]
+ self.group[n+1].H0 = (self.group[n].HEnd[0]*self.group[n].TEnd[0]+self.group[n].HEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0
+ self.group[n+1].Ct0 = (self.group[n].CtEnd[0]*self.group[n].TEnd[0]+self.group[n].CtEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0
+ # get blowing velocity from viscous solver and update inviscid solver
+ for i in range(0, self.group[n].nE):
+ self.group[n].boundary.setU(i, self.group[n].u[i])
+ print('done!')
+ print(' Iter Cd xtr_top xtr_bot Error')
+ print('{0:8d} {1:12.5f} {2:12.5f} {3:12.5f} {4:12.5f}'.format(it, self.vsolver.Cd, self.vsolver.xtr[0], self.vsolver.xtr[1], abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd))
+ print('')
+ # Converged or not
+ if abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd < self.tol:
+ converged = 1
+ else:
+ CdOld = self.vsolver.Cd
+ it += 1
+ self.vsolver.it += 1
+ # save results
+ self.isolver.save(self.writer)
+ self.vsolver.writeFile()
+ print('\n')
diff --git a/dart/viscous/solver.py b/dart/viscous/Steady/StdSolver.py
old mode 100644
new mode 100755
similarity index 98%
rename from dart/viscous/solver.py
rename to dart/viscous/Steady/StdSolver.py
index 9b96228..922bfa3
--- a/dart/viscous/solver.py
+++ b/dart/viscous/Steady/StdSolver.py
@@ -39,9 +39,10 @@
# -- check transonic predictions
# -- implement quasi-2D method for 3D flows
-from dart.viscous.boundary import Boundary
-from dart.viscous.wake import Wake
-import dart.viscous.newton as Newton
+#from dart.viscous.Steady.StdBoundary import Boundary
+#from dart.viscous.Steady.StdWake import Wake
+from matplotlib import pyplot as plt
+import dart.viscous.Steady.Newton as Newton
import numpy as np
from fwk.coloring import ccolors
@@ -78,6 +79,7 @@ class Solver:
wData['delta'] = self.blVectors[9]
wData['xtrTop'] = self.xtr[0]
wData['xtrBot'] = self.xtr[1]
+ self.wData = wData
return wData
def writeFile(self):
@@ -436,6 +438,8 @@ class Solver:
Cd[i-1], Cf[i-1], Hstar[i-1], Hstar2[i-1], Hk[i-1], delta[i-1] = self.__laminarClosure(theta[i-1], H[i-1], Ret[i-1],Me[i-1], group.name)
Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i] = self.__turbulentClosure(theta[i], H[i],Ct[i], Ret[i],Me[i], group.name)
deltaStar[i] = H[i]*theta[i]
+ if group.name=='wake':
+ print('dx = ',xx[i]-xx[i-1])
blwVel[i-1] = (rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1]))
# Normalize the friction and dissipation coefficients by the freestream velocity
Cf = vt**2*Cf
@@ -447,6 +451,7 @@ class Solver:
# Outputs
blScalars = [Cdtot, Cdf, Cdp]
blVectors = [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta]
+ print(len(deltaStar))
return blwVel, xtr, xx, blScalars, blVectors
diff --git a/dart/viscous/Steady/StdWake.py b/dart/viscous/Steady/StdWake.py
new file mode 100755
index 0000000..625b3e4
--- /dev/null
+++ b/dart/viscous/Steady/StdWake.py
@@ -0,0 +1,78 @@
+#!/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.
+
+
+## Wake behind airfoil (around which the boundary layer is computed)
+# Amaury Bilocq
+
+from dart.viscous.Steady.StdBoundary import Boundary
+from dart.viscous.Steady.StdAirfoil import Airfoil
+import numpy as np
+
+class Wake(Boundary):
+ def __init__(self, _boundary):
+ Boundary.__init__(self, _boundary)
+ self.name = 'wake' # type of boundary
+ self.T0 = 0 # initial condition for the momentum thickness
+ self.H0 = 0
+ self.n0 = 9 # wake is always turbulent
+ self.Ct0 = 0
+
+
+ def connectList(self):
+ ''' Sort the value read by the viscous solver/ Create list of connectivity
+ '''
+ N1 = np.zeros(self.nN, dtype=int) # Node number
+ connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
+ connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
+ data = np.zeros((self.boundary.nodes.size(), 10))
+ i = 0
+ for n in self.boundary.nodes:
+ data[i,0] = n.no
+ data[i,1] = n.pos[0]
+ data[i,2] = n.pos[1]
+ data[i,3] = n.pos[2]
+ data[i,4] = self.v[i,0]
+ data[i,5] = self.v[i,1]
+ data[i,6] = self.Me[i]
+ data[i,7] = self.rhoe[i]
+ data[i,8] = self.deltaStar[i]
+ data[i,9] = self.xx[i]
+ i += 1
+ # Table containing the element and its nodes
+ eData = np.zeros((self.nE,3), dtype=int)
+ for i in range(0, self.nE):
+ eData[i,0] = self.boundary.tag.elems[i].no
+ eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
+ eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
+ connectListNodes = data[:,1].argsort()
+ connectListElems[0] = np.where(eData[:,1] == min(data[:,1]))[0]
+ for i in range(1, len(eData[:,0])):
+ connectListElems[i] = np.where(eData[:,1] == eData[connectListElems[i-1],2])[0]
+ data[:,0] = data[connectListNodes,0]
+ data[:,1] = data[connectListNodes,1]
+ data[:,2] = data[connectListNodes,2]
+ data[:,3] = data[connectListNodes,3]
+ data[:,4] = data[connectListNodes,4]
+ data[:,5] = data[connectListNodes,5]
+ data[:,6] = data[connectListNodes,6]
+ data[:,7] = data[connectListNodes,7]
+ data[:,8] = data[connectListNodes,8]
+ data[:,9] = data[connectListNodes,9]
+ # Separated upper and lower part
+ data = np.delete(data,0,1)
+ return connectListNodes, connectListElems, data
diff --git a/dart/viscous/Validation.py b/dart/viscous/Validation.py
new file mode 100755
index 0000000..dc2f13f
--- /dev/null
+++ b/dart/viscous/Validation.py
@@ -0,0 +1,94 @@
+import numpy as np
+import matplotlib.pyplot as plt
+class Validation:
+ def __init__(self,_xfoilFileName = None):
+ self.xfoilFileName=_xfoilFileName
+ def main(self, dartcl, wData):
+ # Load xfoil data
+
+ xFoilIn = 0
+ cl_xfoil = 0
+ cd_xfoil = 0
+
+ if self.xfoilFileName is not None:
+
+ xFoilIn = 1
+ with open('/Users/pauldechamps/lab/dartflo/dart/viscous/XfoilFiles/' + self.xfoilFileName, 'r') as f:
+ cl_xfoil=f.readline()
+ cd_xfoil=f.readline()
+ f.close()
+ # columns : x,ue,theta,Cf,H
+ data=np.loadtxt('/Users/pauldechamps/lab/dartflo/dart/viscous/XfoilFiles/'+self.xfoilFileName,skiprows=3,usecols=(1,3,4,5,6,7))
+
+ # columns : x, delta*, theta, H, Cf, ue
+ data[:,[1,2,3,4,5]]=data[:,[2,3,5,4,2]]
+ data[:,]
+
+ # Find upper/lower side separation
+ for i in range(len(data[:,0])):
+ if data[i,0]<=data[i+1,0]:
+ xStag=i
+ break
+ for i in range(len(wData['x'])):
+ if wData['x'][i]==1:
+ xTeDart=i+1
+ break
+
+ # 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), '|'))
+
+ # Plot variables
+ plot_dStar=plt.figure(1)
+ plt.plot(wData['x'][:xTeDart],wData['delta*'][:xTeDart],'b-',label=r'DARTflo upper')
+ plt.plot(wData['x'][xTeDart:],wData['delta*'][xTeDart:],'r-',label=r'DARTflo lower')
+
+ if xFoilIn:
+ plt.plot(data[:xStag,0],data[:xStag,1],'b--',label=r'Xfoil upper')
+ plt.plot(data[xStag:,0],data[xStag:,1],'r--',label=r'Xfoil lower')
+ #plt.plot(data[xWake:,0],data[xWake:,1],'r--',label=r'Xfoil Wake')
+ plt.legend()
+ plt.xlim([0,1])
+ plt.xlabel('$x/c$')
+ plt.ylabel('$\delta^*$')
+ """ # Theta
+ plot_Theta=plt.figure(2)
+ plt.plot(wData['x'][:xTeDart], wData['theta'][:xTeDart],'b-',label=r'DARTflo upper')
+ plt.plot(wData['x'][xTeDart:], wData['theta'][xTeDart:],'r-',label=r'DARTflo lower')
+ plt.plot(data[:xStag,0],data[:xStag,2],'b--',label=r'Xfoil upper')
+ plt.plot(data[xStag:,0],data[xStag:,2],'r--',label=r'Xfoil lower')
+ #plt.plot(data[xWake:,0],data[xWake:,2],'r--',label=r'Xfoil Wake')
+ plt.legend()
+ plt.xlim([0,1])
+ plt.xlabel('$x/c$')
+ plt.ylabel('$\theta$')"""
+ # H
+ plot_H=plt.figure(3)
+ plt.plot(wData['x'][:xTeDart],wData['H'][:xTeDart],'b-',label=r'DARTflo upper')
+ plt.plot(wData['x'][xTeDart:],wData['H'][xTeDart:],'r-',label=r'DARTflo lower')
+ plt.axvline(x=0.0138)
+ if xFoilIn:
+ plt.plot(data[:xStag,0],data[:xStag,3],'b--',label=r'Xfoil upper')
+ plt.plot(data[xStag:,0],data[xStag:,3],'r--',label=r'Xfoil lower')
+ #plt.plot(data[xWake:,0],data[xWake:,3],'r--',label=r'Xfoil Wake')
+ plt.legend()
+ plt.xlim([0,1])
+ plt.xlabel('$x/c$')
+ plt.ylabel('$H$')
+ # Cf
+ plot_Cf=plt.figure(4)
+ plt.plot(wData['x'][:xTeDart], wData['Cf'][:xTeDart],'b-',label=r'DARTflo upper')
+ plt.plot(wData['x'][xTeDart:], wData['Cf'][xTeDart:],'r-',label=r'DARTflo lower')
+
+ if xFoilIn:
+ plt.plot(data[:xStag,0],data[:xStag,4],'b--',label=r'Xfoil upper')
+ plt.plot(data[xStag:,0],data[xStag:,4],'r--',label=r'Xfoil lower')
+ #plt.plot(data[xWake:,0],data[xWake:,4],'r--',label=r'Xfoil Wake')
+ plt.legend()
+ plt.xlim([0,1])
+ plt.xlabel('$x/c$')
+ plt.ylabel('$C_f$')
+
+ plt.show()
+
diff --git a/dart/viscous/XfoilFiles/BlParM03A0Re6e6.dat b/dart/viscous/XfoilFiles/BlParM03A0Re6e6.dat
new file mode 100644
index 0000000..81e9b54
--- /dev/null
+++ b/dart/viscous/XfoilFiles/BlParM03A0Re6e6.dat
@@ -0,0 +1,186 @@
+0
+0.00666
+# s x y Ue/Vinf Dstar Theta Cf H H* P m K tau Di
+ 0.00000 1.00000 0.00126 0.88651 0.004714 0.003017 0.001226 1.5541 1.6929 0.00237 0.00418 0.00356
+ 0.00840 0.99168 0.00242 0.89902 0.004398 0.002865 0.001335 1.5266 1.7029 0.00232 0.00395 0.00354
+ 0.01982 0.98037 0.00398 0.91583 0.004019 0.002676 0.001486 1.4935 1.7156 0.00224 0.00368 0.00353
+ 0.03304 0.96727 0.00576 0.93146 0.003709 0.002511 0.001625 1.4679 1.7258 0.00218 0.00345 0.00350
+ 0.04772 0.95272 0.00771 0.94569 0.003453 0.002370 0.001752 1.4480 1.7340 0.00212 0.00327 0.00348
+ 0.06337 0.93720 0.00974 0.95828 0.003244 0.002249 0.001864 1.4330 1.7404 0.00207 0.00311 0.00345
+ 0.07959 0.92112 0.01182 0.96934 0.003072 0.002147 0.001961 1.4217 1.7453 0.00202 0.00298 0.00341
+ 0.09610 0.90474 0.01389 0.97909 0.002927 0.002058 0.002047 1.4131 1.7492 0.00197 0.00287 0.00338
+ 0.11275 0.88821 0.01594 0.98779 0.002803 0.001979 0.002123 1.4063 1.7522 0.00193 0.00277 0.00334
+ 0.12949 0.87160 0.01796 0.99564 0.002693 0.001908 0.002191 1.4011 1.7546 0.00189 0.00268 0.00330
+ 0.14626 0.85494 0.01994 1.00282 0.002595 0.001844 0.002253 1.3968 1.7565 0.00185 0.00260 0.00327
+ 0.16305 0.83827 0.02189 1.00945 0.002505 0.001785 0.002310 1.3935 1.7581 0.00182 0.00253 0.00323
+ 0.17984 0.82158 0.02381 1.01564 0.002423 0.001729 0.002364 1.3907 1.7594 0.00178 0.00246 0.00319
+ 0.19665 0.80488 0.02569 1.02146 0.002346 0.001677 0.002414 1.3885 1.7605 0.00175 0.00240 0.00315
+ 0.21346 0.78817 0.02753 1.02699 0.002274 0.001627 0.002462 1.3867 1.7614 0.00172 0.00234 0.00310
+ 0.23027 0.77146 0.02934 1.03228 0.002205 0.001580 0.002508 1.3852 1.7621 0.00168 0.00228 0.00306
+ 0.24707 0.75475 0.03111 1.03736 0.002140 0.001534 0.002553 1.3840 1.7627 0.00165 0.00222 0.00302
+ 0.26388 0.73803 0.03284 1.04228 0.002077 0.001490 0.002596 1.3830 1.7632 0.00162 0.00216 0.00298
+ 0.28068 0.72132 0.03453 1.04706 0.002017 0.001447 0.002639 1.3822 1.7637 0.00159 0.00211 0.00293
+ 0.29748 0.70460 0.03619 1.05172 0.001958 0.001406 0.002681 1.3816 1.7640 0.00155 0.00206 0.00288
+ 0.31427 0.68789 0.03781 1.05630 0.001901 0.001365 0.002722 1.3812 1.7643 0.00152 0.00201 0.00284
+ 0.33105 0.67118 0.03938 1.06080 0.001845 0.001325 0.002763 1.3809 1.7645 0.00149 0.00196 0.00279
+ 0.34783 0.65447 0.04092 1.06523 0.001790 0.001286 0.002804 1.3806 1.7647 0.00146 0.00191 0.00274
+ 0.36460 0.63777 0.04241 1.06962 0.001737 0.001248 0.002846 1.3805 1.7648 0.00143 0.00186 0.00269
+ 0.38135 0.62108 0.04387 1.07398 0.001684 0.001210 0.002888 1.3805 1.7649 0.00140 0.00181 0.00265
+ 0.39809 0.60440 0.04527 1.07831 0.001633 0.001173 0.002930 1.3805 1.7650 0.00136 0.00176 0.00260
+ 0.41482 0.58772 0.04663 1.08262 0.001582 0.001136 0.002972 1.3806 1.7650 0.00133 0.00171 0.00254
+ 0.43154 0.57106 0.04794 1.08691 0.001532 0.001100 0.003015 1.3808 1.7651 0.00130 0.00167 0.00249
+ 0.44823 0.55441 0.04920 1.09120 0.001482 0.001064 0.003059 1.3811 1.7651 0.00127 0.00162 0.00244
+ 0.46491 0.53778 0.05041 1.09548 0.001434 0.001029 0.003104 1.3814 1.7650 0.00123 0.00157 0.00239
+ 0.48157 0.52116 0.05156 1.09976 0.001385 0.000994 0.003150 1.3817 1.7650 0.00120 0.00152 0.00233
+ 0.49820 0.50456 0.05265 1.10404 0.001338 0.000959 0.003197 1.3821 1.7649 0.00117 0.00148 0.00228
+ 0.51481 0.48798 0.05368 1.10833 0.001291 0.000925 0.003245 1.3826 1.7649 0.00114 0.00143 0.00222
+ 0.53140 0.47143 0.05465 1.11261 0.001244 0.000891 0.003294 1.3831 1.7648 0.00110 0.00138 0.00217
+ 0.54796 0.45489 0.05555 1.11689 0.001198 0.000858 0.003345 1.3836 1.7647 0.00107 0.00134 0.00211
+ 0.56448 0.43839 0.05639 1.12116 0.001152 0.000825 0.003397 1.3842 1.7646 0.00104 0.00129 0.00205
+ 0.58098 0.42191 0.05714 1.12543 0.001107 0.000792 0.003450 1.3849 1.7644 0.00100 0.00125 0.00199
+ 0.59744 0.40546 0.05782 1.12969 0.001062 0.000759 0.003505 1.3856 1.7643 0.00097 0.00120 0.00193
+ 0.61386 0.38905 0.05842 1.13393 0.001018 0.000727 0.003562 1.3864 1.7641 0.00094 0.00115 0.00187
+ 0.63024 0.37268 0.05893 1.13815 0.000974 0.000696 0.003621 1.3873 1.7639 0.00090 0.00111 0.00181
+ 0.64657 0.35635 0.05935 1.14233 0.000931 0.000664 0.003681 1.3882 1.7637 0.00087 0.00106 0.00175
+ 0.66286 0.34007 0.05967 1.14647 0.000888 0.000633 0.003743 1.3892 1.7635 0.00083 0.00102 0.00168
+ 0.67910 0.32383 0.05989 1.15056 0.000846 0.000602 0.003808 1.3904 1.7633 0.00080 0.00097 0.00162
+ 0.69527 0.30766 0.06000 1.15457 0.000804 0.000572 0.003875 1.3916 1.7631 0.00076 0.00093 0.00155
+ 0.71139 0.29154 0.06000 1.15851 0.000762 0.000542 0.003944 1.3929 1.7628 0.00073 0.00088 0.00149
+ 0.72743 0.27550 0.05988 1.16234 0.000721 0.000512 0.004016 1.3944 1.7625 0.00069 0.00084 0.00142
+ 0.74340 0.25953 0.05963 1.16605 0.000681 0.000483 0.004090 1.3960 1.7622 0.00066 0.00079 0.00135
+ 0.75928 0.24366 0.05924 1.16961 0.000641 0.000454 0.004168 1.3977 1.7619 0.00062 0.00075 0.00128
+ 0.77506 0.22788 0.05871 1.17300 0.000601 0.000425 0.004249 1.3997 1.7615 0.00059 0.00071 0.00121
+ 0.79074 0.21222 0.05803 1.17617 0.000562 0.000397 0.004333 1.4019 1.7611 0.00055 0.00066 0.00114
+ 0.80628 0.19670 0.05718 1.17908 0.000524 0.000369 0.004421 1.4043 1.7607 0.00051 0.00062 0.00107
+ 0.82167 0.18135 0.05616 1.18169 0.000486 0.000342 0.004513 1.4069 1.7603 0.00048 0.00057 0.00099
+ 0.83688 0.16619 0.05497 1.18392 0.000448 0.000315 0.004610 1.4099 1.7598 0.00044 0.00053 0.00092
+ 0.85186 0.15127 0.05358 1.18570 0.000412 0.000288 0.004711 1.4133 1.7594 0.00041 0.00049 0.00085
+ 0.86655 0.13666 0.05200 1.18693 0.000376 0.000262 0.004818 1.4171 1.7589 0.00037 0.00045 0.00077
+ 0.88087 0.12246 0.05022 1.18747 0.000341 0.000237 0.004931 1.4213 1.7585 0.00033 0.00040 0.00070
+ 0.89470 0.10877 0.04824 1.18717 0.000307 0.000213 0.005049 1.4261 1.7581 0.00030 0.00036 0.00063
+ 0.90789 0.09575 0.04609 1.18586 0.000274 0.000190 0.005171 1.4315 1.7577 0.00027 0.00032 0.00056
+ 0.92028 0.08358 0.04381 1.18334 0.000243 0.000168 0.005297 1.4375 1.7575 0.00023 0.00029 0.00049
+ 0.93169 0.07242 0.04142 1.17942 0.000215 0.000147 0.005422 1.4441 1.7575 0.00021 0.00025 0.00043
+ 0.94200 0.06239 0.03900 1.17393 0.000189 0.000129 0.005546 1.4513 1.7577 0.00018 0.00022 0.00037
+ 0.95118 0.05354 0.03659 1.16673 0.000166 0.000113 0.005664 1.4589 1.7582 0.00015 0.00019 0.00031
+ 0.95926 0.04581 0.03424 1.15767 0.000146 0.000098 0.005776 1.4666 1.7591 0.00013 0.00017 0.00027
+ 0.96635 0.03910 0.03196 1.14656 0.000128 0.000086 0.005888 1.4736 1.7609 0.00011 0.00015 0.00023
+ 0.97255 0.03329 0.02976 1.13312 0.000112 0.000075 0.005999 1.4795 1.7636 0.00010 0.00013 0.00019
+ 0.97802 0.02826 0.02764 1.11701 0.000098 0.000065 0.006068 1.4883 1.7659 0.00008 0.00011 0.00016
+ 0.98286 0.02387 0.02560 1.09802 0.000087 0.000057 0.005960 1.5127 1.7614 0.00007 0.00010 0.00013
+ 0.98718 0.02003 0.02361 1.07580 0.000079 0.000050 0.005938 1.5610 1.7470 0.00006 0.00009 0.00011
+ 0.99107 0.01666 0.02167 1.04912 0.000072 0.000044 0.006114 1.6243 1.7295 0.00005 0.00008 0.00009
+ 0.99460 0.01369 0.01977 1.01676 0.000066 0.000039 0.006275 1.6874 1.7139 0.00004 0.00007 0.00007
+ 0.99781 0.01108 0.01789 0.97760 0.000059 0.000034 0.006368 1.7562 1.6986 0.00003 0.00006 0.00005
+ 1.00078 0.00879 0.01602 0.92996 0.000054 0.000029 0.006411 1.8257 1.6827 0.00003 0.00005 0.00004
+ 1.00352 0.00678 0.01415 0.87187 0.000048 0.000025 0.006354 1.8999 1.6671 0.00002 0.00004 0.00003
+ 1.00607 0.00505 0.01227 0.80117 0.000043 0.000022 0.006113 1.9846 1.6508 0.00001 0.00003 0.00002
+ 1.00847 0.00358 0.01038 0.71544 0.000038 0.000019 0.005681 2.0743 1.6353 0.00001 0.00003 0.00001
+ 1.01072 0.00236 0.00848 0.61305 0.000035 0.000016 0.004920 2.1769 1.6194 0.00001 0.00002 0.00001
+ 1.01286 0.00140 0.00657 0.49349 0.000033 0.000014 0.003973 2.2491 1.6092 0.00000 0.00002 0.00000
+ 1.01490 0.00070 0.00466 0.35814 0.000031 0.000014 0.003087 2.2392 1.6195 0.00000 0.00001 0.00000
+ 1.01685 0.00024 0.00276 0.21164 0.000030 0.000013 0.001909 2.2314 1.6208 0.00000 0.00001 0.00000
+ 1.01872 0.00003 0.00091 0.06128 0.000029 0.000013 0.000562 2.2295 1.6211 0.00000 0.00000 0.00000
+ 1.02053 0.00003 -0.00091 -0.08727 0.000029 0.000013 0.000800 2.2295 1.6211 0.00000 -0.00000 -0.00000
+ 1.02240 0.00024 -0.00276 -0.23741 0.000030 0.000013 0.002139 2.2321 1.6207 0.00000 -0.00001 -0.00000
+ 1.02434 0.00070 -0.00466 -0.38339 0.000031 0.000014 0.003300 2.2398 1.6194 0.00000 -0.00001 -0.00000
+ 1.02638 0.00140 -0.00657 -0.51792 0.000033 0.000014 0.004167 2.2497 1.6179 0.00000 -0.00002 -0.00000
+ 1.02852 0.00236 -0.00848 -0.63666 0.000035 0.000016 0.004691 2.2637 1.6157 0.00001 -0.00002 -0.00001
+ 1.03078 0.00358 -0.01038 -0.73791 0.000038 0.000017 0.004917 2.2791 1.6134 0.00001 -0.00003 -0.00001
+ 1.03317 0.00505 -0.01227 -0.82233 0.000042 0.000018 0.004910 2.2961 1.6108 0.00001 -0.00003 -0.00002
+ 1.03573 0.00678 -0.01415 -0.89181 0.000047 0.000020 0.004754 2.3136 1.6083 0.00002 -0.00004 -0.00002
+ 1.03847 0.00879 -0.01602 -0.94869 0.000051 0.000022 0.004506 2.3313 1.6058 0.00002 -0.00005 -0.00003
+ 1.04143 0.01108 -0.01789 -0.99525 0.000057 0.000024 0.004214 2.3489 1.6034 0.00002 -0.00006 -0.00004
+ 1.04465 0.01369 -0.01977 -1.03346 0.000062 0.000026 0.003906 2.3661 1.6011 0.00003 -0.00006 -0.00005
+ 1.04817 0.01666 -0.02167 -1.06496 0.000069 0.000029 0.003599 2.3829 1.5989 0.00003 -0.00007 -0.00006
+ 1.05206 0.02003 -0.02361 -1.09104 0.000075 0.000031 0.003302 2.3992 1.5968 0.00004 -0.00008 -0.00006
+ 1.05638 0.02387 -0.02560 -1.11274 0.000083 0.000034 0.003021 2.4150 1.5948 0.00004 -0.00009 -0.00007
+ 1.06122 0.02826 -0.02764 -1.13083 0.000091 0.000037 0.002757 2.4303 1.5929 0.00005 -0.00010 -0.00009
+ 1.06669 0.03329 -0.02976 -1.14592 0.000099 0.000040 0.002511 2.4453 1.5911 0.00005 -0.00011 -0.00010
+ 1.07290 0.03910 -0.03196 -1.15845 0.000109 0.000044 0.002283 2.4598 1.5894 0.00006 -0.00013 -0.00011
+ 1.07998 0.04581 -0.03424 -1.16876 0.000119 0.000048 0.002073 2.4740 1.5878 0.00007 -0.00014 -0.00012
+ 1.08806 0.05354 -0.03659 -1.17708 0.000131 0.000052 0.001879 2.4879 1.5862 0.00007 -0.00015 -0.00013
+ 1.09724 0.06239 -0.03900 -1.18358 0.000143 0.000057 0.001703 2.5015 1.5847 0.00008 -0.00017 -0.00015
+ 1.10756 0.07242 -0.04142 -1.18842 0.000156 0.000062 0.001543 2.5149 1.5833 0.00009 -0.00019 -0.00016
+ 1.11897 0.08358 -0.04381 -1.19174 0.000170 0.000067 0.001400 2.5281 1.5819 0.00009 -0.00020 -0.00018
+ 1.13135 0.09575 -0.04609 -1.19372 0.000184 0.000072 0.001274 2.5411 1.5805 0.00010 -0.00022 -0.00019
+ 1.14454 0.10877 -0.04824 -1.19454 0.000199 0.000077 0.001163 2.5540 1.5792 0.00011 -0.00024 -0.00021
+ 1.15837 0.12246 -0.05022 -1.19441 0.000214 0.000083 0.001065 2.5670 1.5779 0.00012 -0.00026 -0.00022
+ 1.17269 0.13666 -0.05200 -1.19348 0.000230 0.000088 0.000979 2.5801 1.5766 0.00013 -0.00027 -0.00024
+ 1.18738 0.15127 -0.05358 -1.19192 0.000245 0.000094 0.000903 2.5934 1.5753 0.00013 -0.00029 -0.00025
+ 1.20236 0.16619 -0.05497 -1.18984 0.000260 0.000099 0.000835 2.6070 1.5740 0.00014 -0.00031 -0.00026
+ 1.21757 0.18135 -0.05616 -1.18734 0.000276 0.000105 0.000774 2.6210 1.5727 0.00015 -0.00033 -0.00028
+ 1.23296 0.19670 -0.05718 -1.18450 0.000291 0.000110 0.000719 2.6354 1.5714 0.00015 -0.00035 -0.00029
+ 1.24851 0.21222 -0.05803 -1.18137 0.000307 0.000115 0.000668 2.6504 1.5701 0.00016 -0.00036 -0.00030
+ 1.26418 0.22788 -0.05871 -1.17801 0.000323 0.000120 0.000622 2.6658 1.5688 0.00017 -0.00038 -0.00031
+ 1.27996 0.24366 -0.05924 -1.17446 0.000339 0.000125 0.000580 2.6819 1.5674 0.00017 -0.00040 -0.00032
+ 1.29584 0.25953 -0.05963 -1.17076 0.000355 0.000131 0.000541 2.6985 1.5661 0.00018 -0.00042 -0.00033
+ 1.31181 0.27550 -0.05988 -1.16693 0.000371 0.000136 0.000505 2.7156 1.5647 0.00018 -0.00043 -0.00034
+ 1.32785 0.29154 -0.06000 -1.16300 0.000387 0.000141 0.000471 2.7333 1.5633 0.00019 -0.00045 -0.00035
+ 1.34397 0.30766 -0.06000 -1.15901 0.000404 0.000146 0.000440 2.7513 1.5620 0.00020 -0.00047 -0.00035
+ 1.36015 0.32383 -0.05989 -1.15499 0.000421 0.000151 0.000411 2.7694 1.5606 0.00020 -0.00049 -0.00036
+ 1.37638 0.34007 -0.05967 -1.15097 0.000437 0.000156 0.000385 2.7869 1.5594 0.00021 -0.00050 -0.00037
+ 1.39267 0.35635 -0.05935 -1.14704 0.000453 0.000161 0.000363 2.8021 1.5583 0.00021 -0.00052 -0.00038
+ 1.40900 0.37268 -0.05893 -1.14339 0.000468 0.000165 0.000347 2.8100 1.5578 0.00022 -0.00054 -0.00039
+ 1.42538 0.38905 -0.05842 -1.14063 0.000477 0.000170 0.000347 2.7936 1.5589 0.00022 -0.00054 -0.00039
+ 1.44181 0.40546 -0.05782 -1.13555 0.000444 0.000176 0.000791 2.5009 1.5351 0.00023 -0.00050 -0.00040
+ 1.45826 0.42191 -0.05714 -1.12271 0.000342 0.000196 0.002821 1.7309 1.6589 0.00025 -0.00038 -0.00046
+ 1.47476 0.43839 -0.05639 -1.11903 0.000330 0.000221 0.004106 1.4834 1.7358 0.00028 -0.00037 -0.00054
+ 1.49129 0.45489 -0.05555 -1.11627 0.000358 0.000251 0.004437 1.4125 1.7622 0.00031 -0.00040 -0.00062
+ 1.50784 0.47143 -0.05465 -1.11313 0.000399 0.000283 0.004375 1.3974 1.7671 0.00035 -0.00044 -0.00069
+ 1.52443 0.48798 -0.05368 -1.10956 0.000444 0.000315 0.004228 1.3963 1.7664 0.00039 -0.00049 -0.00076
+ 1.54104 0.50456 -0.05265 -1.10573 0.000490 0.000347 0.004079 1.3978 1.7647 0.00042 -0.00054 -0.00083
+ 1.55768 0.52116 -0.05156 -1.10175 0.000535 0.000379 0.003946 1.3992 1.7633 0.00046 -0.00059 -0.00089
+ 1.57433 0.53778 -0.05041 -1.09768 0.000580 0.000411 0.003831 1.3998 1.7623 0.00050 -0.00064 -0.00096
+ 1.59101 0.55441 -0.04920 -1.09356 0.000625 0.000443 0.003731 1.3997 1.7617 0.00053 -0.00068 -0.00102
+ 1.60771 0.57106 -0.04794 -1.08939 0.000670 0.000474 0.003641 1.3990 1.7615 0.00056 -0.00073 -0.00108
+ 1.62442 0.58772 -0.04663 -1.08520 0.000714 0.000506 0.003560 1.3981 1.7614 0.00060 -0.00077 -0.00114
+ 1.64115 0.60440 -0.04527 -1.08098 0.000759 0.000538 0.003486 1.3970 1.7615 0.00063 -0.00082 -0.00120
+ 1.65789 0.62108 -0.04387 -1.07673 0.000803 0.000571 0.003416 1.3958 1.7617 0.00066 -0.00086 -0.00125
+ 1.67465 0.63777 -0.04241 -1.07245 0.000848 0.000603 0.003350 1.3946 1.7618 0.00069 -0.00091 -0.00131
+ 1.69141 0.65447 -0.04092 -1.06814 0.000893 0.000636 0.003288 1.3935 1.7620 0.00073 -0.00095 -0.00137
+ 1.70819 0.67118 -0.03938 -1.06377 0.000939 0.000669 0.003228 1.3925 1.7622 0.00076 -0.00100 -0.00142
+ 1.72497 0.68789 -0.03781 -1.05935 0.000986 0.000703 0.003169 1.3916 1.7623 0.00079 -0.00104 -0.00147
+ 1.74177 0.70460 -0.03619 -1.05486 0.001033 0.000737 0.003112 1.3908 1.7624 0.00082 -0.00109 -0.00152
+ 1.75856 0.72132 -0.03453 -1.05028 0.001080 0.000771 0.003056 1.3902 1.7625 0.00085 -0.00113 -0.00157
+ 1.77536 0.73803 -0.03284 -1.04559 0.001129 0.000806 0.003001 1.3898 1.7625 0.00088 -0.00118 -0.00162
+ 1.79217 0.75475 -0.03111 -1.04078 0.001179 0.000842 0.002945 1.3896 1.7624 0.00091 -0.00123 -0.00167
+ 1.80898 0.77146 -0.02934 -1.03580 0.001230 0.000879 0.002890 1.3896 1.7622 0.00094 -0.00127 -0.00172
+ 1.82579 0.78817 -0.02753 -1.03063 0.001283 0.000916 0.002833 1.3898 1.7619 0.00097 -0.00132 -0.00177
+ 1.84259 0.80488 -0.02569 -1.02523 0.001338 0.000955 0.002776 1.3903 1.7615 0.00100 -0.00137 -0.00181
+ 1.85940 0.82158 -0.02381 -1.01955 0.001396 0.000996 0.002717 1.3912 1.7610 0.00103 -0.00142 -0.00186
+ 1.87620 0.83827 -0.02189 -1.01352 0.001456 0.001038 0.002656 1.3925 1.7602 0.00107 -0.00148 -0.00190
+ 1.89299 0.85494 -0.01994 -1.00707 0.001520 0.001082 0.002591 1.3944 1.7593 0.00110 -0.00153 -0.00194
+ 1.90976 0.87160 -0.01796 -1.00010 0.001589 0.001130 0.002522 1.3968 1.7580 0.00113 -0.00159 -0.00199
+ 1.92649 0.88821 -0.01594 -0.99247 0.001664 0.001180 0.002448 1.4002 1.7564 0.00116 -0.00165 -0.00203
+ 1.94315 0.90474 -0.01389 -0.98402 0.001748 0.001236 0.002366 1.4047 1.7543 0.00120 -0.00172 -0.00207
+ 1.95966 0.92112 -0.01182 -0.97452 0.001843 0.001298 0.002275 1.4107 1.7515 0.00123 -0.00180 -0.00210
+ 1.97587 0.93720 -0.00974 -0.96369 0.001953 0.001367 0.002171 1.4190 1.7477 0.00127 -0.00188 -0.00214
+ 1.99152 0.95272 -0.00771 -0.95121 0.002085 0.001448 0.002051 1.4306 1.7427 0.00131 -0.00198 -0.00217
+ 2.00620 0.96727 -0.00576 -0.93676 0.002246 0.001543 0.001912 1.4467 1.7357 0.00135 -0.00210 -0.00220
+ 2.01942 0.98037 -0.00398 -0.92027 0.002445 0.001655 0.001752 1.4691 1.7264 0.00140 -0.00225 -0.00223
+ 2.03084 0.99168 -0.00242 -0.90122 0.002704 0.001791 0.001566 1.5013 1.7136 0.00145 -0.00244 -0.00225
+ 2.03924 1.00000 -0.00126 -0.88651 0.002930 0.001906 0.001427 1.5294 1.7029 0.00150 -0.00260 -0.00226
+ 2.03924 1.00010 -0.00000 0.88651 0.010139 0.004923 0.000000 2.0504
+ 2.04764 1.00850 0.00000 0.86396 0.008930 0.005434 0.000000 1.6354
+ 2.05721 1.01807 0.00000 0.88159 0.007904 0.005054 0.000000 1.5556
+ 2.06811 1.02897 0.00000 0.89695 0.007149 0.004757 0.000000 1.4947
+ 2.08053 1.04139 0.00000 0.91011 0.006578 0.004523 0.000000 1.4460
+ 2.09468 1.05553 0.00000 0.92140 0.006130 0.004336 0.000000 1.4050
+ 2.11079 1.07165 0.00000 0.93111 0.005767 0.004185 0.000000 1.3692
+ 2.12915 1.09000 0.00000 0.93953 0.005465 0.004061 0.000000 1.3370
+ 2.15005 1.11091 0.00000 0.94689 0.005209 0.003958 0.000000 1.3073
+ 2.17387 1.13473 0.00000 0.95337 0.004986 0.003870 0.000000 1.2794
+ 2.20100 1.16186 0.00000 0.95911 0.004790 0.003795 0.000000 1.2532
+ 2.23190 1.19276 0.00000 0.96422 0.004616 0.003731 0.000000 1.2284
+ 2.26711 1.22797 0.00000 0.96878 0.004461 0.003675 0.000000 1.2049
+ 2.30721 1.26807 0.00000 0.97286 0.004322 0.003626 0.000000 1.1829
+ 2.35289 1.31375 0.00000 0.97650 0.004197 0.003583 0.000000 1.1623
+ 2.40493 1.36578 0.00000 0.97975 0.004086 0.003546 0.000000 1.1432
+ 2.46420 1.42506 0.00000 0.98263 0.003988 0.003514 0.000000 1.1258
+ 2.53172 1.49258 0.00000 0.98519 0.003901 0.003486 0.000000 1.1101
+ 2.60863 1.56949 0.00000 0.98743 0.003825 0.003461 0.000000 1.0961
+ 2.69625 1.65711 0.00000 0.98940 0.003759 0.003440 0.000000 1.0837
+ 2.79605 1.75691 0.00000 0.99111 0.003703 0.003422 0.000000 1.0730
+ 2.90974 1.87060 0.00000 0.99258 0.003654 0.003406 0.000000 1.0637
+ 3.03924 2.00010 0.00000 0.99397 0.003612 0.003392 0.000000 1.0557
diff --git a/dart/viscous/XfoilFiles/Case1FreeTrans.dat b/dart/viscous/XfoilFiles/Case1FreeTrans.dat
new file mode 100644
index 0000000..49e1c4f
--- /dev/null
+++ b/dart/viscous/XfoilFiles/Case1FreeTrans.dat
@@ -0,0 +1,186 @@
+0
+0.005
+# s x y Ue/Vinf Dstar Theta Cf H H* P m K tau Di
+ 0.00000 1.00000 0.00126 0.88216 0.002936 0.001903 0.001402 1.5348 1.7010 0.00148 0.00259 0.00222
+ 0.00840 0.99168 0.00242 0.89711 0.002704 0.001786 0.001543 1.5056 1.7120 0.00144 0.00243 0.00221
+ 0.01982 0.98037 0.00398 0.91655 0.002439 0.001647 0.001732 1.4722 1.7252 0.00138 0.00224 0.00219
+ 0.03304 0.96727 0.00576 0.93342 0.002235 0.001533 0.001895 1.4489 1.7348 0.00134 0.00209 0.00216
+ 0.04772 0.95272 0.00771 0.94818 0.002071 0.001437 0.002038 1.4321 1.7420 0.00129 0.00196 0.00213
+ 0.06337 0.93720 0.00974 0.96091 0.001938 0.001355 0.002160 1.4203 1.7473 0.00125 0.00186 0.00210
+ 0.07959 0.92112 0.01182 0.97193 0.001826 0.001285 0.002266 1.4117 1.7511 0.00121 0.00177 0.00207
+ 0.09610 0.90474 0.01389 0.98158 0.001731 0.001223 0.002359 1.4055 1.7540 0.00118 0.00170 0.00203
+ 0.11275 0.88821 0.01594 0.99016 0.001646 0.001167 0.002442 1.4009 1.7561 0.00114 0.00163 0.00199
+ 0.12949 0.87160 0.01796 0.99788 0.001571 0.001116 0.002517 1.3975 1.7578 0.00111 0.00157 0.00195
+ 0.14626 0.85494 0.01994 1.00494 0.001502 0.001069 0.002587 1.3950 1.7591 0.00108 0.00151 0.00191
+ 0.16305 0.83827 0.02189 1.01146 0.001437 0.001024 0.002652 1.3931 1.7600 0.00105 0.00145 0.00187
+ 0.17984 0.82158 0.02381 1.01754 0.001377 0.000982 0.002714 1.3918 1.7608 0.00102 0.00140 0.00182
+ 0.19665 0.80488 0.02569 1.02328 0.001319 0.000941 0.002774 1.3909 1.7613 0.00099 0.00135 0.00178
+ 0.21346 0.78817 0.02753 1.02872 0.001264 0.000903 0.002832 1.3904 1.7617 0.00096 0.00130 0.00173
+ 0.23027 0.77146 0.02934 1.03392 0.001212 0.000865 0.002889 1.3901 1.7620 0.00092 0.00125 0.00168
+ 0.24707 0.75475 0.03111 1.03893 0.001160 0.000828 0.002945 1.3902 1.7622 0.00089 0.00121 0.00164
+ 0.26388 0.73803 0.03284 1.04377 0.001110 0.000792 0.003001 1.3904 1.7623 0.00086 0.00116 0.00159
+ 0.28068 0.72132 0.03453 1.04848 0.001062 0.000757 0.003058 1.3908 1.7623 0.00083 0.00111 0.00154
+ 0.29748 0.70460 0.03619 1.05308 0.001014 0.000723 0.003115 1.3914 1.7623 0.00080 0.00107 0.00149
+ 0.31427 0.68789 0.03781 1.05759 0.000967 0.000689 0.003173 1.3922 1.7622 0.00077 0.00102 0.00144
+ 0.33105 0.67118 0.03938 1.06202 0.000921 0.000656 0.003232 1.3931 1.7620 0.00074 0.00098 0.00138
+ 0.34783 0.65447 0.04092 1.06639 0.000875 0.000623 0.003293 1.3941 1.7619 0.00071 0.00093 0.00133
+ 0.36460 0.63777 0.04241 1.07072 0.000830 0.000590 0.003357 1.3952 1.7617 0.00068 0.00089 0.00128
+ 0.38135 0.62108 0.04387 1.07500 0.000785 0.000557 0.003425 1.3963 1.7616 0.00064 0.00084 0.00122
+ 0.39809 0.60440 0.04527 1.07924 0.000740 0.000525 0.003497 1.3974 1.7615 0.00061 0.00080 0.00116
+ 0.41482 0.58772 0.04663 1.08346 0.000695 0.000493 0.003574 1.3984 1.7615 0.00058 0.00075 0.00110
+ 0.43154 0.57106 0.04794 1.08764 0.000651 0.000461 0.003659 1.3991 1.7617 0.00055 0.00071 0.00105
+ 0.44823 0.55441 0.04920 1.09178 0.000606 0.000429 0.003754 1.3993 1.7622 0.00051 0.00066 0.00098
+ 0.46491 0.53778 0.05041 1.09587 0.000561 0.000397 0.003862 1.3989 1.7630 0.00048 0.00061 0.00092
+ 0.48157 0.52116 0.05156 1.09989 0.000515 0.000366 0.003985 1.3977 1.7643 0.00044 0.00057 0.00086
+ 0.49820 0.50456 0.05265 1.10379 0.000470 0.000334 0.004127 1.3959 1.7660 0.00041 0.00052 0.00079
+ 0.51481 0.48798 0.05368 1.10748 0.000424 0.000301 0.004278 1.3950 1.7674 0.00037 0.00047 0.00072
+ 0.53140 0.47143 0.05465 1.11083 0.000380 0.000269 0.004396 1.4012 1.7661 0.00033 0.00042 0.00065
+ 0.54796 0.45489 0.05555 1.11372 0.000345 0.000238 0.004306 1.4374 1.7527 0.00029 0.00038 0.00058
+ 0.56448 0.43839 0.05639 1.11674 0.000334 0.000210 0.003525 1.5809 1.7023 0.00026 0.00037 0.00050
+ 0.58098 0.42191 0.05714 1.12383 0.000386 0.000188 0.001703 2.0373 1.5929 0.00024 0.00043 0.00043
+ 0.59744 0.40546 0.05782 1.13604 0.000484 0.000173 0.000350 2.7733 1.5604 0.00022 0.00055 0.00040
+ 0.61386 0.38905 0.05842 1.13786 0.000481 0.000170 0.000333 2.8151 1.5575 0.00022 0.00055 0.00039
+ 0.63024 0.37268 0.05893 1.14121 0.000469 0.000165 0.000343 2.8152 1.5575 0.00022 0.00053 0.00038
+ 0.64657 0.35635 0.05935 1.14504 0.000453 0.000161 0.000363 2.8022 1.5583 0.00021 0.00052 0.00038
+ 0.66286 0.34007 0.05967 1.14901 0.000436 0.000156 0.000386 2.7853 1.5595 0.00021 0.00050 0.00037
+ 0.67910 0.32383 0.05989 1.15303 0.000420 0.000151 0.000413 2.7671 1.5608 0.00020 0.00048 0.00036
+ 0.69527 0.30766 0.06000 1.15704 0.000403 0.000146 0.000442 2.7487 1.5622 0.00019 0.00047 0.00035
+ 0.71139 0.29154 0.06000 1.16100 0.000386 0.000141 0.000473 2.7305 1.5635 0.00019 0.00045 0.00034
+ 0.72743 0.27550 0.05988 1.16489 0.000370 0.000135 0.000507 2.7128 1.5649 0.00018 0.00043 0.00034
+ 0.74340 0.25953 0.05963 1.16867 0.000354 0.000130 0.000543 2.6956 1.5663 0.00018 0.00041 0.00033
+ 0.75928 0.24366 0.05924 1.17232 0.000338 0.000125 0.000582 2.6791 1.5677 0.00017 0.00040 0.00032
+ 0.77506 0.22788 0.05871 1.17581 0.000322 0.000120 0.000625 2.6631 1.5690 0.00017 0.00038 0.00031
+ 0.79074 0.21222 0.05803 1.17910 0.000306 0.000115 0.000671 2.6476 1.5704 0.00016 0.00036 0.00030
+ 0.80628 0.19670 0.05718 1.18215 0.000291 0.000110 0.000721 2.6327 1.5717 0.00015 0.00034 0.00028
+ 0.82167 0.18135 0.05616 1.18491 0.000275 0.000104 0.000777 2.6183 1.5730 0.00015 0.00033 0.00027
+ 0.83688 0.16619 0.05497 1.18731 0.000260 0.000099 0.000838 2.6043 1.5743 0.00014 0.00031 0.00026
+ 0.85186 0.15127 0.05358 1.18928 0.000244 0.000094 0.000906 2.5908 1.5756 0.00013 0.00029 0.00025
+ 0.86655 0.13666 0.05200 1.19071 0.000229 0.000088 0.000983 2.5775 1.5769 0.00013 0.00027 0.00023
+ 0.88087 0.12246 0.05022 1.19149 0.000214 0.000083 0.001069 2.5644 1.5781 0.00012 0.00025 0.00022
+ 0.89470 0.10877 0.04824 1.19145 0.000199 0.000077 0.001167 2.5515 1.5794 0.00011 0.00024 0.00021
+ 0.90789 0.09575 0.04609 1.19044 0.000184 0.000072 0.001278 2.5385 1.5808 0.00010 0.00022 0.00019
+ 0.92028 0.08358 0.04381 1.18824 0.000169 0.000067 0.001404 2.5255 1.5821 0.00009 0.00020 0.00018
+ 0.93169 0.07242 0.04142 1.18468 0.000156 0.000061 0.001547 2.5123 1.5835 0.00009 0.00018 0.00016
+ 0.94200 0.06239 0.03900 1.17957 0.000142 0.000057 0.001707 2.4990 1.5850 0.00008 0.00017 0.00015
+ 0.95118 0.05354 0.03659 1.17277 0.000130 0.000052 0.001883 2.4854 1.5865 0.00007 0.00015 0.00013
+ 0.95926 0.04581 0.03424 1.16414 0.000119 0.000048 0.002076 2.4715 1.5881 0.00006 0.00014 0.00012
+ 0.96635 0.03910 0.03196 1.15349 0.000109 0.000044 0.002286 2.4574 1.5897 0.00006 0.00013 0.00011
+ 0.97255 0.03329 0.02976 1.14059 0.000099 0.000040 0.002513 2.4429 1.5914 0.00005 0.00011 0.00010
+ 0.97802 0.02826 0.02764 1.12512 0.000090 0.000037 0.002758 2.4281 1.5932 0.00005 0.00010 0.00008
+ 0.98286 0.02387 0.02560 1.10661 0.000083 0.000034 0.003019 2.4129 1.5950 0.00004 0.00009 0.00007
+ 0.98718 0.02003 0.02361 1.08447 0.000075 0.000031 0.003298 2.3972 1.5970 0.00004 0.00008 0.00006
+ 0.99107 0.01666 0.02167 1.05791 0.000068 0.000029 0.003592 2.3810 1.5991 0.00003 0.00007 0.00005
+ 0.99460 0.01369 0.01977 1.02591 0.000062 0.000026 0.003895 2.3643 1.6013 0.00003 0.00006 0.00005
+ 0.99781 0.01108 0.01789 0.98716 0.000056 0.000024 0.004197 2.3472 1.6036 0.00002 0.00006 0.00004
+ 1.00078 0.00879 0.01602 0.94002 0.000051 0.000022 0.004482 2.3298 1.6060 0.00002 0.00005 0.00003
+ 1.00352 0.00678 0.01415 0.88252 0.000046 0.000020 0.004720 2.3123 1.6085 0.00002 0.00004 0.00002
+ 1.00607 0.00505 0.01227 0.81239 0.000042 0.000018 0.004866 2.2950 1.6110 0.00001 0.00003 0.00002
+ 1.00847 0.00358 0.01038 0.72732 0.000038 0.000017 0.004860 2.2781 1.6135 0.00001 0.00003 0.00001
+ 1.01072 0.00236 0.00848 0.62542 0.000035 0.000015 0.004618 2.2630 1.6158 0.00001 0.00002 0.00001
+ 1.01286 0.00140 0.00657 0.50609 0.000033 0.000014 0.004080 2.2491 1.6180 0.00000 0.00002 0.00000
+ 1.01490 0.00070 0.00466 0.37105 0.000031 0.000014 0.003197 2.2395 1.6195 0.00000 0.00001 0.00000
+ 1.01685 0.00024 0.00276 0.22470 0.000030 0.000013 0.002026 2.2318 1.6207 0.00000 0.00001 0.00000
+ 1.01872 0.00003 0.00091 0.07433 0.000029 0.000013 0.000682 2.2295 1.6211 0.00000 0.00000 0.00000
+ 1.02053 0.00003 -0.00091 -0.07433 0.000029 0.000013 0.000682 2.2295 1.6211 0.00000 -0.00000 -0.00000
+ 1.02240 0.00024 -0.00276 -0.22470 0.000030 0.000013 0.002026 2.2318 1.6207 0.00000 -0.00001 -0.00000
+ 1.02434 0.00070 -0.00466 -0.37105 0.000031 0.000014 0.003197 2.2395 1.6195 0.00000 -0.00001 -0.00000
+ 1.02638 0.00140 -0.00657 -0.50609 0.000033 0.000014 0.004080 2.2491 1.6180 0.00000 -0.00002 -0.00000
+ 1.02852 0.00236 -0.00848 -0.62542 0.000035 0.000015 0.004618 2.2630 1.6158 0.00001 -0.00002 -0.00001
+ 1.03078 0.00358 -0.01038 -0.72732 0.000038 0.000017 0.004860 2.2781 1.6135 0.00001 -0.00003 -0.00001
+ 1.03317 0.00505 -0.01227 -0.81239 0.000042 0.000018 0.004866 2.2950 1.6110 0.00001 -0.00003 -0.00002
+ 1.03573 0.00678 -0.01415 -0.88252 0.000046 0.000020 0.004720 2.3123 1.6085 0.00002 -0.00004 -0.00002
+ 1.03847 0.00879 -0.01602 -0.94002 0.000051 0.000022 0.004482 2.3298 1.6060 0.00002 -0.00005 -0.00003
+ 1.04143 0.01108 -0.01789 -0.98716 0.000056 0.000024 0.004197 2.3472 1.6036 0.00002 -0.00006 -0.00004
+ 1.04465 0.01369 -0.01977 -1.02591 0.000062 0.000026 0.003895 2.3643 1.6013 0.00003 -0.00006 -0.00005
+ 1.04817 0.01666 -0.02167 -1.05791 0.000068 0.000029 0.003592 2.3810 1.5991 0.00003 -0.00007 -0.00005
+ 1.05206 0.02003 -0.02361 -1.08447 0.000075 0.000031 0.003298 2.3972 1.5970 0.00004 -0.00008 -0.00006
+ 1.05638 0.02387 -0.02560 -1.10661 0.000083 0.000034 0.003019 2.4129 1.5950 0.00004 -0.00009 -0.00007
+ 1.06122 0.02826 -0.02764 -1.12512 0.000090 0.000037 0.002758 2.4281 1.5932 0.00005 -0.00010 -0.00008
+ 1.06669 0.03329 -0.02976 -1.14059 0.000099 0.000040 0.002513 2.4429 1.5914 0.00005 -0.00011 -0.00010
+ 1.07290 0.03910 -0.03196 -1.15349 0.000109 0.000044 0.002286 2.4574 1.5897 0.00006 -0.00013 -0.00011
+ 1.07998 0.04581 -0.03424 -1.16414 0.000119 0.000048 0.002076 2.4715 1.5881 0.00006 -0.00014 -0.00012
+ 1.08806 0.05354 -0.03659 -1.17277 0.000130 0.000052 0.001883 2.4854 1.5865 0.00007 -0.00015 -0.00013
+ 1.09724 0.06239 -0.03900 -1.17957 0.000142 0.000057 0.001707 2.4990 1.5850 0.00008 -0.00017 -0.00015
+ 1.10756 0.07242 -0.04142 -1.18468 0.000156 0.000061 0.001547 2.5123 1.5835 0.00009 -0.00018 -0.00016
+ 1.11897 0.08358 -0.04381 -1.18824 0.000169 0.000067 0.001404 2.5255 1.5821 0.00009 -0.00020 -0.00018
+ 1.13135 0.09575 -0.04609 -1.19044 0.000184 0.000072 0.001278 2.5385 1.5808 0.00010 -0.00022 -0.00019
+ 1.14454 0.10877 -0.04824 -1.19145 0.000199 0.000077 0.001167 2.5515 1.5794 0.00011 -0.00024 -0.00021
+ 1.15837 0.12246 -0.05022 -1.19149 0.000214 0.000083 0.001069 2.5644 1.5781 0.00012 -0.00025 -0.00022
+ 1.17269 0.13666 -0.05200 -1.19071 0.000229 0.000088 0.000983 2.5775 1.5769 0.00013 -0.00027 -0.00023
+ 1.18738 0.15127 -0.05358 -1.18928 0.000244 0.000094 0.000906 2.5908 1.5756 0.00013 -0.00029 -0.00025
+ 1.20236 0.16619 -0.05497 -1.18731 0.000260 0.000099 0.000838 2.6043 1.5743 0.00014 -0.00031 -0.00026
+ 1.21757 0.18135 -0.05616 -1.18491 0.000275 0.000104 0.000777 2.6183 1.5730 0.00015 -0.00033 -0.00027
+ 1.23296 0.19670 -0.05718 -1.18215 0.000291 0.000110 0.000721 2.6327 1.5717 0.00015 -0.00034 -0.00028
+ 1.24851 0.21222 -0.05803 -1.17910 0.000306 0.000115 0.000671 2.6476 1.5704 0.00016 -0.00036 -0.00030
+ 1.26418 0.22788 -0.05871 -1.17581 0.000322 0.000120 0.000625 2.6631 1.5690 0.00017 -0.00038 -0.00031
+ 1.27996 0.24366 -0.05924 -1.17232 0.000338 0.000125 0.000582 2.6791 1.5677 0.00017 -0.00040 -0.00032
+ 1.29584 0.25953 -0.05963 -1.16867 0.000354 0.000130 0.000543 2.6956 1.5663 0.00018 -0.00041 -0.00033
+ 1.31181 0.27550 -0.05988 -1.16489 0.000370 0.000135 0.000507 2.7128 1.5649 0.00018 -0.00043 -0.00034
+ 1.32785 0.29154 -0.06000 -1.16100 0.000386 0.000141 0.000473 2.7305 1.5635 0.00019 -0.00045 -0.00034
+ 1.34397 0.30766 -0.06000 -1.15704 0.000403 0.000146 0.000442 2.7487 1.5622 0.00019 -0.00047 -0.00035
+ 1.36015 0.32383 -0.05989 -1.15303 0.000420 0.000151 0.000413 2.7671 1.5608 0.00020 -0.00048 -0.00036
+ 1.37638 0.34007 -0.05967 -1.14901 0.000436 0.000156 0.000386 2.7853 1.5595 0.00021 -0.00050 -0.00037
+ 1.39267 0.35635 -0.05935 -1.14504 0.000453 0.000161 0.000363 2.8022 1.5583 0.00021 -0.00052 -0.00038
+ 1.40900 0.37268 -0.05893 -1.14121 0.000469 0.000165 0.000343 2.8152 1.5575 0.00022 -0.00053 -0.00038
+ 1.42538 0.38905 -0.05842 -1.13786 0.000481 0.000170 0.000333 2.8151 1.5575 0.00022 -0.00055 -0.00039
+ 1.44181 0.40546 -0.05782 -1.13604 0.000484 0.000173 0.000350 2.7733 1.5604 0.00022 -0.00055 -0.00040
+ 1.45826 0.42191 -0.05714 -1.12383 0.000386 0.000188 0.001703 2.0373 1.5929 0.00024 -0.00043 -0.00043
+ 1.47476 0.43839 -0.05639 -1.11674 0.000334 0.000210 0.003525 1.5809 1.7023 0.00026 -0.00037 -0.00050
+ 1.49129 0.45489 -0.05555 -1.11372 0.000345 0.000238 0.004306 1.4374 1.7527 0.00029 -0.00038 -0.00058
+ 1.50784 0.47143 -0.05465 -1.11083 0.000380 0.000269 0.004396 1.4012 1.7661 0.00033 -0.00042 -0.00065
+ 1.52443 0.48798 -0.05368 -1.10748 0.000424 0.000301 0.004278 1.3950 1.7674 0.00037 -0.00047 -0.00072
+ 1.54104 0.50456 -0.05265 -1.10379 0.000470 0.000334 0.004127 1.3959 1.7660 0.00041 -0.00052 -0.00079
+ 1.55768 0.52116 -0.05156 -1.09989 0.000515 0.000366 0.003985 1.3977 1.7643 0.00044 -0.00057 -0.00086
+ 1.57433 0.53778 -0.05041 -1.09587 0.000561 0.000397 0.003862 1.3989 1.7630 0.00048 -0.00061 -0.00092
+ 1.59101 0.55441 -0.04920 -1.09178 0.000606 0.000429 0.003754 1.3993 1.7622 0.00051 -0.00066 -0.00098
+ 1.60771 0.57106 -0.04794 -1.08764 0.000651 0.000461 0.003659 1.3991 1.7617 0.00055 -0.00071 -0.00105
+ 1.62442 0.58772 -0.04663 -1.08346 0.000695 0.000493 0.003574 1.3984 1.7615 0.00058 -0.00075 -0.00110
+ 1.64115 0.60440 -0.04527 -1.07924 0.000740 0.000525 0.003497 1.3974 1.7615 0.00061 -0.00080 -0.00116
+ 1.65789 0.62108 -0.04387 -1.07500 0.000785 0.000557 0.003425 1.3963 1.7616 0.00064 -0.00084 -0.00122
+ 1.67465 0.63777 -0.04241 -1.07072 0.000830 0.000590 0.003357 1.3952 1.7617 0.00068 -0.00089 -0.00128
+ 1.69141 0.65447 -0.04092 -1.06639 0.000875 0.000623 0.003293 1.3941 1.7619 0.00071 -0.00093 -0.00133
+ 1.70819 0.67118 -0.03938 -1.06202 0.000921 0.000656 0.003232 1.3931 1.7620 0.00074 -0.00098 -0.00138
+ 1.72497 0.68789 -0.03781 -1.05759 0.000967 0.000689 0.003173 1.3922 1.7622 0.00077 -0.00102 -0.00144
+ 1.74177 0.70460 -0.03619 -1.05308 0.001014 0.000723 0.003115 1.3914 1.7623 0.00080 -0.00107 -0.00149
+ 1.75856 0.72132 -0.03453 -1.04848 0.001062 0.000757 0.003058 1.3908 1.7623 0.00083 -0.00111 -0.00154
+ 1.77536 0.73803 -0.03284 -1.04377 0.001110 0.000792 0.003001 1.3904 1.7623 0.00086 -0.00116 -0.00159
+ 1.79217 0.75475 -0.03111 -1.03893 0.001160 0.000828 0.002945 1.3902 1.7622 0.00089 -0.00121 -0.00164
+ 1.80898 0.77146 -0.02934 -1.03392 0.001212 0.000865 0.002889 1.3901 1.7620 0.00092 -0.00125 -0.00168
+ 1.82579 0.78817 -0.02753 -1.02872 0.001264 0.000903 0.002832 1.3904 1.7617 0.00096 -0.00130 -0.00173
+ 1.84259 0.80488 -0.02569 -1.02328 0.001319 0.000941 0.002774 1.3909 1.7613 0.00099 -0.00135 -0.00178
+ 1.85940 0.82158 -0.02381 -1.01754 0.001377 0.000982 0.002714 1.3918 1.7608 0.00102 -0.00140 -0.00182
+ 1.87620 0.83827 -0.02189 -1.01146 0.001437 0.001024 0.002652 1.3931 1.7600 0.00105 -0.00145 -0.00187
+ 1.89299 0.85494 -0.01994 -1.00494 0.001502 0.001069 0.002587 1.3950 1.7591 0.00108 -0.00151 -0.00191
+ 1.90976 0.87160 -0.01796 -0.99788 0.001571 0.001116 0.002517 1.3975 1.7578 0.00111 -0.00157 -0.00195
+ 1.92649 0.88821 -0.01594 -0.99016 0.001646 0.001167 0.002442 1.4009 1.7561 0.00114 -0.00163 -0.00199
+ 1.94315 0.90474 -0.01389 -0.98158 0.001731 0.001223 0.002359 1.4055 1.7540 0.00118 -0.00170 -0.00203
+ 1.95966 0.92112 -0.01182 -0.97193 0.001826 0.001285 0.002266 1.4117 1.7511 0.00121 -0.00177 -0.00207
+ 1.97587 0.93720 -0.00974 -0.96091 0.001938 0.001355 0.002160 1.4203 1.7473 0.00125 -0.00186 -0.00210
+ 1.99152 0.95272 -0.00771 -0.94818 0.002071 0.001437 0.002038 1.4321 1.7420 0.00129 -0.00196 -0.00213
+ 2.00620 0.96727 -0.00576 -0.93342 0.002235 0.001533 0.001895 1.4489 1.7348 0.00134 -0.00209 -0.00216
+ 2.01942 0.98037 -0.00398 -0.91655 0.002439 0.001647 0.001732 1.4722 1.7252 0.00138 -0.00224 -0.00219
+ 2.03084 0.99168 -0.00242 -0.89711 0.002704 0.001786 0.001543 1.5056 1.7120 0.00144 -0.00243 -0.00221
+ 2.03924 1.00000 -0.00126 -0.88216 0.002936 0.001903 0.001402 1.5348 1.7010 0.00148 -0.00259 -0.00222
+ 2.03924 1.00010 -0.00000 0.88216 0.008367 0.003806 0.000000 2.1892
+ 2.04764 1.00850 0.00000 0.85638 0.007078 0.004274 0.000000 1.6483
+ 2.05721 1.01807 0.00000 0.87698 0.006118 0.003924 0.000000 1.5511
+ 2.06811 1.02897 0.00000 0.89434 0.005455 0.003664 0.000000 1.4809
+ 2.08053 1.04139 0.00000 0.90867 0.004977 0.003468 0.000000 1.4266
+ 2.09468 1.05553 0.00000 0.92060 0.004613 0.003318 0.000000 1.3818
+ 2.11079 1.07165 0.00000 0.93065 0.004325 0.003200 0.000000 1.3431
+ 2.12915 1.09000 0.00000 0.93924 0.004088 0.003104 0.000000 1.3085
+ 2.15005 1.11091 0.00000 0.94669 0.003888 0.003024 0.000000 1.2770
+ 2.17387 1.13473 0.00000 0.95321 0.003716 0.002957 0.000000 1.2478
+ 2.20100 1.16186 0.00000 0.95899 0.003566 0.002900 0.000000 1.2207
+ 2.23190 1.19276 0.00000 0.96412 0.003434 0.002851 0.000000 1.1956
+ 2.26711 1.22797 0.00000 0.96870 0.003318 0.002809 0.000000 1.1724
+ 2.30721 1.26807 0.00000 0.97279 0.003216 0.002772 0.000000 1.1512
+ 2.35289 1.31375 0.00000 0.97645 0.003126 0.002739 0.000000 1.1320
+ 2.40493 1.36578 0.00000 0.97970 0.003047 0.002711 0.000000 1.1147
+ 2.46420 1.42506 0.00000 0.98260 0.002978 0.002687 0.000000 1.0995
+ 2.53172 1.49258 0.00000 0.98516 0.002919 0.002665 0.000000 1.0862
+ 2.60863 1.56949 0.00000 0.98741 0.002868 0.002647 0.000000 1.0747
+ 2.69625 1.65711 0.00000 0.98938 0.002825 0.002630 0.000000 1.0648
+ 2.79605 1.75691 0.00000 0.99109 0.002788 0.002617 0.000000 1.0565
+ 2.90974 1.87060 0.00000 0.99257 0.002757 0.002605 0.000000 1.0494
+ 3.03924 2.00010 0.00000 0.99392 0.002731 0.002594 0.000000 1.0435
diff --git a/dart/viscous/XfoilFiles/case12Xfoil.dat b/dart/viscous/XfoilFiles/case12Xfoil.dat
new file mode 100644
index 0000000..9fc97d5
--- /dev/null
+++ b/dart/viscous/XfoilFiles/case12Xfoil.dat
@@ -0,0 +1,186 @@
+1.3329
+0.01068
+# s x y Ue/Vinf Dstar Theta Cf H H* P m K tau Di
+ 0.00000 1.00000 0.00126 0.90089 0.014842 0.007386 0.000403 2.0096 1.5738 0.00599 0.01337 0.00850
+ 0.00840 0.99168 0.00242 0.90875 0.013953 0.007132 0.000455 1.9562 1.5837 0.00589 0.01268 0.00848
+ 0.01982 0.98037 0.00398 0.92007 0.012836 0.006791 0.000532 1.8903 1.5972 0.00575 0.01181 0.00845
+ 0.03304 0.96727 0.00576 0.93342 0.011722 0.006419 0.000621 1.8260 1.6117 0.00559 0.01094 0.00841
+ 0.04772 0.95272 0.00771 0.94823 0.010677 0.006042 0.000721 1.7673 1.6263 0.00543 0.01012 0.00838
+ 0.06337 0.93720 0.00974 0.96373 0.009748 0.005680 0.000826 1.7163 1.6400 0.00528 0.00939 0.00834
+ 0.07959 0.92112 0.01182 0.97930 0.008946 0.005346 0.000931 1.6735 1.6524 0.00513 0.00876 0.00830
+ 0.09610 0.90474 0.01389 0.99453 0.008262 0.005044 0.001033 1.6379 1.6634 0.00499 0.00822 0.00825
+ 0.11275 0.88821 0.01594 1.00923 0.007680 0.004775 0.001132 1.6085 1.6728 0.00486 0.00775 0.00821
+ 0.12949 0.87160 0.01796 1.02333 0.007180 0.004533 0.001226 1.5841 1.6811 0.00475 0.00735 0.00817
+ 0.14626 0.85494 0.01994 1.03685 0.006748 0.004315 0.001316 1.5636 1.6882 0.00464 0.00700 0.00812
+ 0.16305 0.83827 0.02189 1.04983 0.006368 0.004118 0.001403 1.5464 1.6943 0.00454 0.00669 0.00807
+ 0.17984 0.82158 0.02381 1.06233 0.006033 0.003939 0.001485 1.5316 1.6997 0.00444 0.00641 0.00803
+ 0.19665 0.80488 0.02569 1.07443 0.005732 0.003774 0.001565 1.5190 1.7044 0.00436 0.00616 0.00798
+ 0.21346 0.78817 0.02753 1.08618 0.005460 0.003621 0.001641 1.5081 1.7086 0.00427 0.00593 0.00793
+ 0.23027 0.77146 0.02934 1.09764 0.005212 0.003478 0.001716 1.4985 1.7123 0.00419 0.00572 0.00788
+ 0.24707 0.75475 0.03111 1.10887 0.004984 0.003345 0.001789 1.4901 1.7156 0.00411 0.00553 0.00782
+ 0.26388 0.73803 0.03284 1.11992 0.004773 0.003219 0.001861 1.4826 1.7186 0.00404 0.00535 0.00777
+ 0.28068 0.72132 0.03453 1.13084 0.004576 0.003100 0.001931 1.4760 1.7212 0.00396 0.00517 0.00772
+ 0.29748 0.70460 0.03619 1.14166 0.004390 0.002987 0.002001 1.4701 1.7236 0.00389 0.00501 0.00766
+ 0.31427 0.68789 0.03781 1.15242 0.004216 0.002878 0.002071 1.4647 1.7258 0.00382 0.00486 0.00760
+ 0.33105 0.67118 0.03938 1.16316 0.004050 0.002774 0.002141 1.4599 1.7278 0.00375 0.00471 0.00754
+ 0.34783 0.65447 0.04092 1.17390 0.003892 0.002674 0.002210 1.4555 1.7296 0.00369 0.00457 0.00748
+ 0.36460 0.63777 0.04241 1.18469 0.003742 0.002578 0.002281 1.4515 1.7313 0.00362 0.00443 0.00742
+ 0.38135 0.62108 0.04387 1.19554 0.003597 0.002485 0.002352 1.4478 1.7329 0.00355 0.00430 0.00736
+ 0.39809 0.60440 0.04527 1.20647 0.003458 0.002394 0.002424 1.4444 1.7343 0.00348 0.00417 0.00729
+ 0.41482 0.58772 0.04663 1.21753 0.003324 0.002306 0.002498 1.4413 1.7357 0.00342 0.00405 0.00722
+ 0.43154 0.57106 0.04794 1.22872 0.003195 0.002221 0.002573 1.4384 1.7369 0.00335 0.00393 0.00716
+ 0.44823 0.55441 0.04920 1.24007 0.003069 0.002138 0.002650 1.4357 1.7381 0.00329 0.00381 0.00709
+ 0.46491 0.53778 0.05041 1.25161 0.002948 0.002057 0.002729 1.4333 1.7392 0.00322 0.00369 0.00701
+ 0.48157 0.52116 0.05156 1.26336 0.002829 0.001977 0.002810 1.4310 1.7402 0.00316 0.00357 0.00694
+ 0.49820 0.50456 0.05265 1.27534 0.002715 0.001900 0.002893 1.4288 1.7411 0.00309 0.00346 0.00686
+ 0.51481 0.48798 0.05368 1.28758 0.002603 0.001824 0.002980 1.4268 1.7420 0.00302 0.00335 0.00678
+ 0.53140 0.47143 0.05465 1.30010 0.002494 0.001750 0.003069 1.4250 1.7428 0.00296 0.00324 0.00670
+ 0.54796 0.45489 0.05555 1.31293 0.002387 0.001677 0.003162 1.4233 1.7436 0.00289 0.00313 0.00662
+ 0.56448 0.43839 0.05639 1.32610 0.002283 0.001606 0.003258 1.4217 1.7443 0.00282 0.00303 0.00653
+ 0.58098 0.42191 0.05714 1.33963 0.002181 0.001536 0.003359 1.4203 1.7450 0.00276 0.00292 0.00644
+ 0.59744 0.40546 0.05782 1.35356 0.002082 0.001467 0.003464 1.4190 1.7456 0.00269 0.00282 0.00635
+ 0.61386 0.38905 0.05842 1.36793 0.001985 0.001400 0.003574 1.4178 1.7462 0.00262 0.00272 0.00626
+ 0.63024 0.37268 0.05893 1.38277 0.001890 0.001334 0.003690 1.4167 1.7467 0.00255 0.00261 0.00616
+ 0.64657 0.35635 0.05935 1.39813 0.001797 0.001269 0.003811 1.4158 1.7472 0.00248 0.00251 0.00606
+ 0.66286 0.34007 0.05967 1.41407 0.001705 0.001205 0.003939 1.4150 1.7476 0.00241 0.00241 0.00595
+ 0.67910 0.32383 0.05989 1.43062 0.001616 0.001142 0.004074 1.4144 1.7479 0.00234 0.00231 0.00585
+ 0.69527 0.30766 0.06000 1.44787 0.001528 0.001081 0.004218 1.4139 1.7482 0.00227 0.00221 0.00573
+ 0.71139 0.29154 0.06000 1.46588 0.001442 0.001020 0.004371 1.4136 1.7484 0.00219 0.00211 0.00562
+ 0.72743 0.27550 0.05988 1.48475 0.001358 0.000960 0.004534 1.4134 1.7486 0.00212 0.00202 0.00550
+ 0.74340 0.25953 0.05963 1.50457 0.001275 0.000902 0.004709 1.4134 1.7487 0.00204 0.00192 0.00537
+ 0.75928 0.24366 0.05924 1.52546 0.001193 0.000844 0.004897 1.4137 1.7487 0.00196 0.00182 0.00524
+ 0.77506 0.22788 0.05871 1.54758 0.001113 0.000787 0.005101 1.4141 1.7486 0.00189 0.00172 0.00510
+ 0.79074 0.21222 0.05803 1.57108 0.001035 0.000732 0.005323 1.4148 1.7484 0.00181 0.00163 0.00496
+ 0.80628 0.19670 0.05718 1.59617 0.000958 0.000677 0.005566 1.4157 1.7482 0.00172 0.00153 0.00481
+ 0.82167 0.18135 0.05616 1.62310 0.000882 0.000623 0.005834 1.4169 1.7478 0.00164 0.00143 0.00465
+ 0.83688 0.16619 0.05497 1.65216 0.000808 0.000569 0.006131 1.4184 1.7474 0.00155 0.00133 0.00449
+ 0.85186 0.15127 0.05358 1.68368 0.000735 0.000517 0.006465 1.4203 1.7468 0.00147 0.00124 0.00431
+ 0.86655 0.13666 0.05200 1.71806 0.000664 0.000467 0.006842 1.4225 1.7461 0.00138 0.00114 0.00413
+ 0.88087 0.12246 0.05022 1.75570 0.000594 0.000417 0.007272 1.4250 1.7453 0.00129 0.00104 0.00394
+ 0.89470 0.10877 0.04824 1.79701 0.000527 0.000369 0.007764 1.4280 1.7444 0.00119 0.00095 0.00374
+ 0.90789 0.09575 0.04609 1.84231 0.000463 0.000324 0.008332 1.4314 1.7435 0.00110 0.00085 0.00353
+ 0.92028 0.08358 0.04381 1.89171 0.000403 0.000281 0.008987 1.4351 1.7425 0.00101 0.00076 0.00331
+ 0.93169 0.07242 0.04142 1.94509 0.000347 0.000241 0.009740 1.4391 1.7415 0.00091 0.00068 0.00309
+ 0.94200 0.06239 0.03900 2.00205 0.000297 0.000206 0.010603 1.4431 1.7406 0.00083 0.00059 0.00288
+ 0.95118 0.05354 0.03659 2.06206 0.000252 0.000174 0.011589 1.4468 1.7400 0.00074 0.00052 0.00266
+ 0.95926 0.04581 0.03424 2.12458 0.000212 0.000147 0.012719 1.4500 1.7398 0.00066 0.00045 0.00245
+ 0.96635 0.03910 0.03196 2.18920 0.000178 0.000123 0.014028 1.4517 1.7405 0.00059 0.00039 0.00224
+ 0.97255 0.03329 0.02976 2.25567 0.000148 0.000102 0.015581 1.4510 1.7424 0.00052 0.00033 0.00203
+ 0.97802 0.02826 0.02764 2.32377 0.000121 0.000084 0.017479 1.4463 1.7462 0.00045 0.00028 0.00183
+ 0.98286 0.02387 0.02560 2.39319 0.000097 0.000068 0.019854 1.4369 1.7525 0.00039 0.00023 0.00163
+ 0.98718 0.02003 0.02361 2.46325 0.000077 0.000054 0.022688 1.4273 1.7594 0.00033 0.00019 0.00142
+ 0.99107 0.01666 0.02167 2.53292 0.000061 0.000042 0.024423 1.4574 1.7516 0.00027 0.00016 0.00120
+ 0.99460 0.01369 0.01977 2.60367 0.000056 0.000033 0.018286 1.7067 1.6761 0.00022 0.00015 0.00097
+ 0.99781 0.01108 0.01789 2.70731 0.000078 0.000026 0.002218 3.0028 1.5159 0.00019 0.00021 0.00078
+ 1.00078 0.00879 0.01602 2.80047 0.000086 0.000021 -0.000333 4.0432 1.5281 0.00017 0.00024 0.00072
+ 1.00352 0.00678 0.01415 2.85102 0.000066 0.000019 0.001214 3.4042 1.5335 0.00016 0.00019 0.00068
+ 1.00607 0.00505 0.01227 2.90955 0.000051 0.000017 0.003812 2.9799 1.5479 0.00014 0.00015 0.00065
+ 1.00847 0.00358 0.01038 2.96013 0.000041 0.000015 0.006780 2.7348 1.5632 0.00013 0.00012 0.00060
+ 1.01072 0.00236 0.00848 2.99057 0.000034 0.000013 0.009966 2.5743 1.5772 0.00012 0.00010 0.00056
+ 1.01286 0.00140 0.00657 2.99110 0.000029 0.000012 0.013204 2.4597 1.5894 0.00011 0.00009 0.00050
+ 1.01490 0.00070 0.00466 2.95345 0.000025 0.000011 0.016155 2.3755 1.5998 0.00009 0.00008 0.00044
+ 1.01685 0.00024 0.00276 2.87337 0.000023 0.000010 0.018443 2.3128 1.6084 0.00008 0.00007 0.00038
+ 1.01872 0.00003 0.00091 2.75432 0.000022 0.000010 0.019442 2.2784 1.6135 0.00007 0.00006 0.00032
+ 1.02053 0.00003 -0.00091 2.60885 0.000021 0.000009 0.019599 2.2556 1.6169 0.00006 0.00005 0.00026
+ 1.02240 0.00024 -0.00276 2.43348 0.000020 0.000009 0.019300 2.2288 1.6212 0.00005 0.00005 0.00021
+ 1.02434 0.00070 -0.00466 2.22710 0.000020 0.000009 0.018107 2.2080 1.6246 0.00005 0.00004 0.00016
+ 1.02638 0.00140 -0.00657 2.00097 0.000020 0.000009 0.016184 2.1935 1.6270 0.00004 0.00004 0.00012
+ 1.02852 0.00236 -0.00848 1.76780 0.000021 0.000010 0.013859 2.1844 1.6285 0.00003 0.00004 0.00009
+ 1.03078 0.00358 -0.01038 1.53897 0.000022 0.000010 0.011450 2.1802 1.6293 0.00002 0.00003 0.00006
+ 1.03317 0.00505 -0.01227 1.32213 0.000024 0.000011 0.009211 2.1786 1.6295 0.00002 0.00003 0.00004
+ 1.03573 0.00678 -0.01415 1.12131 0.000026 0.000012 0.007225 2.1810 1.6291 0.00001 0.00003 0.00003
+ 1.03847 0.00879 -0.01602 0.93765 0.000028 0.000013 0.005581 2.1819 1.6290 0.00001 0.00003 0.00002
+ 1.04143 0.01108 -0.01789 0.77078 0.000031 0.000014 0.004182 2.1903 1.6275 0.00001 0.00002 0.00001
+ 1.04465 0.01369 -0.01977 0.61918 0.000033 0.000015 0.003105 2.1887 1.6278 0.00001 0.00002 0.00001
+ 1.04817 0.01666 -0.02167 0.48137 0.000036 0.000016 0.002169 2.2058 1.6249 0.00000 0.00002 0.00000
+ 1.05206 0.02003 -0.02361 0.35535 0.000039 0.000018 0.001498 2.1957 1.6266 0.00000 0.00001 0.00000
+ 1.05638 0.02387 -0.02560 0.23983 0.000043 0.000019 0.000890 2.2256 1.6217 0.00000 0.00001 0.00000
+ 1.06122 0.02826 -0.02764 0.13308 0.000046 0.000021 0.000472 2.2126 1.6238 0.00000 0.00001 0.00000
+ 1.06669 0.03329 -0.02976 0.03425 0.000053 0.000024 0.000103 2.2295 1.6211 0.00000 0.00000 0.00000
+ 1.07290 0.03910 -0.03196 -0.05772 0.000053 0.000024 0.000173 2.2295 1.6211 0.00000 -0.00000 -0.00000
+ 1.07998 0.04581 -0.03424 -0.14345 0.000062 0.000027 0.000361 2.2596 1.6163 0.00000 -0.00001 -0.00000
+ 1.08806 0.05354 -0.03659 -0.22296 0.000066 0.000029 0.000541 2.2321 1.6206 0.00000 -0.00001 -0.00000
+ 1.09724 0.06239 -0.03900 -0.29650 0.000073 0.000032 0.000622 2.2748 1.6140 0.00000 -0.00002 -0.00000
+ 1.10756 0.07242 -0.04142 -0.36357 0.000079 0.000035 0.000720 2.2571 1.6167 0.00000 -0.00003 -0.00000
+ 1.11897 0.08358 -0.04381 -0.42417 0.000087 0.000038 0.000744 2.2842 1.6126 0.00001 -0.00004 -0.00000
+ 1.13135 0.09575 -0.04609 -0.47802 0.000094 0.000041 0.000780 2.2805 1.6131 0.00001 -0.00004 -0.00001
+ 1.14454 0.10877 -0.04824 -0.52552 0.000101 0.000044 0.000780 2.2951 1.6110 0.00001 -0.00005 -0.00001
+ 1.15837 0.12246 -0.05022 -0.56708 0.000109 0.000047 0.000782 2.2994 1.6104 0.00002 -0.00006 -0.00001
+ 1.17269 0.13666 -0.05200 -0.60340 0.000117 0.000051 0.000770 2.3079 1.6091 0.00002 -0.00007 -0.00002
+ 1.18738 0.15127 -0.05358 -0.63516 0.000124 0.000054 0.000757 2.3142 1.6082 0.00002 -0.00008 -0.00002
+ 1.20236 0.16619 -0.05497 -0.66302 0.000132 0.000057 0.000740 2.3209 1.6073 0.00002 -0.00009 -0.00003
+ 1.21757 0.18135 -0.05616 -0.68756 0.000139 0.000060 0.000722 2.3271 1.6064 0.00003 -0.00010 -0.00003
+ 1.23296 0.19670 -0.05718 -0.70927 0.000146 0.000063 0.000704 2.3331 1.6055 0.00003 -0.00010 -0.00004
+ 1.24851 0.21222 -0.05803 -0.72857 0.000154 0.000066 0.000685 2.3389 1.6047 0.00003 -0.00011 -0.00004
+ 1.26418 0.22788 -0.05871 -0.74580 0.000161 0.000069 0.000666 2.3445 1.6040 0.00004 -0.00012 -0.00005
+ 1.27996 0.24366 -0.05924 -0.76124 0.000168 0.000072 0.000647 2.3500 1.6032 0.00004 -0.00013 -0.00005
+ 1.29584 0.25953 -0.05963 -0.77513 0.000175 0.000074 0.000629 2.3553 1.6025 0.00004 -0.00014 -0.00006
+ 1.31181 0.27550 -0.05988 -0.78767 0.000182 0.000077 0.000612 2.3605 1.6018 0.00005 -0.00014 -0.00006
+ 1.32785 0.29154 -0.06000 -0.79902 0.000189 0.000080 0.000595 2.3655 1.6011 0.00005 -0.00015 -0.00007
+ 1.34397 0.30766 -0.06000 -0.80934 0.000196 0.000083 0.000578 2.3705 1.6005 0.00005 -0.00016 -0.00007
+ 1.36015 0.32383 -0.05989 -0.81873 0.000203 0.000085 0.000563 2.3753 1.5998 0.00006 -0.00017 -0.00008
+ 1.37638 0.34007 -0.05967 -0.82731 0.000210 0.000088 0.000548 2.3799 1.5992 0.00006 -0.00017 -0.00008
+ 1.39267 0.35635 -0.05935 -0.83516 0.000217 0.000091 0.000533 2.3845 1.5986 0.00006 -0.00018 -0.00008
+ 1.40900 0.37268 -0.05893 -0.84236 0.000223 0.000093 0.000519 2.3889 1.5981 0.00007 -0.00019 -0.00009
+ 1.42538 0.38905 -0.05842 -0.84898 0.000230 0.000096 0.000506 2.3932 1.5975 0.00007 -0.00020 -0.00009
+ 1.44181 0.40546 -0.05782 -0.85508 0.000236 0.000099 0.000493 2.3974 1.5970 0.00007 -0.00020 -0.00010
+ 1.45826 0.42191 -0.05714 -0.86072 0.000243 0.000101 0.000481 2.4015 1.5965 0.00007 -0.00021 -0.00010
+ 1.47476 0.43839 -0.05639 -0.86593 0.000250 0.000104 0.000470 2.4054 1.5960 0.00008 -0.00022 -0.00011
+ 1.49129 0.45489 -0.05555 -0.87076 0.000256 0.000106 0.000459 2.4092 1.5955 0.00008 -0.00022 -0.00011
+ 1.50784 0.47143 -0.05465 -0.87524 0.000262 0.000109 0.000448 2.4129 1.5950 0.00008 -0.00023 -0.00012
+ 1.52443 0.48798 -0.05368 -0.87942 0.000269 0.000111 0.000438 2.4164 1.5946 0.00009 -0.00024 -0.00012
+ 1.54104 0.50456 -0.05265 -0.88331 0.000275 0.000114 0.000428 2.4198 1.5942 0.00009 -0.00024 -0.00012
+ 1.55768 0.52116 -0.05156 -0.88694 0.000281 0.000116 0.000419 2.4231 1.5938 0.00009 -0.00025 -0.00013
+ 1.57433 0.53778 -0.05041 -0.89034 0.000287 0.000118 0.000411 2.4263 1.5934 0.00009 -0.00026 -0.00013
+ 1.59101 0.55441 -0.04920 -0.89352 0.000294 0.000121 0.000402 2.4294 1.5930 0.00010 -0.00026 -0.00014
+ 1.60771 0.57106 -0.04794 -0.89651 0.000300 0.000123 0.000394 2.4324 1.5927 0.00010 -0.00027 -0.00014
+ 1.62442 0.58772 -0.04663 -0.89931 0.000306 0.000125 0.000386 2.4353 1.5923 0.00010 -0.00027 -0.00015
+ 1.64115 0.60440 -0.04527 -0.90194 0.000312 0.000128 0.000379 2.4381 1.5920 0.00010 -0.00028 -0.00015
+ 1.65789 0.62108 -0.04387 -0.90442 0.000318 0.000130 0.000372 2.4409 1.5916 0.00011 -0.00029 -0.00015
+ 1.67465 0.63777 -0.04241 -0.90674 0.000323 0.000132 0.000365 2.4437 1.5913 0.00011 -0.00029 -0.00016
+ 1.69141 0.65447 -0.04092 -0.90893 0.000329 0.000135 0.000358 2.4466 1.5910 0.00011 -0.00030 -0.00016
+ 1.70819 0.67118 -0.03938 -0.91097 0.000335 0.000137 0.000352 2.4494 1.5906 0.00011 -0.00031 -0.00016
+ 1.72497 0.68789 -0.03781 -0.91289 0.000341 0.000139 0.000345 2.4524 1.5903 0.00012 -0.00031 -0.00017
+ 1.74177 0.70460 -0.03619 -0.91468 0.000347 0.000141 0.000339 2.4555 1.5899 0.00012 -0.00032 -0.00017
+ 1.75856 0.72132 -0.03453 -0.91634 0.000353 0.000143 0.000333 2.4589 1.5895 0.00012 -0.00032 -0.00018
+ 1.77536 0.73803 -0.03284 -0.91788 0.000359 0.000146 0.000327 2.4626 1.5891 0.00012 -0.00033 -0.00018
+ 1.79217 0.75475 -0.03111 -0.91927 0.000365 0.000148 0.000320 2.4667 1.5886 0.00013 -0.00034 -0.00018
+ 1.80898 0.77146 -0.02934 -0.92053 0.000371 0.000150 0.000314 2.4713 1.5881 0.00013 -0.00034 -0.00019
+ 1.82579 0.78817 -0.02753 -0.92164 0.000377 0.000152 0.000307 2.4766 1.5875 0.00013 -0.00035 -0.00019
+ 1.84259 0.80488 -0.02569 -0.92258 0.000384 0.000155 0.000300 2.4828 1.5868 0.00013 -0.00035 -0.00019
+ 1.85940 0.82158 -0.02381 -0.92334 0.000391 0.000157 0.000293 2.4903 1.5859 0.00013 -0.00036 -0.00020
+ 1.87620 0.83827 -0.02189 -0.92389 0.000399 0.000159 0.000285 2.4993 1.5849 0.00014 -0.00037 -0.00020
+ 1.89299 0.85494 -0.01994 -0.92419 0.000407 0.000162 0.000276 2.5104 1.5837 0.00014 -0.00038 -0.00020
+ 1.90976 0.87160 -0.01796 -0.92421 0.000416 0.000165 0.000266 2.5244 1.5822 0.00014 -0.00038 -0.00021
+ 1.92649 0.88821 -0.01594 -0.92388 0.000426 0.000167 0.000254 2.5424 1.5804 0.00014 -0.00039 -0.00021
+ 1.94315 0.90474 -0.01389 -0.92311 0.000438 0.000171 0.000241 2.5664 1.5779 0.00015 -0.00040 -0.00021
+ 1.95966 0.92112 -0.01182 -0.92178 0.000452 0.000174 0.000224 2.5996 1.5747 0.00015 -0.00042 -0.00021
+ 1.97587 0.93720 -0.00974 -0.91972 0.000471 0.000178 0.000203 2.6478 1.5703 0.00015 -0.00043 -0.00022
+ 1.99152 0.95272 -0.00771 -0.91670 0.000496 0.000182 0.000175 2.7219 1.5642 0.00015 -0.00045 -0.00022
+ 2.00620 0.96727 -0.00576 -0.91250 0.000534 0.000188 0.000138 2.8461 1.5554 0.00016 -0.00049 -0.00022
+ 2.01942 0.98037 -0.00398 -0.90747 0.000592 0.000194 0.000090 3.0548 1.5443 0.00016 -0.00054 -0.00022
+ 2.03084 0.99168 -0.00242 -0.90249 0.000693 0.000200 0.000031 3.4652 1.5323 0.00016 -0.00063 -0.00023
+ 2.03924 1.00000 -0.00126 -0.90089 0.000707 0.000202 0.000027 3.4990 1.5085 0.00016 -0.00064 -0.00022
+ 2.03924 1.00010 -0.00000 0.90089 0.018045 0.007588 0.000000 2.3782
+ 2.04764 1.00850 0.00004 0.89078 0.017063 0.007961 0.000000 2.1432
+ 2.05721 1.01807 0.00020 0.89508 0.015878 0.007806 0.000000 2.0340
+ 2.06811 1.02897 0.00046 0.90013 0.014758 0.007633 0.000000 1.9334
+ 2.08053 1.04138 0.00084 0.90583 0.013713 0.007448 0.000000 1.8411
+ 2.09468 1.05551 0.00135 0.91197 0.012752 0.007259 0.000000 1.7566
+ 2.11079 1.07161 0.00204 0.91837 0.011880 0.007073 0.000000 1.6795
+ 2.12915 1.08995 0.00293 0.92482 0.011095 0.006895 0.000000 1.6092
+ 2.15005 1.11082 0.00406 0.93114 0.010395 0.006729 0.000000 1.5448
+ 2.17387 1.13460 0.00549 0.93722 0.009772 0.006577 0.000000 1.4858
+ 2.20100 1.16167 0.00727 0.94298 0.009218 0.006439 0.000000 1.4316
+ 2.23190 1.19249 0.00949 0.94840 0.008726 0.006315 0.000000 1.3818
+ 2.26711 1.22759 0.01222 0.95345 0.008288 0.006203 0.000000 1.3361
+ 2.30721 1.26755 0.01557 0.95814 0.007899 0.006103 0.000000 1.2943
+ 2.35289 1.31305 0.01967 0.96250 0.007553 0.006013 0.000000 1.2562
+ 2.40493 1.36485 0.02465 0.96652 0.007248 0.005932 0.000000 1.2218
+ 2.46420 1.42381 0.03069 0.97023 0.006978 0.005860 0.000000 1.1909
+ 2.53172 1.49094 0.03797 0.97364 0.006741 0.005795 0.000000 1.1634
+ 2.60863 1.56735 0.04674 0.97677 0.006535 0.005736 0.000000 1.1392
+ 2.69625 1.65433 0.05725 0.97962 0.006356 0.005684 0.000000 1.1182
+ 2.79605 1.75334 0.06981 0.98220 0.006203 0.005638 0.000000 1.1002
+ 2.90974 1.86604 0.08477 0.98451 0.006073 0.005597 0.000000 1.0850
+ 3.03924 1.99432 0.10254 0.98694 0.005955 0.005555 0.000000 1.0721
diff --git a/dart/viscous/XfoilFiles/case15Xfoil.dat b/dart/viscous/XfoilFiles/case15Xfoil.dat
new file mode 100644
index 0000000..109fa31
--- /dev/null
+++ b/dart/viscous/XfoilFiles/case15Xfoil.dat
@@ -0,0 +1,186 @@
+1.3
+0.1
+# s x y Ue/Vinf Dstar Theta Cf H H* P m K tau Di
+ 0.00000 1.00000 0.00126 1.05932 0.098795 0.013714 -0.000005 7.1548 1.5916 0.01539 0.10466 0.02595
+ 0.00840 0.99168 0.00242 1.06069 0.097302 0.013552 -0.000005 7.1304 1.5911 0.01525 0.10321 0.02573
+ 0.01982 0.98037 0.00398 1.06257 0.095216 0.013335 -0.000005 7.0907 1.5904 0.01506 0.10117 0.02544
+ 0.03304 0.96727 0.00576 1.06477 0.092757 0.013088 -0.000005 7.0379 1.5895 0.01484 0.09876 0.02511
+ 0.04772 0.95272 0.00771 1.06724 0.090009 0.012817 -0.000005 6.9732 1.5884 0.01460 0.09606 0.02475
+ 0.06337 0.93720 0.00974 1.06992 0.087069 0.012534 -0.000006 6.8978 1.5872 0.01435 0.09316 0.02436
+ 0.07959 0.92112 0.01182 1.07274 0.084023 0.012245 -0.000006 6.8132 1.5858 0.01409 0.09014 0.02397
+ 0.09610 0.90474 0.01389 1.07566 0.080930 0.011956 -0.000006 6.7204 1.5843 0.01383 0.08705 0.02357
+ 0.11275 0.88821 0.01594 1.07867 0.077822 0.011670 -0.000006 6.6202 1.5826 0.01358 0.08394 0.02318
+ 0.12949 0.87160 0.01796 1.08174 0.074719 0.011389 -0.000006 6.5126 1.5808 0.01333 0.08083 0.02279
+ 0.14626 0.85494 0.01994 1.08489 0.071631 0.011113 -0.000007 6.3979 1.5790 0.01308 0.07771 0.02241
+ 0.16305 0.83827 0.02189 1.08809 0.068565 0.010843 -0.000007 6.2759 1.5770 0.01284 0.07460 0.02203
+ 0.17984 0.82158 0.02381 1.09136 0.065525 0.010580 -0.000007 6.1465 1.5749 0.01260 0.07151 0.02166
+ 0.19665 0.80488 0.02569 1.09468 0.062515 0.010323 -0.000007 6.0095 1.5726 0.01237 0.06843 0.02130
+ 0.21346 0.78817 0.02753 1.09807 0.059538 0.010073 -0.000008 5.8648 1.5702 0.01215 0.06538 0.02094
+ 0.23027 0.77146 0.02934 1.10152 0.056596 0.009830 -0.000008 5.7122 1.5677 0.01193 0.06234 0.02060
+ 0.24707 0.75475 0.03111 1.10503 0.053693 0.009594 -0.000008 5.5516 1.5650 0.01172 0.05933 0.02026
+ 0.26388 0.73803 0.03284 1.10860 0.050830 0.009366 -0.000008 5.3828 1.5621 0.01151 0.05635 0.01993
+ 0.28068 0.72132 0.03453 1.11222 0.048011 0.009146 -0.000008 5.2059 1.5591 0.01131 0.05340 0.01962
+ 0.29748 0.70460 0.03619 1.11590 0.045238 0.008934 -0.000009 5.0207 1.5558 0.01113 0.05048 0.01931
+ 0.31427 0.68789 0.03781 1.11962 0.042515 0.008731 -0.000008 4.8275 1.5522 0.01094 0.04760 0.01902
+ 0.33105 0.67118 0.03938 1.12337 0.039845 0.008537 -0.000008 4.6262 1.5483 0.01077 0.04476 0.01874
+ 0.34783 0.65447 0.04092 1.12714 0.037232 0.008353 -0.000007 4.4172 1.5440 0.01061 0.04197 0.01847
+ 0.36460 0.63777 0.04241 1.13089 0.034681 0.008180 -0.000005 4.2005 1.5392 0.01046 0.03922 0.01821
+ 0.38135 0.62108 0.04387 1.13456 0.032198 0.008020 -0.000003 3.9764 1.5335 0.01032 0.03653 0.01796
+ 0.39809 0.60440 0.04527 1.13806 0.029790 0.007877 0.000002 3.7447 1.5268 0.01020 0.03390 0.01773
+ 0.41482 0.58772 0.04663 1.14117 0.027470 0.007757 0.000009 3.5053 1.5183 0.01010 0.03135 0.01750
+ 0.43154 0.57106 0.04794 1.14374 0.025254 0.007665 0.000018 3.2604 1.5079 0.01003 0.02888 0.01729
+ 0.44823 0.55441 0.04920 1.14753 0.023141 0.007535 0.000101 3.0377 1.5008 0.00992 0.02656 0.01709
+ 0.46491 0.53778 0.05041 1.15493 0.021132 0.007299 0.000188 2.8623 1.5027 0.00974 0.02441 0.01690
+ 0.48157 0.52116 0.05156 1.16404 0.019265 0.007029 0.000284 2.7082 1.5073 0.00952 0.02243 0.01671
+ 0.49820 0.50456 0.05265 1.17482 0.017545 0.006735 0.000392 2.5730 1.5141 0.00930 0.02061 0.01653
+ 0.51481 0.48798 0.05368 1.18716 0.015972 0.006423 0.000512 2.4544 1.5224 0.00905 0.01896 0.01636
+ 0.53140 0.47143 0.05465 1.20096 0.014543 0.006103 0.000646 2.3507 1.5318 0.00880 0.01747 0.01619
+ 0.54796 0.45489 0.05555 1.21608 0.013252 0.005781 0.000792 2.2599 1.5417 0.00855 0.01612 0.01603
+ 0.56448 0.43839 0.05639 1.23239 0.012090 0.005462 0.000951 2.1807 1.5519 0.00830 0.01490 0.01587
+ 0.58098 0.42191 0.05714 1.24974 0.011047 0.005151 0.001121 2.1114 1.5620 0.00805 0.01381 0.01571
+ 0.59744 0.40546 0.05782 1.26803 0.010112 0.004851 0.001300 2.0510 1.5719 0.00780 0.01282 0.01555
+ 0.61386 0.38905 0.05842 1.28714 0.009273 0.004563 0.001489 1.9981 1.5814 0.00756 0.01194 0.01539
+ 0.63024 0.37268 0.05893 1.30700 0.008519 0.004288 0.001686 1.9518 1.5904 0.00732 0.01113 0.01523
+ 0.64657 0.35635 0.05935 1.32756 0.007839 0.004026 0.001892 1.9110 1.5989 0.00710 0.01041 0.01506
+ 0.66286 0.34007 0.05967 1.34880 0.007223 0.003778 0.002105 1.8750 1.6070 0.00687 0.00974 0.01490
+ 0.67910 0.32383 0.05989 1.37073 0.006664 0.003543 0.002327 1.8431 1.6145 0.00666 0.00913 0.01473
+ 0.69527 0.30766 0.06000 1.39338 0.006152 0.003319 0.002558 1.8146 1.6216 0.00644 0.00857 0.01456
+ 0.71139 0.29154 0.06000 1.41682 0.005682 0.003107 0.002800 1.7890 1.6282 0.00624 0.00805 0.01439
+ 0.72743 0.27550 0.05988 1.44114 0.005247 0.002904 0.003055 1.7658 1.6346 0.00603 0.00756 0.01421
+ 0.74340 0.25953 0.05963 1.46644 0.004843 0.002710 0.003324 1.7446 1.6406 0.00583 0.00710 0.01402
+ 0.75928 0.24366 0.05924 1.49288 0.004466 0.002525 0.003612 1.7250 1.6463 0.00563 0.00667 0.01383
+ 0.77506 0.22788 0.05871 1.52061 0.004112 0.002347 0.003921 1.7067 1.6519 0.00543 0.00625 0.01363
+ 0.79074 0.21222 0.05803 1.54986 0.003778 0.002176 0.004257 1.6894 1.6574 0.00523 0.00586 0.01343
+ 0.80628 0.19670 0.05718 1.58087 0.003462 0.002011 0.004626 1.6728 1.6628 0.00503 0.00547 0.01321
+ 0.82167 0.18135 0.05616 1.61393 0.003161 0.001852 0.005034 1.6567 1.6683 0.00482 0.00510 0.01299
+ 0.83688 0.16619 0.05497 1.64941 0.002875 0.001698 0.005492 1.6409 1.6739 0.00462 0.00474 0.01275
+ 0.85186 0.15127 0.05358 1.68771 0.002601 0.001549 0.006012 1.6251 1.6796 0.00441 0.00439 0.01250
+ 0.86655 0.13666 0.05200 1.72930 0.002341 0.001404 0.006607 1.6092 1.6856 0.00420 0.00405 0.01224
+ 0.88087 0.12246 0.05022 1.77470 0.002092 0.001265 0.007296 1.5931 1.6919 0.00399 0.00371 0.01197
+ 0.89470 0.10877 0.04824 1.82439 0.001857 0.001132 0.008097 1.5769 1.6984 0.00377 0.00339 0.01167
+ 0.90789 0.09575 0.04609 1.87878 0.001636 0.001005 0.009027 1.5608 1.7052 0.00355 0.00307 0.01137
+ 0.92028 0.08358 0.04381 1.93804 0.001433 0.000886 0.010101 1.5455 1.7121 0.00333 0.00278 0.01104
+ 0.93169 0.07242 0.04142 2.00202 0.001248 0.000776 0.011315 1.5319 1.7186 0.00311 0.00250 0.01071
+ 0.94200 0.06239 0.03900 2.07027 0.001085 0.000677 0.012643 1.5216 1.7242 0.00290 0.00225 0.01036
+ 0.95118 0.05354 0.03659 2.14208 0.000944 0.000588 0.014020 1.5170 1.7281 0.00270 0.00202 0.01000
+ 0.95926 0.04581 0.03424 2.21666 0.000825 0.000511 0.015316 1.5217 1.7288 0.00251 0.00183 0.00962
+ 0.96635 0.03910 0.03196 2.29321 0.000729 0.000444 0.016281 1.5427 1.7242 0.00233 0.00167 0.00922
+ 0.97255 0.03329 0.02976 2.37100 0.000657 0.000386 0.016447 1.5936 1.7105 0.00217 0.00156 0.00880
+ 0.97802 0.02826 0.02764 2.45021 0.000614 0.000336 0.014948 1.7079 1.6806 0.00202 0.00151 0.00830
+ 0.98286 0.02387 0.02560 2.53708 0.000613 0.000290 0.010589 1.9765 1.6237 0.00187 0.00156 0.00769
+ 0.98718 0.02003 0.02361 2.66354 0.000672 0.000234 0.004055 2.6929 1.5416 0.00166 0.00179 0.00682
+ 0.99107 0.01666 0.02167 2.90004 0.000773 0.000143 -0.001365 5.1073 1.5329 0.00120 0.00224 0.00533
+ 0.99460 0.01369 0.01977 3.11628 0.000790 0.000074 -0.002441 10.1254 1.5774 0.00072 0.00246 0.00352
+ 0.99781 0.01108 0.01789 3.12974 0.000672 0.000071 -0.002607 8.9761 1.5638 0.00069 0.00210 0.00338
+ 1.00078 0.00879 0.01602 3.14572 0.000517 0.000067 -0.002879 7.1865 1.5448 0.00067 0.00163 0.00325
+ 1.00352 0.00678 0.01415 3.17160 0.000367 0.000064 -0.003364 5.3916 1.5310 0.00064 0.00117 0.00310
+ 1.00607 0.00505 0.01227 3.21705 0.000251 0.000058 -0.000944 3.9713 1.5282 0.00060 0.00081 0.00296
+ 1.00847 0.00358 0.01038 3.28586 0.000176 0.000051 0.010426 3.1140 1.5419 0.00056 0.00058 0.00281
+ 1.01072 0.00236 0.00848 3.35403 0.000134 0.000045 0.025711 2.7081 1.5653 0.00050 0.00045 0.00263
+ 1.01286 0.00140 0.00657 3.38367 0.000109 0.000039 0.040838 2.4971 1.5852 0.00045 0.00037 0.00240
+ 1.01490 0.00070 0.00466 3.35732 0.000093 0.000035 0.054172 2.3716 1.6003 0.00039 0.00031 0.00212
+ 1.01685 0.00024 0.00276 3.27022 0.000082 0.000032 0.064209 2.2927 1.6113 0.00034 0.00027 0.00182
+ 1.01872 0.00003 0.00091 3.13037 0.000075 0.000030 0.069038 2.2525 1.6174 0.00030 0.00024 0.00151
+ 1.02053 0.00003 -0.00091 2.95597 0.000071 0.000029 0.070444 2.2275 1.6214 0.00026 0.00021 0.00122
+ 1.02240 0.00024 -0.00276 2.74547 0.000068 0.000029 0.069738 2.2016 1.6256 0.00022 0.00019 0.00096
+ 1.02434 0.00070 -0.00466 2.50115 0.000066 0.000029 0.065539 2.1832 1.6287 0.00018 0.00017 0.00073
+ 1.02638 0.00140 -0.00657 2.23865 0.000067 0.000029 0.058576 2.1719 1.6307 0.00015 0.00015 0.00053
+ 1.02852 0.00236 -0.00848 1.97374 0.000068 0.000030 0.050142 2.1663 1.6317 0.00012 0.00013 0.00038
+ 1.03078 0.00358 -0.01038 1.71902 0.000072 0.000032 0.041499 2.1647 1.6319 0.00010 0.00012 0.00027
+ 1.03317 0.00505 -0.01227 1.48189 0.000076 0.000034 0.033479 2.1661 1.6317 0.00008 0.00011 0.00018
+ 1.03573 0.00678 -0.01415 1.26536 0.000082 0.000037 0.026496 2.1691 1.6312 0.00006 0.00010 0.00012
+ 1.03847 0.00879 -0.01602 1.06957 0.000088 0.000040 0.020608 2.1741 1.6303 0.00005 0.00009 0.00008
+ 1.04143 0.01108 -0.01789 0.89299 0.000095 0.000043 0.015809 2.1784 1.6296 0.00003 0.00009 0.00005
+ 1.04465 0.01369 -0.01977 0.73367 0.000104 0.000047 0.011848 2.1868 1.6281 0.00003 0.00008 0.00003
+ 1.04817 0.01666 -0.02167 0.58909 0.000112 0.000051 0.008754 2.1886 1.6278 0.00002 0.00007 0.00002
+ 1.05206 0.02003 -0.02361 0.45753 0.000123 0.000056 0.006135 2.2038 1.6253 0.00001 0.00006 0.00001
+ 1.05638 0.02387 -0.02560 0.33657 0.000133 0.000060 0.004199 2.1973 1.6264 0.00001 0.00004 0.00000
+ 1.06122 0.02826 -0.02764 0.22525 0.000146 0.000066 0.002489 2.2243 1.6219 0.00000 0.00003 0.00000
+ 1.06669 0.03329 -0.02976 0.12153 0.000155 0.000070 0.001276 2.2155 1.6233 0.00000 0.00002 0.00000
+ 1.07290 0.03910 -0.03196 0.02521 0.000181 0.000081 0.000224 2.2295 1.6211 0.00000 0.00000 0.00000
+ 1.07998 0.04581 -0.03424 -0.06483 0.000181 0.000081 0.000576 2.2295 1.6211 0.00000 -0.00001 -0.00000
+ 1.08806 0.05354 -0.03659 -0.14904 0.000209 0.000092 0.001115 2.2610 1.6161 0.00000 -0.00003 -0.00000
+ 1.09724 0.06239 -0.03900 -0.22671 0.000222 0.000100 0.001632 2.2333 1.6204 0.00001 -0.00005 -0.00000
+ 1.10756 0.07242 -0.04142 -0.29837 0.000247 0.000108 0.001869 2.2729 1.6143 0.00001 -0.00007 -0.00000
+ 1.11897 0.08358 -0.04381 -0.36301 0.000265 0.000117 0.002142 2.2584 1.6165 0.00002 -0.00010 -0.00001
+ 1.13135 0.09575 -0.04609 -0.42114 0.000290 0.000127 0.002230 2.2805 1.6131 0.00002 -0.00012 -0.00002
+ 1.14454 0.10877 -0.04824 -0.47248 0.000312 0.000136 0.002328 2.2798 1.6132 0.00003 -0.00015 -0.00002
+ 1.15837 0.12246 -0.05022 -0.51782 0.000335 0.000146 0.002348 2.2906 1.6116 0.00004 -0.00017 -0.00003
+ 1.17269 0.13666 -0.05200 -0.55763 0.000358 0.000156 0.002356 2.2956 1.6109 0.00005 -0.00020 -0.00004
+ 1.18738 0.15127 -0.05358 -0.59272 0.000381 0.000165 0.002339 2.3022 1.6099 0.00006 -0.00023 -0.00006
+ 1.20236 0.16619 -0.05497 -0.62371 0.000404 0.000174 0.002311 2.3080 1.6091 0.00007 -0.00025 -0.00007
+ 1.21757 0.18135 -0.05616 -0.65121 0.000426 0.000184 0.002275 2.3136 1.6083 0.00008 -0.00028 -0.00008
+ 1.23296 0.19670 -0.05718 -0.67573 0.000449 0.000193 0.002233 2.3189 1.6075 0.00009 -0.00030 -0.00010
+ 1.24851 0.21222 -0.05803 -0.69771 0.000470 0.000201 0.002189 2.3239 1.6068 0.00010 -0.00033 -0.00011
+ 1.26418 0.22788 -0.05871 -0.71749 0.000492 0.000210 0.002143 2.3288 1.6061 0.00011 -0.00035 -0.00012
+ 1.27996 0.24366 -0.05924 -0.73537 0.000513 0.000219 0.002096 2.3334 1.6055 0.00012 -0.00038 -0.00014
+ 1.29584 0.25953 -0.05963 -0.75161 0.000534 0.000227 0.002050 2.3379 1.6049 0.00013 -0.00040 -0.00015
+ 1.31181 0.27550 -0.05988 -0.76642 0.000555 0.000236 0.002004 2.3422 1.6043 0.00014 -0.00043 -0.00017
+ 1.32785 0.29154 -0.06000 -0.77997 0.000575 0.000244 0.001960 2.3463 1.6037 0.00015 -0.00045 -0.00019
+ 1.34397 0.30766 -0.06000 -0.79241 0.000595 0.000252 0.001917 2.3502 1.6032 0.00016 -0.00047 -0.00020
+ 1.36015 0.32383 -0.05989 -0.80387 0.000615 0.000260 0.001875 2.3540 1.6027 0.00017 -0.00049 -0.00022
+ 1.37638 0.34007 -0.05967 -0.81447 0.000635 0.000268 0.001835 2.3575 1.6022 0.00018 -0.00052 -0.00023
+ 1.39267 0.35635 -0.05935 -0.82430 0.000654 0.000275 0.001796 2.3609 1.6017 0.00019 -0.00054 -0.00025
+ 1.40900 0.37268 -0.05893 -0.83344 0.000673 0.000283 0.001759 2.3641 1.6013 0.00020 -0.00056 -0.00026
+ 1.42538 0.38905 -0.05842 -0.84197 0.000692 0.000290 0.001724 2.3672 1.6009 0.00021 -0.00058 -0.00028
+ 1.44181 0.40546 -0.05782 -0.84995 0.000711 0.000298 0.001690 2.3700 1.6005 0.00022 -0.00060 -0.00029
+ 1.45826 0.42191 -0.05714 -0.85745 0.000729 0.000305 0.001658 2.3726 1.6002 0.00022 -0.00062 -0.00031
+ 1.47476 0.43839 -0.05639 -0.86451 0.000747 0.000312 0.001628 2.3750 1.5999 0.00023 -0.00065 -0.00032
+ 1.49129 0.45489 -0.05555 -0.87117 0.000764 0.000319 0.001599 2.3772 1.5996 0.00024 -0.00067 -0.00034
+ 1.50784 0.47143 -0.05465 -0.87749 0.000782 0.000326 0.001572 2.3792 1.5993 0.00025 -0.00069 -0.00035
+ 1.52443 0.48798 -0.05368 -0.88349 0.000799 0.000333 0.001547 2.3809 1.5991 0.00026 -0.00071 -0.00037
+ 1.54104 0.50456 -0.05265 -0.88921 0.000815 0.000340 0.001523 2.3824 1.5989 0.00027 -0.00072 -0.00038
+ 1.55768 0.52116 -0.05156 -0.89469 0.000831 0.000346 0.001501 2.3837 1.5988 0.00028 -0.00074 -0.00040
+ 1.57433 0.53778 -0.05041 -0.89994 0.000847 0.000353 0.001480 2.3847 1.5986 0.00029 -0.00076 -0.00041
+ 1.59101 0.55441 -0.04920 -0.90499 0.000863 0.000359 0.001460 2.3855 1.5985 0.00029 -0.00078 -0.00043
+ 1.60771 0.57106 -0.04794 -0.90987 0.000878 0.000365 0.001443 2.3860 1.5985 0.00030 -0.00080 -0.00044
+ 1.62442 0.58772 -0.04663 -0.91459 0.000892 0.000371 0.001426 2.3862 1.5984 0.00031 -0.00082 -0.00045
+ 1.64115 0.60440 -0.04527 -0.91919 0.000907 0.000377 0.001411 2.3862 1.5984 0.00032 -0.00083 -0.00047
+ 1.65789 0.62108 -0.04387 -0.92366 0.000920 0.000383 0.001397 2.3859 1.5985 0.00033 -0.00085 -0.00048
+ 1.67465 0.63777 -0.04241 -0.92804 0.000934 0.000388 0.001384 2.3853 1.5985 0.00033 -0.00087 -0.00050
+ 1.69141 0.65447 -0.04092 -0.93233 0.000946 0.000394 0.001373 2.3844 1.5987 0.00034 -0.00088 -0.00051
+ 1.70819 0.67118 -0.03938 -0.93656 0.000959 0.000399 0.001363 2.3833 1.5988 0.00035 -0.00090 -0.00052
+ 1.72497 0.68789 -0.03781 -0.94073 0.000971 0.000404 0.001355 2.3818 1.5990 0.00036 -0.00091 -0.00054
+ 1.74177 0.70460 -0.03619 -0.94486 0.000982 0.000409 0.001348 2.3800 1.5992 0.00037 -0.00093 -0.00055
+ 1.75856 0.72132 -0.03453 -0.94897 0.000992 0.000414 0.001342 2.3780 1.5995 0.00037 -0.00094 -0.00057
+ 1.77536 0.73803 -0.03284 -0.95306 0.001003 0.000418 0.001337 2.3756 1.5998 0.00038 -0.00096 -0.00058
+ 1.79217 0.75475 -0.03111 -0.95715 0.001012 0.000423 0.001334 2.3728 1.6002 0.00039 -0.00097 -0.00059
+ 1.80898 0.77146 -0.02934 -0.96126 0.001021 0.000427 0.001332 2.3697 1.6006 0.00039 -0.00098 -0.00061
+ 1.82579 0.78817 -0.02753 -0.96540 0.001029 0.000431 0.001332 2.3661 1.6011 0.00040 -0.00099 -0.00062
+ 1.84259 0.80488 -0.02569 -0.96959 0.001037 0.000435 0.001334 2.3621 1.6016 0.00041 -0.00100 -0.00063
+ 1.85940 0.82158 -0.02381 -0.97384 0.001043 0.000439 0.001337 2.3575 1.6022 0.00042 -0.00102 -0.00065
+ 1.87620 0.83827 -0.02189 -0.97819 0.001049 0.000442 0.001343 2.3523 1.6029 0.00042 -0.00103 -0.00066
+ 1.89299 0.85494 -0.01994 -0.98265 0.001053 0.000445 0.001351 2.3464 1.6037 0.00043 -0.00103 -0.00068
+ 1.90976 0.87160 -0.01796 -0.98727 0.001057 0.000447 0.001362 2.3395 1.6047 0.00044 -0.00104 -0.00069
+ 1.92649 0.88821 -0.01594 -0.99209 0.001058 0.000450 0.001377 2.3314 1.6058 0.00044 -0.00105 -0.00071
+ 1.94315 0.90474 -0.01389 -0.99717 0.001058 0.000452 0.001397 2.3218 1.6071 0.00045 -0.00106 -0.00072
+ 1.95966 0.92112 -0.01182 -1.00258 0.001056 0.000453 0.001424 2.3102 1.6088 0.00045 -0.00106 -0.00073
+ 1.97587 0.93720 -0.00974 -1.00839 0.001050 0.000453 0.001459 2.2959 1.6109 0.00046 -0.00106 -0.00075
+ 1.99152 0.95272 -0.00771 -1.01471 0.001040 0.000452 0.001507 2.2780 1.6135 0.00047 -0.00106 -0.00076
+ 2.00620 0.96727 -0.00576 -1.02155 0.001026 0.000450 0.001570 2.2563 1.6168 0.00047 -0.00105 -0.00078
+ 2.01942 0.98037 -0.00398 -1.02941 0.001002 0.000446 0.001667 2.2246 1.6219 0.00047 -0.00103 -0.00079
+ 2.03084 0.99168 -0.00242 -1.03639 0.000985 0.000443 0.001739 2.2037 1.6253 0.00048 -0.00102 -0.00080
+ 2.03924 1.00000 -0.00126 -1.05932 0.000851 0.000412 0.002361 2.0417 1.6292 0.00046 -0.00090 -0.00080
+ 2.03924 1.00010 -0.00000 1.05932 0.102142 0.014126 0.000000 7.1814
+ 2.04764 1.00850 0.00005 1.05668 0.101482 0.014449 0.000000 6.9751
+ 2.05721 1.01807 0.00025 1.05321 0.100490 0.014875 0.000000 6.7089
+ 2.06811 1.02896 0.00058 1.04914 0.099226 0.015377 0.000000 6.4082
+ 2.08053 1.04137 0.00105 1.04434 0.097649 0.015970 0.000000 6.0720
+ 2.09468 1.05550 0.00171 1.03864 0.095706 0.016674 0.000000 5.6994
+ 2.11079 1.07159 0.00257 1.03182 0.093332 0.017517 0.000000 5.2901
+ 2.12915 1.08991 0.00369 1.02357 0.090453 0.018539 0.000000 4.8441
+ 2.15005 1.11077 0.00512 1.01340 0.086989 0.019801 0.000000 4.3611
+ 2.17387 1.13452 0.00691 1.00032 0.082881 0.021430 0.000000 3.8385
+ 2.20100 1.16156 0.00916 0.98150 0.078185 0.023808 0.000000 3.2587
+ 2.23190 1.19234 0.01194 0.96339 0.072701 0.026128 0.000000 2.7601
+ 2.26711 1.22737 0.01538 0.95197 0.066571 0.027589 0.000000 2.3927
+ 2.30721 1.26725 0.01959 0.94497 0.060272 0.028463 0.000000 2.0987
+ 2.35289 1.31264 0.02473 0.94231 0.054175 0.028784 0.000000 1.8644
+ 2.40493 1.36430 0.03098 0.94340 0.048605 0.028658 0.000000 1.6789
+ 2.46420 1.42309 0.03855 0.94725 0.043774 0.028243 0.000000 1.5333
+ 2.53172 1.48999 0.04769 0.95267 0.039755 0.027691 0.000000 1.4193
+ 2.60863 1.56612 0.05868 0.95863 0.036505 0.027118 0.000000 1.3300
+ 2.69625 1.65274 0.07184 0.96443 0.033918 0.026588 0.000000 1.2596
+ 2.79605 1.75129 0.08757 0.96972 0.031871 0.026125 0.000000 1.2039
+ 2.90974 1.86343 0.10629 0.97438 0.030252 0.025733 0.000000 1.1597
+ 3.03924 1.99101 0.12852 0.98249 0.028567 0.025075 0.000000 1.1232
diff --git a/dart/viscous/__init__.py b/dart/viscous/__init__.py
deleted file mode 100644
index 40a96af..0000000
--- a/dart/viscous/__init__.py
+++ /dev/null
@@ -1 +0,0 @@
-# -*- coding: utf-8 -*-
diff --git a/dart/viscous/airfoil.py b/dart/viscous/airfoil.py
old mode 100644
new mode 100755
index 2bd069d..31a5796
--- a/dart/viscous/airfoil.py
+++ b/dart/viscous/airfoil.py
@@ -19,7 +19,7 @@
## Airfoil around which the boundary layer is computed
# Amaury Bilocq
-from dart.viscous.boundary import Boundary
+from dart.viscous.Boundary import Boundary
import numpy as np
class Airfoil(Boundary):
@@ -34,7 +34,7 @@ class Airfoil(Boundary):
self.HEnd = [0,0]
self.CtEnd = [0,0]
self.nEnd = [0,0]
-
+
def initialConditions(self, Re, dv):
if dv > 0:
self.T0 = np.sqrt(0.075/(Re*dv))
@@ -45,6 +45,7 @@ class Airfoil(Boundary):
self.Ct0 = 0
return self.T0, self.H0, self.n0, self.Ct0
+
def connectList(self):
''' Sort the value read by the viscous solver/ Create list of connectivity
'''
diff --git a/dart/viscous/boundary.py b/dart/viscous/boundary.py
old mode 100644
new mode 100755
diff --git a/dart/viscous/coupler.py b/dart/viscous/coupler.py
old mode 100644
new mode 100755
index 269348d..aa4ce37
--- a/dart/viscous/coupler.py
+++ b/dart/viscous/coupler.py
@@ -21,62 +21,84 @@
import numpy as np
from fwk.coloring import ccolors
-from dart.viscous.airfoil import Airfoil
-from dart.viscous.wake import Wake
+from dart.viscous.Airfoil import Airfoil
+from dart.viscous.Wake import Wake
+import tbox.utils as tboxU
+import dart.utils as floU
class Coupler:
- def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer):
+ def __init__(self, _isolver, _vsolver, _boundaryAirfoil, _boundaryWake, _tol, _writer, _bnd, _timeParam=None):
+ self.bnd=_bnd
self.isolver =_isolver # inviscid solver
self.vsolver = _vsolver # viscous solver
+
self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)] # airfoil and wake python objects
self.tol = _tol # tolerance of the VII
self.writer = _writer
-
+ self.timeParam=_timeParam
def run(self):
''' Perform coupling
'''
- # initialize loop
+ # Initialize loop
it = 0
converged = 0 # temp
CdOld = self.vsolver.Cd
+ print('+---------------------------- VII Solver begin --------------------------------+')
+ print('+------------------------------------------------------------------------------+')
+ print('|{:>79s}'.format('|'))
while converged == 0:
- print(ccolors.ANSI_BLUE + 'Iteration: ', it, ccolors.ANSI_RESET)
+ print('|', ccolors.ANSI_BLUE + 'Iteration: ', it, ccolors.ANSI_RESET, '{:>63s}'.format('|'))
# run inviscid solver
- self.isolver.run()
- print('--- Viscous solver parameters ---')
- print('Tolerance (drag count):', self.tol * 1000)
- print('--- Viscous problem definition ---')
- print('Reynolds number:', self.vsolver.Re)
- print('')
+ print('+----------------------------- Inviscid solver --------------------------------+')
+ self.isolver.run()
for n in range(0, len(self.group)):
- print('Computing for', self.group[n].name, '...', end = ' ')
+ # print('Computing for', self.group[n].name, '...', end = ' ')
for i in range(0, len(self.group[n].boundary.nodes)):
self.group[n].v[i,0] = self.isolver.U[self.group[n].boundary.nodes[i].row][0]
self.group[n].v[i,1] = self.isolver.U[self.group[n].boundary.nodes[i].row][1]
self.group[n].v[i,2] = self.isolver.U[self.group[n].boundary.nodes[i].row][2]
self.group[n].Me[i] = self.isolver.M[self.group[n].boundary.nodes[i].row]
self.group[n].rhoe[i] = self.isolver.rho[self.group[n].boundary.nodes[i].row]
- # run viscous solver
- self.vsolver.run(self.group[n])
- if self.group[n].name == 'airfoil':
- self.group[n+1].T0 = self.group[n].TEnd[0]+self.group[n].TEnd[1]
- self.group[n+1].H0 = (self.group[n].HEnd[0]*self.group[n].TEnd[0]+self.group[n].HEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0
- self.group[n+1].Ct0 = (self.group[n].CtEnd[0]*self.group[n].TEnd[0]+self.group[n].CtEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0
- # get blowing velocity from viscous solver and update inviscid solver
+
+ if it==0:
+ # 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', '')
+ Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
+ tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+ # run viscous solver
+ print('+------------------------------ Viscous solver --------------------------------+')
+ self.vsolver.run(self.group)
+ self.vsolver.writeFile()
+ for n in range(0, len(self.group)):
for i in range(0, self.group[n].nE):
self.group[n].boundary.setU(i, self.group[n].u[i])
- print('done!')
- print(' Iter Cd xtr_top xtr_bot Error')
- print('{0:8d} {1:12.5f} {2:12.5f} {3:12.5f} {4:12.5f}'.format(it, self.vsolver.Cd, self.vsolver.xtr[0], self.vsolver.xtr[1], abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd))
- print('')
+ # Write files
+ '''Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
+ tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+ self.isolver.save(self.writer)'''
+ #self.vsolver.writeFile()
+ # Output coupling statut
+ xtrT_type='forced' if self.vsolver.xtrF[0] is not None else 'free'
+ xtrB_type='forced' if self.vsolver.xtrF[1] is not None else 'free'
+ print('+---------------------------- Postprocessing ----------------------------------+')
+ print('|{0:>6s}| {1:>14s}| {2:>18s}| {3:>18s}| {4:>14s}|'.format('Iter', 'Cd', f'Top xtr ({xtrT_type})', f'Bot xtr ({xtrB_type})','Error'))
+ print('|{0:>6d}| {1:>14f}| {2:>18f}| {3:>18f}| {4:>14f}|'.format(it, self.vsolver.Cd, self.vsolver.xtr[0], self.vsolver.xtr[1], np.log10(abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd)))
+ print('|{:>79}'.format('|'))
+ print('+------------------------------------------------------------------------------+')
+
+ # Write files
# Converged or not
if abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd < self.tol:
converged = 1
else:
CdOld = self.vsolver.Cd
+
it += 1
self.vsolver.it += 1
# save results
self.isolver.save(self.writer)
self.vsolver.writeFile()
+ Cp = floU.extract(self.bnd.groups[0].tag.elems, self.isolver.Cp)
+ tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
print('\n')
diff --git a/dart/viscous/model/BIConditions.py b/dart/viscous/model/BIConditions.py
new file mode 100755
index 0000000..8ecbbff
--- /dev/null
+++ b/dart/viscous/model/BIConditions.py
@@ -0,0 +1,154 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# BIConditions: Initial and boundary conditions #
+# of the pseudo-time marching problem
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+from dart.viscous.model.Closures import Closures
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+class Initial:
+ def __init__(self, BoundaryCond):
+ self.boundCond=BoundaryCond
+ self.mu_air=1.849e-5 # Dynamic viscosity of air at 25°C
+
+ def initialConditions(self, var):
+ """Twaithes integral solution of the (steady) Von Karman momentum
+ integral equation for a laminar flow over a flat plate
+
+ Parameters
+ ----------
+ - var: class
+ Data container: - Solution
+ - Boundary layer parameters
+ - External flow parameters
+ - Mesh
+ - Transition information
+
+ Tasks
+ -----
+ - Compute Twaithes solution of the momentum integral equation for a laminar flow over a flat plate in the
+ regime prescibed by the inviscid solver.
+ - Define boundary conditions at the stagnation point and the first wake point
+ """
+ u0=np.zeros(var.nN*var.nVar)
+
+ for n in range(3): # 3 components (upper, lower sides, wake)
+ if n==2:
+ xx=var.xx[var.bNodes[n]:]
+ vt=var.extPar[var.bNodes[n]*var.nExtPar+2::var.nExtPar]
+ rho=var.extPar[var.bNodes[n]*var.nExtPar+1::var.nExtPar]
+ else:
+ xx=var.xx[var.bNodes[n]:var.bNodes[n+1]]
+ vt=var.extPar[var.bNodes[n]*var.nExtPar+2:var.bNodes[n+1]*var.nExtPar:var.nExtPar]
+ rho=var.extPar[var.bNodes[n]*var.nExtPar+1:var.bNodes[n+1]*var.nExtPar:var.nExtPar]
+
+ Rex=rho*vt*xx/self.mu_air
+ Rex[Rex==0]=1 # Stagnation point
+ ThetaTw=xx/np.sqrt(Rex)
+ # lambda : Dimensionless pressure-gradient parameter
+ lambdaTw=np.zeros(len(ThetaTw))
+ HTw=np.zeros(len(ThetaTw))
+ for idx in range(0,len(lambdaTw)):
+ idxPrev=idx-1
+ lambdaTw[idx]=(rho[idx]*ThetaTw[idx]**2)/self.mu_air*(vt[idx]-vt[idxPrev])#/(xx[idx]-xx[idxPrev])
+ # Empirical correlations for H=f(lambda)
+ if 0<=lambdaTw[idx]<0.1:
+ HTw[idx]=2.61-3.75*lambdaTw[idx]+5.24*lambdaTw[idx]**2
+ elif -0.1<lambdaTw[idx]<0:
+ HTw[idx]=2.088+0.0731/(lambdaTw[idx]+0.14)
+
+ if n==2:
+ u0[var.bNodes[n]*var.nVar+0::var.nVar]=ThetaTw # Theta
+ u0[var.bNodes[n]*var.nVar+1::var.nVar]=HTw # H
+ else:
+ u0[var.bNodes[n]*var.nVar+0:var.bNodes[n+1]*var.nVar:var.nVar]=ThetaTw # Theta
+ u0[var.bNodes[n]*var.nVar+1:var.bNodes[n+1]*var.nVar:var.nVar]=HTw # H
+
+ u0[1::var.nVar]=2 # H
+ u0[2::var.nVar]=0 # N
+ u0[3::var.nVar]=var.extPar[2::var.nExtPar] # Vt
+ u0[4::var.nVar]=0.001 # Ct
+
+ var.u=u0
+ self.boundCond.airfoilBC(var)
+ self.boundCond.wakeBC(var)
+
+
+class Boundaries:
+ """ ---------- Boundary conditions ----------
+
+ Boundary conditions at the stagnation point and on the first point of the wake.
+
+ Attributes
+ ----------
+ Re: Flow Reynolds number.
+ """
+ def __init__(self, _Re):
+ self.Re=_Re
+
+ # Airfoil boundary conditions
+ def airfoilBC(self, var):
+ """Boundary conditions at the stagnation point.
+ Twaithes solution method values used for the SP.
+ """
+ if var.dv[0] > 0:
+ T0 = np.sqrt(0.075/(self.Re*var.dv[0]))
+ else:
+ T0 = 1e-6
+ H0 = 2.23 # One parameter family Falkner Skan
+ n0 = 0
+ ue0=var.extPar[2]
+ Ct0 = 0
+
+ # Stagnation point
+ var.blPar[var.bNodes[0]*var.nBlPar+0] = ue0*T0*var.Re # Ret
+ 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)
+ 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')
+ # Cteq stays = 0
+ #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
+
+ var.u[var.bNodes[0]*var.nVar:(var.bNodes[0]+1)*var.nVar]=[T0, H0, n0, ue0, Ct0]
+ pass
+
+ # Wake boundary conditions
+ def wakeBC(self, var):
+ """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]
+ 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]
+ Ct0=(xEnd_up[0]*xEnd_up[4]+xEnd_lw[0]*xEnd_lw[4])/T0 # ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low)
+
+ var.blPar[var.bNodes[2]*var.nBlPar+0] = ue0*T0*var.Re
+ 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
+
+ var.u[var.bNodes[2]*var.nVar:(var.bNodes[2]+1)*var.nVar]=[T0, H0, n0, ue0, Ct0]
+ pass
diff --git a/dart/viscous/model/CSolver.py b/dart/viscous/model/CSolver.py
new file mode 100644
index 0000000..0ae15f2
--- /dev/null
+++ b/dart/viscous/model/CSolver.py
@@ -0,0 +1,322 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# CSolver : Boundary layer laws and linear system sparse #
+# matrices computation (Residuals, Jacobian) #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.11.11 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+from dart.viscous.model.Closures import Closures
+
+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 == dataCont.bNodes[1] - 1 or iMarker == dataCont.nN-1 or iMarker == dataCont.nN - 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
+ #
+ # 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]
+
+ if np.linalg.det(timeMatrix) == 0:
+ raise RuntimeError('Singular time matrix.')
+ 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[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
+
+
+class sysMatrix(flowEquations):
+ def __init__(self, _nMarker, _nMarkers, setGradients, fouScheme) -> None:
+ flowEquations.__init__(self, setGradients, fouScheme)
+ self._nMarker = _nMarker
+ self._nMarkers = _nMarkers
+ pass
+
+ class Jacobian:
+
+ def _initialize(nVar, nMarkers) -> np.ndarray:
+ return np.zeros((nVar * nMarkers, nVar * nMarkers))
+
+ class linSysRes:
+
+ def _initialize(nVar, nMarkers) -> np.ndarray:
+ return np.zeros(nVar * nMarkers)
+
+
+ def buildJacobian(self, dataCont, marker = -1, order = 1, linSysRes = None, 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.
+
+ - 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.
+
+ - eps (double, optional): Perturbation on the variables. Default : 10^(-8).
+
+ Returns
+ -------
+
+ - JacMatrix (np.ndarray): Jacobian of the system.
+
+ """
+
+ # Build Jacobian for whole domain
+ if marker == -1:
+ sysMatSize = self._nMarker * self._nMarkers
+ Jacobian = np.zeros((sysMatSize, sysMatSize))
+
+ for iMarker in range(1,dataCont.nMarkers):
+ xUp = dataCont.u[iMarker*dataCont.nMarkers : (iMarker+1)*dataCont.nMarkers].copy()
+ xDwn = dataCont.u[iMarker*dataCont.nMarkers : (iMarker+1)*dataCont.nMarkers].copy()
+
+ for iVar in range(dataCont.nVar):
+ xSave = xUp
+ xUp += eps
+ xDwn -= eps
+
+ # Build Jacobian for one point
+ else:
+
+ if order == 1 and linSysRes 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
+
+ JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x=xUp) - linSysRes) / eps
+
+ xUp[iVar] = varSave
+ return JacMatrix
+
+ elif order == 1:
+ raise RuntimeError('First order Jacobian determination requires the linear system residuals')
+
+ 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
+
+ JacMatrix[:,iVar] = (self._blLaws(dataCont, marker, x=xUp) - self._blLaws(dataCont, marker, x=xDwn)) / (2*eps)
+
+ xUp[iVar] = varSave
+ xDwn[iVar] = varSave
+ return JacMatrix
+
+ else:
+ raise RuntimeError('Only first and second orders Jacobian can be computed.')
+
+
+ def buildResiduals(self, dataCont, marker = -1) -> 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.
+ """
+
+ # Build Residual for the whole domain
+ if marker == -1:
+ a=1
+
+ # Build Residual for one point
+ else: return self._blLaws(dataCont, marker, x=None)
\ No newline at end of file
diff --git a/dart/viscous/model/Closures.py b/dart/viscous/model/Closures.py
new file mode 100755
index 0000000..28e4491
--- /dev/null
+++ b/dart/viscous/model/Closures.py
@@ -0,0 +1,234 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# Closures: Closure models to close the #
+# system of PDEs of the boundary layer equations #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+from math import exp
+
+class Closures:
+ def __init__(self) -> None:
+ pass
+
+ def laminarClosure(theta, H, Ret,Me, name):
+ """Laminar closure derived from the Falkner-Skan one-parameter profile family
+ Nishida and Drela (1996)
+ """
+ Hk = (H - 0.29*Me**2)/(1+0.113*Me**2) # Kinematic shape factor
+ if name == "airfoil":
+ Hk = max(Hk, 1.02)
+ Hk = min(Hk,6)
+ else:
+ Hk = max(Hk, 1.00005)
+ Hk = min(Hk,6)
+ Hstar2 = (0.064/(Hk-0.8)+0.251)*Me**2 # Density shape parameter
+ delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta)
+ Hks = Hk - 4.35
+ if Hk <= 4.35:
+ dens = Hk + 1
+ Hstar = 0.0111*Hks**2/dens - 0.0278*Hks**3/dens + 1.528 -0.0002*(Hks*Hk)**2
+ elif Hk > 4.35:
+ Hstar = 1.528 + 0.015*Hks**2/Hk
+ Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction
+ if Hk < 5.5:
+ tmp = (5.5-Hk)**3 / (Hk+1)
+ Cf = (-0.07 + 0.0727*tmp)/Ret
+ elif Hk >= 5.5:
+ tmp = 1 - 1/(Hk-4.5)
+ Cf = (-0.07 + 0.015*tmp**2) /Ret
+ if Hk < 4:
+ Cd = ( 0.00205*(4.0-Hk)**5.5 + 0.207 ) * (Hstar/(2*Ret))
+ elif Hk >= 4:
+ HkCd = (Hk-4)**2
+ denCd = 1+0.02*HkCd
+ Cd = ( -0.0016*HkCd/denCd + 0.207) * (Hstar/(2*Ret))
+ if name == 'wake':
+ Cd = 2*(1.10*(1.0-1.0/Hk)**2/Hk)* (Hstar/(2*Ret))
+ Cf = 0
+ return Cd, Cf, delta, Hstar, Hstar2, Hk
+
+
+ def turbulentClosure(theta, H, Ct, Ret, Me, name):
+ """Turbulent closure derived from Nishida and Drela (1996)
+ """
+ # eliminate absurd transients
+ Ct = min(Ct, 0.30)
+ Ct = max(Ct, 0.0000001)
+ Hk = (H - 0.29*Me**2)/(1+0.113*Me**2)
+ if name == 'wake':
+ Hk = max(Hk, 1.00005)
+ Hk = min(Hk,10)
+ else:
+ Hk = max(Hk, 1.00005)
+ Hk = min(Hk,10)
+ Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me**2
+ gam = 1.4 - 1
+ Fc = np.sqrt(1+0.5*gam*Me**2)
+ if Ret <= 400:
+ H0 = 4
+ elif Ret > 400:
+ H0 = 3 + (400/Ret)
+ if Ret > 200:
+ Ret = Ret
+ elif Ret <= 200:
+ Ret = 200
+ if Hk <= H0:
+ Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))**2)*(1.5/(Hk+0.5))+1.5+4/Ret
+ elif Hk > H0:
+ 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 = np.max((logRt,3))
+ arg = np.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)
+ if name == 'wake':
+ if Us > 0.99995:
+ Us = 0.99995
+ n = 2 # two wake halves
+ Cf = 0 # no friction inside the wake
+ Hkc = Hk - 1
+ 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
+ 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))
+ Hkc = max(Hk-1-18/Ret, 0.01)
+ # dissipation coefficient
+ Hmin = 1 + 2.1/np.log(Ret)
+ Fl = (Hk-1)/(Hmin-1)
+ Dfac = 0.5+0.5*np.tanh(Fl)
+ Cdw = 0.5*(Cf*Us) * Dfac
+ Cdd = (0.995-Us)*Ct**2
+ Cdl = 0.15*(0.995-Us)**2/Ret
+ if n*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H) >= 0:
+ Cteq = np.sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)) #Here it is Cteq^0.5
+ else:
+ Cteq = 0.03
+ Cd = n*(Cdw+ Cdd + Cdl)
+ 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):
+ '''Compute the amplitude of the TS waves envelop for the transition
+ '''
+ Dgr = 0.08
+ Hk1 = 3.5
+ Hk2 = 4
+ Hmi = (1/(Hk-1))
+ logRet_crit = 2.492*(1/(Hk-1))**0.43 + 0.7*(np.tanh(14*Hmi-9.24)+1.0) # value at which the amplification starts to grow
+ if Ret <=0:
+ Ret = 1
+ if np.log10(Ret) < logRet_crit - Dgr :
+ Ax = 0
+ else:
+ Rnorm = (np.log10(Ret) - (logRet_crit-Dgr))/(2*Dgr)
+ if Rnorm <=0:
+ Rfac = 0
+ if Rnorm >= 1:
+ Rfac = 1
+ else:
+ Rfac = 3*Rnorm**2 - 2*Rnorm**3
+ # normal envelope amplification rate
+ m = -0.05+2.7*Hmi-5.5*Hmi**2+3*Hmi**3+0.1*np.exp(-20*Hmi)
+ arg = 3.87*Hmi-2.52
+ dndRet = 0.028*(Hk-1)-0.0345*np.exp(-(arg**2))
+ Ax = (m*dndRet/theta)*Rfac
+ # set correction for separated profiles
+ if Hk > Hk1:
+ Hnorm = (Hk - Hk1)/(Hk2-Hk1)
+ if Hnorm >=1:
+ Hfac = 1
+ else:
+ Hfac = 3*Hnorm**2 - 2*Hnorm**3
+ Ax1 = Ax
+ Gro = 0.3+0.35*np.exp(-0.15*(Hk-5))
+ Tnr = np.tanh(1.2*(np.log10(Ret)-Gro))
+ Ax2 = (0.086*Tnr - 0.25/(Hk-1)**1.5)/theta
+ if Ax2 < 0:
+ Ax2 = 0
+ Ax = Hfac*Ax2 + (1-Hfac)*Ax1
+ if Ax < 0:
+ Ax = 0
+ return Ax
\ No newline at end of file
diff --git a/dart/viscous/model/DataStructures.py b/dart/viscous/model/DataStructures.py
new file mode 100755
index 0000000..6ea77a8
--- /dev/null
+++ b/dart/viscous/model/DataStructures.py
@@ -0,0 +1,217 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# Variables: Solution and external #
+# flow parameters container #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+
+class DGeometry:
+
+ def __init__(self, data, paramDict, _bNodes):
+
+ # Private
+
+ self.__globCoord = data[:,0] # Global coordinates of each marker
+ self.__locCoord = paramDict['xx'] # Local coordinates of each marker
+ self.__nMarkers = len(self.__locCoord) # Number of markers in the computational domain
+
+ self.__boundaryNodes = _bNodes # Boundary nodes of the computational domain [Stagnation point,
+ # First point on lower side,
+ # First wake point]
+
+ # Public
+
+ def GetglobCoord(self): return self.__globCoord
+ def GetlocCoord(self): return self.__locCoord
+ def GetnMarkers(self): return self.__nMarkers
+ def GetglobCoord(self): return self.__bNodes
+ def GetboundaryNodes(self): return self.__boundaryNodes
+
+
+class DSolution:
+
+ def __init__(self, _xtrF):
+
+ # Private
+
+ self.__nVar = 5 # Number of variables involved in the system of PDEs
+ self.__nBlPar = 10 # Number of boundary layer parameters
+
+ self.__xtrF = _xtrF # Forced Transition location on [upper, lower] sides
+ self.__xtr = [1, 1] # Transition location
+
+ # Public
+
+ # Solution vectors
+
+ self.u = [] # Vector containing the solution [theta, H, N, ue, Ctau]
+
+ # Boundary layer parameters
+
+ self.ReTheta = # Reynolds number based on the variable 'theta'.
+ self.disspCoeff = # Dissipation coefficient.
+ self.frictionCoeff = # Wall friction coefficient.
+ self.delta = # Boundary layer thickness.
+ self.Hstar = # Kinetic energy shape parameter.
+ self.Hstar2 = # Density shape parameter.
+ self.Hk = # Kinematic shape factor.
+ self.Cteq = # Equilibrium shear stress coefficient.
+ self.Us = # Equivalent normalized wall slip velocity.
+ self.deltaStar = # Displacement thickness of the boundary layer.
+
+ def GetnVar(self): return self.__nVar
+ def GetnBlPar(self): return self.__nBlPar
+
+ def GetxtrF(self): return self.__xtrF
+ def Getxtr(self): return self.__xtr
+ def Setxtr(self, xtrIn, side):
+ self.__xtr[side] = xtrIn
+ pass
+
+ def initializeSolution(self, DGeometry):
+
+ nMarkers = PGeometry.GetnMarkers()
+ self.u = np.empty(self.__nVar * nMarkers)
+ self.blPar = np.empty(self.__nVar * nMarkers)
+
+ return self.u, self.blPar
+
+class DOuterFlow:
+
+ def __init__(self, data) -> None:
+
+ self.__Mach =
+ self.__density =
+ self.__velocity =
+ self.__dv =
+
+ self.outerDeltaStar =
+ self.outerLocCoord =
+ self.outerRe =
+
+
+
+
+
+class DVariables:
+ """Data container for viscous solver.
+
+ Attributes
+ ----------
+
+ - Ti: float [0,1]
+ - User defined turbulence intensity that impacts the critical amplification ratio Ncrit of the e^N method.
+ - Re: float
+ - Flow Reynolds number.
+
+ - xx: numpy array
+ - Airfoil surfaces and wake mesh.
+
+ - nN: int
+ - Number of cell in the computational domain.
+
+ - bNodes: numpy array
+ - 3 elements array containing the boundary nodes of the computational domain [Stagnation point, First lower side point, First wake point].
+
+ - nVar: int
+ - Number of variables involved in the differential problem.
+
+ - u: numpy array (vector of len(nVar*nN))
+ - Array containing the solution (nVar variables at each of the nN points).
+
+ - blPar: numpy array (vector of len(nBlPar*nN))
+ - Boundary layer parameters.
+
+ - extPar: numpy array (vector of len(nExtPar*nN))
+ - External (inviscid) flow parameters.
+
+ - dv: numpy array (vector of len(nN))
+ - Velocity gradient on each cell.
+
+ - xGlobal: numpy array (vector of len(nN))
+ - Cell position non-dimensionalized by the airfoil chord: x/c
+
+ - xtrF: numpy array of len(2)
+ - User forced transiton on [Upper side, Lower side]
+
+ - xtr: numpy array of len(2)
+ - Flow laminar to turbulent transition location on [Upper side, Lower side]
+
+ - Ncrit: float
+ - Critical amplification ratio for the e^N method (default value= 9.)
+
+ """
+ def __init__(self, data, paramDict, bNodes, xtrF, Ti, Re):
+ self.__Ti=Ti
+ self.Re=Re
+ self.xx=paramDict['xx']
+ self.nN=len(self.xx)
+
+ self.bNodes=bNodes
+ # Solution
+ self.nVar=5
+ self.u=np.zeros(self.nVar*self.nN)
+
+ # Parameters of the boundary layers
+ # [ 0 , 1, 2, 3, 4, 5, 6, 7, 8, 9]
+ # [Ret, Cd, Cf, delta, Hstar, Hstar2, Hk, Cteq, Us, deltaStar]
+ self.nBlPar=10
+ self.blPar=np.zeros(self.nBlPar*self.nN)
+ self.Ret_previous=self.blPar[0::self.nVar]
+
+ self.dueOld_dx=np.zeros(self.nN)
+ self.ddeltaStarOld_dx=np.zeros(self.nN)
+
+ # Parameters of the external (inviscid) flow
+ # [ 0, 1, 2, 3, 4, 5]
+ # [Me, rhoe, vtExt, deltaStarExt, xxExt ReExt]
+ self.nExtPar=6
+ self.extPar=np.zeros(self.nExtPar*self.nN)
+ self.extPar[0::self.nExtPar]=data[:,5]
+ self.extPar[1::self.nExtPar]=data[:,6]
+ self.extPar[2::self.nExtPar]=paramDict['vtExt']
+ self.extPar[3::self.nExtPar]=data[:,7]
+ self.extPar[4::self.nExtPar]=paramDict['xx']
+ self.dv=paramDict['dv']
+ self.xGlobal=data[:,0]
+
+ self.turb=np.zeros(self.nN, dtype=int)
+
+ self.xtrF=xtrF
+ self.xtr=np.zeros(2)
+ if not len(self.xtrF)==2:
+ raise RuntimeError(ccolors.ANSI_RED + "Error in forced transition." + ccolors.ANSI_RESET)
+ for n in range(len(self.xtrF)):
+ if self.xtrF[n] is not None:
+ if 0<=self.xtrF[n]<=1:
+ self.xtr[n]=self.xtrF[n]
+ else:
+ raise ValueError(ccolors.ANSI_RED + "Non-physical transition location." + ccolors.ANSI_RESET)
+ else:
+ self.xtr[n]=1
+
+ self.idxTr=[self.bNodes[1]-1,self.bNodes[2]-1]
+
+ self.Ncrit=self.turbIntensity(Ti)
+ pass
+
+ def turbIntensity(self,Ti):
+ '''Computes the critical amplification ratio Ncrit based on the input
+ free-stream turbulence intensity. Ncrit=9 by default.'''
+ if Ti is not None and Ti<=1:
+ tau=2.7*np.tanh(self.Ti/2.7)
+ Ncrit=-8.43-2.4*np.log(tau/100) # Drela 1998
+ else:
+ Ncrit=9
+ return Ncrit
diff --git a/dart/viscous/model/Discretization.py b/dart/viscous/model/Discretization.py
new file mode 100755
index 0000000..979774b
--- /dev/null
+++ b/dart/viscous/model/Discretization.py
@@ -0,0 +1,105 @@
+################################################################################
+# #
+# Flow viscous module file #
+# Discretization: Computes solution gradient at one point #
+# of the computational domain. Different #
+# finite difference spacial schemes are implemented. #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+
+class Discretization:
+ """Discretize the spacial derivatives required to evaluate the boundary layer laws.
+
+ Attributes
+ ----------
+
+ - xx: Cell coordinates in the (local) x-direction in the boundary layer.
+
+ - xxExt: Cell coordinates in the (local) x-direction for the flow outside the boundary layer.
+
+ - bNodes: Boundary nodes of the computational domain
+
+ Returns
+ -------
+
+ Each method of the class returns the derivatives of the solution and of boundary layer parameters
+
+ - du_dx: Derivative of the solution [dTheta/dx, dH/dx, dN/dx, due/dx, dCt/dx]
+
+ - dHstar_dx: Derivative of Hstar
+
+ - dueExt_dx, ddeltaStarExt_dx: External flow derivatives. Used in the interaction law.
+
+ - c, cExt: Parameters dependent on the cell size 'dx', 'dxExt'. Used in the interation law.
+ """
+ def __init__(self, xx, xxExt, bNodes):
+ dx=np.zeros(len(xx))
+ for i in range(1,len(xx)):
+ dx[i]=xx[i]-xx[0] if i==bNodes[1] else xx[i]-xx[i-1]
+ self.dx=dx
+
+ dxExt=np.zeros(len(xxExt))
+ for i in range(1,len(xx)):
+ 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):
+ """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
+ ddeltaStarExt_dx=(var.extPar[idx*var.nExtPar+3]-var.extPar[idxPrev*var.nExtPar+3])/dxExt
+ c=4/(np.pi*dx)
+ cExt=4/(np.pi*dx)
+ return du_dx, dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt
+
+ def center(self, var, x, idxPrev2, idxPrev,idx, idxNext):
+ """Second order central difference.
+ """
+ #print(idx, idxPrev, idxNext)
+ du= (var.u[idxNext*var.nVar:(idxNext+1)*var.nVar] - var.u[idxPrev*var.nVar:(idxPrev+1)*var.nVar])
+
+ dx = var.xx[idx] - var.xx[idxPrev]
+ dxExt = var.extPar[idx*var.nExtPar + 4] - var.extPar[idxPrev*var.nExtPar + 4]
+
+ du_dx = du / (2*dx)
+ dHstar_dx = (var.blPar[idxNext*var.nBlPar+4]-var.blPar[idxPrev*var.nBlPar+4]) / (2*dx)
+ dueExt_dx = (var.extPar[idxNext*var.nExtPar+2] - var.extPar[idxPrev*var.nExtPar+2]) / (2*dxExt)
+ ddeltaStarExt_dx=(var.extPar[idxNext*var.nExtPar+3]-var.extPar[idxPrev*var.nExtPar+3])/(2*dxExt)
+ c=4/(np.pi*dx)
+ cExt=4/(np.pi*dxExt)
+ return du_dx, dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt
+
+ '''def fod(self, var, x,idxPrev2, idxPrev,idx, idxNext):
+ du=self.u[+var.nVar*idx:+var.nVar*(idx+1)]-x
+ du_dx=du/self.dx[idx]
+ dHstar_dx=(self.Hstar[idxNext]-self.Hstar[idx])/self.dx[idx]
+ return du_dx, dHstar_dx'''
+
+ # Second order backward difference
+ def sou(self, var, x,idxPrev2, idxPrev,idx, idxNext): # Attention if we want gradients at 2nd point after stagnation point on lower side.
+ """Second order upwing difference.
+ """
+ du=-3*x+4*var.u[+var.nVar*idxPrev:+var.nVar*(idxPrev+1)]-var.u[+var.nVar*idxPrev2:+var.nVar*(idxPrev2+1)]
+ du_dx=du/(2*self.dx[idx])
+ dHstar_dx=(3*var.blPar[idx*var.nBlPar+4]-4*var.blPar[idxPrev*var.nBlPar+4]+var.blPar[idxPrev2*var.nBlPar])/(2*self.dx[idx])
+ dueExt_dx=(var.extPar[idxNext*var.nExtPar+2]-var.extPar[idxPrev*var.nExtPar+2])/(2*self.dx[idx])
+ ddeltaStarExt_dx=(var.extPar[idxNext*var.nExtPar+3]-var.extPar[idxPrev*var.nExtPar+3])/(2*self.dx[idx])
+ c=4/(np.pi*self.dx[idx])
+ cExt=4/(np.pi*self.dxExt[idx])
+ return du_dx, dHstar_dx, dueExt_dx, ddeltaStarExt_dx, c, cExt
\ No newline at end of file
diff --git a/dart/viscous/model/TSolver.py b/dart/viscous/model/TSolver.py
new file mode 100644
index 0000000..cf78b7f
--- /dev/null
+++ b/dart/viscous/model/TSolver.py
@@ -0,0 +1,395 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# Time solver : Time discretization #
+# and computation towards steady state #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+from matplotlib import pyplot as plt
+from fwk.coloring import ccolors
+import numpy as np
+import fwk
+from dart.viscous.model.Closures import Closures
+
+class timeConfig:
+ """ ---------- Time integration configuration ----------
+
+ Time integration related computations: CFL adaptation, time step (per cell)
+
+ Attributes
+ ----------
+ cflAdapt: bool
+ 'True' if CFL adaptation is active, 'False' otherwise
+
+ __cflAdaptParam: numpy vect
+ [Lower, Upper] bounds of the CFL
+ """
+ def __init__(self):
+ pass
+
+ def computeTimeStep(self, CFL, dx, invVel):
+ # Speed of sound.
+ freeStreamMach = 0.15
+ gamma = 1.4
+ gradU2 = invVel**2
+ speedOfSound_Squared = 1 / (freeStreamMach * freeStreamMach) + 0.5 * (gamma - 1) * (1 - gradU2)
+ # Velocity of the fastest travelling wave
+ advVel = np.sqrt(speedOfSound_Squared)+invVel
+
+ # Time step computation
+ return CFL*dx/advVel
+
+
+ def adaptCFL(self, resCurr, res0):
+ """CFL number adaptation. The metric driving the adaptation is the residual value.
+ """
+ CFL=self.__CFL0*(res0/resCurr)**0.7
+ CFL=max(CFL,self.__CFL0)
+ return CFL
+
+ def outputState(self, var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, color):
+ if color == 'white':
+ print('| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>8.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'.format(iMarker,
+ innerIters,
+ np.log10(normLinSysRes/normLinSysRes0),
+ CFL,
+ var.xGlobal[iMarker],
+ var.flowRegime[iMarker],
+ var.u[iMarker*var.nVar + 2]))
+ elif color == 'yellow':
+ print(ccolors.ANSI_YELLOW + '| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>8.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'.format(iMarker,
+ innerIters,
+ np.log10(normLinSysRes/normLinSysRes0),
+ CFL,
+ var.xGlobal[iMarker],
+ var.flowRegime[iMarker],
+ var.u[iMarker*var.nVar + 2]),
+ ccolors.ANSI_RESET)
+ elif color == 'green':
+ print(ccolors.ANSI_YELLOW + '| {:>6.0f}| {:>11.0f}| {:>9.4f}| {:>8.2f}| {:>14.6f}| {:>6.0f}| {:>11.3f}|'.format(iMarker,
+ innerIters,
+ np.log10(normLinSysRes/normLinSysRes0),
+ CFL,
+ var.xGlobal[iMarker],
+ var.flowRegime[iMarker],
+ var.u[iMarker*var.nVar + 2]),
+ ccolors.ANSI_RESET)
+
+
+
+class implicitEuler(timeConfig):
+ """ ---------- Implicit Euler time integration ----------
+
+ Solves the unsteady boundary layer equations using the implicit Euler
+ method with Newton method to solve the non-linear system.
+
+ Attributes
+ ----------
+ solver: class
+ Spacial operator and Jacobian computation of the boundary layer laws
+
+ timeParams: dict
+ Contains necessary information for the solver: tolerance, CFL, CFL apdation parameters
+
+ safeguard: int
+ Number of iterations for which the residual has to decrease to unlock CFL adaptation
+ (usually changes after the first coupling iteration)
+ and second order Jacobian evaluation starts.
+
+ safemode: bool
+ Indicates if the solver is in safe mode (1) or not (0). Value switched according to safeguard
+ """
+
+ def __init__(self, _solver, transition, wakeBC, _cfl0):
+ timeConfig.__init__(self)
+ # Initialize time parameters
+
+ self.__CFL0=_cfl0
+ self.__maxIt = 50000
+ self.__tol=1e-6
+ self.solver=_solver
+ self.transSolver = transition
+ self.wakeBC = wakeBC
+ def __str__(self):
+ return 'Implicit Euler'
+
+ 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 timeSolve(self, var):
+ """ ---------- Newton iteration to solve the boundary layer equations ----------
+
+ Damped Newton method in Pseudo-Time. Linear algebraic system (obtained by Taylor linearis.)
+ is solved with a direct method.
+
+ Equations
+ ---------
+ dU/dt = F(U)
+ <=> (U_(n+1)-U_n)/dt = F(U_(n+1)) (Time disc.)
+ <=> (J-I)dU = -F(U_n), J=dF_n/dU_n (Linearis.)
+
+ Parameters
+ ----------
+ - var: class
+ Data container: - Solution
+ - Boundary layer parameters
+ - External flow parameters
+ - Mesh
+ - Transition information
+
+ Tasks
+ ------
+ - Updates the solution 'var.u' through Newton iterations until steady state.
+
+ - Asks 'solver' for spacial operator F and Jacobian J and uses a
+ direct linear solver (LU decomposition).
+
+ - Corrects undesirable results.
+
+ - Asks the solver for transtion position/BC and wake BC to be updated at each iteration.
+ """
+
+ # Flag to diplay the simulation evolution.
+ convergedPoints = np.zeros(var.nN)
+ convergedPoints[0] = 1 # Boundary condition
+
+ # Initialize the flow regime on each cell.
+ self.transSolver.initTransition(var)
+
+ # Transition is locked when it is found on one side. We do not look at the value of N
+ # for the remaining points on this side.
+ lockTrans = 0
+
+ for iMarker in range(1, var.nN):
+ displayState = 0
+
+ if not lockTrans and 1 < iMarker < var.bNodes[2]:
+
+ # Determine if the amplification waves start growing or not.
+ logRet_crit = 2.492*(1/(var.blPar[(iMarker-1)*var.nBlPar + 6]-1))**0.43
+ + 0.7*(np.tanh(14*(1/(var.blPar[(iMarker-1)*var.nBlPar + 6]-1))-9.24)+1.0)
+
+ if var.blPar[(iMarker-1)*var.nBlPar + 0] > 0: # Avoid log of something <= 0
+
+ if np.log10(var.blPar[(iMarker-1)*var.nBlPar + 0]) <= logRet_crit - 0.08:
+
+ var.u[iMarker*var.nVar + 2] = 0
+
+ # Upper side start.
+ if iMarker == var.bNodes[0]+1:
+ print('| - Computing upper side. {:>54}'.format('|'))
+
+ # Lower side start.
+ if iMarker == var.bNodes[1]:
+ lockTrans = 0
+ print('| - Computing lower side. {:>54}'.format('|'))
+
+ # Wake start
+ if iMarker == var.bNodes[2]: # First wake point
+
+ # Wake boundary condition.
+ self.wakeBC(var)
+ convergedPoints[iMarker] = 1
+ lockTrans = 1
+ print('| - Computing wake. {:>60}'.format('|'))
+
+ # Ignore remaining instructions for the first wake point.
+ continue
+
+ if not lockTrans and iMarker-1 < var.bNodes[2] and var.u[(iMarker-1) * var.nVar +2 ] >= var.Ncrit:
+ # Transition
+ print('| Transition located near x/c = {:<.3f}. Marker: {:<.0f} {:>29s}'.format(var.xGlobal[iMarker-1],iMarker-1, '|'))
+
+ # Save laminar solution
+ lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
+
+ # Set up turbulent portion boundary condition
+ avTurb = self.transSolver.solveTransition(var, iMarker - 1)
+
+ # Compute turbulent solution
+ var.flowRegime[iMarker - 1] = 1
+ self.newtonSolver(var, iMarker-1, displayState)
+
+ # Average solutions
+ var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
+
+ # Recompute closures
+ var.blPar[(iMarker-1)*var.nBlPar+1:(iMarker-1)*var.nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
+ var.u[(iMarker-1)*var.nVar+1],
+ var.u[(iMarker-1)*var.nVar+4],
+ var.blPar[(iMarker-1)*var.nBlPar+0],
+ var.extPar[(iMarker-1)*var.nExtPar+0],
+ 'airfoil')
+
+ lockTrans = 1
+
+ # Forced transition
+ if not lockTrans and var.xtrF[0] is not None and (var.flowRegime[iMarker]== 1 and var.flowRegime[iMarker-1] == 0):
+ print('Transition', iMarker)
+ lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
+ avTurb = self.transSolver.transitionBC(var, iMarker - 1)
+ # Compute turbulent solution
+ var.flowRegime[iMarker - 1] = 1
+ self.newtonSolver(var, iMarker-1, displayState)
+ # Average solutions
+ var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
+ # Recompute closures
+ var.blPar[(iMarker-1)*var.nBlPar+1:(iMarker-1)*var.nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
+ var.u[(iMarker-1)*var.nVar+1],
+ var.u[(iMarker-1)*var.nVar+4],
+ var.blPar[(iMarker-1)*var.nBlPar+0],
+ var.extPar[(iMarker-1)*var.nExtPar+0],
+ 'airfoil')
+
+ lockTrans = 1
+ if not lockTrans and var.xtrF[1] is not None and (var.flowRegime[iMarker]== 1 and var.flowRegime[iMarker-1] == 0):
+ print('Transition', iMarker)
+ lamSol = var.u[(iMarker-1)*var.nVar : (iMarker)* var.nVar].copy()
+ avTurb = self.transSolver.transitionBC(var, iMarker - 1)
+ # Compute turbulent solution
+ var.flowRegime[iMarker - 1] = 1
+ self.newtonSolver(var, iMarker-1, displayState)
+ # Average solutions
+ var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar] = (1-avTurb) * lamSol + avTurb * var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
+ # Recompute closures
+ var.blPar[(iMarker-1)*var.nBlPar+1:(iMarker-1)*var.nBlPar+9]=Closures.turbulentClosure(var.u[(iMarker-1)*var.nVar+0],
+ var.u[(iMarker-1)*var.nVar+1],
+ var.u[(iMarker-1)*var.nVar+4],
+ var.blPar[(iMarker-1)*var.nBlPar+0],
+ var.extPar[(iMarker-1)*var.nExtPar+0],
+ 'airfoil')
+
+ # Call Newton solver for the point
+ convergedPoints[iMarker] = self.newtonSolver(var, iMarker, displayState)
+
+ return convergedPoints
+
+
+ def newtonSolver (self, var, iMarker, displayState):
+
+ # Save initial condition in case a point diverges.
+ u0=var.u.copy()
+
+ normLinSysRes0 = np.linalg.norm(self.solver.buildResiduals(var,iMarker))
+
+ # Initial CFL and time step.
+ CFL = self.__CFL0
+ CFLAdapt = 1
+ jacEval = 5
+ dx = var.xx[iMarker] - var.xx[0] if iMarker == var.bNodes[1] else var.xx[iMarker] - var.xx[iMarker - 1]
+ dt = self.computeTimeStep(CFL, dx, var.extPar[iMarker*var.nExtPar + 2])
+ IMatrix=np.identity(var.nVar) / dt
+
+ # Main loop.
+ converged = 0
+ innerIters = 0
+ while innerIters <= self.__maxIt:
+
+ try:
+ # Residuals
+ linSysRes=self.solver.buildResiduals(var,iMarker)
+
+ # Check for convergence
+ normLinSysRes=np.linalg.norm(linSysRes)
+
+ if displayState and innerIters % 200 == 0:
+
+ self.outputState(var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, 'yellow')
+
+ if normLinSysRes <= self.__tol*normLinSysRes0:
+
+ converged = 1
+ if displayState:
+ self.outputState(var, iMarker, innerIters, normLinSysRes, normLinSysRes0, CFL, 'green')
+ break
+
+ if innerIters % jacEval == 0:
+ Jacobian=self.solver.buildJacobian(var, iMarker, linSysRes = linSysRes, eps = 1e-8)
+
+ # Linear solver
+ slnIncr = np.linalg.solve((IMatrix - Jacobian), linSysRes)
+ var.u[iMarker*var.nVar : (iMarker+1)*var.nVar] += slnIncr
+ var.u[iMarker*var.nVar + 1] = max(var.u[iMarker*var.nVar + 1], 1.0005)
+
+ except KeyboardInterrupt:
+ quit()
+ except Exception as e:
+ print('|', ccolors.ANSI_YELLOW + 'Newton solver failed at position {:<.4f}. Marker: {:<.0f} {:>30s}'.format(var.xGlobal[iMarker], iMarker,'|'), ccolors.ANSI_RESET)
+ print(e)
+ print('Current CFL', CFL)
+
+ #var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = u0[iMarker*var.nVar : (iMarker + 1)*var.nVar].copy()
+ if iMarker == var.bNodes[1]:
+ var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[0*var.nVar : 1*var.nVar]
+ elif iMarker != var.transMarkers[0] and iMarker != var.transMarkers[1]:
+ var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = var.u[(iMarker-1)*var.nVar : (iMarker)*var.nVar]
+ else:
+ var.u[iMarker*var.nVar : (iMarker + 1)*var.nVar] = u0[iMarker*var.nVar : (iMarker + 1)*var.nVar].copy()
+
+ CFL = min(max(CFL / 2,0.01),1)
+ dt = self.computeTimeStep(CFL, dx, var.extPar[iMarker*var.nExtPar + 2])
+ IMatrix = np.identity(var.nVar) / dt
+
+ CFLAdapt = 0
+ jacEval = 1
+ displayState = 1
+ print('+------------------------------------------------------------------------------+')
+ print('| {:>6s}| {:>8s}| {:>9s}| {:>8s}| {:>12s}| {:>6s}| {:>8s}|'.format('Point','Inner Iters', 'rms[F]', 'End CFL', 'Position (x/c)', 'Regime', 'Amp. factor'))
+ print('+------------------------------------------------------------------------------+')
+ print('Lowering CFL: {:<.3f}'.format(CFL))
+ continue
+
+ # Correct absurd transcient
+ var.u[iMarker*var.nVar+1] = max(var.u[iMarker*var.nVar+1],1.0005)
+
+
+ # CFL adaptation
+ if CFLAdapt:
+ if normLinSysRes <= normLinSysRes0:
+ CFL = self.__CFL0* (normLinSysRes0/normLinSysRes)**0.7
+ elif normLinSysRes > normLinSysRes0:
+ CFL = self.__CFL0* (normLinSysRes0/normLinSysRes)**0.5
+ CFL = max(CFL, 0.1)
+ dt = self.computeTimeStep(CFL, dx, var.extPar[iMarker*var.nExtPar + 2])
+ IMatrix=np.identity(var.nVar) / dt
+
+ innerIters+=1
+
+ else:
+ # Maximum number of iterations reached.
+ if normLinSysRes >= 1e-3 * normLinSysRes0:
+ print('|', ccolors.ANSI_RED + 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'.format(var.xGlobal[iMarker], np.log10(normLinSysRes/normLinSysRes0),'|'),
+ ccolors.ANSI_RESET)
+ var.u[iMarker*var.nVar:(iMarker+1)*var.nVar] = u0[(iMarker)*var.nVar: (iMarker+1)*var.nVar].copy()
+ name = 'airfoil' if iMarker <= var.bNodes[2] else 'wake'
+
+ if var.flowRegime[iMarker]==0:
+
+ # Laminar closure relations
+ var.blPar[iMarker*var.nBlPar+1:iMarker*var.nBlPar+7]=Closures.laminarClosure(var.u[iMarker*var.nVar + 0], var.u[iMarker*var.nVar + 1], var.blPar[iMarker*var.nBlPar+0], var.extPar[iMarker*var.nExtPar+0], name)
+
+ elif var.flowRegime[iMarker]==1:
+
+ # Turbulent closures relations
+ 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], name)
+
+ else:
+ print('|', ccolors.ANSI_YELLOW+ 'Too many iteration at position {:<.4f}. normLinSysRes = {:<.3f}. {:>30s}'.format(var.xGlobal[iMarker], np.log10(normLinSysRes/normLinSysRes0),'|'),
+ ccolors.ANSI_RESET)
+
+ return converged
\ No newline at end of file
diff --git a/dart/viscous/model/Transition.py b/dart/viscous/model/Transition.py
new file mode 100644
index 0000000..f82f036
--- /dev/null
+++ b/dart/viscous/model/Transition.py
@@ -0,0 +1,87 @@
+import numpy as np
+from dart.viscous.model.Closures import Closures
+
+class Transition:
+ def __init__(self, xtrF):
+ self.upperForced = xtrF[0]
+ self.lowerForced = xtrF[1]
+
+ def initTransition(self, var):
+
+ var.flowRegime=np.zeros(var.nN)
+ var.flowRegime[var.bNodes[2]:]=1
+
+ for n in range(len(var.xtrF)):
+ if var.xtrF[n] is not None:
+ if var.xtrF[n]==1: # Forced laminar flow
+ var.transMarkers[n]=var.bNodes[n+1]
+ var.xtr[n]=1
+ var.flowRegime[var.bNodes[n]:var.bNodes[n+1]]=0
+ elif var.xtrF[n]==0: # Forced turbulent flow
+ var.transMarkers[n]=var.bNodes[n]
+ var.xtr[n]=0
+ var.flowRegime[var.bNodes[n]:var.bNodes[n+1]]=1
+ elif 0<var.xtrF[n]<1:
+ idx0=(abs(var.xGlobal[n*var.bNodes[1]:var.bNodes[1+n]]-0)).argmin()
+ var.transMarkers[n]=(abs(var.xGlobal[n*var.bNodes[1]+idx0:var.bNodes[1+n]]
+ -var.xtrF[n])).argmin()+n*var.bNodes[1]+idx0
+ var.flowRegime[var.transMarkers[n]:var.bNodes[n+1]]=1
+ var.xtr[n]=var.xtrF[n]
+ else:
+ raise ValueError('Non-physical forced transition location')
+
+
+ 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
+
+
+ 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
diff --git a/dart/viscous/model/Variables.py b/dart/viscous/model/Variables.py
new file mode 100755
index 0000000..356e3b6
--- /dev/null
+++ b/dart/viscous/model/Variables.py
@@ -0,0 +1,136 @@
+################################################################################
+# #
+# DARTFlo viscous module file #
+# Variables: Solution and external #
+# flow parameters container #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+from matplotlib import pyplot as plt
+class Variables:
+ """Data container for viscous solver.
+
+ Attributes
+ ----------
+
+ - Ti: float [0,1]
+ - User defined turbulence intensity that impacts the critical amplification ratio Ncrit of the e^N method.
+ - Re: float
+ - Flow Reynolds number.
+
+ - xx: numpy array
+ - Airfoil surfaces and wake mesh.
+
+ - nN: int
+ - Number of cell in the computational domain.
+
+ - bNodes: numpy array
+ - 3 elements array containing the boundary nodes of the computational domain [Stagnation point, First lower side point, First wake point].
+
+ - nVar: int
+ - Number of variables involved in the differential problem.
+
+ - u: numpy array (vector of len(nVar*nN))
+ - Array containing the solution (nVar variables at each of the nN points).
+
+ - blPar: numpy array (vector of len(nBlPar*nN))
+ - Boundary layer parameters.
+
+ - extPar: numpy array (vector of len(nExtPar*nN))
+ - External (inviscid) flow parameters.
+
+ - dv: numpy array (vector of len(nN))
+ - Velocity gradient on each cell.
+
+ - xGlobal: numpy array (vector of len(nN))
+ - Cell position non-dimensionalized by the airfoil chord: x/c
+
+ - xtrF: numpy array of len(2)
+ - User forced transiton on [Upper side, Lower side]
+
+ - xtr: numpy array of len(2)
+ - Flow laminar to turbulent transition location on [Upper side, Lower side]
+
+ - Ncrit: float
+ - Critical amplification ratio for the e^N method (default value= 9.)
+
+ """
+ def __init__(self, data, paramDict, _xx, _xGlobal, bNodes, xtrF, Ti, Re):
+ self.__Ti=Ti
+ self.Re=Re
+ self.xx = _xx
+ self.nN=len(self.xx)
+
+ self.bNodes=bNodes
+ # Solution
+ self.nVar=5
+ self.u=np.zeros(self.nVar*self.nN)
+
+ # Parameters of the boundary layers
+ # [ 0 , 1, 2, 3, 4, 5, 6, 7, 8, 9]
+ # [Ret, Cd, Cf, delta, Hstar, Hstar2, Hk, Cteq, Us, deltaStar]
+ self.nBlPar=10
+ self.blPar=np.zeros(self.nBlPar*self.nN)
+ self.Ret_previous=self.blPar[0::self.nVar]
+
+ self.dueOld_dx=np.zeros(self.nN)
+ self.ddeltaStarOld_dx=np.zeros(self.nN)
+
+ # Parameters of the external (inviscid) flow
+ # [ 0, 1, 2, 3, 4, 5]
+ # [Me, rhoe, vtExt, deltaStarExt, xxExt ReExt]
+ self.nExtPar=6
+ self.extPar=np.zeros(self.nExtPar*self.nN)
+ self.extPar[0::self.nExtPar]=data[:,5]
+ self.extPar[1::self.nExtPar]=data[:,6]
+ self.extPar[2::self.nExtPar]=paramDict['vtExt']
+ self.extPar[3::self.nExtPar]=data[:,7]
+ self.extPar[4::self.nExtPar]=_xx
+ self.dv=paramDict['dv']
+ self.xGlobal = _xGlobal
+
+
+ """plt.plot(self.xGlobal[0:self.bNodes[1]],self.extPar[2:self.bNodes[1]*self.nExtPar:self.nExtPar])
+ plt.plot(self.xGlobal[self.bNodes[1]:self.bNodes[2]],self.extPar[self.bNodes[1]*self.nExtPar + 2 : self.bNodes[2]*self.nExtPar:self.nExtPar])
+ plt.plot(self.xGlobal[self.bNodes[2]:],self.extPar[self.bNodes[2]*self.nExtPar + 2 ::self.nExtPar])
+ plt.axvline(x=0.59194)
+ plt.show()"""
+
+ self.flowRegime = np.zeros(self.nN, dtype=int)
+
+ self.xtrF=xtrF
+ self.xtr=np.zeros(2)
+ if not len(self.xtrF)==2:
+ raise RuntimeError(ccolors.ANSI_RED + "Error in forced transition." + ccolors.ANSI_RESET)
+ for n in range(len(self.xtrF)):
+ if self.xtrF[n] is not None:
+ if 0<=self.xtrF[n]<=1:
+ self.xtr[n]=self.xtrF[n]
+ else:
+ raise ValueError(ccolors.ANSI_RED + "Non-physical transition location." + ccolors.ANSI_RESET)
+ else:
+ self.xtr[n]=1
+
+ self.transMarkers=[self.bNodes[1]-1,self.bNodes[2]-1]
+
+ self.Ncrit=self.turbIntensity(Ti)
+ pass
+
+ def turbIntensity(self,Ti):
+ '''Computes the critical amplification ratio Ncrit based on the input
+ free-stream turbulence intensity. Ncrit=9 by default.'''
+ if Ti is not None and Ti<=1:
+ tau=2.7*np.tanh(self.Ti/2.7)
+ Ncrit=-8.43-2.4*np.log(tau/100) # Drela 1998
+ else:
+ Ncrit=9
+ return Ncrit
diff --git a/dart/viscous/model/vplo.py b/dart/viscous/model/vplo.py
new file mode 100755
index 0000000..d2fdbf5
--- /dev/null
+++ b/dart/viscous/model/vplo.py
@@ -0,0 +1,87 @@
+################################################################################
+# #
+# Flow viscous module file #
+# vplo: Solution plotter #
+# and comparison with XFoil solution #
+# #
+# Author: Paul Dechamps #
+# Université de Liège #
+# Date: 2021.09.10 #
+# #
+################################################################################
+
+
+# ------------------------------------------------------------------------------
+# Imports
+# ------------------------------------------------------------------------------
+import numpy as np
+import matplotlib.pyplot as plt
+
+# ------------------------------------------------------------------------------
+# Plots the solution and compares it to xfoil file given in 'xfoilFiles' folder
+# ------------------------------------------------------------------------------
+class vplo:
+ def __init__(self, case):
+ #self.filename=case
+ self.filename=case
+
+ def showSol(self,pltIn, side, var):
+
+
+ with open('/home/parallels/labLnx/waves/flow/viscous/model/xfoilFiles/'+self.filename,'r') as f:
+ cl_xfoil=f.readline()
+ cd_xfoil=f.readline()
+ f.close()
+ # columns : x,ue,theta,Cf,H
+ data=np.loadtxt('/home/parallels/labLnx/waves/flow/viscous/model/xfoilFiles/'+self.filename,skiprows=3,usecols=(1,3,5,6,7))
+
+ # columns : x, theta, H, ue, Cf
+ data[:,[1,2,3,4]]=data[:,[2,4,1,3]]
+ data[:,]
+
+ for i in range(len(data[:,0])):
+ if data[i,0]==data[i-1,0]:
+ cut=i-1
+ np.insert(data,2,0,axis=1)
+ showXfoil=0
+
+ for i in range(len(pltIn)):
+ if pltIn[i]==1:
+ if i==3:
+ CC=-1
+ else:
+ CC=1
+ if i==0 or i==1:
+ var_xfoil=i+1
+ else:
+ var_xfoil=i
+ if side[0]==1:
+ if showXfoil==1:
+ plt.plot(data[:cut,0],data[:cut,var_xfoil],'r--',lw=0.7)
+ plt.plot(var.xGlobal[var.bNodes[0]:var.bNodes[1]],var.u[i:var.bNodes[1]*var.nVar:var.nVar],'r-',lw=1.5)
+ plt.xlim([0,1])
+ plt.xlabel('x/c')
+
+ if side[1]==1:
+ if showXfoil==1:
+ plt.plot(data[cut:,0],data[cut:,var_xfoil],'g--',lw=0.7)
+ plt.plot(np.insert(var.xGlobal[var.bNodes[1]:var.bNodes[2]],0,var.xGlobal[0]),CC*np.insert(var.u[var.bNodes[1]*var.nVar+i:var.bNodes[2]*var.nVar:var.nVar],0,var.u[i]),'g-',lw=1.5)
+ plt.xlim([0,1])
+ plt.xlabel('x/c')
+ if side[2]==1:
+ plt.plot(var.xGlobal[var.bNodes[2]:],var.u[var.bNodes[2]*var.nVar+i::var.nVar],'r-',lw=1.5)
+ plt.xlim([0,2])
+ plt.xlabel('x/c')
+ plt.xlim([0,2])
+ #plt.axvline(x=var.xtr,ymin=0,'r--',lw=1)
+ if i==0:
+ plt.title('Theta')
+ elif i==1:
+ plt.title('H')
+ elif i==2:
+ plt.title('N')
+ elif i==3:
+ plt.title('Ue')
+ elif i==4:
+ plt.title('Ct')
+ plt.show()
\ No newline at end of file
diff --git a/dart/viscous/wake.py b/dart/viscous/wake.py
old mode 100644
new mode 100755
index 91e1f32..f62982b
--- a/dart/viscous/wake.py
+++ b/dart/viscous/wake.py
@@ -19,8 +19,8 @@
## Wake behind airfoil (around which the boundary layer is computed)
# Amaury Bilocq
-from dart.viscous.boundary import Boundary
-from dart.viscous.airfoil import Airfoil
+from dart.viscous.Boundary import Boundary
+from dart.viscous.Airfoil import Airfoil
import numpy as np
class Wake(Boundary):
--
GitLab