SIMPLEC程序和SIMPLE程序相似。两种算法所使用的表达式唯一的区别就是表面流动速度校正项。和SIMPLE中一样,校正方程可写为:
??Jf?J*f?df?pc0?pc1?
但是系数d_f重新定义为:
??df??A2ap?a??fnb?
nb??可以看出,在压力速度耦合是得到解的主要因素时,使用修改后的校正方程可以加速收
敛。
PISO
压力隐式分裂算子(PISO)的压力速度耦合格式是SIMPLE算法族的一部分,它是基于压力速度校正之间的高度近似关系的一种算法。SIMPLE和SIMPLEC算法的一个限制就是在压力校正方程解出之后新的速度值和相应的流量不满足动量平衡。因此必须重复计算直至平衡得到满足。为了提高该计算的效率,PISO算法执行了两个附加的校正:相邻校正和偏斜校正。
PISO算法的主要思想就是将压力校正方程[69]中解的阶段中的SIMPLE和SIMPLEC算法所需的重复计算移除。经过一个或更多的附加PISO循环,校正的速度会更接近满足连续性和动量方程。这一迭代过程被称为动量校正或者邻近校正。PISO算法在每个迭代中要花费稍多的CPU时间但是极大的减少了达到收敛所需要的迭代次数,尤其是对于过渡问题,这一优点更为明显。
对于具有一些倾斜度的网格,单元表面质量流量校正和邻近单元压力校正差值之间的关系是相当简略的。因为沿着单元表面的压力校正梯度的分量开始是未知的,所以需要进行一个和上面所述的PISO邻近校正中相似的迭代步骤[51]。初始化压力校正方程的解之后,重新计算压力校正梯度然后用重新计算出来的值更新质量流量校正。这个被称为偏斜矫正的过程极大的减少了计算高度扭曲网格所遇到的收敛性困难。PISO偏斜校正可以使我们在基本相同的迭代步中,从高度偏斜的网格上得到和更为正交的网格上不相上下的解。
多相流中强体积力的特定处理
当多项流中存在较大的体积力(如:重力或者表面张力),动量方程中的体积力项和压力梯度项几乎是平衡的,相比较来说,对流项和粘性项的贡献就较小了。除非考虑压力梯度和体积力的局部平衡,否则分离算法的收敛性会很差。FLUENT提供了一种可选的隐式体积力处理,这种处理考虑了上面所说的影响从而使得解更具有鲁棒性。
基本的程序是将包含体积力的校正项增加到表面流动校正方程中(SIMPLE中的方程13)。这样,SIMPLE中的方程9就多了一个额外的体积力校正项,从而使得流动在迭代过程中提早得到真实的压力场。
这一选项只在多项流计算中使用,但是在默认情况下是关闭的。设定多相流计算的说明中包括了打开隐式体积力处理的说明,具体可以在以下几节中找到相关说明:在VOF计算中包括体积力,在气穴计算中包括体积力,在代数滑移混合计算中包括体积力。
除此之外,通过使用体积力的亚松驰因子,FLUENT允许你控制体积力中的变化。
耦合求解器
FLUENT中的耦合求解器同时解连续性、动量、(合适的话)能量和组分输运,并将它们作为一组控制方程或者方程的矢量来处理。随后会按顺序解附加标量的控制方程(也就是说这些附加标量方程相互之间是分离的而且和耦合方程组之间是分离的)。
矢量形式的控制方程
将曲面面积微分dA作为控制体积,积分描述平均流动属性的单组分流体的控制方程系统,相应的笛卡尔坐标形式为:
?WdV???F?G??dA??HdV
V?t?V其中,矢量W,F和G分别定义为:
?0??????v??????u???vu?pi??xi????????????W???v?,F???vv?p?j?,G???yi?
???w???vw?pk????zi???????????E????vE?pv????ijvj?q??矢量H包含了诸如体积力、能量等源项。
在这里,r、v、E和p分别是流体的密度、速度、单位质量的总能量以及压力。T是粘性应力张量,q是热流量。
总能量E和总焓H的关系为:
E?H?p?
其中:
H?h?v22
在方程1中所表示的Navier-Stokes方程在低马赫数下会非常具有刚性,这主要是因为流体速度v和声速c相差太大。对于不可压流来说这种情况尤其真实,不管流动速度有多大,不可压流体中的声速都是无穷大。在这些情况下,方程的数值刚性会导致较差的收敛速度。在FLUENT中,我们采用耦合求解器中的一种被称为(时间导数)预处理[175]的方法克服了这种困难。
预处理
时间导数预处理方法,是用预处理矩阵先乘以矢量形式控制方程中方程1的时间导数项。这一步就重新标度了所解方程系统的声速(特征值),从而减轻了低马赫数和不可压流动中会遇到的数值刚性的影响。
推导预处理矩阵,首先是用微分学中的一个法则(具体看方程就知道了)将守恒量W在控制方程的方程1的因变量形式变形为原始变量Q的形式,结果如下:
?W?QdV???F?G??dA??HdV ?VV?Q?t其中,Q是矢量{p, u, v, w, T }^T,雅克比(Jacobian)矩阵?W/?Q为:
??p??up?W?????pv?Q??w?p???pH其中:
?T???00?Tu??0?0?Tv? 00??Tw???u?v?w?TH??Cp??000?p?????,?T?
?Tp?pT对于理想气体d = 1,对于不可压流体d = 0。
选择原始变量Q作为因变量有几个原因。首先,解不可压流动时它是自然的选择。其次,当我们使用二阶精度时,为了得到粘性流动中更高精度的速度和温度梯度以及无粘流动的压力梯度,我们需要重建Q而不是W。最后,选择压力作为因变量允许系统中的声波的传播被挑选出来[165]。
我们用预处理矩阵G来替换雅克比矩阵?W/?Q(方程3)来实现方程的预处理,这样,预处理系统的守恒形式为:
??QdV???F?G??dA??HdV ?VV?t其中:
0????u???????v0??w0????H?u参数?为: ?1?T???2??U?r?Cp?T??00?Tu???0?Tv? 0??Tw???v?w?TH??Cp??00?? ??方程9中的参考速度U_r在当地选取,从而保证系统的特征值关于对流和耗散时间尺度能够调节的很好[175]。
预处理之后的方程系统(方程5)的特征值为:
u,u,u,u??c?,u??c?
其中:
F?1??Q ?FR?FL??1?A22在这里,d Q是空间差分Q_R - Q_L。流量F_R = F(Q_R)和F_L = F(Q_L)是用表面“右边”
和“左边”的(重建解)矢量Q_R和Q_L计算出来的。矩阵|A(hat)|定义为:
??M?M?1 A其中L是特征值的对角阵,M是G^-1A对角化的模式矩阵,其中A是无粘流量Jacobian矩阵?F/?Q。
对于非预处理系统(和理想气体),方程13化简为在Roe平均用于估计G,|A(hat)|时的Roe流量差分分裂[133]。目前所使用的是Q_R和Q_L的平均。
在当前形式中,方程13可以看成二阶中心差分加上附加的矩阵耗散。附加的矩阵耗散项不仅对对流变量迎风的产生和超音速流动的压力和流量速度有作用,而且还提供了定常状态的压力速度耦合,并对低速和不可压流的稳定性和高效收敛有作用。
定常流的时间步进法
FLUENT中耦合控制方程组(预处理中的方程5)在时间上的离散既可以用于定常计算也可以用于非定常计算。在定常情况下,进行时间步进计算直至达到定常解。耦合方程可以用隐式也可以用显式时间步进来进行时间离散。这两种算法在下面的显式和隐式格式两节中介绍。
显式格式
在显式格式中,预处理的方程5是用多步,时间步进算法[73]来离散时间导数的。从第n步迭代到第n+1步迭代是用m步Runge-Kutta格式来提高解的,m步Runge-Kutta格式为:
R?iNfaces??F?Q??G?Q???A?VH
ii时间步D t是从CFL(Courant-Friedrichs-Lewy)条件计算得到的
?t?CFL?x?max
其中l_max是预处理一节中方程11定义的局部特征值的最大值。
使用全近似多重网格方法(FAS)可以加快显示格式的收敛速度,关于FAS请参阅后面将要介绍的内容。
通过与相邻残差的隐式平均增加格式的支持可以在后面的计算中增加最大时间步长。残差通过Laplacian 光滑算子进行过滤,其中Laplacian 光滑算子为:
Ri?Ri????Rj?Ri?
该方程可以用下面的Jacobi迭代解出:
Riim?Ri???Rjm?11???1(译者注:分母为什么这样写?)
对于e=0.5,两次Jacobi迭代就足可以允许时间步增加一倍。
隐式格式
在显式格式中,控制方程(预处理一节中的方程5)的时间Euler隐式离散和流量的牛顿类型线化界和得到下面delta形式的线化系统[174]:
Nfaces??n?1D?S??Rn ??j,k??Qj??中心D和对角线下面的系数矩阵S_j,k分别为:
VD????tNfacesj?Sj,i
??Fj?GjSj,k????Q??Qk?k??? ?残差矢量R^n和时间步D分别在显式格式中的方程1和方程3中定义。
方程1是用点Gauss-Seidel格式和适应耦合的一组方程的代数多重网格方法解出的。
非定常流的时间离散
对于时间精度计算,可以使用显式和隐式时间步进格式。(隐式方法也被称为双重时间步进法)
显式时间步进
在显式时间步进方法中,使用显式格式一节所介绍的方法,在区域的每一个单元内使用相同的时间步,取消预处理选项。
双重时间步进
为了得到预处理方程的高效的,时间精度解,我们采用双重时间步进法和多步格式。在这里我们在矢量形式的控制方程的方程1中引入伪时间导数:
??WdV???t?V???QdV???F?G??dA??HdV
VV其中t是物理时间,t是时间步进程序中的伪时间。注意:方程1左手边的第二项t??消失了,而重新获得了矢量形式控制方程中的方程1。
方程1中的时间相关项是用一阶或者二阶精度后向时间差分进行隐式离散的。写成半隐式形式如下:
A?e?b?0
其中f_e是精确解。在解收敛之前可能会有一个和近似解f有关的误差d:
A??b?d
我们寻找一个校正f的y,这样,精确解由下式给出:
?e????
将方程7代入到方程3中有: