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
dart/viscous/newton.py deleted 100644 → 0
View file @ 7b79ac85
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
## Newton raphson method coupled with linear solver
# Amaury Bilocq
import numpy as np
from fwk.coloring import ccolors
def newton(f, x, maxIt, tol=1.0e-8):
'''Newton procedure to solve the non linear set of boundary layer equations.
Boundary layer equations are parabolic in x -> downstream marching is used
An implicit marching model is applied (Crank Nicolson) -> matrix inversion is performed
'''
# Numerical Jacobian
def jaco(f, x):
eps = 1.0e-8
n = len(x)
jac = np.zeros((n,n))
xp = x.copy()
xm = x.copy()
for j in range(n):
# perturb solution
xTmp = x[j]
xp[j] += eps
xm[j] -= eps
# compute jacobian
jac[:,j] = (f(xp) - f(xm)) / (2 * eps)
# reset
xp[j] = xTmp
xm[j] = xTmp
return jac
# Backtrack line search
def ls(f, x, dx):
a = 1.0
ires = np.linalg.norm(f(x))
while a > 1/64:
res = np.linalg.norm(f(x + a * dx))
if res < ires:
break
else:
a /= 2
return a
# Solve
for i in range(maxIt):
# solve
jac = jaco(f,x)
dx = np.linalg.solve(jac, -f(x))
# line search
a = ls(f, x, dx)
x += a * dx
fx = f(x)
err = np.linalg.norm(fx)
# safeguard
x[0] = max(x[0],0.)
x[1] = max(x[1],1.0005)
if i != 0 and err < tol:
return x
# Raise error but return solution and error
raise RuntimeError(ccolors.ANSI_YELLOW + 'Newton solver - too many iterations!\n' + ccolors.ANSI_RESET, x, err)
dart/viscous/solver.py deleted 100644 → 0
View file @ 7b79ac85
#!/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.
## Integral boundary layer equations viscous solver
# Amaury Bilocq
#
# TODO list
# - Improve user-input
# -* let the user define the transition location
# -* let the user define the turbulence intensity (which will affect the amplification factor)
# - Improve the physics (robustness)
# -* wake should start turbulent
# -* use empirical correlations to fix Cteq at transition
# -- split the panel where the transition happens instead of averaging
# -* understand why Ct becomes negative in turbulent solver and prevent it
# - Improve the numerics (robustness)
# -* test the convergence by monitoring the change in the interaction velocity (vt) between 2 VVI iterations
# -- implement a more sophisticated line search (3 points quadratic?)
# -* try to use underrelaxation
# -- compute the jacobian analytically
# -* use an inverse coupling procedure when laminar separation occurs
# - Increase range of validity
# -- split the inviscid and viscous mesh
# -- check transonic predictions
# -- implement quasi-2D method for 3D flows
from dart.viscous.boundary import Boundary
from dart.viscous.wake import Wake
import dart.viscous.newton as Newton
import numpy as np
from fwk.coloring import ccolors
class Solver:
def __init__(self, _Re):
'''Initialize the viscous solver
'''
self.Re = _Re # Reynolds number
self.gamma = 1.4 # heat capacity of air
self.Cd = 0 # drag coefficient
self.it = 0 # viscous inviscid iteration for interaction law
self.xtr = [1, 1] # set transition at TE reference coordinates
def dictionary(self):
'''Create dictionary with the coordinates and the boundary layer parameters
'''
wData = {}
wData['x'] = self.data[:,0]
wData['y'] = self.data[:,1]
wData['z'] = self.data[:,2]
wData['x/c'] = self.data[:,0]/max(self.data[:,0])
wData['Cd'] = self.blScalars[0]
wData['Cdf'] = self.blScalars[1]
wData['Cdp'] = self.blScalars[2]
wData['delta*'] = self.blVectors[0]
wData['theta'] = self.blVectors[1]
wData['H'] = self.blVectors[2]
wData['Hk'] = self.blVectors[3]
wData['H*'] = self.blVectors[4]
wData['H**'] = self.blVectors[5]
wData['Cf'] = self.blVectors[6]
wData['Cdissip'] = self.blVectors[7]
wData['Ctau'] = self.blVectors[8]
wData['delta'] = self.blVectors[9]
wData['xtrTop'] = self.xtr[0]
wData['xtrBot'] = self.xtr[1]
return wData
def writeFile(self):
'''Write the results in airfoil_viscous.dat
'''
wData = self.dictionary()
toW = ['delta*', 'H', 'Cf'] # list of variable to be written (should be a user input)
# Write
print('writing file: airfoil_viscous.dat...', end = '')
f = open('airfoil_viscous.dat', 'w+')
f.write('$Aerodynamic coefficients\n')
f.write(' Re Cd Cdp Cdf xtr_top xtr_bot\n')
f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f} {4:15.6f} {5:15.6f}\n'.format(self.Re, wData['Cd'], wData['Cdp'], wData['Cdf'], wData['xtrTop'], wData['xtrBot']))
f.write('$Boundary layer variables\n')
f.write(' x y z x/c')
for s in toW:
f.write(' {0:>15s}'.format(s))
f.write('\n')
for i in range(len(self.data[:,0])):
f.write('{0:15.6f} {1:15.6f} {2:15.6f} {3:15.6f}'.format(wData['x'][i], wData['y'][i], wData['z'][i], wData['x/c'][i]))
for s in toW:
f.write(' {0:15.6f}'.format(wData[s][i]))
f.write('\n')
f.close()
print('done!')
def __getBLcoordinates(self, data):
'''Transform the reference coordinates into airfoil coordinates
'''
nN = len(data[:,0])
nE = nN-1 #nbr of element along the surface
vt = np.zeros(nN)
xx = np.zeros(nN) # position along the surface of the airfoil
dx = np.zeros(nE) # distance along two nodes
dv = np.zeros(nE) # speed difference according to element
alpha = np.zeros(nE) # angle of the element for the rotation matrix
for i in range(0,nE):
alpha[i] = np.arctan2((data[i+1,1]-data[i,1]),(data[i+1,0]-data[i,0]))
dx[i] = np.sqrt((data[i+1,0]-data[i,0])**2+(data[i+1,1]-data[i,1])**2)
xx[i+1] = dx[i]+xx[i]
vt[i] = (data[i,3]*np.cos(alpha[i]) + data[i,4]*np.sin(alpha[i])) # tangential speed at node 1 element i
vt[i+1] = (data[i+1,3]*np.cos(alpha[i]) + data[i+1,4]*np.sin(alpha[i])) # tangential speed at node 2 element i
dv[i] = (vt[i+1] - vt[i])/dx[i] #velocity gradient according to element i. central scheme with half point
return xx, dv, vt, nN, alpha
def __amplRoutine(self, Hk, Ret, theta):
'''Compute the amplitude of the TS waves envelop for the transition
'''
Dgr = 0.08
Hk1 = 3.5
Hk2 = 4
Hmi = (1/(Hk-1))
logRet_crit = 2.492*(1/(Hk-1))**0.43 + 0.7*(np.tanh(14*Hmi-9.24)+1.0) # value at which the amplification starts to grow
if Ret <=0:
Ret = 1
if np.log10(Ret) < logRet_crit - Dgr :
Ax = 0
else:
Rnorm = (np.log10(Ret) - (logRet_crit-Dgr))/(2*Dgr)
if Rnorm <=0:
Rfac = 0
if Rnorm >= 1:
Rfac = 1
else:
Rfac = 3*Rnorm**2 - 2*Rnorm**3
# normal envelope amplification rate
m = -0.05+2.7*Hmi-5.5*Hmi**2+3*Hmi**3+0.1*np.exp(-20*Hmi)
arg = 3.87*Hmi-2.52
dndRet = 0.028*(Hk-1)-0.0345*np.exp(-(arg**2))
Ax = (m*dndRet/theta)*Rfac
# set correction for separated profiles
if Hk > Hk1:
Hnorm = (Hk - Hk1)/(Hk2-Hk1)
if Hnorm >=1:
Hfac = 1
else:
Hfac = 3*Hnorm**2 - 2*Hnorm**3
Ax1 = Ax
Gro = 0.3+0.35*np.exp(-0.15*(Hk-5))
Tnr = np.tanh(1.2*(np.log10(Ret)-Gro))
Ax2 = (0.086*Tnr - 0.25/(Hk-1)**1.5)/theta
if Ax2 < 0:
Ax2 = 0
Ax = Hfac*Ax2 + (1-Hfac)*Ax1
if Ax < 0:
Ax = 0
return Ax
def __laminarClosure(self, theta, H, Ret,Me, name):
''' Laminar closure derived from the Falkner-Skan one-parameter profile family
Nishida and Drela (1996)
'''
Hk = (H - 0.29*Me**2)/(1+0.113*Me**2) # Kinematic shape factor
if name == "airfoil":
Hk = max(Hk, 1.02)
Hk = min(Hk,6)
else:
Hk = max(Hk, 1.00005)
Hk = min(Hk,6)
Hstar2 = (0.064/(Hk-0.8)+0.251)*Me**2 # Density shape parameter
delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta)
Hks = Hk - 4.35
if Hk <= 4.35:
dens = Hk + 1
Hstar = 0.0111*Hks**2/dens - 0.0278*Hks**3/dens + 1.528 -0.0002*(Hks*Hk)**2
elif Hk > 4.35:
Hstar = 1.528 + 0.015*Hks**2/Hk
Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction
if Hk < 5.5:
tmp = (5.5-Hk)**3 / (Hk+1)
Cf = (-0.07 + 0.0727*tmp)/Ret
elif Hk >= 5.5:
tmp = 1 - 1/(Hk-4.5)
Cf = (-0.07 + 0.015*tmp**2) /Ret
if Hk < 4:
Cd = ( 0.00205*(4.0-Hk)**5.5 + 0.207 ) * (Hstar/(2*Ret))
elif Hk >= 4:
HkCd = (Hk-4)**2
denCd = 1+0.02*HkCd
Cd = ( -0.0016*HkCd/denCd + 0.207) * (Hstar/(2*Ret))
if name == 'wake':
Cd = 2*(1.10*(1.0-1.0/Hk)**2/Hk)* (Hstar/(2*Ret))
Cf = 0
return Cd, Cf, Hstar, Hstar2, Hk, delta
def __turbulentClosure(self, theta, H, Ct, Ret, Me, name):
''' Turbulent closure derived from Nishida and Drela (1996)
'''
# eliminate absurd transients
Ct = min(Ct, 0.30)
Ct = max(Ct, 0.0000001)
Hk = (H - 0.29*Me**2)/(1+0.113*Me**2)
if name == 'wake':
Hk = max(Hk, 1.00005)
Hk = min(Hk,10)
else:
Hk = max(Hk, 1.00005)
Hk = min(Hk,10)
Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me**2
gam = self.gamma - 1
Fc = np.sqrt(1+0.5*gam*Me**2)
if Ret <= 400:
H0 = 4
elif Ret > 400:
H0 = 3 + (400/Ret)
if Ret > 200:
Ret = Ret
elif Ret <= 200:
Ret = 200
if Hk <= H0:
Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))**2)*(1.5/(Hk+0.5))+1.5+4/Ret
elif Hk > H0:
Hstar = (Hk-H0)**2*(0.007*np.log(Ret)/(Hk-H0+4/np.log(Ret))**2 + 0.015/Hk) + 1.5 + 4/Ret
Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction
logRt = np.log(Ret/Fc)
logRt = np.max((logRt,3))
arg = np.max((-1.33*Hk, -20))
Us = (Hstar/2)*(1-4*(Hk-1)/(3*H)) # Equivalent normalized wall slip velocity
delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta)
Ctcon = 0.5/(6.7**2*0.75)
if name == 'wake':
if Us > 0.99995:
Us = 0.99995
n = 2 # two wake halves
Cf = 0 # no friction inside the wake
Hkc = Hk - 1
Cdw = 0 # wall contribution
Cdd = (0.995-Us)*Ct**2 # turbulent outer layer contribution
Cdl = 0.15*(0.995-Us)**2/Ret # laminar stress contribution to outer layer
else:
if Us > 0.95:
Us = 0.98
n = 1
Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt/2.3026)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-1))
Hkc = max(Hk-1-18/Ret, 0.01)
# dissipation coefficient
Hmin = 1 + 2.1/np.log(Ret)
Fl = (Hk-1)/(Hmin-1)
Dfac = 0.5+0.5*np.tanh(Fl)
Cdw = 0.5*(Cf*Us) * Dfac
Cdd = (0.995-Us)*Ct**2
Cdl = 0.15*(0.995-Us)**2/Ret
Cteq = np.sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)) #Here it is Cteq^0.5
Cd = n*(Cdw+ Cdd + Cdl)
Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2)
return Cd, Cf, Hstar, Hstar2, Hk, Cteq, delta
def __boundaryLayer(self, data, group, savef=False):
'''Discretize and solve the boundary layer equations for laminar/turbulent flows
'''
xx, dv, vtOld, nN, alpha = self.__getBLcoordinates(data)
Me = data[:,5]
rhoe = data[:,6]
deltaStarOld = data[:,7]
if self.it == 0:
xxOld = xx
else:
xxOld = data[:,8]
#Initialize variables
blwVel = np.zeros(nN-1) # Blowing velocity
vt = np.zeros(nN) # Boundary layer velocity
theta = np.zeros(nN) # Momentum thickness
H = np.zeros(nN) # Shape factor
deltaStar = np.zeros(nN) # Displacement thickness
Hstar = np.zeros(nN) # Kinetic energy shape parameter
Hstar2 = np.zeros(nN) # Density shape parameter
Hk = np.zeros(nN) # Kinematic shape factor
Ret = np.zeros(nN) #Re based on momentum thickness
Cd = np.zeros(nN) # Dissipation coefficient
Cf = np.zeros(nN) # Skin friction
n = np.zeros(nN) # Logarithm of the maximum amplification ratio
Ax = np.zeros(nN) # Amplification factor
Ct = np.zeros(nN) # Shear stress coefficient
Cteq = np.zeros(nN) # Equilibrium shear stress coefficient
delta = np.zeros(nN) # BL thickness
# Boundary conditions
if group.name == 'airfoil':
theta[0], H[0], n[0], Ct[0] = group.initialConditions(self.Re, dv[0])
elif group.name == 'wake':
theta[0] = group.T0
H[0] = group.H0
n[0] = group.n0
Ct[0] = group.Ct0
vt[0] = vtOld[0]
Ret[0] = vt[0]*theta[0]*self.Re
if Ret[0] < 1:
Ret[0] = 1
vt[0] = Ret[0]/(theta[0]*self.Re)
deltaStar[0] = theta[0]*H[0]
# Define transition location ( N = 9) and compute stagnation point
ntr = 9
Cd[0], Cf[0], Hstar[0], Hstar2[0], Hk[0], delta[0] = self.__laminarClosure( theta[0], H[0], Ret[0],Me[0], group.name)
if n[0] < ntr:
turb = 0 # we begin with laminar flow
else :
turb = 1 # typically for the wake
xtr = 0 # if only turbulent flow
itTurb = 0
# Solve the boundary layer equations
for i in range(1,nN):
# x[0] = theta; x[1] = H; x[2] = n for laminar/Ct for turbulent; x[3] = vt from interaction law
def f_lam(x):
f = np.zeros(len(x))
Ret[i] = np.maximum(x[3]*x[0]*self.Re, 1)
Cd[i], Cf[i], Hstar[i], Hstar2[i], Hk[i], delta[i] = self.__laminarClosure(x[0], x[1], Ret[i],Me[i], group.name)
Ta = upw*x[0]+(1-upw)*theta[i-1]
Ha = upw*x[1]+(1-upw)*H[i-1]
uea = upw*x[3]+(1-upw)*vt[i-1]
Reta = upw*Ret[i] + (1-upw)*Ret[i-1]
Mea = upw*Me[i]+(1-upw)*Me[i-1]
Cda, Cfa, Hstara, Hstar2a, Hka, deltaa = self.__laminarClosure(Ta , Ha, Reta, Mea, group.name)
dx = xx[i] - xx[i-1]
Axa = self.__amplRoutine(Hka, Reta, Ta)
dTheta = x[0]-theta[i-1]
due = x[3]-vt[i-1]
dH = x[1]-H[i-1]
dn = x[2] - n[i-1]
dHstar = Hstar[i]-Hstar[i-1]
Hstara = upw*Hstar[i]+(1-upw)*Hstar[i-1]
f[0] = dTheta/dx + (2 + Ha - Mea**2)*(Ta/uea)*due/dx - Cfa/2
f[1] = Ta*(dHstar/dx) + (2*Hstar2a + Hstara*(1-Ha))*Ta/uea * due/dx - 2*Cda + Hstara*Cfa/2
f[2] = dn - dx*Axa
f[3] = uea - 4/(np.pi*dx)*Ta*Ha - (upw*vtOld[i]+(1-upw)*vtOld[i-1]) + 4/(np.pi*(xxOld[i]-xxOld[i-1]))*(upw*deltaStarOld[i]+(1-upw)*deltaStarOld[i-1])
return f
def f_turb(x):
f = np.zeros(len(x))
if group.name == 'wake':
dissipRatio = 0.9 # wall/wake dissipation length ratio
else:
dissipRatio = 1
Ret[i] = x[3]*x[0]*self.Re
Ta = upw*x[0]+(1-upw)*theta[i-1]
xxa = upw*xx[i]+(1-upw)*xx[i-1]
Ha = upw*x[1]+(1-upw)*H[i-1]
Cta = upw*x[2]+(1-upw)*Ct[i-1]
uea = upw*x[3]+(1-upw)*vt[i-1]
Reta = np.maximum(upw*(x[3]*x[0]*self.Re)+(1-upw)*(vt[i-1]*theta[i-1]*self.Re), 1)
Mea = upw*Me[i]+(1-upw)*Me[i-1]
Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i] = self.__turbulentClosure(x[0], x[1],x[2], Ret[i],Me[i], group.name)
Cda, Cfa, Hstara, Hstar2a, Hka,Cteqa, deltaa = self.__turbulentClosure(Ta, Ha,Cta, Reta,Mea, group.name)
dx = xx[i]-xx[i-1]
dTheta = x[0]-theta[i-1]
due = x[3]-vt[i-1]
dH = x[1]-H[i-1]
dHstar = Hstar[i]-Hstar[i-1]
dCt = x[2]-Ct[i-1] # here it is Ct^1/2 which is represented
Ta = upw*x[0]+(1-upw)*theta[i-1]
Cta = upw*x[2]+(1-upw)*Ct[i-1]
f[0] = dTheta/dx + (2 + Ha - Mea**2)*(Ta/uea)*due/dx - Cfa/2
f[1] = Ta*(dHstar/dx) + (2*Hstar2a + Hstara*(1-Ha))*Ta/uea * due/dx - 2*Cda + Hstara*Cfa/2
f[2] = (2*deltaa/Cta)*(dCt/dx) - 5.6*(Cteqa-Cta*dissipRatio)-2*deltaa*(4/(3*Ha*Ta)*(Cfa/2-((Hka-1)/(6.7*Hka*dissipRatio))**2)-1/uea * due/dx)
f[3] = uea - 4/(np.pi*dx)*Ta*Ha - (upw*vtOld[i]+(1-upw)*vtOld[i-1]) + 4/(np.pi*(xxOld[i]-xxOld[i-1]))*(upw*deltaStarOld[i]+(1-upw)*deltaStarOld[i-1])
return f
# Laminar solver
if turb == 0:
if n[i-1] > 8 or i < 5:
upw = 1 # backward Euler
else:
upw = 0.5 # Trapezoidal scheme
x = np.array([theta[i-1],H[i-1], n[i-1], vt[i-1]])
try:
sol = Newton.newton(f_lam, x, 200)
except RuntimeError as e:
sol = e.args[1]
print(str(e.args[0]) + ccolors.ANSI_YELLOW + 'Laminar solver not converged at position: {:.3f}, error: {:.1e}'.format(xx[i], e.args[2]) + ccolors.ANSI_RESET)
# Determine if the amplification waves start to grow or not
logRet_crit = 2.492*(1/(Hk[i]-1))**0.43 + 0.7*(np.tanh(14*(1/(Hk[i]-1))-9.24)+1.0)
if np.log10(Ret[i]) <= logRet_crit - 0.08 :
n[i] = 0
else:
n[i] = sol[2]
xtr = 1
# Turbulent solver
elif turb == 1:
if i <2 and group.name =='wake':
upw = 1 # begin of the wake
elif itTurb <2 and group.name == 'airfoil':
upw = 1
itTurb += 1
else:
upw = 0.5 # trapezoidal scheme
x = np.array([theta[i-1], H[i-1], Ct[i-1], vt[i-1]])
try:
sol = Newton.newton(f_turb, x, 200)
except RuntimeError as e:
sol = e.args[1]
print(str(e.args[0]) + ccolors.ANSI_YELLOW + 'Turbulent solver not converged at position: {:.3f}, error: {:.1e}'.format(xx[i], e.args[2]) + ccolors.ANSI_RESET)
Ct[i] = sol[2]
theta[i] = sol[0]
H[i] = sol[1]
vt[i] = sol[3]
# Transition solver
if n[i] >= ntr and turb ==0:
turb = 1
upw = 1
xtr = (ntr-n[i-1])*((data[i,0]-data[i-1,0])/(n[i]-n[i-1]))+data[i-1,0]
avTurb = (data[i,0]-xtr)/(data[i,0]-data[i-1,0]) # percentage of the subsection corresponding to a turbulent flow
avLam = 1- avTurb # percentage of the subsection corresponding to a laminar flow
# Averaging both solutions
Cteq[i-1]= self.__turbulentClosure(theta[i-1], H[i-1], 0.03, Ret[i-1], Me[i-1], group.name)[5]
Ct[i-1] = avLam*Cteq[i-1] * 0.7
x_turb = np.array([theta[i-1], 1.515, Ct[i-1] , vt[i-1]])
try:
sol_turb = Newton.newton(f_turb, x_turb, 200)
except RuntimeError as e:
sol = e.args[1]
print(str(e.args[0]) + ccolors.ANSI_YELLOW + 'Transition solver not converged at position: {:.3f}, error: {:.1e}'.format(xx[i], e.args[2]) + ccolors.ANSI_RESET)
theta[i] = avLam*sol[0]+avTurb*sol_turb[0]
H[i] = avLam*sol[1]+avTurb*sol_turb[1]
if group.name == 'airfoil':
Ct[i] = max(avTurb*sol_turb[2],0.03)
else:
Ct[i] = max(avTurb*sol_turb[2],0.05)
vt[i] = avLam*sol[3]+avTurb*sol_turb[3]
Cd[i-1], Cf[i-1], Hstar[i-1], Hstar2[i-1], Hk[i-1], delta[i-1] = self.__laminarClosure(theta[i-1], H[i-1], Ret[i-1],Me[i-1], group.name)
Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i] = self.__turbulentClosure(theta[i], H[i],Ct[i], Ret[i],Me[i], group.name)
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]*(xx[i]-xx[i-1]))
# Normalize the friction and dissipation coefficients by the freestream velocity
Cf = vt**2*Cf
Cd = vt**2*Cd
# Compute the various drag coefficients (total, friction, profile)
Cdtot = 2*theta[-1]*(vt[-1]**((H[-1]+5)/2))
Cdf = np.trapz(Cf, data[:,0])
Cdp = Cdtot-Cdf
# Outputs
blScalars = [Cdtot, Cdf, Cdp]
blVectors = [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta]
return blwVel, xtr, xx, blScalars, blVectors
def run(self, group):
''' Run the viscous solver to solve the BL equations and give the outputs to the coupler
'''
group.connectListNodes, group.connectListElems, data = group.connectList()
# Compute BLE for the airfoil (upper section(data[0]) and lower section (data[1]))
if len(data) == 2:
ublwVel, xtrTop, uxx, ublScalars, ublVectors = self.__boundaryLayer(data[0], group, True)
lblwVel, xtrBot, lxx, lblScalars, lblVectors = self.__boundaryLayer(data[1], group)
# Put the upper and lower values together
self.data = np.concatenate((data[0], data[1]))
self.blScalars = np.add(ublScalars, lblScalars)
self.blVectors = np.hstack((ublVectors, lblVectors))
blwVel = np.concatenate((ublwVel, lblwVel))
self.xtr = [xtrTop, xtrBot]
# Boundary conditions for the wake
group.TEnd = [ublVectors[1][-1], lblVectors[1][-1]]
group.HEnd = [ublVectors[2][-1], lblVectors[2][-1]]
group.CtEnd = [ublVectors[8][-1], lblVectors[8][-1]]
# Delete the stagnation point doublon for variables in the VII loop
lblVectors[0] = np.delete(lblVectors[0],0,0) #deltastar
lxx = np.delete(lxx,0,0) # airfoil coordinates
group.deltaStar = np.concatenate((ublVectors[0], lblVectors[0]))
group.xx = np.concatenate((uxx, lxx))
# Compute BLE for the wake
else:
blwVel, _, group.xx, blScalars, blVectors = self.__boundaryLayer(data, group)
group.deltaStar = blVectors[0]
# The drag is measured at the end of the wake
self.Cd = blScalars[0]
self.Cdp = self.Cd-self.blScalars[1] # Cdf comes from the airfoil as there is no friction in the wake
self.blScalars[0] = self.Cd
self.blScalars[2] = self.Cdp
self.Cdf = self.blScalars[1]
# Sort the following results in reference frame
group.deltaStar = group.deltaStar[group.connectListNodes.argsort()]
group.xx = group.xx[group.connectListNodes.argsort()]
group.u = blwVel[group.connectListElems.argsort()]
dart/viscous/wake.py deleted 100644 → 0
View file @ 7b79ac85
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# Copyright 2020 University of Liège
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
## Wake behind airfoil (around which the boundary layer is computed)
# Amaury Bilocq
from dart.viscous.boundary import Boundary
from dart.viscous.airfoil import Airfoil
import numpy as np
class Wake(Boundary):
def __init__(self, _boundary):
Boundary.__init__(self, _boundary)
self.name = 'wake' # type of boundary
self.T0 = 0 # initial condition for the momentum thickness
self.H0 = 0
self.n0 = 9 # wake is always turbulent
self.Ct0 = 0
def connectList(self):
''' Sort the value read by the viscous solver/ Create list of connectivity
'''
N1 = np.zeros(self.nN, dtype=int) # Node number
connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
data = np.zeros((self.boundary.nodes.size(), 10))
i = 0
for n in self.boundary.nodes:
data[i,0] = n.no
data[i,1] = n.pos[0]
data[i,2] = n.pos[1]
data[i,3] = n.pos[2]
data[i,4] = self.v[i,0]
data[i,5] = self.v[i,1]
data[i,6] = self.Me[i]
data[i,7] = self.rhoe[i]
data[i,8] = self.deltaStar[i]
data[i,9] = self.xx[i]
i += 1
# Table containing the element and its nodes
eData = np.zeros((self.nE,3), dtype=int)
for i in range(0, self.nE):
eData[i,0] = self.boundary.tag.elems[i].no
eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
connectListNodes = data[:,1].argsort()
connectListElems[0] = np.where(eData[:,1] == min(data[:,1]))[0]
for i in range(1, len(eData[:,0])):
connectListElems[i] = np.where(eData[:,1] == eData[connectListElems[i-1],2])[0]
data[:,0] = data[connectListNodes,0]
data[:,1] = data[connectListNodes,1]
data[:,2] = data[connectListNodes,2]
data[:,3] = data[connectListNodes,3]
data[:,4] = data[connectListNodes,4]
data[:,5] = data[connectListNodes,5]
data[:,6] = data[connectListNodes,6]
data[:,7] = data[connectListNodes,7]
data[:,8] = data[connectListNodes,8]
data[:,9] = data[connectListNodes,9]
# Separated upper and lower part
data = np.delete(data,0,1)
return connectListNodes, connectListElems, data
amfe @ c11ff6ee
Compare 87212cfb...c11ff6ee
Subproject commit 87212cfb034a92cc83acfeea751f076a84c12779
Subproject commit c11ff6ee646765a19f9d43496d354e939a84b408
scripts/format_code.py deleted 100755 → 0
View file @ 7b79ac85
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
#
# This script runs "clang-format" recursively on all C/C++ files
import sys
import os
import fnmatch
import re
import subprocess
def all_files(root,
patterns='*',
skips='*.git*;*build*',
single_level=False,
yield_folders=False):
# self.checkPath(root)
patterns = patterns.split(';')
skips = skips.split(';')
for path, subdirs, files in os.walk(root):
if yield_folders:
files.extend(subdirs)
files.sort()
for name in files:
for pattern in patterns:
if fnmatch.fnmatch(name, pattern):
fullname = os.path.join(path, name)
ok = True
for skip in skips:
if fnmatch.fnmatch(os.path.relpath(fullname, start=root), skip):
ok = False
if ok:
yield fullname
break
if single_level:
break
def main():
# loop over all files and format them
encs = {}
for f in all_files(os.getcwd(), patterns='*.cpp;*.c;*.h;*.hpp'):
# print(f)
cmd = ['clang-format-10', "-style=file", "-i", f]
retcode = subprocess.call(cmd)
if retcode != 0:
print(f'ERROR: retcode = {retcode}')
break
if __name__ == "__main__":
# print('running format_code.py...')
main()