From 44c0abd084b0f52978125e3d76a43d3532f3b656 Mon Sep 17 00:00:00 2001
From: Bilocq Amaury <amaury.bilocq@student.uliege.be>
Date: Wed, 11 Nov 2020 13:38:57 +0100
Subject: [PATCH] Squashed - Quasi-simultaneous coupler for VII with test cases
---
.gitlab-ci.yml | 2 +-
flow/tests/bli.py | 30 ++++++++++++------------
flow/viscous/airfoil.py | 32 +++++++++++++++-----------
flow/viscous/coupler.py | 2 +-
flow/viscous/solver.py | 51 ++++++++++++++++-------------------------
flow/viscous/wake.py | 4 ++--
6 files changed, 57 insertions(+), 64 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 451ca6ee..f4cf7d46 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -10,7 +10,7 @@ default:
.global_tag: &global_tag_def
tags:
- - mn2l
+# - mn2l
# - warson # you can choose a set of runners here
stages:
diff --git a/flow/tests/bli.py b/flow/tests/bli.py
index 3a89ec2e..c6e8cbc2 100644
--- a/flow/tests/bli.py
+++ b/flow/tests/bli.py
@@ -21,13 +21,13 @@
# Amaury Bilocq
# Test the viscous-inviscid interaction scheme
# Reference to the master's thesis: http://hdl.handle.net/2268/252195
-# Reference test cases with Naca0012:
+# Reference test cases with Naca0012 (different from master's thesis):
# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, p = 2, m = 7, tol = 10^-4, msTE = 0.01, msLE = 0.001
-# -> nIt = 41, Cl = 0.58, Cd = 0.0062, xtrTop = 0.0555, xtrBot = 0.7397
+# -> nIt = 32, Cl = 0.55, Cd = 0.0061, xtrTop = 0.0611, xtrBot = 0.7392
# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, p = 2, m = 7, tol = 10^-4, msTE = 0.01, msLE = 0.001
-# -> nIt = , Cl = 0.69, Cd = 0.0067, xtrTop = 0.0384, xtrBot = 0.7364
-# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, p = 2, m = 11, tol = 5.2*10^-3, msTE = 0.01, msLE = 0.00075
-# -> nIt = , Cl = 1.39 , Cd = 0.011, xtrTop = 0.008, xtrBot = 1
+# -> nIt = 45, Cl = 0.64, Cd = 0.0065, xtrTop = 0.0497, xtrBot = 0.7338
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, p = 2, m = 5, tol = 2.5*10^-3, msTE = 0.01, msLE = 0.00075
+# -> nIt = 39, Cl = 1.30 , Cd = 0.011, xtrTop = 0.010, xtrBot = 1
#
# CAUTION
# This test is provided to ensure that the solver works properly.
@@ -60,7 +60,7 @@ def main():
# define filter parameters and tolerance of the VII
p = 2
- m = 7
+ m = 5
tol = 1e-4
# mesh the geometry
@@ -75,7 +75,7 @@ def main():
pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, M_crit, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True)
tms['pre'].stop()
- # solve inviscid problem
+ # solve viscous problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
isolver = floD.newton(pbl)
@@ -106,22 +106,22 @@ def main():
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
if Re == 1e7 and M_inf == 0 and alpha == 5*math.pi/180:
- tests.add(CTest('Cl', isolver.Cl, 0.58, 5e-2))
+ tests.add(CTest('Cl', isolver.Cl, 0.55, 5e-2)) # Xfoil 0.58
tests.add(CTest('Cd', vsolver.Cd, 0.0062, 0.01))
tests.add(CTest('Cdp', vsolver.Cdp, 0.0018, 0.01))
- tests.add(CTest('xtr_top', vsolver.xtr[0], 0.056, 0.05))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.061, 0.05)) # Xfoil 0.056
tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.740, 0.05))
elif Re == 1e7 and M_inf == 0.5 and alpha == 5*math.pi/180:
- tests.add(CTest('Cl', isolver.Cl, 0.69, 5e-2))
+ tests.add(CTest('Cl', isolver.Cl, 0.65, 5e-2)) # Xfoil 0.69
tests.add(CTest('Cd', vsolver.Cd, 0.0067, 0.01))
tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 0.01))
- tests.add(CTest('xtr_top', vsolver.xtr[0], 0.038, 0.05))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.049, 0.05)) # Xfoil 0.038
tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.736, 0.05))
elif Re == 1e7 and M_inf == 0 and alpha == 12*math.pi/180:
- tests.add(CTest('Cl', isolver.Cl, 1.39, 5e-2))
+ tests.add(CTest('Cl', isolver.Cl, 1.30, 5e-2)) # Xfoil 1.39
tests.add(CTest('Cd', vsolver.Cd, 0.0011, 0.01))
- tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 0.01))
- tests.add(CTest('xtr_top', vsolver.xtr[0], 0.008, 0.05))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0064, 0.01))
+ tests.add(CTest('xtr_top', vsolver.xtr[0], 0.010, 0.05)) # Xfoil 0.008
tests.add(CTest('xtr_bot', vsolver.xtr[1], 1.000, 0.05))
else:
raise Exception('Test not defined for this flow')
@@ -132,4 +132,4 @@ def main():
print('')
if __name__ == "__main__":
- main()
+ main()
\ No newline at end of file
diff --git a/flow/viscous/airfoil.py b/flow/viscous/airfoil.py
index 0dd0ccd1..15fc43c7 100644
--- a/flow/viscous/airfoil.py
+++ b/flow/viscous/airfoil.py
@@ -72,30 +72,34 @@ class Airfoil(Boundary):
eData[i,0] = self.boundary.tag.elems[i].no
eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
- # Defining TE/LE nodes number
+ # Find the stagnation point
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]
+ # 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:
- 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
+ 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
+ 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
+ # 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):
+ if eData[connectListElems[i],2] == upperGlobTE:
upperTE = 1
i += 1
# Sort the pressure side
@@ -106,7 +110,7 @@ class Airfoil(Boundary):
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):
+ if eData[connectListElems[i+1],1] == lowerGlobTE:
lowerTE = 1
i += 1
connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
@@ -125,4 +129,4 @@ class Airfoil(Boundary):
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]
+ return connectListNodes, connectListElems,[uData, lData]
\ No newline at end of file
diff --git a/flow/viscous/coupler.py b/flow/viscous/coupler.py
index 727a9a73..4a672c20 100644
--- a/flow/viscous/coupler.py
+++ b/flow/viscous/coupler.py
@@ -82,4 +82,4 @@ class Coupler:
# save results
self.isolver.save(0, self.writer)
self.vsolver.writeFile()
- print('\n')
+ print('\n')
\ No newline at end of file
diff --git a/flow/viscous/solver.py b/flow/viscous/solver.py
index 6af0f757..54f97a45 100644
--- a/flow/viscous/solver.py
+++ b/flow/viscous/solver.py
@@ -29,7 +29,6 @@
# - Implement quasi-2D method for 3D capability
# - User can define the transition location (l 303)
# - Write the empirical relation between the amplification factor and the turbulence intensity (l 303)
-# - Check why the computation crash if the first point of the wake is computed with the turbulent closure (l 309)
from flow.viscous.boundary import Boundary
from flow.viscous.wake import Wake
@@ -69,6 +68,7 @@ class Solver:
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
@@ -193,8 +193,8 @@ class Solver:
Cd = ( -0.0016*HkCd/denCd + 0.207) * (Hstar/(2*Ret))
if name == 'wake':
Cd = 2*(1.10*(1.0-1.0/Hk)**2/Hk)* (Hstar/(2*Ret))
- Cf = 0
- return Cd, Cf, Hstar, Hstar2, Hk, delta, H
+ 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)
@@ -256,7 +256,7 @@ class Solver:
Cteq = np.sqrt(n*Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*(Hk**2)*H)) #Here it is Cteq^0.5
Cd = n*(Cdw+ Cdd + Cdl)
Ueq = 1.25/Hk*(Cf/2-((Hk-1)/(6.432*Hk))**2)
- return Cd, Cf, Hstar, Hstar2, Hk, Cteq, delta, Ueq, Ct
+ return Cd, Cf, Hstar, Hstar2, Hk, Cteq, delta
def __boundaryLayer(self, data, group):
@@ -269,6 +269,8 @@ class Solver:
#Filter the initial conditions (Me and rho can be filtered too)
if group.name == 'airfoil':
vtOld = savgol_filter(vtOld, self.m, self.p, mode = 'interp')
+ Me = savgol_filter(Me, self.m, self.p, mode = 'interp')
+ rhoe = savgol_filter(rhoe, self.m, self.p, mode = 'interp')
vtInv = vtOld
rhoInv = rhoe
if self.it == 0:
@@ -291,7 +293,6 @@ class Solver:
Ax = np.zeros(nN) # Amplification factor
Ct = np.zeros(nN) # Shear stress coefficient
Cteq = np.zeros(nN) # Equilibrium shear stress coefficient
- Ueq = np.zeros(nN) # Equilibrium velocity gradient
delta = np.zeros(nN) # BL thickness
# Boundary conditions
if group.name == 'airfoil':
@@ -309,27 +310,26 @@ class Solver:
deltaStar[0] = theta[0]*H[0]
# Define transition location ( N = 9) and compute stagnation point
ntr = 9
- if n[0] < 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
- Cd[0], Cf[0], Hstar[0], Hstar2[0], Hk[0], delta[0], H[0] = self.__laminarClosure( theta[0], H[0], Ret[0],Me[0], group.name)
else :
turb = 1 # typically for the wake
- Cd[0], Cf[0], Hstar[0], Hstar2[0], Hk[0], delta[0], H[0] = self.__laminarClosure( theta[0], H[0], Ret[0],Me[0], group.name)
xtr = 0 # if only turbulent flow
itTurb = 0
# Solve the boundary layer equations
for i in range(1,nN):
# x[0] = theta; x[1] = H; x[2] = n for laminar/Ct for turbulent; x[3] = vt from interaction law
def f_lam(x):
- f = np.zeros(len(x))
+ f = np.zeros(len(x))
Ret[i] = np.maximum(x[3]*x[0]*self.Re, 1)
- Cd[i], Cf[i], Hstar[i], Hstar2[i], Hk[i], delta[i], x[1] = self.__laminarClosure(x[0], x[1], Ret[i],Me[i], group.name)
+ 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, Ha = self.__laminarClosure(Ta , Ha, Reta, Mea, group.name)
+ 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]
@@ -357,8 +357,8 @@ class Solver:
uea = upw*x[3]+(1-upw)*vt[i-1]
Reta = np.maximum(upw*(x[3]*x[0]*self.Re)+(1-upw)*(vt[i-1]*theta[i-1]*self.Re), 1)
Mea = upw*Me[i]+(1-upw)*Me[i-1]
- Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i], Ueq[i], x[2] = self.__turbulentClosure(x[0], x[1],x[2], Ret[i],Me[i], group.name)
- Cda, Cfa, Hstara, Hstar2a, Hka,Cteqa, deltaa, Ueqa, Cta = self.__turbulentClosure(Ta, Ha,Cta, Reta,Mea, group.name)
+ 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]
@@ -401,7 +401,7 @@ class Solver:
Ct[i] = sol[2]
theta[i] = sol[0]
H[i] = sol[1]
- vt[i] = sol[3]
+ vt[i] = sol[3]
# Transition solver
if n[i] >= ntr and turb ==0:
turb = 1
@@ -421,10 +421,10 @@ class Solver:
else:
Ct[i] = max(avTurb*sol_turb[2],0.05)
vt[i] = avLam*sol[3]+avTurb*sol_turb[3]
- Cd[i-1], Cf[i-1], Hstar[i-1], Hstar2[i-1], Hk[i-1], delta[i-1], _ = self.__laminarClosure(theta[i-1], H[i-1], Ret[i-1],Me[i-1], group.name)
- Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i], Ueq[i], _ = self.__turbulentClosure(theta[i], H[i],Ct[i], Ret[i],Me[i], group.name)
+ 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] = -1*(rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1])) # TO DO: need to divide by rhoe or not?
+ blwVel[i-1] = -1*(rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*(xx[i]-xx[i-1]))
# Normalize the friction and dissipation coefficients by the freestream velocity
Cf = vt**2*Cf
Cd = vt**2*Cd
@@ -434,7 +434,7 @@ class Solver:
Cdp = Cdtot-Cdf
# Outputs
blScalars = [Cdtot, Cdf, Cdp]
- blVectors = [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct]
+ blVectors = [deltaStar, theta, H, Hk, Hstar, Hstar2, Cf, Cd, Ct, delta]
return blwVel, xtr, xx, blScalars, blVectors
@@ -451,7 +451,7 @@ class Solver:
self.blScalars = np.add(ublScalars, lblScalars)
self.blVectors = np.hstack((ublVectors, lblVectors))
blwVel = np.concatenate((ublwVel, lblwVel))
- blwVel = savgol_filter(blwVel, self.m, self.p, mode = 'interp')
+ blwVel = savgol_filter(blwVel, self.m, self.p, mode = 'interp')
self.xtr = [xtrTop, xtrBot]
# Boundary conditions for the wake
group.TEnd = [ublVectors[1][-1], lblVectors[1][-1]]
@@ -475,15 +475,4 @@ class Solver:
# 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]
-
-
-
-
-
-
-
-
-
-
-
+ group.u = blwVel[group.connectListElems.argsort()]
\ No newline at end of file
diff --git a/flow/viscous/wake.py b/flow/viscous/wake.py
index 0e8da24a..0f3cd5e5 100644
--- a/flow/viscous/wake.py
+++ b/flow/viscous/wake.py
@@ -61,7 +61,7 @@ class Wake(Boundary):
eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
connectListNodes = data[:,1].argsort()
- connectListElems[0] = np.where(eData[:,1] == 1.0)[0]
+ 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]
@@ -76,4 +76,4 @@ class Wake(Boundary):
data[:,9] = data[connectListNodes,9]
# Separated upper and lower part
data = np.delete(data,0,1)
- return connectListNodes, connectListElems, data
+ return connectListNodes, connectListElems, data
\ No newline at end of file
--
GitLab