相场模型介绍

Author : zbzhen,        Modified : Tue Feb 21 20:05:00 2023

1. 个人笔记(川大)

笔记主页
文献概要

主要讲解梯度流或守恒式方程
例如： Allen-Cahn: 或者 Cahn-Hilliard

Navier-stokes方程

相场模型(多相流）

1.1. 模型问题

ϕt+Aϕ+N(ϕ)=0\phi_t + A\phi + N(\phi)=0

方法一：需要解一串方程组
(ϕk)t+Aϕk+N(ϕ1,ϕ2，⋯ )=0(\phi_k)_t + A\phi_k + N(\phi_1, \phi_2，\cdots)=0
这种方法无法面对多维问题。

方法二：对时间离散

ϕn+1−ϕn∂t+Aϕn+1+N(ϕn+1)=0\frac{\phi^{n+1} - \phi^n}{\partial t}+ A\phi^{n+1} + N(\phi^{n+1})=0

1.2. 模型来源

1.2.1. 来源

ϕt=−GδE[ϕ]δϕ\phi_t = -\mathcal G \frac{\delta E[\phi]}{\delta\phi}

根据链式法则
Et=(δE[ϕ]δϕ,ϕt)E_t = (\frac{\delta E[\phi]}{\delta\phi}, \phi_t )

G\mathcal G是一个正定算子,则是梯度流，此时
Et=−(δE[ϕ]δϕ,GδE[ϕ]δϕ)≤0E_t = - (\frac{\delta E[\phi]}{\delta\phi}, \mathcal G\frac{\delta E[\phi]}{\delta\phi} ) \leq 0

G\mathcal G反对称，则是Hamiltomain系统，此时
Et=0E_t = 0

这里变分导数的定义

ddϵE[ϕ+ϵψ]∣ϵ=0=(δEδϕ,ψ)\frac{d}{d \epsilon} E[\phi + \epsilon \psi] |_{\epsilon=0}= (\frac{\delta E}{\delta \phi}, \psi)
若
E(ϕ)=E~(ϕ,∇ϕ)E(\phi) = \tilde E(\phi , \nabla \phi)
则有
Eϕ=−∇⋅E~ϕ+E~∇ϕE_\phi = - \nabla \cdot\tilde E_\phi + \tilde E_{\nabla \phi}

1.2.2. 例子1

根据G\mathcal G的不同选取，会得到不同的方程

G\mathcal G取正定算子
heat equation: E(ϕ)=∫Ω12∣∇ϕ∣2E(\phi)=\int_{\Omega} \frac{1}{2}|\nabla \phi|^{2} and G=I\mathcal{G}=I

Allen-Cahn: E(ϕ)=∫Ω(12∣∇ϕ∣2+14ϵ2(ϕ2−1)2)E(\phi)=\int_{\Omega}\left(\frac{1}{2}|\nabla \phi|^{2}+\frac{1}{4 \epsilon^{2}}\left(\phi^{2}-1\right)^{2}\right) and G=I\mathcal{G}=I

Cahn-Hilliard: E(ϕ)=∫Ω(ϵ2∣∇ϕ∣2+14ϵ(ϕ2−1)2)E(\phi)=\int_{\Omega}\left(\frac{\epsilon}{2}|\nabla \phi|^{2}+\frac{1}{4 \epsilon}\left(\phi^{2}-1\right)^{2}\right) and
G=−Δ\mathcal{G}=-\Delta

Phase-field crystal: E(ϕ)=∫Ω(14ϕ4+α2ϕ2−∣∇ϕ∣2+12∣Δϕ∣2)E(\phi)=\int_{\Omega}\left(\frac{1}{4} \phi^{4}+\frac{\alpha}{2} \phi^{2}-|\nabla \phi|^{2}+\frac{1}{2}|\Delta \phi|^{2}\right)
and G=−Δ\mathcal{G}=-\Delta

L1L^{1} minimization: E(ϕ)=∫Ω∣∇ϕ∣E(\phi)=\int_{\Omega}|\nabla \phi| and G=I\mathcal{G}=I

G\mathcal G取反对称算子
Nonlinear Schrödinger equation:
E(ϕ)=∫Ω(12∣∇ϕ∣2+12F(∣ϕ∣2))E(\phi)=\int_{\Omega}\left(\frac{1}{2}|\nabla \phi|^{2}+\frac{1}{2} F\left(|\phi|^{2}\right)\right) and G=i\mathcal{G}=\mathrm{i}

一维情形
KDV equation: E(ϕ)=∫Ω(12∣∂xϕ∣2+ϕ3),G=∂xE(\phi)=\int_{\Omega}\left(\frac{1}{2}\left|\partial_{x} \phi\right|^{2}+\phi^{3}\right), \mathcal{G}=\partial_{x}

δEδϕ=−ϕxxx+3ϕ2\frac{\delta E}{\delta \phi} = - \phi_{xxx} + 3\phi^2

则有

ϕt=−∂x(−ϕxxx+3ϕ2)=∂xxxϕ+6ϕϕx\phi_t = - \partial_x(- \phi_{xxx} + 3\phi^2 ) = \partial_{xxx} \phi + 6\phi\phi_x

1.3. Allen-Cahn与Cahn-Hilliard能量稳定格式

引入一个新的变量μ\mu
{ϕt=−Gμμ=Eϕ=−Δϕ+F′(ϕ)\begin{cases} {\phi_t} = -\mathcal G\mu\\ \mu =E_{\phi}= - \Delta \phi + F'(\phi) \end{cases}
边界条件∂ϕ∂n∣∂Ω=∂μ∂n∣∂Ω=0\frac{\partial \phi}{\partial \bm n}|_{\partial \Omega}=\frac{\partial \mu}{\partial \bm n}|_{\partial \Omega}= 0或者周期边界
有
Et=−(Gμ,μ)≤0E_t = -(\mathcal G \mu, \mu) \leq 0

证明数值格式的能量稳定性，原理是从原来格式中做相同的推算

1.3.1. 全隐格式(Crank-Nicotson)：

{ϕn+1−ϕnΔt=−Gμn+1+μn2μn+1+μn2=−Δϕn+1+ϕn2+F(ϕn+1)−F(ϕn)ϕn+1−ϕn\begin{cases} \dfrac{\phi^{n+1} - \phi^n}{\Delta t} = -\mathcal G \dfrac{\mu^{n+1}+\mu^n}{2}\\[1.5em] \dfrac{\mu^{n+1}+\mu^n}{2} = -\Delta \dfrac{\phi^{n+1}+\phi^n}{2} +\dfrac{F(\phi^{n+1}) - F(\phi^n)}{\phi^{n+1} - \phi^n} \end{cases}

不难证出上式具有2阶无条件稳定，但是每次需要二阶非线性方程组，并且只能证明出当Δt\Delta t充分小的时候才能存在唯一解，该方法最早是杜强老师大约1991-1992年提出

1.3.2. 凸分裂方法

对非线性项做分解
F(ϕ)=Fc(ϕ)−Fe(ϕ)F(\phi) = F_c(\phi) - F_e(\phi)
并且FcF_c与FeF_e都是凸函数，则可以构造下面的格式
{ϕn+1−ϕnΔt=−Gμn+1μn+1=−Δϕn+1+Fc′(ϕn+1)−Fe′(ϕn)\begin{cases} \dfrac{\phi^{n+1} - \phi^n}{\Delta t} = -\mathcal G \mu^{n+1}\\[1.5em] \mu^{n+1}= -\Delta \phi^{n+1} + F'_c(\phi^{n+1}) - F'_e(\phi^n) \end{cases}

注：一般而言FcF_c是能量耗散部分，FeF_e是能量递增格式

上式子分别与μn+1\mu^{n+1}与Δt(ϕn+1−ϕn)\Delta t(\phi^{n+1} - \phi^n)做内积

利用数学恒等式

2(a−b,a)=(a,a)−(b,b)+(a−b,a−b)2(a-b, a) = (a, a) - (b,b) + (a-b, a-b)

对FcF_c在ϕn+1\phi^{n+1}做泰勒展开

Fc(ϕn)=Fc(ϕn+1)+Fc′(ϕn+1)(ϕn−ϕn+1)+12Fc′′(ξ)(ϕn−ϕn+1)2F_c(\phi^n) = F_c(\phi^{n+1}) + F'_c(\phi^{n+1})(\phi^n - \phi^{n+1}) + \frac{1}{2}F''_c(\xi)(\phi^n - \phi^{n+1}) ^2
同理对FeF_e在ϕn\phi^{n}做泰勒展开

Fe(ϕn+1)=Fe(ϕn)−Fe′(ϕn+1)(ϕn−ϕn+1)+12Fe′′(η)(ϕn−ϕn+1)2F_e(\phi^{n+1}) = F_e(\phi^{n}) - F'_e(\phi^{n+1})(\phi^n - \phi^{n+1}) + \frac{1}{2}F''_e(\eta)(\phi^n - \phi^{n+1}) ^2

两式子合并起来最后可以证明出该格式也是无条件稳定
缺点是： 只有一阶精度，也要解非线性方程组
优点是： 虽然也有非线性，但是非线性项是凸的，因此是存在唯一的，并且唯一解为凸泛函的能量最小解

证明： 只针对简单的情形G=I\mathcal G = I

ϕn+1−ϕnΔt=Δϕn+1−Fc′(ϕn+1)+Fe′(ϕn)\dfrac{\phi^{n+1} - \phi^n}{\Delta t}= \Delta \phi^{n+1} - F'_c(\phi^{n+1}) + F'_e(\phi^n)

定义Q(ϕ)=∫12Δt∣ϕ∣2+12∣∇ϕ∣2+Fc(ϕ)−gn(ϕ)Q(\phi) = \int\frac{1}{2\Delta t} |\phi|^2 + \frac{1}{2}|\nabla \phi|^2 + F_c(\phi) - g^n(\phi)

其中
gn(ϕ)=1Δtϕn+Fe′(ϕn)g^n(\phi) = \frac{1}{\Delta t} \phi^n + F'_e(\phi^n)

凸泛函QQ的最小值就是变分导数为零的解

δQδϕ∣ϕ=ϕn+1=1Δtϕn+1−Δϕn+1+Fc′(ϕn+1)−gn\frac{\delta Q}{\delta \phi}|_{\phi = \phi^{n+1}} = \frac{1}{\Delta t} \phi^{n+1} - \Delta \phi^{n+1} + F'_c(\phi^{n+1}) - g^n

1.3.3. 该方法的特点

该方法最早提出来是1993年

最大缺点：

给一个二阶例子的构造F(ϕ)=(ϕ2−1)2F(\phi) = (\phi^2 -1)^2

取Fc(ϕ)=14(ϕ4+1),Fe(ϕ)=12ϕ2F_c(\phi) = \frac{1}{4}(\phi^4+1), \quad F_e(\phi)=\frac{1}{2}\phi^2

则有Fc′=ϕ3,Fe′=ϕF'_c = \phi^3,\quad F'_e = \phi

{ϕn+1−ϕnΔt=−Gμn+1+μn2μn+1+μn2=−Δϕn+1+ϕn2+A−12(3ϕn−ϕn−1)\begin{cases} \dfrac{\phi^{n+1} - \phi^n}{\Delta t} = -\mathcal G \dfrac{\mu^{n+1}+\mu^n}{2}\\[1.5em] \dfrac{\mu^{n+1}+\mu^n}{2} = -\Delta \dfrac{\phi^{n+1}+\phi^n}{2}+ A - \dfrac{1}{2}(3\phi^n - \phi^{n-1}) \end{cases}

其中
A=(ϕn+1)2−(ϕn)22ϕn+1+ϕn2A = \frac{(\phi^{n+1})^2 - (\phi^n)^2}{2} \frac{\phi^{n+1}+\phi^n}{2}

证明的过程中需要注意下面的拼凑

(3ϕn−ϕn−1,ϕn+1−ϕn)=(ϕn+1+ϕn−(ϕn+1−2ϕn+ϕn−1),ϕn+1−ϕn)=(ϕn+1+ϕn−(ϕn+1−ϕn)+(ϕn−ϕn−1)),ϕn+1−ϕn)\begin{aligned} &(3\phi^n - \phi^{n-1}, \phi^{n+1} - \phi^n)\\ =&(\phi^{n+1}+\phi^n - (\phi^{n+1}-2\phi^n + \phi^{n-1}),\phi^{n+1}-\phi^n)\\ =&(\phi^{n+1}+\phi^n - (\phi^{n+1}-\phi^n) + (\phi^n-\phi^{n-1})),\phi^{n+1}-\phi^n) \end{aligned}

于是得到修正的能量是递减的，其中修正的能量定义为
E~(ϕn)=12∣∣∇ϕ∣∣2+14(∣∣ϕn+1∣∣4−2∣∣ϕn+1∣∣2)+14∣∣ϕn+1−ϕn∣∣\tilde E(\phi^n) = \frac{1}{2}||\nabla \phi||^2 + \frac{1}{4}(||\phi^{n+1}||^4- 2||\phi^{n+1}||^2) + \frac{1}{4}|| \phi^{n+1}-\phi^n||

1.4. 稳定化方法

从半隐格式出发
{ϕn+1−ϕnΔt=−Gμn+1μn+1=−Δϕn+1+F′(ϕn)\begin{cases} \dfrac{\phi^{n+1} - \phi^n}{\Delta t} = -\mathcal G \mu^{n+1}\\[1.5em] \mu^{n+1} = -\Delta \phi^{n+1} + F'(\phi^n) \end{cases}

优点： 计算快
缺点： 不稳定

改进办法: 加上稳定项

{ϕn+1−ϕnΔt=−Gμn+1μn+1=−Δϕn+1+F′(ϕn)\begin{cases} \dfrac{\phi^{n+1} - \phi^n}{\Delta t} = -\mathcal G \mu^{n+1}\\[1.5em] \mu^{n+1} = -\Delta \phi^{n+1} + F'(\phi^n) \end{cases}

能量修正
E(ϕ)=∫12∣∇ϕ∣2+S∣ϕ∣2+F(ϕ)−S∣ϕ∣2dxE(\phi) = \int \frac{1}{2}|\nabla \phi|^2 + S|\phi|^2 + F(\phi) - S|\phi|^2 dx
或者
E(ϕ)=∫12∣(∇ϕ∣2+S∣ϕ∣2)−(S∣ϕ∣2−F(ϕ))dxE(\phi) = \int \frac{1}{2}|(\nabla \phi|^2 + S|\phi|^2 )-( S|\phi|^2-F(\phi))dx
取Fe(ϕ)=S∣ϕ∣2−F(ϕ)F_e(\phi) =S|\phi|^2 - F(\phi)
于是
Fe′′(ϕ)=2S−F′′(ϕ)F''_e(\phi) = 2S - F''(\phi)
如果∣F′′(ϕ)∣≤L|F''(\phi)|\leq L，取S≥L2S\geq \frac{L}{2}，则是特殊凸分裂格式，因此可以无条件稳定

早期能量是下面的格式
F(ϕ)=(1+ϕ)ln⁡(1+ϕ)−(1−ϕ)ln⁡(1−ϕ)+θϕ2F(\phi) = (1+\phi)\ln(1+\phi) - (1-\phi)\ln(1-\phi) + \theta \phi^2
物理上只关注ϕ∈(−1,1)\phi\in(-1,1),

因此可以对能量修正，取
F~(ϕ)={14(ϕ2−1)2,∣ϕ∣≤1F(ϕ),∣ϕ∣>1.\tilde F(\phi)= \begin{cases} \frac{1}{4}(\phi^2 - 1)^2 , \quad & |\phi|\leq 1\\ F(\phi),\quad &| \phi| >1. \end{cases}

但是这种方法无法达到二阶

其实是现有了稳定化方法，再有凸分离方法，本质上这两种方法是有很大联系的

1.5. 拉格朗日乘子法

F(ϕ)=14(ϕ2−1)2F(\phi) = \frac{1}{4}(\phi^2 - 1)^2

引入q=ϕ2q = \phi^2，则F(ϕ)=14q2F(\phi) = \frac{1}{4}q^2

则

{ϕt=−Gμμ=Eϕ=−Δϕ+qϕqt=2ϕϕt\begin{cases} \phi_t = -\mathcal G\mu \\ \mu = E_\phi = -\Delta \phi + q\phi\\ q_t = 2\phi\phi_t \end{cases}

数值格式为
{ϕn+1−ϕn=−ΔtGμn+1μn+1=−Δϕn+1+qn+1ϕnqn+1−qn=2ϕn(ϕn+1−ϕn)\begin{cases} \phi^{n+1} - \phi^n =-\Delta t\mathcal G \mu^{n+1} \\ \mu^{n+1} = -\Delta \phi^{n+1} + q^{n+1}\phi^n\\ q^{n+1} - q^{n} = 2\phi^n(\phi^{n+1} - \phi^n) \end{cases}
上式分别与μn+1,ϕn+1−ϕn,12qn+1\mu^{n+1}, \phi^{n+1}-\phi^n, \frac{1}{2}q^{n+1}做内积再联立可得修正的能量稳定，其中修正的能量定义为

E~(ϕn+1)=12∣∇ϕn+1∣2+12∣∣qn+1∣∣2\tilde E({\phi^{n+1}}) = \frac{1}{2}|\nabla \phi^{n+1}|^2 + \frac{1}{2}||q^{n+1}||^2

优点：

结论：
做守恒型方程，用CN格式比较好，因为它没有多余的耗散项目
做耗散型方程，用BDF格式比较好

二阶格式简记：
对于CN格式：

对ϕt\phi_t项，改写成
ϕn+1−ϕnΔt\frac{\phi^{n+1}- \phi^n}{\Delta t}
对于ϕ\phi项，改写成
3ϕn−ϕn−12\frac{3\phi^{n} - \phi^{n-1}}{2}

对于BDF格式：
对ϕt\phi_t项，改写成
3ϕn+1−4ϕn+ϕn−12Δt\frac{3\phi^{n+1}- 4\phi^n+\phi^{n-1}}{2\Delta t}
对于ϕ\phi项，改写成
2ϕn−ϕn−12\phi^{n} - \phi^{n-1}

不同的地方引入的辅助变量qq，CN格式12(qn+1+qn)\frac{1}{2}(q^{n+1}+q^n)，BDF格式为qn+1q^{n+1}

注：拉格朗日乘子法是从方程出发的，直接对方程进行修正

1.6. 傅里叶谱方法

{∂ϕ−Δϕ=fϕ∣∂Ω=0,or ∂ϕ∂n∣=0,或周期边界\begin{cases} \partial \phi - \Delta \phi = f\\ \phi|_{\partial \Omega}=0, \text{or } \frac{\partial \phi}{\partial n}| = 0, \text{或周期边界} \end{cases}
如果用傅里叶普方法
u=∑ukjeikxeijyu = \sum u_{kj}e^{ikx}e^{ijy}
于是
−Δu=∑(k2+j2)ukjeikxeijy-\Delta u = \sum(k^2 + j^2)u_{kj}e^{ikx}e^{ijy}

1.7. 作业

1.7.1. 作业1

证明拉格朗日乘子法BDF-2格式的稳定性

1.7.2. 作业2

写出Phase-field crystal的凸分裂格式格式

笔记主页