Calculixとは
フリーの解析ソフト。
環境
データの入手
公式サイトからのリンクで行けるGithubのexample。
下記サイトのリンクから、自分のCalculixのバージョンにあったファイルを探す。
https://github.com/mkraska/CalculiX-Examples/releases
eyebar
いきなり難しいサンプルに手を出してしまった。
アルファベット順でexampleを漁ってはだめっぽい。
まずはディレクトリの移動。
cd ~/CalculiX-Examples-2.11/CalculiX-Examples-2.11/Contact/Eyebar
「~」は、サンプルを解凍した場所のアドレス。
モデル形状
下記が解析する形状。赤点が拘束。青点が鏡面対象拘束。赤文字が、Xプラス方向へ2強制変位。緑が接触を考慮する位置。なんでも、荷重で動かすよりも、強制変位で動かした方が、収束しやすいそうな。
ラインL004に属する、エレメントのフェイス（S3エッジ）と、
ラインL00Fに属する、エレメントのフェイス（S4エッジ）の接触が、
フェイストゥーフェイスで考慮されている。
形状でいうと、図で接触しそうなあたり。当たり前ですが。
ちょうど重なって見にくくなってる線部分。
python2 param.py par.eyebar.fbd
と打つと、eyebar.fbdが作成される。
cgxでeyebar.fbd読むと自動的に、全メッシュの出力、接触に必要なフェイスなどを出力はしてくれる。
それらのファイルは、あらかじめ作ってあった.inpが読み込んでくれるので、
ccxでinpを解析すれば結果は出てくる。やったのはこれだけなので、気付いた事を書いていく。
セットを作っていた。
セットを作って、いろいろ指定していた。セットを作って、別ファイルで出力。.inpでそれを読み込んで、拘束や強制変位や接触の条件に使用するノードや要素を指定するために利用していた。
セット名をまとめて作る方法として、以下のようなことをやっていた。
seto set名
〜ノードやメッシュを作成する操作〜
setc set名
この手法で、seto（オープン）とsetc（クローズ）の間の操作で生成されるノードなどを、セットにしていた。
通常はセットを作るときは、setaを使うようです。
seta a n 1
上記は、セット名aに、ノード番号1を追加するやり方です。
セット名をもう一度setaで定義しても、追加されるだけで、
seta a n 2
としても、1は消えない。なので、aには、1と2が入っている状態になる。
ノード1を消したいときは、setrを使う。
setr a n 1
で、セットaから、ノード1が消える。
plot n a
などとして確認できます。
セット名allは、自動的に作られる、全部入りのセット。
plotのオプションを調べる。
cgxで、「エレメントを表示させたい。」「ノードを表示させたい。」「ノード番号を表示させたい。」などというとき、setが用意されているのであれば、plotが一番簡単なのではないかと思う。
plotには、たくさんのオプションが用意されているので、どれを使ったら何が見れるのか分かりにくい。
plot p all
plot l all
plot s all
plot e all
plot n all
plot f all
cgxは、.inpを作るソフトではなかった。
cgxのコマンドで、解析の細かい設定を書き込んだ.inpを作る方法はないみたい。
プリじゃないじゃん。
.inpは手打ちっぽい。
FemapやAbaqus CAEのようなGUIでの設定を期待してはいけないという事ですか。
それはキツイなぁ〜。全部cgxで、それもマウスクリックで設定したかったなぁ〜。
このフォルダの解析を一通り実行するためには、〜\CalculiX-Examples-2.11\CalculiX-Examples-2.11\Scriptsフォルダ内の、
が必要。（monitor.pyは、なくても解析はできる。）
また、post.fbdで、gnuplotを、
monitor.pyでpylabを使用している。
monitor.pyを書き直した
monitor.pyは、cvgファイル（非線形解析をした時のイテレーション回数や残留力などが書かれている）などの情報を図示する。
pylabがうまく動かないので、matplotlib.pyplotに書き換えた。
また、WSL上のmatplotlibは、バックエンドをいじらないと動作しない。
書き直してみた。多分これで動く。
下記コードはhttps://github.com/mkraska/CalculiX-Examplesのscriptを、自分の環境用に書き直したもの。
アドレス書いたし、著作者も特定できるので、問題ないはず。
import sys
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as pltimport numpy
import glob
if len(sys.argv)==1:
print "No jobname given."
files=glob.glob("*.sta")
if len(files)==1:
print "Found", files[0]
job=files[0][:-4]
else:
print "Available .sta files:"
for f in files:
print " ", f
quit()
if len(sys.argv)>1:
print "Jobname:",sys.argv[1]
job = sys.argv[1]
try:
sta=numpy.genfromtxt(job+'.sta',skip_header=2,delimiter=[6,11,7,6,14,14,14])
except:
sta=numpy.array([[ 1. , 1. , 1. , 0 , 0 , 0 , 0 ]])
cvg=numpy.genfromtxt(job+'.cvg',skip_header=4)
if sta.ndim==1:
sta=sta[numpy.newaxis,:]
"""
.sta
0 STEP
1 INC
2 ATT
3 ITRS
4 TOT TIME
5 STEP TIME
6 INC TIME
.cvg
0 STEP
1 INC
2 ATT
3 ITER
4 CONT EL
5 RESID FORCE
6 CORR DISP
7 RESID FLUX
8 CORR TEMP"""
iters=cvg.shape[0] it=range(iters) itinc=cvg.astype(int)[:,1] itstep=cvg.astype(int)[:,0] itdt=numpy.empty([iters]) itsteptime=numpy.empty([iters]) for i in it:
stp = itstep[i]
inc = itinc[i]
if (cvg[i,5]==0.0):
cvg[i,5]=0.5
for j in range(sta.shape[0]):
if (stp==sta.astype(int)[j,0]) and (inc==sta.astype(int)[j,1]):
itdt[i]=sta[j,6]
itsteptime[i]=sta[j,5]
print i, stp, inc, j, itdt[i],itsteptime[i],cvg[i,5]
imax=i
plt.subplot(2,1,1)
plt.title('sta and cvg data of job '+job )
plt.semilogy(it[:imax],itdt[:imax],'-',
it[:imax],cvg[:imax,5],'-',
it[:imax],cvg[:imax,6],'r-')
plt.grid()
plt.legend(['dt','force','disp'],
fontsize='small',framealpha=0.5, loc=2)
sp1=plt.subplot(2,1,2)
plt.plot(it[:imax],itsteptime[:imax],'r-',
it[:imax],itsteptime[:imax],'b-',)
plt.legend(['# cont. el.','step time'],
fontsize='small',framealpha=0.5, loc=2)
plt.ylabel('step time')
plt.xlabel('Iteration')
plt.grid()
sp2=sp1.twinx()
plt.plot(it[:imax],cvg[:imax,4],'r-')
plt.ylabel('# of cont. elements')
plt.savefig(job)
plt.show()
eyebarって何。
.inp
eyebar.inpを見た。
*INCLUDE, INPUT= を使って、対象条件、接触条件、観察したいノード番号、強制変位するノードなどを指定するsetが書かれたファイルを読み込んでいた。
control.sur要らん気がするけど、削って試したわけではない。
板厚情報は＊SOLID SECTIONで指定
*SOLID SECTION で、板厚を設定していた。外枠が1で、内側が2の板厚。
シェル要素なので断面形状はなし。beam要素だったら*BEAM SECTION,ELSET=Eall,MATERIAL=STEEL,SECTION=Rectとか書いて、断面形状を入力していた。
エレメントのセット名と、マテリアル名を指定するのは、beamと同じ。
接触
*SURFACE INTERACTION, NAME=contact
*SURFACE BEHAVIOR, PRESSURE-OVERCLOSURE=LINEAR
1000000,0.1,0.01
*CONTACT PAIR, INTERACTION=contact, TYPE=NODE TO SURFACE
Sseye,Sspin
一行目は、contactと書いてあるけど、この段階ではただの名前（表面相互作用の名前）。
2行目からそのcontactが、どう表面相互作用するか書かれていて、
4行目で、接触のペアを、contactで指定した相互作用させると宣言して
5行目に、スレーブ、マスターの順で接触のペアをセット名で指定していた。
強制変位はSTEP内の*boundaryで指定
強制変位は*boundaryでやるようです。
最初、ファイルを見ていて、何で中心の丸が動いているのか分からなかった。
どこに、荷重とか強制変位書いてあるの？
荷重っぽい、DLOADはアスタリスク二個でコメントアウトしてあるし…
→＊STEPの中の、
*BOUNDARY
Ncontrol,1,1,2
で指定されていた。
解析結果の見方
データセットは、一番数字が大きい25が最終的な結果（2mm変位させたときの結果）。
.cvgファイルの中の、INC(REMENT)の最大値分25が結果として出力された。
このデータを見たい場合は、GUIのウィンドウにマウスカーソルを入れたうえで下記を入力。
ds 122 e 7
25番目の中に入っている、122が応力のデータセット名なので、122。7以外の数値にも、いろいろな結果が入っている。
上記画像のGUIから選ぶと、Freeglutのエラーで止まる。
freeglut (cgx): Menu manipulation not allowed while menus in use.
上記のようなエラーで。
