From fd99999c7f38063c472d0623fb6f4171b34b4e15 Mon Sep 17 00:00:00 2001
From: AmauryBilocq <amaury.bilocq@gmail.com>
Date: Sun, 26 Apr 2020 14:47:01 +0200
Subject: [PATCH] V1.0 Viscous Inviscid interaction with direct solver
---
flow/tests/bliLiftComp.py | 134 ++++++++++
flow/tests/bliLiftIncomp.py | 134 ++++++++++
flow/tests/{bli.py => bliNonLift.py} | 27 +-
flow/viscous/Airfoil.py | 125 ++++++++++
flow/viscous/Boundary.py | 37 +++
flow/viscous/Wake.py | 75 ++++++
flow/viscous/coupler.py | 60 +++--
flow/viscous/newton.py | 69 +++++
flow/viscous/solver.py | 359 +++++++++++++++++++++++++--
9 files changed, 967 insertions(+), 53 deletions(-)
create mode 100644 flow/tests/bliLiftComp.py
create mode 100644 flow/tests/bliLiftIncomp.py
rename flow/tests/{bli.py => bliNonLift.py} (83%)
create mode 100644 flow/viscous/Airfoil.py
create mode 100644 flow/viscous/Boundary.py
create mode 100644 flow/viscous/Wake.py
create mode 100644 flow/viscous/newton.py
diff --git a/flow/tests/bliLiftComp.py b/flow/tests/bliLiftComp.py
new file mode 100644
index 00000000..b924e101
--- /dev/null
+++ b/flow/tests/bliLiftComp.py
@@ -0,0 +1,134 @@
+#!/usr/bin/env python
+# -*- 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.
+
+
+# @brief Compute lifting (linear or nonlinear) viscous flow around a naca0012
+# @authors Adrien Crovato, Amaury Bilocq
+# Test the viscous-inviscid interaction scheme
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+from __future__ import print_function
+import math
+import flow as flo
+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 tbox.gmsh as gmsh
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+def main():
+ # timer
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ # define flow variables
+ Re = 10*10**6
+ alpha = 2*math.pi/180
+ U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+ M_inf = 0.6
+ M_crit = 2 # Squared critical Mach number (above which density is modified)
+ c_ref = 1
+ dim = len(U_inf)
+
+ # mesh the geometry
+ print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+ tms['msh'].start()
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
+ msh = gmsh.MeshLoader("../models/n0012.geo",__file__).execute(**pars)
+ # crack the mesh
+ mshCrck = tbox.MshCrack(msh, dim)
+ mshCrck.setCrack('wake')
+ mshCrck.addBoundaries(['field', 'airfoil', 'downstream'])
+ mshCrck.run()
+ del mshCrck
+ gmshWriter = tbox.GmshExport(msh)
+ gmshWriter.save(msh.name)
+ tms['msh'].stop()
+
+ # set the problem and add medium, initial/boundary conditions
+ tms['pre'].start()
+ pbl = flo.Problem(msh, dim, alpha, M_inf, c_ref, c_ref, 0.25, 0., 0.)
+ if M_inf == 0:
+ pbl.set(flo.Medium(msh, 'field', flo.F0ElRhoL(), flo.F0ElMachL(), flo.F0ElCpL(), flo.F0PsPhiInf(dim, alpha)))
+ else:
+ pbl.set(flo.Medium(msh, 'field', flo.F0ElRho(M_inf, M_crit), flo.F0ElMach(M_inf), flo.F0ElCp(M_inf), flo.F0PsPhiInf(dim, alpha)))
+ pbl.add(flo.Initial(msh, 'field', flo.F0PsPhiInf(dim, alpha)))
+ pbl.add(flo.Dirichlet(msh, 'upstream', flo.F0PsPhiInf(dim, alpha)))
+ pbl.add(flo.Freestream(msh, 'side', tbox.Fct1C(-U_inf[0], -U_inf[1], 0.)))
+ pbl.add(flo.Freestream(msh, 'downstream', tbox.Fct1C(-U_inf[0], -U_inf[1], 0.)))
+ pbl.add(flo.Wake(msh, ['wake', 'wake_', 'field']))
+ pbl.add(flo.Kutta(msh, ['te', 'wake_', 'airfoil', 'field']))
+ bnd = flo.Boundary(msh, ['airfoil', 'field'])
+ pbl.add(bnd)
+ blw = flo.Blowing(msh, ['airfoil', 'field'])
+ blw2 = flo.Blowing(msh, ['wake', 'field'])
+ pbl.add(blw)
+ pbl.add(blw2)
+ tms['pre'].stop()
+
+ # solve inviscid problem
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+ isolver = flo.Newton(pbl)
+ floD.initNewton(isolver)
+ vsolver = floVS.Solver(Re)
+ coupler = floVC.Coupler(isolver, vsolver, gmshWriter,blw,blw2)
+ coupler.run()
+ tms['solver'].stop()
+
+ # extract Cp
+ Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+ tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+ # display results
+ print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+ print(' 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, vsolver.Cd, isolver.Cm))
+
+ # display timers
+ tms['total'].stop()
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+
+ # 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, vsolver.Cd, isolver.Cm, 4), True)
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.0558, 1e-1)) # TODO check value and tolerance
+ tests.add(CTest('Cl', isolver.Cl, 0.292, 5e-2))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0013, 0.01))
+ tests.add(CTest('Cd', vsolver.Cd, 0.0057, 0.01))
+ tests.add(CTest('top_xtr', vsolver.top_xtr, 0.16, 0.5 ))
+ tests.add(CTest('bot_xtr', vsolver.bot_xtr, 0.48, 0.5 ))
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/flow/tests/bliLiftIncomp.py b/flow/tests/bliLiftIncomp.py
new file mode 100644
index 00000000..635cc62a
--- /dev/null
+++ b/flow/tests/bliLiftIncomp.py
@@ -0,0 +1,134 @@
+#!/usr/bin/env python
+# -*- 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.
+
+
+# @brief Compute lifting (linear or nonlinear) viscous flow around a naca0012
+# @authors Adrien Crovato, Amaury Bilocq
+# Test the viscous-inviscid interaction scheme
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+from __future__ import print_function
+import math
+import flow as flo
+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 tbox.gmsh as gmsh
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+def main():
+ # timer
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ # define flow variables
+ Re = 10*10**6
+ alpha = 5*math.pi/180
+ U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+ M_inf = 0.0
+ M_crit = 5 # Squared critical Mach number (above which density is modified)
+ c_ref = 1
+ dim = len(U_inf)
+
+ # mesh the geometry
+ print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+ tms['msh'].start()
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
+ msh = gmsh.MeshLoader("../models/n0012.geo",__file__).execute(**pars)
+ # crack the mesh
+ mshCrck = tbox.MshCrack(msh, dim)
+ mshCrck.setCrack('wake')
+ mshCrck.addBoundaries(['field', 'airfoil', 'downstream'])
+ mshCrck.run()
+ del mshCrck
+ gmshWriter = tbox.GmshExport(msh)
+ gmshWriter.save(msh.name)
+ tms['msh'].stop()
+
+ # set the problem and add medium, initial/boundary conditions
+ tms['pre'].start()
+ pbl = flo.Problem(msh, dim, alpha, M_inf, c_ref, c_ref, 0.25, 0., 0.)
+ if M_inf == 0:
+ pbl.set(flo.Medium(msh, 'field', flo.F0ElRhoL(), flo.F0ElMachL(), flo.F0ElCpL(), flo.F0PsPhiInf(dim, alpha)))
+ else:
+ pbl.set(flo.Medium(msh, 'field', flo.F0ElRho(M_inf, M_crit), flo.F0ElMach(M_inf), flo.F0ElCp(M_inf), flo.F0PsPhiInf(dim, alpha)))
+ pbl.add(flo.Initial(msh, 'field', flo.F0PsPhiInf(dim, alpha)))
+ pbl.add(flo.Dirichlet(msh, 'upstream', flo.F0PsPhiInf(dim, alpha)))
+ pbl.add(flo.Freestream(msh, 'side', tbox.Fct1C(-U_inf[0], -U_inf[1], 0.)))
+ pbl.add(flo.Freestream(msh, 'downstream', tbox.Fct1C(-U_inf[0], -U_inf[1], 0.)))
+ pbl.add(flo.Wake(msh, ['wake', 'wake_', 'field']))
+ pbl.add(flo.Kutta(msh, ['te', 'wake_', 'airfoil', 'field']))
+ bnd = flo.Boundary(msh, ['airfoil', 'field'])
+ pbl.add(bnd)
+ blw = flo.Blowing(msh, ['airfoil', 'field'])
+ blw2 = flo.Blowing(msh, ['wake', 'field'])
+ pbl.add(blw)
+ pbl.add(blw2)
+ tms['pre'].stop()
+
+ # solve inviscid problem
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+ isolver = flo.Newton(pbl)
+ floD.initNewton(isolver)
+ vsolver = floVS.Solver(Re)
+ coupler = floVC.Coupler(isolver, vsolver, gmshWriter,blw,blw2)
+ coupler.run()
+ tms['solver'].stop()
+
+ # extract Cp
+ Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+ tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+ # display results
+ print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+ print(' 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, vsolver.Cd, isolver.Cm))
+
+ # display timers
+ tms['total'].stop()
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+
+ # 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, vsolver.Cd, isolver.Cm, 4), True)
+
+ # check results
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.9531, 1e-1)) # TODO check value and tolerance
+ tests.add(CTest('Cl', isolver.Cl, 0.5658, 5e-2))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0017, 0.01))
+ tests.add(CTest('Cd', vsolver.Cd, 0.0061, 0.01))
+ tests.add(CTest('top_xtr', vsolver.top_xtr, 0.0531, 0.5 ))
+ tests.add(CTest('bot_xtr', vsolver.bot_xtr, 0.7484, 0.5 ))
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/flow/tests/bli.py b/flow/tests/bliNonLift.py
similarity index 83%
rename from flow/tests/bli.py
rename to flow/tests/bliNonLift.py
index 04c3d967..d1f4ef78 100644
--- a/flow/tests/bli.py
+++ b/flow/tests/bliNonLift.py
@@ -36,7 +36,7 @@ import tbox.utils as tboxU
import tbox.gmsh as gmsh
import fwk
from fwk.testing import *
-from fwk.coloring import ccolors
+from fwk.coloring import ccolors
def main():
# timer
@@ -44,6 +44,7 @@ def main():
tms['total'].start()
# define flow variables
+ Re = 5*10**6
alpha = 0*math.pi/180
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
M_inf = 0.0
@@ -54,7 +55,7 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.0075, 'msLe' : 0.0075}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
msh = gmsh.MeshLoader("../models/n0012.geo",__file__).execute(**pars)
# crack the mesh
mshCrck = tbox.MshCrack(msh, dim)
@@ -81,9 +82,10 @@ def main():
pbl.add(flo.Kutta(msh, ['te', 'wake_', 'airfoil', 'field']))
bnd = flo.Boundary(msh, ['airfoil', 'field'])
pbl.add(bnd)
- # TODO also add Blowing on wake...
blw = flo.Blowing(msh, ['airfoil', 'field'])
+ blw2 = flo.Blowing(msh, ['wake', 'field'])
pbl.add(blw)
+ pbl.add(blw2)
tms['pre'].stop()
# solve inviscid problem
@@ -91,8 +93,8 @@ def main():
tms['solver'].start()
isolver = flo.Newton(pbl)
floD.initNewton(isolver)
- vsolver = floVS.Solver(blw)
- coupler = floVC.Coupler(isolver, vsolver, gmshWriter)
+ vsolver = floVS.Solver(Re)
+ coupler = floVC.Coupler(isolver, vsolver, gmshWriter,blw,blw2)
coupler.run()
tms['solver'].stop()
@@ -102,7 +104,7 @@ def main():
# 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('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, isolver.Cm))
# display timers
tms['total'].stop()
@@ -112,16 +114,17 @@ 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))
- else:
- raise Exception('Test not defined for this flow')
+ tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.4132, 1e-1)) # TODO check value and tolerance
+ tests.add(CTest('Cl', isolver.Cl, 0.0, 5e-2))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.00067, 0.01))
+ tests.add(CTest('Cd', vsolver.Cd, 0.0051, 0.01))
+ tests.add(CTest('top_xtr', vsolver.top_xtr, 0.4381, 0.5 ))
+ tests.add(CTest('bot_xtr', vsolver.bot_xtr, 0.4368, 0.5 ))
tests.run()
# eof
diff --git a/flow/viscous/Airfoil.py b/flow/viscous/Airfoil.py
new file mode 100644
index 00000000..3c954dab
--- /dev/null
+++ b/flow/viscous/Airfoil.py
@@ -0,0 +1,125 @@
+#!/usr/bin/env python
+# -*- 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.
+
+
+# @brief Boundary: airfoil
+# @authors 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"
+ 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.H0 = 2.23
+ self.n0 = 0
+ self.Ct0 = 0.03
+ 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(), 9))
+ i = 0
+ for n in self.boundary.nodes:
+ data[i,0] = n.no
+ data[i,1] = n.pos.x[0]
+ data[i,2] = n.pos.x[1]
+ data[i,3] = n.pos.x[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]
+ 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.groups[0].tag.elems[i].no
+ eData[i,1] = self.boundary.groups[0].tag.elems[i].nodes[0].no
+ eData[i,2] = self.boundary.groups[0].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] = self.deltaStar
+ # 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)
+ return connectListElems,[uData, lData]
\ No newline at end of file
diff --git a/flow/viscous/Boundary.py b/flow/viscous/Boundary.py
new file mode 100644
index 00000000..cd9a1c20
--- /dev/null
+++ b/flow/viscous/Boundary.py
@@ -0,0 +1,37 @@
+#!/usr/bin/env python
+# -*- 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.
+
+
+# @brief Mother class of all the physical boundaries
+# @authors Amaury Bilocq
+
+import numpy as np
+
+
+class Boundary(object):
+ def __init__(self, _boundary):
+ self.boundary = _boundary # boundary which is considerated
+ self.nN = len(self.boundary.nodes) # number of nodes
+ self.nE = len(self.boundary.groups[0].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.connectListElems = np.zeros(self.nE) # Connectivity for the blowing velocity
+ self.Cd = 0 # drag coefficient of the physical boundary
+ self.Cdp = 0 # Form drag
+ self.deltaStar = np.zeros(self.nN) # Displacement thickness
\ No newline at end of file
diff --git a/flow/viscous/Wake.py b/flow/viscous/Wake.py
new file mode 100644
index 00000000..667230bd
--- /dev/null
+++ b/flow/viscous/Wake.py
@@ -0,0 +1,75 @@
+#!/usr/bin/env python
+# -*- 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.
+
+
+# @brief Boundary: airfoil
+# @authors 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"
+ self.T0 = 0 # initial condition for the momentum thickness
+ self.H0 = 0
+ self.n0 = 0 # should not look at transition inside the wake
+ 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(), 9))
+ i = 0
+ for n in self.boundary.nodes:
+ data[i,0] = n.no
+ data[i,1] = n.pos.x[0]
+ data[i,2] = n.pos.x[1]
+ data[i,3] = n.pos.x[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]
+ 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.groups[0].tag.elems[i].no
+ eData[i,1] = self.boundary.groups[0].tag.elems[i].nodes[0].no
+ eData[i,2] = self.boundary.groups[0].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] = self.deltaStar
+ # Separated upper and lower part
+ data = np.delete(data,0,1)
+ return connectListElems, data
diff --git a/flow/viscous/coupler.py b/flow/viscous/coupler.py
index 66c10e6e..723f509c 100644
--- a/flow/viscous/coupler.py
+++ b/flow/viscous/coupler.py
@@ -24,42 +24,64 @@ from builtins import range
from builtins import object
import numpy as np
from fwk.coloring import ccolors
+from flow.viscous.Airfoil import Airfoil
+from flow.viscous.Wake import Wake
class Coupler(object):
- def __init__(self, _isolver, _vsolver, _writer):
+ def __init__(self, _isolver, _vsolver, _writer, _boundaryAirfoil, _boundaryWake):
'''...
'''
self.isolver =_isolver # inviscid solver
self.vsolver = _vsolver # viscous solver
self.writer = _writer
+ self.airfoil = Airfoil(_boundaryAirfoil) # create the object airfoil
+ self.wake = Wake(_boundaryWake) # create the object wake
+ self.group = [self.airfoil, self.wake]
def run(self):
''' Perform coupling
'''
# initialize loop
it = 0
- converged = False # temp
- while True:
+ alpha = 1# underrelaxation factor
+ converged = 0 # temp
+ tol = 10**-6
+ CdOld = self.vsolver.Cd
+ # while it < 15:
+ 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].x[0]
- self.vsolver.v[i,1] = self.isolver.U[self.vsolver.boundary.nodes[i].row].x[1]
- self.vsolver.v[i,2] = self.isolver.U[self.vsolver.boundary.nodes[i].row].x[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
+ for n in range(0, len(self.group)):
+ uOld = self.group[n].u
+ 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].x[0]
+ self.group[n].v[i,1] = self.isolver.U[self.group[n].boundary.nodes[i].row].x[1]
+ self.group[n].v[i,2] = self.isolver.U[self.group[n].boundary.nodes[i].row].x[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])
+ self.group[n].u = self.group[n].u[self.group[n].connectListElems]
+ self.group[n].u = uOld + alpha*(self.group[n].u - uOld) # underrelaxation
+ 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
+ self.group[n+1].n0 = self.group[n].nEnd[0]
+ # 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('Computation is done for: ' ,self.group[n].name)
+ print('Number of iterations:' ,+it)
+ print('Cd = ' ,+self.vsolver.Cd)
+ print('Top_xtr = ' ,+self.vsolver.top_xtr)
+ print('Bot_xtr = ' ,+self.vsolver.bot_xtr)
+ if abs(self.vsolver.Cd - CdOld) < tol:
+ converged = 1
else:
- converged = True # perform 2 iterations
- it += 1
-
+ CdOld = self.vsolver.Cd
+ it += 1
# save results
self.isolver.save(0, self.writer)
print('\n')
diff --git a/flow/viscous/newton.py b/flow/viscous/newton.py
new file mode 100644
index 00000000..262818ea
--- /dev/null
+++ b/flow/viscous/newton.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+# -*- 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.
+
+
+# @brief Newton raphson method coupled with linear solver using numpy library
+# @authors Amaury Bilocq
+
+from __future__ import print_function
+from builtins import range
+from builtins import object
+import numpy as np
+import math
+
+
+def newtonRaphson(f, x, maxIt, tol=1.0e-9):
+ ''' Newton procedure to solve the coupled, non linear system. Boundary layer equation are parabolic in x --> Use of downstream marching type of solver
+ As an implicit marching model is applied (Crank Nicolson), a matrix inversion is performed'''
+ # Trial to develop 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
+
+ # Backtrack linesearch
+ def BacktrackingLineSearch(f, x, dx, maxLsIt):
+ a = 1
+ f0 = f(x+dx)
+ fevalIt = 1
+ while fevalIt < maxLsIt:
+ fa = f(x+a*dx)
+ fevalIt += 1
+ print('fa = ', + fa ,'f0 = ', +f0)
+ if abs(fa[1]) > abs(f0[1]): # compare H
+ a = a/2
+ print('a = ' ,+a)
+ else:
+ break
+ return a
+
+ 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
+ # print('New x =',+ x)
+ if np.sqrt(np.dot(dx, dx)) < tol*max(abs(any(x)),1) and x[0]>0 and x[1]>0: return x
+ print("Too many iterations")
+
+
diff --git a/flow/viscous/solver.py b/flow/viscous/solver.py
index af289ab3..a467f217 100644
--- a/flow/viscous/solver.py
+++ b/flow/viscous/solver.py
@@ -22,32 +22,347 @@
from __future__ import print_function
from builtins import range
from builtins import object
+from flow.viscous.Boundary import Boundary
+from flow.viscous.Wake import Wake
+import flow.viscous.newton as newton
import numpy as np
+import math
+import matplotlib.pyplot as plt
+import os.path
+from scipy.io import loadmat
+from scipy.signal import savgol_filter
class Solver(object):
- def __init__(self, _boundary):
+ def __init__(self, _Re):
'''...
'''
- # 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.groups[0].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):
+ self.Re = _Re # Reynolds number
+ self.gamma = 1.4
+ self.Cd = 0
+
+ def run(self, group):
''' Solve BLE
'''
- 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.groups[0].tag.elems[i].nodes[0].pos.x[0]+self.boundary.groups[0].tag.elems[i].nodes[1].pos.x[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
+ def writeFile(data, H, theta, Cf, vt, Cdissip, Ct, Dstar, blwUe):
+ # save_path = 'D:/Amaury/Documents/TFE/xfoil/'
+ name = 'BLparameters'
+ name2 = 'blwU'
+ # complete_name = os.path.join(save_path, name+".dat")
+ # complete_name2 = os.path.join(save_path, name2+".dat")
+ f = open(name, 'w+')
+ for i in range(len(H)):
+ f.write("%s %s %s %s %s %s %s %s\n" % (data[i,0], H[i], theta[i], Cf[i], vt[i], Cdissip[i], Ct[i], Dstar[i]))
+ f.close()
+ f = open(name2, 'w+')
+ for i in range(len(blwUe)):
+ f.write("%s %s\n" % (0, blwUe[i]))
+ f.close()
+
+ def __getBLcoordinates(data):
+ 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] = abs(data[i,3]*np.cos(alpha[i]) + data[i,4]*np.sin(alpha[i])) # tangential speed at node 1 element i
+ vt[i+1] = abs(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, dx, dv, vt, nN
+
+ def boundaryLayer(data):
+ # Get the BL coordinates
+ xx, dx, dv, vt, nN = __getBLcoordinates(data)
+ if group.name == 'airfoil':
+ vt = savgol_filter(vt, 11, 2, mode = 'interp')
+ Me = data[:,5]
+ rhoe = data[:,6]
+ deltaStarOld = data[:,7]
+ #Initialize variable
+ blwVel = np.zeros(nN-1) # Blowing 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
+ dndRet = np.zeros(nN) # Variation of n according to Re_theta
+ m = np.zeros(nN) # Empirical relation
+ l = np.zeros(nN) # Idem
+ Ct = np.zeros(nN) # Shear stress coefficient
+ Cteq = np.zeros(nN) # Equilibrium shear stress coefficient
+ delta = np.zeros(nN) # BL thickness
+ # Initial condition
+ 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
+ idxTr = 0 # To Do: it is the index of the transition, should find a way to delete it for the wake
+ # vt[0] = vt[0]+ 4/(np.pi*dx[i-1])*(theta[0]*H[0]-deltaStarOld[i]) !!!! remains vt[0]?
+ Ret[0] = vt[0]*theta[0]*self.Re
+ deltaStar[0] = theta[0]*H[0]
+ # print('---------------Grid point n 0---------------------')
+ Cd[0], Cf[0], Hstar[0], Hstar2[0], dndRet[0], m[0], l[0], Hk[0], delta[0] = newLaminarClosure( theta[0], H[0], Ret[0],Me[0], group.name)
+ ntr = 9 # amplification ratio at which transition occurs TO DO: params from user (if not, impose a value of 9)
+ if n[0] < 9:
+ turb =0 # we begin with laminar flow
+ else :
+ turb =1 # typically for the wake
+ xtr = 0 # if only turbulent flow
+ # print('------------The solution at grid point ',+0,'=', + theta[0], +H[0], +Ct[0],+vt[0], 'X position:', + data[0,0])
+ 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
+ '''
+ # print('---------------Grid point n',+i,'---------------------')
+ def f_lam(x):
+ f = np.zeros(len(x))
+ Ret[i] = x[3]*x[0]*self.Re
+ Cd[i], Cf[i], Hstar[i], Hstar2[i], dndRet[i], m[i], l[i], Hk[i], delta[i] = newLaminarClosure(x[0], x[1], Ret[i],Me[i], group.name)
+ dTheta = x[0]-theta[i-1]
+ due = x[3]-vt[i-1]
+ dH = x[1]-H[i-1]
+ dHstar = Hstar[i]-Hstar[i-1]
+ Ta = upw*x[0]+(1-upw)*theta[i-1]
+ RhsTheta1 = (2 + x[1] - Me[i]**2) * (x[0]/x[3]) * (due/dx[i-1]) - Cf[i]/2
+ RhsTheta2 = (2 + H[i-1] - Me[i-1]**2) * (theta[i-1]/vt[i-1]) * (due/dx[i-1])- Cf[i-1]/2
+ RhsH1 = (2*Hstar2[i] + Hstar[i]*(1-x[1])) * (x[0]/x[3]) * (due/dx[i-1]) - 2*Cd[i] + 0.5*Cf[i]*Hstar[i]
+ RhsH2 = (2*Hstar2[i-1] + Hstar[i-1]*(1-H[i-1])) * (theta[i-1]/vt[i-1]) * (due/dx[i-1]) - 2*Cd[i-1] + 0.5*Cf[i-1]*Hstar[i-1]
+ RhsN1 = dndRet[i]*((m[i])+1)/2 * l[i]/x[0]
+ RhsN2 = dndRet[i-1]*((m[i-1])+1)/2 * l[i-1]/theta[i-1]
+ f[0] = dTheta/dx[i-1] + ((1-upw)*RhsTheta2+upw*RhsTheta1)
+ f[1] = Ta*(dHstar/dx[i-1]) + ((1-upw)*RhsH2+upw*RhsH1)
+ f[2] = (x[2]-n[i-1])/dx[i-1] - ((1-upw)*RhsN2+upw*RhsN1)
+ # f[3] = x[3]-vt[i]- (4/(np.pi*dx[i-1]))*(x[0]*x[1]-deltaStarOld[i])
+ f[3] = x[3] -vt[i]
+ return f
+ # Central differencing to gain accuracy
+ def f_turb(x):
+ f = np.zeros(len(x))
+ Ret[i] = x[3]*x[0]*self.Re
+ Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i] = newTurbulentClosure(x[0], x[1],x[2], Ret[i],Me[i], group.name)
+ 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]
+ deriv = 1/Cta * dCt/dx[i-1]
+ RhsTheta1 = (2 + x[1] - Me[i]**2) * (x[0]/x[3]) * (due/dx[i-1]) - Cf[i]/2
+ RhsTheta2 = (2 + H[i-1] - Me[i-1]**2) * (theta[i-1]/vt[i-1]) * (due/dx[i-1])- Cf[i-1]/2
+ RhsH1 = (2*Hstar2[i] + Hstar[i]*(1-x[1])) * (x[0]/x[3]) * (due/dx[i-1]) - 2*Cd[i] + 0.5*Cf[i]*Hstar[i]
+ RhsH2 = (2*Hstar2[i-1] + Hstar[i-1]*(1-H[i-1])) * (theta[i-1]/vt[i-1]) * (due/dx[i-1]) - 2*Cd[i-1] + 0.5*Cf[i-1]*Hstar[i-1]
+ delta1 =delta[i]
+ delta2 = delta[i-1]
+ TmpHk1 = 4/(3*x[1]*x[0]) * (0.5*Cf[i] - ((Hk[i]-1)/(6.7*Hk[i]))**2)
+ TmpHk2 = 4/(3*H[i-1]*theta[i-1]) * (0.5*Cf[i-1] - ((Hk[i-1]-1)/(6.7*Hk[i-1]))**2)
+ TmpUe1 = 1/x[3] * due/dx[i-1]
+ TmpUe2 = 1/vt[i-1] * due/dx[i-1]
+ TmpCteq1 = 5.6*(Cteq[i]-x[2])
+ TmpCteq2 = 5.6*(Cteq[i-1]-Ct[i-1])
+ f[0] = dTheta/dx[i-1] + ((1-upw)*RhsTheta2+upw*RhsTheta1)
+ f[1] = Ta*(dHstar/dx[i-1]) + ((1-upw)*RhsH2+upw*RhsH1)
+ f[2] = 2*((upw*delta1+(1-upw)*delta2)*(deriv - (upw*TmpHk1+(1-upw)*TmpHk2) + (upw*TmpUe1+(1-upw)*TmpUe2))) - (upw*TmpCteq1+(1-upw)*TmpCteq2)
+ # f[3] = x[3]-vt[i]- 4/(np.pi*dx[i-1])*(x[0]*x[1]-deltaStarOld[i])
+ f[3] = x[3] -vt[i]
+ return f
+ '''Solver depending on the case '''
+ # Laminar solver
+ if turb == 0:
+ if n[i-1] > 8.5:
+ upw = 0.5 # Trapezoidal scheme
+ elif i < 7:
+ upw = 1 # backward Euler
+ else:
+ upw = 0.5 # Trapezoidal scheme
+ x = np.array([theta[i-1],H[i-1], n[i-1], vt[i]])
+ sol = newton.newtonRaphson(f_lam, x, 120)
+ logRet_crit = 2.492*(1/(H[i-1]-1))**0.43 + 0.7*(np.tanh(14*(1/(H[i-1]-1))-9.24)+1.0) # value at which the amplification starts to growth
+ if np.log10(Ret[i]) <= logRet_crit - 0.08 :
+ n[i] = 0
+ else:
+ n[i] = sol[2]
+ xtr = 1 # no transition if remains laminar
+ # Turbulent solver
+ elif turb == 1:
+ upw = 0.5 # trapezoidal scheme
+ x = np.array([theta[i-1], H[i-1], Ct[i-1], vt[i]])
+ sol = newton.newtonRaphson(f_turb, x, 120)
+ Ct[i] = sol[2]
+ # print('------------The solution at grid point ',+i,'=', + sol,+vt[i], 'X position:', + data[i,0])
+ theta[i] = sol[0]
+ H[i] = sol[1]
+ vt[i] = sol[3]
+ # Transition solver
+ if n[i] >= ntr and turb ==0:
+ # Initial condition
+ turb=1
+ Ct[i-1] = Ct[0]
+ # Interpolation to define position of the transition
+ idxTr = i
+ 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 laminar flow
+ avLam = 1- avTurb # percentage of the subsection corresponding to a turbulent flow
+ # Compute boundary layer parameters with FOU for turbulent case
+ upw = 0
+ x_turbTr = np.array([theta[i-1], 1.515, Ct[i-1], vt[i]])
+ sol_turbTr = newton.newtonRaphson(f_turb, x_turbTr, 30)
+ # # Averaging both solutions
+ theta[i] = avLam*sol[0]+avTurb*sol_turbTr[0]
+ H[i] = avLam*sol[1]+avTurb*sol_turbTr[1]
+ vt[i] = avLam*sol[3]+avTurb*sol_turbTr[3]
+ # theta[i] = avTurb*sol_turbTr[0]
+ # H[i] = avTurb*sol_turbTr[1]
+ Ct[i] = avTurb*sol_turbTr[2]
+ # print('Ok for the transition part')
+ # print('------------The solution at grid point ',+i,'=', +theta[i], +H[i], +Ct[i],+vt[i], 'X position:', + data[i,0])
+ deltaStar[i] = H[i]*theta[i]
+ blwVel[i-1] = (rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*dx[i-1]) # TO DO: need to divide by rhoe or not?
+ Cdf = np.trapz(Cf, data[:,0])
+ Cdtot = 2*theta[-1]*(vt[-1]**((H[-1]+5)/2))
+ blVariables = [blwVel, deltaStar, theta, H, Cf, xtr, Cdf, Cdtot, vt, Cd ,n[idxTr], Ct]
+ return blVariables
+
+ def newLaminarClosure(theta, H, Ret,Me, name):
+ ''' Laminar closure derived from the Falkner-Skan one-parameter profile family
+ Drela(1996)'''
+ Hk = (H - 0.29*Me**2)/(1+0.113*Me**2) # Kinematic shape factor
+ Hstar2 = (0.064/(Hk-0.8)+0.251)*Me**2 # Density shape parameter
+ dndRet = 0.01*np.sqrt((2.4*Hk-3.7+2.5*np.tanh(1.5*Hk-4.65))**2+0.25)
+ l = (6.54*Hk-14.07)/(Hk**2)
+ m = (0.058*(((Hk-4)**2)/(Hk-1))-0.068)/l
+ Hmi = (1/(H-1))
+ # m = -0.05+2.7*Hmi-5.5*Hmi**2+3*Hmi**3+0.1*np.exp(-20*Hmi)
+ # dndRet = 0.028*(Hk-1)-0.0345*np.exp((-3.87*Hmi+2.52)**2)
+ 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:
+ print('Laminar separation')
+ 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))
+ # print(' Ret, Cd, Cf, Hk, H* = ',+ Ret, + Cd, +Cf, +Hk, + Hstar)
+ return Cd, Cf, Hstar, Hstar2, dndRet, m, l, Hk, delta
+
+ def newTurbulentClosure(theta, H, Ct, Ret, Me, name):
+ ''' Turbulent closure derived using the skin friction and velocity profile formulas of Swafford (1983)'''
+ Hk = (H - 0.29*Me**2)/(1+0.113*Me**2)
+ if name == "airfoil":
+ Hk = max(Hk, 1.05)
+ else:
+ Hk = max(Hk, 1.00005)
+ 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:
+ print('The flow is separated here')
+ Hstar = (Hk-H0)**2*(0.007*np.log10(Ret)/(Hk-H0+4/np.log10(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.log10(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
+ Cdi = 0.03*0.75**3*Us**3 # inner shear layer Row 1062 from xfoil
+ else:
+ if Us > 0.95:
+ Us = 0.98
+ n = 1
+ # Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-1))
+ Cf0 = 0.01013/(logRt-1.02) - 0.00075
+ Ho = 1/(1-6.55*np.sqrt(Cf0/2))
+ Cf = (1/Fc)*Cf0*(-0.5 + 0.9/(Hk/Ho-0.4))
+ Hkc = Hk - 1 - 18/Ret
+ Hkc = max(Hkc, 0.01)
+ Cdi = 0
+ Cd = n*((Cf*Us/2)+(0.995-Us)*Ct**2 +0.15*(0.995-Us)**2/Ret)
+ Cteq = np.sqrt(Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*((Hk**2)*H))) #Here it is Cteq^0.5
+ # print('Cteq, Cd, Cf, Hk, H* = ', +Cteq, + Cd, +Cf, +Hk, + Hstar)
+ return Cd, Cf, Hstar, Hstar2, Hk, Cteq, delta
+
+ group.connectListElems, data = group.connectList()
+ if len(data) == 2:
+ ublVariables = boundaryLayer(data[0])
+ lblVariables = boundaryLayer(data[1])
+ blwVel = abs(np.concatenate((ublVariables[0], lblVariables[0])))
+ group.Cd = ublVariables[7] + lblVariables[7]
+ group.Cdp = group.Cd-(ublVariables[6]+lblVariables[6])
+ self.top_xtr = ublVariables[5]
+ self.bot_xtr = lblVariables[5]
+ self.Cdp = group.Cdp
+ group.TEnd = [ublVariables[2][-1], lblVariables[2][-1]]
+ group.HEnd = [ublVariables[3][-1], lblVariables[3][-1]]
+ group.CtEnd = [ublVariables[11][-1], lblVariables[11][-1]]
+ group.nEnd = [ublVariables[10], lblVariables[10]]
+ # Write the results
+ writeFile(np.concatenate((data[0], data[1])), np.concatenate((ublVariables[3], lblVariables[3])), np.concatenate((ublVariables[2], lblVariables[2])),
+ np.concatenate((ublVariables[4], lblVariables[4])), np.concatenate((ublVariables[8], lblVariables[8])),
+ np.concatenate((ublVariables[9], lblVariables[9])), np.concatenate((ublVariables[11], lblVariables[11])),
+ np.concatenate((ublVariables[1], lblVariables[1])), blwVel)
+ lblVariables[1] = np.delete(lblVariables[1],0,0) # issue with the number of point
+ group.deltaStar = np.concatenate((ublVariables[1], lblVariables[1]))
+ else:
+ blVariables = boundaryLayer(data)
+ blwVel = abs(blVariables[0])
+ group.Cd = blVariables[7]
+ group.Cdp = group.Cd-blVariables[6]
+ self.Cd = group.Cd
+ group.deltaStar = blVariables[1]
+
+ group.u = blwVel
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
--
GitLab