Pythonで計算するために多自由度振動系を行列形式にする方法

コンピュータに連立方程式を計算させる時は、線形代数学で学んだ行列形式で方程式を表現する方が効率的です。ここではPythonをはじめ様々なコンピュータシミュレーションのために、多自由度振動系の運動方程式を行列形式へ変換する方法を紹介します。

多自由度振動系モデルの運動方程式立てる

作用している力を考えよう！

今回扱う多自由度振動系とは、以下に示すモデルを指します。このモデルは質点が2つあるので、2自由度を持つと言えます。

1自由度の場合は「Pythonで1自由度非減衰系の自由振動シミュレーション」で説明したように、非常に考えやすい式でした。しかし、このような複数の自由度を持つ場合でも考え方は同じです。

まずは質点に作用している力を考えましょう！

\(m_{1}\)に作用する力から\(m_{1}\)の運動方程式を立てる

運動方程式の基本は\(F=ma\)なので、左辺の力の合計を考えて、右辺の質量と加速度（これも慣性力という力）を考えます。

まずは質点\(m_{1}\)に作用する力ですが、質点自体にも力がかかる可能性があるため、\(f_{1}\)が挙げられますね。図の右方向を正とすればプラスの力です。

次に、\(m_{1}\)の左側にはばねがあるので、フックの法則からばねの伸びに比例して復元力が発生します。ばねの左端は固定点なので変位が発生せず、この復元力は\(x_{1}\)によってのみ発生します。但し、質点を左側に引っ張る方向に作用するので、力の符号は負となり\(-k_{1}x_{1}\)になります。

また、\(m_{1}\)の右側にもばねがあります。このばねに発生する復元力はフックの法則で考えると質点\(m_{1}\)と\(m_{2}\)の両方の相対変位\(x_{2}-x_{1}\)で決まるので\(k_{2}\left ( x_{2}-x_{1} \right )\)になります。

後は右辺の質量と加速度を加えてニュートンの運動方程式を立てると、式(1)が導出されます。

$$f_{1}-k_{1}x_{1}+k_{2}\left ( x_{2}-x_{1} \right )= m_{1}\ddot{x_{1}}　(1)$$

式(1)において、左辺に力\(f_{1}\)だけを置き、その他を全て右辺に移項、\(x_{1}\)と\(x_{2}\)についてまとめると式(2)が得られます。これが質点\(m_{1}\)の運動方程式になります。

$$f_{1}=m_{1}\ddot{x}_{1}+\left(k_{1}+k_{2}\right)x_{1}-k_{2}x_{2}　(2)$$

\(m_{2}\)に作用する力から\(m_{2}\)の運動方程式を立てる

では続いて質点\(m_{2}\)の運動方程式も求めてみましょう！

まずは左辺ですが、先ほどと同様に、\(m_{2}\)自体にも力がかかる可能性があるので\(f_{2}\)を正の方向で置きます。

そして質点\(m_{2}\)の左側にはばねがありますので、方向が負の復元力が働きます。このばねは\(x_{2}-x_{1}\)の相対変位で変位が決まり、ばね定数\(k_{2}\)をかけた\(-k_{2}\left ( x_{2}-x_{1} \right )\)が質点\(m_{2}\)に作用する復元力です。

質点\(m_{2}\)の右側は自由端となりもうばねは無いので、復元力は以上です。

式(1)と同様に右辺に質量と加速度をとると、式(3)が得られます。

$$f_{2}-k_{2}\left ( x_{2}-x_{1} \right )=m_{2}\ddot{x}_{2}　(3)$$

この式(3)も式変形をして整理すると、式(4)を得ることができ、これが質点\(m_{2}\)の運動方程式になります。

$$f_{2}=m_{2}\ddot{x}_{2}-k_{2}x_{1}+k_{2}x_{2}　(4)$$

連立方程式を行列表示させる

各質点における運動方程式が立てられたら、後は簡単です。まず式(2)と式(4)を一組にした連立方程式(5)を立てましょう。

\[
\begin{cases}
f_{1}=m_{1}\ddot{x}_{1}+\left(k_{1}+k_{2}\right)x_{1}-k_{2}x_{2}\\
f_{2}=m_{2}\ddot{x}_{2}-k_{2}x_{1}+k_{2}x_{2}
\end{cases}　(5)
\]

このような連立方程式が出来上がったら、「Pythonで線形代数！連立1次方程式を解く」で実施したように行列形式で書くことができます。もし忘れてしまった方はこちらの記事を参照して頂きたいと思います。

式(6)に式(5)を行列表示させた式を示します。この式は線形代数の行列計算を行うことで式(5)に戻すことも可能です。

\[
\begin{bmatrix}
f_{1}\\
f_{2}\end{bmatrix}
=
\begin{bmatrix}
m_{1} &0 \\
0 &m_{2}
\end{bmatrix}
\begin{bmatrix}
\ddot{x}_{1}\\
\ddot{x}_{2}
\end{bmatrix}
+
\begin{bmatrix}
k_{1}+k_{2} &-k_{2} \\
-k_{2} &k_{2}
\end{bmatrix}
\begin{bmatrix}
x_{1}\\
x_{2}
\end{bmatrix}　(6)
\]

上の式(6)は成分をいちいち書いていると大変なので、式(7)と書くこともあります。ここで、\(\mathbf{F}\)は外力ベクトル、\(\mathbf{M}\)は質量マトリクス、\(\mathbf{K}\)は剛性マトリクス、\(\mathbf{\ddot{x}}, \mathbf{x}\)はそれぞれ加速度ベクトルと変位ベクトルと呼びます。

$$\mathbf{F}=\mathbf{M}\mathbf{\ddot{x}}+\mathbf{K}\mathbf{x}　(7)$$

これで多自由度振動系の運動方程式を行列表示させる方法の説明は終了です。

基本はこの考え方で、ばねが並列になっていようが直列になっていようが、力の作用を考えて運動方程式を立てれば同じように整理することが可能です。

まとめ

本ページではこの後Pythonを使って振動の解析を行うために、まずは基礎知識として運動方程式を立てる方法を学びました。

多自由度の質点モデルは質点に作用する力とその方向を記述すれば運動方程式を作ることができます。

そして、運動方程式が連立する場合は行列表示にすることも、コンピュータで計算を行う場合は有効な手段です。