Juliaでサメを泳がせます。その2
Reference
- Gabriel Weymouth, "Simulation of a swimming dogfish shark," The Julia Language Blog, 2021.
Disclaimer
内容の正確性には可能な限り留意しておりますが、必ずしも正確性を保証するものではありません。 内容に基づく一切の結果について、一切の責任を負いかねます。 予めご了承ください。
はじめに
前回記事ではJuliaCon 2021の講演アーカイブを参考にして、流体シミュレーションのパッケージであるWaterlily.jlを動かしてみるところまでをやった。
このパッケージを使ってサメの泳ぎのシミュレーションができるよ、ということがThe Julia Language Blogの以下の記事にある。記事を書いているのは前述のパッケージの作者である。
この記事は前編に続く中編的な位置づけ、サメを泳がせる前にPluto.jlを練習しておくのと、前回の流体シミュレーションのフォローアップをしておきたい。
WaterLily on Pluto
前述のブログ記事であるSimulation of a swimming dogfish sharkについて、この記事の良いところはPluto.jlというノートブックの仕組みを用いて、記事がそのまま実行できるコードの集合になっていること。 Pluto.jlとは公式の説明によればreactive, lightweight, and simpleを特長としたノートブック。 見た目はJupyter Notebookと似ているものの、とくに対話操作に長けている様子。 Qiitaの紹介記事によれば大学の講義、教育用途での利用もされているようだ。
ブログ記事冒頭にあるPlutoノートブックを開くには、あらかじめJuliaにPlutoパッケージを追加しておく必要がある。久しぶりすぎてJuliaでのパッケージ管理の方法を忘れてたのでおさらいしておく。
- パッケージ管理の全体像は『1から始める Juliaプログラミング』の2.11.1項に詳しい。
]で入れるパッケージ管理モードからadd Plutoを実行してパッケージを追加した。
プロンプトに戻り、using Pluto; Pluto.run()を実行するとブラウザでノートブックが立ち上がる。
用意されたサンプルのノートブックを開くか、新しいノートブックを作るか、もしくはファイルパスを指定してノートブックのファイルを開く。
3つ目のやり方でうまくファイルを指定をするためには、開きたいファイル、今回の例でいえばブログ記事からダウンロードできるPlutoノートブックのファイルがある階層でJuliaを起動しつつusing Pluto; Pluto.run()を実行することで、ファイル名のサジェストが効いてくれるはず。
あれこれ試行錯誤して最終的に次のようなコードを書いた：
n,m = 3*2^6,2^7
center,radius = m/2,m/6
body = AutoBody((x,t)->√sum(abs2,x .- center) - radius)
sim = Simulation((n+2,m+2), [1.,0.], radius; ν=ν, body=body)
sim_step!(sim,10)
# plot the vorcity ω=curl(u) scaled by the body length L and flow speed U
function plot_vorticity(sim)
@inside sim.flow.σ[I] = WaterLily.curl(3,I,sim.flow.u)*sim.L/sim.U
contourf(sim.flow.σ',
color=palette(:BuGn), clims=(-10,10),linewidth=0,
aspect_ratio=:equal,legend=false,border=:none)
end
# make a gif over a swimming cycle
@gif for t ∈ sim_time(sim) .+ range(0.0, stop=20.0, step=0.25)
sim_step!(sim,t)
plot_vorticity(sim)
end
上記コードを3番目のセルに格納した様子。一つのセルに複数行を入れるときにはbegin~endで囲む。
関数plot_vorticityはブログ記事にあるものを流用した。
この関数を使ってgifアニメーションのもととなる渦度分布を作っている（はず）。
セルの右側にあるプルアップメニュー📖はLive docsといって、ここに変数名を入力することでインタラクティブに説明が得られる。
上の図では試しに複合型 (composite types) simを表示させているけどとても便利。
Plutoのインタラクティビティとはあるセルに加えられた変更を他のセルが自動的に検知してくれるということ。
例えば上の図で2番目のセルにある動粘性係数ν = 0.0の値を変更すると3番目のセルが自動で実行されて、対応した流れ場を得ることができる。
Jupyter Notebookではこうした反映は手動で実行していたはずなので、それと比べると楽だし安全になっていてよい。
これは逆にいえば、2番目のセルで動粘性係数の値を変更すると3番目のセルは必ず自動的に実行されてしまう。 特に今回のケースでは3番目のセルの実行には多少時間がかかるので、変数の値を変えるたびに実行してほしくない状況もままあるかもしれない。 そうした場合にはセル右側にあるメニューから "Disable cell" を選択し、セルを無効化することで自動実行を止めることができる。 動粘性係数の値を変えたのちに改めて "Enable cell" とすることで、所望の反映を得ることができる。
なおこの "Disable cell" 機能は最近のアップデートらしく、Pluto.jlのJuliaCon 2021のアーカイブ動画で新機能 (new feature) として紹介されている。
円柱周りの流れ
ところで前回記事では動粘性係数の値を適当に変えることで3種類のそれっぽい流れ場を作っていた。
ソースコードに残されたコメントによれば変数νとはScaled viscosity (Re=UL/ν) であって、たとえばパッケージ内のソースファイルFlow.jlにおいてν*∂(j,CI(I,i),u)みたいに使われている。
簡単な計算により（要出典）予想されるのは、前回記事で用いたν=1.0というのは一様流速度と円柱直径に基づくレイノルズ数が20くらいに相当するということ。
そしてν=0.2とはこの5倍、レイノルズ数が100くらいに対応していると思われる。
前回記事でも引用した円柱周りの流れの論文によれば、レイノルズ数RがR<60～70のときには双子渦領域、60～70<R<150～300では純カルマン渦列領域とあって、それぞれのシミュレーション結果から得られた動画を見ていると、まあ大外しはしてなさそうに思う。
ところでこの引数は呼び出し時に省略することができて、省略した場合にはデフォルト値としてν=0.0が入る。この値が何に対応しているかというと、これはレイノルズ数が無限大、つまり粘性力は何ら働いていないということになる。
円柱周りの流れはレイノルズ数がある値（O(102~3)あたりだったように記憶している）よりも高くなると流れの3次元性の影響が無視できなくなり、2次元のシミュレーションでは原理的に正しい流れを再現できなくなる。
デモンストレーションの動画は一見それっぽく見えても、実は流れの物理としての正しさは保証されていないことには注意しておきたい。
Follow-up
『粘性流体力学 (航空宇宙工学テキストシリーズ)』によれば非圧縮性流体近似においては以下の圧力に関する式を得ることができ、これはポアソン方程式である。「非圧縮性流れにおいては、どこか1か所で何かが起こると、その影響が流れ場全体に及ぶことを意味している。 (p.140)」
（式 (3.10.6) は長いので省略……）
前回記事では圧力の楕円形偏微分方程式に対する解法に言及したけれど、これは解法のテクニカルな観点で楕円形の方程式を扱っているというだけではなく、支配方程式の性質としてそういう側面があるということ。
おわりに
Juliaでサメを泳がせる前のPluto.jlの事前練習と、シミュレーションで注意しておくべき点を補足した。
コンピューターの性能は日進月歩で向上し、解法は年を追うごとに洗練され、またJuliaのような新たなプログラミング言語が登場して、シミュレーションを取り巻く環境は常に変化し続けている。 しかし工学におけるシミュレーションが捉えるべきは物理であって、かつてカルマンが見た100年前の連続体の振る舞いと何ら変わりのないことには留意する必要がある。 一方のルールは物理法則として不変であり、しかしもう一方のルールは変わりつづけていく、そうした中での生存戦略をどうすれば良いか？ 過去の蓄積が効いてくる領域と、今夜勝ちたいバランスをとるには……という感じで、先日のResearchat.fmにあった研究のキャラセレ問題を繰り返し聴きながら考えを巡らせ続けている。■
ライセンス
The Julia Language Blog / MIT License
Copyright (c) 2011-2019: Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and other contributors:
https://github.com/JuliaLang/www.julialang.org/contributors
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Pluto.jl / MIT License
https://github.com/fonsp/Pluto.jl
Copyright (c) 2020 Fons van der Plas
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
WaterLily.jl / MIT "Expat" License
https://github.com/weymouth/WaterLily.jl
Copyright (c) 2020: Gabriel Weymouth.
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
