英语原文共 16 页，剩余内容已隐藏，支付完成后下载完整资料
在OpenFoam中使用对数重构对三维粘弹性方腔顶盖驱动流的数值模拟
Florian Habla, Ming Wei Tan, Johannes Haszlig;lberger, Olaf Hinrichsen
摘要
在这项工作中，我们实现了由Fattal和Kupferman（2004）在基于配置有限体积法（FVM）的开源CFD软件OpenFOAM中提出的粘弹性本构方程的对数重构。这个实现包括一个有效的特征值和特征向量子程序，并且当它与CUBISTA方案等适当的对流方案结合使用时，算法最终在时间和空间上都是二阶精度的。新开发的求解器首先用一个粘弹性流体的启动Poiseuille流的解析解来验证，然后应用于盖驱动的空腔流的三维和瞬态模拟，其中粘弹性流体由Oldroyd- B本构方程。为了检查结果的网格收敛性和四面体网格，显示了一阶迎风格式和CUBISTA格式在四个不同大小的六面体网格上的结果，以显示对非结构化网格的适用性。对于不同的Weissenberg数值得到的结果在主涡流中心的位置，流线模式，速度和应力分布等方面进行了介绍和讨论。我们能够获得Weissenberg数值的足够的网格收敛结果，如果不使用对数构型重构就不可能获得这些结果。 Weissenberg数的稳定性上限没有找到。
1.介绍
对于复杂的粘弹性流体进行高Weissenberg数的数值模拟是一个突出的挑战。幸运的是，在过去几年时间里，数值算法在稳定性和准确性上取得了重要的发展。在学术界用来测试数值算法的基准问题包括收缩流、在球体周围流、方腔流等等。
到目前为止，大部分关于方腔流的研究工作仅仅是理论方面的，这主要是因为这种非常复杂的流动模式在简单的几何图形上得到了发展。牛顿流体在方腔内的数值模拟是直截了当的，关于这方面的文献有和多，例如：sheu和tsai研究了稳定流拓扑结构在具有有限元的三维方腔顶盖驱动中雷诺数Re=400.但是相反，由于预测粘性流体在方腔中流动的复杂性，关于这方面的文献迄今为止只有少数。而对于这方面有这点兴趣部分原因是低Weissenberg数相对来说还是比较容易实现的。比如说：Demir研究了由Crinimale-Ericson-Filbey（CEF）方程控制的粘弹性流体的瞬态流动，对于所有的雷诺数，得到最大的Weissnberg数We max= 0.01。与其类似的，Demir在移动顶盖上施加均匀的速度，这导致了两个上角处的不连续性，这也将可达到的最大的Weissenberg数限制在0.2以下。最近有许多试验加入了正规化的边界条件，去掉两个顶角上的速度和速度梯度，比如Fattal 和 Kupferman的研究。然而，靠近盖子的薄边界薄层和结构上的奇点在顶角的下端处张量还是存在，这也使得这个问题愈发困难，并构成了Weissenberg数的上限。
Fattal和Kupferman提出了所谓的对数化重构（LCR），其中构造张量的对数化演化方程被求解而不是去求解方程的本身。这消除了停滞点处应力（以及结构张量）的指数变化。新的变量（结构张量的对数）可以比结构张量的指数更好的通过基于多项式的插值计算。这种技术的好处就是可以保持结构张量的正确性始终不变。Fattal和Kupferman在FENE-CR和 Oldroyd-B这两种二维粘弹性液体顶盖流动中测试了这个技术，他们表示稳定的模拟可以使得Weissenberg数达到5.动能的振荡表明，当Weissenberg数达到3以上时，收敛时的损失和瞬态流动模式开始发生。Pan和Hao将对数重构技术应用到其有限元代码中，并模拟了一个空腔内粘弹性流体的二维斯托克斯流，达到了3的Weissnberg数。类似于在FENE-CR和Oldroyd-B上做的研究，他们解决了粗糙网格上的构象张量的对数化演化方程。这减少了高频模式，进一步使得对数化方程的解稳定。在这项工作中，一阶迎风差分格式被用于离散本构方程的对流项。然而，当他引入大量的数值扩散时，迎风差分被认为是最不准确的，尽管这有助于稳定解决方案以及显著增加了最大可能实现的Weissenberg数。Oliveira使用有限体积法模拟了具有正则化边界条件的FENE-CR流体的稳态流动和瞬态反冲，其相对较大的延迟率为0.79。在蠕动流条件（Re=0）下，最大的Weissenberg值为10。将本构方程的平流项离散为收敛性和高度准确的CUBISTA方案，其形式上为3阶。Yapici等人使用了有限体积的方法代码，也利用逆风方案在本构方程中处理对流项。他们对在不同雷诺数下的空腔内的一种Oldroyd-B液体的流动进行了模拟，并给出了一种关于蠕变流动情况的Weissenberg数。
在这一工作中，在基于有限体积的开源软件OpenFOAM中实现了对数-构象重组。结果展示了Oldroyd-B流体在启动Poiseuille流动和在蠕变流动条件下的三维空腔的流动，剩下的工作安排如下：在第二章中给出了控制方程，并解释了对数构象方法的理论，在接下来的第三章中，我们描述了OpenFOAM的数值实现，第四章给出了Poiseuille的启动结果，并讨论了三维驱动腔，最后，在第五章中做最后的总结。
2.数学背景
2.1.控制方程
我们考虑由Oldroyd-B本构方程控制的不可压缩和等温粘弹性流体的流动。平衡方程是动量平衡方程和质量平衡方程：
其中U是速度，q是密度，t是时间，p是压力，s是额外的应力张量，可以写成溶剂和聚合物消耗的总和。
（3）
牛顿法则适用于溶剂的消耗:
（4）
eta;是溶剂黏度，tau;是聚合物消耗量，Olroyd-B方程可能会出现，但在这里应该指出许多其他的本构方程，如Giesekus，SPTT或FENE型模型以及其他可能会出现在这里的模式模型。Olroyd-B方程的定义如下：
（5）
lambda;是松弛时间，eta;是溶剂黏度，表示上对流时间的导数。
（6）
延迟率beta;被定义为溶剂黏度和总黏度之比
（7）
Oldroyd-B方程可以用构象张量c来重写
（8）
其中I是单位矩阵，使用方程（8）使得本构方程（5）变为
（9）
聚合物压力的计算不是通过方程（9），而是利用方程（8）得到的。
2.2.对数—构象方法
要使得演化方程（9）具有正确性，构象张量c必须是正定的。在高弹性的流动中，这个属性可能被违反，这往往导致数值计算失败。主要问题表现为一维问题：在高变形率区域，拉伸和松弛项呈现指数增长。对流可以平衡这个增长，然而，由于对流项基于多项式插值，所以对流项不能平衡指数放大，然后导致数值模拟失败。为了应对这种不稳定性，Fattal和Kupferman提出了对数变换的方程式（9），后者被称为“记录 - 构象方法”，将在下文中概述。
由于构造张量c是一个对称正定的矩阵，因此可以根据它对角化：
Lambda;是一个对角矩阵，由c的三个特征值组成，R是一个正交矩阵，由c的三个特征向量组成。任何对角矩阵都可以通过对对角线上的对数元素来进行对数运算。构造张量c的对数被定义为：
因此，现在是我们的新变量。对数构象重组的主要特征是速度梯度的分解。
其中Omega;和N为反对称矩阵，B是对称无痕矩阵，与c交换。X代表旋转，B代表纯扩展，公式（12）的最后一项作为“虚拟部分”。 引入这个分解到等式（9），很容易获得：
为了获得矩阵X和B，速度梯度首先旋转到构象张量c的主轴上:
对于我们随后定义的三维流动问题
矩阵Omega;和B现在可以通过以下公式获得：
其中，在的情况下，可能会导致未定义的零除现象。幸运的是，当c = I时，这只是第一次的情况。在那个情况下，我们定义Omega;=0，和。
最后，对于我们新的变量，我们可以重写公式（13）：
解出这个方程后，可以应用逆变换。新获得的按照以下公式分解：
构象张量可以通过
计算。
这里强调构象张量Lambda;的特征值与工作变量Lambda;w的特征值之间的关系是
而特征向量R对于c和W来说都是相同的。这对于实现一个有效的数值算法是相关的，这将在下一章中描述。
对数构象的一个主要优点是构象张量的严格正定性。 即使在离散水平上也能保证这一点。 因此，高Weissenberg数不稳定将不予考虑。 正定性可以通过检查行列式c 大于等于0来验证。
3. 数值方法
在本节中，我们试图描述在这项工作中开发的数值方法。 上面描述的模型在OpenFoam中实现，发达的求解器可以被认为是Favero等人最初的工作的进一步发展和我们接下来的论文。
3.1. 离散化方案
本构方程中的对流项将由逆风方案或CUBISTA方案离散化。 逆风方案是一阶精确的，虽然它是非常稳定的，但会引入过量的数值扩散，对于大多数复杂的问题来说，它是不够的，通常应该避免。 由Alves等人开发的CUBISTA对流差分方案是专门为粘弹性本构方程的对流项开发的高精度方案。该方案基于三阶QUICK方案，在均匀网格和平滑流程上形式上为三阶。 可以通过
在Sweby图中描述。
其中Phi;是通量限制器，r是内插变量Phi;的连续梯度的比率。 Jasak等人引入了一个修改，以避免需要确定在最远的逆风单元中的值，而这个值是使得非结构化网格不切实际。这个r是从以下公式计算出的:
其中d是节点C和顺风节点D之间的距离矢量。将面值Phi;f与Phi;(r)混合，而Phi;(r)通过迎风差分和由中心差分得到。
通过这种方式，对流项完全隐含地离散化，因为它只有助于矩阵，而不利于源矢量。 这种实现与Alves等人使用的延迟修正方法相反。其中与迎风方案对应的部分是隐含离散的，而其余的高阶部分是离散化的，从而保证了对角线的优势。 延期修正法中的CUBISTA方案最近被Afonso等人应用于配置的有限体积框架内的对数构型对流项。
动量平衡中的对流项在所有的模拟中都是可以忽略的。尽管如此，这个术语仍然保留，用与本构方程中使用的方案相同的方案进行离散化。
出现在动量和对数化本构方程中的时间导数用广义三点方案离散化，也称“齿轮”方案。 这个方案在时间上是二阶准确的，然而，由于需要在步骤n和n-1存储变量，它增加了存储需求。 Xue等人对此方案进行了全面的评估， 对于粘弹性流体的瞬态Poiseuille流。 广义三点方案被发现具有最小的相对误差，并允许在他们的模拟中最大的时间步。然而，在第一个时间步，只有初始条件已知，并且变量n-1还不存在。 因此，隐式欧拉方案被用于第一时间步骤。动量平衡中的扩散项通过二阶精确线性插值和高斯积分来离散化。 在非正交网格的情况下，为了保持二阶精度，正如Jasak所解释的那样使用校正方法。 应用源项的二阶高斯公式（所有项，但方程（17）中的惯性项和对流项）。速度梯度必须在进行分解之前先进行计算。（12）也是使用基于线性插值的二阶精确高斯公式来计算的。
3.2. 动量方程离散化
将式（3）的应力分解与牛顿定律式（4）一起引入动量方程，并将时间导数结果离散化：
在蠕变流动条件（Re趋向于0）和可忽略的溶剂消耗（eta;s趋向于0）的情况下，方程（24）变得奇异。 为了确保方程（24）是可解的并允许速度场的更新，可以引入与eta;p成比例的附加椭圆项：
这被称为双侧扩散（BSD）技术，即使在非蠕变流动条件下和非零溶剂粘度下也常常用于稳定模拟。 然而，BSD并不是短暂问题的适当选择，因为它引入了大量的瞬态扩散，并可能导致虚假响应。 如果l.hs.s上的扩散条件存在问题 和r.hs.s. 不要完全抵消。 在这种情况下，由此产生的误差就像扩散一样。 这可能是因为不同时间步长的速度被用于两个项（U n 1在lhs上的项和U n在rhs上的项），或者尽管对于两个项使用U n 1如公式(25)--在一个时间步骤内不会收敛。 但是，如果达到收敛，则两项完全抵消，避免了BSD的负面影响。 第3.5节进一步讨论。
为了求解动量平衡，我们需要给tau;p加上边界条件。 虽然没有求解PDE（由方程（8）计算），但我们需要施加边界条件，以便能够使用高斯定理计算tau;p的发散。 因此，在坚实的墙壁上，我们将分配一个二阶外推边界条件，如Habla等人所描述的一个零梯度边界条件。 一般使用二阶外推法，因为它更准确（但不太稳定）。 如果稳定性是一个问题，我们将外推的顺序降低到一阶，这是通过忽略与P的第一个空间导数成比例的项来完成的。最终，这个一阶外推等于零梯度边界条件。
3.3. 通量公式和压力方程
l.h.s.上的术语 等式（24）或（25）的隐含条件和r.h.s.明确。压力校正方程是通过将这些方程式转换成半离散化的形式而得到的，其中所有项除了压力梯度以外都是离散的。该时间指示器n 1在下面为了可读性而被放弃。我们可以写出来：
术语包括非对角系数和(显式)源项。除以对角线系数的速度可以写成：
在
全文共23228字，剩余内容已隐藏，支付完成后下载完整资料
资料编号：[11768]，资料为PDF文档或Word文档，PDF文档可免费转换为Word
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
