相场模拟的python实现

空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

相场模拟的python实现

相场模拟

wangchangxin

发布于 2023-08-02

推荐镜像 :Basic Image:bohrium-notebook:2023-03-26

推荐机型 :c2_m4_cpu

相场模拟的python实现

一. 背景

1.1 序参量（相场参量）

1.2 自由能泛函

1.3 核心困难

二. 晶界长大模拟案例

总结与思考：

代码

文本

🎯 本教程旨在快速理解相场模拟及其实现。

在 Bohrium Notebook 界面，你可以点击界面上方蓝色按钮 开始连接，选择基础镜像及任何一款节点配置，稍等片刻即可运行。

一. 背景

代码

文本

相场法(Phase field method) 是以Ginzburg-Landau理论为物理基础，通过微分方程来体现具有特定物理机制的扩散和热力学驱动的综合作用，通过计算机编程求解上述方程，从而获取研究体系（固体界面，固液界面等）在时间和空间上的瞬时状态。

代码

文本

看不懂？没关系。简单来说，相场方法是用来对材料的微观结构进行模拟的经典模型方法，其核心是物理扩散公式，以最小化自由能为原则，实现对相变等物理现象的模拟。相场法研究的更多是微观/介观尺度，而非原子/电子尺度。目前其应用较好的方向有：凝固现象的树枝晶生长，晶粒热处理的长大，结构相变，位错运动等。

代码

文本

1.1 序参量（相场参量）

代码

文本

在相场方法中，对于特定的物理问题，相场模拟需先选择一个变量作为序参量（相场参量），它可以是模型中的具体物理量，也可以是一组变量的函数。该定义好的序参量，可以反应体系的物理性质，比如结晶度（degree of crystallinity）、原子排列（atomic ordering）、晶体对称性（crystal symmetry）、状态（state）、磁化强度（magnetization）、电极化量（electric polarization）、浓度（concentration）、组成（composition）、损伤程度（damage level）等。采用序参量的梯度来追踪定位相界面，便可以观察到界面的演化过程。下面以扩散界面模型作为案例叙述说明。

代码

文本

在传统相变模型中，常因为追踪不同相之间界面的变化而导致大量繁琐的计算。（如下图左图所示），在传统的明锐界面下，界面对能量的贡献与界面能密度和界面拓扑结构有关，而界面形态（interface morphology）轨迹，特别是三维情况下，可能会十分复杂。而相场法假定相之间的界面是扩散性（diffusive）的，介于相邻的相之间的扩散界面由用 $ϕ$ 表示的序参量的连续插值描述，如下图右图所示。

代码

文本

未标题-1.png

代码

文本

扩散模型对比图，纵坐标为相成分，横坐标为位置

代码

文本

简单起见，此时序参量定义为被相组成，红色相为0，蓝色相为1。两者之间为连续梯度变化。（softmax啊有木有！）。在此定义下，微观相分布类似于序变量的一个分布函数 $ϕ (x)$ ，其中 $x$ 是位置向量。而在扩散界面模型下，界面能够通过序参量的梯度进行识别。因此界面形态和界面能量都能够很自然地基于梯度项定义。界面移动，代表两相的比例变化。

代码

文本

有了上述定义，如何进行模拟过程呢？

代码

文本

1.2 自由能泛函

代码

文本

在模拟计算前，我们先根据序参量确定优化的能量目标。这基于两个假定：

代码

文本

假定1： 每个相的自由能密度都依赖于相的热力学状态，因此每个相的体自由能密度（bulk free energy density）都能够由含 $ϕ$ 的公式描述。

假定2： 界面上的自由能密度是不同于相内部，界面位置能够用 $ϕ$ 的梯度来描述。

代码

文本

系统的自由能 $F$ 能够由以下相场自由能泛函给定：

$F (ϕ, \nabla ϕ) = \int_{A} [f (ϕ) + 1/2 k ∣\nabla ϕ ∣^{2}] d V （公式 1 ）$

代码

文本

其中， $F (ϕ)$ 是体自由能密度， $1/2 k ∣\nabla ϕ ∣^{2}$ 是梯度能量密度。（更加精细的话，额外需要考虑∇𝜙的二次项 $α (\nabla ϕ)^{2}$ 。）

代码

文本

在相场理论中，系统的状态和动力学由能量最小原理决定。首先，若系统处于平衡状态时，序参量 $ϕ (x)$ 应当使得自由能 $F$ 取得最小值，即平衡状态下，自由能的一阶变分导数为零，这便是我们的最终目标。

$\frac{\partial F ( ϕ , \nabla ϕ )}{\partial ϕ ( x )} = 0 （公式 2 ）$

代码

文本

其次，对于不受约束的热力学量，微结构演化过程是以整体能量最速下降的方向进行的。在相场理论中，这种关于微结构的整体能量变化由一阶变分导数给定。

代码

文本

$δ F (ϕ, \nabla ϕ) = δ \int_{A} [f (ϕ) + 1/2 k ∣\nabla ϕ ∣^{2}] d V （公式 4 ） = \int_{A} [\frac{\partial f ( ϕ )}{\partial ϕ} - k Δ ϕ] δ ϕ d V + \int_{\partial A} k \nabla δ n ϕ d S$

代码

文本

$k \nabla δ n$ =0,上方程式为：

$\frac{δ F}{δ ϕ} = \frac{\partial f ( ϕ}{\partial ϕ )} - k Δ ϕ$ （公式5)

代码

文本

其中， $Δ$ 为拉普拉斯算子。

代码

文本

1.3 核心困难

代码

文本

"A: 哎？现在只要随着时间或者温度变化，按照 $δ F (ϕ, \nabla ϕ)$ 梯度下降最快的方向，调整 $ϕ$ ,使得F一直下降到最小值，材料组织的的演化过程就得到了哎，问题解决了吗？

B: 然并卵，我们并不知道 $F (ϕ)$ 的形式，这涉及到材料的元素、成分等等复杂的影响，需要考虑各种扩散，边界方程。

A: 不知道形式怎么办？

B: 猜一个吧。

A: 猜的准吗？

B: 不准。

A: 。。。"

代码

文本

自此，广大的材料工作者前赴后继的针对各种场景，提出各种经验模型，希望能够更加符合材料真实的过程。在此不做展开。

代码

文本

二. 晶界长大模拟案例

代码

文本

下面以晶粒长大模型为例子，展示模拟过程。

如下晶粒长大模型使用的能量梯度方程为：

$F (ϕ, \nabla ϕ) = - α ϕ + β ϕ^{3} + 2 γ ϕ (∣ ϕ ∣^{2} - ϕ^{2}) - k Δ ϕ$

代码

文本

模拟代码如下：

代码

文本

[3]

import numpy as np

import matplotlib.pylab as plt

import matplotlib.cm as cm

import random as rd

from matplotlib import animation

def nMP(index, gridsize):

nM1 = index - 1

nM2 = index - 2

nP1 = index + 1

nP2 = index + 2

if index == 0:

nM1 = gridsize - 1

nM2 = gridsize - 2

elif index == 1:

nM2 = gridsize - 1

elif index == gridsize - 1:

nP1 = 0

nP2 = 1

elif index == gridsize - 2:

nP2 = 0

else:

nM1 = index - 1

nM2 = index - 2

nP1 = index + 1

nP2 = index + 2

return nM1, nM2, nP1, nP2

class phaseField():

def __init__(self, gridsize, p, dx, dt, totalTime):

self.gridsize = gridsize

self.p = p

self.dx = dx

self.dt = dt

self.totalTime = totalTime

self.alpha = 1.0

self.beta = 1.0

self.gamma = 1.0

self.kappa = 2.0

self.L = 1.0

self.eta = np.zeros([gridsize, gridsize, p])

self._seeding()

# print self.eta

def _seeding(self):

for i in range(self.p):

x = rd.randint(0, self.gridsize - 1)

y = rd.randint(0, self.gridsize - 1)

p = rd.randint(0, self.p - 1)

self.eta[x, y, p] = 1

def grainGrow(self):

imglist =[]

my_dpi = int(self.gridsize / 5.0)

fig, ax = plt.subplots(1, 1, figsize=(10.5, 8))

for n in range(1, int(self.totalTime / self.dt)):

print('Processing timestep: ', n)

for i in range(self.gridsize):

for j in range(self.gridsize):

iM1, iM2, iP1, iP2 = nMP(i, self.gridsize)

jM1, jM2, jP1, jP2 = nMP(j, self.gridsize)

gradi = 0.5 * (self.eta[iM1, j, :] + self.eta[iP1, j, :] +

self.eta[i, jM1, :] + self.eta[i, jP1, :] -

4.0 * self.eta[i, j, :]) / self.dx ** 2.0 + \

0.25 * ( self.eta[iM2, j, :] + self.eta[iP2, j, :] +

self.eta[i, jM2, :] + self.eta[i, jP2,:] -

4.0 * self.eta[i, j,:]) / self.dx ** 2.0

sumeta = np.sum(self.eta[i, j, :] ** 2.0, axis=-1)

detai = -self.L * (-self.alpha * self.eta[i, j, :] +

self.beta * self.eta[i, j, :] ** 3.0 +

2.0 * self.gamma * self.eta[i, j, :] * (sumeta - self.eta[i, j, :] ** 2) - self.kappa * gradi)

self.eta[i, j, :] = self.eta[i, j, :] + detai * self.dt

grey = np.sum(self.eta ** 2.0, axis=-1)

im = ax.imshow(grey, cmap=cm.Greys_r)

im.axes.get_xaxis().set_visible(False)

im.axes.get_yaxis().set_visible(False)

imglist.append([im])

ani = animation.ArtistAnimation(fig, imglist, interval=200)

ani.save(r'phase2.gif', dpi=my_dpi)

# plt.show()

def main():

p1 = phaseField(gridsize=100, p=400, dx=2, dt = 0.5, totalTime = 40)

p1.grainGrow()

if __name__ == '__main__':

main()

Processing timestep:  1
Processing timestep:  2
Processing timestep:  3
Processing timestep:  4
Processing timestep:  5
Processing timestep:  6
Processing timestep:  7
Processing timestep:  8
Processing timestep:  9
Processing timestep:  10
Processing timestep:  11
Processing timestep:  12
Processing timestep:  13
Processing timestep:  14
Processing timestep:  15
Processing timestep:  16
Processing timestep:  17
Processing timestep:  18
Processing timestep:  19
Processing timestep:  20
Processing timestep:  21
Processing timestep:  22
Processing timestep:  23
Processing timestep:  24
Processing timestep:  25
Processing timestep:  26
Processing timestep:  27
Processing timestep:  28
Processing timestep:  29
Processing timestep:  30
Processing timestep:  31
Processing timestep:  32
Processing timestep:  33
Processing timestep:  34
Processing timestep:  35
Processing timestep:  36
Processing timestep:  37
Processing timestep:  38
Processing timestep:  39
Processing timestep:  40
Processing timestep:  41
Processing timestep:  42
Processing timestep:  43
Processing timestep:  44
Processing timestep:  45
Processing timestep:  46
Processing timestep:  47
Processing timestep:  48
Processing timestep:  49
Processing timestep:  50
Processing timestep:  51
Processing timestep:  52
Processing timestep:  53
Processing timestep:  54
Processing timestep:  55
Processing timestep:  56
Processing timestep:  57
Processing timestep:  58
Processing timestep:  59
Processing timestep:  60
Processing timestep:  61
Processing timestep:  62
Processing timestep:  63
Processing timestep:  64
Processing timestep:  65
Processing timestep:  66
Processing timestep:  67
Processing timestep:  68
Processing timestep:  69
Processing timestep:  70
Processing timestep:  71
Processing timestep:  72
Processing timestep:  73
Processing timestep:  74
Processing timestep:  75
Processing timestep:  76
Processing timestep:  77
Processing timestep:  78
Processing timestep:  79
MovieWriter ffmpeg unavailable; using Pillow instead.

代码

文本

下面为动画展示，可以观察到晶粒形核（随机撒点），晶粒接触形成晶界，晶粒长大晶界合并的过程。

代码

文本

[8]

from IPython.display import Image

from IPython.display import display

display(Image("phase2.gif"))

代码

文本

总结与思考：

相场模拟是一种典型的唯象方法，从现象出发，提炼可能的理论或者经验模型。虽然考虑质量守恒，能量守恒等各种模式，一定程度上能提高这些模型的模拟合理性。但总体上来说，它的精度不高，这也是其发展缓慢的原因。但是，它提供了一种很好的框架，即通过序参数及泛函定义材料结构演化的框架。

对于AI来说，最不怕的就是必需显式的函数形式。通过AI领域的手段，实现泛函的拟合（黑箱模型即可），从而可以使用相场模拟方式，扩展海量的预训练数据。通过这些预训练数据，可以实现对材料组织结构演化的预训练模型。这对于真实数据为小样本的材料领域，可能是一个不错的范式。

代码

文本

相场模拟

已赞2

本文被以下合集收录

机器学习与DFT精华帖

gtang

更新于 2024-09-10

38 篇21 人关注

材料计算

虚白

更新于 2024-08-25

20 篇5 人关注