本ポストはこちらの続きです。
Rで有限要素法:片持ち梁-中間点に集中荷重:たわみ量
const typesetMath = (el) => { if (window.MathJax) { // MathJax Typeset window.MathJax.typeset(); } else if (window.katex...
www.saecanet.com
2026.09.09
本ポストではこちらのポストのコード中の
Rで有限要素法:片持ち梁-自由端集中荷重:たわみ量
const typesetMath = (el) => { if (window.MathJax) { // MathJax Typeset window.MathJax.typeset(); } else if (window.katex...
www.saecanet.com
2026.08.30
# ------------------------------------------------------------
# 3. メッシュ生成と全体行列への組み込み
# ------------------------------------------------------------
n_elem <- 5L
Le <- l / n_elem
n_node <- n_elem + 1L
n_dof <- 2L * n_node
K <- matrix(0, n_dof, n_dof)
F <- numeric(n_dof)
K_e <- beam_element_k(Le) # 等間隔分割なので要素剛性行列は全要素で共通
for (e in seq_len(n_elem)) {
# 要素eの自由度: 節点e, e+1 の (w, θ) 計4自由度
dofs <- c(2 * e - 1, 2 * e, 2 * e + 1, 2 * e + 2)
K[dofs, dofs] <- K[dofs, dofs] + K_e
}のコード各行の意味と全体行列 K について確認します。
コード各行の解説
n_elem <- 5L
# 梁全体をいくつの要素に分割するかを指定する(ここでは5分割)
Le <- l / n_elem
# 支間長さlをn_elem等分した、1要素あたりの長さを計算する(この例では2.0/5=0.4)
n_node <- n_elem + 1L
# 要素数がn_elem個のとき、節点の数はその数より1つ多くなる(要素5個なら節点6個)
n_dof <- 2L * n_node
# 全体の自由度数を計算する。1節点あたり自由度2(たわみw、傾きθ)なので、
# 節点数の2倍が全体の自由度数になる(節点6個なら自由度12個)
K <- matrix(0, n_dof, n_dof)
# 全体剛性行列Kの入れ物を用意する。サイズはn_dof×n_dof(ここでは12×12)で、
# 中身はまだ何も組み込まれていないのですべてゼロで初期化する
F <- numeric(n_dof)
# 全体荷重ベクトルFの入れ物を、n_dof個のゼロで用意する(後の荷重載荷で使う)
K_e <- beam_element_k(Le)
# 要素剛性行列を1回だけ計算しておく。等間隔分割(全要素の長さLeが同じ)なので、
# どの要素についても中身は共通のため、ループの外で1回計算すれば十分
for (e in seq_len(n_elem)) {
# 要素番号eを1からn_elem(この例では1から5)まで1つずつ変えながら、以下の処理を繰り返す
dofs <- c(2 * e - 1, 2 * e, 2 * e + 1, 2 * e + 2)
# 要素eが関係する、全体の自由度番号を求める。
# 要素eは節点e・節点e+1からできており、それぞれが(w, θ)の2自由度を持つため、
# 全体の自由度番号としては [2e-1, 2e, 2e+1, 2e+2] の4つになる
# (例: e=1なら[1,2,3,4]、e=2なら[3,4,5,6]、…というように、隣り合う要素どうしで自由度が2つ重なる)
K[dofs, dofs] <- K[dofs, dofs] + K_e
# 要素eの剛性行列K_e(4×4)を、全体剛性行列Kの該当する位置(dofsで指定した行・列)に足し込む。
# 「上書き」ではなく「足し算」であることがポイントで、隣り合う要素が共有する節点の自由度では、
# 両方の要素からの寄与が重なって加算される
}
cat("--- 1列目から6列目まで ---\n")
K[, 1:6]
cat("\n--- 7列目から12列目まで ---\n")
K[, 7:12]--- 1列目から6列目まで ---
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 164063812 32812762 -164063812 32812762 0 0
[2,] 32812762 8750070 -32812762 4375035 0 0
[3,] -164063812 -32812762 328127625 0 -164063812 32812762
[4,] 32812762 4375035 0 17500140 -32812762 4375035
[5,] 0 0 -164063812 -32812762 328127625 0
[6,] 0 0 32812762 4375035 0 17500140
[7,] 0 0 0 0 -164063812 -32812762
[8,] 0 0 0 0 32812762 4375035
[9,] 0 0 0 0 0 0
[10,] 0 0 0 0 0 0
[11,] 0 0 0 0 0 0
[12,] 0 0 0 0 0 0
--- 7列目から12列目まで ---
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 0 0 0 0
[2,] 0 0 0 0 0 0
[3,] 0 0 0 0 0 0
[4,] 0 0 0 0 0 0
[5,] -164063812 32812762 0 0 0 0
[6,] -32812762 4375035 0 0 0 0
[7,] 328127625 0 -164063812 32812762 0 0
[8,] 0 17500140 -32812762 4375035 0 0
[9,] -164063812 -32812762 328127625 0 -164063812 32812762
[10,] 32812762 4375035 0 17500140 -32812762 4375035
[11,] 0 0 -164063812 -32812762 164063812 -32812762
[12,] 0 0 32812762 4375035 -32812762 8750070全体行列 K の出力の見方
出力された12×12行列は、5つの要素剛性行列(各4×4)を、共有する節点の自由度のところで重ね合わせて組み立てた結果です。
行・列の対応関係
行(列)の番号は、自由度の並び [w1, θ1, w2, θ2, w3, θ3, w4, θ4, w5, θ5, w6, θ6] に対応しています(奇数番目が各節点のたわみ\(w\)、偶数番目が傾き\(\theta\))。
対角線上のブロック
- 1〜2行・列(節点1)と11〜12行・列(節点6):片方の要素だけが関わる自由度なので、値は
K_eの対応する成分(例えば164063812など)がそのまま入っています。 - 3〜4行・列(節点2)、5〜6行・列(節点3)、…、9〜10行・列(節点5):2つの要素が共有する節点なので、両方の要素からの寄与が加算され、値がちょうど2倍(328127625 = 164063812×2 など)になっています。
ゼロになっている部分
例えば [1,5](節点1と節点3の間)のようにゼロになっている成分は、節点1と節点3を直接つなぐ要素が存在しない(間に節点2を挟んでいる)ことを表しています。
梁要素は隣り合う節点としかつながらないため、離れた節点どうしの成分は基本的にゼロのままです。
この「隣接した節点だけが値を持ち、それ以外はゼロ」という構造(帯行列)は、有限要素法の全体剛性行列に共通する特徴です。
対称性
K[3,4]=32812762 と K[4,3]=32812762 のように、行列全体が対称になっていることも確認できます。
これは、剛性行列が「ある自由度に力を加えたときに、他の自由度にどれだけ変位が生じるか」という関係から作られており、剛性行列の \((i, j)\) 成分は、「自由度 \(j\) に単位変位を与えたときに、自由度 \(i\) に生じる力」を表していますので、相反定理により、
- 「自由度 \(j\) を動かして自由度 \(i\) に生じる力」
- 「自由度 \(i\) を動かして自由度 \(j\) に生じる力」
は必ず等しくなるため、\((i,j)\) 成分と \((j,i)\) 成分が常に一致し、行列全体が対称になります。
以上です。
