Scialbを用いた不均一な厚みの板の解析

「テーマ」

 不均一な厚みを持つ板の固有振動数を求める

「方針」

 Scilabを用いた薄板四角形の解析では均一の板厚しか扱えませんでした。ここではFEMの最小要素別に厚みを変えることで不均一な厚みの板の固有振動数を計算してみます。

 結果はFreeCAD計算結果と比較することで検証してみます。


1.解析モデル

1.1 モデル方針

 検討対象は他解析と共通とするため以下で統一しました。

・500×500のアルミ板、厚みは図示のようなおたまじゃくし形状。

・1辺のみ固定

「Detail1」変位法による薄板FEM解析詳細で行っていた均一薄板の解析方法では扱えない形状です。

1.2 Scilaモデル

 以下のモデルを考えて行きます。

・材料:AL(E=7e10[Pa]、ρ=2800[kg/m3]、ポアソン比 0.3)

・寸法:500×500×厚み(厚みは下図)

・分割数:11×11節点10×10要素

・1辺固定:①~⑪節点固定

1.3 不均一な厚みの設定

 ここで変位法による薄板FEM解析において厚みを入れ込む部分を検討します。

もともと薄板におけるFEM解析方法は以下の手順(概略)で計算しました.

詳細は「Detail1」変位法による薄板FEM解析詳細及び参考文献[1]4.2章をご確認ください。

(1)薄板要素の基礎方程式

 D(∂4ω/∂x4+2∂4ω/∂x2∂y2+∂4ω/∂y4)=0

 D=Et3/(12(1-ν2))

   E : 板要素の縦弾性係数

   t:板厚

   ν:板要素のポアソン比

   ω:板のz方向の変位

(2)剛性行列[K]

 固有値計算に用いる剛性行列[K]は以下となります。

   [K] = [C]-T∫[Q]T[D] [Q]dxdy[C]-1

 ここでそれぞれの行列は以下の定義となります。

〇[C]について

ω(x,y)={ 1, x, y, x2, xy, y2, x3, x2y, xy2, y3, x3y, xy3 }{c}

という仮定式から境界条件{δ}と四隅座標データ(x,y)から以下の関係式として計算されます。

{δ}=[C]{c}

  {δ}={ω1, ωx1, ωy1, ω2, ωx2, ωy2, ω3, ωx3, ωy3, ω4, ωx4, ωy4}

  {c}={c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12}T

〇[Q]について

ω(x,y)={ 1, x, y, x2, xy, y2, x3, x2y, xy2, y3, x3y, xy3 }{c}

という仮定式からひずみ{ε}と四隅座標データ(x,y)から以下の関係式として計算されます。

{ε} = [Q]{c}

  {ε} = { -∂2ω/∂x2, -∂2ω/∂y2, 2*∂2ω/∂x∂y}

  {c}={c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12}T

〇[D]について

材料力学より曲げモーメント{Mx,My,Mxy}Tとひずみ{ε}から以下の関係式として計算されます。

{Mx, My, Mxy}T=[D]*{ε}

[D]=Et3/(12*(1-ν2))*

1ν0
ν10
00(1-ν)/2

 以上剛性行列[K]の中で薄板の厚みが関係するのは[D]の行列中のt3のみとなります。

そこで最小要素別にこのtを個別に設定することで不均一な薄板をモデル化してみました。

 かなり粗い近似ですが、まずは単純な計算としました。

(3)整合質量行列[M]

 固有値計算に用いる整合質量行列[M]は以下となります。

[M]=∫{N}T{N}ρdv

以下の2式より変位関数{N(x,y)}={ 1, x, y, x2, xy, y2, x3, x2y, xy2, y3, x3y, xy3 }[C]-1

  ω(x,y)={N(x,y)}{δ}

  ω(x,y)={ 1, x, y, x2, xy, y2, x3, x2y, xy2, y3, x3y, xy3 }{c}

    ={ 1, x, y, x2, xy, y2, x3, x2y, xy2, y3, x3y, xy3 }[C]-1{δ}

 以上整合質量行列[M]の中で薄板の厚みが関係するのは[M]の積分結果から出てくるtのみとなります。そこでこちらも最小要素別にこのtを個別に設定することで不均一な薄板をモデル化してみました。


2.計算結果

 以下にScilabを用いた固有値計算の結果とFreeCADの計算結果を示します。

 結果は固有値も形状も十分な一致を確認出来ました(見やすくするため形状は反転を示しているものもあります)。

〇Scilab 1次mode:4.53Hz
〇FreeCAD 1次mode:4.65Hz
〇Scilab 2次mode:7.82Hz
〇FreeCAD 1次mode:8.88Hz
〇Scilab 3次mode:13.90Hz
〇FreeCAD 1次mode:13.93Hz
〇Scilab 4次mode:21.78Hz
〇FreeCAD 1次mode:21.92Hz
〇Scilab 5次mode:31.80Hz
〇FreeCAD 5次mode:32.01Hz

 尚FreeCADのメッシュ条件は以下としました。

・Netgen(FreeCAD装備)

 [最大要素20mm / 中程度の粗さ]

 板厚を敢えて10倍とし、密度を100倍に修正。

 (FreeCADは薄い板の計算が苦手です。詳しくは薄板三角形の固有振動数及び強度解析step1をご参照ください。)


3.まとめ

 非常に粗い近似ではありますが、FreeCADの計算結果と十分な一致をみることが出来ました。

 FreeCADと違い、Scilab計算結果は動的変化を解析することが出来ます。

 このため本結果をブレードのフラッタ解析に適用することで複雑なブレード形状を扱うことが可能となります。


参考文献

[1]機械構造 振動学 MATALBによる有限要素法と応答解析 小松敬冶 著 森北出版株式会社