ThermoPower.Test.GasComponents.TestCC

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のパスとこのライブラリのパスを設定します。

Macの場合

ライブラリのディレクトリが Docker コンテナから見て。  /home/jmodelica/jmodelica/lib/TermoPower 3.1 になるように配置するものとします。

import os
os.environ['MODELICAPATH'] = '/usr/local/jmodelica/ThirdParty/MSL:/home/jmodelica/jmodelica/lib'

Windowsの場合

ライブラリのディレクトリが  C:\Users\ユーザー名\jmodelica\lib\ThermoPower 3.1 になるように配置します。

import os
os.environ['MODELICAPATH'] = 'C:\\JModelica.org-2.14\\install\\ThirdParty\\MSL;C:\\Users\\ユーザー名\\jmodelica\\lib'
In [1]:
import os
os.environ['MODELICAPATH'] = '/usr/local/jmodelica/ThirdParty/MSL:/home/jmodelica/jmodelica/lib'

モデルをコンパイルして FMU を生成ます。

In [2]:
from pymodelica import compile_fmu
fmu = compile_fmu('ThermoPower.Test.GasComponents.TestCC')

FMU をロードしてシミュレーションを実行します。

In [3]:
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) は以下のようになります。

In [4]:
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 です。

In [5]:
%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]')
Out[5]:
Text(0.5,0,u'time [s]')

コンプレッサーの空気

コンプレッサーの空気の成分の質量分率、質量流量、温度、圧力を抽出します。

空気は、O2, H2O, Ar, N2 の混合ガスです。最初と最後の質量分率は次のようになります。空気は予め加熱され 616.95 K になっています。

In [6]:
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 

空気の質量分率、質量流量、温度をプロットします。

In [7]:
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'])
Out[7]:
<matplotlib.legend.Legend at 0x7f9d7f32e1d0>

排気ガス (Exhaust Gas)

燃焼室内で燃料ガスと空気が混合して、CH4が完全燃焼して発熱します。 排ガス (Exhaust Gas) の成分の質量分率、質量流量、温度、圧力を抽出します。排ガスの質量流量は CombustionChamber1 に接続された PressDrop1 の入口の質量流量を使用します。 0.5 s で燃料の流量が変化するので、排ガスの成分、質量流量、温度も変化します。

In [8]:
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 

排ガスの質量分率、質量流量、温度をプロットします。

In [9]:
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'])
Out[9]:
<matplotlib.legend.Legend at 0x7f9d7f1994d0>

燃料ガス、コンプレッサーの空気、排気ガスの圧力

このモデルは、Compressor Air の質量流量と Fuel Gas の質量流量、バルブ出口の SinkP1の圧力を規定しているので、CombustionChamber1 (燃焼室)の Exhaust Gas の圧力はシミュレーション結果として得られます。Wfuel (燃料源) や Wcompressor (コンプレッサ) の圧力は、接続されている CombustionChamber1 の圧力と等しくなります。

In [10]:
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]')
Out[10]:
Text(0.5,0,u'time [s]')
In [ ]: