Rで有限要素法:荷重と境界条件

コード各行の解説

F[n_dof - 1] <- F[n_dof - 1] + P
# 全体荷重ベクトルFのうち、「最終節点のw(たわみ)自由度」に集中荷重Pを加える。
# 自由度の並びは[w1,θ1,w2,θ2,...]なので、最終節点(n_node番目)のwは
# 全体でn_dof-1番目(θがn_dof番目)にあたる。「上書き」ではなく「加算」にしているのは、
# 他の荷重(等分布荷重など)がすでにFに入っている場合でも、それを消さずに追加できるようにするため

fixed_dofs <- c(1, 2)
# 固定端で拘束する自由度の番号を指定する。1番目は節点1のw、2番目は節点1のθで、
# 「固定端では、たわみも傾きもゼロ」という境界条件(w1=0, θ1=0)に対応する

free_dofs <- setdiff(seq_len(n_dof), fixed_dofs)
# 全自由度の番号(1からn_dofまで)から、fixed_dofsに含まれる番号を取り除き、
# 「拘束されていない(値が未知の)自由度」の番号リストを作る。
# setdiff(A, B)は「Aの中からBに含まれる要素を除いたもの」を返す関数

F の見方

F

 [1]    0    0    0    0    0    0    0    0    0    0 1000    0

要素数5・節点数6の設定なので、F は長さ12(節点数6×自由度2)のベクトルです。

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

• 11番目の成分だけが1000:
  • これは6番目の節点6(=最終節点、自由端)の自由度 \(w_6\) に、集中荷重 \(P=1000\)N が載荷されていることを表しています。
• それ以外がすべてゼロ:
  • 自由端の傾き(\(\theta_6\)、12番目の成分)や、その他すべての節点には、外から直接加わる力・モーメントがないことを表しています。

F は「各自由度に、外部から直接加わっている力(またはモーメント)の大きさ」を並べたベクトルであり、今回はその中で「6番目の節点のたわみ方向に1000Nの力が加わっている」として1箇所だけが非ゼロになっている、と確認できます。

free_dofs の見方

free_dofs

 [1]  3  4  5  6  7  8  9 10 11 12

これは、「まだ値が分かっていない(これから連立方程式を解いて求める)自由度」の番号一覧です。

• 1番と2番が含まれていない:
  • 節点1(固定端)の \(w_1, \theta_1\) は、境界条件により最初から0と分かっているため、「解くべき未知数」から除外されています。
• 3番から12番まで:
  • 節点2〜6(\(w_2,\theta_2,\dots,w_6,\theta_6\))の、合計10個の自由度が「未知数」として残っており、この後 solve(K[free_dofs, free_dofs], F[free_dofs]) によって、これら10個の値がまとめて計算されることになります。

つまり free_dofs は、「全体の連立方程式(12元)から、境界条件によってすでに値の分かっている2つの自由度を除いた、実際に解くべき10元の連立方程式がどの自由度に対応するか」を示すリストです。

以上です。