Add wht case. Updt wbht geo. Minor changes.

flow/cases/wbht.py

@@ -5,7 +5,7 @@
 # @authors Adrien Crovato
 #
 # CAUTION
-# Results might be nonphysical!
+# Gmsh might not generate the wing/tail surfaces properly. Results might not be accurate 
 
 def getParam():
     p = {}
@@ -34,7 +34,7 @@ def getParam():
                  'msFus' : fums, 'msFusCaps' : fucms, 'msFar' : fams} # Parameters for input file model
     p['Dim'] = 3 # Problem dimension
     p['Format'] = 'vtk' # Save format (vtk or gmsh)
-    p['Slice'] = [0.83, 3.0, 4., 6.1, 6.9] # array of (y) coordinates to perform slice along the span (empty if none)
+    p['Slice'] = [0.83, 3.0, 4.0, 6.2, 6.9] # array of (y) coordinates to perform slice along the span (empty if none)
     p['TagId'] = 9 # tag number of physical group to be sliced
     # Markers
     p['Fluid'] = 'field' # Name of physical group containing the fluid
@@ -46,7 +46,7 @@ def getParam():
     p['WakeTips'] = ['wakeTipW', 'wakeTipT'] # LIST of names of physical group containing the edge of the wake
     p['TeTips'] = ['teTipW', 'teTipT'] # LIST of names of physical group containing the edge of the wake and the trailing edge
     # Freestream
-    p['M_inf'] = 0.6 # Freestream Mach number
+    p['M_inf'] = 0.75 # Freestream Mach number
     p['CL'] = 0.4 # Target lift coefficient
     p['AoA'] = 0. # Freestream angle of attack (guess)
     p['dCL'] = 7.5 # Slope of CL w.r.t AoA [1/rad] (guess)

flow/cases/wht.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# @brief Wing-tail configuration file for flow trim script
+# @authors Adrien Crovato
+
+def getParam():
+    p = {}
+    # Specific
+    scl = 2 # scaling factor for lifting surface mesh size
+    wrlems = scl*0.0210 # wing root leading edge mesh size
+    wrtems = scl*0.0210 # wing root trailing edge mesh size
+    wtlems = scl*0.0090 # wing tip leading edge mesh size
+    wttems = scl*0.0090 # wing tip trailing edge mesh size
+    wwlems = scl*0.0042 # wing winglet leading mesh size
+    wwtems = scl*0.0060 # wing winglet trailing mesh size
+    trlems = scl*0.0089 # tail root leading mesh size
+    trtems = scl*0.0089 # tail root trailing mesh size
+    ttlems = scl*0.0058 # tail tip leading mesh size
+    tttems = scl*0.0058 # tail tip trailing mesh size
+    fams = 5. # farfield mesh size
+    # Input/Output
+    p['File'] = '../models/wht.geo' # Input file containing the model
+    p['Pars'] = {'msLeW0' : wrlems, 'msTeW0' : wrtems,
+                 'msLeW1' : wtlems, 'msTeW1' : wttems,
+                 'msLeW2' : wwlems, 'msTeW2' : wwtems,
+                 'msLeT0' : trlems, 'msTeT0' : trtems,
+                 'msLeT1' : ttlems, 'msTeT1' : tttems,
+                 'msFar' : fams} # Parameters for input file model
+    p['Dim'] = 3 # Problem dimension
+    p['Format'] = 'vtk' # Save format (vtk or gmsh)
+    p['Slice'] = [0.01, 3.0, 4.0, 6.2, 6.9] # array of (y) coordinates to perform slice along the span (empty if none)
+    p['TagId'] = 8 # tag number of physical group to be sliced
+    # Markers
+    p['Fluid'] = 'field' # Name of physical group containing the fluid
+    p['Farfield'] = ['farfield', 'downstream'] # LIST of names of physical groups containing the farfield boundaries (downstream should be last element)
+    p['Symmetry'] = 'symmetry' # Name of physical group containing the symmetry boundaries
+    p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface boundary
+    p['Wakes'] = ['wakeW', 'wakeT'] # LIST of names of physical group containing the wake
+    p['WakeTips'] = ['wakeTipW', 'wakeTipT'] # LIST of names of physical group containing the edge of the wake
+    p['TeTips'] = ['teTipW', 'teTipT'] # LIST of names of physical group containing the edge of the wake and the trailing edge
+    # Freestream
+    p['M_inf'] = 0.75 # Freestream Mach number
+    p['CL'] = 0.4 # Target lift coefficient
+    p['AoA'] = 0. # Freestream angle of attack (guess)
+    p['dCL'] = 7.5 # Slope of CL w.r.t AoA [1/rad] (guess)
+    # Geometry
+    p['S_ref'] = 8.9 # Reference surface length
+    p['c_ref'] = 1.7 # Reference chord length
+    p['x_ref'] = 0.0 # Reference point for moment computation (x)
+    p['y_ref'] = 0.0 # Reference point for moment computation (y)
+    p['z_ref'] = 0.0 # Reference point for moment computation (z)
+    # Numerical
+    p['NSolver'] = 'Newton' # Noninear solver type (Picard or Newton)
+    p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
+    p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
+    p['Max_it'] = 15 # Solver maximum number of iterations
+    p['LS_tol'] = 1e-6 # Tolerance on linesearch residual
+    p['Max_it_LS'] = 10 # Linesearch maximum number of iterations
+    p['AV_thrsh'] = 1e-2 # Residual threshold below which the artificial viscosity is decreased
+    p['M_crit'] = 5. # Critical Mach number above which density is damped
+    return p
+
+def main():
+    # Arguments parser
+    import argparse
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-k", help="nb of threads", type=int, default=1)
+    args, extra = parser.parse_known_args()
+
+    # Script call
+    import flow.scripts.trim as t
+    trim = t.Trim(getParam(), args.k)
+    trim.run()
+    trim.disp()
+
+if __name__ == "__main__":
+    main()
\ No newline at end of file

flow/scripts/polar.py

@@ -71,7 +71,7 @@ class Polar(Config):
         print "       M    alpha       Cl       Cd       Cm"
         i = 0
         while i < len(self.alphas):
-            print "{0:8.1f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}".format(self.mach, self.alphas[i], self.Cls[i], self.Cds[i], self.Cms[i])
+            print "{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}".format(self.mach, self.alphas[i], self.Cls[i], self.Cds[i], self.Cms[i])
             i = i+1
         print '\n'
         # display timers

flow/scripts/trim.py

@@ -52,7 +52,7 @@ class Trim(Config):
             if it != 0:
                 self.dCl = (self.solver.Cl - Cl_) / (self.alpha - alpha_)
             # break or store values and update AoA
-            if abs(self.ClTgt - self.solver.Cl) < 0.005:
+            if abs(self.ClTgt - self.solver.Cl) < 0.005 or math.isnan(self.solver.Cl):
                 break
             else:
                 Cl_ = self.solver.Cl
@@ -75,7 +75,7 @@ class Trim(Config):
         # display results
         print ccolors.ANSI_BLUE + 'Trim analysis' + ccolors.ANSI_RESET
         print "       M    alpha      dCl       Cl       Cd       Cm"
-        print "{0:8.1f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}".format(self.mach, self.alpha*180/math.pi, self.dCl, self.solver.Cl, self.solver.Cd, self.solver.Cm)
+        print "{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}".format(self.mach, self.alpha*180/math.pi, self.dCl, self.solver.Cl, self.solver.Cd, self.solver.Cm)
         print '\n'
         # display timers
         self.tms["total"].stop()