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

# ============================================================
# 片持ち梁のたわみ:FEM(オイラー・ベルヌーイ梁要素) vs 厳密解
# ============================================================
#
# 厳密解 (x は自由端Aから測った距離):
#   w(x)     = (P l^3)/(6EI) * { 2 - 3(x/l) + (x/l)^3 }
#   w_max    = (P l^3)/(3EI)          [x = 0, 自由端]
#   θ(x)     = (P l^2)/(2EI) * { 1 - (x/l)^2 }
#   θ_A      = (P l^2)/(2EI)          [x = 0, 自由端]
#
# 変数定義 :
#   l   : 梁の長さ (支間) [m]
#   P   : 自由端にかける集中荷重 [N]
#   E   : ヤング率 [Pa]
#   I   : 断面二次モーメント [m^4]
#   EI  : 曲げ剛性
#   s   : FEM側の座標。固定端(B)を原点 s=0、自由端(A)を s=l とする
#         (厳密解の x は自由端起点のため、x = l - s の関係になる)

library(ggplot2)

# ------------------------------------------------------------
# 1. 梁の諸元と荷重の設定
# ------------------------------------------------------------
l <- 2.0 # 支間長さ [m]
P <- 1000 # 自由端の集中荷重 [N]
E <- 210e9 # ヤング率(鋼材の代表値)[Pa]
I <- 4.1667e-6 # 断面二次モーメント [m^4] (例: 50mm x 100mm矩形断面)
EI <- E * I

# ------------------------------------------------------------
# 2. 梁要素の剛性行列
# ------------------------------------------------------------
# 節点自由度の並び: [w1, θ1, w2, θ2]

beam_element_k <- function(Le) {
  (EI / Le^3) * matrix(c(
    12, 6 * Le, -12, 6 * Le,
    6 * Le, 4 * Le^2, -6 * Le, 2 * Le^2,
    -12, -6 * Le, 12, -6 * Le,
    6 * Le, 2 * Le^2, -6 * Le, 4 * Le^2
  ), nrow = 4, byrow = TRUE)
}

# ------------------------------------------------------------
# 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
}

# ------------------------------------------------------------
# 4. 荷重と境界条件
# ------------------------------------------------------------
# 固定端(節点1, s=0): w=0, θ=0
# 自由端(最後の節点, s=l): 集中荷重 P

F[n_dof - 1] <- F[n_dof - 1] + P # 最終節点のw自由度に荷重を載荷

fixed_dofs <- c(1, 2)
free_dofs <- setdiff(seq_len(n_dof), fixed_dofs)

# ------------------------------------------------------------
# 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 # 厳密解の座標系(自由端起点)に変換

# ------------------------------------------------------------
# 6. 厳密解の計算
# ------------------------------------------------------------
w_exact_fn <- function(x) (P * l^3) / (6 * EI) * (2 - 3 * (x / l) + (x / l)^3)

# ------------------------------------------------------------
# 7. 節点値での比較
# ------------------------------------------------------------
comparison <- data.frame(
  節点位置_自由端起点_x_m = round(x_nodes, 3),
  FEM近似解_mm = round(w_fem * 1000, 8),
  厳密解_mm = round(w_exact_fn(x_nodes) * 1000, 8),
  差_mm = w_fem * 1000 - w_exact_fn(x_nodes) * 1000
)
print(comparison)

cat("\n自由端でのたわみ量:\n")
cat("  FEM近似解:", w_fem[n_node] * 1000, "mm\n")
cat("  厳密解    :", P * l^3 / (3 * EI) * 1000, "mm  (公式: w_max = P l^3 / 3EI)\n")

# ------------------------------------------------------------
# 8. 要素内部も含めた連続的な可視化(3次エルミート補間で描画)
# ------------------------------------------------------------
hermite_w <- function(xi, Le, w1, th1, w2, th2) {
  N1 <- 1 - 3 * xi^2 + 2 * xi^3
  N2 <- Le * (xi - 2 * xi^2 + xi^3)
  N3 <- 3 * xi^2 - 2 * xi^3
  N4 <- Le * (-xi^2 + xi^3)
  N1 * w1 + N2 * th1 + N3 * w2 + N4 * th2
}

fem_curve <- do.call(rbind, lapply(seq_len(n_elem), function(e) {
  xi_seq <- seq(0, 1, length.out = 20)
  w1 <- w_fem[e]
  th1 <- th_fem[e]
  w2 <- w_fem[e + 1]
  th2 <- th_fem[e + 1]
  s_local <- s_nodes[e] + xi_seq * Le
  data.frame(
    x = l - s_local,
    w_mm = hermite_w(xi_seq, Le, w1, th1, w2, th2) * 1000
  )
}))

x_fine <- seq(0, l, length.out = 200)
exact_curve <- data.frame(x = x_fine, w_mm = w_exact_fn(x_fine) * 1000)

p <- ggplot() +
  geom_line(data = exact_curve, aes(x = x, y = w_mm, color = "厳密解"), linewidth = 1.2) +
  geom_line(
    data = fem_curve, aes(x = x, y = w_mm, color = "FEM(要素内補間)"),
    linetype = "dashed", linewidth = 0.8
  ) +
  geom_point(
    data = data.frame(x = x_nodes, w_mm = w_fem * 1000),
    aes(x = x, y = w_mm, color = "FEM(節点)"), size = 2
  ) +
  scale_color_manual(values = c(
    "厳密解" = "gray40", "FEM(要素内補間)" = "steelblue",
    "FEM(節点)" = "firebrick"
  ), name = "") +
  labs(
    title = "片持ち梁のたわみ曲線:FEM(梁要素) vs 厳密解",
    subtitle = paste0("要素数 = ", n_elem, "  (両者はほぼ完全に重なる)"),
    x = "自由端からの距離 x [m]", y = "たわみ量 w [mm]"
  ) +
  theme_minimal(base_size = 12)

print(p)