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データのインポートに挑戦したい
↓関連記事
コメント