Rで有限要素法:連立方程式

本ポストはこちらの続きです。

const typesetMath = (el) => { if (window.MathJax) { // MathJax Typeset window.MathJax.typeset(); } else if (window.katex...

本ポストではこちらのポスト(以降、該当ポスト)のコード中の

Rで有限要素法:片持ち梁-自由端集中荷重:たわみ量

const typesetMath = (el) => { if (window.MathJax) { // MathJax Typeset window.MathJax.typeset(); } else if (window.katex...

# ------------------------------------------------------------
# 5. 連立方程式を解く
# ------------------------------------------------------------
u <- numeric(n_dof)
u[free_dofs] <- solve(K[free_dofs, free_dofs], F[free_dofs])

w_fem <- u[seq(1, n_dof, by = 2)] # 各節点のたわみ
th_fem <- u[seq(2, n_dof, by = 2)] # 各節点の傾き dw/ds (FEM座標s方向)

s_nodes <- seq(0, l, length.out = n_node)
x_nodes <- l - s_nodes # 厳密解の座標系(自由端起点)に変換

のコード各行の意味と、変数について確認します。

コード各行の解説

u <- numeric(n_dof)
# 全自由度分(この例では12個)の解ベクトルuの入れ物を、すべてゼロで用意する。
# 固定端の自由度(u[1], u[2])は境界条件によりゼロのままでよいため、
# この時点で確保しておいた0が、そのまま最終的な答えになる

u[free_dofs] <- solve(K[free_dofs, free_dofs], F[free_dofs])
# 拘束されていない自由度(free_dofs)だけを取り出した部分連立方程式を解く。
# K[free_dofs, free_dofs]は、全体剛性行列Kから固定端に関する行・列を除いた
# 10×10の部分行列、F[free_dofs]は対応する10個分の荷重成分。
# solve(行列, ベクトル)は「行列 × 未知数 = ベクトル」の連立方程式を解くR関数で、
# 求まった10個の値を、uのうちfree_dofsに対応する位置にまとめて代入する

w_fem <- u[seq(1, n_dof, by = 2)]
# uの中から、奇数番目(1,3,5,...)の成分だけを取り出す。
# 自由度の並びが[w1,θ1,w2,θ2,...]なので、奇数番目はすべて「たわみw」の成分にあたり、
# これにより各節点のたわみだけを並べたベクトルができる

th_fem <- u[seq(2, n_dof, by = 2)]
# 同様に、uの偶数番目(2,4,6,...)の成分を取り出す。
# 偶数番目はすべて「傾きθ」の成分にあたるため、各節点の傾きだけを並べたベクトルができる

s_nodes <- seq(0, l, length.out = n_node)
# FEM側の座標s(固定端起点、s=0〜l)で、各節点の位置を等間隔に並べたベクトルを作る

x_nodes <- l - s_nodes
# 厳密解の座標系(自由端起点のx)に変換する。x = l - s の関係を、各節点についてまとめて計算する

出力結果の見方

u(全自由度の解、長さ12)

 [1] 0.0000000000 0.0000000000 0.0001706653 0.0008228506 0.0006338997 0.0014628454 0.0013165609 0.0019199846 0.0021455066 0.0021942682 0.0030475947 0.0022856960

並びは [w1, θ1, w2, θ2, w3, θ3, w4, θ4, w5, θ5, w6, θ6] に対応しています。

最初の2つ(u[1], u[2])が固定端の \(w_1=0, \theta_1=0\) になっている点で、境界条件が正しく反映されていることが確認できます。

残りの10個が、実際に連立方程式を解いて求まった値です。

w_fem(各節点のたわみ、長さ6)

w_fem

[1] 0.0000000000 0.0001706653 0.0006338997 0.0013165609 0.0021455066 0.0030475947

u から奇数番目だけを抜き出したもので、節点1(固定端)から節点6(自由端)まで、たわみが単調に増加していく様子が読み取れます。

最後の値 0.0030475947(m)= 約3.048mmが自由端のたわみで、該当ポストの厳密解 3.047595mm との比較に使われた値です。

th_fem(各節点の傾き、長さ6)

th_fem

[1] 0.0000000000 0.0008228506 0.0014628454 0.0019199846 0.0021942682 0.0022856960

同様に u から偶数番目だけを抜き出したもの。固定端で傾き0から始まり、自由端に向かって傾きが増えていく(たわみ曲線が自由端に近づくほど急になる)様子を表しています。

s_nodes(FEM座標、長さ6)

s_nodes

[1] 0.0 0.4 0.8 1.2 1.6 2.0

n_elem <- 5、支間長さ l <- 2.0 なので、要素長さ Le = 2.0/5 = 0.4 ごとに、固定端(0.0)から自由端(2.0)まで等間隔に並んだ節点位置です。

x_nodes(厳密解の座標、長さ6)

x_nodes

[1] 2.0 1.6 1.2 0.8 0.4 0.0

s_nodes を x = l - s で変換したもので、s_nodes とは逆順になっています(FEMでは固定端が原点、厳密解では自由端が原点のため)。

自由度の種類と対応する力の種類

自由度の種類	対応する力の種類	仕事	仕事の単位
\(w\)(たわみ、単位:m)	力(N)	力 × たわみ	N・m = J
\(\theta\)(傾き、単位:rad、無次元)	モーメント(N・m)	モーメント × 回転角	N・m = J

以上です。