はじめに
２次元FEM応力解析を行ったあと，やはりすぐに結果を見てみたいものである．
このため，昔から図化プログラムを作っており，いつか解説を作ろうと思っていたのだが，めんどうくさくてなかなかできなかった．
今回は，少し時間ができたので，それを試みたものである．
やっていること
私の場合は，２次元FEM応力解析結果に対し，Python matplotlibにより，以下の６個の出力図を作っている．
- 第１主応力コンター
- 第２主応力コンター
- 最大せん断応力コンター
- 第１主応力ベクトル図
- 第２主応力ベクトル図
- 変位モード図
上記のうち，コンターといっても，要素を，その要素の平均応力の値に応じて色分けしただけである．また，最大せん断応力は，第１主応力と第２主応力から定まるモール円の半径と定義している．
解析結果は，以下の記事で紹介したプログラムによるものを対象としている．
- https://qiita.com/damyarou/items/320bad2052bb5fccd75f
- https://qiita.com/damyarou/items/22912863563ba9866b6a
FEM解析結果の出力フォーマットが変われば，図化プログラムの入力部分も変わってくるが，そこらへんは適宜変更ということで．
上述６個の出力図を作ったあと，関連性をチェックしやすくするため，ImageMagickにより１枚の図にまとめている．
まとめ方は以下の通り．画像を作るのが面倒なので，表として書いているがご容赦を．
| 第１主応力 コンター図 |
第１主応力 ベクトル図 |
| 第２主応力 コンター図 |
第２主応力 ベクトル図 |
| 最大せん断応力 コンター図 |
変位 モード図 |
作例はこちら．
環境
- 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解析プログラムに合わせる必要あり．
ここでは，関数の引数に１要素の節点数を示す変数 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')
入力ファイル指定
３節点要素用と４節点要素用の２つの関数を準備した．作図事例は４節点要素のものである．ここで指定しているものは以下の通り．
- １要素節点数
- 入力ファイル名（FEM解析結果ファイル）
- ｘ方向作図範囲
- ｙ方向作図範囲
- 応力区分（応力コンターを作る上での応力区分）
応力区分は，応力コンターを作る上での区分であり，等間隔でなくても構わない．
入力は，カラーバーの凡例に書き込むため，文字列入力としている．
９個の応力値に加え，最初と最後は「カラ」のリストである．
これは，凡例表示したカラーバーを見ればわかると思うが，カラーバー両端の数値は不要であることを意味する．
すなわち，def INP4() の事例でいえば，圧縮側 -0.8MPa より小さな応力は茶色，引張側 +0.8MPa より大きな応力は濃い青で表示されており，カラーバーの両端に数値は不要であるため．
もちろん，この部分は，入力ファイルや作図範囲，応力コンターの区切り値を変えたい場合には書き換える必要がある．
def INP3(): # 関数名（３節点要素用）
nod=3 # １要素節点数
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(): # 関数名（４節点要素用）
nod=4 # １要素節点数
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 を記述している．
やっていることは，以下の通り．
- ３個のコンター図を縦に結合
- ２個のベクトル図と変位モード図を縦に結合
- これら２つの縦長画像を横に結合
- 最後に扱い良さを考慮して，横幅 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
以 上