波が海洋構造物を押す力はどう計算する？モリソン方程式をJavaScriptで実装

「海の波が構造物を押す力って、どうやって計算するんだろう？」

正直、私も最初はそう思ってました。波ってランダムで複雑で、とても1つの式で表せる気がしなかったんです。

でも、ある1つの方程式を知ったとき、「え、これだけでいいの？」と驚きました。それがモリソン方程式です。

しかも、この式のすごいところは、**慣性力（水の加速度による力）と抗力（水の速度による力）**を分けて考えられること。これだけで、波力の本質が見えてくるんです。

今日は、ブラウザで動く海洋構造物 波力シミュレーターを題材に、モリソン方程式の実装までガッツリ解説します。

---

ざっくり本質──「押す力」と「引きずる力」の2種類

まず、直感的に理解しましょう。

波が円柱に当たるとき、力は2種類あります：

1. 慣性力：水が「加速しながらぶつかる」力。水の塊がドンとぶつかるイメージ。
2. 抗力：水の流れが「こすれるように押す」力。川の流れに手を入れたときの感覚。

「波力は、水の『加速』による慣性力と、水の『速度』による抗力の和である」

どちらが支配的かは、波の周期と円柱の太さで決まります。

• 太い柱 × 短い周期 → 慣性力が支配的
• 細い柱 × 長い周期 → 抗力が支配的

この見極めに使うのが、KC数（Keulegan-Carpenter数）。ツールでは自動計算されます。

---

数式で理解する──モリソン方程式

核心の式がこちら：

F = \underbrace{\rho C_m V \frac{du}{dt}}_{\text{慣性力}} + \underbrace{\frac{1}{2}\rho C_d A \, u|u|}_{\text{抗力}}

各項の意味：

• $\rho$：海水の密度（約1025 kg/m³）
• $C_m$：慣性力係数（円柱なら約2.0）
• $V$：部材の体積（$= \frac{\pi D^2}{4} \times L$）
• $\frac{du}{dt}$：水粒子の加速度（エアリー波理論で計算）
• $C_d$：抗力係数（0.6〜1.2が一般的）
• $A$：投影面積（$= D \times L$）
• $u$：水粒子の水平速度
• $u|u|$：速度の2乗（方向を保つため絶対値付き）

エアリー波理論による速度と加速度：

水平速度：

u = \frac{\pi H}{T} \frac{\cosh(k(z+h))}{\sinh(kh)} \cos(\omega t - kx)

水平加速度：

\frac{du}{dt} = -\frac{2\pi^2 H}{T^2} \frac{\cosh(k(z+h))}{\sinh(kh)} \sin(\omega t - kx)

ここで：

• $H$：波高
• $T$：波周期
• $h$：水深
• $k$：波数（$= 2\pi / \lambda$）
• $\omega$：角周波数（$= 2\pi / T$）
• $z$：水深方向の位置（海底が$-h$、静水面が0）

波数$k$は分散関係式から求めます：

\omega^2 = gk \tanh(kh)

これ、解析的には解けないので、数値的に解く必要があります。

---

コードで実装する──JavaScriptで20行

実際にモリソン方程式を計算するコードを書いてみましょう。

// モリソン方程式による波力計算（JavaScript）
function morisonForce(H, T, h, D, L, z, t, Cm, Cd, rho=1025, g=9.81) {
    // 1. 分散関係式をニュートン法で解く（波数kを求める）
    const omega = 2 * Math.PI / T;
    let k = 0.1; // 初期推定値
    for (let i = 0; i < 100; i++) {
        const f = omega*omega - g*k*Math.tanh(k*h);
        const df = -g*Math.tanh(k*h) - g*k*h/Math.cosh(k*h)/Math.cosh(k*h);
        k -= f / df;
        if (Math.abs(f) < 1e-8) break;
    }

    // 2. エアリー波理論：水粒子速度と加速度
    const sinh_kh = Math.sinh(k * h);
    const cosh_kz_plus_h = Math.cosh(k * (z + h));
    const theta = omega * t; // 位相（x=0で計算）

    // 水平速度 u [m/s]
    const u = (Math.PI * H / T) * (cosh_kz_plus_h / sinh_kh) * Math.cos(theta);
    // 水平加速度 du/dt [m/s²]
    const du_dt = -(2 * Math.PI * Math.PI * H / (T * T)) * (cosh_kz_plus_h / sinh_kh) * Math.sin(theta);

    // 3. 幾何パラメータ
    const V = Math.PI * D * D / 4 * L; // 体積 [m³]
    const A = D * L;                   // 投影面積 [m²]

    // 4. モリソン方程式
    const F_inertia = rho * Cm * V * du_dt;       // 慣性力 [N]
    const F_drag = 0.5 * rho * Cd * A * u * Math.abs(u); // 抗力 [N]
    const F_total = F_inertia + F_drag;

    return { F_inertia, F_drag, F_total, u, du_dt, k };
}

// 使用例：波高3m、周期8秒、水深20m、直径1m、長さ10m、海底からの高さ5m
const result = morisonForce(3, 8, 20, 1, 10, -15, 0, 2.0, 0.7);
console.log(`慣性力: ${result.F_inertia.toFixed(1)} N`);
console.log(`抗力: ${result.F_drag.toFixed(1)} N`);
console.log(`合力: ${result.F_total.toFixed(1)} N`);

実行結果の例（t=0秒、波の山の位置）：

慣性力: 0.0 N
抗力: 14892.3 N
合力: 14892.3 N

波の山では加速度が0になるので慣性力がゼロ、速度最大で抗力が支配的になります。

---

数値例で確かめる──現実的なケースを計算

条件：

• 波高 $H = 5$ m
• 周期 $T = 10$ 秒
• 水深 $h = 30$ m
• 円柱径 $D = 2$ m
• 部材長 $L = 10$ m
• 水深位置 $z = -10$ m（水面下10m）
• 慣性力係数 $C_m = 2.0$
• 抗力係数 $C_d = 0.7$

計算ステップ（t=2.5秒、位相90度）：

1. 分散関係式：$\omega = 2\pi/10 = 0.628$ rad/s

   • ニュートン法で $k \approx 0.0403$ rad/m
   • 波長 $\lambda = 2\pi/k \approx 156$ m

2. 速度・加速度：

   • $\sinh(kh) = \sinh(1.209) \approx 1.521$
   • $\cosh(k(z+h)) = \cosh(0.806) \approx 1.343$
   • $u = \frac{\pi \times 5}{10} \times \frac{1.343}{1.521} \times \cos(90^\circ) = 0$ m/s
   • $\frac{du}{dt} = -\frac{2\pi^2 \times 5}{100} \times \frac{1.343}{1.521} \times \sin(90^\circ) \approx -0.871$ m/s²

3. 幾何パラメータ：

   • $V = \pi \times 2^2 / 4 \times 10 \approx 31.42$ m³
   • $A = 2 \times 10 = 20$ m²

4. モリソン方程式：

   • 慣性力：$1025 \times 2.0 \times 31.42 \times (-0.871) \approx -56,082$ N（約-5.7トン）
   • 抗力：$0.5 \times 1025 \times 0.7 \times 20 \times 0 \times 0 = 0$ N
   • 合力：-56,082 N（下向きに約5.7トン）

位相90度では速度が0で抗力ゼロ、加速度最大で慣性力が支配的になることが確認できました。

---

シミュレーターで遊ぶ──3つの実験

実験1：波高を変えてみる（慣性力 vs 抗力の変化）

• 初期条件：波高2m → 波高6m
• 何が起きる？
  • 波高が3倍になると、速度は3倍、加速度も3倍
  • 慣性力は3倍、抗力は9倍（速度の2乗なので）
  • 抗力の方が敏感に反応する

ツールで波高スライダーを動かすと、抗力の波形が急激に大きくなるのがわかります。

実験2：円柱径を変える（太さの影響）

• 直径0.5m → 2.0m（4倍）
• 何が起きる？
  • 体積は$D^2$に比例 → 16倍
  • 投影面積は$D$に比例 → 4倍
  • 慣性力が抗力より支配的になる

太い柱ほど「慣性力支配」になりやすい理由がここにあります。

実験3：潮流速度を加える（現実的な海の再現）

• 潮流速度 $U_c = 0.5$ m/s を追加
• 何が起きる？
  • 抗力の式が $\frac{1}{2}\rho C_d A (u+U_c)|u+U_c|$ に
  • 波の速度に潮流が加算され、抗力が非対称に
  • 波の進行方向と逆方向で力の大きさが変わる

ツールで「潮流速度Uc」を0.5m/sに設定すると、合力の波形が上下非対称になるのが確認できます。

---

現場でハマるポイント──適用限界と落とし穴

1. モリソン方程式は「細長いもの専用」

これはスレンダーボディ理論です。

• 対象：$D/\lambda < 0.2$ 程度（波長に対して十分細い）
• ダメな例：大型タンカー、セミサブ型プラットフォーム全体
• そういう場合はポテンシャル流理論や回折理論が必要

2. 係数は定数じゃない

$C_m$ と $C_d$ は、以下の影響を受けます：

• レイノルズ数（流れの速さと部材径で決まる）
• KC数（渦の放出パターン）
• 表面粗さ（生物付着で変化）
• 群柱効果（隣接する柱の影響）

実務ではDNV規格やISO 19902の推奨値を参照し、実験値で補正します。

3. エアリー波理論の限界

この理論は微小振幅波が前提：

• 使える目安：$H/\lambda < 0.1$（波高/波長）
• 台風時の高波（$H/\lambda > 0.15$）ではストークス波理論などが必要
• 浅海域（$h/\lambda < 0.05$）ではクノイド波理論が必要

ツールの結果を鵜呑みにせず、「入力した波は微小振幅波の仮定を満たしているか？」を常に確認しましょう。

---

まとめ──3つのポイント

1. モリソン方程式は「慣性力 + 抗力」の2項で波力を表現する
   水の加速度と速度を分けて考えることで、波力の本質が見える

2. KC数で支配的な力がわかる
   KC < 5：慣性力支配、KC > 20：抗力支配。ツールで自動計算

3. 実務では係数の設定と理論の限界に注意
   エアリー波理論は微小振幅波が前提。係数は実験値や規格を参照

---

▶ 海洋構造物 波力シミュレーター — ブラウザで即動作、登録不要

実際にスライダーを動かしながら、慣性力と抗力のバランスがどう変わるか、ぜひ体感してみてください。

NovaSolverでは1,600以上の工学シミュレーターを無料公開中 👉 一覧はこちら