Skip to content
Snippets Groups Projects
Commit 016ea870 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

2d and 3d boundary layer examples

parent 8251eabe
No related branches found
No related tags found
No related merge requests found
import gmsh
import sys
import math
import numpy as np
gmsh.initialize(sys.argv)
gmsh.model.add("NACA 0012")
# create the boundary layer by extrusion, or as a mesh constraint?
by_extrusion = False
# rounded trailing edge?
rounded = True
# base mesh sizes
fact = 1;
lc1 = 0.01 * fact
lc2 = 0.3 * fact
# generate curved mesh?
order2 = True
# xy coordinates of top part of NACA 0012 profile
naca = [0.9987518, 0.0014399, 0.9976658, 0.0015870, 0.9947532, 0.0019938,
0.9906850, 0.0025595, 0.9854709, 0.0032804, 0.9791229, 0.0041519,
0.9716559, 0.0051685, 0.9630873, 0.0063238, 0.9534372, 0.0076108,
0.9427280, 0.0090217, 0.9309849, 0.0105485, 0.9182351, 0.0121823,
0.9045085, 0.0139143, 0.8898372, 0.0157351, 0.8742554, 0.0176353,
0.8577995, 0.0196051, 0.8405079, 0.0216347, 0.8224211, 0.0237142,
0.8035813, 0.0258337, 0.7840324, 0.0279828, 0.7638202, 0.0301515,
0.7429917, 0.0323294, 0.7215958, 0.0345058, 0.6996823, 0.0366700,
0.6773025, 0.0388109, 0.6545085, 0.0409174, 0.6313537, 0.0429778,
0.6078921, 0.0449802, 0.5841786, 0.0469124, 0.5602683, 0.0487619,
0.5362174, 0.0505161, 0.5120819, 0.0521620, 0.4879181, 0.0536866,
0.4637826, 0.0550769, 0.4397317, 0.0563200, 0.4158215, 0.0574033,
0.3921079, 0.0583145, 0.3686463, 0.0590419, 0.3454915, 0.0595747,
0.3226976, 0.0599028, 0.3003177, 0.0600172, 0.2784042, 0.0599102,
0.2570083, 0.0595755, 0.2361799, 0.0590081, 0.2159676, 0.0582048,
0.1964187, 0.0571640, 0.1775789, 0.0558856, 0.1594921, 0.0543715,
0.1422005, 0.0526251, 0.1257446, 0.0506513, 0.1101628, 0.0484567,
0.0954915, 0.0460489, 0.0817649, 0.0434371, 0.0690152, 0.0406310,
0.0572720, 0.0376414, 0.0465628, 0.0344792, 0.0369127, 0.0311559,
0.0283441, 0.0276827, 0.0208771, 0.0240706, 0.0145291, 0.0203300,
0.0093149, 0.0164706, 0.0052468, 0.0125011, 0.0023342, 0.0084289,
0.0005839, 0.0042603, 0.0000000, 0.0000000]
# create profile from points
pts = []
l = len(naca)
for i in range(0, l // 2):
pts.append(gmsh.model.occ.addPoint(naca[2 * i], naca[2 * i + 1], 0, lc1))
for i in range(l // 2 - 2, -1, -1):
pts.append(gmsh.model.occ.addPoint(naca[2 * i], -naca[2 * i + 1], 0, lc1))
pts.reverse()
curv = []
curv.append(gmsh.model.occ.addSpline(pts))
if rounded:
# circle as trailing edge
c = gmsh.model.occ.addPoint(0.9985510, 0.0000000, 0, lc1)
curv.append(gmsh.model.occ.addCircleArc(pts[-1], c, pts[0]))
else:
pt = gmsh.model.occ.addPoint(1.0095, 0.0000000, 0, lc1)
curv.append(gmsh.model.occ.addLine(pts[-1], pt))
curv.append(gmsh.model.occ.addLine(pt, pts[0]))
cl = gmsh.model.occ.addCurveLoop(curv)
if by_extrusion:
gmsh.model.occ.synchronize()
# a boundary layer can be created through extrusion using the built-in CAD
# kernel: this creates topological entities that will be filled with a
# discrete geometry (a mesh) during mesh generation
n = np.linspace(1, 1, 7)
d = np.logspace(-4, -2, 7)
extbl = gmsh.model.geo.extrudeBoundaryLayer(gmsh.model.getEntities(1),
n, -d, True)
# create curve loop with "top" curves of the boundary layer
cl2 = gmsh.model.geo.addCurveLoop([c[1] for c in extbl[::2]])
# connect it with the exterior surface
p1 = gmsh.model.geo.addPoint(-1, -1, 0, lc2)
p2 = gmsh.model.geo.addPoint(2, -1, 0, lc2)
p3 = gmsh.model.geo.addPoint(2, 1, 0, lc2)
p4 = gmsh.model.geo.addPoint(-1, 1, 0, lc2)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p1)
cl3 = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
s2 = gmsh.model.geo.addPlaneSurface([cl3, cl2])
gmsh.model.geo.synchronize()
else:
# in 2D, boundary layers can also be specified as a meshing constraint,
# trough the BoundaryLayer field; this is quite a bit more general, as it
# handles intersections between several boundary layers, fans, etc.
s = gmsh.model.occ.addPlaneSurface([cl])
r = gmsh.model.occ.addRectangle(-1, -1, 0, 3, 2)
gmsh.model.occ.cut([(2, r)], [(2, s)])
gmsh.model.occ.synchronize()
f = gmsh.model.mesh.field.add('BoundaryLayer')
gmsh.model.mesh.field.setNumbers(f, 'CurvesList', curv)
gmsh.model.mesh.field.setNumber(f, 'Size', lc1 / 100)
gmsh.model.mesh.field.setNumber(f, 'Ratio', 2.4)
gmsh.model.mesh.field.setNumber(f, 'Quads', 1)
gmsh.model.mesh.field.setNumber(f, 'Thickness', 0.02)
if not rounded:
# create a fan at the trailing edge
gmsh.option.setNumber('Mesh.BoundaryLayerFanElements', 7)
gmsh.model.mesh.field.setNumbers(f, 'FanPointsList', [pt])
gmsh.model.mesh.field.setAsBoundaryLayer(f)
# generate the mesh
gmsh.model.mesh.generate(2)
# 2nd order + fast curving of the boundary layer + optimization
if order2:
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize('HighOrderFastCurving')
gmsh.model.mesh.optimize('HighOrder')
gmsh.write('naca_boundary_layer_2d.msh')
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
gmsh.finalize()
import gmsh
import sys
import math
import numpy as np
gmsh.initialize(sys.argv)
gmsh.model.add("NACA 0012 with a round tip")
# extrusion length of the profle in the z-direction
z = 0.5
# base mesh sizes
fact = 1;
lc1 = 0.01 * fact
lc2 = 0.3 * fact
gmsh.option.setNumber('Mesh.Algorithm3D', 10)
# generate curved mesh?
order2 = False
# xy coordinates of top part of NACA 0012 profile
naca = [0.9987518, 0.0014399, 0.9976658, 0.0015870, 0.9947532, 0.0019938,
0.9906850, 0.0025595, 0.9854709, 0.0032804, 0.9791229, 0.0041519,
0.9716559, 0.0051685, 0.9630873, 0.0063238, 0.9534372, 0.0076108,
0.9427280, 0.0090217, 0.9309849, 0.0105485, 0.9182351, 0.0121823,
0.9045085, 0.0139143, 0.8898372, 0.0157351, 0.8742554, 0.0176353,
0.8577995, 0.0196051, 0.8405079, 0.0216347, 0.8224211, 0.0237142,
0.8035813, 0.0258337, 0.7840324, 0.0279828, 0.7638202, 0.0301515,
0.7429917, 0.0323294, 0.7215958, 0.0345058, 0.6996823, 0.0366700,
0.6773025, 0.0388109, 0.6545085, 0.0409174, 0.6313537, 0.0429778,
0.6078921, 0.0449802, 0.5841786, 0.0469124, 0.5602683, 0.0487619,
0.5362174, 0.0505161, 0.5120819, 0.0521620, 0.4879181, 0.0536866,
0.4637826, 0.0550769, 0.4397317, 0.0563200, 0.4158215, 0.0574033,
0.3921079, 0.0583145, 0.3686463, 0.0590419, 0.3454915, 0.0595747,
0.3226976, 0.0599028, 0.3003177, 0.0600172, 0.2784042, 0.0599102,
0.2570083, 0.0595755, 0.2361799, 0.0590081, 0.2159676, 0.0582048,
0.1964187, 0.0571640, 0.1775789, 0.0558856, 0.1594921, 0.0543715,
0.1422005, 0.0526251, 0.1257446, 0.0506513, 0.1101628, 0.0484567,
0.0954915, 0.0460489, 0.0817649, 0.0434371, 0.0690152, 0.0406310,
0.0572720, 0.0376414, 0.0465628, 0.0344792, 0.0369127, 0.0311559,
0.0283441, 0.0276827, 0.0208771, 0.0240706, 0.0145291, 0.0203300,
0.0093149, 0.0164706, 0.0052468, 0.0125011, 0.0023342, 0.0084289,
0.0005839, 0.0042603, 0.0000000, 0.0000000]
# create profile in z=0 plane from points
pts = []
l = len(naca)
for i in range(0, l // 2):
pts.append(gmsh.model.occ.addPoint(naca[2 * i], naca[2 * i + 1], 0))
for i in range(l // 2 - 2, -1, -1):
pts.append(gmsh.model.occ.addPoint(naca[2 * i], -naca[2 * i + 1], 0))
pts.reverse()
spl = gmsh.model.occ.addSpline(pts)
# circle as trailing edge
c = gmsh.model.occ.addPoint(0.9985510, 0.0000000, 0)
cir = gmsh.model.occ.addCircleArc(pts[-1], c, pts[0])
# extrude the profile along z
ext = gmsh.model.occ.extrude([(1, spl), (1, cir)], 0, 0, z)
# cut the extruded profile at the leading and trailing edge to prepare for the
# rounded wing tip
bb = gmsh.model.occ.getBoundingBox(1, cir)
p1 = gmsh.model.occ.addPoint(0, 0, 0)
p2 = gmsh.model.occ.addPoint(0, 0, z)
c1 = gmsh.model.occ.addLine(p1, p2)
p3 = gmsh.model.occ.addPoint(bb[3], 0, 0)
p4 = gmsh.model.occ.addPoint(bb[3], 0, z)
c2 = gmsh.model.occ.addLine(p3, p4)
gmsh.model.occ.fragment([(1, c1), (1, c2)], gmsh.model.occ.getEntities(2))
# retrieve tip curves that will be revolved
eps = 1e-6
gmsh.option.setNumber('Geometry.OCCBoundsUseStl', 1)
tc = gmsh.model.occ.getEntitiesInBoundingBox(-eps,-eps,z-eps,
bb[3]+eps,1,z+eps, dim=1)
rev = gmsh.model.occ.revolve(tc, 0,0,z, 1,0,0, math.pi/2)
gmsh.model.occ.revolve(rev[0::4], 0,0,z, 1,0,0, math.pi/2)
# glue surfaces, set mesh size and sync to model
gmsh.model.occ.fragment(gmsh.model.occ.getEntities(2), [])
gmsh.model.occ.mesh.setSize(gmsh.model.occ.getEntities(0), lc1)
gmsh.model.occ.synchronize()
# create a boundary layer for all the surfaces through extrusion using the
# built-in CAD kernel: this creates topological entities that will be filled
# with a discrete geometry (a mesh) during mesh generation; more general
# boundary layer meshing constraints are currently only available in 2D, through
# the BoundaryLayer Field - see 'naca_boundary_layer_2d.py'.
n = np.linspace(1, 1, 7)
d = np.logspace(-4, -2, 7)
extbl = gmsh.model.geo.extrudeBoundaryLayer(gmsh.model.getEntities(2),
n, d, True)
# get "top" surfaces of the boundary layer
top = []
for i in range(1, len(extbl)):
if extbl[i][0] == 3:
top.append(extbl[i-1])
# extract the boundary of the top surfaces, which define a hole in the symmetry
# plane, and create the symmetry plane
gmsh.model.geo.synchronize()
bnd = gmsh.model.getBoundary(top)
cl2 = gmsh.model.geo.addCurveLoop([c[1] for c in bnd])
p1 = gmsh.model.geo.addPoint(-1, -1, 0, lc2)
p2 = gmsh.model.geo.addPoint(2, -1, 0, lc2)
p3 = gmsh.model.geo.addPoint(2, 1, 0, lc2)
p4 = gmsh.model.geo.addPoint(-1, 1, 0, lc2)
l1 = gmsh.model.geo.addLine(p1, p2)
l2 = gmsh.model.geo.addLine(p2, p3)
l3 = gmsh.model.geo.addLine(p3, p4)
l4 = gmsh.model.geo.addLine(p4, p1)
cl3 = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
s2 = gmsh.model.geo.addPlaneSurface([cl3, cl2])
# create 3D air box
p11 = gmsh.model.geo.addPoint(-1, -1, 2*z, lc2)
p12 = gmsh.model.geo.addPoint(2, -1, 2*z, lc2)
p13 = gmsh.model.geo.addPoint(2, 1, 2*z, lc2)
p14 = gmsh.model.geo.addPoint(-1, 1, 2*z, lc2)
l11 = gmsh.model.geo.addLine(p11, p12)
l12 = gmsh.model.geo.addLine(p12, p13)
l13 = gmsh.model.geo.addLine(p13, p14)
l14 = gmsh.model.geo.addLine(p14, p11)
l_1_11 = gmsh.model.geo.addLine(p1, p11)
l_2_12 = gmsh.model.geo.addLine(p2, p12)
l_3_13 = gmsh.model.geo.addLine(p3, p13)
l_4_14 = gmsh.model.geo.addLine(p4, p14)
cl3 = gmsh.model.geo.addCurveLoop([l11, l12, l13, l14])
s3 = gmsh.model.geo.addPlaneSurface([cl3])
cl4 = gmsh.model.geo.addCurveLoop([l1, l_2_12, -l11, -l_1_11])
s4 = gmsh.model.geo.addPlaneSurface([cl4])
cl5 = gmsh.model.geo.addCurveLoop([l2, l_3_13, -l12, -l_2_12])
s5 = gmsh.model.geo.addPlaneSurface([cl5])
cl6 = gmsh.model.geo.addCurveLoop([l3, l_4_14, -l13, -l_3_13])
s6 = gmsh.model.geo.addPlaneSurface([cl6])
cl7 = gmsh.model.geo.addCurveLoop([l4, l_1_11, -l14, -l_4_14])
s7 = gmsh.model.geo.addPlaneSurface([cl7])
b = [s[1] for s in top]
b.extend([s2, s3, s4, s5, s6, s7])
sl = gmsh.model.geo.addSurfaceLoop(b)
v = gmsh.model.geo.addVolume([sl])
gmsh.model.geo.synchronize()
# generate the mesh
gmsh.model.mesh.generate(3)
# 2nd order + fast curving of the boundary layer + optimization
if order2:
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize('HighOrderFastCurving')
gmsh.model.mesh.optimize('HighOrder')
gmsh.write('naca_boundary_layer_3d.msh')
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
gmsh.finalize()
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment