高分子混合系の共存線を計算する

説明足りてないところはいつか書く予定です。

まえがき

2成分の高分子混合系が相分離するかどうかは、自由エネルギーできまる。
この現象を端的に表現できるモデルがFlory-Huggins自由エネルギーであり、以下のようにあらわされる。

F = \chi\phi(1-\phi) + \frac{\phi}{N_1}\log\phi + \frac{1-\phi}{N_2}\log(1-\phi)

ここで、$\chi$ は2種類の高分子成分間の相互作用パラメタ、 $\phi$ は第一成分の体積分率、 $N_1, N_2$ はそれぞれ第一、第二成分のセグメント数(≒重合度)である。

簡単にこの式を見てみると、第1項は相互作用項。第1成分と第2成分がコンタクトする確率はそれぞれの体積分率の積に比例し、その比例係数が $\chi$ となっている。$\chi$ が大きいほど相互作用が強く、相分離しやすい。
第2，3項はエントロピー項。高分子を構成するセグメントをランダムに配置するときのエントロピーを計算しているが、セグメントはつながっているという制約から配置に制限がかかり、エントロピーは1/セグメント数に見積もられる。(両辺に無次元量の $N_1$ をかければすぐわかる通り)相互作用を大きくすることとセグメント数を増やすことは等価であり、高分子量体を含む系は相分離しやすいことがわかる。

本記事では高分子混合系の相図、特に2相の共存線の計算を説明する。
共存線は自由エネルギーに共通接線を引くことによって得ることができるが、対称的な場合、つまり $N_1 = N_2$ の場合を除いて解析的に得ることはできない。

共存線の考え方

そのうち書く

対称的な場合の共存線

まずは対称的な場合を確認してみる。

def free_energy(phi, coeffs):
    chi = coeffs[0]
    N1 = coeffs[1]
    N2 = coeffs[2]
    return chi*phi*(1-phi) + phi*np.log(phi)/N1 + (1-phi)*np.log(1-phi)/N2

coeffs = [5.0, 3.0, 2.0] # chi, N1, N2

phi = np.linspace(0, 1, 100)
plt.plot(phi, free_energy(phi, coeffs))

(式を見れば自明だが)自由エネルギーを描いてみればわかる通り、$\phi=0.5$ に対して対称である。
自由エネルギーの形から明らかに、共通接線は最小値の2点で引くことができる。最小値の $\phi$ は自由エネルギーを $\phi$ について微分して0になるところなので、

\frac{\partial F}{\partial \phi} = \chi (1-\phi) + \frac{1}{N_1}\log\left( \frac{\phi}{1-\phi}\right)
= 0

として得られて、結局

\chi = -\frac{1}{N_1(1-2\phi)}\log\left(\frac{\phi}{1-\phi}\right)

が得られるので

coeffs = [1.3, 2.0, 2.0]
phi = np.linspace(0, 1, 100)
plt.xlabel('$\phi$')
plt.ylabel('$\chi$')
chi = coeffs[0]
N = coeffs[1]
plt.plot(phi, -np.log(phi/(1-phi))/(N*(1-2*phi)), color=dora[0])

とすると共存線が得られる。

非対称な場合の共存線

ここからが本題。
結論から言ってしまえば、2本の接線が等しいという条件でNewton法を使う。
共存組成となる体積分率をそれぞれ $\phi_1, \phi_2$ とすると、2本の接線

\begin{cases}
y_1 = F^{\prime}(\phi_1) (\phi-\phi_1) + F(\phi_1), \\\
y_2 = F^{\prime}(\phi_2) (\phi-\phi_2) + F(\phi_2),
\end{cases}

が等しいと言い換えることができる。これを言い換えると、

\begin{cases}
F^{\prime}(\phi_1) - F^{\prime}(\phi_2) = 0, \\\
-\phi_1F^{\prime}(\phi_1) + \phi_2 F^{\prime}(\phi_2) + F(\phi_1) - F(\phi_2) = 0.
\end{cases}

1つ目の式は傾きが等しい条件であり、2つ目の式は任意の $\phi$ について2本の接線の値が等しい条件である。この2つの式を満たす $\{ \phi_1, \phi_2 \}$ が求める共存組成だ。

数値計算のための下準備

ここでは2次元ニュートン法 (2次元ニュートン法による円と円の交点の数値計算 )を使ってみる。
変数はリンク先の記述に準じ、ここでは結果だけを示す。

\begin{cases}
J_{11} = -2\chi + \frac{1}{N_1\phi_1} + \frac{1}{N_2(1-\phi_1)}, \\\
J_{12} = 2\chi + \frac{1}{N_1\phi_2} + \frac{1}{N_2(1-\phi_2)}, \\\
J_{21} = -\frac{N2(1-\phi_1)^2+N_1\phi_1^2}{N_1N_2\phi_1(1-\phi_1)^2}, \\\
J_{22} =  \frac{N2(1-\phi_2)^2+N_1\phi_2^2}{N_1N_2\phi_2(1-\phi_2)^2}. 
\end{cases}

計算結果をそのままプログラムに起こしてもいいが、以下の関数を定義すると見通しが良くなる。

j_1(\phi) = -2\chi + \frac{1}{N_1\phi} + \frac{1}{N_2(1-\phi)}, \\\
j_2(\phi) = -\frac{N_2(1-\phi)^2+N_1\phi^2}{N_1N_2\phi(1-\phi)^2}.

あと、簡単に計算できるが、自由エネルギーの微分とスピノーダルとなる体積分率(自由エネルギーの2階微分を $\phi$ について解く)も示しておく。

F^{\prime}(\phi) = \chi(1-2\phi) +\frac{1}{N_1}\log\phi - \frac{1}{N_2}\log(1-\phi),　\\\
\phi_c = \frac{\alpha - N_1 + N_2 \pm \sqrt{(\alpha+N_2-N_1)^2-4\alpha N_2}}{2\alpha},\\\
\alpha = 2\chi N_1N_2.

補足

あとから見直したら自分でもわからなかったので、もう少し丁寧に。

基本的な考え方としては、接線の「傾き」と

python コード

def free_energy(phi, coeffs):
    chi = coeffs[0]
    N1 = coeffs[1]
    N2 = coeffs[2]
    return chi*phi*(1-phi) + phi*np.log(phi)/N1 + (1-phi)*np.log(1-phi)/N2


def spinodal(coeffs):
    chi = coeffs[0]
    N1 = coeffs[1]
    N2 = coeffs[2]
    alpha = 2*chi*N1*N2
    common = (alpha+N2-N1) / (2*alpha)
    r = np.sqrt( (alpha+N2-N1)**2-4*alpha*N2 )/(2*alpha)
    return common+r, common-r


def dfdphi(phi, coeffs):
    chi = coeffs[0]
    N1 = coeffs[1]
    N2 = coeffs[2]
    return chi*(1-2*phi) + (1+np.log(phi))/N1 - (1+np.log(1-phi))/N2


def fe_tangent(phi, p0, coeffs):
    chiN = dfdphi(p0, coeffs)*(phi-p0)+free_energy(p0, coeffs)
    return chiN


def j1(phi, coeffs):
    chi = coeffs[0]
    N1 = coeffs[1]
    N2 = coeffs[2]
    return -2*chi + 1/(N1*phi) + 1/(N2*phi)


def j2(phi, coeffs):
    chi = coeffs[0]
    N1 = coeffs[1]
    N2 = coeffs[2]
    return - (N2*(1-phi)**2 + N1*phi**2)/(N1*N2*phi*(1-phi)**2)


def fx(psi, coeffs):
    return dfdphi(psi[0], coeffs) - dfdphi(psi[1], coeffs)


def fy(psi, coeffs):
    return free_energy(psi[0], coeffs) - free_energy(psi[1], coeffs) - psi[0]*dfdphi(psi[0], coeffs) + psi[1]*dfdphi(psi[1], coeffs)


def fxfy(psi, coeffs):
    x = fx(psi, coeffs)
    y = fy(psi, coeffs)
    return np.array([x, y])


def jacobian(psi, coeffs):
    p1, p2 = psi
    J11 = j1(p1, coeffs)
    J12 = -j1(p2, coeffs)
    J21 = j2(p1, coeffs)
    J22 = -j2(p2, coeffs)
    return np.array([[J11, J12], [J21, J22]])

やりたいことの本体。

def binodal(coeffs, psi=np.array([-1, -1]), verbose=False, eps=1e-10, max_loop=1e9):
    if psi[0] < 0.0:
        p1, p2 = spinodal(coeffs)
        psi = np.array([p1, p2])
    delta = np.array([1e1, 1e1])
    k = 0
    flag = True
    while (abs(delta[0]) > eps and abs(delta[1]) > eps):
        try:
            J = jacobian(psi, coeffs)
            invJ = np.linalg.inv(J)
            f = fxfy(psi, coeffs)
            delta = np.dot(invJ, f)
            while(psi[0]-delta[0] > 1.0 or psi[1]-delta[1] < 0.0):
                # 0<phi<1に収まらないとlog(phi), log(1-phi)がnanになるのでdeltaを適度に小さくする
                delta *= 0.5
            psi = psi - delta
        except:
            print('could not calculate inverse J')
            print('J = ')
            print(J)
            print('delta = ')
            print(delta)
            print('psi = ')
            print(psi)
            flag = False
            break

        if (verbose and k % 10000 == 0):
            print('k = %s:' % k)
            print('dphi1 = {:>12.5e}'.format(
                delta[0]), ',\t dphi2 = {:>12.5e}'.format(delta[1]))
            print('phi1  = {:>12.5e}'.format(
                psi[0]),   ',\t phi2  = {:>12.5e}'.format(psi[1]))
        k += 1
        if (k > max_loop):
            print('could not convergence')
            flag = False
            break
    return psi, flag

共通接線が得られているかどうか、確認

coeffs = [5.0, 3.0, 2.0]
psi, flag = binodal(coeffs, eps=1e-15)
phi = np.linspace(0, 1, 100)
plt.plot(phi, free_energy(phi, coeffs))
plt.plot(phi, fe_tangent(phi, psi[0], coeffs))
plt.plot(phi, fe_tangent(phi, psi[1], coeffs))

こんな感じで共通接線を見つけることができる。相図を書きたければ $\chi$ の値だけ適当な刻みで大きくしていけば書くことができる。

相図作成

def phase_diagram(max_chi, N1, N2, dev=50, eps=1e-15):
    # 計算を効率良くするために臨界点を計算しておく
    phi_c = np.sqrt(N2)/(np.sqrt(N1)+np.sqrt(N2))
    min_chi = (np.sqrt(N1)+np.sqrt(N2))**2/(2*N1*N2)
    
    # chiの刻み幅
    delta_chi = (max_chi-min_chi) / dev
    chi = min_chi

    binodal_line = []
    spinodal_line = []

    binodal_line.append([phi_c, chi])
    spinodal_line.append([phi_c, chi])
    
    # 最初の共存線を得るための初期値
    psi = spinodal([chi+delta_chi, N1, N2])

    while(chi < max_chi):
        chi = chi + delta_chi
        coeffs = [chi, N1, N2]
        psi, flag = binodal(coeffs, psi, eps=eps)
        if(flag): # binodal()計算が収束したときのみデータを追加
            binodal_line.append([psi[0], chi])
            binodal_line.append([psi[1], chi])
            phi_sp = spinodal([chi, N1, N2])
            spinodal_line.append([phi_sp[0], chi])
            spinodal_line.append([phi_sp[1], chi])
        else: # 収束しない場合さらに大きなchiでも収束しないので打ち切る
            break
    return np.array(sorted(spinodal_line)), np.array(sorted(binodal_line)) # phi についてソートしておくと共存線書きやすい

# 実施例
max_chi = 8
N1 = 3
N2 = 1

sp_line, bi_line = phase_diagram(max_chi, N1, N2)
plt.xlabel('$\phi$')
plt.ylabel('$\chi$')
plt.xlim([0, 1])
plt.plot(sp_line[:, 0], sp_line[:, 1], color=dora[0], label='spinodal')
plt.plot(bi_line[:, 0], bi_line[:, 1], color=dora[1], label='binodal')
plt.legend()