LoginSignup
4

More than 5 years have passed since last update.

Python matplotlibによる2次元FEM応力解析結果の図化

Last updated at Posted at 2018-02-05

はじめに

2次元FEM応力解析を行ったあと,やはりすぐに結果を見てみたいものである.
このため,昔から図化プログラムを作っており,いつか解説を作ろうと思っていたのだが,めんどうくさくてなかなかできなかった.
今回は,少し時間ができたので,それを試みたものである.

やっていること

私の場合は,2次元FEM応力解析結果に対し,Python matplotlibにより,以下の6個の出力図を作っている.

  • 第1主応力コンター
  • 第2主応力コンター
  • 最大せん断応力コンター
  • 第1主応力ベクトル図
  • 第2主応力ベクトル図
  • 変位モード図

上記のうち,コンターといっても,要素を,その要素の平均応力の値に応じて色分けしただけである.また,最大せん断応力は,第1主応力と第2主応力から定まるモール円の半径と定義している.

解析結果は,以下の記事で紹介したプログラムによるものを対象としている.

FEM解析結果の出力フォーマットが変われば,図化プログラムの入力部分も変わってくるが,そこらへんは適宜変更ということで.

上述6個の出力図を作ったあと,関連性をチェックしやすくするため,ImageMagickにより1枚の図にまとめている.
まとめ方は以下の通り.画像を作るのが面倒なので,表として書いているがご容赦を.

 第1主応力 
 コンター図
 第1主応力 
 ベクトル図 
 第2主応力 
 コンター図
 第2主応力 
 ベクトル図 
最大せん断応力
 コンター図
  変位
 モード図

作例はこちら.

環境

  • MacBook Pro (Retina, 13-inch, Mid 2014)
  • macOS High Sierra
  • Python 3.6.4
  • ImageMagick 7.0.7

このプログラムの問題点,それはとても遅いことである.

GMTで作図したら,早くなる.
GMTでの作図例(https://qiita.com/damyarou/items/23607fb0a252b53ebebb)

作図プログラム解説

私は,個人的にpyplotを直接呼び込むのが好きなので,そのようにプログラムしている.
以下,プログラム全文を区切って説明する.少し長いがご勘弁を.

モジュール読み込み

  • sysは,コマンドライン入力を使うためにインポートしているが,ここでは使っていない.
  • フォント指定を行うために,FontProperties を読み込んでいる.桁数を揃えて見栄えを良くするため,一部等幅フォントを使っている.
#***************************************
# Diagrams for 2D FEM result
#***************************************
import matplotlib.pyplot as plt
import numpy as np
import sys
from matplotlib.font_manager import FontProperties

データ読み込み

この部分は,各人の持っているFEM解析プログラムに合わせる必要あり.
ここでは,関数の引数に1要素の節点数を示す変数 nod を入れており,三角形要素・四角形要素双方に対応させている.

def DATAINP(nod,fnameR):
    # Data input
    fin=open(fnameR,'r')
    text=fin.readline()
    text=fin.readline()
    text=text.strip()
    text=text.split()
    npoin=int(text[0])
    nele =int(text[1])
    nsec =int(text[2])
    npfix=int(text[3])
    nlod =int(text[4])
    NSTR =int(text[5])

    node=np.zeros([nod+1,nele],dtype=np.int)
    x   =np.zeros(npoin,dtype=np.float64)
    y   =np.zeros(npoin,dtype=np.float64)
    disx=np.zeros(npoin,dtype=np.float64)
    disy=np.zeros(npoin,dtype=np.float64)
    sig1=np.zeros(nele,dtype=np.float64)
    sig2=np.zeros(nele,dtype=np.float64)
    angl=np.zeros(nele,dtype=np.float64)
    taum=np.zeros(nele,dtype=np.float64)

    text=fin.readline()
    for i in range(0,nsec):
        text=fin.readline()

    text=fin.readline()
    for i in range(0,npoin):
        text=fin.readline()
        text=text.strip()
        text=text.split()
        x[i]=float(text[1])
        y[i]=float(text[2])

    text=fin.readline()
    for i in range(0,npfix):
        text=fin.readline()

    text=fin.readline()
    for i in range(0,nele):
        text=fin.readline()
        text=text.strip()
        text=text.split()
        for j in range(0,nod):
            node[j,i]=int(text[j+1]) #node,sectuon number

    text=fin.readline()
    for i in range(0,npoin):
        text=fin.readline()
        text=text.strip()
        text=text.split()
        disx[i]=float(text[1]) # displacement in x-direction
        disy[i]=float(text[2]) # displacement in y-direction

    text=fin.readline()
    for i in range(0,nele):
        text=fin.readline()
        text=text.strip()
        text=text.split()
        sig1[i]=float(text[4])
        sig2[i]=float(text[5])
        angl[i]=float(text[6])
        taum[i]=0.5*abs(sig1[i]-sig2[i])
    fin.close()
    return npoin,nele,node,x,y,disx,disy,sig1,sig2,angl,taum

応力コンター作図(要素の色塗り)

関数 SIGCOL

  • ここでは応力値を10段階に分けて色塗りするようにしている.
  • しきい値配列 sg は,解析モデルやケースにより異なるため,あとで出てくる「入力ファイル指定」の際,指定するようにしている.

関数 CON

  • 要素に色塗りを行う関数.
  • 塗り分ける色は10色としている.圧縮側(マイナス応力)は茶色から薄い黄色,引張側(プラス応力)は水色から濃い青として色を指定した.完全に私の趣味である.
  • 色塗りは,plt.fill(xa,ya,facecolor=col[icol],edgecolor='#666666',lw=0.1) のように,facecolor で応力値に応じた色で塗りつぶし,edgecolor で要素境界を灰色で描くようにしている.
  • 応力の最大・最小を示すテキストは,plt.text(0.9,0.9,sinfo,va='top',ha='right',fontproperties=fp,transform=plt.gca().transAxes) の部分で作図領域の右上に表示している.フォントは fontproperties=fp で指定し(ここでは等幅のMenlo.ttc),transform=plt.gca().transAxes でグラフ領域の相対座標指定(左下[0,0],右上[1,1])している.
  • 凡例のカラーバーは自作した.グラフ領域の中に plt.axes([xcb,ycb,wcb,hcb]) で小さな領域を作り,そこを plt.axvspan で区切って塗りつぶしている.
def SIGCOL(sig,sg):
    if sg[8]<sig: icol=9
    if sg[7]<sig  and  sg[8]>=sig: icol=8
    if sg[6]<sig  and  sg[7]>=sig: icol=7
    if sg[5]<sig  and  sg[6]>=sig: icol=6
    if sg[4]<sig  and  sg[5]>=sig: icol=5
    if sg[4]>=sig and  sg[3]<=sig: icol=4
    if sg[3]>sig  and  sg[2]<=sig: icol=3
    if sg[2]>sig  and  sg[1]<=sig: icol=2
    if sg[1]>sig  and  sg[0]<=sig: icol=1
    if sg[0]>sig: icol=0
    return icol

def CON(nod,nnn,_sig,node,x,y,nele,ss,fp):
    _s=[]
    for i in range(1,len(ss)-1):
        _s=_s+[float(ss[i])]
    sg=np.array(_s)
    col=['#cc0000','#ff0066','#ff9966','#ffcc66','#ffffaa','#00ffff','#00ccff','#0099ff','#0000ff','#000088']
    if nnn==0:
        # first principal stress
        ls0='* First principal stress'
        ls1='sig1_max={0:8.3f} MPa'.format(np.max(_sig))
        ls2='sig1_min={0:8.3f} MPa'.format(np.min(_sig))
    if nnn==1:
        # second principal stress
        ls0='* Second principal stress'
        ls1='sig2_max={0:8.3f} MPa'.format(np.max(_sig))
        ls2='sig2_min={0:8.3f} MPa'.format(np.min(_sig))
    if nnn==2:
        # maximum shearing stress
        ls0='* Maximum shearing stress'
        ls1='tau_max={0:8.3f} MPa'.format(np.max(_sig))
        ls2='tau_min={0:8.3f} MPa'.format(np.min(_sig))
    for ne in range(0,nele):
        sig=_sig[ne]
        icol=SIGCOL(sig,sg)
        nn=node[0:nod,ne]-1
        xa=x[nn]
        ya=y[nn]
        plt.fill(xa,ya,facecolor=col[icol],edgecolor='#666666',lw=0.1)
    print(ls0)
    sinfo=ls0+'\n'+ls1+'\n'+ls2
    plt.text(0.9,0.9,sinfo,va='top',ha='right',fontproperties=fp,transform=plt.gca().transAxes)
    # Legend (color bar)
    xcb=0.2; ycb=0.8; wcb=0.2; hcb=0.02
    plt.axes([xcb,ycb,wcb,hcb])
    plt.xlim([0,10])
    plt.ylim([0,1])
    plt.axis('off')
    for i in range(0,10):
        x1=float(i)
        x2=float(i+1)
        plt.axvspan(x1,x2, facecolor=col[i], alpha=1)
    for i in range(0,11):
        plt.text(float(i),-0.2,ss[i],fontsize=12,va='top',ha='center')
    plt.text(5,1.1,'Comp.  (MPa)  Tens.',fontsize=10,va='bottom',ha='center')

変位モード作図

  • 変位の最大値が,図上で 0.5m (変数 scl_dis=0.5)となるように変位を変換している.
  • 変位後の形状は plt.fill(xa,ya,facecolor='#ffffcc',edgecolor='#aaaaaa',lw=0.5) で要素を薄い黄色で,境界を薄い灰色で着色している.
  • 変位前の形状は plt.fill(xa,ya,facecolor='none',edgecolor='#aaaaaa',lw=0.2) で要素境界のみを薄い灰色で着色している.facecolor='none' になかなか気が付かなかった.
def DIS(nod,dmax,disx,disy,x,y,node,nele,fp):
    # displacement
    scl_dis=0.5
    ls0='* Displacement mode'
    ls1='disx_max={0:8.3f} mm'.format(np.max(disx))
    ls2='disx_min={0:8.3f} mm'.format(np.min(disx))
    ls3='disy_max={0:8.3f} mm'.format(np.max(disy))
    ls4='disy_min={0:8.3f} mm'.format(np.min(disy))
    rx=x+disx/dmax*scl_dis
    ry=y+disy/dmax*scl_dis
    for ne in range(0,nele):
        nn=node[0:nod,ne]-1
        xa=rx[nn]
        ya=ry[nn]
        plt.fill(xa,ya,facecolor='#ffffcc',edgecolor='#aaaaaa',lw=0.5)
    for ne in range(0,nele):
        nn=node[0:nod,ne]-1
        xa=x[nn]
        ya=y[nn]
        plt.fill(xa,ya,facecolor='none',edgecolor='#aaaaaa',lw=0.2)
    print(ls0)
    sinfo=ls0+'\n'+ls1+'\n'+ls2+'\n'+ls3+'\n'+ls4
    plt.text(0.98,0.98,sinfo,va='top',ha='right',fontproperties=fp,transform=plt.gca().transAxes)

応力ベクトル作図

  • 矢印を書くのは面倒なので,引張応力は青線,圧縮応力は赤線で表示するようにした.
  • 応力の最大値が図上で 1.5m (変数 scl_vec=1.5)となるように線の長さを変換している.
  • 凡例は,引張応力を示す青線と,圧縮応力を示す赤線を長さゼロのダミーで描き,凡例化し左上に表示している.
def VEC(nod,nnn,vmax,_sig,angl,x,y,node,nele,fp):
    scl_vec=1.5
    if nnn==4:
        # first principal stress
        _ang=angl
        ls0='* First principal stress'
        ls1='sig1_max={0:8.3f} MPa'.format(np.max(_sig))
        ls2='sig1_min={0:8.3f} MPa'.format(np.min(_sig))
    if nnn==5:
        # second principal stress
        _ang=angl+90.0
        ls0='* Second principal stress'
        ls1='sig2_max={0:8.3f} MPa'.format(np.max(_sig))
        ls2='sig2_min={0:8.3f} MPa'.format(np.min(_sig))
    xg=np.zeros(nele)
    yg=np.zeros(nele)
    for ne in range(0,nele):
        nn=node[0:nod,ne]-1
        xa=x[nn]
        ya=y[nn]
        plt.fill(xa,ya,facecolor='#ffffcc',edgecolor='#aaaaaa',lw=0.2)
        xg[ne]=np.average(xa)
        yg[ne]=np.average(ya)
    al=abs(_sig)/vmax*scl_vec
    x1=xg-0.5*al*np.cos(np.radians(_ang))
    y1=yg-0.5*al*np.sin(np.radians(_ang))
    x2=xg+0.5*al*np.cos(np.radians(_ang))
    y2=yg+0.5*al*np.sin(np.radians(_ang))
    col=['#ff0000','#0000ff']
    for ne in range(0,nele):
        if _sig[ne]<0.0:
            icol=0
        else:
            icol=1
        plt.plot([x1[ne],x2[ne]],[y1[ne],y2[ne]],color=col[icol],lw=0.5)
    print(ls0)
    sinfo=ls0+'\n'+ls1+'\n'+ls2
    plt.text(0.98,0.98,sinfo,va='top',ha='right',fontproperties=fp,transform=plt.gca().transAxes)
    # Legend
    plt.plot([xmin],[ymin],'-',color=col[1],lw=1.5,label='Tension')
    plt.plot([xmin],[ymin],'-',color=col[0],lw=1.5,label='Compression')
    plt.legend(loc='upper left')        

入力ファイル指定

3節点要素用と4節点要素用の2つの関数を準備した.作図事例は4節点要素のものである.ここで指定しているものは以下の通り.

  • 1要素節点数
  • 入力ファイル名(FEM解析結果ファイル)
  • x方向作図範囲
  • y方向作図範囲
  • 応力区分(応力コンターを作る上での応力区分)

応力区分は,応力コンターを作る上での区分であり,等間隔でなくても構わない.
入力は,カラーバーの凡例に書き込むため,文字列入力としている.
9個の応力値に加え,最初と最後は「カラ」のリストである.
これは,凡例表示したカラーバーを見ればわかると思うが,カラーバー両端の数値は不要であることを意味する.
すなわち,def INP4() の事例でいえば,圧縮側 -0.8MPa より小さな応力は茶色,引張側 +0.8MPa より大きな応力は濃い青で表示されており,カラーバーの両端に数値は不要であるため.

もちろん,この部分は,入力ファイルや作図範囲,応力コンターの区切り値を変えたい場合には書き換える必要がある.

def INP3():                   # 関数名(3節点要素用)
    nod=3                     # 1要素節点数
    fnameR='out_arch3.txt'    # 入力ファイル名(FEM解析結果ファイル)
    xmin=-7 ; xmax=7  ; dx=1  # x方向作図領域
    ymin=64 ; ymax=77 ; dy=1  # y方向作図領域
    ss=['','-10','-6',-4,'-2','0','1','2','3','5',''] # 応力区分
    return nod,fnameR,xmin,xmax,dx,ymin,ymax,dy,ss

def INP4():                   # 関数名(4節点要素用)
    nod=4                     # 1要素節点数
    fnameR='out_arch4.txt'    # 入力ファイル名(FEM解析結果ファイル)
    xmin=-6 ; xmax=6  ; dx=1  # x方向作図領域
    ymin=175; ymax=185; dy=1  # y方向作図領域
    ss=['','-.8','-.6','-.4','-.2','0','.2','.4','.6','.8',''] # 応力区分
    return nod,fnameR,xmin,xmax,dx,ymin,ymax,dy,ss

メインルーチン(データ読み込みと作図準備)

  • データを読み込んだあと,作図のためのいくつかの処理を行う.
  • ここでは,細かいことは言わず,1MPa = 100t/m2 である.
  • vmax(応力の最大値)と dmax(変位の最大値)は,作図上,応力ベクトル線長あるいは変位を制御するために用いている.
#===================
# Main routine
#===================
#nod,fnameR,xmin,xmax,dx,ymin,ymax,dy,ss=INP3()
nod,fnameR,xmin,xmax,dx,ymin,ymax,dy,ss=INP4()
npoin,nele,node,x,y,disx,disy,sig1,sig2,angl,taum=DATAINP(nod,fnameR)

print('xmin, xmax ={0:10.3f} {1:10.3f}'.format(np.min(x),np.max(x)))
print('ymin, ymax ={0:10.3f} {1:10.3f}'.format(np.min(y),np.max(y)))
#stry=input('Continue? y/n ')
#if stry=='n': sys.exit()

# t/m2 -> MPa
sig1=sig1*0.01
sig2=sig2*0.01
taum=taum*0.01
vec1_max=max(abs(np.max(sig1)),abs(np.min(sig1)))
vec2_max=max(abs(np.max(sig2)),abs(np.min(sig2)))
vmax=np.max([vec1_max,vec2_max])
# m -> mm
disx=disx*1000.0
disy=disy*1000.0
dmax=np.max([np.max(np.abs(disx)),np.max(np.abs(disy))])

メインルーチン(作図部分)

  • 等幅フォントで表示したい部分に対し,fp でフォントとそのサイズを指定している.
  • 出力画像ファイル名はリストで指定し,順次作図・保存していく.
fsz=14
fp = FontProperties(fname=r'/System/Library/Fonts/Menlo.ttc', size=fsz)
if nod==3: flist=['_fig_pl3_con1.png','_fig_pl3_con2.png','_fig_pl3_con3.png',
                  '_fig_pl3_disp.png','_fig_pl3_vec1.png','_fig_pl3_vec2.png']
if nod==4: flist=['_fig_pl4_con1.png','_fig_pl4_con2.png','_fig_pl4_con3.png',
                  '_fig_pl4_disp.png','_fig_pl4_vec1.png','_fig_pl4_vec2.png']
for nnn,fnameF in enumerate(flist):
    ixx=10
    fig = plt.figure(figsize=(ixx,ixx/(xmax-xmin)*(ymax-ymin)),facecolor='w')
    plt.rcParams['font.family'] = 'Arial Narrow'
    plt.rcParams["font.size"] = fsz
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.xlabel('x-direction (m)')
    plt.ylabel('y-direction (m)')
    plt.gca().spines['right'].set_visible(False)
    plt.gca().spines['top'].set_visible(False)
    plt.gca().yaxis.set_ticks_position('left')
    plt.gca().xaxis.set_ticks_position('bottom')
    plt.gca().set_aspect('equal', adjustable='box')

    if nnn==0: CON(nod,nnn,sig1,node,x,y,nele,ss,fp)
    if nnn==1: CON(nod,nnn,sig2,node,x,y,nele,ss,fp)
    if nnn==2: CON(nod,nnn,taum,node,x,y,nele,ss,fp)
    if nnn==3: DIS(nod,dmax,disx,disy,x,y,node,nele,fp)
    if nnn==4: VEC(nod,nnn,vmax,sig1,angl,x,y,node,nele,fp)
    if nnn==5: VEC(nod,nnn,vmax,sig2,angl,x,y,node,nele,fp)

    plt.savefig(fnameF, dpi=300, bbox_inches="tight", pad_inches=0.2)
    plt.clf()

ImageMagickによる処理

Jupyter notebook を使っているので,シェルスクリプトを実行するため,先頭に %%bash を記述している.
やっていることは,以下の通り.

  • 3個のコンター図を縦に結合
  • 2個のベクトル図と変位モード図を縦に結合
  • これら2つの縦長画像を横に結合
  • 最後に扱い良さを考慮して,横幅 1000px にリサイズしている.
%%bash

convert -append _fig_pl4_con1.png _fig_pl4_con2.png _fig_pl4_con3.png _1.png
convert -append _fig_pl4_vec1.png _fig_pl4_vec2.png _fig_pl4_disp.png _2.png
convert +append _1.png _2.png fig_pl4_cv.png
rm _1.png _2.png

convert -resize 1000x fig_pl4_cv.png fig_pl4_cvs.png

以 上

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
4