ABSTRACT

In this article, basic governing equations of the liquid crystalline micro flow under the electric field were established according to Lesile-Ericksen theory, and the programming and simulation process was assisted by MATLAB software. In order to improve the efficiency of the simulation, the optimization was done according to the characteristics of the program. Since the time and space scale was small, and the number of cycles was large, the key subprogram was optimized with serial computation. In addition, allocating memory for all variables in advance and optimization single program sentence were used to reduce the operating time. The results of optimal program are in accord with the primary programs. In the simulation process of one second, the computing time was reduced from about 912.83 hours to about 47.37 hours, and the calculating efficiency was improved nearly 20 times.

Keywords:Backflow Effect, MATLAB Software, Micro Flow, Program Optimization, Serial Computation

电场作用下液晶微流动数值计算程序 编写与优化

林绿叶，刘春波，王安标

河南工业大学机电学院，河南 郑州

收稿日期：2017年10月5日；录用日期：2017年10月18日；发布日期：2017年10月24日

摘 要

本文是根据液晶引流效应，以小分子液晶Leslie-Ericksen理论为基础，建立了模拟电场作用下产生的液晶微流动基本控制方程组，并运用MATLAB软件辅助编程计算。为提高求解效率，针对程序时间、空间尺度小，循环次数较多的特点，重点对关键子程序段进行了优化，得到了串行计算(for循环)的最佳计算方法。此外，通过预先对所有的变量分配内存和优化单行程序的方法，有效减少了程序运行时间。执行结果显示：优化后的程序运行结果与原程序结果吻合，模拟计算1秒时长时，计算时间由原来约912.83个小时减少到约47.37个小时，计算效率提高近20倍。

关键词 :液晶引流效应，MATLAB软件，微流动，程序优化，串行计算

Copyright © 2017 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

微流体驱动和控制技术作为微机械电子系统(Microelectronics mechanical systems, MEMS)的一个重要分支，是MEMS发展需要解决的关键技术之一。随着科学技术的发展，各种新型微流体驱动技术不断产生，电渗驱动、压电驱动等都已得到实际应用 [1] 。但是每种微流体驱动技术都伴随着难以克服的缺点，传统技术微型化难，电渗驱动对驱动材料敏感，压电驱动驱动电压过高，驱动力小等 [2] [3] 。因此，开发全新的微流体驱动技术仍然是这一领域亟待解决的问题。

近年来，出现了一种利用特殊流体或流体的特殊性质实现微流体驱动与控制的驱动方式，这种驱动方式对所驱动的粒子没有限制，响应快，驱动器的设计加工简单，易于实现微型化。其中液晶引流驱动得到了广泛关注 [4] [5] [6] ，驱动原理主要是利用非牛顿流体液晶同时具有各向异性与流动性的特点，在电磁或温度场的作用下，液晶的指向矢会发生转动，进而引起微流动。这种现象被称为液晶的引流效应，且这种微流动发生于分子尺度，响应快，精度高，且可实现低电压下的精确驱动，对驱动材料无要求。此外，液晶具有流动性，采用这种驱动方式可以简化微驱动器的结构设计，易于制造 [7] 。

液晶作为一种典型的粘弹性非牛顿流体，其特殊的力学性质导致流动状态与普通流体差别很大，加上微尺度的原因，液晶引流驱动相关的数值模拟较之普通流体复杂很多，针对普通流体的商业软件无法进行液晶微流动的模拟计算，需要针对不同的流动状态人工编程计算。目前关于粘弹性流体微流动的数值计算研究较少，使用较多的数值研究方法主要有玻尔兹曼方程、N-S运动方程和分子动力学模型等 [8] 。如栗雪娟等用WCCBS_SU方法模拟微可压缩粘弹流体 [9] ，郭春海等用差分格式和乘方格式的有限元模拟微混合器的流场 [8] ，付一志、焦群英运用龙格–库塔法求解薄壁细胞收微吸管和探针共同作用时内部流场等 [10] 。但是针对液晶的微流动计算的方法研究还没有。与普通非牛顿流体相比，液晶微流动计算的难度主要表现在粘弹性系数较多，本构方程复杂，导致最终的方程组推导难度较大。其次，数值计算时间与空间步长很小，计算量很大，程序执行时间特别长，因此，针对程序特点的优化过程必不可少。

本文将针对液晶微流动的特点，以电场作用下引起液晶微流动为例进行数值模拟。为了提高程序的编写及执行速度，将模拟程序编制过程从纯人工改为基于MATLAB软件进行编制，并依托MATLAB软件进行程序优化，提高运算速率。

2. 数学模型

基于液晶引流效应的微流动数值计算模型的构建基础为小分子液晶的Lesile_Ericksen理论(简称L-E理论) [11] 。L-E理论为计算小分子液晶流动的权威理论，适用于各种粘弹性流体。选用的液晶为棒状小分子液晶5CB，基本的控制方程组包括连续性方程、纳维–斯托克斯方程(N-S方程)、电场引起的外力方程、本构方程、亥姆霍兹自由能密度方程以及液晶微流动计算特有的角运动方程 [12] 。

2.1. 基础控制方程组

连续性方程：

$\nabla \cdot v=0$ (1)

运动方程：

$\rho \frac{Dv}{Dt}=F+\nabla \cdot \left(-pI+\tau \right)$ (2)

本构方程：

$\tau ={\alpha }_{1}nnn\cdot A\cdot n+{\alpha }_{2}nN+{\alpha }_{3}Nn+{\alpha }_{4}A+{\alpha }_{5}nn\cdot A+{\alpha }_{6}A\cdot nn-\frac{\partial F}{\partial \nabla n}\cdot {\left(\nabla n\right)}^{\text{T}}$ (3)

在液晶中施加强度为E的电场，所引起的电场力F大小为：

$F=\left[\left({\epsilon }_{\perp }E+{\epsilon }_a\left(n\cdot E\right)n\right)\cdot \nabla \right]E$ (4)

自由能密度方程：

$2F=K_1{\left(\nabla \cdot n\right)}^2+K_2{\left(n\cdot \nabla \times n\right)}^2+K_3{\left|n\times \nabla \times n\right|}^2$ (5)

角运动方程：

$n\times \left\{G+\frac{\partial F}{\partial n}-\nabla \cdot \left(\frac{\partial F}{\partial \nabla n}\right)+{\lambda }_{1}N+{\lambda }_{2}A\cdot n\right\}=0$ (6)

2.2. 方程组化简展开与初始边界条件

本计算针对的物理模型是常见的液晶盒，两ITO玻璃板(y方向)间隔为110 μm，与另外两个方向的尺寸相比可以忽略，因此计算中不用考虑厚度方向的速度，直接按简单的三维模型，如图1所示的坐标系进行数值计算。其中，单位向量n代表液晶的指向矢。

图1. 计算坐标系

对于液晶盒内微流动的计算，主动力为电场力，阻力主要为粘弹性力，重力影响可忽略，故动量方程可化简为：

y轴方向剪切力为主动力，可认为x轴方向加速度仅由y方向切应力所致。故n_x、n_y、n_z对x、y的偏导数为0，即：

$\nabla \cdot {\tau }_1=\frac{\partial {\tau }_{1,1}}{\partial x}+\frac{\partial {\tau }_{1,2}}{\partial y}+\frac{\partial {\tau }_{1,3}}{\partial z}=\frac{\partial {\tau }_{1,2}}{\partial y}$

将其代入运动方程并展开整理得：

$\begin{array}{c}\rho \frac{\partial u}{\partial t}={\alpha }_1\left\{\frac{\partial }{\partial y}\left(n_x^2{n}_y^2\right)\frac{\partial u}{\partial y}+\frac{\partial }{\partial y}\left(n_x{n}_y^2{n}_z\right)\frac{\partial w}{\partial y}+n_z{n}_y^2\left(n_x\frac{{\partial }^{2}u}{\partial y^2}+n_z\frac{{\partial }^{2}w}{\partial y^2}\right)\right\}+\frac{{\alpha }_4}{2}\frac{{\partial }^{2}w}{\partial y^2}\\ +{\alpha }_2\frac{\partial }{\partial y}\left(n_y\frac{\partial n_x}{\partial t}\right)+{\alpha }_3\frac{\partial }{\partial y}\left(n_x\frac{\partial n_y}{\partial t}\right)-\frac{{\alpha }_2-{\alpha }_5}{2}\left(2n_y\frac{\partial n_y}{\partial y}\frac{\partial u}{\partial y}+n_y^2\frac{{\partial }^{2}u}{\partial y^2}\right)\\ +\frac{{\alpha }_3+{\alpha }_6}{2}\left\{\frac{\partial }{\partial y}\left(n_x{n}_z\right)\frac{\partial w}{\partial y}+2n_x\frac{\partial n_x}{\partial y}\frac{\partial u}{\partial y}+n_x{n}_z\frac{{\partial }^{2}w}{\partial y^2}+n_x^2\frac{{\partial }^{2}u}{\partial y^2}\right\}\end{array}$ (7a)

$\begin{array}{c}\rho \frac{\partial w}{\partial t}={\alpha }_1\left\{\frac{\partial }{\partial y}\left(n_x^2{n}_y^2\right)\frac{\partial w}{\partial y}+\frac{\partial }{\partial y}\left(n_x{n}_y^2{n}_z\right)\frac{\partial u}{\partial y}+n_z{n}_y^2\left(n_x\frac{{\partial }^{2}u}{\partial y^2}+n_z\frac{{\partial }^{2}w}{\partial y^2}\right)\right\}+\frac{{\alpha }_4}{2}\frac{{\partial }^{2}w}{\partial y^2}\\ +{\alpha }_2\frac{\partial }{\partial y}\left(n_y\frac{\partial n_z}{\partial t}\right)+{\alpha }_3\frac{\partial }{\partial y}\left(n_x\frac{\partial n_y}{\partial t}\right)-\frac{{\alpha }_2-{\alpha }_5}{2}\left(2n_y\frac{\partial n_y}{\partial y}\frac{\partial w}{\partial y}+n_y^2\frac{{\partial }^{2}w}{\partial y^2}\right)\\ +\frac{{\alpha }_3+{\alpha }_6}{2}\left\{\frac{\partial }{\partial y}\left(n_x{n}_z\right)\frac{\partial u}{\partial y}+2n_z\frac{\partial n_x}{\partial y}\frac{\partial w}{\partial y}+n_x{n}_z\frac{{\partial }^{2}u}{\partial y^2}+n_z^2\frac{{\partial }^{2}w}{\partial y^2}\right\}\end{array}$ (7b)

同理化简展开角运动方程可得：

${\gamma }_1\frac{\partial n_x}{\partial t}=n_y\left\{2\left(k_3-k_2\right)\frac{\partial n_x}{\partial y}\frac{\partial n_y}{\partial y}-{\alpha }_2\frac{\partial u}{\partial y}\right\}+\frac{{\partial }^2{n}_x}{\partial y^2}\left\{k_2+\left(k_3-k_2\right)n_y^2\right\}$ (7c)

$\begin{array}{c}{\gamma }_1\frac{\partial n_y}{\partial t}=n_y\Delta \epsilon E^2-{\alpha }_3\left(n_x\frac{\partial u}{\partial y}+n_y\frac{\partial w}{\partial y}\right)-\left(k_3-k_2\right)n_y\left\{{\left(\frac{\partial n_x}{\partial y}\right)}^2+{\left(\frac{\partial n_y}{\partial y}\right)}^2-{\left(\frac{\partial n_z}{\partial y}\right)}^2\right\}\\ +\frac{{\partial }^2{n}_y}{\partial y^2}\left\{k_1+\left(k_3-k_2\right)n_y^2\right\}\end{array}$ (7d)

${\gamma }_1\frac{\partial n_z}{\partial t}=n_y\left\{2\left(k_3-k_2\right)\frac{\partial n_z}{\partial y}\frac{\partial n_y}{\partial y}-{\alpha }_2\frac{\partial w}{\partial y}\right\}+\frac{{\partial }^2{n}_z}{\partial y^2}\left\{k_2+\left(k_3-k_2\right)n_y^2\right\}$ (7e)

在此坐标系下，沿着y轴方向施加电场，则微流动微流动速度v、分子指向矢n与电场E可分别表示为：

$v=\left(u,0,w\right)$

$n=\left(\cos \varphi \sin \theta ,\sin \varphi ,\cos \varphi \cos \theta \right)$

$E=\left(0,E,0\right)$

依据L-E理论，液晶分子指向矢n为单位矢量，因此有：

$n_x^2+n_y^2+n_z^2=1$ (8)

初值设定为：液晶盒两ITO板间距离为110 μm，两平板固定，初始状态静止。整个盒内液晶分子的配向场初始值为扭转角 $\varphi =0^{\circ }$ ，倾斜角 $\theta =5^{\circ }$ ，初始速度为0。沿y轴方向施加的连续方波电场参数为幅值为5 V，频率为1 Hz，占空比为20%。

3. 数值计算

3.1. 数值计算方法及参数

对方程组进行离散化处理后进行数值计算，时间方向上采用二阶龙格库塔法，计算间隔为10⁻⁹；空间上采用中心差分法，空间网格大小为 $\Delta y=H/100$ ，H为液晶盒厚度。若采用单步校正的修恩法 [13] ，编程计算的算子为

$y_{i+1}=y_i+\frac{1}{2}\left(k_1+k_2\right)h$ (9)

式中， $y_i,y_{i+1}$ 为某物理量(如液晶微流动的速度，指向矢)旧值、新值， $h$ 为前进步长， $k_1$ 为按 $\left(x_i,y_i\right)$ 计算出的斜率， $k_2$ 为按 $\left(x_{i+1},y_{i+1}\right)$ 计算出的斜率。

3.2. MATLAB程序编写

根据上述初始条件和数学模型，输入Leslie粘性系数( ${\alpha }_1~{\alpha }_6$ )、弹性系数( $k_1~k_3$ )、液晶密度 $\rho$ 、电场强度E等常数及液晶初始速度v和指向矢n。其中，指向矢的数值由液晶分子倾斜角 $\theta$ 和扭转角 $\varphi$ 按线性分布计算得到，如图1。将液晶盒沿 $y$ 轴方向100等分，按线性分布计算这101个节点处的倾斜角和扭转角大小，进而计算出指向矢n大小。

按照编程计算算子用Matlab软件进行编程，步骤如下：

1) 计算指向矢n、速度v对时间和空间的偏导数

根据中心差分法，用差商代替微商，计算 $n,v$ 对 $t,y$ 的一阶、二阶及混合偏导数，并将计算各值代入式7c~式7e计算 $n$ 对时间 $t$ 的偏导数，再用式(8)进行微调。将算得的 $n$ 对时间 $t$ 的偏导数代入式7a和式7b计算 $v$ 对时间 $t$ 的偏导数，并保存。根据上两步计算的结果作为k₁、k₂的值.

2) 计算下一时刻n、v值

将n、v对时间、空间的偏导数和n、v旧值以及k₁、k₂的值代入式9计算下一个节点n、v的值，并替换原来的值。

3) 以0.01 s为单位保存n、v值，并作图

首先判断现在n、v值是否需要输出保存，若程序总共要运行1 s，那么每隔0.01 s保存一个数据即满足作图要求。对多周期计算，同样需01 s保存100个值。然后，判断程序是否执行完，若已执行完，用所得数据作图。所得图形可明显看到不同波形、不同盒厚、不同场强、电场不同频率、电场不同占空比对液晶微流动速度的影响。

4. 程序优化

整个程序每执行一次，需计算10⁹次循环体(计算1 s时)，用4核4G的计算机初次运行程序，大约需要约912.83个小时。根据程序特点着重从以下几个方面对程序进行优化。

4.1. 优化变量定义

matlab可以自动调整变量的大小，且对向量化的指令非常敏感，因此采用预先为变量分配内存空间的方法可使调用内存的指令次数降为1，大大减少计算的时间。如运用以下程序段对变量nx, ny, nz, vx, vz, dn1dt定义大小和分配内存。

程序段1：n1=zeros(m+1,2,'double');

v1=zeros(m+1,2,'double');

dn1dt=zeros(1,m+1,'double');

dv1dt=zeros(1,m+1,'double');

pddn1dt=zeros(m+1,2,'double');

pddv1dt=zeros(m+1,2,'double');

nx=zeros(m+1,kmax/kstep,'double');

u=zeros(m+1,kmax/kstep,'double');

通过对变量n1、n2、n3、v1、v2、dn1dt、dn2dt、dn3dt、dv1dt、dv2dt、pddn1dt、pddn2dt、pddn3dt、pddv1dt、pddv2dt、nx、ny、nz、u、w预定义大小，可大大降低运算时间，提高程序段的运算速度。

4.2. 优化程序段

按照上述方法，先定义程序中所有变量的大小，再计算每个程序段执行1次的时间。可得从第1步到第4步每步运行一次时间分别为0.004305S、0.071581S、0.000153S、0.000222S。由此可见，第2步执行一次所需时间最长。

由于第1步程序段在整个程序执行1次中，仅执行1次，而循环体内部子程序执行10⁹次，因此整个程序的优化重点应该在循环体部分，即程序编写过程中的第2、3、4步。从上述每个程序段的执行一次所需时间可知：第2步程序运行时间最长，是第3、4步运行之和的190多倍，第2步是一个二重循环(2阶龙格-库塔法)，所以第2步应该是整个程序优化重点。

第2步中需要计算两次的子程序包括4个部分，分别为“计算 $n$ 对时间 $t$ 的偏导数”子程序、“ $n$ 对时间 $t$ 的偏导数微调”子程序、“计算 $v$ 对时间 $t$ 的偏导数”子程序和“n、v中间值计算”子程序。每个子程序执行一次所需时间分别为：0.002455S、0.002756S、0.015191S和0.003936S。由此可见四个子程序执行时间相差不大，应该分别对其进行优化。注意，在MATLAB进行计算时，同一台电脑执行同一段程序运行时间不完全相同，但相差不大，这是因为某两个时刻电脑执行的任务不尽相同，但是差别不会太大，因此可以把运行时间作为优化程序参数。

4.2.1. 优化程序中的子程序

本程序由于循环体必须按照顺序执行，不能采用并行计算，但是对于循环体内部的子函数可以分块优化。对于子程序的优化，主要是改变程序的编写方式，也就是改变计算方法，常用的编程方法主要有“串行计算(for循环)”、“并行计算(parfor循环)”、“向量计算” [14] 。

以循环次数为横坐标，计算程序段的运行时间为纵坐标，分别对并行计算，for循环以及向量计算进行绘图。对于“计算 $n$ 对时间 $t$ 的偏导数”子程序，计算结果如图2所示，可见对于此段子程序，向量计算最优；图3为对“ $n$ 对dnidt(i = 1, 2, 3)的偏导数微调”子程序计算结果，可见对于本子程序for循环计算最优；图4为对“计算 $v$ 对时间 $t$ 的偏导数”子程序计算结果，可见对本子程序for循环计算最优。

在并行计算中，对循环体变量有一定的限制。并行计算的变量分为：循环体变量、Sliced变量、临时变量等五类。本文中用到的是循环体变量。在使用循环体变量时，要求循环体必须具有一定的独立性、循环体之间也不相关，在“n、v中间值计算”的子程序中的循环变量不符合并行计算对变量的要求，因此对次子程序的优化不考虑并行计算，重点针对for循环计算和向量计算方法。如图5所示，此为“n中间值计算”子程序计算结果，可见在循环次数非常大时，for循环计算最优。如图6为“v中间值计算”子程序计算结果，同样是for循环计算时间较短。

图2. “计算n对时间t的偏导数”子程序优化

图3. n对dnidt(i = 1, 2, 3)微调子程序优化

图4. “计算v对时间t的偏导数微调”子程序优化

图5. “n中间值计算”子程序优化

图6. “v中间值计算”子程序优化

4.2.2. 优化程序段中单的行程序

在数值计算中，优化每一行的程序同样会减少整个程序运行时间，达到非常好的优化效果。如本程序中，在“计算v对t的偏导数”的子程序中对 $n_y^2$ 共有8次计算(由式7a和式7b可知)。若在进行计算式7a和式7b前，就已计算并保存了对 $n_y^2$ 的计算结果，则单次循环可以少计算乘法7次。对10⁹次的整个循环体，相当于少计算乘法1.4 × 10¹⁰次(2阶龙格–库塔法计算两次子程序)。

4.3. 优化结果

经过上述程序优化过程，在同样的硬件条件下(普通4核4G)将原程序计算1 s所需时间由原来的约912.83个小时，减少到约47.37小时，计算效率提高近20倍。

5. 结论

针对液晶引流驱动中的微流动引流的计算问题，采用MATLAB软件辅助编程及数值计算方法，并针对程序特点分别使用了优化变量定义，程序段及单行程序段方法对所编写程序进行了优化处理。通过优化计算发现for循环的计算方法适合进行本程序的计算，相同的计算硬件及时长条件下，计算时间由原来的约912.83小时优化到约47.37个小时，计算效率提高近20倍。

基金项目

国家自然科学基金资助项目(11372003)；河南省高校科技创新人才支持计划(14HASTIT007)。

文章引用

林绿叶,刘春波,王安标. 电场作用下液晶微流动数值计算程序编写与优化
The Written and Optimization of the Liquid Crystalline Micro Flow Simulation Program under the Electric Field Effect[J]. 计算机科学与应用, 2017, 07(10): 974-983. http://dx.doi.org/10.12677/CSA.2017.710110

参考文献 (References)