第八章 边界条件

第八章边界条件任何数值模拟都可以认为仅仅是在物理区域或系统的一部分中进行的。

区域的断层产生了人工边界，在这个断层中有我们处理的物理量。

此外，还有暴露在流体中的自然边界。

边界条件的数值处理需要特别注意。

在实际的系统中处理不当模拟就会出现偏差。

与此同时，稳定性和求解方案中的合成速度同样对数值模拟有消极的影响。

下边边界条件的类型是我们在欧拉方程和N-S方程中数值计算最常见的几种：·固体壁面·外表面的远场和流体内部流出或流出的表面·对称面·平整切割和周期性边界。

·平板间的边界这些边界条件的处理问题在以后几节中会进行详细的介绍。

对于文献中进一步涉及的边界条件，比如壁面上的热辐射或者是自由表面上的（热辐射），读者可以在3.4节中了解。

8.1 虚拟单元的概念在我们讨论边界条件时，我们需要提到虚拟单元（也可以被称作虚拟点）这个概念。

在规则的网格中这种方法非常的流行。

然而，在不规则网格中，虚拟单元仍然有很多的优点。

虚拟单元是在物理区域外部附加层上的一些网格点。

这个可以由图8.1中的二维规则网格中看到。

正如我们看到的，整个计算区域被两层虚拟单元包围着（由虚线标出），虚拟单元（点）通常不会像区域内的网格一样产生（除过多平板的网格）。

尽管它仍然有几何形状，比如体积或者表面的矢量，但是它仅仅是虚拟的。

利用虚拟单元可以简化计算沿边界的通量，梯度，散度等等。

这是由于在边界上可以将空间离散的模型进行扩展。

正如图8.1中我们看到的，在物理区域内同样可以进行离散。

因此，我们可以在所有的“物理”网格点中求解控制方程。

这种方法可以使离散工作非常简单。

此外，所有规则的网格点可以存在在一个单独的区域内，这在矢量计算中非常很有用。

虚拟单元不但包含有守恒变量，同时也有几何量。

很明显的是，虚拟单元层必须完全覆盖物理区域外。

几何量通常由边界的控制体积来求得。

在多网格平板中（3.1节），所有的流体变量和几何变量可以从相邻的平板求得。

图8.1中灰色的虚拟单元产生了一个问题，由于不清楚怎么设置它的值（如果没有相邻网格板），通过环绕形的离散模型，我们不需要知道它的值。

但是，在计算梯度（粘性通量—见4.4节）或者在多网格中计算转移量时就很有比必要知道它的值。

通常，如图8.1中用箭头表示的一样，相邻的“规则的”虚拟单元的平均值是很有用的。

8．2 固体壁面8.2.1非粘性流动在非粘性流动中，流体流过表面，由于没有摩擦力，速度矢量与壁面相切，与就是说在壁面的法线方向没有作用力。

即：0=⋅n v 在壁面上， （8.1）n 表示壁面上的单位法相矢量。

因此，相对应的速度V 等于0。

即，对流通量的矢量简化成只有压力项：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=0 0 )(w z w y w x c p n p n p n F ω （8.2）其中Pw 为壁面压力。

规则的单元中心方案在以单元中心的方案中，在单元的质心可以求出压力。

但是，在等式（8.2）中，边界表面单元的Pw 需要求解出。

我们可以很容易的通过内部区域推导出壁面的压力。

考虑到图8.2，我们可以简单的设Pw=P2。

通过任一个两两相连点，我们可以取得较高的精确度。

)3(2132p p p w -=（8.3）或者一个三点相连的推导公式是)31015(81432p p p p w +-=（8.4）为了解释网格的延伸，到壁面的距离可以替代常数系数。

【1】上式的推导公式8.3和8.4没有涉及到网格和表面的几何形状。

一个可供选择的方法—Rizzi 发明的被称作为法向动量的关系。

它是基于在非粘性流动中，壁面是流线型的。

在式8.1中沿着流线型壁面的法线流动方向的导数为0。

并且动量方程可以替换成如下的公式：p n n v v ∇⋅=∇⋅⋅)(ρ （8.5）等式（8.5）包括密度，速度和壁面法线方向的压力。

法向动量在【1】中给出了精确的结果。

但是在复杂的几何表面，法向动量的数值求解存在着问题。

推导中详细的描述和精确的比较可以在文献【1】中找到。

虚拟单元中守恒变量的值在壁面内部能够被线性的推导出：320321232W W W W W W -=-= （8.6）公式（8.6）中的参数和图8.2中的相对应。

如果虚拟单元在空间上离散后可以被利用，那么对流通量的计算就可以和公式（8.2）一致。

规则的单元顶点的方案通过重叠的控制体积（见4.2.2），对单元顶点的离散可以直接由方程（8.1）在边界条件上应用。

壁面的对流通量可以通过公式（8.2）来求解。

壁面的压力Pw 可以通过计算平均节点的值来表示，二维公式节点的值可以通过（4.26）计算，三维公式节点的值可以通过（4.27）计算。

分类公式（二维公式中的（4.30））现在仅仅解释了两个单元（在三维中是4个）。

在重叠的控制体积内，单元顶点的值可以通过许多不同的方法求解。

一种方法是通过公式（8.2）将壁面上每个表面的控制体积分离。

因此，通过图8.3，我们可以写出：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=+++------0 )()()()()()(0 0 )()()()()()(0 )(2,4/12,2,4/12,2,4/12,2,4/12,12,4/12,12,4/12,12,,i w i z i w i y i w i x i w i z i w i y i w i x i w c p n p n p n p n p n p n F （8.7）公式（8.7）中的压力2,4/1)(-i w p 和2,4/1)(+i w p 可以通过线性的插值表示：[]2,12,2,4/1)()(341)(+++=i w i wi w p pp （8.8）对应的三维公式可以表示未划分的网格。

另外一种可能的方法是在各自的壁面节点利用公式（8.2）来进行计算。

壁面的压力Pw 可以简单的设它等于P2。

单位法向矢量可以通过节点2上所有的面的法向矢量的平均值来表示。

这种方法要求对速度矢量进行修正。

通过时间步进的求解，壁面上的速度矢量被投影到切平面。

【3】【4】2,2,2,2,2,)(])([i av i av i i corri n n v v v ⋅⋅-=）（ （8.9） av n 是平均单位法向矢量。

通过这种方法，流动将会变成与壁面相切。

为了对虚拟点进行赋值，在内部流场中，利用与公式（8.6）相关的公式推导守恒变量是非常有效的。

不规则的以单元为中心的方案公式（8.1）中壁面的边界条件也同样适用于未构建以单元为中心的方案。

如果边界单元是四面体、六面体或者是棱形（壁面上有三角形的表面），压力可以通过公式（8.3）来推导。

相邻的单元（图8.2中的3号）可以通过5.2.1中表面的数值来代替。

对于三角形或者四面体单元，在文献【5】【6】中通过一层虚拟单元来表示。

壁面上边界单元的速度矢量可以表示虚拟单元中速度分量。

比如：在图8.2中的虚拟单元1速度可以如下来表示：n V v v 2212-= （8.10）在这里z y x n w n v n u V 2222++=是推导的速度，Tz y x n n n n ],,[=表示的是壁面的单位法向矢量。

我们假设虚拟单元中的压力和密度与边界单元的值相等。

（即就是2p p w =）。

规则的双重中线的方案公式（8.1）中的边界条件需要更多的注意离散化的双重中线问题。

图8.4和图8.5分别表示了二维和三维的情况。

公式（8.2）中的对流通量可以在壁面上的每一个控制体积上的表面分开来计算。

这种方法与第一种构建单元顶点的方案相同。

对于一个四面体单元（就像图8.4中的1-3-4-5），压力可以通过公式8.8来代替。

对于六面体，棱形或者棱锥形，控制体积的表面是四边形（就像图8.5的面1-4-5-6），推导的公式可以如下来表示：)339(1615641int p p p p p +++=（8.11）如果边界元素是四面体（或者二维中的三角形），壁面的压力应该通过有限体积法来计算。

【7】在壁面的1-2部分，例如，*21-表面的压力可以通过如下的公式计算：)5(6121int p p p +=（8.12）对于四面体，例如图8.5中的壁面1-2-3，压力可以如下表示：)6(81321int p p p p ++=（8.13）8.2.2 粘性流动对于经过固体壁面的粘性流体，介于表面和流体的中垂直于表面的为零。

因此，我们我们可以称之为非滑动的边界条件。

在一个固定的壁面，速度分量在表面变为：0===w v u （8.14）对于非滑动的边界条件，有两个基本的结果。

第一点，我们不需要解壁面上的动量方程，这个在单元顶点方案中已经应用。

第二点，对流通量通过非滑动的壁面的公式已经在公式（8.2）中给出，并且在公式（2.24）中的项已经简化成T k ∇=Θ，因此，对流通量中的壁面压力同样可以通过非粘性流动来求得。

但是，虚拟的单元处理的方法不同。

以单元为中心的方案公式（8.14）中非滑动边界条件的处理可以通过利用虚拟单元来简化。

在一个绝热的壁面（没有热通量通过壁面），我们可以取（如图8.2）2121212121 w w v v u u E E -=-=-===，，，ρρ （8.15）同样的对于单元0和3。

方法同样适用于构建的以单元为中心的方案的和非构建的以单元为中心的方案（参考文献【6】）。

如果壁面的温度已知，速度分量仍然像方程（8.15）中一样是相反的。

利用壁面的温度可以将周围的温度线性的表示。

由于压力的梯度垂直于壁面并且为零，边界元素的压力同样可以用虚拟单元来表示（例如210 p p p ==）。

虚拟单元中的密度和总能量可以通过推导的值来计算。

单元顶点的方案由于动量方程不需要求解，壁面上的对流通量没有作用。

公式（2.23）中的粘性通量仅仅对能量方程中垂直于壁面的温度梯度有作用。

对于一个绝热的壁面，n T w ⋅∇为零。

因此，我们没有必要求解任何壁面上的对流或者粘性通量。

为了防止避免节点上非零速度分量的产生，动量方程中剩余的项都为零。

在已知壁面温度的情况下，我们可以直接设壁面（例如图8.3中的节点)2,(i ）的总能量为（假设为理想气体）w i pi T c E 2,2,)(ργρ=(8.16)在这里，w T 表示给定的壁面温度。

剩余的动量和能量方程都设为零。

对于未建立一单元为中心的方案同样适用。

另外一种方法，对于一些应用来说比较简略，并且根本没有解决壁面的控制方程。

同样的，密度和能量可以直接简化为RT p w i i 3,2,=ρ 和γρ3,2,)(i i p E =（8.17）方程（8.17）中假设垂直于壁面没有压力梯度（因此32p p =）。

由于所有的守恒变量都是一定的，剩余的所有的方程都等于零。

这个方法同样适用于未构建的网格。