液晶的研究对于基础科学研究和工业应用都有重要意义. 许多专家和学者先后从微观分子到宏观连续建立了不同层次的数学理论^{[1, 2]}. 目前研究主要集中在用以描述液晶中的缺陷、分子分布规律以及动力学规律的Ericksen-Leslie动力学方程组^[3-6].

简化的Ericksen-Leslie模型既保持了原模型的数学性质, 又降低了研究的难度. 求解该方程最大的困难在于指向矢的约束在数值模拟时难以实现, 为此提出一些相对松弛的约束条件, 如加罚方法和鞍点方法. 加罚方法研究由[7]提出, 之后被广泛研究应用, 可参考[8-12]. 使用加罚方法的不足之处在于构造的数值格式稳定性过于依赖. 在文献[13]中将加罚方法和极限方法置于统一的框架下进行研究和探索, 该方法在克服指向矢的约束的同时引入的拉格朗日乘子还可以克服非线性带来的困难. 值得注意的是鞍点方法的研究较少, 新的实验工作和理论推广需要大量的建模和分析工作, 从分析的角度来看, 现有的工作远远不够. 本文使用增量投影算法结合显隐处理构造了一个一阶精度的、线性的、无条件稳定格式, 以期为相关研究提供参考.

在$ Q=(0, T)\times\Omega , \Omega \in R^N(N=2, 3) $上, 简化的Ericksen-Leslie方程如下:

$ \begin{equation} \left\{ \begin{aligned} \partial_{t}\textbf{d}+(\textbf{u}\cdot\bigtriangledown)\textbf{d} -\gamma(|\nabla d|^2\textbf{d}-\bigtriangleup \textbf{d})=0, \\ |\textbf{d}|=1, \\ \partial_{t}\textbf{u}+(\textbf{u}\cdot\bigtriangledown)\textbf{u}-\nu\triangle \textbf{u}+\nabla p+\lambda\bigtriangledown\cdot\left((\bigtriangledown \textbf{d})^t\bigtriangledown \textbf{d}\right)=0, \\ \bigtriangledown\cdot \textbf{u}=0. \end{aligned} \right. \end{equation} $

(2.1)

其中, $ \textbf{d} $: $ Q\rightarrow {R}^{N} $表示向列液晶流分子的局部指向矢, $ \textbf{u} $: $ Q\rightarrow {R}^{N} $表示流体流动速度, $ p $: $ Q\rightarrow {R}^{N} $表示压力; $ \lambda>0 $是弹性常数; $ \gamma>0 $是松弛时间常数; $ \nu>0 $是流体粘度常数. $ (\nabla \textbf{d})^T $是$ \nabla \textbf{d}=(\partial_j d_i)_{i, j} $的转置, $ |\textbf{d}|=|d(x, t)| $表示$ \textbf{d} $的欧几里得范数, $ (\textbf{u}\cdot \bigtriangledown)\textbf{d}=\sum^M_{i=1}u_i \partial_i \textbf{d} $, $ \triangle \textbf{d}=\sum^M_{i=1}\partial_{ii} \textbf{d} $.

上述方程中梯度的平方增长$ |\nabla d|^2d $以及弹性应力张量$ \bigtriangledown\cdot\left(\left(\bigtriangledown \textbf{d}\right)^T\bigtriangledown \textbf{d}\right) $这两个复杂的非线性项以及非线性代数约束$ |d|=1 $给我们的理论分析及数值模拟带来极大的困难, 特别是调和映照热流方程中的$ |\nabla d|^2d $在高维情况必定会造成奇性.

简化的Ericksen-Leslie方程的鞍点结构如下:

$ \begin{equation} \left\{ \begin{aligned} \partial_{t}\textbf{d}+\left(\textbf{u}\cdot\bigtriangledown\right)\textbf{d}-\gamma\left(q\textbf{d}-\bigtriangleup \textbf{d}\right)=0, \\ \textbf{d}\cdot \textbf{d}-\varepsilon^2q=1, \\ \partial_{t}\textbf{u}+\left(\textbf{u}\cdot\bigtriangledown\right)\textbf{u}-\nu\triangle \textbf{u}+\nabla p+\lambda\bigtriangledown\cdot\left(\left(\bigtriangledown \textbf{d}\right)^T\bigtriangledown \textbf{d}\right)=0, \\ \bigtriangledown\cdot \textbf{u}=0, \end{aligned} \right. \end{equation} $

(2.2)

其中, $ q=\frac{1}{\varepsilon^2}\left(\textbf{d} \cdot \textbf{d}-1\right) $. 对比方程(2.1) 可以发现, 使用鞍点模型引入乘子$ q $, 既克服了$ |\nabla d|^2d $非线性项带来的困难, 又松弛了$ |d|=1 $的限制.

根据(2.2) 第二式可得

$ \begin{equation} q\left(\left(\nabla \textbf{d}\right)^T\nabla \textbf{d}\right)=\frac{1}{2}q\nabla\left(|\textbf{d}|^2\right)=\frac{1}{2}q\nabla\left(\varepsilon^2q\right)=\frac{1}{4\varepsilon^2}\nabla q^2. \end{equation} $

(2.3)

现在对非线性项$ \bigtriangledown\cdot\left(\left(\bigtriangledown \textbf{d}\right)^T\bigtriangledown \textbf{d}\right) $进行如下处理^[13]

$ \begin{equation} \bigtriangledown\cdot\left(\left(\bigtriangledown \textbf{d}\right)^T\bigtriangledown \textbf{d}\right)=\left(\bigtriangledown \textbf{d}\right)^T\triangle d+\frac{1}{2}\nabla\left(|\nabla \textbf{d}|^2\right). \end{equation} $

(2.4)

结合(2.3) 和(2.4), 得

$ \begin{align} \nabla\cdot\left(\left(\nabla \textbf{d}\right)^T\nabla \textbf{d}\right) &=\left(\nabla \textbf{d}\right)^T\triangle \textbf{d}-q\left(\nabla \textbf{d}\right)^T \textbf{d}+\frac{1}{2}\nabla\left(|\nabla \textbf{d}|^2\right)+\frac{1}{4\varepsilon^2}\nabla q^2\\ &=\left(\nabla \textbf{d}\right)^T\left(\triangle \textbf{d}-q\textbf{d}\right)+\frac{1}{2}\nabla\left(|\nabla \textbf{d}|^2+\frac{1}{2\varepsilon^2}q^2\right)\\ &=\left(\nabla \textbf{d}\right)^T\left(\frac{1}{\gamma}\partial_{t}\textbf{d}+\left(\textbf{u}\cdot\bigtriangledown\right)\textbf{d}\right)+\frac{1}{2}\nabla\left(|\nabla \textbf{d}|^2+\frac{1}{2\varepsilon^2}q^2\right). \end{align} $

(2.5)

对(2.2)第二式求导并将压力调整为: $ {\widetilde{p}=p+\frac{\lambda}{2}\left(|\nabla \textbf{d}|^2\right)+\frac{\lambda}{2\varepsilon^2}q^2} $, 为了简便仍将调整后的压力记为$ p $, 得

$ \begin{equation} \left\{ \begin{aligned} \partial_{t}\textbf{d}+\left(\textbf{u}\cdot\bigtriangledown\right)\textbf{d}-\gamma\left(q\textbf{d}-\bigtriangleup \textbf{d}\right)=0, \\ q_{t}=\frac{2}{\varepsilon^2}dd_{t}, \\ \partial_{t}\textbf{u}+\left(\textbf{u}\cdot\bigtriangledown\right)\textbf{u}-\nu\triangle \textbf{u}+\nabla p+\frac{1}{\gamma}\left(\nabla \textbf{d}\right)^T\left(\partial_{t}\textbf{d}+\left(\textbf{u}\cdot\bigtriangledown\right)\textbf{d}\right)=0, \\ \bigtriangledown\cdot \textbf{u}=0. \end{aligned} \right. \end{equation} $

(2.6)

此外, 我们对速度场和方向场分别添加诺伊曼边界条件和齐次狄利克雷边界条件

$ \begin{equation} \partial_n \textbf{d}(x, t)=0 \textbf{u}(x, t)=0, (x, t)\in\partial\Omega , \end{equation} $

(2.7)

给定初始条件为

$ \begin{equation} \textbf{d}(x, t)=\textbf{d}^0(x), \textbf{u}(x, t)=\textbf{u}^0(x, t), (x, t)\in\Omega , \end{equation} $

(2.8)

其中, $ \textbf{u}^0:\Omega\rightarrow \mathcal{R}^M $, $ \textbf{u}^0 \in \textbf{H} $, $ \textbf{d}^0:\Omega\rightarrow \mathcal{R}^M $, $ \textbf{d}^0 \in \textbf{H}^1 $, $ \ |\textbf{d}^0|=1 \in \Omega $.

增量投影方法分为使压力显式化的粘性步和矫正压力的投影步^{[14, 15]}.

粘性步:

$ \begin{equation} \left\{ \begin{aligned} \frac{ \widetilde{u}^{n+1}-u^n}{\delta t}+\left(u^n\cdot\bigtriangledown\right)\widetilde{u}^{n+1}-\nu\Delta \widetilde{u}^{n+1} +\frac{\lambda}{\gamma}\left(\nabla d^n\right)^T\left(\frac{d^{n+1}-d^n}{\delta t}+\widetilde{u}^{n+1}\cdot\bigtriangledown d^n\right) +\nabla p^n=0, \\ \widetilde{u}^{n+1}|_{\partial\Omega}=0. \end{aligned} \right. \end{equation} $

(3.1)

投影步:

$ \begin{equation} \left\{ \begin{aligned} \frac{u^{n+1}-\widetilde{u}^{n+1}}{\delta t}+\nabla\left(p^{n+1}-p^{n}\right)=0, \\ \nabla\cdot u^{n+1}=0, \\ u^{n+1}\cdot n|_{\partial \Omega}=0. \end{aligned} \right. \end{equation} $

(3.2)

把时间区间[0, T] 均匀的分成$ N $个长度为$ \delta t=T/N $的小区间. $ t^n=n\delta t $, $ 0\leq n\leq N $. 给定初始条件$ {d^{0}, u^{0}, q^0=\frac{\left(d^0\right)^2-1}{\varepsilon^2}} $. 结合增量投影算法, 给出(2.6) 的时间离散格式

Step 1: 求解$ d^{n+1} $和中间变量$ \widetilde{u}^{n+1} $

$ \begin{equation} \left\{ \begin{aligned} W+\gamma q^{n+1} d^n-\gamma\triangle d^{n+1}+\frac{\gamma}{\varepsilon^2}(d^{n+1}-d^{n})=0, \\ \begin{split} \frac{ \widetilde{u}^{n+1}-u^n}{\delta t}-\nu\Delta \widetilde{u}^{n+1}+ \frac{\lambda}{\gamma}\left(\nabla d^n\right)^TW+B(u^n, \widetilde{u}^{n+1}) +\nabla p^n=0, \\ \end{split}\\ \partial_n d^{n+1}|_{\partial\Omega}=0, \\ \widetilde{u}^{n+1}|_{\partial\Omega}=0. \end{aligned} \right. \end{equation} $

(3.3)

其中, $ W=\frac{d^{n+1}-d^n}{\delta t}+\widetilde{u}^{n+1}\cdot \bigtriangledown d^n $, $ B\left(u, v\right):=:\left(u\cdot\bigtriangledown\right)v+\frac{1}{2}\left(\nabla\cdot u\right)v $.

Step 2: 求解$ p^{n+1} $

$ \begin{equation} \left\{ \begin{aligned} \bigtriangledown\cdot\tilde{u}^{n+1}+\bigtriangleup\left(p^{n+1}-p^n\right)=0, \\ \partial_n (p^{n+1}-p^n)|_{\partial\Omega}=0. \end{aligned} \right. \end{equation} $

(3.4)

Step 3: 求解$ u^{n+1} $

$ \begin{equation} \left\{ \begin{aligned} \frac{u^{n+1}-\widetilde{u}^{n+1}}{\delta t}+\nabla\left(p^{n+1}-p^{n}\right)=0, \\ \nabla\cdot u^{n+1}=0, \\ u^{n+1}\cdot n|_{\partial \Omega}=0. \end{aligned} \right. \end{equation} $

(3.5)

Step 4: 求解$ q^{n+1} $

$ \begin{equation} \left\{ \begin{aligned} \frac{\varepsilon^2}{2}\frac{q^{n+1}-q^n}{\delta t}=d^n\frac{d^{n+1}-d^n}{\delta t}, \\ \partial_n q^{n+1}|_{\partial\Omega}=0. \end{aligned} \right. \end{equation} $

(3.6)

在方程(3.3) 中增加稳定化项$ \frac{\gamma}{\varepsilon^2}(d^{n+1}-d^{n}) $是为了平衡使用拉格朗日乘子简化非线性项带来的误差. 事实上, $ q^{n+1} $的计算可以从Step 1中解耦出来. 由(3.6) 得知$ q^{n+1}=q^n+\frac{2}{\varepsilon^2}d^{n}\left(d^{n+1}-d^{n}\right) $带入(3.3) 得

$ \begin{equation} \begin{array}{lcl} \frac{d^{n+1}-d^n}{\delta t} +\widetilde{u}^{n+1}\cdot\bigtriangledown d^n+\gamma q^nd^n+\gamma\left(\frac{2}{\varepsilon^2} d^{n}\left(d^{n+1}-d^{n}\right)\right) d^n-\gamma\triangle d^{n+1}\frac{\gamma}{\varepsilon^2}(d^{n+1}-d^{n})=0. \end{array} \nonumber \end{equation} $

另外, 从这四个分步可以看出, 只有在Step 1中$ d^{n+1} $和$ \widetilde{u}^{n+1} $有弱耦合, 如果我们令(3.5) 第一式的第二项为$ \widetilde{u}^{n+1}\cdot\bigtriangledown d^n $, 则我们的半离散格式在求解所有未知数的时候就是全解耦的.

定理   若$ \left(d^{n+1}_{h}, q^{n+1}_{h}, \widetilde{u}^{n+1}_{h}, p^{n+1}_{h}, u^{n+1}_{h}\right) $为方程(3.3-3.6) 的解, 则满足如下能量定律:

$ \begin{eqnarray} \begin{array}{lcl} \left(E\left(u^{n+1}, d^{n+1}, q^{n+1}\right)+\frac{\delta t}{2}\|\nabla p^{n+1}\|^2\right) -\left(E\left(u^n, d^n, q^n\right)+\frac{\delta t}{2}\|\nabla p^n\|\right) \leq0 \end{array} \end{eqnarray} $

(3.7)

其中, $ E\left(u, d, q\right)=\int_{\Omega}\frac{1}{2}|\textbf{u}|^2+\frac{\lambda}{2}|\nabla \textbf{d}|^2+\frac{\lambda\varepsilon^2}{4}q^2 $.

证   将(3.3) 第一式与$ \frac{\lambda}{\gamma}\frac{d^{n+1}-d^n}{\delta t} $做內积得

$ \begin{equation} \begin{array}{lcl} \frac{\lambda}{\gamma}\|W\|^2-\frac{\lambda}{\gamma}\left(W, \widetilde{u}^{n+1}\cdot \bigtriangledown d^n\right)+\left(\lambda q^{n+1}d^n, \frac{d^{n+1}-d^n}{\delta t}\right)+\frac{\lambda}{\varepsilon^2}\|d^{n+1}-d^n\|^2\\ +\frac{\lambda}{2\delta t}\left(\|\nabla d^{n+1}\|^2-\|\nabla d^{n}\|^2 +\|\nabla (d^{n+1}-d^n)\|^2\right)=0. \end{array} \end{equation} $

(3.8)

令(3.3) 第二式与$ \widetilde{u}^{n+1} $做內积, 得

$ \begin{eqnarray} \begin{array}{lcl} \frac{1}{2\delta t}\left(\|\widetilde{u}^{n+1}\|^2-\|u^n\|^2+\|\widetilde{u}^{n+1}-u^n\|^2\right)+\nu\|\nabla \widetilde{u}^{n+1}\|^2+\left(\nabla p^n, \widetilde{u}^{n+1}\right)\\ + \frac{\lambda}{\gamma}\left(\left(\bigtriangledown d^n\right)^TW, \widetilde{u}^{n+1}\right)=0. \end{array} \end{eqnarray} $

(3.9)

令(3.4) 和$ \delta t p^n $做內积可以得到

$ \begin{equation} \frac{\delta t}{2}\left(\|\nabla p^{n+1}\|^2-\|\nabla p^n\|^2-\|\nabla\left(p^{n+1}-p^n\right)\|^2\right)-\left(\nabla p^n, \widetilde{u}^{n+1}\right)=0. \end{equation} $

(3.10)

若(3.5) 第一式分别与$ u^{n+1} $和$ \frac{\delta t}{2} $做內积, 得

$ \begin{eqnarray} \frac{1}{2\delta t}\left(\|u^{n+1}\|^2-\|\widetilde{u}^{n+1}\|^2+\|u^{n+1}-\widetilde{u}^{n+1}\|^2\right)=0. \end{eqnarray} $

(3.11)

$ \begin{equation} \frac{1}{2\delta t}\|u^{n+1}-\widetilde{u}^{n+1}\|^2 =\frac{\delta t}{2}\|\nabla\left(p^{n+1}-p^n\right)\|^2. \end{equation} $

(3.12)

(3.6) 和$ \lambda q^{n+1} $做內积, 得

$ \begin{eqnarray} \frac{\lambda\varepsilon^2}{4\delta t}\left(\|q^{n+1}\|^2-\|q^n\|^2+\|q^{n+1}-q^n\|^2\right) =\lambda\left(d^n\frac{\left(d^{n+1}-d^n\right)}{\delta t}, q^{n+1}\right). \end{eqnarray} $

(3.13)

整合(3.8–3.13), 得

$ \begin{eqnarray} \begin{array}{lcl} \frac{\lambda}{\gamma}\|W\|^2+\frac{\lambda}{2\delta t}\left(\|\nabla d^{n+1}\|^2-\|\nabla d^{n}\|^2+\|\nabla\left(d^{n+1}-d^n\right)\|^2\right)+\frac{\lambda\varepsilon^2}{4\delta t}\left(\|q^{n+1}\|^2-\|q^n\|^2\right.\\ \left.+\|q^{n+1}-q^n\|^2\right)+\frac{1}{2\delta t}\left(\|u^{n+1}\|^2-\|u^n\|^2+\|\widetilde{u}^{n+1}-u^n\|^2\right) +\nu\|\nabla\left(\widetilde{u}^{n+1}\right)\|^2 \\+\frac{\delta t}{2}\left(\|\nabla p^{n+1}\|^2-\|\nabla p^{n}\|^2\right)+\frac{\lambda}{\varepsilon^2}\|d^{n+1}-d^n\|^2=0.\\ \end{array} \end{eqnarray} $

(3.14)

该定理得证.

令$ T_h $是在$ \Omega $上网格大小为$ h $的拟一致剖分, $ X_h $和$ Y_h $分别是$ H^1_{0}(\Omega) $和$ H^1(\Omega) $的有限元近似, 记$ M_{h}=Y_{h}\bigcap L^2_{0}(\Omega):=\{q_{h}\in Y_{h}; \int_{\Omega}q_{h}dx=0\} $, 若$ X_h $和$ M_h $分别是速度场和压力场的稳定近似空间, 则一定存在一个常数$ C $使得

$ \begin{equation} \sup \limits_{u_h \in X_h}\frac{\left(\nabla\cdot u_h, q_h\right)}{\|u_h\|^2}\geq C\|q_h\|^2_{L^2}, \forall q_h\in M_h. \end{equation} $

(3.15)

现给出(3.3–3.6)的全离散格式, 对任意$ (\psi_h, v_{h}, g_{h}, m_h, \varphi_{h}, )\in Y_h\times X_h\times X_{h}\times Y_{h}\times Y_{h}, $求解$ (d^{n+1}_{h}, q^{n+1}_{h}, \widetilde{u}^{n+1}_{h}, p^{n+1}_{h}, u^{n+1}_{h})\in Y_{h}\times Y_{h}\times X_{h}\times M_{h}\times X_{h} $, 使得

Step 1:

$ \begin{equation} \left\{ \begin{aligned} \left(W_h, \psi_{h}\right)+\gamma\left(\nabla d^{n+1}_h, \nabla\psi_{h}\right)+\gamma\left(q^{n+1}_{h}d^n_h, \psi_{h}\right)+\left(\frac{\gamma}{\varepsilon^2}(d^{n+1}_h-d^{n}_h), \psi_{h}\right)=0, \\ \begin{split} \left(\frac{\widetilde{u}^{n+1}_{h}-u^{n}_{h}}{\delta t}, v_{h}\right)+\nu\left(\nabla \widetilde{u}^{n+1}_h, \nabla v_{h}\right) +\left(B\left(u^{n}_{h}, \widetilde{u}^{n+1}_{h}\right), v_{h}\right) +\left(\nabla p^n_h, v_{h}\right)\\ +\frac{\gamma}{\lambda}\left(\left(\nabla d^{n}_{h}\right)^TW_h, v_{h}\right)=0, \end{split} \end{aligned} \right. \end{equation} $

(3.16)

其中, $ W_h=\frac{d^{n+1}_h-d^n_h}{\delta t}+\widetilde{u}^{n+1}_h\cdot \bigtriangledown d^n_h $.

Step 2:

$ \begin{equation} \left(\bigtriangledown\cdot\tilde{u}^{n+1}_h, g_h\right)+\left(\bigtriangleup\left(p^{n+1}_h-p^n_h\right), g_h\right)=0. \end{equation} $

(3.17)

Step 3:

$ \begin{equation} \left\{ \begin{aligned} \left({u^{n+1}_{h}-\tilde{u}^{n+1}_{h}}, g_h\right)+\delta t\left(\nabla\left(p^{n+1}_h-p^n_h\right), g_h\right)=0, \\ \left(\nabla\cdot u^{n+1}_h, m_h\right)=0. \end{aligned} \right. \end{equation} $

(3.18)

Step 4:

$ \begin{equation} \frac{\varepsilon^2}{2}\left(\frac{q^{n+1}_{h}-q^{n}_{h}}{\delta t}, \varphi_h\right) =\left(d^{n}_{h}\frac{d^{n+1}_{h}-d^{n}_{h}}{\delta t}, \varphi_h\right). \end{equation} $

(3.19)

全离散格式保持半离散格式的能量稳定, 不再详细赘述.

本小节中, 借助FreeFem进行数值模拟来验证本文所提离散格式的有效性. 相关程序代码可参考文献[16-18]. 使用$ P_1 $函数空间近似$ Y_h $, 对于$ Y_h\times X_h $取$ P1\times P1b $混合有限元空间近.

首先, 在单位圆$ \Omega=\{(x, y):x^2+y^2<1\} $上计算$ T=0.1 $时$ L^2 $范数意义下的误差和收敛率, 计算公式如下

$ \begin{equation*} E_{u_h}=\|u_{h_i}^t-u^{t}_{h_{i+1}}\|_{L^2}, \quad E_{u^t}=\|u_h^{t_{i}}-u^{t_{i+1}}_{h}\|_{L^2}, \end{equation*} $

$ \begin{equation*} R_{u_h}=\frac{log\left({E(u_{h_i})}/{E(u_{h_{i+1}})}\right)}{log\left({h_i/h_{i+1}}\right)}, \qquad R_{u^t}=\frac{log\left({E(u^{t_{i}})}/{E(u^{t_{i+1}})}\right)}{log\left({t_i/t_{i+1}}\right)}, \end{equation*} $

其中, $ E_{u_h} $代表计算变量$ u $时, 当时间网格固定时由空间网格取值不同造成的空间误差(可简记为空间误差); $ E_{u^t} $代表计算变量$ u $时, 当空间网格固定时由时间网格取值不同造成的时间误差(简记为时间误差). $ R_{u_h} $代表变量$ u $的空间误差收敛速率; $ R^{u_t} $代表变量$ u $的空间误差收敛速率. 时间区间限定为[0, 0.1], 初值为

$ \begin{equation*} u_0=0, \quad q_{0}=\frac{\left(d^0\right)^2-1}{\varepsilon^2}, \quad p_0=\frac{\lambda}{2}|\nabla d_{0}|^2+\frac{\lambda \varepsilon^2}{4}q^2_0, \end{equation*} $

$ \begin{equation*} d_0=\left(sin\left(a\right), cos\left(a\right)\right), \quad a=\pi\left(x^2+y^2\right)^2. \end{equation*} $

其他参数设定为$ \gamma=\lambda=\nu=1 $, $ \varepsilon=0.05 $.

表 1表明本文所构建的数值格式求解所得的$ d, u $和$ p $时间误差在$ L^2 $范数意义下具有一阶精度. 由表 2可以得出该格式得出的$ d, u $的空间误差在$ L^2 $范数意义下的收敛速度是二阶, $ p $的收敛率为1.5.

其次, 在正方形区域$ \Omega=[-1, 1]\times[-1, 1] $上验证所提出格式的能量稳定性以及描述液晶中的缺陷、分子分布规律的准确性.定义Ericksen-Leslie方程的离散能量为:

$ \begin{equation} E=\frac{1}{2}\|\textbf{u}_{h}(t)\|^2+\frac{\lambda}{2}\|\nabla \textbf{d}_{h}(t)\|^2+\frac{\varepsilon^2}{4}q_h(t)^2. \end{equation} $

(4.1)

$ E $代表模型中涉及到的总能量, 其中, $ E_u=\frac{1}{2}\|\textbf{u}_{h}(t)\|^2 $表示动能, $ E_d=\frac{\lambda}{2}\|\nabla \textbf{d}_{h}(t)\|^2 $表示弹性能. 修正后的能量为:

$ \begin{equation} E_{app}=E\left(u^{n+1}, d^{n+1}, q^{n+1}\right)+\frac{\delta t}{2}\|\nabla p^{n+1}\|^2 \end{equation} $

(4.2)

时间为$ T=1 $, 参数取$ \lambda=\gamma=\nu=1 $, $ \varepsilon=0.05, $ $ h=2^{-5}, \delta t=0.001 $. 初始条件取:

$ \begin{equation*} \textbf{u}_0=0, \quad q_0=\frac{\left(d^0\right)^2-1}{\varepsilon^2}, \quad p_0=\frac{\lambda}{2}|\nabla d_0|^2+\frac{\lambda \varepsilon^2}{4}q^2_0, \end{equation*} $

$ \begin{equation*} \textbf{d}_0=\widetilde{\textbf{d}}/\sqrt{|\widetilde{\textbf{d}}|^2+\varepsilon^2}, \widetilde{\textbf{d}}=\left(x^2+y^2-0.25, y\right)^T. \end{equation*} $

图 1详细反映了各能量在湮灭过程中的变化, 其中(a) 验证了所提格式的能量稳定性. 图 2和图 3分别描述了在奇异点湮灭过程中的方向场和速度场的变化, 说明所提格式可以准确描述液晶中的缺陷、分子分布规律. 结合图 1、图 2以及图 3, 可以看出弹性能和总能量在湮灭的时候骤减, 而动能在湮灭时刻达到最大.

本文提出了求解Ericksen-Leslie模型的一阶精度、线性稳定的数值格式, 分析并验证了该格式的稳定性及有效性. 在未来的研究中, 可以结合多重网格或自适应网格进一步提高计算效率、精度.