鏡板の無い円筒圧力容器の解析をしてみます。単位系はN-mm系を想定しています。
円筒の内半径は200mm、外半径は201mm。つまり肉厚は1mmです。高さは250mmとしました。円を正36角形で近似しました。
sheet 4 200 0 0 201 0 0 201 0 250 200 0 250
revolve 0 0 0 0 0 10 36
基本節点間隔は、25mmです。節点数1,728、表面三角形要素数3,456です。寸法の範囲は、(-201, -201,0)-(201,201,250)。基本節点間隔次第では一様なメッシュが出来ます。
図-1 表面パッチ(Meshman_ParticleViewer_HPCで表示)
そのままでメッシュ分割をすると要素が細かく成り過ぎるので、メッシュ作成時に「表面形状を補正する」にチェックをせずにメッシュを切ります。補正を止めるとアルゴリズムの頑健性が落ちますが元の形状や表面のメッシュの模様も維持されます。メッシュ生成に失敗する時には、チェックを入れると良いでしょう。
要素数5,472、節点数10,656の四面体二次要素が作成されました。
ヤング率は200,000MPa、ポアソン比は0.3としました。
ここ迄来て、対称性を利用して解けば良かったと気付きましたが、構わず進めます。境界条件の設定が少し難しくなります。
自由に半径方向に伸び縮み出来るようにする為、拘束を緩めに設定します。
図-2 節点8(0,-201,0)を選択中
拘束節点を選択する時に節点座標値が分からない為意図した節点がどうか分からない問題が有りました。AdvOnWinに機能追加をしました。図-2で見て分かるようにステータスバーを新設して、ピックした節点と面グループの情報を表示するようにしました。これで意図した節点をピックしやすくなります。Ver. 0.44bとしてリリースしようと思います。
図-3 節点8の拘束条件(XYZ方向)
これで3つの並進剛体モードを拘束しました。
図-4 節点26(0,201,0)を選択中
図-5 節点26の拘束条件(XZ方向)
これで、X軸回りとZ軸回りを拘束しました。
図-6 節点35(201,0,0)を選択中
図-7 節点35の拘束条件(Z方向)
図-8 面グループ0を選択中
剛体モードを拘束する拘束条件が分かり難いので、別に説明図を作成しました。
図-9 全ての拘束条件
節点8のZ方向と節点26のZ方向でX軸回りの回転を拘束します。節点8のX方向と節点26のX方向でZ軸回りの回転を拘束します。
節点8のZ方向と節点35のZ方向でY軸回りの回転を拘束します。
勿論剛体モードを止めるたけの拘束方法は1通りだけではありません。
荷重は内圧としました。
図-9 内圧条件(0.2MPa)
全ての境界条件をテキストで表示します。境界条件設定用の3D画面(BCAgent Ver.0.43bと言うタイトルの画面です)において情報(I)>境界条件(B)を選択します。
図-10 全ての境界条件をテキスト表示するメニュー項目
図-11が表示されます。これは実はcndという拡張子で保存される境界条件ファイルをテキスト表示しただけです。
図-11 全境界条件のテキスト表示
内圧p[MPa]を受ける薄肉円筒の周方向応力の理論解。
$$\sigma_\theta =\frac{a}{t}p$$
但し、a:内半径、t:板厚。
実際の値を入力すると、
$$\sigma_\theta =\frac{200}{1}0.2=40.$$
出典:日本機械学会計算力学技術者認定事業標準問題集1 計算力学技術者2級第9版(2016/07/31)付録6.2節。
図-12 最大主応力(周方向応力)37.5MPa～43.5MPa
誤差は-6.3%～8.8%です。少し大きいですね。図-12にZの端に近い方に赤や青の集中部が見えてますが、本来これは無い方が良いですね。図-1の表面パッチと見比べるとメッシュが一様分割では無い事が原因のようです。
円に内接する正36角形で円を近似しているので、周の実際の長さLは
$$L =2nasin(\frac{\pi}{n})$$
ここでn:角数、a:内半径。
真の周は1257mm、36角形の周は1255mm。よって圧力が0.1%本当より高く加わるので、正確には割り引く必要が有ります。差が少ないのでここでは割引計算は省きます。
図-13 X方向変位の分布(-4.08e-2～4.21e-2)
内圧p[MPa]を受ける薄肉円筒の内半径の増加量の理論解。
$$\delta a =\frac{a^2p}{tE}$$
但し、a:内半径[mm]、t:板厚[mm]、E:ヤング率[MPa]、p:内圧[MPa]。
実際の値を代入すると、4.e-2[mm]です。誤差は-2%から+5%です。
図-14 Y方向変位の分布(-1.58e-3～8.15e-2)
Y方向の場合は、理論解は同じ筈ですが、計算値はX方向変位と随分値が違いますね。理由は自分で考えてみて下さい。
次はZ方向歪です。
図-15 歪zz(-6.13e-5～-5.81e-5)
内圧p[MPa]を受ける薄肉円筒の軸方向歪の理論解。
$$\epsilon_{zz} =-\nu\frac{ap}{tE}$$
実際の値を代入すると、6.e-5[mm]です。誤差は-2.1%から+3.2%です。
