第八章 边界条件.doc

1、第八章 边界条件任何数值模拟都可以认为仅仅是在物理区域或系统的一部分中进行的。区域的断层产生了人工边界，在这个断层中有我们处理的物理量。此外，还有暴露在流体中的自然边界。边界条件的数值处理需要特别注意。在实际的系统中处理不当模拟就会出现偏差。与此同时，稳定性和求解方案中的合成速度同样对数值模拟有消极的影响。下边边界条件的类型是我们在欧拉方程和 N-S 方程中数值计算最常见的几种：固体壁面外表面的远场和流体内部流出或流出的表面对称面平整切割和周期性边界。平板间的边界这些边界条件的处理问题在以后几节中会进行详细的介绍。对于文献中进一步涉及的边界条件，比如壁面上的热辐射或者是自由表面上的（热辐射）

2、，读者可以在 3.4 节中了解。8.1 虚拟单元的概念在我们讨论边界条件时，我们需要提到虚拟单元（也可以被称作虚拟点）这个概念。在规则的网格中这种方法非常的流行。然而，在不规则网格中，虚拟单元仍然有很多的优点。虚拟单元是在物理区域外部附加层上的一些网格点。这个可以由图 8.1 中的二维规则网格中看到。正如我们看到的，整个计算区域被两层虚拟单元包围着（由虚线标出） ，虚拟单元（点）通常不会像区域内的网格一样产生（除过多平板的网格） 。尽管它仍然有几何形状，比如体积或者表面的矢量，但是它仅仅是虚拟的。利用虚拟单元可以简化计算沿边界的通量，梯度，散度等等。这是由于在边界上可以将空间离散的模型进行扩展

3、。正如图 8.1 中我们看到的，在物理区域内同样可以进行离散。因此，我们可以在所有的“物理”网格点中求解控制方程。这种方法可以使离散工作非常简单。此外，所有规则的网格点可以存在在一个单独的区域内，这在矢量计算中非常很有用。虚拟单元不但包含有守恒变量，同时也有几何量。很明显的是，虚拟单元层必须完全覆盖物理区域外。几何量通常由边界的控制体积来求得。在多网格平板中（3.1 节） ，所有的流体变量和几何变量可以从相邻的平板求得。图 8.1 中灰色的虚拟单元产生了一个问题，由于不清楚怎么设置它的值（如果没有相邻网格板） ，通过环绕形的离散模型，我们不需要知道它的值。但是，在计算梯度（粘性通量见 4.4

4、节）或者在多网格中计算转移量时就很有比必要知道它的值。通常，如图 8.1 中用箭头表示的一样，相邻的“规则的”虚拟单元的平均值是很有用的。82 固体壁面8.2.1 非粘性流动在非粘性流动中，流体流过表面，由于没有摩擦力，速度矢量与壁面相切，与就是说在壁面的法线方向没有作用力。即：在壁面上， （8.1）0nvn 表示壁面上的单位法相矢量。因此，相对应的速度 V 等于 0。即，对流通量的矢量简化成只有压力项：（8.2）0 )(wzyxcpnF其中 Pw 为壁面压力。规则的单元中心方案在以单元中心的方案中，在单元的质心可以求出压力。但是，在等式（8.2）中，边界表面单元的 Pw 需要求解出。我们可以

5、很容易的通过内部区域推导出壁面的压力。考虑到图 8.2，我们可以简单的设 Pw=P2。通过任一个两两相连点，我们可以取得较高的精确度。（8.3）)3(21ppw或者一个三点相连的推导公式是（8.4）)05(8432为了解释网格的延伸，到壁面的距离可以替代常数系数。 【1】上式的推导公式 8.3 和 8.4 没有涉及到网格和表面的几何形状。一个可供选择的方法Rizzi 发明的被称作为法向动量的关系。它是基于在非粘性流动中，壁面是流线型的。在式8.1 中沿着流线型壁面的法线流动方向的导数为 0。并且动量方程可以替换成如下的公式：（8.5）pnv)(等式（8.5）包括密度，速度和壁面法线方向的压力。

6、法向动量在【1】中给出了精确的结果。但是在复杂的几何表面，法向动量的数值求解存在着问题。推导中详细的描述和精确的比较可以在文献【1】中找到。虚拟单元中守恒变量的值在壁面内部能够被线性的推导出：（8.6）3201W公式（8.6）中的参数和图 8.2 中的相对应。如果虚拟单元在空间上离散后可以被利用，那么对流通量的计算就可以和公式（8.2）一致。规则的单元顶点的方案通过重叠的控制体积（见 4.2.2） ，对单元顶点的离散可以直接由方程（8.1）在边界条件上应用。壁面的对流通量可以通过公式（8.2）来求解。壁面的压力 Pw 可以通过计算平均节点的值来表示，二维公式节点的值可以通过（4.26）计算，三

7、维公式节点的值可以通过（4.27）计算。分类公式（二维公式中的（4.30） ）现在仅仅解释了两个单元（在三维中是 4 个） 。在重叠的控制体积内，单元顶点的值可以通过许多不同的方法求解。一种方法是通过公式（8.2）将壁面上每个表面的控制体积分离。因此，通过图 8.3，我们可以写出：（8.7）0 )(0 0 )(0 )( 2,4/12,2,4/12,2,4/12,1,2,4/12,12, iwizyiixiwizyiixiwc pnpnF公式（8.7）中的压力 和 可以通过线性的插值表示：2,4/1)(iwp2,4/1)(iw（8.8）2,12,/ )()(3iwii p对应的三维公式可以表示未

8、划分的网格。另外一种可能的方法是在各自的壁面节点利用公式（8.2）来进行计算。壁面的压力Pw 可以简单的设它等于 P2。单位法向矢量可以通过节点 2 上所有的面的法向矢量的平均值来表示。这种方法要求对速度矢量进行修正。通过时间步进的求解，壁面上的速度矢量被投影到切平面。 【3】 【4】（8.9）2,2,2,2, )(iaviaviicori nv）（是平均单位法向矢量。通过这种方法，流动将会变成与壁面相切。avn为了对虚拟点进行赋值，在内部流场中，利用与公式（8.6）相关的公式推导守恒变量是非常有效的。不规则的以单元为中心的方案公式（8.1）中壁面的边界条件也同样适用于未构建以单元为中心的方案

9、。如果边界单元是四面体、六面体或者是棱形（壁面上有三角形的表面） ，压力可以通过公式（8.3）来推导。相邻的单元（图 8.2 中的 3 号）可以通过 5.2.1 中表面的数值来代替。对于三角形或者四面体单元，在文献【5】 【6】中通过一层虚拟单元来表示。壁面上边界单元的速度矢量可以表示虚拟单元中速度分量。比如：在图 8.2 中的虚拟单元 1 速度可以如下来表示：（8.10）nVv21在这里 是推导的速度， 表示的是壁面的单位法zyxnwvnuV22 Tzyx,向矢量。我们假设虚拟单元中的压力和密度与边界单元的值相等。 （即就是 ） 。2pw规则的双重中线的方案公式（8.1）中的边界条件需要更多

10、的注意离散化的双重中线问题。图 8.4 和图 8.5 分别表示了二维和三维的情况。公式（8.2）中的对流通量可以在壁面上的每一个控制体积上的表面分开来计算。这种方法与第一种构建单元顶点的方案相同。对于一个四面体单元（就像图 8.4 中的 1-3-4-5） ，压力可以通过公式 8.8 来代替。对于六面体，棱形或者棱锥形，控制体积的表面是四边形（就像图 8.5 的面 1-4-5-6） ，推导的公式可以如下来表示：（8.11）)39(16564int pp如果边界元素是四面体（或者二维中的三角形） ，壁面的压力应该通过有限体积法来计算。 【7】在壁面的 1-2 部分，例如， 表面的压力可以通过如下的

11、公式计算：*2（8.12）)(21int对于四面体，例如图 8.5 中的壁面 1-2-3，压力可以如下表示：（8.13）)6(8321int pp8.2.2 粘性流动对于经过固体壁面的粘性流体，介于表面和流体的中垂直于表面的为零。因此，我们我们可以称之为非滑动的边界条件。在一个固定的壁面，速度分量在表面变为：（8.14）0wvu对于非滑动的边界条件，有两个基本的结果。第一点，我们不需要解壁面上的动量方程，这个在单元顶点方案中已经应用。第二点，对流通量通过非滑动的壁面的公式已经在公式（8.2）中给出，并且在公式（2.24）中的项已经简化成 ，因此，Tk对流通量中的壁面压力同样可以通过非粘性流动来

12、求得。但是，虚拟的单元处理的方法不同。以单元为中心的方案公式（8.14）中非滑动边界条件的处理可以通过利用虚拟单元来简化。在一个绝热的壁面（没有热通量通过壁面） ，我们可以取（如图 8.2）（8.15）212121 wvuE，同样的对于单元 0 和 3。方法同样适用于构建的以单元为中心的方案的和非构建的以单元为中心的方案（参考文献【6】 ） 。如果壁面的温度已知，速度分量仍然像方程（8.15）中一样是相反的。利用壁面的温度可以将周围的温度线性的表示。由于压力的梯度垂直于壁面并且为零，边界元素的压力同样可以用虚拟单元来表示（例如 ） 。虚拟单元中的密度和总210 p能量可以通过推导的值来计算。单

13、元顶点的方案由于动量方程不需要求解，壁面上的对流通量没有作用。公式（2.23）中的粘性通量仅仅对能量方程中垂直于壁面的温度梯度有作用。对于一个绝热的壁面，为零。因此，我们没有必要求解任何壁面上的对流或者粘性通量。为了防止避nTw免节点上非零速度分量的产生，动量方程中剩余的项都为零。在已知壁面温度的情况下，我们可以直接设壁面（例如图 8.3 中的节点 ）的)2,(i总能量为（假设为理想气体）(8.16)wipiTcE2,2,)(在这里， 表示给定的壁面温度。剩余的动量和能量方程都设为零。对于未建立wT一单元为中心的方案同样适用。另外一种方法，对于一些应用来说比较简略，并且根本没有解决壁面的控制方

14、程。同样的，密度和能量可以直接简化为和 （8.17）RTpwii3,2,3,2,)(iipE方程（8.17）中假设垂直于壁面没有压力梯度（因此 ） 。由于所有的守恒32p变量都是一定的，剩余的所有的方程都等于零。这个方法同样适用于未构建的网格。但是，压力的推导在三角形或四边形网格中需要一些附加条件。如果壁面是绝缘的，虚拟单元中的值如下所示：（8.18） 3,1,3,1,3,1, , iiiiii wvuE，以上同样适用于节点 0 和 4.如果壁面的温度已知，虚拟单元中的温度可以从内部推导出：（8.19）3,0,3,1, 2 2iwiiwi TT和按照图 8.3 所示，速度分量和公式（8.18）

15、中的同样是相反的。密度和能量可以通过内部的温度值和压力 来算出。3,ip8.3 远场外部流体通过机翼，翅膀，汽车和其它另外的结构数值模拟通常是在一个有限的区域内传导。正是如此，人为加上的无限远的边界条件就很有必要。无限远边界条件的数值表示通常要满足两个基本的要求。第一，与无限的区域相比较，无限远处的横截面应该对流体的求解没有显著的影响。第二，任何流出的扰动不能对原来的流场有影响。由于它们具有抛物线的特点，低于或接近于音速的流体通常对远场边界条件比较敏感。处理不当的话就会降低稳定状态流动下的集中。此外，解的精确度也会受到影响。在人为规定的边界上求解吸收流出波浪的方法有很多。 【10】【15】 。

16、我们在文献【16】中可以找到以前所学的非反映条件的边界。在以下两个部分中，我们将讨论由 Whitfield 和 Janus 提出的特征变量的概念，我们将在远场边界条件下讨论升力体。8.3.1 特征变量的概念通过雅克比对流通量的特征值，在特征变量的条件下，信息通常在计算区域内外进行传递。比如，在亚音速的来流条件下有四个流进的特性（在三维情况下）和一个流出的特性。对于亚音速去流时情况正好相反。根据 Kress 的一维理论，边界外边的特性应该等于内部的内部的特性。剩余的条件应该决定于区域内部的解。Whitfield 和 Janus12的方法是基于垂直于壁面的一维欧拉方程的特性。这种方法可以很好的处理

17、大量已构建的网格和未构建的网格问题。它不但适用于远场条件，也同样适用于非粘性固体壁面。在远场条件下两个基本的流动情形可以在图 8.6 中看到。流体可以自由的进出所限定的区域。因此，根据当地的马赫数，可以得到四种不同的远场条件：1.超音速来流2.超音速去流3.亚音速来流4.亚音速去流超音速来流对于超音速来流，所有的特征值都有相同的符号。由于流体进入一个物理区域，边界上的守恒变量仅仅取决于自由液面。因此，(8.20)abW的值主要与已知的马赫数 和两个流体的角度有关。 （攻角和边切角）aWM超音速去流在这种情况下，所有的特征值同样有相同的符号。但是，流体流出物理区域并且所有的守恒变量决定于区域内的

18、解。这里可以设：(8.21)dbW亚音速来流这里，物理区域内有四个特征变量区域外有一个。因此，根据自由液面的值来计算四个特征变量。区域外的一个特征变量可以通过物理区域的内部四个变量进行推导。边界条件如下所示：(8.22)/(/)( )()()(210020cpnwvu wnvuncpbazabybaxab dazdaydaxdab 这里， 代表初始状态，初始状态通常等于内部的点（如图 8.6 中的 d 点） 。A0c和点的值取决于自由液面的状态。亚音速去流在亚音速去流的情况下，四个流体变量（密度和三个速度分量）可以通过物理区域内部进行推导。剩余的第五个变量（压力）可以通过外部得出。远场边界的初

19、始变量可以通过文献【12】得到：(8.23)/(/)(002cpnwvupbdzdbybdxdbabPa 是静态压力。虚拟单元中的物理特性可以从状态 b 和 d 线性的推导。8.3.2 升力体的修正上面在远场边界条件下假定是没有循环的，这个对于升力体在亚音速和接近声速的情况下是不正确的。正式由于这样的原因，远场边界将不得不远离物体。否则，流体的求解将不太准确。如果自由液面的流动包括一个单独的点的影响（三维中的马蹄形的顶点） ，那么到远场的距离就会大大的缩短（一个数量级） 。这个顶点可以设为升力体的中心。接下来，我们将会在二维和三维中对顶点进行修正。二维中顶点的修正我们在这里讲的方法是参考的文献

20、【18】 。自由液面的速度修正分量可以用以下的表达式表示（假设流体可以压缩）：（8.24）cos)(sin12i)(i2* 22* Mdvu是环量， 是极坐标中远场的点 是攻角， 是当地马赫数，特别的，),(环量可以有下式求得：(8.25)LaCv21通过利用库塔条件，在方程（8.25） ， 表示机翼的弦长， 是升力系数，可以通过L表面压力求得。公式（8.24）中的极坐标可以表示如下：(8.26)refrefrefxydtan)()(22这里 和 是坐标的初始点（落在顶点上，比如在弦的 1/4 处） 。refxrfy修正的自由液面的压力表示为：(8.27)1/(/12*/)1(* pvp其中

21、，修正的自由液面的密度可以通过状态方程来表示：2*2*vuv(8.28)/1*p修正量 带替公式（8.22）或者（8.23）中的 。*,和pvu aapvu和,以上的顶点修正公式（8.24）（8.25）只有在亚音速的条件下才能使用。但是，对于自由液面条件下的修正同样适用于接近音速的情况下。这种情况在图 8.7 中可以说明。在图中，升力系数到远场边界的距离是已知的。远场的半径可以设为 5%，20% ，50%和 99%的弦长。正如我们看到的，没有用顶点修正的仿真和远场的距离相对独立。相反的是，在距离为 20 时，计算的结果比较的精确。这样导致了网格单元/点的减少。这一点在文献【19】中利用高阶项来

22、进行顶点修正来说明，远场边界能够仅仅允许 5%弦长的精度。三维中的顶点修正机翼中关于远场边界的影响可以近似的认为是马蹄形的顶点。修订的自由液面的速度分量可以从文献【20】 【21】得到：(8.29) CylzBylzwAxlzlzvAu22* 222*2 )()(2其中 表示环量， 表示笛卡尔坐标系中的远场点，l 表示一半的跨度。此外，),(yx在公式（8.29）中，假设流体在 x 轴的正方向，机翼沿着 z 轴放置。公式（8.29）中的A，B 和 C 如下所示：(8.30)xCBlzA1简写成如下的式子：(8.31)2221Mylzx其中 是当地的马赫数。环量可以通过式（8.25）来计算。A

23、表示实际的弦线。修M正值压力和密度可以分别从公式（8.27）和（8.28）来求得。公式（8.22）或（8.23）中的量 可以被它们的修正值代替。公式（8.29）中修正的速度分量的表达式aapwvu和,变成了无限的，其中顶点的线通过流出的边界。这些点可以表示为： ,0 , ,0 , ylzylz并且 。为了避免出现数值奇点，在文献【21】中对值进行了限制farx 22)()(ylzylz和在公式（8.29）中是对于四分之一的翼展，例如 ，这种测量在沿着顶点线/时降低了距离 的修正速度 的精度。lzl 和 /l wv 和文献【21】中的数值结果表示的是按照公式（8.29）中提供的方法，升力和拉力系

24、数精度都会降低。通过实验发现距离远场边界 7*l 时两个系数会得到精确的值。8.4 入口和出口的边界基于 N-S 方程下有很多的方法处理出口和入口的边界条件。在这里，我们主要关注在涡轮机械应用中发展起来的方法。适合于非反映的进口和入口边界条件在文献【27】-【30】中有所描述。Giles 的文章【31】和 Hirsch 和 Verhoff32提出了欧拉方程形式的非反映边界条件，这个边界条件是介于主体和进口或者出口平面的一段小小的距离。在确定的条件下，进口和出口的边界是周期性的，它与速度、压力以及温度的梯度有关。这种类型的流体偶尔可以碰到，比如，在热交换的仿真中。周期性的进口和出口的边界条件的问

25、题在文献【34】- 【36】里关于不同频道的 LES（发射脱离系统）中可以找到。亚音速的进口一个普通的公式包括总压力，总温度和两个流动角。一个特征变量从流体区域内可以推导出来。一个可以通过黎曼不变式解出，其定义如下：(8.32)12ddcnvR参数 d 表示区域内部的状态（如图 8.6a） 。黎曼不变式用来决定是绝对速度还是边界上的速度。通过实验可以发现，选择声音的速度可以得到一个较稳定结果，特别是在低马赫数下，因此，我们设(8.33) 21)(1cos(cos12cos)( 20 RRb是与边界相关的流动角， 是声音的滞止速度，因此：0(8.34)2cosdvn并且(8.35)2201ddc

26、这里 表示的是在内部点 d 的合速度（图 8.6） 。公式（8.34）中的单位法向矢量 垂2dv n直区域向外。一些静态的值比如温度，压力，密度或者边界上的绝对速度如下所示：(8.36)(20)1/(2020bpdbbbbbTcvRTpc和 是已知的总的温度和压力， 和 分别代表特定的气体常数和在常压下的热0Tp系数。进口的速度分量可以将 沿着两个流体角度进行分解。 （在二维中是一个）2dv亚音速的出口在涡轮机械中，静态压力通常用来描述出口。亚音速下的出口边界可以由公式（8.23）中流出的情况做相同的处理。仅仅是环境压力 Pa 换成了静态出口压力。在虚拟单元中的变量可以通过边界线性推导以及内部

27、的 d 点。8.5 对称平面如果流体是线性或者平面对称，第一种情况下肯定是没有通量穿过边界。也就是说垂直于对称边界的速度为零。此外，下面的梯度消失：垂直于边界的标量的梯度垂直于边界的切向速度的梯度沿着边界的法向速度的梯度（由于 ）0nv我们可以把上面的条件写成如下式所示：（8.37）0)(nvttU式中，U 表示的是标量的变量， 表示的切于对称边界的矢量。以单元为中心的方案对称边界条件的处理可以通过虚拟单元来进行大大的简化。虚拟单元中的流体变量可以通过反映的单元来表示。这个就意味着虚拟单元中例如密度和压力的标量和相反的内部单元是相同的。例如：（8.38）3021 U和对应的符号参考图 8.2。

28、速度分量如公式（8.10）中的边界上表示的一样。虚拟单元中的法向速度的法向梯度与对称的单元的法向梯度相等，但是符号相反。单元顶点的方案（双重控制体）两种不同的方法可以解决这个问题，一种可能性是通过反射边界的网格来建立缺失的半个控制体积。通量和梯度通过在内部利用反映出的流体变量来计算（同上） 。第二种方法是计算半个控制体积的通量（但是不是沿着边界） 。剩余的垂直于对称平面的分量都用零来表示。同样对于贴于边界的控制体积其表面的法向矢量可以进行修正。 （如图 8.4 中的 点）*2修正包括去点所有表面矢量的分量，其分量垂直于对称平面。梯度同样可以由公式（8.37）进行修正。8.6 坐标切割这种情况仅

29、仅是在建成的网格中遇到。坐标切割表示的是人为的，非物理的边界，它是一条加在网格点上的线（三维中是个平面）它有不同的计算坐标，但是有相同的物理方位，这就意味着这个网格与自己重叠。正如我们将在 11.1.1 中看到的，坐标切割在被称作拓扑网格中出现。流体变量和它们的梯度将通过切割连续。OFigC或 者)5.1(最好的处理切割边界条件的方法是利用虚拟单元。这种情况在图 8.8 进行了描述。正如我们看到的，这里的虚拟层不是虚拟的，但是它们与切面的反面的网格是一致的，或者在虚拟的点可以直接从反向的单元得到。在以单元为中心的方案中，通过边界单元表面的通量可以像内部流场一样精确的计算。切面边界由两种不同的以

30、单元为中心的方案来解决。一种是在切面产生一个完整的控制体积。 （第二部分可以通过图 8.8 的虚线来表示） 。利用虚拟点，通量可以当作区域内的点进行求解。第二种方法是组合每一个分开的半个控制体积。在图 8.8b 点 2 和点 5 剩余的控制体是加上的。图中点 2 和点 5 中很好的组合起来非常重要。8.7 周期性边界有一些关于流场对一个或者多个坐标方向是周期性变化的实际应用。在这种情况下，只要仅仅模拟出重复地区中的一个就可以很好的解决问题。剩余的物理区域通过修正使其符合周期性的边界条件。我们可以区分两种基本的周期性的边界。第一个是周期性的移动。这个意味着通过坐标的移动一个周期性的边界可以转移到

31、另一个边界。第二种类型是周期性的边界，它是由于坐标的旋转。因此我们可以称之为旋转的周期性。在下面，我们将会求解以单元为中心的周期性的边界条件以及和单元顶点方案的周期性边界条件。我们同样需要考虑旋转的周期性。进一步关于周期性边界的细节可以在文献【37】 【38】中找到。以单元为中心的方案利用虚拟单元的概念可以对周期性边界条件进行简化。我们可以从表示涡轮机械的图8.9 为例。在垂直方向，涡轮机械的结构是周期性的。阴影中的单元 1 和 2 分别放在较低和较高的周期性边界。由于周期性的原因，第一个虚拟单元层与相反的周期边界的边界单元相一致。第二层虚拟单元与第二层的物理单元相一致等等等等。因此，所有在虚

32、拟单元中的标量值（密度，压力）都对应了相应的物理单元。例如：(8.39) 21 U和矢量值（速度，梯度）同样在移动性的周期性变化中具有相同的关系。旋转周期边界需要一个矢量变量的修正。这个将在后面进行讨论。单元顶点的方案（双重控制体积）这种情况在图 8.10 中进行了描述。一种处理周期性边界的方法是由围在阴影控制体积的表面组成的合成通量来解决。剩余的在点 1 和点 2 合起来组成完整的净通量，因此：(8.40) 1,2,1 RRsumsum和点 1 和点 2 的部分控制体积也应该加上。在移动性的周期变化边界下，反相边界剩余的控制体没有变化。例如： 于是就有 ，旋转性的周期21 和 sumsu,2

33、,1 变化边界在利用公式（8.40）之前需要对动量方程进行变化。旋转性周期旋转性周期的情况是基于旋转的坐标系统。因此，所有的矢量项，例如速度或者标量的梯度都会随着一同转换。标量项例如压力或者密度对于坐标的旋转是不变的。如果我们假设旋转坐标和 x 轴平行，旋转矩阵变为：(8.41)cosin0i01R介于周期性边界 A 和 B 的角 顺时针是正的。因此，比如，速度矢量从边界 A 移动到 B 可以写成（图 8.9 中的单元 和图 8.10 中的点 ）2,141(8.42)ABvR很容易看出 x 轴分量 通过旋转没有变化。因此 。所有的流体量的梯),.(Aueiv ABu度以相同的方式转移。正如上面

34、陈述的，在以单元为顶点的方案中，剩余的动量方程在使用公式（8.40）必须要进行修正。式（8.41）中的旋转矩阵的应用可以变成如下的公式：(8.43)wvuAvuBwvusmRR,图（8.43）的上标 表示的是三个动量方程。wvu,8.8 介于平板网格的界面在讨论 3.1 节中构建的网格进行空间离散时，很明显通常它可能在复杂的几何区域内产生一个单独的网格。 （如图 3.4） 。我们提出两种可能的方案来解决这个问题。第一种方法是多平板方法，第二种是 技术。在下文，我们将主要阐述多平板方法。关于多平Chimera板进一步的讨论，读者可以参考文献【39】-【45】 。文献 【46】对多平板方法很好的进

35、行了介绍。更多的关于 技术的可以在文献【47】- 【51】 ，这里我们不做讨论。在多平板方法中，物理区域被分割成确定数量的虚拟部分。因此，计算区域同样分成了相同的平板数。总体上，求解一个特定的平板的物理量需要流场在相邻的一个或多个平板。因此，我们需要提供一个可以有效的在平板之间进行信息交换的数据结构。如果不同的处理器来解决平板的控制方程的话，这个结构也要求其可以交流。第一部分的数据结构包括平板边界的编号。图 8.12 表示了一种特殊的编号方法。这种编号的思想总结起来如下所示：boundary 1: =IBEGiboundary 2: =IENDboundary 3: j =JBEGbounda

36、ry 4: j =JENDboundary 5: k =KBEGboundary 6: k =KEND所有的平板用相同的编号方案是非常重要的。计算区域中的网格点指数 i，j ，k 定义如下：IBEG i IENDJBEG j JENDKBEG k KEND单元的指数 IJK 需要以相同的方法定义为以单元为中心的方案。由于多平板方法经常通过虚拟的单元来解决，因此物理单元在开始或者结束时将会有一定的偏移（见图 8.1）每一个平板的边界被分成一些不重叠的小块。在相同的边界条件下，这种分类可以规范不同的边界条件。图 8.13 对此进行了描述。唯一区别每一个小块的方法是相应的平板数量以及平板边界的数量。

37、此外，小块的原点，高度和宽度同样需要存储。正是由于这个原因，图 8.13 用坐标 L1BEG , L1END ,L2BEG 和 L2END 来表示。小块坐标轴按照环形方向来确定的。与就是说，如果我们设定了 i 轴，j 和 k 就将分别是第一和第二个环形方向。在j 轴上，环形的方向将会分别变成 k 和 i。因此，由于图 8.13 中的小块是在 j =JBEG 的边界，在 k 方向设定了 轴，在 i 方向设定了 轴。这种环形方向唯一定义了小块的位置。1l2l这些表示介于平板间数据结构的剩余部分确保了其可以在这些小块中互相交换。由于这个原因，我们需要扩展上述通过相邻的平板和小块的数据结构。流场中介于

38、两个平板间的一些量的变化在图 8.14 中进行了描述。这个过程主要包括两步。第一步，区域里通过相邻小块虚拟层重叠的变量可以写成单独的虚拟单元或者可以暂时的存储（图 8.14 中的 ）所有的平板都要这样做。第二步，平板间的数据和 BA进行交换。这就意味着 可以写成 B 板的虚拟层， 可以写成 A 板的虚拟层。和 BA 如果两个小块有不同的位置，数据必须相应的进行转换。如果网格线和平板的表面不匹配，可以参考文献【52】 【53】 。8.9 不规则网格边界的流体梯度5.3.4 我们已经讲过在双重中线的方案中需要注意一些流体梯度的值。如果在三角形或者四边形网格中计算梯度，通过式子（5.46）我们可以用

39、 Green-Gauss 方法来计算。由式（8.12）可以计算网格中的边界（不包括对称边界和周期性边界）或者用（8.13）来代替平均算法。否则，梯度的计算精度将会不准确。考虑到图 8.4 中的项，边界节点 1 的贡献可以写成： 25611SnU）（这里 是边界表面节点 1 和 2 之间的长度（二等分） 。图 8.5 中节点 1 三角形表面12S1-2-3 的贡献由对应的公式（8.13）可以表示为： 3681221Sn）（是三角形 1-2-3 中灰色的区域。在混合的网格中，通常比较适合用小方格来表示312S虚拟的边界。 【8】 （见图 5.15）以单元为中心的方案对对称和周期性边界没有特殊的要求。通量求解的方法与 8.5 和 8.7节相同。这个同样适用于双重中线方案，如果梯度用最小方格法，只需要在 8.5 节中将某些梯度设定为零即可。 （公式 8.37）如果在双重中线方案中用 Green-Gauss 方法，则需要对贴于边界的控制体积表面的法向矢量进行修正。 （像图 8.4 中的 点。 ）如果垂直于对称平面，则可以设表面矢量的分量*2为零即可。最后，梯度可以由 8.5 节中的方法进行修正。在周期性的边界中，边界两面的梯度和体积需要以当前的通量相加（如等式（8.40） ） 。在旋转性周期的情况下，梯度需要通过公式（8.41）提供的旋转矩阵来进行转换。