diff --git a/dart/api/core.py b/dart/api/core.py index 7da8be7f62f8f1ae029cd8e2210b67f6c210c023..e729b8ddbbc215b546e8e42a4b8e7a26737aef3e 100644 --- a/dart/api/core.py +++ b/dart/api/core.py @@ -136,7 +136,7 @@ def initDart(params, scenario='aerodynamic', task='analysis'): if 'Symmetry' in params: mshCrck.addBoundaries([params['Symmetry']]) if 'WakeTips' in params: - mshCrck.setExcluded(params['WakeTips'][i]) # 2.5D computations + mshCrck.setExcluded(params['WakeTips'][i]) # 3D computations (not excluded for 2.5D) mshCrck.run() # save the updated mesh in native (gmsh) format and then set the writer to the requested format tbox.GmshExport(_msh).save(_msh.name) @@ -191,12 +191,19 @@ def initDart(params, scenario='aerodynamic', task='analysis'): bnd = dart.Body(_msh, [params['Fuselage'], params['Fluid']]) _pbl.add(bnd) # add Wake/Kutta boundary conditions + # TODO refactor Wake "exclusion" for 3D? if _dim == 2: _pbl.add(dart.Wake(_msh, [params['Wake'], params['Wake']+'_', params['Fluid']])) _pbl.add(dart.Kutta(_msh, [params['Te'], params['Wake']+'_', params['Wing'], params['Fluid']])) else: for i in range(len(params['Wakes'])): - _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['TeTips'][i]])) + if 'WakeExs' in params: + _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['WakeExs'][i]])) # 3D + fuselage + elif 'WakeTips' in params: + _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['WakeTips'][i]])) # 3D + else: + _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid']])) # 2.5D + _pbl.add(dart.Kutta(_msh, [params['Tes'][i], params['Wakes'][i]+'_', params['Wings'][i], params['Fluid']])) # Direct (forward) solver creation # NB: more solvers/options are available but we restrict the user's choice to the most efficient ones diff --git a/dart/api/internal/trim.py b/dart/api/internal/trim.py index 2555218013c5e20a211ed44fe66de850b985924a..be1c330c4935622140a81068ffbb83b6efc182f7 100644 --- a/dart/api/internal/trim.py +++ b/dart/api/internal/trim.py @@ -111,7 +111,7 @@ class Trim: self.sol.save(self.wrtr) # extract Cp if self.dim == 2: - Cp = dU.extract(self.bnd.groups[0].tag.elems, self.solver.Cp) + Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp) tU.write(Cp, "Cp_airfoil.dat", "%1.5e", ", ", "x, y, z, Cp", "") elif self.dim == 3 and self.format == 'vtk' and self.slice: dU.writeSlices(self.msh.name, self.slice, self.tag) diff --git a/dart/benchmark/onera.py b/dart/benchmark/onera.py index 08e2f85e3d40dcfe2b3e21696f871645d320cf8f..75209e232fcbc4c17ffc551834b2590f87be87cf 100644 --- a/dart/benchmark/onera.py +++ b/dart/benchmark/onera.py @@ -66,7 +66,7 @@ def main(): # set the problem and add fluid, initial/boundary conditions tms['pre'].start() - pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') + pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip') tms['pre'].stop() # solve problem diff --git a/dart/benchmark/onera_lfs.msh b/dart/benchmark/onera_lfs.msh index bad883483fd52e1bc2012df01111620350a1d674..b69a99789cfecd7d3db417b14288b4271dc85713 100644 --- a/dart/benchmark/onera_lfs.msh +++ b/dart/benchmark/onera_lfs.msh @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:da782b51f3be6ff2b8d222ed9ad1aba46849b2d018df294ac5f780a0258a74c8 -size 39450759 +oid sha256:c235746b196f5b428da900f56a024eb69c5d1baa2766409c0e1413109c2cd7fb +size 36951138 diff --git a/dart/cases/coyote.py b/dart/cases/coyote.py index 6dda2e0d08de2e3f96136091afd2e83cb368361d..fc493e30324286c0bd5cf268fcb1b797bd2bda72 100644 --- a/dart/cases/coyote.py +++ b/dart/cases/coyote.py @@ -44,7 +44,7 @@ def getParam(): p['Dim'] = 3 # Problem dimension p['Format'] = 'vtk' # Save format (vtk or gmsh) p['Slice'] = [1.8, 3.5, 5.3] # array of (y) coordinates to perform slice along the span (empty if none) - p['TagId'] = 9 # tag number of physical group to be sliced + p['TagId'] = 11 # tag number of physical group to be sliced # Markers p['Fluid'] = 'field' # Name of physical group containing the fluid p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element) @@ -52,7 +52,8 @@ def getParam(): p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface body 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 + p['WakeExs'] = ['wakeExW', 'wakeExT'] # LIST of names of physical group containing the edge of the wake and the intersection of lifting surface with fuselage (to be excluded from Wake B.C.) + p['Tes'] = ['teW', 'teT'] # LIST of names of physical group containing the trailing edge # Freestream p['M_inf'] = 0.8 # Freestream Mach number p['AoA_begin'] = 1 # Freestream angle of attack (begin) diff --git a/dart/cases/n0012_3.py b/dart/cases/n0012_3.py deleted file mode 100644 index 57700c2e6cbe4b92065fcc678300cef2369fb776..0000000000000000000000000000000000000000 --- a/dart/cases/n0012_3.py +++ /dev/null @@ -1,73 +0,0 @@ -#!/usr/bin/env python3 -# -*- 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. - - -## NACA 0012 wing configuration file for flow polar script -# Adrien Crovato - -def getParam(): - import os.path - from fwk.wutils import parseargs - p = {} - # Specific - span = 1. # span length - # Arguments - args = parseargs() - p['Threads'] = args.k - p['Verb'] = args.verb - # Input/Output - p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n0012_3.geo') # Input file containing the model - p['Pars'] = {'spn' : span, 'lgt' : 3., 'wdt' : 3., 'hgt' : 3., 'msLe' : 0.02, 'msTe' : 0.02, 'msF' : 1.} # Parameters for input file model - p['Dim'] = 3 # Problem dimension - p['Format'] = 'vtk' # Save format (vtk or gmsh) - p['Slice'] = [0.01, 0.25, 0.5, 0.75, 0.95] # array of (y) coordinates to perform slice along the span (empty if none) - p['TagId'] = 4 # tag number of physical group to be sliced - # Markers - p['Fluid'] = 'field' # Name of physical group containing the fluid - p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of name of physical groups containing the farfield boundaries (upstream/downstream should be first/last element) - p['Symmetry'] = 'symmetry' # Name of physical group containing the symmetry boundaries - p['Wings'] = ['wing'] # LIST of names of physical groups containing the lifting surface body - p['Wakes'] = ['wake'] # LIST of names of physical group containing the wake - p['WakeTips'] = ['wakeTip'] # LIST of names of physical group containing the edge of the wake - p['TeTips'] = ['teTip'] # LIST of names of physical group containing the edge of the wake and the trailing edge - # Freestream - p['M_inf'] = 0.3 # Freestream Mach number - p['AoA_begin'] = -3. # Freestream angle of attack (begin) - p['AoA_end'] = 3. # Freestream angle of attack (end) - p['AoA_step'] = 1. # Freestream angle of attack (step) - # Geometry - p['S_ref'] = 1.*span # Reference surface length - p['c_ref'] = 1. # Reference chord length - p['x_ref'] = 0.25 # Reference point for moment computation (x) - p['y_ref'] = 0. # Reference point for moment computation (y) - p['z_ref'] = 0. # Reference point for moment computation (z) - # Numerical - p['LSolver'] = 'MUMPS' # Linear solver (PARDISO, GMRES, or MUMPS) - p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual - p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual - p['Max_it'] = 10 # Solver maximum number of iterations - return p - -def main(): - # Script call - import dart.api.internal.polar as p - polar = p.Polar(getParam()) - polar.run() - polar.disp() - -if __name__ == "__main__": - main() diff --git a/dart/cases/wbht.py b/dart/cases/wbht.py index 4744d14e7fac451e54428a842c6e4b6e37acc4ee..73a8f1daaa682aca82263d5d7027baf1f2ecb3d8 100644 --- a/dart/cases/wbht.py +++ b/dart/cases/wbht.py @@ -56,7 +56,7 @@ def getParam(): p['Dim'] = 3 # Problem dimension p['Format'] = 'vtk' # Save format (vtk or gmsh) 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'] = 10 # tag number of physical group to be sliced + p['TagId'] = 12 # tag number of physical group to be sliced # Markers p['Fluid'] = 'field' # Name of physical group containing the fluid p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element) @@ -65,9 +65,10 @@ def getParam(): p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface body 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 + p['WakeExs'] = ['wakeExW', 'wakeExT'] # LIST of names of physical group containing the edge of the wake and the intersection of lifting surface with fuselage (to be excluded from Wake B.C.) + p['Tes'] = ['teW', 'teT'] # LIST of names of physical group containing the trailing edge # Freestream - p['M_inf'] = 0.75 # Freestream Mach number + p['M_inf'] = 0.72 # 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) diff --git a/dart/cases/wht.py b/dart/cases/wht.py index ce7301700cb9d5967ab2f500b323740fe02a7d6d..dd0fc7d10496d2bdafe0eacdc246bbcf7445c00e 100644 --- a/dart/cases/wht.py +++ b/dart/cases/wht.py @@ -59,9 +59,9 @@ def getParam(): p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface body 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 + p['Tes'] = ['teW', 'teT'] # LIST of names of physical group containing the trailing edge # Freestream - p['M_inf'] = 0.75 # Freestream Mach number + p['M_inf'] = 0.72 # 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) diff --git a/dart/default.py b/dart/default.py index 3a03df5626170c58ef3be646579d518954072a69..93cb86c6e6acdee97831fc74d0af51b802155767 100644 --- a/dart/default.py +++ b/dart/default.py @@ -57,7 +57,7 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None): gmshWriter.save(msh.name) return msh, gmshWriter -def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld = 'field', upstr = 'upstream', ff = 'farfield', dnstr = 'downstream', wk = 'wake', te = None, tp = None, vsc = False, dbc = False): +def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld = 'field', upstr = 'upstream', ff = 'farfield', dnstr = 'downstream', wk = 'wake', te = 'te', tp = None, vsc = False, dbc = False): """Initialize problem Parameters @@ -97,7 +97,7 @@ def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld te : str (optional) name of physical tag contaning the trailing edge tp : str (optional) - name of physical tag contaning the wake tip (and the trailing edge for 3D) + name of physical tag contaning the wake tip (only for 3D) vsc : str (optional) name of physical tag contaning the body to apply transpiration B.C. dbc : bool (optional) @@ -121,15 +121,12 @@ def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld pbl.add(dart.Freestream(msh, ff, dim, alpha, beta)) pbl.add(dart.Freestream(msh, dnstr, dim, alpha, beta)) # add wake and Kutta conditions - if dim == 2 and te is not None: - wake = dart.Wake(msh, [wk, wk+'_', fld]) - pbl.add(wake) - pbl.add(dart.Kutta(msh, [te, wk+'_', sur, fld])) - elif dim == 3: - wake = dart.Wake(msh, [wk, wk+'_', fld, tp]) - pbl.add(wake) + if tp is None: + wake = dart.Wake(msh, [wk, wk+'_', fld]) # 2 and 2.5D else: - wake = None + wake = dart.Wake(msh, [wk, wk+'_', fld, tp]) # 3D + pbl.add(wake) + pbl.add(dart.Kutta(msh, [te, wk+'_', sur, fld])) # add body of interest bnd = dart.Body(msh, [sur, fld]) pbl.add(bnd) diff --git a/dart/models/agard445.geo b/dart/models/agard445.geo index c7007513110879f3906df863a6c26e1f16dcb322..d26866470a8cb6b02b872b5cb6494022d109d0ff 100644 --- a/dart/models/agard445.geo +++ b/dart/models/agard445.geo @@ -354,7 +354,7 @@ Mesh.SmoothNormals = 1; // Trailing edge and wake tip Physical Line("wakeTip") = {162}; -Physical Line("teTip") = {61, 162}; +Physical Line("te") = {61}; // Internal Field: Physical Volume("field") = {1}; diff --git a/dart/models/coyote.geo b/dart/models/coyote.geo index c9e3615f9aa4efa89c7409c58c96e60ef3679703..bdf9a0cc9fe4ec222421b044ec6bd1a3b1ee3d53 100644 --- a/dart/models/coyote.geo +++ b/dart/models/coyote.geo @@ -61,9 +61,11 @@ Coherence; /* Physical */ // Tips -Physical Line("teTipW") = {155:158, 160:161, 164:167, 169:170}; +Physical Line("wakeExW") = {156:158, 160, 165:167, 169}; +Physical Line("teW") = {155, 161, 164, 170}; Physical Line("wakeTipW") = {160, 169}; -Physical Line("teTipT") = {173:176, 178:181}; +Physical Line("wakeExT") = {173, 175, 176, 178, 180, 181}; +Physical Line("teT") = {174, 179}; Physical Line("wakeTipT") = {173, 178}; // Downstream diff --git a/dart/models/cylinder_2D5.geo b/dart/models/cylinder_2D5.geo deleted file mode 100644 index cc70813611ef7e12c5d07f74371e69da126c9cf4..0000000000000000000000000000000000000000 --- a/dart/models/cylinder_2D5.geo +++ /dev/null @@ -1,161 +0,0 @@ -/* 2D5 Cylinder */ - -// Parameters -// domain -DefineConstant[ dmt = { 1.0, Name "Cylinder diameter" } ]; -DefineConstant[ lgt = { 3.0, Name "Channel length" } ]; -DefineConstant[ wdt = { 2.0, Name "Channel width" } ]; -DefineConstant[ hgt = { 3.0, Name "Channel height" } ]; -DefineConstant[ msN = { 0.1, Min 0.1, Max 1, Step 0.1, Name "Nearfield mesh size" } ]; -DefineConstant[ msF = { 1.0, Min 0.1, Max 1, Step 0.1, Name "Farfield mesh size" } ]; - -// mesh algo -msh3D = 2; -msh2D = 5; - -// Geometry -// cylinder -Point(0) = {0, 0, 0}; -Point(1) = {dmt/2, 0, 0, msN}; -Point(2) = {0, 0, dmt/2, msN}; -Point(3) = {-dmt/2, 0, 0, msN}; -Point(4) = {0, 0, -dmt/2, msN}; -Point(10) = {0, wdt, 0}; -Point(11) = {dmt/2, wdt, 0, msN}; -Point(12) = {0, wdt, dmt/2, msN}; -Point(13) = {-dmt/2, wdt, 0, msN}; -Point(14) = {0, wdt, -dmt/2, msN}; - -// box -Point(31) = {lgt/2, 0, hgt/2, msF}; -Point(32) = {-lgt/2, 0, hgt/2, msF}; -Point(33) = {-lgt/2, 0, -hgt/2, msF}; -Point(34) = {lgt/2, 0, -hgt/2, msF}; -Point(35) = {lgt/2, wdt, hgt/2, msF}; -Point(36) = {-lgt/2, wdt, hgt/2, msF}; -Point(37) = {-lgt/2, wdt, -hgt/2, msF}; -Point(38) = {lgt/2, wdt, -hgt/2, msF}; - -// midplane -Point(41) = {lgt/2, 0, 0, msF}; -Point(42) = {lgt/2, wdt, 0, msF}; -Point(43) = {-lgt/2, wdt, 0, msF}; -Point(44) = {-lgt/2, 0, 0, msF}; - -// wing -Circle(1) = {1, 0, 2}; -Circle(2) = {2, 0, 3}; -Circle(3) = {3, 0, 4}; -Circle(4) = {4, 0, 1}; -Circle(5) = {11, 10, 12}; -Circle(6) = {12, 10, 13}; -Circle(7) = {13, 10, 14}; -Circle(8) = {14, 10, 11}; -Line(9) = {1, 11}; -Line(10) = {2, 12}; -Line(11) = {3, 13}; -Line(12) = {4, 14}; - -// box -Line(21) = {41, 31}; -Line(22) = {31, 32}; -Line(23) = {32, 44}; -Line(24) = {44, 33}; -Line(25) = {33, 34}; -Line(26) = {34, 41}; -Line(27) = {42, 35}; -Line(28) = {35, 36}; -Line(29) = {36, 43}; -Line(30) = {43, 37}; -Line(31) = {37, 38}; -Line(32) = {38, 42}; -Line(33) = {31, 35}; -Line(34) = {32, 36}; -Line(35) = {33, 37}; -Line(36) = {34, 38}; - -// midplane -Line(41) = {1, 41}; -Line(42) = {11, 42}; -Line(43) = {13, 43}; -Line(44) = {3, 44}; -Line(45) = {41, 42}; -Line(46) = {42, 11}; -Line(47) = {13, 43}; -Line(48) = {43, 44}; - -// wing -Line Loop(1) = {1, 10, -5, -9}; -Surface(1) = {-1}; -Line Loop(2) = {2, 11, -6, -10}; -Surface(2) = {-2}; -Line Loop(3) = {3, 12, -7, -11}; -Surface(3) = {-3}; -Line Loop(4) = {4, 9, -8, -12}; -Surface(4) = {-4}; - -// symmetry -Line Loop(5) = {41, 21, 22, 23, -44, -2, -1}; -Plane Surface(5) = {-5}; -Line Loop(12) = {42, 27, 28, 29, -43, -6, -5}; -Plane Surface(12) = {12}; -Line Loop(6) = {44, 24, 25, 26, -41, -4, -3}; -Plane Surface(6) = {-6}; -Line Loop(13) = {43, 30, 31, 32, -42, -8, -7}; -Plane Surface(13) = {13}; - -// upstream -Line Loop(7) = {48, -23, 34, 29}; -Plane Surface(7) = {-7}; -Line Loop(8) = {48, 24, 35, -30}; -Plane Surface(8) = {8}; - -// downstream -Line Loop(9) = {45, 27, -33, -21}; -Plane Surface(9) = {-9}; -Line Loop(10) = {45, -32, -36, 26}; -Plane Surface(10) = {10}; - -// farfield -Line Loop(11) = {22, 34, -28, -33}; -Plane Surface(11) = {11}; -Line Loop(14) = {25, 36, -31, -35}; -Plane Surface(14) = {14}; - -// wake -Line Loop(15) = {41, 45, -42, -9}; -Plane Surface(15) = {15}; - -// internal -Line Loop(16) = {11, 43, 48, -44}; -Plane Surface(16) = {16}; - -// field -Surface Loop(1) = {7, 5, 9, 12, 11, 2, 1, 16, 15}; -Volume(1) = {1}; -Surface Loop(2) = {6, 8, 14, 10, 13, 4, 3, 16, 15}; -Volume(2) = {2}; - -// Mesh -// 2D -MeshAlgorithm Surface{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} = msh2D; -// 3D -Mesh.Algorithm3D = msh3D; -Mesh.OptimizeNetgen = 1; -Mesh.Smoothing = 10; -Mesh.SmoothNormals = 1; - -// Physical -Physical Line("te") = {9}; -Physical Surface("cylinder") = {1, 2}; -Physical Surface("cylinder_") = {3, 4}; -Physical Surface("symmetry") = {5, 12}; -Physical Surface("symmetry_") = {6, 13}; -Physical Surface("upstream") = {7, 8}; -Physical Surface("downstream") = {9}; -Physical Surface("downstream_") = {10}; -Physical Surface("farfield") = {11, 14}; -Physical Surface("wake") = {15}; -Physical Volume("field") = {1}; -Physical Volume("field_") = {2}; - diff --git a/dart/models/cylinder_3.geo b/dart/models/cylinder_3.geo deleted file mode 100644 index f5c9602503fe6c07e9c567eee05ed1ff85f5b568..0000000000000000000000000000000000000000 --- a/dart/models/cylinder_3.geo +++ /dev/null @@ -1,203 +0,0 @@ -/* 3D Cylinder */ - -// Parameters -// domain -DefineConstant[ dmt = { 1.0, Name "Cylinder diameter" } ]; -DefineConstant[ spn = { 2.5, Name "Cylinder span" } ]; -DefineConstant[ tpt = { .05, Name "Cylinder tip thickness" } ]; // percentage of span -DefineConstant[ lgt = { 8.0, Name "Channel length" } ]; -DefineConstant[ wdt = { 4.0, Name "Channel width" } ]; -DefineConstant[ hgt = { 6.0, Name "Channel height" } ]; -DefineConstant[ msN = { 0.1, Min 0.1, Max 1, Step 0.1, Name "Nearfield mesh size" } ]; -DefineConstant[ msF = { 1.0, Min 0.1, Max 1, Step 0.1, Name "Farfield mesh size" } ]; - -// mesh algo -msh3D = 2; -msh2D = 5; - -// Geometry -// cylinder -Point(0) = {0, 0, 0}; -Point(1) = {dmt/2, 0, 0, msN}; -Point(2) = {0, 0, dmt/2, msN}; -Point(3) = {-dmt/2, 0, 0, msN}; -Point(4) = {0, 0, -dmt/2, msN}; -Point(10) = {0, spn, 0}; -Point(11) = {dmt/2, spn, 0, msN}; -Point(12) = {0, spn, dmt/2, msN}; -Point(13) = {-dmt/2, spn, 0, msN}; -Point(14) = {0, spn, -dmt/2, msN}; - -// tip -Point(21) = {0, spn+tpt*spn, 0, msN}; - -// box -Point(31) = {lgt/2, 0, hgt/2, msF}; -Point(32) = {-lgt/2, 0, hgt/2, msF}; -Point(33) = {-lgt/2, 0, -hgt/2, msF}; -Point(34) = {lgt/2, 0, -hgt/2, msF}; -Point(35) = {lgt/2, wdt, hgt/2, msF}; -Point(36) = {-lgt/2, wdt, hgt/2, msF}; -Point(37) = {-lgt/2, wdt, -hgt/2, msF}; -Point(38) = {lgt/2, wdt, -hgt/2, msF}; - -// midplane -Point(41) = {lgt/2, 0, 0, msF}; -Point(42) = {lgt/2, spn, 0, msF}; -Point(43) = {lgt/2, wdt, 0, msF}; -Point(44) = {dmt/2, wdt, 0, msF}; -Point(45) = {0, wdt, 0, msF}; -Point(46) = {-dmt/2, wdt, 0, msF}; -Point(47) = {-lgt/2, wdt, 0, msF}; -Point(48) = {-lgt/2, spn, 0, msF}; -Point(49) = {-lgt/2, 0, 0, msF}; - -// wing -Circle(1) = {1, 0, 2}; -Circle(2) = {2, 0, 3}; -Circle(3) = {3, 0, 4}; -Circle(4) = {4, 0, 1}; -Circle(11) = {11, 10, 12}; -Circle(12) = {12, 10, 13}; -Circle(13) = {13, 10, 14}; -Circle(14) = {14, 10, 11}; -Line(21) = {1, 11}; -Line(22) = {2, 12}; -Line(23) = {3, 13}; -Line(24) = {4, 14}; - -// tip -Ellipse(31) = {11, 10, 10, 21}; -Ellipse(32) = {12, 10, 10, 21}; -Ellipse(33) = {13, 10, 10, 21}; -Ellipse(34) = {14, 10, 10, 21}; - -// box -Line(41) = {41, 31}; -Line(42) = {31, 32}; -Line(43) = {32, 49}; -Line(44) = {49, 33}; -Line(45) = {33, 34}; -Line(46) = {34, 41}; -Line(47) = {43, 35}; -Line(48) = {35, 36}; -Line(49) = {36, 47}; -Line(50) = {47, 37}; -Line(51) = {37, 38}; -Line(52) = {38, 43}; -Line(53) = {31, 35}; -Line(54) = {32, 36}; -Line(55) = {33, 37}; -Line(56) = {34, 38}; - -// midplane -Line(61) = {1, 41}; -Line(62) = {11, 42}; -Line(63) = {11, 44}; -Line(64) = {21, 45}; -Line(65) = {13, 46}; -Line(66) = {13, 48}; -Line(67) = {3, 49}; -Line(68) = {41, 42}; -Line(69) = {42, 43}; -Line(70) = {43, 44}; -Line(71) = {44, 45}; -Line(72) = {45, 46}; -Line(73) = {46, 47}; -Line(74) = {47, 48}; -Line(75) = {48, 49}; - -// cylinder -Line Loop(1) = {1, 22, -11, -21}; -Surface(1) = {-1}; -Line Loop(2) = {2, 23, -12, -22}; -Surface(2) = {-2}; -Line Loop(3) = {3, 24, -13, -23}; -Surface(3) = {-3}; -Line Loop(4) = {4, 21, -14, -24}; -Surface(4) = {-4}; -//tip -Line Loop(5) = {11, 32, -31}; -Surface(5) = {-5}; -Line Loop(6) = {12, 33, -32}; -Surface(6) = {-6}; -Line Loop(7) = {13, 34, -33}; -Surface(7) = {-7}; -Line Loop(8) = {14, 31, -34}; -Surface(8) = {-8}; - -// symmetry -Line Loop(9) = {61, 41, 42, 43, -67, -2, -1}; -Plane Surface(9) = {-9}; -Line Loop(10) = {4, 61, -46, -45, -44, -67, 3}; -Plane Surface(10) = {10}; - -// upstream -Line Loop(11) = {43, -75, -74, -49, -54}; -Plane Surface(11) = {11}; -Line Loop(12) = {75, 44, 55, -50, 74}; -Plane Surface(12) = {12}; - -// downstream -Line Loop(13) = {41, 53, -47, -69, -68}; -Plane Surface(13) = {13}; -Line Loop(14) = {68, 69, -52, -56, 46}; -Plane Surface(14) = {14}; - -// farfield -Line Loop(15) = {42, 54, -48, -53}; -Plane Surface(15) = {15}; -Line Loop(16) = {47, 48, 49, -73, -72, -71, -70}; -Plane Surface(16) = {16}; -Line Loop(17) = {51, 52, 70, 71, 72, 73, 50}; -Plane Surface(17) = {17}; -Line Loop(18) = {45, 56, -51, -55}; -Plane Surface(18) = {18}; - -// wake -Line Loop(19) = {61, 68, -62, -21}; -Surface(19) = {19}; - -// internal -Line Loop(20) = {62, 69, 70, -63}; -Surface(20) = {20}; -Line Loop(21) = {64, -71, -63, 31}; -Surface(21) = {-21}; -Line Loop(22) = {65, -72, -64, -33}; -Surface(22) = {-22}; -Line Loop(23) = {74, -66, 65, 73}; -Surface(23) = {23}; -Line Loop(24) = {75, -67, 23, 66}; -Surface(24) = {24}; - -// field -Surface Loop(1) = {15, 9, 13, 16, 11, 2, 6, 5, 1, 24, 23, 22, 21, 20, 19}; -Volume(1) = {1}; -Surface Loop(2) = {10, 4, 8, 7, 3, 14, 17, 18, 12, 24, 23, 22, 21, 20, 19}; -Volume(2) = {2}; - -// Mesh -// cylinder, farfield, symmetry, midplane -MeshAlgorithm Surface{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24} = msh2D; -// 3D -Mesh.Algorithm3D = msh3D; -Mesh.OptimizeNetgen = 1; -Mesh.Smoothing = 10; -Mesh.SmoothNormals = 1; - -// Physical -Physical Line("teTip") = {21, 62}; -Physical Line("wakeTip") = {62}; -Physical Surface("cylinder") = {1, 2}; -Physical Surface("cylinder_") = {3, 4}; -Physical Surface("tip") = {5, 6, 7, 8}; -Physical Surface("symmetry") = {9}; -Physical Surface("symmetry_") = {10}; -Physical Surface("upstream") = {11, 12}; -Physical Surface("downstream") = {13}; -Physical Surface("downstream_") = {14}; -Physical Surface("farfield") = {15, 16, 17, 18}; -Physical Surface("wake") = {19}; -Physical Volume("field") = {1}; -Physical Volume("field_") = {2}; - diff --git a/dart/models/n0012_3.geo b/dart/models/n0012_3.geo deleted file mode 100644 index c8b02ca41949324232d062a0b5289585c7198591..0000000000000000000000000000000000000000 --- a/dart/models/n0012_3.geo +++ /dev/null @@ -1,442 +0,0 @@ -/* Rectangular NACA0012 wing */ -// Initially generated by unsRgridWingGen.m -// For gmsh 4, use line 334 instead of line 335 (Bezier <- BSpline) - -// Parameters -// domain and mesh -DefineConstant[ spn = { 2.0, Name "Wing span" } ]; -DefineConstant[ lgt = { 3.0, Name "Channel length" } ]; -DefineConstant[ wdt = { 3.0, Name "Channel width" } ]; -DefineConstant[ hgt = { 3.0, Name "Channel height" } ]; -DefineConstant[ msLe = { 0.05, Name "Leading edge mesh size" } ]; -DefineConstant[ msTe = { 0.05, Name "Trailing edge mesh size" } ]; -DefineConstant[ msF = { 1.0, Name "Farfield mesh size" } ]; - -//// GEOMETRY - - -/// Points - -// Airfoil 1: naca0012, 101 points -Point(1) = {1.000000,0.000000,0.000000,msTe}; -Point(2) = {0.999013,0.000000,0.000143}; -Point(3) = {0.996057,0.000000,0.000572}; -Point(4) = {0.991144,0.000000,0.001280}; -Point(5) = {0.984292,0.000000,0.002260}; -Point(6) = {0.975528,0.000000,0.003501}; -Point(7) = {0.964888,0.000000,0.004990}; -Point(8) = {0.952414,0.000000,0.006710}; -Point(9) = {0.938153,0.000000,0.008643}; -Point(10) = {0.922164,0.000000,0.010770}; -Point(11) = {0.904508,0.000000,0.013071,msTe}; -Point(12) = {0.885257,0.000000,0.015523}; -Point(13) = {0.864484,0.000000,0.018106}; -Point(14) = {0.842274,0.000000,0.020795}; -Point(15) = {0.818712,0.000000,0.023569}; -Point(16) = {0.793893,0.000000,0.026405}; -Point(17) = {0.767913,0.000000,0.029279}; -Point(18) = {0.740877,0.000000,0.032168}; -Point(19) = {0.712890,0.000000,0.035048}; -Point(20) = {0.684062,0.000000,0.037896}; -Point(21) = {0.654508,0.000000,0.040686}; -Point(22) = {0.624345,0.000000,0.043394}; -Point(23) = {0.593691,0.000000,0.045992}; -Point(24) = {0.562667,0.000000,0.048455}; -Point(25) = {0.531395,0.000000,0.050754}; -Point(26) = {0.500000,0.000000,0.052862}; -Point(27) = {0.468605,0.000000,0.054749}; -Point(28) = {0.437333,0.000000,0.056390}; -Point(29) = {0.406309,0.000000,0.057755}; -Point(30) = {0.375655,0.000000,0.058819}; -Point(31) = {0.345492,0.000000,0.059557}; -Point(32) = {0.315938,0.000000,0.059947}; -Point(33) = {0.287110,0.000000,0.059971,msTe}; -Point(34) = {0.259123,0.000000,0.059614}; -Point(35) = {0.232087,0.000000,0.058863}; -Point(36) = {0.206107,0.000000,0.057712}; -Point(37) = {0.181288,0.000000,0.056159}; -Point(38) = {0.157726,0.000000,0.054206}; -Point(39) = {0.135516,0.000000,0.051862}; -Point(40) = {0.114743,0.000000,0.049138}; -Point(41) = {0.095492,0.000000,0.046049}; -Point(42) = {0.077836,0.000000,0.042615}; -Point(43) = {0.061847,0.000000,0.038859}; -Point(44) = {0.047586,0.000000,0.034803}; -Point(45) = {0.035112,0.000000,0.030473}; -Point(46) = {0.024472,0.000000,0.025893}; -Point(47) = {0.015708,0.000000,0.021088}; -Point(48) = {0.008856,0.000000,0.016078}; -Point(49) = {0.003943,0.000000,0.010884}; -Point(50) = {0.000987,0.000000,0.005521}; -Point(51) = {0.000000,0.000000,0.000000,msLe}; -Point(52) = {0.000987,0.000000,-0.005521}; -Point(53) = {0.003943,0.000000,-0.010884}; -Point(54) = {0.008856,0.000000,-0.016078}; -Point(55) = {0.015708,0.000000,-0.021088}; -Point(56) = {0.024472,0.000000,-0.025893}; -Point(57) = {0.035112,0.000000,-0.030473}; -Point(58) = {0.047586,0.000000,-0.034803}; -Point(59) = {0.061847,0.000000,-0.038859}; -Point(60) = {0.077836,0.000000,-0.042615}; -Point(61) = {0.095492,0.000000,-0.046049}; -Point(62) = {0.114743,0.000000,-0.049138}; -Point(63) = {0.135516,0.000000,-0.051862}; -Point(64) = {0.157726,0.000000,-0.054206}; -Point(65) = {0.181288,0.000000,-0.056159}; -Point(66) = {0.206107,0.000000,-0.057712}; -Point(67) = {0.232087,0.000000,-0.058863}; -Point(68) = {0.259123,0.000000,-0.059614}; -Point(69) = {0.287110,0.000000,-0.059971,msTe}; -Point(70) = {0.315938,0.000000,-0.059947}; -Point(71) = {0.345492,0.000000,-0.059557}; -Point(72) = {0.375655,0.000000,-0.058819}; -Point(73) = {0.406309,0.000000,-0.057755}; -Point(74) = {0.437333,0.000000,-0.056390}; -Point(75) = {0.468605,0.000000,-0.054749}; -Point(76) = {0.500000,0.000000,-0.052862}; -Point(77) = {0.531395,0.000000,-0.050754}; -Point(78) = {0.562667,0.000000,-0.048455}; -Point(79) = {0.593691,0.000000,-0.045992}; -Point(80) = {0.624345,0.000000,-0.043394}; -Point(81) = {0.654508,0.000000,-0.040686}; -Point(82) = {0.684062,0.000000,-0.037896}; -Point(83) = {0.712890,0.000000,-0.035048}; -Point(84) = {0.740877,0.000000,-0.032168}; -Point(85) = {0.767913,0.000000,-0.029279}; -Point(86) = {0.793893,0.000000,-0.026405}; -Point(87) = {0.818712,0.000000,-0.023569}; -Point(88) = {0.842274,0.000000,-0.020795}; -Point(89) = {0.864484,0.000000,-0.018106}; -Point(90) = {0.885257,0.000000,-0.015523}; -Point(91) = {0.904508,0.000000,-0.013071,msTe}; -Point(92) = {0.922164,0.000000,-0.010770}; -Point(93) = {0.938153,0.000000,-0.008643}; -Point(94) = {0.952414,0.000000,-0.006710}; -Point(95) = {0.964888,0.000000,-0.004990}; -Point(96) = {0.975528,0.000000,-0.003501}; -Point(97) = {0.984292,0.000000,-0.002260}; -Point(98) = {0.991144,0.000000,-0.001280}; -Point(99) = {0.996057,0.000000,-0.000572}; -Point(100) = {0.999013,0.000000,-0.000143,msTe}; - -// Airfoil 2: naca0012, 101 points -Point(102) = {1.000000,spn,0.000000,msTe}; -Point(103) = {0.999013,spn,0.000143}; -Point(104) = {0.996057,spn,0.000572}; -Point(105) = {0.991144,spn,0.001280}; -Point(106) = {0.984292,spn,0.002260}; -Point(107) = {0.975528,spn,0.003501}; -Point(108) = {0.964888,spn,0.004990}; -Point(109) = {0.952414,spn,0.006710}; -Point(110) = {0.938153,spn,0.008643}; -Point(111) = {0.922164,spn,0.010770}; -Point(112) = {0.904508,spn,0.013071,msTe}; -Point(113) = {0.885257,spn,0.015523}; -Point(114) = {0.864484,spn,0.018106}; -Point(115) = {0.842274,spn,0.020795}; -Point(116) = {0.818712,spn,0.023569}; -Point(117) = {0.793893,spn,0.026405}; -Point(118) = {0.767913,spn,0.029279}; -Point(119) = {0.740877,spn,0.032168}; -Point(120) = {0.712890,spn,0.035048}; -Point(121) = {0.684062,spn,0.037896}; -Point(122) = {0.654508,spn,0.040686}; -Point(123) = {0.624345,spn,0.043394}; -Point(124) = {0.593691,spn,0.045992}; -Point(125) = {0.562667,spn,0.048455}; -Point(126) = {0.531395,spn,0.050754}; -Point(127) = {0.500000,spn,0.052862}; -Point(128) = {0.468605,spn,0.054749}; -Point(129) = {0.437333,spn,0.056390}; -Point(130) = {0.406309,spn,0.057755}; -Point(131) = {0.375655,spn,0.058819}; -Point(132) = {0.345492,spn,0.059557}; -Point(133) = {0.315938,spn,0.059947}; -Point(134) = {0.287110,spn,0.059971,msTe}; -Point(135) = {0.259123,spn,0.059614}; -Point(136) = {0.232087,spn,0.058863}; -Point(137) = {0.206107,spn,0.057712}; -Point(138) = {0.181288,spn,0.056159}; -Point(139) = {0.157726,spn,0.054206}; -Point(140) = {0.135516,spn,0.051862}; -Point(141) = {0.114743,spn,0.049138}; -Point(142) = {0.095492,spn,0.046049}; -Point(143) = {0.077836,spn,0.042615}; -Point(144) = {0.061847,spn,0.038859}; -Point(145) = {0.047586,spn,0.034803}; -Point(146) = {0.035112,spn,0.030473}; -Point(147) = {0.024472,spn,0.025893}; -Point(148) = {0.015708,spn,0.021088}; -Point(149) = {0.008856,spn,0.016078}; -Point(150) = {0.003943,spn,0.010884}; -Point(151) = {0.000987,spn,0.005521}; -Point(152) = {0.000000,spn,0.000000,msLe}; -Point(153) = {0.000987,spn,-0.005521}; -Point(154) = {0.003943,spn,-0.010884}; -Point(155) = {0.008856,spn,-0.016078}; -Point(156) = {0.015708,spn,-0.021088}; -Point(157) = {0.024472,spn,-0.025893}; -Point(158) = {0.035112,spn,-0.030473}; -Point(159) = {0.047586,spn,-0.034803}; -Point(160) = {0.061847,spn,-0.038859}; -Point(161) = {0.077836,spn,-0.042615}; -Point(162) = {0.095492,spn,-0.046049}; -Point(163) = {0.114743,spn,-0.049138}; -Point(164) = {0.135516,spn,-0.051862}; -Point(165) = {0.157726,spn,-0.054206}; -Point(166) = {0.181288,spn,-0.056159}; -Point(167) = {0.206107,spn,-0.057712}; -Point(168) = {0.232087,spn,-0.058863}; -Point(169) = {0.259123,spn,-0.059614}; -Point(170) = {0.287110,spn,-0.059971,msTe}; -Point(171) = {0.315938,spn,-0.059947}; -Point(172) = {0.345492,spn,-0.059557}; -Point(173) = {0.375655,spn,-0.058819}; -Point(174) = {0.406309,spn,-0.057755}; -Point(175) = {0.437333,spn,-0.056390}; -Point(176) = {0.468605,spn,-0.054749}; -Point(177) = {0.500000,spn,-0.052862}; -Point(178) = {0.531395,spn,-0.050754}; -Point(179) = {0.562667,spn,-0.048455}; -Point(180) = {0.593691,spn,-0.045992}; -Point(181) = {0.624345,spn,-0.043394}; -Point(182) = {0.654508,spn,-0.040686}; -Point(183) = {0.684062,spn,-0.037896}; -Point(184) = {0.712890,spn,-0.035048}; -Point(185) = {0.740877,spn,-0.032168}; -Point(186) = {0.767913,spn,-0.029279}; -Point(187) = {0.793893,spn,-0.026405}; -Point(188) = {0.818712,spn,-0.023569}; -Point(189) = {0.842274,spn,-0.020795}; -Point(190) = {0.864484,spn,-0.018106}; -Point(191) = {0.885257,spn,-0.015523}; -Point(192) = {0.904508,spn,-0.013071,msTe}; -Point(193) = {0.922164,spn,-0.010770}; -Point(194) = {0.938153,spn,-0.008643}; -Point(195) = {0.952414,spn,-0.006710}; -Point(196) = {0.964888,spn,-0.004990}; -Point(197) = {0.975528,spn,-0.003501}; -Point(198) = {0.984292,spn,-0.002260}; -Point(199) = {0.991144,spn,-0.001280}; -Point(200) = {0.996057,spn,-0.000572}; -Point(201) = {0.999013,spn,-0.000143,msTe}; - -// Tip: -Point(5101) = {0.999013,spn,0.000000}; -Point(5102) = {0.996057,spn+0.000125,0.000000}; -Point(5103) = {0.991144,spn+0.000331,0.000000}; -Point(5104) = {0.984292,spn+0.000620,0.000000}; -Point(5105) = {0.975528,spn+0.000989,0.000000}; -Point(5106) = {0.964888,spn+0.001437,0.000000}; -Point(5107) = {0.952414,spn+0.001963,0.000000}; -Point(5108) = {0.938153,spn+0.002563,0.000000}; -Point(5109) = {0.922164,spn+0.003237,0.000000}; -Point(5110) = {0.904508,spn+0.003981,0.000000,msTe}; -Point(5111) = {0.885257,spn+0.004791,0.000000}; -Point(5112) = {0.864484,spn+0.005666,0.000000}; -Point(5113) = {0.842274,spn+0.006602,0.000000}; -Point(5114) = {0.818712,spn+0.007594,0.000000}; -Point(5115) = {0.793893,spn+0.008640,0.000000}; -Point(5116) = {0.767913,spn+0.009734,0.000000}; -Point(5117) = {0.740877,spn+0.010873,0.000000}; -Point(5118) = {0.712890,spn+0.012052,0.000000}; -Point(5119) = {0.684062,spn+0.013266,0.000000}; -Point(5120) = {0.654508,spn+0.014511,0.000000}; -Point(5121) = {0.624345,spn+0.015781,0.000000}; -Point(5122) = {0.593691,spn+0.017072,0.000000}; -Point(5123) = {0.562667,spn+0.018379,0.000000}; -Point(5124) = {0.531395,spn+0.019696,0.000000}; -Point(5125) = {0.500000,spn+0.021019,0.000000}; -Point(5126) = {0.468605,spn+0.022341,0.000000}; -Point(5127) = {0.437333,spn+0.023658,0.000000}; -Point(5128) = {0.406309,spn+0.024965,0.000000}; -Point(5129) = {0.375655,spn+0.026256,0.000000}; -Point(5130) = {0.345492,spn+0.027526,0.000000}; -Point(5131) = {0.315938,spn+0.028771,0.000000}; -Point(5132) = {0.287110,spn+0.029985,0.000000,msTe}; -Point(5133) = {0.019492,spn+0.041801,0.000000}; - -// Dummy tip center: -Point(5349) = {0.904508,spn,0.000000}; -Point(5350) = {0.287110,spn,0.000000}; - -// Box -Point(10001) = {1+lgt, 0, 0, msF}; -Point(10002) = {1+lgt, 0, hgt, msF}; -Point(10003) = {-lgt, 0, hgt, msF}; -Point(10004) = {-lgt, 0, -hgt, msF}; -Point(10005) = {1+lgt, 0, -hgt, msF}; -Point(10011) = {1+lgt, wdt, hgt, msF}; -Point(10012) = {-lgt, wdt, hgt, msF}; -Point(10013) = {-lgt, wdt, -hgt, msF}; -Point(10014) = {1+lgt, wdt, -hgt, msF}; - -// Wake -Point(10021) = {1+lgt, spn, 0, msF}; - -/// Lines - -// Airfoil 1: -Spline(1) = {1,2,3,4,5,6,7,8,9,10,11}; -Spline(2) = {11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33}; -Spline(3) = {33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51}; -Spline(4) = {51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69}; -Spline(5) = {69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91}; -Spline(6) = {91,92,93,94,95,96,97,98,99,100,1}; - -// Airfoil 2: -Spline(7) = {102,103,104,105,106,107,108,109,110,111,112}; -Spline(8) = {112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134}; -Spline(9) = {134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152}; -Spline(10) = {152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170}; -Spline(11) = {170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192}; -Spline(12) = {192,193,194,195,196,197,198,199,200,201,102}; - -// Airfoil 1 to airfoil 2: -Line(61) = {1,102}; -Line(62) = {11,112}; -Line(63) = {33,134}; -Line(64) = {51,152}; -Line(65) = {69,170}; -Line(66) = {91,192}; - -// Tip: -Spline(121) = {102,5101,5102,5103,5104,5105,5106,5107,5108,5109,5110}; -Spline(122) = {5110,5111,5112,5113,5114,5115,5116,5117,5118,5119,5120,5121,5122,5123,5124,5125,5126,5127,5128,5129,5130,5131,5132}; -If(GMSH_MAJOR_VERSION >= 4) - Bezier(123) = {5132,5133,152}; -Else - BSpline(123) = {5132,5133,152}; -EndIf -Ellipse(124) = {112,5349,112,5110}; -Ellipse(125) = {134,5350,134,5132}; -Ellipse(126) = {170,5350,170,5132}; -Ellipse(127) = {192,5349,192,5110}; - -// Box -Line(201) = {10001,10002}; -Line(202) = {10002,10003}; -Line(203) = {10003,10004}; -Line(204) = {10004,10005}; -Line(205) = {10005,10001}; -Line(211) = {10011,10012}; -Line(212) = {10012,10013}; -Line(213) = {10013,10014}; -Line(214) = {10014,10011}; -Line(221) = {1, 10001}; -Line(222) = {102, 10021}; -Line(231) = {10002,10011}; -Line(232) = {10003,10012}; -Line(233) = {10004,10013}; -Line(234) = {10005,10014}; - -// Wake -Line(241) = {10001, 10021}; - -/// Line loops & Surfaces - -// Wing 1: -Line Loop(1) = {1,62,-7,-61}; -Line Loop(2) = {2,63,-8,-62}; -Line Loop(3) = {3,64,-9,-63}; -Line Loop(4) = {4,65,-10,-64}; -Line Loop(5) = {5,66,-11,-65}; -Line Loop(6) = {6,61,-12,-66}; -Surface(1) = {-1}; -Surface(2) = {-2}; -Surface(3) = {-3}; -Surface(4) = {-4}; -Surface(5) = {-5}; -Surface(6) = {-6}; - -// Wingtip: -Line Loop(11) = {7,124,-121}; -Line Loop(12) = {8,125,-122,-124}; -Line Loop(13) = {9,-123,-125}; -Line Loop(14) = {10,126,123}; -Line Loop(15) = {11,127,122,-126}; -Line Loop(16) = {12,121,-127}; -Surface(11) = {-11}; -Surface(12) = {-12}; -Surface(13) = {-13}; -Surface(14) = {-14}; -Surface(15) = {-15}; -Surface(16) = {-16}; - -// Symmetry -Line Loop(21) = {201, 202, 203, 204, 205}; -Line Loop(22) = {1, 2, 3, 4, 5, 6}; -Plane Surface(21) = {-21, 22}; - -// Downsteam -Line Loop(31) = {201, 231, -214, -234, 205}; -Plane Surface(31) = {31}; - -// Farfield -Line Loop(41) = {233, -212, -232, 203}; -Plane Surface(41) = {41}; -Line Loop(42) = {204, 234, -213, -233}; -Plane Surface(42) = {42}; -Line Loop(43) = {202, 232, -211, -231}; -Plane Surface(43) = {43}; -Line Loop(44) = {213, 214, 211, 212}; -Plane Surface(44) = {44}; - -// Wake -Line Loop(51) = {221, 241, -222, -61}; -Surface(51) = {51}; - -/// Volume -Surface Loop(1) = {1,2,3,4,5,6,11,12,13,14,15,16,21,31,41,42,43,44}; -Volume(1) = {1}; - -/// Embbeded -Line{221} In Surface{21}; -Line{241} In Surface{31}; -Surface{51} In Volume{1}; - -//// MESHING ALGORITHM - -/// 2D: - -///Wing, farfield and symmetry plane: -MeshAlgorithm Surface{1,2,3,4,5,6,21,31,41,42,43,44,51} = 5; - -///Tip: -MeshAlgorithm Surface{11,12,13,14,15,16} = 1; - -/// 3D: - -Mesh.Algorithm3D = 2; -//Mesh.OptimizeNetgen = 1; -Mesh.Optimize = 1; -Mesh.Smoothing = 10; -Mesh.SmoothNormals = 1; - - - -//// PHYSICAL GROUPS - -// Trailing edge and wake tip -Physical Line("wakeTip") = {222}; -Physical Line("teTip") = {61, 222}; - -// Internal Field: -Physical Volume("field") = {1}; - -// Wing: -Physical Surface("wing") = {1,2,3,11,12,13}; -Physical Surface("wing_") = {4,5,6,14,15,16}; - -// Symmetry: -Physical Surface("symmetry") = {21}; - -// Downstream: -Physical Surface("downstream") = {31}; - -// Farfield: -Physical Surface("upstream") = {41}; -Physical Surface("farfield") = {42,43,44}; - -// Wake: -Physical Surface("wake") = {51}; diff --git a/dart/models/oneraM6.geo b/dart/models/oneraM6.geo index 6c14627beea6c5588ab593509a2b95cf405377c6..25b777bf212557af77618346ac2a17f4690748bd 100644 --- a/dart/models/oneraM6.geo +++ b/dart/models/oneraM6.geo @@ -559,7 +559,7 @@ Mesh.SmoothNormals = 1; // Trailing edge and wake tip Physical Line("wakeTip") = {162}; -Physical Line("teTip") = {61, 162}; +Physical Line("te") = {61}; // Internal Field: Physical Volume("field") = {1}; diff --git a/dart/models/rae_25.geo b/dart/models/rae_25.geo new file mode 100644 index 0000000000000000000000000000000000000000..1deee3cd4c47b1440c640a1be8d462993215b809 --- /dev/null +++ b/dart/models/rae_25.geo @@ -0,0 +1,400 @@ +/* 2.5D extruded RAE 2822 wing */ + +// Parameters +// domain +DefineConstant[ spn = { 1.0, Name "Wing span" }, + lgt = { 6.0, Name "Channel length" }, + hgt = { 6.0, Name "Channel height" }, + msN = { 0.05, Min 0.1, Max 1, Step 0.1, Name "Nearfield mesh size" }, + msF = { 1.0, Min 0.1, Max 1, Step 0.1, Name "Farfield mesh size" } ]; + +// Geometry +// wing +Point(1) = {1.000000, -spn/2, 0.000000, msN}; +Point(2) = {0.999398, -spn/2, 0.000128, msN}; +Point(3) = {0.997592, -spn/2, 0.000510, msN}; +Point(4) = {0.994588, -spn/2, 0.001137, msN}; +Point(5) = {0.990393, -spn/2, 0.002001, msN}; +Point(6) = {0.985016, -spn/2, 0.003092, msN}; +Point(7) = {0.978470, -spn/2, 0.004401, msN}; +Point(8) = {0.970772, -spn/2, 0.005915, msN}; +Point(9) = {0.961940, -spn/2, 0.007622, msN}; +Point(10) = {0.951995, -spn/2, 0.009508, msN}; +Point(11) = {0.940961, -spn/2, 0.011562, msN}; +Point(12) = {0.928864, -spn/2, 0.013769, msN}; +Point(13) = {0.915735, -spn/2, 0.016113, msN}; +Point(14) = {0.901604, -spn/2, 0.018580, msN}; +Point(15) = {0.886505, -spn/2, 0.021153, msN}; +Point(16) = {0.870476, -spn/2, 0.023817, msN}; +Point(17) = {0.853553, -spn/2, 0.026554, msN}; +Point(18) = {0.835779, -spn/2, 0.029347, msN}; +Point(19) = {0.817197, -spn/2, 0.032176, msN}; +Point(20) = {0.797850, -spn/2, 0.035017, msN}; +Point(21) = {0.777785, -spn/2, 0.037847, msN}; +Point(22) = {0.757051, -spn/2, 0.040641, msN}; +Point(23) = {0.735698, -spn/2, 0.043377, msN}; +Point(24) = {0.713778, -spn/2, 0.046029, msN}; +Point(25) = {0.691342, -spn/2, 0.048575, msN}; +Point(26) = {0.668445, -spn/2, 0.050993, msN}; +Point(27) = {0.645142, -spn/2, 0.053258, msN}; +Point(28) = {0.621490, -spn/2, 0.055344, msN}; +Point(29) = {0.597545, -spn/2, 0.057218, msN}; +Point(30) = {0.573365, -spn/2, 0.058845, msN}; +Point(31) = {0.549009, -spn/2, 0.060194, msN}; +Point(32) = {0.524534, -spn/2, 0.061254, msN}; +Point(33) = {0.500000, -spn/2, 0.062029, msN}; +Point(34) = {0.475466, -spn/2, 0.062530, msN}; +Point(35) = {0.450991, -spn/2, 0.062774, msN}; +Point(36) = {0.426635, -spn/2, 0.062779, msN}; +Point(37) = {0.402455, -spn/2, 0.062562, msN}; +Point(38) = {0.378510, -spn/2, 0.062133, msN}; +Point(39) = {0.354858, -spn/2, 0.061497, msN}; +Point(40) = {0.331555, -spn/2, 0.060660, msN}; +Point(41) = {0.308658, -spn/2, 0.059629, msN}; +Point(42) = {0.286222, -spn/2, 0.058414, msN}; +Point(43) = {0.264302, -spn/2, 0.057026, msN}; +Point(44) = {0.242949, -spn/2, 0.055470, msN}; +Point(45) = {0.222215, -spn/2, 0.053753, msN}; +Point(46) = {0.202150, -spn/2, 0.051885, msN}; +Point(47) = {0.182803, -spn/2, 0.049874, msN}; +Point(48) = {0.164221, -spn/2, 0.047729, msN}; +Point(49) = {0.146447, -spn/2, 0.045457, msN}; +Point(50) = {0.129524, -spn/2, 0.043071, msN}; +Point(51) = {0.113495, -spn/2, 0.040585, msN}; +Point(52) = {0.098396, -spn/2, 0.038011, msN}; +Point(53) = {0.084265, -spn/2, 0.035360, msN}; +Point(54) = {0.071136, -spn/2, 0.032644, msN}; +Point(55) = {0.059039, -spn/2, 0.029874, msN}; +Point(56) = {0.048005, -spn/2, 0.027062, msN}; +Point(57) = {0.038060, -spn/2, 0.024219, msN}; +Point(58) = {0.029228, -spn/2, 0.021348, msN}; +Point(59) = {0.021530, -spn/2, 0.018441, msN}; +Point(60) = {0.014984, -spn/2, 0.015489, msN}; +Point(61) = {0.009607, -spn/2, 0.012480, msN}; +Point(62) = {0.005412, -spn/2, 0.009416, msN}; +Point(63) = {0.002408, -spn/2, 0.006306, msN}; +Point(64) = {0.000602, -spn/2, 0.003165, msN}; +Point(65) = {0.000000, -spn/2, 0.000000, msN}; +Point(66) = {0.000602, -spn/2, -0.003160, msN}; +Point(67) = {0.002408, -spn/2, -0.006308, msN}; +Point(68) = {0.005412, -spn/2, -0.009443, msN}; +Point(69) = {0.009607, -spn/2, -0.012559, msN}; +Point(70) = {0.014984, -spn/2, -0.015649, msN}; +Point(71) = {0.021530, -spn/2, -0.018707, msN}; +Point(72) = {0.029228, -spn/2, -0.021722, msN}; +Point(73) = {0.038060, -spn/2, -0.024685, msN}; +Point(74) = {0.048005, -spn/2, -0.027586, msN}; +Point(75) = {0.059039, -spn/2, -0.030416, msN}; +Point(76) = {0.071136, -spn/2, -0.033170, msN}; +Point(77) = {0.084265, -spn/2, -0.035843, msN}; +Point(78) = {0.098396, -spn/2, -0.038431, msN}; +Point(79) = {0.113495, -spn/2, -0.040929, msN}; +Point(80) = {0.129524, -spn/2, -0.043326, msN}; +Point(81) = {0.146447, -spn/2, -0.045610, msN}; +Point(82) = {0.164221, -spn/2, -0.047773, msN}; +Point(83) = {0.182803, -spn/2, -0.049805, msN}; +Point(84) = {0.202150, -spn/2, -0.051694, msN}; +Point(85) = {0.222215, -spn/2, -0.053427, msN}; +Point(86) = {0.242949, -spn/2, -0.054994, msN}; +Point(87) = {0.264302, -spn/2, -0.056376, msN}; +Point(88) = {0.286222, -spn/2, -0.057547, msN}; +Point(89) = {0.308658, -spn/2, -0.058459, msN}; +Point(90) = {0.331555, -spn/2, -0.059046, msN}; +Point(91) = {0.354858, -spn/2, -0.059236, msN}; +Point(92) = {0.378510, -spn/2, -0.058974, msN}; +Point(93) = {0.402455, -spn/2, -0.058224, msN}; +Point(94) = {0.426635, -spn/2, -0.056979, msN}; +Point(95) = {0.450991, -spn/2, -0.055257, msN}; +Point(96) = {0.475466, -spn/2, -0.053099, msN}; +Point(97) = {0.500000, -spn/2, -0.050563, msN}; +Point(98) = {0.524534, -spn/2, -0.047719, msN}; +Point(99) = {0.549009, -spn/2, -0.044642, msN}; +Point(100) = {0.573365, -spn/2, -0.041397, msN}; +Point(101) = {0.597545, -spn/2, -0.038043, msN}; +Point(102) = {0.621490, -spn/2, -0.034631, msN}; +Point(103) = {0.645142, -spn/2, -0.031207, msN}; +Point(104) = {0.668445, -spn/2, -0.027814, msN}; +Point(105) = {0.691342, -spn/2, -0.024495, msN}; +Point(106) = {0.713778, -spn/2, -0.021289, msN}; +Point(107) = {0.735698, -spn/2, -0.018232, msN}; +Point(108) = {0.757051, -spn/2, -0.015357, msN}; +Point(109) = {0.777785, -spn/2, -0.012690, msN}; +Point(110) = {0.797850, -spn/2, -0.010244, msN}; +Point(111) = {0.817197, -spn/2, -0.008027, msN}; +Point(112) = {0.835779, -spn/2, -0.006048, msN}; +Point(113) = {0.853553, -spn/2, -0.004314, msN}; +Point(114) = {0.870476, -spn/2, -0.002829, msN}; +Point(115) = {0.886505, -spn/2, -0.001592, msN}; +Point(116) = {0.901604, -spn/2, -0.000600, msN}; +Point(117) = {0.915735, -spn/2, 0.000157, msN}; +Point(118) = {0.928864, -spn/2, 0.000694, msN}; +Point(119) = {0.940961, -spn/2, 0.001033, msN}; +Point(120) = {0.951995, -spn/2, 0.001197, msN}; +Point(121) = {0.961940, -spn/2, 0.001212, msN}; +Point(122) = {0.970772, -spn/2, 0.001112, msN}; +Point(123) = {0.978470, -spn/2, 0.000935, msN}; +Point(124) = {0.985016, -spn/2, 0.000719, msN}; +Point(125) = {0.990393, -spn/2, 0.000497, msN}; +Point(126) = {0.994588, -spn/2, 0.000296, msN}; +Point(127) = {0.997592, -spn/2, 0.000137, msN}; +Point(128) = {0.999398, -spn/2, 0.000035, msN}; + +Point(301) = {1.000000, spn/2, 0.000000, msN}; +Point(302) = {0.999398, spn/2, 0.000128, msN}; +Point(303) = {0.997592, spn/2, 0.000510, msN}; +Point(304) = {0.994588, spn/2, 0.001137, msN}; +Point(305) = {0.990393, spn/2, 0.002001, msN}; +Point(306) = {0.985016, spn/2, 0.003092, msN}; +Point(307) = {0.978470, spn/2, 0.004401, msN}; +Point(308) = {0.970772, spn/2, 0.005915, msN}; +Point(309) = {0.961940, spn/2, 0.007622, msN}; +Point(310) = {0.951995, spn/2, 0.009508, msN}; +Point(311) = {0.940961, spn/2, 0.011562, msN}; +Point(312) = {0.928864, spn/2, 0.013769, msN}; +Point(313) = {0.915735, spn/2, 0.016113, msN}; +Point(314) = {0.901604, spn/2, 0.018580, msN}; +Point(315) = {0.886505, spn/2, 0.021153, msN}; +Point(316) = {0.870476, spn/2, 0.023817, msN}; +Point(317) = {0.853553, spn/2, 0.026554, msN}; +Point(318) = {0.835779, spn/2, 0.029347, msN}; +Point(319) = {0.817197, spn/2, 0.032176, msN}; +Point(320) = {0.797850, spn/2, 0.035017, msN}; +Point(321) = {0.777785, spn/2, 0.037847, msN}; +Point(322) = {0.757051, spn/2, 0.040641, msN}; +Point(323) = {0.735698, spn/2, 0.043377, msN}; +Point(324) = {0.713778, spn/2, 0.046029, msN}; +Point(325) = {0.691342, spn/2, 0.048575, msN}; +Point(326) = {0.668445, spn/2, 0.050993, msN}; +Point(327) = {0.645142, spn/2, 0.053258, msN}; +Point(328) = {0.621490, spn/2, 0.055344, msN}; +Point(329) = {0.597545, spn/2, 0.057218, msN}; +Point(330) = {0.573365, spn/2, 0.058845, msN}; +Point(331) = {0.549009, spn/2, 0.060194, msN}; +Point(332) = {0.524534, spn/2, 0.061254, msN}; +Point(333) = {0.500000, spn/2, 0.062029, msN}; +Point(334) = {0.475466, spn/2, 0.062530, msN}; +Point(335) = {0.450991, spn/2, 0.062774, msN}; +Point(336) = {0.426635, spn/2, 0.062779, msN}; +Point(337) = {0.402455, spn/2, 0.062562, msN}; +Point(338) = {0.378510, spn/2, 0.062133, msN}; +Point(339) = {0.354858, spn/2, 0.061497, msN}; +Point(340) = {0.331555, spn/2, 0.060660, msN}; +Point(341) = {0.308658, spn/2, 0.059629, msN}; +Point(342) = {0.286222, spn/2, 0.058414, msN}; +Point(343) = {0.264302, spn/2, 0.057026, msN}; +Point(344) = {0.242949, spn/2, 0.055470, msN}; +Point(345) = {0.222215, spn/2, 0.053753, msN}; +Point(346) = {0.202150, spn/2, 0.051885, msN}; +Point(347) = {0.182803, spn/2, 0.049874, msN}; +Point(348) = {0.164221, spn/2, 0.047729, msN}; +Point(349) = {0.146447, spn/2, 0.045457, msN}; +Point(350) = {0.129524, spn/2, 0.043071, msN}; +Point(351) = {0.113495, spn/2, 0.040585, msN}; +Point(352) = {0.098396, spn/2, 0.038011, msN}; +Point(353) = {0.084265, spn/2, 0.035360, msN}; +Point(354) = {0.071136, spn/2, 0.032644, msN}; +Point(355) = {0.059039, spn/2, 0.029874, msN}; +Point(356) = {0.048005, spn/2, 0.027062, msN}; +Point(357) = {0.038060, spn/2, 0.024219, msN}; +Point(358) = {0.029228, spn/2, 0.021348, msN}; +Point(359) = {0.021530, spn/2, 0.018441, msN}; +Point(360) = {0.014984, spn/2, 0.015489, msN}; +Point(361) = {0.009607, spn/2, 0.012480, msN}; +Point(362) = {0.005412, spn/2, 0.009416, msN}; +Point(363) = {0.002408, spn/2, 0.006306, msN}; +Point(364) = {0.000602, spn/2, 0.003165, msN}; +Point(365) = {0.000000, spn/2, 0.000000, msN}; +Point(366) = {0.000602, spn/2, -0.003160, msN}; +Point(367) = {0.002408, spn/2, -0.006308, msN}; +Point(368) = {0.005412, spn/2, -0.009443, msN}; +Point(369) = {0.009607, spn/2, -0.012559, msN}; +Point(370) = {0.014984, spn/2, -0.015649, msN}; +Point(371) = {0.021530, spn/2, -0.018707, msN}; +Point(372) = {0.029228, spn/2, -0.021722, msN}; +Point(373) = {0.038060, spn/2, -0.024685, msN}; +Point(374) = {0.048005, spn/2, -0.027586, msN}; +Point(375) = {0.059039, spn/2, -0.030416, msN}; +Point(376) = {0.071136, spn/2, -0.033170, msN}; +Point(377) = {0.084265, spn/2, -0.035843, msN}; +Point(378) = {0.098396, spn/2, -0.038431, msN}; +Point(379) = {0.113495, spn/2, -0.040929, msN}; +Point(380) = {0.129524, spn/2, -0.043326, msN}; +Point(381) = {0.146447, spn/2, -0.045610, msN}; +Point(382) = {0.164221, spn/2, -0.047773, msN}; +Point(383) = {0.182803, spn/2, -0.049805, msN}; +Point(384) = {0.202150, spn/2, -0.051694, msN}; +Point(385) = {0.222215, spn/2, -0.053427, msN}; +Point(386) = {0.242949, spn/2, -0.054994, msN}; +Point(387) = {0.264302, spn/2, -0.056376, msN}; +Point(388) = {0.286222, spn/2, -0.057547, msN}; +Point(389) = {0.308658, spn/2, -0.058459, msN}; +Point(390) = {0.331555, spn/2, -0.059046, msN}; +Point(391) = {0.354858, spn/2, -0.059236, msN}; +Point(392) = {0.378510, spn/2, -0.058974, msN}; +Point(393) = {0.402455, spn/2, -0.058224, msN}; +Point(394) = {0.426635, spn/2, -0.056979, msN}; +Point(395) = {0.450991, spn/2, -0.055257, msN}; +Point(396) = {0.475466, spn/2, -0.053099, msN}; +Point(397) = {0.500000, spn/2, -0.050563, msN}; +Point(398) = {0.524534, spn/2, -0.047719, msN}; +Point(399) = {0.549009, spn/2, -0.044642, msN}; +Point(400) = {0.573365, spn/2, -0.041397, msN}; +Point(401) = {0.597545, spn/2, -0.038043, msN}; +Point(402) = {0.621490, spn/2, -0.034631, msN}; +Point(403) = {0.645142, spn/2, -0.031207, msN}; +Point(404) = {0.668445, spn/2, -0.027814, msN}; +Point(405) = {0.691342, spn/2, -0.024495, msN}; +Point(406) = {0.713778, spn/2, -0.021289, msN}; +Point(407) = {0.735698, spn/2, -0.018232, msN}; +Point(408) = {0.757051, spn/2, -0.015357, msN}; +Point(409) = {0.777785, spn/2, -0.012690, msN}; +Point(410) = {0.797850, spn/2, -0.010244, msN}; +Point(411) = {0.817197, spn/2, -0.008027, msN}; +Point(412) = {0.835779, spn/2, -0.006048, msN}; +Point(413) = {0.853553, spn/2, -0.004314, msN}; +Point(414) = {0.870476, spn/2, -0.002829, msN}; +Point(415) = {0.886505, spn/2, -0.001592, msN}; +Point(416) = {0.901604, spn/2, -0.000600, msN}; +Point(417) = {0.915735, spn/2, 0.000157, msN}; +Point(418) = {0.928864, spn/2, 0.000694, msN}; +Point(419) = {0.940961, spn/2, 0.001033, msN}; +Point(420) = {0.951995, spn/2, 0.001197, msN}; +Point(421) = {0.961940, spn/2, 0.001212, msN}; +Point(422) = {0.970772, spn/2, 0.001112, msN}; +Point(423) = {0.978470, spn/2, 0.000935, msN}; +Point(424) = {0.985016, spn/2, 0.000719, msN}; +Point(425) = {0.990393, spn/2, 0.000497, msN}; +Point(426) = {0.994588, spn/2, 0.000296, msN}; +Point(427) = {0.997592, spn/2, 0.000137, msN}; +Point(428) = {0.999398, spn/2, 0.000035, msN}; + +// box +Point(631) = {1+lgt/2, -spn/2, hgt/2, msF}; +Point(632) = {-lgt/2, -spn/2, hgt/2, msF}; +Point(633) = {-lgt/2, -spn/2, -hgt/2, msF}; +Point(634) = {1+lgt/2, -spn/2, -hgt/2, msF}; +Point(635) = {1+lgt/2, spn/2, hgt/2, msF}; +Point(636) = {-lgt/2, spn/2, hgt/2, msF}; +Point(637) = {-lgt/2, spn/2, -hgt/2, msF}; +Point(638) = {1+lgt/2, spn/2, -hgt/2, msF}; + +// midplane +Point(641) = {1+lgt/2, -spn/2, 0, msF}; +Point(642) = {1+lgt/2, spn/2, 0, msF}; +Point(643) = {-lgt/2, spn/2, 0, msF}; +Point(644) = {-lgt/2, -spn/2, 0, msF}; + +// wing +Spline(1) = {1:44}; +Spline(2) = {44:65}; +Spline(3) = {65:86}; +Spline(4) = {86:128,1}; +Spline(5) = {301:344}; +Spline(6) = {344:365}; +Spline(7) = {365:386}; +Spline(8) = {386:428,301}; +Line(9) = {1, 301}; +Line(10) = {44, 344}; +Line(11) = {65, 365}; +Line(12) = {86, 386}; + +// box +Line(21) = {641, 631}; +Line(22) = {631, 632}; +Line(23) = {632, 644}; +Line(24) = {644, 633}; +Line(25) = {633, 634}; +Line(26) = {634, 641}; +Line(27) = {642, 635}; +Line(28) = {635, 636}; +Line(29) = {636, 643}; +Line(30) = {643, 637}; +Line(31) = {637, 638}; +Line(32) = {638, 642}; +Line(33) = {631, 635}; +Line(34) = {632, 636}; +Line(35) = {633, 637}; +Line(36) = {634, 638}; + +// midplane +Line(41) = {1, 641}; +Line(42) = {301, 642}; +Line(43) = {365, 643}; +Line(44) = {65, 644}; +Line(45) = {641, 642}; +Line(46) = {643, 644}; + +// wing +Line Loop(1) = {1, 10, -5, -9}; +Surface(1) = {-1}; +Line Loop(2) = {2, 11, -6, -10}; +Surface(2) = {-2}; +Line Loop(3) = {3, 12, -7, -11}; +Surface(3) = {-3}; +Line Loop(4) = {4, 9, -8, -12}; +Surface(4) = {-4}; + +// symmetry +Line Loop(5) = {41, 21, 22, 23, -44, -2, -1}; +Plane Surface(5) = {-5}; +Line Loop(12) = {42, 27, 28, 29, -43, -5, -6}; +Plane Surface(12) = {12}; +Line Loop(6) = {44, 24, 25, 26, -41, -4, -3}; +Plane Surface(6) = {-6}; +Line Loop(13) = {43, 30, 31, 32, -42, -8, -7}; +Plane Surface(13) = {13}; + +// upstream +Line Loop(7) = {46, -23, 34, 29}; +Plane Surface(7) = {-7}; +Line Loop(8) = {46, 24, 35, -30}; +Plane Surface(8) = {8}; + +// downstream +Line Loop(9) = {45, 27, -33, -21}; +Plane Surface(9) = {-9}; +Line Loop(10) = {45, -32, -36, 26}; +Plane Surface(10) = {10}; + +// farfield +Line Loop(11) = {22, 34, -28, -33}; +Plane Surface(11) = {11}; +Line Loop(14) = {25, 36, -31, -35}; +Plane Surface(14) = {14}; + +// wake +Line Loop(15) = {41, 45, -42, -9}; +Plane Surface(15) = {15}; + +// internal +Line Loop(16) = {11, 43, 46, -44}; +Plane Surface(16) = {16}; + +// field +Surface Loop(1) = {7, 5, 9, 12, 11, 1, 2, 16, 15}; +Volume(1) = {1}; +Surface Loop(2) = {6, 8, 14, 10, 13, 3, 4, 16, 15}; +Volume(2) = {2}; + +// Mesh +// 2D +MeshAlgorithm Surface{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16} = 5; +// 3D +Mesh.Algorithm3D = 2; +Mesh.OptimizeNetgen = 1; +Mesh.Smoothing = 10; +Mesh.SmoothNormals = 1; + +// Physical +Physical Line("te") = {9}; +Physical Surface("wing") = {1, 2}; +Physical Surface("wing_") = {3, 4}; +Physical Surface("symmetry") = {5, 12}; +Physical Surface("symmetry_") = {6, 13}; +Physical Surface("upstream") = {7, 8}; +Physical Surface("downstream") = {9}; +Physical Surface("downstream_") = {10}; +Physical Surface("farfield") = {11, 14}; +Physical Surface("wake") = {15}; +Physical Volume("field") = {1}; +Physical Volume("field_") = {2}; diff --git a/dart/models/rae_3.geo b/dart/models/rae_3.geo new file mode 100644 index 0000000000000000000000000000000000000000..2c76facfc85fe7c42ea4173652badb3f87d5b6d7 --- /dev/null +++ b/dart/models/rae_3.geo @@ -0,0 +1,430 @@ +/* RAE 2822 rectangular wing */ + +// Parameters +// domain and mesh +DefineConstant[ spn = { 1.0, Name "Wing span" }, + lgt = { 6.0, Name "Channel length" }, + wdt = { 3.0, Name "Channel width" }, + hgt = { 6.0, Name "Channel height" }, + msN = { 0.05, Name "Leading edge mesh size" }, + msF = { 1.0, Name "Farfield mesh size" } ]; + +//// GEOMETRY + + +/// Points + +// Airfoil 1 +Point(1) = {1.000000, 0.0000, 0.000000, msN}; +Point(2) = {0.999398, 0.0000, 0.000128, msN}; +Point(3) = {0.997592, 0.0000, 0.000510, msN}; +Point(4) = {0.994588, 0.0000, 0.001137, msN}; +Point(5) = {0.990393, 0.0000, 0.002001, msN}; +Point(6) = {0.985016, 0.0000, 0.003092, msN}; +Point(7) = {0.978470, 0.0000, 0.004401, msN}; +Point(8) = {0.970772, 0.0000, 0.005915, msN}; +Point(9) = {0.961940, 0.0000, 0.007622, msN}; +Point(10) = {0.951995, 0.0000, 0.009508, msN}; +Point(11) = {0.940961, 0.0000, 0.011562, msN}; +Point(12) = {0.928864, 0.0000, 0.013769, msN}; +Point(13) = {0.915735, 0.0000, 0.016113, msN}; +Point(14) = {0.901604, 0.0000, 0.018580, msN}; +Point(15) = {0.886505, 0.0000, 0.021153, msN}; +Point(16) = {0.870476, 0.0000, 0.023817, msN}; +Point(17) = {0.853553, 0.0000, 0.026554, msN}; +Point(18) = {0.835779, 0.0000, 0.029347, msN}; +Point(19) = {0.817197, 0.0000, 0.032176, msN}; +Point(20) = {0.797850, 0.0000, 0.035017, msN}; +Point(21) = {0.777785, 0.0000, 0.037847, msN}; +Point(22) = {0.757051, 0.0000, 0.040641, msN}; +Point(23) = {0.735698, 0.0000, 0.043377, msN}; +Point(24) = {0.713778, 0.0000, 0.046029, msN}; +Point(25) = {0.691342, 0.0000, 0.048575, msN}; +Point(26) = {0.668445, 0.0000, 0.050993, msN}; +Point(27) = {0.645142, 0.0000, 0.053258, msN}; +Point(28) = {0.621490, 0.0000, 0.055344, msN}; +Point(29) = {0.597545, 0.0000, 0.057218, msN}; +Point(30) = {0.573365, 0.0000, 0.058845, msN}; +Point(31) = {0.549009, 0.0000, 0.060194, msN}; +Point(32) = {0.524534, 0.0000, 0.061254, msN}; +Point(33) = {0.500000, 0.0000, 0.062029, msN}; +Point(34) = {0.475466, 0.0000, 0.062530, msN}; +Point(35) = {0.450991, 0.0000, 0.062774, msN}; +Point(36) = {0.426635, 0.0000, 0.062779, msN}; +Point(37) = {0.402455, 0.0000, 0.062562, msN}; +Point(38) = {0.378510, 0.0000, 0.062133, msN}; +Point(39) = {0.354858, 0.0000, 0.061497, msN}; +Point(40) = {0.331555, 0.0000, 0.060660, msN}; +Point(41) = {0.308658, 0.0000, 0.059629, msN}; +Point(42) = {0.286222, 0.0000, 0.058414, msN}; +Point(43) = {0.264302, 0.0000, 0.057026, msN}; +Point(44) = {0.242949, 0.0000, 0.055470, msN}; +Point(45) = {0.222215, 0.0000, 0.053753, msN}; +Point(46) = {0.202150, 0.0000, 0.051885, msN}; +Point(47) = {0.182803, 0.0000, 0.049874, msN}; +Point(48) = {0.164221, 0.0000, 0.047729, msN}; +Point(49) = {0.146447, 0.0000, 0.045457, msN}; +Point(50) = {0.129524, 0.0000, 0.043071, msN}; +Point(51) = {0.113495, 0.0000, 0.040585, msN}; +Point(52) = {0.098396, 0.0000, 0.038011, msN}; +Point(53) = {0.084265, 0.0000, 0.035360, msN}; +Point(54) = {0.071136, 0.0000, 0.032644, msN}; +Point(55) = {0.059039, 0.0000, 0.029874, msN}; +Point(56) = {0.048005, 0.0000, 0.027062, msN}; +Point(57) = {0.038060, 0.0000, 0.024219, msN}; +Point(58) = {0.029228, 0.0000, 0.021348, msN}; +Point(59) = {0.021530, 0.0000, 0.018441, msN}; +Point(60) = {0.014984, 0.0000, 0.015489, msN}; +Point(61) = {0.009607, 0.0000, 0.012480, msN}; +Point(62) = {0.005412, 0.0000, 0.009416, msN}; +Point(63) = {0.002408, 0.0000, 0.006306, msN}; +Point(64) = {0.000602, 0.0000, 0.003165, msN}; +Point(65) = {0.000000, 0.0000, 0.000000, msN}; +Point(66) = {0.000602, 0.0000, -0.003160, msN}; +Point(67) = {0.002408, 0.0000, -0.006308, msN}; +Point(68) = {0.005412, 0.0000, -0.009443, msN}; +Point(69) = {0.009607, 0.0000, -0.012559, msN}; +Point(70) = {0.014984, 0.0000, -0.015649, msN}; +Point(71) = {0.021530, 0.0000, -0.018707, msN}; +Point(72) = {0.029228, 0.0000, -0.021722, msN}; +Point(73) = {0.038060, 0.0000, -0.024685, msN}; +Point(74) = {0.048005, 0.0000, -0.027586, msN}; +Point(75) = {0.059039, 0.0000, -0.030416, msN}; +Point(76) = {0.071136, 0.0000, -0.033170, msN}; +Point(77) = {0.084265, 0.0000, -0.035843, msN}; +Point(78) = {0.098396, 0.0000, -0.038431, msN}; +Point(79) = {0.113495, 0.0000, -0.040929, msN}; +Point(80) = {0.129524, 0.0000, -0.043326, msN}; +Point(81) = {0.146447, 0.0000, -0.045610, msN}; +Point(82) = {0.164221, 0.0000, -0.047773, msN}; +Point(83) = {0.182803, 0.0000, -0.049805, msN}; +Point(84) = {0.202150, 0.0000, -0.051694, msN}; +Point(85) = {0.222215, 0.0000, -0.053427, msN}; +Point(86) = {0.242949, 0.0000, -0.054994, msN}; +Point(87) = {0.264302, 0.0000, -0.056376, msN}; +Point(88) = {0.286222, 0.0000, -0.057547, msN}; +Point(89) = {0.308658, 0.0000, -0.058459, msN}; +Point(90) = {0.331555, 0.0000, -0.059046, msN}; +Point(91) = {0.354858, 0.0000, -0.059236, msN}; +Point(92) = {0.378510, 0.0000, -0.058974, msN}; +Point(93) = {0.402455, 0.0000, -0.058224, msN}; +Point(94) = {0.426635, 0.0000, -0.056979, msN}; +Point(95) = {0.450991, 0.0000, -0.055257, msN}; +Point(96) = {0.475466, 0.0000, -0.053099, msN}; +Point(97) = {0.500000, 0.0000, -0.050563, msN}; +Point(98) = {0.524534, 0.0000, -0.047719, msN}; +Point(99) = {0.549009, 0.0000, -0.044642, msN}; +Point(100) = {0.573365, 0.0000, -0.041397, msN}; +Point(101) = {0.597545, 0.0000, -0.038043, msN}; +Point(102) = {0.621490, 0.0000, -0.034631, msN}; +Point(103) = {0.645142, 0.0000, -0.031207, msN}; +Point(104) = {0.668445, 0.0000, -0.027814, msN}; +Point(105) = {0.691342, 0.0000, -0.024495, msN}; +Point(106) = {0.713778, 0.0000, -0.021289, msN}; +Point(107) = {0.735698, 0.0000, -0.018232, msN}; +Point(108) = {0.757051, 0.0000, -0.015357, msN}; +Point(109) = {0.777785, 0.0000, -0.012690, msN}; +Point(110) = {0.797850, 0.0000, -0.010244, msN}; +Point(111) = {0.817197, 0.0000, -0.008027, msN}; +Point(112) = {0.835779, 0.0000, -0.006048, msN}; +Point(113) = {0.853553, 0.0000, -0.004314, msN}; +Point(114) = {0.870476, 0.0000, -0.002829, msN}; +Point(115) = {0.886505, 0.0000, -0.001592, msN}; +Point(116) = {0.901604, 0.0000, -0.000600, msN}; +Point(117) = {0.915735, 0.0000, 0.000157, msN}; +Point(118) = {0.928864, 0.0000, 0.000694, msN}; +Point(119) = {0.940961, 0.0000, 0.001033, msN}; +Point(120) = {0.951995, 0.0000, 0.001197, msN}; +Point(121) = {0.961940, 0.0000, 0.001212, msN}; +Point(122) = {0.970772, 0.0000, 0.001112, msN}; +Point(123) = {0.978470, 0.0000, 0.000935, msN}; +Point(124) = {0.985016, 0.0000, 0.000719, msN}; +Point(125) = {0.990393, 0.0000, 0.000497, msN}; +Point(126) = {0.994588, 0.0000, 0.000296, msN}; +Point(127) = {0.997592, 0.0000, 0.000137, msN}; +Point(128) = {0.999398, 0.0000, 0.000035, msN}; + +// Airfoil 2 +Point(301) = {1.000000, spn, 0.000000, msN}; +Point(302) = {0.999398, spn, 0.000128, msN}; +Point(303) = {0.997592, spn, 0.000510, msN}; +Point(304) = {0.994588, spn, 0.001137, msN}; +Point(305) = {0.990393, spn, 0.002001, msN}; +Point(306) = {0.985016, spn, 0.003092, msN}; +Point(307) = {0.978470, spn, 0.004401, msN}; +Point(308) = {0.970772, spn, 0.005915, msN}; +Point(309) = {0.961940, spn, 0.007622, msN}; +Point(310) = {0.951995, spn, 0.009508, msN}; +Point(311) = {0.940961, spn, 0.011562, msN}; +Point(312) = {0.928864, spn, 0.013769, msN}; +Point(313) = {0.915735, spn, 0.016113, msN}; +Point(314) = {0.901604, spn, 0.018580, msN}; +Point(315) = {0.886505, spn, 0.021153, msN}; +Point(316) = {0.870476, spn, 0.023817, msN}; +Point(317) = {0.853553, spn, 0.026554, msN}; +Point(318) = {0.835779, spn, 0.029347, msN}; +Point(319) = {0.817197, spn, 0.032176, msN}; +Point(320) = {0.797850, spn, 0.035017, msN}; +Point(321) = {0.777785, spn, 0.037847, msN}; +Point(322) = {0.757051, spn, 0.040641, msN}; +Point(323) = {0.735698, spn, 0.043377, msN}; +Point(324) = {0.713778, spn, 0.046029, msN}; +Point(325) = {0.691342, spn, 0.048575, msN}; +Point(326) = {0.668445, spn, 0.050993, msN}; +Point(327) = {0.645142, spn, 0.053258, msN}; +Point(328) = {0.621490, spn, 0.055344, msN}; +Point(329) = {0.597545, spn, 0.057218, msN}; +Point(330) = {0.573365, spn, 0.058845, msN}; +Point(331) = {0.549009, spn, 0.060194, msN}; +Point(332) = {0.524534, spn, 0.061254, msN}; +Point(333) = {0.500000, spn, 0.062029, msN}; +Point(334) = {0.475466, spn, 0.062530, msN}; +Point(335) = {0.450991, spn, 0.062774, msN}; +Point(336) = {0.426635, spn, 0.062779, msN}; +Point(337) = {0.402455, spn, 0.062562, msN}; +Point(338) = {0.378510, spn, 0.062133, msN}; +Point(339) = {0.354858, spn, 0.061497, msN}; +Point(340) = {0.331555, spn, 0.060660, msN}; +Point(341) = {0.308658, spn, 0.059629, msN}; +Point(342) = {0.286222, spn, 0.058414, msN}; +Point(343) = {0.264302, spn, 0.057026, msN}; +Point(344) = {0.242949, spn, 0.055470, msN}; +Point(345) = {0.222215, spn, 0.053753, msN}; +Point(346) = {0.202150, spn, 0.051885, msN}; +Point(347) = {0.182803, spn, 0.049874, msN}; +Point(348) = {0.164221, spn, 0.047729, msN}; +Point(349) = {0.146447, spn, 0.045457, msN}; +Point(350) = {0.129524, spn, 0.043071, msN}; +Point(351) = {0.113495, spn, 0.040585, msN}; +Point(352) = {0.098396, spn, 0.038011, msN}; +Point(353) = {0.084265, spn, 0.035360, msN}; +Point(354) = {0.071136, spn, 0.032644, msN}; +Point(355) = {0.059039, spn, 0.029874, msN}; +Point(356) = {0.048005, spn, 0.027062, msN}; +Point(357) = {0.038060, spn, 0.024219, msN}; +Point(358) = {0.029228, spn, 0.021348, msN}; +Point(359) = {0.021530, spn, 0.018441, msN}; +Point(360) = {0.014984, spn, 0.015489, msN}; +Point(361) = {0.009607, spn, 0.012480, msN}; +Point(362) = {0.005412, spn, 0.009416, msN}; +Point(363) = {0.002408, spn, 0.006306, msN}; +Point(364) = {0.000602, spn, 0.003165, msN}; +Point(365) = {0.000000, spn, 0.000000, msN}; +Point(366) = {0.000602, spn, -0.003160, msN}; +Point(367) = {0.002408, spn, -0.006308, msN}; +Point(368) = {0.005412, spn, -0.009443, msN}; +Point(369) = {0.009607, spn, -0.012559, msN}; +Point(370) = {0.014984, spn, -0.015649, msN}; +Point(371) = {0.021530, spn, -0.018707, msN}; +Point(372) = {0.029228, spn, -0.021722, msN}; +Point(373) = {0.038060, spn, -0.024685, msN}; +Point(374) = {0.048005, spn, -0.027586, msN}; +Point(375) = {0.059039, spn, -0.030416, msN}; +Point(376) = {0.071136, spn, -0.033170, msN}; +Point(377) = {0.084265, spn, -0.035843, msN}; +Point(378) = {0.098396, spn, -0.038431, msN}; +Point(379) = {0.113495, spn, -0.040929, msN}; +Point(380) = {0.129524, spn, -0.043326, msN}; +Point(381) = {0.146447, spn, -0.045610, msN}; +Point(382) = {0.164221, spn, -0.047773, msN}; +Point(383) = {0.182803, spn, -0.049805, msN}; +Point(384) = {0.202150, spn, -0.051694, msN}; +Point(385) = {0.222215, spn, -0.053427, msN}; +Point(386) = {0.242949, spn, -0.054994, msN}; +Point(387) = {0.264302, spn, -0.056376, msN}; +Point(388) = {0.286222, spn, -0.057547, msN}; +Point(389) = {0.308658, spn, -0.058459, msN}; +Point(390) = {0.331555, spn, -0.059046, msN}; +Point(391) = {0.354858, spn, -0.059236, msN}; +Point(392) = {0.378510, spn, -0.058974, msN}; +Point(393) = {0.402455, spn, -0.058224, msN}; +Point(394) = {0.426635, spn, -0.056979, msN}; +Point(395) = {0.450991, spn, -0.055257, msN}; +Point(396) = {0.475466, spn, -0.053099, msN}; +Point(397) = {0.500000, spn, -0.050563, msN}; +Point(398) = {0.524534, spn, -0.047719, msN}; +Point(399) = {0.549009, spn, -0.044642, msN}; +Point(400) = {0.573365, spn, -0.041397, msN}; +Point(401) = {0.597545, spn, -0.038043, msN}; +Point(402) = {0.621490, spn, -0.034631, msN}; +Point(403) = {0.645142, spn, -0.031207, msN}; +Point(404) = {0.668445, spn, -0.027814, msN}; +Point(405) = {0.691342, spn, -0.024495, msN}; +Point(406) = {0.713778, spn, -0.021289, msN}; +Point(407) = {0.735698, spn, -0.018232, msN}; +Point(408) = {0.757051, spn, -0.015357, msN}; +Point(409) = {0.777785, spn, -0.012690, msN}; +Point(410) = {0.797850, spn, -0.010244, msN}; +Point(411) = {0.817197, spn, -0.008027, msN}; +Point(412) = {0.835779, spn, -0.006048, msN}; +Point(413) = {0.853553, spn, -0.004314, msN}; +Point(414) = {0.870476, spn, -0.002829, msN}; +Point(415) = {0.886505, spn, -0.001592, msN}; +Point(416) = {0.901604, spn, -0.000600, msN}; +Point(417) = {0.915735, spn, 0.000157, msN}; +Point(418) = {0.928864, spn, 0.000694, msN}; +Point(419) = {0.940961, spn, 0.001033, msN}; +Point(420) = {0.951995, spn, 0.001197, msN}; +Point(421) = {0.961940, spn, 0.001212, msN}; +Point(422) = {0.970772, spn, 0.001112, msN}; +Point(423) = {0.978470, spn, 0.000935, msN}; +Point(424) = {0.985016, spn, 0.000719, msN}; +Point(425) = {0.990393, spn, 0.000497, msN}; +Point(426) = {0.994588, spn, 0.000296, msN}; +Point(427) = {0.997592, spn, 0.000137, msN}; +Point(428) = {0.999398, spn, 0.000035, msN}; + +// Box +Point(10001) = {1+lgt/2, 0, 0, msF}; +Point(10002) = {1+lgt/2, 0, hgt/2, msF}; +Point(10003) = {-lgt/2, 0, hgt/2, msF}; +Point(10004) = {-lgt/2, 0, -hgt/2, msF}; +Point(10005) = {1+lgt/2, 0, -hgt/2, msF}; +Point(10011) = {1+lgt/2, wdt, hgt/2, msF}; +Point(10012) = {-lgt/2, wdt, hgt/2, msF}; +Point(10013) = {-lgt/2, wdt, -hgt/2, msF}; +Point(10014) = {1+lgt/2, wdt, -hgt/2, msF}; + +// Wake +Point(10021) = {1+lgt/2, spn, 0, msF}; + +/// Lines + +// Airfoil 1: +Spline(1) = {1:14}; +Spline(2) = {14:46}; +Spline(3) = {46:65}; +Spline(4) = {65:84}; +Spline(5) = {84:116}; +Spline(6) = {116:128,1}; + +// Airfoil 2: +Spline(7) = {301:314}; +Spline(8) = {314:346}; +Spline(9) = {346:365}; +Spline(10) = {365:384}; +Spline(11) = {384:416}; +Spline(12) = {416:428,301}; + +// Airfoil 1 to airfoil 2: +Line(61) = {1,301}; +Line(62) = {14,314}; +Line(63) = {46,346}; +Line(64) = {65,365}; +Line(65) = {84,384}; +Line(66) = {116,416}; + +// Box +Line(201) = {10001,10002}; +Line(202) = {10002,10003}; +Line(203) = {10003,10004}; +Line(204) = {10004,10005}; +Line(205) = {10005,10001}; +Line(211) = {10011,10012}; +Line(212) = {10012,10013}; +Line(213) = {10013,10014}; +Line(214) = {10014,10011}; +Line(221) = {1, 10001}; +Line(222) = {301, 10021}; +Line(231) = {10002,10011}; +Line(232) = {10003,10012}; +Line(233) = {10004,10013}; +Line(234) = {10005,10014}; + +// Wake +Line(241) = {10001, 10021}; + +/// Line loops & Surfaces + +// Wing 1: +Line Loop(1) = {1,62,-7,-61}; +Line Loop(2) = {2,63,-8,-62}; +Line Loop(3) = {3,64,-9,-63}; +Line Loop(4) = {4,65,-10,-64}; +Line Loop(5) = {5,66,-11,-65}; +Line Loop(6) = {6,61,-12,-66}; +Surface(1) = {-1}; +Surface(2) = {-2}; +Surface(3) = {-3}; +Surface(4) = {-4}; +Surface(5) = {-5}; +Surface(6) = {-6}; + +// Wingtip: +Line Loop(11) = {7:12}; +Plane Surface(11) = {-11}; + +// Symmetry +Line Loop(21) = {201, 202, 203, 204, 205}; +Line Loop(22) = {1, 2, 3, 4, 5, 6}; +Plane Surface(21) = {-21, 22}; + +// Downsteam +Line Loop(31) = {201, 231, -214, -234, 205}; +Plane Surface(31) = {31}; + +// Farfield +Line Loop(41) = {233, -212, -232, 203}; +Plane Surface(41) = {41}; +Line Loop(42) = {204, 234, -213, -233}; +Plane Surface(42) = {42}; +Line Loop(43) = {202, 232, -211, -231}; +Plane Surface(43) = {43}; +Line Loop(44) = {213, 214, 211, 212}; +Plane Surface(44) = {44}; + +// Wake +Line Loop(51) = {221, 241, -222, -61}; +Surface(51) = {51}; + +/// Volume +Surface Loop(1) = {1,2,3,4,5,6,11,21,31,41,42,43,44}; +Volume(1) = {1}; + +/// Embbeded +Line{221} In Surface{21}; +Line{241} In Surface{31}; +Surface{51} In Volume{1}; + +//// MESHING ALGORITHM + +/// 2D: + +///Wing, farfield and symmetry plane: +MeshAlgorithm Surface{1,2,3,4,5,6,11,21,31,41,42,43,44,51} = 5; + +/// 3D: + +Mesh.Algorithm3D = 2; +//Mesh.OptimizeNetgen = 1; +Mesh.Optimize = 1; +Mesh.Smoothing = 10; +Mesh.SmoothNormals = 1; + + + +//// PHYSICAL GROUPS + +// Trailing edge and wake tip +Physical Line("wakeTip") = {222}; +Physical Line("te") = {61}; + +// Internal Field: +Physical Volume("field") = {1}; + +// Wing: +Physical Surface("wing") = {1,2,3,11}; +Physical Surface("wing_") = {4,5,6}; + +// Symmetry: +Physical Surface("symmetry") = {21}; + +// Downstream: +Physical Surface("downstream") = {31}; + +// Farfield: +Physical Surface("upstream") = {41}; +Physical Surface("farfield") = {42,43,44}; + +// Wake: +Physical Surface("wake") = {51}; diff --git a/dart/models/wbht.geo b/dart/models/wbht.geo index 0e79a8e9f93db547430c2e713af28971313e8708..9f252d80dc453d9c0f44d661fd8b2facf94c3dd5 100644 --- a/dart/models/wbht.geo +++ b/dart/models/wbht.geo @@ -950,9 +950,9 @@ Line(3026) = {3049, 3149}; /// Wake // Wing -Point(4001) = {25, 0, -2.26}; -Point(4002) = {25, 6.1, -0.66}; -Point(4003) = {25, 7, -0.41}; +Point(4001) = {25, 0, -0.534}; +Point(4002) = {25, 6.1, 0.01642}; +Point(4003) = {25, 7, 0.2577}; Line(4001) = {2001, 4001}; Line(4002) = {2301, 4002}; @@ -1213,10 +1213,12 @@ Mesh.SmoothNormals = 1; /* Physical */ // Tips: -Physical Line("teTipW") = {4092,4091,4094,2041,4003}; +Physical Line("teW") = {4094,2041}; Physical Line("wakeTipW") = {4003}; -Physical Line("teTipT") = {4103,4102,4101,4100,4105,4106}; +Physical Line("wakeExW") = {4092, 4091, 4003}; +Physical Line("teT") = {4105}; Physical Line("wakeTipT") = {4106}; +Physical Line("wakeExT") = {4102,4101,4100, 4103, 4106}; // Symmetry: Physical Surface("symmetry") = {1}; diff --git a/dart/models/wht.geo b/dart/models/wht.geo index 53a64ea4da10f4f750a2dcfc12cdd35e66ceee8a..4fb94de3675ce307bef1b05c6334ba0bce12b3b5 100644 --- a/dart/models/wht.geo +++ b/dart/models/wht.geo @@ -30,7 +30,7 @@ zOfst = 0.56; // Domain Point(1) = {-10, 0, -10}; Point(2) = {17, 0, -10}; -Point(4001) = {17, 0, -1.39}; +Point(4001) = {17, 0, -0.03763}; Point(4011) = {17, 0, 1.06}; Point(3) = {17, 0, 10}; Point(4) = {-10, 0, 10}; @@ -833,8 +833,8 @@ Line(3026) = {3049, 3149}; /// Wake // Wing -Point(4002) = {17, 6.1, 0.09}; -Point(4003) = {17, 7, 0.29}; +Point(4002) = {17, 6.1, 0.5764}; +Point(4003) = {17, 7, 0.8177}; Line(4001) = {2001, 4001}; Line(4002) = {2301, 4002}; @@ -988,9 +988,9 @@ Mesh.SmoothNormals = 1; /* Physical */ // Tips: -Physical Line("teTipW") = {2031, 2041, 4003}; +Physical Line("teW") = {2031, 2041}; Physical Line("wakeTipW") = {4003}; -Physical Line("teTipT") = {3021, 4012}; +Physical Line("teT") = {3021}; Physical Line("wakeTipT") = {4012}; // Symmetry: diff --git a/dart/src/wAdjoint.cpp b/dart/src/wAdjoint.cpp index 00e62aa3463dcd1ffe109e2bd7c81b572e91399d..23fe93fbda6af18f0484d7a5c37276ab51d8fff1 100644 --- a/dart/src/wAdjoint.cpp +++ b/dart/src/wAdjoint.cpp @@ -117,8 +117,8 @@ void Adjoint::run() dCf_X[1][i] -= dCd_XU[i]; // dCd/dX - dRu/dX^T * lambdaCdU } // compute mesh adjoint (totals) - std::vector<double> lamClMsh = this->solveMesh(dCf_X[0]); // lambdaClX = dRx/dX^T * (dCl/dX - dRu/dX^T * lambdaClU) - std::vector<double> lamCdMsh = this->solveMesh(dCf_X[1]); // lambdaCdX = dRx/dX^T * (dCd/dX - dRu/dX^T * lambdaCdU) + std::vector<double> lamClMsh = this->solveMesh(dCf_X[0]); // lambdaClX = dRx/dX^-T * (dCl/dX - dRu/dX^T * lambdaClU) + std::vector<double> lamCdMsh = this->solveMesh(dCf_X[1]); // lambdaCdX = dRx/dX^-T * (dCd/dX - dRu/dX^T * lambdaCdU) // store solution for (auto n : sol->pbl->msh->nodes) { diff --git a/dart/src/wKuttaElement.cpp b/dart/src/wKuttaElement.cpp index c66a8de5ddd78f4c6dc934f6c9af502b473dae9a..360efa5eb6a1e738558a3818c7024e5a20ede6d3 100644 --- a/dart/src/wKuttaElement.cpp +++ b/dart/src/wKuttaElement.cpp @@ -28,26 +28,57 @@ KuttaElement::KuttaElement(size_t _no, Element *&_surE, Element *&_volE, // Sanity checks if (surE->type() == ELTYPE::LINE2) nDim = 2; + else if (surE->type() == ELTYPE::TRI3) + nDim = 3; else { std::stringstream err; err << "KuttaElement only implemented for surface elements of type " - << ELTYPE::LINE2 << " (" << surE->type() << " was given)!\n"; + << ELTYPE::LINE2 << " or " << ELTYPE::TRI3 + << " (" << surE->type() << " was given)!\n"; throw std::runtime_error(err.str()); } - if (volE->type() == ELTYPE::TRI3) + if (volE->type() == ELTYPE::TRI3 || volE->type() == ELTYPE::TETRA4) nCol = volE->nodes.size(); else { std::stringstream err; err << "KuttaElement only implemented for volume elements of type " - << ELTYPE::TRI3 << " (" << volE->type() << " was given)!\n"; + << ELTYPE::TRI3 << " or " << ELTYPE::TETRA4 + << " (" << volE->type() << " was given)!\n"; throw std::runtime_error(err.str()); } // Count TE nodes + mapNodes(); nRow = teN.size(); } +/** + * @brief Map the wake nodes to their local volume node indices + */ +void KuttaElement::mapNodes() +{ + // Find local volume index of contributing surface nodes + for (auto p : teN) + { + for (size_t i = 0; i < volE->nodes.size(); ++i) + { + if (surE->nodes[p.first] == volE->nodes[i]) + { + vsMap[p.first] = i; + break; + } + } + } + // Check that all contributing surface nodes have been mapped + if (vsMap.size() != teN.size()) + { + std::stringstream err; + err << "dart::KuttaElement: could not map some contributing nodes from surface element " << *surE << " to " << *volE << " (number of mapped nodes: " << vsMap.size() << ")!\n"; + throw std::runtime_error(err.str()); + } +} + void KuttaElement::write(std::ostream &out) const { out << "KuttaElement #" << no diff --git a/dart/src/wKuttaElement.h b/dart/src/wKuttaElement.h index b48c48c952139c0e5e722dcab6c24ec9b50fc5c6..f3fe8d8e27f873b399378dc0095cfe878b8873a5 100644 --- a/dart/src/wKuttaElement.h +++ b/dart/src/wKuttaElement.h @@ -21,6 +21,7 @@ #include "wObject.h" #include "wElement.h" #include <vector> +#include <map> namespace dart { @@ -29,16 +30,12 @@ namespace dart * @brief Super-element to apply Kutta condition * * Made up of one surface element at the trailing edge and the connected volume - * element. Only Line2 (surface) and Tri3 (volume) elements are curently - * implemented. - * Since Tri3 are used, the evaluation of the volume integrals is performed at - * dummy Gauss point #0. If volume elements with non-constant shape function - * derivatives are used, those derivates should be recomputed at the true - * surface element Gauss points. - * Kutta works well in 2D, but fails in 3D unless the geometry follows strict - * restrictions (thin and sharp TE, simple planform, same mesh elements on TE - * suction and pressure side, ...). Kutta contribution is therefore disabled for - * 3D cases. + * element. Only Line2/Tri3 (surface) and Tri3/Tetra4 (volume) elements are + * currently implemented. + * Since Tri3 or Tetra4 are used, the evalutaion of the volume integrals is + * performed at dummy Gauss point #0. If volume elements with non-constant shape + * function derivatives are used, those derivates should be recomputed at the + * true surface element Gauss points. * @authors Adrien Crovato */ class DART_API KuttaElement : public fwk::wObject @@ -54,12 +51,16 @@ public: size_t nDim; ///< dimension of volume element std::vector<std::pair<size_t, tbox::Node *>> teN; ///< Map of local node indices and trailing edge nodes + std::map<size_t, size_t> vsMap; ///< map between surface node and volume node indices int sign; ///< +1 if upper element, -1 if lower element KuttaElement(size_t _no, tbox::Element *&_surE, tbox::Element *&_volE, std::vector<std::pair<size_t, tbox::Node *>> &_teN, int _sign); virtual void write(std::ostream &out) const override; #endif + +private: + void mapNodes(); }; } // namespace dart diff --git a/dart/src/wKuttaResidual.cpp b/dart/src/wKuttaResidual.cpp index cbca6a696e549756a1c58bb41224f1222691120e..6d0b9e3b806bd3ca058e16499b9adc56284465db 100644 --- a/dart/src/wKuttaResidual.cpp +++ b/dart/src/wKuttaResidual.cpp @@ -27,18 +27,30 @@ using namespace dart; /** * @brief Build the residual matrix for a fixed-point iteration, on one Kutta element * - * A = \int psi * v_u*grad_phi dS + * A = 1/S * \int (psi + grad_psi*v_inf) * v*grad_phi dS + * NB: v_inf = [1, 0, 0] */ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<double> const &phi) { - // Get shape functions, Jacobian and Gauss points + // Get shape functions and Gauss points Cache &sCache = ke.surE->getVCache(); Gauss &sGauss = sCache.getVGauss(); + Eigen::MatrixXd const &dSfV = ke.volE->getVCache().getDsf(0); + // Get Jacobian + Eigen::MatrixXd const &iJ = ke.volE->getJinv(0); // Get size size_t nRow = ke.nRow; + size_t nDim = ke.nDim; + // Stabilization + Eigen::MatrixXd dSf(nRow, nDim); + for (size_t i = 0; i < nRow; ++i) + for (size_t j = 0; j < nDim; ++j) + dSf(i, j) = dSfV(j, ke.vsMap.at(ke.teN[i].first)); + Eigen::VectorXd dSfp(nRow); + dSfp = 0.5 * sqrt(ke.surE->getVol()) * dSf * iJ.transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim); // Velocity - Eigen::RowVectorXd Aup = ke.volE->computeGradient(phi, 0).transpose() * ke.volE->getJinv(0) * ke.volE->getVCache().getDsf(0); + Eigen::RowVectorXd Aup = 1 / ke.surE->getVol() * ke.volE->computeGradient(phi, 0).transpose() * iJ * dSfV; // Build Eigen::MatrixXd A = Eigen::MatrixXd::Zero(nRow, ke.nCol); @@ -48,7 +60,7 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do Eigen::VectorXd const &sf = sCache.getSf(k); Eigen::VectorXd sfp(nRow); for (size_t i = 0; i < nRow; ++i) - sfp(i) = sf(ke.teN[i].first); + sfp(i) = sf(ke.teN[i].first) + dSfp(i); // wake contribution A += ke.sign * sfp * Aup * sGauss.getW(k) * ke.surE->getDetJ(k); } @@ -58,19 +70,30 @@ Eigen::MatrixXd KuttaResidual::buildFixed(KuttaElement const &ke, std::vector<do /** * @brief Build the residual vector, on one Kutta element * - * b = \int psi * (grad_phi)^2 dS + * b = 1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 dS + * NB: v_inf = [1, 0, 0] */ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double> const &phi) { - // Get shape functions, Jacobian and Gauss points + // Get shape functions and Gauss points Cache &sCache = ke.surE->getVCache(); Gauss &sGauss = sCache.getVGauss(); + Eigen::MatrixXd const &dSfV = ke.volE->getVCache().getDsf(0); + // Get Jacobian + Eigen::MatrixXd const &iJ = ke.volE->getJinv(0); // Get size size_t nRow = ke.nRow; + size_t nDim = ke.nDim; + // Stabilization + Eigen::MatrixXd dSf(nRow, nDim); + for (size_t i = 0; i < nRow; ++i) + for (size_t j = 0; j < nDim; ++j) + dSf(i, j) = dSfV(j, ke.vsMap.at(ke.teN[i].first)); + Eigen::VectorXd dSfp(nRow); + dSfp = 0.5 * sqrt(ke.surE->getVol()) * dSf * iJ.transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim); // Velocity - Eigen::VectorXd dPhi = ke.volE->computeGradient(phi, 0); - double dPhi2 = dPhi.squaredNorm(); + double dPhi2 = 1 / ke.surE->getVol() * ke.volE->computeGradient(phi, 0).squaredNorm(); // Build Eigen::VectorXd b = Eigen::VectorXd::Zero(nRow); @@ -80,7 +103,7 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double> Eigen::VectorXd const &sf = sCache.getSf(k); Eigen::VectorXd sfp(nRow); for (size_t i = 0; i < nRow; ++i) - sfp(i) = sf(ke.teN[i].first); + sfp(i) = sf(ke.teN[i].first) + dSfp(i); // wake contribution b += ke.sign * sfp * dPhi2 * sGauss.getW(k) * ke.surE->getDetJ(k); } @@ -90,7 +113,8 @@ Eigen::VectorXd KuttaResidual::build(KuttaElement const &ke, std::vector<double> /** * @brief Build the gradient of the residual vector with respect to the flow variable (jacobian), on one Kutta element * - * A = \int psi * 2*(grad_phi_u)*dgrad_phi_u dS + * A = 1/S * \int (psi + grad_psi*v_inf) * 2*(grad_phi_u)*dgrad_phi_u dS + * NB: v_inf = [1, 0, 0] */ Eigen::MatrixXd KuttaResidual::buildGradientFlow(KuttaElement const &ke, std::vector<double> const &phi) { @@ -101,27 +125,62 @@ Eigen::MatrixXd KuttaResidual::buildGradientFlow(KuttaElement const &ke, std::ve /** * @brief Build the gradient of the residual vector with respect to the nodes, on one Kutta element - * - * A = d( \int psi * (grad_phi)^2 dS ) - * = \int psi * d(grad_phi)^2 dS - * + \int psi * (grad_phi)^2 ddS + * + * A = d( 1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 dS ) + * = d1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 dS + * + 1/S * \int dgrad_psi*v_inf * (grad_phi)^2 dS + * + 1/S * \int (psi + grad_psi*v_inf) * d(grad_phi)^2 dS + * + 1/S * \int (psi + grad_psi*v_inf) * (grad_phi)^2 ddS + * NB: v_inf = [1, 0, 0] */ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(KuttaElement const &ke, std::vector<double> const &phi) { // Get shape functions and Gauss points Cache &surCache = ke.surE->getVCache(); Gauss &surGauss = surCache.getVGauss(); + Eigen::MatrixXd const &dSfV = ke.volE->getVCache().getDsf(0); // Get Jacobian Eigen::MatrixXd const &iJ = ke.volE->getJinv(0); std::vector<Eigen::MatrixXd> const &dJ = ke.volE->getGradJ(0); + // Get surface and gradient + double surf = ke.surE->getVol(); + std::vector<double> const &dSurf = ke.surE->getGradVol(); // Get size size_t nRow = ke.nRow; size_t nDim = ke.nDim; size_t nCol = ke.nCol; size_t nSur = ke.surE->nodes.size(); - // velocity + // Stabilization term + Eigen::VectorXd vi = Eigen::Vector3d(1, 0, 0).middleRows(0, nDim); + double sqSurf = sqrt(surf); + Eigen::MatrixXd dSf(nRow, nDim); + for (size_t i = 0; i < nRow; ++i) + for (size_t m = 0; m < nDim; ++m) + dSf(i, m) = dSfV(m, ke.vsMap.at(ke.teN[i].first)); + Eigen::VectorXd gradSfp = 0.5 * sqSurf * dSf * iJ.transpose() * vi; + // Gradient of stabilization term (volume and surface) + std::vector<Eigen::VectorXd> dGradSfV(nDim * nCol); + for (size_t i = 0; i < nCol; ++i) + { + for (size_t m = 0; m < nDim; ++m) + { + size_t idx = i * nDim + m; + dGradSfV[idx] = 0.5 * sqSurf * dSf * (-iJ * dJ[idx] * iJ).transpose() * vi; // S * dgrad_psi + } + } + std::vector<Eigen::VectorXd> dGradSfS(nDim * nSur); + for (size_t i = 0; i < nSur; ++i) + { + for (size_t m = 0; m < nDim; ++m) + { + size_t idx = i * nDim + m; + dGradSfS[idx] = 0.5 * 0.5 / sqSurf * dSurf[idx] * dSf * iJ.transpose() * vi; // dS * grad_psi + } + } + // Velocity and normalization Eigen::VectorXd dPhi = ke.volE->computeGradient(phi, 0); + double nrm = 1 / surf; // Build Eigen::MatrixXd Av = Eigen::MatrixXd::Zero(nRow, nDim * nCol); @@ -136,7 +195,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku Eigen::VectorXd const &sf = surCache.getSf(k); Eigen::VectorXd sfp(nRow); for (size_t i = 0; i < nRow; ++i) - sfp(i) = sf(ke.teN[i].first); + sfp(i) = sf(ke.teN[i].first) + gradSfp(i); // Volume for (size_t i = 0; i < nCol; ++i) @@ -145,7 +204,9 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku { size_t idx = i * nDim + m; // gradient of velocity jump - Av.col(idx) += ke.sign * w * sfp * 2 * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // psi * dgrad_phi^2 + Av.col(idx) += ke.sign * nrm * w * sfp * 2 * dPhi.transpose() * (-iJ * dJ[idx] * dPhi) * detJ; // 1/S * (psi + grad_psi) * dgrad_phi^2 * detJ + // gradient of stabilization term (volume) + Av.col(idx) += ke.sign * nrm * w * dGradSfV[idx] * dPhi.dot(dPhi) * detJ; // 1/S * dgrad_psi * grad_phi^2 * detJ } } // Surface @@ -154,8 +215,12 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> KuttaResidual::buildGradientMesh(Ku for (size_t m = 0; m < nDim; ++m) { size_t idx = i * nDim + m; + // gradient of normalization + As.col(idx) -= ke.sign * nrm * nrm * dSurf[idx] * w * sfp * dPhi.dot(dPhi) * detJ; // d1/S * (psi + grad_psi) * grad_phi^2 * detJ + // gradient of stabilization term (surface) + As.col(idx) += ke.sign * nrm * w * dGradSfS[idx] * dPhi.dot(dPhi) * detJ; // 1/S * dgrad_psi * grad_phi^2 * detJ // gradient of jacobian determinant - As.col(idx) += ke.sign * w * sfp * dPhi.dot(dPhi) * dDetJ[idx]; // psi * grad_phi^2 * ddetJ + As.col(idx) += ke.sign * nrm * w * sfp * dPhi.dot(dPhi) * dDetJ[idx]; // 1/S * (psi + grad_psi) * grad_phi^2 * ddetJ } } } diff --git a/dart/src/wProblem.cpp b/dart/src/wProblem.cpp index 73d6f912234591de19674f8ef6cd03e9275d95ee..06d219dc192c5e9d5e4ac9b96f53c1c26a7b9618 100644 --- a/dart/src/wProblem.cpp +++ b/dart/src/wProblem.cpp @@ -194,63 +194,42 @@ void Problem::check() const // Sanity checks if (fluid == NULL) throw std::runtime_error("No Fluid provided!\n"); + if (bodies.empty()) + throw std::runtime_error("No Body provided!\n"); if (iIC == NULL) throw std::runtime_error("No Initial Condition provided!\n"); if (fBCs.empty()) throw std::runtime_error("No Freestream B.C. provided!\n"); if (wBCs.empty()) throw std::runtime_error("No Wake B.C. provided!\n"); + if (kSCs.empty()) + throw std::runtime_error("No Kutta condition provided!\n"); // Three-dimension problem if (nDim == 3) { // Volume element type for (auto e : fluid->tag->elems) - { if (e->type() != ELTYPE::TETRA4) - { - std::stringstream err; - err << "3D solver is only implemented for volume elements of type " - << ELTYPE::TETRA4 << " (" << e->type() << " was given)!\n"; - throw std::runtime_error(err.str()); - } - } + throwUnsupElType(e, "3", "volume", ELTYPE::TETRA4); + // Surface element type (Body) + for (auto surf : bodies) + for (auto e : surf->groups[0]->tag->elems) + if (e->type() != ELTYPE::TRI3) + throwUnsupElType(e, "3", "surface", ELTYPE::TRI3); // Surface element type (Freestream B.C.) for (auto surf : fBCs) - { for (auto e : surf->tag->elems) - { if (e->type() != ELTYPE::TRI3) - { - std::stringstream err; - err << "3D solver is only implemented for surface elements of type " - << ELTYPE::TRI3 << " (" << e->type() << " was given)!\n"; - throw std::runtime_error(err.str()); - } - } - } + throwUnsupElType(e, "3", "surface", ELTYPE::TRI3); // Surface element type (Wake B.C.) for (auto surf : wBCs) { for (auto e : surf->groups[0]->tag->elems) - { if (e->type() != ELTYPE::TRI3) - { - std::stringstream err; - err << "3D solver is only implemented for surface elements of type " - << ELTYPE::TRI3 << " (" << e->type() << " was given)!\n"; - throw std::runtime_error(err.str()); - } - } + throwUnsupElType(e, "3", "surface", ELTYPE::TRI3); for (auto e : surf->groups[1]->tag->elems) - { if (e->type() != ELTYPE::TRI3) - { - std::stringstream err; - err << "3D solver is only implemented for surface elements of type " - << ELTYPE::TRI3 << " (" << e->type() << " was given)!\n"; - throw std::runtime_error(err.str()); - } - } + throwUnsupElType(e, "3", "surface", ELTYPE::TRI3); } // Blowing B.C. if (!bBCs.empty()) @@ -261,52 +240,27 @@ void Problem::check() const { // Volume element type for (auto e : fluid->tag->elems) - { if (e->type() != ELTYPE::TRI3) - { - std::stringstream err; - err << "2D solver is only implemented for volume elements of type " - << ELTYPE::TRI3 << " (" << e->type() << " was given)!\n"; - throw std::runtime_error(err.str()); - } - } + throwUnsupElType(e, "2", "volume", ELTYPE::TRI3); + // Surface element type (Body) + for (auto surf : bodies) + for (auto e : surf->groups[0]->tag->elems) + if (e->type() != ELTYPE::LINE2) + throwUnsupElType(e, "2", "surface", ELTYPE::LINE2); // Surface element type (Freestream B.C.) for (auto surf : fBCs) - { for (auto e : surf->tag->elems) - { if (e->type() != ELTYPE::LINE2) - { - std::stringstream err; - err << "2D solver is only implemented for surface elements of type " - << ELTYPE::LINE2 << " (" << e->type() << " was given)!\n"; - throw std::runtime_error(err.str()); - } - } - } + throwUnsupElType(e, "2", "surface", ELTYPE::LINE2); // Surface element type (Wake B.C.) for (auto surf : wBCs) { for (auto e : surf->groups[0]->tag->elems) - { if (e->type() != ELTYPE::LINE2) - { - std::stringstream err; - err << "2D solver is only implemented for surface elements of type " - << ELTYPE::LINE2 << " (" << e->type() << " was given)!\n"; - throw std::runtime_error(err.str()); - } - } + throwUnsupElType(e, "2", "surface", ELTYPE::LINE2); for (auto e : surf->groups[1]->tag->elems) - { if (e->type() != ELTYPE::LINE2) - { - std::stringstream err; - err << "2D solver is only implemented for surface elements of type " - << ELTYPE::LINE2 << " (" << e->type() << " was given)!\n"; - throw std::runtime_error(err.str()); - } - } + throwUnsupElType(e, "2", "surface", ELTYPE::LINE2); } } else @@ -318,6 +272,15 @@ void Problem::check() const } } +void Problem::throwUnsupElType(Element const *e, std::string const &pdim, std::string const &edim, ELTYPE type) const +{ + std::stringstream err; + err << pdim << "D solver is only implemented for " + << edim << " elements of type " << type + << " (" << e->type() << " was given)!\n"; + throw std::runtime_error(err.str()); +} + void Problem::write(std::ostream &out) const { out << "dart::Problem parameters" diff --git a/dart/src/wProblem.h b/dart/src/wProblem.h index 3ba02003b099911f8da965d56dd85c4401f50e17..d59911e6863d6db28eafc223371731f2095517b5 100644 --- a/dart/src/wProblem.h +++ b/dart/src/wProblem.h @@ -19,6 +19,7 @@ #include "dart.h" #include "wObject.h" +#include "wElement.h" #include "wF1Ct.h" #include <memory> #include <iostream> @@ -87,6 +88,9 @@ public: void initGradElems(); virtual void write(std::ostream &out) const override; #endif + +private: + void throwUnsupElType(tbox::Element const *e, std::string const &pdim, std::string const &edim, tbox::ELTYPE type) const; }; } // namespace dart diff --git a/dart/src/wSolver.cpp b/dart/src/wSolver.cpp index 7f69f841eb1f6c7357f6c3ef3a2c0f2be733ec97..3b11893a558e00df78cb679bf60e2dc48fc7535e 100644 --- a/dart/src/wSolver.cpp +++ b/dart/src/wSolver.cpp @@ -21,6 +21,7 @@ #include "wAssign.h" #include "wBody.h" #include "wWake.h" +#include "wKutta.h" #include "wLoadFunctional.h" #include "wMshData.h" @@ -110,6 +111,8 @@ Solver::Solver(std::shared_ptr<Problem> _pbl, std::cout << "Dirichlet on " << *bnd->tag << "\n"; for (auto bnd : pbl->wBCs) std::cout << "Wake on " << *bnd->groups[0]->tag << "\n"; + for (auto bnd : pbl->kSCs) + std::cout << "Kutta on " << *bnd->groups[0]->tag << "\n"; for (auto bnd : pbl->bodies) std::cout << "Body on " << *bnd->groups[0]->tag << "\n"; std::cout << "--- Freestream conditions ---\n" diff --git a/dart/src/wWakeElement.h b/dart/src/wWakeElement.h index 0cf0a30cb76206fad57a05296069a4fa4ff18b2b..43ead69dca1bbbef7c1d4651addf02801baf0d72 100644 --- a/dart/src/wWakeElement.h +++ b/dart/src/wWakeElement.h @@ -30,8 +30,8 @@ namespace dart * @brief Super-element to apply wake boundary conditions * * Made up of two surface elements on the wake and the connected volume - * elements. Only Line2 (surface) and Tri3 (volume) elements are curently - * implemented. + * elements. Only Line2/Tri3 (surface) and Tri3/Tetra4 (volume) elements are + * currently implemented. * Since Tri3 or Tetra4 are used, the evalutaion of the volume integrals is * performed at dummy Gauss point #0. If volume elements with non-constant shape * function derivatives are used, those derivates should be recomputed at the diff --git a/dart/src/wWakeResidual.cpp b/dart/src/wWakeResidual.cpp index 671eb02674d119592d428eb2959ba23f74edd740..4af31768ca2f18c7370533247264bdba82f190e2 100644 --- a/dart/src/wWakeResidual.cpp +++ b/dart/src/wWakeResidual.cpp @@ -66,7 +66,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> WakeResidual::buildFixed(WakeElemen Eigen::VectorXd const &sf = sCache.getSf(k); Eigen::VectorXd sfp(nRow); for (size_t i = 0; i < nRow; ++i) - sfp(i) = dSfp(i) + sf(we.wkN[i].first); + sfp(i) = sf(we.wkN[i].first) + dSfp(i); // wake contribution Au += sfp * Aup * sGauss.getW(k) * we.surUpE->getDetJ(k); Al += sfp * Alw * sGauss.getW(k) * we.surUpE->getDetJ(k); @@ -99,10 +99,8 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con Eigen::VectorXd dSfp(nRow); dSfp = 0.5 * sqrt(we.surUpE->getVol()) * dSf * we.volUpE->getJinv(0).transpose() * Eigen::Vector3d(1, 0, 0).middleRows(0, nDim); // Velocity squared - Eigen::VectorXd dPhiU = we.volUpE->computeGradient(phi, 0); - Eigen::VectorXd dPhiL = we.volLwE->computeGradient(phi, 0); - double dPhi2u = dPhiU.squaredNorm(); - double dPhi2l = dPhiL.squaredNorm(); + double dPhi2u = we.volUpE->computeGradient(phi, 0).squaredNorm(); + double dPhi2l = we.volLwE->computeGradient(phi, 0).squaredNorm(); // Build Eigen::VectorXd bu = Eigen::VectorXd::Zero(nRow); @@ -113,7 +111,7 @@ std::tuple<Eigen::VectorXd, Eigen::VectorXd> WakeResidual::build(WakeElement con Eigen::VectorXd const &sf = sCache.getSf(k); Eigen::VectorXd sfp(nRow); for (size_t i = 0; i < nRow; ++i) - sfp(i) = dSfp(i) + sf(we.wkN[i].first); + sfp(i) = sf(we.wkN[i].first) + dSfp(i); // wake contribution bu += sfp * dPhi2u * sGauss.getW(k) * we.surUpE->getDetJ(k); bl += sfp * dPhi2l * sGauss.getW(k) * we.surUpE->getDetJ(k); diff --git a/dart/tests/adjoint.py b/dart/tests/adjoint.py index 3d3932a85195fe4d8199e67ff40562b8f3ef51ef..6d35b848abe2b85bef09d01625878dbfce1115f8 100644 --- a/dart/tests/adjoint.py +++ b/dart/tests/adjoint.py @@ -55,7 +55,7 @@ def main(): # set the problem and add fluid, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil') tms['pre'].stop() # solve problem @@ -114,17 +114,17 @@ def main(): if M_inf == 0. and alpha == 0*np.pi/180: tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 6.9, 1e-2)) tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.0, 1e-3)) - tests.add(CTest('dCl_dMsh', dClX, 55.439, 1e-3)) - tests.add(CTest('dCd_dMsh', dCdX, 0.482, 1e-3)) - tests.add(CTest('dCl_dAoA (msh)', dClAoA, 6.9, 1e-2)) + tests.add(CTest('dCl_dMsh', dClX, 55.1, 1e-2)) + tests.add(CTest('dCd_dMsh', dCdX, 0.482, 1e-2)) + tests.add(CTest('dCl_dAoA (msh)', dClAoA, 7.1, 1e-2)) # FD 7.1364 (step = 1e-6) tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.0, 1e-3)) elif M_inf == 0.7 and alpha == 2*np.pi/180: tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 11.8, 1e-2)) # FD 11.835 (step = 1e-5) tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.127, 1e-2)) # 0.12692 (step = 1e-5) - tests.add(CTest('dCl_dMsh', dClX, 82.67, 1e-3)) - tests.add(CTest('dCd_dMsh', dCdX, 0.908, 1e-3)) - tests.add(CTest('dCl_dAoA (msh)', dClAoA, 11.8, 1e-2)) - tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.127, 1e-2)) + tests.add(CTest('dCl_dMsh', dClX, 83.6, 1e-2)) + tests.add(CTest('dCd_dMsh', dCdX, 0.904, 1e-2)) + tests.add(CTest('dCl_dAoA (msh)', dClAoA, 12.1, 1e-2)) # FD 12.1443 (step = 1e-6) + tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.129, 1e-2)) # FD 0.1289 (step = 1e-6) else: raise Exception('Test not defined for this flow') tests.run() diff --git a/dart/tests/bli.py b/dart/tests/bli.py index fd2f186f340c56719e664fcdca3dc8edb89fed6b..9014bcc0ab9feac43ed1dcb451d597950fd2cf60 100644 --- a/dart/tests/bli.py +++ b/dart/tests/bli.py @@ -57,7 +57,7 @@ def main(): M_inf = 0.5 c_ref = 1 dim = 2 - tol = 1e-4 # tolerance for the VII (usually one drag count) + tol = 1e-3 # tolerance for the VII (usually one drag count) #TODO tolerance has been decreased with new Kutta, hopefully will increase with Paul's dev # mesh the geometry print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) @@ -68,7 +68,7 @@ def main(): # set the problem and add fluid, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True) + pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', vsc = True, dbc=True) tms['pre'].stop() # solve viscous problem diff --git a/dart/tests/cylinder.py b/dart/tests/cylinder.py index 88dc36af29ef3ae7e7955fb12e128d880625c8d8..2d5f07835f07e2607f88aa357f6a434ee4292b7b 100644 --- a/dart/tests/cylinder.py +++ b/dart/tests/cylinder.py @@ -59,7 +59,7 @@ def main(): # set the problem and add fluid, initial/boundary conditions tms['pre'].start() - pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, l_ref, l_ref, 0., 0., 0., 'cylinder', te = 'te', dbc=True) + pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, l_ref, l_ref, 0., 0., 0., 'cylinder', dbc=True) tms['pre'].stop() # solve problem diff --git a/dart/tests/cylinder2D5.py b/dart/tests/cylinder2D5.py deleted file mode 100644 index 107e6a4c8867a1ce083a9b647300ba4250ff18cb..0000000000000000000000000000000000000000 --- a/dart/tests/cylinder2D5.py +++ /dev/null @@ -1,134 +0,0 @@ -#!/usr/bin/env python3 -# -*- 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. - - -## Compute the (non)linear flow around extruded 2D cylinder, with various B.C. -# Adrien Crovato -# -# Test solver convergence and boundary conditions -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import math -import dart.default as floD -import fwk -from fwk.testing import * -from fwk.coloring import ccolors - -def main(): - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - alpha = 3*math.pi/180 - U_inf = [math.cos(alpha), 0., math.sin(alpha)] # norm should be 1 - M_inf = 0.1 - dim = 3 - - # define dimension and mesh size - dmt = 1.0 # cylinder diameter - lgt = 8.0 # channel length - hgt = 6.0 # channel height - wdt = 3.0 # channel width - l_ref = dmt # reference length - S_ref = dmt*wdt # reference area - fms = 0.5 # farfield mesh size - nms = 0.1 # nearfield mesh size - - # mesh the geometry - print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) - tms["msh"].start() - pars = {'dmt' : dmt, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms} - msh, gmshWriter = floD.mesh(dim, 'models/cylinder_2D5.geo', pars, ['field', 'cylinder', 'symmetry', 'downstream']) - tms["msh"].stop() - - # set the problem and add fluid, initial/boundary conditions - tms['pre'].start() - pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'te', dbc=True) - tms['pre'].stop() - - # solve problem - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - solver = floD.newton(pbl) - cnvrgd = solver.run() - solver.save(gmshWriter) - tms['solver'].stop() - - # compute mean error in Body, Neumann boundary condition (works only for Tri3 elements) - gradn = [] - for e in bnd.groups[0].tag.elems: - gradx = 0 - grady = 0 - gradz = 0 - ux = e.nodes[1].pos[0] - e.nodes[0].pos[0] - uy = e.nodes[1].pos[1] - e.nodes[0].pos[1] - uz = e.nodes[1].pos[2] - e.nodes[0].pos[2] - vx = e.nodes[2].pos[0] - e.nodes[0].pos[0] - vy = e.nodes[2].pos[1] - e.nodes[0].pos[1] - vz = e.nodes[2].pos[2] - e.nodes[0].pos[2] - nx = uy*vz - uz*vy - ny = uz*vx - ux*vz - nz = ux*vy - uy*vx - for n in e.nodes: - gradx += solver.U[n.row][0] / e.nodes.size() - grady += solver.U[n.row][1] / e.nodes.size() - gradz += solver.U[n.row][2] / e.nodes.size() - gradn_ = gradx * nx + grady * ny + gradz * nz - gradn.append(gradn_) - errNeu = sum(gradn) / len(gradn) - # compute largest error in Dirichlet boundary condition - errDir = 0. - for n in dirichlet.nodes: - errCur = solver.phi[n.row] - (n.pos[0]*U_inf[0] + n.pos[1]*U_inf[1] + n.pos[2]*U_inf[2]) - if (abs(errCur) > errDir): - errDir = abs(errCur) - # compute largest error in Wake boundary condition - errWak = 0. - for n in wake.nodMap: - errCur = solver.Cp[n.row] - solver.Cp[wake.nodMap[n].row] - if (abs(errCur) > errWak): - errWak = abs(errCur) - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - - # visualize solution - floD.initViewer(pbl) - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - if (not cnvrgd): - raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET) - tests = CTests() - tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True)) - tests.add(CTest('Neumann B.C. mean error', errNeu, 0., 1e-2)) - tests.add(CTest('Dirichlet B.C. max error', errDir, 0., 1e-12)) - tests.add(CTest('Wake/Kutta B.C. max error', errWak, 0., 1e-1)) - tests.run() - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/tests/cylinder3.py b/dart/tests/cylinder3.py deleted file mode 100644 index 869eca38e59673a42dd7db4e238425c57fc6dcdf..0000000000000000000000000000000000000000 --- a/dart/tests/cylinder3.py +++ /dev/null @@ -1,137 +0,0 @@ -#!/usr/bin/env python3 -# -*- 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. - - -## Compute the (non)linear flow around 3D cylinder, with various B.C. -# Adrien Crovato -# -# Test solver convergence and boundary conditions -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import math -import dart.utils as floU -import dart.default as floD -import fwk -from fwk.testing import * -from fwk.coloring import ccolors - -def main(): - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - alpha = 3*math.pi/180 - U_inf = [math.cos(alpha), 0., math.sin(alpha)] # norm should be 1 - M_inf = 0.1 - dim = 3 - - # define dimension and mesh size - dmt = 1.0 # cylinder diameter - spn = 2.0 # cylinder span - tpt = 5e-2 # cylinder tip thickness (ratio wrt span) - lgt = 8.0 # channel length - hgt = 6.0 # channel height - wdt = 3.0 # channel width - l_ref = dmt # reference length - S_ref = dmt*spn # reference area - fms = 0.5 # farfield mesh size - nms = 0.1 # nearfield mesh size - - # mesh the geometry - print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) - tms['msh'].start() - pars = {'dmt' : dmt, 'spn' : spn, 'tpt' : tpt, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms} - msh, gmshWriter = floD.mesh(dim, 'models/cylinder_3.geo', pars, ['field', 'cylinder', 'symmetry', 'downstream'], wktp = 'wakeTip') - tms['msh'].stop() - - # set the problem and add fluid, initial/boundary conditions - tms['pre'].start() - pbl, dirichlet, wake, bnd,_ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'teTip', dbc=True) - tms['pre'].stop() - - # solve problem - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - solver = floD.newton(pbl) - cnvrgd = solver.run() - solver.save(gmshWriter) - tms['solver'].stop() - - # compute mean error in Body, Neumann boundary condition (works only for Tri3 elements) - gradn = [] - for e in bnd.groups[0].tag.elems: - gradx = 0 - grady = 0 - gradz = 0 - ux = e.nodes[1].pos[0] - e.nodes[0].pos[0] - uy = e.nodes[1].pos[1] - e.nodes[0].pos[1] - uz = e.nodes[1].pos[2] - e.nodes[0].pos[2] - vx = e.nodes[2].pos[0] - e.nodes[0].pos[0] - vy = e.nodes[2].pos[1] - e.nodes[0].pos[1] - vz = e.nodes[2].pos[2] - e.nodes[0].pos[2] - nx = uy*vz - uz*vy - ny = uz*vx - ux*vz - nz = ux*vy - uy*vx - for n in e.nodes: - gradx += solver.U[n.row][0] / e.nodes.size() - grady += solver.U[n.row][1] / e.nodes.size() - gradz += solver.U[n.row][2] / e.nodes.size() - gradn_ = gradx * nx + grady * ny + gradz * nz - gradn.append(gradn_) - errNeu = sum(gradn) / len(gradn) - # compute largest error in Dirichlet boundary condition - errDir = 0. - for n in dirichlet.nodes: - errCur = solver.phi[n.row] - (n.pos[0]*U_inf[0] + n.pos[1]*U_inf[1] + n.pos[2]*U_inf[2]) - if (abs(errCur) > errDir): - errDir = abs(errCur) - # compute largest error in Wake boundary condition - errWak = 0. - for n in wake.nodMap: - errCur = solver.Cp[n.row] - solver.Cp[wake.nodMap[n].row] - if (abs(errCur) > errWak): - errWak = abs(errCur) - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - - # visualize solution - floD.initViewer(pbl) - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - if (not cnvrgd): - raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET) - tests = CTests() - tests.add(CTest('iteration count', solver.nIt, 4, 1, forceabs=True)) - tests.add(CTest('Neumann B.C. mean error', errNeu, 0., 1e-2)) - tests.add(CTest('Dirichlet B.C. max error', errDir, 0., 1e-12)) - tests.add(CTest('Wake/Kutta B.C. max error', errWak, 0., 1e-1)) - tests.run() - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/tests/lift.py b/dart/tests/lift.py index 54f58b33149fe1038c74525592615c90a119789b..1c8ee82f1b62023c79550761527bd1288585ea9b 100644 --- a/dart/tests/lift.py +++ b/dart/tests/lift.py @@ -53,7 +53,7 @@ def main(): # set the problem and add fluid, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil') tms['pre'].stop() # solve problem @@ -97,6 +97,7 @@ def main(): tests.add(CTest('iteration count', solver.nIt, 12, 3, forceabs=True)) tests.add(CTest('min(Cp)', min(Cp[:,3]), -1.28, 5e-2)) tests.add(CTest('Cl', solver.Cl, 0.388, 5e-2)) + tests.add(CTest('Cd', solver.Cd, 0.0013, 5e-4, forceabs=True)) tests.add(CTest('Cm', solver.Cm, 0.0, 1e-2)) else: raise Exception('Test not defined for this flow') diff --git a/dart/tests/lift3.py b/dart/tests/lift3.py deleted file mode 100644 index 67c2ad3d74fdf53c8ed07166138bf2363f95ac80..0000000000000000000000000000000000000000 --- a/dart/tests/lift3.py +++ /dev/null @@ -1,101 +0,0 @@ -#!/usr/bin/env python3 -# -*- 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. - - -## Compute lifting (linear or nonlinear) potential flow around a NACA 0012 rectangular wing -# Adrien Crovato -# -# Test the nonlinear shock-capturing capability and the Kutta condition -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import math -import dart.utils as floU -import dart.default as floD -import fwk -from fwk.testing import * -from fwk.coloring import ccolors - -def main(): - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - alpha = 3*math.pi/180 - M_inf = 0.3 - dim = 3 - - # define dimension and mesh size - spn = 1.0 # wing span - lgt = 3.0 # channel half length - hgt = 3.0 # channel half height - wdt = 3.0 # channel width - S_ref = 1.*spn # reference area - c_ref = 1 # reference chord - fms = 1.0 # farfield mesh size - nms = 0.02 # nearfield mesh size - - # mesh the geometry - print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) - tms['msh'].start() - pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms} - msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip') - tms['msh'].stop() - - # set the problem and add fluid, initial/boundary conditions - tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') - tms['pre'].stop() - - # solve problem - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver'].start() - solver = floD.newton(pbl) - solver.run() - solver.save(gmshWriter) - tms['solver'].stop() - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - - # visualize solution - floD.initViewer(pbl) - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - tests = CTests() - if alpha == 3*math.pi/180 and M_inf == 0.3 and spn == 1.0: - tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True)) - tests.add(CTest('CL', solver.Cl, 0.135, 5e-2)) - tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035 - tests.add(CTest('CS', solver.Cs, 0.0121, 5e-2)) - tests.add(CTest('CM', solver.Cm, -0.0278, 1e-2)) - else: - raise Exception('Test not defined for this flow') - tests.run() - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/tests/meshDef3.py b/dart/tests/meshDef3.py deleted file mode 100644 index d571ee61ffbcd3b53f3dce91e427d21e49a2e55a..0000000000000000000000000000000000000000 --- a/dart/tests/meshDef3.py +++ /dev/null @@ -1,134 +0,0 @@ -#!/usr/bin/env python3 -# -*- 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. - - -## Compute the flow around a NACA 0012 wing on a deformed mesh -# Adrien Crovato -# -# Test the mesh deformation process -# -# CAUTION -# This test is provided to ensure that the solver works properly. -# Mesh refinement may have to be performed to obtain physical results. - -import numpy as np -import dart.utils as floU -import dart.default as floD -import fwk -from fwk.testing import * -from fwk.coloring import ccolors - -def main(): - # timer - tms = fwk.Timers() - tms['total'].start() - - # define flow variables - alpha0 = 0 # must be zero for this testfile - M_inf = 0.3 - dim = 3 - - # define dimension and mesh size - spn = 1.0 # wing span - lgt = 3.0 # channel half length - hgt = 3.0 # channel half height - wdt = 3.0 # channel width - S_ref = 1.*spn # reference area - c_ref = 1 # reference chord - fms = 1.0 # farfield mesh size - nms = 0.05 # nearfield mesh size - - # parameters for mesh deformation - alfa = 3*np.pi/180 - xc = 0.25 - dz_max = 0.2*spn - - # mesh an airfoil - print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) - tms['msh'].start() - pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms} - msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip') - tms['msh'].stop() - # create mesh deformation handler - mshDef = floD.morpher(msh, dim, ['wing']) - - # set the problem and add fluid, initial/boundary conditions - tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') - tms['pre'].stop() - - # solver problem for 0°AOA - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver0'].start() - solver = floD.newton(pbl) - solver.run() - tms['solver0'].stop() - - # deform the mesh (dummy rotation + bending) - print(ccolors.ANSI_BLUE + 'PyDeforming...' + ccolors.ANSI_RESET) - tms['deform'].start() - mshDef.savePos() - rMat = np.matrix([[np.cos(alfa), np.sin(alfa)], [-np.sin(alfa), np.cos(alfa)]]) - vNods = {} - for e in bnd.groups[0].tag.elems: - for n in e.nodes: - if n.no not in vNods: - vNods[n.no] = n - for no in vNods: - nod = np.matrix([[vNods[no].pos[0]-xc], [vNods[no].pos[2]]]) - new = rMat * nod - vNods[no].pos[0] = new[0,0] - vNods[no].pos[2] = new[1,0] + dz_max/(spn*spn)*(vNods[no].pos[1]*vNods[no].pos[1]) - # deform the mesh - mshDef.deform() - gmshWriter.save(msh.name+"_def") - tms['deform'].stop() - - # solve problem for alfa AOA - print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) - tms['solver1'].start() - solver.run() - solver.save(gmshWriter) - tms['solver1'].stop() - - # display results - print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) - print(' case M alpha Cl Cd Cm') - print(' true {0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}'.format(M_inf, alfa*180/np.pi, solver.Cl, solver.Cd, solver.Cm)) - - # display timers - tms['total'].stop() - print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) - print('CPU statistics') - print(tms) - - # check results - print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) - tests = CTests() - if alfa == 3*np.pi/180 and M_inf == 0.3 and spn == 1.0: - tests.add(CTest('iteration count', solver.nIt, 3, 1, forceabs=True)) - tests.add(CTest('CL', solver.Cl, 0.135, 5e-1)) - tests.add(CTest('CD', solver.Cd, 0.0062, 1e-2)) # Tranair (NF=0.0062, FF=0.0030), Panair 0.0035 - tests.add(CTest('CS', solver.Cs, -0.011, 5e-2)) - tests.add(CTest('CM', solver.Cm, 0.0061, 1e-2)) - tests.run() - - # eof - print('') - -if __name__ == "__main__": - main() diff --git a/dart/tests/meshDef.py b/dart/tests/morpher.py similarity index 96% rename from dart/tests/meshDef.py rename to dart/tests/morpher.py index c50e5af78d270351db7117ae710d3cb81246839f..c82f27719227c722d27d4db4d80d0d17d546d4d3 100644 --- a/dart/tests/meshDef.py +++ b/dart/tests/morpher.py @@ -55,7 +55,7 @@ def main(): pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01} msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream']) # create mesh deformation handler - mshDef = floD.morpher(msh, dim, ['airfoil']) + morpher = floD.morpher(msh, dim, ['airfoil']) # mesh the reference geometry pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01, 'angle' : alfa, 'xRot' : xc} msh_ref, gmshWriter_ref = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream']) @@ -63,9 +63,9 @@ def main(): # set the problem and add fluid, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, c_ref, c_ref, xc+dx, dy, 0., 'airfoil', te = 'te') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, c_ref, c_ref, xc+dx, dy, 0., 'airfoil') # set the refrence problem and add fluid, initial/boundary conditions - pbl_ref, _, _, bnd_ref, _ = floD.problem(msh_ref, dim, alpha0, 0., M_inf, c_ref, c_ref, xc, 0., 0., 'airfoil', te = 'te') + pbl_ref, _, _, bnd_ref, _ = floD.problem(msh_ref, dim, alpha0, 0., M_inf, c_ref, c_ref, xc, 0., 0., 'airfoil') tms['pre'].stop() # solver problem for 0°AOA @@ -78,7 +78,7 @@ def main(): # deform the mesh (dummy rotation + translation) print(ccolors.ANSI_BLUE + 'PyDeforming...' + ccolors.ANSI_RESET) tms['deform'].start() - mshDef.savePos() + morpher.savePos() rMat = np.matrix([[np.cos(alfa), np.sin(alfa)], [-np.sin(alfa), np.cos(alfa)]]) vNods = {} for e in bnd.groups[0].tag.elems: @@ -91,7 +91,7 @@ def main(): vNods[no].pos[0] = new[0,0] + dx + xc vNods[no].pos[1] = new[1,0] + dy # deform the mesh - mshDef.deform() + morpher.deform() gmshWriter.save("n0012_def") tms['deform'].stop() diff --git a/dart/tests/nonlift.py b/dart/tests/nonlift.py index 06136600610981edca7e4fa6350b10172aafdb0e..2710fcbc16126072583f14e4a928b3301b71d615 100644 --- a/dart/tests/nonlift.py +++ b/dart/tests/nonlift.py @@ -54,7 +54,7 @@ def main(): # set the problem and add fluid, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil') tms['pre'].stop() # solve problem @@ -93,6 +93,7 @@ def main(): elif M_inf == 0.8: tests.add(CTest('iteration count', solver.nIt, 15, 3, forceabs=True)) tests.add(CTest('min(Cp)', min(Cp[:,3]), -0.89, 5e-2)) + tests.add(CTest('Cd', solver.Cd, 0.0049, 5e-4, forceabs=True)) else: raise Exception('Test not defined for this flow') tests.add(CTest('Cl', solver.Cl, 0., 1e-2)) diff --git a/dart/tests/adjoint3.py b/dart/tests/rae_25.py similarity index 62% rename from dart/tests/adjoint3.py rename to dart/tests/rae_25.py index a9bb5a768a233508a06f76d13f278033695c75e6..16473232cd6f7ca7b8c8ef751b82e70efbabf53e 100644 --- a/dart/tests/adjoint3.py +++ b/dart/tests/rae_25.py @@ -2,13 +2,13 @@ # -*- 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. @@ -16,22 +16,20 @@ # limitations under the License. -## Compute adjoint solution of lifting (linear or nonlinear) potential flow around a NACA 0012 +## Compute the nonlinear flow around extruded 2D rae2822 # Adrien Crovato # -# Test the adjoint solver +# Test transonic capabilities and kutta condition # # CAUTION # This test is provided to ensure that the solver works properly. # Mesh refinement may have to be performed to obtain physical results. import numpy as np -import dart.utils as floU import dart.default as floD -import tbox.utils as tboxU import fwk from fwk.testing import * -from fwk.coloring import ccolors +from fwk.coloring import ccolors def main(): # timer @@ -39,31 +37,38 @@ def main(): tms['total'].start() # define flow variables - alpha = 3*np.pi/180 - M_inf = 0.3 - c_ref = 1 - span = 1 + alpha = 1*np.pi/180 + M_inf = 0.72 dim = 3 + # define dimension and mesh size + spn = 0.1 # wing span + lgt = 6.0 # channel length + hgt = 6.0 # channel height + wdt = 3.0 # channel width + c_ref = 1.0 # reference length + S_ref = spn # reference area + fms = 1.0 # farfield mesh size + nms = 0.01 # nearfield mesh size + # mesh the geometry print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) - tms['msh'].start() - pars = {'spn' : span, 'lgt' : 3, 'wdt' : 3, 'hgt' : 3, 'msLe' : 0.02, 'msTe' : 0.02, 'msF' : 1.} - msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip') - # create mesh deformation handler + tms["msh"].start() + pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms} + msh, gmshWriter = floD.mesh(dim, 'models/rae_25.geo', pars, ['field', 'wing', 'symmetry', 'downstream']) morpher = floD.morpher(msh, dim, ['wing']) - tms['msh'].stop() + tms["msh"].stop() # set the problem and add fluid, initial/boundary conditions tms['pre'].start() - pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, span*c_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0.25, 0., 0., 'wing') tms['pre'].stop() # solve problem print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) tms['solver'].start() solver = floD.newton(pbl) - solver.run() + cnvrgd = solver.run() solver.save(gmshWriter) tms['solver'].stop() @@ -83,7 +88,6 @@ def main(): dCdX += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]).dot(np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]])) dClX = np.sqrt(dClX) dCdX = np.sqrt(dCdX) - # recover gradient wrt to AoA from mesh gradient drot = np.array([[-np.sin(alpha), 0, np.cos(alpha)], [0, 0, 0], [-np.cos(alpha), 0, -np.sin(alpha)]]) dClAoA = 0 @@ -95,8 +99,8 @@ def main(): # display results print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) - print(' M alpha Cl Cd dCl_da dCd_da') - print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}'.format(M_inf, alpha*180/np.pi, solver.Cl, solver.Cd, adjoint.dClAoa, adjoint.dCdAoa)) + 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/np.pi, solver.Cl, solver.Cd, solver.Cm)) # display timers tms['total'].stop() @@ -104,25 +108,25 @@ def main(): print('CPU statistics') print(tms) + # visualize solution + floD.initViewer(pbl) + # check results print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) + if (not cnvrgd): + raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET) tests = CTests() - if M_inf == 0.3 and alpha == 3*np.pi/180 and span == 1: - tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 2.5, 5e-2)) # FD 2.526 (step 1e-5) - tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.146, 5e-2)) # FD 0.1460 (step 1e-5) - tests.add(CTest('dCl_dMsh', dClX, 2.016, 5e-2)) - tests.add(CTest('dCd_dMsh', dCdX, 0.133, 5e-2)) - tests.add(CTest('dCl/dAoA (msh)', dClAoA, 2.5, 5e-2)) - tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.146, 5e-2)) - elif M_inf == 0.7 and alpha == 5*np.pi/180 and span == 2: - tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 4.9, 5e-2)) - tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.454, 5e-2)) - tests.add(CTest('dCl_dMsh', dClX, 2.571, 5e-2)) - tests.add(CTest('dCd_dMsh', dCdX, 0.204, 5e-2)) - tests.add(CTest('dCl/dAoA (msh)', dClAoA, 4.9, 5e-2)) - tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.454, 5e-2)) - else: - raise Exception('Test not defined for this flow') + tests.add(CTest('iteration count', solver.nIt, 6, 1, forceabs=True)) + tests.add(CTest('CL', solver.Cl, 0.61, 5e-2)) # 2D, 0.62 + tests.add(CTest('CD', solver.Cd, 0.0016, 1e-2)) # 2D 0.0009 (3D refined 0.0006) + tests.add(CTest('CS', solver.Cs, 0.0000, 1e-3, forceabs=True)) + tests.add(CTest('CM', solver.Cm, -0.116, 5e-2)) # 2D -0.117 + tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 12.7, 5e-2)) # 2D 12.7, FD 12.65 (1e-6) + tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.073, 5e-2)) # 2D 0.051, FD, 0.073 (1e-6) + tests.add(CTest('dCl_dMsh', dClX, 38.0, 5e-2)) + tests.add(CTest('dCd_dMsh', dCdX, 2.0, 5e-2)) + tests.add(CTest('dCl/dAoA (msh)', dClAoA, 12.8, 5e-2)) # 2D 12.9 + tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.073, 5e-2)) # 2D 0.051 tests.run() # eof diff --git a/dart/tests/rae_3.py b/dart/tests/rae_3.py new file mode 100644 index 0000000000000000000000000000000000000000..74a4120e1e80c55e3464698f631a6d02023bcdde --- /dev/null +++ b/dart/tests/rae_3.py @@ -0,0 +1,161 @@ +#!/usr/bin/env python3 +# -*- 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. + + +## Compute the nonlinear flow around rectangular 3D rae2822 +# Adrien Crovato +# +# Test the direct and adjoint solvers, and the mesh morpher in 3D +# +# CAUTION +# This test is provided to ensure that the solver works properly. +# Mesh refinement may have to be performed to obtain physical results. + +import numpy as np +import dart.default as floD +import fwk +from fwk.testing import * +from fwk.coloring import ccolors + +def main(): + # timer + tms = fwk.Timers() + tms['total'].start() + + # define flow variables + alpha = 0 + M_inf = 0.8 + dim = 3 + + # define dimension and mesh size + spn = 1.0 # wing span + lgt = 6.0 # channel length + hgt = 6.0 # channel height + wdt = 3.0 # channel width + c_ref = 1.0 # reference length + S_ref = spn # reference area + fms = 1.0 # farfield mesh size + nms = 0.02 # nearfield mesh size + + # parameters for mesh deformation + alfa = 1*np.pi/180 + xc = 0.25 + dz_max = 0.2*spn + + # mesh the geometry and create mesh deformation handler + print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET) + tms["msh"].start() + pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msN' : nms, 'msF' : fms} + msh, gmshWriter = floD.mesh(dim, 'models/rae_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip') + morpher = floD.morpher(msh, dim, ['wing']) + tms["msh"].stop() + + # set the problem and add fluid, initial/boundary conditions + tms['pre'].start() + pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip') + solver = floD.newton(pbl) # create the solver now so element memory is up to date + tms['pre'].stop() + + # deform the mesh (dummy rotation + bending) + print(ccolors.ANSI_BLUE + 'PyDeforming...' + ccolors.ANSI_RESET) + tms['deform'].start() + morpher.savePos() + rMat = np.matrix([[np.cos(alfa), np.sin(alfa)], [-np.sin(alfa), np.cos(alfa)]]) + vNods = {} + for e in bnd.groups[0].tag.elems: + for n in e.nodes: + if n.no not in vNods: + vNods[n.no] = n + for no in vNods: + nod = np.matrix([[vNods[no].pos[0]-xc], [vNods[no].pos[2]]]) + new = rMat * nod + vNods[no].pos[0] = new[0,0] + vNods[no].pos[2] = new[1,0] + dz_max/(spn*spn)*(vNods[no].pos[1]*vNods[no].pos[1]) + # deform the mesh + morpher.deform() + gmshWriter.save(msh.name+"_def") + tms['deform'].stop() + + # solve problem + print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET) + tms['solver'].start() + cnvrgd = solver.run() + solver.save(gmshWriter) + tms['solver'].stop() + + # solve adjoint problem + print(ccolors.ANSI_BLUE + 'PyGradient...' + ccolors.ANSI_RESET) + tms['adjoint'].start() + adjoint = floD.adjoint(solver, morpher) + adjoint.run() + adjoint.save(gmshWriter) + tms['adjoint'].stop() + + # compute norm of gradient wrt to mesh + dClX = 0 + dCdX = 0 + for n in msh.nodes: + dClX += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1], adjoint.dClMsh[n.row][2]]).dot(np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1], adjoint.dClMsh[n.row][2]])) + dCdX += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]).dot(np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]])) + dClX = np.sqrt(dClX) + dCdX = np.sqrt(dCdX) + # recover gradient wrt to AoA from mesh gradient + drot = np.array([[-np.sin(alfa), 0, np.cos(alfa)], [0, 0, 0], [-np.cos(alfa), 0, -np.sin(alfa)]]) + dClAoA = 0 + dCdAoA = 0 + for n in bnd.nodes: + dx = drot.dot(np.array([n.pos[0] - 0.25, n.pos[1], n.pos[2]])) + dClAoA += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1], adjoint.dClMsh[n.row][2]]).dot(dx) + dCdAoA += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]).dot(dx) + + # display results + print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET) + print(' M alpha Cl Cd dCl_da dCd_da') + print('{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}'.format(M_inf, alfa*180/np.pi, solver.Cl, solver.Cd, adjoint.dClAoa, adjoint.dCdAoa)) + + # display timers + tms['total'].stop() + print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET) + print('CPU statistics') + print(tms) + + # visualize solution + floD.initViewer(pbl) + + # check results + print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET) + if (not cnvrgd): + raise Exception(ccolors.ANSI_RED + 'Flow solver failed to converge!' + ccolors.ANSI_RESET) + tests = CTests() + tests.add(CTest('iteration count', solver.nIt, 7, 1, forceabs=True)) + tests.add(CTest('CL', solver.Cl, 0.21, 5e-2)) + tests.add(CTest('CD', solver.Cd, 0.008, 1e-2)) + tests.add(CTest('CS', solver.Cs, -0.020, 5e-2)) + tests.add(CTest('CM', solver.Cm, -0.080, 1e-2)) + tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 3.0, 5e-2)) # FD 3.0144 (1e-6) + tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.2, 5e-2)) # FD 0.2007 (1e-6) + tests.add(CTest('dCl_dMsh', dClX, 2.3, 5e-2)) + tests.add(CTest('dCd_dMsh', dCdX, 0.2, 5e-2)) + tests.add(CTest('dCl/dAoA (msh)', dClAoA, 3.0, 5e-2)) + tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.2, 5e-2)) + tests.run() + + # eof + print('') + +if __name__ == "__main__": + main() diff --git a/dart/validation/agard.py b/dart/validation/agard.py index 3e048c0585be28f3cc6d914508771f2d17e46ca4..86d5b9b767a88ab9107adba31809b85365a4d772 100644 --- a/dart/validation/agard.py +++ b/dart/validation/agard.py @@ -69,7 +69,7 @@ def main(): # set the problem and add fluid, initial/boundary conditions tms['pre'].start() - pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') + pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip') tms['pre'].stop() # solve problem diff --git a/dart/validation/onera.py b/dart/validation/onera.py index 5739b25dceff469f07aecee6a60399e5ab805ad3..6957094a6f26428f4bc3814c8ca757320c9cc1eb 100644 --- a/dart/validation/onera.py +++ b/dart/validation/onera.py @@ -67,7 +67,7 @@ def main(): # set the problem and add fluid, initial/boundary conditions tms['pre'].start() - pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip') + pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip') tms['pre'].stop() # solve problem diff --git a/ext/amfe b/ext/amfe index 87212cfb034a92cc83acfeea751f076a84c12779..054ba00c327e4f2745ff9c559832247f90c3cec9 160000 --- a/ext/amfe +++ b/ext/amfe @@ -1 +1 @@ -Subproject commit 87212cfb034a92cc83acfeea751f076a84c12779 +Subproject commit 054ba00c327e4f2745ff9c559832247f90c3cec9