振動の運動方程式をエクセルVBAを使って解く

資料請求番号:SH43 TS41

Runge-Kutta法を使った振動計算

常微分方程式を使うと、様々な自然現象を記述することができます。
その常微分方程式の数値解法としてルンゲクッタ(Runge-kutta)法があります。

※Runge-kutta法の詳しい解説はこちら


今回は、振動運動の運動方程式をRunge-kutta法で解くためのプログラムをExcelのVBAで作成しました。
その備忘録として、本記事にソースコードを書き残し、計算結果を実現象と照らし合わせてみました。

※ソースコードはこちら

振動の運動方程式

ストーク

今回、Runge-kutta法で解く方程式は以下の通りだ。振動の運動方程式。

シママ

xを変位とすればxの2回微分で加速度になるから、この式はma=Fの運動方程式の形をしているのね。

ストーク

せや。右辺第一項はフックの法則を表している。バネを引っ張れば、それに対して戻ろうとする力が働くだろう?
だから-kxになってるたい。

シママ

第2項はこれ・・・抵抗力かしら?なんでxの微分がついているのかしら?

ストーク

確かに抵抗力なんやが、これは・・・そうやな・・・
水の中で手をゆっくり動かしたときと早く動かしたときとどっちが手を前に進めにくいか?で考えれば
抵抗力は運動の速度に比例するってことが分かると思うたい。

シママ

なるほど・・・。

ストーク

第3項は強制振動による力。外部からmF0cosωtだけ力が加わっていると考える。力を加えるんだから符号はプラスやな。

ストーク

それで、両辺mで割って・・・

ストーク

と、おけば、振動に関する2階の微分方程式が完成する。これをRunge-kutta法で解くたい。

振動の運動方程式(2階微分方程式)をRunge-Kutta法で解く

2階微分方程式→2個の連立方程式

シママ

あれ?でも・・・Runge-kuttaで解く方程式って、こんな形の1階常微分方程式じゃなかったかしら?

シママ

この方程式、Runge-kuttaで解けるの?

ストーク

解ける。振動の運動方程式は、以下の2つの常微分方程式の連立と見ればいいだろう?
それで、2つの連立微分方程式はオマエ、こないだ実装したやないか。

シママ

あ、ホントだ。これならできそうね。
2つの方程式をt, x, vの関数としてみて・・・

シママ

こんな風にRunge-Kutta法を適用すれば出来そう!

ストーク

おっし。じゃあ、振動の運動方程式をRunge-kuttaで解くマクロを作ってみろよ。

シママ

うん。すぐにできるわよ。前回のソースコードをこの方程式に変えればいいだけでしょう?

・・・・

シママ

できた。

ストーク

おお、早いな。じゃあ、検証していこうか。

ソースコードと検証

ソースコード

シママ

こんな感じのワークシートで・・・

シママ

このソースコードを走らせると・・・

/* –vibration formula rk-method ——– */

Public omega0, gamma, force, omega As Double

Sub rkvib()
Dim x0, v0, init, ed, h, h2, kx(4), kv(4) As Double
Dim i, j As Integer

‘Calculate propertv
init = Cells(1, 3)
ed = Cells(2, 3)
h = Cells(3, 3)
h2 = Cells(4, 3)

‘Constants
omega0 = Cells(5, 3)
gamma = Cells(6, 3)
force = Cells(7, 3)
omega = Cells(8, 3)

‘Initial Condition
x0 = Cells(4, 7)
v0 = Cells(4, 8)
x = x0
v = v0
t = init

‘Runge-kutta roop

For i = 0 To ((ed – init) / h)

j = 6 + i / h2

Cells(1, 6) = i
Cells(2, 6) = j – 6

If i Mod h2 = 0 Then

Cells(j, 6) = t
Cells(j, 7) = x
Cells(j, 8) = v

End If

kx(1) = h * F1(t, x, v)
kv(1) = h * F2(t, x, v)
kx(2) = h * F1(t + h / 2, x + kx(1) / 2, v + kv(1) / 2)
kv(2) = h * F2(t + h / 2, x + kx(1) / 2, v + kv(1) / 2)
kx(3) = h * F1(t + h / 2, x + kx(2) / 2, v + kv(2) / 2)
kv(3) = h * F2(t + h / 2, x + kx(2) / 2, v + kv(2) / 2)
kx(4) = h * F1(t + h, x + kx(3), v + kv(3))
kv(4) = h * F2(t + h, x + kx(3), v + kv(3))
nx = x + (kx(1) + 2 * kx(2) + 2 * kx(3) + kx(4)) / 6
nv = v + (kv(1) + 2 * kv(2) + 2 * kv(3) + kv(4)) / 6
nt = t + h

t = nt
x = nx
v = nv

Next i

End Sub

Function F1(ByVal t As Double, ByVal x As Double, ByVal v As Double) As Double
F1 = v
End Function

Function F2(ByVal t As Double, ByVal x As Double, ByVal v As Double) As Double
F2 = -omega0 ^ 2 * x – 2 * gamma * v + force * Cos(omega * t)
End Function

/* –vibration formula rk-method End——– */

シママ

ジャジャ~ン!!

ストーク

(・・・かわいいな。)

シママ

とりあえずは、単振動、つまりω=γ=F0=0でやってみたよ。その方が検証しやすいでしょ?

ストーク

ああ、そうだな。

単振動

上記の振動の運動方程式のうち、ω=γ=F0=0とし、抵抗力も外力もないとした場合、つまり

で表される運動を単振動と言います。

ストーク

初期条件はx(0)=1, v(0)=0でやったのか。

シママ

うん。バネを引っ張った状態ね。そこから手を離したら、バネは余弦波の振動をするでしょ?

ストーク

せや。確かにコサインの波が出来てるな。

シママ

それで、ω0を大きくしたら振動数は大きくなって波長が小さくなるところとか、
ω0をπにしたら、この曲線はcos(πt)になって1秒で半周、2秒で1周するところとかちゃんと再現できてるでしょ?

ストーク

確かに。

シママ

ついでだからω0を0.1刻みで0~1まで変化させたときの動画を作ってみたよ!

ストーク

へぇ~。こんなことができるんやなぁ~。現象の表現がより視覚的になってええな!

シママ

この動画の作り方のソースコードはまた今度まとめるね!

減衰振動

γに正の値を与えると、減衰しながら振動します。振動の減衰度合いが強すぎると振動しないまま運動が止まります。
減衰しながら振動する運動を減衰振動、振動せずに停止する運動を過減衰といいます。
減衰振動と過減衰の境目を臨界減衰といいます。臨界減衰を呈しているときγω0を満たします。

シママ

ω0を2にして、減衰係数を上げていく。γの値がω0を超えたら過減衰になるのよ。
この性質もちゃんと表現できてるでしょ?

ストーク

確かにな。解析的に解くと、減衰振動・臨界減衰・過減衰とでは解曲線が変わってきてしまうからなぁ~。
うっかり、関数を切り替え忘れてて間違ったデータを出すこともあるし、そもそも切り替えがめんどくさかったりするから、数値解析は便利やな。

シママ

そうなのよ~。解析的に解けるとは言っても、その解析解が場合分けされているような現象だと、サッと答えを出しにくいのよね。
そういう意味でも、支配方程式を丸ごと数値解析してしまった方が早い場合もあるわよね。もちろん、精度が良ければの話だけどね。

シママ

こちらもγを0~3まで変化させたときの動画を作ってみたよ!

ストーク

おおっ。減衰振動・臨界減衰・過減衰の違いがよく分かってええな!

減衰振動+強制振動

本来は減衰して振動が止まってしまうところを、何らかの外力(強制振動)を加えることによって、振動を保たせ続けることができます。

シママ

本来、γ=0.1の減衰を受けて止まってしまう振動運動に対して・・・

シママ

例えば、F0=0.5 ω=1のF0cosωtという外力を加えれば、振動運動を続けることができるのよ。

ストーク

なるほどな。最初は安定しないかもしれないけど、最終的に減衰と外力が釣り合って、平衡状態になるんやな。

ストーク

これ、設計を上手くやれば最初から安定した振動が得られるかもな。

シママ

う~ん。言われてみればそうね。例えば、F0を強くしたりして0から4に変更すれば・・・・

シママ

あ、ホラ、F0=3がちょうどいいね!

ストーク

おおっ~。これなら外部環境に合わせてどんな外力を加えればいいのかわかるから、振動を使った機械系の設計ができるな!

共振

ストーク

振動を使った機械系で怖いのが共振だが・・・共振は再現できるのか?

シママ

もちろん!今、ω0=2, γ=0.1, F0=3, ω=1で安定している系があるじゃない?
共振は自己の周波数ω0と外部周波数ωが一緒になった時に起きるのよね。

ストーク

ふぅ~ん。良く知ってんじゃねぇか。

シママ

もう共振という危険な現象は、安全会議で耳タコになるくらい出てくるからね~・・・。実際、富士事業所にあるようなデッカイ機械が大きく揺れたらホント怖いし・・。

シママ

じゃあ、共振をシミュレーションしてみるね。ωω0に近づけていくのよ。

ストーク

おおっ、振幅が7倍になりよった!!

ストーク

ダッシュポット(減衰力)と外力の比によってはもっと大変なことになりそうたいね。

シママ

例えば、γが1/10になってγ=0.01になったとしたら・・・


※縦軸のスケールに注意

シママ

・・・25倍か~。怖いね。

ストーク

怖い。ウチにはタービンとか攪拌とか回るものが特に多いから気をつけないといかんたいね・・・。

まとめ

今回は振動の運動方程式をrunge-kutta法で解くための、エクセルのVBAプログラムを作成し、公開しました。
このプログラムを使用して、実際に単振動、減衰振動、強制振動の系について計算し、振動の運動の性質について説明しました。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です