こんにちは（@t_kun_kamakiri）(‘◇’)ゞ
最近、モデリカのためのオープンソースであるOpenModelicaとか使って遊んでいます。
本記事の本筋は以下の内容です。
OpenModelicaを使ってバネ＋ダンパーによる減衰振動のモデル作りから計算までを体験してみましょう(^^)/
OpenModelicaについては前の記事でちょろっと説明しました。
OpenModelicaのバージョンは、「1.13.2 (64-bit)」を使っています。
はじめてのModelicaプログラミング -1日で読める わかる Modelica入門- (MBD Lab Series)
Modelicaによるシステムシミュレーション入門 MBD Lab Series
バネ＋ダンパーによる減衰振動をOpenModelicaで体験してみよう
バネの振動における理論を説明した後、
OpenModelicaを使ってバネの振動の解析をしたいと思います(^^)/
バネとダンパーによる減衰振動の理論は？
ダンパーの振動の運動方程式
と、このように書くことができます。
- バネ定数：\(k=10\)[N/m]
- 減衰定数：\(k=10\)[N\dot s/m]
- 質量：\(m=10\)[kg]
とします。
(1)式は以下のように置くことで、見通しの良い形に変形することができます。
- \(\gamma=\frac{D}{2m}\)
- \(\omega_{0}=\sqrt{\frac{k}{m}}\)
ダンパーの振動の運動方程式
(2)式の一般解、\(x(t)=e^{\lambda t}\)と置き(2)式に代入して以下の特性方程式より\(\lambda\)を決めることで求めることができます。
「減衰定数から来る項\(\gamma\)」と「振動から来る項\(\omega_{0}\)」の大小関係によって、
以下の4つの解に分類することができます。
※2回の微分方程式なので、(1)から求めた\(x\)では「積分定数(\(A,B\))」が2つあります。
今回、OpenModelicaで解析しようとする条件は、
- バネ定数：\(k=10\)[N/m]
- 減衰定数：\(k=10\)[N・s/m]
- 質量：\(m=10\)[kg]
\(\gamma=0.5\)、\(\omega_{0}=1\)なので、
No.2（減衰振動）の解析をすることになります。
この「積分定数」は、初期条件を定めることで決定されます。
初期条件を以下のように決めます。
- 質点の初期位置：\(x_{0}=0[m]\)
- 質点の初速度：\(v_{0}=10[m/s]\)
そうすると、以下のように解が定まります。
(1)式の特殊解
※(4)式から初期条件を定めて得られた解なので、特殊解と呼ばれています。
※\(\omega=\sqrt{\omega_{0}^2-\gamma^2}\)
具体的に数字を入れると、
となります。
バネの振動をOpenModelicaでモデル作成
では、OpenModlicaを使ってバネのモデル構築をしてみましょう。
操作方法は以下の動画を参考にしてください。
OpenModelicaで得られた解と理論解(5)式とを比較
OpenModelicaの結果を「csvファイル」として出力ができるので、csvファイルにして理論解の波形と重ねてみます。
OpenModelicaのcsvファイル出力
※出力方法は以下の記事の後半を参考にして下さい。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 | import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt #理論解 m = 10 # 質量[kg] k = 10 # 剛性[N/m] state0 = [0.0, 10.0] # 初期値[x0, v0] omega = np.sqrt(k/m) #固有角振動数 Damp = 10 #減衰項 gamma=Damp/2/m omega_gamma=np.sqrt(omega**2-gamma**2) t0 = 0 # 初期時間[s] tf = 10 # 終了時間[s] dt = 0.005 # 時間刻み[s] t = np.arange(t0, tf+dt, dt) # 時間軸配列 x_thory = np.exp(-gamma*t)*((state0[1]+gamma*state0[0])/omega_gamma * np.sin(omega_gamma * t)+state0[0]*np.cos(omega_gamma*t)) #配列要素として読み込む array=np.genfromtxt('./Damp_Spring.csv',delimiter=',',dtype=float,skip_header=2) array_t=np.transpose(array) #リスト型に変換 time_ana=list(array_t[0]) x_ana=list(array_t[1]) #グラフの作成 plt.plot(time_ana,x_ana,label="OpenModelica") plt.plot(t,x_thory,label="theory") ##フォントサイズの指定 plt.xticks(fontsize=7) plt.yticks(fontsize=7) ##グリッドの有無 plt.grid(True) ##x軸とy軸のラベル plt.xlabel("time[sec]") plt.ylabel("x[m]") #凡例の位置 plt.legend(loc='best') ##グラフを描く plt.show() |
グラフの比較
一致しているので、グラフが重なっています(^^)/
Follow @t_kun_kamakiri