このシミュレーターが計算していること、そして合わなかったこと

ルアーの動きシミュレーターの中身です。

結論を先に書きます。振動周波数は実測に近い値が出ました。振れ角は、実測の1桁下しか出ませんでした。この記事の後半は、なぜ足りないのかの話です。

計算したのは私です。論文の再現ではありません。式の組み立ても、合わなかったことの説明も、私が考えたものです。

何を解いているか

上から見た平面の中だけを扱います。

ルアーを剛体とみなし、自由度を2つに取ります。

y — ルアーの横位置
θ — ヨー角。頭が左右どちらを向いているか

鉛直方向、つまり潜っていく向きの動きは解いていません。理由は後述します。

流体力の出し方

ルアーを長さ方向に30個の区画に切り、区画ごとに流体力を計算して足します。細長い物体を扱うときの標準的なやり方です。

各区画は、その位置での平面図上の幅だけの面積を持ちます。リップは長方形、ボディは楕円としています。

区画 i が受ける横向きの力は、次のようになります。

F_i = q × 面積_i × Cn(α_i)
q   = ½ × 水の密度 × 巻き速度²

Cn は法線力係数で、迎角に対して Cn(α)= 1.8 × sin(2α) としました。45度で最大になり、90度で0に戻る形です。臺田ら2002が測った潜行力係数が40〜45度で最大になっていることに合わせてあります。

α_i は区画ごとの局所的な迎角です。ここが要点で、ルアーが回転していると、区画ごとに迎角が違います。

α_i = θ − (ẏ + ζ_i × θ̇)/ 巻き速度

ζ_i はラインアイから測った区画の位置で、頭に向かう向きを正としています。頭側の区画と尾側の区画では、回転による迎角の変化が逆向きになります。これが姿勢を元に戻す働きを生みます。

ラインアイが基準点になっています

モーメントはラインアイのまわりで計算しています。臺田ら1999が「ルアーはタイドアイの位置を張力の作用点として潜行する」と書いているのに合わせました。ラインの張力はこの点に働くので、この点まわりのモーメントには現れません。

したがって、ラインアイの位置を動かすと ζ_i が全部ずれ、姿勢を戻す力の強さが変わります。スライダーで Ct を動かしたときに挙動が変わるのは、これが理由です。

ラインそのものは、横方向のばねとして扱っています。ばね定数は抗力を糸の長さで割った値です。

時間を進める方法:ルンゲ・クッタ法

運動方程式が出す答えは「いまの状態から、次の一瞬でどれだけ変化するか」です。位置と速度が分かれば加速度が計算できる。これを少しずつ積み上げて、時間を先へ進めます。

いちばん素朴なやり方は、いまの加速度がしばらく続くと思って、そのぶんだけ進めるというものです。オイラー法といいます。

次の速度 = いまの速度 + いまの加速度 × 刻み幅
次の位置 = いまの位置 + いまの速度   × 刻み幅

簡単ですが、曲がる運動に弱いのが難点です。円を描くように動くものを計算すると、外へ外へとずれていきます。加速度は刻みの間にも変わっているのに、入口の値で代表させているからです。

ルンゲ・クッタ法は、刻みの中で加速度を何度か測り直します。4次のものは、1回の刻みで4回測ります。

区間の入口で加速度を測る
その値で区間の真ん中まで仮に進めて、そこで測り直す
2で得た値を使って、もう一度区間の真ん中へ進めて測り直す
3で得た値で区間の出口まで進めて、そこで測る

そして4つを重みつきで平均します。真ん中で測った2つを重く見ます。

変化量 = 刻み幅 ÷ 6 × (1回目 + 2×2回目 + 2×3回目 + 4回目)

坂道を降りるときに、出発点の傾きだけで一歩の大きさを決めるのではなく、途中の傾きも確かめてから一歩を決める、という違いです。計算量は4倍になりますが、刻み幅を粗くしても同じ精度が出るので、結局は得になります。

「4次」というのは、刻み幅を半分にすると誤差が16分の1になるという意味です。オイラー法は、半分にしても誤差は半分にしかなりません。

なぜこの方法にしたか

このサイトの飛跡シミュレーターでは、別の方法(速度ベルレ法)を使いました。あちらは空気抵抗を除けばエネルギーが保存する運動で、ベルレ法にはエネルギーを長時間ほとんど失わないという性質があります。何秒計算しても放物線が崩れないのはそのためです。

今回はエネルギーが保存しません。渦がエネルギーを入れ、流体力がエネルギーを奪う。入る量と出る量がつり合ったところで振幅が決まる、という運動です。ベルレ法の利点が効かないので、単純に精度の高い方法を選びました。

刻み幅は0.35ミリ秒です。4.5秒ぶんを計算して、最後の1.2秒だけを使って周波数と振れ角を出しています。始めの数秒は振れ方が落ち着くまでの過渡的な状態なので、捨てています。

刻み幅を細かくしても結果が変わらないことは確かめました。0.25ミリ秒で計算し直すと、周波数のずれは0.2%以内、振れ角のずれは1.1%以内でした。

準定常の流体力だけでは、振動が止まります

ここが作っていていちばん意外だったところです。

上の式だけで解くと、振動は減衰して止まります。ルアーは少し傾いた姿勢で落ち着き、そのまま真っ直ぐ進みます。振れません。

理由は式を見れば分かります。回転による迎角の項 ζ_i × θ̇ / U は、頭側と尾側の両方で回転を妨げる向きに働きます。つまり減衰項です。流体力は、姿勢を戻すだけでなく、揺れを止める働きもします。

実物が振れ続けているということは、揺らし続けている何かがあるということです。

論文はそれをカルマン渦だと説明しています。流れの中の物体の後方では両側から交互に渦が離れ、そのたびに横向きの力が周期的に変わります。三木ら2001も臺田ら2001も、この説明を採っています。

そこで、渦の力を別の式で足しました。後流振動子と呼ばれる、渦の放出をひとつの振動子で表すモデルです。ファン・デル・ポール型の非線形振動子を1つ用意し、ルアーの横加速度と結合させます。

Q̈ + ε ω_s (Q² − 1)Q̇ + ω_s² Q = A × (ボディの横加速度)/ 体高
渦の力 = q × ボディの平面面積 × CL × Q
ω_s = 2π × ストローハル数 × 巻き速度 ÷ 体高

渦が離れる周波数は、巻き速度に比例し、体高に反比例します。ストローハル数は0.25としました。

この結合があると、渦の放出とルアーの揺れが引き込み合います。実物のルアーで起きていることも、おそらくこれです。

迎角は、計算をあきらめました

鉛直方向、つまり深さの向きを解いていない理由がここです。

はじめは、鉛直面内でのモーメントの釣り合いから迎角を出そうとしました。ルアーを横から見た図の中で、流体力と浮力がつり合う姿勢を探す、という計算です。ラインの張力はラインアイに働くのでモーメントを持たず、残るのは流体力と浮力だけになります。これなら解けるはずでした。

出てきた迎角は、1度前後でした。実測は6〜57度です。桁が違います。

原因は分かります。浮力が小さすぎるのです。比重0.94のルアーだと、浮力から重さを引いた残りは0.01ニュートン程度しかありません。同じ条件で流体力は0.7ニュートン近くあります。浮力では、流体力に釣り合う姿勢を作れません。

つまり、実物の迎角を決めているのは、私が式に入れたどれでもありません。リップの取り付け角、ラインアイのボディ軸からの上下のずれ、ボディの上下非対称——どれが効いているのかを特定できませんでした。

そこで、迎角は論文の記述から与えることにしました。

臺田ら2001: 流速40〜80 cm/s において迎角は6度から57度の範囲。Ct とともに増大する
三木ら2001: Ct = 0.25 で迎角は35〜40度。流速が上がると小さくなる

この2つを通る単調な曲線を引いて、Ct と巻き速度から迎角を返す関数にしました。計算ではなく、実測の読み替えです。シミュレーターの表示にも「実測から与えた値」と書いてあります。

そして、この迎角が渦の強さを決めます。迎角が大きいほど流れが大きく剥がれ、渦が強くなると考えて、渦の力を迎角の正弦に比例させました。ラインアイを前に付けすぎると振れなくなるのは、この経路です。

合ったこと

周波数は巻き速度に比例します

巻き速度	計算した振動周波数	実測 (臺田ら2001、図からの読み取り)
40 cm/s	4.1 Hz	5〜6 Hz
60 cm/s	6.2 Hz	7〜9 Hz
80 cm/s	8.3 Hz	10〜15 Hz

比例関係が出ました。三浦ら1985が「流速と振動数の間には非常に高い直線性が認められる」と報告しているとおりです。

絶対値は2〜4割低く出ています。ストローハル数を0.3程度に上げれば合いますが、合うように選んだ数字になってしまうので、標準的な0.25のままにしてあります。

体高を高くすると、周波数が下がります

体高比 κ (= 体高 ÷ ボディ長)	計算した振動周波数 (巻き速度80 cm/s)
0.25	11.4 Hz
0.41	7.0 Hz
0.62	4.5 Hz
1.00	2.8 Hz

三木ら2001は、κ = 0.25 付近で周波数が最大となり、κ の増大に伴って小さくなったと報告しています。向きは合っています。

ただし、実測の変化はもっと緩やかです。図から読み取ると、実測は約22 Hz から約11 Hz へ、半分になる程度です。計算は4分の1になっています。体高への依存が急すぎます。

ラインアイを後ろへ移すと、振れが大きくなります

ラインアイの位置 Ct	計算した振れ角	計算した迎角
0.10	0.3°	10°
0.15	0.7°	19°
0.20	1.0°	27°
0.25	1.4°	34°
0.30	1.8°	41°
0.35	2.1°	48°

臺田ら2001は Ct = 0.15 と 0.18 では振動が認められなかったと報告しています。向きは合っています。

合わなかったこと

振れ角が1桁足りません

これがいちばん大きな食い違いです。

体高比 κ	計算した振れ角	実測 (三木ら2001、図からの読み取り)
0.25	1.0°	約0°
0.41	1.6°	約11°
0.62	2.3°	約20°
1.00	2.5°	約35°

体高が小さい範囲では近い値です。κ = 0.25 では実測がほぼ0度、計算が1度。同じくらいです。

体高が大きくなると、まったく足りません。κ = 1.00 で実測35度に対して計算2.5度。14分の1です。

理由は、モデルの構造にあると思います。ここからは私の推測です。

計算しているのは、渦の力に対する強制的な応答です。渦が押す、ルアーが少し振れる、という関係です。渦の放出周波数は、ルアー自身が姿勢を戻そうとする固有の周波数より高いところにあり、共振していません。共振していない領域では、力を加えても振幅は大きくなりません。

実物で起きているのは、おそらくこれとは別のことです。大きく振れるルアーでは、ルアーの揺れそのものが渦の放出を引きずり、両者が完全に同期していると考えられます。この状態になると振幅は一気に大きくなります。後流振動子はこの状態を表せるモデルなのですが、私の設定では引き込みが起きるところまで行きませんでした。

もうひとつ、三木ら2001が挙げている説明も効いていそうです。体高が大きいと、流れに対してボディが下流方向に長くなり、渦がボディの側面付近で形成されるようになる。渦ができる場所そのものが移ります。区画ごとに面積をかけるだけのモデルでは、渦の位置は表せません。

潜行抗力比が、実測の範囲を超えます

シミュレーターは潜行抗力比も表示していますが、ラインアイを前寄りにすると3.4まで上がります。臺田ら2002の実測は最大2.54です。

抗力係数の式を 0.10 + 1.1 sin²α という簡単な形で置いたためで、迎角が小さい領域で抗力を低く見積もりすぎています。ここは合わせ込んでいません。

潜行深度は出していません

計算すればそれらしい数字は出せますが、ラインの抗力とたわみを入れていないので、絶対値がまったく合いません。試したところ、20mのコースで10メートル以上潜るという結果になりました。臺田ら1999の実測は2〜4メートルです。

出さないことにしました。かわりに、臺田ら2002と直接比べられる潜行抗力比を表示しています。

このモデルで言えないこと

シンキングペンシルには使えません。リップがないルアーは、そもそも定常的な姿勢を取りません。後編の記事に書いたとおり、回流水槽の計測でも同じ理由で数値が取れていません。

フックが入っていません。三木ら2001が課題に挙げているとおり、フックは重心と浮心の位置を変え、姿勢と振動に影響するはずです。

渦を解いていません。渦の放出は1つの振動子で代表させているだけで、流れ場は計算していません。東京電機大学の講究録が2015年に「張力との釣り合いも考慮したうえで周期運動について考察する必要がある」と未解決のまま結んでいる問題です。私が片手間に解けるものではありません。

振れ角が水中でどれだけの水の動きを生むかも、計算していません。これはベイトを感じ取る仕組みの記事で扱った側線の話につながりますが、そもそも測ったデータが見つかっていません。

そして、どう泳ぐと釣れるかは、このモデルからは何も出ません。

続けて読むなら

ウォブリングとローリングは、別の動きではありませんでした— このモデルを作る前に調べた、呼び分けの話
リップの大きさの影響とルアーを比べる物差し— リップの形が何を変えるかと、比べるための物差し
ウォブリングはスレやすいのか— 動きの違いが釣果に効くのかどうか
大きいルアーは大きい魚を釣るのか— 動きではなくサイズを変えたときに何が起きるか
「良く釣れるルアー」の正体— 振れ角の差が釣果に出るかどうかは、個人の釣行では確かめられません