TestCC では、燃料ガスの低位発熱量(Lower Heating Value of fuel)をパラメータ HH に設定しました。このモデルに含まれる CombustionChamberモデルは、以下のように設定することによって、反応物と生成物の生成エンタルピー(enthalpy of formation)の差から発熱量を計算することができるようになります。
TestCC にこのような改造を施したモデル TestCC_includeEOF を作成し、シミュレーションを実行します。
環境変数 MODELICAPATH に MSL のパスと ThermoPower ライブラリのパスを設定します。
Windows の場合は、
import os
os.environ['MODELICAPATH'] = 'C:\\JModelica.org-2.14\\install\\ThirdParty\\MSL;C:\\Users\\ユーザー名\\jmodelica\\lib
のように設定します。
以下は Mac の場合です。
import os
os.environ['MODELICAPATH'] = '/usr/local/jmodelica/ThirdParty/MSL:/home/jmodelica/jmodelica/lib'
TestCC のソースコードを改編して TestCC_includeEOF を作成します。主な改編点は次のようになります。
テキストベースで作成したので、コンポーネントの位置などに関する情報は記述してません。
%%writefile TestCC_includeEOF.mo
model TestCC_includeEOF
import SI = Modelica.SIunits;
replaceable package Air = ThermoPower.Media.Air(excludeEnthalpyOfFormation = false);
replaceable package Fuel = ThermoPower.Media.NaturalGas(excludeEnthalpyOfFormation = false);
replaceable package Exhaust = ThermoPower.Media.FlueGas(excludeEnthalpyOfFormation = false);
ThermoPower.Gas.SourceMassFlow Wcompressor(redeclare package Medium = Air, w0=158, T=616.95);
ThermoPower.Gas.CombustionChamber CombustionChamber1(
redeclare package Air = Air, redeclare package Fuel = Fuel, redeclare package Exhaust = Exhaust,
initOpt=ThermoPower.Choices.Init.Options.steadyState, HH= 0.0, pstart=11.2e5, V=0.1, S=0.1);
ThermoPower.Gas.SourceMassFlow Wfuel(redeclare package Medium = Fuel, use_in_w0=true);
ThermoPower.Gas.PressDrop PressDrop1(redeclare package Medium = Exhaust,
FFtype=ThermoPower.Choices.PressDrop.FFtypes.OpPoint,
rhonom=3.3, wnom=158.9, pstart=11.2e5, dpnom=0.426e5);
ThermoPower.Gas.SensT SensT1(redeclare package Medium = Exhaust);
Modelica.Blocks.Sources.Step Step1(startTime=0.5, height=-0.3, offset=3.1);
ThermoPower.Gas.ValveLin ValveLin1(redeclare package Medium = Exhaust, Kv=161.1/9.77e5);
ThermoPower.Gas.SinkPressure SinkP1(redeclare package Medium = Exhaust);
Modelica.Blocks.Sources.Constant Constant1(k=1);
inner ThermoPower.System system;
equation
connect(Wfuel.flange, CombustionChamber1.inf);
connect(Wcompressor.flange, CombustionChamber1.ina);
connect(CombustionChamber1.out, PressDrop1.inlet);
connect(PressDrop1.outlet, SensT1.inlet);
connect(Step1.y, Wfuel.in_w0);
connect(ValveLin1.outlet, SinkP1.flange);
connect(SensT1.outlet, ValveLin1.inlet);
connect(Constant1.y, ValveLin1.cmd);
end TestCC_includeEOF;
Overwriting TestCC_includeEOF.mo
モデルをコンパイルして FMU を作成します。
from pymodelica import compile_fmu
fmu = compile_fmu('TestCC_includeEOF', 'TestCC_includeEOF.mo')
FMU をロードしてシミュレーションを実行します。
from pyfmi import load_fmu
model = load_fmu(fmu)
res = model.simulate(final_time=5)
Final Run Statistics: --- Number of steps : 87 Number of function evaluations : 107 Number of Jacobian evaluations : 3 Number of function eval. due to Jacobian eval. : 21 Number of error test failures : 0 Number of nonlinear iterations : 99 Number of nonlinear convergence failures : 0 Number of state function evaluations : 89 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.0161700248718 seconds.
シミュレーション結果から、燃焼室温度、燃料ガス温度、コンプレッサー空気温度を抽出します。
t = res['time']
Tc = res['CombustionChamber1.fluegas.T']
Tf = res['Wfuel.gas.T']
Ta = res['Wcompressor.gas.T']
比較するために、DyMat を使って TestCC の燃焼室温度も抽出します。
import DyMat
d = DyMat.DyMatFile('ThermoPower_Test_GasComponents_TestCC_result.mat')
t1 = d.abscissa('CombustionChamber1.fluegas.T', valuesOnly=True)
Tc1 = d.data('CombustionChamber1.fluegas.T')
TestCC の燃焼室温度、 TestCC_includeEOF の燃焼室温度、燃料ガス温度、コンプレッサ空気温度をプロットします。
%matplotlib notebook
import matplotlib.pyplot as plt
plt.figure(1)
plt.plot(t1,Tc1,t,Tc,t,Tf, t,Ta)
plt.legend(['Exhaust Gas, HH = 41.6e6 J/kg', 'Exhaust Gas, include enthalpy of formation', 'Fuel Gas', 'Compressor Air'],bbox_to_anchor=(0.1,0.8), loc='upper left')
plt.ylim(0,1400)
plt.xlim(0,5)
plt.ylabel('Temperature [k]')
plt.xlabel('Time [s]')
Text(0.5,0,u'Time [s]')