引用元：土井健司、高木力、 光永靖、 鳥澤眞介 (2021) 数値流体力学的手法による平行泳法の魚の流体力学的効果. PLoS ONE 16(5): e0250837 https://doi.org/10.1371/journal.pone.0250837
( PLoS ONE は、Creative Commons, CC BY 4.0を受けています。）
編集者：ロイ・グルカ、コースタルカロライナ大学、アメリカ合衆国
受理：2020年9月1日 2021年4月14日受理 2021年5月3日公開
著作権： © 2021 Doi et al. この記事は、原著者と出典を明記することを条件に、あらゆる媒体での無制限の使用、配布、複製を許可するクリエイティブ・コモンズ表示ライセンスの条件の下で配布されるオープンアクセス論文です。
データの利用：すべての関連データは原稿の中にあります。
資金提供：日本学術振興会の科学研究費補助金（C）［助成番号20K04355］。
利害関係： 著者らは、利害関係がないことを宣言している。
魚類が群れを形成する理由は様々である。しかし、魚の群れのエネルギー効率を向上させる流体力学的なメカニズムは、まだ解明されていない。また、実際の遊泳運動に基づいたシミュレーションによる魚類モデルの例は少なく、既存のモデルでは動きが単純である。
そこで、本研究では、ビワマス（サケ科魚類）の遊泳挙動を画像解析により解析し、その遊泳運動を定式化した。さらに、定式化した遊泳運動を用いて数値流体力学解析を行い、実際の魚の遊泳運動で魚体モデルに作用する流体力を求めた。計算された周囲の流体力から平行遊泳時の魚体モデルの遊泳効率を求め、近接距離を変化させて比較した。
また、魚体モデル周囲の流れ場についても検討した。2匹の魚モデルが平行に泳ぐ場合の遊泳効率は、モデルの全長をLとしたとき、0.4Lの距離で離れると約10%向上することがわかった。さらに、魚体後方の流れ場を同位相と逆位相の条件下で、個体間距離0.8Lと1.2Lで調査した。2個体の鼻の中点から0.5-2.0Lの距離範囲における見かけの流速は，遊泳速度よりも小さかった。魚体模型の圧力分布では、尾びれで圧力が上昇していた。興味深いことに、尾柄に似た圧力分布マップが得られた。負の推力を避けるためには、今回得られた圧力分布マップのように、体の後部を細くする必要がある。
多くの魚は、捕食者の回避、高い摂餌効率、繁殖機会の向上などの理由から、群れを形成する。魚類の群れ行動や個体間の関係については，様々な研究が行われている。流体力学的な観点からは，群れで泳ぐと個体で泳ぐよりもエネルギー効率が向上する。
コンピュータ技術の発展に伴い、水生生物の遊泳メカニズムの研究に数値流体力学（CFD）を応用することが多くなっている。この方法は，空間を多数の小さな要素に分割し，各要素の基本流体方程式を数値的に解いて近似解を得るものである。
生きた魚の群れについては，ヨーロッパスズキ（Dicentrarchus labrax）は，単独遊泳時よりも群れでいる時の方がエネルギー消費量が少ない。しかし，生きた魚の実験では，魚体に作用する力を詳細に分離して測定することは不可能である。
本研究では、サンプル魚からサケ科魚類の遊泳運動を導出し、その運動を3次元矩形板魚体モデルに適用して、体側面における圧力分布を評価した。魚体モデルの翼断面には、NACA（National Advisory Committee for Aeronautics）4桁シリーズを使用した。これは上下対称の翼型であり，これまでの研究でも魚体モデルとして使用されている。ウナギは体全体が波打つように動くanguilliform，マグロは尾だけが振動するtunniform，サケ科は体の半分が波打つように動くsubcarangiformである。
本研究では，流体力学の観点から CFD 解析を行い，魚体周りの流れ場と遊泳効率を評価することで，魚群における並列泳の効果を明らかにすることを目的とした。平石らは，ニジマス（Oncorhynchus mykiss，サケ科）の遊泳運動を解析し，数式で表現している。Videler and Hessは，セイス（Pollachius virens）とサバ（Scomber scombrus）の遊泳運動について検討し，その結果を発表している。我々は、サケ科魚類の遊泳運動の画像を解析した。泳ぎを定式化し，CFD 解析により魚体モデルに作用する流体力 を算出した。体側波動関数を CFD に適用した例はいくつかあるが，その横方向の振れ幅は単純である。使用した体表波動関数は，尾端に近づくにつれて単純包絡線が大きくなるような振幅を有していた。平行移動が遊泳効率に及ぼす影響を確認するため，魚体モデルに作用する抗力と推力を抽出して遊泳効率を計算し，近接距離を変えて比較するとともに，魚体周りの流れ場についても検討した。
この研究は、近畿大学動物研究・生命倫理委員会（許可番号：KAAG-30-001）の特別な承認を受け、「研究動物の世話と使用に関する指導原則」に厳密に則って実施されたものです。
遊泳運動のシミュレーション
滋賀県の琵琶湖でビワマス（Oncorhynchus sp.）を捕獲し、1匹の試験魚として研究に使用した（表1）。試験魚は、循環式水槽（PT-110（改良型）、容積：0.396 m3、観察部：1箇所）で1点泳がせた。幅 1.1m，深さ 0.3m，高さ 0.2m 株式会社西日本流体工学研究所）内で泳がせ，水槽上部からデジタルビデオカメラ（HDC-TM750，パナソニック株式会社）により 60fps で記録した。
表1. 試験魚のパラメータ
撮影したビデオデータをもとに、魚体上の鼻先から尾鰭の先端までの13点の時系列2次元位置座標を求めた。座標の抽出には、研究室で開発したソフトウェアを使用した。体軸上の点は、映像から位置が特定できる魚体上の点から設定した。各点は、鼻先（a）、目（b、c）、胸鰭基部（d、e）、胸鰭後端（f、g）、背鰭付け根（h）、背鰭後端（i）、脂鰭付け根（j）、尾鰭付け根（k）、尾鰭後端（l、m）に割り振られた。さらに、b-c, d-e, f-g, l-mの中点を体軸上の点、すなわちそれぞれn, o, p, qとみなした（図1）。尾びれ後端については、尾びれの上端と下端の中点を尾びれ先端の位置とし、尾びれ後端は、尾びれの上端と下端の中点を尾びれ先端の位置とした。
図1. 体軸上の点
体軸上の各点の2次元位置座標から体軸に沿った振幅を計算し、体軸上の進行波の周波数と波長を求めた。
ここで、h（x，t）は魚体断面の進行方向に垂直な方向の変位、x（BL）は魚の進行方向と平行で反対側にある体長上の鼻からの距離、t（s）は時間、fは周波数、λ（BL）は進行波の波長である。また、a、b、c、dは遊泳運動の形状を表すパラメータである。体軸上の位置xと振幅の関係は、aebxとcedxという2つの指数関数の和で表すことができる。遊泳の推力は主にcedxで表される部分である尾の振動によって発生し、尾は頭部付近の振幅の影響をわずかに受ける一方、頭部付近の振幅は尾の振幅の影響を受けるため、cedxのxの値から頭部付近の値を代入して得られる値を引いたものがaebxである。式2の未知パラメータcとdを求めるため、x＞0.4の振幅値を対数変換し、最小二乗法で線形近似した。x < 0.25の部分では，x < 0.25の点でのcedxの値を求め，その値を振幅値から引き，同様に近似することにより，式2の未知定数a，bを推定した。
ソリタリースイミングのCFD解析
CFD 解析に用いたモデルは，左右対称の翼を持つ長方形 の板状モデルである。魚類モデルの長さ（L）、高さ、最大翼厚は、ビワマスの実測データと等しくなるように、NACA4桁シリーズの翼（図2）をベースに、それぞれ43.9、8.93、4.51cmとした。モデル作成には3次元コンピュータ支援設計（CAD）ソフト（AutoCAD 2017、Autodesk）を使用した。
図2. CFD解析に使用したモデル
これらの方程式の解は、連続式（式4、式5）を満たすことで得ることができる。式中、Uは速度ベクトル、tは時間、Fは単位体積あたりの力ベクトル、ρは流体密度、pは圧力、μは粘性係数である。
本研究では、3次元熱流体解析システム（SCRYU/Tetra V14, Software CRADLE CO., LTD）を用いて、Navier-Stokes方程式と連続方程式を有限体積法で数値的に解き、物理量を算出した。
魚類モデルにおける要素の移動条件として、実際のビワマスの遊泳運動から導き出された遊泳運動関数を使用した。要素移動条件とは、移動する物体や境界形状が変化する物体の流体解析を行う際に、任意ラグランジュオイラー（ALE）法を用いて、要素を構成する接点を体積領域に対して移動させる条件である。本解析では，導出した遊泳運動関数を用いて，遊泳運動中の魚体モデル の各要素の節点を移動させ，CFD 解析を実施した。
独泳のCFD解析では，以下の条件を設定した。解析領域は，魚体モデルの左右の形状変化，背水の 視点，壁の影響を考慮して，軸方向 5.0L，左右方向 4.0L，上下方向 2.0L とした。Lは魚体模型の全長になる。
魚類模型の前面から5段階の流れの定常流を流入し、後面から圧力と流れ方向を一定にして流出させた。流入速度は生魚実験で使用した速度と同じで、生魚実験のそれぞれで記録された各速度での遊泳運動から遊泳運動関数を算出した（図3）。
図3. ソリタリースイミングの CFD 解析のための解析領域。(a) 上面図、(b) 側面図。
壁面条件は，解析領域の流入面，流出面を除く全ての面において，フリースリップで摩擦応力の影響を受けないものと仮定した。流体温度は 20℃，流体密度は 998kgm-3，粘性係数は 1.02×10-3 (Pa・s)，熱伝導率は 0 と仮定した。3 次元魚体モデルのレイノルズ数およびストローハル数は，それぞれ 2.60×105 および 0.34 とした。乱流の影響を考慮するために，低レイノルズ数から高レイノルズ数までの非定常乱流の解析に適した AKN k-ε モデルを使用した。
CFD 解析では、解析領域を小さな要素に分割し、それぞれ の流体方程式を解いて解を求める。図4に今回の CFD 解析におけるメッシュ分割を示す。テトラ要素サイズは 2～64mm の範囲とし，魚体モ デルに近づくにつれて細かくなるようにテトラ要素サイズが 2mm となるようにモデル最近傍要素のサイズを分割している。解析領域の総要素数は 1.20×106 であった。実際の流れでは、壁からの摩擦応力により流体と壁の間に流速勾配が生じ、境界層と呼ばれる薄い流れの層が形成される。この境界層における流速の勾配を表す境界層要素（プリズム層）と呼ばれる魚の表面に平行な層を挿入することで、乱流解析の精度を高めることができる。境界層要素は、魚のモデルの表面を3層で覆うように、変化率1で挿入し、1層目のプリズム層の厚さは0.4mmとした。
図4. 孤立遊泳の CFD 解析におけるメッシュ分割
CFD 解析では，モデル表面上の各要素に作用する流体力 を任意の成分に変換することができる。そのため，魚体モデル表面に作用する流体力のうち，前方 向の成分をスラスト力，逆方向の成分を泳ぐ魚体モデルに作用 する抗力として計算することが可能であり，その結果，魚体モデ ルに作用する流体力を算出することができる。本解析では、魚体モデルに作用する体軸方向の流体力のうち、圧力のスラスト方向成分をスラスト力と定義した。抗力は，圧力の抗力方向成分と摩擦抗力の和と定義した。
パラレルスイミングのCFD解析
解析エリア内に2匹の並列魚モデルを設置し、並行して回遊する際の背水の変化を観察した。2個体の運動時のメッシュの変形による計算の発散を避けるため、重合格子法を採用した。重合格子法では、独立領域と呼ばれるメッシュに従属領域と呼ばれるメッシュを重ね合わせることで、一つの解析空間を表現する。この方法の利点は、領域ごとに運動条件を設定し、複雑な物体の移動や衝突を解析することができる点である。
パラレルスイミングに使用したモデルは、ソロスイミングに使用したモデルと同じ左右対称の翼の断面形状を持つ長方形のプレートモデルである。遊泳運動は、単独遊泳に適用したものと同じ遊泳運動関数を適用した。
CFD解析の条件は以下の通りである。解析領域は，魚体モデルの左右の動きを考慮した大きさ とした（図5）。解析に用いた流入速度は，単独泳法の CFD 解析で求めたものである。このとき，各解析条件で推力と抗力がほぼ均衡するように尾びれ振動の周波数を変化させた。壁面条件は，解析領域内の流入面，流出面を除く全ての面において，フリースリップで摩擦応力の影響を受けないものとした。流体温度は 20℃，流体密度は 998kg/m3，粘性係数は 1.02×10-3（Pa・s），熱伝導率は 0 とし た。3 次元魚体モデルのレイノルズ数は 2.60×105 であった。3次元魚体モデルのストロハル数は0.42〜0.44であった。乱流解析モデルには，低レイノルズ数から高レイノルズ数までの非定常乱流の解析に適した AKN k-ε モデルを使用した。
図5.パラレルスイミングのCFD解析のための解析領域。(a) 平面図，(b) 側面図。
近接から単独遊泳まで同じ解析条件で並行遊泳モデルの個体間距離を比較するため、距離を0.4, 0.6, 0.8, 1.2, 1.6, 2.0L とし、振幅運動の位相が同相と逆相の場合をそれぞれ解析した。入口流速は1.52BL s-1とした。
図6に本解析におけるメッシュ分割を示す。テトラ要素のサイズは2～32mmの範囲とし、魚体モデルに近づいたときに細かくなるようにテトラ要素のサイズが4mmとなるように最も近いモデルの要素のサイズを分割した。境界層要素は、魚体モデルの表面を覆うように8層、変化率1.1で挿入し、第1プリズム層の厚さは0.8mmとした。
図6. パラレルスイミングのCFD解析におけるメッシュ分割
遊泳の効率
遊泳時の効率を比較するための指標を算出することを試みた。魚類の遊泳効率を評価するために，しばしばFroude効率と呼ばれる指標が用いられる。Hemelrijkらは，Froude効率から群泳の流体力学的な優位性を報告している。
モデル表面に作用する力を構成要素に分解することで、推力と横力を算出することができる。横力のパワーは、尾振り運動の速度の横成分と、魚体モデルの表面に作用する圧力の横成分の積で求めることができる。横方向成分とは、泳ぐ方向（Y方向）に垂直な方向を表す。本研究では、魚体モデルの表面要素に着目した。表面要素に作用する流体力から総力に対する有用力の比を算出した。また、異なる条件を比較するための指標として、遊泳効率を用いた。
図7. 魚体表面の要素にかかる力Fiの分解
Uは流入流速、V0は要素の速度、uiは要素の進行方向速度、viは要素の横方向速度、FTiは要素に作用する推力、FSiは要素に作用する横力である。
Borazjani と Sotiropoulosが示唆したように，推力と抗力の値が釣り合った後に，遊泳効率を評価した。
遊泳運動の定式化
式2の未知パラメータa、b、c、dは、体軸上の各点における振幅の値を用いて指数近似により求めた。表2に式2の未知定数と、各流速における遊泳運動関数を求めるために用いた式3のλ、fの値を示す。なお、遊泳速度のBLは、供試魚の体長を基準とした。魚の体軸上の各点の時系列データを、実測データと、遊泳速度1.52 BL s-1で遊泳運動関数から求めたデータとを比較したところ、図8（R²＝0.96）のようになった。
図8. 魚の体軸上の各点におけるデータの時系列
遊泳速度1.52BL s-1における実測データ（実線）と遊泳運動関数から求めた結果（破線）を比較したものである。
表2. 式2の定数a, b, c, dと式3の尾びれの波長(λ)と周波数（f）
ソリッドスイミングのCFD解析
図9は、遊泳速度1.52BL s-1で4ストローク尾打ちを行った際に、魚体モデルの表面に作用するスラスト力と抗力を示している。推力と抗力の変化には周期性があることが分かる。また、尾びれの先端で振幅が最大となる点と、スラスト力が最大となる点の間に位相差があることがわかる。
図9. 魚のモデル表面に作用する推力と抗力
これらの力は、遊泳速度1.52BL s-1で尾びれを3ストローク分行った際に作用したものである。
推力と抗力の値を比較するために、推力と抗力の時間平均値をそれぞれFT、FDと定義し、各泳速度で比較した（表3）。今回の解析で用いた遊泳運動関数は、すべて循環水槽内を一定速度で定常的に泳ぐ生物実験により決定したものであり、この状態で推力と抗力が釣り合うようになっている。CFD 計算で解析した 5 つの遊泳速度において，FT/FD 比は遊泳速度 1.52 BL s-1 で最も 1 に近い値となった（表 3）．図10は、遊泳速度1.52BL s-1で尾びれの尾端が最大速度になったときのモデル表面の圧力分布である。
図10. モデルの面圧分布
遊泳速度1.52BL s-1における最大速度（a）時の尾びれ尾端部の右側（1）と左側（2）を示している
表3. 推力と抗力（それぞれFTとFDとする）の時間平均値を各水泳速度で比較した
パラレルスイミングのCFD解析
図11は、同位相および逆位相の場合において、個々の距離が0.4Lと0.8Lの平行移動中の魚の周りの流れを示している。
図11.モデル周辺の速度分布：(a) 0.4Lと同位相の個別距離、(b) 0.4Lと逆位相の個別距離、(c) 0.8Lと同位相の個別距離、 (d) 0.8Lと逆位相の個別距離、 (e) 1.2Lと同位相の個別距離、(f)1.2Lと逆位相の個別距離。
図12では、横軸は、個体間距離0.8L、逆位相下での2個体の鼻の中点を原点とし、x軸に平行な点線上の距離を表している。図11、図12に示すように、同相、逆相のいずれの条件においても、個体間距離0.4Lでは、2個体の鼻の中点から0.5～2.0Lの距離範囲における見かけの流速は入口流速より大きい。一方、個体間距離0.8Lと1.2Lの両フェーズでは、前述の距離範囲での見かけの流速は流入流速より低くなる。
図12. 速度分布の各点における速度値
横軸は個体間距離0.8L、逆位相の条件で2個体の鼻の中点を原点としてx軸に平行な点線上の距離を表している。縦軸は点線上の速度値を表す。長破線は入口流速(1.52 BL s-1)を示す。
同位相と逆位相の条件下で分離した2個体の遊泳効率を比較した（図13）。同位相条件下では大きな変化は見られない。しかし、逆位相の条件下（破線）、距離0.4Lでは、遊泳効率は約10%向上した（図13）。
図13.同位相（実線）と逆位相（破線）の条件下で2個体が並行して泳ぐ場合の遊泳効率
図10に示した魚モデルの圧力分布では、平面2の尾びれの圧力が上昇している。興味深いことに、平面1では、尾柄と同様の等値線マップが観察される。Wuは，尾びれは力を発生させるために推進するのであり，負の推力を発生させないためには，体尾部を細くする必要があると報告している。
魚体後流では、図12に示すように、個体間距離が0.8Lのときに低速度分布域が発生する。最大限の節約を実現するためには，魚は逆位相の尾部運動で連続的に泳ぐことが望ましいと報告されている。このため、2匹の魚に追従する個体では、航跡上の見かけの流速が実際の遊泳速度より低くなる可能性がある。このことは、後続の魚がある距離を移動した後に、エネルギー節約を実現することを意味する。
抗力は遊泳速度の2乗に比例し、パワーは遊泳速度の3乗に比例する。個体間の横方向の距離が0.8Lで逆位相の場合、2個体間の正中線上の鼻から1〜2Lの距離の流速は、尾びれの1周期で約0.969倍となる。したがって、遊泳に必要な動力は約0.910倍となり、省エネ効果があると推定される。魚体後方の流速は0.974倍で、同相状態では鼻からの距離は1〜2Lの範囲にある。
魚の鼻から1〜2Lの距離に個体がいて、逆位相から同位相に変わったとしても、魚はある程度得をすることになる。一方、反位相の状態でも鼻先から2〜3Lの距離では見かけの流速は1.03倍であった。群れとしては、1-2Lの距離範囲に後方遊泳の個体がいたほうが、魚のエネルギーを節約できる。個体間距離が1.2Lの場合、実際の後方遊泳速度は約0.98倍となる。この場合、個体間距離が0.8Lのときほどの省エネ効果は得られない。このように、個体間距離が0.8Lの場合、後続の遊泳者はエネルギーを得ていると考えられる。個体間距離が0.4Lに近づくと、見かけの流速が単独泳者より大きくなり、負荷が大きくなる可能性がある。したがって、省エネにはならない。
Weihs は、魚群の菱形配列が理論上後続個体の省エネ効果を高めるのに有効であると報告しているが、本研究での横泳ぎ下の値は節約率約9％であり、Weihsの報告ほど高い利得は得られていない。これは、Weihsがポテンシャル理論を適用して渦列の単純な重ね合わせから計算を行ったため、実際の粘性流体の効果が十分に反映されていないためと思われる。また，3次元物体モデルの形状をコピーしていないため，その形状に起因している可能性もある。Vermaらは、一直線に泳ぐのではなく、中心から外れたところを泳ぐと効果的であると報告している。
遊泳効率は、個体間距離0.4L、逆位相の状態で最も高くなった。これは、図.14に示すように、尾びれの振幅運動の最大速度で表面圧力の等値線が描かれるため、推力が高くなったものと思われる。Baoらは、逆位相で運動する2枚のフラッタリングフォイルが後流の干渉によりベナール・フォン・カールマン渦の通りを誘発し、推進性能を向上させることを報告している。同様の現象が本研究の結果に影響を与えた可能性がある。また、個体によって尾鰭振動の位相が常に保たれているわけではない。また、尾鰭振動の位相は常に一定ではないため、高い利得が得られる期間と得られない期間が繰り返される可能性がある。しかし、トータルバランスで考えると、遊泳中はプラスのエネルギーゲインが得られることになる。
図14.尾びれの振幅運動の最大速度における表面圧力の等高線プロット：(a) 個体間距離0.4L、(b) 個体間距離2.0Lの場合。
生きた魚の実験では，高速遊泳時に群れが凝集し，他の個体の動きと同期する傾向が見られた。この報告は，より横方向にコンパクトな魚の群れがより多くのエネルギーを獲得できるという我々の結果を支持すると考えられる。
魚群の中で泳ぐことだけが省エネではない。バースト＆グライド泳法と各個体の省エネを組み合わせることで、全体の省エネをさらに高めることができる。
また、本研究では、個体の加速度、減速度を考慮していない。今後，CFD 解析により，個体の加減速を考慮した効率計算を行う必要がある。
魚群の3次元構造をさらに解析することで、魚の群れにおける個体間の流体力学的な相互作用をより正確にシミュレーションすることができるようになる。
