求解极坐标系下Allen-Cahn方程的有限差分方法

A Finite Difference Method for the Allen-Cahn Equation in Polar Coordinate System

摘要: 非线性Allen-Cahn方程是材料学中相场模拟模型的一类重要方程，用于描述二元合金在一定温度下进行相位分离的过程，在许多科学领域中有广泛的应用。以往的工作主要在矩形区域上考虑求解，本文研究有效的数值方法在圆形区域上求解该方程。首先将Allen-Cahn方程中拉普拉斯算子写为极坐标的形式，然后采用中心差分方法在空间r方向和θ方向分别进行空间离散，得到非线性常微分方程组，并将网格上的数值解以矩阵形式表示。在时间离散过程中，采用积分因子法结合Krylov子空间的方法进行求解。最后给出数值试验。

Abstract: The nonlinear Allen-Cahn equation is an important phase field simulation model in materials science. It is used to describe the phase separation process of binary alloys at a certain temperature and has been widely used in many scientific fields. The previous work mainly considers the solution on the rectangular region. In this paper, the effective numerical method is used to solve the equation on a circular region. First, the Laplacian operator in Allen-Cahn equation is written as the form of polar coordinates. Then, the central difference method is used to carry out spatial discretization in the space direction r and θ respectively, and the numerical solution on the grid is expressed in matrix form by using large sparse linear system. In the implementation of time discretization, the integral factor method combined with Krylov subspace method was used to solve the problem. Finally, numerical experiments are given.

1. 引言

Allen-Cahn方程最初是作为研究结晶固体反相位边界运动的模型而被提出的 [1]，该方程及其变形在图像处理 [2]、平均曲率运动 [3]、晶体的生长 [4]、和材料科学等实际问题中都发挥着极为重要的作用。

由于此类相场模型没有真解，因此选取有效的数值方法进行模拟尤为重要。Allen-Cahn方程的数值逼近方法有多种，包括半解析傅里叶谱法 [5]，谱方法 [6]，间断有限元方法 [7]，算子分裂方法 [8] 以及求解高维Allen-Cahn方程的紧致差分交替方向隐式方法 [9] 等。近年来，也有学者采用自适应有限元方法 [10]，无网格方法 [11] 对Allen-Cahn方程进行数值求解，上述工作均在直角坐标系下完成。但是，在许多实际问题中，需要考虑圆形区域，这方面的研究工作还不是很多 [12]。本文考虑在圆盘区域上求解Allen-Cahn方程，首先将拉普拉斯算子表示为极坐标中的形式，沿r方向和 $\theta$ 方向进行网格剖分，然后采用中心差分方法在网格点上对Allen-Cahn方程进行空间离散，并将网格点上的数值解以矩阵形式表示。接下来发展积分因子方法进行时间离散，在时间离散过程中，利用梯形求积公式进行近似，使其精度为二阶精度。在实现过程中求解指数矩阵与向量的乘积时，采用Krylov子空间的方法来近似，而不单独求解指数矩阵。

2. 数值方法

考虑计算区域为 $\Omega =\left\{\left(x,y\right):x^2+y^2<\text{1}\right\}$ 的二维非线性Allen-Cahn方程：

$\frac{\text{d}u}{\text{d}t}=\gamma \Delta u+f\left(u\right),\left(x,y,t\right)\in \Omega \times \left(0,T\right],$ (1)

其中， $\gamma$ 为扩散系数； $\Delta u=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}$ 为拉普拉斯算子； $f\left(u\right)=\frac{1}{{\epsilon }^2}\left(u^3-u\right)$ ； $\epsilon$ 是一个给定的正参数且 $\epsilon \ll 1$，

设该方程在r方向边界条件为齐次Neumann边界，在 $\theta$ 方向的边界条件为周期边界：

$\frac{\partial u\left(\text{0},\theta \right)}{\partial r}=\frac{\partial u\left(\text{1},\theta \right)}{\partial r}=0;u\left(r,0\right)=u\left(r,2\pi \right).$ (2)

下面将方程(1)写成极坐标系下的形式：

$\frac{\text{d}u}{\text{d}t}=\gamma \left(\frac{{\partial }^{2}u}{\partial r^2}+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}+\frac{1}{r}\frac{\partial u}{\partial r}\right)+f\left(u\right),\left(r,\theta \right)\in \Omega .$ (3)

将计算区域 $\left[\text{0},\text{1}\right]\times \left[0,2\pi \right]$ 分别沿r和 $\theta$ 方向划分成 $N_r\times N_{\theta }$ 个网格，网格节点表示为

$\left(r_i,{\theta }_j\right)=\left(\left(i-\frac{1}{2}\right)h_r,j{h}_{\theta }\right)$ $i=1,2,\cdots ,N_r;j=0,1,2,\cdots ,N_{\theta }$。在r方向的网格步长为 $h_r=\text{1}/N_r$，在 $\theta$ 方向的网

格步长为 $h_{\theta }=2\pi /N_{\theta }$。时间方向剖分为 $t_n=n\Delta t,n=1,2,\cdots$。其中 $\Delta t$ 为时间步长。设u在网格节点 $\left(r_i,{\theta }_j\right)$ 处的数值解为 $u_{i,j}$，结合泰勒展开式，在在网格点处对(3)式中的项进行离散，略去余项后，得到(3)式的二阶中心差分格式：

$\frac{\text{d}{u}_{i,j}}{\text{d}t}=\gamma \left(\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_r{}^2}+\frac{u_{i+1,j}-u_{i-1,j}}{2r_i{h}_r}+\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{r_i{}^2{h}_{\theta }{}^2}\right)+f\left(u_{i,j}\right).$ (4)

下面对边界上的差分格式进行修正使其达到二阶局部截断误差。考虑到方程(1)在r方向具有齐次

Neumann边界条件，当 $i=1$ 时 $r_1=\frac{\Delta r}{2}$，且 ${\left(\frac{\partial u}{\partial r}\right)}_{1,j}=0$，当 $i=N_r$ 时， ${\left(\frac{\partial u}{\partial r}\right)}_{N_r,j}=0$，将(4)式在r方向边

界的差分格式写为

$\{\begin{array}{l}\frac{\text{d}{u}_{1,j}}{\text{d}t}=\gamma \left(\frac{u_{2,j}-2u_{1,j}}{h_r{}^2}+\frac{u_{2,j}}{2r_1{h}_r}+\frac{u_{1,j+1}-2u_{1,j}+u_{1,j-1}}{r_1{}^2{h}_{\theta }{}^2}\right)+f_{1,j},i=1,\\ \frac{\text{d}{u}_{N_r,j}}{\text{d}t}=\gamma \left(\frac{u_{N_r-1,j}-u_{N_r,j}}{h_r{}^2}+\frac{u_{N_r,j}-u_{N_r-1,j}}{2r_{N_r}{h}_r}+\frac{u_{N_r,j+1}-2u_{N_r,j}+u_{N_r,j-1}}{r_{N_r}{}^2{h}_{\theta }{}^2}\right)+f_{N_r,j},i=N_r.\end{array}$ (5)

考虑到u在 $\theta$ 方向是以 $2\pi$ 为周期的，因此边界值为 $u_{i,0}=u_{i,N_{\theta }},u_{i,1}=u_{i,N_{\theta }+1}$，进而将(4)式在 $\theta$ 方向边界的差分格式写为

$\{\begin{array}{l}\frac{\text{d}{u}_{i,1}}{\text{d}t}=\gamma \left(\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_r{}^2}+\frac{u_{i+1,j}-u_{i-1,j}}{2r_i{h}_r}+\frac{u_{i,N_{\theta }}-2u_{i,1}+u_{i,2}}{r_i{}^2{h}_{\theta }{}^2}\right)+f_{i,1},j=1,\\ \frac{\text{d}{u}_{i,N_{\theta }}}{\text{d}t}=\gamma \left(\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_r{}^2}+\frac{u_{i+1,j}-u_{i-1,j}}{2r_i{h}_r}+\frac{u_{i,N_{\theta }-1}-2u_{i,N_{\theta }}+u_{i,1}}{r_i{}^2{h}_{\theta }{}^2}\right)+f_{i,N_{\theta }},j=N_{\theta }.\end{array}$ (6)

定义离散解 $u_{i,j},1\le i\le N_r,1\le j\le N_{\theta }$ 为 $N_r\times N_{\theta }$ 阶的矩阵U，将非线性项记为F。将矩阵U和矩阵F按列重新排列成向量的形式，记 $\mathcal{U}=\text{vec}\left(U\right)$， $\mathcal{F}\left(\mathcal{U}\right)=\text{vec}\left(F\right)$，结合(5)式与(6)式，定义矩阵A为

$A=\left(\begin{array}{ccccc}T-2D& D& & & D\\ D& T-2D& D& & \\ & \ddots & \ddots & \ddots & \\ & & D& T-2D& D\\ D& & & D& T-2D\end{array}\right),$ (7)

其中矩阵 $D=\text{diag}\left({\beta }_1,{\beta }_2,\cdots ,{\beta }_{N_r}\right),{\beta }_i=\frac{1}{{\left(i-1/2\right)}^2{h}_{\theta }^2{h}_r^2},1\le i\le N_r$。

矩阵 $T=\frac{1}{h_r^2}\left(\begin{array}{ccccc}-2& 1+{\lambda }_1& & & \\ 1-{\lambda }_2& -2& 1+{\lambda }_2& & \\ & \ddots & \ddots & \ddots & \\ & & 1-{\lambda }_{N_r-1}& -2& 1+{\lambda }_{N_r-1}\\ & & & 1-{\lambda }_{N_r}& -1+{\lambda }_{N_r}\end{array}\right)$， ${\lambda }_i=\frac{1}{2\left(i-1/2\right)},1\le i\le N_r$。

由上述矩阵的定义，将(1)式经过空间离散后，得到一组非线性常微分方程组：

$\frac{\text{d}u}{\text{d}t}=\gamma A\mathcal{U}+\mathcal{F}\left(\mathcal{U}\right),\left(r,\theta \right)\in \Omega .$ (8)

下面采用积分因子方法求解方程(8)。两边同乘积分因子 ${\text{e}}^{-At}$，然后在一个时间步长内做积分得到：

${\mathcal{U}}_{n+1}={\text{e}}^{A\Delta t}{\mathcal{U}}_n+{\text{e}}^{A\Delta t}{\displaystyle {\int }_0^{\Delta t}{\text{e}}^{-A\tau }\mathcal{F}\left(\mathcal{U}\right)\text{d}\tau }.$ (9)

对(9)式应用梯形积分公式近似得到：

${\mathcal{U}}_{n+1}={\text{e}}^{A\Delta t}\left({\mathcal{U}}_n+\frac{\Delta t}{2}\mathcal{F}\left({\mathcal{U}}_n\right)\right)+\frac{\Delta t}{2}\mathcal{F}\left({\mathcal{U}}_{n+1}\right).$ (10)

令上式中的 ${\mathcal{U}}_n+\frac{\Delta t}{2}\mathcal{F}\left({\mathcal{U}}_n\right)=v$，则(11)式可写为

${\mathcal{U}}_{n+1}={\text{e}}^{A\Delta t}v+\frac{\Delta t}{2}\mathcal{F}\left({\mathcal{U}}_{n+1}\right).$ (11)

虽然A为大型稀疏矩阵，而 ${\text{e}}^{A\Delta t}$ 为稠密矩阵，因此直接求解指数矩阵困难较大，本文利用Krylov子空间 [13] $K_m\left(\Delta t,v\right)=span\left\{v,\left(A\Delta t\right)v,\cdots ,{\left(A\Delta t\right)}^{m-1}v\right\}$ 的元素来近似指数矩阵与向量的乘积，即 ${\text{e}}^{A\Delta t}v$，而不是单独求解指数矩阵。采用不动点迭代法在每个单元上求解非线性方程组(11)，其中迭代初值选取为 ${\mathcal{U}}_{n+1}^0={\mathcal{U}}_n$。

3. 数值实验

应用本文提出的有限差分法和积分因子法求解下面的数值算例。考虑二维计算区域 $\Omega =\left\{\left(x,y\right):x^2+y^2\le \text{1}\text{.5}\right\}$，针对具有 $\epsilon \ll 1$ 方向具有齐次Neumann边界， $\theta$ 方向具有周期边界条件的方程(1)进行求解。初值选取形状像哑铃的函数，见图1(a)。

$u_0\left(x,y\right)=\{\begin{array}{l}\tanh \left(\frac{3}{\epsilon }\left({\left(x-0.5\right)}^2+y^2-{\left(0.39\right)}^2\right)\right),x>0.14,\\ \tanh \left(\frac{3}{\epsilon }\left(y^2-{\left(0.15\right)}^2\right)\right),-0.3\le x\le 0.14,\\ \tanh \left(\frac{3}{\epsilon }\left({\left(x+0.5\right)}^2+y^2-{\left(0.25\right)}^2\right)\right),x<-0.3,\end{array}$

选取参数 $\epsilon =0.05$，并在每个空间方向上取128个等分点，时间步长 $\Delta t=5\times {10}^{-5}$。选取 $t_1=3.53\times {10}^{-2},$ $t_2=1.18\times {10}^{-1},t_3=2\times {10}^{-1}$ 这3个时间点，其数值结果如图1所示。

(a) (b) (c) (d)

Figure 1. Numerical solution of the case $t_1=3.53\times {10}^{-2},t_2=1.18\times {10}^{-1},t_3=2\times {10}^{-1}$

图1. 算例在 $t_1=3.53\times {10}^{-2},t_2=1.18\times {10}^{-1},t_3=2\times {10}^{-1}$ 时的数值解

上图显示了在不同时间点上Allen-Cahn方程解的数值结果。结果表明，随着时间增大，其形状由初始时间 $t_0$ 的哑铃状不断向圆形聚拢，最终消失。数值结果与文献 [7] 中的结果吻合。

基金项目

辽宁省自然科学基金(20180550996)资助的课题。