Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • Paul.Dechamps/dartflo
1 result
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 1591 deletions
vii/__init__.py deleted 100644 → 0
View file @ 3a33458d
# -*- coding: utf-8 -*-
# dart MODULE initialization file
# 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.
import fwk
import tbox
from viiw import *
vii/_src/CMakeLists.txt deleted 100644 → 0
View file @ 3a33458d
# 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.
# CMake input file of the SWIG wrapper around "viiw.so"
INCLUDE(${SWIG_USE_FILE})
FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
FILE(GLOB ISRCS *.i)
SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
SET(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_FLAGS} "-interface" "_viiw") # avoids "import _module_d" with MSVC/Debug
SWIG_ADD_LIBRARY(viiw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
SET_PROPERTY(TARGET viiw PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON)
MACRO_DebugPostfix(viiw)
TARGET_INCLUDE_DIRECTORIES(viiw PRIVATE ${PROJECT_SOURCE_DIR}/vii/_src
${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
${PYTHON_INCLUDE_PATH})
TARGET_LINK_LIBRARIES(viiw PRIVATE vii tbox fwk ${PYTHON_LIBRARIES})
INSTALL(FILES ${EXECUTABLE_OUTPUT_PATH}/\${BUILD_TYPE}/viiw.py DESTINATION ${CMAKE_INSTALL_PREFIX})
INSTALL(TARGETS viiw DESTINATION ${CMAKE_INSTALL_PREFIX})
vii/_src/viiw.h deleted 100644 → 0
View file @ 3a33458d
/*
* 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.
*/
#include "DDriver.h"
\ No newline at end of file
vii/_src/viiw.i deleted 100644 → 0
View file @ 3a33458d
/*
* 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.
*/
// SWIG input file of the 'dart' module
%feature("autodoc","1");
%module(docstring=
"'viiw' module: Viscous inviscid interaction for dart 2021-2022
(c) ULiege - A&M",
directors="0",
threads="1"
) viiw
%{
#include "vii.h"
#include "fwkw.h"
#include "tboxw.h"
#include "viiw.h"
%}
%include "fwkw.swg"
// ----------- MODULES UTILISES ------------
%import "tboxw.i"
// ----------- VII CLASSES ----------------
%include "vii.h"
%shared_ptr(vii::Driver);
%immutable vii::Driver::Re; // read-only variable (no setter)
%immutable vii::Driver::Cdt;
%immutable vii::Driver::Cdt_sec;
%immutable vii::Driver::Cdf;
%immutable vii::Driver::Cdf_sec;
%immutable vii::Driver::Cdp;
%immutable vii::Driver::Cdp_sec;
%immutable vii::Driver::CFL0;
%immutable vii::Driver::verbose;
%include "DDriver.h"
\ No newline at end of file
vii/api/__init__.py deleted 100644 → 0
View file @ 3a33458d
# -*- coding: utf-8 -*-
vii/api/core.py deleted 100644 → 0
View file @ 3a33458d
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2022 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.
## Initialize VII computation
# Paul Dechamps
def initVII(cfg, icfg, iSolverName='DART'):
"""
Inputs
------
- cfg: Dict with options
- Re (float): Flow Reynolds number.
- Minf (float): Freestream Mach number (only used for time step computation).
- CFL0 (float): Initial CFL number.
- Verb (int): Verbosity level of the viscous solver.
- couplIter (int): Maximum number of coupling iterations.
- couplTol (float): Tolerance to end computation (drag count between 2 iterations).
- iterPrint (int): Number of iterations between which the coupler outputs its state.
- resetInv (bool): Resets the inviscid solver at each iteration (True), reuses previous
state as initial condition (False).
Note
All inputs have default values except Re.
- icg: Dict with options. Inviscid solver configuration.
- iSolverName (string): Name of the inviscid solver.
Outputs
-------
dict with keys
- coupler: Coupler python object to run the viscous-inviscid algorithm.
- solversAPI: Interface between the solver objects to ensure proper communication.
- vSolver: vii::Driver object to run the viscous computation.
- iSolverObjects: Dict containing.
- All inviscid solver dependencies (depend on iSolverName).
@todo: - Correctly handle 3D computation parameters.
"""
# Imports
from vii.coupler import Coupler
import vii
if iSolverName == 'DART':
from vii.interfaces.dart.DartInterface import DartInterface as interface
from dart.api.core import initDart
iSolverObjects = initDart(icfg, scenario=icfg['scenario'], task=icfg['task'], viscous = 1)
# Check viscous solver parameters.
if 'Re' in cfg and cfg['Re'] > 0:
_Re = cfg['Re']
else:
raise RuntimeError('Missing or invalid Reynolds number')
if 'Minf' in cfg and cfg['Minf'] > 0:
_Minf = cfg['Minf']
else:
_Minf = 0.1
if 'CFL0' in cfg and cfg['CFL0'] > 0:
_CFL0 = cfg['CFL0']
else:
_CFL0 = 1
if 'Verb' in cfg and 0<= cfg['Verb'] <= 3:
_verbose = cfg['Verb']
else:
_verbose = 1
_span = 0
_nSections = 1
# Check coupler parameters.
if 'couplIter' in cfg and cfg['couplIter'] > 0:
__couplIter = cfg['couplIter']
else:
__couplIter = 150
if 'couplTol' in cfg:
__couplTol = cfg['couplTol']
else:
__couplTol = 1e-4
if 'iterPrint' in cfg:
__iterPrint = cfg['iterPrint']
else:
__iterPrint = 1
if 'resetInv' in cfg:
__resetInv = cfg['resetInv']
else:
__resetInv = False
# Viscous solver object.
vSolver = vii.Driver(_Re, _Minf, _CFL0, _nSections, _span, verbose =_verbose)
# Solvers interface.
solversAPI = interface(iSolverObjects['sol'], vSolver, [iSolverObjects['blwb'], iSolverObjects['blww']])
# Coupler
coupler = Coupler(solversAPI, vSolver, _maxCouplIter = __couplIter, _couplTol = __couplTol, _iterPrint = __iterPrint, _resetInv = __resetInv)
return {'coupler': coupler,
'solversAPI': solversAPI,
'vSolver': vSolver,
'iSolverObjects': iSolverObjects}
vii/coupler.py deleted 100644 → 0
View file @ 3a33458d
#!/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.
# VII Coupler
# Paul Dechamps
import fwk
from fwk.coloring import ccolors
import math
import numpy as np
class Coupler:
def __init__(self, iSolverAPI, vSolver, _maxCouplIter=150, _couplTol=1e-4, _iterPrint=1, _resetInv=False):
self.iSolverAPI = iSolverAPI
self.vSolver = vSolver
self.maxCouplIter = _maxCouplIter
self.couplTol = _couplTol
self.resetInv = _resetInv
self.iterPrint = _iterPrint if self.iSolverAPI.getVerbose() == 0 and self.vSolver.verbose == 0 else 1
self.tms = fwk.Timers()
def run(self):
# Aerodynamic coefficients.
aeroCoeffs = np.zeros((0, 2))
# Convergence parameters.
couplIter = 0
cdPrev = 0.0
while couplIter < self.maxCouplIter:
# Run inviscid solver.
self.tms['inviscid'].start()
if self.resetInv:
self.iSolverAPI.reset()
self.iSolverAPI.run()
self.tms['inviscid'].stop()
# Write inviscid Cp distribution.
if couplIter == 0:
self.iSolverAPI.writeCp()
# Impose inviscid boundary in the viscous solver.
self.tms['processing'].start()
self.iSolverAPI.imposeInvBC(couplIter)
self.tms['processing'].stop()
# Run viscous solver.
self.tms['viscous'].start()
exitCode = self.vSolver.run(couplIter)
self.tms['viscous'].stop()
# Impose blowing boundary condition in the inviscid solver.
self.tms['processing'].start()
self.iSolverAPI.imposeBlwVel()
self.tms['processing'].stop()
aeroCoeffs = np.vstack((aeroCoeffs, [self.iSolverAPI.getCl(), self.vSolver.Cdt]))
# Check convergence status.
error = abs((self.vSolver.Cdt - cdPrev) / self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
if error <= self.couplTol:
print(ccolors.ANSI_GREEN, '{:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.2f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.getxtr(0, 0), self.vSolver.getxtr(0, 1), np.log10(error)), ccolors.ANSI_RESET)
return aeroCoeffs
cdPrev = self.vSolver.Cdt
if couplIter == 0:
print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'.format(self.iSolverAPI.getAoA()*180/math.pi, self.iSolverAPI.getMinf(), self.vSolver.getRe()/1e6))
print('{:>5s}| {:>10s} {:>10s} | {:>10s} {:>8s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
print(' ----------------------------------------------------------')
if couplIter % self.iterPrint == 0:
print(' {:>4.0f}| {:>10.4f} {:>10.4f} | {:>10.4f} {:>8.4f} {:>8.3f}'.format(couplIter, self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.getxtr(0, 0), self.vSolver.getxtr(0, 1), np.log10(error)))
if self.iSolverAPI.getVerbose() != 0 or self.vSolver.verbose != 0:
print('')
couplIter += 1
self.tms['processing'].stop()
else:
print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
print('')
print('{:>5s}| {:>10s} {:>10s} | {:>10s} {:>8s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
print(' ----------------------------------------------------------')
print(ccolors.ANSI_RED, '{:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.2f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.getxtr(0, 0), self.vSolver.getxtr(0, 1), np.log10(error)), ccolors.ANSI_RESET)
return aeroCoeffs
\ No newline at end of file
vii/interfaces/dart/DartInterface.py deleted 100644 → 0
View file @ 3a33458d
#!/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.
# Dart interface
# Paul Dechamps
import numpy as np
import dart.utils as floU
import tbox.utils as tboxU
from ..viscous.DAirfoil import Airfoil
from ..viscous.DWake import Wake
from ..viscous import DGroupBuilder
class DartInterface:
def __init__(self, dartSolver, vSolver, blw):
self.solver = dartSolver
self.vSolver = vSolver
self.createProblem(blw[0], blw[1])
self.nElmUpperPrev = 0
def createProblem(self, _boundaryAirfoil, _boundaryWake):
"""Create data structure for information transfer.
Args
----
- _boundaryAirfoil: dart::Blowing data structure for the airfoil.
- _boundaryWake: dart::Blowing data structure for the wake.
"""
self.group = [Airfoil(_boundaryAirfoil), Wake(_boundaryWake)]
self.blwVel = [None, None, None]
self.deltaStar = [None, None, None]
self.xx = [None, None, None]
def run(self):
return self.solver.run()
def writeCp(self):
""" Write Cp distribution around the airfoil on file.
"""
Cp = floU.extract(self.group[0].boundary.tag.elems, self.solver.Cp)
tboxU.write(Cp, 'Cp_Inviscid.dat', '%1.5e',', ', 'x, y, z, Cp', '')
def save(self, writer):
self.solver.save(writer)
def getAoA(self):
return self.solver.pbl.alpha
def getMinf(self):
return self.solver.pbl.M_inf
def getCl(self):
return self.solver.Cl
def getCm(self):
return self.solver.Cm
def getVerbose(self):
return self.solver.verbose
def reset(self):
self.solver.reset()
def imposeInvBC(self, couplIter):
""" Extract the inviscid paramaters at the edge of the boundary layer
and use it as a boundary condition in the viscous solver.
Params
------
- couplIter (int): Coupling iteration.
"""
# Retreive parameters.
for n in range(0, len(self.group)):
for i in range(0, len(self.group[n].boundary.nodes)):
self.group[n].v[i,0] = self.solver.U[self.group[n].boundary.nodes[i].row][0]
self.group[n].v[i,1] = self.solver.U[self.group[n].boundary.nodes[i].row][1]
self.group[n].v[i,2] = self.solver.U[self.group[n].boundary.nodes[i].row][2]
self.group[n].Me[i] = self.solver.M[self.group[n].boundary.nodes[i].row]
self.group[n].rhoe[i] = self.solver.rho[self.group[n].boundary.nodes[i].row]
dataIn = [None, None]
for n in range(0, len(self.group)):
self.group[n].connectListNodes, self.group[n].connectListElems, dataIn[n] = self.group[n].connectList(couplIter)
self.fMeshDict, _, _ = DGroupBuilder.mergeVectors(dataIn, 0, 1)
fMeshDict = self.fMeshDict.copy()
if ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) != self.nElmUpperPrev) or couplIter == 0):
# Initialize mesh
self.vSolver.setGlobMesh(0, 0, fMeshDict['globCoord'][:fMeshDict['bNodes'][1]], fMeshDict['yGlobCoord'][:fMeshDict['bNodes'][1]], fMeshDict['zGlobCoord'][:fMeshDict['bNodes'][1]])
self.vSolver.setGlobMesh(0, 1, fMeshDict['globCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['yGlobCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['zGlobCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.setGlobMesh(0, 2, fMeshDict['globCoord'][fMeshDict['bNodes'][2]:], fMeshDict['yGlobCoord'][fMeshDict['bNodes'][2]:], fMeshDict['zGlobCoord'][fMeshDict['bNodes'][2]:])
# Initialize parameters at the edge of the boundary layer.
self.vSolver.setxxExt(0, 0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]])
self.vSolver.setxxExt(0, 1, fMeshDict['locCoord'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.setxxExt(0, 2, fMeshDict['locCoord'][fMeshDict['bNodes'][2]:])
self.vSolver.setDeltaStarExt(0, 0, fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
self.vSolver.setDeltaStarExt(0, 1, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.setDeltaStarExt(0, 2, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
else:
# Parameters at the edge of the boundary layer only.
self.vSolver.setxxExt(0, 0, fMeshDict['xxExt'][:fMeshDict['bNodes'][1]])
self.vSolver.setxxExt(0, 1, fMeshDict['xxExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.setxxExt(0, 2, fMeshDict['xxExt'][fMeshDict['bNodes'][2]:])
self.vSolver.setDeltaStarExt(0, 0, fMeshDict['deltaStarExt'][:fMeshDict['bNodes'][1]])
self.vSolver.setDeltaStarExt(0, 1, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.setDeltaStarExt(0, 2, fMeshDict['deltaStarExt'][fMeshDict['bNodes'][2]:])
# Inviscid boundary condition.
self.nElmUpperPrev = len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]])
self.vSolver.setVelocities(0, 0, fMeshDict['vx'][:fMeshDict['bNodes'][1]], fMeshDict['vy'][:fMeshDict['bNodes'][1]], fMeshDict['vz'][:fMeshDict['bNodes'][1]])
self.vSolver.setVelocities(0, 1, fMeshDict['vx'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vy'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['vz'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.setVelocities(0, 2, fMeshDict['vx'][fMeshDict['bNodes'][2]:], fMeshDict['vy'][fMeshDict['bNodes'][2]:], fMeshDict['vz'][fMeshDict['bNodes'][2]:])
self.vSolver.setMe(0, 0, fMeshDict['Me'][:fMeshDict['bNodes'][1]])
self.vSolver.setMe(0, 1, fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.setMe(0, 2, fMeshDict['Me'][fMeshDict['bNodes'][2]:])
self.vSolver.setRhoe(0, 0, fMeshDict['rho'][:fMeshDict['bNodes'][1]])
self.vSolver.setRhoe(0, 1, fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
self.vSolver.setRhoe(0, 2, fMeshDict['rho'][fMeshDict['bNodes'][2]:])
def imposeBlwVel(self):
""" Extract the solution of the viscous calculation (blowing velocity)
and use it as a boundary condition in the inviscid solver.
"""
# Retreive and set blowing velocities.
for iRegion in range(3):
self.blwVel[iRegion] = np.asarray(self.vSolver.extractBlowingVel(0, iRegion))
self.deltaStar[iRegion] = np.asarray(self.vSolver.extractDeltaStar(0, iRegion))
self.xx[iRegion] = np.asarray(self.vSolver.extractxx(0, iRegion))
blwVelAF = np.concatenate((self.blwVel[0], self.blwVel[1]))
# Remove stagnation point doublon.
deltaStarAF = np.concatenate((self.deltaStar[0], self.deltaStar[1][1:]))
xxAF = np.concatenate((self.xx[0], self.xx[1][1:]))
# Blowing velocity, displacement thickness and mesh (local coordinates) distributions in the wake.
blwVelWK = self.blwVel[2]
deltaStarWK = self.deltaStar[2]
xxWK = self.xx[2]
# Update data structures for information transfer.
self.group[0].u = blwVelAF[self.group[0].connectListElems.argsort()]
self.group[1].u = blwVelWK[self.group[1].connectListElems.argsort()]
self.group[0].deltaStar = deltaStarAF[self.group[0].connectListNodes.argsort()]
self.group[1].deltaStar = deltaStarWK[self.group[1].connectListNodes.argsort()]
self.group[0].xx = xxAF[self.group[0].connectListNodes.argsort()]
self.group[1].xx = xxWK[self.group[1].connectListNodes.argsort()]
# Impose blowing velocities.
for n in range(0, len(self.group)):
for i in range(0, self.group[n].nE):
self.group[n].boundary.setU(i, self.group[n].u[i])
\ No newline at end of file
vii/interfaces/dart/__init__.py deleted 100644 → 0
View file @ 3a33458d
# -*- coding: utf-8 -*-
\ No newline at end of file
vii/interfaces/viscous/DAirfoil.py deleted 100755 → 0
View file @ 3a33458d
#!/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 .DBoundary 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, couplIter):
""" Sort the value read by the viscous solver/ Create list of connectivity
"""
N1 = np.zeros(self.nN, dtype=int) # Node number
connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
data = np.zeros((self.boundary.nodes.size(), 10))
i = 0
for n in self.boundary.nodes:
data[i,0] = n.no
data[i,1] = n.pos[0]
data[i,2] = n.pos[1]
data[i,3] = n.pos[2]
data[i,4] = self.v[i,0]
data[i,5] = self.v[i,1]
data[i,6] = self.Me[i]
data[i,7] = self.rhoe[i]
data[i,8] = self.deltaStar[i]
data[i,9] = self.xx[i]
i += 1
# Table containing the element and its nodes
eData = np.zeros((self.nE,3), dtype=int)
for i in range(0, self.nE):
eData[i,0] = self.boundary.tag.elems[i].no
eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
# Find the stagnation point
if couplIter == 0:
self.idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
idxStag = self.idxStag
globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
# Find TE nodes (assuming suction side element is above pressure side element in the y direction)
idxTE = np.where(data[:,1] == np.max(data[:,1]))[0] # TE nodes
idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], eData[:,2] == data[idxTE[0],0]))[0] # TE element 1
idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], eData[:,2] == data[idxTE[1],0]))[0] # TE element 2
ye0 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe0[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe0[0],2])[0],2]) # y coordinates of TE element 1 CG
ye1 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe1[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe1[0],2])[0],2]) # y coordinates of TE element 2 CG
if ye0 - ye1 > 0: # element 1 containing TE node 1 is above element 2
upperGlobTE = int(data[idxTE[0],0])
lowerGlobTE = int(data[idxTE[1],0])
else:
upperGlobTE = int(data[idxTE[1],0])
lowerGlobTE = int(data[idxTE[0],0])
# Connectivity
connectListElems[0] = np.where(eData[:,1] == globStag)[0]
N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes -> globStag
connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
i = 1
upperTE = 0
lowerTE = 0
# Sort the suction part
while upperTE == 0:
N1[i] = eData[connectListElems[i-1],2] # Second node of the element
connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes
connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
if eData[connectListElems[i],2] == upperGlobTE:
upperTE = 1
i += 1
# Sort the pressure side
connectListElems[i] = np.where(eData[:,2] == globStag)[0]
connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
N1[i] = eData[connectListElems[i],2]
while lowerTE == 0:
N1[i+1] = eData[connectListElems[i],1] # First node of the element
connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes
connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
if eData[connectListElems[i+1],1] == lowerGlobTE:
lowerTE = 1
i += 1
connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
data[:,0] = data[connectListNodes,0]
data[:,1] = data[connectListNodes,1]
data[:,2] = data[connectListNodes,2]
data[:,3] = data[connectListNodes,3]
data[:,4] = data[connectListNodes,4]
data[:,5] = data[connectListNodes,5]
data[:,6] = data[connectListNodes,6]
data[:,7] = data[connectListNodes,7]
data[:,8] = data[connectListNodes,8]
data[:,9] = data[connectListNodes,9]
# Separated upper and lower part
data = np.delete(data,0,1)
uData = data[0:np.argmax(data[:,0])+1]
lData = data[np.argmax(data[:,0])+1:None]
lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
return connectListNodes, connectListElems,[uData, lData]
vii/interfaces/viscous/DBoundary.py deleted 100755 → 0
View file @ 3a33458d
#!/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
vii/interfaces/viscous/DGroupBuilder.py deleted 100755 → 0
View file @ 3a33458d
import numpy as np
def mergeVectors(dataIn, mapMesh, nFactor):
"""Merges 3 groups upper/lower sides and wake.
"""
for k in range(len(dataIn)):
if len(dataIn[k]) == 2: # Airfoil case.
xx_up, dv_up, vtExt_up, _, alpha_up = getBLcoordinates(dataIn[k][0])
xx_lw, dv_lw, vtExt_lw, _, alpha_lw = getBLcoordinates(dataIn[k][1])
else: # Wake case
xx_wk, dv_wk, vtExt_wk, _, alpha_wk = getBLcoordinates(dataIn[k])
if mapMesh:
# Save initial mesh.
cMesh = np.concatenate((dataIn[0][0][:, 0], dataIn[0][1][1:, 0], dataIn[1][:, 0]))
cMeshbNodes = [0, len(dataIn[0][0][:, 0]), len(dataIn[0][0][:, 0]) + len(dataIn[0][1][1:, 0])]
cxx = np.concatenate((xx_up, xx_lw[1:], xx_wk))
cvtExt = np.concatenate((vtExt_up, vtExt_lw[1:], vtExt_wk))
cxxExt = np.concatenate((dataIn[0][0][:, 8], dataIn[0][1][1:, 8], dataIn[1][:, 8]))
# Create fine mesh.
fMeshUp = createFineMesh(dataIn[0][0][:, 0], nFactor)
fMeshLw = createFineMesh(dataIn[0][1][:, 0], nFactor)
fMeshWk = createFineMesh(dataIn[1][:, 0], nFactor)
fMesh = np.concatenate((fMeshUp, fMeshLw[1:], fMeshWk))
fMeshbNodes = [0, len(fMeshUp), len(fMeshUp) + len(fMeshLw[1:])]
fxx_up = interpolateToFineMesh(xx_up, fMeshUp, nFactor)
fxx_lw = interpolateToFineMesh(xx_lw, fMeshLw, nFactor)
fxx_wk = interpolateToFineMesh(xx_wk, fMeshWk, nFactor)
fxx = np.concatenate((fxx_up, fxx_lw[1:], fxx_wk))
fvtExt_up = interpolateToFineMesh(vtExt_up, fMeshUp, nFactor)
fvtExt_lw = interpolateToFineMesh(vtExt_lw, fMeshLw, nFactor)
fvtExt_wk = interpolateToFineMesh(vtExt_wk, fMeshWk, nFactor)
fvtExt = np.concatenate((fvtExt_up, fvtExt_lw[1:], fvtExt_wk))
fMe_up = interpolateToFineMesh(dataIn[0][0][:, 5], fMeshUp, nFactor)
fMe_lw = interpolateToFineMesh(dataIn[0][1][:, 5], fMeshLw, nFactor)
fMe_wk = interpolateToFineMesh(dataIn[1][:, 5], fMeshWk, nFactor)
fMe = np.concatenate((fMe_up, fMe_lw[1:], fMe_wk))
frho_up = interpolateToFineMesh(dataIn[0][0][:, 6], fMeshUp, nFactor)
frho_lw = interpolateToFineMesh(dataIn[0][1][:, 6], fMeshLw, nFactor)
frho_wk = interpolateToFineMesh(dataIn[1][:, 6], fMeshWk, nFactor)
frho = np.concatenate((frho_up, frho_lw[1:], frho_wk))
fdStarExt_up = interpolateToFineMesh(dataIn[0][0][:, 7], fMeshUp, nFactor)
fdStarExt_lw = interpolateToFineMesh(dataIn[0][1][:, 7], fMeshLw, nFactor)
fdStarExt_wk = interpolateToFineMesh(dataIn[1][:, 7], fMeshWk, nFactor)
fxxExt_up = interpolateToFineMesh(dataIn[0][0][:, 8], fMeshUp, nFactor)
fxxExt_lw = interpolateToFineMesh(dataIn[0][1][:, 8], fMeshLw, nFactor)
fxxExt_wk = interpolateToFineMesh(dataIn[1][:, 8], fMeshWk, nFactor)
fdStarext = np.concatenate((fdStarExt_up, fdStarExt_lw[1:], fdStarExt_wk))
fxxext = np.concatenate((fxxExt_up, fxxExt_lw[1:], fxxExt_wk))
fdv = np.concatenate((dv_up, dv_lw[1:], dv_wk))
# Create dictionnaries for the fine and coarse meshes.
fMeshDict = {'globCoord': fMesh, 'bNodes': fMeshbNodes, 'locCoord': fxx,
'vtExt': fvtExt, 'Me': fMe, 'rho': frho, 'deltaStarExt': fdStarext, 'xxExt': fxxext, 'dv': fdv}
cMe = np.concatenate((dataIn[0][0][:, 5], dataIn[0][1][1:, 5], dataIn[1][:, 5]))
cRho = np.concatenate((dataIn[0][0][:, 6], dataIn[0][1][1:, 6], dataIn[1][:, 6]))
cdStarExt = np.concatenate((dataIn[0][0][:, 7], dataIn[0][1][1:, 7], dataIn[1][:, 7]))
cMeshDict = {'globCoord': cMesh, 'bNodes': cMeshbNodes, 'locCoord': cxx,
'vtExt': cvtExt, 'Me': cMe, 'rho': cRho, 'deltaStarExt': cdStarExt, 'xxExt': cxxExt, 'dv': fdv}
dataUpper = np.zeros((len(fMeshUp), 0))
dataLower = np.zeros((len(fMeshLw), 0))
dataWake = np.zeros((len(fMeshWk), 0))
for iParam in range(len(dataIn[0][0][0, :])):
interpParamUp = interpolateToFineMesh(dataIn[0][0][:, iParam], fMeshUp, nFactor)
dataUpper = np.column_stack((dataUpper, interpParamUp))
interpParamLw = interpolateToFineMesh(dataIn[0][1][:, iParam], fMeshLw, nFactor)
dataLower = np.column_stack((dataLower, interpParamLw))
interpParamWk = interpolateToFineMesh(dataIn[1][:, iParam], fMeshWk, nFactor)
dataWake = np.column_stack((dataWake, interpParamWk))
# Remove stagnation point doublon.
dataLower = np.delete(dataLower, 0, axis=0)
else:
Mesh = np.concatenate((dataIn[0][0][:, 0], dataIn[0][1][:, 0], dataIn[1][:, 0]))
Meshy = np.concatenate((dataIn[0][0][:, 1], dataIn[0][1][:, 1], dataIn[1][:, 1]))
Meshz = np.concatenate((dataIn[0][0][:, 2], dataIn[0][1][:, 2], dataIn[1][:, 2]))
MeshbNodes = [0, len(dataIn[0][0][:, 0]), len(dataIn[0][0][:, 0]) + len(dataIn[0][1][:, 0]), len(dataIn[0][0][:, 0]) + len(dataIn[0][1][:, 0]) + len(dataIn[1][:, 0])]
# Concatenate all parameters (upper side, lower side, wake)
xx = np.concatenate((xx_up, xx_lw, xx_wk))
vtExt = np.concatenate((vtExt_up, vtExt_lw, vtExt_wk))
vx = np.concatenate((dataIn[0][0][:, 3], dataIn[0][1][:, 3], dataIn[1][:, 3]))
vy = np.concatenate((dataIn[0][0][:, 4], dataIn[0][1][:, 4], dataIn[1][:, 4]))
vz = np.concatenate((dataIn[0][0][:, 5], dataIn[0][1][:, 5], dataIn[1][:, 5]))
alpha = np.concatenate((alpha_up, alpha_lw, alpha_wk))
dv = np.concatenate((dv_up, dv_lw, dv_wk))
Me = np.concatenate((dataIn[0][0][:, 5], dataIn[0][1][:, 5], dataIn[1][:, 5]))
rho = np.concatenate((dataIn[0][0][:, 6], dataIn[0][1][:, 6], dataIn[1][:, 6]))
dStarext = np.concatenate((dataIn[0][0][:, 7], dataIn[0][1][:, 7], dataIn[1][:, 7]))
xxExt = np.concatenate((dataIn[0][0][:, 8], dataIn[0][1][:, 8], dataIn[1][:, 8]))
# Create dictionnaries for the fine and coarse meshes which are equivalent if meshMap is disabled.
fMeshDict = {'globCoord': Mesh, 'yGlobCoord': Meshy, 'zGlobCoord': Meshz, 'bNodes': MeshbNodes, 'locCoord': xx,
'vtExt': vtExt, 'vx': vx, 'vy': vy, 'vz': vz, 'Me': Me, 'rho': rho, 'deltaStarExt': dStarext, 'xxExt': xxExt, 'dv': dv}
cMeshDict = fMeshDict.copy()
dataUpper = dataIn[0][0]
dataLower = dataIn[0][1][1:, :]
dataWake = dataIn[1]
data = np.concatenate((dataUpper, dataLower, dataWake), axis=0)
return fMeshDict, cMeshDict, data
def getBLcoordinates(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) # Outer flow velocity.
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]
# Tangential speed at node 1 element i
vt[i] = (data[i, 3] * np.cos(alpha[i]) + data[i, 4] * np.sin(alpha[i]))
# Tangential speed at node 2 element i
vt[i+1] = (data[i+1, 3] * np.cos(alpha[i]) + data[i+1, 4] * np.sin(alpha[i]))
# Velocity gradient according to element i.
dv[i] = (vt[i+1] - vt[i]) / dx[i]
return xx, dv, vt, nN, alpha
def interpolateToFineMesh(cVector, fMesh, nFactor):
""" Interpolate variables of cVector on the fine mesh fMesh.
Params
------
- cVector (numpy.ndarray): Distribution to be interpolated.
- fMesh (numpy.ndarray): Local tangential coordinate (xi) of the fine mesh.
- nFactor (int): Number of times each cell was split when creating the fine mesh.
"""
fVector = np.zeros(len(fMesh)-1)
for cPoint in range(len(cVector) - 1):
dv = cVector[cPoint + 1] - cVector[cPoint]
for fPoint in range(nFactor):
fVector[nFactor * cPoint + fPoint] = cVector[cPoint] + \
fPoint * dv / nFactor
fVector = np.append(fVector, cVector[-1])
return fVector
def createFineMesh(cMesh, nFactor):
""" Create fine mesh distribution in the local frame of reference.
Params
------
- cMesh (numpy.ndarray): Initial mesh (used for the inviscid calculations.
- nFactor (int): Number of times each cell will be split to create the fine mesh.
"""
fMesh = []
for iPoint in range(len(cMesh) - 1):
dx = cMesh[iPoint + 1] - cMesh[iPoint]
for iRefinedPoint in range(nFactor):
fMesh = np.append(
fMesh, cMesh[iPoint] + iRefinedPoint * dx / nFactor)
fMesh = np.append(fMesh, cMesh[-1])
return fMesh
vii/interfaces/viscous/DWake.py deleted 100755 → 0
View file @ 3a33458d
#!/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 & Paul Dechamps
from .DBoundary import Boundary
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, couplIter):
""" Sort the value read by the viscous solver/ Create list of connectivity
"""
N1 = np.zeros(self.nN, dtype=int) # Node number
connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
data = np.zeros((self.boundary.nodes.size(), 10))
i = 0
for n in self.boundary.nodes:
data[i,0] = n.no
data[i,1] = n.pos[0]
data[i,2] = n.pos[1]
data[i,3] = n.pos[2]
data[i,4] = self.v[i,0]
data[i,5] = self.v[i,1]
data[i,6] = self.Me[i]
data[i,7] = self.rhoe[i]
data[i,8] = self.deltaStar[i]
data[i,9] = self.xx[i]
i += 1
# Table containing the element and its nodes
eData = np.zeros((self.nE,3), dtype=int)
for i in range(0, self.nE):
eData[i,0] = self.boundary.tag.elems[i].no
eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
connectListNodes = data[:,1].argsort()
# connectListElems[0] = np.where(eData[:,1] == min(data[:,1]))[0]
connectListElems[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
vii/src/CMakeLists.txt deleted 100644 → 0
View file @ 3a33458d
# 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.
# CMake input file of dart.so
FILE(GLOB SRCS *.h *.cpp *.inl *.hpp)
ADD_LIBRARY(vii SHARED ${SRCS})
MACRO_DebugPostfix(vii)
TARGET_INCLUDE_DIRECTORIES(vii PUBLIC ${PROJECT_SOURCE_DIR}/vii/src)
TARGET_LINK_LIBRARIES(vii tbox)
INSTALL(TARGETS vii DESTINATION ${CMAKE_INSTALL_PREFIX})
SOURCE_GROUP(base REGULAR_EXPRESSION ".*\\.(cpp|inl|hpp|h)")
vii/src/DBoundaryLayer.cpp deleted 100644 → 0
View file @ 3a33458d
#include "DBoundaryLayer.h"
#include "DClosures.h"
#include "DDiscretization.h"
#include "DEdge.h"
#include <cmath>
#include <iomanip>
using namespace vii;
/**
* @brief Construct a new BoundaryLayer::BoundaryLayer object.
*
* @param _xtrF Forced transition location (default=-1; free transition).
*/
BoundaryLayer::BoundaryLayer(double _xtrF)
{
if (_xtrF < 0. && _xtrF != -1.)
throw std::runtime_error("Fatal error: Invalid transition location.\n");
xtrF = _xtrF;
xtr = 1.0;
closSolver = new Closures();
mesh = new Discretization();
blEdge = new Edge();
transMarker = 0;
}
BoundaryLayer::~BoundaryLayer()
{
delete closSolver;
delete mesh;
delete blEdge;
// std::cout << "~vii::BoundaryLayer()\n";
}
/**
* @brief Set the mesh and resize all boundary layer quantities accordingly.
*
* @param x Position x in the global frame of reference.
* @param y Position y in the global frame of reference.
* @param z Position z in the global frame of reference.
*/
void BoundaryLayer::setMesh(std::vector<double> x,
std::vector<double> y,
std::vector<double> z)
{
size_t nMarkers = x.size();
mesh->setGlob(x, y, z);
// Resize and reset all variables if needed.
if (mesh->getDiffnElm() != 0)
{
u.resize(nVar * nMarkers, 0.);
Ret.resize(nMarkers, 0.);
cd.resize(nMarkers, 0.);
cf.resize(nMarkers, 0.);
Hstar.resize(nMarkers, 0.);
Hstar2.resize(nMarkers, 0.);
Hk.resize(nMarkers, 0.);
ctEq.resize(nMarkers, 0.);
us.resize(nMarkers, 0.);
deltaStar.resize(nMarkers, 0.);
delta.resize(nMarkers, 0.);
regime.resize(nMarkers, 0);
}
}
/**
* @brief Set the boundary condition at the stagnation point.
*
* @param Re Freestream Reynolds number.
*/
void BoundaryLayer::setStagBC(double Re)
{
double dv0 = (blEdge->getVt(1) - blEdge->getVt(0)) / (mesh->getLoc(1) - mesh->getLoc(0));
if (dv0 > 0)
u[0] = sqrt(0.075 / (Re * dv0)); // Theta
else
u[0] = 1e-6; // Theta
u[1] = 2.23; // H
u[2] = 0.; // N
u[3] = blEdge->getVt(0); // ue
u[4] = 0.; // Ctau
Ret[0] = u[3] * u[0] * Re;
if (Ret[0] < 1)
{
Ret[0] = 1.;
u[3] = Ret[0] / (u[0] * Re);
}
deltaStar[0] = u[0] * u[1];
std::vector<double> lamParam(8, 0.);
closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], blEdge->getMe(0), name);
Hstar[0] = lamParam[0];
Hstar2[0] = lamParam[1];
Hk[0] = lamParam[2];
cf[0] = lamParam[3];
cd[0] = lamParam[4];
delta[0] = lamParam[5];
ctEq[0] = lamParam[6];
us[0] = lamParam[7];
}
/**
* @brief Set the boundary condition at the first wake point.
*
* @param Re Freestream Reynolds number
* @param UpperCond Variables at the last point on the upper side.
* @param LowerCond Variables at the last point on the lower side.
*/
void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond, std::vector<double> LowerCond)
{
if (name != "wake")
std::cout << "Warning: Imposing wake boundary condition on " << name << std::endl;
u[0] = UpperCond[0] + LowerCond[0]; // (Theta_up+Theta_low).
u[1] = (UpperCond[0] * UpperCond[1] + LowerCond[0] * LowerCond[1]) / u[0]; // ((Theta*H)_up+(Theta*H)_low)/(Theta_up+Theta_low).
u[2] = 9.; // Amplification factor.
u[3] = blEdge->getVt(0); // Inviscid velocity.
u[4] = (UpperCond[0] * UpperCond[4] + LowerCond[0] * LowerCond[4]) / u[0]; // ((Ct*Theta)_up+(Theta)_low)/(Ct*Theta_up+Theta_low).
Ret[0] = u[3] * u[0] * Re; // Reynolds number based on the momentum thickness.
if (Ret[0] < 1)
{
Ret[0] = 1;
u[3] = Ret[0] / (u[0] * Re);
}
deltaStar[0] = u[0] * u[1];
// Laminar closures.
std::vector<double> lamParam(8, 0.);
closSolver->laminarClosures(lamParam, u[0], u[1], Ret[0], blEdge->getMe(0), name);
Hstar[0] = lamParam[0];
Hstar2[0] = lamParam[1];
Hk[0] = lamParam[2];
cf[0] = lamParam[3];
cd[0] = lamParam[4];
delta[0] = lamParam[5];
ctEq[0] = lamParam[6];
us[0] = lamParam[7];
for (size_t k = 0; k < mesh->getnMarkers(); ++k)
regime[k] = 1;
}
/**
* @brief Compute the blowing velocity in each station of the region.
*
*/
void BoundaryLayer::computeBlwVel()
{
deltaStar[0] = u[0] * u[1];
blowingVelocity.resize(mesh->getnMarkers() - 1, 0.);
for (size_t iPoint = 1; iPoint < mesh->getnMarkers(); ++iPoint)
{
deltaStar[iPoint] = u[iPoint * nVar + 0] * u[iPoint * nVar + 1];
blowingVelocity[iPoint - 1] = (blEdge->getRhoe(iPoint) * blEdge->getVt(iPoint) * deltaStar[iPoint] - blEdge->getRhoe(iPoint - 1) * blEdge->getVt(iPoint - 1) * deltaStar[iPoint - 1]) / (blEdge->getRhoe(iPoint) * (mesh->getLoc(iPoint) - mesh->getLoc(iPoint - 1)));
}
}
/**
* @brief Print solution at a given station.
*
* @param iPoint Marker id.
* @note Function used for debugging.
*/
void BoundaryLayer::printSolution(size_t iPoint) const
{
std::cout << "Pt " << iPoint << "Reg " << regime[iPoint] << std::endl;
std::cout << "x " << mesh->getx(iPoint) << "xx " << mesh->getLoc(iPoint) << "xxExt " << mesh->getExt(iPoint) << std::endl;
std::cout << std::scientific << std::setprecision(15);
std::cout << std::setw(10) << "T " << u[iPoint * nVar + 0] << "\n"
<< std::setw(10) << "H " << u[iPoint * nVar + 1] << "\n"
<< std::setw(10) << "N " << u[iPoint * nVar + 2] << "\n"
<< std::setw(10) << "ue " << u[iPoint * nVar + 3] << "\n"
<< std::setw(10) << "Ct " << u[iPoint * nVar + 4] << "\n"
<< std::setw(10) << "dStar (BL) " << deltaStar[iPoint] << "\n"
<< std::setw(10) << "dStar (old) " << blEdge->getDeltaStar(iPoint) << "\n"
<< std::setw(10) << "vt " << blEdge->getVt(iPoint) << "\n"
<< std::setw(10) << "Me " << blEdge->getMe(iPoint) << "\n"
<< std::setw(10) << "rhoe " << blEdge->getRhoe(iPoint) << std::endl;
}
\ No newline at end of file
vii/src/DBoundaryLayer.h deleted 100644 → 0
View file @ 3a33458d
#ifndef DBOUNDARYLAYER_H
#define DBOUNDARYLAYER_H
#include "vii.h"
#include "DClosures.h"
#include "DDiscretization.h"
#include "DEdge.h"
namespace vii
{
/**
* @brief Boundary layer region upper/lower side or wake.
* @author Paul Dechamps
*/
class VII_API BoundaryLayer
{
private:
unsigned int const nVar = 5; ///< Number of variables of the partial differential problem.
double nCrit = 9.0; ///< Critical amplification factor.
public:
Discretization *mesh; ///< 1D surface boundary layer mesh.
Edge *blEdge; ///< Quantites at the inviscid boundary (edge of the boundary layer).
Closures *closSolver; ///< Closure relations class.
std::string name; ///< Name of the region.
// Boundary layer variables.
std::vector<double> u; ///< Solution vector (theta, H, N, ue, Ct).
std::vector<double> Ret; ///< Reynolds number based on the momentum thickness (theta).
std::vector<double> cd; ///< Local dissipation coefficient.
std::vector<double> cf; ///< Local friction coefficient.
std::vector<double> Hstar; ///< Kinetic energy shape parameter (thetaStar/theta).
std::vector<double> Hstar2; ///< Density shape parameter (deltaStarStar/theta).
std::vector<double> Hk; ///< Kinematic shape parameter (int(1-u/ue d_eta)).
std::vector<double> ctEq; ///< Equilibrium shear stress coefficient (turbulent BL).
std::vector<double> us; ///< Equivalent normalized wall slip velocity.
std::vector<double> delta; ///< Boundary layer thickness.
std::vector<double> deltaStar; ///< Dispacement thickness (int(1-rho*u/rhoe*ue d_eta)).
// Transition related variables.
std::vector<int> regime; ///< Laminar (0) or turbulent (1) regime.
std::vector<double> blowingVelocity; ///< Blowing velocity.
size_t transMarker; ///< Marker id of the transition location.
double xtrF; ///< Forced transition location in the global frame of reference.
double xtr; ///< Transition location in the global frame of reference.
BoundaryLayer(double _xtrF = -1.0);
~BoundaryLayer();
// Boundary conditions.
void setStagBC(double Re);
void setWakeBC(double Re, std::vector<double> upperConditions, std::vector<double> lowerConditions);
void setMesh(std::vector<double> x,
std::vector<double> y,
std::vector<double> z);
// Getters.
size_t getnVar() const { return nVar; };
double getnCrit() const { return nCrit; };
// Others
void printSolution(size_t iPoint) const;
void computeBlwVel();
};
} // namespace vii
#endif // DBOUNDARYLAYER_H
vii/src/DClosures.cpp deleted 100644 → 0
View file @ 3a33458d
#include "DClosures.h"
#include <cmath>
#include <iostream>
using namespace vii;
/**
* @brief Laminar closure relations at a given station.
*
* @param closureVars Vector where the computed variables will be stored.
* @param theta Momentum thickness.
* @param H Shape factor of the boundary layer.
* @param Ret Reynolds number based on the momentum thickness.
* @param Me Local Mach number.
* @param name Name of the region.
*/
void Closures::laminarClosures(std::vector<double> &closureVars, double theta,
double H, double Ret, const double Me,
const std::string name)
{
H = std::max(H, 1.00005);
// Kinematic shape factor H*.
double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
if (name == "wake")
Hk = std::max(Hk, 1.00005);
else
Hk = std::max(Hk, 1.05000);
// Density shape parameter.
double Hstar2 = (0.064 / (Hk - 0.8) + 0.251) * Me * Me;
// Boundary layer thickness.
double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta);
double Hstar = 0.;
double Hks = Hk - 4.35;
if (Hk <= 4.35)
Hstar = 0.0111 * Hks * Hks / (Hk + 1) - 0.0278 * Hks * Hks * Hks / (Hk + 1) + 1.528 - 0.0002 * (Hks * Hk) * (Hks * Hk);
else
Hstar = 1.528 + 0.015 * Hks * Hks / Hk;
// Whitfield's minor additional compressibility correction.
Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me);
// Friction coefficient.
double tmp = 0.;
double cf = 0.;
if (Hk < 5.5)
{
tmp = (5.5 - Hk) * (5.5 - Hk) * (5.5 - Hk) / (Hk + 1.0);
cf = (-0.07 + 0.0727 * tmp) / Ret;
}
else if (Hk >= 5.5)
{
tmp = 1.0 - 1.0 / (Hk - 4.5);
cf = (-0.07 + 0.015 * tmp * tmp) / Ret;
}
// Dissipation coefficient.
double Cd = 0.;
if (Hk < 4)
Cd = (0.00205 * std::pow(4.0 - Hk, 5.5) + 0.207) * (Hstar / (2 * Ret));
else
{
double HkCd = (Hk - 4) * (Hk - 4);
double denCd = 1 + 0.02 * HkCd;
Cd = (-0.0016 * HkCd / denCd + 0.207) * (Hstar / (2 * Ret));
}
// Wake relations.
if (name == "wake")
{
Cd = 2 * (1.10 * (1.0 - 1.0 / Hk) * (1.0 - 1.0 / Hk) / Hk) * (Hstar / (2 * Ret));
cf = 0.;
}
double us = 0.;
double ctEq = 0.;
closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us};
}
/**
* @brief Turbulent closure relations at a given station.
*
* @param closureVars Vector where the computed variables will be stored.
* @param theta Momentum thickness.
* @param H Shape factor of the boundary layer.
* @param Ct Shear stress coefficient.
* @param Ret Reynolds number based on the momentum thickness.
* @param Me Local Mach number.
* @param name Name of the region.
*/
void Closures::turbulentClosures(std::vector<double> &closureVars, double theta, double H, double Ct, double Ret, const double Me, const std::string name)
{
H = std::max(H, 1.00005);
Ct = std::min(Ct, 0.3);
// Ct = std::max(std::min(Ct, 0.30), 0.0000001);
double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
if (name == "wake")
Hk = std::max(Hk, 1.00005);
else
Hk = std::max(Hk, 1.05000);
double Hstar2 = ((0.064 / (Hk - 0.8)) + 0.251) * Me * Me;
double gamma = 1.4 - 1.;
double Fc = std::sqrt(1 + 0.5 * gamma * Me * Me);
double H0 = 0.;
if (Ret <= 400)
H0 = 4.0;
else
H0 = 3 + (400 / Ret);
if (Ret <= 200)
Ret = 200;
double Hstar = 0.;
if (Hk <= H0)
Hstar = ((0.5 - 4 / Ret) * ((H0 - Hk) / (H0 - 1)) * ((H0 - Hk) / (H0 - 1))) * (1.5 / (Hk + 0.5)) + 1.5 + 4 / Ret;
else
Hstar = (Hk - H0) * (Hk - H0) * (0.007 * std::log(Ret) / ((Hk - H0 + 4 / std::log(Ret)) * (Hk - H0 + 4 / std::log(Ret))) + 0.015 / Hk) + 1.5 + 4 / Ret;
// Whitfield's minor additional compressibility correction.
Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me);
double logRt = std::log(Ret / Fc);
logRt = std::max(logRt, 3.0);
double arg = std::max(-1.33 * Hk, -20.0);
// Equivalent normalized wall slip velocity.
double us = (Hstar / 2) * (1 - 4 * (Hk - 1) / (3 * H));
// Boundary layer thickness.
double delta = std::min((3.15 + H + (1.72 / (Hk - 1))) * theta, 12 * theta);
double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
double Hkc = 0.;
double Cdw = 0.;
double Cdd = 0.;
double Cdl = 0.;
double Hmin = 0.;
double Fl = 0.;
double Dfac = 0.;
double cf = 0.;
double Cd = 0.;
double n = 1.0;
double ctEq = 0.;
if (name == "wake")
{
if (us > 0.99995)
us = 0.99995;
n = 2.0; // two wake halves
cf = 0.0; // no friction inside the wake
Hkc = Hk - 1;
Cdw = 0.0; // Wall contribution.
Cdd = (0.995 - us) * Ct * Ct; // Turbulent outer layer contribution.
Cdl = 0.15 * (0.995 - us) * (0.995 - us) / Ret; // Laminar stress contribution to outer layer.
ctEq = std::sqrt(4 * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5.
}
else
{
if (us > 0.95)
us = 0.98;
n = 1.0;
cf = (1 / Fc) * (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) + 0.00011 * (std::tanh(4 - (Hk / 0.875)) - 1));
Hkc = std::max(Hk - 1 - 18 / Ret, 0.01);
// Dissipation coefficient.
Hmin = 1 + 2.1 / std::log(Ret);
Fl = (Hk - 1) / (Hmin - 1);
Dfac = 0.5 + 0.5 * std::tanh(Fl);
Cdw = 0.5 * (cf * us) * Dfac;
Cdd = (0.995 - us) * Ct * Ct;
Cdl = 0.15 * (0.995 - us) * (0.995 - us) / Ret;
ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
}
if (n * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H) < 0)
std::cout << "Negative sqrt encountered " << std::endl;
// Dissipation coefficient.
Cd = n * (Cdw + Cdd + Cdl);
closureVars = {Hstar, Hstar2, Hk, cf, Cd, delta, ctEq, us};
}
/**
* @brief Turbulent closure relations at a given station.
*
* @param closureVars Computed equilibrium shear stress coefficient.
* @param theta Momentum thickness.
* @param H Shape factor of the boundary layer.
* @param Ct Shear stress coefficient.
* @param Ret Reynolds number based on the momentum thickness.
* @param Me Local Mach number.
* @param name Name of the region.
*/
void Closures::turbulentClosures(double &closureVars, double theta, double H, double Ct, double Ret, const double Me, const std::string name)
{
H = std::max(H, 1.00005);
Ct = std::min(Ct, 0.3);
double Hk = (H - 0.29 * Me * Me) / (1 + 0.113 * Me * Me);
if (name == "wake")
Hk = std::max(Hk, 1.00005);
else
Hk = std::max(Hk, 1.05000);
double H0 = 0.;
if (Ret <= 400)
H0 = 4.0;
else
H0 = 3 + (400 / Ret);
if (Ret <= 200)
Ret = 200;
double Hstar = 0.;
if (Hk <= H0)
Hstar = ((0.5 - 4 / Ret) * ((H0 - Hk) / (H0 - 1)) * ((H0 - Hk) / (H0 - 1))) * (1.5 / (Hk + 0.5)) + 1.5 + 4 / Ret;
else
Hstar = (Hk - H0) * (Hk - H0) * (0.007 * std::log(Ret) / ((Hk - H0 + 4 / std::log(Ret)) * (Hk - H0 + 4 / std::log(Ret))) + 0.015 / Hk) + 1.5 + 4 / Ret;
// Whitfield's minor additional compressibility correction.
Hstar = (Hstar + 0.028 * Me * Me) / (1 + 0.014 * Me * Me);
// Equivalent normalized wall slip velocity.
double us = (Hstar / 2) * (1 - 4 * (Hk - 1) / (3 * H));
double Ctcon = 0.5 / (6.7 * 6.7 * 0.75);
double Hkc = 0.;
double ctEq = 0.;
if (name == "wake")
{
if (us > 0.99995)
us = 0.99995;
Hkc = Hk - 1;
ctEq = std::sqrt(4 * Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
}
else
{
if (us > 0.95)
us = 0.98;
Hkc = std::max(Hk - 1 - 18 / Ret, 0.01);
ctEq = std::sqrt(Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H)); // Here it is ctEq^0.5
}
if (Hstar * Ctcon * ((Hk - 1) * Hkc * Hkc) / ((1 - us) * (Hk * Hk) * H) < 0)
std::cout << "Negative sqrt encountered " << std::endl;
closureVars = ctEq;
}
\ No newline at end of file
vii/src/DClosures.h deleted 100644 → 0
View file @ 3a33458d
#ifndef DCLOSURES_H
#define DCLOSURES_H
#include "vii.h"
#include <vector>
#include <string>
namespace vii
{
/**
* @brief Boundary layer closure relations.
* @author Paul Dechamps
*/
class VII_API Closures
{
public:
static void laminarClosures(std::vector<double> &closureVars, double theta,
double H, double Ret, const double Me,
const std::string name);
static void turbulentClosures(std::vector<double> &closureVars, double theta,
double H, double Ct, double Ret, double Me,
const std::string name);
static void turbulentClosures(double &closureVars, double theta,
double H, double Ct, double Ret, const double Me,
const std::string name);
};
} // namespace vii
#endif // DCLOSURES_H
vii/src/DDiscretization.cpp deleted 100644 → 0
View file @ 3a33458d
#include "DDiscretization.h"
#include <cmath>
using namespace vii;
Discretization::Discretization()
{
nMarkers = 0;
nElm = 0;
prevnMarkers = 0;
}
Discretization::~Discretization() {}
void Discretization::setGlob(std::vector<double> &_x, std::vector<double> &_y, std::vector<double> &_z)
{
if (_x.size() != _y.size() || _y.size() != _z.size() || _x.size() != _z.size())
throw std::runtime_error("vii::Discretization Wrong mesh sizes.\n");
nMarkers = _x.size();
x.resize(nMarkers, 0.);
y.resize(nMarkers, 0.);
z.resize(nMarkers, 0.);
nElm = nMarkers - 1;
x = _x;
y = _y;
z = _z;
computeBLcoord();
}
/**
* @brief Compute nodes coordinates in the local frame of reference.
*/
void Discretization::computeBLcoord()
{
loc.resize(nMarkers, 0.);
alpha.resize(nElm, 0.);
dx.resize(nElm, 0.);
for (size_t iPoint = 0; iPoint < nElm; ++iPoint)
{
alpha[iPoint] = std::atan2(y[iPoint + 1] - y[iPoint], x[iPoint + 1] - x[iPoint]);
// dx = sqrt(delta x ^2 + delta y ^2).
dx[iPoint] = std::sqrt((x[iPoint + 1] - x[iPoint]) * (x[iPoint + 1] - x[iPoint]) + (y[iPoint + 1] - y[iPoint]) * (y[iPoint + 1] - y[iPoint]));
loc[iPoint + 1] = loc[iPoint] + dx[iPoint];
}
}
void Discretization::setExt(std::vector<double> extM)
{
ext.resize(extM.size(), 0.);
ext = extM;
};
vii/src/DDiscretization.h deleted 100644 → 0
View file @ 3a33458d
#ifndef DDISCRETIZATION_H
#define DDISCRETIZATION_H
#include "vii.h"
#include "DBoundaryLayer.h"
namespace vii
{
/**
* @brief Boundary layer mesh.
* @author Paul Dechamps
*/
class VII_API Discretization
{
private:
std::vector<double> x; ///< x coordinate in the global frame of reference.
std::vector<double> y; ///< y coordinate in the global frame of reference.
std::vector<double> z; ///< z coordinate in the global frame of reference.
std::vector<double> loc; ///< xi coordinate in the local frame of reference.
std::vector<double> ext; ///< xi coordinate in the local frame of reference, at the edge of the boundary layer.
std::vector<double> dx; ///< Cell size.
std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference.
size_t nMarkers; ///< Number of points in the domain.
size_t nElm; ///< Number of cells in the domain.
public:
size_t prevnMarkers; ///< Number of points in the domain at the previous iteration.
Discretization();
~Discretization();
// Getters.
size_t getnMarkers() const { return nMarkers; };
double getLoc(size_t iPoint) const { return loc[iPoint]; };
double getx(size_t iPoint) const { return x[iPoint]; };
double gety(size_t iPoint) const { return y[iPoint]; };
double getz(size_t iPoint) const { return z[iPoint]; };
double getExt(size_t iPoint) const { return ext[iPoint]; };
double getAlpha(size_t iElm) const { return alpha[iElm]; };
size_t getDiffnElm() const { return nMarkers - prevnMarkers; };
// Setters.
void setGlob(std::vector<double> &_x, std::vector<double> &_y, std::vector<double> &_z);
void setExt(std::vector<double> extM);
void computeBLcoord();
};
} // namespace vii
#endif // WDISCRETIZATION_H