テラメカニクス入門 — 不整地走行の力学を理解する

月面ローバーが柔らかいレゴリスの上を走ろうとしたとき、車輪はどこまで沈むのでしょうか。火星探査ローバー Spirit は2009年、砂地にはまり込んで動けなくなり、二度と脱出できませんでした。舗装路を走る自動車では考えられない事態ですが、地球外の天体では「地面そのもの」が最大の敵になり得るのです。

砂浜を歩いたことがある人なら、足が砂にめり込んで歩きにくい経験をしたことがあるでしょう。靴底が砂に沈むほど、足を引き抜くのに余分な力が必要になります。これと全く同じことが、月面や火星の地表を走る車輪にも起こります。では、車輪がどれだけ沈むのか、どれだけの力で前に進めるのかを、定量的に予測するにはどうすればよいでしょうか。その答えを与えるのがテラメカニクス（Terramechanics）です。

テラメカニクスは「土壌と車両の相互作用の力学」であり、以下のような幅広い分野で不可欠な知識です。

• 惑星探査ローバー: 月面・火星表面のレゴリス上での走行性能を予測し、車輪の設計やミッション計画に活用します。Spirit のスタック事故のような事態を事前にシミュレーションで回避する技術の基盤です
• 軍用車両: 泥濘地や砂漠での走破性能の評価に、Bekker 理論が1950年代から使われてきました。テラメカニクスはそもそも軍事技術として誕生しました
• 農業機械: トラクターが畑を走る際の沈下や走行抵抗の予測は、エネルギー効率の最適化に直結します
• 災害対応ロボット: 瓦礫や軟弱地盤の上で活動するレスキューロボットの移動能力の評価に用いられます

本記事の内容

• テラメカニクスの基本概念 — 土壌と車輪の相互作用とは何か
• Bekker の圧力-沈下モデル — 車輪がどれだけ沈むかを予測する
• Janosi-Hanamoto のせん断応力モデル — 車輪が地面を掴む力を予測する
• 沈下と走行抵抗の関係 — なぜ沈むと走りにくくなるのか
• スリップ率と牽引力 — 車輪の空転と推進力の定量的関係
• 月面・火星レゴリスの土壌パラメータ
• Spirit 砂地スタック事例の力学的分析
• Python による車輪-土壌相互作用のシミュレーション

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 惑星探査ローバーの移動機構 — ローバーの車輪構造やロッカーボギー機構の基礎

また、以下の数学的知識があるとスムーズに読み進められます。

• 積分の基礎（定積分の計算）
• 三角関数の基本的な性質

テラメカニクスとは — 土壌力学と車両力学の融合

「固い地面」の常識が通用しない世界

私たちが日常で使う自動車は、アスファルトやコンクリートといった「変形しない地面」の上を走ることを前提に設計されています。タイヤと路面の間の摩擦力で加速・減速し、路面はほとんど変形しません。この場合、タイヤの力学だけを考えればよく、地面側の力学は無視できます。

しかし、月面のレゴリス、火星の砂地、地球上でも泥濘地や砂浜のような軟弱地盤の上では話が全く異なります。車輪は地面にめり込み、地面は車輪の荷重で圧縮され、車輪が回転すると地面のせん断変形（ずれ）が発生します。車輪と地面が互いに影響し合う「相互作用」が走行性能を支配するのです。

この相互作用を体系的に扱う学問がテラメカニクス（Terramechanics）です。ラテン語の “terra”（大地）と “mechanics”（力学）を組み合わせた造語で、1950年代に M.G. Bekker がその基礎を確立しました。

テラメカニクスの3つの柱

テラメカニクスは、大きく3つの要素から構成されます。

第1の柱: 土壌の力学特性。地面がどのような力学的性質を持つかを定量化します。具体的には、荷重をかけたときにどれだけ沈むか（圧力-沈下特性）、横方向の力をかけたときにどれだけずれるか（せん断特性）を実験的に測定し、パラメータとして表現します。

第2の柱: 車輪-土壌の接触力学。車輪が土壌に接触している領域で、圧力（法線応力）とせん断応力がどのように分布するかを解析します。車輪の形状、荷重、回転速度、スリップ率などが入力となり、接触面での応力分布が出力されます。

第3の柱: 車両全体の走行性能。個々の車輪の力学を統合して、車両全体の牽引力、走行抵抗、登坂能力などを予測します。ローバーの設計では、この段階でミッション要求（「傾斜15度の砂地を登れるか」など）を満たすかどうかを判定します。

本記事では、第1と第2の柱に焦点を当て、車輪1つと土壌の相互作用を詳しく解析します。まずは第1の柱、すなわち土壌の圧力-沈下特性を定量化する Bekker モデルから始めましょう。

Bekker の圧力-沈下モデル — 車輪はどこまで沈むのか

板が砂に沈む実験をイメージする

Bekker モデルの本質を理解するために、簡単な実験をイメージしてみましょう。砂場に平らな板を置き、上からゆっくりと荷重を増やしていきます。最初は板はほとんど沈みませんが、荷重を増やすにつれて沈下量が増えていきます。この「圧力と沈下の関係」をグラフに描くと、指数関数的な曲線が得られます。砂を押し固めるほど抵抗が大きくなるため、圧力の増加に対して沈下が頭打ちになっていくのです。

Bekker は、この圧力-沈下関係を以下のシンプルな経験式で表現しました。

Bekker の圧力-沈下式

板の幅を $b$、沈下量を $z$ としたとき、圧力 $p$ は次式で与えられます。

$$ p = \left(\frac{k_c}{b} + k_\phi\right) z^n $$

各パラメータの意味は以下の通りです。

この式の形を直感的に理解しましょう。まず $z^n$ の部分は、沈下量が大きくなるほど圧力が増加することを表しています。指数 $n$ が 1 なら線形、1 より大きければ沈下に対して圧力が急激に増加する（硬い土壌）、1 より小さければ圧力の増加が緩やか（柔らかい土壌）です。

次に $k_c/b + k_\phi$ の部分は、板の幅 $b$ への依存性を表しています。幅の広い板ほど $k_c/b$ の寄与が小さくなり、$k_\phi$ が支配的になります。これは物理的に、幅広の接地面では土壌の内部摩擦角による抵抗が主要因になることに対応しています。逆に幅の狭い接地面では、土壌の粘着力 $k_c$ の寄与が相対的に大きくなります。

沈下パラメータの物理的意味

3つのパラメータ $k_c$、$k_\phi$、$n$ は、土壌の種類によって大きく異なります。いくつかの代表的な土壌のパラメータ値を以下に示します。

乾燥砂では $k_c$ が非常に小さく $k_\phi$ が大きいことから、砂の抵抗は主に粒子間の摩擦によるものであることが読み取れます。一方、粘土質土壌では $k_c$ が大きく、粘着力の寄与が無視できないことがわかります。

剛体車輪の沈下量の導出

ここからが本題です。Bekker の圧力-沈下式を使って、剛体の円筒車輪が軟弱地盤の上に置かれたときの沈下量を求めましょう。

車輪の半径を $r$、幅を $b$、車輪にかかる垂直荷重を $W$ とします。車輪が土壌に沈下量 $z_0$ だけめり込んだとき、車輪と土壌の接触面での幾何学的関係を考えます。

車輪の最下点からの角度を $\theta$ とすると、角度 $\theta$ における局所的な沈下量 $z(\theta)$ は幾何学的に次のように表されます。

$$ z(\theta) = r(\cos\theta – \cos\theta_1) $$

ここで $\theta_1$ は車輪と土壌の接触が始まる角度（接触角の最大値）で、車輪の最大沈下量 $z_0$ と以下の関係があります。

$$ z_0 = r(1 – \cos\theta_1) $$

この関係を導きましょう。車輪の中心を原点にとると、車輪の最下点は中心から距離 $r$ だけ下にあります。角度 $\theta$ の位置における車輪表面の高さは $-r\cos\theta$ です。一方、変形前の地表面は車輪最下点から $z_0$ だけ上にあり、高さは $-(r – z_0)$ です。したがって、角度 $\theta$ における局所沈下量は以下のようになります。

地表面の高さから車輪表面の高さを引くと、

$$ z(\theta) = -(r – z_0) – (-r\cos\theta) = r\cos\theta – r + z_0 $$

$z_0 = r(1 – \cos\theta_1)$ を代入すると、

$$ z(\theta) = r\cos\theta – r + r(1 – \cos\theta_1) = r(\cos\theta – \cos\theta_1) $$

が得られます。接触が始まる角度 $\theta_1$ では $z(\theta_1) = 0$ となることも確認できます。

次に、車輪にかかる垂直荷重 $W$ と接触面の圧力分布の関係を立てます。角度 $\theta$ における圧力 $p(\theta)$ の垂直方向成分を、接触領域 $(-\theta_1 \leq \theta \leq \theta_1)$ にわたって積分すると荷重 $W$ に等しくなります。

$$ W = rb \int_{-\theta_1}^{\theta_1} p(\theta) \cos\theta \, d\theta $$

車輪の対称性から、この積分は次のように書き換えられます。

$$ W = 2rb \int_{0}^{\theta_1} p(\theta) \cos\theta \, d\theta $$

ここで Bekker の圧力-沈下式を代入します。$p(\theta) = (k_c/b + k_\phi) z(\theta)^n$ より、

$$ W = 2rb \left(\frac{k_c}{b} + k_\phi\right) \int_{0}^{\theta_1} [r(\cos\theta – \cos\theta_1)]^n \cos\theta \, d\theta $$

整理すると、

$$ W = 2r^{n+1}b \left(\frac{k_c}{b} + k_\phi\right) \int_{0}^{\theta_1} (\cos\theta – \cos\theta_1)^n \cos\theta \, d\theta $$

この積分は一般的な $n$ に対しては解析的に解けないため、数値積分で求めます。この式は、与えられた荷重 $W$ に対して接触角 $\theta_1$（したがって沈下量 $z_0$）を決定する暗黙的な方程式です。数値的に $\theta_1$ を求め、$z_0 = r(1 – \cos\theta_1)$ で沈下量を算出します。

ただし、沈下量が車輪半径に比べて小さい場合（$z_0 \ll r$、つまり $\theta_1 \ll 1$）には、小角度近似 $\cos\theta \approx 1 – \theta^2/2$ を使って解析的な近似解を得ることができます。Bekker はこの近似の下で次の簡易式を導きました。

$$ z_0 \approx \left(\frac{3W}{b(2r)^{1/2}\left(\frac{k_c}{b} + k_\phi\right)}\right)^{\frac{2}{2n+1}} $$

この式は「車輪の沈下は荷重に対して $2/(2n+1)$ 乗で増加する」ことを示しています。$n = 1$ のとき $z_0 \propto W^{2/3}$ となり、荷重を2倍にしても沈下は $2^{2/3} \approx 1.59$ 倍にしかなりません。これは、沈下が進むにつれて接触面積が増え、圧力の集中が緩和されるためです。

では、この圧力-沈下モデルを Python で実装し、車輪パラメータと土壌パラメータを変えたときの沈下量の変化を可視化してみましょう。

Python で沈下量を計算する

以下のコードでは、数値積分を用いて荷重と沈下量の関係を正確に計算し、Bekker の近似解と比較します。

import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq
import matplotlib.pyplot as plt

# --- 土壌パラメータ（乾燥砂）---
kc = 0.99e3      # N/m^(n+1)
kphi = 1528.4e3   # N/m^(n+2)
n = 1.1

# --- 車輪パラメータ ---
r = 0.15   # 車輪半径 [m]（火星ローバー想定）
b = 0.16   # 車輪幅 [m]

# Bekker 係数（幅を考慮）
k_eq = kc / b + kphi

def vertical_load(theta1, r, b, k_eq, n):
    """接触角 theta1 に対する垂直荷重を数値積分で計算"""
    def integrand(theta):
        return (np.cos(theta) - np.cos(theta1))**n * np.cos(theta)
    integral, _ = quad(integrand, 0, theta1)
    return 2 * r**(n + 1) * b * k_eq * integral

def find_theta1(W, r, b, k_eq, n, theta_max=np.pi/3):
    """荷重 W に対する接触角 theta1 を数値的に求める"""
    func = lambda th: vertical_load(th, r, b, k_eq, n) - W
    # theta1=0 では荷重0、theta1=theta_max で十分大きな荷重を仮定
    return brentq(func, 1e-6, theta_max)

def bekker_approx_sinkage(W, r, b, k_eq, n):
    """Bekker の近似式による沈下量"""
    return (3 * W / (b * np.sqrt(2 * r) * k_eq))**(2 / (2 * n + 1))

# --- 荷重を変化させて沈下量を計算 ---
W_array = np.linspace(10, 300, 100)  # 荷重 [N]
z0_numerical = []
z0_approx = []

for W in W_array:
    theta1 = find_theta1(W, r, b, k_eq, n)
    z0 = r * (1 - np.cos(theta1))
    z0_numerical.append(z0 * 1000)  # mm に変換
    z0_approx.append(bekker_approx_sinkage(W, r, b, k_eq, n) * 1000)

# --- 可視化 ---
plt.figure(figsize=(9, 5))
plt.plot(W_array, z0_numerical, 'b-', linewidth=2, label='Numerical (exact)')
plt.plot(W_array, z0_approx, 'r--', linewidth=2, label='Bekker approx.')
plt.xlabel('Vertical Load W [N]', fontsize=12)
plt.ylabel('Sinkage $z_0$ [mm]', fontsize=12)
plt.title('Wheel Sinkage vs. Vertical Load (Dry Sand)', fontsize=13)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

上のグラフから、以下の点が読み取れます。

1. 荷重の増加に対して沈下量は非線形に増加する — 荷重を2倍にしても沈下は2倍にはなりません。これは接触面積の増大による圧力の再分配効果を反映しています。$n = 1.1$ の乾燥砂では、沈下は荷重の約 $0.625$ 乗（$= 2/(2 \times 1.1 + 1)$）に比例します。
2. Bekker の近似式は小沈下量では精度が良い — 荷重が小さく沈下が浅い領域（グラフの左側）では、近似解と数値解がほぼ一致しています。一方、荷重が大きくなると近似解と数値解の差が開いていきます。これは小角度近似の妥当性が崩れるためです。
3. 火星ローバー程度の荷重（50〜100N）では沈下は数mm程度 — Curiosity ローバーの各車輪荷重は約145N ですが、乾燥砂上の沈下量は10mm 未満と予測されます。ただしこれは静的な沈下であり、走行時のスリップによる追加沈下は含まれていません。

Bekker モデルにより「どれだけ沈むか」がわかるようになりました。しかし、車輪が前に進むためには地面を「掴んで蹴る」必要があります。次は、車輪と土壌の間のせん断力 — つまり推進力の源を定量化する Janosi-Hanamoto モデルを見ていきましょう。

Janosi-Hanamoto のせん断応力モデル — 地面を掴む力

車輪はなぜ前に進めるのか

舗装路の上では、タイヤのゴムと路面の間の摩擦力が車を前に進ませます。しかし砂地や泥の上では、車輪が土壌にめり込み、土壌そのものをせん断（ずらす）ことで推進力を得ます。直感的には、スコップで砂を掘るときに砂が崩れるのと同じ現象です。車輪の回転によって車輪-土壌境界で土壌がせん断変形を受け、その抵抗力が車輪を前に押す力になります。

このせん断抵抗力は無限に大きくなるわけではありません。砂を横に押し続けると、最初は抵抗が増えますが、やがて砂が完全にずれて（破壊されて）抵抗は一定値に落ち着きます。この「せん断変位とせん断応力の関係」をモデル化したのが Janosi-Hanamoto の式です。

Mohr-Coulomb の破壊基準

せん断応力モデルの前提として、土壌のせん断強度を規定する Mohr-Coulomb の破壊基準 を確認しておきましょう。土壌が破壊（せん断破壊）するときの最大せん断応力 $\tau_{\max}$ は、法線応力 $\sigma$ に対して以下の線形関係で与えられます。

$$ \tau_{\max} = c + \sigma \tan\phi $$

ここで $c$ は土壌の粘着力（cohesion）、$\phi$ は内部摩擦角（internal friction angle）です。

この式は直感的に理解できます。粘着力 $c$ は、土壌粒子同士をくっつけている力（湿った砂が形を保つのはこの力のおかげ）であり、荷重がゼロでも一定のせん断抵抗を発揮します。一方、$\sigma \tan\phi$ の項は、荷重が大きいほど粒子同士が強く噛み合い、ずれにくくなる効果を表しています。乾燥砂では $c \approx 0$ で、摩擦角 $\phi$ のみがせん断強度を決定します。

Janosi-Hanamoto の式

Mohr-Coulomb は「最終的な強度」を与えますが、車輪の走行では土壌が「徐々にずれていく過程」が重要です。Janosi と Hanamoto は、せん断変位 $j$ に対するせん断応力 $\tau$ の発達過程を以下の指数関数モデルで表現しました。

$$ \tau = \tau_{\max}\left(1 – e^{-j/K}\right) = (c + \sigma \tan\phi)\left(1 – e^{-j/K}\right) $$

ここで $K$ はせん断変形係数（shear deformation modulus）と呼ばれ、単位は長さ [m] です。

この式が表す物理を理解しましょう。せん断変位 $j$ がゼロのとき $\tau = 0$（まだ土壌はずれていない）、$j$ が増加するにつれて $\tau$ は指数関数的に増加し、$j \gg K$ で $\tau \to \tau_{\max}$ に漸近します。パラメータ $K$ は「せん断応力が最大値の $(1 – 1/e) \approx 63\%$ に達するまでに必要な変位量」を意味します。

$K$ が小さい土壌（例えば密に締め固められた砂）では、わずかなずれで最大強度に達します。$K$ が大きい土壌（例えば緩い砂やローバーが踏み荒らした砂）では、大きなずれが必要です。

代表的な土壌のせん断パラメータを以下に示します。

月面レゴリスの粘着力が非常に小さいのは、大気がなく水分も存在しないため、粒子間の粘着力がほぼゼロに近いことを反映しています。一方、内部摩擦角は35度とやや大きく、これは角張ったレゴリス粒子同士の噛み合わせ効果によるものです。

車輪回転におけるせん断変位

車輪が回転しながら前進するとき、車輪-土壌接触面の各点でのせん断変位 $j$ を計算する必要があります。ここでスリップ率の概念が登場しますが、詳しくは次のセクションで扱います。ここでは、角度 $\theta$ における局所せん断変位 $j(\theta)$ が以下で与えられることを示しておきます。

$$ j(\theta) = r[\theta_1 – \theta – (1 – s)(\sin\theta_1 – \sin\theta)] $$

ここで $s$ はスリップ率（次セクションで定義）です。

この式の導出を簡潔に示します。車輪上の一点が接触領域に入る瞬間（角度 $\theta_1$）から角度 $\theta$ まで移動する間に、車輪表面の移動距離と車輪中心の前進距離に差が生じます。車輪表面の移動距離は $r(\theta_1 – \theta)$（弧の長さ）、車輪中心の水平移動距離は $r(1-s)(\sin\theta_1 – \sin\theta)$ です。両者の差がせん断変位となります。

$$ j(\theta) = r(\theta_1 – \theta) – r(1 – s)(\sin\theta_1 – \sin\theta) $$

スリップがゼロ ($s = 0$) の場合、車輪の周速度と前進速度が一致し、$j(\theta)$ は小さな値になります。スリップが大きい ($s \to 1$) 場合、車輪は空転に近い状態となり、$j(\theta) \approx r(\theta_1 – \theta)$ と大きなせん断変位が発生します。

Janosi-Hanamoto モデルにより、車輪の各接触点でのせん断応力が計算できるようになりました。このせん断応力が推進力の源ですが、同時に車輪の沈下は走行抵抗を生みます。次は、沈下と走行抵抗、そしてスリップ率と牽引力の関係を定量的に見ていきましょう。

スリップ率の定義と物理的意味

スリップとは何か

自動車を運転していて、凍結した道路でアクセルを踏み込んだとき、タイヤが空転して車がなかなか前に進まない経験をしたことはないでしょうか。車輪が1回転する間に、硬い路面なら $2\pi r$ だけ前進するはずですが、滑りやすい路面ではそれより短い距離しか進みません。この「理想的な前進距離と実際の前進距離のずれ」を定量化するのがスリップ率です。

スリップ率の数学的定義

スリップ率 $s$ は、車輪の角速度を $\omega$、車輪の半径を $r$、車両の実際の前進速度を $V$ として次のように定義されます。

$$ s = 1 – \frac{V}{r\omega} $$

この定義を直感的に理解しましょう。

• $s = 0$（スリップなし）: $V = r\omega$ で、車輪の周速度と前進速度が完全に一致します。硬い路面での理想的な走行状態です
• $0 < s < 1$（正のスリップ）: $V < r\omega$ で、車輪が前進距離に対して「余分に」回転しています。砂地での駆動時に典型的です
• $s = 1$（完全空転）: $V = 0$ で、車輪は回転しているが車両は全く前進しません。最悪のスタック状態です
• $s < 0$（負のスリップ）: $V > r\omega$ で、車輪が前進距離に対して「足りない」回転をしています。制動（ブレーキ）時に生じます

惑星探査ローバーの運用では、スリップ率は極めて重要な指標です。Spirit ローバーでは車輪オドメトリとビジュアルオドメトリの比較からスリップ率をリアルタイムで推定し、スリップ率が一定値（例えば40%）を超えるとコマンドで走行を停止させるルールが運用されていました。にもかかわらず Spirit がスタックしたのは、後述するように複合的な要因が重なったためです。

スリップ率と牽引力の関係

車輪が駆動されてスリップが生じると、車輪-土壌接触面でせん断変位が発生し、Janosi-Hanamoto モデルに基づいてせん断応力（推進力の源）が生まれます。スリップ率が大きいほどせん断変位が大きくなり、当然ながら牽引力も変化します。

車輪が発揮する総牽引力（Gross Traction Force）$F$ は、接触面のせん断応力を水平方向に積分して求められます。

$$ F = rb \int_{-\theta_1}^{\theta_1} \tau(\theta) \, d\theta = 2rb \int_{0}^{\theta_1} \tau(\theta) \, d\theta $$

ここで $\tau(\theta) = (c + p(\theta)\tan\phi)(1 – e^{-j(\theta)/K})$ です。

$p(\theta)$ は Bekker モデルから計算される局所圧力、$j(\theta)$ は前セクションで導出したせん断変位です。この積分は解析的に求められないため、数値積分を使います。

スリップ率と牽引力の典型的な関係は以下の特徴を持ちます。

1. $s = 0$ では $j \approx 0$ なので $\tau \approx 0$、牽引力はほぼゼロ
2. スリップ率の増加とともに牽引力は急激に増加
3. $s$ が $0.15 \sim 0.3$ 程度でピークに達する（最適スリップ率）
4. それ以上のスリップでは土壌の過度な変形により牽引力が飽和または緩やかに減少

この特徴は「最適スリップ率」の存在を意味しています。ローバーの制御においては、不必要に大きなスリップを生じさせず、最適スリップ率付近で走行させることがエネルギー効率の面で重要です。

車輪の推進力の仕組みがわかったところで、次は推進を妨げる力 — 走行抵抗について考えましょう。

沈下と走行抵抗 — 沈むほど走りにくくなる力学

なぜ砂浜では歩きにくいのか

砂浜を歩くと、硬い道路を歩くよりもずっと疲れます。これは足が沈み込むたびに、砂を押しのけて圧縮するためにエネルギーを消費しているからです。車輪の場合も同様に、土壌に沈下する際に圧密仕事（compaction work）を行い、これが走行抵抗の主要な成分になります。

圧密抵抗（Compaction Resistance）

車輪が前進するとき、車輪の前方で土壌が圧縮され、後方に轍（わだち）が残ります。車輪の沈下量 $z_0$ の轍を単位距離だけ残すために必要なエネルギーは、Bekker の圧力-沈下式を沈下量について積分することで得られます。

まず、幅 $b$ の車輪が単位前進距離あたりに行う圧密仕事を求めます。車輪が $dx$ だけ前進するとき、幅 $b$ の帯状の土壌が深さ $z_0$ まで圧縮されます。この圧密仕事は次のように計算できます。

$$ \frac{dE}{dx} = b \int_{0}^{z_0} p(z) \, dz $$

Bekker の式 $p(z) = (k_c/b + k_\phi) z^n$ を代入して積分を実行すると、

$$ \frac{dE}{dx} = b \left(\frac{k_c}{b} + k_\phi\right) \frac{z_0^{n+1}}{n+1} $$

この単位距離あたりの仕事が、走行抵抗力に等しくなります（仕事 = 力 $\times$ 距離 の関係から）。したがって圧密抵抗 $R_c$ は次式で与えられます。

$$ R_c = \frac{b\left(\frac{k_c}{b} + k_\phi\right) z_0^{n+1}}{n+1} $$

この式から重要な知見が得られます。走行抵抗は沈下量の $(n+1)$ 乗に比例するため、沈下量が2倍になると走行抵抗は $2^{n+1}$ 倍に急増します。$n = 1.1$ の乾燥砂では $2^{2.1} \approx 4.3$ 倍です。砂地ではわずかな沈下増加が走行抵抗を劇的に悪化させることがわかります。

ブルドージング抵抗

車輪がスリップしながら前進すると、車輪の前方に土壌が盛り上がるブルドージング効果が発生します。これは車輪が沈み込んだ状態で回転する際に、前方の土壌を押し上げる必要があるためです。ブルドージング抵抗 $R_b$ は、スリップ率が高いほど大きくなり、以下の近似式で推定されます。

$$ R_b \approx b \cdot z_0 \cdot (c \cot\phi + \frac{1}{2} \gamma z_0) \cdot \frac{s}{1-s} \cdot \frac{1}{\cos^2(45° – \phi/2)} $$

ここで $\gamma$ は土壌の単位体積重量です。この式は厳密ではありませんが、ブルドージング抵抗がスリップ率に大きく依存することを示しています。$s$ が 1 に近づくとブルドージング抵抗が発散するため、高スリップ状態でのスタックの一因となります。

正味牽引力（Net Drawbar Pull）

車輪が実際に車両を前に引っ張る力は、総牽引力から走行抵抗を差し引いた正味牽引力（Net Drawbar Pull）$DP$ です。

$$ DP = F – R_c – R_b $$

$DP > 0$ であれば車両は前進でき、$DP = 0$ では走行限界、$DP < 0$ では車輪がスタックすることを意味します。

正味牽引力がゼロになる条件は、テラメカニクスにおいて最も実用的な情報の一つです。ローバーの設計では、想定される最悪の土壌条件（最も柔らかいレゴリス）と最悪の地形条件（最大傾斜角）において $DP > 0$ が保証されることを確認します。

ここまでの力学モデルを統合して、スリップ率と正味牽引力の関係を Python でシミュレーションしてみましょう。

Python による車輪-土壌相互作用シミュレーション

統合モデルの実装

以下のコードでは、Bekker の圧力-沈下モデルと Janosi-Hanamoto のせん断モデルを組み合わせ、スリップ率に対する牽引力と走行抵抗を計算します。

import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq
import matplotlib.pyplot as plt

# ===== 土壌パラメータ（乾燥砂）=====
kc = 0.99e3        # N/m^(n+1)
kphi = 1528.4e3     # N/m^(n+2)
n_soil = 1.1        # 沈下指数
c = 1.0e3           # 粘着力 [Pa]
phi = np.radians(30)  # 内部摩擦角 [rad]
K_shear = 0.025     # せん断変形係数 [m]

# ===== 車輪パラメータ =====
r = 0.15      # 車輪半径 [m]
b = 0.16      # 車輪幅 [m]
W = 100.0     # 垂直荷重 [N]（火星重力を考慮した値）

k_eq = kc / b + kphi  # 等価沈下係数

def compute_theta1(W, r, b, k_eq, n_soil):
    """荷重 W に対する接触角 theta1 を数値的に求める"""
    def load_eq(theta1):
        def integrand(theta):
            return (np.cos(theta) - np.cos(theta1))**n_soil * np.cos(theta)
        integral, _ = quad(integrand, 0, theta1)
        return 2 * r**(n_soil + 1) * b * k_eq * integral - W
    return brentq(load_eq, 1e-6, np.pi / 3)

def compute_forces(slip, theta1, r, b, k_eq, n_soil, c, phi, K_shear):
    """
    スリップ率 slip に対する牽引力と圧密抵抗を計算
    """
    # --- 総牽引力（せん断応力の積分）---
    def traction_integrand(theta):
        # 局所圧力
        z_local = r * (np.cos(theta) - np.cos(theta1))
        p_local = k_eq * z_local**n_soil
        # 局所せん断変位
        j_local = r * (theta1 - theta - (1 - slip) * (np.sin(theta1) - np.sin(theta)))
        j_local = max(j_local, 0)
        # Janosi-Hanamoto のせん断応力
        tau_max = c + p_local * np.tan(phi)
        tau = tau_max * (1 - np.exp(-j_local / K_shear))
        return tau

    F_traction, _ = quad(traction_integrand, 0, theta1)
    F_traction *= 2 * r * b

    # --- 圧密抵抗 ---
    z0 = r * (1 - np.cos(theta1))
    R_c = b * k_eq * z0**(n_soil + 1) / (n_soil + 1)

    return F_traction, R_c

# --- 接触角の計算 ---
theta1 = compute_theta1(W, r, b, k_eq, n_soil)
z0 = r * (1 - np.cos(theta1))
print(f"接触角: {np.degrees(theta1):.2f} deg")
print(f"沈下量: {z0 * 1000:.2f} mm")

# --- スリップ率を変化させて力を計算 ---
slips = np.linspace(0.01, 0.95, 80)
F_list = []
Rc_list = []
DP_list = []

for s in slips:
    F, Rc = compute_forces(s, theta1, r, b, k_eq, n_soil, c, phi, K_shear)
    F_list.append(F)
    Rc_list.append(Rc)
    DP_list.append(F - Rc)

F_arr = np.array(F_list)
Rc_arr = np.array(Rc_list)
DP_arr = np.array(DP_list)

# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左: 力の内訳
axes[0].plot(slips * 100, F_arr, 'b-', linewidth=2, label='Gross Traction $F$')
axes[0].plot(slips * 100, Rc_arr * np.ones_like(slips), 'r--', linewidth=2,
             label='Compaction Resistance $R_c$')
axes[0].plot(slips * 100, DP_arr, 'g-', linewidth=2.5, label='Net Drawbar Pull $DP$')
axes[0].axhline(y=0, color='k', linestyle=':', alpha=0.5)
axes[0].set_xlabel('Slip Ratio [%]', fontsize=12)
axes[0].set_ylabel('Force [N]', fontsize=12)
axes[0].set_title('Traction Force vs. Slip Ratio (Dry Sand)', fontsize=13)
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

# 右: 正味牽引力のみ拡大
axes[1].plot(slips * 100, DP_arr, 'g-', linewidth=2.5)
axes[1].axhline(y=0, color='r', linestyle='--', alpha=0.7, label='DP = 0 (stall)')
# 最適スリップ率の表示
idx_max = np.argmax(DP_arr)
axes[1].axvline(x=slips[idx_max] * 100, color='purple', linestyle=':', alpha=0.7,
                label=f'Optimal slip = {slips[idx_max]*100:.0f}%')
axes[1].scatter([slips[idx_max] * 100], [DP_arr[idx_max]], color='purple',
                s=80, zorder=5)
axes[1].set_xlabel('Slip Ratio [%]', fontsize=12)
axes[1].set_ylabel('Net Drawbar Pull [N]', fontsize=12)
axes[1].set_title('Net Drawbar Pull vs. Slip Ratio', fontsize=13)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\n最適スリップ率: {slips[idx_max]*100:.1f}%")
print(f"最大正味牽引力: {DP_arr[idx_max]:.2f} N")

このシミュレーション結果から、以下の重要な知見が読み取れます。

1. 牽引力はスリップ率の増加とともに飽和する — 左図で青い線（総牽引力）はスリップ率10〜20%で急激に立ち上がり、その後は緩やかな増加に転じます。これは Janosi-Hanamoto モデルの指数関数的飽和（$1 – e^{-j/K}$）を反映しています。スリップが大きくなればせん断変位は増えますが、せん断応力は最大値 $\tau_{\max}$ で頭打ちになります。
2. 走行抵抗は静的沈下のみを考慮した場合は一定 — 赤い破線の圧密抵抗はスリップ率によらず一定です。これは簡略化のため動的沈下（スリップによる追加沈下）を無視しているためで、実際にはスリップが増えると沈下も増加し、走行抵抗はさらに大きくなります。
3. 正味牽引力にはピーク（最適スリップ率）が存在する — 右図で緑の曲線にピークが確認できます。このスリップ率より小さいと牽引力が不足し、大きいと走行抵抗の増加やエネルギーの浪費が生じます。ローバーの走行制御ではこの最適点付近に制御することが望ましいです。

複数の土壌タイプの比較

次に、異なる土壌タイプ（乾燥砂、砂質ローム、月面レゴリス推定値）での正味牽引力の比較を行います。

import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq
import matplotlib.pyplot as plt

# --- 車輪パラメータ ---
r = 0.15   # 車輪半径 [m]
b = 0.16   # 車輪幅 [m]
W = 100.0  # 垂直荷重 [N]

# --- 土壌パラメータセット ---
soils = {
    'Dry Sand': {
        'kc': 0.99e3, 'kphi': 1528.4e3, 'n': 1.1,
        'c': 1.0e3, 'phi': np.radians(30), 'K': 0.025
    },
    'Sandy Loam': {
        'kc': 5.27e3, 'kphi': 1515.0e3, 'n': 0.7,
        'c': 1.7e3, 'phi': np.radians(29), 'K': 0.025
    },
    'Lunar Regolith (est.)': {
        'kc': 1.4e3, 'kphi': 820.0e3, 'n': 1.0,
        'c': 0.17e3, 'phi': np.radians(35), 'K': 0.018
    },
}

def compute_all(W, r, b, soil_params):
    """土壌パラメータに対する接触角と力を計算"""
    kc = soil_params['kc']
    kphi = soil_params['kphi']
    n_s = soil_params['n']
    c_s = soil_params['c']
    phi_s = soil_params['phi']
    K_s = soil_params['K']
    k_eq = kc / b + kphi

    # 接触角
    def load_eq(theta1):
        def integrand(theta):
            return (np.cos(theta) - np.cos(theta1))**n_s * np.cos(theta)
        integ, _ = quad(integrand, 0, theta1)
        return 2 * r**(n_s + 1) * b * k_eq * integ - W
    theta1 = brentq(load_eq, 1e-6, np.pi / 3)
    z0 = r * (1 - np.cos(theta1))

    # スリップごとの力
    slips = np.linspace(0.01, 0.90, 60)
    DP_vals = []
    for s in slips:
        def trac_int(theta):
            z_l = r * (np.cos(theta) - np.cos(theta1))
            p_l = k_eq * z_l**n_s
            j_l = r * (theta1 - theta - (1 - s) * (np.sin(theta1) - np.sin(theta)))
            j_l = max(j_l, 0)
            tau = (c_s + p_l * np.tan(phi_s)) * (1 - np.exp(-j_l / K_s))
            return tau
        F, _ = quad(trac_int, 0, theta1)
        F *= 2 * r * b
        R_c = b * k_eq * z0**(n_s + 1) / (n_s + 1)
        DP_vals.append(F - R_c)
    return slips, np.array(DP_vals), z0

# --- 計算と可視化 ---
plt.figure(figsize=(10, 6))
colors = ['#1f77b4', '#ff7f0e', '#2ca02c']

for (name, params), color in zip(soils.items(), colors):
    slips, DP, z0 = compute_all(W, r, b, params)
    plt.plot(slips * 100, DP, linewidth=2.5, color=color,
             label=f'{name} (sinkage={z0*1000:.1f} mm)')

plt.axhline(y=0, color='k', linestyle=':', alpha=0.5)
plt.xlabel('Slip Ratio [%]', fontsize=12)
plt.ylabel('Net Drawbar Pull [N]', fontsize=12)
plt.title('Net Drawbar Pull Comparison across Soil Types\n'
          f'(Wheel: r={r*100:.0f}cm, b={b*100:.0f}cm, W={W:.0f}N)',
          fontsize=13)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

このグラフから、土壌タイプによる走行性能の違いが明確に読み取れます。

1. 砂質ロームが最も大きな正味牽引力を発揮する — 粘着力 $c$ と摩擦角 $\phi$ の両方がバランスよく高く、かつ沈下指数 $n = 0.7$ が小さいため沈下も深いですが、せん断強度が高い土壌です。沈下は大きいものの、走行抵抗を上回る十分な牽引力が得られます。
2. 月面レゴリスは粘着力がほぼゼロだが摩擦角が大きい — 粘着力 $c = 0.17$ kPa は他の土壌に比べて一桁以上小さいですが、内部摩擦角 $\phi = 35°$ が大きいため、荷重がかかった状態ではそれなりの牽引力を発揮します。ただし、低荷重（月面の低重力環境）では摩擦の寄与も小さくなるため注意が必要です。
3. すべての土壌で15〜30%のスリップ率がピーク領域 — 最適スリップ率はおおよそ15〜30%の範囲にあります。これはローバーの走行制御におけるスリップ管理の目標値を決める重要な情報です。

土壌タイプによって牽引力が大きく異なることが確認できました。では、実際の惑星表面 — 月面と火星の表面 — はどのような土壌特性を持っているのでしょうか。次のセクションでレゴリスの物理特性を詳しく見ていきましょう。

月面・火星レゴリスの土壌特性

レゴリスとは

惑星や衛星の表面を覆う、岩石が風化・粉砕されてできた粒子状の層をレゴリス（Regolith）と呼びます。月面のレゴリスは数十億年にわたる微小隕石の衝突（インパクトガーデニング）によって生成された細かい粒子の層で、厚さは数メートルから十数メートルに及びます。火星のレゴリスも同様に存在しますが、大気の存在による風化プロセスも加わり、月面とは異なる特性を示します。

月面レゴリスの特性

Apollo ミッション（1969-1972年）では、宇宙飛行士が月面で直接土壌の力学試験を行い、月面車（Lunar Roving Vehicle, LRV）の走行データも取得されました。これらの貴重なデータから、月面レゴリスの土壌力学パラメータが決定されています。

粒径分布: 月面レゴリスの中央粒径は約60〜80 $\mu$m で、地球上の細かい砂に相当します。しかし粒子の形状は大きく異なります。月面では大気や水による風化がないため、粒子は角張った不規則な形状をしています。この角張りが粒子間の噛み合わせ（インターロッキング）を生み、比較的大きな内部摩擦角の原因となっています。

密度と空隙率: 表面付近のかさ密度は約1.5 g/cm$^3$、深さ30cmでは約1.75 g/cm$^3$ に増加します。真密度（粒子そのもの）は約3.1 g/cm$^3$ なので、空隙率は約40〜52%です。表面ほど緩く、深部ほど締まっているという勾配が存在します。

せん断パラメータ: Apollo の測定によると、粘着力 $c = 0.1 \sim 1.0$ kPa、内部摩擦角 $\phi = 25° \sim 50°$ と幅広い値が報告されています。この幅は、密度の勾配（表面は緩く深部は密）に対応しています。テラメカニクスのシミュレーションでは、表面付近の値として $c = 0.17$ kPa、$\phi = 35°$ がよく使われます。

重力の影響: 月面の重力は地球の約 $1/6$ です。車輪にかかる垂直荷重が小さくなるため、沈下量は地球上よりも小さくなります。しかし同時に、摩擦による牽引力の上限 $\sigma \tan\phi$ も小さくなるため、牽引力の確保が課題になります。結果的に、月面でのスリップ率は地球上よりも高くなりがちです。

火星レゴリスの特性

火星の地表はローバーミッション（Sojourner、Spirit、Opportunity、Curiosity、Perseverance）によって直接調査されてきましたが、Apollo のように人間が土壌力学試験を行ったわけではないため、パラメータの不確実性は月面より大きくなっています。

粒径と組成: 火星のレゴリスは、玄武岩質の砂（粒径約100〜250 $\mu$m）、細粒のダスト（粒径数 $\mu$m）、および様々なサイズの岩石が混在しています。風によるダストの堆積があり、一部の地域では非常に細かく柔らかいダスト層が存在します。

密度: 推定かさ密度は約1.1〜1.6 g/cm$^3$ で、場所による変動が大きいとされています。特にダストが堆積した場所では密度が低く、極めて柔らかい表面となります。

せん断パラメータ: Spirit と Opportunity の車輪沈下データからの推定では、粘着力 $c = 0 \sim 2$ kPa、内部摩擦角 $\phi = 15° \sim 35°$ とされています。場所による変動が月面以上に大きく、同じローバーのミッション中でも走行場所によって特性が劇的に変化することがあります。

重力の影響: 火星の重力は地球の約 $0.38$ 倍です。月面より重力は大きいため、月面ほどではありませんが、やはり牽引力の確保は重要な課題です。

地球との比較

テラメカニクスパラメータの比較を以下の表にまとめます。

低重力環境では荷重が小さくなるため沈下は減りますが、牽引力も同時に減少します。ローバーの設計では、この沈下と牽引力のバランスを適切に評価する必要があります。

月面・火星のレゴリス特性を把握したところで、テラメカニクスの実際の重要性を実感できる事例として、2009年の Spirit ローバーのスタック事故を力学的に分析してみましょう。

Spirit 砂地スタック事例 — テラメカニクスの教訓

事故の概要

2009年4月、NASA の Mars Exploration Rover（MER）Spirit は、火星の Home Plate と呼ばれる地質構造の南側にある “Troy” と名付けられた砂地で走行不能（スタック）に陥りました。約10ヶ月にわたる脱出の試みが行われましたが、最終的に Spirit は自力脱出を果たせず、2010年3月を最後に通信が途絶しました。

Spirit は2004年の着陸から5年以上にわたって火星表面を探査し、数々の科学的発見をもたらした優れたローバーでした。当初の設計寿命は90 sol（火星日）でしたが、予想をはるかに超えて運用が続けられていました。それだけに、砂地スタックによるミッション終了は痛恨の出来事でした。

Troy の地形と土壌特性

Troy は表面が一見すると硬い地殻（クラスト）で覆われていましたが、その下は非常に柔らかい硫酸塩を多く含む細粒土壌でした。Spirit がこの地殻を破って進入したとき、車輪は数センチ以上沈み込みました。

後の分析から、Troy の土壌は以下のような極端に不利な特性を持っていたことが推定されています。

• 低粘着力: $c \approx 0.1 \sim 0.3$ kPa — 非常に弱い粒子間結合
• 低摩擦角: $\phi \approx 15° \sim 20°$ — 通常の火星レゴリスの半分程度
• 高い沈下性: 表面クラストが崩れた後の下層土壌は、極めて圧縮しやすい

さらに状況を悪化させたのは、Spirit の6つの車輪のうち右前輪がすでに故障（2006年にモーター破損）しており、5輪走行を強いられていたことです。故障した右前輪はロックされたまま引きずられる状態であり、追加の走行抵抗となっていました。

テラメカニクスによる分析

Spirit のスタックをテラメカニクスの観点から分析しましょう。スタックが起きるのは正味牽引力 $DP$ がゼロ以下になる場合です。Troy の土壌が通常のレゴリスと比べてどれほど不利だったかを定量的に見ます。

まず、内部摩擦角が $30°$ から $18°$ に低下した場合の影響を考えます。せん断強度は $\tau_{\max} = c + \sigma\tan\phi$ で与えられるため、$\tan 30° = 0.577$ から $\tan 18° = 0.325$ への低下は、摩擦による牽引力が約44%も減少することを意味します。

次に、Spirit の5輪走行の影響です。右前輪が引きずられることで追加の抵抗が発生し、残り5つの車輪の正味牽引力が通常の $5/6$ に減少するだけでなく、ロックされた車輪の引きずり抵抗分がさらに差し引かれます。

これらの複合要因をシミュレーションで確認してみましょう。

import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq
import matplotlib.pyplot as plt

# --- 共通の車輪パラメータ（Spirit MER）---
r = 0.13     # 車輪半径 [m]（MER 車輪）
b = 0.16     # 車輪幅 [m]
g_mars = 3.72  # 火星重力加速度 [m/s^2]
m_rover = 185  # ローバー質量 [kg]

# 通常の6輪走行時: 各車輪荷重
W_normal = m_rover * g_mars / 6  # ≈ 114.7 N

# 5輪走行時: 機能する車輪の荷重増加（簡略化）
W_5wheel = m_rover * g_mars / 5  # ≈ 137.6 N

# --- 2つの土壌シナリオ ---
scenarios = {
    'Normal Mars Regolith\n(6 wheels)': {
        'kc': 1.0e3, 'kphi': 900.0e3, 'n': 1.0,
        'c': 0.6e3, 'phi': np.radians(30), 'K': 0.02,
        'W': W_normal, 'wheels': 6, 'drag': 0
    },
    'Troy Soil (Spirit)\n(5 wheels + drag)': {
        'kc': 0.8e3, 'kphi': 600.0e3, 'n': 0.8,
        'c': 0.2e3, 'phi': np.radians(18), 'K': 0.03,
        'W': W_5wheel, 'wheels': 5, 'drag': 20  # 故障車輪の引きずり抵抗 [N]
    }
}

def sim_scenario(params, r, b):
    kc = params['kc']
    kphi = params['kphi']
    n_s = params['n']
    c_s = params['c']
    phi_s = params['phi']
    K_s = params['K']
    W = params['W']
    n_wheels = params['wheels']
    drag = params['drag']
    k_eq = kc / b + kphi

    # 接触角
    def load_eq(th1):
        def integ(th):
            return (np.cos(th) - np.cos(th1))**n_s * np.cos(th)
        val, _ = quad(integ, 0, th1)
        return 2 * r**(n_s + 1) * b * k_eq * val - W
    theta1 = brentq(load_eq, 1e-6, np.pi / 2.5)
    z0 = r * (1 - np.cos(theta1))

    slips = np.linspace(0.01, 0.95, 80)
    DP_total = []
    for s in slips:
        def trac_int(theta):
            z_l = r * (np.cos(theta) - np.cos(theta1))
            p_l = k_eq * z_l**n_s
            j_l = r * (theta1 - theta - (1 - s) * (np.sin(theta1) - np.sin(theta)))
            j_l = max(j_l, 0)
            tau = (c_s + p_l * np.tan(phi_s)) * (1 - np.exp(-j_l / K_s))
            return tau
        F, _ = quad(trac_int, 0, theta1)
        F *= 2 * r * b
        R_c = b * k_eq * z0**(n_s + 1) / (n_s + 1)
        # 車両全体の正味牽引力
        DP_vehicle = n_wheels * (F - R_c) - drag
        DP_total.append(DP_vehicle)
    return slips, np.array(DP_total), z0

# --- 計算と可視化 ---
plt.figure(figsize=(10, 6))
colors = ['#2ca02c', '#d62728']
linestyles = ['-', '-']

for (name, params), color, ls in zip(scenarios.items(), colors, linestyles):
    slips, DP, z0 = sim_scenario(params, r, b)
    plt.plot(slips * 100, DP, color=color, linewidth=2.5, linestyle=ls,
             label=f'{name}\n  sinkage={z0*1000:.1f}mm/wheel')

plt.axhline(y=0, color='k', linestyle='--', linewidth=1.5, alpha=0.7,
            label='DP = 0 (stall condition)')
plt.fill_between(np.linspace(1, 95, 80), 0, -300, alpha=0.1, color='red')
plt.xlabel('Slip Ratio [%]', fontsize=12)
plt.ylabel('Total Net Drawbar Pull [N]', fontsize=12)
plt.title("Spirit Rover: Normal Regolith vs. Troy Soil Scenario", fontsize=13)
plt.legend(fontsize=9, loc='best')
plt.grid(True, alpha=0.3)
plt.ylim(-150, None)
plt.tight_layout()
plt.show()

このシミュレーション結果は、Spirit がなぜ Troy でスタックしたかを力学的に説明しています。

1. 通常のレゴリスでは十分な正味牽引力が得られる — 緑の曲線は広いスリップ率範囲で正の牽引力を示しており、通常の火星レゴリス上では6輪走行の Spirit は十分な走行能力を持っていたことがわかります。
2. Troy の条件では牽引力が劇的に低下する — 赤い曲線は正味牽引力が大幅に低下しており、高スリップ率でかろうじて正になるか、全域で負になる可能性を示しています。低い摩擦角、低い粘着力、5輪走行、引きずり抵抗という複合要因が重なった結果です。
3. 一度沈下が進むと脱出困難に — このシミュレーションは静的沈下のみを考慮していますが、実際にはスリップによる動的沈下が追加されます。車輪が空転するほど沈下が深まり、沈下が深まるほど走行抵抗が増え、さらに空転する — という正のフィードバックループに陥ります。

Spirit の教訓とテラメカニクスの役割

Spirit の事故は、テラメカニクスの重要性を痛感させる出来事でした。この教訓を受け、以降のミッションでは以下のような対策が講じられています。

車輪設計の改良: Curiosity ローバーの車輪は直径50cmと Spirit の26cmの約2倍に大型化され、接地面積を増やして沈下を低減しています。Bekker モデルの沈下式 $z_0 \propto W^{2/(2n+1)} / r^{1/(2n+1)}$ から、車輪半径を大きくすることが沈下低減に有効であることがわかります。

スリップ監視と自動停止: ビジュアルオドメトリによるスリップ率のリアルタイム推定が改良され、閾値超過時の自動停止がより確実に機能するようになりました。Perseverance では、地形の危険度を事前に評価する AI システムも搭載されています。

地面特性の事前評価: 軌道上からの地表観測データ（HiRISE 画像、熱慣性マップなど）を活用して、走行経路の土壌特性を事前に推定する手法が発展しています。熱慣性は地表粒子サイズと相関があり、細粒のダスト堆積地域を事前に特定するのに役立ちます。

Spirit の事例から学んだ教訓は、テラメカニクスの予測能力をミッション運用に活かすことの重要性です。では最後に、ここまでに学んだすべてのモデルを組み合わせた包括的なシミュレーションとして、傾斜地でのローバーの登坂限界を計算してみましょう。

傾斜地での登坂限界シミュレーション

傾斜面でのテラメカニクス

これまでは水平面での走行を考えてきましたが、実際の惑星表面には傾斜があります。傾斜角 $\alpha$ の斜面を登る場合、以下の修正が必要になります。

まず、車輪にかかる垂直荷重は斜面に垂直な成分のみになるため、

$$ W_\perp = \frac{mg\cos\alpha}{N_w} $$

となります。ここで $m$ はローバー質量、$g$ は重力加速度、$N_w$ は車輪数です。荷重が減少するため沈下は減りますが、同時に牽引力の上限も下がります。

一方、斜面に沿った方向に重力の成分が追加の走行抵抗として作用します。

$$ R_{\text{slope}} = mg\sin\alpha $$

したがって、傾斜面での正味牽引力の条件は次のようになります。

$$ DP_{\text{slope}} = \sum_{i=1}^{N_w} (F_i – R_{c,i}) – mg\sin\alpha > 0 $$

全車輪が同一条件であると仮定すれば、

$$ DP_{\text{slope}} = N_w(F – R_c) – mg\sin\alpha $$

この式から、登坂可能な最大傾斜角を「$DP_{\text{slope}} = 0$ となる $\alpha$」として求めることができます。以下のコードで、スリップ率と傾斜角の両方を変化させて、登坂可能領域を可視化します。

import numpy as np
from scipy.integrate import quad
from scipy.optimize import brentq
import matplotlib.pyplot as plt

# --- パラメータ（Curiosity級ローバー + 火星レゴリス）---
r = 0.25       # 車輪半径 [m]（Curiosity: 25cm）
b = 0.40       # 車輪幅 [m]（Curiosity: 40cm）
m = 899        # ローバー質量 [kg]（Curiosity: 899kg）
g = 3.72       # 火星重力加速度 [m/s^2]
Nw = 6         # 車輪数

# 火星レゴリス（中程度の硬さ）
kc = 1.0e3
kphi = 900.0e3
n_s = 1.0
c_s = 0.6e3
phi_s = np.radians(30)
K_s = 0.02
k_eq = kc / b + kphi

# --- 傾斜角とスリップ率のグリッド ---
alphas_deg = np.linspace(0, 30, 50)
slips = np.linspace(0.01, 0.70, 50)

DP_map = np.zeros((len(alphas_deg), len(slips)))

for i, alpha_deg in enumerate(alphas_deg):
    alpha = np.radians(alpha_deg)
    W_perp = m * g * np.cos(alpha) / Nw  # 各車輪の垂直荷重

    # 接触角
    try:
        def load_eq(th1):
            def integ(th):
                return (np.cos(th) - np.cos(th1))**n_s * np.cos(th)
            val, _ = quad(integ, 0, th1)
            return 2 * r**(n_s + 1) * b * k_eq * val - W_perp
        theta1 = brentq(load_eq, 1e-6, np.pi / 2.5)
    except ValueError:
        DP_map[i, :] = np.nan
        continue

    z0 = r * (1 - np.cos(theta1))
    R_c = b * k_eq * z0**(n_s + 1) / (n_s + 1)
    R_slope = m * g * np.sin(alpha)

    for j_idx, s in enumerate(slips):
        def trac_int(theta):
            z_l = r * (np.cos(theta) - np.cos(theta1))
            p_l = k_eq * z_l**n_s
            j_l = r * (theta1 - theta - (1 - s) * (np.sin(theta1) - np.sin(theta)))
            j_l = max(j_l, 0)
            tau = (c_s + p_l * np.tan(phi_s)) * (1 - np.exp(-j_l / K_s))
            return tau
        F, _ = quad(trac_int, 0, theta1)
        F *= 2 * r * b
        DP_map[i, j_idx] = Nw * (F - R_c) - R_slope

# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(15, 6))

# 左: ヒートマップ
S, A = np.meshgrid(slips * 100, alphas_deg)
cf = axes[0].contourf(S, A, DP_map, levels=30, cmap='RdYlGn')
axes[0].contour(S, A, DP_map, levels=[0], colors='black', linewidths=2.5)
plt.colorbar(cf, ax=axes[0], label='Net Drawbar Pull [N]')
axes[0].set_xlabel('Slip Ratio [%]', fontsize=12)
axes[0].set_ylabel('Slope Angle [deg]', fontsize=12)
axes[0].set_title('Drawbar Pull Map\n(Green=Go, Red=Stuck)', fontsize=13)

# 右: 特定スリップ率での牽引力 vs 傾斜角
slip_targets = [0.10, 0.20, 0.30, 0.50]
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']
for s_target, color in zip(slip_targets, colors):
    idx = np.argmin(np.abs(slips - s_target))
    axes[1].plot(alphas_deg, DP_map[:, idx], linewidth=2, color=color,
                 label=f'Slip = {s_target*100:.0f}%')
axes[1].axhline(y=0, color='k', linestyle='--', linewidth=1.5, alpha=0.7)
axes[1].set_xlabel('Slope Angle [deg]', fontsize=12)
axes[1].set_ylabel('Net Drawbar Pull [N]', fontsize=12)
axes[1].set_title('Drawbar Pull vs. Slope at Fixed Slip', fontsize=13)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

このシミュレーション結果から、傾斜地走行に関する以下の重要な知見が得られます。

1. 登坂限界は傾斜角20度前後 — 左図のヒートマップで、黒い実線（$DP = 0$ の等高線）が登坂限界を示しています。火星レゴリス（中程度の硬さ）で Curiosity 級のローバーの場合、適切なスリップ率（20〜30%）であれば傾斜20度程度まで登坂可能であることが予測されます。これは NASA が公表している Curiosity の登坂能力（最大約20度）とおおむね一致しています。
2. スリップ率が小さすぎても大きすぎても登坂能力は下がる — 右図でスリップ率10%の曲線は、30%の曲線よりも早くゼロを下回ります。スリップが小さいと十分な牽引力が発生しません。一方、スリップ率50%では牽引力は大きいですが、動的沈下（このモデルでは未考慮）によって実際にはさらに性能が悪化します。
3. 急斜面では垂直荷重の減少と重力抵抗の増加のダブルパンチ — 傾斜角が大きくなると、垂直荷重 $W_\perp = mg\cos\alpha/N_w$ が減少して牽引力が下がり、同時に重力抵抗 $mg\sin\alpha$ が増加します。これが正味牽引力を急速に低下させる原因です。

テラメカニクスの発展と課題

古典モデルの限界

ここまで解説してきた Bekker モデルと Janosi-Hanamoto モデルは、テラメカニクスの基礎として広く使われていますが、いくつかの重要な限界があります。

静的モデルの限界: Bekker の圧力-沈下モデルは本来、板の静的な押し込み試験のデータをフィッティングするために開発されたものです。車輪の動的な走行では、スリップに伴う追加沈下、土壌の粒子流動、車輪前方の土壌の盛り上がり（ブルドージング）など、静的モデルでは捉えきれない現象が生じます。

パラメータの場所依存性: 同じ惑星表面でも、場所によって土壌パラメータは大きく変動します。Spirit の Troy 事故が示したように、わずか数メートルの距離で土壌特性が劇的に変化することがあります。事前に全走行経路の土壌パラメータを正確に知ることは不可能であり、テラメカニクスの予測精度はパラメータの不確実性に強く依存します。

重力・大気の影響: 低重力環境では土壌粒子の挙動が変化し、地球上の実験で得たパラメータをそのまま適用できない可能性があります。大気がない月面では真空中での粒子挙動（静電付着、凝集など）も考慮が必要です。

近年の発展

テラメカニクスの研究は、古典モデルの限界を克服するために多くの方向に発展しています。

離散要素法（DEM: Discrete Element Method）: 土壌を個々の粒子の集合としてモデル化し、粒子間の接触力を直接シミュレーションする手法です。粒子の形状、摩擦、粘着力を個別に設定でき、Bekker モデルでは捉えられない微視的な現象（粒子の流動、チェーンフォースなど）を再現できます。計算コストが高いですが、GPU 計算の進展により大規模シミュレーションが可能になりつつあります。

有限要素法（FEM）との結合: 車輪を変形体として扱い（Curiosity の車輪はアルミニウム製で実際に変形する）、土壌との相互作用を FEM で解くアプローチです。車輪の変形が接触面積や圧力分布に影響を与える効果を考慮できます。

機械学習の活用: ローバーの走行データ（車輪電流、車輪速度、スリップ率、沈下量など）から、テラメカニクスパラメータをリアルタイムで推定する試みも進んでいます。事前に想定されたパラメータセットに依存するのではなく、走行中に得られるデータから土壌特性を学習し、走行戦略を適応的に変更するアプローチです。

Resistivity terrain model: NASA JPL が Curiosity の運用で開発した経験的モデルで、車輪モーター電流から直接走行抵抗を推定します。物理ベースのテラメカニクスモデルと経験的モデルを組み合わせることで、より信頼性の高い走行可能判定が可能になっています。

これらの発展は、テラメカニクスが単なる理論体系ではなく、現在も活発に研究が続いている実践的な工学分野であることを示しています。

まとめ

本記事では、テラメカニクス — 土壌と車輪の相互作用の力学 — について、基礎理論から Python シミュレーション、実際のミッション事例までを解説しました。

• Bekker の圧力-沈下モデル $p = (k_c/b + k_\phi)z^n$ は、車輪の沈下量を予測する基礎式です。土壌パラメータ $(k_c, k_\phi, n)$ と車輪の形状から、与えられた荷重に対する沈下量を定量的に求めることができます
• Janosi-Hanamoto のせん断モデル $\tau = (c + \sigma\tan\phi)(1 – e^{-j/K})$ は、車輪が地面から得る牽引力の源を記述します。Mohr-Coulomb の破壊基準に指数関数的なせん断変位依存性を加えたモデルです
• スリップ率と正味牽引力の関係には最適値が存在し、この最適スリップ率付近で走行させることがエネルギー効率と安全性の面で重要です
• 月面・火星レゴリスの特性は地球上の土壌とは大きく異なり、低重力と低粘着力の組み合わせが牽引力確保の課題を生みます
• Spirit のスタック事故は、テラメカニクスの予測能力をミッション運用に統合することの重要性を実証しました。低摩擦角、5輪走行、引きずり抵抗の複合要因が走行不能を引き起こしました
• 傾斜地での登坂限界は、垂直荷重の減少と重力抵抗の増加のダブルパンチにより、火星レゴリスでは概ね20度前後となります

テラメカニクスは、惑星探査ローバーの車輪設計、走行制御、ミッション計画のすべてにおいて基盤となる知識です。本記事で扱った古典的な Bekker / Janosi-Hanamoto モデルは、近年では DEM シミュレーションや機械学習との融合によってさらに発展を続けています。

次のステップとして、以下の記事も参考にしてください。

• 惑星探査ローバーの移動機構 — ローバーの車輪構造やロッカーボギー機構の基礎
• ローバーの自律ナビゲーション — SLAMと経路計画 — テラメカニクスと組み合わせた自律走行制御