噴煙柱モデルをRunge–Kutta法で解く手順(Woods, 1988)
Y. Mitsui
2026-05-13
1 目的
このノートでは、噴煙柱モデルを高度 \(z\) に関する常微分方程式として書き、4次のRunge–Kutta法で数値的に解く手順を説明する。
対象とする微分方程式は一般に
\[ \frac{d\mathbf{y}}{dz} = \mathbf{f}(z,\mathbf{y}) \]
という形である。
ここでは、状態変数として
\[ \mathbf{y} = \begin{pmatrix} M \\ J \\ H \end{pmatrix} \]
を用いる。
それぞれ、
\[ M = \beta U L^2 \]
\[ J = \beta U^2 L^2 \]
\[ H = C_p \theta \beta U L^2 \]
である。
ここで、
- \(M\)：質量フラックスに比例する量
- \(J\)：運動量フラックスに比例する量
- \(H\)：エンタルピーフラックスに比例する量
- \(U\)：噴煙柱の鉛直速度
- \(L\)：噴煙柱半径
- \(\beta\)：噴煙柱のバルク密度
- \(\theta\)：噴煙柱温度
- \(\alpha\)：周囲大気密度
- \(T\)：周囲大気温度
である。
Woods (1988) では、噴煙柱を下部の gas-thrust region と上部の buoyancy-driven convective region に分けて扱い、数値計算には4次Runge–Kutta法と5 m刻みが用いられている。Figure 2 は、\(U_0=300\) m/s、\(\theta_0=1000\) K、\(n_0=0.03\) として、火口半径 \(L_0=20,50,100,150,200\) m を変えたときの速度・温度・半径・密度プロファイルを示している。
このR Markdownでは、Runge–Kutta法の具体的な手順を確認したうえで、Figure 2に近い構成の4パネル図を作成する。
2 Runge–Kutta法でこのモデルを解くとは何をすることか
このモデルでは、噴煙柱の性質が高度 \(z\) とともにどのように変化するかを求めたい。
つまり、時間発展問題ではなく、高度方向の発展問題として、
\[ \frac{d\mathbf{y}}{dz} = \mathbf{f}(z,\mathbf{y}) \]
を解く。
ここで、\(\mathbf{y}\) は噴煙柱の状態を表すベクトルであり、本ノートでは
\[ \mathbf{y} = (M,J,H) \]
とする。
Runge–Kutta法が直接計算するのは、速度 \(U\)、半径 \(L\)、密度 \(\beta\)、温度 \(\theta\) そのものではない。Runge–Kutta法が更新するのは、あくまで状態変数
\[ M,\quad J,\quad H \]
である。
ただし、微分方程式の右辺 \(\mathbf{f}(z,\mathbf{y})\) を計算するには、その時点での
\[ U,\quad L,\quad \beta,\quad \theta \]
が必要になる。
したがって、1ステップの計算では次の操作を行う。
- 現在の高度 \(z\) と状態変数 \((M,J,H)\) を用意する。
- \((M,J,H)\) から \(U,L,\beta,\theta\) を復元する。
- 復元した物理量を使って、\((dM/dz,dJ/dz,dH/dz)\) を計算する。
- Runge–Kutta法により、次の高度 \(z+\Delta z\) における \((M,J,H)\) を求める。
- 新しい \((M,J,H)\) から、再び \(U,L,\beta,\theta\) を復元する。
この繰り返しにより、火口から上方へ向かって、噴煙柱の速度、半径、温度、密度の鉛直プロファイルを得る。
重要なのは、Runge–Kutta法自体は噴煙柱の物理を知っているわけではない、という点である。Runge–Kutta法は、与えられた右辺
\[ \mathbf{f}(z,\mathbf{y}) \]
を何度も評価し、その傾きを組み合わせて次の状態を高精度に推定する数値的な手続きである。噴煙柱の物理は、右辺 \(\mathbf{f}(z,\mathbf{y})\) の中に入っている。
3 4次Runge–Kutta法の基本
Euler法であれば、現在地の傾きだけを使って、
\[ \mathbf{y}_{i+1} = \mathbf{y}_i + \Delta z\mathbf{f}(z_i,\mathbf{y}_i) \]
と進める。
しかし、これでは1ステップの間で傾きが変化する効果を十分に考慮できない。
4次Runge–Kutta法では、1ステップの中で右辺を4回評価する。
まず、現在地での傾きを
\[ \mathbf{k}_1 = \mathbf{f}(z_i,\mathbf{y}_i) \]
として求める。
次に、\(\mathbf{k}_1\) を使って半ステップ進んだと仮定し、その場所での傾きを
\[ \mathbf{k}_2 = \mathbf{f} \left( z_i+\frac{\Delta z}{2}, \mathbf{y}_i+\frac{\Delta z}{2}\mathbf{k}_1 \right) \]
として求める。
さらに、\(\mathbf{k}_2\) を使ってもう一度半ステップ先の傾きを
\[ \mathbf{k}_3 = \mathbf{f} \left( z_i+\frac{\Delta z}{2}, \mathbf{y}_i+\frac{\Delta z}{2}\mathbf{k}_2 \right) \]
として求める。
最後に、\(\mathbf{k}_3\) を使って1ステップ先まで進んだと仮定し、
\[ \mathbf{k}_4 = \mathbf{f} \left( z_i+\Delta z, \mathbf{y}_i+\Delta z\mathbf{k}_3 \right) \]
を求める。
そして、これら4つの傾きを
\[ 1:2:2:1 \]
の重みで平均し、
\[ \mathbf{y}_{i+1} = \mathbf{y}_i + \frac{\Delta z}{6} \left( \mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4 \right) \]
として次の高度の状態を求める。
本モデルに即して言えば、Runge–Kutta法の1ステップは次のように解釈できる。
現在の高度で、まず噴煙柱の速度、半径、温度、密度を求める。次に、その状態でどれだけ空気を取り込むか、運動量がどれだけ変化するか、エンタルピーがどれだけ変化するかを計算する。ただし、その傾きをそのまま使うのではなく、半ステップ先、さらに半ステップ先、1ステップ先でも同様の評価を行う。最後に、それらを重み付き平均して、次の高度の状態を決める。
4 モデル式の概要
4.1 状態変数
状態変数は
\[ \mathbf{y} = (M,J,H) \]
とする。
速度は
\[ U = \frac{J}{M} \]
で復元できる。
半径は
\[ L = \sqrt{\frac{M}{\beta U}} \]
で復元できる。
温度は
\[ \theta = \frac{H}{C_p M} \]
で復元できる。
4.2 gas-thrust region
gas-thrust region では、Woods (1988) の式に対応する形として、速度変化を
\[ U \frac{dU}{dz} = -\frac{U^2}{8L} \sqrt{\frac{\alpha}{\beta}} + \frac{g(\alpha-\beta)}{\beta} \]
とする。
したがって、
\[ \frac{dU}{dz} = -\frac{U}{8L} \sqrt{\frac{\alpha}{\beta}} + \frac{g(\alpha-\beta)}{\beta U} \]
である。
また、運動量フラックスは
\[ \frac{dJ}{dz} = g(\alpha-\beta)L^2 \]
である。
ここで、
\[ J = MU \]
なので、
\[ \frac{dJ}{dz} = U\frac{dM}{dz} + M\frac{dU}{dz} \]
である。したがって、
\[ \frac{dM}{dz} = \frac{1}{U} \left( \frac{dJ}{dz} - M\frac{dU}{dz} \right) \]
として、gas-thrust region における質量フラックスの増加率を求める。
この形にすると、簡略的には
\[ \frac{dM}{dz} \sim \frac{UL}{8}\sqrt{\alpha\beta} \]
に対応する。 \(L\) が抜けると、エントレインメント量が大きく異なり、Figure 2 の形から大きく外れる。
4.3 convective region
噴煙柱が周囲大気より軽くなり、
\[ \beta < \alpha \]
となったら、convective region に入ったとみなす。
convective region では、
\[ \frac{dM}{dz} = 2kUL\alpha \]
とする。
運動量フラックスは gas-thrust region と同様に、
\[ \frac{dJ}{dz} = g(\alpha-\beta)L^2 \]
である。
4.4 エンタルピーフラックス
両領域で、エンタルピーフラックスは
\[ \frac{dH}{dz} = C_aT\frac{dM}{dz} + \frac{U^2}{2}\frac{dM}{dz} - \alpha U L^2 g \]
とする。
右辺第1項は取り込まれる大気のエンタルピー、第2項は運動エネルギーの寄与、第3項は位置エネルギーの増加に対応する。
5 数値計算の実装
5.1 計算の準備
rm(list = ls())
# -----------------------------
# 物理定数
# -----------------------------
g <- 9.81 # 重力加速度 [m/s^2]
Ra <- 285 # 乾燥空気の気体定数 [J/(kg K)]
Rm <- 462 # 火山ガスの気体定数 [J/(kg K)]
Ca <- 998 # 空気の定圧比熱 [J/(kg K)]
Cm <- 1617 # 火山ガスの定圧比熱 [J/(kg K)]
Cs <- 1617 # 固体火山砕屑物の比熱 [J/(kg K)]
rho_s <- 2500 # 固体粒子密度 [kg/m^3]
k_ent <- 0.09 # convective region のエントレインメント係数
# -----------------------------
# 大気モデルのパラメータ
# -----------------------------
T0 <- 273 # 地表大気温度 [K]
P0 <- 101325 # 地表気圧 [Pa]
mu <- 6.5e-3 # 対流圏の気温減率 [K/m]
omega <- 2.0e-3 # 成層圏での気温増加率 [K/m]
H1 <- 11000 # 対流圏界面高度 [m]
H2 <- 20000 # 成層圏下部の代表高度 [m]5.2 大気温度・気圧・密度
ambient_temperature <- function(z) {
if (z < H1) {
T <- T0 - mu * z
} else if (z < H2) {
T <- T0 - mu * H1
} else {
T <- T0 - mu * H1 + omega * (z - H2)
}
return(T)
}
ambient_pressure <- function(z) {
if (z < H1) {
Tz <- T0 - mu * z
P <- P0 * (Tz / T0)^(g / (Ra * mu))
} else if (z < H2) {
T1 <- T0 - mu * H1
P1 <- P0 * (T1 / T0)^(g / (Ra * mu))
P <- P1 * exp(-g * (z - H1) / (Ra * T1))
} else {
T1 <- T0 - mu * H1
P1 <- P0 * (T1 / T0)^(g / (Ra * mu))
P2 <- P1 * exp(-g * (H2 - H1) / (Ra * T1))
Tz <- T1 + omega * (z - H2)
P <- P2 * (Tz / T1)^(-g / (Ra * omega))
}
return(P)
}
ambient_density <- function(z) {
T <- ambient_temperature(z)
P <- ambient_pressure(z)
alpha <- P / (Ra * T)
return(alpha)
}5.3 1ケースを計算する関数
run_column_case <- function(
U0 = 300,
L0 = 100,
theta0 = 1000,
n0 = 0.03,
dz = 5,
zmax = 50000,
u_stop = 1e-3
) {
# -----------------------------
# 初期条件
# -----------------------------
P_vent <- ambient_pressure(0)
beta0 <- 1 / ((1 - n0) / rho_s + n0 * Rm * theta0 / P_vent)
M0_local <- beta0 * U0 * L0^2
J0_local <- beta0 * U0^2 * L0^2
Cp0 <- n0 * Cm + (1 - n0) * Cs
H0_flux <- Cp0 * theta0 * M0_local
y0_local <- c(M = M0_local, J = J0_local, H = H0_flux)
# -----------------------------
# 状態変数から物理量を復元
# -----------------------------
recover_state_local <- function(z, y) {
M <- as.numeric(y["M"])
J <- as.numeric(y["J"])
H <- as.numeric(y["H"])
if (!is.finite(M) || !is.finite(J) || !is.finite(H)) {
return(NULL)
}
if (M <= 0) {
return(NULL)
}
U <- J / M
if (!is.finite(U) || U <= 0) {
return(NULL)
}
# 数値誤差でMがごくわずかにM0_localを下回る場合は補正
if (M < M0_local && abs(M - M0_local) < 1e-8 * M0_local) {
M <- M0_local
}
if (M < M0_local) {
return(NULL)
}
# 初期に存在する固体と火山ガスは保存されると考える
mass_solid <- (1 - n0) * M0_local
mass_volc_gas <- n0 * M0_local
# 増加分を取り込まれた大気とみなす
mass_air <- M - M0_local
if (mass_air < 0 && abs(mass_air) < 1e-8 * M0_local) {
mass_air <- 0
}
if (mass_air < 0) {
return(NULL)
}
mass_gas <- mass_air + mass_volc_gas
if (mass_gas <= 0) {
return(NULL)
}
n <- mass_gas / M
Rg <- (mass_air * Ra + mass_volc_gas * Rm) / mass_gas
Cp <- (mass_air * Ca + mass_volc_gas * Cm + mass_solid * Cs) / M
theta <- H / (Cp * M)
if (!is.finite(theta) || theta <= 0) {
return(NULL)
}
P <- ambient_pressure(z)
beta <- 1 / ((1 - n) / rho_s + n * Rg * theta / P)
if (!is.finite(beta) || beta <= 0) {
return(NULL)
}
L <- sqrt(M / (beta * U))
if (!is.finite(L) || L <= 0) {
return(NULL)
}
alpha <- ambient_density(z)
T <- ambient_temperature(z)
if (!is.finite(alpha) || alpha <= 0 || !is.finite(T) || T <= 0) {
return(NULL)
}
return(list(
z = z,
M = M,
J = J,
H = H,
U = U,
L = L,
beta = beta,
theta = theta,
n = n,
Rg = Rg,
Cp = Cp,
alpha = alpha,
T = T,
P = P
))
}
# -----------------------------
# 微分方程式の右辺
# -----------------------------
rhs_local <- function(z, y, region = "gas_thrust") {
s <- recover_state_local(z, y)
if (is.null(s)) {
return(c(M = NA, J = NA, H = NA))
}
M <- s$M
U <- s$U
L <- s$L
beta <- s$beta
alpha <- s$alpha
T <- s$T
if (region == "gas_thrust") {
# Woods (1988) model A に対応する速度方程式
dUdz <- - (U / (8 * L)) * sqrt(alpha / beta) +
g * (alpha - beta) / (beta * U)
# 運動量保存
dJdz <- g * (alpha - beta) * L^2
# J = M U より dM/dz を求める
dMdz <- (dJdz - M * dUdz) / U
} else if (region == "convective") {
# Morton型のエントレインメント
dMdz <- 2 * k_ent * U * L * alpha
# 運動量保存
dJdz <- g * (alpha - beta) * L^2
} else {
stop("region must be 'gas_thrust' or 'convective'")
}
# エンタルピーフラックス
dHdz <- Ca * T * dMdz + 0.5 * U^2 * dMdz - alpha * U * L^2 * g
return(c(M = dMdz, J = dJdz, H = dHdz))
}
# -----------------------------
# 4次Runge--Kutta法の1ステップ
# -----------------------------
rk4_local <- function(z, y, dz, region) {
k1 <- rhs_local(z, y, region)
if (any(!is.finite(k1))) return(rep(NA, length(y)))
k2 <- rhs_local(z + dz / 2, y + dz * k1 / 2, region)
if (any(!is.finite(k2))) return(rep(NA, length(y)))
k3 <- rhs_local(z + dz / 2, y + dz * k2 / 2, region)
if (any(!is.finite(k3))) return(rep(NA, length(y)))
k4 <- rhs_local(z + dz, y + dz * k3, region)
if (any(!is.finite(k4))) return(rep(NA, length(y)))
y_next <- y + dz * (k1 + 2 * k2 + 2 * k3 + k4) / 6
names(y_next) <- names(y)
return(y_next)
}
# -----------------------------
# 1ステップ実演用の値を保存
# -----------------------------
z_demo <- 0
y_demo <- y0_local
region_demo <- "gas_thrust"
s_demo0 <- recover_state_local(z_demo, y_demo)
k1_demo <- rhs_local(z_demo, y_demo, region_demo)
k2_demo <- rhs_local(z_demo + dz / 2, y_demo + dz * k1_demo / 2, region_demo)
k3_demo <- rhs_local(z_demo + dz / 2, y_demo + dz * k2_demo / 2, region_demo)
k4_demo <- rhs_local(z_demo + dz, y_demo + dz * k3_demo, region_demo)
y_demo_next <- y_demo + dz * (k1_demo + 2 * k2_demo + 2 * k3_demo + k4_demo) / 6
names(y_demo_next) <- names(y_demo)
s_demo1 <- recover_state_local(z_demo + dz, y_demo_next)
# -----------------------------
# 高度方向に積分
# -----------------------------
z <- 0
y <- y0_local
region <- "gas_thrust"
gas_thrust_height <- NA_real_
output <- data.frame()
for (i in 1:ceiling(zmax / dz)) {
s <- recover_state_local(z, y)
if (is.null(s)) {
break
}
output <- rbind(output, data.frame(
step = i,
z = z,
region = region,
L0 = L0,
U0 = U0,
theta0 = theta0,
n0 = n0,
M = s$M,
J = s$J,
H = s$H,
U = s$U,
L = s$L,
beta = s$beta,
alpha = s$alpha,
theta = s$theta,
T = s$T,
n = s$n,
gas_thrust_height = gas_thrust_height
))
# gas-thrust region から convective region への切り替え
if (region == "gas_thrust" && s$beta < s$alpha) {
region <- "convective"
gas_thrust_height <- z
}
# Figure 2 に近づけるため、中立浮力高度では止めない
# 速度がほぼゼロになったら最大高度付近とみなして止める
if (s$U <= u_stop) {
break
}
y_next <- rk4_local(z, y, dz, region)
if (any(!is.finite(y_next))) {
break
}
y <- y_next
z <- z + dz
}
attr(output, "demo") <- list(
s0 = s_demo0,
k1 = k1_demo,
k2 = k2_demo,
k3 = k3_demo,
k4 = k4_demo,
y0 = y_demo,
y_next = y_demo_next,
s1 = s_demo1
)
return(output)
}5.4 4次Runge–Kutta法の1ステップ実演
代表例として、
\[ U_0 = 300 \ {\rm m/s},\quad L_0 = 100 \ {\rm m},\quad \theta_0 = 1000 \ {\rm K},\quad n_0 = 0.03 \]
を用いる。
base_result <- run_column_case(
U0 = 300,
L0 = 100,
theta0 = 1000,
n0 = 0.03,
dz = 5,
zmax = 50000
)
demo <- attr(base_result, "demo")
cat("Initial state vector y0:\n")## Initial state vector y0:print(demo$y0)## M J H
## 2.186978e+07 6.560935e+09 3.536344e+13cat("\nInitial physical variables at z = 0 m:\n")##
## Initial physical variables at z = 0 m:cat("U =", demo$s0$U, "\n")## U = 300cat("L =", demo$s0$L, "\n")## L = 100cat("beta =", demo$s0$beta, "\n")## beta = 7.289928cat("theta =", demo$s0$theta, "\n")## theta = 1000cat("\nk1:\n")##
## k1:print(demo$k1)## M J H
## 11554.4 -587386.9 3629664348.6cat("\nk2:\n")##
## k2:print(demo$k2)## M J H
## 11551.63 -585840.54 3625976730.65cat("\nk3:\n")##
## k3:print(demo$k3)## M J H
## 11551.63 -585840.63 3625979000.62cat("\nk4:\n")##
## k4:print(demo$k4)## M J H
## 11548.87 -584290.55 3622302278.10cat("\nState vector at z = 5 m:\n")##
## State vector at z = 5 m:print(demo$y_next)## M J H
## 2.192754e+07 6.558006e+09 3.538157e+13cat("\nPhysical variables at z = 5 m:\n")##
## Physical variables at z = 5 m:cat("U =", demo$s1$U, "\n")## U = 299.0762cat("L =", demo$s1$L, "\n")## L = 102.8055cat("beta =", demo$s1$beta, "\n")## beta = 6.937061cat("theta =", demo$s1$theta, "\n")## theta = 998.8845この出力は、Runge–Kutta法の1ステップで、\(k_1,k_2,k_3,k_4\) がどのような値になり、それらから次の高度の状態変数がどのように得られるかを確認するためのものである。
5.5 代表ケースの高度方向積分
head(base_result)## step z region L0 U0 theta0 n0 M J H
## 1 1 0 gas_thrust 100 300 1000 0.03 21869784 6560935262 3.536344e+13
## 2 2 5 gas_thrust 100 300 1000 0.03 21927542 6558006062 3.538157e+13
## 3 3 10 gas_thrust 100 300 1000 0.03 21985273 6555092382 3.539966e+13
## 4 4 15 gas_thrust 100 300 1000 0.03 22042976 6552194294 3.541772e+13
## 5 5 20 gas_thrust 100 300 1000 0.03 22100651 6549311874 3.543574e+13
## 6 6 25 gas_thrust 100 300 1000 0.03 22158300 6546445195 3.545373e+13
## U L beta alpha theta T n
## 1 300.0000 100.0000 7.289928 1.302294 1000.0000 273.0000 0.03000000
## 2 299.0762 102.8055 6.937061 1.301628 998.8845 272.9675 0.03255503
## 3 298.1583 105.5485 6.618828 1.300963 997.7723 272.9350 0.03509541
## 4 297.2464 108.2337 6.330369 1.300297 996.6635 272.9025 0.03762129
## 5 296.3402 110.8653 6.067692 1.299632 995.5581 272.8700 0.04013279
## 6 295.4399 113.4470 5.827489 1.298968 994.4559 272.8375 0.04263002
## gas_thrust_height
## 1 NA
## 2 NA
## 3 NA
## 4 NA
## 5 NA
## 6 NAtail(base_result)## step z region L0 U0 theta0 n0 M J
## 4378 4378 21885 convective 100 300 1000 0.03 435550009 5079163858
## 4379 4379 21890 convective 100 300 1000 0.03 435562703 4613556825
## 4380 4380 21895 convective 100 300 1000 0.03 435574731 4094869843
## 4381 4381 21900 convective 100 300 1000 0.03 435585960 3499578786
## 4382 4382 21905 convective 100 300 1000 0.03 435596174 2778864319
## 4383 4383 21910 convective 100 300 1000 0.03 435604917 1786568517
## H U L beta alpha theta T
## 4378 7.320009e+13 11.661494 22851.76 0.07152283 0.05418581 163.3142 205.27
## 4379 7.318650e+13 10.592176 23985.37 0.07147806 0.05413776 163.2793 205.28
## 4380 7.317279e+13 9.401073 25467.84 0.07143336 0.05408976 163.2443 205.29
## 4381 7.315891e+13 8.034186 27558.22 0.07138873 0.05404180 163.2093 205.30
## 4382 7.314483e+13 6.379451 30936.50 0.07134418 0.05399389 163.1742 205.31
## 4383 7.313046e+13 4.101351 38595.71 0.07129973 0.05394602 163.1389 205.32
## n gas_thrust_height
## 4378 0.9512945 740
## 4379 0.9512959 740
## 4380 0.9512972 740
## 4381 0.9512985 740
## 4382 0.9512996 740
## 4383 0.9513006 740if (nrow(base_result) == 0) {
stop("base_result が空です。初期条件またはモデル式を確認してください。")
}
cat("最大到達高度 [m]:", max(base_result$z), "\n")## 最大到達高度 [m]: 21910cat("最大到達高度 [km]:", max(base_result$z) / 1000, "\n\n")## 最大到達高度 [km]: 21.91summary(base_result[, c("z", "U", "L", "beta", "alpha", "theta", "T")])## z U L beta
## Min. : 0 Min. : 4.101 Min. : 100 Min. :0.0713
## 1st Qu.: 5478 1st Qu.:134.916 1st Qu.: 1369 1st Qu.:0.1450
## Median :10955 Median :144.339 Median : 2480 Median :0.2878
## Mean :10955 Mean :138.538 Mean : 3312 Mean :0.3954
## 3rd Qu.:16433 3rd Qu.:151.761 3rd Qu.: 4411 3rd Qu.:0.4963
## Max. :21910 Max. :300.000 Max. :38596 Max. :7.2899
## alpha theta T
## Min. :0.05395 Min. : 163.1 Min. :201.5
## 1st Qu.:0.13969 1st Qu.: 205.3 1st Qu.:201.5
## Median :0.35554 Median : 269.9 Median :204.5
## Mean :0.45678 Mean : 341.0 Mean :219.6
## 3rd Qu.:0.71453 3rd Qu.: 405.2 3rd Qu.:237.4
## Max. :1.30229 Max. :1000.0 Max. :273.05.6 作図用の安全な関数
safe_plot_profile <- function(x, z, xlab, main) {
ok <- is.finite(x) & is.finite(z)
if (sum(ok) < 2) {
message(paste0(main, ": 有限な値が不足しているため作図できません。"))
return(invisible(NULL))
}
plot(
x[ok],
z[ok] / 1000,
type = "l",
xlab = xlab,
ylab = "Height [km]",
main = main
)
}5.7 代表ケースのプロファイル
op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
safe_plot_profile(
base_result$U,
base_result$z,
xlab = "Velocity U [m/s]",
main = "Velocity profile"
)
abline(v = 0, lty = 2)
safe_plot_profile(
base_result$L,
base_result$z,
xlab = "Column radius L [m]",
main = "Radius profile"
)
ok <- is.finite(base_result$theta) & is.finite(base_result$T) & is.finite(base_result$z)
if (sum(ok) >= 2) {
plot(
base_result$theta[ok],
base_result$z[ok] / 1000,
type = "l",
xlab = "Temperature [K]",
ylab = "Height [km]",
main = "Temperature profile"
)
lines(base_result$T[ok], base_result$z[ok] / 1000, lty = 2)
legend("topright", legend = c("Column", "Ambient"), lty = c(1, 2), bty = "n")
}
ok <- is.finite(base_result$beta) & is.finite(base_result$alpha) & is.finite(base_result$z)
if (sum(ok) >= 2) {
plot(
base_result$beta[ok],
base_result$z[ok] / 1000,
type = "l",
xlab = "Density [kg/m^3]",
ylab = "Height [km]",
main = "Density profile"
)
lines(base_result$alpha[ok], base_result$z[ok] / 1000, lty = 2)
legend("topright", legend = c("Column", "Ambient"), lty = c(1, 2), bty = "n")
}
par(op)6 Figure 2 型プロット
6.1 Figure 2 型プロットの計算
Woods (1988) の Figure 2 と同様に、
\[ U_0 = 300 \ {\rm m/s} \]
\[ \theta_0 = 1000 \ {\rm K} \]
\[ n_0 = 0.03 \]
として、火口半径
\[ L_0 = 20,\ 50,\ 100,\ 150,\ 200 \ {\rm m} \]
を変える。
L0_values <- c(20, 50, 100, 150, 200)
fig2_result <- do.call(
rbind,
lapply(L0_values, function(LL) {
run_column_case(
U0 = 300,
L0 = LL,
theta0 = 1000,
n0 = 0.03,
dz = 5,
zmax = 50000
)
})
)
if (is.null(fig2_result) || nrow(fig2_result) == 0) {
stop("fig2_result が空です。Figure 2型プロットを作成できません。")
}
head(fig2_result)## step z region L0 U0 theta0 n0 M J H
## 1 1 0 gas_thrust 20 300 1000 0.03 874791.4 262437410 1.414538e+12
## 2 2 5 gas_thrust 20 300 1000 0.03 886343.0 262321587 1.418189e+12
## 3 3 10 gas_thrust 20 300 1000 0.03 897889.2 262209109 1.421822e+12
## 4 4 15 gas_thrust 20 300 1000 0.03 909430.0 262100033 1.425439e+12
## 5 5 20 gas_thrust 20 300 1000 0.03 920965.4 261994414 1.429039e+12
## 6 6 25 gas_thrust 20 300 1000 0.03 932495.7 261892307 1.432622e+12
## U L beta alpha theta T n
## 1 300.0000 20.00000 7.289928 1.302294 1000.0000 273.0000 0.03000000
## 2 295.9594 22.64349 5.840940 1.301628 994.4762 272.9675 0.04264194
## 3 292.0284 25.04487 4.901852 1.300963 989.0324 272.9350 0.05495284
## 4 288.2025 27.26778 4.243967 1.300297 983.6670 272.9025 0.06694561
## 5 284.4780 29.35244 3.757563 1.299632 978.3785 272.8700 0.07863250
## 6 280.8510 31.32615 3.383428 1.298968 973.1652 272.8375 0.09002514
## gas_thrust_height
## 1 NA
## 2 NA
## 3 NA
## 4 NA
## 5 NA
## 6 NAtail(fig2_result)## step z region L0 U0 theta0 n0 M J
## 21409 6011 30050 convective 200 300 1000 0.03 897212103 12811046231
## 21410 6012 30055 convective 200 300 1000 0.03 897221395 11674107631
## 21411 6013 30060 convective 200 300 1000 0.03 897230216 10413179408
## 21412 6014 30065 convective 200 300 1000 0.03 897238478 8976156226
## 21413 6015 30070 convective 200 300 1000 0.03 897246041 7259171492
## 21414 6016 30075 convective 200 300 1000 0.03 897252623 4980203076
## H U L beta alpha theta T
## 21409 1.501943e+14 14.278726 55017.17 0.02075913 0.01344280 158.1713 221.60
## 21410 1.501679e+14 13.011401 57651.65 0.02074684 0.01343176 158.1419 221.61
## 21411 1.501413e+14 11.605917 61061.11 0.02073457 0.01342073 158.1125 221.62
## 21412 1.501147e+14 10.004203 65787.52 0.02072232 0.01340971 158.0831 221.63
## 21413 1.500879e+14 8.090503 73177.43 0.02071008 0.01339869 158.0536 221.64
## 21414 1.500609e+14 5.550503 88374.82 0.02069785 0.01338769 158.0241 221.65
## n gas_thrust_height
## 21409 0.9054240 1485
## 21410 0.9054249 1485
## 21411 0.9054259 1485
## 21412 0.9054267 1485
## 21413 0.9054275 1485
## 21414 0.9054282 1485cat("各ケースの最大高度:\n")## 各ケースの最大高度:print(aggregate(z ~ L0, data = fig2_result, FUN = max))## L0 z
## 1 20 12250
## 2 50 16585
## 3 100 21910
## 4 150 26225
## 5 200 300756.2 Figure 2 型の4パネル図
op <- par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
line_types <- c(1, 2, 3, 4, 5)
# -----------------------------
# (a) Velocity
# -----------------------------
plot(
NA,
xlim = range(fig2_result$U, finite = TRUE),
ylim = range(fig2_result$z / 1000, finite = TRUE),
xlab = "Velocity [m/s]",
ylab = "Height [km]",
main = "(a) Velocity"
)
for (i in seq_along(L0_values)) {
LL <- L0_values[i]
tmp <- fig2_result[fig2_result$L0 == LL, ]
ok <- is.finite(tmp$U) & is.finite(tmp$z)
lines(tmp$U[ok], tmp$z[ok] / 1000, lty = line_types[i])
}
legend(
"topright",
legend = paste0("L0 = ", L0_values, " m"),
lty = line_types,
bty = "n",
cex = 0.8
)
# -----------------------------
# (b) Temperature
# -----------------------------
plot(
NA,
xlim = range(fig2_result$theta, finite = TRUE),
ylim = range(fig2_result$z / 1000, finite = TRUE),
xlab = "Temperature [K]",
ylab = "Height [km]",
main = "(b) Temperature"
)
for (i in seq_along(L0_values)) {
LL <- L0_values[i]
tmp <- fig2_result[fig2_result$L0 == LL, ]
ok <- is.finite(tmp$theta) & is.finite(tmp$z)
lines(tmp$theta[ok], tmp$z[ok] / 1000, lty = line_types[i])
}
# -----------------------------
# (c) Radius
# -----------------------------
plot(
NA,
xlim = range(fig2_result$L, finite = TRUE),
ylim = range(fig2_result$z / 1000, finite = TRUE),
xlab = "Radius [m]",
ylab = "Height [km]",
main = "(c) Radius"
)
for (i in seq_along(L0_values)) {
LL <- L0_values[i]
tmp <- fig2_result[fig2_result$L0 == LL, ]
ok <- is.finite(tmp$L) & is.finite(tmp$z)
lines(tmp$L[ok], tmp$z[ok] / 1000, lty = line_types[i])
}
# -----------------------------
# (d) Density
# -----------------------------
plot(
NA,
xlim = range(c(fig2_result$beta, fig2_result$alpha), finite = TRUE),
ylim = range(fig2_result$z / 1000, finite = TRUE),
xlab = "Density [kg/m^3]",
ylab = "Height [km]",
main = "(d) Density"
)
for (i in seq_along(L0_values)) {
LL <- L0_values[i]
tmp <- fig2_result[fig2_result$L0 == LL, ]
ok <- is.finite(tmp$beta) & is.finite(tmp$z)
lines(tmp$beta[ok], tmp$z[ok] / 1000, lty = line_types[i])
}
# 大気密度を太い破線で追加
tmp_ref <- fig2_result[fig2_result$L0 == L0_values[1], ]
ok <- is.finite(tmp_ref$alpha) & is.finite(tmp_ref$z)
lines(
tmp_ref$alpha[ok],
tmp_ref$z[ok] / 1000,
lty = 2,
lwd = 2
)
legend(
"topright",
legend = c(paste0("L0 = ", L0_values, " m"), "Ambient"),
lty = c(line_types, 2),
lwd = c(rep(1, length(L0_values)), 2),
bty = "n",
cex = 0.75
)
par(op)7 Figure 2との違いについて
下記は重要な点である。
- gas-thrust region で \(dM/dz\) を直接簡略式で置かず、速度方程式と運動量保存式から整合的に求めるようにした。
- 中立浮力高度、すなわち \(\beta\) が再び \(\alpha\) に近づくところで計算を止めず、速度がほぼ0になるまで積分するようにした。
- Figure 2と同じく、火口半径 \(L_0=20,50,100,150,200\) m の5ケースを比較するようにした。
なお、このR MarkdownはRunge–Kutta法の理解とFigure 2型の再現を目的にした整理版であり、元論文の全てのモデル実装上の細部、たとえば傘型雲領域の扱い、半径定義の係数、初期圧力条件、上部の停止条件、粒子脱落の無視の仕方、比熱と気体定数の更新の細部などを完全には再現していない。
ただし、Figure 2の主要な特徴、すなわち
- 下部 gas-thrust region で半径が急激に増加すること
- 火口半径が大きいほど高く上がりやすいこと
- 大きな火口半径では速度が途中で再増加しうること
- 温度が高度とともに低下すること
- 密度が周囲大気密度より小さくなることで convective region に移行すること
を確認するためのコードとして使える。
8 まとめ
Runge–Kutta法で行っていることは、単に「式を解く」ことではなく、各高度での物理状態を復元しながら、微分方程式の右辺を何度も評価し、その平均的な傾きによって噴煙柱の状態を少しずつ上方へ進めることである。
この意味で、Runge–Kutta法は、噴煙柱モデルの物理過程そのものを与える方法ではない。物理過程は、質量保存、運動量保存、エンタルピー収支、エントレインメントの式としてモデルに入っている。Runge–Kutta法は、それらの式から得られる高度方向の変化率を使って、数値的にプロファイルを構成するための方法である。
Woods (1988) のような噴煙柱モデルでは、解析解を直接書くことは難しい。そのため、火口での初期条件を与え、そこから5 m程度の小さな刻みで状態を更新していく。このようにして、速度、温度、半径、密度が高度とともにどのように変化するかを計算する。
今回のR Markdownでは、まずRunge–Kutta法が何を更新しているのかを整理し、その後にモデル式を示し、最後にRコードで実際に高度方向へ積分した。したがって、計算結果の図を見るときには、各線が単に描かれているのではなく、各高度で\((M,J,H)\)を更新し、そこから\(U,L,\beta,\theta\)を復元する作業の繰り返しによって得られている、と理解すればよい。