摘要： 为比较最小二乘有限元法（Least Square Finite Element Method，LSFEM）和有限体积法在CFD应用中的优劣，采用最小二乘法离散不可压N-S方程的有限元模型，得到正定对称线性系统，采用高效的预处理共轭梯度法求解方程组；利用LSFEM和基于有限体积法的FLUENT分别计算Kovasznay流动、定常二维和三维后台阶流动以及非定常圆柱绕流等4个实例并比较计算结果. 结果表明，LSFEM比有限体积法的收敛性和精确性更好，在CFD领域的应用价值很高.
关键词：
最小二乘有限元法； 有限体积法； CFD； 预处理共轭梯度法； Kovasznay流动； 后台阶�流动； 圆柱绕流
中图分类号： O35；TB115.1
文献标志码： �A
Application comparison of least square finite element method and finite volume method in CFD
TAO Sha�1, YANG Zhigang�1, JIANG Bonan�2, GU Wenjun�1
(1. Shanghai Automotive Wind Tunnel Center, Tongji University, Shanghai 201804, China;
2. Department of Mathematics and Statistics, Oakland University, Miami 48309, USA)
Abstract： To compare the advantages and disadvantages of Least Square Finite Element Method(LSFEM) and finite volume method in actural applications, a least square method is used to discrete the finite element model of incompressible N-S equation, the positive definite symmetric linear system is obtained, and the equations are solved by efficient preconditioned conjugate gradient method; the least square finite element method and FLUENT which is based on finite volume method are separately used to calculate the examples such as Kovasznay flow, steady 2D and 3D backward-facing step flows, and unsteady flow around cylinder. The comparison results indicate that the least square finite element method performs better convergence rate and accuracy than finite volume method and shows good application value in CFD.
Key words：
least square finite element method; finite volume method; CFD; preconditioned conjugate gradient method; Kovasznay flow; backward-facing step flow; flow around cylinder
0 引 言
近30年来，最小二乘有限元法（Least Square Finite Element Method，LSFEM）在流体动力学中的应用得到越来越多的关注，其引入残量最小化概念求解N-S方程，具有很强的稳定性，不需要验证LBB条件［1］，能更灵活地选择有限元计算域.此外，LSFEM求解各种条件下N-S方程的离散形式都正定对称，有利于选择高效、统一的求解方法而不需要任何迎风格式和调节参数，使得基于LSFEM的计算代码具有很强的通用性.［2］
目前，大多数针对LSFEM的研究主要侧重于可行性探索和论证［2-8］，尚未见对其实际应用价值的研究.本文用基于LSFEM的计算代码模拟各种条件下的流体算例，与现今广泛应用的基于有限体积法的FLUENT软件进行比较分析，探索LSFEM在CFD领域推广应用的价值和途径.
1 不可压N-S方程的LSFEM简介
1.1 不可压N-S方程1阶偏微分形式
不可压N-S方程可写成无量纲形式
式中：U为来流特征速度；L为流场特征长度；ν为运动黏度.
由于式(1)中含有速度的2阶导数，直接对其使用�LSFEM�需要采用连续可微的形函数.为避免该问题，可引入独立未知量的涡量ω[WTBX]=（ω�x,ω�y,ω�z）=��Δ�×u[WTBX]，得不可压�N-S�方程的“速度-涡量-压力”1阶偏微分形式［2］
由于式(3)中的对流项u・�Δ�u为非线性，需用牛顿法进行线性化处理，得
u�0・�Δ�u+u・�Δ�u�0+�Δ�p+
1Re�Δ�×ω=
f+u�0・�Δ�u�0
(4)
将上述方法运用于笛卡尔坐标系下的二维问题，可得1阶偏微分方程组
��
u�0�u�x+u�u�0�x+v�0�u�y+v�u�0�y+�p�x+1Re�ω�y=f�x+u�0�u�0�x+v�0�u�0�y
u�0�v�x+u�v�0�x+v�0�v�y+v�v�0�y+�p�y－1Re�ω�x=f�y+u�0�v�0�x+v�0�v�0�y
ω+�u�y－�v�x=0
�u�x+�v�y=0(5)
1.2 基于1阶偏微分方程组的LSFEM
用最小二乘法求解1阶偏微分方程组的基本原理为残差的L2范数最小化.考虑如下1阶偏微分静态边值问题［3］
Au=f， x∈Ω
Bu=g， x∈Γ[JB)][JY](6)
式中：Ω为有界实数空间；Γ为Ω的分段连续边界；f和g为给定向量；B为边界代数算子；A为1阶偏微分算子，
Au=�n��d�i=1A�i�u�x�i+A�0[WTHX]u[JY] (7)
式中：n��d�为空间维度；[WTHX]u = (u��1,u��2,…,u��m)�T�代表m个[WTHX]x=(x�1,x�2,…,x��n��d�)的未知函数；[WTHX]A�i和[WTHX]A�0为随[WTHX]x变化的N���eq�×m矩阵，N���eq�为方程组中的方程数量.
设式(6)解的试探函数为ν，ν属于希尔伯特空间L�2(Ω)且满足边界条件
Bν=g， x∈Γ (8)
并定义残差函数
R=Aν－f
(9)
求解式(6)可得适当的ν，使残差的�L2�范数‖[WTHX]R‖�2取最小值.
令
运用变分法求解
��lim�t→0�d��d�tI([WTHX]u +[WTHX]vt)≡2∫Ω ([WTHX]Av)�T�([WTHX]Au-[WTHX]f)�d�Ω=0 (11)
得
B(u,v)=F(v)
B(u,v)≡(Au,Av)
F(v)≡(f,Av)[JB)][JY](12)
将式(5)所示的二维�N-S�方程组写成如式(6)所示的1阶偏微分系统，得
A�0=0000�u�0�x�u�0�y00�v�0�x�v�0�y000001
A�1=1000u�00100u�001Re0－100
A�2=0100v�0001Re0v�0101000
f=0f�x+u�0�u�0�x+v�0�u�0�yf�y+u�0�v�0�x+v�0�v�0�y0
u=uvpω[JB)][JY](13)
利用式(12)所示方法求解式(13).式(12)为对称双线性系统，适于求解偏微分方程的各类边值问题，因此LSFEM在工程中各个领域的应用具有统一的计算格式，不需另加数值差分格式.当遇到不同的问题时，只需重新输入相对应的A�0，A�i和f值即可，给编写通用的求解代码带来很大的便利.
引入有限元分析，将原定义域分解为一系列子元素，在每个子元素内使用适当的形函数
因此，LSFEM对任意适定的式(6)进行离散化，都可得到正定对称矩阵K[WTBZ]，有利于选用高效、统一的迭代方法求解方程组.本文采用预处理共轭梯度法［9］，该法在求解大型正定对称矩阵系统时具有很高的迭代率和收敛性.
2 计算实例和对比结果
2.1 Kovasznay流动
考虑具有精确解的Kovasznay流动［10］，验证LSFEM离散定常不可压N-S方程的收敛性.KOVASZNAY［10］给出如下满足N-S方程的精确解
u(x,y)=1－�e�λx�cos�(2�π�y)
v(x,y)=λ2�π��e�λx�sin�(2�π�y)
p(x,y)=p�0－12�e�2λx
λ=Re2－(Re�24+4�π��2)1/2(19)
式中：�p��0为参考压力(任意常数).Kovasznay流动的结果与二维周期圆柱列后的尾迹流动非常相似，图1和2分别为�Re�=40时精确解的流函数分布和速度�u�的分布.
根据离散方法收敛率的定义可知，随着网格尺寸的不断减小，计算结果的L2误差也逐渐减小，并且二者在对数坐标系中呈线性关系，其斜率即为收敛率.本文在［-0.5,1.5］×［-0.5,1.5］的计算域内分别划分8×8，16×16，32×32和64×64个网格进行数值计算，计算域边界上采用精确解的边界值作为第一类边界条件.分别采用LSFEM代码与FLUENT计算得到速度�u和v�残差的L2范数收敛率，见图3.
由图3可知，2阶精度LSFEM的收敛率达到4，表现出超收敛性；而2阶精度FLUENT的收敛率只有2，因此LSFEM的收敛性远优于FLUENT.
2.2 定常二维后台阶流动
后台阶流动是工程中常见的流动现象之一，其流动现象较复杂，本文仅对其层流范畴内的流动进行分析.ARMALY等［11］对后台阶分离流进行较为系统和全面的实验研究，本文参考的数据来源于他们的实验结果.根据ARMALY等的实验设置设定后台阶几何模型，见图4，几何模型的尺寸见表1，各边界条件和待测参数设置见图5.
用LSFEM和FLUENT模拟二维后台阶流动所得值与实验值对比见图7，可知，当�Re�小于400时，LSFEM得到的数值结果比FLUENT更接近于实验结果；当�Re�大于400时，由于实验存在的三维效应，2种数值解与实验值都有一定偏差；此外，FLUENT在�Re�大于600后，结果出现较大振荡，无法收敛，而LSFEM仍然得到稳定解，也说明LSFEM的稳定性优于FLUENT.
2.3 定常三维后台阶流动
为进一步验证LSFEM对于三维流动计算的适用性，同时也为排除实验结果中由于三维效应影响所产生的误差，本文在上述二维后台阶流动算例基础上进行三维后台阶流动数值模拟.
LSFEM和FLUENT模拟三维后台阶流动所得�x��1值与实验值对比见图8，可知，当�Re�较大时，三维计算结果比二维计算结果更接近实验结果，实验值确实受到三维效应的影响；LSFEM的三维计算结果明显比FLUENT更接近实验结果，进一步证明LSFEM在求解定常N-S方程时具有比FLUENT更高的精确性.此外，LSFEM可采用更高阶的精度以得到更精确的数值结果.因此，在对计算精度要求较高时，LSFEM更胜任.
2.4 非定常圆柱绕流
圆柱绕流不仅是经典的流体力学问题，而且在工程实际中也具有重要意义.当流场中�Re�大于40时，圆柱上的附着涡瓦解，下游流场不再定常并在后缘上、下两侧出现周期性轮流脱落的涡，形成卡门涡街.WILLIAMSON［13］曾对该流动情况下�Re�与斯特劳哈数的关系进行过详尽研究.本文以WILLIAMSON的研究结果为参照，考察2种数值方法求解非定常N-S方程的准确性.
圆柱绕流几何模型见图9，定常速度入口和�Re�分别取100和200.分析通过2种数值方法得到的斯特劳哈数，LSFEM和FLUENT模拟非定常圆柱绕流所得到的斯特劳哈数与参照值对比见表2，可知，LSFEM在非定常数值计算中也具有较高的准确性.
3 结束语
简介LSFEM求解不可压N-S方程的过程，由此具体说明LSFEM的优点.通过计算实例证明，相对于目前广泛应用于CFD的FLUENT，LSFEM具有更高的收敛率和精确性.此外，不可压N-S方程的LSFEM求解代码不需要另加任何迎风格式和调节参数，代码简单、通用，避免人为设置参数导致的结果偏差问题.从理论上讲，LSFEM继承FLUENT的高阶特性，可采用高阶网格达到更高精度，具有很高的实际推广价值.
参考文献：
�［1］ PEHLIVANOV A I, CAREY G F, LAZAROV R D. Least-squares mixed finite elements for second-order elliptic problems［J］. SIAM J Numer Anal, 1994, 31(5): 1368-1377.
�［2］ JIANG Bonan, CHANG C L. Least-squares finite elements for the Stokes problem［J］. Comput Methods Appl Mech & Eng, 1990, 78(3): �297-311.
�［3］ JIANG Bonan. The least-squares finite element method: theory and applications in computational fluid dynamics and electromagnetic (scientific computation)［M］. Berlin: Springer, 1998: 3-68.
�［4］ BOCHEV P B. Analysis of least-squares finite element methods for the Navier-Stokes equations［J］. SIAM J Numer Anal, 1997, 34(5): �1817-1844.
�［5］ BOCHEV P B, GUNZBURGER M D. Least-squares finite element methods［M］. Berlin: Springer, 2009: 69-364.
�［6］ PONTAZA J P, REDDY J N. Least-squares finite element formulations for viscous incompressible and compressible fluid flows［J］. Comput Methods Appl Mech & Eng, 2006, 195(19-22): 2454�2494.
�［7］ PONTAZA J P, REDDY J N. Spectral/hp least-squares finite element formulation for the Navier-Stokes equations［J］. J Comput Phys, 2003, 190(2): 523-549.
�［8］ KWON K C, PARK S H, JIANG B N, �et al�. The least-squares meshfree method for solving linear elastic problems［J］. Comput Mech, 2002, �30(3): 196-211.
�［9］ CAREY G F, JIANG B N. Least-squares finite element method and preconditioned conjugate gradient solution［J］. Int J Numer Methods Eng, 1987, 24(7): 1283-1296.
�［10］ KOVASZNAY L I. Laminar flow behind a two-dimensional grid［J］. Math Proc Cambridge Philos Soc, 1948, 44(1): 58-62.
�［11］ ARMALY B F, DURST F, PEREIRA J C F, �et al�. Experimental and theoretical investigation of backward-facing step flow［J］. J Fluid Mech, 1983（127）: 473-496.
�［12］ ERTURK E. Numerical solutions of 2D steady incompressible flow over backward-facing step, part I: high Reynolds number solution［J］. �Comput & Fluids, 2008, 37(6): 633-655.
�［13］ WILLIAMSON C H K. A series in 1/�Re� to represent the Strouhal-Reynolds number relationship of the cylinder wake［J］. J Fluids & Struct, 1998, 12(8): 1073-1085.
(编辑 陈锋杰)