Ericksen-Leslie模型的一阶能量稳定数值格式

刘 盟，贾宏恩*

(太原理工大学数学科学学院，山西太原 030024)

液晶是物质介于液态和晶态之间的中间状态，是一种中间相.通常按照液晶形成的条件可分为热致型、溶致型和压致型.热致型液晶按照对称性又分为近晶型、向列型和胆甾型，其中向列相液晶流的研究和应用较为广泛.20世纪60年代，液晶的流体动力学理论已经出现[1-2]，特别是Ericksen-Leslie模型广为人知.然而，由于整个系统非常复杂，因此大多数研究都是对模型进行简化和近似后进行的[3].

当流体不可压缩时，得到了简化的Ericksen-Leslie模型[4]，该模型主要由Navier-Stokes方程耦合各向异性弹性应力张量和一个映到球面的对流调和映射热流方程组成.模型中的非线性项给理论分析和数值算例的设计带来了极大的困难，尤其是必须处处满足指向矢量的限制，在数值模拟时是非常难以实现的.为了对指向矢量的约束条件进行松弛，人们提出了各种方法研究液晶流问题，如加罚法[5-8]、鞍点法[9]、谱方法[10-11]等.

针对具有拉伸效应的Ericksen-Leslie加罚模型[6,7,12],文中通过非增量压力校正投影法[13-14],提出了一种无条件能量稳定的解耦的半离散格式，并通过算例验证了其数值精度.

1 模型介绍

本节对具有拉伸效应的Ericksen-Leslie加罚模型进行介绍.

给定有界开区域Ω⊂RM(M=2,3)，其边界∂Ω光滑.定义ΩT=Ω×(0,T],T>0是一个固定的时间点，那么对于任意的(x,t)∈ΩT，有

(1)

(2)

其中，ε>0是加罚参数.函数fε(d)是标量势函数Fε(d)的梯度，即

(3)

对于模型(1)，我们将给出其在速度场上的齐次Dirichlet边界条件和在指向场上的齐次Neumann边界条件：

u(x,t)=0,∂nd(x,t)=0,x∈∂Ω,

(4)

以及初始条件：

u(x,0)=u0(x),d(x,0)=d0(x),x∈Ω.

(5)

为了更好地理解解耦格式，我们引入关于模型(1)的能量耗散定律[7,8,15]：

其中

一般地，文中Lp(Ω)(p≥1)表示定义在Ω上的p阶可积函数空间，且范数定义为

如果p=2,则||·||L2简写为||·||, 记(·,·)为空间L2上的内积，即对u(x),v(x)∈L2,

假设m是非负整数，经典的Sobolev空间

Hm(Ω)={v∈L2:∂kv∈L2,∀|k|≤m},

对应的范数为

2 半离散格式

本节使用非增量压力校正投影法给出模型(1)的一个一阶解耦的半离散格式，并给出该格式无条件能量稳定的证明.

令(d0,u0,p0)是(d0,u0,p0)的一组合适的近似值，且d0=d0,u0=u0,p0=0.我们用等时间步长离散模型(1)，即0=t0

1)求解dn+1.

(6)

其中

同时，HF>0是Fε(d)的Hessie矩阵的L∞范数的界[5],例如

(9)

其中M是空间维数.

2)使用非增量压力校正投影方法，求解un+1.

3)求解pn+1.

-ΔtΔpn+1+·un+1=0.

(11)

因为格式(6)～(11)是线性的，使用类似文献[5-6]的技巧，很容易证明数值解的存在唯一性.下面给出格式(6)～(11)的能量估计.

首先，给出一个在能量耗散证明中需要用到的不等式，其证明详见文献[5-6].

引理1[15]给定d,u,ω∈H1(Ω)且u|∂Ω=0，则下列等式成立：

其中

证明(13)式左边改写为

同理，(13)式右边改写为

综合以上两式，引理得证. 】

定理1格式(6)～(11)的数值解(dn+1,un+1,pn+1)满足能量耗散定律

其中

证明首先，方程(6)的第一式与λΔtωn+1作内积，得

其次，方程(6)的第二式与λ(dn+1-dn)作内积，得

同时，将上式与Δtun+1做内积，则有

此外，(11)式与pn+1做内积，得到

Δt(pn+1,pn+1)-(un+1,pn+1)=0,

即

(15)式两边平方，结合(21)式可得

结合引理1，得到

最后，将(16),(17),(19),(20)和(22)式相加，并利用(23)式，得到

显然(14)式成立.】

3 数值实例

本节首先讨论奇异点的湮灭在不同时刻的状态，随后对解耦格式的收敛阶进行验证.所有的数值算例均使用有限元法进行逼近，指向、速度、压力和辅助变量(d,u,p,ω)分别使用有限元空间(P2,P2,P1,P0)近似，其中P2表示标准分片二次元，P1表示分片线性元，P0表示分片常数元.数值算例均由FreeFem++[16]和Matlab编程实现.

3.1 奇异点的湮灭

例1[5-9]考虑初始条件

选择计算区域Ω=(-1,1)×(-1,1)，时间步长Δt=0.001，对空间使用32×32的标准三角形网格剖分；选取参数ε=0.05,a=1,ν=1,λ=1,γ=1,M=2.图1,2分别展示了指向场和速度场在t=0.001,0.4,0.6,0.8时的快照，可以看到2个奇异点通过形成4个旋涡的速度场被带到原点.

图1 两奇异点湮灭时指向场的演变

3.2 收敛阶

例2考虑初始条件

选择计算区域Ω=(0,1)×(-1/2,1/2)，时间步长Δt=1.25×10-3,6.25×10-4,3.125×10-4,1.5625×10-4,7.8125×10-5,对空间使用64×64的标准三角形网格进行剖分；选取参数ε=0.05,a=1,ν=1,λ=1,γ=1,M=2.图3,图4和表1分别给出了指向、速度和压力在T=0.1时误差e的L2范数、H1范数和收敛阶.可以看出，时间步长越小，指向场数值解的收敛阶越接近1，而速度场和压力场数值解的收敛阶在Δt=3.125×10-4以后便趋于稳定，收敛情况好于指向场.

图2 两奇异点湮灭时速度场的演变

图3 指向d,速度u和压力p误差e的L2范数

图4 指向d,速度u和压力p误差e的H1范数

表1 指向、速度和压力的时间收敛阶

4 结束语

利用非增量压力校正投影法提出了一个新的解耦的数值格式求解带有拉伸效应的向列相液晶流问题的数值解，这是一个无条件能量稳定数值格式.数值实验首先给出了2个奇异点的湮灭过程,t=0.612时，两奇异点湮灭.其次，验证了格式的数值精度，时间步长Δt越小，数值精度越好.