【gmsh v4.6.0】pythonを使って二次元翼型のメッシュを切る

gmshで二次元翼型解析のためのメッシュを作成する

作成するメッシュは以下の2種類

  • 構造格子
  • 非構造格子+境界層メッシュ

また,mshファイルはpythonを使って作成する

解析に使ったファイルはGithubにアップしてある.
mtkbirdman.com/2D_airfoil

はじめに

OpenFOAMの2次元翼型解析のチュートリアル(≫OpenFOAM_ User Guide_ Turbulent flow over NACA0012 airfoil (2D))を参考に,2次元翼型用のメッシュを作成してみる

イメージはこんな感じ

そして実際に完成したメッシュがこれ

非構造格子+境界層メッシュ
構造格子

実行環境

実行環境を以下に示す

OSWindows10 Home 64bit
CPUIntel Core i5-8250U CPU@ 1.60GHz 1.80GHz
Ubuntu20.04 LTS
python3.8.0
gmsh4.6.0

参考

gmshのインストールについてはこちら
【gmsh v4.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

それでは実際にメッシュを作成していく

全体の流れ

最初にプログラム全体のフローチャートとソースコードを示す

フローチャートの点線部分が構造格子と非構造格子で違う部分で,それ以外のプログラムは共通である

プログラムにおいては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データのインポートに挑戦したい

Sponsored Link

コメント

タイトルとURLをコピーしました