数学物理学报, 2023, 43(4): 1255-1268

Cahn-Hilliard方程的自适应间断有限体积元法

曾纪尧, 李剑,*

陕西科技大学 数学与数据科学学院 西安 710021

An Adaptive Discontinuous Finite Volume Element Method for the Cahn-Hilliard Equation

Zeng Jiyao, Li Jian,*

School of Mathematics and Data Science, Shaanxi University of Science and Technology, Xi'an 710021

通讯作者: *李剑, E-mail: jianli@sust.edu.cn

收稿日期: 2022-04-6   修回日期: 2022-09-13  

基金资助: 国家自然科学基金(11771259)
陕西省人工智能联合实验室(2022JC-SYS-05)
陕西省教育厅创新团队项目(21JP013)
部分资助和陕西省自然科学基础研究计划重点项目(2023-JC-ZD-02)

Received: 2022-04-6   Revised: 2022-09-13  

Fund supported: NSF of China(11771259)
Shaanxi Provincial Joint Laboratory of Artificial Intelligence(2022JC-SYS-05)
Innovative Team Project of Shaanxi Provincial Department of Education(21JP013)
Shaanxi Province Natural Science Basic Research Program Key Project(2023-JC-ZD-02)

摘要

Cahn-Hilliard 方程是一类重要的四阶非线性扩散方程, 具有丰富的物理背景和深刻的研究价值, 但方程的非线性势项 $ f(u) $ 的存在以及小参数 $ \epsilon $ 导致的强刚性会给数值模拟带来诸多挑战, 因此设计高效、准确的数值方案满足方程离散能量定律是非常重要的. 间断有限体积元方法 (DFVEM) 采用低阶元, 具有精度高、操作简单、适合工程应用的网格自适应等优点. 该文对Cahn-Hilliard 方程利用 DFVEM 结合全隐格式进行求解, 证明了全离散格式质量守恒和能量耗散的重要理论结果. 数值实验提出一种自适应时间步进策略, 验证了方法的有效性.

关键词： Cahn-Hilliard 方程; 间断有限体积元法; 离散能量耗散

Abstract

The Cahn-Hilliard equation is an important class of fourth-order nonlinear diffusion equations with rich physical background and profound research value. In the numerical simulation, the existence of the nonlinear potential term $ f(u) $ of the equation and the strong rigidity caused by small parameter $ \epsilon $ will bring many challenges, so it is important to design efficient and accurate numerical schemes to satisfy the discrete energy law of the equations. In this paper, the Cahn-Hilliard equation is solved by the discontinuous finite volume element method (DFVEM) combined with the fully implicit scheme, and the important theoretical results of mass conservation and energy dissipation in the fully discrete scheme are proved. At the same time, the semi-discrete format error estimates are given. Finally, numerical experiments propose an adaptive time stepping strategy and verify the effectiveness of the method.

Keywords： The Cahn-Hilliard equation; Discontinuous finite volume element method; Discrete energy dissipation

PDF (2200KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

曾纪尧, 李剑. Cahn-Hilliard方程的自适应间断有限体积元法[J]. 数学物理学报, 2023, 43(4): 1255-1268

Zeng Jiyao, Li Jian. An Adaptive Discontinuous Finite Volume Element Method for the Cahn-Hilliard Equation[J]. Acta Mathematica Scientia, 2023, 43(4): 1255-1268

1 引言

Cahn-Hilliard 方程由 Cahn 和 Hilliard[2]在 1958 年研究热力学两相物质之间的相互扩散现象时提出. 本文考虑如下的 Cahn-Hilliard 方程

$\begin{matrix} \left\{ \begin{array}{l} \frac{{\partial u}}{{\partial t}} = \Delta w,uad (x,t) \in \Omega \times (0,T),\\ w = f\left( u \right)- {{{\epsilon ^{{\rm 2}}}}}\Delta u,\\ {u\left( {x,0} \right) = {u_0}\left( x \right),uad x \in \Omega,}\\ {{\partial _n}u\left| {_{\partial \Omega }} \right. = 0,uad {\partial _n}w\left| {_{\partial \Omega }} \right. = 0.} \end{array} \right. \end{matrix}$

其中 $ \Omega \subset {R^2} $ 是具有 Lipschitz 边界的有界凸域, $ n $表示 $ \partial \Omega $ 外法向量. $ T $ 是正实数, 参数 $ \epsilon>0 $ 代表界面宽度, 而 $ u(x,t) $ 是相位函数. $ f\left( u \right) = F'\left( u \right),F\left( u \right) = \frac{1}{4}{\left( {{u^2} - 1} \right)^2} $ 是双阱势函数. 在本文中, $ f(u) $ 满足以下条件[10]: 存在一个常数 $ C_0 $ 使得

$\begin{matrix} \mathop {\max }\limits_{u \in R} \left| {f'\left( u \right)} \right| \le {C_0}. \end{matrix} $

作为描述相场模型的最基本方程之一, Cahn-Hilliard 方程具有重要的研究意义和应用价值, 如材料科学[11],图像分析[6]和相变[4]等领域. 国内外学 者对Cahn-Hilliard 方程提出许多有效的数值解法, 如Zhang 等[14]构造无条件 能量稳定的有限差分格式, 提出了一种自适应时间步进策略, 这种策略避免了在数值模拟中使用大时间步长时可能导致的精度损失. 叶[13]利用 Legendre 拟谱法对Cahn-Hilliard 方程进行数值求解, 证明了离散解的存在唯一性, 并给出了最佳误差估计. Feng等[5]提 出全离散有限元格式, 并建立先验误差估计, 文中给出的全离散格式是最优阶的. 文献[13]介绍了一种求解 Cahn-Hilliard 方程的间断有限元方法, 这种方法在不诉诸任何迭代的情况下进行有效地求解, 并在效率、准确性和保持所需解的属性方面具有良好性能. 文献[8]利用间断有限体积元方法求解耦合 Navier-Stokes-Cahn-Hilliard 模型, 并对其质量守恒和能量耗散特性进行了分析, 该数值方法高效且易于处理具有不同类型边界条件的复杂几何形状. 在对 Cahn-Hilliard 方程进行数值模拟时, 间断有限体积元方法相较于其他方法,具有空间结构简单, 在时间相关问题中使用块对角质量矩阵等优点, 同时使用到的低阶元, 更适用于做网格灵活的自适应算法, 时间步长根据能量的演变而变化, 这样更好的解决了分割元素的自由度过大导致的计算量问题. 为了方法具有更好的稳定性, 时间离散选择全隐格式.

本文结构如下:第2节介绍 Cahn-Hilliard 方程间断有限体积元空间离散格式; 第3节结合全隐格式进行时间离散得到全离散格式, 并严格证明了质量守恒性质和离散能量耗散定律; 第5节通过数值算例验证所给格式的可行性和有效性. 最后, 第6节给出了结论.

2 间断有限体积元法

在区域 $ \Omega $ 上引进原始剖分和对偶剖分. 令原始剖分 $ {\mathbb{R}e _h} $ 是 一个三角剖分, 且 $ h =\max \limits_{K \in {\mathbb{R}e _h}} {h_K} $, 其中 $ h_K $ 是单元 $ K $ 的直径. 每个元素 $ K $ 在 $ {\mathbb{R}e _h} $ 通过将三角形的重心连接到每个顶点被分成三个子三角形. 这样的三角形 $ T_j $ 构成了原始剖分 $ {\mathbb{R}e _h} $ 的对偶剖分 $ {{\cal T}_h} $ (见图1).

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   左:原始剖分(实线部分)和对偶剖分(虚线部分); 右:元素$ K $及其对偶元素$ T_j $


定义与相场变量 $ u $ 相关的有限维试探函数空间 $ V_h $ 和有限维检验函数空间 $ V^*_h $

$ \begin{eqnarray*}&& {V_h} = {{\rm }}\left\{ {u \in {L^2}\left( \Omega \right):u\left| {_K} \right. \in {{\cal P}_1}\left( K \right),uad \forall K \in {\mathbb{R}e _h}} \right\}, \\ && {V^*_h} = \left\{ {u \in {L^2}\left( \Omega \right):u\left| {_T} \right. \in {{\cal P}_0}\left( T \right),uad \forall T \in {{\cal T}_h}} \right\}. \end{eqnarray*}$

定义与化学势变量 $ w $ 相关的有限维试探函数空间 $ P_h $ 和有限维检验函数空间 $ P^*_h $

$ \begin{eqnarray*} && {P_h} = {{\rm }}\left\{ {w \in {L^2}\left( \Omega \right):w\left| {_K} \right. \in {{\cal P}_1}\left( K \right),uad \forall K \in {\mathbb{R}e _h}} \right\},\\&& {P^*_h} = \left\{ {w \in {L^2}\left( \Omega \right):w\left| {_T} \right. \in {{\cal P}_0}\left( T \right),uad \forall T \in {{\cal T}_h}} \right\}, \end{eqnarray*}$

式中, $ {{\cal P}_l}\left( T \right) $ 由 $ T $ 上所有次数小于或等于$ l $ 的多项式组成.

设 $ e $ 为 $ {\mathbb{R}e _h} $ 中两个元素 $ K_1 $ 和 $ K_2 $ 的公共内边, 设 $ \textbf{n}_1 $ 和 $ \textbf{n}_2 $ 分别为指向 $ K_1 $ 和 $ K_2 $ 的单位外法向量. 在 $ e $ 上定义标量 $ q $ 的均值 $ \left\{ \cdot \right\} $ 和跳跃 $ [ \cdot] $

$\left\{ q \right\} = \frac{1}{2}\left( {q\left| {_{\partial {K_1}}} \right. + q\left| {_{\partial {K_2}}} \right.} \right),uad \left[ q \right] = q\left| {_{\partial {K_1}}} \right.{{\bf{n}}_1} + q\left| {_{\partial {K_2}}} \right.{{\bf{n}}_2}.$

如果 $ e $ 是 $ \Omega $ 的边界, 则定义

$\left\{q \right\} = q,uad \left[ q \right] = q \cdot \textbf{n}.$

令 $ V(h) = {V_h} + {H^2}\left( \Omega \right) $, $ P(h) = {P_h} + {H^2}\left( \Omega \right) $, 定义一个映射 $ \gamma: V\left( h \right) \to {V^*_h}, P\left( h \right) \to {P^*_h} $

$\gamma u\left| {_T} \right. = \frac{1}{{{h_e}}}\int_e {u\left| {_T} \right.{{\rm d}}s, uad \gamma w\left| {_T} \right. = \frac{1}{{{h_e}}}\int_e {w\left| {_T} \right.{\rm d}s,} uad {{\rm }}T \in {{\cal T}_h}},$

其中 $ h_e $ 是边 $ e $ 的长度.

引理 2.1 对于 $ {u_h} \in {V_h}, K \in {\mathbb{R}e _h} $, 映射 $ \gamma $ 具有以下属性

$ \begin{matrix} \begin{array}{ll} \int_K ( {u_h} - \gamma {u_h}){{\rm d}}x = 0, & \int_e ( {u_h} - \gamma {u_h}){{\rm d}s = 0,} \\[3mm] {\left\| {\gamma {u_h} - {u_h}} \right\|_K} \le C{h_K}{\left| {{u_h}} \right|_{1,T}}, & \left[ {\gamma {u_h}} \right] = \frac{1}{{{h_e}}}\int_e {\left[ {{u_h}} \right]} {{\rm d}s,} \\[2mm] {\left\| {\left[ {\gamma {u_h}} \right]} \right\|_e} \le {\left\| {\left[ {{u_h}} \right]} \right\|_e},& {\left\| {\gamma {u_h}} \right\|_0} = {\left\| {{u_h}} \right\|_0}. \end{array} \end{matrix}$

对于 $ w_h \in {P_h} $, 这些结论也是成立的.

算子 $ \gamma $ 相对于 $ L^2 $ 内积是自伴的, 即

$ \begin{matrix} \left( {{u_h},\gamma {\psi _h}} \right) = \left( {{\psi _h}, \gamma {u_h}} \right),uad \forall {u_h},{\psi _h} \in {V_h}. \end{matrix}$

定义分析中使用的 $ V(h) $ 和 $ P(h) $ 的离散规范

$\begin{eqnarray*}&& \left\| {\left| u \right|} \right\|_h^2 = \sum\limits_{K \in {\mathbb{R}e _h}} {\left| u \right|_{1,K}^2 + } \sum\limits_{e \in \Gamma } {h_e^{ - \beta }} \left\| {\left[ u \right]} \right\|_e^2, \\&& \left\| {\left| w \right|} \right\|_h^2 = \sum\limits_{K \in {\mathbb{R}e _h}} {\left| w \right|_{1,K}^2 + } \sum\limits_{e \in \Gamma } {h_e^{ - \beta }} \left\| {\left[ w \right]} \right\|_e^2. \end{eqnarray*}$

定义 $ {\left\| {\left| {{v_h}} \right|} \right\|_0} = \left( {{v_h},\gamma {v_h}} \right) $, 则 $ {\left\| \cdot \right\|_0} $ 和 $ {\left\| {\left| \cdot \right|} \right\|_0} $ 是等价的, 即存在与网格尺寸无关的正常数 $ C_1 $ 和 $ C_2 $, 使得

${C_1}{\left\| {\left| {{v_h}} \right|} \right\|_0} \le {\left\| {{v_h}} \right\|_0} \le {C_2}{\left\| {\left| {{v_h}} \right|} \right\|_0},uad\forall {v_h} \in {V_h}.$

将(1.1)式第一个式子两端同时乘上 $ \gamma {\psi _h} \in V_h^* $, 第二个式子两端同时乘上 $ \gamma {\chi _h} \in P_h^* $, 并通过格林公式我们的到

$ \begin{equation} \left( {\frac{{\partial u}}{{\partial t}},\gamma {\psi _h}} \right) - \sum\limits_{T \in {{\cal T}_h}} {\int_{\partial T} {\nabla w\cdot{\bf{n}}\gamma {\psi _h}} {{\rm d}s}} {{\rm = 0,}} \end{equation}$
$\begin{equation} - {{{\epsilon ^{{\rm 2}}}}}\sum\limits_{T \in {{\cal T}_h}} {\int_{\partial T} {\nabla u\cdot{\bf{n}}\gamma {\chi _h}} {{\rm d}}s} {{\rm + }}\left( {f\left( u \right),\gamma {\chi _h}} \right) = \sum\limits_{T \in {{\cal T}_h}} {\int_{T} {w\cdot\gamma {\chi _h}} {{\rm d}}x}. \end{equation}$

设 $ \Gamma $ 表示三角形 $ K $ 边界的并集同时令 $ {\Gamma _0}: = \Gamma \backslash \partial \Omega $. 然后, 我们得到

$\begin{matrix} \sum\limits_{T \in {{\cal T}_h}} {\int_{\partial T} {\nabla u \cdot {\bf{n}}\gamma {v_h}{\rm d}s} } &= &\sum\limits_{K \in {\mathbb{R}e _h}} {\sum\limits_{i = 1}^3 {\int_{\partial {T_i}} {\nabla u \cdot {\bf{n}}\gamma {v_h}{\rm d}s} } } \\ &= &\sum\limits_{K \in {\mathbb{R}e _h}} {\sum\limits_{i = 1}^3 {\int_{{J_{i + 1}}Q{J_i}} {\nabla u \cdot {\bf{n}}\gamma {v_h}{\rm d}s} } } + \sum\limits_{K \in {\mathbb{R}e _h}} {\int_{\partial K} {\nabla u \cdot {\bf{n}}\gamma {v_h}{\rm d}s} }, \end{matrix}$

在方程中 $ J_4=J_1 $. 然后通过简单的计算给出

$\begin{matrix} \sum\limits_{K \in {\mathbb{R}e _h}} {\int_{\partial K} {qv \cdot {\bf{n}}{{\rm d}s = }\sum\limits_{e \in \Gamma } {\int_e {\left[ q \right]} } } } \cdot \left\{ v \right\}{{\rm d}s + }\sum\limits_{e \in {\Gamma _{{\rm 0}}}} {\int_e {\left\{ q \right\}} } \cdot \left[ v \right]{{\rm d}s}. \end{matrix}$

通过使用(2.6)式并结合 $ {\left[ {\nabla u} \right]_e} = 0, {\left[ {\nabla w} \right]_e} = 0,$$ \forall e \in {\Gamma _0} $, 可以得出

$\sum_{T \in T_{h}} \int_{\partial T} \nabla u \cdot \mathbf{n} \gamma \chi_{h} \mathrm{~d} s=\sum_{K \in \Re_{h}} \sum_{i=1}^{3} \int_{J_{i+1} Q J_{i}} \nabla u \cdot \mathbf{n} \gamma \chi_{h} \mathrm{~d} s+\sum_{\mathrm{e} \in \Gamma} \int_{e}\left[\gamma \chi_{h}\right]\{\nabla u\} \mathrm{d} s$
$\sum_{T \in T_{h}} \int_{\partial T} \nabla w \cdot \mathbf{n} \gamma \psi_{h} \mathrm{~d} s=\sum_{K \in \Re_{h}} \sum_{i=1}^{3} \int_{J_{i+1} Q J_{i}} \nabla w \cdot \mathbf{n} \gamma \psi_{h} \mathrm{~d} s+\sum_{\mathrm{e} \in \Gamma} \int_{e}\left[\gamma \psi_{h}\right]\{\nabla w\} \mathrm{d} s$.

那么(1.1)式的弱格式为对 $ \forall \left( {{\psi _h},{\chi _h}} \right) \in \left( {{V_h} \times {P_h}} \right) $, 我们有

$\begin{equation} \left( {\frac{{\partial u}}{{\partial t}},\gamma {\psi _h}} \right) - \sum\limits_{K \in {\mathbb{R}e _{{ h}}}} {\sum\limits_{i = 1}^{{\rm 3}} {\int_{{J_{i + 1}}Q{J_{i }}} {\nabla w\cdot{\bf{n}}\gamma {\psi _h}{{\rm d}s} - \sum\limits_{{{\rm e}} \in \Gamma } {\int_e {\left[ {\gamma {\psi _h}} \right]} } } } } \left\{ {\nabla w} \right\}{{\rm d}s = 0,} \end{equation}$
$\begin{equation} -{{{\epsilon ^2}}} \sum\limits_{K \in {\mathbb{R}e _{{ h}}}} {\sum\limits_{i = 1}^{{\rm 3}} {\int_{{J_{i + 1}}Q{J_{i }}} {\nabla u\cdot{\bf{n}}\gamma {\chi _h}{{\rm d}s}} - {{{\epsilon ^2}}}\sum\limits_{{{\rm e}} \in \Gamma } {\int_e {\left[ {\gamma {\chi _h}} \right]} } } } \left\{ {\nabla u} \right\}{{\rm d}s + }\left( {f\left( u \right),\gamma {\chi _h}} \right){{\rm = }}\left( {w,\gamma {\chi _h}} \right). \end{equation}$

结合(2.3)-(2.10)式, Cahn-Hilliard 方程间断有限体积元方法的空间离散格式为: 对 $ \forall \left( {{\psi _h},{\chi _h}} \right) \in \left( {{V_h} \times {P_h}} \right) $, 使得

$ \begin{equation}\label{DFVEM1} \left( {\frac{{\partial {u_h}}}{{\partial t}},\gamma {\psi _h}} \right) + {A}\left( {{w_h},\gamma {\psi _h}} \right) = 0, \end{equation}$
$\begin{equation}\label{DFVEM2} {{{\epsilon ^2}}}{A}\left( {{u_h},\gamma {\chi _h}} \right) + \left( {f\left( u \right),\gamma {\chi _h}} \right){{\rm = }}\left( {w,\gamma {\chi _h}} \right), \end{equation}$

其中

$\begin{matrix} A\left( {{u_h},\gamma {v_h}} \right)& = &- \sum\limits_{K \in {\mathbb{R}e _{{ h}}}} {\sum\limits_{i = 1}^3 {{\int _{{J_{i + 1}}Q{J_{i }}}}\nabla {u_h}\cdot{\bf{n}}\gamma {v_h}{{\rm d}s}} } - \sum\limits_{{{\rm e}} \in \Gamma } {\int_e {\left\{ {\nabla {u_h}} \right\}\left[ {\gamma {v_h}} \right]} } {{\rm d}s} \\ &&+ \theta \sum\limits_{e \in \Gamma } {\int_e {\left\{ {\nabla {v_h}} \right\} \left[ {\gamma {u_h}} \right]} } {{\rm d}}s + \sum\limits_{e \in \Gamma } {\frac{\alpha }{{h_e^\beta }}\int_e {\left[ {{u_h}} \right]\left[ {{v_h}} \right]} }{\rm d}s, \end{matrix}$

对于 $ A(\cdot,\cdot) $ 项, 参数 $ \theta=-1,0,1 $ 将导致对称、不完全和非对称内部惩罚方法, 分别用 SIPG、IIPG 和 NIPG 表示. 需要强调的是, 接下来给出的 Cahn-Hilliard 方程的全离散格式, 我们考虑 $ \theta=-1 $ 的情况.

假设 $ u $ 是(1.1)式的解, 则由 $ {\left[ {\gamma u} \right]_e} = 0 $, 可知 $ u $ 满足

$\begin{equation} \begin{array}{l} \left( {\frac{{\partial u}}{{\partial t}},\gamma {\psi _h}} \right) + A\left( {w,\gamma {\psi _h}} \right) = 0, \\[2mm] {\epsilon ^2}A\left( {u,\gamma {\chi _h}} \right) + \left( {f\left( u \right),\gamma {\chi _h}} \right) - \left( {w,\gamma {\chi _h}} \right) = 0. \\ \end{array} \end{equation}$

3 全离散格式

接下来, 我们再结合后向欧拉格式进行时间离散, 得到全离散格式. 令 $ u_h^j = {u_h}\left( {{t^j}} \right)$, $ {t^j} = j\Delta t $, $ \Delta t = \frac{T}{N} $, 则 Cahn-Hilliard 方程的全离散间断有限体积元格式为

$\begin{equation} \begin{array}{l} \bigg (\frac{{u_h^j - u_h^{j - 1}}}{{\Delta t}},\gamma {\psi _h}\bigg) + {A}\left( {w_h^j,\gamma {\psi _h}} \right) = 0,uad \forall {\psi _h} \in {V_h},uad \\[3mm] {{{\epsilon ^2}}}{A}\left( {u_h^j,\gamma {\chi _h}} \right) +\left( {f\left( {u_h^j} \right)\gamma {\chi _h}} \right) = \left( {w_h^j,\gamma {\chi _h}} \right),uad \forall {\chi _h} \in {P_h}. \end{array} \end{equation}$

初始条件为 $ u_h^0 = {u_{0,h}} $, 并且 $ f_h^j = {\left( {u_h^j} \right)^3} - {u_h^{j-1}} $. 我们定义向后差分算子

$ \begin{equation}\label{BDQ} {d_t}u_h^j=\frac{{u_h^j - u_h^{j - 1}}}{{\Delta t}}. \end{equation}$

为了证明全离散格式的能量稳定性, 首先给出下面的等价关系, 更多细节可以在文献[1,7,12]中得到.

定理 3.1$ \forall u,\psi \in V\left( h \right)\ $, $ \forall w,\chi \in P\left( h \right) $, 我们有

$\begin{aligned}-\sum_{K \in \Re_{h}} \sum_{i=1}^{3} \int_{J_{i+1} Q J_{i}} \nabla w \cdot \mathbf{n} \gamma \psi \mathrm{d} s= & (\nabla w, \nabla \psi)+\sum_{K \in \Re_{h}} \int_{\partial K} \nabla w \cdot \mathbf{n}(\gamma \psi-\psi) \mathrm{d} s \\ & -\sum_{K \in \Re_{h}}(\nabla \cdot \nabla w, \gamma \psi-\psi)_{K},\end{aligned}$
$\begin{aligned}-\sum_{K \in \Re_{h}} \sum_{i=1}^{3} \int_{J_{i+1} Q J_{i}} \nabla u \cdot \mathbf{n} \gamma \chi \mathrm{d} s= & (\nabla u, \nabla \chi)+\sum_{K \in \Re_{h}} \int_{\partial K} \nabla u \cdot \mathbf{n}(\gamma \chi-\chi) \mathrm{d} s \\ & -\sum_{K \in \Re_{h}}(\nabla \cdot \nabla u, \gamma \chi-\chi)_{K} \cdot\end{aligned}$

如果 $ \forall u_h,{\psi}_h \in V\left( h \right) $, $ \forall w_h,{\chi}_h \in P\left( h \right) $, 我们得到

$-\sum_{K \in \Re_{h}} \sum_{i=1}^{3} \int_{J_{i+1} Q J_{i+1}} \nabla w_{h} \cdot \mathbf{n} \gamma \psi_{h} \mathrm{~d} s=\left(\nabla w_{h}, \nabla \psi_{h}\right)$
$-\sum_{K \in \Re_{h}} \sum_{i=1}^{3} \int_{J_{i+1} Q J_{i+1}} \nabla u_{h} \cdot \mathbf{n} \gamma \chi_{h} \mathrm{~d} s=\left(\nabla u_{h}, \nabla \chi_{h}\right)$.

接下来将证明时间离散质量守恒.

定理 3.2 令$ \left\{ {\big({u_h^j},{w_h^j} \big)} \right\}_{j = 1}^N\ $ 是离散系统(3.1)式的解, 那么这个系统有下面的质量守恒定律

$ \begin{equation} \int_\Omega {u_h^j} {{\rm d}}x = \int_\Omega {{u_0}} {{\rm d}}x,uad j \ge 1. \end{equation}$

证 令(3.1)式中的 $ {\psi _h}=1 $, 然后使用引理2.1和定理3.1,我们得到

$\begin{equation} \int_\Omega {u_h^j} {{\rm d}}x = \int_\Omega {u_h^{j - 1}} {{\rm d}}x,uad j \ge 1, \end{equation}$

这样就完成了定理的证明.证毕.

下面, 将验证 Cahn-Hilliard 系统的能量耗散性质. 首先定义相应的离散能量泛函

$\begin{equation}E\left( {{t^j};u_h^j} \right) = \frac{{{\epsilon ^2}}}{2}\left| {\left\| {u_h^j} \right\|} \right|_h^2 + \left( {F\left( {u_h^j} \right),1} \right).\end{equation}$

定理 3.3 令 $ \left\{ {\big( {u_h^j},{w_h^j} \big)} \right\}_{j = 1}^N\ $ 是离散系统(3.1)式的解, 那么这个系统符合以下基本能量定律

$ \begin{equation} E\left( {{t^j};u_h^j} \right) \le E\left( {{t^{j-1}};u_h^{j-1}} \right). \end{equation}$

证 令(3.1)式中的 $ {\psi _h} = \Delta tw_h^j $,(3.1)式可以改写成

$\begin{equation}\bigg( {\frac{{u_h^j - u_h^{j - 1}}}{{\Delta t}},\Delta t\gamma w_h^j} \bigg) + {A}\left( {w_h^j,\Delta t\gamma w_h^j} \right) = 0.\end{equation}$

对于双线性形式 $ {A}\left( { \cdot, \cdot } \right) $ 可以得到

$\begin{equation}{A}\left( {w_h^j, \Delta t\gamma w_h^j} \right) \ge \Delta tC\left| {\left\| {w_h^j} \right\|} \right|_h^2.\end{equation}$

同理令(3.1)式中的 $ {\chi _h} = {u_h^j - u_h^{j - 1}} $, 则(3.1)式可以改写成

$\begin{matrix}{{{\epsilon ^2}}}{A}\left( {u_h^j,\gamma \left( {u_h^j - u_h^{j - 1}} \right)} \right) + \left( {f\left( {u_h^j} \right),\gamma \left( {u_h^j - u_h^{j - 1}} \right)} \right) = \left( {w_h^j,\gamma \left( {u_h^j - u_h^{j - 1}} \right)} \right).\end{matrix}$

为了更好的分析, 这里定义辅助双线性形式 $ {A^s}\left( { \cdot, \cdot } \right) $

$\begin{equation}{A^s}\left( {u,\chi} \right) = \left( {\nabla u,\nabla \chi} \right) - \sum\limits_{e \in \Gamma } {\int_e {\left\{ {\nabla u} \right\} \cdot \left[ \chi \right]} } {{\rm d}s} -\sum\limits_{e \in \Gamma } {\int_e {\left\{ {\nabla \chi} \right\} \cdot \left[ u \right]} } {{\rm d}s + }\sum\limits_{{{\rm e}} \in \Gamma } {\frac{\alpha }{{{h}_{{\rm e}}^\beta }}\int_{{\rm e}} {\left[ u \right]\left[ \chi \right]} } {{\rm d}s}.\end{equation}$

通过定理3.1得到

$\begin{equation}A\left( {u,\gamma \chi} \right) - {A^s}\left( {u,\chi} \right)= - \sum\limits_{e \in \Gamma } {\int_e {\left\{ {\nabla u} \right\} \cdot \left[ {\gamma \chi - \chi} \right]} } {{\rm d}}s -\sum\limits_{e \in \Gamma } {\int_e {\left\{ {\nabla \chi} \right\} \cdot \left[ {\gamma u - u} \right]} } {{\rm d}}s = 0.\end{equation}$

(3.13)式的第一项计算可得

$\begin{equation}A\left( {u_h^j,\gamma \left( {u_h^j - u_h^{j - 1}} \right)} \right) = {A^s}\left( {u_h^j,u_h^j - u_h^{j - 1}} \right)= \Delta t{A^s}\left( {u_h^j,{d_t}u_h^j} \right).\end{equation}$

利用代数恒等式 $ a\left( {a - b} \right) = \frac{1}{2}\left( {{a^2} - {b^2}} \right) + \frac{1}{2}{\left( {a - b} \right)^2} \ge \frac{1}{2}\left( {{a^2} - {b^2}} \right) $ 我们可以得到

$\begin{equation}\Delta tA^s\left( {u_h^j,{d_t}u_h^j} \right) \ge \frac{1}{2}\left| {\left\| {u_h^j} \right\|} \right|_h^2 - \frac{1}{2}\left| {\left\| {u_h^{j - 1}} \right\|} \right|_h^2.\end{equation}$

根据 $ \left| {f\left( u \right)} \right| \le C $, 则有

$\begin{matrix}\left( {f\left( {u_h^j} \right),\gamma \left( {u_h^j - u_h^{j - 1}} \right)} \right)& =&\left( {f\left( {u_h^j} \right),u_h^j - u_h^{j - 1}} \right) \\&&+ \left( {f\left( {u_h^j} \right),\gamma u_h^j - u_h^j - \left( {\gamma u_h^{j - 1} - u_h^{j - 1}} \right)} \right) \\&\ge& \Delta t\left( {f\left( {u_h^j} \right),{d_t}u_h^j} \right).\end{matrix}$

对 $ f ({u_h^j}) $ 进行改写

$\begin{matrix} f\left( {u_h^j} \right) &= &{\left( {u_h^j} \right)^3} - u_h^{j - 1} \\ &= &u_h^j\left( {{{\left| {u_h^j} \right|}^2} - 1} \right) + \Delta t{d_t}u_h^j \\ &= &\frac{1}{2}\left[ {\left( {u_h^j + u_h^{j - 1}} \right) + \Delta t{d_t}u_h^j} \right]\left( {{{\left| {u_h^j} \right|}^2} - 1} \right) + \Delta t{d_t}u_h^j.\end{matrix}$

因此, 方程(3.13)左侧的第二项可写成以下形式

$\begin{matrix} \Delta t\left( {f\left( {u_h^j} \right),{d_t}u_h^j} \right) &\ge& \frac{{\Delta t}}{4}{d_t}{\left\| {\left( {{{\left| {u_h^j} \right|}^2} - 1} \right)} \right\|^2} + \frac{{\Delta {t^2}}}{4}{\left\| {{d_t}\left( {{{\left| {u_h^j} \right|}^2} - 1} \right)} \right\|^2} + \frac{{\Delta {t^2}}}{2}{\left\| {{d_t}u_h^j} \right\|^2} \\ &=& \left( {F\left( {u_h^j} \right) - F\left( {u_h^{j - 1}} \right),1} \right) + \frac{{\Delta {t^2}}}{4}{\left\| {{d_t}\left( {{{\left| {u_h^j} \right|}^2} - 1} \right)} \right\|^2} + \frac{{\Delta {t^2}}}{2}{\left\| {{d_t}u_h^j} \right\|^2} \\ &\ge& \left( {F\left( {u_h^j} \right) - F\left( {u_h^{j - 1}} \right),1} \right).\end{matrix}$

可以直接得到

$\begin{equation}\left( {w_h^j,\gamma \left( {u_h^j - u_h^{j - 1}} \right)} \right) = \left( {\frac{{u_h^j - u_h^{j - 1}}}{{\Delta t}},\Delta t\gamma w_h^j} \right).\end{equation}$

最后, 将(3.1)-(3.21)式整理得到

$\begin{equation}\frac{{{\epsilon ^2}}}{2}\left| {\left\| {u_h^j} \right\|} \right|_h^2 - \frac{{{\epsilon ^2}}}{2}\left| {\left\| {u_h^{j - 1}} \right\|} \right|_h^2 + \left( {F\left( {u_h^j} \right) - F\left( {u_h^{j - 1}} \right),1} \right) \le - \Delta tC\left| {\left\| {w_h^j} \right\|} \right|_h^2 \le 0.\end{equation}$

这样我们就证明得出(3.1)式是无条件能量稳定的.证毕.

4 误差估计

在本节中, 我们将为半离散间断有限体积元方法开发范数 $ {\left\| {\left| \cdot \right|} \right\|_h} $ 中的误差估计. 双线性形式 $ A\left( { \cdot,\gamma \cdot } \right) $ 的强制性和有界性在文献[3,12]中得到证明.

引理 4.1(强制性) 存在一个与 $ h $ 无关的正常数 $ C $, 使得

$\begin{equation}A\left( {{u_h},\gamma {u_h}} \right) \ge C\left\| {\left| {{u_h}} \right|} \right\|_h^2,uad \forall {u_h} \in {V_h}.\end{equation}$

引理 4.2(有界性) 对于 $ \forall {u_h},{v_h} \in {V_h} $ 存在一个与 $ h $ 无关的正常数 $ C $, 使得

$\begin{equation} A\left({{u_h},\gamma {v_h}} \right) \le C{\left\| {\left| {{u_h}} \right|} \right\|_h}{\left\| {\left| {{v_h}} \right|} \right\|_h}. \end{equation}$

引理 4.3假设 $ u\in H_0^1\left( \Omega \right) \cap {H^2}\left( \Omega \right) $ 且 $ {u_h}, v_h \in {V_h} $ 是问题的解,

$ \begin{equation} A\left( {u,\gamma {v_h}} \right) = \left( {g,\gamma {v_h}} \right), uad A\left( {u_h,\gamma {v_h}} \right) = \left( {g,\gamma {v_h}} \right), \end{equation}$

那么存在一个独立于 $ h $ 的正常数 $ C $ 使得

$ \begin{equation} {\left\| {\left| {u - {u_h}} \right|} \right\|_h} \le Ch{\left| u \right|_2}. \end{equation}$

上述引理使用文献[12] 中的方法可以得到. 方程(4.3)的一阶 $ L^2 $ - 误差估计, 将在后面的分析中使用,

$\begin{equation} \left\| {u - {u_h}} \right\| \le Ch{\left| u \right|_2}. \end{equation}$

我们引入一个椭圆投影算子 $ {P_h}:H_0^1\left( \Omega \right) \cap {H^2}\left( \Omega \right) \to {V_h} $ 被定义为

$ \begin{equation} A\left( {Pu,\gamma {v_h}} \right) = A\left( {u,\gamma {v_h}} \right),\forall {v_h} \in {V_h}. \end{equation}$

那么通过引理4.3我们有

$ \begin{equation} \begin{array}{l} {\left\| {\left| {u - Pu} \right|} \right\|_h} \le Ch{\left| u \right|_2},\forall u \in H_0^1\left( \Omega \right) \cap {H^2}\left( \Omega \right), \\ \left\| {u - Pu} \right\| \le Ch{\left| u \right|_2},\forall u \in H_0^1\left( \Omega \right) \cap {H^2}\left( \Omega \right). \\ \end{array} \end{equation}$

注 4.1 上述引理对 $ w \in {P_h} $ 也适用.

定理 4.1 令 $ u $ 和 $ u_h $ 分别是(2.14)式和(2.13)式的解. 若 $ {u_t} \in {L^2}\left( {0,L;{H^2}\left( \Omega \right)} \right) $, 那么存在一个独立于 $ h $ 的正常数 $ C $ 使得

$ \begin{matrix} \int_0^t {{{\left\| {\left| {w - {w_w}} \right|} \right\|}_h}} {\rm d}\tau + {\epsilon ^2}{\left\| {\left| {u - {u_u}} \right|} \right\|_h} \le Ch\left( {\int_0^t {\left( {{{\left| u \right|}_2} + {{\left| {{u_t}} \right|}_2} + {{\left| w \right|}_2}} \right)} {\rm d}\tau + {{\left| u \right|}_2}} \right). \end{matrix}$

证 令 $ u - {u_h} = u - Pu + Pu - {u_h} = {\rho _u} + {\theta _u} $ 和 $ w - {w_h} = w - Pw + Pw - {w_h} = {\rho _w} + {\theta _w} $. 然后从(2.14)式中减去(2.13)式得到

$ \begin{equation} \begin{array}{l} \left( {{{\left( {u - {u_h}} \right)}_t},\gamma {\psi _h}} \right) + A\left( {\left( {w - {w_h}} \right),\gamma {\psi _h}} \right) = 0, \\ {\epsilon ^2}A\left( {u - {u_h},\gamma {\chi _h}} \right) + \left( {f\left( u \right) - f\left( {{u_h}} \right),\gamma {\chi _h}} \right) - \left( {w - {w_h},\gamma {\chi _h}} \right) = 0. \\ \end{array} \end{equation}$

令 $ {\psi _h} = {\rho _w} + {\theta _w},{\chi _h} = {\rho _{ut}} + {\theta _{ut}} $, 将上式相加得到

$ \begin{equation} A\left( {{\rho _w} + {\theta _w},\gamma \left( {{\rho _w} + {\theta _w}} \right)} \right) + {\epsilon ^2}A\left( {{\rho _u} + {\theta _u},\gamma \left( {{\rho _{ut}} + {\theta _{ut}}} \right)} \right) + \left( {f\left( u \right) - f\left( {{u_h}} \right),\gamma \left( {{\rho _{ut}} + {\theta _{ut}}} \right)} \right) = 0. \end{equation}$

通过(4.6)式可将简化为

$ \begin{matrix} A\left( {{\theta _w},\gamma {\theta _w}} \right) + {\epsilon ^2}A\left( {{\theta _u},\gamma {\theta _{ut}}} \right) & =& - \left( {f\left( u \right) - f\left( {Pu} \right),\gamma {\theta _{ut}}} \right) - \left( {f\left( u \right) - f\left( {Pu} \right),\gamma {\rho _{ut}}} \right) \\&& - \left( {f\left( {Pu} \right) - f\left( {{u_h}} \right),\gamma {\theta _{ut}}} \right) - \left( {f\left( {Pu} \right) - f\left( {{u_h}} \right),\gamma {\rho _{ut}}} \right) \\ &=& {{\rm I}_1} + {{\rm I}_2} + {{\rm I}_3} + {{\rm I}_4}. \end{matrix}$

我们通过 Young-不等式得到

$\begin{equation} \begin{array}{l} \left| {{{\rm I}_1}} \right| = \left| {\left( {f'\left( \xi \right){\rho _u},\gamma {\rho _{ut}}} \right)} \right| \le \left\| {f'} \right\|\left| {\left( {{\rho _u},\gamma {\rho _{ut}}} \right)} \right| \le \frac{{{\alpha _1}}}{2}{\left\| {{\rho _{ut}}} \right\|^2} + \frac{{{C_0}}}{{2{\alpha _1}}}{\left\| {{\rho _u}} \right\|^2}, \\[3mm] \left| {{{\rm I}_2}} \right| = \left| {\left( {f'\left( \xi \right){\rho _u},\gamma {\theta _{ut}}} \right)} \right| \le \left\| {f'} \right\|\left| {\left( {{\rho _u},\gamma {\theta _{ut}}} \right)} \right| \le \frac{{{\alpha _2}}}{2}\left\| {{\theta _{ut}}} \right\|_0^2 + \frac{{{C_0}}}{{2{\alpha _2}}}{\left\| {{\rho _u}} \right\|^2}. \end{array} \end{equation}$

对于 $ f(u)=u^3-u $ 的情况, 一个简单的计算给出

$ \begin{equation} f\left( u \right) = f\left( v \right) = {\left( {u - v} \right)^3} + f'\left( u \right)\left( {u - v} \right) - 3u{\left( {u - v} \right)^2}. \end{equation}$

可以将 $ {{\rm I}_3} $ 和 $ {{\rm I}_4} $ 重写为

$ {{\rm I}_3} = {{{\rm I}'}_3} + {{{\rm I}''}_3} + {{{\rm I}'''}_3},uad {{\rm I}_4} = {{{\rm I}'}_4} + {{{\rm I}''}_4} + {{{\rm I}'''}_4}.$

其中

$ \begin{matrix}&& \left| {{{{\rm I}'}_3}} \right| = \left| {\left( {{{\left( {{\theta _u}} \right)}^3},\gamma {\rho _{ut}}} \right)} \right| \le C{\left\| {\left| {{{\left( {{\theta _u}} \right)}^3}} \right|} \right\|_0}{\left\| {\left| {{\rho _{ut}}} \right|} \right\|_0} \le C{\left\| {{\rho _{ut}}} \right\|^2} + C\left\| {{{\left( {{\theta _u}} \right)}^3}} \right\|_0^2, \\ && \left| {{{{\rm I}''}_3}} \right| = \left| {\left( {f'\left( {Pu} \right){\theta _u},\gamma {\rho _{ut}}} \right)} \right| \le \frac{{{\alpha _3}}}{2}\left\| {{\rho _{ut}}} \right\|^2 + \frac{{{C_0}}}{{2{\alpha _3}}}\left\| {{\theta _u}} \right\|_0^2, \\ && \left| {{{{\rm I}'''}_3}} \right| = \left| { - \left( {3\left( {Pu} \right){(\theta _u)^2},\gamma {\rho _{ut}}} \right)} \right| \le \frac{{{\alpha _4}}}{2}\left\| { {\rho _{ut}}} \right\|^2 + \frac{{{C_1}}}{{2{\alpha _4}}}\left\| {{{\left( {{\theta _u}} \right)}^2}} \right\|_0^2. \end{matrix}$

通过(4.13)式同理可得

$\begin{matrix} && {{{\rm I}'}_4} = -\left( {{{\left( {{\theta _u}} \right)}^3},\gamma {\theta_{ut}}} \right) = -\frac{\rm d}{{{\rm d}t}}\left\| {{{\left( {{\theta _u}} \right)}^2}} \right\|_0^2, \\&& \left| {{{{\rm I}''}_4}} \right| = \left| {\left( -{f'\left( {Pu} \right){\theta _u},\gamma {\theta _{ut}}} \right)} \right| \le \frac{{{\alpha _5}}}{2}\left\| {{\theta _{ut}}} \right\|_0^2 + \frac{{{C_0}}}{{2{\alpha _5}}}\left\| {{\theta _u}} \right\|_0^2, \\ && \left| {{{{\rm I}'''}_4}} \right| = \left| { - \left( {3\left( {Pu} \right){{\left( {{\theta _u}} \right)}^2},\gamma {\theta _{ut}}} \right)} \right| \le \frac{{{\alpha _6}}}{2}\left\| {{\theta _{ut}}} \right\|_0^2 + \frac{{{C_1}}}{{2{\alpha _6}}}\left\| {{{\left( {{\theta _u}} \right)}^2}} \right\|_0^2. \end{matrix}$

将(4.12)-(4.15)式整理得到

$ \begin{matrix}&&\left\| {\left| {{\theta _w}} \right|} \right\|_h^2 + \frac{\rm d}{{{\rm d}t}}\left( {{\epsilon ^2}A\left( {{\theta _u},\gamma {\theta _u}} \right) + \left\| {{{\left( {{\theta _u}} \right)}^2}} \right\|_0^2} \right) \\ &\le& \left( {\frac{{{\alpha _1} + {\alpha _3} + {\alpha _4}}}{2} + C} \right){\left\| {{\rho _{ut}}} \right\|^2} + \left( {\frac{{{C_0}}}{{2{\alpha _1}}} + \frac{{{C_0}}}{{2{\alpha _2}}}} \right){\left\| {{\rho _u}} \right\|^2} \frac{{{\alpha _2} + {\alpha _5} + {\alpha _6}}}{2}\left\| {{\theta _{ut}}} \right\|_0^2 \\ &&+ \left( {\frac{{{C_0}}}{{2{\alpha _3}}} + \frac{{{C_0}}}{{2{\alpha _5}}}} \right)\left\| {{\theta _u}} \right\|_0^2 + \left( {\frac{{{C_1}}}{{2{\alpha _4}}} + \frac{{{C_1}}}{{2{\alpha _6}}}} \right)\left\| {{{\left( {{\theta _u}} \right)}^2}} \right\|_0^2 + C\left\| {{{\left( {{\theta _u}} \right)}^3}} \right\|_0^2. \end{matrix}$

根据引理4.3得到

$ \begin{equation} \int_0^t {{{\left\| {\left| {{\theta _w}} \right|} \right\|}_h}} {\rm d}\tau + {\epsilon ^2}{\left\| {\left| {{\theta _u}} \right|} \right\|_h} \le Ch\int_0^t {\left( {{{\left| u \right|}_2} + {{\left| {{u_t}} \right|}_2}} \right)} {\rm d}\tau. \end{equation}$

最后通过三角不等式完成证明.

注 4.2 Cahn-Hilliard 方程全离散间断有限体积元方法的误差估计类似于定理 4.1 的证明.

注 4.3 其中 $ {\alpha _i},i = 1,2, \cdots,6 $ 是一个任意正常数, 为了简便, 到最后都选择为 $ C $.

5 数值算例

本节中, 我们将通过数值算例来验证所得数值格式的有效性. 算例 1, 我们使用一致网格剖分来验证 $ L^2 $ -误差和 $ H^1 $ -误差, 并验证其收敛阶, 结果表明了程序的正确性. 算例2, 结合网格自适应验证方法的有效性. 使用方形流体气泡的演变来测试相场模型的表面张力效应, 并且我们增大时间步长, 数值格式仍保持稳定. 在所有的数值例子中, 我们取惩罚参数 $ \alpha= 1000 $和$ \beta = 1 $.

算例1 我们给出了所提出方法的数值误差和收敛速度. 选择 Cahn-Hilliard 方程的精确解为

$\begin{equation} u = \sin \left( {\pi x} \right)\sin \left( {\pi y} \right)\sin(t). \end{equation}$

定义误差

$\begin{equation} {e_u} = u - {u_h},{e_w} = w - {w_h}. \end{equation}$

在这个算例中, 设计算域为 $ \Omega {{\rm = }}\left[ {0,1} \right] \times \left[ {0,1} \right] $, $ \epsilon = 1 $, $ \Delta t= 4h^2 $和$ T= 1 $. 在表1和表2我们列出了 $ \theta = -1, 0, 1 $ 时的数值误差和收敛阶.

表1   $ \theta = -1,0,1 $时 $ u $ 的数值误差和收敛阶

新窗口打开| 下载CSV


表2   $ \theta = -1,0,1 $时$ w $的数值误差和收敛阶

新窗口打开| 下载CSV


算例 2 在 $ \Omega {{\rm = }}\left[ {0,1} \right] \times \left[ {0,1} \right] $ 的平方域内模拟了方形流体气泡的演化, 选择 $ \epsilon=0.01,$$\Delta t=0.001 $, 初始相位参数如下所示

$\begin{equation}{u_0} = - \tanh \left( {\frac{{\left| {x + y - 1} \right| + \left| {x - y} \right| - 0.5}}{{\sqrt 2 \epsilon }}} \right).\end{equation}$

图2,4,6和8显示了气泡在表面张力作用下的动态演化, 气泡在表面张力作用下变成了一个圆. 图3,5,7和9显示了对应的网格结果表明网格能够准确地捕捉相界面的演变. 同时我们发现方形气泡迅速变形为圆形气泡并且总能量是单调递减的, 并趋向于稳定. 图10展示的是离散相对质量和总能量的演变, 说明数值格式遵循质量守恒性质和离散能量耗散.

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2   $ \Delta t = 0.001 $时不同时刻的初始值图像及数值解


图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   $ \Delta t = 0.001 $时不同时刻对应的自适应网格


图4

新窗口打开| 下载原图ZIP| 生成PPT

图4   $ \Delta t = 0.005 $时不同时刻的初始值图像及数值解


图5

新窗口打开| 下载原图ZIP| 生成PPT

图5   $ \Delta t = 0.005 $时不同时刻对应的自适应网格


图6

新窗口打开| 下载原图ZIP| 生成PPT

图6   $ \Delta t = 0.01 $时不同时刻的初始值图像及数值解


图7

新窗口打开| 下载原图ZIP| 生成PPT

图7   $ \Delta t = 0.01 $时不同时刻对应的自适应网格


图8

新窗口打开| 下载原图ZIP| 生成PPT

图8   $ \Delta t = 0.05 $时不同时刻的初始值图像及数值解


图9

新窗口打开| 下载原图ZIP| 生成PPT

图9   $ \Delta t = 0.05 $时不同时刻对应的自适应网格


图10

新窗口打开| 下载原图ZIP| 生成PPT

图10   左:能量定律;右:质量守恒


6 总结

本文对 Cahn-Hilliard 方程的间断有限体积元方法进行了研究.空间离散利用间断有限体积元方法, 时间离散采用向后欧拉方法, 得到 Cahn-Hilliard 方程全离散格式, 验证了该格式质量守恒性质和离散能量耗散.最后通过数值实验提出了一种自适应时间步进策略计数数值格式的误差和收敛阶, 证明了该格式的有效性.

参考文献

Bi C J, Liu M M.

A discontinuous finite volume element method for second-order elliptic problems

Numerical Methods for Partial Differential Equations, 2012, 28(2): 425-440

DOI:10.1002/num.v28.2      URL     [本文引用: 1]

Cahn J W, Hilliard J E.

Free energy of a nonuniform system I: Interfacial free energy

The Journal of Chemical Physics, 1958, 28(2): 258-267

DOI:10.1063/1.1744102      URL     [本文引用: 1]

It is shown that the free energy of a volume V of an isotropic system of nonuniform composition or density is given by : NV∫V [f0(c)+κ(▿c)2]dV, where NV is the number of molecules per unit volume, ▿c the composition or density gradient, f0 the free energy per molecule of a homogeneous system, and κ a parameter which, in general, may be dependent on c and temperature, but for a regular solution is a constant which can be evaluated. This expression is used to determine the properties of a flat interface between two coexisting phases. In particular, we find that the thickness of the interface increases with increasing temperature and becomes infinite at the critical temperature Tc, and that at a temperature T just below Tc the interfacial free energy σ is proportional to (Tc−T)32.

Chou H S, Ye X.

Unified analysis of finite volume methods for second order elliptic problems

Siam Journal on Numerical Analysis, 2007, 45: 1639-1653

DOI:10.1137/050643994      URL     [本文引用: 1]

Feng X, Karakashian O A.

Fully discrete dynamic mesh discontinuous Galerkin methods for the Cahn-Hilliard equation of phase transition

Mathematics of computation, 2007, 76(259): 1093-1117

DOI:10.1090/S0025-5718-07-01985-0      URL     [本文引用: 1]

Feng X B, Prohl A.

Analysis of a fully discrete finite element method for the phase field model and approximation of its sharp interface limits

Mathematics of Computation, 2004, 73(246):541-567

DOI:10.1090/mcom/2004-73-246      URL     [本文引用: 1]

Jia Z.

Discovering phase field models from image data with the pseudo-spectral physics informed neural networks

Communications on Applied Mathematics and Computation, 2021, 3(2): 357-369

DOI:10.1007/s42967-020-00105-2      [本文引用: 1]

Kumar S, Ruiz-Baier R.

Equal order discontinuous finite volume element methods for the Stokes problem

Journal of Scientific Computing, 2015, 65(3): 956-978

DOI:10.1007/s10915-015-9993-7      URL     [本文引用: 1]

Li R, Gao Y L, Chen J, et al.

Discontinuous finite volume element method for a coupled Navier-Stokes-Cahn-Hilliard phase field model

Advances in Computational Mathematics, 2020, 46(2): 1-35

DOI:10.1007/s10444-020-09758-2      [本文引用: 1]

Liu H L, Yin P M.

Unconditionally energy stable discontinuous Galerkin schemes for the Cahn-Hilliard equation

Journal of Computational and Applied Mathematics, 2021, 390: 113375

DOI:10.1016/j.cam.2020.113375      URL    

Shen J, Yang X F.

Numerical approximations of Allen-Cahn and Cahn-Hilliard equations

Discrete and Continuous Dynamical Systems, 2010, 28: 1669-1691

DOI:10.3934/dcds.2010.28.1669      URL     [本文引用: 1]

Shen J, Yang X F, Yu H J.

Efficient energy stable numerical schemes for a phase field moving contact line model

Journal of Computational Physics, 2015, 284: 617-630

DOI:10.1016/j.jcp.2014.12.046      URL     [本文引用: 1]

Ye X.

A new discontinuous finite volume method for elliptic problems

SIAM Journal on Numerical Analysis, 2004, 42(3): 1062-1072

DOI:10.1137/S0036142902417042      URL     [本文引用: 3]

叶兴德, 程晓良.

Cahn-Hilliard 方程的拟谱逼近

数学物理学报, 2002, 22A(2): 270-280

[本文引用: 2]

Ye X D, Cheng X L.

Legendre collocation approximation for Cahn-Hilliard equation

Acta Mathematica Scientia, 2002, 22A(2): 270-280

[本文引用: 2]

Zhang Z R, Qiao Z H.

An adaptive time-stepping strategy for the Cahn-Hilliard equation

Communications in Computational Physics, 2012, 11(4): 1261-1278

DOI:10.4208/cicp.300810.140411s      URL     [本文引用: 1]

This paper studies the numerical simulations for the Cahn-Hilliard equation which describes a phase separation phenomenon. The numerical simulation of the Cahn-Hilliard model needs very long time to reach the steady state, and therefore large time-stepping methods become useful. The main objective of this work is to construct the unconditionally energy stable finite difference scheme so that the large time steps can be used in the numerical simulations. The equation is discretized by the central difference scheme in space and fully implicit second-order scheme in time. The proposed scheme is proved to be unconditionally energy stable and mass-conservative. An error estimate for the numerical solution is also obtained with second order in both space and time. By using this energy stable scheme, an adaptive time-stepping strategy is proposed, which selects time steps adaptively based on the variation of the free energy against time. The numerical experiments are presented to demonstrate the effectiveness of the adaptive time-stepping approach.

/

〈 〉

版权所有 © 《数学物理学报》杂志社
地址： 武汉市71010号信箱(430071)　电话: 027-87199206(中文版)　027-87199087(英文版)
E-mail: actams@apm.ac.cn
本系统由北京玛格泰克科技发展有限公司设计开发