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

@@ -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