搜索
一、求解N-S方程的数值方法求解N-S方程的数值方法可分为有网格和无网格两类。有网格方法的求解思路是将连续空间划分为离散网格(或单元)，并在网格点(或单元)上离散控制方程，从而寻求控制方程的数值解；无网格方法则将着眼点回归到流体运动的基本单元质点上，通过追踪质点运动性质的变化过程来获得流体解答。网格数值模拟方法发展较早，求解方法较成熟，应用也较广；而无网格方法则刚刚起步，具有很大的发展空间。下面分别就两类中较有代表性的方法予以介绍。1.1网格数值模拟方法数值计算网格可分为不随流体运动而改变形状的欧拉网格和随流体运动改变形状的拉格朗日网格两类。欧拉网格划分较简单，但在用于具有自由表面流体运动数值计算中需要对自由表面位置进行适时追踪；拉格朗日网格虽然无需进行自由表面追踪，但需要在计算过程中随时调整计算网格，以适应流体自由表面形状，额外增加了计算量，并且在流体变形剧烈处还会出现网格失效现象。有限差分法(FiniteDifferenceMethod,FDM)和边界元法(BoundaryElementMethod,BEM)都是欧拉网格方法中常用的方法。在拉格朗日网格法中，有限元法和边界元法可直接应用，而有限差分法一般需对坐标做适体变换。用欧拉网格法求解含有自由表面的水波动力问题难点在于追踪自由表面。目前比较有代表性的自由表面追踪方法有PIC(Particle-In-Cell)法、、MAC(Marker-And-Cell)法和VOF(Volume－Of－Fluid)法等。以下分别予以介绍。1.1.1质点网格法(Particle-in-Cell，简称PIC法)PIC法是Harlow提出的一种用欧拉矩形网格计算多种介质流体运动的方法。该方法把流体既视为连续介质，又视为带有一定质量的质点，然后研究质点在经过固定欧拉网格上的运动性质。PIC法具有计算多相流和处理三维自由表面的能力，曾经成功运用于模拟二维流体发生剧烈变形的情况。其缺点是占用较多计算机内存；模式只能给出自由表面单元位置，而不能给出自由表面精确位置；另外在流体变形剧烈处易出现较大插值误差。1.1.2MAC法(Marker-And-Cell)针对PIC法的上述缺点，20世纪60年代中期，Harlow,Welch提出了一种改进方法，称为MAC法。最初的MAC法在所有单元内部都布满无质量的标记点，通过跟踪这些点即可判断自由表面的位置。MAC法在求解时以速度分量和压力为因变量，在固定网格上采用有限差分法离散控制方程。在MAC法之后，人们又先后提出了SMAC法、ABAMC法、SUMAC法和TUMMAC法［11］系列。新版本模型改为只在自由表面处设标记点，对这些点的连续跟踪可得到自由表面的准确位置，且所耗机时大为减少。Hirt和Nichlos在总结当时的自由表面跟踪方法的基础上还提出了“线段法”(LineSegmentMethod)的概念，可近似确定自由表面的位置，但这种方法推广到三维的情况较困难。MAC法于20世纪60年代提出后先是用于涌潮的传播研究，后自80年代始大量用于波浪研究。如，Miyata运用一种二维MAC格式研究了波浪破碎；Sakai等也提出一种基于早期MAC法的数值模型，来模拟水波泼溅以及波浪破碎过程中二阶乃至三阶漩涡的生成。Gao和Zhao运用二维MAC模型，研究了波浪、建筑物及沙滩的交互作用。在MAC法之后，学者们开发了更具通用性和计算效率的模型。如，Chorin提出的投影法(ProjectionMethod)就是以其数学上的精确、简明而著称的一个。以MAC法为基础，Hirt等和Nichols等开发了易于操作的SOLA(solution-algorithm)模型。Kothe等利用不完全Cholesky共轭梯度法，将表面张力模拟为一种体积力，开发出了RIPPLE模型。1.1.3VOF法（Volume-Of-Fluid）1.1.3.1算法思想Nichols等是最早用利用VOF法的研究群体。VOF法假设在整个流动区域内流体密度为一常数0，空气密度设为0，引进流体体积参数(或者是VOF函数)0/F,从连续方程出发建立关于F的对流方程：.0FuFtF在流体单元内为1，在空单元中则为0，而在自由表面单元中，F值介于0和1之间。给定速度场就可以通过上述对流方程确定F值，进而确定任一时刻的自由表面单元位置。1.1.3.2方法介绍供体－受体法(Donor-AcceptorMethod)：由于VOF函数在自由表面附近(1-1)有很大梯度，用一般的激波捕捉方法虽可求解上述VOF对流方程，但由于存在严重数值耗散，难于精确确定自由表面位置。为此，Hirt和Nichols提出了供体－受体法(Donor-AcceptorMethod)。该法首先依据流速在x,y方向的梯度大小确定自由水面的走向，并根据流速的方向判断自由水面附近的单元格为供体格还是受体格，并进一步在运动方程求解中决定采用逆风或顺风格式。这就有效地控制了数值耗散。有关数值实现的详细信息可参阅Hirt和Nichlos、Kothe、Liu和Lin等相关文献。倾斜界面法(SlopingInterfaceMethod)：上述供体－受体方法对自由表面几何形状进行了简化处理，导致供体－受体法精度不高。Youngs引入了倾斜界面法，较之于供体－受体法，此类方法可更精确确定自由表面的位置。但目前在算法实现以及向三维扩展时有一定困难，也就妨碍了倾斜界面法的实际应用。层面设置法(LevelSetMethod)：该法由Sussman等于1994年提出，是早期借助移动网格跟踪界面方法的改进版本。该方法很容易确定分界面并且数值耗散较小。但由于这种格式没有明确利用质量守恒方程，因此复杂情况下容易造成质量亏损。1.1.3.3应用情况Delft水力学研究小组自20世纪90年代早期就利用VOF方法进行波浪动力学研究。VanderMeer等提出一种模型，研究作用在沿岸建筑物上的卷破波。Iwata等研究了由于水下建筑物所致的波浪破碎变形。Lin和Liuk-e利用紊流模型封闭N-S方程，开发了一种基于VOF法的模型，对倾斜海滩上孤立波破碎等诸多水波动问题进行了研究。Sussman等耦合VOF法和层面设置法，用于三维情况的不可压缩流体的研究。上述方法虽然解决了自由表面的追踪问题，但当自由表面有剧烈变形时，各类方法在处理对流项时，数值耗散问题仍然都是不可避免的。1.2无网格数值模拟方法无网格数值方法也称为粒子方法，该方法完全摒弃了数值网格的概念，由于运动方程中不直接出现对流项，自然也就不会出现网格方法中的数值耗散现象。目前文献上报道的粒子方法有PAF法(Particle-And-Force)、无网格Euler/N-S方程解法、自由元(FreeElement)Galekin法、离散元(DiscreteElement)法、SPH(Smoothed-Particle-Hydrodynamics)法以及MPS法(MovingParticleSemi-implicitMethod)等。在水动力学研究中开展较多的方法主要有MPS和SPH两种，以下对其予以介绍。1.2.1SPH方法SPH方法的基础是插值理论。SPH通过某个插值函数在质点上给出变量的近似值，从而将以偏微分方程形式给出的连续流体的守恒方程转换成积分方程。控制方程中的压力梯度项、扩散项等空间导数项也通过核子函数表示。Cummins和Rudman提出了一种SPH投影法(SPHprojectionmethod)，通过求解压力泊松方程确定压力，从而确立了具有严格算法的不可压缩模型。SPH方法源于天体物理学，后扩展至水动力学研究。随后，Edmond等用不可压缩SPH模型求解Lagrange形式的N-S方程，并将模型与大涡数值模拟相耦合，以模拟近岸孤立波的作用机理。继而，Shao等将不可压缩SPH模型用于模拟带有自由表面的牛顿和非牛顿流体的研究。1.2.2MPS方法Koshizuka等提出了一种粒子方法，用于计算不可压缩黏性流体，并将之命名为移动粒子半隐式方法(Moving-ParticleSemi-implicitmethod)，简称为MPS法。该方法中，粒子间相互作用通过一个核子函数来表征，梯度和扩散算子运算都通过核子函数表达。流体的不可压缩性是通过控制粒子密度数不变实现的。计算过程分为显式(ExplicitStage)和半隐式两个阶段。显式阶段仅考虑重力和黏滞项，半隐式阶段用隐式方法求解压力泊松方程，但速度的二次校正仍然为显式。目前，该方法已用于波浪破碎机理和两相流的研究中。另外，Gotoh还将亚格子尺度紊流模型引入MPS模型，获得了N-S方程的Lagrange解。总体来看，水波动力学中的无网格数值方法还处于初期探索阶段，运用于实际工程中的例子还少见报道。但由于其独有的可很方便地追踪自由表面的特性，决定了其在某些特定的工程实际问题中有网格方法无可比拟的优越性，因而也具有很好的发展前景。二、湍流及其数值模拟方法2.1湍流湍流是流体的一种流动状态。当流速很小时，流体分层流动，互不混合，称为层流，也称为稳流或片流；逐渐增加流速，流体的流线开始出现波浪状的摆动，摆动的频率及振幅随流速的增加而增加，此种流况称为过渡流；当流速增加到很大时，流线不再清楚可辨，流场中有许多小漩涡，层流被破坏，相邻流层间不但有滑动，还有混合。这时的流体作不规则运动，有垂直于流管轴线方向的分速度产生，这种运动称为湍流，又称为乱流、扰流或紊流。2.2湍流的数值模拟工程中绝大多数流体都是湍流。深层次地揭示流体运动的紊动结构及能量耗散过程是计算流体力学应用研究的极其重要的课题。出于不同的研究目的，湍流数值模拟有以下三个不同的层次：直接数值模拟(DirectNumericalSimulation，简称DNS)、雷诺平均数值模拟(ReynoldsAveragedNavier-Stokes，简称RANS)和大涡数值模拟(LargeEddySimulation，简称LES)。以下分别简要介绍。2.2.1直接数值模拟由于N-S方程本身是封闭的，故从原则上讲可以求解所有湍流问题。用DNS法直接求解N-S方程，能获得最精细的流场信息。但在目前计算机发展水平下，进行DNS的应用研究是不现实的，目前只限于小规模的低雷诺数简单湍流物理机制的研究。DNS常用的数值方法是谱方法或伪谱法。Orszag等最早用DNS计算了各向同性湍流。2.2.2雷诺平均数值模拟由于雷诺平均方程是不封闭的，故RANS的核心思想是建立雷诺应力封闭模型，使得平均运动方程可解。目前比较常用的模型有零方程模型、一方程模型(k方程模型)、二方程模型(涡黏性模型，k-e模型)、代数应力模型(k-e-A模型)。其中k-e模型是目前应用最广泛的湍流模型。在可以预见的将来，即使DNS求解复杂湍流成为现实(耗时很长)，能够快速地算出满足一定工程精度要求的湍流统计模式仍会受到工程师们的青睐。事实上，RANS是目前工程界处理湍流问题的唯一方法。2.2.3大涡数值模拟LES的基本思想是：把湍流瞬时运动通过某种滤波方法分解为大尺度运动和小尺度运动两部分。大涡运动通过直接求解N-S方程计算；小涡运动的影响概化为亚格子雷诺应力，需通过建立模型求解，模型称为亚格子尺度模型(SubgridScaleModel)。气象学家Deardoff首次把大涡模拟用于有工程意义的槽道中的流体运动数值计算中。1972年起，Stanford大学的Ferziger和Reynolds领导的集体开始对LES做深入系统的研究。苏铭德1982年以来提出了一种代数应力模型，计算了槽道中的流体流动。Watanabe用大涡模拟计算了三维波浪破碎。LES至今在气体动力学中的研究开展较多，在水动力学中的研究还十分有限，最主要的困难还不在于计算机的限制，而在其方法本身，如现有的亚格子尺度模型仍很不完善以及近壁模型的入流出流边界等问题。三、水力机械空化数值模拟问题描述：高速离心泵在工作时，容易在叶轮叶片处产生汽蚀，会严重影响离心泵的性能。在叶轮入口处配置诱导轮,可以提高离心