見出し画像

Pythonライブラリ(科学技術計算ツール):Cantera


1.概要

 化学熱力学、化学反応、輸送現象などの問題を扱うオープンソースであるCanteraを紹介します。

1-1.Canteraの特徴

 Canteraはオープンソースであり、化学反応や熱力学の計算に使用できます。学習することで化学系での高度なシミュレーションができたり、プロセスシミュレーションを実施することが可能です。

 その他ソフトウェアとの比較は下図の通りです。

  • CHEMKIN (Reaction Design)

    • 世界的に使われている反応解析ソフトウエアパッケージです。Chemkin-II まではSandia National Lab.(USA)が配布していましたが、その後商用化されて現在に至っています。最新のCHEMKIN-PROにはReaction Path Analyzer や Particle Tracking Module なども備えられ、またMulti-Zone Engine Modelにも対応します。また、各種CFDソフトウエアとも対応しています。(記:2010/09/10)

  • Cantera

    • 当初からオープン・ソースで開発が進められている、フリーの反応解析ソフトウエアです。現在はコミュニティベースで開発が続けられており,2012年8月に各種機能が充実した最新版としてVer.2.0 がリリースされています。(記:2012/08/31)

  • Kintecus

    • Chemkin-II 形式の反応機構入力ファイルを変換して読み込むことができる、反応解析ソフトウエアパッケージです。アカデミックユーザーは登録すれば無料で使用できます。(記:2010/09/10)

  • Combustion Simulations (ELTE, Hungary / mirror site: U. Leeds, UK)

    • T. Turanyi, M. J. Pilling らが開発した下記の反応モデル解析ツールが公開されていますが、2002年2月で更新が止まっています。

      • KINALC: CHEMKIN-II の計算結果を処理するpost processor

      • FluxViewer: Java で書かれた 反応経路解析ツール

      • KINAL: 反応解析に用いる計算パッケージ

  • OpenFOAM:流体解析用のオープンソース。化学機構の計算も可能(OpenFOAM初級入門

1-2.Canteraでできること

 できることの一部は下記の通り。その他は公式Docsをご確認ください。

  1. 混合物の熱力学的特性と輸送特性を評価

  2. 化学平衡を計算

  3. 化学種の生成速度の評価

  4. 大規模反応機構による速度論シミュレーションの実施

  5. 一次元火炎のシミュレーション

  6. 反応経路解析

  7. 攪拌反応器のネットワークを使用したプロセスシミュレーションの作成

  8. 非理想流体のモデル化

1-3.用語集

 Canteraでの用語集は下記の通り

  • CK: Chemkin

  • CT: Cantera

  • CTI: Cantera Input (Removed in Cantera 3.0)

  • CTML: Cantera Markup Language (Removed in Cantera 3.0)

  • HKFT: Helgeson-Kirkham-Flowers-Tanger

  • HMW: Harvie, Møller, and Weare

  • IAPWS: International Association for the Properties of Water and Steam

  • MFTP: Mixture Fugacity ThermoPhase

  • PDSS: Pressure-dependent standard state

  • RT: Product of the gas constant (R) and the temperature

  • SHE: Single half-electrode

  • SP: Surface Problem

  • SS: Standard state

  • SSTP: SingleSpeciesTP (ThermoPhase)

  • STIT: SpeciesThermoInterpType

  • VCS: Villars Cruise Smith

  • VPSS: Variable Pressure Standard State

  • VPSSTP: Variable Pressure Standard State ThermoPhase

  • wrt: with respect to

2.環境構築

2-1.インストール

 公式Docs「Installing & Compiling」から要領を確認します。

https://cantera.org/

 Pythonで実施する場合は下記がありますが、私はPipで実施しました。
※Cantera 3.0.0はPython versions 3.8, 3.9, 3.10, and 3.11.

  • If you don't already have Python installed (or already use Conda), the easiest way to install the Cantera Python interface on all operating systems is by using Conda.

  • If you already have a different Python installation, Cantera can be installed using Pip.

[Terminal]
pip install cantera

3.QuickStart1

 まずは公式「Python Tutorial」に従って動作確認します。
※これだけだと理解が難しいですが、とりあえず触ってみます

https://cantera.org/tutorials/index.html#cantera-next-steps

3-1.Solutionクラス作成/モデル読み込み

 Canteraでは扱う物質をSolutionクラスで表現します。
 このクラスに相、化学種、反応などを定義したinput fileを渡すことでインスタンス化します。インスタンスをメソッドで実行(gas1())すると現在の状態(下記参照)を表示できます。

  • 温度 (Temperature)

  • 圧力 (Pressure)

  • 密度 (Density)

  • 相平均分子量

  • エンタルピー、内部エネルギー、エントロピー、ギブズエネルギー

  • 比熱容量 (c_p, c_v)

  • 種ごとのモル分率 (X)、質量分率 (Y)

[IN]
import cantera as ct
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

#DataFrameを並列で表示するためのクラス
class HorizontalDisplay:
    def __init__(self, *args):
        self.args = args

    def _repr_html_(self):
        template = '<div style="float: left; padding: 10px;">{0}</div>'
        return "\n".join(template.format(arg._repr_html_())
                         for arg in self.args)
        

print(f'Version: {ct.__version__}') #バージョンの確認

#物体表現(Object Representation)
gas1 = ct.Solution('gri30.yaml') #GRI- Mech 3.0の読み込み
print(gas1)
gas1()
[OUT]
Version: 3.0.1
<cantera.composite.Solution object at 0x000001CD358089E0>

  gri30:

       temperature   300 K
          pressure   1.0133e+05 Pa
           density   0.081894 kg/m^3
  mean mol. weight   2.016 kg/kmol
   phase of matter   gas

                          1 kg             1 kmol     
                     ---------------   ---------------
          enthalpy             26469             53361  J
   internal energy       -1.2108e+06        -2.441e+06  J
           entropy             64910        1.3086e+05  J/K
    Gibbs function       -1.9447e+07       -3.9204e+07  J
 heat capacity c_p             14311             28851  J/K
 heat capacity c_v             10187             20536  J/K

                      mass frac. Y      mole frac. X     chem. pot. / RT
                     ---------------   ---------------   ---------------
                H2                 1                 1           -15.717
     [  +52 minor]                 0                 0  

 3-1ー1.Classの簡易紹介

 Canteraには3つのクラスがあります。

  • Solution: 1つの「相(phase)」を表すコンテナ
    例: 気相、液相、固相、混合気体、溶液など

  • Interface: 界面や表面反応相を表すコンテナ
    例: 触媒表面、粒子表面、バルク相同士の接触面

  • SolutionArray: 複数の状態(state)情報を一括管理する配列コンテナ
    例: 計算過程で得られた状態を時系列で格納、空間的な分布を表すための多数のSolution状態の集合

 3-1ー2.Solutionクラスの簡易紹介

 SolutionクラスのAPIは下記の通り。一般的なインスタンス方法は、定義時にInput fileを渡します。
 Solutionでは下記メソッドで化学特性を計算できます。

  • thermo() :熱力学特性(エンタルピー, エントロピー, 内部エネルギー, 密度など)を計算

  • kinetics():反応速度論(反応機構に基づく生成消費速度、前進・逆反応速度定数など)にアクセス

  • transport():輸送特性(拡散係数、粘性、熱伝導率など)を取得

  • adjacent() 隣接相(他のSolution)との関係を管理

[API]
class cantera.Solution(infile='', name='', *, origin=None, 
                       yaml=None, thermo=None, species=(), 
                       kinetics=None, reactions=())

 3-1ー3.gri30.yaml

 CanteraではInput fileとして反応機構をYAMLファイルで表現します。今回使用したgri30.yamlの詳細は下記の通りです。

高温酸化反応の制御に関する化学反応論的考察* (高温酸化反応に対する NO 共存の影響)
https://chem.tut.ac.jp/keel/umc/modelsite.html

 参考としてGitHub公開のgri30.yamlの一部は下記の通りです。

[Breadcrumbscantera/data/gri30.yaml]
description: |-
  GRI-Mech Version 3.0 7/30/99  CHEMKIN-II format
  See README30 file at anonymous FTP site unix.sri.com, directory gri;
  WorldWideWeb home page http://www.me.berkeley.edu/gri_mech/ or
  through http://www.gri.org , under 'Basic  Research',
  for additional information, contacts, and disclaimer

  Updated webpage at http://combustion.berkeley.edu/gri-mech/version30/text30.html

generator: ck2yaml
input-files: [gri30.inp, gri30_thermo.dat, gri30_tran.dat]
cantera-version: 2.5.0
date: Wed, 11 Dec 2019 16:59:02 -0500

units: {length: cm, time: s, quantity: mol, activation-energy: cal/mol}

phases:
- name: gri30
  thermo: ideal-gas
  elements: [O, H, C, N, Ar]
  species: [H2, H, O, O2, OH, H2O, HO2, H2O2, C, CH, CH2, CH2(S), CH3, CH4,
    CO, CO2, HCO, CH2O, CH2OH, CH3O, CH3OH, C2H, C2H2, C2H3, C2H4, C2H5,
    C2H6, HCCO, CH2CO, HCCOH, N, NH, NH2, NH3, NNH, NO, NO2, N2O, HNO, CN,
    HCN, H2CN, HCNN, HCNO, HOCN, HNCO, NCO, N2, AR, C3H7, C3H8, CH2CHO,
    CH3CHO]
  kinetics: gas
  transport: mixture-averaged
  state: {T: 300.0, P: 1 atm}

species:
- name: H2
  composition: {H: 2}
  thermo:
    model: NASA7
    temperature-ranges: [200.0, 1000.0, 3500.0]
    data:
    - [2.34433112, 7.98052075e-03, -1.9478151e-05, 2.01572094e-08, -7.37611761e-12,
      -917.935173, 0.683010238]
    - [3.3372792, -4.94024731e-05, 4.99456778e-07, -1.79566394e-10, 2.00255376e-14,
      -950.158922, -3.20502331]
    note: TPIS78
  transport:
    model: gas
    geometry: linear
    well-depth: 38.0
    diameter: 2.92
    polarizability: 0.79
    rotational-relaxation: 280.0
- name: H
  composition: {H: 1}
  thermo:
    model: NASA7
    temperature-ranges: [200.0, 1000.0, 3500.0]
    data:
    - [2.5, 7.05332819e-13, -1.99591964e-15, 2.30081632e-18, -9.27732332e-22,
      2.54736599e+04, -0.446682853]
    - [2.50000001, -2.30842973e-11, 1.61561948e-14, -4.73515235e-18, 4.98197357e-22,
      2.54736599e+04, -0.446682914]
    note: L7/88
  transport:
    model: gas
    geometry: atom
    well-depth: 145.0
    diameter: 2.05
- name: O
~以下省略~

3-2.objectへのRead/Write

  作成したSolutionオブジェクトにある情報の読み込み/書込みが可能です。公式Docsに属性・メソッドの説明はありますが一覧表が無いため、組み込み関数のdir()で簡単に一覧表作成しました。

[IN]
import string
index_str = string.ascii_uppercase
index_list = [i for i in index_str]
#属性のリストを作成
_names = dir(gas1)
names = [i for i in _names if not i.startswith('_')]
#アルファベットの辞書を作成
names_dict = {c:[] for c in index_list}
for name in names:
    names_dict[name[0].upper()].append(name)
#最大数のリストを確認
numlist_max = max([len(v) for v in names_dict.values()])
#DataFrameの作成
columns = [i for i in range(numlist_max)]
df = pd.DataFrame(index=index_list, columns=columns)
for i, (k, v) in enumerate(names_dict.items()):
    if len(v) == 0:
        continue
    else:
        for i, n in enumerate(v):
            df.loc[k, i] = n

#NaNを空白に変換
df = df.fillna('')
df.to_excel('cantera_object_representation.xlsx')

[OUT]

 3-2-1.データ取得(読み込み)

 Solutionの属性からデータを取得することが出来ます。属性の一部は下記の通りです。

  • property CK_mode:Boolean to indicate if the chemkin interpretation is used.

  • property DP:Get/Set density [kg/m^3] and pressure [Pa].

  • property DPX:Get/Set density [kg/m^3], pressure [Pa], and mole fractions.

  • property DPY:Get/Set density [kg/m^3], pressure [Pa], and mass fractions.

  • property G:Get the total Gibbs free energy [J] represented by the Quantity.

  • property H:Get the total enthalpy [J] represented by the Quantity.

  • property HP:Get/Set enthalpy [J/kg or J/kmol] and pressure [Pa].

  • property HPX:Get/Set enthalpy [J/kg or J/kmol], pressure [Pa] and mole fractions.

  • property HPY:Get/Set enthalpy [J/kg or J/kmol], pressure [Pa] and mass fractions.

  • property P:Pressure [Pa].

  • property P_sat:Saturation pressure [Pa] at the current temperature.

  • property Pe:Get electron Pressure [Pa].

  • property S:Get the total entropy [J/K] represented by the Quantity.

  • property SP:Get/Set entropy [J/kg/K or J/kmol/K] and pressure [Pa].

  • property SPX:Get/Set entropy [J/kg/K or J/kmol/K], pressure [Pa], and mole fractions.

  • property SPY:Get/Set entropy [J/kg/K or J/kmol/K], pressure [Pa], and mass fractions.

  • property SV:Get/Set entropy [J/kg/K or J/kmol/K] and specific volume [m^3/kg or m^3/kmol].

  • property SVX:Get/Set entropy [J/kg/K or J/kmol/K], specific volume [m^3/kg or m^3/kmol], and mole fractions.

  • property SVY:Get/Set entropy [J/kg/K or J/kmol/K], specific volume [m^3/kg or m^3/kmol], and mass fractions.

  • property T:Temperature [K].

  • property TD:Get/Set temperature [K] and density [kg/m^3 or kmol/m^3].

  • property TDX:Get/Set temperature [K], density [kg/m^3 or kmol/m^3], and mole fractions.

  • property TDY:Get/Set temperature [K] and density [kg/m^3 or kmol/m^3], and mass fractions.

  • property TP:Get/Set temperature [K] and pressure [Pa].

  • property TPX:Get/Set temperature [K], pressure [Pa], and mole fractions.

  • property TPY:Get/Set temperature [K], pressure [Pa], and mass fractions.

  • property T_sat:Saturation temperature [K] at the current pressure.

  • property Te:Get/Set electron Temperature [K].

  • property U:Get the total internal energy [J] represented by the Quantity.

  • property UV:Get/Set internal energy [J/kg or J/kmol] and specific volume [m^3/kg or m^3/kmol].

  • property UVX:Get/Set internal energy [J/kg or J/kmol], specific volume [m^3/kg or m^3/kmol], and mole fractions.

  • property UVY:Get/Set internal energy [J/kg or J/kmol], specific volume [m^3/kg or m^3/kmol], and mass fractions.

  • property V:Get the total volume [m^3] represented by the Quantity.

  • property X:Get/Set the species mole fractions. Can be set as an array, as a dictionary, or as a string. Always returns an array:

  • property Y:Get/Set the species mass fractions. Can be set as an array, as a dictionary, or as a string. Always returns an array:

 実行時にエラーになったものはコメントアウトしました。

[IN]
print(gas1.CK_mode) #Bool型:Chemkinモードの確認
print(gas1.DP) #密度/Density[kg/m^3]と圧力/Pressure[Pa]の取得
print(gas1.DPX) #密度/Density[kg/m^3]と圧力/Pressure[Pa]とモル分率の取得
print(gas1.DPY) #密度/Density[kg/m^3]と圧力/Pressure[Pa]と質量分率の取得
# print(gas1.G) #ギブスエネルギー/Gibbs free energy[J]の取得
# print(gas1.H) #エンタルピー/Enthalpy[J]の取得

print('')
print(gas1.HP) #エンタルピー/Enthalpy[J/kg or J/kmol]と圧力/Pressure[Pa]の取得
print(gas1.HPX) #エンタルピー/Enthalpy[J/kg or J/kmol]と圧力/Pressure[Pa]とモル分率の取得
print(gas1.HPY) #エンタルピー/Enthalpy[J/kg or J/kmol]と圧力/Pressure[Pa]と質量分率の取得
print(gas1.P) #圧力/Pressure[Pa]の取得
# print(gas1.P_sat) #飽和圧/Saturation pressure[Pa]の取得
# print(gas1.Pe) #電子圧力/Electron Pressure[Pa]の取得
# print(gas1.S) #エントロピー/Entropy[J/K]の取得


print('')
print(gas1.SP) #エントロピー/Entropy[J/kg/K or J/kmol/K]と圧力/Pressure[Pa]の取得
print(gas1.SPX) #エントロピー/Entropy[J/kg/K or J/kmol/K]と圧力/Pressure[Pa]とモル分率の取得
print(gas1.SPY) #エントロピー/Entropy[J/kg/K or J/kmol/K]と圧力/Pressure[Pa]と質量分率の取得
print(gas1.SV) #エントロピー/Entropy[J/kg/K or J/kmol/K]と比体積/Specific volume[m^3/kg or m^3/kmol]の取得
print(gas1.SVX) #エントロピー/Entropy[J/kg/K or J/kmol/K]と比体積/Specific volume[m^3/kg or m^3/kmol]とモル分率の取得
print(gas1.SVY) #エントロピー/Entropy[J/kg/K or J/kmol/K]と比体積/Specific volume[m^3/kg or m^3/kmol]と質量分率の取得


print('')
print(gas1.T) #温度/Temperature[K]の取得
print(gas1.TD) #温度/Temperature[K]と密度/Density[kg/m^3 or kmol/m^3]の取得
print(gas1.TDX) #温度/Temperature[K]と密度/Density[kg/m^3 or kmol/m^3]とモル分率の取得
print(gas1.TDY) #温度/Temperature[K]と密度/Density[kg/m^3 or kmol/m^3]と質量分率の取得
print(gas1.TP) #温度/Temperature[K]と圧力/Pressure[Pa]の取得
print(gas1.TPX) #温度/Temperature[K]と圧力/Pressure[Pa]とモル分率の取得
print(gas1.TPY) #温度/Temperature[K]と圧力/Pressure[Pa]と質量分率の取得


print('')
# print(gas1.T_sat) #飽和温度/Saturation temperature[K]の取得 
print(gas1.Te) #電子温度/Electron Temperature[K]の取得
# print(gas1.U) #内部エネルギー/Internal energy[J]の取得
print(gas1.UV) #内部エネルギー/Internal energy[J/kg or J/kmol]と比体積/Specific volume[m^3/kg or m^3/kmol]の取得
print(gas1.UVX) #内部エネルギー/Internal energy[J/kg or J/kmol]と比体積/Specific volume[m^3/kg or m^3/kmol]とモル分率の取得
print(gas1.UVY) #内部エネルギー/Internal energy[J/kg or J/kmol]と比体積/Specific volume[m^3/kg or m^3/kmol]と質量分率の取得
# print(gas1.V) #体積/Volume[m^3]の取得
print(gas1.X) #モル分率の取得
print(gas1.Y) #質量分率の取得
[OUT]
False
(0.08189392763801234, 101325.00000000001)
(0.08189392763801234, 101325.00000000001, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
(0.08189392763801234, 101325.00000000001, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))

(26468.50456294105, 101325.00000000001)
(26468.50456294105, 101325.00000000001, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
(26468.50456294105, 101325.00000000001, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
101325.00000000001

(64910.063854689455, 101325.00000000001)
(64910.063854689455, 101325.00000000001, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
(64910.063854689455, 101325.00000000001, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
(64910.063854689455, 12.21091757157138)
(64910.063854689455, 12.21091757157138, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
(64910.063854689455, 12.21091757157138, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))

300.0
(300.0, 0.08189392763801234)
(300.0, 0.08189392763801234, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
(300.0, 0.08189392763801234, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
(300.0, 101325.00000000001)
(300.0, 101325.00000000001, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
(300.0, 101325.00000000001, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))

300.0
(-1210802.7183765294, 12.21091757157138)
(-1210802.7183765294, 12.21091757157138, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
(-1210802.7183765294, 12.21091757157138, array([1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0.]))
[1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0.]
[1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0.]

 特定の成分だけ抽出したい場合、Solutionオブジェクトはスライスのような機能があり、これを利用して値の抽出が可能です。

[IN]
Xmajor = gas1['CH4', 'O2', 'CO2', 'H2O', 'N2'].X #主要成分のモル分率の取得
print(Xmajor)
major = gas1['CH4', 'O2', 'CO2', 'H2O', 'N2'] #主要成分の質量分率の取得
print(major)
cp_major = major.partial_molar_cp #主要成分の部分モル定圧比熱の取得
print(cp_major)
wdot_major = major.net_production_rates #主要成分の生成速度の取得
print(wdot_major)
[OUT]
Xmajor = gas1['CH4', 'O2', 'CO2', 'H2O', 'N2'].X #主要成分のモル分率の取得
print(Xmajor)
major = gas1['CH4', 'O2', 'CO2', 'H2O', 'N2'] #主要成分の質量分率の取得
print(major)
cp_major = major.partial_molar_cp #主要成分の部分モル定圧比熱の取得
print(cp_major)
wdot_major = major.net_production_rates #主要成分の生成速度の取得
print(wdot_major)

 3-2-2.データの書き込み

 Solutionオブジェクトへの書き込みは属性に合わせてリストなどを渡すことで実装できます。

 下記事例ではSolution.TPで(温度, 圧力)の情報を更新しました。

[IN]
gas1.TP = 1200, 101325 #温度、圧力を書き換え
gas1()
[OUT]
  gri30:

       temperature   1200 K
          pressure   1.0133e+05 Pa
           density   0.020473 kg/m^3
  mean mol. weight   2.016 kg/kmol
   phase of matter   gas

                          1 kg             1 kmol     
                     ---------------   ---------------
          enthalpy        1.3295e+07        2.6802e+07  J
   internal energy        8.3457e+06        1.6825e+07  J
           entropy             85222        1.7181e+05  J/K
    Gibbs function       -8.8972e+07       -1.7937e+08  J
 heat capacity c_p             15377             31000  J/K
 heat capacity c_v             11253             22686  J/K

                      mass frac. Y      mole frac. X     chem. pot. / RT
                     ---------------   ---------------   ---------------
                H2                 1                 1           -17.978
     [  +52 minor]                 0                 0  

【モル分率Xの書き換え】
 次にGas成分を書き換えます。上記より元々は水素(H2)のみでしたが、組成とモル分率をXで書き換えました。
※後述の通りString(文字列)でも辞書型(dict)でもどちらも可能です

[IN]
print(gas1.X)

gas1.X = 'CH4:1, O2:2, N2:7.52' #ガス組成の設定 
print(gas1.X)

[OUT]
[1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0.]

[0.         0.         0.         0.19011407 0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.09505703 0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.
 0.         0.         0.         0.         0.         0.7148289
 0.         0.         0.         0.         0.        ]

 X属性の出力は各成分(53種類)のモル分率となります。下記に可視化した結果を示します。

[IN]
names_species = gas1.species_names #53種類の化学種の名前
df_species = pd.DataFrame({'化学種':names_species, 'モル分率X':gas1.X}) #DataFrameの作成
print(f'Xの合計:{df_species["モル分率X"].sum()}')
display(HorizontalDisplay(df_species.iloc[:18], df_species.iloc[18:36], df_species.iloc[36::]))

[OUT]
Xの合計:1.0

 モル分率の書き換えは辞書型でも対応できます。

[IN]
phi = 0.8 #余剰空気比の設定
gas1.X = {'CH4':1, 'O2':2/phi, 'N2':2*3.76/phi} #ガス組成の設定:3.76はN2のモル比(79/21=3.76)
print(gas1.X)

[OUT]

 Solution.Xに全て1を渡せば全成分のモル分率が等分され、Solution.Yに全て1を渡せば全成分の質量分率が等分されます。

[IN]
gas1.X = np.ones(53) #全ての成分のモル分率を1に設定
df_species = pd.DataFrame({'化学種':names_species, 'モル分率X':gas1.X, '質量分率Y':gas1.Y}) #DataFrameの作成
display(HorizontalDisplay(df_species.iloc[:18], df_species.iloc[18:36], df_species.iloc[36::]))

gas1.Y = np.ones(53) #全ての成分の質量分率を1に設定
df_species = pd.DataFrame({'化学種':names_species, 'モル分率X':gas1.X, '質量分率Y':gas1.Y}) #DataFrameの作成
display(HorizontalDisplay(df_species.iloc[:18], df_species.iloc[18:36], df_species.iloc[36::]))

[OUT]

4.Next Step:事例実施

 Canteraでは次のステップとして事例集があり、次のレベルを学習できます。

https://cantera.org/examples/index.html

 色々いけそうだったら時間かけて学習したいのですが、とりあえず今回は紹介だけとします。

5.Documents

 Canteraでは学習用に「Canteraライブラリ」と「科学関係の理論」に関するドキュメントが提供されています。

5-1.Canteraライブラリの詳細

 2024年12月現在でのCanteraに関するドキュメントは下記の通り

5-2.Input File

 Canteraを扱うには相や解きたい問題を表現したInput fileが必要です。Input Fileの作成方法に関するドキュメントは下記の通りです。

  • YAML Input File Users' Guide:Input fileのフォーマットに関するガイド

    • Phases and their InterfacesFor each phase that appears in a problem, a corresponding entry should be present in the input file(s). We'll start by describing the entries for phases of various types, and then look at how to define interfaces between phases.

    • Elements and SpeciesFor each species declared as part of a phase description, both the species and the elements that it is comprised of must be defined. Here, we describe how both are defined.

    • ReactionsCantera supports a number of different types of reactions, including several types of homogeneous reactions, surface reactions, and electrochemical reactions. For each, there is a corresponding entry type. Here, we describe how to declare each type of reaction and provide the necessary parameters to calculate the reaction rate for each.

  • YAML Input File API Reference:参考API

  • YAML Format TutorialYAMLファイルの構文や構造を説明

 新しいInput FIlesの変換/作成方法は下記の通りです。

5-3.科学の理論式

 Canteraではライブラリ内で使用される理論式やモデルの紹介をしています。

https://cantera.org/science/index.html

 各説明は下記の通りです。

  • 化学反応速度論(Chemical Kinetic Theory)

    1. ThermodynamicsThe theory behind how Cantera calculates species and phase thermodynamic properties.

    2. Kinetics and Reaction RatesThe models and equations that Cantera uses to calculate chemical reaction rates.

    3. TransportThe models that Cantera uses to calculate transport properties and rates.

  • 反応器およびフレームモデル(Cantera Reactor and Flame Models)

    • ReactorsCantera provides a range of generalized zero-dimensional models that can be given a range of initial and boundary conditions and can also be linked to form reactor networks.

    • FlamesCantera includes a set of models for representing steady-state, quasi-one-dimensional reacting flows, which can be used to simulate a number of common flames.

 Reactorsでは更に詳細に分類されています。

https://cantera.org/science/reactors/reactors.html


参考資料

Cantera公式

技術ブログ

https://www.rccm.co.jp/icem/rccm_fluid/posts/cantera_surfacechemistry1/

あとがき

 プロセスシミュレーションをPythonで実施する最適なライブラリを探しているけど、これがベストかわからないのでいったん中止。
これがよさそうなら追って追記

いいなと思ったら応援しよう!