Calculix 20本ノック:11本目: exampleの実行 Contact/Hertz2D

10本完了したので、目標を20本に変えてみた。
その11本目。

Calculixとは

フリーの解析ソフト。

• ccx　ソルバー
• cgx　プリポスト

環境

• WSL

• Xming

• Calculix 2.11

• gmsh

  gmshの機能を使ってメッシングすることがある。

• imagimagick

  hcpyで画像作成するとき、convertするので使う。

• pip2

• matplotlib

• gnuplot

ノック11本目: exampleの実行 Contact/Hertz2D

データの入手

公式サイトからのリンクで行けるGithubのexample。

下記サイトのリンクから、自分のCalculixのバージョンにあったファイルを探す。

https://github.com/mkraska/CalculiX-Examples/releases

Contact

Hertz2D

平面ひずみ要素で、四角と1/4円をぶつけて接触させる。その時の応力分布をみる。
強制変位と荷重の2パターンで1/4円を移動させている。

ccx2.11では正しい結果がでない？

ccx2.11では収束しなかった。（頑張って条件をいじって結果をだしても、READMEのコンターの分布と同じにならなかった。）
ccx2.16のwindows版を使ったら、一発で通った。結果はREADMEの結果と同等のものが作れた。

READMEにあるmonitor.pyの結果と同じように収束しなかった。
このEXAMPLE全体に言えるが、初期設定では収束しないやつとかがある。
これもccxのバージョンのせいだと思っている。

内容

Hertzの接触を、FEMで求める。
球体とか円柱が、他の面等に接触した時の応力分布等を表したHertzの式というのがあるらしい。
周波数のHertz（Hz）と同じ人が考えたらしい。

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

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

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

サンプルの実行

cgx -b pre.fbd

をすると、メッシュ、境界条件などのファイルが出力される。
これらのファイルは、.inpで読み込まれて使用する。

ccx Hertz

とすれば、同じ結果が出てくる…はずだった。

収束しない

STEP1は強制変位で動かしている。
STEP2は荷重変位で動かしている。
強制変位で動かしたほうが、収束はしやすいそうな。

STEP1は収束した。
ただ、STEP１はREADMEに書いてあるようなイテレーション回数で、収束しない。
STEP2の荷重で動かしているほうは収束しないで終わる。

なんで？

ccx2.11では結果が違う。

終息したSTEP1の結果だけを見てみる。

READMEにあるような、対称な応力分布になってくれない。片側の応力が出てこない。
接触面で応力が出ているけれども、サンプル画像では、内側で応力が発生していた。

ccx2.16(windows版)を使ってみる。

READMEと同じ結果が出てくる。
きっとこれが正しい結果なのだろう。同じ結果なので載せない。

ccxのバージョンは新しいのに限る。
以下の「余計なこと」は、すべてccx2.16でやった。

いろいろやってみる。

master slaveを入れ替えてみた

下記のようにマスタースレーブをひっくり返した。

*CONTACT PAIR, INTERACTION=contact, TYPE=SURFACE TO SURFACE
Ssblk,Sscyl

結果は変わらなかった。

平面応力要素でやった（デフォルトは平面ひずみ要素）

平面ひずみだから意味のある解析なのですが、あえて平面応力でやってみました。

ファイルall.mshの

*ELEMENT, TYPE=CPE8, ELSET=Eall

を

*ELEMENT, TYPE=CPS8, ELSET=Eall

に書き換えた。

cgx Windows版をメインで使わない理由

windowsで最新版(2.16)が使えるのに、わざわざlinux上の2.11を使う理由は、

• cgxのGUIのメニューが一部使用できない箇所があった
  （常にではなく、frdのロード時（？）など限られたタイミングで起こる。）。
  下記画像、（左）WSL（Xming）のcgx2.11の右クリックメニュー（右）windows版のcgx2.16の右クリックメニュー
  windows版ではdatasetsが出ない。載せていないが、矢印を展開した先でも似たような現象が散見された。

だけど、

• ccxは、windows版でも可。
  ccxはむしろ新しい方がいい。

Calculix 10本ノック: 10本目: exampleの実行 Kasten

Calculixとは

フリーの解析ソフト。

• ccx　ソルバー

• cgx　プリポスト

自分の環境

• WSL

• Xming

• Calculix 2.11

• imagimagick

データの入手

公式サイトのリンクから行けるGithubのページCalculixのexample。

下記サイトのリンクから、使用しているCalculixのバージョンに合ったファイルを探す。

https://github.com/mkraska/CalculiX-Examples/releases

ノック10本目: exampleの実行 Kasten

Kasten

Kastenは、ドイツ語で"箱"。

大変形、塑性、不整のある座屈（小さい内圧）、シェル要素とreadmeの書き出しにあった。Buckling with imperfectionは不整のある座屈と翻訳するらしい。

箱が座屈するモデル。 S8 (8-節点 四角形シェル要素)

 cd ~/CalculiX-Examples-2.11/CalculiX-Examples-2.11/Kasten

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

exampleには存在しないけれど

1/4モデルを作成して、

2. 材料非線形解析なしの場合

3. 非線形解析なしの場合。

を用意した。

1/4モデルの作成

1/4モデル作成のため、pre.fbdを下のように書き直した。

 pnt p1 0 0 0
 pnt p2 50 0 0
 pnt p3 0 50 0
 line l1 p1 p2 12
 line l2 p1 p3 12
 ​
 seta lower all
 swep all upper tra 0 0 150 20
 flip A002
 ​
 seta miryz l L004
 seta mirxz l L005
 ​
 plot la all
 plus pa all
 plus sa all
 rot y
 rot d 10
 rot r 10
 frame
 hcpy png
 #sys mv hcpy_1.png geo.png
 ​
 div all auto 3
 elty all qu8
 mesh all
 send all abq
 ​
 send all abq pres -0.01
 send lower abq nam
 send upper abq nam
 send mirxz abq spc 246
 send miryz abq spc 156

1/4モデルの為に、Kasten.inpの*boundary部分を下記のように書き直した。

 *boundary
 Nlower,1,6
 Nupper,1,6
 *include,input=mirxz_246.bou
 *include,input=miryz_156.bou

1/4モデルの為に、post.fbdをミラーコピーするように下記のように書き直した。

 copy all base1 mir x
 move base1 tra 100 0 0
 copy all base2 mir y
 move base2 tra 0 100 0

1/4モデルの解析結果

非線形解析なしの場合

非線形がない場合を見るために、Kasten.inpのStepのnlgeomを削った。

材料非線形がない場合

材料非線形がない場合を見るために、Kasten.inpの*plasticと、その下3行を削った。

あと、上記だけだと収束しなくなったので、*controlを追加して、収束条件を甘くしてしまった。

 *CONTROLS,PARAMETERS=FIELD
 5.0e-1,1e-1,1e-2,,2e-8,1.e-5,1.e-3,1.e-8

10本完成してしまった。次は25本くらいを目標に挙げておく。

やはり、example頼りになりそう。

Calculix 10本ノック: 9本目: exampleの実行 Linear/Mesh1

Calculixとは

フリーの解析ソフト。

• ccx　ソルバー
• cgx　プリポスト

環境

• WSL

• Xming

• Calculix 2.11

• gmsh

• imagimagick

• pip2

• matplotlib

データの入手

公式サイトからのリンクで行けるGithubにあるCalculixのexample。

下記サイトのリンクから、自分のCalculixのバージョンにあったファイルを探す。

https://github.com/mkraska/CalculiX-Examples/releases

ノック9本目: exampleの実行 Linear/

Mesh1

アルミ製皿の固有値解析。
ポストのshapes.fblを使うと、複数の固有値のアニメーションを一気に作れる。
要素はIncompatible mode eight-node brick element (C3D8I)らしい。

ディレクトリの移動

cd ~/CalculiX-Examples-2.11/CalculiX-Examples-2.11/Linear/Mesh1

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

実行

cgx -b pre.fbl
ccx modal
cgx -b shapes.fbl

pre.fblの内容

手打ちで作ってみたら、途中、GUIからframeを使わないと、全体が見えなくなった。

下記で円弧が出現してた。作成途中でできたD002、D003を通る、中心点D001の円弧の命令？（ 14）で書いたline。

line ! D002 D003 D001 14

1/4を作って、コピーで1/1にする。コピーのコマンドは下記。
最後のxは、yz面に対してミラーコピーの意味。yだったらx-z面

copy all new mir x
copy all new mir y

すべてコピーなので、不要な重複が出る。マージで消す。
マージは下記。

merg p all 0.001
merg l all 0.001
merg s all 0.001

all.mshの内容

先ほどのpre.fblの結果、all.mshが出力される。

all.mshには、節点座標とエレメントが使用する節点の情報が書いてある。
エレメントはIncompatible mode eight-node brick element (C3D8I)を使っていた。

modal.inpの内容

メッシュを読んでエレメントセットを用意して、

マテリアルを指定して、

上記マテリアルとエレメントセットをプロパティに設定して、

下記で、固有値計算になる。

*STEP*
frequency
〜省略〜
*END STEP

このmodal.inpが、ccxの解析で使う入力ファイルになる。
ccxで解析された結果は、modal.frdとmodal.datに出力される。

プロパティ設定の下の1は、3次元要素なので要らない気がする。あっても動くが。

shapes.fbl

出力した12個の固有値のアニメーションを自動で作るcgxスクリプト。

はじめに、wihleで使用する変数をvaluで用意している。

次に、readでmodal.frdを読み込んでいる。

次に、rotで角度を整えている。

最後に、whileで、すべてのデータセットのアニメーションを作成していた。
Whileについては以前やったので割愛。
whileはendwhileで終わるが、endwhileがテキストエディタの最終行に来たらエラーで止まるのは要注意。

cgxで最短でアニメを描く方法

多分、gifアニメは3行で出力できる。

アニメーションは、moviでフレームを設定したら、dsのaオプションで出力される。
多分imagimagickか何かが必要になると思う。出力されるファイルは、リピートされないgifアニメ。

read modal.frd new
movi frames 30
ds 1 a

Calculix 10本ノック: 8本目:Contact/Shell1

Calculix 10本ノック: 8本目:Contact/Shell1

Calculixとは

フリーの解析ソフト。

• ccx　ソルバー
• cgx　プリポスト

環境

• WSL

• Xming

• Calculix 2.11

• gmsh

• imagimagick

• pip2

• matplotlib

データの入手

公式サイトからのリンクで行けるGithubにあるCalculixのexample。

下記サイトのリンクから、自分のCalculixのバージョンにあったファイルを探す。

https://github.com/mkraska/CalculiX-Examples/releases

ディレクトリの移動。

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

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

内容

下のような8節点シェルで作ったモデルの接触解析。
ハニカム（ハチの巣、ハニーコム）構造のサンドイッチの最小構成みたいな形の上に半球。
半球に対して、z方向に-2の強制変位をかけている。

初期設定でうまく流れなかったが、contact pairのmaster,slaveを入れ替えたらうまくいった。

Too many cutbacks

初期.inpのまま解析実行したら、下記エラーで止まった。

too slow convergence; the increment size is decreased to 3.125000e-02
the increment is reattempted

*ERROR: too many cutbacks
best solution and residuals are in the frd file

cutbacksの回数を増やしたら、止まらないんじゃないかと思い、cutbacksを増やす方法を検索。
controlをいじるといいらしい。 controlには引数が多すぎる。cutbackの回数以外は、初期設定にしておいた。
もしかすると、省略したら（何も書かずにコンマで区切ったら）よかったのかも。

下記*controlを*step(stepの外でも内でもいい？)に追加したら、このエラーは避けられたが…

** cutback 15
*CONTROLS, PARAMETERS=TIME INCREMENTATION
4,8,9,16,10,4,,15,,
0.25,0.5,0.75,0.85,,,1.5,

結局別のエラーが出た。
下記は結局出てきたエラー。

too slow convergence; the increment size is decreased to 7.629395e-06
the increment is reattempted

*ERROR: increment size smaller than minimum
best solution and residuals are in the frd file

これを解消する気にはならなかった。
*controlは消して、別の方法を探った。

contact pairがおかしかった

下記の2行目で、slave,masterを決めていた。

*contact pair, interaction=tool, type=surface to surface
Stopneg,Sind

dependent (slave), independent(master)の順番。

これを、下記のように逆さにしたら収束した。

*contact pair, interaction=tool, type=surface to surface
Sind, Stopneg

ちなみにsindは下記。

stopnegは下記。

 

結果

表示しているのはds 2 e 23で出せる、最悪主応力(worst Principal Stress)
READMEにある結果と最大値最小値一緒だし、同じ結果でしょ。多分。

Calculix 10本ノック: 7本目:exampleの実行 Contact/Shellassembly

Calculixとは

フリーの解析ソフト。

• ccx　ソルバー
• cgx　プリポスト

環境

• WSL

• Xming

• Calculix 2.11

• gmsh

  gmshの機能を使ってメッシングすることがある。

• imagimagick

  hcpyで画像作成するとき、convertするので使う。

• pip2

• matplotlib

データの入手

公式サイトからのリンクで行けるGithubのexample。

下記サイトのリンクから、自分のCalculixのバージョンにあったファイルを探す。

https://github.com/mkraska/CalculiX-Examples/releases

ノック7本目: exampleの実行 Contact/Shellassembly

Contact

shell_assembly

4種類の.inpが入っていた。

• ペナルティ法でnode to surfaceの接触
• ペナルティ法でsurface to surfaceの接触
• *tieを使ったMPC接触
• *equationを使ったMPC接触

固有値計算してる。

他のサイトでは、「接触を考慮した固有値解析はできない」と書いてある。 )
不安になるが、まぁ、できるんならいいんじゃない…？
いや、良くないのかなぁ。「接触 固有値」でググると、皆、いろいろ工夫してるみたいだし。

どうやら、.inpファイルのstep, perturbationに鍵があるらしい。

READMEざっくり訳

*tieとペナルティ法を使った、MPC接触は、シェルのエッジとシェルの表面の間の、硬い接続（？）になります。 *tieを使った、Node to surfaceペナルティ法やMPC接触は似た結果になります。 ペナルティ接触の固有値の周波数は、ふつう、MPC接触を使った場合よりもきわめて低く、もっともらしいです。ペナルティ接触はcomplianceを足しているからです。

Surface to surface接触は、MPC接触とは全く異なった固有周波数になります。0付近に5つの周波数しかありません（6つの剛体モードがあるべき）。

（surface to surfaceもnode to surfaceのどちらも）ペナルティ接触は、モード解析が必要です。モード解析は、静的なstepに先立って摂動ステップとして行われます。 摂動解析の最初の結果のインクリメントはモードシェイプではないことに気を付けてください。

*equationを使ったMPC接触は、シェルとフェイス接触を、ヒンジ接続にします。（7剛体モード）

う〜ん分からない。訳が合ってるのかも不安。

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

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

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

最新版だとShell0に入っているサンプル。

操作

cgx -b pre.fbd
ccx tie
cgx tie.frd

こんな感じでtieの部分をequ、pc-ss、pc-nsに替えて実行する。

pre.fbd

seto〜setcを使ってセットを作っていた。

上側のプレート1をswepで作って、

下側のプレート2をswepで作って、

qu4（inpでいうs4（4節点シェルエレメント））を指定してメッシュ作成して(all.mshにS4と記入される)、

プレート1のサーフェイスを出力して(S1.sur)、

プレート2のエッジ（L007 ）を出力して(k1.sur)、

すべてのメッシュを出力してして(all.msh)、

プレート1のサーフェイスとプレート2のエッジ(L007)を123456のMPCとして出力していた(k1.equ)。

PgDn, PgUpでplotするsetの切り替えができた！

PgDn、PgUpでセット間を行き来できるのは知らなかった。

plot fa all

と打った後に、PgDn、PgUpすると、作成した各セットを行き来できるので、いちいちplotコマンドを打たなくていい。作ってる最中に便利かも。

.inpを見る

tie.inp

*tie,name=t1,position tolerance=0.1
Sk1,Ss1

とあった

pc-ss.inp

stepの部分

*step,nlgeom
*static
*end step
*STEP,perturbation
*frequency
~~~
~~~
*END STEP

と書いてあった。
「PERTUBATIONってなんだ？」となったので、取説を検索すると、

The parameter PERTURBATION is allowed for FREQUENCY and BUCKLE steps only. If it is specified, the last *STATIC step is taken as reference state and used to calculate the stiffness matrix.

パラメーターPERTURBATION は*FREQUENCY(固有値解析) *BUCKLE(座屈解析)にだけ使用できます。もし、これが指定されると、最後の *STATIC stepの剛性マトリックスが参照され、計算に使用されます。

たしかに、直前に*staticしてるstepがある。その結果を読み込んで、固有値計算をしているらしい。

接触の部分

*surface interaction, name=Klebung
*surface behavior, pressure-overclosure=tied
1.E7
*friction
1,1.E7
*contact pair, interaction=Klebung, type=surface to surface, adjust=1
Sk1,Ss1

pc-ss.inpとは、pressure-overclosure=tiedと、type=surface to surfaceが違う。

*frictionの下の数値が、摩擦係数と粘性slopeらしい。取説には下の様に書いてある。

For face-to-face penalty contact with PRESSURE-OVERCLOSURE=TIED the value of the friction coefficient is irrelevant.

PRESSURE-OVERCLOSURE=TIEDを使ったface to faceペナルティ接触にとっては、摩擦係数は無関係です。

とあるので、1っていう数字は要らないんじゃないかと思ったが、なきゃ無いでエラーが起こるのでそのままにした。

pc-ss.inp

接触の部分

*surface interaction, name=Klebung
*surface behavior, pressure-overclosure=linear
1.E7,10000,10000
*friction
1,1.E7
*contact pair, interaction=Klebung, type=node to surface, adjust=1
Sk1,Ss1

pc-ns.inpとは、pressure-overclosure=linearと、type=node to surfaceが違う。

equ.inp

これは、下記で読み込んでいる、k1.equにの条件が書かれているところに特徴があった。

*include, input=k1.equ

*equation を使って、MPCの式が書かれていた。

*EQUATION4
4
6,1,-1.000000000000
2,1,0.800000000000
11,1,0.200000000000
12,1,0.000000000000

4というのが、式の項の数で、
その下に、ノード番号、拘束方向、変数の係数が書かれている。変数は、ノードの拘束方向の変位。

k1.equを見ると、近い4つのノード同士で式を作っていた。ノードがどうやって選ばれたのかは知らない。

結果

固有値の結果は、.datファイルに出力される。

cgxのwindowの「Time:」って書いてある横の数値は、そのまま周波数として読めそう（CYCLES/TIME）。

*el fileS

PgDn, PgUpでplotするdatasetの切り替えができた！

plotコマンドのところでも書いたが、dsコマンドでも、PgDn, PgUpデータセットの切り替えができる。
下記コマンド等、適当に入力後、PgDn, PgUpで、データセットを切り替えられる。

ds 1 e 1

問題点

cgxのオプション

cgxを使って結果の読取りをする際、cgxの起動時に -bオプションををつけると結果が表示されない。何もつけない（自動）か、-readで出てくると思う。

cgxのdataset

cgxのGUIからdatasetを選択すると、アプリケーションが勝手に終了して下記の表示が出た。

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

表示させるためには、コマンドから下記の様に入力。

ds 1 e 4

データセット1の4番目の結果を表示。何番に何が入っているかは、出力要求によって違う。

windows版ならGUIは問題なく動くと思う。