• {{article.mytitle}}

    • 相场模型介绍

    1. 个人笔记(川大)

    笔记主页
    文献概要

    • NS方程
    • 相场模型
    • 梯度流

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

    Navier-stokes方程

    相场模型(多相流）

    1.2. 模型问题

    ϕ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.3. 模型来源

    1.3.1. 来源

    • 物理模型
    • 自由能或Hamiltomain
      例如
      E[ϕ]=∫Ω12∣∇ϕ∣2+F[ϕ]dΩE[\phi] = \int_{\Omega} \frac{1}{2}|\nabla \phi|^2 + F[\phi] d\Omega

    ϕ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.3.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.4. 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.4.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.4.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.4.3. 该方法的特点

    该方法最早提出来是1993年

    最大缺点：

    • 非线性
    • 很难达到高阶，连二阶都困难(依赖F(ϕ)F(\phi))

    给一个二阶例子的构造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.5. 稳定化方法

    从半隐格式出发
    {ϕ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.6. 拉格朗日乘子法

    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.7. 傅里叶谱方法

    {∂ϕ−Δϕ=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.8. 作业

    1.8.1. 作业1

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

    1.8.2. 作业2

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

    笔记主页