⽔・蒸気性質シンポジウム2023

2023.11.29
慶應義塾⼤学 三⽥キャンパス

分⼦シミュレーションで貢献する地球
温暖化係数が極めて低い冷媒（熱エネ
ルギー輸送媒体）の開発● 近藤 智恵⼦
今井 友暁
川原 豪真
奥村 哲也

⻑崎⼤学
⻑崎⼤学
⻑崎⼤学
⻑崎⼤学

1

冷媒とは．．．熱エネルギーを輸送する有機物質
OutdoorCompressor

Heat
Exchanger

Cool air

15 C

25 C Heat
Exchanger
Heat
エアコンに使われている
R32の分⼦構造

80C

IndoorExpansion Valve45C

⽣体試料⽤など －70C

【冷媒 Refrigerant】
熱を輸送するガス

漁獲物⽤など －60C

蒸気⇔液に変化させて潜熱でエネルギーを輸送

半導体製造など －80C

出典︓朝⽇新聞
特集築地

出典︓NHK

出典︓Forbes India
2021 June. 282冷媒選択が難しい理由
GWP︓地球温暖化係数(Global Warming Potential）⼤気中に放出された単位質量当たりの当該物質が地球温暖化に与える影響を、CO2が１として相対的に表した係数（積分年数︓20年、100年、500年の値が公開されているが，通常GWP100を⽤いる）

単位質量当たり放射強制⼒（⾚外光吸収率）と、分解され減衰してゆく⼤気中濃度の積算
TH RFd
RF: Radiative forcing due to a pulse release

GWP  TH

0

 RF

CO 2

d

of 1 kg trace compound
TH: Chosen time horizon (100yr)

0

GWP of HFCs
HFC134a

1300

HFC410A

2088

HFC404A

3922

HFC23

14600

HFC32

675

CO2濃度上昇のため IPCC 6ARでは値が修正された．

オゾン層破壊係数（塩素を使⽤しない）
消費電⼒
懸案 ユニット体積
事 燃焼性
項
急性毒性
PFAS

(COPを下げない)(体積能⼒を下げない)

PFAS/Regulation (EC) No 1907/2006
CF3基を持つ分⼦は、⼀部分解してTFAとなり、⽔に⾮常に溶けやすく、発がん性がある。TFAは⾃然界に存在する物質であり、冷媒がどの程度の影響を与えるかを正しく⾒積もることは難しい。

(微燃は許容される)(⼼筋の炎症など)

TFA (trifluoro acetic acid)︓
C2HF3O2 or CF3COOH（TFAを⽣成しない）

3

適⽤できる温度範囲は物質ごとに異なる（超低温） 冷凍機器
空調機器（-100）‐80 〜 -40℃ -5〜35℃

⾃動⾞A/C-5〜45℃

⼩規模発電⽤ORC
100〜200℃

現在使⽤中の

HFCs
HFC23

HFC32

HFC134a

HFC245fa

GWP=14600 *Lifetime 228yr

GWP=771 *Lifetime 5.2yr

GWP=1300 *Lifetime 13.4yr

GWP=858 *Lifetime 7.7yr

HFO1123

HFO1234yf

HFO1336mzz(Z)GWP＝0.005 *Lifetime 0.4days

GWP＝0.501 *Lifetime 12days

新冷媒候補

HFOs
Lifetime is very short
C=C reacts with OHin atmosphere.

＊IPCC_AR6

HFO1132a

GWP ＝0.252*
Lifetime 4.7days

超低温域熱物性？

安定な物質との混合？

GWP ＝2.08 *Lifetime 27days

PFAS規制により欧州
では使⽤できない．

4

ひとつの物質の調査に必要な時間は・・・2〜5年︖
表⾯張⼒測定臨界パラメータ
飽和密度・圧⼒
気液平衡
サンプル供給0.5〜1.5年

⾼精度測定（国内国際連携で分担）

1〜2年
低温の表⾯張⼒が測定できるけど，低温の飽和密度データがまだ無いのよね

⾼精度状態⽅程式、
最適化された混合モデルで信頼性の⾼い熱物性値

冷却能力・消費電力を計算．
間接的CO２排出量低減は必ず検討に⼊れる必要がある。

0.5年（随時）

0.5〜2年

状態⽅程式開発に必要な物性情報
 理想気体⽐熱 第2，第3ビリアル係数
 飽和密度・飽和圧⼒
 気液相平衡（混合系）

2050年
カーボンニュートラル実現（⽇本政府表明）

さらに，
混合パラメータの調査は膨⼤に

それでも遅すぎるという研究者の意⾒があるなか．．．

臨界付近の測定値を参照した分⼦シミュレーションにより，暫定的でもいい
ので低温域の物性を迅速に得たい．網羅が難しい混合系の調査も加速したい．

分⼦シミュレーション

5と ひと⼝で⾔っても．．．

量⼦化学計算（Quantum Mechanics Simulation）
⼿法︓ シュレディンガー⽅程式を解く分⼦軌道法（ab initio)
コーンシャム⽅程式を解く密度汎関数法(DFT)
・構造最適化・スキャニング・エネルギー計算
・振動計算

分⼦動⼒学計算（Molecular Dynamics Simulation）
⼿法︓(古典⼒学に沿って) 時間発展的に計算を⾏う
・液・蒸気単相 （⾮理想気体分）
・気液平衡系(2成分，3成分)
・液滴などの挙動
・⾼分⼦やたんぱく質の構造

地球温暖化係数理想気体⽐熱分⼦⼒場構築

⾮理想気体分⽐熱
飽和圧⼒・密度
気液平衡データ

その他，
モンテカルロ法 (Monte Carlo Simulation)など
QM-MM法6

QM（１）︓地球温暖化係数
TH RFd

GWP  TH 0
 RFCO 2 d0

HFC23
GWP=14600 *Lifetime 228yr

RF: Radiative forcing due to a pulse releaseof 1 kg trace compound
量⼦⼒学計算（QM）を⾏って推定が可能．
ヘテロクラスターでなくても，そこそこ
実験値と⼀致するようだ．

TH: Chosen time horizon (100yr)
HFC32

⼤気圏あるいは成層圏での分解速度が必要．
現在はチャンバー試験が主流︖
簡易推算式が提案されているが，QMでより正確な推算も可能となるかもしれない.

GWP=771 *Lifetime 5.2yr

図R23,R32モノマーの振動計算結果
CCSD/DGDZVPレベル
7

QM（１）︓地球温暖化係数
Nimbus4 で観測されたサハラ上空の地球放射強度＊．
Hanel et al., Appl. Opt., 10, 1376 (1971)
地球→宇宙⽅向の放射．これに対し
太陽から地球に届く放射を太陽放射
という．

Atm
Window『 Atm Window * 』と冷媒の⾚外吸収スペクトルが⼀致．
＊⽔蒸気やCO2に邪魔されない
地球放射強度域
これがGWPが⾼い要因の⼀つ

候補物質[R1123，R1132a，R1132(E)]ではR1132aのGWPが⾼そうだ．
やはり，成層圏での分解速度が分からなければ予測は難しい．．．

図CCSD/DGDZVPレベルの振動計算結果

8

QM（2）︓理想気体⽐熱 cv0
HFO1132a
GWP ＝0.252*
Lifetime 4.7days

HFC32
GWP=771 *Lifetime 5.2yr
2

Akasaka et al., I. J. Refrig. (2020)

15R

EOS(Tillner-Roth and Yokozeki, 1997)
B3LYP/6-311G(d,p) Ultra Fine
Scaling Factor0.989
Scaling Factor1.0
Scaling Factor1.1
Scaling Factor0.9
measured(Hozumi et al. 1994)

R1132aUltra Fine
1.5 B3LYP/6-311G(d,p)

Scaling Factor0.989
Scaling Factor1.0
Scaling Factor1.1
Scaling Factor0.9-1-1

Cv [kJ kg K ]

R32

1

1

1

0

Cv0 [kJ kg-1 K-1]

2

measured (Demirdesen et al., 2023)
Plank-Einstein formula

0.8
0.6

0.5

3R

200
50

100

500

1000

5000

10000

50

temperature[K]

100

500

300

1000

400 500
5000

10000

temperature [K]R32では，Scaling Factor 0.989の計算結果が最も⼀致した．
R1132aでは，Scaling Factorは1.0〜1.1であったことから，実⽤温度範囲での推算精度はおおよそ10％のようだ．
R32
15Rに漸近しなかった．分⼦振動が激しい⾼温域では、⾮調和性に起因する計算誤差が⽣じたと考えられる．
R1132a⼆重結合を有するHFO系物質は（ほぼ“ねじり”が無い）、⾼温域でも⽐較的⾮調和振動による影響が出にくいと考えられるが、理由は明確ではない

9

QM（2）︓⼒場構築（分⼦の設計図 結合⻑など）
1
1
2
2
U bonded   kr  r  r0    k    0    k 1  cos  n  0  
2
2
U non-bonded    12   6 ijij
  4 ij   
 
4 0 rijr rij 
ij  


結合⻑

-7.2588

平衡結合⻑
とバネ定数

qi q j-7.2590結合⾓-7.25930

r0(c) R1132a Bonding-7.25932

C=C0(d) R1132a Angle∠F-C-C-7.25934-7.2594-7.2590-7.2591-5-7.2592
MO Method(CCSD/DGTZVP)Fitted curve
AMBER03-7.2593-7.2594

C-H-7.2591Potential energy ×10

Potential energy ×10

C-Fr-7.25936-7.25933-1[kJ mol ]-1[kJ mol ]構造最適化↓
スキャニング-5-7.2592∠H-C-C-7.25934-7.25935-7.25936-7.25930∠F-C-F-7.25932
MO Method(CCSD/DGTZVP)Fitted curve
AMBER03-7.25934-7.25936-7.2592-7.25934-7.2593-7.25935∠H-C-H-7.2594
0.09

0.1

0.11

0.12

0.13

Bond length, r [nm]

図

0.14

0.15

-7.25936
105

110

115

120

125

130

Angle,  [deg ]CCSD/DGTZVPレベルの結合エネルギ―計算結果10

QM（2）︓⼒場構築（分⼦の設計図 点電荷）
1
1
2
2
U bonded   kr  r  r0    k    0    k 1  cos  n  0  
22
  12   6 ijij
  4 ij 
U non-bonded    
 
4 0 rijr rij 
ij  qi q j

クーロン⼒
点電荷の割当てR1123

エネルギー計算（ESPデータ取得）↓
RESP電荷の取得（ChelpG電荷）

1C0.4358e2C-0.2155e3H0.2358e

4F-0.1495e

5F-0.1872e

6F-0.1193eμ1.68

F: 電⼦吸引基 (Negative)
H: 電⼦供与基 (Positive)

双極⼦モーメント
測定値 (Kano, 2020)
1.33 Debye

R1132a1C-0.9421e2C

0.6942e3H

0.3389e4H

0.3389e5F-0.2150e

6F-0.2150e

μ1.6515液相での，誘起双極⼦を考慮に⼊れると，理想気体状態の永久双極⼦より20〜30％程⾼い⽅が良いとされている．
しかし，その程度を⾒積もることは難しい．
11

QM（2）︓⼒場構築（分⼦の設計図 LJパラメータ）
10

1
1
2
2
U bonded   kr  r  r0    k    0    k 1  cos  n  0  
2
2

(a)
5

  12   6 ij
ij
U non-bonded    4 ij   
 
4 0 rijr rij 
ij  vdW 項qi q jヘテロダイマーでエネルギー計算↓
修正が必要な原⼦と⽅向性を特定し，後はMDで微調整【General FF】

AMBER

CHARMM
CROMOSOPLSCOMPASS

SPC, SPC/E, TIP3P, TIP4P【Database】
LigParGen(proteins*)(gas-liquid, proteins)(gas-liquid, proteins)(gas-liquid)(gas-liquid)

TraPPE（water）-5-1

Inter molecular energy, U [kJ mol ]Lennard-Jonesパラメータ   の調整

0

-10
10
OPLS-AA/L
tuned in MD
MO method (BSSE corrected)
MP2/aug-cc-pVDZ
MO method (BSSE corrected)
CCSD/DGDZVP

(b)
5
0
-5
-10
10

(c)
5
0
-5
-10
0.3

0.7

0.8

distance between centers of gravity [nm]

0.4

0.5

0.6

12

MD（1）︓気液平衡モデル

 Zero-gravitational condition
 ”Periodic boundary” isapplied to each surface【Cal. conditions】
Software【To get “Sat. Densities”】

Gromacs2016
Double-precision

Density distribution averaged over 1-4ns.

Number of molecules 2000~6000
LJ and Coulomb
cut off / switch

1500

230 K
280 K
330 K
360 K

1.5 nm / 1.4 nm

Control method

density [kg m-3]

Temp. coupling
Nose-Hoover
Coulomb cal. method Particle Mesh Ewald
NVT 4~6 ns

Time integration

Leapfrog

Time step

0.25 ~ 0.5 fs

Sampling time1 ~ 2 ns

1000

Fitting correlation: ( z) 


500

0

5
10
height [nm]    2

    2  z  z0  tanh 
 D 

2D︓Interfacial thickness
Defined as 7 - 93% of the sat.
Density difference.
13

MD（1）︓気液共存曲線
R1123， R1132a

150
100
0

HFO1123
EOS (Akasaka, 2020)
GEMC (Raabe, 2016)
MD (tuned LJ parameter)
measured (Fukushima et al. , 2015)
measured (Higashi and Akasaka, 2016)
HFO1132a
EOS (Low, 2018)
Proposed corrrelation
MD (tuned LJ parameter)
measured (Otto and Thomas, 1964)
measured (Low, 2018)measured (Higashi and Perera, 2022)

500

300
250
200
150

1000
-3

density [kg m ]1500

Propanea
Tcrit=369.82 Ka
Ttp = 85.46 KaLemmon, 2009

350

temperature [K]fine-tuning
region

250
200

400

R1123 Ttrp = 105 K
R1132a Ttrp = 110 K

300

temperature[K]Propane1000

EOS (Lemmon, 2009)
MD simulation (202309)
measured data
Maass, 1921
Luo,1981
Kahre, 1973
Diller, 1982
Rodosevich, 1973 Kaminishi, 1988
McClune, 1976
Abdulagatov, 1995
Haynes, 1977
Glos, 2004

200

400

-3

density [kg m ]

600

800

LJパラメータの最終調整．
OPLS-AA/Lから出発し他の物質で提案されている⽂献値を参照（例えばRaabeら）．
最終的に気液共存曲線が測定値を⼀致するよう決定． （飽和液密度の偏差 <2％）
新冷媒の測定値は，Propaneよりも随分と少ない．特に超低温域（〜150 K）は実⽤範囲であるにもかかわらず密度の測定値が無い．状態⽅程式開発が難しい理由の⼀つ．

14

MD（2）︓飽和圧⼒

表⾯張⼒

飽和圧⼒ (R1123,R1132a)

表⾯張⼒: R1132(E)saturation pressure [MPa]4.0
3.5 HFO1123
3.0
2.5
2.0
1.5

EOS (Akasaka, 2021)
MD (tuned LJ parameter)
measured (Higashi and Akasaka, 2016)

HFO1132a
EOS (Low, 2018)
MD (tuned LJ parameter)
measured (Higashi, 2022)
measured (Tommasetti et al, 2021)
measured (Kayukawa, 2022)

1.0
0.5
0

150

200

250

temperature[K]300

ただし，cutoff以遠補正無し．Irving-Kirkwood
補正を⼊れると約10％⾼い値に．

R1132aは130 Kまで測定値が存在する（Kayukawa, Tomassettiら）

   Pzz Pxx  Pyy  Lz
   tail22


気液共存曲線が⼀致すると．．．飽和圧⼒，表⾯張⼒も良く⼀致する．
今のところ．偶然かもしれないが．

15

MD（3）︓ビリアル係数(圧⼒の) 第2 第3ビリアル係数-1

B2V×10 [m kg ]

1.0

0

3

0
-0.5
-1.5

250 K

-2.0
-2.5

0.5

0.6

0.7

0.8

center of mass, r [nm]0.9

稀薄蒸気のPvT計算を⾏い動径分布関数 g(r)からB2vを算出

B2V 

R1123
EOS (Akasaka et al.)Measurement (Kano)
MD
R1132aMD-10

R1132a
NVT 8~10 ns

200 K

0.4

-5

3

1000 K
600 K
400 K
300 K

-1.0


Z
 2 N A  1  g  r  r 2 drr1   0

1.0

C3V×106 [m6 kg-2][1- g (r)]

0.5

200

400

600

800

1000

200

400

600

800

1000

5

0

-5

temperature [K]
気液共存曲線，飽和圧⼒が⼀致すると，第2ビリアル係数も⼀致する．．．今のところ．
ところが，第3ビリアル係数はまだ上⼿くいっていない（測定値からの算出も難しい︖)

16

MD（4）︓混合系への応⽤

VLE & PvT
【 PvT properties of R1132a/CO2】【VLE pressure of R1132a/CO2】
PR EOS with an assumption of kij 0.02

3.0

2.7

測定値

2.5
2.0
1.5

liq. vap.

MDmeasured (Tomassetti and Di Nicola, 2023)
VTPR-EOS
PR-EOS with an assumption of kij=0.02

263.15 K

1.7

1.2

0.7

1.0
238.15 K

0.5

MD

0.2
220

270

223.15 K

0(Tomassettiら)

2.2

pressure, MPa

equilibrium pressure [MPa]R1132a/CO 2

0.2

0.4

0.6

0.8

-1

320

370

420

temperature, K

1.0mole fraction of R1132a [mol mol ]【VLE pressure of R32/CF3I】
2.0

vapor liquid
Li et al. (2021)
Sakoda-Higashi (2021)
1.6
MD
PR-EOS (Sakoda, 2021)
1.4
XR0 model

equilibrium pressure [MPa]

1.8

293 K

1.2

単⼀物質で微調整した⼒場を⽤いると，混合系でも良く⼀致する．
今のところ．偶然かもしれないが．．．

1.0
0.8

263 K

0.6
0.4
0.2
0

R32/CF3I

0.2
0.4
0.6
0.8
R32 mole fraction [mol mol -1]

1.0

17

MD（5）︓測定の（スーパー）アシスト
表⾯張⼒測定分⼦シミュレーション
enrichment, E imax( ρi ( z ) )max( ρi , ρi )Vapor-liquidinterfacerelativeabsorption
Bulk
Liquid

⽰差⽑管上昇法（Okada, Higashiら）
0.15 ｍNm-1の精度で-45℃まで測定可能（今週には-80℃まで達成︖）

x, yz

⽰差⽑管上昇法による測定には飽和密度の値が必要

臨界点付近の⾼精度測定値を参照し，
低温域の飽和密度を分⼦シミュレーショ
ンで計算．
Parachor法で密度との相関を検証しつつ．

Parachor法による混合系表⾯張⼒の推算には飽和密度と組成の値が必要

低温域のVLEを分⼦シミュレーションで計算．Parachor法で検証しつつ．
気液界⾯を分⼦レベルで観察し，混合系の予測モデルの妥当性を考察．

18

MDの飽和密度と表⾯張⼒測定値の相関

280

Parachor法[Macleod, Transactions of the Faraday Society. 19, 38-41 (1923)]

260

バルクの液および蒸気のモル飽和密度，モル分率から表⾯張⼒を推算する⽅法

  [P]       : 液 の モル飽和密度  mol cm -3 Mixture  : 蒸 気のモ ル 飽和密度  mol cm -3 

x : 液のモル分 率  mol mol -1 y : 蒸気 の モル分率  mol mol -1 

4N   Pi    xi    yi   i 1[PR1336mzz(Z)]=232[PR1336mzz(E)]=221

220[PR1224yd(Z)]=207[ P ]: Parachor
4240

Parachor, [P]Single[PRE356mmz]=259

200[PR1233zd(E)]=202[PR245fa]=193[PR1234ze(Z)]=178180[PR1243zf]=168160[PR1234yf]=172 [PR1234ze(E)]=172

140[PEthanol]=127120[PR1123]=123[PR1132(E)]=114

100※後に，Fowler(1937)によって統計的に検証された

80[PR32]=92

250

300
350
400
Temperature, T [ K ]450 パラコール[P]は，多くの有機物質において換算温度約0.9以下で温度に依存せず， ⼀定の値をとる事が経験的に知ら
れている

19

Parachor法

は あらゆる

混合系に利⽤可能か︖

Parachor法 [Macleod, Transactions of the Faraday Society. 19, 38-41 (1923)]
バルクの液および蒸気のモル飽和密度，モル分率から表⾯張⼒を推算する⽅法
Single[ P ]: Parachor  [P]       : 液 の モル飽和密度  mol cm -3 Mixture  : 蒸 気のモ ル 飽和密度  mol cm -3 

4

4

N

   Pi    xi    yi   i 1


x : 液のモル分 率  mol mol -1 y : 蒸気 の モル分率  mol mol -1 ※後に，Fowler(1937)によって統計的に検証された

 HFC，HFOなどのフロン類では測定値とパラコールの推算値は⽐較的よく⼀致する
 ⽔と有機化合物の混合物では，測定値とパラコールの推算値は⼤きく乖離し，推奨されない．
[Log et al., Int. J Thermophys. 44, 110 (2023)]20表⾯活性とは濃縮度, E
密度プロファイルの⾮単調性を定量化した指標

Ei max( ρi ( z ) )max( ρi , ρi )Ei  α0  fa (xi) fb (VLE)(経験モデル)[Stephan & Hasse, Int. Rev. Phys. Chem. 39, 319-349 (2020)]

enrichment, E 

相対吸着量, Γi

界⾯における単位⾯積当たりの吸着分⼦数を定量化した指標1  σ 
Γ2(1)  
( RT )   ln a 2 Tmax( ρi ( z ) )max( ρi , ρi )Vapor-liquidinterfacerelativeabsorption(ギブスの吸着式)

Bulk
Liquid ( ρ1( Z )  ρ1 ) ( ρ2( Z )  ρ2 ) 
Γ2(1)  ( ρ2  ρ2 )  dz( ρ1  ρ1)( ρ2  ρ2 ) 
 

x, y
z

21

濃縮度 ︓ 局所成分密度のprofile
濃縮度
Water/Ethanol

Water: SPC/E
Ethanol: CCSD/DGDZVP +RESP + AMBER3

Ei 

Vapor-liquidinterfacemax( ρi ( z ) )max( ρi , ρi )Bulk
Liquidx, yxethanol=0.1 [mol/mol]

xethanol=0.2 [mol/mol]z

R32/CF3I

relativeabsorption

xR32=0.1 [mol/mol]

xR32=0.33 [mol/mol]

xR32=0.67 [mol/mol]22相対吸着量ethanol
water

1

2

3

界⾯
バルク

相対吸着量, Γ
界⾯における単位⾯積当たりの吸着分⼦数を定量化した指標(ギブスの吸着式)

R32/CF3I (at 263K)
MD
fitting STN curve(Gibbs adsorption equation)

6

2

Γ

1  σ ( RT )   ln a 2 T

-6Adsorption [mol/m ](1)
2[10 ] 8 ( ρ1( Z )  ρ1 ) ( ρ2( Z )  ρ2 ) Γ  ( ρ2  ρ2 )  dz( ρ1  ρ1)( ρ2  ρ2 ) 
 
(1)
2

 MDとギブスの吸着式の結果は，異なる．
しかし，xR32 = 0.4 [mol/mol]付近で吸着量が最⼤となる似た傾向を⽰した．

ethanol/water (at 298K)
MD
Yano. (2005)

4

2

1 2

3

0
0

0.2

 R32/CF3Iの相対吸着量は，ethanol/waterの相対吸着量に⽐べて，1/8程しかないことが分かる．

0.4

0.6

0.8

1.0

Mole fraction of R32, ethanol [mol/mol]23

相対吸着量とパラコールの予測精度

1  σ ( RT )   ln a 2 T

(ギブスの吸着式)

 ( ρ1( Z )  ρ1 ) ( ρ2( Z )  ρ2 ) 
Γ2(1)  ( ρ2  ρ2 )  dz( ρ1  ρ1)( ρ2  ρ2 )  
 R32/CF3Iの表⾯張⼒は，パラコールの予測値とよく⼀致している ethanol/waterの表⾯張⼒は，パラコール法の予測と⼤きく乖離している パラコール法は，界⾯での表⾯過剰量が⼤きい混合物では表⾯張⼒の予測に失敗するが，ほとんどの冷媒の混合物で適切に予測ができている．

ethanol/water (at 298K)
MD
Strey et al. (1999)
parachor method

4020[10-6] 8
R32/CF3I (at 263K)
MD
fitting STN curve(Gibbs adsorption equation)

62

Γ2(1)  

R32/CF3I (at 263K)
MD
Kawahara et al. (2023)
fitting STN curve
parachor method

60

Adsorption [mol/m ]界⾯における単位⾯積当たりの吸着分⼦数を定量化した指標

Surface tension [mN/m]

80

相対吸着量, Γethanol/water (at 298K)
MD
Yano. (2005)

4

2

0
0

0.2

0.4

0.6

0.8

Mole fraction of R32, ethanol [mol/mol]

1.0

24

分⼦配向︓バルク液と界⾯での違い

r : radiusa : width of red area

 定めたベクトルとｚ軸とがなす⾓θを算出し、分⼦配向の確⽴密度分布を算出

 液バルクの分⼦配向はR32，
CF3Iともにsinθ型の関数に従い，完全にランダムであ
る．

0.025
frequency distribution [ - ]S  2 r sin   a

0.020
0.015
0.010
0.005

bulk (sin)

0
0

20

40

60

80 100 120 140 160 180
degree[° ] R32は界⾯とバルクでわずかな乖離が⾒られ，CF3Iは60度付近に明らかなピークが⾒られる。

25

分⼦の極性と分⼦配向の関係
R32

CF3I

ethanol

5.0e-2-5.0e-2図

H2 O

等電荷密度⾯

0.0004 eÅ-3 上に表⽰した静電ポテンシャル

 分⼦は，分⼦間⼒の強い部分（⽔素結合，⼤きな静電引⼒が作⽤する部分）を液体内部に残す配向状態をとる R32と⽔は，丸みを帯びた⼩さな分⼦軌道となっており，局所的な電荷の偏りを持たない．
→分⼦配向による界⾯⾃由エネルギーへの影響は⼩さく、配向傾向はわずか
 CF3Iとethanolは，局所的に電荷の偏った部分を有しておりその部分を液体内部に向けた配向状態をとった．
26

結⾔
地球温暖化に影響が⼤きいとされる冷媒の開発のため，分⼦シミュレーションを導⼊した．
計算精度が不⼗分なものもあるが，以下のことが可能となった．

量⼦化学計算放射強制⼒の計算（地球温暖化係数の予測）
理想気体⽐熱の計算，
分⼦動⼒学計算に必要な分⼦⼒場の構築

分⼦動⼒学計算
気液飽和曲線（超低温域への外挿）
飽和圧⼒（超低温域への外挿）

ビリアル係数 (超低温，超⾼温域への外挿)
混合系のVLE（低温域への外挿，より多くの組成・組合せ）
表⾯張⼒（測定に必要な

密度・組成の計算．

界⾯の濃縮度・相対吸吸着量・分⼦配向の計算）

27

最後に分⼦シミュレーションは測定にとって代わるものではありません．精度が悪いというご批判はごもっともです．
とはいえ，物性物理，先端計算化学の分野も⽇々進歩しており，より正確で，計算コストを抑えた⼿法が考案されています．それら分野の努⼒・成果をエンジニアリング分野が享受して，温暖化緩和へ向け迅速に対応することが⼤切ではないでしょうか．
未だ分⼦シミュレーションには，満⾜な精度を得るため⾼精度な測定データが必要です．
測定値を参照して検証が出来た場合，測定が難しい低温・⾼温域や，沢⼭の混合組合せについて計算が可能となります．
そして⼀度に複数の物性を取得し，物性間の相関関係を確認することが出来ます．
冷凍空調分野でこの利点を⽣かすべきではないでしょうか．
分⼦シミュレーションによるミクロスコピックな観察は，様々な物理的背景の理解を助け，さらに精度の良い物性計算へ発展させるために有⽤ではないでしょうか．
謝辞
貴重な測定データや状態⽅程式をご提供くださった，九州⼤学
粥川先⽣，九産⼤ ⾚坂先⽣ に⼼より御礼申し上げます．

東先⽣ 迫⽥先⽣，産総研

狩野先⽣

分⼦動⼒学計算はJST 創発的研究⽀援事業（吉⽥パネル）JPMJFR2172の，量⼦化学計算は科研費基盤（C）22K03946の⽀援を受け実施しました．ここに記し謝意を表します．
28