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)
Final Run Statistics: --- Number of steps : 84 Number of function evaluations : 103 Number of Jacobian evaluations : 3 Number of function eval. due to Jacobian eval. : 21 Number of error test failures : 0 Number of nonlinear iterations : 95 Number of nonlinear convergence failures : 0 Number of state function evaluations : 86 Number of time events : 1 Solver options: Solver : CVode Linear multistep method : BDF Nonlinear solver : Newton Linear solver type : DENSE Maximal order : 5 Tolerances (absolute) : [ 1.00000000e+00 5.00000000e-04 1.00000000e-07 1.00000000e-07 1.00000000e-07 1.00000000e-07 1.00000000e-03] Tolerances (relative) : 0.0001 Simulation interval : 0.0 - 5.0 seconds. Elapsed simulation time: 0.0150210857391 seconds.
シミュレーション結果を確認していきます。
燃料ガスの成分、流量、温度、圧力を抽出します。燃料ガス 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]))
Mass Fraction 0.02000000 0.01200000 0.96800000 0.02000000 0.01200000 0.96800000 Mass Flow Rate 3.10000000 2.80000000 Temperature 300.00000000 300.00000000
燃料ガスの質量分率、質量流量、温度をプロットします。流量は 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]')
Text(0.5,0,u'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]))
Mass Fraction 0.23000000 0.01500000 0.00500000 0.75000000 0.23000000 0.01500000 0.00500000 0.75000000 Mass Flow Rate 158.00000000 158.00000000 Temperature 616.95000000 616.95000000
空気の質量分率、質量流量、温度をプロットします。
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'])
<matplotlib.legend.Legend at 0x7f9d7f32e1d0>
燃焼室内で燃料ガスと空気が混合して、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]))
Mass Fraction 0.15126641 0.00490379 0.05654653 0.05133045 0.73595282 0.15875311 0.00491294 0.05259591 0.04644949 0.73728856 Mass Flow Rate 161.10000000 160.80000075 Temperature 1266.89896971 1210.67809888
排ガスの質量分率、質量流量、温度をプロットします。
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'])
<matplotlib.legend.Legend at 0x7f9d7f1994d0>
このモデルは、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]')
Text(0.5,0,u'time [s]')