1 このファイルの目的

3.9「平板の二次元曲げ（たわみ）」 は、式が多く見えますが、中心はかなり単純です。

薄い板に荷重がかかると、板は曲がる。
その曲がりにくさは 曲げ剛性 \(D\) で決まり、
たわみ \(w(x)\) は 4階微分方程式で表される。

最終的に重要なのは、次の式です。

\[ D\frac{d^4 w}{dx^4} = q(x)-P\frac{d^2w}{dx^2}. \]

ここで、

です。

まずは水平方向の力 \(P\) を無視して、

\[ D\frac{d^4 w}{dx^4}=q(x) \]

を理解します。

2 何を表している式か

2.1 たわみ \(w(x)\)

板を横から見て、もとの水平な状態からどれだけ下に変位したかを \(w(x)\) とします。

この図では、点線がもとの板、曲線がたわんだ板です。

2.2 曲げ剛性 \(D\)

曲げ剛性は

\[ D=\frac{Eh^3}{12(1-\nu^2)} \]

です。

ここで、

  • \(E\)：ヤング率
  • \(h\)：板の厚さ
  • \(\nu\)：ポアソン比

です。

特に重要なのは

\[ D \propto h^3 \]

という点です。

つまり、板の厚さが2倍になると、曲げ剛性は8倍になります。厚い板は非常に曲がりにくくなります。

E <- 100e9      # Pa
nu <- 0.25
h <- seq(1, 10, by = 0.1)  # m
D <- E * h^3 / (12 * (1 - nu^2))

plot(h, D, type = "l", lwd = 2,
     xlab = "plate thickness h (m)",
     ylab = "flexural rigidity D (N m)",
     main = "曲げ剛性 D は厚さ h の3乗に比例する")

3 式の直感

\[ D\frac{d^4 w}{dx^4}=q(x) \]

は、少し見慣れない式ですが、意味は次のように考えられます。

地球科学では、リソスフェアを薄い弾性板とみなし、火山島・海山・堆積物などの荷重によるたわみを考えるときに使います。

4 簡単なシミュレーション

ここでは、長さ \(L\) の板の両端を単純支持とします。

単純支持とは、

\[ w(0)=w(L)=0 \]

かつ、端で曲げモーメントがゼロ、すなわち

\[ w''(0)=w''(L)=0 \]

という条件です。

この条件では、たわみを正弦波の和で表すと便利です。

\[ w(x)=\sum_{n=1}^{N} a_n \sin\left(\frac{n\pi x}{L}\right). \]

同様に荷重も

\[ q(x)=\sum_{n=1}^{N} q_n \sin\left(\frac{n\pi x}{L}\right) \]

と展開すると、

\[ a_n= \frac{q_n}{D(n\pi/L)^4} \]

となります。

つまり、荷重を波に分け、それぞれの波長に対してたわみを計算して足し合わせます。

5 Rで使う関数

flexural_rigidity <- function(E, h, nu) {
  E * h^3 / (12 * (1 - nu^2))
}

beam_deflection <- function(x, q, L, D, N = 100, P = 0) {
  # Simply supported beam/plate strip.
  # Solve: D * fourth_derivative(w) = q - P * second_derivative(w)
  # Equivalent modal relation:
  # a_n = q_n / (D k_n^4 - P k_n^2)

  dx <- x[2] - x[1]
  w <- rep(0, length(x))

  for (n in 1:N) {
    k <- n * pi / L
    basis <- sin(k * x)

    # sine coefficient of q(x)
    qn <- (2 / L) * sum(q * basis) * dx

    denom <- D * k^4 - P * k^2

    # avoid numerical instability near buckling
    if (abs(denom) < 1e-30) next

    an <- qn / denom
    w <- w + an * basis
  }

  w
}

6 例1：中央付近に荷重がある場合

点荷重そのものは数学的にはデルタ関数になりますが、ここでは扱いやすいように、中央付近に幅をもったガウス型の荷重を置きます。

L <- 100e3       # 100 km
E <- 70e9        # 70 GPa
nu <- 0.25
h <- 10e3        # 10 km
D <- flexural_rigidity(E, h, nu)

x <- seq(0, L, length.out = 1001)

q0 <- 1e7        # N/m^2 相当の荷重スケール
sigma <- 8e3     # 荷重の幅
q <- q0 * exp(-(x - L/2)^2 / (2 * sigma^2))

w <- beam_deflection(x, q, L, D, N = 150)

plot(x / 1000, q, type = "l", lwd = 2,
     xlab = "x (km)", ylab = "load q(x)",
     main = "中央付近に置いた荷重")

plot(x / 1000, w, type = "l", lwd = 2,
     xlab = "x (km)", ylab = "deflection w (m)",
     main = "荷重によるたわみ")
abline(h = 0, lty = 2)

この例では、中央付近に荷重があるため、板は中央付近で大きくたわみます。

7 例2：板の厚さを変える

次に、板の厚さ \(h\) を変えてみます。

曲げ剛性は \(h^3\) に比例するため、厚さを少し変えるだけでも、たわみは大きく変わります。

hs <- c(5e3, 10e3, 20e3)  # 5, 10, 20 km
w_list <- list()

for (hi in hs) {
  Di <- flexural_rigidity(E, hi, nu)
  w_list[[as.character(hi)]] <- beam_deflection(x, q, L, Di, N = 150)
}

plot(x / 1000, w_list[[1]], type = "l", lwd = 2,
     xlab = "x (km)", ylab = "deflection w (m)",
     main = "板の厚さによるたわみの違い",
     ylim = range(unlist(w_list)))

for (i in 2:length(hs)) {
  lines(x / 1000, w_list[[i]], lwd = 2, lty = i)
}

legend("topright",
       legend = paste0("h = ", hs / 1000, " km"),
       lwd = 2,
       lty = 1:length(hs),
       bty = "n")

厚さが 5 km から 10 km になると、曲げ剛性は \(2^3=8\) 倍になります。したがって、同じ荷重に対するたわみはかなり小さくなります。

8 例3：水平方向圧縮力 \(P\) の効果

添付資料の式では、水平方向の力 \(P\) も入っています。

\[ D\frac{d^4 w}{dx^4} = q(x)-P\frac{d^2w}{dx^2}. \]

これは

\[ D\frac{d^4 w}{dx^4} + P\frac{d^2w}{dx^2} = q(x) \]

と書いても同じです。

圧縮力 \(P\) があると、板は座屈しやすくなり、同じ荷重でもたわみが大きくなります。

単純支持の板では、最初の座屈荷重はおおよそ

\[ P_{\rm cr}=D\left(\frac{\pi}{L}\right)^2 \]

です。

ここでは \(P\) をその一部として与えてみます。

D <- flexural_rigidity(E, h, nu)
Pcr <- D * (pi / L)^2

P_values <- c(0, 0.3 * Pcr, 0.6 * Pcr, 0.85 * Pcr)
wP_list <- list()

for (Pi in P_values) {
  wP_list[[as.character(Pi)]] <- beam_deflection(x, q, L, D, N = 150, P = Pi)
}

plot(x / 1000, wP_list[[1]], type = "l", lwd = 2,
     xlab = "x (km)", ylab = "deflection w (m)",
     main = "水平圧縮力 P によるたわみの増大",
     ylim = range(unlist(wP_list)))

for (i in 2:length(P_values)) {
  lines(x / 1000, wP_list[[i]], lwd = 2, lty = i)
}

legend("topright",
       legend = paste0("P/Pcr = ", round(P_values / Pcr, 2)),
       lwd = 2,
       lty = 1:length(P_values),
       bty = "n")

\(P\) が大きくなると、分母

\[ Dk^4-Pk^2 \]

が小さくなるため、たわみが大きくなります。

これは、圧縮されている板がより曲がりやすくなることを表しています。

9 ここまでのまとめ

3.9節のポイントは、次のように整理できます。

  1. リソスフェアのような薄い弾性板は、荷重を受けるとたわむ。
  2. たわみは \(w(x)\) で表す。
  3. 板の曲がりにくさは曲げ剛性 \[ D=\frac{Eh^3}{12(1-\nu^2)} \] で決まる。
  4. 厚さ \(h\) が重要で、\(D\propto h^3\) である。
  5. 基本方程式は \[ D\frac{d^4 w}{dx^4}=q(x)-P\frac{d^2w}{dx^2} \] である。
  6. \(P=0\) なら、荷重 \(q(x)\) に対して \[ D\frac{d^4w}{dx^4}=q(x) \] となる。
  7. 圧縮力 \(P\) があると、板はよりたわみやすくなる。

10 地球科学的な読み替え

この式を地球科学に対応させると、次のようになります。

数学の記号 地球科学での意味
\(w(x)\) リソスフェアの沈み込み・隆起
\(q(x)\) 火山島、海山、堆積物などの荷重
\(D\) リソスフェアの曲げ剛性
\(h\) 有効弾性厚
\(P\) 水平方向のテクトニック圧縮力

したがって、観測されたたわみの形を説明できる \(D\) や \(h\) を探すことで、リソスフェアの有効弾性厚を推定できます。

11 おまけ：荷重を変えて遊ぶ

最後に、荷重の幅を変えるとたわみがどう変わるかを見ます。

sigmas <- c(4e3, 8e3, 16e3)
wS_list <- list()

for (si in sigmas) {
  qi <- q0 * exp(-(x - L/2)^2 / (2 * si^2))
  wS_list[[as.character(si)]] <- beam_deflection(x, qi, L, D, N = 150)
}

plot(x / 1000, wS_list[[1]], type = "l", lwd = 2,
     xlab = "x (km)", ylab = "deflection w (m)",
     main = "荷重の幅によるたわみの違い",
     ylim = range(unlist(wS_list)))

for (i in 2:length(sigmas)) {
  lines(x / 1000, wS_list[[i]], lwd = 2, lty = i)
}

legend("topright",
       legend = paste0("sigma = ", sigmas / 1000, " km"),
       lwd = 2,
       lty = 1:length(sigmas),
       bty = "n")

同じ荷重スケールでも、荷重の幅が広いほど、より広い範囲でたわみが生じます。

12 最小限の理解

3.9節は、細部の導出を追うと難しく見えますが、最小限では次の理解で十分です。

\[ \boxed{ D\frac{d^4 w}{dx^4}=q(x)-P\frac{d^2w}{dx^2} } \]

これは、

荷重 \(q\) が板を曲げ、曲げ剛性 \(D\) がそれに抵抗する。
水平圧縮力 \(P\) があると、曲がりやすくなる。

という式です。

地球科学では、リソスフェアを薄い弾性板とみなして、火山島・海山・堆積盆地・沈み込み帯などのたわみを説明するために使います。