基于双调和修正方法的块共聚物方程能量稳定谱方法

《Computational Materials Science》：An energy-stable spectral method for the block copolymer equation via the biharmonic modified approach

　　•首次将BHM方案应用于非局部嵌段共聚物方程。•线性半隐式方案，每步仅需一次求解，无需辅助变量。•在M1、α和Δt的温和条件下可应用离散能量耗散定理。•在N=2563的大规模3D模拟中，质量守恒精度可达O(10?13)。•在笔记本电脑硬件上开展大规模3D模拟时，CPU处理速度相比GPU可提升40倍。

引言
相场模型是一类用于描述材料形成过程及微观结构演变的重要偏微分方程。这类模型中的基础方程为Cahn–Hilliard方程，该方程被广泛用于研究二元混合物及材料中的相分离现象[1]。针对Cahn–Hilliard方程的数值方法已得到深入研究，其中在[2]、[3]中提出了无条件稳定的数值方案。该方程还衍生出了众多变体，用以描述不同的物理过程和建模场景[4]，其应用范围从图像修复到生物学领域[5]、[6]、[7]。

将Cahn–Hilliard框架扩展到聚合物体系后，便得到了嵌段共聚物方程，这是一种描述由短程界面相互作用与长程非局部相互作用共同驱动的微观相分离现象的基础相场模型[8]、[9]、[10]、[11]。当温度降至临界点以下时，共聚物熔体会出现微观相分离，从而形成包括层片、圆柱体以及螺旋形结构在内的丰富形态结构[12]、[13]、[14]。诸如Cahn–Hilliard方程以及相场晶体方程[15]、[16]、[17]这类相场模型与嵌段共聚物方程在结构上存在相似性，为这些模型开发的数值方法往往也能直接应用于嵌段共聚物体系[18]、[19]。

近期，Yang等人[20]提出了一种针对嵌段共聚物方程的Saul’yev型显式方案。虽然这类基于扫描的方案计算简单，但其本质上是顺序式的，因此不太适合现代GPU架构，也会限制其计算规模扩展能力。正是这些问题促使我们开展了本研究。

我们提出了双调和修正方法（BHM），这是一种半隐式谱方法，最初由Bertozzi、Ju和Lu[18]为具有可变迁移率的四阶扩散方程提出，此处我们将其适配用于嵌段共聚物方程。该方程的形式为：?ut=??2Δ2u+Δ(u3?u)+α(u??u)，其中?>0为界面宽度，α>0用于描述长程聚合物拉伸效应，u?则表示u的空间平均值。BHM方案无需进行能量分解，也不需要辅助变量或扩大系统规模，每个时间步仅需进行一次线性求解。其线性结构使其非常适合用于GPU加速的谱方法。基于我们之前在谱方法以及GPU加速相场计算方面的研究成果[19]、[21]、[22]，我们证明了BHM能够为捕捉二维和三维嵌段共聚物体系的丰富形态动态提供简单高效的框架，且所有相关代码均公开可用，便于他人复现实验结果。

本研究的主要贡献如下：我们提出了适用于嵌段共聚物方程的BHM方案，并建立了离散能量耗散定理，该定理证明了在稳定参数M1和非局部系数α的温和实际条件下，该方案具有能量稳定性。该方案在傅里叶空间中的半隐式结构使其能够通过GPU加速，在消费级硬件上的运算速度相比CPU可提升至40倍。通过在随机初始条件下系统改变平均成分u?，我们在二维和三维环境中进行的数值实验成功再现了包括层片结构、螺旋形结构以及球形微域在内的所有嵌段共聚物形态。

本文的其余部分结构如下：第2节介绍了Ohta–Kawasaki能量泛函、H?1梯度流方程、BHM半隐式离散化方法以及离散能量耗散定理；第3节展示了数值实验结果，包括二维形态验证、GPU性能测试以及不同平均成分条件下的全三维模拟结果；第5节则对研究内容进行了总结，并指出了未来的研究方向。

片段摘录
嵌段共聚物能量泛函
我们采用[8]、[12]中给出的嵌段共聚物能量泛函的Ohta–Kawasaki形式，并遵循[20]中的尺度设定：?E[u]=∫ΩF(u)+?22|?u|2dx+α2∫Ω∫ΩG(x?y)(u(x)?u?)(u(y)?u?)dydx，其中F(u)=14u4?12u2，u?=1|Ω|∫Ωudx，G则是具有周期边界条件的?Δ方程的格林函数。通过引入辅助势ψ使得?Δψ=u?u?，非局部项可简化为α2∫Ω|?ψ|2dx，该项可在每个时间步通过快速傅里叶变换高效求解。

控制方程H?1梯度流方程为：?ut=??2Δ2u+Δ(u3?u)+α(u??u)。

物理验证：Cahn–Hilliard极限情况（α=0）
为了在已有理论体系的基础上验证BHM方案的合理性，我们首先考虑α=0的极限情况，此时该方程可还原为经典的Cahn–Hilliard方程。所有模拟均从随机初始条件U0=0.01rand(N,N)+u?开始，计算区域为Ω=[0,2π]2，网格点数量为N=2562，且采用周期边界条件。需说明的是，这种情况并不符合定理1中的假设条件（即要求α≥L），但即便如此，BHM方案依然能够产生预期的旋涡分解以及域粗化动态，具体结果如图所示。

讨论：迈向二阶时间精度
在本研究中，我们的目标是构建一种简单、能量稳定的一阶BHM方案，用于完整的非局部嵌段共聚物模型，同时证明该方案具备在大型二维和三维模拟中进行GPU并行计算的潜力。我们做出这样的选择是有意为之：BHM方案的单线性求解结构以及其对角线型的傅里叶结构，使其既易于实现，又在GPU硬件上具有极高的计算效率，同时避免了IEQ/SAV类方法所需的辅助变量以及系统规模扩大的问题（详见2.2–2.3节）。

结论
我们提出了用于嵌段共聚物方程的BHM方案，这是一种能量稳定的单求解谱方法，无需进行能量分解，也不需要辅助变量或扩大系统规模。与我们之前针对Cahn–Hilliard方程所提出的BHM方案[19]不同，那次研究并未给出稳定性证明，而本研究则为完整的非局部嵌段共聚物模型建立了严格的离散能量耗散定理。其中Ohta–Kawasaki长程项带来了截然不同的稳定性结构，而非仅仅是Cahn–Hilliard方程的延伸版本。

代码可用性
本文中介绍的所有二维和三维模拟对应的MATLAB源代码，在论文被接收后将会发布在https://github.com/sauloorizaga/BHM-BCP网站上。通过稍作参数调整，即可使用这些代码复现本文中的所有图表，而且该实现方案还可以被扩展应用于其他相场模型。

关于 manuscript 编写过程中生成式AI及AI辅助技术的声明
在撰写本文的过程中，作者使用了Claude（Anthropic）和Gemini工具，以协助语言编辑、LaTeX格式排版以及提升写作效率。在使用这些工具之后，作者对所有内容进行了必要的审核和修改，并对最终发表文章的内容承担全部责任。所有的数学模型、算法、数值方案以及计算结果均完全由作者独立完成。

利益冲突声明
作者声明自己不存在任何可能影响本文研究结果的已知财务利益或个人关系。

致谢
本研究得到了美国国立卫生研究院国家普通医学科学研究所颁发的机构发展奖（IDeA）的支持，资助编号为P20GM103451。作者感谢新墨西哥大学（UNM）提供了用于本研究的MATLAB软件以及并行计算工具箱的许可。同时，作者也感谢匿名审稿人提供的宝贵反馈，这些反馈有助于提升本文的质量。

Saulo Orizaga