5.6.4 浸入边界法
Peskin在1972年提出浸入边界(Immersed Boundary，IB)法，并将IB应用于血液在可收缩的心脏瓣膜中流动的数值模拟。IB的产生为解决复杂边界情况下的流固耦合提供了新的思路和方法。在流体运动方程中，IB添加了体积力项来描述物体边界和流场之间的相互作用，并在数值计算中使用两套网格体系。一方面，使用固定不变的欧拉网格计算流场；另一方面，使用一系列拉格朗日点来表示物体的边界。然后通过两套网格的信息交换来演化固体边界和流场之间的相互作用。当前，IB在模拟血液流动、湍流的直接模拟以及多相流等方面取得了巨大的成功，在流体力学领域有着很好的发展前景。
图5.22 浸入边界及节点之间的关系[11]
X—浸入边界点的位置；x—流场节点的位置；h—欧拉点的节点间距
接下来，详细介绍IB的基本原理及程序实现。浸入边界与流场节点之间的关系如图5.22所示。图中的空心点称为欧拉点(在LBM中的网格节点)，该点本身不发生移动，只描述某个固定的位置，在这些位置上的流体速度和压力由Navier-Stokes方程描述；图中的实心点称为拉格朗日点，该点可在流体作用下发生移动，并且相邻点之间存在内力(拉伸、压缩、弯曲、剪切等)约束，其运动规律由流体和各相邻点之间的内力共同控制，由这些点组成的边界即浸入边界。
浸入边界法的核心公式为
式中，δ(·)——Dirac delta函 数，简 称Delta函数。
Delta函数的形式不止一种，本质上是以一种以空间距离为依据的加权函数。一种常用的格式为
式中，r——浸入边界点与欧拉点在坐标轴方向的距离。
浸入边界法有两个核心步骤。
步骤1，将浸入边界上的内力(体积力)分配到周围，即实施公式f(x，t)=∫F(s，t)δ(x-X)d s的描述，如图5.23所示。
图5.23 拉格朗日点上的力分配至附近的欧拉点示意(书后附彩插)
对图5.23左图中间的灰黑圆点A，假如其位置坐标是(5.62，7.55)，该点上x方向的内力分量是0.002、y方向的内力分量是0.001。要将这两个力分量体现在对流场的影响上，就需要分配和转移其到流场节点x(欧拉点)上。图5.23右图是转移结果示意。其转移规则用Delta函数控制。式(5.22)中的r是节点X与周围欧拉点x在x方向或y方向的距离，
的最大值是2；若超过2，则δ(r)=0。因此，点A实际上只与周围的16个欧拉点存在关系，即图5.23中红色虚线所圈定的范围。
浸入边界上的节点距离也是受到限制的，一般情况下，相邻点距离d s在0.4～0.8个网格距离。如果同时考虑浸入边界上所有点上内力对流场节点的影响，那么分配到流体上的力则分布在浸入边界两侧，如图5.24所示。图中，红色的“*”代表浸入边界上的点，蓝色的点则是浸入边界节点上的力所影响的周围流体节点范围。
下面列出IB步骤1的核心代码，其功能是将浸入边界节点的内力转移分配到流场节点上。
图5.24 浸入边界(拉格朗日点)与附近流体节点(欧拉点)之间的关系(书后附彩插)
其中，第2、3行中的floor()函数表示向下取整，如floor(2.3)=2、floor(5.89)=5。因此，对于本例点A的坐标为(5.62，7.55)的情况，“ii=floor(IBx(1，i))-1：floor(IBx(1，i))+2”返回的结果为4、5、6、7，而“jj=floor(IBy(1，i))-1：floor(IBy(1，i))+2”返回的结果为6、7、8、9。这分别是点A周围16个欧拉点的x坐标和y坐标。第6～9行，计算Delta函数，也就是式(5.22)的代码形式，这里有
δ(xi)=1。第10、11行，将分配得到的拉格朗日力累加到欧拉点上，是公式f(x，t)=∫F(s，t)δ(x-X)d s的程序表述。
步骤2，将点A附近16个欧拉点上的速度，采用δ(·)函数加权后，累加到点A上，从而使得点A也具有一个运动速度——式(5.21)中的U(s，t)。这一过程如图5.25所示。
图5.25 欧拉点的速度加权合成拉格朗日点的速度示意(书后附彩插)
从图5.25中可以看到，速度的传递方向正好和步骤1中浸入边界点内力的传递方向相反。此外，如果将这个运动速度与时间步长相乘，就可以更新点A的位置，对所有浸入边界节点进行位置更新，就可以实现运动边界的模拟，即公式
=U(s，t)的表述。
对照步骤1的代码，步骤2的MATLAB代码如下：
容易看出，该代码与步骤1代码的差异在于第10、11行，这里的10、11行是公式U(s，t)=∫u(x，t)δ(x-X)d x的代码表述。