20200201

DualSPHysicsとBlenderで動画を作った

DesignSPHysics(FreeCADのアドオンで、中にDualSPHysicsが入っている)とBlender(+VisualSPHysics)を使って動画を作成した。

試した

粒子で豆腐っぽい形の液体を作って、それを下に落とす。boundaryは描いてない。妥当な結果か分からない。


以下、やり方等のメモ。

FreeCAD+DesignSPHysicsを使ってSPH計算

インストール

  1. FreeCADをインストールする。

  2. DesignSPHysicsをInstallation instructionを参考にインストールする。

    a. GitHubからダウンロードする。

    b. フォルダの名前をDesignSPHysics に変更して、 Windowsなら、%appdata%/FreeCAD/Mod Linux なら ~/.FreeCAD/Mod に入れる。

    c. 同フォルダ内のDesignSPHysics.FCMacro を、macroフォルダに入れる (Windowsなら%appdata%/FreeCAD/Macro Linuxなら; ~/.FreeCAD/Macro 。)

    d. FreeCADを立ち上げて、マクロ(M)>マクロ…で、DesignSPHysics.FCMacroと表示されていれば正しくインストールされている。

09.png

DesignSPHysics.FCMacroの使用方法

公式の動画がYouTubeに上がっていた。

始める前に見たほうがいいかもしれない。

バグ?FreeCAD0.18、DesignSPHysics v0.6.0.1912-12-3

円柱をx500mm,y500mmに置いて、Run GenCaseした形状をParaViewで確認したところ、おかしな形状になった。
xmlを見ると、下のポイントのx,y座標が、(500,500)ではなく(0,0)のままだった。
xmlを手動で修正して、手動でGenCaseする必要がある。この場合、FreeCAD上のRun GenCaseではだめ。
cmdで手動で”~~\DualSPHysics_v4.4\bin\windows\GenCase4_win64.exe”を実行。
やり方は、RUN_DIRECTORYフォルダにあるexampleの.batファイルを参考に、自分でcmdで操作。
出力ファイルを作ったら、FreeCADのRUNでも解析できる。
(FreeCADに頼らない解析実行の勉強にはなる。)

より難解な形状を作ったり解析したりするには、DesignSPHysics使うしかない?バグあったけど、XMLの手入力はさすがに厳しい。形状が複雑になるにつれ、手作業でxmlを書くのは難しくなる。

ParaViewで解析結果の確認

上の動画はBlenderで作ったけど、結果の確認は基本ParaViewで行う。
ParaViewにvtkをロードして結果が確認できる。

10.png

vtkの準備

まずはvtkファイルの準備。これがないと結果が見れない。

  • 粒子で出力したいなら

    • PartVTKOut4でvtkを出力→点群が出てくる。
      • 連番のvtkが出てくる。
  • 液体っぽい表現で出力したいなら

    • IsoSurface4でvtkを出力→流体の表面を模したメッシュが出てくる。
      • Blenderのpythonスクリプトでmarching cube(後述)する必要はない。

ParaView

  1. ParaViewを起動

  2. VTKをすべて選択して、ParaViewにドラッグドロップする。

  3. 左側のメニューの目のマークを押して表示。
  4. 緑色の三角ボタン(Play)を押すとアニメーションされる。

paraviewは、 科学技術可視化 のためのソフトなので、レンダリングは充実していない。

paraview作ってるロスアラモス国立研究所って、アメリカのニューメキシコ州にあるマンハッタン計画で有名な研究所らしい。知らんかった。

FreeCADを使わずに計算を実行する例

DesignSPHysicsのフォルダ内に、examplesが入っている。←FreeCADを使わずに計算を実行している。

DualSPHysicsは、コマンドラインから実行するソフト。
examplesでは、コマンドライン作業をバッチファイルを使って実行させている。

examplesのフォルダに入っている.batファイルがバッチファイル。バッチファイルは、cmdの操作を記述できる。

バッチファイルで作業を簡略化

バッチファイルはダブルクリックして使う。テキストエディタで編集して、処理の順序を変えたりもできる。面倒なコマンドプロンプトの作業もダブルクリックだけで終わる。バッチファイルとは。 )

examplesの.batファイルは、プログラムファイルの場所が正しく書かれていないため、自分で書き直す必要がある。

.batファイルで、isosurfaceを作るには

例えば、isosurfaceの作成を、.batファイルで書くと

@echo off

rem "name" and "dirout" are named according to the testcase

set name=1
set dirout=%name%_out

set isosurface="X:\app\SPH\DualSPHysics_v4.4\bin\windows\IsoSurface4_win64.exe"

%isosurface% -dirin %dirout% -saveiso %dirout%/Surface

5行目のnameの部分は、解析ファイルごとに書き換えないといけない。

Blender+VisualSPHysicsを使った流体のレンダリング

VisualSPHysicsを使うと、連番vtkファイルを、アニメーションのシーケンスにしてくれる。

isosurfaceなら、そのままアニメーションレンダリングが可能。

点群のままならば、レンダリングされないので注意。


11.png12.gif

必要なもの

情報は随時更新されるはずなので、最新のものに従っていただくとして、このサイトを参考に必要なものをインストールする。

  • Blender(2.8以上 2.7なら導入に相当苦労する)
  • Python3.?(cmdでpythonと打った時に2.7が出る状態だと、導入できなかった。Blenderの中のPythonにPATHを通したらできた。)
  • Visual Studio 2015 redistributable package

  • visualSPHysics (エラーが起きて導入難しい時は、Forkされた古いほうを使う。)

    • Blenderが2.7 なら、最新版は無理。(Pythonが2.7でも無理かも?)
    • OpenGLのバージョン問題でBlender2.8が導入できなかった。

visualSPHysicsの使い方

Blenderのアドオンを無事導入できたら、

3D Viewの中にマウスカーソルを入れて、Shift+A>Mesh>DualSPHysics objectを押すとファイルダイアログが出現。

  • 左下に、操作できるGUIある。
  • 連番ファイルは、どれか一つを選ぶだけで全部選択してくれる。

アプライすると、アニメーションのシーケンスが下に出ているかと思います。

繰り返すけど、inportされたのが点群データだった場合、Blenderでレンダリングされない。(isosurfaceを出力しよう!)

Blenderのvertexとして結果をinportし、marching cubeでメッシュしてレンダリング

!苦行!なんでそんな事やるの!?isosurface4でmarching cube法で処理された結果出力できますけど?(マニュアル(DualSPHysics_v4.0_GUIDE.pdf)の11.6 Surface representationに書いてありました。)

それでもやりたい人向けに、書き残しておく。

marching cube法を試すには

マーチングキューブ法は、YouTubeの解説が充実している。

githubからmarching cubeのコードをダウンロードして使う。使い方ものってた。
しかし、これはmain関数を書き直す必要があった。理由は、点群データをメッシュする用途で書かれていなかったから。

書き直した部分だけ抜粋。

import mathutils
def main(a):
print("start calculation of isosurface")
bpy.context.scene.objects.active = bpy.data.objects[a]
obj = context.object
mesh = obj.data
size = len(mesh.vertices)
kd = mathutils.kdtree.KDTree(size)
for i, v in enumerate(mesh.vertices):
kd.insert(v.co, i)
kd.balance()
def scalarfield(pos):
co_find = (pos[0],pos[1],pos[2])
co, index, dist = kd.find(co_find)
return dist

p0=-0.1,-0.1,-0.1
p1=1.1,1.1,1
res=77#200
resolution=(res,res,res)
isolevel=0.02
start = time.time()
isosurface(p0,p1,resolution,isolevel,scalarfield)
elapsed = time.time()-start
print("end test %r"%elapsed)
del obj, kd


###########間は省略###############


if
__name__=="__main__":
# __
# ||
# ||
#_||_
#\ /
# \/ENTER OBJECT NAME.
a="ENTER_OBJECT_NAME"
bpy.ops.object.select_all(action='DESELECT')
bpy.data.objects[a].select = True
main(a)

かなり遅い。

結論

Blenderで液体っぽく表現したいなら、IsoSurface4を使いましょう。そうすれば、isosurfaceのvtkを出力してくれます。BlenderのPythonでやる必要はありません。これに気付かず、正月休みを無駄にしました。

(正月に記す。出すのは2月だが。)

posted by yuchan at 07:00 | Comment(0) | SPH

20200202

Calculix 10本ノック: 1本目

タグ:Windows WSL

Windows上で動かすLinux(WSL)で、Calculixを使えるようにした。
WSLでは基本GUIがないから、ひと手間必要だった。

Calculixとは

フリーの解析ソフト。
OSは基本的にLinux(Windows版もあるが、なぜかcgxの起動が遅い(僕だけ?))。
Abaqusと書式が似ているらしい。

  • ccx ソルバー
  • cgx プリポスト(プリとしてあまり使われない。プリにはgmeshかFreeCADがスタンダード?)

〇本ノックとは

例題を繰り返すやつ。Qiitaとかでよく見る。


ノック1本目: ccx, cgxのインストール

windows上でLinuxが動く状態にする

windows上でUbuntu(Linuxの有名ディストリビューション)が動く(いわゆる Windows Subsystem for Linux)状態にする。

ubuntu

とすると、

yusuke@yusuke-PC-2:~$

と表示されて、ubuntuが起動した状態になればOK。


ccx、cgxをインストールする

この動画を参考にインストールできる。要するに、

sudo apt-get install calculix-ccx
sudo apt-get install calculix-cgx

をする。

インストールが終わったら、

cgx -b a.fbd

上記コマンドで実行できると思いきや、

yusuke@yusuke-PC-2:~$ cgx -b aaa
on a Linux machine, nodename yusuke-PC-2, release 4.4.0-18362-Microsoft, version #476-Microsoft Fri Nov 01 16:53:00 PST 2019, machine x86_64
parameters:3 arguments:2
freeglut (cgx): failed to open display ‘’

エラーで何も動かない。

エラーの解決策

このエラーの解決策は、このサイトのウィンドウが出ないの箇所 )に載っていて、

  • Xming-6-9-0-31-setup.exe
  • Xming-fonts-x-x-x-x-setup.exe(x-x-x-xはバージョン)

をダウンロードしてインストール 後に、下記のコマンドで環境変数の設定して、

export DISPLAY=localhost:0.0

再度cgx

cgx -b aaa.fbd

とすると無事に、cgxが立ち上がる。(ちなみに、Windowsを起動するたびに、Xmingを立ち上げる必要がある。Xmingなしだとやはり同じエラー。後述のbash.bashrcへの記入(Linuxのviとかで)も必要。

1.png

yusuke@yusuke-PC-2:~$ cgx -b aaa.fbd
on a Linux machine, nodename yusuke-PC-2, release 4.4.0-18362-Microsoft, version #476-Microsoft Fri Nov 01 16:53:00 PST 2019, machine x86_64
parameters:3 arguments:2
GL_MAX_EVAL_ORDER:30
ERROR: The input file “aaa.fbd” could not be opened.

↑こんな表示が出る


(重要)exportは、永続性がない。

先ほどexportコマンドで、環境変数を設定した。
しかし、このままではexportは毎回入力しなければならない。
それを回避したいなら、

/etc/bash.bashrc

に、export文を追記する。 /etc/bash.bashrc は、毎回読み込まれるらしい。


サクラエディタで /etc/bash.bashrcを触ってはいけない

/etc/bash.bashrcの編集は、必ずLinuxのコマンドラインからviなどを使って行う。

上記フォルダを、windowsGUIからサクラエディタを使って編集してはいけない。

Windows上での場所が分かったからって、このLinuxファイルをWindowsのツールでいじってはいけない。

サクラエディタとかwindowsから編集すると、ubuntu起動時に、

-bash: /etc/bash.bashrc: Permission denied

と表示されて、bashが読み込まれなくなる。

ユーザー名の色も消える。

解決策は、

> ubuntu$ chmod 644 /etc/bash.bashrc

ここに載ってた。 https://tutorialmore.com/questions-114467.htm

 
Linuxのフォルダの場所

https://qiita.com/kalafinalice/items/70a76d35398ab11af778 にあったのだが、

/etc/bash.bashrcは、Windows上では

“C:\Users\yusuke\AppData\Local\Packages\CanonicalGroupLimited.UbuntuonWindows_79rhkp1fndgsc\LocalState\rootfs\etc\bash.bashrc”

にあった。yusukeの部分は各ユーザー名により違う。

 
ccx

ccxはGUI無いから、余計な事しないでも立ち上がる。

Windows上のファイルにアクセスするためには

例えばデスクトップにファイルを置いて

ccx -i "/mnt/c/Users/yusuke/Desktop/eyebar"

c:¥にアクセスするために、

/mnt/c

と書き換えなくてはならないのが結構面倒くさい。

解析実行に関しての詳細は、ノック2本目で書く。

タグ:Windows WSL
posted by yuchan at 19:00 | Comment(0) | Calculix

20200205

Calculix 10本ノック: 2本目

Calculixとは

フリーの解析ソフト。Abaqusと書式が似ているらしい。

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

ノック2本目: 解析の実行

データの入手

下記サイトのリンクから、beam.inpを取得する。

https://www.xsim.info/articles/CalculiX/Ex-Cantilever-Beam.html


WSLに入る

ubuntuを起動するために、windowsのcmdで

ubuntu

と入力。


解析の実行

当然windowsのフォルダに置いただろうから、デスクトップに置いたら

ccx -i "/mnt/c/Users/yusuke/Desktop/beam"

アドレスは”C:\Users\yusuke\Desktop\beam.inp”となる(yusukeの部分は人によって違う。私のPCの名前はyusukeなのでyusuke。)が、実行時は下記のように書き直す。

  1. c:¥ を /mnt/c/ に直す。
  2. 拡張子.inpを削る。

解析中

CalculiX Version 2.11, Copyright(C) 1998-2015 Guido Dhondt
CalculiX comes with ABSOLUTELY NO WARRANTY. This is free
software, and you are welcome to redistribute it under
certain conditions, see gpl.htm


You are using an executable made on So 31. Jul 13:26:31 CEST 2016

The numbers below are estimated upper bounds

number of:

nodes: 261

#

省略

#

Job finished

 
Linuxのフォルダの場所

https://qiita.com/kalafinalice/items/70a76d35398ab11af778

ccx -i "/mnt/c/Users/yusuke/Desktop/beam"

cgx -gx -read “/mnt/c/Users/yusuke/Desktop/beam.frd”

cd /mnt/c/Users/yusuke/Desktop

cgx -b “/mnt/c/Users/yusuke/Desktop/beam.frd”


解析結果を見る2つのスタイル

cgxで解析結果を見るために、二つのスタイルがある。

  1. 解析結果の.frdファイルをそのまま読む

    cgx -b "/mnt/c/Users/yusuke/Desktop/beam.frd"

    ccxで出力された解析結果ファイルをそのまま読む方法。

  1. 解析操作のスクリプト.fbdファイルから.frdファイルを読む。

    cgx -b "/mnt/c/Users/yusuke/Desktop/kaiseki2.fbd"

    cgxでどんな描画をするかをスクリプトを書いて指定する方法。

    スクリプト内で、1のfrdファイルを読み込んでいるので、どっちにしろ.frdファイルは必要になる。

    デスクトップに置いた、ミーゼスを見るためのfbdの例。

    read /mnt/c/Users/yusuke/Desktop/beam.frd new
    ds 2 e 7
    seta base f all
    view elem
    plot f all n
    plus fv base

    一行目にあるように、ファイルのアドレスは大事。

    ファイル名しか書かないと、カレントディレクトリしか探さない。

    まだそれほど手の込んだことを考えていないので、ミーゼスを表示できれば良しとする。

beam.gif
posted by yuchan at 07:00 | Comment(0) | Calculix

20200207

Calculix 10本ノック: 3本目: exampleの実行 CAD OnshapeTutorial

Calculixとは

フリーの解析ソフト。Abaqusと書式が似ているらしい。

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

環境

  • WSL

  • Xming

  • Calculix 2.11

  • gmsh

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

  • imagimagick

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

    sudo apt install imagemagick
  • pip2

  • matplotlib

ノック3本目: exampleの実行 CAD OnshapeTutorial

データの入手

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

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

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

WSLでLinuxのフォルダを使うというのに、windowsのデスクトップに解凍してしまった。

cd "/mnt/c/Users/yusuke/Desktop/CalculiX-Examples-2.11"

でyusukeの部分を各ユーザーで移動されたし。

デスクトップにcalculixサンプルのフォルダを置く.png

pythonのバージョンは2.7

pythonで実行させる場合もある。その際は、pythonのバージョンに気を付けられたし。

WSLならpython2もpython3も入ってる。
pythonもしくはpython2と打てば、python2が起動したし、
python3と打てば、python3が起動した。

.pyファイルのprintが

print ~

とスペース区切りで書いてあったら、python2。

CAD

OnshapeTutorial

ごく小さい段差がジオメトリ上にあって、オートメッシュで切ると細かいメッシュがたくさんできる。だけど、そんな所細かく切っても意味ない。そういう形状を、calculixのジオメトリ編集機能を使ってキレイに切りなおす。(自分でやったわけではない)

拘束は、足元の変位を拘束(1,2,3)。荷重は頭っぽい位置に面荷重。

hcpy_4.png

hcpy_5.png

gmshと連携して、STEPファイル→inpファイルにできる(inpファイルはcgxが読める)。みたいな理解。

cd ./CalculiX-Examples-2.11/CAD/OnshapeTutorial
cgx -b run.fbd

run.fbdの実行→gmshの終了ですべて終わる。run.fbdの内容は、

  • gmshでノード、エレメント、荷重位置、拘束位置を出力(解析に関する命令は入っていない。)
  • 出力した.inpを読んで、loadには、圧力1を設定して出力、supportと指定したところには何も指定しないでノード番号を出力
  • あらかじめ用意されていた.inpに、上記ファイルが読み込まれる設定になっており、この.inpが実行される。
  • 結果を読み込んで画像を保存して終わり。(要imagimagick)

.inpを作りこむことが勉強になる気がするんですけど、.inpはあらかじめ作ってある。
それどころか、run.fbdを実行すれば、何を書いてあるんだか分からなくても、絵が出てきてしまうという。

cgxでの荷重、拘束の表示の仕方

cgxでは、plot機能を使えば、fbdでLOADなどに指定したsetをplotする事で表示できる。
毎回plotしてたら、前回plotした情報が消えてしまうので、plotの次は、plusを使う。

FEMAPの様な、グラフィカルな表示は期待できない。
そもそも、荷重や拘束の設定は、.inp(もしくは.dlo)で行っているから、cgxの段階では、エクスポート前のsetに過ぎない。
このsetを表示する事で、荷重、拘束の表示とか言う事はできる。

オートメッシュにあまり興味がそそられないので、適当に飛ばす。
「後で必要になったら、.fbdを参考にすればいいや」としか思えない。

animationの機能に頼らずに、結果を動画にしてみた。うつらうつらしてるみたいになってしまったのは、基本ベクトルの何倍の係数をかけるかという設定のせいだと思う。

新規.gif



WSL上だと、datasetを選択すると

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

とか表示されてウィンドウが消えてしまう。
解決方法は、freeglutのバージョンを下げればいいらしいけど、なんか知らんが下げられなかった。

おとなしくコマンドから、

ds 2 e 7

(データセット 2 のエンティティ7を表示)とかして頑張って目的のデータを探した。

posted by yuchan at 00:00 | Comment(0) | Calculix

20200209

Calculix 10本ノック: 4本目:cgxで四則演算、while

Calculixとは

フリーの解析ソフト。

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

環境

  • WSL

  • Xming

  • Calculix 2.11

ノック4本目:cgxで四則演算、while

cgxのを使って、たし算ひき算わり算かけ算をやる。

cgxのコマンド実行には、

  • fbdファイルに下記を書き込んで実行
  • cgxのウィンドウをアクティブにして、マウスカーソルをウィンドウ内に入れてからコマンド入力

の二つの方法がある。

Calculixのcgx内で四則演算するには、valuを使って変数に入れてから、再びvaluを使って計算する。

valuを使った計算は、値をいったん変数に入っていなければ受け付けられなかった。
値をそのまま入力できなかった。
例えば、5 + 10 = 15は、

valu v1 5
valu v2 10

prnt v

valu v3 + v1 v2

prnt v v3

同様にして、

5-10=

valu v4 - v1 v2
valu v5 * v1 v2
valu v6 / v1 v2

While構文

valu a1 0
valu a2 10
valu stop 50

while a1 < stop
valu a1 + a1 a2
prnt v a1
endwhile

.fbdファイルに保存して、実行する場合の注意点
endwhileはテキストファイルの最終行になってはいけない。さもなければ、
key:ENDWHILE from string endwhile not known
と表示されて、whileされない。

prnt v

とすると、変数名すべてを表示する。

prnt v v1のやり方だと、正しく表示してくれないような…気のせいか?

でも、これで簡単な制御ができますね。

posted by yuchan at 07:00 | Comment(0) | Calculix

Calculix 10本ノック: 5本目: cgxを使ってモデル化→解析

Calculixとは

フリーの解析ソフト。

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

ノック5本目: cgxでモデル作成:梁

梁モデルを、cgxだけでモデル化する。コマンドは手打ちする。

今回はコマンドを頻繁に使う。
(重要)必ず、マウスポインターはGUIの中に入れてから、コマンド入力する。
そうでない場合、コマンドラインからの入力を受け付けてくれない。

テキストエディタにコマンドを書いておいて、

cgx -b ファイル名

をすれば、同じ内容になる。

Calculix(だけで行う解析)は、

  1. cgxによるモデル作成

  2. テキストエディタで.inpの作成

の二つが重要になってくると思ったので、やる。

点を打って、z10の所までスイープ→線になる。

pnt ! 0 0 0
seta center p all
swep center py tra 0 10 0

この状態では、何も表示されないので、

​ 点を表示したければ

plot p all

線を表示したければ

plot l all

線の名前も表示したければ、

plot la all

line_name.png

線の再定義

line L001 D001 D002 20

上記の代わりに、前述したswepの所で、

swep center py tra 0 10 0 20
と打っても同じ要素ができた。

要素タイプの選択

elty all be2

CalculiX GraphiX, Version 2.16のUser Manualに載っている要素を表にしてみた。

                 
梁beambe2be2rbe2fbe2dbe3be3rbe3f
平面三角triangletr3tr3utr3etr3str3ctr6tr6utr6etr6str6c
平面四角quadqu4qu4equ4squ4cqu4rqu4erqu4srqu4crqu8qu8equ8squ8cqu8rqu8erqu8srqu8cr
八面体hexahe8he8fhe8ihe8rhe20he20r
三角柱pentape6pe6fpe15pe15r
四面体tetrate4te4fte10te10mte10t

c: axisymmetric軸対象
e: plain strain平面ひずみ
s: plain stress平面応力
u: unstructured mesh訳がすぐに出てこなかった。(これ?),
r: reduced integration低減積分
i: incompatible modes
f: fluid element for ccx,
t: initial temperatures are interpolated linearly within the tet element (ccx:C3D10T)).

この状態では、まだエレメントは作成されていない。

mesh all

と書いて初めてエレメントができる。

ノード、エレメントの表示

エレメントができたので、ノードもできたはず。
いちいちコマンドを入力しなければ、表示されない。

plotコマンドで表示しようとすると、以前のplot結果が消えてしまうため、plusコマンドを使う。

plus n all
plus e all

node.png

出力

.inpに読ませるメッシュファイル(all.mshファイル)の出力のやり方は、

send all abq

これで、all.mshが、カレントフォルダに出力される。

all.mshには下記が書かれている。

  • ノード
    • ノードのセット名
    • ノード座標
  • エレメント
    • エレメントのエレメントタイプ
    • エレメントセット名
    • エレメントが使用しているノード番号
*NODE, NSET=Nall
1,0.000000000000e+00,0.000000000000e+00,0.000000000000e+00
2,0.000000000000e+00,0.000000000000e+00,1.000000000000e+01
3,0.000000000000e+00,0.000000000000e+00,2.000000000000e+01
4,0.000000000000e+00,0.000000000000e+00,3.000000000000e+01
5,0.000000000000e+00,0.000000000000e+00,4.000000000000e+01
6,0.000000000000e+00,0.000000000000e+00,5.000000000000e+01
7,0.000000000000e+00,0.000000000000e+00,6.000000000000e+01
8,0.000000000000e+00,0.000000000000e+00,7.000000000000e+01
9,0.000000000000e+00,0.000000000000e+00,8.000000000000e+01
10,0.000000000000e+00,0.000000000000e+00,9.000000000000e+01
11,0.000000000000e+00,0.000000000000e+00,1.000000000000e+02
*ELEMENT, TYPE=B31, ELSET=Eall
1, 1, 2
2, 2, 3
3, 3, 4
4, 4, 5
5, 5, 6
6, 6, 7
7, 7, 8
8, 8, 9
9, 9, 10
10, 10, 11



全部まとめると下記のようになった。

pnt ! 0 0 0
seta center p all

swep center py tra 0 0 10 10

##line L001 D001 D002 10

elty all be2

mesh all
send all abq

inpファイルの作成。

メッシュができたので、このまま.inpの編集をする。
全部手打ち。
サンプルファイルのコピペもありだと思う。

先ほどのcgx作業で完成したall.mshを読み込むために、下記を記入。

*INCLUDE, INPUT=all.msh

カレントフォルダになければ、アドレスを入力するしかない。

境界条件は下記で、「ノード1を1から6まで完全固定」になると思う。

*BOUNDARY1 ,1, 6

材質は、名前がSTEELで、ヤング率210000、ポアソン比0.3。

*MATERIAL,NAME=STEEL*ELASTIC210000,0.3

梁要素の断面特性は、断面を適用する要素のセット名をEall(all.mshファイル内で、定義されているエレメントのセット名)、材料は上記で設定したSTEEL、断面形状は矩形。

*BEAM SECTION,ELSET=Eall,MATERIAL=STEEL,SECTION=Rect15,151, 0, 0

解析条件は、

*step

*end step

の間に書く。

静解析の場合、書くことは

  • 静解析の指定
  • 荷重

  • 出力要求(ノード)

  • 出力要求(エレメント)

静解析の指定は下記。

*static

集中荷重は下記で、「ノード11を、1方向(X方向)へ、1000の荷重を付加」

*CLOAD11, 1, 1000.0

出力ファイルへの出力要求は、ノードの結果の場合下記で、「U (変位)」

*NODE FILEU

出力ファイルへの出力要求は、エレメントの結果の場合下記で、「S (応力)」

*el file, output=1DS

まとめるとこうなった。

*INCLUDE, INPUT=all.msh

*BOUNDARY
1 ,1,6

*MATERIAL,NAME=STEEL
*ELASTIC
210000,0.3

*BEAM SECTION,ELSET=Eall,MATERIAL=STEEL,SECTION=Rect
15,15
1., 0, 0

*step

*static

*CLOAD
11,1,1000.0

*node file
U
*el file, output=1D
S

*end step

上記をsolve.inpと名前を付けたら、下記コマンドで解析実行。

ccx solve

job finishedと出たら、

cgx -read solv.frd

マウスポインターをGUIの中に入れてから、
全変位を見たければ、下記コマンド入力する。

ds 1 e 4

もしくは、ミーゼス応力を見たければ

ds 2 e 7

hcpy_1.png

う〜ん不便。cgxでのメッシュ作成、ccx用の.inpファイル作成、全部不便。

GUIって便利だなと改めて思った。

posted by yuchan at 19:00 | Comment(0) | Calculix

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.

上記のようなエラーで。

posted by yuchan at 19:00 | Comment(0) | Calculix

20200216

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

line_name.png

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

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

plot fa all

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

12.gif

.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つのノード同士で式を作っていた。ノードがどうやって選ばれたのかは知らない。

node_all.png

結果

固有値の結果は、.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は問題なく動くと思う。

posted by yuchan at 07:00 | Comment(0) | Calculix

20200220

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を入れ替えたらうまくいった。

model.png

model2.png

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は下記。

ind.png


stopnegは下記。

 topneg.png

結果

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

result.png

posted by yuchan at 07:00 | Comment(0) | Calculix

20200223

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を使わないと、全体が見えなくなった。

frame.png

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

line ! D002 D003 D001 14

1/4.png

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)を使っていた。

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

マテリアルを指定して、

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

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

*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については以前やったので割愛。
whileendwhileで終わるが、endwhileがテキストエディタの最終行に来たらエラーで止まるのは要注意。

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

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

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

read modal.frd new
movi frames 30
ds 1 a
posted by yuchan at 07:00 | Comment(0) | Calculix