一种流体模拟方法、装置、系统及存储介质与流程

1.本技术涉及电润湿技术领域，尤其涉及一种流体模拟方法、装置、系统及存储介质。


背景技术：

2.相关技术中，电润湿器件，也称为电湿润器件或电流体器件，因为具有低能耗、高效率的特点，在微流控以及电子显示领域应用广泛。电润湿器件在研究开发的过程中，为了提高效率，通常会对电润湿器件进行流体的模拟仿真。而流体的模拟仿真，会应用到数学、流体力学以及计算机，因此，计算流体力学，作为介于数学、流体力学和计算机之间的交叉学科，在电润湿器件的模拟仿真过程中得到了广泛应用。
3.目前，通过计算流体力学对电润湿器件进行模拟仿真，通常采用有限元方法以求解流体的运动状态方程，如以有限元方法求解纳维斯托克斯方程，会采用cfd软件(仿真计算流体力学软件)进行有限元法的相关数值计算。但是，有限元方法是一种对流体运动方程进行细致微分处理的方法，此种模拟方法对于计算机的性能要求较高，因此会产生较高的模拟仿真成本，同时，微分处理过程的计算量较大，计算机需要花费较长的时间进行处理，会造成建模时间较长，降低了模拟仿真的效率；此外，此种同时模拟方法对于模型设计者的经验要求较高，对于初次使用者需要较长的时间进行磨合，增加了时间上的成本，也进一步降低了模拟仿真的效率。


技术实现要素：

4.本技术旨在至少解决现有技术中存在的技术问题之一。为此，本技术提出一种流体模拟方法、装置、系统及存储介质，能够更加准确的模拟流体的电润湿现象，模拟效果更加接近实际。
5.根据本技术第一方面实施例的流体模拟方法，包括：
6.获取模拟数据集；其中，所述模拟数据集包括流体的第一模拟物理量和待电湿润的器件的多个第二模拟物理量；
7.根据第一物理模拟量和第二物理模拟量定义器件的初始物理场分布；其中，物理场为二维网络模型，所述二维网格模型中的每一格点对应不同的坐标；
8.对所述模拟数据集执行循环过程直至满足预设条件，并将满足预设条件时的当前物理场分布确定为稳态分布状态或目标状态，所述循环过程包括：
9.根据所述第一模拟物理量、所述第二模拟物理量、此前物理场分布中的此前相场分布计算的相场、所述此前物理场分布中的此前密度场分布计算的密度场，计算得到化学势和自由能密度函数；
10.根据所述化学势、所述自由能密度函数、所述相场以及所述密度场，计算得到所述流体在每一格点的第一压强；
11.根据所述二维网格模型上三相接触区域对应的相场、所述此前物理场分布中的此
前电势分布，修正所述三相接触区域对应的第一电场强度，以确定每一格点的第二电场强度；
12.将所述第一压强、所述第二电场强度进行耦合，得到每一格点的第二压强；
13.根据所述第二压强，计算得到所述流体的当前密度平衡分布函数、当前相平衡分布函数；
14.根据所述当前相平衡分布函数、所述当前密度平衡分布函数，计算得到每一格点上所述流体的所述当前物理场分布；其中，所述当前物理场分布包括当前相场分布、当前密度场分布、当前电势分布；
15.将当前相场变化量与预设变化量比较或将当前循环次数与预设循环次数比较；
16.其中，所述预设条件为所述当前相场变化量小于所述预设变化量或所述当前循环次数等于所述预设循环次数。
17.根据本技术的一些实施例，所述获取模拟数据集，包括以下步骤：
18.获取所述流体的多个第一实际物理量、所述器件的多个第二实际物理量；其中，多个第二实际物理量包括器件横向长度l
x
、器件纵向长度ly、固体介电层厚度ds、电极位置、格点大小δ
x
、时间步长δ
t
、电极的位置和电势v
si
；
19.对每一所述第一实际物理量、每一所述第二实际物理量进行无量纲归一化处理，得到与所述第一实际物理量对应的所述第一模拟物理量、与所述第二实际物理量对应的第二模拟物理量；
20.根据每一所述第一模拟物理量、每一所述第二模拟物理量，建立所述模拟数据集。
21.根据本技术的一些实施例，根据第一物理模拟量和第二物理模拟量定义的初始物理场分布，包括：
22.定义初始速度、初始相场、初始电势。
23.根据本技术的一些实施例，所述根据所述化学势、所述自由能密度函数，计算得到所述流体在每一格点的第一压强，包括：
24.根据预设的润湿势、所述相场对应的自由能密度、所述自由能密度函数、所述化学势、所述第一模拟物理量、所述第二模拟物理量、所述密度场，计算得到所述流体在每一格点的所述第一压强。
25.根据本技术的一些实施例，所述二维网格模型上三相接触线区域对应的相场、所述此前物理场分布中的此前电势分布，修正所述三相接触线区域对应的第一电场强度，以确定每一格点的第二电场强度，包括：
26.当三相接触线区域对应的相场落入预设范围，根据所述此前电势分布、预设的增强公式对所述相场对应的所述第一电场强度进行增强处理，得到每一格点的所述第二电场强度；其中，所述增强公式为：
27.e(φ)＝e(φ)*b；其中，e表示电场强度，φ表示相场，e(φ)表示所述相场与所述电场强度的关系式，b表示增强系数，θ0表示流体在固体界面稳定时的接触角，δ
x
表示每一格点的尺寸。
28.根据本技术的一些实施例，所述将所述第一压强、所述第二电场强度进行耦合，得到每一格点的第二压强，包括：
29.根据公式获取联合介电常数；其中，ε表示联合介电常数，φ表示相场，ε1表示第一介电常数，ε2表示第二介电常数；
30.根据所述联合介电常数、所述第二电场强度、所述第一压强，计算得到所述流体在每一格点的所述第二压强。
31.根据本技术的一些实施例，在所述对所述模拟数据集执行循环过程直至满足预设条件，并将满足预设条件时的当前物理场分布确定为稳态分布状态或目标状态之前，包括：
32.当所述格点为边界点，对密度分布函数中的速度参数进行反弹算法处理，以修正所述边界点的所述此前密度场分布。
33.根据本技术的一些实施例，在所述对所述模拟数据集执行循环过程直至满足预设条件，并将满足预设条件时的当前的物理场分布确定为稳态分布状态或目标状态之前，还包括：
34.当所述格点为边界点，对相分布函数中的速度参数进行反弹算法处理，以修正所述边界点的所述此前相场分布。
35.根据本技术第三方面实施例的流体模拟系统，包括：
36.至少一个存储器；
37.至少一个处理器；
38.至少一个程序；
39.所述程序被存储在所述存储器中，所述处理器执行至少一个所述程序以实现如第一方面实施例所述的流体模拟方法。
40.根据本技术第四方面实施例的计算机可读存储介质，所述计算机可读存储介质存储有计算机可执行指令，所述计算机可执行指令用于使计算机执行如第一方面实施例所述的流体模拟方法。
41.根据本技术实施例的流体模拟方法，至少具有如下有益效果：首先，获取模拟数据集；其中，模拟数据集包括流体的多个第一模拟物理量和待电湿润的器件的多个第二模拟物理量；根据第一物理模拟量和第二物理模拟量定义器件的初始物理场分布；其中，物理场为二维网络模型，二维网格模型中的每一格点对应不同的坐标；之后，对模拟数据集执行循环过程直至满足预设条件，并将满足预设条件时的当前物理场分布确定为稳态分布状态或目标状态，循环过程包括：根据第一模拟物理量、第二模拟物理量、此前物理场分布中的此前相场分布计算的相场、此前物理场分布中的此前密度场分布计算的密度场，计算得到化学势和自由能密度函数；根据化学势、自由能密度函数、相场以及密度场，计算得到流体在每一格点的第一压强；根据二维网格模型上三相接触区域对应的相场、此前物理场分布中的此前电势分布，修正三相接触区域对应的第一电场强度，以确定每一格点的第二电场强度；将第一压强、第二电场强度进行耦合，得到每一格点的第二压强；根据第二压强，计算得到流体的当前密度平衡分布函数、当前相平衡分布函数；根据当前相平衡分布函数、当前密度平衡分布函数，计算得到每一格点上流体的当前物理场分布；其中，当前物理场分布包括当前相场分布、当前密度场分布、当前电势分布；将当前相场变化量与预设变化量比较或将当前循环次数与预设循环次数比较；其中，预设条件为当前相场变化量小于预设变化量或当前循环次数等于预设循环次数。本技术的流体模拟方法，通过将第二电场强度和第一压
强进行耦合，得到第二压强，然后就可以通过第二压强计算得到密度分布；一方面，较大程度的简化了计算过程，通过计算量的简化，进而降低了对计算机的高性能要求；另一方面，操作者只需对相应的实际物理量进行设置就可以进行仿真模拟，较大程度的降低了操作者的要求，节省了学习的时间成本；最后一方面，通过修正电场强度后再对电场与流体场进行耦合，模拟仿真得到的参数更加准确，使得模拟流体在器件内的电润湿现象更加接近实际。因此，本技术的流体模拟方法，能够更加快速而准确地模拟流体的电润湿现象，模拟效果更加接近实际。
42.本技术的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本技术的实践了解到。
附图说明
43.下面结合附图和实施例对本技术做进一步的说明，其中：
44.图1为本技术一个实施例所提供的流体模拟系统的结构示意图；
45.图2为本技术实施例所提供的流体模拟方法的流程示意图；
46.图3为本技术实施例所提供的电场强度修正前后的电压和湿润角的变化曲线对比示意图；
47.图4为本技术实施例所提供的介电常数分布与相场关系图；
48.图5为本技术另一实施例所提供的流体模拟系统的结构示意图。
49.附图标记：
50.导电流体100、介电流体110、固体介电层120、存储器200、处理器300。
具体实施方式
51.下面详细描述本技术的实施例，实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本技术，而不能理解为对本技术的限制。
52.在本技术的描述中，需要理解的是，涉及到方位描述，例如上、下、前、后、左、右等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本技术和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本技术的限制。
53.在本技术的描述中，若干的含义是一个以上，多个的含义是两个以上，大于、小于、超过等理解为不包括本数，以上、以下、以内等理解为包括本数。如果有描述到第一、第二只是用于区分技术特征为目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量或者隐含指明所指示的技术特征的先后关系。
54.本技术的描述中，初始的含义为模拟计算开始的时刻，此前的含义是上个循环结束时的时刻，当前的含义是指当前循环所对应的时刻。
55.本技术的描述中，除非另有明确的限定，设置、安装、连接等词语应做广义理解，所属技术领域技术人员可以结合技术方案的具体内容合理确定上述词语在本技术中的具体含义。
56.本技术的描述中，参考术语“一个实施例”、“一些实施例”、“示意性实施例”、“示
例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本技术的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
57.下面对本技术中的名词进行解释。
58.电润湿：
59.是指通过改变液体与绝缘基板之间的电压，来改变液滴在基板上的润湿性(即改变接触角)，使液滴发生形变和运动的现象。电润湿显示器件最早由飞利浦门下的专业显示器厂商liquavista公司研制出来的，其研究成果发表在2003年《nature》杂志上。电润湿显示器件通常由上下两个结构层组成，上结构层由上基板、导电层、支撑柱和密封胶框组成，下结构层由下基板、导电层、疏水绝缘层、弱亲水像素墙组成。像素墙的图案限定了显示器件的像素，像素墙之间的区域为显示区域，像素中填充有油墨和电解质溶液两种互不相溶的流体，通过施加电压来使两种流体的运动来产生显示效果。具体地，当没有施加电压的时候，油墨铺展在憎水层的表面，显示油墨的颜色；当施加电压时，油墨收缩，显示下基板的颜色。
60.格子玻尔兹曼：
61.格子玻尔兹曼方法(lattice boltzmann method,lbm)是一种基于介观模拟尺度的计算流体力学方法。该方法相比于其他传统cfd计算方法，具有介于微观模型和宏观连续模型的介观模型特点，因此具备流体间相互作用描述简单、边界条件易于设置、可并行计算等优势。
62.自由能模型：
63.自由能模型基于相场理论，其中自由能泛函用于解释界面张力效应，并以热力学一致的方式描述界面动力学的演化。
64.young-lippmann方程：
65.描述了在施加电势前后平行基板上的导电液滴的润湿角变化，其中θ和θ0分别代表液滴存在电势差后的稳定润湿角以及没有电势存在时的平衡润湿角。此方程中，ε0和εs分别表示真空中介电常数以及固体介电层的相对介电常数，v是正负电极的电势差，γ为两流体界面的表面张力，ds为固体介电层的厚度。
66.以格子玻尔兹曼理论为基础的数值模拟方法用于模拟流体流动和流体物理建模，与传统的基于求解宏观连续流体方程的数值方法不同，其基本思想是建立简化的流体动力学模型，使得计算速度有很大的提升。
67.多物理场耦合模型的cfd软件(仿真计算流体力学软件)如comsol、ansys，或者开源的openfoam等。以有限元方法求解流体的运动状态方程如纳维斯托克斯方程(n-s方程)的偏微分解。以comsol为例，有限元数值计算的基本流程为：首先，将求解域划分成若干有限元网格；然后，运用插值基函数构造待求解量的近似解；最后，运用加权余量法或变分原理将偏微分方程进行转化，使其成为离散形式的代数方程。应用有限元法进行数值求解时，运用节点插值可得到求解域内任意点的场量。此种对流体运动方程进行细致微分处理的方法，其模拟方式对于计算机性能要求高，计算效率低，同时建模时网格划分对于模型设计者
要求较高的应用经验。
68.下面参照图1描述根据本技术实施例的流体模拟系统。
69.可以理解的是，如图1所示，流体模拟系统，包括：
70.待电湿润的器件，器件设置有固体介电层120；
71.流体，流体包括导电流体100和介电流体110；
72.固体介电层120、导电流体100和介电流体110共同构成了一个流体场。
73.在这个流体场中，格点大小δ
x
、时间步长δ
t
、电极的位置、电势v
si
、器件的横向长度l
x
、器件纵向长度ly、固体介电层厚度ds作为第二实际物理量，需要在这里进行定义并赋值，具体参数通过器件的初始技术方案进行参考，也可以理解为进行模拟的初步方案；同理，在这个流体场中，流体的流体密度ρ
si
、流体接触界面迁移率m
si
、两个流体接触面的表面张力γ
si
以及运动粘度v
si
，两流体接触面的界面厚度l
si
，流体在固体界面稳定接触角θ0、导电流体的介电常数ε1和介电流体的介电常数ε2、流体的尺寸r、以及固体介电层的介电常数εs作为流体的第一实际物理量，都需要在这里进行定义并赋值。值得注意的是，上述物理量均采用国际标准单位制定义。
74.下面参照图2描述根据本技术实施例的流体模拟方法。
75.可以理解的是，如图2所示，流体模拟方法，包括：
76.步骤s100，获取模拟数据集；其中，模拟数据集包括流体的第一模拟物理量和待电湿润的器件的多个第二模拟物理量。
77.可以理解的是，获取模拟数据集，包括以下步骤：
78.获取流体的多个第一实际物理量、器件的多个第二实际物理量；其中，如图2所示，多个第二实际物理量包括器件横向长度l
x
、器件纵向长度ly、固体介电层厚度ds、格点大小δ
x
、时间步长δ
t
、电极的位置、电势v
si
；
79.对每一第一实际物理量、每一第二实际物理量进行无量纲归一化处理，得到与第一实际物理量对应的第一模拟物理量、与第二实际物理量对应的第二模拟物理量；
80.根据每一第一模拟物理量、每一第二模拟物理量，建立模拟数据集。
81.需要说明的是，格子玻尔兹曼的模拟参量均为无量纲数，为将实际器件单位与模拟参量相结合，需引入参考单位制。对每个模拟的基本物理量设置对应的归一化参考值，则模拟物理量等于实际物理量除以归一化参考值所得到的无量纲数。设置物理量的归一化参考值ρ
ref
、迁移率参考值m
ref
、运动粘度参考值ν
ref
＝(δ
x
)2/δ
t
、表面张力参考值γ
ref
＝ρ
ref
*(δ
x
)3/(δ
t
)2以及电压参考值v
ref
，参考值是为了对模拟物理量进行归一化的处理，选取的参考值的具体大小原则上应保证并不影响实际模拟效果。
82.需要说明的是，对第一模拟物理量、第二模拟物理量进行定义：
83.第一模拟物理量包括：流体密度ρ(x,t)、运动粘度ν、流体的介电常数ε、流体接触界面迁移率m、表面张力γ、两流体接触面的界面厚度l、弛豫时间τ；各个固体边界对流体的润湿性θ0、流体的尺寸r；其中，τ表示粒子分布由非平衡态达到平衡态的时间。
84.第二模拟物理量包括：器件横向长度l
x
、器件纵向长度ly、固体介电层厚度ds、电极的位置、电势v
app
；其中，电极包括正电极和负电极。
85.对所有第一实际物理量、第二实际物理量进行无量纲归一化处理，得到第一实际物理量对应的第一模拟物理量、第二实际物理量对应的第二模拟物理量，无量纲的物理量
等于对应的实际器件的物理量除以归一化参考值，物理量a＝a
si
/a
ref
。例如，运动粘度ν＝ν
si
/ν
ref
；迁移率m＝m
si
/m
ref
，密度ρ(x,t)＝ρ
si
/ρ
ref
，电极电势v
app
＝v
si
/v
ref
，依次无量纲化处理后，获得以上各个模拟物理量的无量纲值。例如，实际中第一流体的密度为ρ
si
＝800kg/m3，则密度参考值设置为ρ
ref
＝1000kg/m3，则数值模拟时第一流体的密度ρ为0.8。
86.步骤s110，根据第一物理模拟量和第二物理模拟量定义器件的初始物理场分布；其中，物理场为二维网络模型，二维网格模型中的每一格点对应不同的坐标。
87.需要说明的是，通过步骤s100，已经初步搭建了流体模拟系统的二维模型，此时，对二维模型进行均一化网格处理，得到二维网络模型，根据单个网格尺寸δ
x
、基本时间步δ
t
、矩形模型的边界位置、两种流体在模型中的初始位置、固体材料组分的形状与位置、电极位置。获得模拟网格分布，将网格中格点定义为离散化运动的粒子位置。其中，x、y分别表示模拟网格i(x,y)上网格点的横坐标和纵坐标，x∈(1,l
x
/δ
x
，y∈(1,ly/δ
x
)。下面的表述中用x表示位置坐标(x,y)。
88.值得注意的是，初始物理场分布包括初始密度场分布、初始相场分布以及初始电势分布，物理场包括密度场、相场、速度场、电势；密度场、相场、速度场、电势在设定的区域内分布，则称之为物理场分布，进一步地，在设定的区域内建立网络，就可以得到二维网络模型。
89.可以理解的是，根据第一物理模拟量和第二物理模拟量定义器件的初始物理场分布，包括：
90.定义初始速度、初始相场、初始电势。
91.值得注意的是，通过定义初始速度、初始相场、初始电势，就可以得到初始密度场分布、初始相场分布以及初始电势分布。
92.需要说明的是，在模拟之前，需要对参数的初始值进行定义，定义相场的初始状态，初始相场：
93.其中x0＝(l
x
/2,r0)是中心液滴的初始圆心位置，r0为液滴半径。
94.定义初始速度u(x,t＝0)＝0，初始的电极电势v
app
(x,t＝0)＝0，初始的电极电势v
app
随时间t变化规律如下：
[0095][0096]
为了保证时电场迭代的正确性，初始的电极电势v
app
逐步增大到v
max
，且所用的时间为t
vsetp
。
[0097]
需要说明的是，v
max
为每次模拟时预设的初始的电势v
app
的实际最大值的无量纲化后的数值。
[0098]
对模拟数据集执行循环过程直至满足预设条件，并将满足预设条件时的当前物理场分布确定为稳态分布状态或目标状态，循环过程包括：
[0099]
步骤s120，根据第一模拟物理量、第二模拟物理量、此前物理场分布中的此前相场分布计算的相场、此前物理场分布中的此前密度场分布计算的密度场，计算得到化学势和自由能密度函数。
[0100]
需要说明的是，此前物理场分布，在第一次进行循环时，可以理解为初始物理场分布。
[0101]
需要说明的是，使用格子-玻尔兹曼理论中的自由能模型处理两种非相溶流体，为对两种流体进行区分，定义导电流体和介电流体的相场φ分别为φ＝1、φ＝-1，以在进行数值模拟时快速区分两种流体。-1《φ《1时，为两流体的混合态，主要集中在两相流体的接触界面区域。
[0102]
需要说明的是，使用d2q9模型(二维、9个离散的速度方向)定义流体的离散物理量，流体运动速度为ei，i为粒子速度的方向序列(i∈[0,8]的9个整数)，表示粒子运动的9种速度方向，进而得到如下速度公式：
[0103][0104]
ωi表示各个方向对应的权重系数。依照热力学守恒ω0＝4/9，ω1＝ω2＝ω3＝ω4＝1/9，
[0105]
ω5＝ω6＝ω7＝ω8＝1/36。设置格子声速cs：
[0106]
需要说明的是，定义初始物理场分布时，定义的是相场、密度场、速度场、电场的初始值，也就是初始相场、初始密度场、初始速度场、初始电场，因而需要根据相场计算相分布函数gi(x,t)，相分布函数gi(x,t)的含义，可以理解为表示在各时间点各个格点x处存在的流体的相场φ的概率值；并且，相平衡分布函数g
ieq
(x,t)的含义，可以理解为表示在系统处于平衡态时各个格点对应的位置x处流体的相场φ的值；此外，相平衡分布函数也是在定义相场后计算出来的，需要根据化学势、相场等相关参数进行计算。根据相分布函数和相平衡分布函数，得到的相场演化公式如下，也称之为相场分布：
[0107][0108]
通常定义相场弛豫时间：τg＝1；并且，在相场演化公式中，为碰撞过程，gi(x+eiδ
t
,t+δ
t
)-gi(x,t)为迁移过程。
[0109]
宏观层面，离散的相场分布与宏观量之间遵循如下限制：
[0110][0111][0112][0113]
在公式中，ei和u均为矩阵表示，则和u
t
均为矩阵运算中的转置，i表示模拟网络，φ为相场，μ为化学势，m'为黏度项，因此，相平衡分布函数根据公式以及公式
中的参数，例如相场φ、化学势μ，进行综合计算得到。
[0114]
相场的梯度和拉普拉斯函数定义了相场与时间t和位置x的关系，表达为：
[0115][0116][0117]
需要说明的是，定义u为流体的速度场u(x,t)，表示各个格点的速度分布；
[0118]
需要说明的是，计算速度场需要根据相场平衡分布函数与宏观相函数的限制关系式或者密度场平衡分布函数与宏观密度的限制关系式。
[0119]
值得注意的是，相平衡分布函数是流体状态达到稳定平衡态时的理想值。流体粒子通过不断的碰撞和迁移运动，从当前状态趋向于平衡状态。也就是说，当前的相分布函数的值通过不断的演化趋向于相平衡分布函数的值，最终达到稳定状态。也就是下文中提到的，当前相场变化量小于预设变化量(微分阈值)时，可以认为此时的相分布函数和相平衡分布函数无限接近了，即流体达到稳定状态。
[0120]
需要说明的是，同理，在定义初始密度场后，需要根据密度场分别计算密度分布函数fi(x,t)、密度平衡分布函数f
ieq
(x,t)；对于密度分布函数fi(x,t)的含义，可以理解为表示流体粒子在各个时间各格点对应的位置x处的密度ρ分布概率值；对于密度平衡分布函数f
ieq
(x,t)的含义，可以理解为表示流体粒子在各个时间各格点对应的位置x处的密度ρ的值；密度分布函数fi(x,t)以及密度平衡分布函数f
ieq
(x,t)需要满足与宏观密度ρ(x,t)之间的限制关系：
[0121][0122][0123][0124]
格点中的粒子运动分为碰撞和迁移两部分，碰撞和迁移不断进行，则流体场随着时间不断演化运动直至达到粒子分布的平衡态，流体场表现为宏观稳定状态，表示碰撞迁移运动的碰撞迁移公式，也就是密度场分布，如下所示：
[0125][0126]
其中fi(x,t)为t时刻粒子在x处的第i个方向的分布状态，fi(x+eiδ
t
,t+δ
t
)为δ
t
时刻后，粒子经碰撞迁移到x+eiδ
t
位置i方向的分布函数，流体的运动粘度满足表示碰撞过程，fi(x+eiδ
t
,t+δ
t
)-fi(x,t)表示迁移过程。
[0127]
需要说明的是，fi(x,t)在初始t＝0时等于f
ieq
(x,t)，则此刻fi(x+eiδ
t
,t+δ
t
)下一个时刻(t＝1)的fi(x,t＝1)，则由碰撞迁移公式计算得到，依次类推fi(x,t)和fi(x+eiδ
t
,t+
δ
t
)均由此依次迭代计算得到，上述gi(x,t)具有相同的计算方法。
[0128]
值得注意的是，密度平衡分布函数是流体状态达到稳定平衡态时的理想值。流体粒子通过不断的碰撞和迁移运动，从当前状态趋向于平衡状态。也就是说，当前的密度分布函数的值通过不断的演化趋向于密度平衡分布函数的值，最终达到稳定状态。也就是下文中提到的，当前相场变化量小于预设变化量(微分阈值)时，可以认为此时的密度分布函数和密度平衡分布函数无限接近了，即流体达到稳定状态。
[0129]
需要说明的是，上文中已经介绍了通过此前相场分布计算相场，通过此前密度场分布计算密度场，对于化学势和自由能密度函数，计算过程如下：
[0130]
根据流体的相场、第一模拟物理量、第二模拟物理量，计算得到的化学势；化学势函数如下：
[0131]
化学势：
[0132]
需要说明的是，根据流体的相场、第一模拟物理量、第二模拟物理量、流体的密度场，计算得到自由能密度函数；其中，密度场通过压强和自由能产生间接关系，因此，需要考虑到密度场对自由能的影响；自由能密度函数如下：
[0133][0134]
根据流体的相场、化学势、第一模拟物理量、第二模拟物理量、自由能密度函数，就可以计算得到无电场时，流体场的自由能f：
[0135][0136]
其中，ω为流体区域体积，s为流体区域表面积。
[0137]
润湿势：
[0138][0139]
决定了流体在固体介电层表面的润湿能力的大小。
[0140]
步骤s121，根据化学势、自由能密度函数、相场以及密度场，计算得到流体在每一格点的第一压强。
[0141]
需要说明的是，根据相场、密度场、第一模拟物理量、第二模拟物理量、化学势、自由能密度函数，就可以计算得到流体在每一格点的第一压强。
[0142]
流体的压强张量p是影响流体运动的力学因素，是能量的体密度，在不加电场时，为第一压强，表达式为：
[0143][0144]
步骤s122，根据二维网格模型上三相接触区域对应的相场、此前物理场分布中的此前电势分布，修正三相接触区域对应的第一电场强度，以确定每一格点的第二电场强度。
[0145]
需要说明的是，引入电势分布函数hi(x,y,t)，以表示电势在t时刻在以x表示(x,y)位置且为i方向的电势变化概率：
[0146][0147]
为平衡态时分布函数，也可以称之为电势平衡分布函数：τh＝1。
[0148]
设置初态时的电势v(x,t＝0)，则h
ieq
(x,t＝0)，设置初态时hi(x,t＝0)＝h
ieq
(x,t＝0)。
[0149]
需要说明的是，设置导电流体相与介电流体相的相场阈值φ0，实际情况下，两相流体电导率相差较大时，不能以φ＝0的平均中间态，作为区分两相流体混合区域的流体导电性的参考。根据实际电润湿器件中两相流体的电导率，决定相对应的相场阈值φ0。进一步地，两相流体分别为相对介电常数为78的水和相对介电常数为2的介电油墨，设置阈值φ0＝0.45，当φ≥0.45时，认为流体为介电流体，φ《0.45时，认为流体为导电流体。
[0150]
对于介电流体相的电势∑
ihi
＝v，介电流体遵循介电流体中电势v的分布遵循：
[0151]
表示不同的v对应的φ的阶段函数。
[0152]
定义导电流体是完美导体，对于导电流体相的电势v等于电极电势v
app
，则当定义导电流体是完美导体时，导电流体相的电场强度根据计算，得到导电流体相的电场强度e＝0。
[0153]
需要说明的是，介电流体中双电层的德拜长度远小于格子长度，电势差u为器件中电势分布不均一造成的区域差异。
[0154]
需要说明的是，电势平衡分布函数表示电势在t时刻在以x表示(x,y)位置且为i方向的电势；通过电势平衡分布函数来表示t时刻的电势分布，电势平衡分布函数。
[0155]
可以理解的是，根据二维网格模型上三相接触区域对应的相场、此前物理场分布中的此前电势分布，修正三相接触区域对应的第一电场强度，以确定每一格点的第二电场强度，包括：
[0156]
当三相接触区域对应的相场落入预设范围，根据此前电势分布、预设的增强公式对相场对应的第一电场强度进行增强处理，得到每一格点的第二电场强度；其中，增强公式为：e(φ)＝e(φ)*b；其中，e表示电场强度，φ表示相场，e(φ)表示所述相场与所述电场强度的关系式，b表示增强系数，θ0表示流体在固体界面稳定时的接触角，δ
x
表示每一格点的尺寸。
[0157]
可以理解的是，第一电场强度按照如下公式计算：
[0158][0159]
可以理解为根据电势计算电势差u，并对电势差进行微分计算，得到第一电场强度。
[0160]
需要说明的是，相场φ的预设范围为φ0∈(0.45,0.8)，φ0称之为阈值，当相场落入预设范围后，通过结合相场φ的大小以及电势分布，就可以对对三相接触区域的第一电场强度做增强处理，得到第二电场强度。
[0161]
需要说明的是，对于流体中不在三相接触线附近的区域，即(0.45,0.8)的部分，其电场强度e(φ)不做上述调整，依旧遵循原数值解。
[0162]
值得注意的是，通过上述修正后得到第二电场强度，当然，未做修正的第一电场强度也作为第二电场强度计算第二压强，修正前后的效果如图3所示。
[0163]
步骤s123，将第一压强、第二电场强度进行耦合，得到每一格点的第二压强。
[0164]
可以理解的是，将第一压强、第二电场强度进行耦合，得到每一格点的第二压强，包括：
[0165]
根据公式获取联合介电常数；其中，ε表示联合介电常数，φ表示相场，ε1表示第一介电常数，ε2表示第二介电常数；
[0166]
根据联合介电常数、第二电场强度、第一压强，计算得到流体在每一格点的第二压强。需要说明的是，不同于以往两相流体中的介电常数关于相场的的线性分布关系本技术对介电常数与相场的关系做了调整，调整后的介电常数与流体相的关系如图4所示。
[0167]
需要说明的是，引入电场后，流体的第二压强为描述流体受到的应力，表达式为：
[0168][0169]
步骤s124，根据第二压强，计算得到流体的当前密度平衡分布函数、当前相平衡分布函数。
[0170]
需要说明的是，根据第二压强计算得到流体的当前密度平衡分布函数的计算表达式如下：
[0171][0172][0173]
此处，p为第二压强，且为二阶张量，
[0174]
值得注意的是，相平衡分布函数并不直接由第二压强决定，而是间接受到第二压强影响，因压强决定了流体的应力从而间接影响了两相流体的运动趋势即分布趋势。相平衡分布函数的表达式为：
[0175][0176]
其中，
[0177]
需要说明的是，当前密度平衡分布函数和当前相平衡分布函数表示的是依据当前压强状态下计算得到的密度场和相场稳定分布时的理想数值解，并非当前的密度场和相场的分布状态也不是计算达到稳态分布时的解。其作用在于下一时间步的密度场和相场的分布状态由此时刻的实际分布状态经碰撞和迁移后决定，而其碰撞和迁移的计算过程中将使用当前的密度场和相场的平衡状态。
[0178]
步骤s125，根据当前相平衡分布函数、当前密度平衡分布函数，计算得到每一格点上流体的当前物理场分布；其中，当前物理场分布包括当前相场分布、当前密度场分布、当前电势分布。
[0179]
需要说明的是，当前物理场分布还包括当前速度分布，同时，根据上述公式可知，当前速度分布u(x,t)分别和当前相场分布、当前密度场分布密切相关，而当前相场分布又会影响到当前电势分布；而前文中有提到计算速度场需要根据相场平衡分布函数与宏观相函数的限制关系式或者密度场平衡分布函数与宏观密度的限制关系式，因此，当前速度分布与当前相场分布、当前密度场分布相互影响。
[0180]
步骤s126，将当前相场变化量与预设变化量比较或将当前循环次数与预设循环次数比较。
[0181]
需要说明的是，需要预设判定的条件如下：
①
设置最大迭代次数t
max
，也就是预设循环次数，当迭代次数，也就是当前循环次数，达到t
max
时结束循环运算。
②
设置相场的微分阈值errphi，也就是预设变化量，当整个流体系统中相场的变化值δφ，也就是当前相场变化量，均小于微分阈值errphi时，判定流体系统达到稳定状态，结束循环运算。
[0182]
需要说明的是，当当前相场变化量大于预设变化量的阈值且当前循环次数小于预设循环次数时，计算将进入下一个循环状态，时间进入下一个时间步，计算从步骤s120重新开始，则当前物理场分布在进入下一个循环后，对于新循环而言，当前物理场分布将变为新循环的上一个时间步的物理场分布，称之为此前的物理场分布，以此进行循环运算。
[0183]
需要说明的是，当前相场变化量和当前循环次数每一次循环过程都会重新计算一次。
[0184]
可以理解的是，在执行循环过程直至满足预设条件之前，在每个循环过程中，还包括：
[0185]
当格点为边界点，对密度分布函数中的速度参数进行反弹算法处理，以修正边界点的此前密度场分布。
[0186]
可以理解的是，在执行循环过程直至满足预设条件之前，在每个循环过程中，还包括：
[0187]
当格点为边界点，对相分布函数中的速度参数进行反弹算法处理，以修正边界点
的此前相场分布。
[0188]
需要说明的是，流体场中流体的运动，以及施加电压时所产生的电场，需要在二维网络模型中定义边界条件：流体场和电场整体采用无渗透无滑移边界条件。在液体和固体边界，对流体的密度分布函数fi(x,t)和相分布函数gi(x,t)做反弹算法处理，首要改变其中的速度ei，具体反弹方法如下：
[0189]
与边界碰撞前
→
边界碰撞后
[0190]
e0
→
e0；
[0191]
e1
→
e3；e2
→
e4；e3
→
e1；e4
→
e2；
[0192]
e5
→
e7；e6
→
e8；e7
→
e5；e8
→
e6；
[0193]
界面处(x＝xb)的润湿的相场及速度边界条件如下：
[0194][0195]
u(xb)＝0。
[0196]
其中为界面单位法向量。
[0197]
流体边界处电势v(xb)＝vb,vb为边界处电势，电极正极边界处vb＝v
app
,电极负极边界处vb＝0，流体宏观速度u(xb)＝0。
[0198]
步骤s127，将当前物理场分布确定为稳态分布状态或目标状态；其中，预设条件为当前相场变化量小于预设变化量或当前循环次数等于预设循环次数。
[0199]
需要说明的是，当相场变化量小于预设变化量的阈值时将当前物理场分布确定为稳态分布状态，或，当前循环次数等于预设循环次数将当前物理场分布确定为目标状态。
[0200]
需要说明的是，当前物理场分布包含模拟中所有物理量的分布，包括：
[0201]
相场、密度场、压强、电场、速度场、第一物理模拟量、第二物理模拟量。
[0202]
可以理解的是，循环计算结束后需要将数据进行保存，并输出结果。
[0203]
根据本技术实施例的流体模拟方法，至少具有如下有益效果：首先，获取模拟数据集；其中，模拟数据集包括流体的多个第一模拟物理量和待电湿润的器件的多个第二模拟物理量；根据第一物理模拟量和第二物理模拟量定义器件的初始物理场分布；其中，物理场为二维网络模型，二维网格模型中的每一格点对应不同的坐标；之后，对模拟数据集执行循环过程直至满足预设条件，并将满足预设条件时的当前物理场分布确定为稳态分布状态或目标状态，循环过程包括：根据第一模拟物理量、第二模拟物理量、此前物理场分布中的此前相场分布计算的相场、此前物理场分布中的此前密度场分布计算的密度场，计算得到化学势和自由能密度函数；根据化学势、自由能密度函数、相场以及密度场，计算得到流体在每一格点的第一压强；根据二维网格模型上三相接触线区域对应的相场、此前物理场分布中的此前电势分布，修正三相接触区域对应的第一电场强度，以确定每一格点的第二电场强度；将第一压强、第二电场强度进行耦合，得到每一格点的第二压强；根据第二压强，计算得到流体的当前密度平衡分布函数、当前相平衡分布函数；根据当前相平衡分布函数、当前密度平衡分布函数，计算得到每一格点上流体的当前物理场分布；其中，当前物理场分布包括当前相场分布、当前密度场分布、当前电势分布；将当前相场变化量与预设变化量比较或将当前循环次数与预设循环次数比较；其中，预设条件为当前相场变化量小于预设变化量
或当前循环次数等于预设循环次数。本技术的流体模拟方法，通过将第二电场强度和第一压强进行耦合，得到第二压强，然后就可以通过第二压强计算得到密度分布；一方面，较大程度的简化了计算过程，通过计算量的简化，进而降低了对计算机的高性能要求；另一方面，操作者只需对相应的实际物理量进行设置就可以进行仿真模拟，较大程度的降低了操作者的要求，节省了学习的时间成本；最后一方面，通过修正电场强度后再对电场与流体场进行耦合，模拟仿真得到的参数更加准确，使得模拟流体在器件内的电润湿现象更加接近实际。因此，本技术的流体模拟方法，能够更加快速而准确地模拟流体的电润湿现象，模拟效果更加接近实际。。
[0204]
需要说明的是，本技术的流体模拟方法所得到的湿润角θ0的仿真值和young-lippmann方程的理论值如图3所示，采用本技术的流体模拟方法，所得到的湿润角θ0的仿真值更加接近实际的理论值，进而可使得在电润湿器件的设计开发过程中，一方面，能够更好的观察电润湿器件设计出来的参数，在发生电润湿现象时的变化过程；另一方面，通关观察变化过程以及电润湿器件施加电压时得到的稳态以得到模拟数据，能够根据模拟数据进行设计参考，例如，对参数进行优化，为电润湿器件的工程设计提供更加精确的参考依据。
[0205]
下面参照图5描述根据本技术实施例的流体模拟系统。
[0206]
如图5所示，本技术实施例的流体模拟系统，可以是任意类型的智能终端，如手机、平板电脑、个人计算机等。
[0207]
具体地，流体模拟系统，包括：
[0208]
至少一个存储器200；
[0209]
至少一个处理器300；
[0210]
至少一个程序；
[0211]
程序被存储在存储器200中，处理器300执行至少一个程序以实现上述的流体模拟方法。图5以一个处理器300为例。
[0212]
处理器300和存储器200可以通过总线或其他方式连接，图5以通过总线连接为例。
[0213]
存储器200作为一种非暂态计算机可读存储介质，可用于存储非暂态软件程序、非暂态性计算机可执行程序以及信号，如本技术实施例中的流体模拟系统对应的程序指令/信号。处理器300通过运行存储在存储器200中的非暂态软件程序、指令以及信号，从而执行各种功能应用以及数据处理，即实现上述方法实施例的流体模拟方法。
[0214]
存储器200可以包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需要的应用程序；存储数据区可存储上述流体模拟方法的相关数据等。此外，存储器200可以包括高速随机存取存储器，还可以包括非暂态存储器，例如至少一个磁盘存储器件、闪存器件、或其他非暂态固态存储器件。在一些实施例中，存储器200可选包括相对于处理器300远程设置的存储器，这些远程存储器可以通过网络连接至该流体模拟系统。上述网络的实例包括但不限于互联网、企业内部网、局域网、移动通信网及其组合。
[0215]
一个或者多个信号存储在存储器200中，当被一个或者多个处理器300执行时，执行上述任意方法实施例中的流体模拟方法。例如，执行以上描述的图2中的方法步骤s100至s127。
[0216]
下面参照图5描述根据本技术实施例的计算机可读存储介质。
[0217]
如图5所示，计算机可读存储介质存储有计算机可执行指令，该计算机可执行指令
被一个或多个处理器300执行，例如，被图5中的一个处理器300执行，可使得上述一个或多个处理器300执行上述方法实施例中的流体模拟方法。例如，执行以上描述的图2中的方法步骤s100至s127。
[0218]
以上所描述的装置实施例仅仅是示意性的，其中作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
[0219]
通过以上的实施方式的描述，本领域普通技术人员可以理解，上文中所公开方法中的全部或某些步骤、系统可以被实施为软件、固件、硬件及其适当的组合。某些物理组件或所有物理组件可以被实施为由处理器，如中央处理器、数字信号处理器或微处理器执行的软件，或者被实施为硬件，或者被实施为集成电路，如专用集成电路。这样的软件可以分布在计算机可读介质上，计算机可读介质可以包括计算机存储介质和通信介质。如本领域普通技术人员公知的，术语计算机存储介质包括在用于存储信息的任何方法或技术中实施的易失性和非易失性、可移除和不可移除介质。计算机存储介质包括但不限于ram、rom、eeprom、闪存或其他存储器技术、cd-rom、数字多功能盘或其他光盘存储、磁盒、磁带、磁盘存储或其他磁存储装置、或者可以用于存储期望的信息并且可以被计算机访问的任何其他的介质。此外，本领域普通技术人员公知的是，通信介质通常包含计算机可读信号、数据结构、程序模块或者诸如载波或其他传输机制之类的调制数据信号中的其他数据，并且可包括任何信息递送介质。
[0220]
上面结合附图对本技术实施例作了详细说明，但是本技术不限于上述实施例，在所属技术领域普通技术人员所具备的知识范围内，还可以在不脱离本技术宗旨的前提下作出各种变化。此外，在不冲突的情况下，本技术的实施例及实施例中的特征可以相互组合。