追記8/15
この計算をNVIDIA製GPU積んだゲーミングPCで高速化した記事を投稿しました．
追記9/11
ルンゲクッタ法の更新式で両方の変数を動かす必要があるところ、片方しか動かしておらず間違えていました。全体の整合を取りながらうまく修正する術を思いつかなかったので、ルンゲクッタ法以外の時間発展アルゴリズムの追加やポテンシャルの追加を行いながら修正した記事を投稿しました。
環境
- Mac mini(M1) Big Sur
- Python 3.9.6(arm64)
- matplotlib 3.4.3
- numpy 1.21.1
- scipy 1.7.1
目標
2次元に配置した原子の格子振動を計算して熱伝導の過程を可視化したい．
背景
唐突ですが皆さんのスマホは今日もカイロになってますか？スマホに限らず身の回りの電子製品は高機能化に伴って発熱量が増す一方です．半導体は昇温すると電気抵抗が低下するために，一度発熱すると発熱スパイラルに陥ります．そうすると微細な回路が溶けたり基板が焼けたりして壊れますし，リチウムイオン電池が繋がれている場合は最悪，電解液が沸騰して吹き出して引火して燃えます．
こんな事態を未然に防ぐために発熱量の大きな電子機器には放熱機構が取り付けられます．分かりやすいのはPCのCPUにあるヒートシンクですね．空冷用のファンなんかもそうです．
ところでCPUとヒートシンクの間に何か塗りませんか？そう熱伝導グリースとかいうやつです．厚塗りするとこいつ自体の熱抵抗で逆に放熱性が落ちますが，全く塗らないと全然ヒートシンクに熱が逃がせない大事なやつですね．今回は発熱するCPUから熱伝導グリースへと熱が伝わる過程を分子レベルでシミュレーション・・・・・・するために，まずは単純な2次元格子点の熱運動を可視化してみます．
本当はLAMMPSでも使って分子動力学の計算をするのが正しいのですが，素人が独力で触るにはポテンシャル選択の敷居が高すぎるので，バネの運動でお茶を濁します．
追記8/16
盆休みも終わりなのでしばらく追究しませんが、やりたかったのは基準振動数が大きく異なる（はず）の金属ヒートシンクと有機グリースの間の熱抵抗は0にはできませんよということです。0ではないはずだけどどれくらいだ？無視できる程度か？というのを概算したかったです。国立大学工学部が共同運営している？こちらの記事の動画2dなんかいい感じですね。
問題設定
原子配置
以下の図の様に5*10個の質量$m$の質点を配置します．面倒なので明示していませんが，隣り合う格子点はバネ定数$k$，長さ$l$のバネで接続してあります．境界条件として$x$軸方向の両端は強制的に単振動させます．温度勾配は両端の振動数$\omega$と振幅$a$に差を与えることで生じさせます．ここで質点の質量$m, k, l, \omega, a$はいずれも原子のスケールに合わせます．
運動方程式
各原子をバネで繋いているので以下の式が成り立ちます．
m \frac{d^2X_i}{dt^2} = \sum_{j=0}^{n-1} -k_{ij}\left(\left|X_i - X_j\right|-l_{ij}\right)\frac{X_i-X_j}{|X_i-X_j|}
ここで大文字はベクトル，$k_{ij}, l_{ij}$はそれぞれ原子$k, l$間のバネ定数とバネの長さを表します．原子に付与する番号は左上から順につけることにします．すると2次元配置（$a行b列$）あっても1次元（$i=na+b$)の形で扱うことができます（numpy.ndarray.flattenする）．というか3次元でも同じです．バネの連結状態を与える時などにもこの2次元配置で考えてから1次元化をおこないます．
解法
ルンゲクッタ法を用いて微分方程式の数値解を求めます．要するに$y_{n+1} = y_n + hdy/dt$とするよりも微小変化分の精度あげましょうということです．ルンゲクッタ法は1階常微分方程式しか解けないので，1階微分を別の変数と見なすことで2階常微分方程式を1階常微分方程式の2元連立方程式に書き換えます．
\left\{
\begin{align}
m\frac{dV_i}{dt} &= f(t, V, X) = \sum_{j=0}^{n-1} -k_{ij}\left(\left|X_i - X_j\right|-l_{ij}\right)\frac{X_i-X_j}{|X_i-X_j|} \\
V_i &= g(t, V, X) = \frac{dX_i}{dt}
\end{align}
\right.
こうすることで$V, X$それぞれで差分が計算できるので，後はルンゲクッタ法に従った差分を足し合わせていけば運動が更新できます．下の拙コードでは4次のルンゲクッタ法を自前で書いてしましましたが，scipyの関数を使えばよかった様です（使い方のリンクはsolve_ivpの前身のodeintですがドキュメント読む限り使い方は同じです）．動作確認とアニメーションが欲しかったので2次元格子点しか試していませんが，3次元格子点に配置しても動作するはずです．その場合は初期化部分に新規軸の項を書き足して下さい．格子点じゃない場合は頑張って下さい．
from pathlib import Path
from io import StringIO
import numpy as np
from scipy.sparse import csr_matrix, coo_matrix
from matplotlib import pyplot as plt
from matplotlib import animation
class RungeKutta(object):
def __init__(self, dt, grad):
"""
y_n+1 = y_n + h/6*(k1 + 2*k2 * 2*k3 + k4)
t_n+1 = t_n + h
k1 = f(t, y) = dy/dt|(t,y)
k2 = f[t + h/2, y + h/2*k1)
Paramaters
---------
dt: interval time
grad: gradient
"""
self.dt = dt
self.grad = grad
def calc_k(self, t, y, **kwargs):
k1 = self.grad(t, y, **kwargs)
k2 = self.grad(t + self.dt/2, y + self.dt/2*k1, **kwargs)
k3 = self.grad(t + self.dt/2, y + self.dt/2*k2, **kwargs)
k4 = self.grad(t + self.dt, y + self.dt*k3, **kwargs)
return k1, k2, k3, k4
def calc(self, t, y, **kwargs):
k1, k2, k3, k4 = self.calc_k(t, y, **kwargs)
y_new = y + self.dt/6*(k1 + 2*k2 + 2*k3 + k4)
t_new = t + self.dt
return t_new, y_new
def hooke(t, v, **kwargs):
"""
|F| = k|(|r1-r2|)-l|
= k|((x1-x2)**2 + (y1-y2)**2)**(1/2) - l|
Parameters
----------
X: (n_dims, nx, ny, nz)
t: time
K: constant tensor of hooke's law (n_dims, nx, ny, nz)
L: length tensor of springs
Return
------
force vector: X.shape
"""
X = kwargs['X']
n = X.shape[0]
R = []#scipy.sparse CANNOT handle tensor.
r = 0
for i in range(n):
#force to 0 for NOT connected items
xi = kwargs['C'].multiply(X[i].reshape(-1,1))
xj = kwargs['C'].multiply(X[i].reshape(1,-1))
"""
xi = [[x0, x0, ..., x0],
[x1, x1, ..., x1],
...
[xn, xn, ..., xn]]
xj = [[x0, x1, ..., xn],
[x0, x1, ..., xn],
...
[x0, x1, ..., xn]]
"""
R.append(xi-xj)
r += (xi-xj).power(2)
L = kwargs['C'].multiply(kwargs['L'])
r = r.sqrt()
f = kwargs['K'].multiply(r-L)
F = []
for i in range(n):
F.append(
np.nan_to_num(f.multiply(R[i])/r)
.sum(axis=-1)
.reshape(X[i].shape)
/kwargs['M']
)
F = -np.array(F)
return F
def velosity(t, x, **kwargs):
return kwargs['V']
def boundary(t, phi, a, k, m):
"""
d2x/dt2 = -k/m*x
x = A*exp(sqrt(k/m)t*i)
"""
omega = np.sqrt(k/m)
x = np.cos(omega*t + phi)*a*np.array([[0.5, 1]])
#y = np.sin(omega*t + phi)
vx = -omega*np.sin(omega*t + phi)*a*np.array([[0.5, 1]])
#vy = -omega*np.cos(omega*t + phi)
return x, vx
def main(Nx=5, Ny=3, dt=1.e-16, upto=1.e-12,
k=1.e3, l=1.e-10, m=1.e-25, a=3.e-12):
t = 0
# Define K
I = np.ones(Nx*Ny).reshape(Ny,Nx)
E = (np.eye(Nx*Ny))
# Connectivity to upper atom.
# np.roll(Iu*E, -Nx, axis=1) affords
# atom i(row i) is connected(1)
# or not(0) to atom j(column j).
Iu = I.copy()
#Atoms in row 0 are not connected to last row.
Iu[0,:]=0
Iu = Iu.reshape(-1,1)
# Connectivity to lower atom.
Il = I.copy()
Il[:,0] = 0
Il = Il.reshape(-1,1)
# Connectivity to right atom.
Ir = I.copy()
Ir[:,-1] = 0
Ir = Ir.reshape(-1,1)
# Connectivity to left atom.
Id = I.copy()
Id[-1,:] = 0
Id = Id.reshape(-1,1)
C = (np.roll(E*Iu, -Nx, axis=1)
+ np.roll(E*Il, -1, axis=1)
+ np.roll(E*Ir, 1, axis=1)
+ np.roll(E*Id, Nx, axis=1)
)
C = csr_matrix(C)
K = C*k
# Define L
L = C*l
# Define M
M = np.ones(Nx*Ny).reshape(Ny,Nx)*m
# Define x0
# shape: (2, Nx, Ny)
x = np.arange(1,Nx+1)*l
y = np.arange(Ny,0,-1)*l
x, y = np.meshgrid(x,y)
# use for boundary
xb0 = x.copy()
yb0 = y.copy()
# 99% data are in u +- 3s
dx = np.random.randn(Ny,Nx)*l/300
dy = np.random.randn(Ny,Nx)*l/300
dx[:,0] = 0
dx[:,-1] = 0
x = x + dx
y = y + dy
phi = np.random.rand(Ny*2).reshape(Ny, 2)*2*np.pi
#phi = np.zeros(Ny*2).reshape(Ny,2)
xb, vb = boundary(0, phi, a, k, m)
x[:,0] += xb[:,0]
x[:, -1] += xb[:,1]
x = np.stack([x, y])
v = np.zeros_like(x)
v[0,:,0] += vb[:,0]
v[0,:,-1] += vb[:,1]
V = RungeKutta(dt, hooke)
X = RungeKutta(dt, velosity)
folder = Path(
f'k{k}_l{l}_m{m}_a{a}_dt{dt}'
)
folder.mkdir(exist_ok=True)
wd = folder/f'Nx{Nx}_Ny{Ny}_time{upto}'
wd.mkdir(exist_ok=True)
res_t = open(wd/f'time.txt', 'w')
res_vx = open(wd/f'velosity_x.txt', 'w')
res_vy = open(wd/f'velosity_y.txt', 'w')
res_x = open(wd/f'x.txt', 'w')
res_y = open(wd/f'y.txt', 'w')
fig, ax = plt.subplots()
ax.scatter(x[0], x[1], c='k')
fig.savefig(wd/'Init.png', bbox_inches='tight')
np.savetxt(res_t, np.array([t]))
np.savetxt(res_vx, v[0])
np.savetxt(res_vy, v[1])
np.savetxt(res_x, x[0])
np.savetxt(res_y, x[1])
for txt in [res_t, res_vx, res_vy, res_x, res_y]:
txt.write('\n')
while t<=upto:
t_new, v_new = V.calc(
t, v, X=x, K=K,
L=L, M=M, C=C
)
xb ,vb = boundary(t_new, phi, a, k, m)
xb[:,0] += xb0[:,0]
xb[:,-1] += xb0[:, -1]
v_new[0, :, 0] = vb[:,0]
v_new[0, :, -1] = vb[:,1]
v_new[1, :, 0] = 0
v_new[1, :, -1] = 0
t_new, x_new = X.calc(t, x, V=v)
x_new[0, :, 0] = xb[:,0]
x_new[0, :, -1] = xb[:,1]
x_new[1, :, 0] = yb0[:, 0]
x_new[1, :, -1] = yb0[:, -1]
t = t_new
v = v_new
x = x_new
np.savetxt(res_t, np.array([t]))
np.savetxt(res_vx, v[0])
np.savetxt(res_vy, v[1])
np.savetxt(res_x, x[0])
np.savetxt(res_y, x[1])
for txt in [res_t, res_vx, res_vy, res_x, res_y]:
txt.write('\n')
res_t.close()
res_vx.close()
res_vy.close()
res_x.close()
res_y.close()
return wd
def plot(wd, Nx, Ny, dt, interval, time):
res = {}
for txt in [f'time.txt',
f'x.txt',
f'y.txt',
f'velosity_x.txt',
f'velosity_y.txt',
]:
with open(wd/txt) as f:
snap = f.read().split('\n\n')
block = []
for item in snap:
array = np.loadtxt(StringIO(item))
block.append(array)
block = np.array(block[:-1])
res[txt[:-len(f'.txt')]] = block
fig, ax = plt.subplots()
fig.set_size_inches(15,8)
ims = []
#bottle neck
for i in range(len(res['time'])):
axes = ax.scatter(
res['x'][i], res['y'][i],
c='k')
title = ax.text(
0.5, 1.01,
f'time = {res["time"][i]}',
ha='center', va='bottom',
transform=ax.transAxes,
fontsize='large'
)
ims.append([axes, title])
fig.set_size_inches(15,8)
ani = animation.ArtistAnimation(
fig, ims, interval=interval,
blit=True, repeat_delay=1000)
ani.save(wd/f'harmonics_2d.gif')
if __name__=='__main__':
"""キッテル「固体物理学入門 上」丸善出版より引用
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Elastic Stiffness/(*E12 dyn/cm^2)
Crystal C11 C12 C44 Temperature/K density/(g/cm^3)
--------------------------------------------------
Cu 1.684 1.214 0.754 300 9.018(0 K)
Ag 1.240 0.937 0.461 300 10.635(0 K)
Au 1.923 1.631 0.420 300 19.488(0 K)
Al 1.068 0.607 0.282 300 2.733(0 K)
Atomic radii/A
covalent bonding(tetrahedron) Octet ion
--------------------------------------------
Cu 1.35 -
Ag 1.52 1.26
Au - 1.37
Al 1.26 0.50
引用ここまで
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1 dyn = 1e-5 N
1e12 dyn/cm^2 = 1e7 N/cm^2
= 1e3 N/m^2
= 1 kPa
stress = C*epsilon
force = area*stress
k*(x1-x2-l) ~ l**2 * C*epsilon
= l**2 * C*(x1-x2-l)/l
= C*l*(x1-x2-l)
k = C*l
k/(1e-7 Pa)
----------
Cu 2
Ag 2
Au 2.5
Al 1
Atomic weight mass/(1e-26 kg)
-------------------
Cu 63.5 10.5
Ag 107 17.7
Au 196 32.5
Al 26.9 4.46
"""
Nx = 10
Ny = 5
dt = 1.e-16
upto = 1.e-14
interval=1
k = 1.e2
l = 1.e-10
m = 3.e-25
a = l*0.03
wd = main(Nx, Ny, dt, upto, k, l, m, a)
plot(wd, Nx, Ny, dt, interval, upto)
原子の質量$m$が$10^{-25}$ kgに対してバネ定数$k$が$10^2$ N/mと非常に大きいです．このため時間刻みを非常に小さくしないと，原子がとんでもない勢いで吹っ飛んでいく計算になりました．上の例では0.1 fs刻みで10 fsまで100ステップ計算してます．経過時間が短いので振動の様子がほとんど分かりませんね．
1000 fs = 1 psまで10,000ステップ計算してもいいのですが，私の環境で40秒で計算が終わるのに対してアニメーションの描画に20分以上もかかります．またアニメーションの更新間隔を1 msにしているつもりですがそうなってくれないため，10,000ステップ描画させると1ループ表示させるのにも物凄いスロー再生です（元々フェムト秒スケールなのでスロー再生は当たり前ですが）．
結論
とりあえずバネで繋がれた格子点の運動はシミュレーションできるようになったが，アニメーションがスロー再生すぎて正常な運動になっているか確認が難しい．そのアニメーションも描画に時間がかかり過ぎて作業効率が悪い．
今後
- PyQtGraphを使って描画を高速化することで作業効率を上げる
- 平衡状態に到達したことを確認する手法を検討する
- 温度勾配の再現に最適な振幅と振動数を検討する
- 質量やバネ定数，結合数，対称性が異なる2相を連結して熱伝導の様子を可視化する
先は長い・・・
追記8/13
マッカーリ・サイモン「物理化学（上）」東京化学同人によると二原子分子のバネ定数は$10^2$ N/m程度なのに対して，キッテルに記載されている弾性率から雑に計算した金属のバネ定数は$10^{-7}$ N/mと全然違ってしまっています．ひとまずコード中の$k$の概算コメントは無視して下さい．