シミュレーションのトビラ

シミュレーション技術に関する個人的メモ

シミュレーションのトビラのトビラ

こんにちは。猫の奴隷Bです。
このブログは、コンピュータを使ったシミュレーションに関する私的備忘録です。
順次追加されていく(はずです)ので、良かったらみていってください。

2023/10/27 私の会社のサイトが立ち上がったので、そちらにも転記しています。こちらもどうぞご覧ください。
https://www.neconote-engineer.co.jp/2023/10/27/neco%e3%83%bbno%e3%83%bbmemo/

EXCELで有限要素法の強度解析(6)

全体剛性マトリクス

 前回節点ごとの要素剛性マトリクスを求めましたので、これをモデル全体について作成します。
 これまでのところをおさらいすると、
 要素剛性マトリクス [  K_e] は、

 [ K_e ] = t \Delta [ D ]^T [ B ]^T  [ B ]

でした。
 要素剛性マトリクスは、要素ごとに存在しますので、イメージ的には下図のようになります。
[Ke]のイメージ
 要素剛性マトリクス [  K_e] にはBマトリクスが含まれていますが、Bマトリクスは節点ごとの成分を含んでいます(後述します)。そのため、要素剛性マトリクス [  K_e] も節点ごとの成分を持ちます。三角形一次要素ですので、一つの要素に3つの接点があり、それぞれx座標とy座標がありますので、ひとつの要素ごとに6つの成分を含むことになります。

Dマトリクスのイメージ

 要素剛性マトリクス [  K_e ] の式の中の、Dマトリクス [ D  ] は、

 [ D ] = \dfrac{E}{(1+\nu)(1-2\nu)} \left[ \begin{array}{cccccc} 1-\nu & \nu & \nu & 0 & 0 & 0 \\ \nu & 1-\nu & \nu & 0 & 0 & 0 \\ \nu & \nu & 1-\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{2} ( 1-2 \nu) & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{2} ( 1-2 \nu)  & 0\\ 0 & 0 & 0 & 0 & 0 & \frac{1}{2} ( 1-2 \nu)  \\ \end{array} \right]

であり、材料の物性によって決まるマトリクスで、形状には依存しませんでしたね。
 今回は平面ひずみモデルを考えますので、z座標を省略すると、Dマトリクス [ D  ] は下記のように簡略化できます。
 [ D ] = \dfrac{E}{(1+\nu)(1-2\nu)} \left[ \begin{array}{ccc} 1-\nu & \nu & 0 \\ \nu & 1-\nu  & 0 \\  0 & 0 & \frac{1}{2} ( 1-2 \nu)  \end{array} \right]

 応力とひずみの関係式は、Dマトリクスを用いて\{\sigma\}=[D]\{\varepsilon\}という形でした。つまり、Dマトリクスは応力とひずみとを関連付けるマトリクスであり、下図のようなイメージとなります。
[D]のイメージ
 よって、平面ひずみモデルの場合には、Dマトリクスは3行3列のマトリクスとなります。

Bマトリクスのイメージ

 要素剛性マトリクス [ K_e ] の式の中の、Bマトリクス [ B  ] は、

 \begin{align} [ B  ] &= \left[ \begin{array}{cccccc} \dfrac{\partial N_1}{\partial x} & 0 & \dfrac{\partial N_2}{\partial x} & 0 & \dfrac{\partial N_3}{\partial x} & 0 \\ 0 & \dfrac{\partial N_1}{\partial y} & 0 & \dfrac{\partial N_2}{\partial y} & 0 & \dfrac{\partial N_3}{\partial y} \\ \dfrac{\partial N_1}{\partial y} & \dfrac{\partial N_1}{\partial x} & \dfrac{\partial N_2}{\partial y} & \dfrac{\partial N_2}{\partial x} & \dfrac{\partial N_3}{\partial y} & \dfrac{\partial N_3}{\partial x} \end{array} \right] \end{align}

で、これは形状によって変わる形状関数N_nを含んでいます。
 ちなみに形状関数N_nは、下記のような式でした。
\begin{cases} N_1(x,y)=\dfrac{1}{2\Delta} \{ (y_2-y_3)x + (x_3-x_2)y +(x_2y_3-x_3y_2) \} \\ N_2(x,y)=\dfrac{1}{2\Delta} \{ (y_3-y_1)x + (x_1-x_3)y +(x_3y_1-x_1y_3) \} \\ N_3(x,y)=\dfrac{1}{2\Delta} \{ (y_1-y_2)x + (x_2-x_1)y +(x_1y_2-x_2y_1) \} \end{cases}

 つまり、ひとつの三角形一次要素を構成している3つの接点の情報を持つマトリクスとなります。
 ひずみと各節点の変位の関係式は、Bマトリクスを用いて\{\varepsilon\}=[B]\{u\}という形でした。つまり、Bマトリクスはひずみと節点の変位(座標)とを関連付けるマトリクスであり、下図のようになります。
[B]のイメージ
 よって、平面ひずみモデルの場合には、Bマトリクスは3行6列のマトリクスとなります。

要素剛性マトリクスの中身のイメージ

 DマトリクスとBマトリクスのイメージから、要素剛性マトリクス [  K_e ] の具体的なイメージを考えてみます。
 まず、Dマトリクスの転置行列とBマトリクスの転置行列を掛けていますので、イメージ的には下図のようになります。

[D]t[B]tのイメージ
 Bマトリクスが「節点別の方向ごとのひずみ成分」だったのに対して、Dマトリクスを乗じたことで、「節点別の方向ごとの応力成分」になりました。
 次に、これに転置ではないBマトリクスを乗じるわけですから、イメージ的には下図のようになります。
[Ke]の中身のイメージ
 よって、平面ひずみモデルの場合には、Keマトリクスは6行6列のマトリクスとなります。
 この要素剛性マトリクス [K_e]が要素の数だけ作成されることになります。

全体剛性マトリクス

 要素剛性マトリクス [K_e]のイメージがつかめたところで、いよいよモデル全体に展開します。
 要素剛性マトリクス [K_e]の各節点番号は、前述の説明では便宜上①~③として記載しましたが、実際にはモデル全体における節点番号で記載されます。今回の片持ち梁の例では、要素は1~6の三角形一次要素と、①~⑧の接点により構成されています。
 例えば、「要素1」は節点①②⑦の3つによって構成されていますので、要素剛性マトリクス [K_1]は①x①y②x②y⑦x⑦yの6行6列のマトリクスとなります。
 同様に、「要素2」は節点②③⑥の3つによって構成されていますので、要素剛性マトリクス [K_2]は②x②y③x③y⑥x⑥yの6行6列のマトリクスとなります。その他の要素も同様です。

 また、モデル全体の全体剛性マトリクス [K]は、8個の節点のxy座標のマトリクスとなり、合計で16列16行のマトリクスとなるわけです。
 要素ごとの要素剛性マトリクス [K_e]を全体剛性マトリクス[K]に当てはめてみます。分かりやすいように、要素1と要素2だけ表示すると、下図のようになります。

[K]のイメージ
 ここで、節点②に着目すると、要素1と要素2で重複しています。モデルをみても、節点②を共有していますね。この場合、節点②については、「要素1」の要素剛性マトリクス [K_1]と、「要素2」の要素剛性マトリクス [K_2]を足し合わせる形となります。
 この計算を6個の要素全てについて行い、モデル全体の全体剛性マトリクス [K]を完成させます。

 次は、境界条件を設定して、マトリクス(連立方程式)を解いていきます。

前へ 5.要素剛性マトリクス