2.1.0 - Feature viscous-inviscid interaction

Waves v2.1.0 (flow-v1.7) - Viscous-inviscid interaction

This MR adds a viscous-inviscid interaction capability to Flow.

Changes

• Add VII feature, by @Amaury.Bilocq
• Minor refactoring and bugfixing

Viscous-inviscid interaction

Method
The viscous-inviscid interaction (VII) method allows to correct an inviscid flow computation by modeling the boundary layer, thus taking into account the viscous effects only where they are relevant. The VII method implemented in Flow relies on a quasi-simultaneous coupling of the full potential equation to an integral form of the boundary layer equations. More details can be found in Amaury Bilocq's master thesis.
Implementation
The code related to the VII method is written in python and located under viscous. A unit test, consisting of three test cases, is also available. Note that these tests slightly differ from those described in the master thesis, as a bug was found and corrected after the publication of the thesis.
Further work
Although the present code works and produce satisfactory results, it is limited to 2D subsonic flows and not very robust, as a filter is required to smoothen the oscillations present in the inviscid variables passed to the viscous solver. Furthermore, the code could be refactored for more modularity, or even translated to C++. Suggestions for future work are listed in the viscous solver.

Tests

Tests are passing on ubuntu20.04 (python 3.8.2), and msys2 (python 3.8.1).

flow/tests/bli.py

@@ -16,9 +16,18 @@
 # limitations under the License.
 
 
-# @brief Compute lifting (linear or nonlinear) viscous flow around a naca0012
-# @authors Adrien Crovato, Amaury Bilocq
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+#
+# 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
+#    -> nIt = 41, 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
+#    -> 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
 #
 # CAUTION
 # This test is provided to ensure that the solver works properly.
@@ -29,10 +38,11 @@ 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 fwk
 from fwk.testing import *
-from fwk.coloring import ccolors 
+from fwk.coloring import ccolors
 
 def main():
     # timer
@@ -40,16 +50,23 @@ def main():
     tms['total'].start()
 
     # define flow variables
-    alpha = 0*math.pi/180
-    M_inf = 0.0
+    Re = 1e7
+    alpha = 5*math.pi/180
+    U_inf = [math.cos(alpha), math.sin(alpha)] # norm should be 1
+    M_inf = 0.5
     M_crit = 5 # Squared critical Mach number (above which density is modified)
     c_ref = 1
     dim = 2
 
+    # define filter parameters and tolerance of the VII 
+    p = 2
+    m = 7
+    tol = 1e-4
+
     # mesh the geometry
     print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
     tms['msh'].start()
-    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.0075, 'msLe' : 0.0075}
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
@@ -62,8 +79,8 @@ def main():
     print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
     tms['solver'].start()
     isolver = floD.newton(pbl)
-    vsolver = floVS.Solver(blw)
-    coupler = floVC.Coupler(isolver, vsolver, gmshWriter)
+    vsolver = floVS.Solver(Re, p, m)
+    coupler = floVC.Coupler(isolver, vsolver, blw[0], blw[1], tol, gmshWriter)
     coupler.run()
     tms['solver'].stop()
 
@@ -72,8 +89,8 @@ def main():
     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, isolver.Cd, isolver.Cm))
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cd, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
 
     # display timers
     tms['total'].stop()
@@ -83,16 +100,32 @@ def main():
 
     # 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, isolver.Cd, isolver.Cm, 4), True)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cd, isolver.Cm, 4), True)
 
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
-    if M_inf == 0 and alpha == 0*math.pi/180:
-        tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.41, 1e-1)) # TODO check value and tolerance
-        tests.add(CTest('Cl', isolver.Cl, 0., 5e-2))
+    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('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_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('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_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('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('xtr_bot', vsolver.xtr[1], 1.000, 0.05))
     else:
         raise Exception('Test not defined for this flow')
+    
     tests.run()
 
     # eof