Skip to content
Snippets Groups Projects
Commit fd99999c authored by AmauryBilocq's avatar AmauryBilocq
Browse files

V1.0 Viscous Inviscid interaction with direct solver

parent bee73258
No related branches found
No related tags found
No related merge requests found
#!/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()
#!/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()
...@@ -36,7 +36,7 @@ import tbox.utils as tboxU ...@@ -36,7 +36,7 @@ import tbox.utils as tboxU
import tbox.gmsh as gmsh import tbox.gmsh as gmsh
import fwk import fwk
from fwk.testing import * from fwk.testing import *
from fwk.coloring import ccolors from fwk.coloring import ccolors
def main(): def main():
# timer # timer
...@@ -44,6 +44,7 @@ def main(): ...@@ -44,6 +44,7 @@ def main():
tms['total'].start() tms['total'].start()
# define flow variables # define flow variables
Re = 5*10**6
alpha = 0*math.pi/180 alpha = 0*math.pi/180
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1 U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
M_inf = 0.0 M_inf = 0.0
...@@ -54,7 +55,7 @@ def main(): ...@@ -54,7 +55,7 @@ def main():
# mesh the geometry # mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start() 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) msh = gmsh.MeshLoader("../models/n0012.geo",__file__).execute(**pars)
# crack the mesh # crack the mesh
mshCrck = tbox.MshCrack(msh, dim) mshCrck = tbox.MshCrack(msh, dim)
...@@ -81,9 +82,10 @@ def main(): ...@@ -81,9 +82,10 @@ def main():
pbl.add(flo.Kutta(msh, ['te', 'wake_', 'airfoil', 'field'])) pbl.add(flo.Kutta(msh, ['te', 'wake_', 'airfoil', 'field']))
bnd = flo.Boundary(msh, ['airfoil', 'field']) bnd = flo.Boundary(msh, ['airfoil', 'field'])
pbl.add(bnd) pbl.add(bnd)
# TODO also add Blowing on wake...
blw = flo.Blowing(msh, ['airfoil', 'field']) blw = flo.Blowing(msh, ['airfoil', 'field'])
blw2 = flo.Blowing(msh, ['wake', 'field'])
pbl.add(blw) pbl.add(blw)
pbl.add(blw2)
tms['pre'].stop() tms['pre'].stop()
# solve inviscid problem # solve inviscid problem
...@@ -91,8 +93,8 @@ def main(): ...@@ -91,8 +93,8 @@ def main():
tms['solver'].start() tms['solver'].start()
isolver = flo.Newton(pbl) isolver = flo.Newton(pbl)
floD.initNewton(isolver) floD.initNewton(isolver)
vsolver = floVS.Solver(blw) vsolver = floVS.Solver(Re)
coupler = floVC.Coupler(isolver, vsolver, gmshWriter) coupler = floVC.Coupler(isolver, vsolver, gmshWriter,blw,blw2)
coupler.run() coupler.run()
tms['solver'].stop() tms['solver'].stop()
...@@ -102,7 +104,7 @@ def main(): ...@@ -102,7 +104,7 @@ def main():
# display results # display results
print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
print(' M alpha Cl Cd Cm') 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 # display timers
tms['total'].stop() tms['total'].stop()
...@@ -112,16 +114,17 @@ def main(): ...@@ -112,16 +114,17 @@ def main():
# visualize solution and plot results # visualize solution and plot results
floD.initViewer(pbl) 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 # check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests() tests = CTests()
if M_inf == 0 and alpha == 0*math.pi/180: tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.4132, 1e-1)) # TODO check value and tolerance
tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.41, 1e-1)) # TODO check value and tolerance tests.add(CTest('Cl', isolver.Cl, 0.0, 5e-2))
tests.add(CTest('Cl', isolver.Cl, 0., 5e-2)) tests.add(CTest('Cdp', vsolver.Cdp, 0.00067, 0.01))
else: tests.add(CTest('Cd', vsolver.Cd, 0.0051, 0.01))
raise Exception('Test not defined for this flow') 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() tests.run()
# eof # eof
......
#!/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
#!/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
#!/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
...@@ -24,42 +24,64 @@ from builtins import range ...@@ -24,42 +24,64 @@ from builtins import range
from builtins import object from builtins import object
import numpy as np import numpy as np
from fwk.coloring import ccolors from fwk.coloring import ccolors
from flow.viscous.Airfoil import Airfoil
from flow.viscous.Wake import Wake
class Coupler(object): class Coupler(object):
def __init__(self, _isolver, _vsolver, _writer): def __init__(self, _isolver, _vsolver, _writer, _boundaryAirfoil, _boundaryWake):
'''... '''...
''' '''
self.isolver =_isolver # inviscid solver self.isolver =_isolver # inviscid solver
self.vsolver = _vsolver # viscous solver self.vsolver = _vsolver # viscous solver
self.writer = _writer 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): def run(self):
''' Perform coupling ''' Perform coupling
''' '''
# initialize loop # initialize loop
it = 0 it = 0
converged = False # temp alpha = 1# underrelaxation factor
while True: 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) print(ccolors.ANSI_BLUE + 'Iteration: ', it, ccolors.ANSI_RESET)
# run inviscid solver # run inviscid solver
self.isolver.run() self.isolver.run()
# get velocity at edge of BL from inviscid solver and update viscous solver for n in range(0, len(self.group)):
for i in range(0, len(self.vsolver.boundary.nodes)): uOld = self.group[n].u
self.vsolver.v[i,0] = self.isolver.U[self.vsolver.boundary.nodes[i].row].x[0] for i in range(0, len(self.group[n].boundary.nodes)):
self.vsolver.v[i,1] = self.isolver.U[self.vsolver.boundary.nodes[i].row].x[1] self.group[n].v[i,0] = self.isolver.U[self.group[n].boundary.nodes[i].row].x[0]
self.vsolver.v[i,2] = self.isolver.U[self.vsolver.boundary.nodes[i].row].x[2] self.group[n].v[i,1] = self.isolver.U[self.group[n].boundary.nodes[i].row].x[1]
# run viscous solver self.group[n].v[i,2] = self.isolver.U[self.group[n].boundary.nodes[i].row].x[2]
self.vsolver.run() self.group[n].Me[i] = self.isolver.M[self.group[n].boundary.nodes[i].row]
# get blowing velocity from viscous solver and update inviscid solver self.group[n].rhoe[i] = self.isolver.rho[self.group[n].boundary.nodes[i].row]
for i in range(0, self.vsolver.nE): # run viscous solver
self.vsolver.boundary.setU(i, self.vsolver.u[i]) self.vsolver.run(self.group[n])
# convergence test self.group[n].u = self.group[n].u[self.group[n].connectListElems]
if converged: self.group[n].u = uOld + alpha*(self.group[n].u - uOld) # underrelaxation
break 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: else:
converged = True # perform 2 iterations CdOld = self.vsolver.Cd
it += 1 it += 1
# save results # save results
self.isolver.save(0, self.writer) self.isolver.save(0, self.writer)
print('\n') print('\n')
#!/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")
...@@ -22,32 +22,347 @@ ...@@ -22,32 +22,347 @@
from __future__ import print_function from __future__ import print_function
from builtins import range from builtins import range
from builtins import object 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 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): class Solver(object):
def __init__(self, _boundary): def __init__(self, _Re):
'''... '''...
''' '''
# TODO generalize to handle several boundaries (use python array) self.Re = _Re # Reynolds number
self.boundary = _boundary # boundary on which BLE are sovled self.gamma = 1.4
self.nN = len(self.boundary.nodes) # number of nodes self.Cd = 0
self.nE = len(self.boundary.groups[0].tag.elems) # number of elements
self.v = np.zeros((self.nN, 3)) # velocity at edge of BL def run(self, group):
self.u = np.zeros(self.nE) # blowing velocity
def run(self):
''' Solve BLE ''' Solve BLE
''' '''
print('--- I am a fake BL solver and I am setting dummy blowing velocities! ---') def writeFile(data, H, theta, Cf, vt, Cdissip, Ct, Dstar, blwUe):
# 6th order fitting of delta_star for naca 0012 at alpha = 0, from XFoil # save_path = 'D:/Amaury/Documents/TFE/xfoil/'
a = 0.000054925976379640116 name = 'BLparameters'
b = 0.007407077078697043 name2 = 'blwU'
c = -0.043918316829113194 # complete_name = os.path.join(save_path, name+".dat")
d = 0.1487463222137225 # complete_name2 = os.path.join(save_path, name2+".dat")
e = -0.20162146328539832 f = open(name, 'w+')
f = 0.0880698637668166 for i in range(len(H)):
g = 0.004436812366681563 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]))
for i in range(0, self.nE): f.close()
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]) f = open(name2, 'w+')
# CAUTION: a positive velocity will result in a suction! for i in range(len(blwUe)):
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 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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment