有限元法等传统基于网格的算法在针对一些情况复杂的问题，比如移动边界、界面变化、大变形(流固耦合、高速撞击等)、自由表面流等问题时，由于计算网格可能会出现畸变、重构等缺陷，从而限制了该类算法在数值模拟中计算效率及精度。而无网格方法^[1]是将求解域划分为一系列节点，不需要对求解域进行网格划分^[2]。没有了网格也就没有了网格之间拓扑关系的制约，能够部分或者彻底地摆脱计算过程对于网格的依赖^[3]，这样就可以应用于结构大变形、动态裂纹扩展、结构优化等领域。无网格方法最早提出主要是用来解决天体问题，在1977年Lucy提出光滑粒子法(SPH)并且成功应用于天体物理学领域中^[4]。后来Johnson等针对轴对称几何问题提出了归一化平滑函数方法(NSF)^[5]，通过2个圆柱体碰撞数值验证表明该方法提高了SPH法的计算精度。在1981年Lancaster等在传统的最小二乘法的基础上提出移动最小二乘法(MLS)^[6]，称为构造无网格法的形函数的新方法，该方法能够解决传统的最小二乘法离散数据量大等缺点。后来Belytschko将MLS应用到无网格方法中，最先提出无网格伽辽金方法(element-free Galerkin method)方法，使用高斯积分方法对求解域中的背景网格进行积分^[7-9]，并将EFG方法应用于弹性问题、热传导问题^[10-12]，计算结果表明在收敛性上及精度上EFG法比有限元法具有很明显的优势^[13-17]。

1 EFG法理论与基础

无网格伽辽金算法采用移动最小二乘技术来构造其形函数，假定二维流场中的常函数在节点x邻域Ω_x内局部近似为u(x)

(1)

式中, u^h(x)为函数u(x)的近似表达式, p_j(x)为基函数, 本文中基函数主要取线性基函数, 即p_j(x)=[1, x, y]^T, m为基函数的个数; b_j(x)为系数向量, b_i(x)为待定系数。b_i(x)可以通过求近似函数u^h(x)在点x的邻域Ω_x内误差加权平方和的最小值来确定。

(2)

(3)

求解J(x)关于b(x)的极值, 最终得到的场函数近似表达式如下所示, 其中Φ(x)为形函数。

(4)

(5)

计算得到的A(x)和B(x)的表达式如下所示, ω_i(x)为权函数。

(6)

(7)

由于使用EFG方法构造的形函数不满足Kronecker-δ性质, 因而不能像有限元法那样直接施加本质边界条件, 无网格方法中对于边界条件的施加主要有罚函数法、拉格朗日法、配点法等。在本文中主要使用罚函数法(penalty function method, PM)施加本质边界条件, PM法主要是在本质边界条件中引入一个罚因子α, 使用该方法可以简化离散方程中的系数矩阵, 施加后系数矩阵依然具备带宽及正定的性质, 而且未引进多余的变量数。其中对任意表达式ξ的PM法实施过程如下所示

(8)

式中, u代表某一变量, u_i表示该变量的取值。

本文中权函数选用三次样条权函数, 其形式如下所示

(9)

在EFG法中, 节点处的影响域形状一般有矩形和圆形, 本文中为了简化计算选取矩形影响域, 节点i处的矩形支持域尺寸表达式如下

(10)

式中, d_xi和d_yi分别为矩形支持域的长度和宽度, 其中d_max是控制节点影响域大小的比例参数, 因此d_max的取值直接关系到数值结果的计算精度。d_max的大小一般要求满足以下2个条件:①要保证所有的节点权函数的影响域的并集能覆盖整个求解域; ②保证局部影响域包含的节点数不能少于基函数的项数。为了提高精度和保证矩阵可逆运算, 节点权函数的影响域应尽可能大些, 但为了保持逼近函数的局部特性, 节点权函数的影响域也不能太大, 否则会无谓增加计算量, r_x和r_y分别为水平方向和竖直方向上相邻两点之间的间距。则可以得到权函数的具体表达式为

(11)

由于使用MLS近似方法构造EFG方法的近似函数会需要进行矩阵求逆, 求逆过程计算量较大, 不同于FEM方法, EFG法中的近似函数大都不属于多项式, 无法直接求出其积分表达式, 只能够对其进行数值积分, 因此EFG法需要在求解域中设置一定数量的背景网格, 并且在背景网格中进一步设定一系列有规律的积分节点, 这些节点就是高斯积分点, 然后对每一个背景网格中高斯点进行高阶积分来保证求解精度。背景网格及高斯积分点划分模型如图 1所示。

2 控制方程离散形式求解 2.1 计算流体力学基本数学方程

对于某一个具体的流动问题, 首先确定能够表达其问题的微分方程和边界条件。对于一般的计算流体力学问题, 假设该流体是牛顿流体, 且不考虑热效应的影响, 则考虑其惯性项和时间项的二维不可压缩黏性流动问题数学表达式及其边界条件为:连续性方程

运动性方程

边界条件

2.2 加权余量方程

流场内任意一个节点的速度可以表示为该节点处的插值函数与该节点邻域内所包含系列节点的速度乘积之和, 同样也适用于节点处压力的求解, 可以近似表达为:

(12)

对控制方程的离散推导使用Galerkin离散方法, 将连续性方程和运动性方程与权函数相乘并在计算域内进行积分, 权函数与插值函数保持一致, 得出其加权积分形式如下所示:

(13)

将(13)式中的时间项和对流项中的密度提出, 其余项在计算域内应用分部积分得到其在自然边界条件下的弱形式积分表达式如下所示:

(14)

2.3 压力项处理

对于公式(14), 引入足够大的系数罚因子λ, 则压力项可以用速度项表示, 这样可以将连续性方程引入到运动方程中, 从而简化求解方程的形式, 使用速度和罚因子表达的压力方程如下所示:

(15)

将(15)式带入弱形式积分表达式中, 得到:

(16)

以及

(17)

2.4 施加本质边界条件

采用罚函数法施加本质边界条件, 引入罚因子α, 施加边界条件后的方程如下所示:

(18)

以及

(19)

2.5 惯性项离散

对于惯性项的离散, 一般有速度提出法和直接推导法2种。速度项提出法是将节点邻域内所包含的系列节点速度向量提出, 其实施过程如下所示:

(20)

使用直接推导法离散惯性项是将节点的速度直接使用形函数表示, 其实施过程如下所示:

(21)

将施加边界条件后的方程写成如下所示的矩阵方程形式, 其中使用速度项提出法离散惯性项后所对应的矩阵方程如下所示:

(22)

使用直接推导法离散惯性项后所对应的矩阵方程如下所示:

(23)

式中，和分别代表x方向和y方向的速度的一阶导数, u, v分别代表x方向和y方向的速度, 上述方程中各个系数矩阵表示如下所示:

(24)

(25)

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

3 非线性流动问题研究 3.1 定常非线性流动计算 3.1.1 数学模型

对于定常流动问题, 在求解时需要将N-S方程中的时间项去掉。考虑惯性项的情况下, 方程成为非线性方程, 对于非线性项的处理为Newton-Raphson迭代法和直接线性交替迭代法。使用N-R迭代方法求解使用速度项提出法离散惯性项后所对应的矩阵方程, 得到矩阵方程如下式所示。

(34)

式中

(35)

使用直接线性交替迭代方法处理使用直接推导法离散惯性项后所对应的矩阵方程, 可以得到矩阵方程如下式所示。

(36)

3.1.2 算例分析

首先采用EFG法求解平板相向运动的流动问题, 示意图如图 2所示。该模型可以作为对不可压流场模拟结果的验证。

矩形域主要由3块平板组成, 坐标系如图所示; 流场域的长度w取值为1 m, 宽度H取值为1 m; 流场内部的流体密度取值ρ=1 000 kg/s², 流体黏度取值μ=0.001 N·m/s²。

通过推导, 得出该流动问题的解析解表达式如下所示:

(37)

程序中的参数设置如下:EFG法中的流场域内离散结点均匀分布, 节点分布为441(21×21);压力罚因子λ取值为10¹⁰, 速度罚参数α取值为10⁶; 速度采取线性插值; 支持域的缩放系数d_max取1.2^[18]。解析解中节点设置与EFG法相同。

本算例中将使用直接线性交替迭代法处理非线性问题并进行数值模拟, 并作出流场稳定后速度分布和流线图, 并与解析解结果进行对比分析, 图 4和图 5展示了EFG法与解析解的速度云图, 对比可以发现, 2种算法的速度云图趋势一致。观察EFG解的速度场云图可以发现, 流体满足不可压缩性。说明使用EFG模拟可以得到符合实际情况的流场分布。

在左端附近和右端出口附近分别选取x=0.2 m和x=0.8 m的2条竖线, 并分别在这2条垂直线上选取一系列节点, 求解这些节点处的速度数值解以及EFG解和解析解之间的相对误差, 其中的速度单位为m/s, 速度及相对误差结果展现在表 1和表 2中。从表中可以看出对于所选择的系列节点, EFG解与解析解的相对误差最大不超过3.66%, 从而可以验证EFG算法在此算例中的精度比较高。

3.2 非定常非线性流动计算 3.2.1 数学模型

对于非定常问题的求解, 需要在定常问题的迭代求解过程中外加一个时间迭代层, 当每一迭代时刻计算N(N的选取与所计算的问题有关)个迭代步后, 即可进入下一个时间迭代层的计算。对于时间项的离散, 采用向前差分法处理N-S方程中的时间导数项, 时间项可以表示如下:

(38)

采取θ加权法, 引入θ因子, 将方程中的速度项可以表示成当前时间步长n+1与上一时间步长n的速度分布结果之和, 其表示式如所示

(39)

式中, θ在0到1之间取值, 将θ取值为0则上式为显示格式, 此时每求解每一个时间步时不需要求解非线性方程, 缺点是该方法极不容易收敛; θ取值为2/3时则为伽辽金格式; θ取值为0则为隐式格式, 此时在每一个时间步长内仍需求解非线性方程, 不过此时计算收敛速度较快。

时间项采用全隐式离散形式, 即将三次样条权函数中的速度表达成为隐式格式, 此时速度表达式如下所示:

(40)

将隐式形式的速度的表达式带入三次样条权函数方程得到总体N-S方程的隐式向后差分格式如下式所示:

(41)

将θ取值为0即此时时间项采取全显式离散格式, 将三次样条权函数方程中的速度均表达成为隐式格式, 此时速度表达式如下:

带入三次样条权函数方程得到下式, 其中Δt为时间步长。

(42)

式中，f_u与f_v的表达式如下所示:

(43)

这个算例是对二维方柱进行绕流计算分析。流场域结构如图 6所示, 该流场域的长度为8 m, 宽度为4 m, 在流场区域中节点按照均匀布置, 均匀分布5 290个节点改变不同的来流速度研究速度对方柱绕流流场的影响。

3.2.2 算例分析

针对一系列不同雷诺数的方柱绕流问题, 基于EFG法的非定常算法进行数值模拟, 计算时间步长Δt取1 s。

当来流速度较小时, 即雷诺数Re小于50时, 流场的初始条件按照还未发生改变的流动状态施加, 采用直接线性交替迭代法求解, 计算得到稳定后方柱周围的速度云图以及流线图如图 8所示。

当来流速度继续增加时, 非线性作用更加强烈, 此时需要使用Newton-Raphson迭代法和直接线性交替迭代法来计算, 得出的方柱周围流场内的法向速度云图以及流线图如图 9所示。

图 10展示了Re为100时分别使用EFG法和FEM法计算得到的竖线x=3.5 m上的水平方向速度拟合图。从图中看出2种方法计算得到的速度拟合图拥有很好的吻合度。再将稳定后的流场结构与使用EFG法稳态求解结果对比发现两图流动趋势一致, 并且回流区长度也一致, 由此可见, 在适当的时间步长下, 使用非定常EFG求解法同样能够细致入微地描述方柱周围流体的流动过程。

4 结论

本文针对典型的非线性流动问题进行数值研究, 主要工作内容及结论如下:

1) 借鉴FEM算法中的Garlerkin离散方法, 对N-S方程进行EFG法离散, 使用罚函数法施加压力边界条件和速度边界条件; 分别采取速度项提出和直接推导法离散方程中的惯性项。

2) 针对定常非线性流动问题, 首先推导了使用Newton-Raphson迭代法和直接线性交替迭代法求解定常非线性N-S方程的矩阵形式。并以矩形域上下平板相向运动流动问题为例进行数值模拟, 验证了EFG法的可靠性。

3) 针对非定常非线性流动问题, 使用θ加权法对N-S方程中的时间项离散, 并推导了EFG法非定常求解公式; 以一系列不同雷诺数的方柱绕流问题为例, 使用EFG法进行数值实验。