calculate, equilibrium を使って熱力学計算¶
前章では熱力学データベースを使って状態図を計算する方法について学びました。このときにはまず、熱力学パラメータから正則溶体近似に基づいてギブスエネルギーを計算します。そして全ての相のギブスエネルギーを計算したのちに、系の平衡状態（無限に時間をおいた場合に到達する最もエネルギーが低い状態）における安定相を計算しています。
ここで Calculate を使うと、指定した相の物性値（ギブスエネルギーを含む）を計算することができます。
また equilibrium を使うと、複数の相が関わる平衡状態における物性値を計算することができます。
つまり binplot は Calculate と equilibrium をさまざまな温度、組成で行って状態図の計算を行っています。
ここでは、 Al-Zn ２元系を用いて Calculate と equilibrium の計算結果について説明します。
Al-Zn (S. Mey, 1993)¶
Al-Zn 状態図には、FCC相の溶解度ギャップが含まれています。
参考文献：S. an Mey, Zeitschrift für Metallkunde 84(7) (1993) 451-455.
必要なライブラリとデータベースを読み込みます。
[1]:
%matplotlib inline
from pycalphad import Database, calculate, equilibrium, variables as v
dbf = Database('alzn_mey.tdb')
comps = ['AL', 'ZN', 'VA']
次に考慮する相を選択します。 alzn_mey.tdbにどのような相が含まれているか確かめるには、メモ帳などのテキストエディタで alzn_mey.tdbを開き、PHASE で始まっている行を確認してください。LinuxかMacを使っている場合には、以下のようにしても相の情報が得られます。
[2]:
!grep PHASE alzn_mey.tdb
$ PARAMETERS FOR LIQUID PHASE
PHASE LIQUID % 1 1.0 !
PHASE FCC_A1 % 1 1.0 !
PHASE HCP_A3 % 1 1.0 !
alzn_mey.tdbにはLIQUID、FCC_A1、HCP_A3相が含まれています。これらの相をphases_alznとして取り込むには以下のようにします。
[3]:
phases_alzn = ['LIQUID', 'FCC_A1', 'HCP_A3']
xarray データセット¶
pycalphadにおいて、 calculate または equilibrium で行った計算結果は、xarray データセットとして返されます。xarray データセットは N 次元の表形式のデータ構造です。 これは、Pandas DataFrame の N 次元版とよく似ています。
このノートブックでは、pycalphad の xarray データセットの構造とその基本的な使用法について説明します。
xarray データセットの詳細なドキュメントとチュートリアルは、xarray ドキュメント を参照してください。
データ構造¶
各データセットには、計算条件と物性値が格納されます。 3 つの重要な用語があります。
Dimensions: 圧力(P)、温度(T)などがDimensions（次元）として用いられます。Coordinates:Dimensionsのラベルです。Data variables: ギブスエネルギー、混合エネルギー、組成など、pycalphad によって計算される物性値です。
calculate()¶
calculate は、指定した相の物性値を計算するために使用します。実行する際、4つの引数を指定する必要があります。
calculate(database, components, phases, Pressure, Temperature)`
例えば１気圧（1013125 Pa）、\(T\) = 500 K ~ 1100 Kの間を100 Kステップで計算するには以下のようにします。
calculate(dbf, comps, phase, P=101325, T=(500, 1100, 100))
計算時に5つのdimensions/coordinatesをもつxarrayデータセットがつくられます。
component: 構成元素N: 全原子数 (mol).（1に固定されるため数えない）P: 圧力 (Pa).T: 温度 (K).points: 計算点。各内部自由度のすべてのインデックスを表します。calculateはデフォルトですべての内部自由度にわたって計算点をサンプリングします。pointsには物理的な意味や順序はありません。internal_dof: 内部自由度 (inner_dof) は、任意の相のサイトフラクション配列のインデックスです。FCC_A1相の副格子モデル(AL、ZN)の場合には、internal_dof は整数 0 と 1（ALサイトは0、ZN サイトは1)です。
Calculate によって出力される計算結果には 4 つのデータがあります。
Phase:相。必ず含まれます。X: モル分率で表した各成分の組成。温度、圧力、ポイントの関数。Y:inner_dof配列内の各インデックスのサイトフラクション。温度、圧力、ポイントのインデックスの関数。output: 計算される物性値。 デフォルトはモルギブズエネルギーGM。
Calculate を使って、FCC_A1相の物性値を計算するには以下のようにします。
[4]:
calc_result = calculate(dbf, comps, 'FCC_A1', P=101325, T=(500, 1100, 100))
print(calc_result)
<xarray.Dataset>
Dimensions: (N: 1, P: 1, T: 6, points: 4001, component: 2, internal_dof: 2)
Coordinates:
* component (component) <U2 'AL' 'ZN'
* N (N) float64 1.0
* P (P) float64 1.013e+05
* T (T) float64 500.0 600.0 700.0 800.0 900.0 1e+03
Dimensions without coordinates: points, internal_dof
Data variables:
X (N, P, T, points, component) float64 1.0 1e-14 ... 0.7439 0.2561
Phase (N, P, T, points) <U6 'FCC_A1' 'FCC_A1' ... 'FCC_A1' 'FCC_A1'
Y (N, P, T, points, internal_dof) float64 1.0 1e-14 ... 0.2561
GM (N, P, T, points) float64 -1.559e+04 -2.01e+04 ... -4.81e+04
Attributes:
phase_indices: {'FCC_A1': slice(0, 4001, None)}
/Users/iikubo/miniconda3/envs/pycalphad/lib/python3.9/site-packages/pycalphad/core/utils.py:54: RuntimeWarning: invalid value encountered in divide
pts[:, cur_idx:end_idx] /= pts[:, cur_idx:end_idx].sum(axis=1)[:, None]
計算結果を操作する¶
計算結果のxarrayデータセット（ここでは calc_result ）を操作するには、Pandas 配列と同様に、selを使用して座標値でデータを選択するか、iselを使用して座標のインデックスでデータを選択します。適切にデータ変数を選択することで、プロットや計算の実行に役立つ xarray データセットを得ることができます (DataArrays vs Datasets を参照)。
例えば、calc_result における50 番目のポイント (インデックスで選択) の 500 K (値で選択) における ZN (internal_dof = 1) のサイトフラクションを取得するには以下のようにします。
[5]:
print(calc_result.Y.isel(internal_dof=1, points=49).sel(T=500))
<xarray.DataArray 'Y' (N: 1, P: 1)>
array([[0.97648824]])
Coordinates:
* N (N) float64 1.0
* P (P) float64 1.013e+05
T float64 500.0
DataArray の values 属性にアクセスすると、多次元 NumPy 配列が返されます
[6]:
print(calc_result.Y.isel(internal_dof=1, points=49).sel(T=500).values)
[[0.97648824]]
[7]:
print(calc_result.X.sel(component='ZN').sel(T=500).values)
[[[1.00000000e-14 1.00000000e+00 1.00000000e+00 ... 8.43004350e-01
8.87927218e-01 2.56066359e-01]]]
[8]:
print(calc_result.GM.sel(T=500).values)
[[[-15589.38872794 -20099.46240961 -20099.46240961 ... -20683.03658045
-20722.47948049 -17338.16556129]]]
計算によって得られたギブスエネルギーを図示します。
[14]:
import matplotlib.pyplot as plt
plt.scatter(calc_result.X.sel(component='ZN'), calc_result.GM, marker='.')
[14]:
<matplotlib.collections.PathCollection at 0x14ee18280>
軸タイトルや凡例は以下のようにつけます。
[17]:
import matplotlib.pyplot as plt
from pycalphad.plot.utils import phase_legend
# Get the colors that map phase names to colors in the legend
legend_handles, color_dict = phase_legend(phases_alzn)
fig = plt.figure(figsize=(6,4))
ax = fig.gca()
ax.scatter(calc_result.X.sel(component='ZN'), calc_result.GM, marker='.')
# Format the plot
ax.set_xlabel('X(Zn)')
ax.set_ylabel('GM (J/mol)')
ax.set_xlim((0, 1))
ax.set_ylim((-60000, 0))
ax.legend(handles=legend_handles, loc='upper right')
plt.show()
T = 500 Kで３つの相（LIQUID、FCC_A1、HCP_A3）のギブスエネルギーを比較するには以下のようにします。
[21]:
import matplotlib.pyplot as plt
from pycalphad.plot.utils import phase_legend
# Get the colors that map phase names to colors in the legend
legend_handles, color_dict = phase_legend(phases_alzn)
fig = plt.figure(figsize=(6,4))
ax = fig.gca()
# Loop over phases, calculate the Gibbs energy, and scatter plot GM vs. X(RE)
for phase_name in phases_alzn:
result = calculate(dbf, comps, phase_name, P=101325, T=500, output='GM')
ax.scatter(result.X.sel(component='ZN'), result.GM, marker='.', s=5, color=color_dict[phase_name])
# Format the plot
ax.set_xlabel('X(Zn)')
ax.set_ylabel('GM (J/mol)')
ax.set_xlim((0, 1))
ax.legend(handles=legend_handles, loc='upper right')
plt.show()
/Users/iikubo/miniconda3/envs/pycalphad/lib/python3.9/site-packages/pycalphad/core/utils.py:54: RuntimeWarning: invalid value encountered in divide
pts[:, cur_idx:end_idx] /= pts[:, cur_idx:end_idx].sum(axis=1)[:, None]
/Users/iikubo/miniconda3/envs/pycalphad/lib/python3.9/site-packages/pycalphad/core/utils.py:54: RuntimeWarning: invalid value encountered in divide
pts[:, cur_idx:end_idx] /= pts[:, cur_idx:end_idx].sum(axis=1)[:, None]
/Users/iikubo/miniconda3/envs/pycalphad/lib/python3.9/site-packages/pycalphad/core/utils.py:54: RuntimeWarning: invalid value encountered in divide
pts[:, cur_idx:end_idx] /= pts[:, cur_idx:end_idx].sum(axis=1)[:, None]
equilibrium()¶
また equilibrium を使うと、複数の相が関わる平衡状態における物性値を計算することができます。
equilibrium は、複数の相が関わる平衡状態における物性値を計算するために使用します。実行する際、4つの引数を指定する必要があります。
equilibrium(database, components, phases, conditions, output)`
例えば１気圧（1013125 Pa）、\(T\) = 500 K ~ 1000 Kの間を100 KステップでHM（モルエンタルピー）を計算するには以下のようにします。
equilibrium(dbf, comps , phases, {v.X('ZN'):(0,1,0.05), v.T: (500, 1000, 100), v.P:101325}, output='HM')
[22]:
phases = ['LIQUID', 'FCC_A1', 'HCP_A3']
eq_result = equilibrium(dbf, comps , phases, {v.X('ZN'):(0,1,0.05), v.T: (500, 1000, 100), v.P:101325}, output='HM')
print(eq_result)
<xarray.Dataset>
Dimensions: (N: 1, P: 1, T: 5, X_ZN: 20, vertex: 3, component: 2,
internal_dof: 2)
Coordinates:
* N (N) float64 1.0
* P (P) float64 1.013e+05
* T (T) float64 500.0 600.0 700.0 800.0 900.0
* X_ZN (X_ZN) float64 1e-10 0.05 0.1 0.15 0.2 ... 0.75 0.8 0.85 0.9 0.95
* vertex (vertex) int64 0 1 2
* component (component) <U2 'AL' 'ZN'
Dimensions without coordinates: internal_dof
Data variables:
NP (N, P, T, X_ZN, vertex) float64 1.0 nan nan 1.0 ... 1.0 nan nan
GM (N, P, T, X_ZN) float64 -1.559e+04 -1.615e+04 ... -5.068e+04
MU (N, P, T, X_ZN, component) float64 -1.559e+04 ... -5.065e+04
X (N, P, T, X_ZN, vertex, component) float64 1.0 1e-10 ... nan nan
Y (N, P, T, X_ZN, vertex, internal_dof) float64 1.0 1e-10 ... nan
Phase (N, P, T, X_ZN, vertex) <U6 'FCC_A1' '' '' ... 'LIQUID' '' ''
HM (N, P, T, X_ZN) float64 5.194e+03 5.859e+03 ... 2.528e+04
Attributes:
engine: pycalphad 0.10.2
created: 2024-05-21T06:28:59.183045
equilibrium() データセットには、6 つの dimensions/coordinates があります。
N: 全原子数 (mol).（1に固定されるため数えない）P: 圧力 (Pa)。T: 温度 (K)。X_ZN: これは上のconditionsに渡された元素の組成です。v.X('ZN')としたのでX_ZNになります。vertex:vertex(頂点)は、平衡状態に含まれる相のインデックスです。vertex自体には固有の物理的な意味はありません。平衡状態に含まれる相の数を記述するのに十分な数が自動的生成されます。これはギブスの相律に反するほど大きくなることはありません。component: 系に含まれる元素。internal_dof:inner_dof(内部自由度) は、任意の相のサイトフラクション配列内のインデックスです。FCC_A1相の副格子モデル(AL、ZN)の場合には、internal_dof は整数 0 と 1（ALサイト(0)とZN サイト(1)です。
equilibrium() によって返されるデータセットは calculate 計算に非常に似ており、少なくとも 6 つのデータ変数があります。
Phase: 与えられた条件における平衡状態の相を表す文字列。Phaseの数はlen(vertex)個あります。vertexで記述されるインデックスよりも平衡状態にある相の数が少ない場合は常に、相の値は「』』」で埋められます。 FCC_A1 の単相領域の場合、相の値は['FCC_A1', '']になります。複数のフェーズが存在する場合、それらは必ずしもソートされていないことに注意してください。NP: 平衡状態にある各相の相分率。他に平衡相がない場合 (例: 単相['FCC_A1', ''])、相が存在しないためNPの値は 0 ではなくnanになります。MU: 計算された条件における各成分の化学ポテンシャル。X: 計算された条件における各相の各元素の平衡組成。Y: 計算された条件のにおける各相の各サイトの占有率。GM:calculateのoutputと同じです。outputの値に関係なく、常に出力されます。output: (optional)outputは、equilibriumに渡された出力キーワードによって常に計算される平衡状態の物性値（プロパティ）です。calculateとは異なり、GMは常に報告されるため、outputはGMに追加されます。
例として、\(T\) = 800 K におけるHCP_A3 相の相分率を組成の関数として求めてみましょう。
相分率の値にアクセスするには、coodinates（座標）のインデックス、または値が必要になります。
[28]:
print(eq_result.NP.sel(Phase='FCC_A1', T=800))
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[28], line 1
----> 1 print(eq_result.NP.sel(Phase='FCC_A1', T=800))
File ~/miniconda3/envs/pycalphad/lib/python3.9/site-packages/xarray/core/dataarray.py:1583, in DataArray.sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
1473 def sel(
1474 self,
1475 indexers: Mapping[Any, Any] | None = None,
(...)
1479 **indexers_kwargs: Any,
1480 ) -> Self:
1481 """Return a new DataArray whose data is given by selecting index
1482 labels along the specified dimension(s).
1483
(...)
1581 Dimensions without coordinates: points
1582 """
-> 1583 ds = self._to_temp_dataset().sel(
1584 indexers=indexers,
1585 drop=drop,
1586 method=method,
1587 tolerance=tolerance,
1588 **indexers_kwargs,
1589 )
1590 return self._from_temp_dataset(ds)
File ~/miniconda3/envs/pycalphad/lib/python3.9/site-packages/xarray/core/dataset.py:3033, in Dataset.sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
2972 """Returns a new dataset with each array indexed by tick labels
2973 along the specified dimension(s).
2974
(...)
3030 DataArray.sel
3031 """
3032 indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel")
-> 3033 query_results = map_index_queries(
3034 self, indexers=indexers, method=method, tolerance=tolerance
3035 )
3037 if drop:
3038 no_scalar_variables = {}
File ~/miniconda3/envs/pycalphad/lib/python3.9/site-packages/xarray/core/indexing.py:185, in map_index_queries(obj, indexers, method, tolerance, **indexers_kwargs)
182 options = {"method": method, "tolerance": tolerance}
184 indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "map_index_queries")
--> 185 grouped_indexers = group_indexers_by_index(obj, indexers, options)
187 results = []
188 for index, labels in grouped_indexers:
File ~/miniconda3/envs/pycalphad/lib/python3.9/site-packages/xarray/core/indexing.py:146, in group_indexers_by_index(obj, indexers, options)
144 raise KeyError(f"no index found for coordinate {key!r}")
145 elif key not in obj.dims:
--> 146 raise KeyError(
147 f"{key!r} is not a valid dimension or coordinate for "
148 f"{obj.__class__.__name__} with dimensions {obj.dims!r}"
149 )
150 elif len(options):
151 raise ValueError(
152 f"cannot supply selection options {options!r} for dimension {key!r}"
153 "that has no associated coordinate or index"
154 )
KeyError: "'Phase' is not a valid dimension or coordinate for Dataset with dimensions Frozen({'N': 1, 'P': 1, 'T': 5, 'X_ZN': 20, 'vertex': 3})"
Phaseは dimension/coordinate ではなくデータ変数であるため、上のような指定の仕方をするとエラーが起きます。
sel または isel コマンドを使用するには、HCP_A3 相がデータ変数として格納されているセルのインデックスを事前に調べておく必要があります。
これをひとつひとつ調べることは面倒な作業であるが、xarrayに備わっているマスキング masking を利用すれば、条件に一致するデータ値を見つけることができます。
以下の例では、データ変数にFCC_A1相が格納されているセルの相分率を表示しています。
[29]:
print(eq_result.NP.where(eq_result.Phase=='FCC_A1'))
<xarray.DataArray 'NP' (N: 1, P: 1, T: 5, X_ZN: 20, vertex: 3)>
array([[[[[1. , nan, nan],
[1. , nan, nan],
[0.97607896, nan, nan],
[0.92129859, nan, nan],
[0.86651821, nan, nan],
[0.81173783, nan, nan],
[0.75695745, nan, nan],
[0.70217707, nan, nan],
[0.64739669, nan, nan],
[0.59261631, nan, nan],
[0.53783594, nan, nan],
[0.48305556, nan, nan],
[0.42827518, nan, nan],
[0.3734948 , nan, nan],
[0.31871442, nan, nan],
[0.26393404, nan, nan],
[0.20915367, nan, nan],
[0.15437329, nan, nan],
[0.09959291, nan, nan],
[0.04481253, nan, nan]],
...
[[1. , nan, nan],
[0.75935221, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan]]]]])
Coordinates:
* N (N) float64 1.0
* P (P) float64 1.013e+05
* T (T) float64 500.0 600.0 700.0 800.0 900.0
* X_ZN (X_ZN) float64 1e-10 0.05 0.1 0.15 0.2 ... 0.75 0.8 0.85 0.9 0.95
* vertex (vertex) int64 0 1 2
[23]:
print(eq_result.NP.where(eq_result.Phase=='FCC_A1').sel(P=101325, T=800))
<xarray.DataArray 'NP' (N: 1, X_ZN: 20, vertex: 3)>
array([[[1. , nan, nan],
[1. , nan, nan],
[1. , nan, nan],
[1. , nan, nan],
[0.89739922, nan, nan],
[0.71825009, nan, nan],
[0.53910096, nan, nan],
[0.35995184, nan, nan],
[0.18080271, nan, nan],
[ nan, 0.00165358, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan],
[ nan, nan, nan]]])
Coordinates:
* N (N) float64 1.0
P float64 1.013e+05
T float64 800.0
* X_ZN (X_ZN) float64 1e-10 0.05 0.1 0.15 0.2 ... 0.75 0.8 0.85 0.9 0.95
* vertex (vertex) int64 0 1 2
\(T\) = 800 Kにおいて、各相の相分率がZn濃度（X(ZN)）に対して変化する図を作ってください。
[9]:
%matplotlib inline
import matplotlib.pyplot as plt
from pycalphad import Database, calculate, equilibrium, variables as v
from pycalphad.plot.utils import phase_legend
dbf = Database('alzn_mey.tdb')
comps = ['AL', 'ZN', 'VA']
phases = ['LIQUID', 'FCC_A1', 'HCP_A3']
eq_result = equilibrium(dbf, comps , phases, {v.X('ZN'):(0,1,0.05), v.T: (500, 1000, 100), v.P:101325}, output='HM')
# Get the colors that map phase names to colors in the legend
legend_handles, color_dict = phase_legend(phases)
fig = plt.figure(figsize=(6,4))
ax = fig.gca()
ax.scatter(eq_result.X.sel(component='ZN').sel(P=101325, T=800), eq_result.NP.where(eq_result.Phase=='LIQUID').sel(P=101325, T=800), marker='.')
ax.scatter(eq_result.X.sel(component='ZN').sel(P=101325, T=800), eq_result.NP.where(eq_result.Phase=='FCC_A1').sel(P=101325, T=800), marker='.')
ax.scatter(eq_result.X.sel(component='ZN').sel(P=101325, T=800), eq_result.NP.where(eq_result.Phase=='HCP_A3').sel(P=101325, T=800), marker='.')
#plt.scatter(eq_result.X.sel(component='ZN').sel(P=101325, T=800), eq_result.NP.where(eq_result.Phase=='FCC_A1').sel(P=101325, T=800), marker='.')
#plt.scatter(eq_result.X.sel(component='ZN').sel(P=101325, T=800), eq_result.NP.where(eq_result.Phase=='LIQUID').sel(P=101325, T=800), marker='.')
#plt.scatter(eq_result.X.sel(component='ZN').sel(P=101325, T=800), eq_result.NP.where(eq_result.Phase=='HCP_A3').sel(P=101325, T=800), marker='.')
# Format the plot
ax.set_xlabel('X(Zn)')
ax.set_ylabel('Phase fraction')
ax.set_xlim((0, 1))
#ax.set_ylim((-60000, 0))
ax.legend(handles=legend_handles, loc='lower right')
plt.show()
課題¶
Moodle上で提出する。
JupyterNotebook形式(*.ipynb)で提出する。
学籍番号、氏名を忘れずに記載する。
コードを実装するにあたって考えた過程や工夫したポイント、動作確認の方法等を記述してあれば評価時に考慮します。
moodleで配布した講義資料内にはAl-Mg_Zhong.tdb、NI_AL_DUPIN_2001.TDB、alfe_sei.TDB、nbre_liu.tdbが含まれている。
Al-Mg_Zhong.tdb：Al-Mg２元系の熱力学データベース
NI_AL_DUPIN_2001.TDB：Ni-Al２元系の熱力学データベース
alfe_sei.TDB：Al-Fe２元系の熱力学データベース
nbre_liu.tdb：Nb-Re２元系の熱力学データベース
Al-Mg_Zhong.tdbを読み込み、\(T\) = 800 Kにおいて、各相の相分率がMg濃度（X(MG)）に対して変化する図を作ってください。得られた計算結果が状態図と対応していることを説明してください。
[ ]: