慣性振動の図示 - てれじょんのメモ帳

これはなに

　地獄の前期を終えて夏休みに入ったのでぼちぼちfortranとか総観気象学入門とかの進捗を生んでいるのですが、AtCoderの提出コードのように自分が書いたコードがすぐ見れないのがネックでした。そこでこのブログを思い出したのでなんかやったらちょこちょこ書いていきたいと思います。発見したことを発信するというよりは自分の備忘録として使いたいので、そんなに濃い内容はないかと思いますが生暖かい目で見てくれれば幸いです。
　今回は慣性振動をfortranとgnuplotで図示したのでそれについて書いていきます。地衡風が生じる理由として「どんな方向にものをぶん投げてもコリオリ力で(北半球なら)進行方向の右に曲げられて、そのうち定常状態に落ち着いて地衡風になるよね」っていうよくある説明が厳密には正確じゃないな～ってことが分かると思います。

本題

　出発点はおなじみ運動方程式です。空気塊の速度成分を $(u,v)$ とすれば

$\frac{du}{dt}-fv=-\frac{1}{\rho}\frac{\partial p}{\partial x}-\gamma u \\ \frac{dv}{dt}+fu=-\frac{1}{\rho}\frac{\partial p}{\partial y}-\gamma v$

で与えられます。ただし摩擦項は速度の大きさに比例するとしました。(厳密には $\nu \nabla^2 u_i$ です。)
ここで気圧傾度力は時間・空間的に一定として

$-\frac{1}{\rho}\frac{\partial p}{\partial x}=G_x,~~-\frac{1}{\rho}\frac{\partial p}{\partial y}=G_y$

とすれば見慣れた一階連立常微分方程式になるので

$W:=u+vi,~~G:=-G_x+G_yi$

とおいてやれば

$\frac{dW}{dt}=-(\gamma+fi)W+G$

になります。これを初期条件 $W(0)=W_0$ の下で解いて

$W=\left(W_0-\frac{G}{\gamma+fi}\right)e^{-(\gamma+fi)t}+\frac{G}{\gamma+fi}$

を得ます。これをいい感じに積分すれば変位になるので、初期条件 $x(0)+y(0)i=0$ のもとで積分して

$x+yi=\frac{1}{\delta}\left(W_0-\frac{G}{\delta}\right)(1-e^{-\delta t})+\frac{G}{\delta}t$

となります(ただし $\delta=\gamma+fi$ )。微分方程式の問題集に出てきそう。

　解が得られたので、後はfortranでコードを書いて解の様子を出力し、gnuplotで書くだけです。ソースコードは以下のようになりました。

module globals  !初期条件や定数を格納するモジュール
   real(8),parameter :: pi=3.141592653589793d0
   complex(8) :: G=(-1.0d0,1.0d0),W=(1.3d0,3.0d0)
   real(8),parameter :: f=0.7d0,gm=0.1d0
end module globals

module subporog !地衡風・地上風成分を出力するモジュール
   use globals
   implicit none
contains
   subroutine geo(t,n)
      integer :: i,n,fo=12
      real(8) t,dt,nt
      complex(8) :: gt,gma=(gm,f)
      dt=t/dble(n)
      open(fo,file='geostrophic(frictional)')
      do i=0,n
         nt=dble(i)*dt
         gt=G*nt/gma
         write(fo,*)real(gt),aimag(gt)
      enddo
      close(fo)
   end subroutine geo
         
end module subporog

program main !メイン部分
   use globals
   use subporog
   implicit none
   real(8) t,dt,nt
   integer :: n,i,fo=11
   complex(8) :: ft,gma = (gm,f)
   write(*,*)"input parameter t'(0<t<t')"
   read(*,*)t
   write(*,*)"input step n(integer)"
   read(*,*)n
   dt=t/dble(n)
   open(fo,file="f=07(frictional).d")
   do i=0,n
      nt=dble(i)*dt
      ft=(W-G/gma)*(1.0d0-exp(-gma*nt))/gma+G*nt/gma
      write(fo,*)real(ft),aimag(ft)
   enddo
   close(fo)

   call geo(t,n)

end program main

fortranは複素数型が定義されているんですね。メチャメチャ便利です。
パラメーター $t$ を $0\leq t\leq 55$ 、はじめにぶん投げる方向を $(1.3,3.0)$ 、気圧傾度力の方向を $(-1.0,1.0)$ として図示したものが以下になります。ステップ数は2000にしました。

f:id:jhonson1415:20200812203251p:plain — 図1:出力した画像

　摩擦に引きずられて地上風が地衡風の左側にズレていることや、十分時間が過ぎた後では摩擦を考慮したほうは地上風に近づいている一方、摩擦を無視した方ではいつまで経っても振動しているのが見えますね。風が地上風に落ち着くには摩擦の効果がデカいというのがよくわかります。
　地上風は等圧線と概ね30度くらいで交わるって言われてるのでパラメーターの初期値を適当に設定しすぎた感は否めませんが、fortranに慣れ親しむのが目的だったのでこれで良しとします。暇な人居たらいい感じにパラメーター設定してみてください。

最後に

　fortranは神。息をするように扱える様になりたいです。