1 このノートの目的

このノートでは、 Love の応力関数 と Hankel 変換 を、茂木モデルそのものではなく、数学的な道具として丁寧に整理します。

目的は次の3点です。

  1. なぜ Love の応力関数を導入するのか。
  2. なぜ軸対称問題では Hankel 変換が自然に出てくるのか。
  3. Love の応力関数と Hankel 変換を組み合わせると、何が簡単になるのか。

最初に結論を書くと、

\[ \boxed{ \text{Love の応力関数は、軸対称弾性問題を1つのスカラー関数にまとめる道具} } \]

であり、

\[ \boxed{ \text{Hankel 変換は、円筒座標の半径方向微分を波数空間で簡単にする道具} } \]

です。

2 軸対称弾性問題とは何か

円筒座標 \((s,\varphi,z)\) を使います。

軸対称とは、すべての量が \(\varphi\) に依存しない、ということです。

\[ \frac{\partial}{\partial \varphi}=0 \]

また、ねじれがない問題では、

\[ u_\varphi=0 \]

とできます。

したがって変位は

\[ \mathbf{u}=(u_s,0,u_z) \]

です。

このように、軸対称問題では本質的に \(s\)-\(z\) 平面の2次元問題になります。

3 直接解くと何が大変か

外力がない静的な線形弾性体では、変位は

\[ \nabla^2 \mathbf{u} + \frac{1}{1-2\nu} \nabla(\nabla\cdot\mathbf{u})=0 \]

を満たします。

これはベクトル方程式なので、本来は

\[ u_s(s,z), \quad u_z(s,z) \]

を連立して解く必要があります。

しかも、応力境界条件は変位の微分で書かれます。

たとえば地表 \(z=0\) で自由表面なら、

\[ \sigma_{zz}=0,\quad \sigma_{zs}=0 \]

です。

つまり、

という状況になります。

そこで、変位と応力を同時に整理できるスカラー関数を導入します。それが Love の応力関数 です。

4 Love の応力関数の基本形

軸対称でねじれのない問題では、Love の応力関数 \(\phi(s,z)\) を用いて、

\[ u_s = -\frac{1}{2\mu} \frac{\partial^2 \phi}{\partial s \partial z} \]

\[ u_z = \frac{1}{2\mu} \left[ 2(1-\nu)\nabla^2\phi - \frac{\partial^2 \phi}{\partial z^2} \right] \]

と書けます。

ここで、円筒座標の軸対称 Laplacian は

\[ \nabla^2\phi = \frac{1}{s} \frac{\partial}{\partial s} \left( s\frac{\partial \phi}{\partial s} \right) + \frac{\partial^2 \phi}{\partial z^2} \]

です。

この形で変位を定義すると、\(\phi\) は

\[ \nabla^2\nabla^2\phi=0 \]

を満たします。

これは 重調和方程式 です。

つまり、

\[ \boxed{ (u_s,u_z)\text{ を直接解く問題} \quad\Longrightarrow\quad \phi\text{ を解く問題} } \]

に変わります。

5 応力も \(\phi\) から出る

Love の応力関数が便利なのは、変位だけでなく応力も同じ \(\phi\) から出せることです。

代表的に、地表境界条件で重要な2つを書くと、

\[ \sigma_{zz} = \frac{\partial}{\partial z} \left[ (2-\nu)\nabla^2\phi - \frac{\partial^2\phi}{\partial z^2} \right] \]

\[ \sigma_{zs} = \frac{\partial}{\partial s} \left[ (1-\nu)\nabla^2\phi - \frac{\partial^2\phi}{\partial z^2} \right] \]

です。

したがって、地表で

\[ \sigma_{zz}=p(s),\quad \sigma_{zs}=q(s) \]

のような軸対称荷重が与えられたとき、それを \(\phi\) の条件に変換して解くことができます。

6 ここまでのイメージ

Love の応力関数の役割は、図式的には次の通りです。

\[ \phi(s,z) \quad\Longrightarrow\quad u_s,\ u_z \quad\Longrightarrow\quad \sigma_{ss},\sigma_{\varphi\varphi},\sigma_{zz},\sigma_{zs} \]

つまり、\(\phi\) を1つ決めれば、変位と応力がまとめて決まります。

7 しかし、まだ偏微分方程式である

ここまでで、ベクトル方程式はスカラー方程式になりました。

しかし、

\[ \nabla^2\nabla^2\phi=0 \]

はまだ \(s\) と \(z\) に関する偏微分方程式です。

これをさらに簡単にするために使うのが Hankel 変換 です。

8 Fourier 変換との比較

まず、Fourier 変換を思い出します。

デカルト座標で

\[ f(x) \]

を扱うとき、Fourier 変換では

\[ e^{ikx} \]

の重ね合わせで表します。

このとき、

\[ \frac{d^2}{dx^2}e^{ikx}=-k^2e^{ikx} \]

なので、微分が単なる \(-k^2\) の掛け算になります。

これが Fourier 変換の大きな利点です。

9 円筒座標では Bessel 関数が Fourier 波の役割をする

軸対称円筒座標では、半径方向の Laplacian は

\[ \frac{1}{s} \frac{\partial}{\partial s} \left( s\frac{\partial}{\partial s} \right) \]

です。

この演算子に対して、Fourier 波 \(e^{ikx}\) と同じ役割をするのが Bessel 関数 \(J_0(\xi s)\) です。

実際に、

\[ \frac{1}{s} \frac{\partial}{\partial s} \left[ s\frac{\partial}{\partial s}J_0(\xi s) \right] = -\xi^2 J_0(\xi s) \]

が成り立ちます。

つまり、

\[ J_0(\xi s) \]

は、円筒座標の半径方向 Laplacian に対する固有関数です。

これが Hankel 変換で Bessel 関数が出てくる理由です。

10 Rで Bessel 関数の形を見る

s <- seq(0, 40, length.out = 1000)
xi <- 0.4

J0 <- besselJ(xi * s, 0)
J1 <- besselJ(xi * s, 1)

plot(s, J0, type = "l", lwd = 2,
     xlab = "s", ylab = "Bessel function",
     main = expression(J[0](xi*s) ~ "と" ~ J[1](xi*s)))
lines(s, J1, lwd = 2, lty = 2)
abline(h = 0, lty = 3)
legend("topright",
       legend = c("J0", "J1"),
       lwd = 2,
       lty = c(1, 2),
       bty = "n")

\(J_0\) は中心 \(s=0\) で最大になり、振動しながら減衰するような形をしています。

軸対称な荷重や変位の「半径方向の波」として使われます。

11 Hankel 変換の定義

0次の Hankel 変換では、関数 \(\phi(s,z)\) を

\[ \phi(s,z) = \int_0^\infty \Phi(\xi,z)J_0(\xi s)\xi d\xi \]

と表します。

逆変換は

\[ \Phi(\xi,z) = \int_0^\infty \phi(s,z)J_0(\xi s)s ds \]

です。

Fourier 変換と対応させると、

デカルト座標 軸対称円筒座標
\(x\) \(s\)
波数 \(k\) 波数 \(\xi\)
\(e^{ikx}\) \(J_0(\xi s)\)
Fourier 変換 Hankel 変換

という対応です。

12 Hankel 変換で何が簡単になるか

\[ \phi(s,z) = \int_0^\infty \Phi(\xi,z)J_0(\xi s)\xi d\xi \]

と書くと、Laplacian は

\[ \nabla^2\phi \quad\Longleftrightarrow\quad \left( \frac{\partial^2}{\partial z^2} - \xi^2 \right)\Phi \]

に変わります。

したがって、重調和方程式

\[ \nabla^2\nabla^2\phi=0 \]

は、Hankel 変換後には

\[ \left( \frac{\partial^2}{\partial z^2} - \xi^2 \right)^2 \Phi(\xi,z)=0 \]

になります。

これは \(z\) だけの常微分方程式です。

つまり、

\[ \boxed{ s,z\text{ の偏微分方程式} \quad\Longrightarrow\quad \xi\text{ ごとの }z\text{ の常微分方程式} } \]

になります。

これが Hankel 変換の最大の意味です。

13 変換後の一般解

\[ \left( \frac{d^2}{dz^2} - \xi^2 \right)^2 \Phi=0 \]

の解は、指数関数と \(z\) をかけた項の組み合わせになります。

半無限媒質を \(z<0\) とし、深部 \(z\to-\infty\) で変位が消える解だけを残すと、

\[ \Phi(\xi,z)=[C(\xi)+D(\xi)z]e^{\xi z} \]

となります。

ここで、\(C(\xi)\), \(D(\xi)\) は境界条件から決まる関数です。

この形が重要です。

14 短波長ほど浅いところで消える

\(e^{\xi z}\) の性質を図示します。

ここでは \(z<0\) とします。

z <- seq(-30, 0, length.out = 500)
xis <- c(0.05, 0.1, 0.2, 0.4)

plot(z, exp(xis[1] * z), type = "l", lwd = 2,
     xlab = "z", ylab = expression(exp(xi*z)),
     main = "短波長成分ほど深さ方向に早く減衰する",
     ylim = c(0, 1))

for (i in 2:length(xis)) {
  lines(z, exp(xis[i] * z), lwd = 2, lty = i)
}

legend("topleft",
       legend = paste0("xi = ", xis),
       lwd = 2,
       lty = 1:length(xis),
       bty = "n")

\(\xi\) が大きいほど、\(z<0\) で急速に小さくなります。

これは地球物理的には、

短波長の荷重や変形は浅部に局在し、
長波長の変形ほど深くまで影響する

という直感につながります。

15 地表応力と \(C(\xi),D(\xi)\)

地表 \(z=0\) に軸対称な応力を与えます。

\[ \sigma_{zz}(s,0)=p(s) \]

\[ \sigma_{zs}(s,0)=q(s) \]

これらを Hankel 変換して、

\[ p(s)= \int_0^\infty P(\xi)J_0(\xi s)\xi d\xi \]

\[ q(s)= \int_0^\infty Q(\xi)J_1(\xi s)\xi d\xi \]

と書きます。

このとき、境界条件から \(C(\xi)\), \(D(\xi)\) は

\[ C(\xi) = \xi^{-3} \left[ -2\nu P(\xi)+(1-2\nu)Q(\xi) \right] \]

\[ D(\xi) = \xi^{-2} \left[ P(\xi)+Q(\xi) \right] \]

と決まります。

つまり、

\[ \boxed{ \text{地表応力 }p,q \rightarrow \text{Hankel 係数 }P,Q \rightarrow C,D \rightarrow \Phi \rightarrow \phi \rightarrow u,\sigma } \]

という流れです。

16 地表変位の式

地表 \(z=0\) での変位は、

\[ u_s(s,0)= \frac{1}{2\mu} \int_0^\infty \left[ (1-2\nu)P(\xi)+2(1-\nu)Q(\xi) \right] J_1(\xi s)d\xi \]

\[ u_z(s,0)= \frac{1}{2\mu} \int_0^\infty \left[ 2(1-\nu)P(\xi)+(1-2\nu)Q(\xi) \right] J_0(\xi s)d\xi \]

となります。

この式が意味していることは明快です。

地表応力の波数成分 \(P(\xi),Q(\xi)\) を、
Bessel 関数で重ね合わせると、
地表変位が得られる。

17 簡単な数値実験1：Hankel空間で荷重を与える

ここでは、実空間の荷重 \(p(s)\) を直接与えるのではなく、Hankel 空間で

\[ P(\xi)=P_0e^{-d\xi},\quad Q(\xi)=0 \]

と与えます。

これは厳密な火山モデルではなく、Hankel 空間の成分を積分して地表変位が出ることを見るための練習です。

mu <- 30e9
nu <- 0.25
P0 <- 1e6
d <- 5e3

xi <- seq(1e-6, 0.01, length.out = 6000)
dxi <- xi[2] - xi[1]

P_xi <- P0 * exp(-d * xi)
Q_xi <- 0 * xi

surface_us <- function(s) {
  integrand <- ((1 - 2 * nu) * P_xi + 2 * (1 - nu) * Q_xi) *
    besselJ(s * xi, 1)
  sum(integrand) * dxi / (2 * mu)
}

surface_uz <- function(s) {
  integrand <- (2 * (1 - nu) * P_xi + (1 - 2 * nu) * Q_xi) *
    besselJ(s * xi, 0)
  sum(integrand) * dxi / (2 * mu)
}

s <- seq(0, 30e3, length.out = 300)

us <- sapply(s, surface_us)
uz <- sapply(s, surface_uz)

plot(s / 1000, uz, type = "l", lwd = 2,
     xlab = "s (km)", ylab = "surface displacement",
     main = "Hankel空間の鉛直荷重 P(xi) から計算した地表変位")
lines(s / 1000, us, lwd = 2, lty = 2)
abline(h = 0, lty = 3)
legend("topright",
       legend = c("uz", "us"),
       lwd = 2,
       lty = c(1, 2),
       bty = "n")

ここでは \(Q=0\) としたので、鉛直荷重だけによる変形を見ています。

18 簡単な数値実験2：\(d\) を変える

上の例で

\[ P(\xi)=P_0e^{-d\xi} \]

としました。

\(d\) が大きいほど、高波数成分、つまり短波長成分が強く抑えられます。

そのため、地表変位はなめらかで広がった形になります。

d_values <- c(2e3, 5e3, 10e3)
s <- seq(0, 40e3, length.out = 300)

uz_list <- list()

for (dd in d_values) {
  P_xi_tmp <- P0 * exp(-dd * xi)
  uz_tmp <- sapply(s, function(ss) {
    integrand <- 2 * (1 - nu) * P_xi_tmp * besselJ(ss * xi, 0)
    sum(integrand) * dxi / (2 * mu)
  })
  uz_list[[as.character(dd)]] <- uz_tmp
}

plot(s / 1000, uz_list[[1]], type = "l", lwd = 2,
     xlab = "s (km)", ylab = "uz",
     main = "P(xi)=P0 exp(-d xi)：d による広がりの違い",
     ylim = range(unlist(uz_list)))

for (i in 2:length(d_values)) {
  lines(s / 1000, uz_list[[i]], lwd = 2, lty = i)
}

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

この例では、\(d\) が大きいほど変位のピークは低くなり、空間的には広がります。

これは、深い圧力源ほど地表変位が広い範囲に分布する、という直感にもつながります。

19 簡単な数値実験3：Bessel関数の重ね合わせを見る

Hankel 変換は、関数を \(J_0(\xi s)\) の重ね合わせで表すものです。

ここでは、いくつかの \(\xi\) 成分を重ね合わせると、局在した形が作れることを見ます。

s <- seq(0, 50e3, length.out = 800)
xi_comp <- c(1/4e3, 1/8e3, 1/16e3, 1/25e3)
amp <- c(1.0, 0.8, 0.6, 0.4)

components <- sapply(seq_along(xi_comp), function(i) {
  amp[i] * besselJ(xi_comp[i] * s, 0)
})

combined <- rowSums(components)

plot(s / 1000, components[,1], type = "l", lwd = 2,
     xlab = "s (km)", ylab = "amplitude",
     main = "J0 成分の重ね合わせ",
     ylim = range(c(components, combined)))

for (i in 2:ncol(components)) {
  lines(s / 1000, components[,i], lwd = 1, lty = i)
}

lines(s / 1000, combined, lwd = 3)

legend("topright",
       legend = c(paste0("component ", 1:4), "sum"),
       lwd = c(rep(1, 4), 3),
       lty = c(1:4, 1),
       bty = "n")

実際の Hankel 変換では、有限個ではなく、連続的な \(\xi\) の成分を積分して重ね合わせます。

20 Love応力関数とHankel変換の全体像

ここまでの流れをまとめると、次のようになります。

20.1 Step 1：軸対称性を使う

\[ \frac{\partial}{\partial\varphi}=0,\quad u_\varphi=0 \]

なので、未知量は \(u_s,u_z\) になります。

20.2 Step 2：Love の応力関数を導入する

\[ u_s = -\frac{1}{2\mu} \frac{\partial^2 \phi}{\partial s \partial z} \]

\[ u_z = \frac{1}{2\mu} \left[ 2(1-\nu)\nabla^2\phi - \frac{\partial^2 \phi}{\partial z^2} \right] \]

これにより、未知量は \(\phi\) 1つになります。

20.3 Step 3：\(\phi\) は重調和方程式を満たす

\[ \nabla^2\nabla^2\phi=0 \]

20.4 Step 4：Hankel 変換する

\[ \phi(s,z) = \int_0^\infty \Phi(\xi,z)J_0(\xi s)\xi d\xi \]

すると、

\[ \nabla^2 \rightarrow \frac{d^2}{dz^2}-\xi^2 \]

なので、

\[ \left( \frac{d^2}{dz^2}-\xi^2 \right)^2 \Phi=0 \]

になります。

20.5 Step 5：半無限媒質の深部で消える解を選ぶ

\[ \Phi(\xi,z)=[C(\xi)+D(\xi)z]e^{\xi z} \]

20.6 Step 6：地表応力から \(C,D\) を決める

\[ \sigma_{zz}(s,0)=p(s),\quad \sigma_{zs}(s,0)=q(s) \]

を Hankel 変換して \(P(\xi),Q(\xi)\) を得ます。

そこから \(C(\xi),D(\xi)\) が決まり、変位も決まります。

21 最小限の理解

Love の応力関数と Hankel 変換は、別々の役割を持っています。

道具 役割
Love の応力関数 ベクトル変位問題をスカラー関数の問題にする
Hankel 変換 円筒座標の半径方向微分を波数空間で簡単にする
Bessel 関数 軸対称問題で Fourier 波の代わりになる
\(e^{\xi z}\) 半無限媒質の深さ方向減衰を表す

一言でいうと、

\[ \boxed{ \text{Love応力関数で未知量を1つにし、Hankel変換で半径方向を波数分解する} } \]

ということです。

22 付録：Rで使った主な式

このノートで数値計算に使った地表変位の式は、

\[ u_s(s,0)= \frac{1}{2\mu} \int_0^\infty \left[ (1-2\nu)P(\xi)+2(1-\nu)Q(\xi) \right] J_1(\xi s)d\xi \]

\[ u_z(s,0)= \frac{1}{2\mu} \int_0^\infty \left[ 2(1-\nu)P(\xi)+(1-2\nu)Q(\xi) \right] J_0(\xi s)d\xi \]

です。

ここで、

です。