Rで有限要素法:等価節点荷重

# ------------------------------------------------------------
# 1. 梁の諸元と荷重の設定
# ------------------------------------------------------------
l <- 2.0
a <- 0.4
b <- 0.8
c <- l - a - b # 区間③の長さ(自動計算)
p <- 600 # 荷重区間での分布荷重の強さ [N/m]
E <- 210e9
I <- 4.1667e-6
EI <- E * I

# 荷重区間のFEM座標(固定端起点)での範囲
s_load_lo <- l - (a + b) # = c
s_load_hi <- l - a # = b + c

# 分布荷重をFEM座標sの関数として定義。区間[s_load_lo, s_load_hi]でのみp、他は0。
q_dist <- function(s) ifelse(s >= s_load_lo & s <= s_load_hi, p, 0)

# ------------------------------------------------------------
# 2. 形状関数と要素剛性行列
# ------------------------------------------------------------
shape_N <- function(xi, Le) {
  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)
  c(N1, N2, N3, N4)
}

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. 等価節点荷重をガウス求積で数値的に計算する(5点ガウス求積)
# ------------------------------------------------------------
gauss_pts <- c(
  0,
  -1 / 3 * sqrt(5 - 2 * sqrt(10 / 7)),
  1 / 3 * sqrt(5 - 2 * sqrt(10 / 7)),
  -1 / 3 * sqrt(5 + 2 * sqrt(10 / 7)),
  1 / 3 * sqrt(5 + 2 * sqrt(10 / 7))
)
gauss_wts <- c(
  128 / 225,
  (322 + 13 * sqrt(70)) / 900,
  (322 + 13 * sqrt(70)) / 900,
  (322 - 13 * sqrt(70)) / 900,
  (322 - 13 * sqrt(70)) / 900
)

beam_element_f <- function(s0, Le) {
  f_e <- numeric(4)
  for (k in seq_along(gauss_pts)) {
    xi <- (gauss_pts[k] + 1) / 2
    jac <- Le / 2
    s <- s0 + xi * Le
    q <- q_dist(s)
    f_e <- f_e + shape_N(xi, Le) * q * jac * gauss_wts[k]
  }
  f_e
}

# ------------------------------------------------------------
# 4. メッシュ生成(荷重区間の境界にちょうど節点が来るように要素数を選ぶ)
# ------------------------------------------------------------
n_elem <- 10L # Le=0.2、s_load_lo=0.8, s_load_hi=1.6 はともにLeの倍数
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)) {
  s0 <- (e - 1) * Le
  f_e <- beam_element_f(s0, Le)
  dofs <- c(2 * e - 1, 2 * e, 2 * e + 1, 2 * e + 2)
  K[dofs, dofs] <- K[dofs, dofs] + K_e
  F[dofs] <- F[dofs] + f_e
}