04. 결정핵생성기구 구현편

앞에서 시스템의 임계반경 이상의 결정핵이 생기기 위한 수식에 대하여 알아봤고, 이번편은 그러한 수식을 파이썬으로 구현하는것이다. 우선 제일 먼저 구현할 코드는 임계반경을 구하는 코드이다.

임계반경 구하는 식

수식에 대한 자세한 설명은 앞선편에서 했기때문에 앞으로도 생략하겠다.

 

우선, 물질의 고유한 특성과 연관된 표면자유에너지, 부피자유에너지, 용융점 이 3가지는 철을 기준으로 잡고서 해당 값들을 입력해주고, 계산해보면 다음과 같다.

def critical_radious(mt, t, me, gamma): # mt = 용융점(k) , t = 현재온도(k) , me = 용융잠열(J/m^3) , gamma = 표면에너지( J/m^2 )

    r = (-2 * gamma * mt) / ((mt - t) *  me)

    return r

melting_tempearture = 1811 #K

temperature = 1537 #K

melting_energy = -1.85 * 10**9 #J/m^3

surface_enrgy = 0.204 #J/m^2

critical_r = critical_radious(melting_tempearture , temperature, melting_energy , surface_enrgy)

print(f"{critical_r}m^2")

 

출력값 : 1.4576602880252514e-09m^2 (약 1.5나노미터제곱)

 

임계반경에 대한 계산이 끝났으면, 다음은 임계에너지에 대한 계산을 해주면 된다. 식은 아래와 같다.

 

def critical_energy(mt, t, me, gamma): # mt = 용융점(k) , t = 현재온도(k) , me = 용융잠열(J/m^3) , gamma = 표면에너지( J/m^2 )

    E = ((16*np.pi*(gamma**3))/(3*(me**2))) * (mt**2 / (mt-t)**2)

    return E

melting_tempearture = 1846 #K

temperature = 1546 #K

melting_energy = -1.85 * 10**9 #J/m^3

surface_enrgy = 0.204 #J/m^2

critical_E = critical_energy(melting_tempearture , temperature, melting_energy , surface_enrgy)

print(f"{critical_E}J")

출력값 : 1.573680379254878e-18J

그리고 최종적으로 이를 볼츠만 분포의 활성화에너지 항에 넣으면 완성이다. 수식은 아래와 같다.

 

완성한 전체 프로그램은 아래와 같다.

 

 

import numpy as np

def critical_radious(mt, t, me, gamma): # mt = 용융점(k) , t = 현재온도(k) , me = 용융잠열(J/g) , gamma = 표면에너지(J)

    r = (-2 * gamma * mt) / ((mt-t) *  me)

    return r

def critical_energy(mt, t, me, gamma): # mt = 용융점(k) , t = 현재온도(k) , me = 용융잠열(J/m^3) , gamma = 표면에너지( J/m^2 )

    E = ((16*np.pi*(gamma**3))/(3*(me**2))) * (mt**2 / (mt-t)**2)

    return E

def boltzmann_distribution(critical_energy, t):

    distribution = np.exp(-(critical_energy)/(boltzmann_constance*t))

    return distribution

boltzmann_constance = 1.380 * 10**-23 #J/K

melting_tempearture = 1846 #K

temperature = 1446 #K

melting_energy = -1.85 * 10**9 #J/m^3

surface_enrgy = 0.204 #J/m^2

N = 6.02 * 10**22 # /개

critical_r = critical_radious(melting_tempearture , temperature, melting_energy , surface_enrgy) #임계반경 계산

critical_E = critical_energy(melting_tempearture , temperature, melting_energy , surface_enrgy) #임계에너지 계산

num_of_nuclei = N * boltzmann_distribution(critical_E, temperature)

#계산값 출력구문 (확인용)

print(f"{critical_r}m^2")

print(f"{critical_E}J")

print(f"{num_of_nuclei}")

N은 계산 편의상 아보가르드수로 넣었으며, 동작시키면 계산값들이 정상적으로 나온다. 하지만 이러한 계산값들이 올바르게 동작하는지 알 필요가 있기때문에 온도차이Vs결정핵의 생성갯수로 그래프를 플롯팅하여, 원자핵의 생성의 추세를 확인하는 작업이 필요하다.

 

#빈배열 생성

x = []

y = []

#용융점에서 10도씩 낮추는 배열생성

for i in range(10, 1010 ,10):

    x.append((melting_tempearture - i))

x.reverse()

#위의 값을 이용하여 100개의 온도

for i in range(0,100):

    m = critical_energy(melting_tempearture , x[i], melting_energy , surface_enrgy)

    n = 100000 * boltzmann_distribution(m, temperature)

    y.append(n)

plt.plot(x,y)

plt.grid()

plt.show()

1800~800도 사이의 원자핵 생성분포

전체 계의 원자수를 만개라고 가정하면, 800도에서는 약 275개가량의 원자핵을 가지게됨을 알수있다. 실제보다 훨씬 많은 핵생성으로 보여지기때문에, 상수들의 수정이 필요해보인다.

 

프로그램은 성공적으로 작동한다. 하지만 이 프로그램은 아래와 같은 한계점을 보인다.

 

1. 균일 핵생성만을 가정했기때문에 비균일 핵생성은 모델링하지 못한다.

2. 핵생성에 필요한 에너지를 부피와 계면에 따른 자유에너지변화에 의한 임계반경에너지의 볼츠만 분포로만 모델링했지만, 실제로는 원자의 확산도 필요하기때문에 이에 대한 항이 추가되어야한다.

 

이 두개의 문제점은 좀 더 시간을 들인다면 해결가능하다. 하지만 일단 이번회차는 여기까지하고, 다음주부터는 위의 두가지 문제를 고려하여 해결해볼지, 아니면 실제 핵생성과 성장을 모델링할지 고민해보도록하겠다.

 

 

------------------------

2024.04.04

시험기간이라 바쁘긴하지만 아무리 생각해도 결정핵생성을 모델링하는데 임계반경에 의한 임계에너지를 볼츠만분포에 의하여 계산하는건 부적합한것 같아서 급하게 수정한다. 가장 큰 문제는 급냉이 커질수록(온도가 낮아질수록), 핵의 갯수가 기하급수적으로 무한정 증가한다는것인데, 이는 사실상 불가능하다. 실제로는 온도가 낮아질수록 확산에 필요한 에너지가 부족해지기 때문에, 실제로는 특정온도에서 최고점을 가져야하기때문에, 두개의 EXP함수로 나타내는것이 바람직하다. 실제로는 아래의 식이 좀 더 실제와 가깝다.

 

G*는 임계반경에 의한 에너지,  -Qd는 확산에 의한 항이다.

 

철의 원자핵 생성에 필요한 Self Diffusion energys는 철의 자기확산계수로 가정하여 진행한다. 도저히 찾을수가없어, 임의의 값을 넣었다.

 

완성된 전체프로그램과 실행결과는 밑에.

import numpy as np

import matplotlib.pyplot as plt

def critical_radious(mt, t, me, gamma): # mt = 용융점(k) , t = 현재온도(k) , me = 용융잠열(J/g) , gamma = 표면에너지(J/m^2)

    r = (-2 * gamma * mt) / ((mt-t) *  me)

    return r

def critical_energy(mt, t, me, gamma): # mt = 용융점(k) , t = 현재온도(k) , me = 용융잠열(J/m^3) , gamma = 표면에너지( J/m^2 )

    E = ((16*np.pi*(gamma**3))/(3*(me**2))) * ((mt**2) / (mt-t)**2)

    return E

def boltzmann_distribution(critical_energy, self_diffusion_energy, t):

    distribution = np.exp(-(critical_energy)/(boltzmann_constance*t)) * np.exp(-(self_diffusion_energy)/(boltzmann_constance*t))

    return distribution

boltzmann_constance = 1.380 * 10**-23 #J/K

melting_tempearture = 1800#K

temperature = 1700 #K

melting_energy = -1.85 * 10**9 #J/m^3

surface_energy = 0.204 #J/m^2

N = 6.02 * 10**22 # /개

SDE = 1.021 * (10**-18)#자기확산에너지

critical_r = critical_radious(melting_tempearture , temperature, melting_energy , surface_energy) #임계반경 계산

critical_E = critical_energy(melting_tempearture , temperature, melting_energy , surface_energy) #임계에너지 계산

num_of_nuclei =  N * boltzmann_distribution(critical_E, SDE, temperature)

print(np.exp(-(SDE)/(boltzmann_constance* temperature)))

print(np.exp(-(critical_E)/(boltzmann_constance*temperature)))

#계산값 출력구문 (확인용)

print(f"{critical_r}m^2")

print(f"{critical_E}J")

print(f"{num_of_nuclei}개")

print(f"{SDE}J is SDE")

# 온도 배열 생성

temperatures = np.linspace(100, melting_tempearture-1, 500)  # 100K부터 용융점까지

# 각 온도에 대한 임계에너지와 볼츠만 분포 계산

critical_energies = critical_energy(melting_tempearture, temperatures, melting_energy, surface_energy)

boltzmann_distributions = boltzmann_distribution(critical_energies, SDE, temperatures)

# 그래프 그리기

plt.figure(figsize=(10, 6))

plt.plot(temperatures, boltzmann_distributions, label='Boltzmann Distribution')

plt.xlabel('Temperature (K)')

plt.ylabel('Distribution')

plt.grid(True)

plt.legend()

plt.show()

결과물이 만족스럽지는 않다. 그래프가 너무 좁게 나타났다. 아마 계수를 적절하게 만진다면 해결할수있을것같으니 물리적 상수들에 대해서 좀더 고민해봐야겠다.

 

ps. 임계반경은 구해놓고 사용하지않았는데, 이는 이후에 모델링을 실제로 작성할때 필요해서 구해놓은 값이다.