Modelica Standard Library の MixtureGasNasa を使った例題として、ThermoPower ライブラリの CombustionChamber (燃焼室) のテストモデル TestCC のシミュレーションを実行します。
このモデルは、コンプレッサから出た空気(Compressor Air)と燃料ガス(Fuel Gas)が燃焼室で混合し、燃料ガス内のCH4が完全燃焼します。そして、加熱された排ガス(Exhaust Gas) が排出されます。燃焼室モデルのパラメータとして燃料ガスの低位発熱量(Lower Heating value of fuel) HH = 41.6e6 J/kg を与えています。 0.5秒に燃料ガスの質量流量を 3.1 kg/s から 2.8 kg/s に減少させます。モデルの概要は、
に解説を書きました。
ThermoPower ライブラリを GitHub mirror や svn://svn.code.sf.net/p/thermopower/svn/trunk から入手します。 そして、環境変数 MODELICAPATH に、MSLのパスとこのライブラリのパスを設定します。
ライブラリのディレクトリが Docker コンテナから見て。 /home/jmodelica/jmodelica/lib/TermoPower 3.1 になるように配置するものとします。
import os
os.environ['MODELICAPATH'] = '/usr/local/jmodelica/ThirdParty/MSL:/home/jmodelica/jmodelica/lib'
ライブラリのディレクトリが C:\Users\ユーザー名\jmodelica\lib\ThermoPower 3.1 になるように配置します。
import os
os.environ['MODELICAPATH'] = 'C:\\JModelica.org-2.14\\install\\ThirdParty\\MSL;C:\\Users\\ユーザー名\\jmodelica\\lib'
import os
os.environ['MODELICAPATH'] = '/usr/local/jmodelica/ThirdParty/MSL:/home/jmodelica/jmodelica/lib'
モデルをコンパイルして FMU を生成ます。
from pymodelica import compile_fmu
fmu = compile_fmu('ThermoPower.Test.GasComponents.TestCC')
FMU をロードしてシミュレーションを実行します。
from pyfmi import load_fmu
model = load_fmu(fmu)
res = model.simulate(final_time=5)
シミュレーション結果を確認していきます。
燃料ガスの成分、流量、温度、圧力を抽出します。燃料ガス N2、CO2、CH4 の混合ガスで、最初と最後の質量分率 (Mass Fraction) と質量流量 (Mass Flow Rate)、温度 (Temperature) は以下のようになります。
t = res['time']
Wf_N2 = res['Wfuel.gas.X[1]']
Wf_CO2 = res['Wfuel.gas.X[2]']
Wf_CH4 = res['Wfuel.gas.X[3]']
Wf_w = res['Wfuel.w']
Wf_T = res['Wfuel.gas.T']
Wf_p = res['Wfuel.gas.p']
print('Mass Fraction')
print(3*'%10.8f ' %(Wf_N2[0], Wf_CO2[0], Wf_CH4[0]))
n=len(t)-1
print(3*'%10.8f ' %(Wf_N2[n], Wf_CO2[n], Wf_CH4[n]))
print('Mass Flow Rate '+2*'%10.8f ' %(Wf_w[0], Wf_w[n]))
print('Temperature '+2*'%10.8f ' %(Wf_T[0], Wf_T[n]))
燃料ガスの質量分率、質量流量、温度をプロットします。流量は t=0.5 秒で最初 3.1 kg/s から 2.8 kg/s に変化します。 温度は 300 K です。
%matplotlib notebook
import matplotlib.pyplot as plt
plt.figure(1)
plt.subplot(3,1,1)
plt.title('Fuel Gas')
plt.plot(t,Wf_N2,t,Wf_CO2,t,Wf_CH4)
plt.ylim(0,1)
plt.xlim(0,5)
plt.xticks(color="None")
plt.legend(['N2','CO2','CH4'], loc='center right', ncol=3, title='Mass Fraction')
plt.ylabel('[kg/kg]')
plt.subplot(3,1,2)
plt.plot(t, Wf_w)
plt.xlim(0,5)
plt.xticks(color="None")
plt.ylabel('[kg/s]')
plt.legend(['Mass Flow Rate'])
plt.subplot(3,1,3)
plt.plot(t,Wf_T)
plt.xlim(0,5)
plt.ylabel('[K]')
plt.legend(['Temperature'])
plt.xlabel('time [s]')
コンプレッサーの空気の成分の質量分率、質量流量、温度、圧力を抽出します。
空気は、O2, H2O, Ar, N2 の混合ガスです。最初と最後の質量分率は次のようになります。空気は予め加熱され 616.95 K になっています。
Wc_O2 = res['Wcompressor.gas.X[1]']
Wc_H2O = res['Wcompressor.gas.X[2]']
Wc_Ar = res['Wcompressor.gas.X[3]']
Wc_N2 = res['Wcompressor.gas.X[4]']
Wc_w = res['Wcompressor.w']
Wc_T = res['Wcompressor.gas.T']
Wc_p = res['Wcompressor.gas.p']
print('Mass Fraction')
print(4*'%10.8f ' %(Wc_O2[0], Wc_H2O[0], Wc_Ar[0], Wc_N2[0]))
print(4*'%10.8f ' %(Wc_O2[n], Wc_H2O[n], Wc_Ar[n], Wc_N2[n]))
print('Mass Flow Rate '+2*'%10.8f ' %(Wc_w[0], Wc_w[n]))
print('Temperature '+2*'%10.8f ' %(Wc_T[0], Wc_T[n]))
空気の質量分率、質量流量、温度をプロットします。
plt.figure(2)
plt.subplot(3,1,1)
plt.title('Compressor Air')
plt.plot(t,Wc_O2,t,Wc_H2O,t,Wc_Ar,t,Wc_N2)
plt.ylim(0,1)
plt.xlim(0,5)
plt.xticks(color="None")
plt.legend(['O2','H2O','Ar','N2'], loc='center right', ncol=4, title='Mass Fraction')
plt.ylabel('[kg/kg]')
plt.subplot(3,1,2)
plt.plot(t, Wc_w)
plt.xlim(0,5)
plt.xticks(color="None")
plt.ylabel('[kg/s]')
plt.legend(['Mass Flow Rate'])
plt.subplot(3,1,3)
plt.plot(t,Wc_T)
plt.xlim(0,5)
plt.ylabel('[K]')
plt.xlabel('time [s]')
plt.legend(['Temperature'])
燃焼室内で燃料ガスと空気が混合して、CH4が完全燃焼して発熱します。 排ガス (Exhaust Gas) の成分の質量分率、質量流量、温度、圧力を抽出します。排ガスの質量流量は CombustionChamber1 に接続された PressDrop1 の入口の質量流量を使用します。 0.5 s で燃料の流量が変化するので、排ガスの成分、質量流量、温度も変化します。
Cc_O2 = res['CombustionChamber1.fluegas.X[1]']
Cc_Ar = res['CombustionChamber1.fluegas.X[2]']
Cc_H2O = res['CombustionChamber1.fluegas.X[3]']
Cc_CO2 = res['CombustionChamber1.fluegas.X[4]']
Cc_N2 = res['CombustionChamber1.fluegas.X[5]']
Cc_w = res['PressDrop1.w']
Cc_T = res['CombustionChamber1.fluegas.T']
Cc_p = res['CombustionChamber1.fluegas.p']
print('Mass Fraction')
print(5*'%10.8f ' %(Cc_O2[0], Cc_Ar[0], Cc_H2O[0], Cc_CO2[0], Cc_N2[0]))
print(5*'%10.8f ' %(Cc_O2[n], Cc_Ar[n], Cc_H2O[n], Cc_CO2[n], Cc_N2[n]))
print('Mass Flow Rate '+2*'%10.8f ' %(Cc_w[0], Cc_w[n]))
print('Temperature '+2*'%10.8f ' %(Cc_T[0], Cc_T[n]))
排ガスの質量分率、質量流量、温度をプロットします。
plt.figure(3)
plt.subplot(3,1,1)
plt.title('Exhaust Gas')
plt.plot(t,Cc_O2,t,Cc_Ar,t,Cc_H2O,t,Cc_CO2,t,Cc_N2)
plt.ylim(0,1)
plt.xlim(0,5)
plt.xticks(color="None")
plt.legend(['O2','Ar','H2O','CO2','N2'], loc='center right', ncol=5, title='Mass Fraction')
plt.ylabel('[kg/kg]')
plt.subplot(3,1,2)
plt.plot(t, Cc_w)
plt.xlim(0,5)
plt.xticks(color="None")
plt.ylabel('[kg/s]')
plt.legend(['Mass Flow Rate'])
plt.subplot(3,1,3)
plt.plot(t,Cc_T)
plt.xlim(0,5)
plt.ylabel('[K]')
plt.xlabel('time [s]')
plt.legend(['Temperature'])
このモデルは、Compressor Air の質量流量と Fuel Gas の質量流量、バルブ出口の SinkP1の圧力を規定しているので、CombustionChamber1 (燃焼室)の Exhaust Gas の圧力はシミュレーション結果として得られます。Wfuel (燃料源) や Wcompressor (コンプレッサ) の圧力は、接続されている CombustionChamber1 の圧力と等しくなります。
plt.figure(4)
plt.title('Pressure')
plt.plot(t, Wf_p, t, Wc_p, t, Cc_p)
plt.ylim(10e5,12e5)
plt.xlim(0,5)
plt.xlabel('time [s]')