#!/usr/bin/env python
# coding: utf-8

# ### TestCC_includeEOF - 反応物と生成物の生成エンタルピーの差から発熱量を計算する
# TestCC　では、燃料ガスの低位発熱量(Lower Heating Value of fuel)をパラメータ HH に設定しました。このモデルに含まれる　CombustionChamberモデルは、以下のように設定することによって、反応物と生成物の生成エンタルピー(enthalpy of formation)の差から発熱量を計算することができるようになります。
# 
# * すべての Media （物性モデル) について、excludeEnthalpyOfFormation = false にする。
# * HH = 0 にする。
# 
# 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 の場合です。

# In[1]:


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


# TestCC のソースコードを改編して TestCC_includeEOF を作成します。主な改編点は次のようになります。
# 
# * 全体モデル TestCC_include の replaceable package として Air、Fuel、Exhaust を宣言し、各コンポーネントはこれらを redeclare するようにしました。
# * Air, Fuel, Exhaust の宣言で、 excludeEnthalpyOfFormation = false を設定しました。
# * CombustionChamber1 のパラメータで、HH=0 としました。
# 
# テキストベースで作成したので、コンポーネントの位置などに関する情報は記述してません。

# In[2]:


get_ipython().run_cell_magic('writefile', 'TestCC_includeEOF.mo', 'model TestCC_includeEOF\n   import SI = Modelica.SIunits;\n     \n  replaceable package Air = ThermoPower.Media.Air(excludeEnthalpyOfFormation = false);\n  replaceable package Fuel = ThermoPower.Media.NaturalGas(excludeEnthalpyOfFormation = false);\n  replaceable package Exhaust = ThermoPower.Media.FlueGas(excludeEnthalpyOfFormation = false);\n   \n  ThermoPower.Gas.SourceMassFlow Wcompressor(redeclare package Medium = Air, w0=158, T=616.95);\n  \n  ThermoPower.Gas.CombustionChamber CombustionChamber1(\n      redeclare package Air = Air, redeclare package Fuel = Fuel, redeclare package Exhaust = Exhaust,\n     initOpt=ThermoPower.Choices.Init.Options.steadyState, HH= 0.0, pstart=11.2e5,  V=0.1, S=0.1);\n  \n  ThermoPower.Gas.SourceMassFlow Wfuel(redeclare package Medium = Fuel, use_in_w0=true);  \n  ThermoPower.Gas.PressDrop PressDrop1(redeclare package Medium = Exhaust,\n      FFtype=ThermoPower.Choices.PressDrop.FFtypes.OpPoint,\n      rhonom=3.3, wnom=158.9, pstart=11.2e5, dpnom=0.426e5);  \n  ThermoPower.Gas.SensT SensT1(redeclare package Medium = Exhaust);  \n  Modelica.Blocks.Sources.Step Step1(startTime=0.5, height=-0.3, offset=3.1);  \n  ThermoPower.Gas.ValveLin ValveLin1(redeclare package Medium = Exhaust, Kv=161.1/9.77e5);\n  ThermoPower.Gas.SinkPressure SinkP1(redeclare package Medium = Exhaust);\n  Modelica.Blocks.Sources.Constant Constant1(k=1);\n  inner ThermoPower.System system;\n\nequation\n  connect(Wfuel.flange, CombustionChamber1.inf);\n  connect(Wcompressor.flange, CombustionChamber1.ina);\n  connect(CombustionChamber1.out, PressDrop1.inlet);\n  connect(PressDrop1.outlet, SensT1.inlet);\n  connect(Step1.y, Wfuel.in_w0);\n  connect(ValveLin1.outlet, SinkP1.flange);\n  connect(SensT1.outlet, ValveLin1.inlet);\n  connect(Constant1.y, ValveLin1.cmd);\nend TestCC_includeEOF;\n')


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

# In[3]:


from pymodelica import compile_fmu
fmu = compile_fmu('TestCC_includeEOF', 'TestCC_includeEOF.mo')


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

# In[4]:


from pyfmi import load_fmu
model = load_fmu(fmu)
res = model.simulate(final_time=5)


# シミュレーション結果から、燃焼室温度、燃料ガス温度、コンプレッサー空気温度を抽出します。

# In[5]:


t = res['time']
Tc = res['CombustionChamber1.fluegas.T']
Tf = res['Wfuel.gas.T']
Ta = res['Wcompressor.gas.T']


# 比較するために、DyMat を使って TestCC の燃焼室温度も抽出します。

# In[6]:


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 の燃焼室温度、燃料ガス温度、コンプレッサ空気温度をプロットします。

# In[7]:


get_ipython().run_line_magic('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]')


# In[ ]: