英语原文共 33 页

使用OpenFOAM电驱动流动的数据模拟

Francisco Pimenta^a和Manuel A. Alves^b

CEFT，波尔图大学工程学院化学工程系，Rua Dr. Roberto Frias，4200-465波尔图，葡萄牙

在这项工作中，我们研究了在有限体积框架中电驱动流（EDF）模型中安装使用OpenFOAM。泊松能斯特-普朗克模型用于带电物质的传输，并且它与纳维-斯托克斯方程耦合，控制流体流动。此外，还执行了泊松-玻尔兹曼和德拜-休克尔模型。仔细处理离散化方法和边界条件，以确保通用网格中离子种类的保守性和空间和时间的二阶精度。所开发的求解器的适用性通过两个相关的EDF中说明：导电圆柱体周围的感应电荷电渗和跨离子选择性膜的电荷传输，在当前工作中开发的作为开源免费提供的求解器，可以成为研究通用电驱动流的有价值且通用的工具。

一、导言

微米和纳米尺度上的流体和电场之间的相互作用是几种现象的起源，并且越来越受关注。这些互动经常发生分为两类，电流体动力学（EHD）和电动力学（EK）^1,2我们在整个工作中统称为电驱动流（EDF）。 而电流体动力学其特征在于界面处净电荷的积累，区别电动力学的特征是形成中性双电层（EDLs）¹。 电动力学和电流体动力学在燃料电池³中发挥重要作用，电解液⁴，分离技术^5,6，液体在微流体系统中泵送^7,8和混合^8,9，以及其他应用。

通过理论模型理解电动力学和电流体动力学现象已经成为过去几十年来所追求的挑战。然而，仍有一些问题在最基本的理论^1,10中和其他物理领域一样仍然没有答案。这是由于理论在攻克所涉及的所有事件方面的失败，由于简化或仅仅因为它们是未知的，以及很难从可用模型中提取解决方案，关于这个最后一个问题，封闭式解决方案通常难以获得常在电驱动流中出现的即使是简单的几何形状的方程组。匹配渐近线^11,12和数值的方法^13-17是在解决由此产生的方程组经常采用的方法，其中后者特别适合模拟复杂几何形状的电驱动流。

有限差异^13,18-25，有限元^13,26-31，有限体积^{16,32,33,34-43}和玻尔兹曼网格^15,44-46是已经成功应用模拟电驱动流的数据方法。在可用的软件包中，基于有限元方法的一直是很受欢迎的选择。 在开源软件包的范围内，一些使用可以在文献中找到的Gerris^32,33和OpenFOAM^34-39,47-49，但这些软件包中没有一个能够始终如一地为EDF提供各种即用型模型，尽管两者都允许用户构建自己的模型。在这项工作中，我们使用OpenFOAM安装包由于以下特点：缺乏付费许可，处理通用多面体网格，并行计算能力，创建/改变源代码的可能性，模块之间易于集成（流变模型，多相模型等）并得到学术界和工业界的广泛认可。

关于使用OpenFOAM模拟电驱动流，Nandi和Aluru³⁷将泊松能斯特-普朗克方程与纳维-斯托克斯方程耦合以研究在微纳米通道中的I-V曲线。然后推导出面积平均的多离子传输模型来研究类似的问题^38,39。与斯托克斯方程耦合的泊松能斯特-普朗克模型也已用于研究稳定条件^47-49下纳米孔内的流动。 Zografos等人³⁶使用德拜-休克尔模型优化了收缩-扩展微通道以产生均匀的拉伸电渗流。在两相电驱动流的范围内，利马和阿维拉³⁴实施了泄漏电介质模型，以研究在电场作用下粘弹性液滴的变形。Roghair等人³⁵使用类似的模型来分析牛顿流体的电润湿，其中流体和固体区域以耦合方式进行数值模拟。据我们所知，只有Roghair等人³⁵的工作中使用的求解器已经提供给公众，然而，这只能实现牛顿流体的两相电驱动流的泄漏介电模型。

这项工作报告了模型的数值实现，使用OpenFOAM工具箱模拟牛顿流体的单相电驱动流。这些模型在rheoTool的顶部实现，rheoTool是一个开源工具箱，用于使用OpenFOAM模拟广义牛顿/粘弹性流体流，现在能够模拟压力和电驱动流（单独或混合）。我们使用泊松能斯特-普朗克，泊松-玻尔兹曼和德拜-休克尔模型将这项工作局限于动电问题的分析，但其他模型（滑动模型，欧姆模型等）和其他案例可以在rheoTool中找到。在两个特定电驱动流的模拟中进一步证明了求解器的适用性：围绕导电圆柱体的感应电荷电渗（ICEO）和穿过离子选择性膜的电荷传输。这些情况也用于评估数值算法的准确性和稳定性，并且所获得的结果可用于基准目的。总的来说，这项工作是为了增加通用开源解算器的可用性，内置电驱动流模型，用于一般流变学的流体。

本文的其余部分安排如下：第二部分介绍了控制方程，第三部分描述了在OpenFOAM有限体积框架中数值实现的细节。 第四节讨论了应用案例中数值方法的性能。 最后，第五部分介绍了未来工作的主要结论和观点。

二、管理方程式

在电驱动流的模拟中，可以识别两个不同但耦合的分量：流体动力学分量，由连续性和动量方程表示; 电子元件，通常由泊松能斯特-普朗克方程表示，用于在稀释电解质中传输带电物质。这两个组件将在下一节中介绍。

A.流体动力学：连续性和动量方程

考虑在电力的作用下牛顿流体的瞬态，不可压缩，等温，单相，层流。连续性（方程1）和动量（方程2）方程式

其中u是速度矢量，t是时间，p是压力，f E代表电气每单位体积的力，rho;是流体密度，eta;是粘度。

通过条件f E确保流体动力学与电力之间的耦合。识别磁效应，该条件可以从静电麦克斯韦应力张量得出。

其中ε是流体的电容率，E是电场，rho;E是电荷密度。 在该工作中解决的所有问题中，流体的电容率被认为是恒定的，在那种情况下fE=rho;EE，虽然电场通常由电势的负梯度定义，但rho;E的计算依赖于模型，如下所示。

B.泊松 - 能斯特 - 普朗克方程

在电场作用下，低离子强度电解质中带电物质（离子）的传输可用能斯特 - 普朗克方程描述

其中c i是摩尔浓度，D i是扩散系数，mu;i是电迁移率，Psi;是电势，在静电中具有E=-▽Psi;。指数i代表电解质中的每个单独物质。方程（4）左边的项代表ci的材料导数，右边的第一项是由于扩散引起的ci通量，最后一项代表ci通过动作的运输电场，也称为电迁移。这种电迁移术语形式上类似于对流术语，其中可以识别电迁移速度（或离子化后的通量），uMI=mu;i▽Psi;。在这项工作中考虑的所有电解质都遵循能斯特-爱因斯坦关系，mu;i=Di ezi/kT，其中zi是电荷效率，e是基本电荷，k是玻尔兹曼常数，T是绝对温度。但是，为了通用性，我们将保留通用u的符号，因为奈斯特-爱因斯坦关系不适用于几种聚电解质，例如DNA分子⁵⁰。给定域中的电位分布可以用高斯法（忽略极化）计算。

电荷密度定义为

其中F代表法拉第常数，m是不同离子种类的数量。方程（1）-（6）形成了用于模拟大多数EDF的基本理论。然而，这组偏微分方程的电分量经常被简化或变换以避免数值问题，例如，由可以共存的显着不同的长度和时间尺度产生的数值刚度，或者仅仅是为了能够推导出封闭形式的解。接下来介绍其中两个简化模型。此后，由等式（4）-（6）定义的Poisson-Nernst-Planck模型缩写为PNP。

C. 泊松 - 玻尔兹曼模型

泊松 - 玻尔兹曼模型是广泛采用的PNP模型的简化，并且基于离子物质遵循玻尔兹曼分布的假设^44,51,52

其中c i，0是i的参考浓度，其中Psi;=Psi;0, Psi;代表固有电位，如后面3.1节所述。方程（7）也可以看出，作为将能斯特-普朗克方程（方程4）与参考点（其中c i = c i，Psi;=Psi;0）相结合的结果，假设零材料导数（通过对流的稳态和可忽略的离子传输）。变量Psi;已被Psi;取代，以便在外部施加的电势与固有电位共存的一般情况下简化Psi;0的定义。与通用PNP模型相比，该简化模型具有有限的适用性，主要是由于假设玻尔兹曼平衡到固定参考点。关于这些限制的更多细节可以在其他地方⁵²找到，但粗略地说，泊松-玻尔兹曼模型基本上对于具有可忽略的电荷传输，非重叠中性双电层和低内在电位的稳态问题是有效的。

我们将假设A，这相当于考虑参考点中的电中性，通常是批量解。除了前面提到的限制之外，该假设不会产生进一步的限制，并且仅用于从等式中移除变量A. 当在方程（6）中插入方程（7）时，获得所谓的泊松 - 玻尔兹曼模型的电荷密度。

方程（5）和（8）定义泊松-玻尔兹曼模型，此后缩写为PB，其中电气元件现在与流体动力学相分离（但不是相对）。 另外，关于系统的电气元件的唯一未知的是电势，其由等式（5）和（8）计算。

D. 德拜-休克尔模型

在低电势极限下，PB模型可以进一步简化为所谓的德拜-休克尔模型。在此近似下，电荷密度可以通过将泰勒级数中的方程（8）的指数项扩展到第二期。

方程（5）和（9）定义了Debye-Huuml;ckel模型，此后缩写为DH。同样，电气元件也与流体动力学分离。

三、数值方法

有限体积框架中的控制方程的数值离散遵循OpenFOAM中可用的标准程序，其他地方描述了相关的偏微分方程，例如Moukalled等人⁵³和Pimenta和Alves⁵⁴。 该方法适用于并置网格的一致运算符。 为简明起见，此处未介绍标准运算符的离散化过程（时间导数，对流项等）。 在本节中，我们将分析限制在EDF的特殊性。

A.将电势分配到外部和内部电位。

当用简化模型（PB和DH）解决EDF问题时，通常可以方便地分解两个变量中的电位，Psi;=psi; phi;Ext，其中phi;Ext是外部施加的电势，psi;是固有地存在于通道中的电势，通常与EDL形成⁵¹有关。 在这种方法下，高斯定律也被分解为两个方程。

独立地和顺序地解决每个变量。遵循这种策略，通常的做法是计算时仅考虑phi;Ext在动量方程（方程3）的电力项的电场中的作用，因此fE=-rho;Ephi;Ext，正如Kyoungjin等人⁵¹所解释的。这主要是为了避免在墙壁附近产生大的压力梯度，这是由于这些区域中Psi;的正常梯度-这种贡献可以通过考虑它由相反的压力梯度平衡而在动量方程中抵消，假设它不会影响流量⁵¹。 在粘弹性流体的情况下，这种转变可能是特别重要的（有利的），其中大的不平衡压力梯度的发展可能人为地触发流动不稳定性。

B. EDF中的壁面边界条件

在微米级和纳米级的EDF中使用的几种材料是不可渗透的和电绝缘的，例如玻璃和聚二甲基硅氧烷（PDMS）。 与电解质接触时，这些材料在固/液界面处形成双电层。 考虑到PNP模型，通常分配给这些表面的边界条件是零速度，离子的无通量条件和指定的电位或指定的表面电荷（表面上的常数或可变）。

墙上的无通量条件导致（假设墙壁没有穿透）

其中S f是壁法向矢量，其大小代表墙壁上的面f的面积。 方程（11）形式上为c i的Robin型边界条件。 为了便于实现，有几种方法可以将方程（11）转换为Neumann或Dirichlet边界条件，其中三种方法如下

用ci，P表示在拥有计算ci，f的边界面的单元中心处计算的ci的值（下标P和f分别代表单元和面），n=Sf/|Sf|是单位矢量，正常面对f。 方法（I）直接从方程（11）得到，只需在方程的不同侧隔离每个项。在离散边界上的梯度隔离所有c i，f项之后，可从方法（I）获得方法（II）。 方法（III）基于细胞和面部之间的方程（11）的分析积分导出c i，f的表达式，其自然地导致边界附近的c i的指数变化。当使用ci，f的解析表达式计算ci的梯度时，保留这种指数行为，而不是取▽ci|f*n=(ci,f-ci,p)/ delta;，delta;是中心和面f之间的距离。这最后一种方法通常将指数变化线性化，导致粗网格中的误差较大。

上述三种非常相似的方法在强加F i,f = 0时很有效。实际上，它们都在离散水平上满足方程（11），其表明扩散和电迁移通量在边界处相互抵消。这可以通过忽略扩散和电迁移通量到源向量的作用（在无通量边界处）从而在矩阵构建阶段完成，避免在壁面上计算ci或其梯度的需要。实际上，上述三种方法只能根据当该变量用于后处理目的时检索的ci值的准确性来区分，或者当它用于某些高阶方法时，可以提高非正交网格的精度。此外，当使用高分辨率方案用于方程（4）的对流项时，还需要c i,f的值来计算限制器。基于这些标准，在非正交网格中，方法（III）被认为更准确并且更少依赖网格。因此，在整个这项工作中使用方法（III）。应当注意，当使用该方法来评估c i,f时，考虑来自前一时间步长/迭代的c i,P的值，导致明确的边界条件（隐式实现也是可能的）。

关于壁上的电位，如上所述，两种常见的选择是施加电势，或者定义表面电荷密度（sigma;），这可能与电势有关，▽Psi;|f*n=sigma;/ε。对于PNP模型，在壁上施加电势是不太自然的边界条件，并且如果使用单个电势变量，则实施在不规则形状的壁上也更复杂。这个问题的一个常见解决方法是分解PNP模型中的电位。

当使用简化的PB和DH模型时，只需要两个电位的边界条件，因为它们是唯一被解决的与电相关的变量。考虑到前一节讨论的电位分解，外部施加的电场与绝缘壁相切，而固有电位可能与zeta;电位（xi;）或

资料编号：[3687]