diff --git a/flow/default.py b/flow/default.py
index 9af0ed3a25150675978222c7cc0367f0278804af..ebc7aeecf54bf2949fbc38fa207d613cee110d8d 100644
--- a/flow/default.py
+++ b/flow/default.py
@@ -65,11 +65,14 @@ def problem(msh, dim, alpha, beta, minf, mcrit, sref, lref, xref, yref, zref, su
bnd = flo.Boundary(msh, [sur, fld])
pbl.add(bnd)
if vsc:
- blw = flo.Blowing(msh, 'airfoil')
+ blw = flo.Blowing(msh, sur)
+ blww = flo.Blowing(msh, wk)
pbl.add(blw)
+ pbl.add(blww)
else:
blw = None
- return pbl, dirichlet, wake, bnd, blw
+ blww = None
+ return pbl, dirichlet, wake, bnd, [blw, blww]
## Initialize Picard solver
def picard(pbl):
diff --git a/flow/tests/bli.py b/flow/tests/bli.py
index 4a6d2262f90fded20f9d07b88d19bbc84a2904fb..3a89ec2ed433952e7ad907ffd6dc443c63311c38 100644
--- a/flow/tests/bli.py
+++ b/flow/tests/bli.py
@@ -16,9 +16,18 @@
# limitations under the License.
-# @brief Compute lifting (linear or nonlinear) viscous flow around a naca0012
-# @authors Adrien Crovato, Amaury Bilocq
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+#
+# Amaury Bilocq
# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012:
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, p = 2, m = 7, tol = 10^-4, msTE = 0.01, msLE = 0.001
+# -> nIt = 41, Cl = 0.58, Cd = 0.0062, xtrTop = 0.0555, xtrBot = 0.7397
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, p = 2, m = 7, tol = 10^-4, msTE = 0.01, msLE = 0.001
+# -> nIt = , Cl = 0.69, Cd = 0.0067, xtrTop = 0.0384, xtrBot = 0.7364
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, p = 2, m = 11, tol = 5.2*10^-3, msTE = 0.01, msLE = 0.00075
+# -> nIt = , Cl = 1.39 , Cd = 0.011, xtrTop = 0.008, xtrBot = 1
#
# CAUTION
# This test is provided to ensure that the solver works properly.
@@ -29,10 +38,11 @@ import flow.utils as floU
import flow.default as floD
import flow.viscous.solver as floVS
import flow.viscous.coupler as floVC
+import tbox
import tbox.utils as tboxU
import fwk
from fwk.testing import *
-from fwk.coloring import ccolors
+from fwk.coloring import ccolors
def main():
# timer
@@ -40,16 +50,23 @@ def main():
tms['total'].start()
# define flow variables
- alpha = 0*math.pi/180
- M_inf = 0.0
+ Re = 1e7
+ alpha = 5*math.pi/180
+ U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+ M_inf = 0.5
M_crit = 5 # Squared critical Mach number (above which density is modified)
c_ref = 1
dim = 2
+ # define filter parameters and tolerance of the VII
+ p = 2
+ m = 7
+ tol = 1e-4
+
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.0075, 'msLe' : 0.0075}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
@@ -62,8 +79,8 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
isolver = floD.newton(pbl)
- vsolver = floVS.Solver(blw)
- coupler = floVC.Coupler(isolver, vsolver, gmshWriter)
+ vsolver = floVS.Solver(Re, p, m)
+ coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
coupler.run()
tms['solver'].stop()
@@ -72,8 +89,8 @@ def main():
tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
# display results
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
- print(' M alpha Cl Cd Cm')
- print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alpha*180/math.pi, isolver.Cl, isolver.Cd, isolver.Cm))
+ print(' Re M alpha Cl Cd Cdp Cdf Cm')
+ print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
# display timers
tms['total'].stop()
@@ -83,16 +100,32 @@ def main():
# visualize solution and plot results
floD.initViewer(pbl)
- tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, isolver.Cd, isolver.Cm, 4), True)
+ tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- if M_inf == 0 and alpha == 0*math.pi/180:
- tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.41, 1e-1)) # TODO check value and tolerance
- tests.add(CTest('Cl', isolver.Cl, 0., 5e-2))
+ if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.58, 5e-2))
+ tests.add(CTest('Cd', vsolver.Cd, 0.0062, 0.01))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 0.01))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.056, 0.05))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.05))
+ elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 0.69, 5e-2))
+ tests.add(CTest('Cd', vsolver.Cd, 0.0067, 0.01))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 0.01))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.038, 0.05))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.05))
+ elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180:
+ tests.add(CTest('Cl', isolver.Cl, 1.39, 5e-2))
+ tests.add(CTest('Cd', vsolver.Cd, 0.0011, 0.01))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 0.01))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.008, 0.05))
+ tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.05))
else:
raise Exception('Test not defined for this flow')
+
tests.run()
# eof
diff --git a/flow/viscous/airfoil.py b/flow/viscous/airfoil.py
new file mode 100644
index 0000000000000000000000000000000000000000..0dd0ccd147eda95c7f3952eae6078e1afe5df29b
--- /dev/null
+++ b/flow/viscous/airfoil.py
@@ -0,0 +1,128 @@
+#!/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 flow.viscous.boundary 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
+ # Defining TE/LE nodes number
+ idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
+ globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
+ idxTE = np.where(data[:,1] == 1.0)[0]
+ if idxTE[0] < idxTE[1]:
+ upperIdxTE = idxTE[0]
+ lowerIdxTE = idxTE[1]
+ else:
+ upperIdxTE = idxTE[1]
+ lowerIdxTE = idxTE[0]
+ upperGlobTE = data[upperIdxTE,0] # Number of the upper TE node
+ lowerGlobTE = data[lowerIdxTE,0] # Number of the lower TE node
+ connectListElems[0] = np.where(eData[:,1] == globStag)[0]
+ N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes
+ 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] == int(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] == int(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/flow/viscous/boundary.py b/flow/viscous/boundary.py
new file mode 100644
index 0000000000000000000000000000000000000000..3a8e862bb038441ce3b6d8398dc22a13f4247a59
--- /dev/null
+++ b/flow/viscous/boundary.py
@@ -0,0 +1,37 @@
+#!/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/flow/viscous/coupler.py b/flow/viscous/coupler.py
index a75855a54c97586bee7c76e7ea5e3e94825944a2..727a9a73ebd9693fa132fda9d0cb81c4741da85d 100644
--- a/flow/viscous/coupler.py
+++ b/flow/viscous/coupler.py
@@ -16,18 +16,21 @@
# limitations under the License.
-# @brief Viscous-inviscid coupler
-# @authors Adrien Crovato, Amaury Bilocq
+## Viscous-inviscid coupler (quasi-simultaneous coupling)
+#
+# Amaury Bilocq
import numpy as np
from fwk.coloring import ccolors
+from flow.viscous.airfoil import Airfoil
+from flow.viscous.wake import Wake
class Coupler:
- def __init__(self, _isolver, _vsolver, _writer):
- '''...
- '''
+ 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):
@@ -35,28 +38,48 @@ class Coupler:
'''
# initialize loop
it = 0
- converged = False # temp
- while True:
+ 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()
- # get velocity at edge of BL from inviscid solver and update viscous solver
- for i in range(0, len(self.vsolver.boundary.nodes)):
- self.vsolver.v[i,0] = self.isolver.U[self.vsolver.boundary.nodes[i].row][0]
- self.vsolver.v[i,1] = self.isolver.U[self.vsolver.boundary.nodes[i].row][1]
- self.vsolver.v[i,2] = self.isolver.U[self.vsolver.boundary.nodes[i].row][2]
- # run viscous solver
- self.vsolver.run()
- # get blowing velocity from viscous solver and update inviscid solver
- for i in range(0, self.vsolver.nE):
- self.vsolver.boundary.setU(i, self.vsolver.u[i])
- # convergence test
- if converged:
- break
+ print('--- Viscous solver parameters ---')
+ print('Filter windows length:', self.vsolver.m)
+ print('Filter polynomial order:', self.vsolver.p)
+ print('Tolerance:', self.tol)
+ 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:
- converged = True # perform 2 iterations
- it += 1
-
+ CdOld = self.vsolver.Cd
+ it += 1
+ self.vsolver.it += 1
# save results
self.isolver.save(0, self.writer)
+ self.vsolver.writeFile()
print('\n')
diff --git a/flow/viscous/newton.py b/flow/viscous/newton.py
new file mode 100644
index 0000000000000000000000000000000000000000..bed39ec64ca0a43df345cce476da8c4c5e8da6c2
--- /dev/null
+++ b/flow/viscous/newton.py
@@ -0,0 +1,51 @@
+#!/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.
+
+
+## Newton raphson method coupled with linear solver
+#
+# Amaury Bilocq
+
+import numpy as np
+from fwk.coloring import ccolors
+
+def newton(f, x, maxIt, tol=1.0e-9):
+ '''Newton procedure to solve the non linear set of boundary layer equations.
+ Boundary layer equations are parabolic in x -> downstream marching is used
+ An implicit marching model is applied (Crank Nicolson) -> matrix inversion is performed
+ '''
+ # Compute Jacobian
+ def jacobian(f, x):
+ dx = 1.0e-8
+ n = len(x)
+ jac = np.zeros((n,n))
+ f0 = f(x)
+ for j in range(n):
+ Dxj = (abs(x[j])*dx if x[j] != 0 else dx)
+ x_plus = [(xi if k != j else xi+Dxj) for k, xi in enumerate(x)]
+ jac[:,j] = (f(x_plus)-f0)/Dxj # numerical solution of the jacobian
+ return jac, f0
+ # Solve
+ for i in range(maxIt):
+ jac, f0 = jacobian(f,x)
+ if np.sqrt(np.dot(f0,f0)/len(x)) < tol and x[0]>0 and x[1]>0: return x
+ dx = np.linalg.solve(jac, -f0)
+ x = x + dx
+ x[1] = max(x[1],1.0005) # to avoid too low shape parameter. Solver must be redone to be more robust
+ if np.sqrt(np.dot(dx, dx)) < tol*max(abs(any(x)),1) and x[0]>0 and x[1]>0: return x
+ # Error
+ raise RuntimeError(ccolors.ANSI_RED + 'Viscous Newton solver - too many iterations, aborting!\n' + ccolors.ANSI_RESET)
diff --git a/flow/viscous/solver.py b/flow/viscous/solver.py
index 3b52be11afd0d7e6a9f90bea6c8f29338e5bdec6..6af0f757cd87ee6b61dcdc004de6fea9c2c6942f 100644
--- a/flow/viscous/solver.py
+++ b/flow/viscous/solver.py
@@ -16,35 +16,474 @@
# limitations under the License.
-# @brief Boundary Layer Equations viscous solver
-# @authors Adrien Crovato, Amaury Bilocq
+## Integral boundary layer equations viscous solver
+#
+# Amaury Bilocq
+#
+# TODO list
+# - Improve the filter in boundaryLayer (l 263)
+# - Inverse procedure at the singularity location (such as in Xfoil)
+# - Create a new mesh for the viscous solver
+# - Improve the transition solver by refining the grid to split the transitional element into a laminar element and a turbulent element.
+# Here it is an average between laminar and turbulent solution (l 397-416)
+# - Implement quasi-2D method for 3D capability
+# - User can define the transition location (l 303)
+# - Write the empirical relation between the amplification factor and the turbulence intensity (l 303)
+# - Check why the computation crash if the first point of the wake is computed with the turbulent closure (l 309)
+from flow.viscous.boundary import Boundary
+from flow.viscous.wake import Wake
+import flow.viscous.newton as Newton
import numpy as np
+from scipy.signal import savgol_filter
-class Solver:
- def __init__(self, _boundary):
- '''...
+class Solver:
+ def __init__(self, _Re, _p, _m):
+ '''Initialize the viscous solver
'''
- # TODO generalize to handle several boundaries (use python array)
- self.boundary = _boundary # boundary on which BLE are sovled
- 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
-
- def run(self):
- ''' Solve BLE
+ self.Re = _Re # Reynolds number
+ self.p = _p # Filtering order
+ self.m = _m # Filtering bandwith
+ self.gamma = 1.4 # heat capacity of air
+ self.Cd = 0 # drag coefficient
+ self.it = 0 # viscous inviscid iteration for interaction law
+ self.xtr = [1, 1] # set transition at TE reference coordinates
+
+ 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['xtrTop'] = self.xtr[0]
+ wData['xtrBot'] = self.xtr[1]
+ return wData
+
+ def writeFile(self):
+ '''Write the results in airfoil_viscous.dat
+ '''
+ wData = self.dictionary()
+ toW = ['delta*', 'H', 'Cf'] # 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 __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 __amplRoutine(self, Hk, Ret, theta):
+ '''Compute the amplitude of the TS waves envelop for the transition
'''
- print('--- I am a fake BL solver and I am setting dummy blowing velocities! ---')
- # 6th order fitting of delta_star for naca 0012 at alpha = 0, from XFoil
- a = 0.000054925976379640116
- b = 0.007407077078697043
- c = -0.043918316829113194
- d = 0.1487463222137225
- e = -0.20162146328539832
- f = 0.0880698637668166
- g = 0.004436812366681563
- for i in range(0, self.nE):
- x = 0.5*(self.boundary.tag.elems[i].nodes[0].pos[0]+self.boundary.tag.elems[i].nodes[1].pos[0])
- # CAUTION: a positive velocity will result in a suction!
- self.u[i] = -(b + 2*c*x + 3*d*x*x* + 4*e*x*x*x + 5*f*x*x*x*x + 6*g*x*x*x*x*x) # ue = d(delta_star) / dx
+ 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
+
+ def __laminarClosure(self, 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, Hstar, Hstar2, Hk, delta, H
+
+ def __turbulentClosure(self, 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 = self.gamma - 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
+ Cteq = np.sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)) #Here it is Cteq^0.5
+ Cd = n*(Cdw+ Cdd + Cdl)
+ Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2)
+ return Cd, Cf, Hstar, Hstar2, Hk, Cteq, delta, Ueq, Ct
+
+
+ def __boundaryLayer(self, data, group):
+ '''Discretize and solve the boundary layer equations for laminar/turbulent flows
+ '''
+ xx, dv, vtOld, nN, alpha = self.__getBLcoordinates(data)
+ Me = data[:,5]
+ rhoe = data[:,6]
+ deltaStarOld = data[:,7]
+ #Filter the initial conditions (Me and rho can be filtered too)
+ if group.name == 'airfoil':
+ vtOld = savgol_filter(vtOld, self.m, self.p, mode = 'interp')
+ vtInv = vtOld
+ rhoInv = rhoe
+ if self.it == 0:
+ xxOld = xx
+ else:
+ xxOld = data[:,8]
+ #Initialize variables
+ blwVel = np.zeros(nN-1) # Blowing velocity
+ vt = np.zeros(nN) # Boundary layer velocity
+ theta = np.zeros(nN) # Momentum thickness
+ H = np.zeros(nN) # Shape factor
+ deltaStar = np.zeros(nN) # Displacement thickness
+ Hstar = np.zeros(nN) # Kinetic energy shape parameter
+ Hstar2 = np.zeros(nN) # Density shape parameter
+ Hk = np.zeros(nN) # Kinematic shape factor
+ Ret = np.zeros(nN) #Re based on momentum thickness
+ Cd = np.zeros(nN) # Dissipation coefficient
+ Cf = np.zeros(nN) # Skin friction
+ n = np.zeros(nN) # Logarithm of the maximum amplification ratio
+ Ax = np.zeros(nN) # Amplification factor
+ Ct = np.zeros(nN) # Shear stress coefficient
+ Cteq = np.zeros(nN) # Equilibrium shear stress coefficient
+ Ueq = np.zeros(nN) # Equilibrium velocity gradient
+ delta = np.zeros(nN) # BL thickness
+ # Boundary conditions
+ if group.name == 'airfoil':
+ theta[0], H[0], n[0], Ct[0] = group.initialConditions(self.Re, dv[0])
+ elif group.name == 'wake':
+ theta[0] = group.T0
+ H[0] = group.H0
+ n[0] = group.n0
+ Ct[0] = group.Ct0
+ vt[0] = vtOld[0]
+ Ret[0] = vt[0]*theta[0]*self.Re
+ if Ret[0] < 1:
+ Ret[0] = 1
+ vt[0] = Ret[0]/(theta[0]*self.Re)
+ deltaStar[0] = theta[0]*H[0]
+ # Define transition location ( N = 9) and compute stagnation point
+ ntr = 9
+ if n[0] < 9:
+ turb =0 # we begin with laminar flow
+ Cd[0], Cf[0], Hstar[0], Hstar2[0], Hk[0], delta[0], H[0] = self.__laminarClosure( theta[0], H[0], Ret[0],Me[0], group.name)
+ else :
+ turb = 1 # typically for the wake
+ Cd[0], Cf[0], Hstar[0], Hstar2[0], Hk[0], delta[0], H[0] = self.__laminarClosure( theta[0], H[0], Ret[0],Me[0], group.name)
+ xtr = 0 # if only turbulent flow
+ itTurb = 0
+ # Solve the boundary layer equations
+ for i in range(1,nN):
+ # x[0] = theta; x[1] = H; x[2] = n for laminar/Ct for turbulent; x[3] = vt from interaction law
+ def f_lam(x):
+ f = np.zeros(len(x))
+ Ret[i] = np.maximum(x[3]*x[0]*self.Re, 1)
+ Cd[i], Cf[i], Hstar[i], Hstar2[i], Hk[i], delta[i], x[1] = self.__laminarClosure(x[0], x[1], Ret[i],Me[i], group.name)
+ Ta = upw*x[0]+(1-upw)*theta[i-1]
+ Ha = upw*x[1]+(1-upw)*H[i-1]
+ uea = upw*x[3]+(1-upw)*vt[i-1]
+ Reta = upw*Ret[i] + (1-upw)*Ret[i-1]
+ Mea = upw*Me[i]+(1-upw)*Me[i-1]
+ Cda, Cfa, Hstara, Hstar2a, Hka, deltaa, Ha = self.__laminarClosure(Ta , Ha, Reta, Mea, group.name)
+ dx = xx[i] - xx[i-1]
+ Axa = self.__amplRoutine(Hka, Reta, Ta)
+ dTheta = x[0]-theta[i-1]
+ due = x[3]-vt[i-1]
+ dH = x[1]-H[i-1]
+ dn = x[2] - n[i-1]
+ dHstar = Hstar[i]-Hstar[i-1]
+ Hstara = upw*Hstar[i]+(1-upw)*Hstar[i-1]
+ f[0] = dTheta/dx + (2 + Ha - Mea**2)*(Ta/uea)*due/dx - Cfa/2
+ f[1] = Ta*(dHstar/dx) + (2*Hstar2a + Hstara*(1-Ha))*Ta/uea * due/dx - 2*Cda + Hstara*Cfa/2
+ f[2] = dn - dx*Axa
+ f[3] = uea - 4/(np.pi*dx)*Ta*Ha - (upw*vtOld[i]+(1-upw)*vtOld[i-1]) + 4/(np.pi*(xxOld[i]-xxOld[i-1]))*(upw*deltaStarOld[i]+(1-upw)*deltaStarOld[i-1])
+ return f
+ def f_turb(x):
+ f = np.zeros(len(x))
+ if group.name == 'wake':
+ dissipRatio = 0.9 # wall/wake dissipation length ratio
+ else:
+ dissipRatio = 1
+ Ret[i] = x[3]*x[0]*self.Re
+ Ta = upw*x[0]+(1-upw)*theta[i-1]
+ xxa = upw*xx[i]+(1-upw)*xx[i-1]
+ Ha = upw*x[1]+(1-upw)*H[i-1]
+ Cta = upw*x[2]+(1-upw)*Ct[i-1]
+ uea = upw*x[3]+(1-upw)*vt[i-1]
+ Reta = np.maximum(upw*(x[3]*x[0]*self.Re)+(1-upw)*(vt[i-1]*theta[i-1]*self.Re), 1)
+ Mea = upw*Me[i]+(1-upw)*Me[i-1]
+ Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i], Ueq[i], x[2] = self.__turbulentClosure(x[0], x[1],x[2], Ret[i],Me[i], group.name)
+ Cda, Cfa, Hstara, Hstar2a, Hka,Cteqa, deltaa, Ueqa, Cta = self.__turbulentClosure(Ta, Ha,Cta, Reta,Mea, group.name)
+ dx = xx[i]-xx[i-1]
+ dTheta = x[0]-theta[i-1]
+ due = x[3]-vt[i-1]
+ dH = x[1]-H[i-1]
+ dHstar = Hstar[i]-Hstar[i-1]
+ dCt = x[2]-Ct[i-1] # here it is Ct^1/2 which is represented
+ Ta = upw*x[0]+(1-upw)*theta[i-1]
+ Cta = upw*x[2]+(1-upw)*Ct[i-1]
+ f[0] = dTheta/dx + (2 + Ha - Mea**2)*(Ta/uea)*due/dx - Cfa/2
+ f[1] = Ta*(dHstar/dx) + (2*Hstar2a + Hstara*(1-Ha))*Ta/uea * due/dx - 2*Cda + Hstara*Cfa/2
+ f[2] = (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)
+ f[3] = uea - 4/(np.pi*dx)*Ta*Ha - (upw*vtOld[i]+(1-upw)*vtOld[i-1]) + 4/(np.pi*(xxOld[i]-xxOld[i-1]))*(upw*deltaStarOld[i]+(1-upw)*deltaStarOld[i-1])
+ return f
+ # Laminar solver
+ if turb == 0:
+ if n[i-1] > 8 or i < 5:
+ upw = 1 # backward Euler
+ else:
+ upw = 0.5 # Trapezoidal scheme
+ x = np.array([theta[i-1],H[i-1], n[i-1], vt[i-1]])
+ sol = Newton.newton(f_lam, x, 200)
+ # Determine if the amplification waves start to grow or not
+ logRet_crit = 2.492*(1/(Hk[i]-1))**0.43 + 0.7*(np.tanh(14*(1/(Hk[i]-1))-9.24)+1.0)
+ if np.log10(Ret[i]) <= logRet_crit - 0.08 :
+ n[i] = 0
+ else:
+ n[i] = sol[2]
+ xtr = 1
+ # Turbulent solver
+ elif turb == 1:
+ if i <2 and group.name =='wake':
+ upw = 1 # begin of the wake
+ elif itTurb <2 and group.name == 'airfoil':
+ upw = 1
+ itTurb += 1
+ else:
+ upw = 0.5 # trapezoidal scheme
+ x = np.array([theta[i-1], H[i-1], Ct[i-1], vt[i-1]])
+ sol = Newton.newton(f_turb, x, 200)
+ Ct[i] = sol[2]
+ theta[i] = sol[0]
+ H[i] = sol[1]
+ vt[i] = sol[3]
+ # Transition solver
+ if n[i] >= ntr and turb ==0:
+ turb = 1
+ upw = 1
+ xtr = (ntr-n[i-1])*((data[i,0]-data[i-1,0])/(n[i]-n[i-1]))+data[i-1,0]
+ avTurb = (data[i,0]-xtr)/(data[i,0]-data[i-1,0]) # percentage of the subsection corresponding to a turbulent flow
+ avLam = 1- avTurb # percentage of the subsection corresponding to a laminar flow
+ # Averaging both solutions
+ Cteq[i-1]= self.__turbulentClosure(theta[i-1], H[i-1], 0.03, Ret[i-1], Me[i-1], group.name)[5]
+ Ct[i-1] = avLam*Cteq[i-1] * 0.7
+ x_turb = np.array([theta[i-1], 1.515, Ct[i-1] , vt[i-1]])
+ sol_turb = Newton.newton(f_turb, x_turb, 200)
+ theta[i] = avLam*sol[0]+avTurb*sol_turb[0]
+ H[i] = avLam*sol[1]+avTurb*sol_turb[1]
+ if group.name == 'airfoil':
+ Ct[i] = max(avTurb*sol_turb[2],0.03)
+ else:
+ Ct[i] = max(avTurb*sol_turb[2],0.05)
+ vt[i] = avLam*sol[3]+avTurb*sol_turb[3]
+ 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], Ueq[i], _ = self.__turbulentClosure(theta[i], H[i],Ct[i], Ret[i],Me[i], group.name)
+ deltaStar[i] = H[i]*theta[i]
+ blwVel[i-1] = -1*(rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1])) # TO DO: need to divide by rhoe or not?
+ # Normalize the friction and dissipation coefficients by the freestream velocity
+ Cf = vt**2*Cf
+ Cd = vt**2*Cd
+ # Compute the various drag coefficients (total, friction, profile)
+ Cdtot = 2*theta[-1]*(vt[-1]**((H[-1]+5)/2))
+ Cdf = np.trapz(Cf, data[:,0])
+ Cdp = Cdtot-Cdf
+ # Outputs
+ blScalars = [Cdtot, Cdf, Cdp]
+ blVectors = [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct]
+ return blwVel, xtr, xx, blScalars, blVectors
+
+
+ def run(self, group):
+ ''' Run the viscous solver to solve the BL equations and give the outputs to the coupler
+ '''
+ group.connectListNodes, group.connectListElems, data = group.connectList()
+ # Compute BLE for the airfoil (upper section(data[0]) and lower section (data[1]))
+ if len(data) == 2:
+ ublwVel, xtrTop, uxx, ublScalars, ublVectors = self.__boundaryLayer(data[0], group)
+ lblwVel, xtrBot, lxx, lblScalars, lblVectors = self.__boundaryLayer(data[1], group)
+ # Put the upper and lower values together
+ self.data = np.concatenate((data[0], data[1]))
+ self.blScalars = np.add(ublScalars, lblScalars)
+ self.blVectors = np.hstack((ublVectors, lblVectors))
+ blwVel = np.concatenate((ublwVel, lblwVel))
+ blwVel = savgol_filter(blwVel, self.m, self.p, mode = 'interp')
+ self.xtr = [xtrTop, xtrBot]
+ # Boundary conditions for the wake
+ group.TEnd = [ublVectors[1][-1], lblVectors[1][-1]]
+ group.HEnd = [ublVectors[2][-1], lblVectors[2][-1]]
+ group.CtEnd = [ublVectors[8][-1], lblVectors[8][-1]]
+ # Delete the stagnation point doublon for variables in the VII loop
+ lblVectors[0] = np.delete(lblVectors[0],0,0) #deltastar
+ lxx = np.delete(lxx,0,0) # airfoil coordinates
+ group.deltaStar = np.concatenate((ublVectors[0], lblVectors[0]))
+ group.xx = np.concatenate((uxx, lxx))
+ # Compute BLE for the wake
+ else:
+ blwVel, _, group.xx, blScalars, blVectors = self.__boundaryLayer(data, group)
+ group.deltaStar = blVectors[0]
+ # The drag is measured at the end of the wake
+ self.Cd = blScalars[0]
+ self.Cdp = self.Cd-self.blScalars[1] # Cdf comes from the airfoil as there is no friction in the wake
+ self.blScalars[0] = self.Cd
+ self.blScalars[2] = self.Cdp
+ self.Cdf = self.blScalars[1]
+ # Sort the following results in reference frame
+ group.deltaStar = group.deltaStar[group.connectListNodes.argsort()]
+ group.xx = group.xx[group.connectListNodes.argsort()]
+ group.u = blwVel[group.connectListElems]
+
+
+
+
+
+
+
+
+
+
+
diff --git a/flow/viscous/wake.py b/flow/viscous/wake.py
new file mode 100644
index 0000000000000000000000000000000000000000..0e8da24ab423f1f4d34bb89a34780bb27095aced
--- /dev/null
+++ b/flow/viscous/wake.py
@@ -0,0 +1,79 @@
+#!/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 flow.viscous.boundary import Boundary
+from flow.viscous.airfoil 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] == 1.0)[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