gmshで二次元翼型解析のためのメッシュを作成する
作成するメッシュは以下の2種類
- 構造格子
- 非構造格子+境界層メッシュ
また，mshファイルはpythonを使って作成する
解析に使ったファイルはGithubにアップしてある．≫https://github.com/mtkbirdman/mtkbirdman.com/tree/master/OpenFOAM/2D_airfoil_ver2
はじめに
OpenFOAMの2次元翼型解析のチュートリアル(≫OpenFOAM_ User Guide_ Turbulent flow over NACA0012 airfoil (2D))を参考に，2次元翼型用のメッシュを作成してみる
イメージはこんな感じ
そして実際に完成したメッシュがこれ
実行環境
実行環境を以下に示す
| OS | Windows10 Home 64bit |
| CPU | Intel Core i5-8250U CPU@ 1.60GHz 1.80GHz |
| Ubuntu | 20.04 LTS |
| python | 3.8.0 |
| gmsh | 4.6.0 |
参考
gmshのインストールについてはこちら
gmshの基本的な使い方についてはこちら
≫自作ソルバにGmshを！ - Qiita
≫Gmsh の使い方 - PENGUINITIS
≫Gmsh reference manual (stable release)
pythonのモジュールであるgmshについてのドキュメントはこちら
≫the Gmsh Python API (v4.6)
≫the Gmsh Python tutorials
↓Pythonのインストールについてはこちら
それでは実際にメッシュを作成していく
全体の流れ
最初にプログラム全体のフローチャートとソースコードを示す
フローチャートの点線部分が構造格子と非構造格子で違う部分で，それ以外のプログラムは共通である
プログラムにおいてはstructure=True/Falseで切り替えができるようにしてある
def genMesh(airfoilFile,structure=True):
#read airfoil file
ar = np.loadtxt(airfoilFile, skiprows=1)
#removing duplicate end point
if np.max(np.abs(ar[0] - ar[(ar.shape[0]-1)]))<1e-6:
ar = ar[:-1]
#calculate TE angle
TE_angle_U=math.atan((ar[1][1]-ar[0][1])/(ar[1][0]-ar[0][0]))
TE_angle_D=math.atan((ar[ar.shape[0]-1][1]-ar[0][1])/(ar[ar.shape[0]-1][0]-ar[0][0]))
#initialize gmsh and add model
gmsh.initialize()
#gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add("airfoil")
#set point index
TE_pointIndex = 1000
LE_pointIndex = TE_pointIndex+np.argmin(ar,axis=0)[0]
MAX_pointIndex = TE_pointIndex+ar.shape[0]-1
#add airfoil points
pointIndex = TE_pointIndex
for n in range(ar.shape[0]):
lc = 1e-2 if n != np.argmin(ar,axis=0)[0] else 1e-3
gmsh.model.geo.addPoint(ar[n][0],ar[n][1],0.,lc,pointIndex)
pointIndex += 1
#add poiuts
lc = 2.
R=50.
gmsh.model.geo.addPoint(1.+R*math.cos(math.pi/2+TE_angle_U),R*math.sin(math.pi/2+TE_angle_U),0.,lc,1)
gmsh.model.geo.addPoint(1.-R,0.,0.,lc,2)
gmsh.model.geo.addPoint(1.+R*math.cos(math.pi/2+TE_angle_D),-R*math.sin(math.pi/2+TE_angle_D),0.,lc,3)
gmsh.model.geo.addPoint(R,-R*math.sin(math.pi/2+TE_angle_D),0.,lc,4)
gmsh.model.geo.addPoint(R,0.,0.,lc,5)
gmsh.model.geo.addPoint(R,R*math.sin(math.pi/2+TE_angle_U),0.,lc,6)
#add Circle and Line
gmsh.model.geo.addCircleArc(1,TE_pointIndex,2,1)
gmsh.model.geo.addCircleArc(2,TE_pointIndex,3,2)
gmsh.model.geo.addLine(3,4,3)
gmsh.model.geo.addLine(4,5,4)
gmsh.model.geo.addLine(5,6,5)
gmsh.model.geo.addLine(6,1,6)
gmsh.model.geo.addLine(TE_pointIndex,1,7)
gmsh.model.geo.addLine(LE_pointIndex,2,8)
gmsh.model.geo.addLine(TE_pointIndex,3,9)
gmsh.model.geo.addLine(TE_pointIndex,5,10)
#add airfoil spline (upper surface and under surface)
list_index=[pointIndex for pointIndex in range(TE_pointIndex,LE_pointIndex+1)]
gmsh.model.geo.addSpline(list_index,11)
list_index=[pointIndex for pointIndex in range(LE_pointIndex,MAX_pointIndex+1)]
list_index.append(TE_pointIndex)
gmsh.model.geo.addSpline(list_index,12)
if structure:
#add CurveLoop and PlenSurface
gmsh.model.geo.addCurveLoop([-11,7,1,-8], 1)
gmsh.model.geo.addCurveLoop([8,2,-9,-12], 2)
gmsh.model.geo.addCurveLoop([9,3,4,-10], 3)
gmsh.model.geo.addCurveLoop([10,5,6,-7], 4)
gmsh.model.geo.addPlaneSurface([1],1)
gmsh.model.geo.addPlaneSurface([2],2)
gmsh.model.geo.addPlaneSurface([3],3)
gmsh.model.geo.addPlaneSurface([4],4)
#set TransfiniteCurve
NH=128;NW=128;NA=128
list_index=[1,2]
for CurveIndex in [1,2]:
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NA)
prog=1.05
for (CurveIndex,direction) in zip([4,5,7,9],[False,True,True,True]):
prog_tmp = prog if direction else -prog
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NH,meshType="Progression",coef=prog_tmp)
prog=1.04
CurveIndex=10
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NW,meshType="Progression",coef=prog_tmp)
for CurveIndex in [3,6]:
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NW)
prog=1.07
CurveIndex=8
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NW,meshType="Progression",coef=prog)
prog=1.03
for (CurveIndex,direction) in zip([11,12],[False,True]):
prog_tmp = prog if direction else -prog
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NA,meshType="Progression",coef=prog_tmp)
for SurfaceIndex in [1,2,3,4]:
gmsh.model.geo.mesh.setTransfiniteSurface(SurfaceIndex)
gmsh.model.geo.mesh.setRecombine(2,SurfaceIndex)
#extrude
gmsh.model.geo.extrude([(2,SurfaceIndex)],0.,0.,1.,[1],[1.],recombine=True)
#add PhysicalGroup and set Name
gmsh.model.addPhysicalGroup(2,[1,2,3,4],1) ;gmsh.model.setPhysicalName(2,1,"front")
gmsh.model.addPhysicalGroup(2,[34,56,78,100],2) ;gmsh.model.setPhysicalName(2,2,"back")
gmsh.model.addPhysicalGroup(2,[29,47],3) ;gmsh.model.setPhysicalName(2,3,"inlet")
gmsh.model.addPhysicalGroup(2,[73,91],4) ;gmsh.model.setPhysicalName(2,4,"exit")
gmsh.model.addPhysicalGroup(2,[95],5) ;gmsh.model.setPhysicalName(2,5,"top")
gmsh.model.addPhysicalGroup(2,[69],6) ;gmsh.model.setPhysicalName(2,6,"bottom")
gmsh.model.addPhysicalGroup(2,[21,55],7) ;gmsh.model.setPhysicalName(2,7,"aerofoil")
gmsh.model.addPhysicalGroup(3,[1,2,3,4],8) ;gmsh.model.setPhysicalName(3,8,"internal")
else:
#add CurveLoop and PlenSurface
gmsh.model.geo.addCurveLoop([1,2,3,4,5,6], 1)
gmsh.model.geo.addCurveLoop([11,12], 2)
gmsh.model.geo.addPlaneSurface([1,2],1)
#set BoundaryLayer field
gmsh.model.mesh.field.add("BoundaryLayer",1)
gmsh.model.mesh.field.setNumbers(1,"EdgesList",[11,12])
gmsh.model.mesh.field.setNumber(1,"Quads",1)
gmsh.model.mesh.field.setNumber(1,"hwall_n",1e-3)
gmsh.model.mesh.field.setNumber(1,"thickness",1e-2)
gmsh.model.mesh.field.setAsBoundaryLayer(1)
#extrude
gmsh.model.geo.extrude([(2,1)],0.,0.,1.,[1],[1.],recombine=True)
#add PhysicalGroup and set Name
gmsh.model.addPhysicalGroup(2,[1],1) ;gmsh.model.setPhysicalName(2,1,"front")
gmsh.model.addPhysicalGroup(2,[54],2) ;gmsh.model.setPhysicalName(2,2,"back")
gmsh.model.addPhysicalGroup(2,[25,29],3);gmsh.model.setPhysicalName(2,3,"inlet")
gmsh.model.addPhysicalGroup(2,[37,41],4);gmsh.model.setPhysicalName(2,4,"exit")
gmsh.model.addPhysicalGroup(2,[45],5) ;gmsh.model.setPhysicalName(2,5,"top")
gmsh.model.addPhysicalGroup(2,[33],6) ;gmsh.model.setPhysicalName(2,6,"bottom")
gmsh.model.addPhysicalGroup(2,[49,53],7);gmsh.model.setPhysicalName(2,7,"aerofoil")
gmsh.model.addPhysicalGroup(3,[1],8) ;gmsh.model.setPhysicalName(3,8,"internal")
#generate mesh and finalize gmsh
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(3)
gmsh.write('airfoil.msh')
gmsh.finalize()
print("gmsh > done!")
if os.system("gmshToFoam airfoil.msh > /dev/null") != 0:
print("error during conversion to OpenFoam mesh!")
return(-1)
print("gmshToFOAM > done!")
#set boundary
with open("constant/polyMesh/boundary", "rt") as inFile:
with open("constant/polyMesh/boundaryTemp", "wt") as outFile:
inBlock = False
inAerofoil = False
for line in inFile:
if "front" in line or "back" in line:
inBlock = True
elif "aerofoil" in line:
inAerofoil = True
if inBlock and "type" in line:
line = line.replace("patch", "empty")
inBlock = False
if inAerofoil and "type" in line:
line = line.replace("patch", "wall")
inAerofoil = False
outFile.write(line)
os.rename("constant/polyMesh/boundaryTemp","constant/polyMesh/boundary")
return(0)1つずつ説明していく
翼型座標の読み込み~後縁角の計算
該当部分のコードを以下に示す
後縁に2つ点がある場合は1つに減らし，後縁角はコードラインを挟んで上側と下側に分けて計算している
#read airfoil file
ar = np.loadtxt(airfoilFile, skiprows=1)
#removing duplicate end point
if np.max(np.abs(ar[0] - ar[(ar.shape[0]-1)]))<1e-6:
ar = ar[:-1]
#calculate TE angle
TE_angle_U=math.atan((ar[1][1]-ar[0][1])/(ar[1][0]-ar[0][0]))
TE_angle_D=math.atan((ar[ar.shape[0]-1][1]-ar[0][1])/(ar[ar.shape[0]-1][0]-ar[0][0]))initialize~モデルの追加
gmshを始めるときは，gmsh.initialize()を行い，モデルを追加する必要がある
≫tutorial_python_t1.py - GitLab
#initialize gmsh and add model
gmsh.initialize()
#gmsh.option.setNumber("General.Terminal", 1)
gmsh.model.add("airfoil")gmsh.option.setNumber("General.Terminal", 1)を加えると，gmshが実行されているときのログを画面に出力することができる．デバッグのときなどはエラーを見つけるのに役に立つ
Pointの追加
gmsh.model.geo.addPointメソッドを使ってPointを追加していく．
使い方はこんな感じ
gmsh.model.geo.addPoint(x座標, y座標, z座標, メッシュの大きさ, Pointのタグ)x座標，y座標，z座標，メッシュの大きさの単位はすべて[m]で，メッシュの大きさはその点に隣接するメッシュのみに適用される
実際に加えたPointとソースコードを以下に示す
#set point index
TE_pointIndex = 1000
LE_pointIndex = TE_pointIndex+np.argmin(ar,axis=0)[0]
MAX_pointIndex = TE_pointIndex+ar.shape[0]-1
#add airfoil points
pointIndex = TE_pointIndex
for n in range(ar.shape[0]):
lc = 1e-2 if n != np.argmin(ar,axis=0)[0] else 1e-3
gmsh.model.geo.addPoint(ar[n][0],ar[n][1],0.,lc,pointIndex)
pointIndex += 1
#add poiuts
lc = 2.
R=50.
gmsh.model.geo.addPoint(1.+R*math.cos(math.pi/2+TE_angle_U),R*math.sin(math.pi/2+TE_angle_U),0.,lc,1)
gmsh.model.geo.addPoint(1.-R,0.,0.,lc,2)
gmsh.model.geo.addPoint(1.+R*math.cos(math.pi/2+TE_angle_D),-R*math.sin(math.pi/2+TE_angle_D),0.,lc,3)
gmsh.model.geo.addPoint(R,-R*math.sin(math.pi/2+TE_angle_D),0.,lc,4)
gmsh.model.geo.addPoint(R,0.,0.,lc,5)
gmsh.model.geo.addPoint(R,R*math.sin(math.pi/2+TE_angle_U),0.,lc,6)np.argminを使って，x座標が最小の点のインデックスを取得し，前縁点のpointIndexを計算している
理由は後述するが，翼型の点を追加するときに前縁点のメッシュサイズのみ1e-2に設定してある
Point1~3は後縁点(1，0)を中心に半径Rとなる円周上に配置しており，原点とPoint1，Point3をつなぐ線はそれぞれ翼後縁と直角になるようにしている
Curveの追加
先ほどと追加した点をつないで線を追加する
使うメソッドは以下の通り
#円弧を追加
gmsh.model.geo.addCircleArc(始点のタグ,原点のタグ,終点のタグ,Circleのタグ)
#直線を追加
gmsh.model.geo.addLine(始点のタグ,終点のタグ,Lineのタグ)
#スプラインを追加
gmsh.model.geo.addSpline([点のタグのlist],Splineのタグ)円弧は視点から終点に向かって反時計回りに描かれる
実際に加えたCurveとソースコードを以下に示す
※Nurb11とNurb12が重なっている
≫NURBSとは？ - Qiita
#add Circle and Line
gmsh.model.geo.addCircleArc(1,TE_pointIndex,2,1)
gmsh.model.geo.addCircleArc(2,TE_pointIndex,3,2)
gmsh.model.geo.addLine(3,4,3)
gmsh.model.geo.addLine(4,5,4)
gmsh.model.geo.addLine(5,6,5)
gmsh.model.geo.addLine(6,1,6)
gmsh.model.geo.addLine(TE_pointIndex,1,7)
gmsh.model.geo.addLine(LE_pointIndex,2,8)
gmsh.model.geo.addLine(TE_pointIndex,3,9)
gmsh.model.geo.addLine(TE_pointIndex,5,10)
#add airfoil spline (upper surface and under surface)
list_index=[pointIndex for pointIndex in range(TE_pointIndex,LE_pointIndex+1)]
gmsh.model.geo.addSpline(list_index,11)
list_index=[pointIndex for pointIndex in range(LE_pointIndex,MAX_pointIndex+1)]
list_index.append(TE_pointIndex)
gmsh.model.geo.addSpline(list_index,12)Line7～10は領域を4つに分けるための直線で，構造格子を作るときに使用する．(非構造格子では使わない)
翼型のスプラインは上面と下面に分ける．(理由は後述)
翼下面のスプラインを構成する点は，前縁点~後縁の1つ前の点~後縁点であることに注意する
Surfaceの追加
さきほど追加した線を組み合わせて閉曲線を追加し，その閉曲線をもとに平面を追加する
#閉曲線を追加
gmsh.model.geo.addCurveLoop([曲線のタグのlist],CurveLoopのタグ)
#平面を追加
gmsh.model.geo.addPlaneSurface([閉曲線のタグのlist],PlaneSurfaceのタグ)曲線には始点→終点という向きがあり，閉曲線を追加するときにはすべての線の向きがきちんと1周できるようにそろっている必要がある．曲線の向きを逆にしたいときは，タグに"-"をつければいい
また，addPlaneSurfaceメソッドの閉曲線のlistは，[一番外側の閉曲線，くり抜く閉曲線1，くり抜く閉曲線，・・・]と指定する
ソースコードを以下に示す
if structure:
#add CurveLoop and PlenSurface
gmsh.model.geo.addCurveLoop([-11,7,1,-8], 1)
gmsh.model.geo.addCurveLoop([8,2,-9,-12], 2)
gmsh.model.geo.addCurveLoop([9,3,4,-10], 3)
gmsh.model.geo.addCurveLoop([10,5,6,-7], 4)
gmsh.model.geo.addPlaneSurface([1],1)
gmsh.model.geo.addPlaneSurface([2],2)
gmsh.model.geo.addPlaneSurface([3],3)
gmsh.model.geo.addPlaneSurface([4],4)
else:
#add CurveLoop and PlenSurface
gmsh.model.geo.addCurveLoop([1,2,3,4,5,6], 1)
gmsh.model.geo.addCurveLoop([11,12], 2)
gmsh.model.geo.addPlaneSurface([1,2],1)構造格子の場合は領域を4つの平面に分けているのに対して，非構造格子では外側の閉曲線から翼型の閉曲線をくり抜いていて1つの平面にしている
メッシュの定義
構造格子と非構造格子でかなり違うので，それぞれ説明していく
非構造格子
gmshでは，デフォルトで非構造格子が選択される
今回はメッシュサイズの変更と境界層メッシュの導入を行う
メッシュサイズの調整
先ほど説明したように，点を加えるときの引数によってその点に接するメッシュのサイズを指定することができ，点と点の間のメッシュサイズはいい感じに決定される
曲線を追加したときの曲線上のメッシュサイズは，曲線の始点と終点で指定したメッシュサイズをもとに良い感じに決定される．
下の画像は一辺10 [m]の正方形で，4つの点のメッシュサイズを10，1，0.1，0.01と指定したものである
翼型のスプライン曲線を追加したときに翼の上面と下面でスプライン曲線を分けて前縁点のみメッシュサイズを小さく指定したのは，翼の後縁から前縁にかけて細かくなるようなメッシュを指定するためである
境界層メッシュの定義
境界層メッシュはmesh size fieldsを使って定義できる
≫[Gmsh] Boundary Layer Field
≫Gmsh reference manual (stable release) - 6.3.1 Specifying mesh element sizes
≫tutorial_python_t10.py - GitLab
addメソッドで"BoundaryLayer"を指定してfieldを追加し，setNumbersメソッドやsetNumberメソッドでオプションの値を指定する
最後にsetAsBoundaryLayerメソッドで，追加したfieldをBoundaryLayerとして設定する
if structure:
(略)
else:
#set BoundaryLayer field
gmsh.model.mesh.field.add("BoundaryLayer",1)
gmsh.model.mesh.field.setNumbers(1,"EdgesList",[11,12])
gmsh.model.mesh.field.setNumber(1,"Quads",1)
gmsh.model.mesh.field.setNumber(1,"hwall_n",1e-3)
gmsh.model.mesh.field.setNumber(1,"thickness",1e-2)
gmsh.model.mesh.field.setAsBoundaryLayer(1)それぞれのオプションの意味を以下に示す
| Option | 説明 | default value |
| EdgesList | 境界層が必要な幾何学モデルの曲線のタグ | - |
| Quads | 境界層でのrecombined要素の生成 | 0 |
| hwall_n | 壁垂直方向のメッシュサイズ | 0.1 |
| thickness | 境界層の最大厚さ | 0.01 |
EdgeListで境界層メッシュを加える曲線を入力し，thicknessで境界層メッシュの厚さ，hwall_nでメッシュサイズを決定する．Quadsを1にすると，可能な限り三角形メッシュを四角形に変換してくれる
もちろんhwall_n<thicknessである必要がある
構造格子
構造格子では，Transfinite meshesを使ってメッシュサイズを決定する
≫Gmsh reference manual (stable release) - 6.3.2 Structured grids
≫tutorial_python_t6.py - GitLab
gmsh.model.geo.mesh.setTransfiniteCurve(Curveのタグ,分割数,分割タイプ,係数)
gmsh.model.geo.mesh.setTransfiniteSurface(Surfaceのタグ)
gmsh.model.geo.mesh.setRecombine(次元数(=2),Surfaceのタグ)setTransfiniteCurveメソッドでは分割のタイプとして"Progression"と"Bump"を選択することができる
"Progression"を選択すると，指定した係数(>1)倍だけ曲線の向きにメッシュが大きくなっていく．曲線の逆側からメッシュを大きくしたい場合は負の値を入力すればいい
"Bump"を選択すると，指定した係数倍だけ曲線の両端のメッシュの大きさが変化する．(係数>1ならメッシュが大きくなる)
デフォルトの引数としてmeshType="Progression"，coef=1.が設定されているので，曲線のタグと分割数のみ指定した場合は同じ大きさのメッシュで等分される．
setTransfiniteSurfaceメソッドで指定するSurfaceは次の条件を満たしている必要がある
- 点の数が4つ(もしくは3つ)
- 全ての辺が
setTransfiniteCurveメソッドによって分割されている - 向かい合う辺の分割数が同じ (係数や分割方法は違ってもよい)
これによって，Surface内のメッシュがきれいに整えられる
最後にsetRecombineによって三角形メッシュが四角形に直され，きれいな構造格子が完成する
if structure:
#set TransfiniteCurve
NH=128;NW=128;NA=128
list_index=[1,2]
for CurveIndex in [1,2]:
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NA)
prog=1.05
for (CurveIndex,direction) in zip([4,5,7,9],[False,True,True,True]):
prog_tmp = prog if direction else -prog
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NH,meshType="Progression",coef=prog_tmp)
prog=1.04
CurveIndex=10
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NW,meshType="Progression",coef=prog_tmp)
for CurveIndex in [3,6]:
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NW)
prog=1.07
CurveIndex=8
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NW,meshType="Progression",coef=prog)
prog=1.03
for (CurveIndex,direction) in zip([11,12],[False,True]):
prog_tmp = prog if direction else -prog
gmsh.model.geo.mesh.setTransfiniteCurve(CurveIndex,NA,meshType="Progression",coef=prog_tmp)
for SurfaceIndex in [1,2,3,4]:
gmsh.model.geo.mesh.setTransfiniteSurface(SurfaceIndex)
gmsh.model.geo.mesh.setRecombine(2,SurfaceIndex)
else:
(略)
Extrude
ここまで作成してきた2次元のSurfaceをextrudeメソッドを使って押し出して3次元のVolumeにする
≫tutorial_python_t3.py - GitLab
gmsh.model.geo.extrude([(次元数,タグ)のlist],nx,ny,nz,[要素数のlist],[要素厚さのlist],recombineのタグ)1つ目の引数では押し出したいPointやCurve，Surfaceの次元数とそのタグを指定する
2~4つ目の引数では押し出しの方向ベクトルのx，y，z成分を入力する．今回はxy平面をz方向に1 [m]押し出すので0.，0.，1.である
5つ目の引数では各Layerの要素数を指定し，6つ目の引数で各レイヤーの高さを指定する．今回はLayerは1つしかないので[1]，[1.]である
7つ目の引数でrecombine=Trueとすれば，押し出した方向のメッシュができる限り四角に変換される
if structure:
for SurfaceIndex in [1,2,3,4]:
gmsh.model.geo.extrude([(2,SurfaceIndex)],0.,0.,1.,[1],[1.],recombine=True)
else:
gmsh.model.geo.extrude([(2,1)],0.,0.,1.,[1],[1.],recombine=True)PhysicalGroupの追加
addPhysicalGroupメソッドを使ってこれまでに作ったSurfaceやVolumeをPhysicalGroupに追加し，setPhysicalNameメソッドで名前を付ける
このPhysicalGroupはOpenFOAMで解析するときの境界条件の設定などに使われる
gmsh.model.addPhysicalGroup(次元数,[タグのlist],PhisicalGroupのタグ)
gmsh.model.setPhysicalName(次元数,PhisicalGroupのタグ,"名前")extrudeメソッドで新たに追加されたSurfaceやVolumeのタグは自動的に振られるので，現時点ではどのSurfaceに何番のタグが振られているのかわからない
そのため，先にこれ以降のプログラムを完成させて実行し，作成されたmshファイルをGUIアプリで開いてタグを確認する必要がある
左上のメニューからTools→Visibilityでウィンドウを開き，Treeのタブでそれぞれのタグが確認できるので，これを参考にプログラムに入力していく
if structure:
gmsh.model.addPhysicalGroup(2,[1,2,3,4],1) ;gmsh.model.setPhysicalName(2,1,"front")
gmsh.model.addPhysicalGroup(2,[34,56,78,100],2) ;gmsh.model.setPhysicalName(2,2,"back")
gmsh.model.addPhysicalGroup(2,[29,47],3) ;gmsh.model.setPhysicalName(2,3,"inlet")
gmsh.model.addPhysicalGroup(2,[73,91],4) ;gmsh.model.setPhysicalName(2,4,"exit")
gmsh.model.addPhysicalGroup(2,[95],5) ;gmsh.model.setPhysicalName(2,5,"top")
gmsh.model.addPhysicalGroup(2,[69],6) ;gmsh.model.setPhysicalName(2,6,"bottom")
gmsh.model.addPhysicalGroup(2,[21,55],7) ;gmsh.model.setPhysicalName(2,7,"aerofoil")
gmsh.model.addPhysicalGroup(3,[1,2,3,4],8) ;gmsh.model.setPhysicalName(3,8,"internal")
else:
gmsh.model.addPhysicalGroup(2,[1],1) ;gmsh.model.setPhysicalName(2,1,"front")
gmsh.model.addPhysicalGroup(2,[54],2) ;gmsh.model.setPhysicalName(2,2,"back")
gmsh.model.addPhysicalGroup(2,[25,29],3);gmsh.model.setPhysicalName(2,3,"inlet")
gmsh.model.addPhysicalGroup(2,[37,41],4);gmsh.model.setPhysicalName(2,4,"exit")
gmsh.model.addPhysicalGroup(2,[45],5) ;gmsh.model.setPhysicalName(2,5,"top")
gmsh.model.addPhysicalGroup(2,[33],6) ;gmsh.model.setPhysicalName(2,6,"bottom")
gmsh.model.addPhysicalGroup(2,[49,53],7);gmsh.model.setPhysicalName(2,7,"aerofoil")
gmsh.model.addPhysicalGroup(3,[1],8) ;gmsh.model.setPhysicalName(3,8,"internal")以上でモデル作成は終了である
synclonize~finalize
モデルの作成が終わったので，synchronizeメソッドで内部のCADデータをgmshの形式に変換し，generateメソッドでメッシュを切り，writeメソッドでファイルに書き出したあとにfinalizeメソッドでgmshを終了する
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(3)
gmsh.write('airfoil.msh')
gmsh.finalize()generateメソッドの引数は切りたいメッシュの次元数で，今回は3次元メッシュを切るよう設定している
writeメソッドではファイルの形式を指定することでmshファイル以外の形式のメッシュファイルにも出力できる
gmshToFOAMの実行~boundaryの書き換え
これについては次の記事を参考にしてほしい
≫【OpenFOAM v2006】二次元翼型解析
まとめ
pythonでgmshのモジュールを使ってメッシュファイルの作成を行った
次はCADデータのインポートに挑戦したい
↓関連記事
コメント