1. 引言
文献 [1] 中提出了一种最小二乘等几何方法，从最小二乘泛函导出变分方程，采用NURBS基函数构造有限元空间，用该方法求解了线性的Stokes流动，文献 [2] 用牛顿法线性化对流项，将最小二乘等几何方法推广到非线性N-S方程的求解，研究了稳态流动的模拟，文献 [3] 对时间项采用隐式差分格式进行离散，研究了该方法在非稳态流动模拟中的应用。研究表明最小二乘等几何方法兼有最小二乘有限元 [4] 与等几何分析 [5] 的优点，而且可以引入多重网格方法以提高计算效率 [6] ，对于雷诺数较小的流动问题的能给出较好的数值结果。
虽然文献 [1] [2] [3] 通过构造解析解或求解基准流动问题初步验证了最小二乘等几何方法的适用性，但是对较复杂流动的适用性需进一步研究，本文以突然启动圆柱绕流问题为对象，研究该方法求解较复杂流动的特性。直径为D的圆柱位于静止的流体中，从零时刻开始突然以恒定的速度U匀速运动。圆柱突然启动后的瞬态流动是典型的非定常分离流动，高雷诺数下流动的复杂性给数值模拟带来了挑战，研究人员用不同方法对不同雷诺数下的流动进行了模拟，包括有限差分法 [7] [8] 、涡量法 [9] [10] 、LBM法 [11] [12] [13] 、小波配位法 [14] 等。其中涡量法是求解这个问题的较为有效的方法，但圆柱表面的无滑移边界条件处理困难。另外，由于圆柱从静止突然加速到U，在初始时刻存在时间上的奇异性，这也给数值模拟带来了困难。Bouard [15] 等对突然启动圆柱绕流进行了实验研究，所考虑的雷诺数(
)最大为9500，实验发现，随着雷诺数的增加，流动中会产生鼓包现象、孤立的二次涡、α-现象和β-现象等不同流动特征。
本文对雷诺数为550~9500的流动演化过程进行了数值模拟，数值解准确捕捉到了不同雷诺数下的产生的流动特征，计算结果与文献中的数值及实验结果进行了对比。
2. 控制方程离散
最小二乘等几何方法用控制方程的残差定义最小二乘形式的“能量”泛函，由该泛函导出变分问题，再由NURBS构造有限维空间进行控制方程离散。该方法对不具有自伴随性质的流动问题也能得到对称、正定的系数矩阵；计算网格在分析中自动生成，不需要人工划分网格；除了采用加密网格、提高多项式次数来提高数值解精度，还可以通过升阶算法提高NURBS基函数的光滑性来提高精度。
非定常粘性不可压流动由N-S方程描述，其控制方程的无量纲形式可写成：
(1)
式中：u为流动速度，p为压力，f为流体所受到的体积力，Re为雷诺数，Ω为流动区域，Γ是Ω的边界，
为时间区间，u0为初始时刻的速度分布，us为边界上的速度分布。初始速度应满足不可压条件以满足适定性要求。前面两式分别反映了流体流动中的动量守恒和质量守恒，后面两式为初始条件和边界条件。
先对控制方程中的时间项用稳定性好的全隐式格式进行离散，然后在每一个时间步上用牛顿法进行线性化，最后用最小二乘等几何方法迭代求解。时间离散精度在程序中通过向后多步差分的步数进行控制，空间离散精度在程序中通过NURBS升阶进行控制，具体过程见 [2] [3] 。
3. 计算模型
本文采用Matlab实现了求解二维非定常粘性不可压流动最小二乘等几何算法，虽然在高雷诺数下圆柱绕流最终发展成三维流动，但是实验表明在流动的初期仍然可以近似看成二维的，流场分布具有对称性。计算区域如图1所示，计算中以圆柱为参考系，圆柱表面的速度边界为u = v = 0，外部的边界上u = U，v = 0，初始时刻的流场均匀分布u = U，v = 0，p = 0，采用向后三阶差分进行时间离散，NURBS基函数的可导阶数k = 6，Newton迭代的收敛条件是相邻两次迭代的误差小于10−5。分别计算了雷诺数为550、3000和9500的流动，Re = 550、3000的时间步长为0.01，Re = 9500的时间步长为0.005，网格随雷诺数的增大而加密，三种情况下网格的总自由度数分别为28080、82908和127368，其中Re = 9500的网格如图2所示。
Figure 1. Domain of transient flow around an impulsively started circular cylinder
图1. 瞬时启动圆柱绕流计算区域示意图
4. 计算结果及分析
4.1. Re = 550
流动初始时刻是无旋的，由于圆柱的运动，在其表面对称地产生涡，t = 1时在离后滞止点40˚左右次涡鼓起并将主涡推离圆柱。图3是t = 2.5时流线分布情况，圆柱后存在一个较大的主涡，在离后滞止点40˚左右存在一个孤立的二次涡，受到二次涡的挤压，主涡的流线发生了变形，流线分布与文献 [15] 中实验的流线分布一致。
图4是Re = 550时的流动参数变化情况，并与文献中的数值和试验结果进行了比较。Koumoutsakos [9] 和Ploumhans [10] 采用的是涡量法，Li [11] 采用的是LBM法，不同方法的阻力系数是一致的。由于初始时刻在时间上的奇异性，阻力系数刚开始很大并快速下降，达到最小点后，随着次涡的产生和逐渐增大，圆柱前后压降增大使阻力升高，在t = 1.5左右达到最大值，此后由于次涡的强度逐渐减弱，主涡远离圆柱，阻力慢慢衰减。图4(b)中离散点为Bouard [15] 等的实验结果，可以看到最小二乘等几何方法的计算结果与实验值一致。
4.2. Re = 3000
图5是t = 2.5时的流线分布情况，与实验的流线分布一致，从图中可以看到Bouard [15] 等所谓的α-现象。图6是阻力随时间的变化情况，并与Anderson [7] 等、George [8] 等以及Koumoutsakos [9] 等的计算结果进行了对比，从图中可以看出不同方法的结果具有一致性。阻力曲线经历了初始时刻的快速下降后，在t = 0.25左右从底部逐渐回升，在t = 1.8左右达到最大值，然后逐渐降低。与Re = 550不同是在t = 1.1到1.5之间存在一个平台期，阻力在这一段时间内增加缓慢，这是因为在这一段时间内圆柱表面的压力分布变化较小，而粘性阻力已经达到稳定。图6是圆柱后对称轴线上的速度分布，与实验结果一致。
Figure 3. Streamline for Re = 550 (t = 2.5)
图3. Re = 550流线分布(t = 2.5)
(a) 阻力系数随时间变化 (b) 对称轴线上的速度分布
Figure 4. Streamline for Re = 550 (t = 2.5)
图4. Re = 550流线分布(t = 2.5)
Figure 5. Streamline for Re = 3000 (t = 2.5)
图5. Re = 3000流线分布(t = 2.5)
(a) 阻力系数随时间变化 (b) 对称轴线上的速度分布
Figure 6. Flow parameter for Re = 550
图6. Re = 3000流线分布(t = 2.5)
4.2. Re = 9500
图7中对比了计算和实验的流线分布，不同时刻下两者一致，t = 1.0时可以看到Bouard [15] 等所谓的β-现象。图8是圆柱后对称轴线上的速度分布，分别与实验结果以及George [8] 等的计算值进行了对比。与实验值之间存在一定的差异，但基本上是一致的。George等采用有限差分和投影方法，计算中采用半交错网格，共327,000个节点，总自由度数约100万，从图中可以看出本文的曲线与George等的曲线几乎完全重合，而本文的自由度数只是其13%左右，这表明本文的方法具有较高的精度。阻力系数随时间的变化曲线如图9所示，图中与Anderson [7] 等、George [8] 等以及Koumoutsakos [9] 等的计算结果进
(a) 与实验结果对比 (b) 与有限差分法结果对比
Figure 8. Velocity along the axis of symmetry for Re = 9500
图8. Re = 9500对称轴线上的速度分布
Figure 9. Variation of drag coefficient with time for Re = 9500
图9. Re = 9500阻力系数随时间变化
行了对比，不同的方法计算的结果基本一致。
5. 结论
本文应用最小二乘等几何方法模拟了静止流场中圆柱突然启动后初始阶段的流动演化过程。在圆柱后存在大小不同的涡，涡的结构随时间的演化过程和雷诺数变化而表现出不同形式，存在鼓包、孤立涡、α-现象和β-现象等流动特征，最小二乘等几何方法准确捕捉到了这些流动特征，与实验结果一致，表明该方法可以与复杂流动模拟。应用NURBS升阶算法提高基函数的光滑性后，最小二乘等几何方法具有高的精度，与有限差分法相比降低了准确模拟流动所需的自由度数。