一种新方法来模拟粘弹性流体自由表面流动

文摘

自由表面流出现在各种各样的工程应用。预测的动态特性问题,具体数值方法需要准确地捕获自由表面的形状。提出了一种新的方法,它结合了场(ALE)技术与有限体积法(有限体积法)来模拟含粘弹性自由表面流动。基于一个开源CFD叫做OpenFOAM工具箱,我们设计了一个ALE-FVM自由表面仿真平台。同时,膨胀流动一直与我们的调查提出平台自由表面现象的进一步分析。结果验证了该方法的正确性和有效性,自由表面模拟牛顿流体和粘弹性流体。

1。介绍

近年来,越来越多的兴趣自由表面模拟(如膨胀和泡沫增长)。膨胀,作为一种典型的粘弹性流体现象,被广泛研究。数值模拟挤压需要解决自由表面流动问题。数值方案免费表面处理可以分为固定网格(接口捕捉)和移动网格(界面跟踪)方法。固定网格方法捕捉接口在欧拉描述,包括流体的体积(受到)方法(1),水平集方法(2),标记和细胞(MAC)方法(3),通过这些方法中的数值方案界面重构。虽然这些方法可以解决尖锐的接口和大变形问题,自由表面的位置通常是不够精确的。相反,移动网格方法明确跟踪接口的位置的计算网格,保证接口的精确位置,但可能减少网格的质量在网在大变形情况下移动。场(ALE)方法(4移动网格方法之一。它结合了利用拉格朗日方法跟踪接口和欧拉方法求解控制方程的效率。因此,它已广泛应用于波自由表面模拟。

ALE方法通常是结合有限元法(FEM) (5),界面上的牵引边界条件可以在一个整洁的方式通过集成力动量方程的部分。然而,随着发展的有限体积法(有限体积法),有限体积法在流体模拟系统更具吸引力在当地保护好财产。因此,结合ALE法与有限体积法能够代表自由表面的位置直接计算网格,同时受益于良好的保护财产和有限体积法。Tuković和Jasak6)开发了一个平台来模拟牛顿流体自由表面问题通过使用ALE方法在有限体积法系统。为了得到准确的速度从表面上看,一个新的表面网格加上传统的网格,但它增加了计算成本和复杂性。此外,它不能解决与粘弹性流体。这样做最大的挑战是在自由表面上设置适当的边界条件。

本文提出了一种新的数值方法对二维粘弹性挤压流动模拟。最好的作者的知识,本文是第一个实现结合啤酒技术和有限体积法在研究时间粘弹性挤压模拟,它扩展了ALE方法的应用领域。本文的主要贡献可以总结如下:(我)我们解决的挑战与有限体积法耦合ALE方法模拟时间粘弹性挤压流动。设置适当的自由表面边界条件的障碍以及解决通过推导狄利克雷压力边界条件,引入平衡力对自由表面。该方法能够处理牛顿流体和粘弹性流体。(2)我们实现了一个耦合ALE-FVM仿真平台基于OpenFOAM(开源领域操作和操作)7]。小说全面验证了数值方法比较与先前的文献。相比Tuković和Jasak的工作(6),我们的方法大大降低了存储和计算成本放弃额外的表面网格。此外,基于OpenFOAM,我们的自由表面数值解算器可以很容易地扩展为大规模并行模拟大量的处理器核心。

剩下的纸是组织如下。节中给出了该方法的细节2。报告的数值结果与讨论部分3。结论部分4。

2。方法

在本节中,Oldroyd-B模型的数学公式和数值算法在本文将使用。

2.1。数学模型

等温和不可压缩粘弹性流体流动的控制方程包括质量守恒方程(连续性),动量方程和本构方程。这些Oldroyd-B模型方程的形式给出 在哪里是流体密度,是流体的速度,是细胞的速度,是压力,是应力张量的聚合物的贡献。和分别是溶剂粘度和聚合物粘度。特征松弛时间,是变形速率张量,它可以表示为 聚合物的upper-convected时间导数应力张量的定义是

为方便与他人的作品,以上变量通常转化为无量纲的形式。新的无量纲变量雷诺数和韦森伯数介绍了通过扩展以下变量 在这特征长度(的一半宽度死在这工作),特点是速度(平均速度在这个工作),然后呢的总和吗和;也就是说,总粘度。是粘度比参数。

2.2。数值算法

2.2.1。自由表面边界条件的实现

膨胀模拟中最困难的问题是自由表面边界条件的实现。边界条件对自由表面可以分为两类,运动学边界条件和动态边界条件。运动学边界条件要求流体速度的法向分量等于界面网格速度的法向分量,用 上面的方程意味着流体之间的相对速度和网格垂直于界面是零。在OpenFOAM,流体的速度是存储在单元中心,这意味着很难直接得到界面网格的速度。然而,流体的流量通过每个单元可以很容易获得。因此,界面的速度可以根据计算净表面质量流量。详细描述关于运动学边界条件的实现可以在[6,8]。接口界面节点移动后,内部节点需要根据移动解决拉普拉斯方程的解决方案: 在哪里扩散系数。解算器选择移动节点在这个工作是拉普拉斯脸分解(9),可以保持良好的网在大变形情况下的质量。选择和扩散系数成反比自由表面的距离的平方。再生内部网格以这种方式形成的平等的潜力和减少非正交性进而增加网格的质量。

动态边界条件也被称为牵引边界条件。忽略表面张力,这意味着牵引施加流体的气氛应该等于大气力量作用于流体。不失一般性,大气压力设置为。边界条件是给定的 在哪里是总应力张量,即柯西应力张量在自由表面,然后呢是额外的应力张量。在OpenFOAM求解方程时,庇索算法需要这个词和这个词单独存在于这个方程。然而,这个术语是一个组合,,和。因此,该边界条件很难满足直接。大家都知道这一点,在有限体积法系统中,每一项的方程应集成在细胞体积。因此,右侧(1)是解决。根据高斯定理,空间导数项可以转化为积分在细胞表面边界体积。然后,我们有 在哪里面对区域向量和吗是脸的单位法向量。和代表脸周围细胞和自由表面的脸,分别。

动态边界条件要求;然而,它就不会是这样了如果任何边界条件是不恰当的。注意的是边界条件只会影响细胞拥有边界的计算,我们介绍一种平衡力量在自由表面(它是零别处包括内部字段和其他边界字段),以补偿零边界力的影响。通过添加这个词到原始的动量方程,右侧(1)可以解决 可以看出,总应力的影响对自由表面以这种方式抵消,这意味着动力边界条件同样是满意。动量方程实际上是解决的形式: 在哪里表明这个词只是应用于自由表面边界。从简单的迭代算法解线性动量方程和连续性方程,需要适当的压力边界条件在求解泊松方程的步骤。根据力平衡在界面的法线方向,狄利克雷边界条件的压力 重要的是要注意,这里只是一个导出量的压力条件,只出现了,因为需要解决的泊松方程在隔离。牵引边界条件时真正的物理量,它必须满足流体。前提是和共同解决,不需要压力边界条件了(10]。

2.2.2。开发者的方法数值稳定

为了提高稳定性的线性动量方程,特别是对于高数字流,离散弹性分割应力(开发者)数值Matallah提出的策略等。11使用)。动量方程可以写成 在哪里和添加人工粘度比。这个词作为源项。一个典型的选择是。这种策略大大增加数值稳定,可以应用于模拟各种本构模型。

2.2.3。仿真过程

除了自由表面边界条件的处理,仿真过程类似于(6如下所示:(1)为一个新的时间步,速度场,压力场聚合物应力场,流体质量流量从以前的时间步。最初的净质量通量通过界面的面孔根据减去前面的界面位置计算体积的流体质量通量脸通量。(2)移动界面网格点,这样净质量流量通过接口可以补偿和质量守恒定律是满意。(3)根据界面网格点的位移,内部网格点的位置是通过求解拉普拉斯方程获得。(4)外层迭代循环的开始。(一)新的估计聚合物应力张量通过求解指定的本构方程与先前的速度场;(b)设置边界条件和估计的压力平衡的力量;(c)计算最新的速度场和压力简单的算法;(d)移动网格点的接口类似于步骤(2),直到达到规定的精度。(5)如果该项目没有到达最后的时间,返回步骤(1)。

再次执行步骤(4)(d)为了避免的操作步骤(3)中,这是一个耗时的步骤。这样的过程是合理的位移提供接口在一个时间步长小于最近的界面单元的厚度。这条件很容易满足减少的时间间隔。

2.2.4。其他通用数值方案用于仿真

除了上述方法外,还用于低松弛算法提高数值模拟的稳定性,松弛参数和分别是0.3和0.5。

向后欧拉方法和二阶时间导数。线性系统的解算器,PCG(预处理共轭梯度法)与AMG(代数Multi-Grid方法)预处理用于压力,而BiCGstab (Bi-Conjugate梯度稳定方法)与对角线不完整图片右下方(DILU)预处理用于速度和压力(12,13]。绝对的公差是的速度和分别对压力和压力。

3所示。结果与讨论

3.1。基本的几何形状和边界条件

素描造型的几何形状和边界条件的牛顿和Oldroyd-B液体如图1。的一半宽度死去选为特征长度。上游和下游的长度设置为和,分别。他们应该足够长的时间以确保下游流量入口的模具出口前和出口完全开发。摘要长度设置和。

压力,速度,和聚合物压力耦合的方式解决了在我们的系统中,我们通常给所有正确的边界条件或任何其中一个边界。为了避免不必要的失败数值在入口附近,充分发展流动的条件强加于入口边界。的速度在入口方向,而速度方向是零。因此,流入的平均速度和墙附近的剪切速率。聚合物在进口Oldroyd-B模型应力张量是根据分析的解决方案。无滑动边界条件设置在墙上。中心线的对称边界条件设置所有字段。对称边界条件由OpenFoam本身提供。这意味着有飞机通量对流和扩散通量;也就是说,垂直速度和剪切应力分量为零(和)。我们假设的流出口已经充分发展;因此,出口压力设置为零。注意,这里的压力只是一个参考压力用于解决控制方程。的自由表面,压力最初设置为固定值为零,和它将改变动态模拟(见(11))。没有特别声明,零梯度边界条件设置为其余边界。

3.2。验证牛顿流体的计算模型

通过设置,,,牛顿流体,Oldroyd-B模型降解。在等温条件下,忽略重力和表面张力的影响,肿胀率,定义为挤出物的最大直径除以死,的直径从我们的模拟。是一个很好的协议的结果Mitsoulis et al。14)和乔治奥Boudouvis (15),这是。

3.3。验证计算模型使用Oldroyd-B与粘弹性流体模型

基于K-BKZ模型,坦纳16]推导分析公式考虑粘性和弹性对膨胀率的贡献 在哪里第一法向应力差吗剪切应力。他们都是评估在墙上。引入的参数可恢复的剪切,肿胀比率可以估计

Oldroyd-B模型,充分发展的可恢复的剪切泊肃叶流等于

因为在我们的案例中和,可恢复的剪切。

网格收敛性进行了测试有七个网格,如表所示1。和的相对最小长度吗设在和分别设在。这些细胞出现出口附近的死亡。通过所有这七个网格的比较,发现的价值在确定湿润比中扮演主导的角色。此外,如果和、溶胀比的相对误差低于0.2%。该解决方案通过使用M6达到收敛。时间间隔是M7和在其他情况下(M1除外)。

奇点附近的压力值和聚合物应力如图2,在那里在的范围以及不同的网格密度(M4, M6)。所示,两条曲线非常接近,这表明给定的网格密度是充分的。的趋势,压力逐渐下降,然后经历大幅下跌之后,轻微的海拔在退出之前。随后,从正到负显著下降,后来趋于0。如数据所示2 (b)和2 (c),有一个最高点和最低点退出,这两点。此外,的曲线比一个更清晰的。作为增加的价值大幅减少,达到底部,长到一个极点,才逐渐下降至零。的概要文件的趋势以及奇点附近的最大和最小值,该方法的结果是接近的谱单元法采用三阶多项式老人工作(17),再次表明,当前网能够提供可靠的结果。小偏差的最大值和最小值的位置是归因于这样一个事实:速度是存储在本文和细胞的细胞中心顶点(17),分别。

一般来说,增加了、聚合物的组件的压力应该是类似但更大的价值。然而,这个数字远不同于直觉。图3显示的最低点和哪个是负峰这是积极的。力以及方向由聚合物引起的压力组件是,在那里垂直分量的标准界面是正常的。因此,更有可能正在接近一个(对应于一个小溶胀比温柔的曲线),大的影响是什么。随着力的增加,挤压越来越大。

膨胀率在图给出4与他人的作品。注意,在他们的作品中,相对参数的定义是不同的。例如,Russo和菲利普斯18和多美et al。19完整的宽度作为特征长度。我们只是改变他们在我们的工作是相同的一个。此外,钩针和Keunings [20.和多美et al。19雷诺数为0和0.5,分别。Russo和菲利普斯18)测试两种情况和结果揭示的事实没有重大变化时通过不同雷诺数足够小。所有的这些三篇论文类似的模拟。我们的模拟是一种类似于钩针和Keunings20.和多美et al。19那些不同的弛豫时间在不同的情况下;因此,曲线的趋势是相同的。然而,Russo和菲利普斯(18)改变了输入速度,也就是说,两壁附近的剪切速率和韦森伯数都发生了变化。在他的情况下,恢复的剪切速率远远不同于我们模拟一个即使是相同的。这是曲线的趋势的原因似乎与众不同。所有这些曲线显示的溶胀比韦森伯数的增加而增加。从图4上面,我们的结果是钩针和Keunings20.和多美et al。19]但低于Russo和菲利普斯[18]。的增加数,弹性效应决定的膨胀起着主导的作用。此外,轻微的肿胀比率较低数量见图4是由于负面的价值和也被观察到20.]。

稳态自由表面在不同的配置文件数字显示在图5。可以看出,一个更大的值将导致增加曲线的斜率。这是因为第一法向应力大的退出渠道大于较小的数量。这种现象也表明肿胀的主要因素是出口附近的第一法向应力值的死亡。在高数量,六面体单元变形大大增加了网非正交性。最大的网非正交性已达到高达40和阈值OpenFOAM通常设置为70。尽管非正交性值可以减少网格细化,会增加计算时间。从我们获得的最大溶胀比模拟大约是2.26使用网格M1。

速度等高线组件为,0.5,0.8,1.0为,显示在图1.06。通过增加数,速度的水平分量增加在的方向。的轮廓线扩大向上和向前。垂直速度的组件也会增加由于较大的弹性效应导致一个更大的溶胀比。

第一法向应力差的形象(在之前退出,退出死的大范围的如图7。这些线是采样设在以固定。显然,每一行的最大价值增加变得更大,这是一致的(15)。所有这些四行从0开始在中间的死()。至于,它逐渐上升,直到接近墙,在那里成长瞬间峰值。随后,它将显著下降。其他的资料值小于1.0的曲线转移向右向下。在细节,没有峰值的情况。一般来说,一个更大的导致更高的峰值,这表明第一法向应力占主导地位的挤压过程。

降低粘度比从1到0意味着逐渐增加的贡献弹性应力和减少纯粘性应力的影响。当达到零,模型降低upper-convected麦克斯韦(UCM)模型。为了调查的作用在膨胀模拟,模拟下和正在进行。图8说明了膨胀率将增加的趋势数量增加,而不管的价值。此外,结果离开制革厂商理论明显。考虑到一个极端的例子,方法1、溶胀比几乎仍是一个不断即使对于一个较大的值根据制革厂商理论。然而,这并不符合实验结果。的情况下,这个案子,,可恢复的剪切分别是3和2.7。根据(15),的溶胀比前一个应该大于后一个,而事实是相反的。Russo和菲利普斯18]属性差异,制革厂商公式来源于KBK-Z积分本构方程只考虑单一弛豫时间;然而,Oldroyd-B模型拥有一个滞后时间除了弛豫时间。这个结果可能与这一事实有关伸长的反应远比一个吗如图9。因此,具体的反应也在膨胀过程中扮演着重要的角色时,应该考虑预测最后的溶胀比。的仿真结果21]表明,流动具有足够高拉伸粘度能膨胀到一个明显大型国有即使在流率相对较低,忽略他们的表现在简单剪切流动。然而,降低粘度比会导致增加可采剪切,反过来,增加肿胀配给的价值根据(15)。通过比较数据的价值4和8合理地说,制革厂商理论仍具有一定的意义。除此之外,仔细比较,发现,当粘度比是非常大的(的情况下),这是一种非常稀液、肿胀率显著下降。相比之下,当相对较小(的情况下)、溶胀比不太敏感。这种现象也被观察到22]。注意,因为不同的粘度比的定义,本文是等价的在他们的作品。

我们可以看到从(13),剪切应力的值不会改变一旦剪切速率决定。因此,第一法向应力差是决定了溶胀比唯一的变量。数据7和10证明了高比例的聚合压力贡献将导致一个更大的价值膨胀,这将提供更大的力量。根据Oldroyd-B模型的解析解,小将导致更高价值的第一法向应力。因此,第一法向应力起着至关重要的作用在决定挤压流动的行为,可以用来预测溶胀比的最终结果。

4所示。结论

在这工作,跟踪新方法提出了粘弹性流体的自由表面。设置一个合适的自由表面边界条件的障碍以及解决了推导狄利克雷压力边界条件,引入平衡力。提出的数值方法是实现一OpenFOAM,广泛用于CFD模拟(23,24]。与先前的文献,比较结果表明,我们的解决方案是接近与其他方法(如有限元法、谱元法和MAC),这验证了该数值方法的正确性和有效性。此外,结果表明,膨胀率和粘度比的增加逐步下降。与此同时,简要分析粘度比的影响。此外,提供适当的本构方程,为新型高分子材料膨胀行为的调查可以执行基于我们的新平台。此外,有限体积法系统适合并行化以来,自由表面模拟在我们平台将很容易扩展到大量处理器并行计算。

利益冲突

作者宣称没有利益冲突有关的出版。

确认

我们应感谢曼彻斯特大学提供的计算资源,部分资金从中国的国家自然科学基金资助。61221491,没有。61303071,没有。61120106005,和开放基金(没有。201503 - 01和没有。201503 - 02)高性能计算的国家重点实验室。

引用

1. a·r·西奇,j·l·Favero n·s·m·卡多佐和h . Jasak“粘弹性流体分析内部和在自由表面流使用软件OpenFOAM”电脑和化学工程,34卷,不。12日,第1993 - 1984页,2010年。视图:出版商的网站|谷歌学术搜索
2. Osher和j . a . Sethian”方面与curvature-dependent传播速度:基于哈密顿雅可比算法公式,“计算物理学杂志,卷79,不。1,12-49,1988页。视图:出版商的网站|谷歌学术搜索|MathSciNet
3. 麦基,m·f·本诉g·费雷拉et al .,“MAC方法,”电脑&液体,37卷,不。8,907 - 930年,2008页。视图:出版商的网站|谷歌学术搜索|MathSciNet
4. c·w·赫特,a . a . Amsden和j·l·库克“任意lagrangianeulerian所有流速度计算方法(计算物理学杂志转载,14卷,pg 227 - 253, 1974),“计算物理学杂志,卷135,不。2、203 - 216年,1997页。视图:谷歌学术搜索
5. m·福丁和d . Esselaoui粘弹性有限元程序流。”国际期刊液体中的数值方法,7卷,不。10日,1035 - 1052年,1987页。视图:出版商的网站|谷歌学术搜索
6. z Tuković和h . Jasak移动网格有限体积界面追踪方法,表面张力界面流体为主,“电脑&液体,55卷,第84 - 70页,2012年。视图:出版商的网站|谷歌学术搜索|MathSciNet
7. h . Jasak“OpenFOAM:开源CFD研究和产业”国际造船和海洋工程》杂志上,1卷,不。2、89 - 94年,2009页。视图:出版商的网站|谷歌学术搜索
8. 美国Muzaferija和m . Perić”计算自由表面流动采用有限体积法和移动网格,”数值传热,B部分:基本面,32卷,不。4、369 - 384年,1997页。视图:出版商的网站|谷歌学术搜索
9. h . JasakŽ。Tuković”自动网格运动的非结构化有限体积方法,”交易的Famena,30卷,不。2,页1 - 2004。视图:谷歌学术搜索
10. n Phanthien和r·坦纳“粘弹性有限体积方法,”流变学系列,8卷,第359 - 331页,1999年。视图:出版商的网站|谷歌学术搜索
11. h . Matallah·汤森·m·f·韦伯斯特,“粘弹流,复苏和stress-splitting方案”非牛顿流体力学杂志》上,卷75,不。2 - 3、139 - 166年,1998页。视图:出版商的网站|谷歌学术搜索
12. m·A . Ajiz和A·詹宁斯“一个健壮的完整Choleski-conjugate梯度算法,”国际期刊工程中的数值方法,20卷,不。5,949 - 966年,1984页。视图:出版商的网站|谷歌学术搜索|MathSciNet
13. j·j·李,和c c。陆”,不完全Lu预处理对于大规模密集复杂的线性系统从电磁波散射问题,“计算物理学杂志,卷185,不。1,第175 - 158页,2003。视图:出版商的网站|谷歌学术搜索|Zentralblatt数学
14. e . Mitsoulis g·c·乔治奥,z Kountouriotis,”牛顿挤出物膨胀,各种影响因素研究”电脑&液体57卷,第207 - 195页,2012年。视图:出版商的网站|谷歌学术搜索|MathSciNet
15. g·c·乔治奥和a·g·Boudouvis聚合牛顿extrudate-swell问题的解决方案。”国际期刊液体中的数值方法卷,29号3、363 - 371年,1999页。视图:出版商的网站|谷歌学术搜索|Zentralblatt数学|MathSciNet
16. r·坦纳,工程流变学英国牛津,牛津大学出版社,2002年。
17. 美国老人,复杂的粘弹性数值模拟流使用间断伽辽金光谱/惠普元素方法[博士。论文),2013年。
18. 鲁索g和t . n .菲利普斯“光谱元素预测Oldroyd-B液体的膨胀,“电脑&液体,43卷,不。1,第118 - 107页,2011。视图:出版商的网站|谷歌学术搜索|MathSciNet
19. m . f .多美,n . Mangiavacchi j . A . Cuminato A . Castelo和s·麦基,“一个有限差分方法模拟非定常粘弹性自由表面流动,”非牛顿流体力学杂志》上,卷106,不。2 - 3、61 - 106年,2002页。视图:出版商的网站|谷歌学术搜索|Zentralblatt数学
20. m . j .钩针和r . Keunings有限元分析的死高弹性流体的膨胀,“非牛顿流体力学杂志》上,10卷,不。3 - 4、339 - 356年,1982页。视图:出版商的网站|谷歌学术搜索
21. a . Al-Muslimawi h . r . Tamaddon-Jahromi, m·f·韦伯斯特”模拟粘弹性和viscoelastoplastic膨胀流动,”非牛顿流体力学杂志》上卷。191年,45-56,2013页。视图:出版商的网站|谷歌学术搜索
22. J.-R。克莱蒙特和m . Normandin”数值模拟的挤出物膨胀Oldroyd-B液体使用流管分析和简化近似,”非牛顿流体力学杂志》上,50卷,不。2 - 3、193 - 215年,1993页。视图:出版商的网站|谷歌学术搜索|Zentralblatt数学
23. X.-W。杨郭,邹,x, x。元,m .王”界面不稳定和混乱的流变响应在二元聚合物混合物在剪切流,”RSC的进步,4卷,不。105年,第61177 - 61167页,2014年。视图:出版商的网站|谷歌学术搜索
24. X.-W。郭,W.-G。杨,X.-H。曹,y, X.-J。杨,“非平衡稳定状态的纠缠聚合物混合物在剪切流,”机械工程的发展,7卷,不。6、1 - 10,2015页。视图:出版商的网站|谷歌学术搜索

版权

版权©2015余曹等。这是一个开放的分布式下文章知识共享归属许可,它允许无限制的使用、分配和复制在任何媒介,提供最初的工作是正确引用。