Rで有限要素法:4辺単純支持コンクリートスラブ:たわみ量

# ============================================================
# 4辺単純支持・中央集中荷重を受ける矩形板のたわみ:厳密解(Navierの解)
# ============================================================
#
# 支配方程式 :
#   D ∇^4 w = q(x, y),   D = E h^3 / (12 (1 - ν^2))
#
# 境界条件 :
#   w = 0,  ∂^2w/∂n^2 = 0   (境界上のすべての点で)
#
# 厳密解 :
#   w(x,y) = (4P)/(π^4 D a b) *
#            Σ_m Σ_n [ sin(mπ/2) sin(nπ/2) / (m^2/a^2 + n^2/b^2)^2 ] *
#                     sin(mπx/a) sin(nπy/b)
#
# 変数定義 :
#   E   : ヤング率  [Pa]
#   h   : 板(床)の厚さ  [m]
#   nu  : ポアソン比
#   D   : 板の曲げ剛性  [N・m]
#   a,b : x方向・y方向の板の辺の長さ  [m]
#   P   : 中央にかける集中荷重  [N]
#   w   : たわみ量  [m]

library(ggplot2)

# ------------------------------------------------------------
# 1. コンクリート床の物性値・寸法・荷重の設定
# ------------------------------------------------------------
E <- 25e9 # ヤング率: コンクリートの代表値 25 GPa
h <- 0.2 # 床の厚さ: 0.2 m
nu <- 0.2 # ポアソン比: コンクリートの代表値
a <- 5 # 床のスパン(x方向): 5 m
b <- 5 # 床のスパン(y方向): 5 m
P <- 50000 # 中央集中荷重: 50 kN

D <- E * h^3 / (12 * (1 - nu^2))
cat("曲げ剛性 D =", D, "N・m\n")

# ------------------------------------------------------------
# 2. Navierの二重sin級数によるたわみの計算
# ------------------------------------------------------------
# 級数は m, n が奇数の項のみ寄与する(sin(mπ/2)が偶数で0になるため)。
# n_terms 項(奇数番号のみ)まで足し合わせて近似する。

n_terms <- 40L # 級数の打ち切り項数(奇数のみ n_terms 個使用)
m_vals <- seq(1, by = 2, length.out = n_terms) # 1, 3, 5, ...
n_vals <- seq(1, by = 2, length.out = n_terms)

deflection_navier <- function(x, y) {
  total <- 0
  for (m in m_vals) {
    for (n in n_vals) {
      # sin(mπ/2), sin(nπ/2) は奇数の m, n に対して常に ±1
      sm <- sin(m * pi / 2)
      sn <- sin(n * pi / 2)
      denom <- (m^2 / a^2 + n^2 / b^2)^2
      total <- total + sm * sn / denom * sin(m * pi * x / a) * sin(n * pi * y / b)
    }
  }
  (4 * P) / (pi^4 * D * a * b) * total
}

# ------------------------------------------------------------
# 3. 格子点上でたわみを評価
# ------------------------------------------------------------
n_grid <- 21L
x_seq <- seq(0, a, length.out = n_grid)
y_seq <- seq(0, b, length.out = n_grid)

grid_df <- expand.grid(x = x_seq, y = y_seq)
grid_df$w <- mapply(deflection_navier, grid_df$x, grid_df$y)

w_center <- deflection_navier(a / 2, b / 2)
cat("中央でのたわみ量 (厳密解): ", w_center * 1000, "mm\n")

# ------------------------------------------------------------
# 4. 可視化
# ------------------------------------------------------------
p <- ggplot(grid_df, aes(x = x, y = y, fill = w * 1000)) +
  geom_tile() +
  scale_fill_viridis_c(name = "たわみ量\n[mm]") +
  coord_equal() +
  labs(
    title = "4辺単純支持コンクリート床のたわみ(厳密解:Navierの級数解)",
    subtitle = paste0("中央集中荷重 P = ", P / 1000, " kN, 床厚 h = ", h, " m"),
    x = "x [m]", y = "y [m]"
  ) +
  theme_minimal(base_size = 12)

print(p)

# 後続のFEMスクリプトとの比較用に、厳密解の中央たわみ量を保存しておく
saveRDS(
  list(w_center = w_center, D = D, a = a, b = b, P = P, E = E, h = h, nu = nu),
  "plate_exact_reference.rds"
)

# ============================================================
# 4辺単純支持・中央集中荷重を受ける矩形板のたわみ:FEM近似解
# (ACM: Adini-Clough-Melosh 非適合長方形板要素)
# ============================================================
#
# 各節点の自由度 : [w, ∂w/∂x, ∂w/∂y]  (計3自由度)
# 要素内のたわみは12項の3次多項式で近似する:
#   w(x,y) = a1 + a2 x + a3 y + a4 x^2 + a5 xy + a6 y^2
#            + a7 x^3 + a8 x^2 y + a9 x y^2 + a10 y^3
#            + a11 x^3 y + a12 x y^3
# (4節点 × 3自由度 = 12 個の未知係数 a1〜a12 に対応する)
#
# ひずみエネルギー :
#   U = (D/2) ∫∫ [ w_xx^2 + 2ν w_xx w_yy + w_yy^2 + 2(1-ν) w_xy^2 ] dA
# ここで w_xx = ∂^2w/∂x^2 等。この2次形式から要素剛性行列を導く。

# 前スクリプトで保存した厳密解と物性値を読み込む
ref <- readRDS("plate_exact_reference.rds")
E <- ref$E
h <- ref$h
nu <- ref$nu
a_len <- ref$a
b_len <- ref$b
P <- ref$P
D <- ref$D

# ------------------------------------------------------------
# 1. 構造格子メッシュの生成
# ------------------------------------------------------------
nx <- 21L # x方向節点数(奇数分割になるよう nx-1 を偶数にし、中心に節点を置く)
ny <- 21L
hx <- a_len / (nx - 1)
hy <- b_len / (ny - 1)
ae <- hx / 2 # 要素の半幅 (element half-width, x方向)
be <- hy / 2 # 要素の半幅 (element half-width, y方向)

nodes <- expand.grid(i = 0:(nx - 1), j = 0:(ny - 1)) |>
  transform(x = i * hx, y = j * hy) |>
  subset(select = c(x, y))
nodes <- as.matrix(nodes)
n_node <- nrow(nodes)

node_id <- function(i, j) j * nx + i + 1L

# ------------------------------------------------------------
# 2. 単項式基底とその微分
# ------------------------------------------------------------
# 要素中心を原点とするローカル座標 (x, y) について定義する。

P_vec <- function(x, y) c(1, x, y, x^2, x * y, y^2, x^3, x^2 * y, x * y^2, y^3, x^3 * y, x * y^3)
Px_vec <- function(x, y) c(0, 1, 0, 2 * x, y, 0, 3 * x^2, 2 * x * y, y^2, 0, 3 * x^2 * y, y^3)
Py_vec <- function(x, y) c(0, 0, 1, 0, x, 2 * y, 0, x^2, 2 * x * y, 3 * y^2, x^3, 3 * x * y^2)
Pxx_vec <- function(x, y) c(0, 0, 0, 2, 0, 0, 6 * x, 2 * y, 0, 0, 6 * x * y, 0)
Pyy_vec <- function(x, y) c(0, 0, 0, 0, 0, 2, 0, 0, 2 * x, 6 * y, 0, 6 * x * y)
Pxy_vec <- function(x, y) c(0, 0, 0, 0, 1, 0, 0, 2 * x, 2 * y, 0, 3 * x^2, 3 * y^2)

# 要素の4隅のローカル座標(要素中心が原点)
corners <- list(c(-ae, -be), c(ae, -be), c(ae, be), c(-ae, be))

# ------------------------------------------------------------
# 3. 形状関数を「節点値に一致する」条件から構成する
# ------------------------------------------------------------
# 12x12行列 C の各行は、各節点での w, ∂w/∂x, ∂w/∂y の値を
# 単項式係数 a1〜a12 で表したもの。
# 「多項式の係数 a = C^{-1} × 節点自由度ベクトル d」という関係により、
# 形状関数はハードコードせず、行列演算から自動的に求まる。

C <- matrix(0, 12, 12)
for (k in seq_along(corners)) {
  cxy <- corners[[k]]
  C[3 * (k - 1) + 1, ] <- P_vec(cxy[1], cxy[2])
  C[3 * (k - 1) + 2, ] <- Px_vec(cxy[1], cxy[2])
  C[3 * (k - 1) + 3, ] <- Py_vec(cxy[1], cxy[2])
}
C_inv <- solve(C)

# 曲げ剛性マトリクス D_b
D_b <- D * matrix(c(1, nu, 0, nu, 1, 0, 0, 0, 2 * (1 - nu)), nrow = 3, byrow = TRUE)

# ------------------------------------------------------------
# 4. 要素剛性行列の計算(3点ガウス求積)
# ------------------------------------------------------------
gauss_pts <- c(-sqrt(3 / 5), 0, sqrt(3 / 5))
gauss_wts <- c(5 / 9, 8 / 9, 5 / 9)

element_stiffness <- function() {
  Ke <- matrix(0, 12, 12)
  for (gi in seq_along(gauss_pts)) {
    for (gj in seq_along(gauss_pts)) {
      x <- gauss_pts[gi] * ae
      y <- gauss_pts[gj] * be
      jac <- ae * be
      B <- rbind(Pxx_vec(x, y), Pyy_vec(x, y), Pxy_vec(x, y)) %*% C_inv # 3x12
      Ke <- Ke + t(B) %*% D_b %*% B * gauss_wts[gi] * gauss_wts[gj] * jac
    }
  }
  Ke
}
K_e <- element_stiffness() # 全要素で形状・寸法が同じなので、1回だけ計算すればよい

# ------------------------------------------------------------
# 5. 全体行列への組み込み
# ------------------------------------------------------------
n_dof <- 3L * n_node
K <- matrix(0, n_dof, n_dof)
F <- numeric(n_dof)

dof_of <- function(node) c(3 * (node - 1) + 1, 3 * (node - 1) + 2, 3 * (node - 1) + 3)

for (j in 0:(ny - 2)) {
  for (i in 0:(nx - 2)) {
    n1 <- node_id(i, j)
    n2 <- node_id(i + 1, j)
    n3 <- node_id(i + 1, j + 1)
    n4 <- node_id(i, j + 1)
    gdofs <- c(dof_of(n1), dof_of(n2), dof_of(n3), dof_of(n4))
    K[gdofs, gdofs] <- K[gdofs, gdofs] + K_e
  }
}

# ------------------------------------------------------------
# 6. 荷重の載荷:中央節点への集中荷重
# ------------------------------------------------------------
ic <- (nx - 1) %/% 2
jc <- (ny - 1) %/% 2
center_node <- node_id(ic, jc)
F[3 * (center_node - 1) + 1] <- F[3 * (center_node - 1) + 1] + P

# ------------------------------------------------------------
# 7. 境界条件:単純支持(縁でたわみ w = 0、回転は自由)
# ------------------------------------------------------------
is_boundary <- nodes[, 1] == 0 | nodes[, 1] == a_len | nodes[, 2] == 0 | nodes[, 2] == b_len
boundary_nodes <- which(is_boundary)
fixed_dofs <- 3 * (boundary_nodes - 1) + 1 # 各境界節点の w-自由度のみ拘束

free_dofs <- setdiff(seq_len(n_dof), fixed_dofs)

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

w_h <- u[seq(1, n_dof, by = 3)] # 各節点のたわみ自由度のみ取り出す

w_center_fem <- w_h[center_node]
cat("中央でのたわみ量 (FEM近似解):", w_center_fem * 1000, "mm\n")
cat("中央でのたわみ量 (厳密解)  :", ref$w_center * 1000, "mm\n")
cat("相対誤差:", (w_center_fem - ref$w_center) / ref$w_center * 100, "%\n")

# ------------------------------------------------------------
# 9. 可視化と比較
# ------------------------------------------------------------
plot_df <- data.frame(x = nodes[, 1], y = nodes[, 2], w_mm = w_h * 1000)

p <- ggplot(plot_df, aes(x = x, y = y, fill = w_mm)) +
  geom_tile() +
  scale_fill_viridis_c(name = "たわみ量\n[mm]") +
  coord_equal() +
  labs(
    title = "4辺単純支持コンクリート床のたわみ(FEM近似解:ACM板要素)",
    subtitle = paste0(
      "節点数: ", nx, "×", ny, ", 中央たわみ誤差: ",
      round((w_center_fem - ref$w_center) / ref$w_center * 100, 2), "%"
    ),
    x = "x [m]", y = "y [m]"
  ) +
  theme_minimal(base_size = 12)

print(p)

# ============================================================
# ACM板要素FEMの収束確認
# ============================================================
#
# 4辺単純支持・中央集中荷重を受けるコンクリート床について、
# メッシュを細かくしていったときにFEM近似解が厳密解(Navierの級数解)に
# どのように収束していくかを確認する。
#
# 変数定義 :
#   nx_list : 検証する節点数(x方向 = y方向)のベクトル
#   w_center_fem : 各メッシュでの中央たわみ量(FEM近似解)
#   w_center_exact : 中央たわみ量の厳密解(Navierの級数解)
#   rel_error : 相対誤差 (w_fem - w_exact) / w_exact

# ------------------------------------------------------------
# 1. 物性値・寸法・荷重の設定(前回と同じコンクリート床)
# ------------------------------------------------------------
E <- 25e9
h <- 0.2
nu <- 0.2
a_len <- 5
b_len <- 5
P <- 50000
D <- E * h^3 / (12 * (1 - nu^2))

# ------------------------------------------------------------
# 2. 厳密解(中央たわみ)の計算 — 前回スクリプトと同じ級数
# ------------------------------------------------------------
n_terms <- 40L
m_vals <- seq(1, by = 2, length.out = n_terms)
n_vals <- seq(1, by = 2, length.out = n_terms)

deflection_navier <- function(x, y) {
  total <- 0
  for (m in m_vals) {
    for (n in n_vals) {
      sm <- sin(m * pi / 2)
      sn <- sin(n * pi / 2)
      denom <- (m^2 / a_len^2 + n^2 / b_len^2)^2
      total <- total + sm * sn / denom * sin(m * pi * x / a_len) * sin(n * pi * y / b_len)
    }
  }
  (4 * P) / (pi^4 * D * a_len * b_len) * total
}
w_center_exact <- deflection_navier(a_len / 2, b_len / 2)
cat("厳密解(中央たわみ):", w_center_exact * 1000, "mm\n")
cat("\n一辺あたりの節点数と厳密解との誤差\n")

# ------------------------------------------------------------
# 3. ACM要素FEMを1関数にまとめる(メッシュ密度を引数化)
# ------------------------------------------------------------
solve_plate_fem <- function(nx, ny) {
  hx <- a_len / (nx - 1)
  hy <- b_len / (ny - 1)
  ae <- hx / 2
  be <- hy / 2

  nodes <- expand.grid(i = 0:(nx - 1), j = 0:(ny - 1)) |>
    transform(x = i * hx, y = j * hy) |>
    subset(select = c(x, y))
  nodes <- as.matrix(nodes)
  n_node <- nrow(nodes)
  node_id <- function(i, j) j * nx + i + 1L

  P_vec <- function(x, y) c(1, x, y, x^2, x * y, y^2, x^3, x^2 * y, x * y^2, y^3, x^3 * y, x * y^3)
  Px_vec <- function(x, y) c(0, 1, 0, 2 * x, y, 0, 3 * x^2, 2 * x * y, y^2, 0, 3 * x^2 * y, y^3)
  Py_vec <- function(x, y) c(0, 0, 1, 0, x, 2 * y, 0, x^2, 2 * x * y, 3 * y^2, x^3, 3 * x * y^2)
  Pxx_vec <- function(x, y) c(0, 0, 0, 2, 0, 0, 6 * x, 2 * y, 0, 0, 6 * x * y, 0)
  Pyy_vec <- function(x, y) c(0, 0, 0, 0, 0, 2, 0, 0, 2 * x, 6 * y, 0, 6 * x * y)
  Pxy_vec <- function(x, y) c(0, 0, 0, 0, 1, 0, 0, 2 * x, 2 * y, 0, 3 * x^2, 3 * y^2)

  corners <- list(c(-ae, -be), c(ae, -be), c(ae, be), c(-ae, be))
  C <- matrix(0, 12, 12)
  for (k in seq_along(corners)) {
    cxy <- corners[[k]]
    C[3 * (k - 1) + 1, ] <- P_vec(cxy[1], cxy[2])
    C[3 * (k - 1) + 2, ] <- Px_vec(cxy[1], cxy[2])
    C[3 * (k - 1) + 3, ] <- Py_vec(cxy[1], cxy[2])
  }
  C_inv <- solve(C)
  D_b <- D * matrix(c(1, nu, 0, nu, 1, 0, 0, 0, 2 * (1 - nu)), nrow = 3, byrow = TRUE)

  gauss_pts <- c(-sqrt(3 / 5), 0, sqrt(3 / 5))
  gauss_wts <- c(5 / 9, 8 / 9, 5 / 9)

  K_e <- matrix(0, 12, 12)
  for (gi in seq_along(gauss_pts)) {
    for (gj in seq_along(gauss_pts)) {
      x <- gauss_pts[gi] * ae
      y <- gauss_pts[gj] * be
      jac <- ae * be
      B <- rbind(Pxx_vec(x, y), Pyy_vec(x, y), Pxy_vec(x, y)) %*% C_inv
      K_e <- K_e + t(B) %*% D_b %*% B * gauss_wts[gi] * gauss_wts[gj] * jac
    }
  }

  n_dof <- 3L * n_node
  K <- matrix(0, n_dof, n_dof)
  F <- numeric(n_dof)
  dof_of <- function(node) c(3 * (node - 1) + 1, 3 * (node - 1) + 2, 3 * (node - 1) + 3)

  for (j in 0:(ny - 2)) {
    for (i in 0:(nx - 2)) {
      n1 <- node_id(i, j)
      n2 <- node_id(i + 1, j)
      n3 <- node_id(i + 1, j + 1)
      n4 <- node_id(i, j + 1)
      gdofs <- c(dof_of(n1), dof_of(n2), dof_of(n3), dof_of(n4))
      K[gdofs, gdofs] <- K[gdofs, gdofs] + K_e
    }
  }

  ic <- (nx - 1) %/% 2
  jc <- (ny - 1) %/% 2
  center_node <- node_id(ic, jc)
  F[3 * (center_node - 1) + 1] <- F[3 * (center_node - 1) + 1] + P

  is_boundary <- nodes[, 1] == 0 | nodes[, 1] == a_len | nodes[, 2] == 0 | nodes[, 2] == b_len
  boundary_nodes <- which(is_boundary)
  fixed_dofs <- 3 * (boundary_nodes - 1) + 1
  free_dofs <- setdiff(seq_len(n_dof), fixed_dofs)

  u <- numeric(n_dof)
  u[free_dofs] <- solve(K[free_dofs, free_dofs], F[free_dofs])

  u[3 * (center_node - 1) + 1] # 中央節点のたわみ自由度を返す
}

# ------------------------------------------------------------
# 4. 複数のメッシュ密度で実行し、誤差の推移を記録
# ------------------------------------------------------------
# 注意: nx は奇数にすること(中心に節点を置くため、nx-1が偶数である必要がある)

nx_list <- c(5L, 9L, 15L, 21L, 27L, 33L, 41L)

results <- data.frame(nx = nx_list, w_center_fem = NA_real_)
for (k in seq_along(nx_list)) {
  nx <- nx_list[k]
  results$w_center_fem[k] <- solve_plate_fem(nx, nx)
  cat("nx =", nx, " -> w_center =", results$w_center_fem[k] * 1000, "mm\n")
}

results <- results |>
  transform(
    h_mesh = a_len / (nx - 1), # メッシュの代表寸法(要素1辺の長さ)
    rel_error_pct = (w_center_fem - w_center_exact) / w_center_exact * 100
  )

cat("\n要素の代表寸法(h)と厳密解との誤差\n")
print(results)

# ------------------------------------------------------------
# 5. 可視化(相対誤差 vs メッシュサイズ、両対数プロット)
# ------------------------------------------------------------
p1 <- ggplot(results, aes(x = h_mesh, y = abs(rel_error_pct))) +
  geom_line(color = "steelblue") +
  geom_point(size = 2, color = "steelblue") +
  scale_x_log10() +
  scale_y_log10() +
  labs(
    title = "ACM板要素FEMの収束確認",
    subtitle = "要素サイズを小さくするほど、厳密解への誤差が減少する",
    x = "要素の代表寸法 h [m](対数軸)",
    y = "中央たわみの相対誤差の絶対値 [%](対数軸)"
  ) +
  theme_minimal(base_size = 12)

p2 <- ggplot(results, aes(x = nx, y = w_center_fem * 1000)) +
  geom_line(color = "darkorange") +
  geom_point(size = 2, color = "darkorange") +
  geom_hline(yintercept = w_center_exact * 1000, linetype = "dashed", color = "gray30") +
  annotate("text",
    x = max(nx_list) * 0.1, y = w_center_exact * 1000,
    label = "厳密解", vjust = -1, color = "gray30"
  ) +
  labs(
    title = "中央たわみ量の収束(FEM vs 厳密解)",
    x = "一辺あたりの節点数 nx",
    y = "中央たわみ量 [mm]"
  ) +
  theme_minimal(base_size = 12)

print(p1)
print(p2)