CalculiX CrunchiX USER'S MANUAL version 2.5 日本語訳
ver1.0 今井康貴, 有尾元基（佐賀大学海洋エネルギー研究センター）
5章

簡単な例題

5.1 片持ち梁（beamp.inp）
この節では、自由端で荷重される片持ち梁を分析する。
片持ち梁の形状、荷重および境界条件を図 1 に示す。梁の寸法は 1×1× 8m である。 9 ×10 N36の点荷重を受ける。左端は全方向に固定する。長さと力の単位は 1m および 1MN をとる。図 2
は CalculiXGraphiX(cgx) で梁の形状を作りメッシュ分割した結果である。図では要素ラベルだけを表示している。
CalculiX 入力は、基本的に、形状と境界条件を記述するモデル定義部、荷重を定義する 1 つ
以上のステップから構成される。
モデル定義部はファイルのはじめから最初の*STEP までである。全ての入力は * からはじまるキーワードにより行われる。キーワードは入力データの種類を示す。* STEP はキーワードの一種である。ほとんどのキーワードは、モデル定義（最初の *STEP 以前）もしくはステップ（*STEP と *END STEP の間）のどちらかで使用される。両方で使用されるものは少数である。
図 4 の例では、最初のキーワードは *HEADING であり、その後に問題内容が記述される。これは出力に対する影響はなく、単に識別に役立つ。次に *NODE と座標系が記述される。1 行のデータはコンマで区切る。132 文字を超えないこと。キーワードは必要に応じて繰り返すことができる。例えば、*NODE キーワードは各節点の前に書ける。
次に、*ELEMENT によって形状が定義される。形状の定義は要素のリストであり、要素タイ
プ、要素に属する節点、要素オーダーが記述される。要素タイプはパラメタである。
この梁のケースでは 20 節点の立方体要素（簡易積分型）C3D20R を使用した。

図 1: 梁問題における形状と境界条件1図 2: 梁の要素分割

Model Definition

Material description

*STEP
Step1
*END STEP

2

*STEP
Step2*END STEP

・
・
・
・
*STEP
Step n*END STEP

図 3: CalculiX 入力構造3

図 4: 入力ファイル

ELSET=Eall によって、*ELEMENT に続く要素はすべて、Eall セットに格納される。
このセットは後で材質定義に使用される。それぞれの要素は、20 個の節点番号とともに記述される。*NODE および *ELEMENT が形状が記述の最重要部である。残ったモデル定義は幾何学
的境界条件および材質の記述である。
この梁問題の幾何学的境界条件は z=0 の固定である。この目的のため、z=0 にある節点を*NSET によって節点群 FIX に格納する。群に属する節点は、*NSET の次の行に書く。
*BOUNDARY によって FIX に属する節点は 1，2，3 方向（x,y,z に相当する）に固定される。
図 4 の 3 つの*BOUNDARY は一緒に書いても良い。
4

*BOUNDARY
FIX,1
FIX,2
FIX,3

あるいはさらに短くできる*BOUNDARY* 境界
FIX,1,3
これは自由度 1～3 の固定を示す（ゼロがセットされる）。

次は材質の記述である。この節は特別である。というのは、同じ材質はグループ化される必要があり、記述は*STEP 前ならどこでも良いためである。
材質の記述は常に*MATERIAL で始まる。パラメタ NAME により材質の名前を定義できる。
2

ここでは材質は線形弾性、ヤング率 210, 000.0 MN / m 、ポアソン比 0.3 の鋼とする。これら特性は *ELASTIC の次の行に記述する。これで材質の定義は終了。次に、*SOLID SECTION により材質が要素群 Eall に割り当てられる。
モデル定義部の最終部 *NSET,NSET=LOAD で荷重をかける節点群 LOAD を定義する。モデル定義部と*STEP 間の** から始まる行はコメントである。コメントはどんな場所にも導入できる。CalculiX はコメントを無視する。コメントは入力ファイルを分かりやすくするのに役立つ。
この問題では、1 ステップだけで良い。ステップは常に *STEP で始まり、*END STEP で終わる。*STATIC は手続きを定義する。*STATIC は、準静的に荷重が載荷されることを示す(質量慣性を考えない)。他には *FREQUENCY、*BUCKLE、*MODAL DYNAMIC、*STEADYSTATE DYNAMICS および *DYNAMIC がある。次に、*CLOAD で節点群 LOAD に集中荷重が適用される。y 方向に大きさ 1 の荷重が作用する。全荷重は 9 である。
最後に、ユーザー指定項目の表示およびファイル記録を行う。*NODE PRINT と *EL PRINT
は .dat ファイルにアスキー形式で記録される PRINT が無い場合、.dat ファイルは作成されない。*NODE PRINT と *EL PRINT の後には出力する節点群および要素群を記述する。要素情報は積分点に格納される。
一方、*NODE FILE および *EL FILE は、.frd アスキーファイル出力を管理する。T このフ
ァイルにある結果は CalculiX GraphiX (cgx) で見ることができる。*NODE FILE および *EL
FILE で選択された項目は、完全なモデルのために常に格納される。要素量は節点に外挿され、同一節点への複数の寄与は平均される。*NODE PRINT、*EL PRINT、*NODE FILE、*EL
FILE の項目は文字によって選択される。例えば、U は変位、S は(コーシー)応力である。
5

入力は *END STEP で終わる。
梁（beam）問題の出力は beam.dat および beam.frd である。beam.dat は節点群 Nall の変位
および Eall 積分点の応力が記録される。beam.frd は全節点の変位および外挿された応力を記録する。これは可視化プログラム CalculiXGraphiX(cgx) の入力になる。cgx 結果を図 5、6、7 に示す。

図 5: 梁の変形6図 6:梁の軸垂直応力7Figure 7: 梁中の VonMises 応力
図 5 は荷重下における梁の変形を示す。予想通りに、梁は端部の横力により曲がる。図 6 は軸方向の垂直応力を示す。曲げモーメントにより、梁高さ方向の線形分布が示される。図 7 は梁内部の VonMises 応力を示す。

5.2 圧縮力よる梁の振動周波数（beamf2.inp）
前節からの梁を考慮し、その固有振動数と固有振動モードを決定する。横方向に違った振動周
波数を持たせるため、断面形状を 1× 1 から 1× 1.5 に変える。全長は変わらず、 8 のままとする。
入力を図 8 に示す。前節の入力とほぼ同じである。入力ファイルは beamf2.inp である。
違いはステップにある。第一ステップにおいて 梁の端部に圧縮力が適用される。節点群 LASTに属する節点には、z 方向に −48.155 の荷重が適用される。これは z 負方向に 48.155 の力がかかることに等しい。第二ステップは周波数である。PURTABATION により、前段階の静的ステップの変形および応力を周波数計算に取り入れる。*FREQUENCY およびそれに続く行が、モード解析ステップであることを示し、固有周波数の低い方から 10 個を決定する。固有周波数8

は .dat ファイルに保存される。表 2 は、前負荷の有無による固有周波数（cycle/s）の比較を示す。また ABAQUS の結果も示す。前負荷無しの入力は test フォルダ中の beamf.inp である。
前負荷により固有周波数が低下する。これは、低周波数において顕著である。実際には、最低の曲げ振動数は非常に低いため、座屈が生じるだろう。事実、座屈を決定する１つの方法は、最低の固有周波数がゼロになる点まで圧縮荷重を増加させることである。この例では、座屈荷重は21× 48.155 = 1011.3 であることを意味する。係数 21 は同じ荷重が 21 節点に作用することを表す。もう一つの座屈荷重を決定する方法は *BUCKLE の使用である。この入力ファイルは
beamb.inp である。

図 8: 入力ファイル

9

表 2:
Without preload
CalculiXABAQUS
13,096.
13,096.
19,320.
19,319.
76,840.
76,834.
86,955.
86,954.
105,964.
105,956.
162,999.
162,998.
197,645.
197,540.
256,161.
256,029.
261,140.
261,086.
351,862.
351,197.

With preload
CalculiXABAQUS
705.
1,780.
14,614.
14,822.
69,731.
70,411.
86,544.
86,870.
101,291.
102,148.
162,209.
163,668.
191,581.
193,065.
251,858.
253,603.
259,905.
260,837.
345,729.
347,688.

慣性短軸に対する第 2 曲げモードを図 9、第 1 ねじれモードを図 10 に示す。

図 9: 慣性短軸に対する第 2 曲げモード10

図 10: 慣性短軸に対する第 1 ねじれモード
5.3 細い軸をもつ回転円板の振動計算（rotor.inp）
これは複雑な振動計算の例である。外径 10 、内径 2 、厚さ 0.25 の円板が、外径 2 、内径 1 の中空軸に取り付けられる。入力ファイルは rotor.inp である。軸の中央に円板を取付ける。軸の両端は完全固定する。円板から軸端までの長さは 50 である。この例の入力ファイルを図 11 に示す。
ファイルには節点と要素が定義される。節点群 Nfix は軸端の節点であり全方向に固定される。
材質は一般鋼である。遠心荷重計算のために密度が必要であることに注目する。
円板が回転しているため、遠心力として予荷重がかかる。したがって、第一段階は、この荷重による変形および応力を計算するための非線形幾何静的なステップとである。frequency ステッ
プ中のパラメタ perturbation により、この予荷重は剛性行列の周波数計算に取り込まれる。結果の固有周波数は rotor.dat ファイルのはじめに保存される。図 12 は 9000 回転の場合である。
*FREQUENCY ステップでは固有値問題が解かれ、固有値は構造物の固有周波数の 2 乗である。
図 12 の第 1 列が固有値、2～4 列が固有周波数である。固有値が負の場合、虚数の固有周波数となる。 9000 rad / s で回転するローターの２つの最低固有値の場合である。約 6000 rad / s 以下の回転速度では、全ての固有振動数が実数である。 18000 rad / s までの回転速度の関数としての11

最低の固有振動数が図 13 に示される(+と x のカーブ)。
図 11: ローターの入力（省略。rotor.inp 参照）
EIGENVALUE OUTPUT
MODE NOEIGENVALUEFREQUENCY

REAL PART(RAD/TIME)IMAGINARY PART(CYCLES/TIME(RAD/TIME)

1 -0.2744426E+09 0.0000000E+00 0.0000000E+00 0.1656631E+05
2 -0.2744426E+09 0.0000000E+00 0.0000000E+00 0.1656631E+05
3 0.1775557E+06 0.4213736E+03 0.6706369E+02 0.0000000E+00
4 0.1775558E+06 0.4213737E+03 0.6706370E+02 0.0000000E+00
5 0.7206542E+09 0.2684500E+05 0.4272515E+04 0.0000000E+00
6 0.7206542E+09 0.2684500E+05 0.4272515E+04 0.0000000E+00
7 0.2028445E+10 0.4503826E+05 0.7168061E+04 0.0000000E+00
8 0.2740625E+10 0.5235098E+05 0.8331917E+04 0.0000000E+00
9 0.2740625E+10 0.5235098E+05 0.8331917E+04 0.0000000E+00
10 0.5597315E+10 0.7481521E+05 0.1190721E+05 0.0000000E+00
EIGENVALUE OUTPUT
MODE NOFREQUENCY
REAL PARTIMAGINARY PART

(RAD/TIME) (CYCLES/TIME)(RAD/TIME)

1 0.5089832E+01 0.8100720E+00 -0.1026829E-10
2 0.1173764E+05 0.1868104E+04 0.2672643E-03
3 0.1453582E+05 0.2313447E+04 0.8041058E-05
4 0.2338142E+05 0.3721269E+04 -0.1226430E-04
5 0.3485137E+05 0.5546768E+04 -0.1510576E-03
6 0.3836838E+05 0.6106517E+04 -0.1763646E-04
7 0.4503827E+05 0.7168064E+04 0.7438635E-04
8 0.4957786E+05 0.7890561E+04 -0.1156950E-02
9 0.6937749E+05 0.1104177E+05 -0.1043579E-03
10 0.7710084E+05 0.1227098E+05 -0.2964534E-03
図 12: ローターの固有振動数（rotor.dat 参照）12

図 13: 回転速度と固有振動数

虚数の固有振動数の物理的な意味は何か？振動計算に起因する固有モードは e

iω t

を含む。

固有振動数 ω が実数の場合、サインもしくはコサインである。 ω が虚数の場合、増加もしくは減少する指数関数である。[16] したがって虚数の固有振動数は、反応はもはや振動しない。
それは無限に増加し、システムは破壊される。図 13 を見ると、最低次の固有周波数は減少し、回転速度 6000 rad / s 近くでゼロになる。その点において固有振動数は虚数になり、ローターは破損する。これは長い間エンジニアを当惑させた。というのは、実システムは破損なしに臨界超過の速度に達するためである。
本質的な点は、計算が（軸に固定された）回転移動座標系で行われることである。ニュートン則は加速する座標系において有効ではない。そして回転座標系は加速している。ニュートン則の修正項はコリオリ力として与えられる。コリオリの力の導入は複素非線形固有値系になる。これは*COMPLEX FREQUENCY で解かれる（6.8.3 節参照）。結果の固有周波数は実数であるが、固有振動モードは通常複素数であることを示すことは可能である。これは回転固有振動モードになる。

13

図 14: 2 節の固有モード*COMPLEX FREQUENCY を使用するためには、*FREQUENCY ステップにおいて Coriolis
無しの固有振動モードを計算し保存しておく必要がある(STORAGE=YES)。図 11 参照のこと。
複素周波数応答はこれら固有振動モードの線形結合で表される。*COMPLEX FREQUENCY ス
テップで要求される固有振動数の個数は、それ以前の*FREQUENCY ステップにおける固有振動数の個数を超過するべきでない。固有振動モードが複素数であるため、それらは振幅と位相で保存される。*NODE FILE の下に PU と記述すること。
回転軸の正確な固有値は図 13 の直線になる。それぞれの線は固有モードを表す。最下端の減少する線は、-Z 方向を見た時の 2 節反時計（CCW）固有モードである。最上端の減少する線は
3 節反時計固有モード、最下端および最上端の増加する線は 2 節時計固有モードである。節は半径方向の運動がゼロになる位置である。図 14 は 2 節固有モードを示す。図 15 は 3 節固有モードを示す。図 13 において、x と y のスケールが同一であれば、線の傾きは 45 度になる。
両方の増加する直線が全く同一の固有モードに相当することは驚きである。例えば、回転速度

5816 rad / s では、同一の固有モードが回転速度ゼロおよび 11632 rad / s の固有振動数において生じる。しかし、固有振動モードが回転系で計算されることを思い出すこと。つまり、軸とともに回転する観察者によって観察されたものである。固定観察者の振動を得るために、周波数は原
点 を 通 り グ ラ フ を 二 分 す る 45 度 の 直 線 に 対 し て 考 慮 さ れ な け れ ば な ら な い 。 観 察 者 は5816 rad / s および −5816 rad / s で同一の固有モード（cw と ccw）を見る
14

最後に、コリオリ力は必ずしも適切だとは限らない。一般に、細長い回転構造物（大きなブレード...）はコリオリによる重要な振動数シフトを示す。

図 15: 3 節の固有モード

5.4 炉の熱計算（furnace.inp）
この問題は、図 16 に描かれた炉の熱計算を扱う。炉は温度 Tb の底板から構成され規定される。
非常に短時間内に 300K から 1000K まで直線的に変わる、その後は一定のままである。炉の側壁は外界から分離されているが、炉の他の壁の放射によって熱交換を行う。側壁と底の放射率は

ε = 1 である。炉の天板は放射によって熱交換する。周囲温度は 300K である。
天板の放射は ε = 0.8 である。更に、天板は一定速度 0.001 kg / s で移動する流体(空気)は対流によって熱交換が行われる。上右角の流体温度は 300K である。オーブンの壁は 10 cm の鋼で作られ て い る 。 鉄 鋼 株 用 の 材 質 定 数 は 次 の と お り で あ る 。 熱 伝 導 度 k = 50 W / mK 、 比 熱
c = 446 W / kgK お よ び 密 度 ρ = 7800 kg / m3 。 空 気 の 材 質 定 数 は 次 の と お り で あ る 。 熱C p = 1000 W / kgK お よ び密 度 ρ = 1 kg / m3 。 対流 係数 h = 25W / m2 K で あ る。 炉の 寸法 は0.3 × 0.3 × 0.3m3 (立方体)である。 t = 0 では全部 T = 300 K である。位置 A～E における温度を
15

時間の関数として求める。
入力を図 17 に示す。これは節点の定義から始まる。構造で最も大きな節点番号は 602 である。
節点 603～608 は流体の節点である。つまり流体節点が定義された。z=0.3 は炉の天井、z=0 は底である。節点 603 は流体温度が 300 K の点（入口）、節点 608 は出口に相当する。他の点はそれらの中間に位置する。

図 16: 炉の記述

流体の節点の座標は実際の計算に入らない。*FILM の対流の定義だけが炉と流体の間の熱交換を支配する。*ELEMENT を用いて炉壁を 6 節点シェル要素で定義する。更に、流体の節点も要素(要素タイプ D)、いわゆるネットワーク要素に割り当てる。これらの要素は流体に材質特性を割り当てるために必要である。通常、材質特性は節点ではなく要素に割り当てられる。ネット
ワーク要素は 2 つの端部節点から構成される。温度は未知数である。中央要素で要素を通過する
質量流量を定義する。603 から 613 の流体の節点は要素 301 から 305 のネットワーク要素に割り当てられる。
次に、2 つの節点群を定義する。GAS は全ての流体節点、Ndown は、炉の底の全ての節点である。*PHYSICAL CONSTANTS は、放射の分析に必要である。それは絶対零度を定義する、ケルビンとステファンボルツマン定数のためゼロ。入力は SI 単位を使用する。
次に、STEEL の材質定数を定義する。熱分析であるから熱伝導率、比熱および密度を定義する。*SHELL SECTION は 要素群 FURNACE に鋼材質を割り当てる。要素群 FURNACE は以前に *ELEMENT によって定義される。それは、炉に属する全ての要素を含む。厚さ 0.01m
16が割り当てられる。

17

図 17: 炉の入力ファイル
18

GAS の材質定数は密度と比熱である。これらは流体用の定数である。流体中の伝導は考慮されない。すべてのネットワーク要素を含む要素群 EGAS に材質 GAS を割り当てる。*INITIAL
CONDITIONS は全節点（炉の節点および流体の節点）に初期温度 300 K を割り当てる。
*AMPLITUDE はランプ関数（傾斜関数）を定義する。 0.0 で 0.3 で開始し、 1.0 において 1.0 となるまで直線的に増加する。これは炉の底における温度境界条件の定義に使用される。これでモデル定義は終了である。
最初のステップは t = 0 から t = 1 の温度境界条件の線形増加を記述する。*STEP の INC = 100
はこのステップにおいて 100 段階の増分を考える。その手続きは *HEAT TRANSFER である、
つまり、我々は純粋な熱分析を行う。未知数は温度である。また、機械的な未知数(例えば変位)はない。ステップ時間は 1 である。最初の増分量は 0.1 である。両方を *HEAT TRANSFER の次の行に記述する。*HEAT TRANSFER にはパラメタ STEADY STATE が無く、これが非定常解析であることを示す。
次に、温度境界条件が来る。炉の底は 1000K を維持するが、振幅 A1 で変更される。温度境
界条件は 0.3 × 1000 = 300 K で開始し、 t = 1 において 1000 K になるまで直線的に増加する。第2 の境界条件は、節点 603 の(流体)温度が 300 K を維持することを指定する。これは入口温度である。
「11」が温度自由度であることに注目する。
流体の質量流量は *BOUNDARY で定義される。ネットワーク要素の中央節点の最初の自由度に適用される。第 1 行は節点 609 の質量流量が 0.001 であることを示す。節点 609 はネットワー
ク要素 301 の中央節点である。質量流量が正であるため、流体は節点 603 から節点 604 へ流れる。
つまり、ネットワー要素 301 の第 1 節点から第 3 節点に流れる。ユーザーは質量保存を確認する。これはプログラムでもチェックされる。
放射境界条件の最初のセットは炉底の上面は放射率 1 および環境温度 1000K で空洞放射を指定する。ある位置での到達率が総計 1 にならない場合、空洞放射は環境温度が使用される。1 に不足している部分は周囲環境へ放射される。各行の 1 番目は要素番号、2 番目は放射面である。
一般に、これらの行は、cgx (CalculiXGraphiX) で自動的に生成される。
2 番目、3 番目の部分は炉の天井と側壁の内部放射を定義する。4 番目は炉の天井の上面から環境への放射を定義する。 300 K を維持する。天井の放射率は 0.8 である。
次にフィルム条件が記述される。炉の天井の上面における強制対流を対流係数 h = 25W / mKで定義する。*FILM の次行は、要素 51 の第 2 面が節点 604 を通じて強制対流することを示す。

19

図 18: t=3001s の温度分布

図 19: 選択された位置における温度20

行の最後は対流係数である。したがって、流体と熱交換する面は、熱交換が行われる適切な流体
節点を指定する必要がある。
最後に、*NODE FILE は温度を .frd ファイルに保存し、*NODE PRINT は流体温度を .dat
ファイルに保存する。
入力ファイルは test フォルダの furnace.inp である。この分析では、第 2 ステップを追加した。
残り 3000 秒間炉底温度を一定に維持する。
計算の間に何が起こるか？炉の壁面と天井は壁中の熱伝導および底面からの放射により加熱される。炉の天井は周囲への放射および流体との対流によって熱を失う。流体との熱交換により温度分布は非対称である。入口において流体は低温であるため、炉は出口より入口で多くの熱を失う。出口において流体温度は高くなり、流体と炉の温度差は小さくなる。対流により、入口から出口に向かって温度が上昇すると予想される。熱伝導により、天井の中央において温度が最低になると予想される。両者の効果が重なる。 t = 3001 s の温度分布を図 18 に示す。炉の底から天井に向かって温度変化がある。天井では温度は対称ではない。これは図 19 にも見られる。A～Eの温度が時間の関数で表示される。
定常状態条件に到達していないことに注目する。また、シェル要素のような 2D 要素は、自動的に適正な厚さをもつ 3D 要素に拡張されることに注意する。したがって、絵はは 3D 要素を示す。これは CalculiX GraphiX で描かれた。

5.5 ダム下の浸出 （dam.inp）
この節では、ダム下の地下水流れを分析する。図 20 にダムの形状を示す。これは文献[22]の 1
章の練習問題 30 である。長さの単位はフィート（ 1 ft = 0.3048 m ）である。ダムの水位は上流側 20 フ ィー ト、下 流側 5 フィ ートで ある。 ダム下の 土は異 方性で ある。上 流の浸 透率 は

k1 = 4k2 = 10 −2 cm / s である。下流は 25k3 = 100k4 = 10 −2 cm / s である。興味あるのは水力勾配である。これはパイピングが発生するか否かを予測するためである。パイピングは、土が地下水によって持ち去られることであり、不安定な状態を構成する。（通常ダム下流側に発生する）経験則として、水力勾配が 1 の場合にパイピングが生じる。
6.8.14 節は静止している地下水流れの支配方程式が熱伝導方程式と同じことを示す。水のヘッドに相当する変数は温度である。また、流速は熱の流れに相当する。有限要素分析法では SI 単位を採用した。フィートをメートルに変換した。

21

図 20: ダムの形状
更に、上流および無限下流に水を通さない垂直壁を想定した。（実際にはダム中点から 30 m 上流
および 30 m 下流）

境界条件は以下のとおりである。
1.ダム上流および下流の左側および右側の鉛直境界、底面は不透生性である。これら境界において法線方向の水流速がゼロである。これは熱フラックスがゼロに相当する。

2. ヘッドを定義するため、ｚ軸を底面から上向きにとる（ヘッドの定義は式 108 参照）。大気圧

P0 はゼロと仮定する。上流側の全ヘッドは 28 フィート、下流側の全ヘッドは 13 フィートである。
熱問題では、これは温度境界条件に相当する。

入力ファイルの要約を図 21 に示す。完全な入力ファイルは dam.inp である。問題は 2 次元であり、CalculiX GraphiX によるメッシュでは qu8 要素が用いられた。ダムに近接する部分を高解像度でメッシュ分割するため、バイアスが使用された。図 22 にメッシュ分割を示す。
入力ファイルの最初に節点および要素を定義する。CalculiX GraphiX の qu8 要素は、send コ
マンドにより、デフォルトで CalculiX CrunchiX の S8 要素に変換される。しかし、ここでは平面要素がより適切である。問題の計算が熱であり機械ではないため、平面歪み要素(CPE8)ある
いは平面応力要素(CPS8)をとるかは些細な問題である。*ELSET で 2 種類の異なる土の要素セットを定義する。全ヘッド一定が適用される節点は *NSET で定義される。土の浸透性は熱問題の熱伝導係数に相当する。浸透性は *CONDUCTIVITY(TYPE=ORTHO) を使用して直交的と定22

義される。次行は、x、y および z 方向の浸透性である(SI 単位: m/s)。z-方向の値は些細である。
というのはその方向に勾配が期待されないためである。SOLID SECTION で土領域に材質特性を割り当てる。計算は静的なため、*INITIAL CONDITIONS は不要である。しかし CalculiX
CrunchiX は形式的に伝熱計算でそれを必要とする。

図 21: 入力23

図 22: 全揚程24

図 23: y 方向の浸出速度*HEAT TRANSFER ステップにおいて STEADY STATE 計算を指定する。追加のステップ情報は無しである。これはステップ長のおよび初期増分にデフォルト値(1)を使用することを意味
する*BOUNDARY は上流および下流の全ヘッドを定義する。熱の自由度は 11 である。最後に*NODE PRINT, *NODE FILE, *EL FILE で出力を指定する。NT は温度を表し、この問題ではそれに等価な全ヘッド（図 22）である。HFL は熱フラックスを表し、この問題ではそれに等価な地下水速度である。図 23 の y 方向速度成分である。
上流は浸透性が高いため全ヘッド勾配は小さい。下流は逆である。流速は特に下流に重要であ−4る 。 流 速 は図 23 の 赤い 部 分 に おい て 2.25 ×10 m / s に 達 する 。 下 流の y 方 向 浸透性 は

4 ×10−4 m / s であるため、これは水力勾配 0.56 に相当する。これは 1 未満であるため、パイピ
ングが生じる。ダムのすぐ隣りにおいて流速が最大である。
この例では、CalculiX GraphiX 伝熱能力を用いて浸出問題を解決できることを示した。同じことはラプラス型方程式で支配された他の現象に当てはまる。

5.6 水配管システム（pipe.inp）
CalculiX では、定常状態の水力学（油圧）・空気力学ネットワーク計算が可能である。単独でも、あるいは構造の機械、熱の計算と一緒に実施可能である。ここでは[10]にある配管を分析する。入力ファイルは pipe.inp である。
ネットワークを図 24 に示す。ネットワークは以下で構成される。
•

表面水位 14.5 m の上流側の貯水槽•

A: 収縮 0.8 の入口。

•

AB: 直径 0.2 m 、全長 5 m のパイプ。

•

B: 半径 0.3 m 、 45 度のベンド•

BC: 直径 0.2 m 、長さ 5 m のパイプ。

•

CD: 直径 0.3 m 、長さ 5 m のパイプ。

•

DE: 直径 0.15 m 、長さ 2.5 m のパイプ。

•

E: α = 0.5 のゲートバルブ。

•

直径 0.15 m 、長さ 1.56 m のパイプ。

•

水位 6.5 m の貯水槽

25

図 24:
全てのパイプは摩擦係数 n = 0.015 である。入力を図 25 に示す。
CalculiX では、線ネットワークは 3 節点ネットワーク要素(D タイプ要素)でモデル化される。
要素の角節点において温度と圧力は未知である。それらはそれぞれ自由度 0 および 2 に割り当てられる。中央節点では、質量フラックスは未知数で、自由度 1 に割り当てられる。ネットワーク
要素の特性は *FLUID SECTION で定義される。それらは 6.3 節（気体）および 6.4 節（液体）で扱われる。このネットワークでは以下が必要である。

•

液体を表現するダミーのネットワーク入口要素はネットワーク(element1)を入力している。それは節点番号 0 として第 1 節点と見なされる。

•

地点 A における PIPE ENTRANCE ネットワーク要素。この要素は水深を考慮に入れる。
特別の貯水槽要素がないことに注目すること。水位の差は要素の角節点に適切な座標を割り付けることにより、任意の要素タイプ中に考慮することができる。

•

A-B 間を表す PIPE MANNING 型要素 (要素 3)•

地点 B を表す PIPE BEND 型要素（要素 4）•

B-C 間を表す PIPE MANNING 型要素（要素 5）•

地点 C で直径増加を表す PIPE ENLARGEMENT 型要素（要素 6）•

C-D 間を表す PIPE MANNING 型要素（要素 7）•

地点 D で直径減少を表す PIPE CONTRACTION 型要素（要素 8）•

D-E 間を表す PIPE MANNING 型要素（要素 9）26

•

地点 E でバルブを表す PIPE GATE VALVE 型要素（要素 10）•

位置 E(要素 10)のバルブ用のタイプ PIPE GATE VALVE 型要素。

•

E-F 間を表す PIPE MANNING 型要素（要素 11）•

貯槽への出口を表す PIPE ENLARGEMENT 型要素（要素 12）貯槽入口要素は無いことに注目。貯槽入口は直径増加で表される。液体がネットワークを離れることを表す•

ダミーのネットワーク出口要素(要素 13)

図 25: 入力ファイル（省略。pipe.inp 参照）
入力ファイルでは、これらの要素はすべて D タイプ要素で定義される。節点は正確な座標で記述する。*FLUID SECTION によって、要素をそれぞれ適切に記述する。ダミーのネットワー
ク入口および出口要素は *FLUID SECTION で定義される。
ネットワークの材質特性は*DENSITY に記述される密度、*FLUID SECTION に記録される比熱および動粘度である。比熱は伝熱がモデル化されている場合のみ必要である。ここでは必要ない。水の動粘性係数は 1750xE-6 Ns/m[29]である。境界条件は節点 3 および 25 の大気圧である。
両方とも貯槽の液体表面を表す。圧力はネットワーク要素の角要素において自由度 2 である。
ネ ッ ト ワ ー ク 分 析 は *COUPLED TEMPERATURE-DISPLACEMENT ま た は *HEATTRANSFER で行われる。ここで構造を考慮に入れない。そのため伝熱解析を行う。最後に、重力荷重が指定される。これは水力ネットワークに不可欠である。結節点の出力では、NT は自由度 0 が必要である。
一方、U は自由度 1 から 3 である。今回は質量フラックス（中央節点の DOF1）および圧力（角要素の DOF2）に興味があるため、*NODE PRINT の次行に U を記述する。公式に U は変位であり、そのため.dat ファイルには変位と記入される。
.dat ファイル中の結果は以下のようになる。
パイプ中の質量フラックス(中央節点の DOF1、1 番目)は一定値を 89.592 kg / s である。これは[10]の結果 89.4 l / s によく合う。全ての節点および要素の定義が図 25 に記述されていない。
出力は、地点 A が節点 5、地点 B が要素 7-9、地点 C が 11-13、地点 D が 15-17、地点 E が 1921、地点 F が節点 23 に相当する。出力ファイルの 2 番目は圧力である。ベンド、バルブおよび収縮は圧力を低下させる
一方拡大は圧力を増加させる（速度が低下する）。
ネットワークの構造（パイプ壁など）もモデル化されれば、流体圧力は自動的に要素表面に割り当てられる。これは*DLOAD 中の PxNP ラベルによって行われる。

5.7 蓋で駆動されるキャビティ （liquid1.inp）
蓋で駆動されるキャビティは、非圧縮粘性流れの有名なベンチマーク問題である[69]。問題の形状を図 26 に示す。ここでは正方形のキャビティを考える。三方を滑り無し固定壁で囲まれ、27

天井が接線方向に単位長さで動く。下端左角を基準静圧力ゼロとおく。本問題では レイノルズ
数 400 の流速および圧力分布を考える。
入力ファイル（liquid1.inp）を図 27 に示す。問題は 2 次元であるが、単位厚さをもつ 3 次元問題としてモデル化される。というのは CalculiX は 3 次元流れしか扱えないためである。図 28に 2D 射影されたメッシュを示す。6 節点くさび要素が使用される。厚さ方向は 1 層である。結果が層厚方向において同じため、これは十分である。入力ファイルは節点座標および要素で始まる。流体の要素タイプは構造要素と同じであり、C が F で置換された F3D6 を用いる。蓋をの節点および滑り無し壁に属する節点は、それぞれ節点群 Nin および Nwall に集められる。節点
群 N1 は印刷のために作成される。これは蓋に近い節点の集合である。
*BOUNDARY の次行に斉次境界条件（値ゼロ）が記述される。壁における速度はゼロ（滑り無し）である。法線方向および厚さ方向も同様にゼロである。さらに、キャビティ下端左角の参照点（節点 1 および節点 2376）は圧力ゼロである。材質定義は、密度、熱容量、動粘性係数から構成される。密度は 1 にセットされる。熱容量と動粘性係数は *FLUID CONSTANTS の次行に記述される。計算は定常であり熱容量は必要ないため、この値は無関係である。レイノルズ数が 400 になるよう動粘性係数の値を選んだ。レイノルズ数は代表速度と代表長さの積を動粘性率で割って定義される。蓋の速度は 1 である。その長さは 1 である。また、密度が 1 であるため、粘性係数および動粘性係数は同じである。従って、動粘性係数は値 1/ 400 をとる。 *SOLIDSECTION で材質特性が要素に割り当てられる。

28

図 26: 蓋で駆動されるキャビティの形状
No Data
図 27: 蓋で駆動されるキャビティの入力データ29

図 28: 蓋で駆動されるキャビティのメッシュ分割

問題の未知数は流速と静圧である。温度境界条件は提供されない。したがって温度は無関係である。
*INITIAL

CONDITIONS,TYPE=FLUIDVELOCITY,*INITIAL

CONDITIONS,TYPE=PRESSURE
によって未知数の初期値をゼロにする。流速はそれぞれの自由度について初期条件を指定する必要がある。
ステップは *STEP で始まる。しかし、繰り返しの最大数は流体計算向けであり、パラメタ
INCF で管理される。定常状態の計算は *STATIC が使用される。この次行の値は流体計算には関係ない。計算が定常的に行われるように増分が自動的に選択される。蓋のゼロでない接線速度は *BOUNDARY の次行に記述される。非斉次（非ゼロ）境界条件は 1 つのステップ内で定義されなければならない。ステップは節点の流速 V および静圧 PS の印刷リクエストで終わる。印刷の周期は FREQUENCYF パラメタによって 200 で定義される。これは 200 増分ごとに結果が保存することを意味する。
x 方向(蓋の接線方向)の速度分布は、図 29 に示される。最も小さな値( −0.33 )およびその位置
は文献[69]に非常によく一致する。図 30 は速度ベクトルを示す。蓋の近く、下端左角および下端右角に大きな変化がある。図 31 の圧力分布は渦中心および上端左角の低圧域を示す。上端右角は、渦の x 成分のよどみ点であり、圧力上昇によって特徴づけられる。
30

図 29: 蓋で駆動されるキャビティの流速の x 成分31

図 30: 蓋で駆動されるキャビティの流速分布32

図 31: 蓋で駆動されるキャビティの圧力分布

5.8 跳水のある水路（channel2.inp）
開きょ流れが一次元のネットワークでモデル化できることはあまり知られていない。支配方程式
である Bresse 式（6.8.18 節参照）および可能な流体タイプが 6.5 節に記述される。
この問題の入力を図 32 に示す。それは CalculiX test の 1 つである。水路は 1 行に並んだ 3
節点ネットワーク要素(タイプ D)6 つから構成される。節点は架空の座標系をもつ。それらは計算を入力しない。しかし、それらは .frd ファイルに保存される。CalculiX GraphiX による可視化は、正確な座標を使用すると有利な場合がある。ネットワークでは通常エントリーの最終節点および出口要素にはラベルゼロがある。材質は水で、密度、熱容量および動粘性係数が特徴である。次に *ELSET で要素が適切な要素群に保存される。要素群は *FLUID SECTION で参照される。
*FLUID SECTION を分析すると水路の構造が明らかになる。上流に水門があり、下流には大きな貯槽がある。また両者は直線の水路によって接続される。ほとんどの上流要素のように、水門は 2 つの要素から構成される。水門要素および水門開放要素である。水門はその下端で水深を固定するが、この水深は下流の水位によって制御される背水曲線によって超えられる。

Available downloads for the solver (CalculiX CrunchiX: ccx) : test example 参照33

図 32: 跳水のある水路の入力

図 33: 跳水のある水路の水深
水門は幅 10 m （水路にわたって一定）、傾斜 0.005 、ゲート高 0.8 m である。さらに、下流の水門開放要素は(3)のように表される。水門開始要素は同じ幅および傾斜を持っている、長さ 0.1mである。長さが負の場合、真の長さは節点座標から計算される。角 θ はゼロである。それは断面が台形ではなく長方形であることを意味する。*FLUID SECTION に MANNING オプションがあるため、その次のパラメタ( 0.01m−1 3s )は Manning 係数である。上流の水門要素のラベルは(2)で与えられる。水路要素の定数は節 6.5 に記述される。ここで重要なのは長さ 49.8 m である。
最後の要素(貯槽)は再び非常に短い要素(長さ 0.1m )である。水門開口部や貯槽のような要素（実際に物理的な長さを持っていない要素）の長さは、小さくしておかれるべきである。
次に、境界条件を定義する。貯槽の水深は 2.7m である、一方、質量流量は 60000 kg / s である 。 CalculiX の ネ ッ ト ワ ー ク 計 算 は 定 常 伝 熱 計 算 の 特 殊 な ケ ー ス で あ る た め 、 *HEATTRANSFER,STEADY STATE が使用される。駆動力は重力である。
CalculiX を実行すると、跳水が要素 4（直線水路要素）の相対的位置 0.67 にあるというメッセージが現われる。これは図 33 において明らかである。水路が描かれる。水門は x = 5 m に位置する。 x = 55 m から貯槽になる。チャネルの底は黒色である。水門後方の水位は規定されない。計算の結果 3.667 m であった。水門における水位は 0.8 m に制御される。前部水曲線（上流条件で制御される曲線）は下流に伝達し、 x = 38.5 m にある跳水によって後部水曲線（下流条件で制御される曲線）につながる。言い換えると、跳水は上流の臨界流れと下流の未臨界流れを接34続する。臨界水深は図の点線で示される。これはフルード数が 1 になる水深（臨界流れ）である。
水路流れでは、機械的変位の自由度が、質量流量、水深および臨界深度のために使用される。
したがって、*NODE PRINT 次行の U はこの情報を .dat ファイルに保存する。同じ情報は*NODE FILE 次行に MF、DEPT、HCRI を記述することにより、 .frd ファイルに保存される。

5.9 梁要素を用いた片持ち梁
以前に、厚い片持ち梁は体積要素でモデル化された。前節では似た例題（6.2.28）で 4 面体梁要素が使用された。梁要素は定義が簡単である。それらは線上の 3 節点から構成される。内部的に体積要素へ拡張される。2 つのタイプの梁要素がある。B32 要素（C3D20 要素に内部変換される）および B32R(低減積分)要素（C3D20R 要素へ内部変換される）。本節の結果に基づいて、
B32R 要素が高く推奨される。断面力が必要な場合は特に、B32 要素は回避されるべきである。
片持ち梁は長さ 100 mm 、 2 × 2 mm の正方形断面をもつ。梁の軸はワールド座標の z 方向にとる。梁は単に 1 つの要素でモデル化され端点に x 方向に単位荷重が載荷される（図 34）。我々は積分点節点 a および b における応力、固定点における断面力、自由端の変位に興味がある。積分
点 A の位置は x = −1

3 , y = −1(

3 , z = 50 1 + 1)

3 である。B の節点座標は x = −1 , y = 1 ,z = 100 である。材質は等方性線形弾性をもつ。ヤング率 100, 000 MPa 、ポアソン比 0.3 である。
この入力ファイルは simplebeam.inp に似ている（図 35）。
積 分 点 の 応 力 は *EL PRINT に よ っ て 得 ら れ る 。 節 点 に お け る 応 力 は *EL FILE にOUTPUT=3D オプションを付けて得られる（OUTPUT=3D オプションはデフォルトで付く）。
一方、断面力は SECTION FORCES オプションを付けて得られる（*EL FILE, SECTIONFORCES ）。このオプションは OUTPUT=3D と排他である。変位は拡張していない視界（つまり OUTPUT=2D オプション）が最もよい。これは本例では、2 回計算することを示す。
OUTPUT=3D オプションで 1 回、SECTION FORCES オプションで 1 回である。

図 34: 梁の形状35

図 35: 梁の入力
結果を表 3 に示す。{mm,N,s,K}単位系を使用した。reference は単純梁理論を用いた解析値[51]である。結果は理論と全体的に一致している。積分点における応力は正確に一致する。また、外挿された節点に作用する法線応力も一致する。せん断力は特別の対処を必要とする。梁では、せん断力は、断面を横切って放物線状に変化する。二次の体積要素は断面を横切る線形の応力変化だけを計算できる。したがって、放物線的な変化は近似される断面を横切る一定せん断応力によって近似される。低減積分点（ ± 1

3 ）が放物線的な応力変化が平均値に達する点であるため、
積分点の値は正確である!節点への外挿値は同じ一定値をとり、当然違う。角の正確な値は 0 であるため。
断面力は以下の手順で得られる。
1. 積分点(積分点 a のような要素内部)で応力を計算する。
2. それらの応力から(節点 b のような)角節点へ外挿する
3. 隣接する角節点間の補間によって中央節点における応力を計算する
4. 断面にある全節点の応力を用いて面内の低減積分点(積分点 c のような面の形状関数を使用する点)における応力を内挿する
36

5. これら応力を数値的に積分する。
表 3 に示されるように、この手続きは正方形梁の断面力を正確に評価する。
梁先端の変位(表 3 の ux)は計算と理論が 10%違う。端点におけるせん断力による梁の変形は 3
次オーダーである。しかし、C3D20R 要素は 2 次までしかシミュレートできない。5 要素を用いると、偏差は 2.4% まで減少する(表 4)。積分点が固定点に近づいたことに注意する。1 要素の場
合と同じ位置であるが、要素が固定点の隣りにきた？
同じ梁に自由端で 1Nmm のトルクが作用場合を考える。結果を表 5 に示す。
トルクは完全に一致する。端点におけるねじれは理論値と 15%違う。
（uy が節点 b に相当する点における y 方向変位である）節点 b におけるせん断応力は明らかに違う（角節点ではせん断応
力は発生しない）。しかし、積分点面内の節点の内挿値を用いた積分値は正確なトルクになる。
より多くの要素を使用しても表 5 の値は変わらない。
同じ例題を円形断面で計算してみる。(半径= 1mm 、同じ長さ 100 mm 、矩形断面と同じ境界条
件および材質)円形断面では、要素の頂点節点を x, y = ±0.7071, ±0.7071 に置く。一方、中央節
点を x, y = 0, ±1 、 x, y = ±1, 0 に置く。積分点を x, y = ±0.5210 に置く。1 要素の結果を表 6 に、
5 要素の結果を表 7 に示す。
表 3: 正方形断面梁の曲げ計算結果（１要素）
resultvalue

Reference

34.15134.151σ xz ( a )

−0.25

−0.25

Fxx
M yy

−1.

−1.

100.

100.

σ zz ( b )

75.

75.

σ xz ( b )

−0.25

0.

ux

2.25

2.50

σ zz ( a )表 4: 正方形断面梁の曲げ計算結果（5 要素）
resultvalue

Reference

41.47141.471σ xz ( a )

−0.25

−0.25

Fxx
M yy

−1.

−1.

100.

100.

σ zz ( b )

75.

75.

σ xz ( b )

−0.25

0.

ux

2.44

2.50

σ zz ( a )表 5: 正方形断面梁のねじれ計算結果（１要素）
resultvalue
37

Referenceσ xz ( a )−0.21651−

σ yz ( a )−0.21651

−

M zz

1.

1.

−0.3750

σ yz ( b )−0.3750uy9.75 × 10−4

1.1525 × 10 −3

σ xz ( b )表 6: 円形断面梁の曲げ計算結果（1 要素）
resultvalue

Reference

34.00

52.26

σ xz ( a )−0.322−0.318Fxx
M yy−0.99996

−1.

58.7

100.

σ zz ( b )

62.8

90.03

σ xz ( b )−0.322−0.318

ux

2.91

4.24

σ zz ( a )表 7: 円形断面梁の曲げ計算結果（5 要素）
resultvalue

Reference

59.77

63.41

σ xz ( a )−0.322−0.318Fxx
M yy−0.99996

−1.

102.

100.

σ zz ( b )

109.

90.03

σ xz ( b )−0.322−0.318

ux

3.86

4.24

σ zz ( a )1 要素の場合、せん断応力が理論値に近い。せん断力はもっと近い。これは円形周囲を区分的に二次近似した断面積が小さく、わずかに高いせん断応力を補償するという事から説明される。同様の結果はトルクの場合も発生する。しかし、積分点における法線応力は理論値とかなり違う。
曲げモーメントがかなり小さな節点の場合と同様である。同じことは x 方向の変形にも当てはまる。5 要素を使用すると著しく改善する。曲げモーメントは理論値とわずか 2%の差である。自由端変位は 9%の差である。ここで再び変位が 3 次オーダーであることを議論する。二次要素は二次オーダーまでしか計算できない。より多くの要素の使用が結果を改善する。
形断面梁にトルクを作用させた結果を表 8 に示す。
（1 要素の場合。5 要素でも結果は同じ）
トルク（Mzz）は理論値と一致するが、積分点におけるせん断応力は理論値と 6%の差がある。
このため頂点におけるせん断値には理論値と 19%の差が発生する。
38

表 8: 円形断面梁の曲げ計算結果（1 要素）
resultvalue

Referenceσ yz ( a )−0.309−0.331

M zz

0.999994

1.

−0.535−0.450σ yz ( b )−0.535−0.450uy

1.54 × 10−3

1.66 × 10−3

σ xz ( a )−0.309σ xz ( b )−0.331

積分点への内挿はせん断力を −0.305 MPa と評価する。これら応力の積分が最終的に完全なトル
クとなる。梁端点のねじれ角は異論値と 7%違う。
以上をまとめると、C3D20R 要素を用いると良い結果が得られる。

矩形断面では:•

断面力は正確である。

•

積分点における応力は正確である。

•

十分な要素が使用されれば、曲げの変位は正確である。

•

ねじれ角は理論と多少離れる(15%)。

円形断面の場合:•

せん断力およびトルク断面力は正確である。

•

十分な要素が使用される場合、曲げモーメントは正確である。

•

十分な要素が使用されれば、曲げの変位は正確である。

•

ねじれ角は理論と多少離れる(7%)。

一般に、断面力から応力を計算することが推奨される。ただ一つの欠点は C3D20R 要素がアワーグラス問題、つまり不思議な変位に結びつく場合があることである。しかし、断面図を横切る
平均変位は十分理論と一致する。加えて、このアワーグラス問題の場合、非線形計算が収束しない問題がある。これは、CalculiX では拡張節点座標をわずかに乱す(0.1%) ことにより改善された。
同様計算が B32 要素を用いて行なわれたが、結果は思わしくなかった。特に曲げにおいて断面力は理論値から離れた。

図 36: 具体的な片持ち梁の入力（省略）5.10 鉄筋コンクリート片持ち梁（concretebeam.inp）
本節の目的は、鉄筋コンクリートの自重による応力の計算である。本問題の特徴は、複合材構造、
39

圧縮のみの材質（コンクリート）の扱いである。

梁は断面 1×1 m 、長さ 10 m である。コンクリートの密度は 2350 kg / m である。鋼の密度は
237800 kg / m3 である。それぞれヤング率は 14000 MPa および 210000 MPa である。鋼は梁上面（引張側）のみ、表面から 9.5cm に位置する。層厚は 1cm である。（実際には鋼はコンクリート
内部に棒状に存在する。薄い層は近似モデル化である。棒の断面積が層の断面積に等しいことを確認すること。）シェル構造を使用すると複合材質を簡単に表現できる。この機能は梁要素ではまだ使用できない。従って、梁は幅 1m および長さ 10 m の板としてモデル化した。図 36 の*SHELL SECTION の次行に、層厚とそれらの材質を記述する。梁の上から順に記述する。方向（上から下）はシェル要素の法線方向により制御される。*ELEMENT の下に記述される要素の節点の順番で制御される複合材シェルの場合、それぞれの層内に積分点が 2 つある。S8R 要素を使用する。中立軸の位置を求めるため、コンクリート部は複数の層でモデル化された。コンク
リート部に 10 層、鋼部に 1 層使用した。

コンクリートは張力を保持できないが、圧縮下では線形弾性的に変形する。この性質は
COMPRESSION ONLY でモデル化する。CalculiX では、これはユーザー材質の例である。ユーザー材質の名前は固定文字で表す。今回は COMPRESSION ONLY とした。材質は 80 文字まで記述できるので、残り 64 文字は自由に選択可能である。今回は他に文字を入れない。複数の圧縮材が使用され、それらを区別する場合には追加の文字を入れる。COMPRESSION ONLY
材質は 2 つの定数で特徴づけられる。第一はヤング率、第二はユーザーが許可する引張応力である。
今回は 0.1MPa とした。SI 単位を使用した。
単純梁理論では、鋼の引張応力 152.3 MPa 、コンクリート下端で最大圧縮応力 7.77 MPa となる。有限要素計算はそれぞれ 152 MPa 、 7.38 MPa を予言する(図 37)。これはほぼ理論値である。
CalculiX では、複合構造の画像は常に三次元に拡張される。図 38 から複合材の正確な寸法および鋼部の高い引張応力が分かる。

40

図 37:固定端における梁の高さ方向の軸方向応力41図 38:固定端における梁の高さ方向の軸方向応力

図 39: 変形したシート内の最大主応力5.11 薄いシートのしわ（leifer1.inp、leifer2.inp）
本問題の入力ファイルは leifer1.inp と leifer2.inp である。これは J.Leifer によって 2003 年に最初に考案された。構造は、1 辺 229 mm の正方形、厚さ 0.0762 mm の薄いシートである。1
辺が固定され、これに並行に反対側に向かって 1 ㎜動かす。ヤング率とポアソン比はそれぞれ3790 MPa および 0.38 である。実験結果はせん断変形によるしわの生成を確認する。
この実験を 2 つのアプローチでシミュレートする。両方において、シートは二次シェル要素でシミュレートされる。1 番目のシミュレーション(leifer1)では、材質は等方的な線形弾性材質であり、しわはシートの座屈により発生する。この座屈を拡張するため、シートに直行する座標系（本シミュレーションでは z 方向）をわずかに偶然性の方法で乱す。入力ファイルの座標系を確認すること。さらに、シミュレーションは非常に小さな時間ステップで動的に行なわれる。図
39 に、変形したシートの中で最大主応力を示す。( x = 0 端が固定された。 x = 229 の端は y の負方向に 1mm 移動した)しわが発生する。最も小さな主応力箇所では、シートに仮想的に圧力が作用しない。それらは座屈により削除された。この種のシミュレーションは非常に長い時間がかかる。(1 つの時間ステップに 336 回計算が必要)。
42

図 40: 等方シミュレーション中のせん断力43

図 41: 引張シミュレーション中のせん断力
圧力が無いため、正確な圧力を求める 2 番目の方法を試みる。材質を等方性としてシミュレー
トする代わりに引張のみの材質モデルを使用する。これは収束がはるかに早い長所(計算時間が短い)がある。図 40 および 41 は、両方のシミュレーションのせん断力を比較する。それらはよく一致する(しわが発生しない等方性シミュレーション中のせん断力分布は全く異なる)同じこと
が他の応力成分にも当てはまる。しかし、テンションのみの材質の使用は平面変形に結びつかない。最も小さな主歪みを見ることによりしわを間接的に導出する。(図 42)大きな負の値はしわの存在を指す。

44

図 42: 引張力のみのシミュレーションにおける最小主歪45