<< Calculix 10本ノック: 5本目: cgxを使ってモデル化→解析 | TOP | Calculix 10本ノック: 7本目:exampleの実行 Contact/Shellassembly >>

20200211

Calculix 10本ノック: 6本目: exampleの実行 Contact/eyebar

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

いきなり難しいサンプルに手を出してしまった。
アルファベット順でexampleを漁ってはだめっぽい。

まずはディレクトリの移動。

cd ~/CalculiX-Examples-2.11/CalculiX-Examples-2.11/Contact/Eyebar

「~」は、サンプルを解凍した場所のアドレス。

モデル形状

下記が解析する形状。赤点が拘束。青点が鏡面対象拘束。赤文字が、Xプラス方向へ2強制変位。緑が接触を考慮する位置。なんでも、荷重で動かすよりも、強制変位で動かした方が、収束しやすいそうな。

model.png

ラインL004に属する、エレメントのフェイス（S3エッジ）と、
ラインL00Fに属する、エレメントのフェイス（S4エッジ）の接触が、
フェイストゥーフェイスで考慮されている。

line.png

形状でいうと、図で接触しそうなあたり。当たり前ですが。
ちょうど重なって見にくくなってる線部分。

contact parts.png

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 p all.png

plot l all

plot l all.png

plot s all

plot s all.png

plot e all

plot e all.png

plot n all

plot n all.png

plot f all

plot f all.png

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()

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が結果として出力された。

which.png

このデータを見たい場合は、GUIのウィンドウにマウスカーソルを入れたうえで下記を入力。

ds 122 e 7

25番目の中に入っている、122が応力のデータセット名なので、122。7以外の数値にも、いろいろな結果が入っている。

上記画像のGUIから選ぶと、Freeglutのエラーで止まる。

freeglut (cgx): Menu manipulation not allowed while menus in use.

上記のようなエラーで。

【Calculixの最新記事】
posted by yuchan at 19:00 | Comment(1) | Calculix