diff --git a/flow/models/n0012.db b/flow/models/n0012.db
new file mode 100644
index 0000000000000000000000000000000000000000..b8f485733b337b3280214d824d6f6b73da657805
Binary files /dev/null and b/flow/models/n0012.db differ
diff --git a/flow/models/n0012s.geo b/flow/models/n0012s.geo
new file mode 100644
index 0000000000000000000000000000000000000000..d42b35316e2b1c2431b8678081465b0f9bdca1b3
--- /dev/null
+++ b/flow/models/n0012s.geo
@@ -0,0 +1,298 @@
+/* Naca0012 */
+
+// Geometry
+DefineConstant[ xLgt = { 5.0, Name "Domain length (x-dir)" } ];
+DefineConstant[ yLgt = { 5.0, Name "Domain length (y-dir)" } ];
+
+// Mesh
+nChord = 200;
+pChord = 0.3;
+nLe = 115;
+pLe = 1.015;
+nPerp = 25;
+pPerp = 1.25;
+
+/**************
+ Geometry
+ **************/
+
+//Point(201) = {1.000000,0.000000,0};
+Point(200) = {0.999753,-0.000035,0};
+Point(199) = {0.999013,-0.000141,0};
+Point(198) = {0.997781,-0.000317,0};
+Point(197) = {0.996057,-0.000562,0};
+Point(196) = {0.993844,-0.000876,0};
+Point(195) = {0.991144,-0.001258,0};
+Point(194) = {0.987958,-0.001707,0};
+Point(193) = {0.984292,-0.002222,0};
+Point(192) = {0.980147,-0.002801,0};
+Point(191) = {0.975528,-0.003443,0};
+Point(190) = {0.970440,-0.004147,0};
+Point(189) = {0.964888,-0.004909,0};
+Point(188) = {0.958877,-0.005729,0};
+Point(187) = {0.952414,-0.006603,0};
+Point(186) = {0.945503,-0.007531,0};
+Point(185) = {0.938153,-0.008510,0};
+Point(184) = {0.930371,-0.009537,0};
+Point(183) = {0.922164,-0.010610,0};
+Point(182) = {0.913540,-0.011726,0};
+Point(181) = {0.904508,-0.012883,0};
+Point(180) = {0.895078,-0.014079,0};
+Point(179) = {0.885257,-0.015310,0};
+Point(178) = {0.875056,-0.016574,0};
+Point(177) = {0.864484,-0.017868,0};
+Point(176) = {0.853553,-0.019189,0};
+Point(175) = {0.842274,-0.020535,0};
+Point(174) = {0.830656,-0.021904,0};
+Point(173) = {0.818712,-0.023291,0};
+Point(172) = {0.806454,-0.024694,0};
+Point(171) = {0.793893,-0.026111,0};
+Point(170) = {0.781042,-0.027539,0};
+Point(169) = {0.767913,-0.028974,0};
+Point(168) = {0.754521,-0.030414,0};
+Point(167) = {0.740877,-0.031856,0};
+Point(166) = {0.726995,-0.033296,0};
+Point(165) = {0.712890,-0.034733,0};
+Point(164) = {0.698574,-0.036163,0};
+Point(163) = {0.684062,-0.037582,0};
+Point(162) = {0.669369,-0.038988,0};
+Point(161) = {0.654508,-0.040378,0};
+Point(160) = {0.639496,-0.041747,0};
+Point(159) = {0.624345,-0.043094,0};
+Point(158) = {0.609072,-0.044414,0};
+Point(157) = {0.593691,-0.045705,0};
+Point(156) = {0.578217,-0.046962,0};
+Point(155) = {0.562667,-0.048182,0};
+Point(154) = {0.547054,-0.049362,0};
+Point(153) = {0.531395,-0.050499,0};
+Point(152) = {0.515705,-0.051587,0};
+Point(151) = {0.500000,-0.052625,0};
+Point(150) = {0.484295,-0.053608,0};
+Point(149) = {0.468605,-0.054534,0};
+Point(148) = {0.452946,-0.055397,0};
+Point(147) = {0.437333,-0.056195,0};
+Point(146) = {0.421783,-0.056924,0};
+Point(145) = {0.406309,-0.057581,0};
+Point(144) = {0.390928,-0.058163,0};
+Point(143) = {0.375655,-0.058666,0};
+Point(142) = {0.360504,-0.059087,0};
+Point(141) = {0.345492,-0.059424,0};
+Point(140) = {0.330631,-0.059674,0};
+Point(139) = {0.315938,-0.059834,0};
+Point(138) = {0.301426,-0.059902,0};
+Point(137) = {0.287110,-0.059876,0};
+Point(136) = {0.273005,-0.059754,0};
+Point(135) = {0.259123,-0.059535,0};
+Point(134) = {0.245479,-0.059217,0};
+Point(133) = {0.232087,-0.058799,0};
+Point(132) = {0.218958,-0.058280,0};
+Point(131) = {0.206107,-0.057661,0};
+Point(130) = {0.193546,-0.056940,0};
+Point(129) = {0.181288,-0.056119,0};
+Point(128) = {0.169344,-0.055197,0};
+Point(127) = {0.157726,-0.054176,0};
+Point(126) = {0.146447,-0.053056,0};
+Point(125) = {0.135516,-0.051839,0};
+Point(124) = {0.124944,-0.050527,0};
+Point(123) = {0.114743,-0.049121,0};
+Point(122) = {0.104922,-0.047624,0};
+Point(121) = {0.095492,-0.046037,0};
+Point(120) = {0.086460,-0.044364,0};
+Point(119) = {0.077836,-0.042608,0};
+Point(118) = {0.069629,-0.040770,0};
+Point(117) = {0.061847,-0.038854,0};
+Point(116) = {0.054497,-0.036863,0};
+Point(115) = {0.047586,-0.034800,0};
+Point(114) = {0.041123,-0.032668,0};
+Point(113) = {0.035112,-0.030471,0};
+Point(112) = {0.029560,-0.028212,0};
+Point(111) = {0.024472,-0.025893,0};
+Point(110) = {0.019853,-0.023517,0};
+Point(109) = {0.015708,-0.021088,0};
+Point(108) = {0.012042,-0.018607,0};
+Point(107) = {0.008856,-0.016078,0};
+Point(106) = {0.006156,-0.013503,0};
+Point(105) = {0.003943,-0.010884,0};
+Point(104) = {0.002219,-0.008223,0};
+Point(103) = {0.000987,-0.005521,0};
+Point(102) = {0.000247,-0.002779,0};
+Point(101) = {0.000000,0.000000,0};
+Point(100) = {0.000247,0.002779,0};
+Point(99) = {0.000987,0.005521,0};
+Point(98) = {0.002219,0.008223,0};
+Point(97) = {0.003943,0.010884,0};
+Point(96) = {0.006156,0.013503,0};
+Point(95) = {0.008856,0.016078,0};
+Point(94) = {0.012042,0.018607,0};
+Point(93) = {0.015708,0.021088,0};
+Point(92) = {0.019853,0.023517,0};
+Point(91) = {0.024472,0.025893,0};
+Point(90) = {0.029560,0.028212,0};
+Point(89) = {0.035112,0.030471,0};
+Point(88) = {0.041123,0.032668,0};
+Point(87) = {0.047586,0.034800,0};
+Point(86) = {0.054497,0.036863,0};
+Point(85) = {0.061847,0.038854,0};
+Point(84) = {0.069629,0.040770,0};
+Point(83) = {0.077836,0.042608,0};
+Point(82) = {0.086460,0.044364,0};
+Point(81) = {0.095492,0.046037,0};
+Point(80) = {0.104922,0.047624,0};
+Point(79) = {0.114743,0.049121,0};
+Point(78) = {0.124944,0.050527,0};
+Point(77) = {0.135516,0.051839,0};
+Point(76) = {0.146447,0.053056,0};
+Point(75) = {0.157726,0.054176,0};
+Point(74) = {0.169344,0.055197,0};
+Point(73) = {0.181288,0.056119,0};
+Point(72) = {0.193546,0.056940,0};
+Point(71) = {0.206107,0.057661,0};
+Point(70) = {0.218958,0.058280,0};
+Point(69) = {0.232087,0.058799,0};
+Point(68) = {0.245479,0.059217,0};
+Point(67) = {0.259123,0.059535,0};
+Point(66) = {0.273005,0.059754,0};
+Point(65) = {0.287110,0.059876,0};
+Point(64) = {0.301426,0.059902,0};
+Point(63) = {0.315938,0.059834,0};
+Point(62) = {0.330631,0.059674,0};
+Point(61) = {0.345492,0.059424,0};
+Point(60) = {0.360504,0.059087,0};
+Point(59) = {0.375655,0.058666,0};
+Point(58) = {0.390928,0.058163,0};
+Point(57) = {0.406309,0.057581,0};
+Point(56) = {0.421783,0.056924,0};
+Point(55) = {0.437333,0.056195,0};
+Point(54) = {0.452946,0.055397,0};
+Point(53) = {0.468605,0.054534,0};
+Point(52) = {0.484295,0.053608,0};
+Point(51) = {0.500000,0.052625,0};
+Point(50) = {0.515705,0.051587,0};
+Point(49) = {0.531395,0.050499,0};
+Point(48) = {0.547054,0.049362,0};
+Point(47) = {0.562667,0.048182,0};
+Point(46) = {0.578217,0.046962,0};
+Point(45) = {0.593691,0.045705,0};
+Point(44) = {0.609072,0.044414,0};
+Point(43) = {0.624345,0.043094,0};
+Point(42) = {0.639496,0.041747,0};
+Point(41) = {0.654508,0.040378,0};
+Point(40) = {0.669369,0.038988,0};
+Point(39) = {0.684062,0.037582,0};
+Point(38) = {0.698574,0.036163,0};
+Point(37) = {0.712890,0.034733,0};
+Point(36) = {0.726995,0.033296,0};
+Point(35) = {0.740877,0.031856,0};
+Point(34) = {0.754521,0.030414,0};
+Point(33) = {0.767913,0.028974,0};
+Point(32) = {0.781042,0.027539,0};
+Point(31) = {0.793893,0.026111,0};
+Point(30) = {0.806454,0.024694,0};
+Point(29) = {0.818712,0.023291,0};
+Point(28) = {0.830656,0.021904,0};
+Point(27) = {0.842274,0.020535,0};
+Point(26) = {0.853553,0.019189,0};
+Point(25) = {0.864484,0.017868,0};
+Point(24) = {0.875056,0.016574,0};
+Point(23) = {0.885257,0.015310,0};
+Point(22) = {0.895078,0.014079,0};
+Point(21) = {0.904508,0.012883,0};
+Point(20) = {0.913540,0.011726,0};
+Point(19) = {0.922164,0.010610,0};
+Point(18) = {0.930371,0.009537,0};
+Point(17) = {0.938153,0.008510,0};
+Point(16) = {0.945503,0.007531,0};
+Point(15) = {0.952414,0.006603,0};
+Point(14) = {0.958877,0.005729,0};
+Point(13) = {0.964888,0.004909,0};
+Point(12) = {0.970440,0.004147,0};
+Point(11) = {0.975528,0.003443,0};
+Point(10) = {0.980147,0.002801,0};
+Point(9) = {0.984292,0.002222,0};
+Point(8) = {0.987958,0.001707,0};
+Point(7) = {0.991144,0.001258,0};
+Point(6) = {0.993844,0.000876,0};
+Point(5) = {0.996057,0.000562,0};
+Point(4) = {0.997781,0.000317,0};
+Point(3) = {0.999013,0.000141,0};
+Point(2) = {0.999753,0.000035,0};
+Point(1) = {1.000000,0.000000,0};
+
+Spline(1) = {80,79,78,77,76,75,74,73,72,71,70,69,68,67,66,65,64,63,62,61,60,59,58,57,56,55,54,53,52,51,50,49,48,47,46,45,44,43,42,41,40,39,38,37,36,35,34,33,32,31,30,29,28,27,26,25,24,23,22,21,20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1};
+Spline(2) = {101,100,99,98,97,96,95,94,93,92,91,90,89,88,87,86,85,84,83,82,81,80};
+Spline(3) = {122,121,120,119,118,117,116,115,114,113,112,111,110,109,108,107,106,105,104,103,102,101};
+Spline(4) = {1,200,199,198,197,196,195,194,193,192,191,190,189,188,187,186,185,184,183,182,181,180,179,178,177,176,175,174,173,172,171,170,169,168,167,166,165,164,163,162,161,160,159,158,157,156,155,154,153,152,151,150,149,148,147,146,145,144,143,142,141,140,139,138,137,136,135,134,133,132,131,130,129,128,127,126,125,124,123,122};
+
+// Farfield
+Point(10001) = {1+xLgt, 0, 0};
+Point(10002) = {1+xLgt, yLgt, 0};
+Point(10003) = {1, yLgt, 0};
+Point(10004) = {-xLgt, yLgt, 0};
+Point(10005) = {-xLgt, 0, 0};
+Point(10006) = {-xLgt,-yLgt, 0};
+Point(10007) = {1, -yLgt, 0};
+Point(10008) = {1+xLgt, -yLgt, 0};
+
+Line(10001) = {10001, 10002};
+Line(10002) = {10002, 10003};
+Line(10003) = {10003, 10004};
+Line(10004) = {10004, 10005};
+Line(10005) = {10005, 10006};
+Line(10006) = {10006, 10007};
+Line(10007) = {10007, 10008};
+Line(10008) = {10008, 10001};
+
+// Perpendicular lines
+Line(10011) = {1, 10001};
+Line(10012) = {1, 10003};
+Line(10013) = {80, 10004};
+Line(10014) = {101, 10005};
+Line(10015) = {122, 10006};
+Line(10016) = {1, 10007};
+
+// Internal field
+Line Loop(10001) = {10011, 10001, 10002, -10012};
+Line Loop(10002) = {1, 10012, 10003, -10013};
+Line Loop(10003) = {2, 10013, 10004, -10014};
+Line Loop(10004) = {3, 10014, 10005, -10015};
+Line Loop(10005) = {4, 10015, 10006, -10016};
+Line Loop(10006) = {-10011, 10016, 10007, 10008};
+
+Plane Surface(10001) = {10001};
+Plane Surface(10002) = {10002};
+Plane Surface(10003) = {10003};
+Plane Surface(10004) = {10004};
+Plane Surface(10005) = {10005};
+Plane Surface(10006) = {10006};
+
+/*************************
+ Meshing
+ *************************/
+
+Transfinite Line {1, 4, 10003, 10006} = nChord Using Bump pChord;
+Transfinite Line {2, 10004} = nLe Using Progression 1/pLe;
+Transfinite Line {3, 10005} = nLe Using Progression pLe;
+Transfinite Line {10011, 10012, 10013, 10014, 10015, 10016, 10001, 10007} = nPerp Using Progression pPerp;
+Transfinite Line {10002, 10008} = nPerp Using Progression 1/pPerp;
+
+Transfinite Surface {10001} = {1, 10001, 10002, 10003};
+Transfinite Surface {10002} = {1, 10003, 10004, 80};
+Transfinite Surface {10003} = {80, 10004, 10005, 101};
+Transfinite Surface {10004} = {122, 10006, 10005, 101};
+Transfinite Surface {10005} = {1, 10007, 10006, 122};
+Transfinite Surface {10006} = {1, 10007, 10008, 10001};
+//Recombine Surface {10001, 10002, 10003, 10004, 10005, 10006};
+
+/*************************
+ Physical Groups
+ *************************/
+
+Physical Point("te") = {1};
+Physical Line("upstream") = {10004, 10005};
+Physical Line("side") = {10002, 10003, 10006, 10007};
+Physical Line("downstream") = {10001, 10008};
+Physical Line("airfoil") = {1, 2};
+Physical Line("airfoil_") = {3, 4};
+Physical Line("wake") = {10011};
+Physical Surface("field") = {10001, 10002, 10003, 10004, 10005, 10006};
+
diff --git a/flow/models/rae2822.db b/flow/models/rae2822.db
new file mode 100644
index 0000000000000000000000000000000000000000..a340d2d615e7884045fd00448803c3550a5e15bd
Binary files /dev/null and b/flow/models/rae2822.db differ
diff --git a/flow/tests/bliLiftComp.py b/flow/tests/bliLiftComp.py
index b924e10179e743eaa0803439cab1f66884d2dae3..20d8e2e3dc4e1ede70084c5eef5275662b8d5d64 100644
--- a/flow/tests/bliLiftComp.py
+++ b/flow/tests/bliLiftComp.py
@@ -19,7 +19,14 @@
# @brief Compute lifting (linear or nonlinear) viscous flow around a naca0012
# @authors Adrien Crovato, 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:
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, p = 2, m = 7, tol = 10^-4, msTE = 0.01, msLE = 0.001
+# Results: Cl = 0.58, Cd = 0.0062, xtrTop = 0.0555, xtrBot = 0.7397
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, p = 2, m = 7, tol = 10^-4, msTE = 0.01, msLE = 0.001
+# Results: 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.15*10^-3, msTE = 0.01, msLE = 0.00075
+# Results: Cl = 1.39 , Cd = 0.011, xtrTop = 0.008, xtrBot = 1
# CAUTION
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
@@ -29,8 +36,8 @@ import math
import flow as flo
import flow.utils as floU
import flow.default as floD
-import flow.viscous.solver as floVS
-import flow.viscous.coupler as floVC
+import flow.viscous.Solver as floVS
+import flow.viscous.Coupler as floVC
import tbox
import tbox.utils as tboxU
import tbox.gmsh as gmsh
@@ -45,13 +52,18 @@ def main():
# define flow variables
Re = 10*10**6
- alpha = 2*math.pi/180
+ alpha = 5*math.pi/180
U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
- M_inf = 0.6
+ M_inf = 0.5
M_crit = 2 # Squared critical Mach number (above which density is modified)
c_ref = 1
dim = len(U_inf)
+ # define filter parameters and tolerance of the VII
+ p = 2
+ m = 7
+ tol = 10**-4
+
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
@@ -93,8 +105,8 @@ def main():
tms['solver'].start()
isolver = flo.Newton(pbl)
floD.initNewton(isolver)
- vsolver = floVS.Solver(Re)
- coupler = floVC.Coupler(isolver, vsolver, gmshWriter,blw,blw2)
+ vsolver = floVS.Solver(Re, p, m)
+ coupler = floVC.Coupler(isolver, vsolver, gmshWriter,blw,blw2, tol)
coupler.run()
tms['solver'].stop()
@@ -119,12 +131,11 @@ def main():
# check results
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.0558, 1e-1)) # TODO check value and tolerance
- tests.add(CTest('Cl', isolver.Cl, 0.292, 5e-2))
- tests.add(CTest('Cdp', vsolver.Cdp, 0.0013, 0.01))
- tests.add(CTest('Cd', vsolver.Cd, 0.0057, 0.01))
- tests.add(CTest('top_xtr', vsolver.top_xtr, 0.16, 0.5 ))
- tests.add(CTest('bot_xtr', vsolver.bot_xtr, 0.48, 0.5 ))
+ tests.add(CTest('Cl', isolver.Cl, 0.68, 5e-2))
+ tests.add(CTest('Cd', vsolver.Cd, 0.0067, 0.01))
+ tests.add(CTest('Cdp', vsolver.Cdp, 0.0025, 0.01))
+ tests.add(CTest('top_xtr', vsolver.xtr[0], 0.0374, 0.05 ))
+ tests.add(CTest('bot_xtr', vsolver.xtr[1], 0.7391, 0.05 ))
tests.run()
# eof
diff --git a/flow/tests/bliLiftIncomp.py b/flow/tests/bliLiftIncomp.py
deleted file mode 100644
index 635cc62a246c378bd221fe8508a49d8358b21d0a..0000000000000000000000000000000000000000
--- a/flow/tests/bliLiftIncomp.py
+++ /dev/null
@@ -1,134 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-# @brief Compute lifting (linear or nonlinear) viscous flow around a naca0012
-# @authors Adrien Crovato, Amaury Bilocq
-# Test the viscous-inviscid interaction scheme
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-from __future__ import print_function
-import math
-import flow as flo
-import flow.utils as floU
-import flow.default as floD
-import flow.viscous.solver as floVS
-import flow.viscous.coupler as floVC
-import tbox
-import tbox.utils as tboxU
-import tbox.gmsh as gmsh
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-def main():
- # timer
- tms = fwk.Timers()
- tms['total'].start()
-
- # define flow variables
- Re = 10*10**6
- alpha = 5*math.pi/180
- U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
- M_inf = 0.0
- M_crit = 5 # Squared critical Mach number (above which density is modified)
- c_ref = 1
- dim = len(U_inf)
-
- # mesh the geometry
- print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
- tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
- msh = gmsh.MeshLoader("../models/n0012.geo",__file__).execute(**pars)
- # crack the mesh
- mshCrck = tbox.MshCrack(msh, dim)
- mshCrck.setCrack('wake')
- mshCrck.addBoundaries(['field', 'airfoil', 'downstream'])
- mshCrck.run()
- del mshCrck
- gmshWriter = tbox.GmshExport(msh)
- gmshWriter.save(msh.name)
- tms['msh'].stop()
-
- # set the problem and add medium, initial/boundary conditions
- tms['pre'].start()
- pbl = flo.Problem(msh, dim, alpha, M_inf, c_ref, c_ref, 0.25, 0., 0.)
- if M_inf == 0:
- pbl.set(flo.Medium(msh, 'field', flo.F0ElRhoL(), flo.F0ElMachL(), flo.F0ElCpL(), flo.F0PsPhiInf(dim, alpha)))
- else:
- pbl.set(flo.Medium(msh, 'field', flo.F0ElRho(M_inf, M_crit), flo.F0ElMach(M_inf), flo.F0ElCp(M_inf), flo.F0PsPhiInf(dim, alpha)))
- pbl.add(flo.Initial(msh, 'field', flo.F0PsPhiInf(dim, alpha)))
- pbl.add(flo.Dirichlet(msh, 'upstream', flo.F0PsPhiInf(dim, alpha)))
- pbl.add(flo.Freestream(msh, 'side', tbox.Fct1C(-U_inf[0], -U_inf[1], 0.)))
- pbl.add(flo.Freestream(msh, 'downstream', tbox.Fct1C(-U_inf[0], -U_inf[1], 0.)))
- pbl.add(flo.Wake(msh, ['wake', 'wake_', 'field']))
- pbl.add(flo.Kutta(msh, ['te', 'wake_', 'airfoil', 'field']))
- bnd = flo.Boundary(msh, ['airfoil', 'field'])
- pbl.add(bnd)
- blw = flo.Blowing(msh, ['airfoil', 'field'])
- blw2 = flo.Blowing(msh, ['wake', 'field'])
- pbl.add(blw)
- pbl.add(blw2)
- tms['pre'].stop()
-
- # solve inviscid problem
- print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
- tms['solver'].start()
- isolver = flo.Newton(pbl)
- floD.initNewton(isolver)
- vsolver = floVS.Solver(Re)
- coupler = floVC.Coupler(isolver, vsolver, gmshWriter,blw,blw2)
- coupler.run()
- tms['solver'].stop()
-
- # extract Cp
- Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
- tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
- # display results
- print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
- print(' M alpha Cl Cd Cm')
- print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, isolver.Cm))
-
- # display timers
- tms['total'].stop()
- print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
- print('CPU statistics')
- print(tms)
-
- # visualize solution and plot results
- floD.initViewer(pbl)
- tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
-
- # check results
- print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- tests = CTests()
- tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.9531, 1e-1)) # TODO check value and tolerance
- tests.add(CTest('Cl', isolver.Cl, 0.5658, 5e-2))
- tests.add(CTest('Cdp', vsolver.Cdp, 0.0017, 0.01))
- tests.add(CTest('Cd', vsolver.Cd, 0.0061, 0.01))
- tests.add(CTest('top_xtr', vsolver.top_xtr, 0.0531, 0.5 ))
- tests.add(CTest('bot_xtr', vsolver.bot_xtr, 0.7484, 0.5 ))
- tests.run()
-
- # eof
- print('')
-
-if __name__ == "__main__":
- main()
diff --git a/flow/tests/bliNonLift.py b/flow/tests/bliNonLift.py
deleted file mode 100644
index d1f4ef78781a29457645418bbd38e7a2a63d1cf7..0000000000000000000000000000000000000000
--- a/flow/tests/bliNonLift.py
+++ /dev/null
@@ -1,134 +0,0 @@
-#!/usr/bin/env python
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-# http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-
-# @brief Compute lifting (linear or nonlinear) viscous flow around a naca0012
-# @authors Adrien Crovato, Amaury Bilocq
-# Test the viscous-inviscid interaction scheme
-#
-# CAUTION
-# This test is provided to ensure that the solver works properly.
-# Mesh refinement may have to be performed to obtain physical results.
-
-from __future__ import print_function
-import math
-import flow as flo
-import flow.utils as floU
-import flow.default as floD
-import flow.viscous.solver as floVS
-import flow.viscous.coupler as floVC
-import tbox
-import tbox.utils as tboxU
-import tbox.gmsh as gmsh
-import fwk
-from fwk.testing import *
-from fwk.coloring import ccolors
-
-def main():
- # timer
- tms = fwk.Timers()
- tms['total'].start()
-
- # define flow variables
- Re = 5*10**6
- alpha = 0*math.pi/180
- U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
- M_inf = 0.0
- M_crit = 5 # Squared critical Mach number (above which density is modified)
- c_ref = 1
- dim = len(U_inf)
-
- # mesh the geometry
- print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
- tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
- msh = gmsh.MeshLoader("../models/n0012.geo",__file__).execute(**pars)
- # crack the mesh
- mshCrck = tbox.MshCrack(msh, dim)
- mshCrck.setCrack('wake')
- mshCrck.addBoundaries(['field', 'airfoil', 'downstream'])
- mshCrck.run()
- del mshCrck
- gmshWriter = tbox.GmshExport(msh)
- gmshWriter.save(msh.name)
- tms['msh'].stop()
-
- # set the problem and add medium, initial/boundary conditions
- tms['pre'].start()
- pbl = flo.Problem(msh, dim, alpha, M_inf, c_ref, c_ref, 0.25, 0., 0.)
- if M_inf == 0:
- pbl.set(flo.Medium(msh, 'field', flo.F0ElRhoL(), flo.F0ElMachL(), flo.F0ElCpL(), flo.F0PsPhiInf(dim, alpha)))
- else:
- pbl.set(flo.Medium(msh, 'field', flo.F0ElRho(M_inf, M_crit), flo.F0ElMach(M_inf), flo.F0ElCp(M_inf), flo.F0PsPhiInf(dim, alpha)))
- pbl.add(flo.Initial(msh, 'field', flo.F0PsPhiInf(dim, alpha)))
- pbl.add(flo.Dirichlet(msh, 'upstream', flo.F0PsPhiInf(dim, alpha)))
- pbl.add(flo.Freestream(msh, 'side', tbox.Fct1C(-U_inf[0], -U_inf[1], 0.)))
- pbl.add(flo.Freestream(msh, 'downstream', tbox.Fct1C(-U_inf[0], -U_inf[1], 0.)))
- pbl.add(flo.Wake(msh, ['wake', 'wake_', 'field']))
- pbl.add(flo.Kutta(msh, ['te', 'wake_', 'airfoil', 'field']))
- bnd = flo.Boundary(msh, ['airfoil', 'field'])
- pbl.add(bnd)
- blw = flo.Blowing(msh, ['airfoil', 'field'])
- blw2 = flo.Blowing(msh, ['wake', 'field'])
- pbl.add(blw)
- pbl.add(blw2)
- tms['pre'].stop()
-
- # solve inviscid problem
- print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
- tms['solver'].start()
- isolver = flo.Newton(pbl)
- floD.initNewton(isolver)
- vsolver = floVS.Solver(Re)
- coupler = floVC.Coupler(isolver, vsolver, gmshWriter,blw,blw2)
- coupler.run()
- tms['solver'].stop()
-
- # extract Cp
- Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
- tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
- # display results
- print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
- print(' M alpha Cl Cd Cm')
- print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, isolver.Cm))
-
- # display timers
- tms['total'].stop()
- print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
- print('CPU statistics')
- print(tms)
-
- # visualize solution and plot results
- floD.initViewer(pbl)
- tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
-
- # check results
- print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
- tests = CTests()
- tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.4132, 1e-1)) # TODO check value and tolerance
- tests.add(CTest('Cl', isolver.Cl, 0.0, 5e-2))
- tests.add(CTest('Cdp', vsolver.Cdp, 0.00067, 0.01))
- tests.add(CTest('Cd', vsolver.Cd, 0.0051, 0.01))
- tests.add(CTest('top_xtr', vsolver.top_xtr, 0.4381, 0.5 ))
- tests.add(CTest('bot_xtr', vsolver.bot_xtr, 0.4368, 0.5 ))
- tests.run()
-
- # eof
- print('')
-
-if __name__ == "__main__":
- main()
diff --git a/flow/viscous/Airfoil.py b/flow/viscous/Airfoil.py
index 3c954dab33752fbb1e9ac7498933ec78622107b8..bb6d11f7f4f3f53ecd3bef743b99c98368d6e9a9 100644
--- a/flow/viscous/Airfoil.py
+++ b/flow/viscous/Airfoil.py
@@ -43,7 +43,7 @@ class Airfoil(Boundary):
self.H0 = 2.23 # One parameter family Falkner Skan
# self.H0 = 2.23
self.n0 = 0
- self.Ct0 = 0.03
+ self.Ct0 = 0
return self.T0, self.H0, self.n0, self.Ct0
def connectList(self):
@@ -52,7 +52,7 @@ class Airfoil(Boundary):
N1 = np.zeros(self.nN, dtype=int) # Node number
connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
- data = np.zeros((self.boundary.nodes.size(), 9))
+ data = np.zeros((self.boundary.nodes.size(), 10))
i = 0
for n in self.boundary.nodes:
data[i,0] = n.no
@@ -63,6 +63,8 @@ class Airfoil(Boundary):
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)
@@ -116,10 +118,11 @@ class Airfoil(Boundary):
data[:,5] = data[connectListNodes,5]
data[:,6] = data[connectListNodes,6]
data[:,7] = data[connectListNodes,7]
- data[:,8] = self.deltaStar
+ 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)
- return connectListElems,[uData, lData]
\ No newline at end of file
+ lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
+ return connectListNodes, connectListElems,[uData, lData]
\ No newline at end of file
diff --git a/flow/viscous/Boundary.py b/flow/viscous/Boundary.py
index cd9a1c20e6582a79bd616e492ada6290b60bcaa7..9526a829e42dee90e8126f2b552e291c66aaaae3 100644
--- a/flow/viscous/Boundary.py
+++ b/flow/viscous/Boundary.py
@@ -31,7 +31,7 @@ class Boundary(object):
self.u = np.zeros(self.nE) # blowing velocity
self.Me = np.zeros(self.nN) # Mach number at the edge of the boundary layer
self.rhoe = np.zeros(self.nN) # Density at the edge of the boundary layer
- self.connectListElems = np.zeros(self.nE) # Connectivity for the blowing velocity
- self.Cd = 0 # drag coefficient of the physical boundary
- self.Cdp = 0 # Form drag
- self.deltaStar = np.zeros(self.nN) # Displacement thickness
\ No newline at end of file
+ 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
\ No newline at end of file
diff --git a/flow/viscous/Wake.py b/flow/viscous/Wake.py
index 667230bd11ae4c3f1bc64c055655059493a1c05f..bae06e5331c69e7dce6cd92af0e7641433ac3fc5 100644
--- a/flow/viscous/Wake.py
+++ b/flow/viscous/Wake.py
@@ -29,7 +29,7 @@ class Wake(Boundary):
self.name = "wake"
self.T0 = 0 # initial condition for the momentum thickness
self.H0 = 0
- self.n0 = 0 # should not look at transition inside the wake
+ self.n0 = 9 # wake is always turbulent
self.Ct0 = 0
@@ -39,7 +39,7 @@ class Wake(Boundary):
N1 = np.zeros(self.nN, dtype=int) # Node number
connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
- data = np.zeros((self.boundary.nodes.size(), 9))
+ data = np.zeros((self.boundary.nodes.size(), 10))
i = 0
for n in self.boundary.nodes:
data[i,0] = n.no
@@ -50,6 +50,8 @@ class Wake(Boundary):
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)
@@ -69,7 +71,8 @@ class Wake(Boundary):
data[:,5] = data[connectListNodes,5]
data[:,6] = data[connectListNodes,6]
data[:,7] = data[connectListNodes,7]
- data[:,8] = self.deltaStar
+ data[:,8] = data[connectListNodes,8]
+ data[:,9] = data[connectListNodes,9]
# Separated upper and lower part
data = np.delete(data,0,1)
- return connectListElems, data
+ return connectListNodes, connectListElems, data
diff --git a/flow/viscous/coupler.py b/flow/viscous/coupler.py
index 723f509cb5359a5f5185411f7682a8b7848eb013..c2aabe26e7df30cd95143211a95c4aa461c61202 100644
--- a/flow/viscous/coupler.py
+++ b/flow/viscous/coupler.py
@@ -28,7 +28,7 @@ from flow.viscous.Airfoil import Airfoil
from flow.viscous.Wake import Wake
class Coupler(object):
- def __init__(self, _isolver, _vsolver, _writer, _boundaryAirfoil, _boundaryWake):
+ def __init__(self, _isolver, _vsolver, _writer, _boundaryAirfoil, _boundaryWake, _tol):
'''...
'''
self.isolver =_isolver # inviscid solver
@@ -37,23 +37,20 @@ class Coupler(object):
self.airfoil = Airfoil(_boundaryAirfoil) # create the object airfoil
self.wake = Wake(_boundaryWake) # create the object wake
self.group = [self.airfoil, self.wake]
+ self.tol = _tol # tolerance of the VII
def run(self):
''' Perform coupling
'''
# initialize loop
it = 0
- alpha = 1# underrelaxation factor
converged = 0 # temp
- tol = 10**-6
CdOld = self.vsolver.Cd
- # while it < 15:
while converged == 0:
print(ccolors.ANSI_BLUE + 'Iteration: ', it, ccolors.ANSI_RESET)
# run inviscid solver
self.isolver.run()
for n in range(0, len(self.group)):
- uOld = self.group[n].u
for i in range(0, len(self.group[n].boundary.nodes)):
self.group[n].v[i,0] = self.isolver.U[self.group[n].boundary.nodes[i].row].x[0]
self.group[n].v[i,1] = self.isolver.U[self.group[n].boundary.nodes[i].row].x[1]
@@ -62,26 +59,27 @@ class Coupler(object):
self.group[n].rhoe[i] = self.isolver.rho[self.group[n].boundary.nodes[i].row]
# run viscous solver
self.vsolver.run(self.group[n])
- self.group[n].u = self.group[n].u[self.group[n].connectListElems]
- self.group[n].u = uOld + alpha*(self.group[n].u - uOld) # underrelaxation
if self.group[n].name == "airfoil":
self.group[n+1].T0 = self.group[n].TEnd[0]+self.group[n].TEnd[1]
self.group[n+1].H0 = (self.group[n].HEnd[0]*self.group[n].TEnd[0]+self.group[n].HEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0
self.group[n+1].Ct0 = (self.group[n].CtEnd[0]*self.group[n].TEnd[0]+self.group[n].CtEnd[1]*self.group[n].TEnd[1])/self.group[n+1].T0
- self.group[n+1].n0 = self.group[n].nEnd[0]
# get blowing velocity from viscous solver and update inviscid solver
for i in range(0, self.group[n].nE):
self.group[n].boundary.setU(i, self.group[n].u[i])
print('Computation is done for: ' ,self.group[n].name)
print('Number of iterations:' ,+it)
print('Cd = ' ,+self.vsolver.Cd)
- print('Top_xtr = ' ,+self.vsolver.top_xtr)
- print('Bot_xtr = ' ,+self.vsolver.bot_xtr)
- if abs(self.vsolver.Cd - CdOld) < tol:
+ print('Top_xtr = ' ,+self.vsolver.xtr[0])
+ print('Bot_xtr = ' ,+self.vsolver.xtr[1])
+ print('Error = ' ,+abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd)
+ # Converged or not
+ if abs(self.vsolver.Cd - CdOld)/self.vsolver.Cd < self.tol:
converged = 1
else:
CdOld = self.vsolver.Cd
- it += 1
+ it += 1
+ self.vsolver.it += 1
# save results
self.isolver.save(0, self.writer)
+ self.vsolver.writeFile()
print('\n')
diff --git a/flow/viscous/newton.py b/flow/viscous/newton.py
index 262818eafa2114ad4cb2d8b6e599bbc0a7bb0650..e6130db646fcbc51eab04e9b11354ec64f82f4b8 100644
--- a/flow/viscous/newton.py
+++ b/flow/viscous/newton.py
@@ -27,43 +27,45 @@ import math
def newtonRaphson(f, x, maxIt, tol=1.0e-9):
- ''' Newton procedure to solve the coupled, non linear system. Boundary layer equation are parabolic in x --> Use of downstream marching type of solver
- As an implicit marching model is applied (Crank Nicolson), a matrix inversion is performed'''
- # Trial to develop Jacobian
- def jacobian(f, x):
- dx = 1.0e-8
- n = len(x)
- jac = np.zeros((n,n))
- f0 = f(x)
- for j in range(n):
- Dxj = (abs(x[j])*dx if x[j] != 0 else dx)
- x_plus = [(xi if k != j else xi+Dxj) for k, xi in enumerate(x)]
- jac[:,j] = (f(x_plus)-f0)/Dxj # numerical solution of the jacobian
- return jac, f0
+ ''' Newton procedure to solve the coupled, non linear system. Boundary layer equation are parabolic in x --> Use of downstream marching type of solver
+ As an implicit marching model is applied (Crank Nicolson), a matrix inversion is performed'''
+ # Trial to develop Jacobian
+ def jacobian(f, x):
+ dx = 1.0e-8
+ n = len(x)
+ jac = np.zeros((n,n))
+ f0 = f(x)
+ for j in range(n):
+ Dxj = (abs(x[j])*dx if x[j] != 0 else dx)
+ x_plus = [(xi if k != j else xi+Dxj) for k, xi in enumerate(x)]
+ jac[:,j] = (f(x_plus)-f0)/Dxj # numerical solution of the jacobian
+ return jac, f0
# Backtrack linesearch
- def BacktrackingLineSearch(f, x, dx, maxLsIt):
- a = 1
- f0 = f(x+dx)
- fevalIt = 1
- while fevalIt < maxLsIt:
- fa = f(x+a*dx)
- fevalIt += 1
- print('fa = ', + fa ,'f0 = ', +f0)
- if abs(fa[1]) > abs(f0[1]): # compare H
- a = a/2
- print('a = ' ,+a)
- else:
- break
- return a
+ def BacktrackingLineSearch(f, x, dx, maxLsIt):
+ a = 0.5
+ f0 = f(x+dx)
+ fevalIt = 1
+ while fevalIt < maxLsIt:
+ fa = f(x+a*dx)
+ fevalIt += 1
+ # print('fa = ', + fa ,'f0 = ', +f0)
+ if abs(any(fa)) > abs(any(f0)): # compare H
+ a = a/2
+ print('a = ' ,+a)
+ else:
+ break
+ return a
- for i in range(maxIt):
- jac, f0 = jacobian(f,x)
- if np.sqrt(np.dot(f0,f0)/len(x)) < tol and x[0]>0 and x[1]>0: return x
- dx = np.linalg.solve(jac, -f0)
- x = x + dx
+ for i in range(maxIt):
+ jac, f0 = jacobian(f,x)
+ if np.sqrt(np.dot(f0,f0)/len(x)) < tol and x[0]>0 and x[1]>0: return x
+ dx = np.linalg.solve(jac, -f0)
+ # a = BacktrackingLineSearch(f, x, dx, maxIt)
+ x = x + dx
+ x[1] = max(x[1],1.0005) # to avoid too low shape parameter. Solver must be redone to be more robust
# print('New x =',+ x)
- if np.sqrt(np.dot(dx, dx)) < tol*max(abs(any(x)),1) and x[0]>0 and x[1]>0: return x
- print("Too many iterations")
+ if np.sqrt(np.dot(dx, dx)) < tol*max(abs(any(x)),1) and x[0]>0 and x[1]>0: return x
+ print("Too many iterations")
diff --git a/flow/viscous/solver.py b/flow/viscous/solver.py
index e74482bf90fd94b39fd0b78ad0d58475ae0bcf37..a60fabacc99a39a1fe91e6cf574a45232bf3a297 100644
--- a/flow/viscous/solver.py
+++ b/flow/viscous/solver.py
@@ -24,7 +24,7 @@ from builtins import range
from builtins import object
from flow.viscous.Boundary import Boundary
from flow.viscous.Wake import Wake
-import flow.viscous.newton as newton
+import flow.viscous.Newton as Newton
import numpy as np
import math
import matplotlib.pyplot as plt
@@ -33,328 +33,440 @@ from scipy.io import loadmat
from scipy.signal import savgol_filter
class Solver(object):
- def __init__(self, _Re):
- '''...
+ '''To do:
+ - Improve the filter in boundaryLayer (l 263);
+ - Inverse procedure at the singularity location ( such as in Xfoil);
+ - Create a new mesh for the viscous solver;
+ - Improve the transition solver by refining the grid to split the transitional element into a laminar element and a turbulent element.
+ Here it is an average between laminar and turbulent solution (l 397-416);
+ - 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);
+ - Verify why the computation crash if the first point of the wake is computed with the turbulent closure (l 309).
+ '''
+
+ def __init__(self, _Re, _p, _m):
+ '''Initialize the viscous solver
'''
self.Re = _Re # Reynolds number
- self.gamma = 1.4
- self.Cd = 0
+ self.p = _p # Filtering order
+ self.m = _m # Filtering bandwith
+ 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 run(self, group):
- ''' Solve BLE
- '''
- def writeFile(data, H, theta, Cf, vt, Cdissip, Ct, Dstar, blwUe):
- # save_path = 'D:/Amaury/Documents/TFE/xfoil/'
- name = 'BLparameters'
- name2 = 'blwU'
- # complete_name = os.path.join(save_path, name+".dat")
- # complete_name2 = os.path.join(save_path, name2+".dat")
- f = open(name, 'w+')
- for i in range(len(H)):
- f.write("%s %s %s %s %s %s %s %s\n" % (data[i,0], H[i], theta[i], Cf[i], vt[i], Cdissip[i], Ct[i], Dstar[i]))
- f.close()
- f = open(name2, 'w+')
- for i in range(len(blwUe)):
- f.write("%s %s\n" % (0, blwUe[i]))
- f.close()
+ def dictionary(self):
+ '''Create dictionary with the needed 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])
+ wBlVariables = {}
+ wBlVariables['Cd'] = self.blScalars[0]
+ wBlVariables['Cdf'] = self.blScalars[1]
+ wBlVariables['Cdp'] = self.blScalars[2]
+ wBlVariables['delta*'] = self.blVectors[0]
+ wBlVariables['theta'] = self.blVectors[1]
+ wBlVariables['H'] = self.blVectors[2]
+ wBlVariables['Hk'] = self.blVectors[3]
+ wBlVariables['H*'] = self.blVectors[4]
+ wBlVariables['H**'] = self.blVectors[5]
+ wBlVariables['Cf'] = self.blVectors[6]
+ wBlVariables['Cdissip'] = self.blVectors[7]
+ wBlVariables['Ctau'] = self.blVectors[8]
+ wBlVariables['xtrTop'] = self.xtr[0]
+ wBlVariables['xtrBot'] = self.xtr[1]
+ return wData, wBlVariables
- def __getBLcoordinates(data):
- nN = len(data[:,0])
- nE = nN-1 #nbr of element along the surface
- vt = np.zeros(nN)
- xx = np.zeros(nN) # position along the surface of the airfoil
- dx = np.zeros(nE) # distance along two nodes
- dv = np.zeros(nE) # speed difference according to element
- alpha = np.zeros(nE) # angle of the element for the rotation matrix
- for i in range(0,nE):
- alpha[i] = np.arctan2((data[i+1,1]-data[i,1]),(data[i+1,0]-data[i,0]))
- dx[i] = np.sqrt((data[i+1,0]-data[i,0])**2+(data[i+1,1]-data[i,1])**2)
- xx[i+1] = dx[i]+xx[i]
- vt[i] = abs(data[i,3]*np.cos(alpha[i]) + data[i,4]*np.sin(alpha[i])) # tangential speed at node 1 element i
- vt[i+1] = abs(data[i+1,3]*np.cos(alpha[i]) + data[i+1,4]*np.sin(alpha[i])) # tangential speed at node 2 element i
- dv[i] = (vt[i+1] - vt[i])/dx[i] #velocity gradient according to element i. central scheme with half point
- return xx, dx, dv, vt, nN
+ def writeFile(self):
+ '''Write the results in airfoil_viscous.dat'''
+ wData, wBlVariables = self.dictionary()
+ f = open('airfoil_viscous.dat', 'w+')
+ f.write('$Aerodynamic coefficients \n')
+ f.write("%s, %s, %s, %s, %s, %s \n" % ('Re','Cd', 'Cdf', 'Cdp', 'xtrTop', 'xtrBot'))
+ f.write("%f, %f, %f, %f, %f, %f \n" %(self.Re, wBlVariables['Cd'], wBlVariables['Cdf'], wBlVariables['Cdp'], wBlVariables['xtrTop'], wBlVariables['xtrBot'] ))
+ f.write('$Boundary layer variables \n')
+ f.write("%s, %s, %s, %s, %s, %s, %s \n" % ('x','y', 'z', 'x/c', 'delta*', 'H', 'Cf'))
+ for i in range(len(self.data[:,0])):
+ f.write("%e, %e, %e, %e, %e, %e, %e \n" % (wData['x'][i], wData['y'][i], wData['z'][i], wData['x/c'][i], wBlVariables['delta*'][i], wBlVariables['H'][i], wBlVariables['Cf'][i]))
+ f.close()
+ print('writing file: airfoil_viscous.dat... done!')
- def boundaryLayer(data):
- # Get the BL coordinates
- xx, dx, dv, vt, nN = __getBLcoordinates(data)
- if group.name == 'airfoil':
- vt = savgol_filter(vt, 11, 2, mode = 'interp')
- Me = data[:,5]
- rhoe = data[:,6]
- deltaStarOld = data[:,7]
- #Initialize variable
- blwVel = np.zeros(nN-1) # Blowing velocity
- theta = np.zeros(nN) # Momentum thickness
- H = np.zeros(nN) # Shape factor
- deltaStar = np.zeros(nN) # Displacement thickness
- Hstar = np.zeros(nN) # Kinetic energy shape parameter
- Hstar2 = np.zeros(nN) # Density shape parameter
- Hk = np.zeros(nN) # Kinematic shape factor
- Ret = np.zeros(nN) #Re based on momentum thickness
- Cd = np.zeros(nN) # Dissipation coefficient
- Cf = np.zeros(nN) # Skin friction
- n = np.zeros(nN) # Logarithm of the maximum amplification ratio
- dndRet = np.zeros(nN) # Variation of n according to Re_theta
- m = np.zeros(nN) # Empirical relation
- l = np.zeros(nN) # Idem
- Ct = np.zeros(nN) # Shear stress coefficient
- Cteq = np.zeros(nN) # Equilibrium shear stress coefficient
- delta = np.zeros(nN) # BL thickness
- # Initial condition
- if group.name == "airfoil":
- theta[0], H[0], n[0], Ct[0] = group.initialConditions(self.Re, dv[0])
- elif group.name == "wake":
- theta[0] = group.T0
- H[0] = group.H0
- n[0] = group.n0
- Ct[0] = group.Ct0
- idxTr = 0 # To Do: it is the index of the transition, should find a way to delete it for the wake
- # vt[0] = vt[0]+ 4/(np.pi*dx[i-1])*(theta[0]*H[0]-deltaStarOld[i]) !!!! remains vt[0]?
- Ret[0] = vt[0]*theta[0]*self.Re
- deltaStar[0] = theta[0]*H[0]
- # print('---------------Grid point n 0---------------------')
- Cd[0], Cf[0], Hstar[0], Hstar2[0], dndRet[0], m[0], l[0], Hk[0], delta[0] = newLaminarClosure( theta[0], H[0], Ret[0],Me[0], group.name)
- ntr = 9 # amplification ratio at which transition occurs TO DO: params from user (if not, impose a value of 9)
- if n[0] < 9:
- turb =0 # we begin with laminar flow
- else :
- turb =1 # typically for the wake
- xtr = 0 # if only turbulent flow
- # print('------------The solution at grid point ',+0,'=', + theta[0], +H[0], +Ct[0],+vt[0], 'X position:', + data[0,0])
- for i in range(1,nN):
- '''x[0] = theta; x[1] = H; x[2] = n for laminar/Ct for turbulent; x[3] = vt from interaction law
- '''
- # print('---------------Grid point n',+i,'---------------------')
- def f_lam(x):
- f = np.zeros(len(x))
- Ret[i] = x[3]*x[0]*self.Re
- Cd[i], Cf[i], Hstar[i], Hstar2[i], dndRet[i], m[i], l[i], Hk[i], delta[i] = newLaminarClosure(x[0], x[1], Ret[i],Me[i], group.name)
- dTheta = x[0]-theta[i-1]
- due = x[3]-vt[i-1]
- dH = x[1]-H[i-1]
- dHstar = Hstar[i]-Hstar[i-1]
- Ta = upw*x[0]+(1-upw)*theta[i-1]
- RhsTheta1 = (2 + x[1] - Me[i]**2) * (x[0]/x[3]) * (due/dx[i-1]) - Cf[i]/2
- RhsTheta2 = (2 + H[i-1] - Me[i-1]**2) * (theta[i-1]/vt[i-1]) * (due/dx[i-1])- Cf[i-1]/2
- RhsH1 = (2*Hstar2[i] + Hstar[i]*(1-x[1])) * (x[0]/x[3]) * (due/dx[i-1]) - 2*Cd[i] + 0.5*Cf[i]*Hstar[i]
- RhsH2 = (2*Hstar2[i-1] + Hstar[i-1]*(1-H[i-1])) * (theta[i-1]/vt[i-1]) * (due/dx[i-1]) - 2*Cd[i-1] + 0.5*Cf[i-1]*Hstar[i-1]
- RhsN1 = dndRet[i]*((m[i])+1)/2 * l[i]/x[0]
- RhsN2 = dndRet[i-1]*((m[i-1])+1)/2 * l[i-1]/theta[i-1]
- f[0] = dTheta/dx[i-1] + ((1-upw)*RhsTheta2+upw*RhsTheta1)
- f[1] = Ta*(dHstar/dx[i-1]) + ((1-upw)*RhsH2+upw*RhsH1)
- f[2] = (x[2]-n[i-1])/dx[i-1] - ((1-upw)*RhsN2+upw*RhsN1)
- # f[3] = x[3]-vt[i]- (4/(np.pi*dx[i-1]))*(x[0]*x[1]-deltaStarOld[i])
- f[3] = x[3] -vt[i]
- return f
- # Central differencing to gain accuracy
- def f_turb(x):
- f = np.zeros(len(x))
- Ret[i] = x[3]*x[0]*self.Re
- Cd[i], Cf[i], Hstar[i], Hstar2[i],Hk[i],Cteq[i], delta[i] = newTurbulentClosure(x[0], x[1],x[2], Ret[i],Me[i], group.name)
- dTheta = x[0]-theta[i-1]
- due = x[3]-vt[i-1]
- dH = x[1]-H[i-1]
- dHstar = Hstar[i]-Hstar[i-1]
- dCt = x[2]-Ct[i-1] # here it is Ct^1/2 which is represented
- Ta = upw*x[0]+(1-upw)*theta[i-1]
- Cta = upw*x[2]+(1-upw)*Ct[i-1]
- deriv = 1/Cta * dCt/dx[i-1]
- RhsTheta1 = (2 + x[1] - Me[i]**2) * (x[0]/x[3]) * (due/dx[i-1]) - Cf[i]/2
- RhsTheta2 = (2 + H[i-1] - Me[i-1]**2) * (theta[i-1]/vt[i-1]) * (due/dx[i-1])- Cf[i-1]/2
- RhsH1 = (2*Hstar2[i] + Hstar[i]*(1-x[1])) * (x[0]/x[3]) * (due/dx[i-1]) - 2*Cd[i] + 0.5*Cf[i]*Hstar[i]
- RhsH2 = (2*Hstar2[i-1] + Hstar[i-1]*(1-H[i-1])) * (theta[i-1]/vt[i-1]) * (due/dx[i-1]) - 2*Cd[i-1] + 0.5*Cf[i-1]*Hstar[i-1]
- delta1 =delta[i]
- delta2 = delta[i-1]
- TmpHk1 = 4/(3*x[1]*x[0]) * (0.5*Cf[i] - ((Hk[i]-1)/(6.7*Hk[i]))**2)
- TmpHk2 = 4/(3*H[i-1]*theta[i-1]) * (0.5*Cf[i-1] - ((Hk[i-1]-1)/(6.7*Hk[i-1]))**2)
- TmpUe1 = 1/x[3] * due/dx[i-1]
- TmpUe2 = 1/vt[i-1] * due/dx[i-1]
- TmpCteq1 = 5.6*(Cteq[i]-x[2])
- TmpCteq2 = 5.6*(Cteq[i-1]-Ct[i-1])
- f[0] = dTheta/dx[i-1] + ((1-upw)*RhsTheta2+upw*RhsTheta1)
- f[1] = Ta*(dHstar/dx[i-1]) + ((1-upw)*RhsH2+upw*RhsH1)
- f[2] = 2*((upw*delta1+(1-upw)*delta2)*(deriv - (upw*TmpHk1+(1-upw)*TmpHk2) + (upw*TmpUe1+(1-upw)*TmpUe2))) - (upw*TmpCteq1+(1-upw)*TmpCteq2)
- # f[3] = x[3]-vt[i]- 4/(np.pi*dx[i-1])*(x[0]*x[1]-deltaStarOld[i])
- f[3] = x[3] -vt[i]
- return f
- '''Solver depending on the case '''
- # Laminar solver
- if turb == 0:
- if n[i-1] > 8.5:
- upw = 0.5 # Trapezoidal scheme
- elif i < 7:
- upw = 1 # backward Euler
- else:
- upw = 0.5 # Trapezoidal scheme
- x = np.array([theta[i-1],H[i-1], n[i-1], vt[i]])
- sol = newton.newtonRaphson(f_lam, x, 120)
- logRet_crit = 2.492*(1/(H[i-1]-1))**0.43 + 0.7*(np.tanh(14*(1/(H[i-1]-1))-9.24)+1.0) # value at which the amplification starts to growth
- if np.log10(Ret[i]) <= logRet_crit - 0.08 :
- n[i] = 0
- else:
- n[i] = sol[2]
- xtr = 1 # no transition if remains laminar
- # Turbulent solver
- elif turb == 1:
- upw = 0.5 # trapezoidal scheme
- x = np.array([theta[i-1], H[i-1], Ct[i-1], vt[i]])
- sol = newton.newtonRaphson(f_turb, x, 120)
- Ct[i] = sol[2]
- # print('------------The solution at grid point ',+i,'=', + sol,+vt[i], 'X position:', + data[i,0])
- theta[i] = sol[0]
- H[i] = sol[1]
- vt[i] = sol[3]
- # Transition solver
- if n[i] >= ntr and turb ==0:
- # Initial condition
- turb=1
- Ct[i-1] = Ct[0]
- # Interpolation to define position of the transition
- idxTr = i
- xtr = (ntr-n[i-1])*((data[i,0]-data[i-1,0])/(n[i]-n[i-1]))+data[i-1,0]
- avTurb = (data[i,0]-xtr)/(data[i,0]-data[i-1,0]) # percentage of the subsection corresponding to a laminar flow
- avLam = 1- avTurb # percentage of the subsection corresponding to a turbulent flow
- # Compute boundary layer parameters with FOU for turbulent case
- upw = 0
- x_turbTr = np.array([theta[i-1], 1.515, Ct[i-1], vt[i]])
- sol_turbTr = newton.newtonRaphson(f_turb, x_turbTr, 30)
- # # Averaging both solutions
- theta[i] = avLam*sol[0]+avTurb*sol_turbTr[0]
- H[i] = avLam*sol[1]+avTurb*sol_turbTr[1]
- vt[i] = avLam*sol[3]+avTurb*sol_turbTr[3]
- # theta[i] = avTurb*sol_turbTr[0]
- # H[i] = avTurb*sol_turbTr[1]
- Ct[i] = avTurb*sol_turbTr[2]
- # print('Ok for the transition part')
- # print('------------The solution at grid point ',+i,'=', +theta[i], +H[i], +Ct[i],+vt[i], 'X position:', + data[i,0])
- deltaStar[i] = H[i]*theta[i]
- blwVel[i-1] = (rhoe[i]*vt[i]*deltaStar[i]-rhoe[i-1]*vt[i-1]*deltaStar[i-1])/(rhoe[i]*dx[i-1]) # TO DO: need to divide by rhoe or not?
- Cdf = np.trapz(Cf, data[:,0])
- Cdtot = 2*theta[-1]*(vt[-1]**((H[-1]+5)/2))
- blVariables = [blwVel, deltaStar, theta, H, Cf, xtr, Cdf, Cdtot, vt, Cd ,n[idxTr], Ct]
- return blVariables
-
- def newLaminarClosure(theta, H, Ret,Me, name):
- ''' Laminar closure derived from the Falkner-Skan one-parameter profile family
- Drela(1996)'''
- Hk = (H - 0.29*Me**2)/(1+0.113*Me**2) # Kinematic shape factor
- Hstar2 = (0.064/(Hk-0.8)+0.251)*Me**2 # Density shape parameter
- dndRet = 0.01*np.sqrt((2.4*Hk-3.7+2.5*np.tanh(1.5*Hk-4.65))**2+0.25)
- l = (6.54*Hk-14.07)/(Hk**2)
- m = (0.058*(((Hk-4)**2)/(Hk-1))-0.068)/l
- Hmi = (1/(H-1))
- # m = -0.05+2.7*Hmi-5.5*Hmi**2+3*Hmi**3+0.1*np.exp(-20*Hmi)
- # dndRet = 0.028*(Hk-1)-0.0345*np.exp((-3.87*Hmi+2.52)**2)
- delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta)
- Hks = Hk - 4.35
- if Hk <= 4.35:
- dens = Hk + 1
- Hstar = 0.0111*Hks**2/dens - 0.0278*Hks**3/dens + 1.528 -0.0002*(Hks*Hk)**2
- elif Hk > 4.35:
- Hstar = 1.528 + 0.015*Hks**2/Hk
- Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction
- if Hk <= 5.5:
- tmp = (5.5-Hk)**3 / (Hk+1)
- Cf = (-0.07 + 0.0727*tmp)/ Ret
- elif Hk > 5.5:
- print('Laminar separation')
- tmp = 1 - 1/(Hk-4.5)
- Cf = (-0.07 + 0.015*tmp**2) / Ret
- if Hk <= 4:
- Cd = ( 0.00205*(4.0-Hk)**5.5 + 0.207 ) * (Hstar/(2*Ret))
- elif Hk > 4:
- HkCd = (Hk-4)**2
- denCd = 1+0.02*HkCd
- Cd = ( -0.0016*HkCd/denCd + 0.207) * (Hstar/(2*Ret))
- if name == 'wake':
- Cd = 2*(1.10*(1.0-1.0/Hk)**2/Hk)*(Hstar/(2*Ret))
- # print(' Ret, Cd, Cf, Hk, H* = ',+ Ret, + Cd, +Cf, +Hk, + Hstar)
- return Cd, Cf, Hstar, Hstar2, dndRet, m, l, Hk, delta
+ def __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 newTurbulentClosure(theta, H, Ct, Ret, Me, name):
- ''' Turbulent closure derived using the skin friction and velocity profile formulas of Swafford (1983)'''
- Hk = (H - 0.29*Me**2)/(1+0.113*Me**2)
- if name == "airfoil":
- Hk = max(Hk, 1.05)
+ 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:
- Hk = max(Hk, 1.00005)
- Hstar2 = ((0.064/(Hk-0.8))+0.251)*Me**2
- gam = self.gamma - 1
- Fc = np.sqrt(1+0.5*gam*Me**2)
- if Ret <= 400:
- H0 = 4
- elif Ret > 400:
- H0 = 3 + (400/Ret)
- if Ret > 200:
- Ret = Ret
- elif Ret <= 200:
- Ret = 200
- if Hk <= H0:
- Hstar = ((0.5-4/Ret)*((H0-Hk)/(H0-1))**2)*(1.5/(Hk+0.5))+1.5+4/Ret
- elif Hk > H0:
- print('The flow is separated here')
- Hstar = (Hk-H0)**2*(0.007*np.log10(Ret)/(Hk-H0+4/np.log10(Ret))**2 + 0.015/Hk) + 1.5 + 4/Ret
- Hstar = (Hstar + 0.028*Me**2)/(1+0.014*Me**2) # Whitfield's minor additional compressibility correction
- logRt = np.log10(Ret/Fc)
- logRt = np.max((logRt,3))
- arg = np.max((-1.33*Hk, -20))
- Us = (Hstar/2)*(1-4*(Hk-1)/(3*H)) # Equivalent normalized wall slip velocity
- delta = min((3.15+ H + (1.72/(Hk-1)))*theta, 12*theta)
- Ctcon = 0.5/(6.7**2*0.75)
- if name == "wake":
- if Us > 0.99995:
- Us = 0.99995
- n = 2 # two wake halves
- Cf = 0 # no friction inside the wake
- Hkc = Hk - 1
- Cdi = 0.03*0.75**3*Us**3 # inner shear layer Row 1062 from xfoil
- else:
- if Us > 0.95:
- Us = 0.98
- n = 1
- # Cf = (1/Fc)*(0.3*np.exp(arg)*(logRt)**(-1.74-0.31*Hk)+0.00011*(np.tanh(4-(Hk/0.875))-1))
- Cf0 = 0.01013/(logRt-1.02) - 0.00075
- Ho = 1/(1-6.55*np.sqrt(Cf0/2))
- Cf = (1/Fc)*Cf0*(-0.5 + 0.9/(Hk/Ho-0.4))
- Hkc = Hk - 1 - 18/Ret
- Hkc = max(Hkc, 0.01)
- Cdi = 0
- Cd = n*((Cf*Us/2)+(0.995-Us)*Ct**2 +0.15*(0.995-Us)**2/Ret)
- Cteq = np.sqrt(Hstar*Ctcon*((Hk-1)*Hkc**2)/((1-Us)*((Hk**2)*H))) #Here it is Cteq^0.5
- # print('Cteq, Cd, Cf, Hk, H* = ', +Cteq, + Cd, +Cf, +Hk, + Hstar)
- return Cd, Cf, Hstar, Hstar2, Hk, Cteq, delta
+ 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, H
+
+ 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, Ueq, Ct
- group.connectListElems, data = group.connectList()
+
+ def __boundaryLayer(self, data, group):
+ '''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]
+ #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')
+ vtInv = vtOld
+ rhoInv = rhoe
+ 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
+ Ueq = np.zeros(nN) # Equilibrium velocity gradient
+ 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
+ if n[0] < 9:
+ 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))
+ 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)
+ 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)
+ 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], 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)
+ 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]])
+ sol = Newton.newtonRaphson(f_lam, x, 200)
+ # 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]])
+ sol = Newton.newtonRaphson(f_turb, x, 200)
+ 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]])
+ sol_turb = Newton.newtonRaphson(f_turb, x_turb, 200)
+ 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], Ueq[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?
+ # 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]
+ 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:
- ublVariables = boundaryLayer(data[0])
- lblVariables = boundaryLayer(data[1])
- blwVel = abs(np.concatenate((ublVariables[0], lblVariables[0])))
- group.Cd = ublVariables[7] + lblVariables[7]
- group.Cdp = group.Cd-(ublVariables[6]+lblVariables[6])
- self.top_xtr = ublVariables[5]
- self.bot_xtr = lblVariables[5]
- self.Cdp = group.Cdp
- group.TEnd = [ublVariables[2][-1], lblVariables[2][-1]]
- group.HEnd = [ublVariables[3][-1], lblVariables[3][-1]]
- group.CtEnd = [ublVariables[11][-1], lblVariables[11][-1]]
- group.nEnd = [ublVariables[10], lblVariables[10]]
- # Write the results
- writeFile(np.concatenate((data[0], data[1])), np.concatenate((ublVariables[3], lblVariables[3])), np.concatenate((ublVariables[2], lblVariables[2])),
- np.concatenate((ublVariables[4], lblVariables[4])), np.concatenate((ublVariables[8], lblVariables[8])),
- np.concatenate((ublVariables[9], lblVariables[9])), np.concatenate((ublVariables[11], lblVariables[11])),
- np.concatenate((ublVariables[1], lblVariables[1])), blwVel)
- lblVariables[1] = np.delete(lblVariables[1],0,0) # issue with the number of point
- group.deltaStar = np.concatenate((ublVariables[1], lblVariables[1]))
+ ublwVel, xtrTop, uxx, ublScalars, ublVectors = self.__boundaryLayer(data[0], group)
+ 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))
+ 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]]
+ 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:
- blVariables = boundaryLayer(data)
- blwVel = abs(blVariables[0])
- group.Cd = blVariables[7]
- group.Cdp = group.Cd-blVariables[6]
- self.Cd = group.Cd
- group.deltaStar = blVariables[1]
-
- group.u = blwVel
+ 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
+ # 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]