Calculixとは
フリーの解析ソフト。
- ccx ソルバー
- cgx プリポスト
環境
WSL
Xming
Calculix 2.11
gmsh
gmshの機能を使ってメッシングすることがある。
imagimagick
hcpyで画像作成するとき、内部でconvertしているので使う。
sudo apt install imagemagickで入手。
pip2
matplotlib
ノック6本目: exampleの実行 Contact/eyebar
データの入手
公式サイトからのリンクで行けるGithubのexample。
下記サイトのリンクから、自分のCalculixのバージョンにあったファイルを探す。
https://github.com/mkraska/CalculiX-Examples/releases
Contact
eyebar
まずはディレクトリの移動。
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フォルダ内の、
- param.py
- monitor.py
が必要。（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 plt#import pylab
import numpy
import glob
# print sys.argv[1]
# job = sys.argv[1]
# processing command line arguments, get the
# jobname
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:
# before the first increment is completed, the sta-file is empty, so we
# fake some contents to not crash the
sta=numpy.array([[ 1. , 1. , 1. , 0 , 0 , 0 , 0 ]])
cvg=numpy.genfromtxt(job+'.cvg',skip_header=4)
# ensure 2dim array in case of single data line
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"""
# iteration counter
iters=cvg.shape[0] # number of iterations
it=range(iters) # range for loop over iterations
# increment number for iteration
itinc=cvg.astype(int)[:,1] # increment number for the iterations
itstep=cvg.astype(int)[:,0] # step number for the iteration
itdt=numpy.empty([iters]) # time step for the iteration
itsteptime=numpy.empty([iters]) # Step time for iterations
for i in it:
stp = itstep[i]
inc = itinc[i]
# ccx writes values below 1e-6 as zero
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
## Plot force residuals, disp correction and time step
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)
# step time and number of contact elements
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().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 725番目の中に入っている、122が応力のデータセット名なので、122。7以外の数値にも、いろいろな結果が入っている。
上記画像のGUIから選ぶと、Freeglutのエラーで止まる。
freeglut (cgx): Menu manipulation not allowed while menus in use.上記のようなエラーで。
