在IBM中，将结构对流体的作用等效为体积力。传统的LBM方程不包含这个体积力。如果需要考虑这个体积力的话，其相应的控制方程也需要发送变化。

体积力对方程的影响主要集中在粒子的速度上。

1 控制方程

在LBE中考虑作用力的影响，直观的方法是在演化方程中增加一个作用力项：

$$
f_{i}\left(x+c_{i} \delta t, t+\delta t\right)-f_{i}(x, t)=-\frac{1}{\tau}\left[f_{i}(x, t)-f_{i}^{e q}(x, t)\right]+\delta t F_{i}
$$

作用力项中的不同形式，导致不同的作用力模型。这里考虑的是郭照力老师提出来的CZS模型。

在LBGK模型中，发送一次碰撞的时间为：

$$
\tau \delta
$$

由作用力引起的流体动量变化为：

$$
\mathbf{F} \tau \delta
$$

因此，可以看到，体积力的影响主要是影响粒子的宏观速度。通过改变平衡态分布函数，来考虑作用力的影响。而平衡态分布函数，是密度与宏观速度的函数。即：

$$
f_{i}^{e q}=E_{i}\left(\rho, \mathbf{u}^{*}\right)
$$

其中带星号的u表示的是平衡态速度。其值与普通的LBGK方程中的平衡态速度不一样，普通的平衡态速度：

$$
\rho \mathbf{u}=\sum_{i} \mathbf{e}{i} f{i}
$$

修正后的平衡态速度：

$$
\rho \mathbf{u}^{*}=\sum_{i} \mathbf{e}{i} f{i}+\frac{\delta t}{2} \mathbf{F}
$$

修正后的平衡态速度与标准的平衡态速度相比，多了关于体积力的一项。这实际上是对碰撞前后的速度的一个平均。

GZS模型的流体宏观速度和修正的平衡态速度一样：

$$
\rho \mathbf{u}=\sum_{i} \mathbf{e}{i} f{i}+\frac{\delta t}{2} \mathbf{F}
$$

力源项的计算：

$$
F_{i}=\left(1-\frac{1}{2 \tau}\right) \omega_{i}\left[\frac{\mathbf{e}{i}-\mathbf{v}}{c{s}^{2}}+\frac{\left(\mathbf{e}{i} \cdot \mathbf{v}\right)}{c{s}^{4}} \mathbf{e}_{i}\right] \cdot \mathbf{F}
$$

总结来看，其控制方程主要是4个：

  1. 包含外力的LBE方程:

$$
f_{i}\left(\mathbf{x}+\mathbf{e}{i} \Delta t, t+\Delta t\right)-f{i}(\mathbf{x}, t)=-\frac{1}{\tau}\left(f_{i}(\mathbf{x}, t)-f_{i}^{e q}(\mathbf{x}, t)\right)+F_{i} \Delta t
\tag{1.1}
$$

  1. 修正的平衡态速度方程:

$$
\rho \mathbf{u}^{*}=\sum_{i} \mathbf{e}{i} f{i}+\frac{\delta t}{2} \mathbf{F}\tag{1.2}
$$

  1. 修正的流体速度方程：

$$
\rho \mathbf{u}=\sum_{i} \mathbf{e}{i} f{i}+\frac{\delta t}{2} \mathbf{F}\tag{1.3}
$$

  1. 体积力计算方程：

$$
F_{i}=\left(1-\frac{1}{2 \tau}\right) \omega_{i}\left[\frac{\mathbf{e}{i}-\mathbf{v}}{c{s}^{2}}+\frac{\left(\mathbf{e}{i} \cdot \mathbf{v}\right)}{c{s}^{4}} \mathbf{e}_{i}\right] \cdot \mathbf{F}\tag{1.4}
$$

  1. 平衡态分布函数：

$$
f_{i}^{e q}=E_{i}\left(\rho, \mathbf{u}^{*}\right)\tag{1.5}
$$

2 离散方程

离散主要是针对包含外力的LBE方程，即对以下方程进行离散：

$$
f_{i}\left(\mathbf{x}+\mathbf{e}{i} \Delta t, t+\Delta t\right)-f{i}(\mathbf{x}, t)=-\frac{1}{\tau}\left(f_{i}(\mathbf{x}, t)-f_{i}^{e q}(\mathbf{x}, t)\right)+F_{i} \Delta t \tag{2.1}
$$

  • 碰撞：碰撞和此前的LBGK变化不大，只是增加了体积力而已：

$$
f_{i}^{\star}(\mathbf{x}, t)=f_{i}(\mathbf{x}, t)-\frac{\Delta t}{\tau}\left[f_{i}(\mathbf{x}, t)-f_{i}^{e q}(\mathbf{x}, t)\right]+F_{i} \Delta t \tag{2.2}
$$

带星号的分布函数表示的是碰撞之后的分布函数。此时并没有进行时间推移

  • 迁移：

$$
f_{i}\left(x+c_{i} \Delta t, t+\Delta t\right)=f_{i}^{\star}(x, t) \tag{2.3}
$$

3 关于与IBM的结合

包含外力项的流体方程，其最大的作用在于实现IBM，即浸入边界法。因为浸入边界的核心就是把边界对流体的作用等效为一个力源项。

因此，我们来看，在和IBM结合的时候，如何来使用。实际上，不同的边界力计算模型，会采用不同的算法。这里我们首先只考虑一个简单的情况，即边界力采用罚方法，或者解析力法从结构域直接求得。

首先我们假设边界力是已经求出来了。即F我们已知。

  1. 这样，我们就要分别通过根据式1.2，1.3计算出平衡态速度u*以及宏观流体速度u；
  2. 根据计算得到的平衡态速度u*，根据式1.5计算得到平衡态分布函数feq；
  3. 根据宏观流体速度u，通过式1.4计算得到离散的体积力Fi；
  4. 把离散的体积力Fi和平衡态分布函数feq带入到式2.2，执行碰撞操作；
  5. 执行迁移操作，至此完成一次循环计算。

4 Palabos里的外力项模型

4.1 包含外力项模型的动态

在Palabos里面，包含外力项模型的动态为：ExternalForceDynamics

其继承关系图如下图：

ExternalForceDynamics的继承关系
ExternalForceDynamics的继承关系

我们来看看，ExternalForceDynamics下的计算宏观速度的两个方法：computeVelocity和computeVelocityExternal。先看第一个： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
template<typename T, template<typename U> class Descriptor>
void ‪ExternalForceDynamics::computeVelocity (
        ‪Cell const& cell, ‪Array::d>& u ) const
{
    T rhoBar;
    ‪Array::d> force, j;
    ‪momentTemplates::get_rhoBar_j(cell,rhoBar, j);
    force.‪from_cArray(cell.‪getExternal(Descriptor::ExternalField::forceBeginsAt));
    
    T invRho = Descriptor::invRho(rhoBar);
    for (‪plint iD = 0; iD < Descriptor::d; ++iD)
    {
        u[iD] = j[iD]*invRho + force[iD]/(T)2;
    }
}

这个是ExternalForceDynamics下面的计算宏观速度的其中一个方法，输入是单元的引用，及一个速度矢量。通过该方法，可以计算单个单元的宏观速度。可以看到，其宏观速度的计算方法如式1.3.

其中：momentTemplates::get_rhoBar_j(cell,rhoBar, j);从cell里面取得rhoBar和j，用于后续的计算。force.‪from_cArray(cell.‪getExternal(Descriptor::ExternalField::forceBeginsAt));从cell里面取得外力，用于后续的计算。

第二个宏观速度的计算方法computeVelocityExternal ：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
template<typename T, template<typename U> class Descriptor>
void ‪ExternalForceDynamics::computeVelocityExternal (
        ‪Cell const& cell, T rhoBar, ‪Array::d> const& j,
        ‪Array::d>& u ) const
{
    ‪Array::d> force;
    force.‪from_cArray(cell.‪getExternal(Descriptor::ExternalField::forceBeginsAt));
    
    T invRho = Descriptor::invRho(rhoBar);
    for (‪plint iD = 0; iD < Descriptor::d; ++iD)
    {
        u[iD] = j[iD]*invRho + force[iD]/(T)2;
    }
}

computeVelocityExternal 这个方法和上面的几乎一样，不一样的在于computeVelocityExternal中的rhobar，j，都是从外部给与，而不是从cell中获得。

4.2 体积力的计算

接下来，我们来看看体积力的计算，即式1.4的计算。体积力的计算主要源自于addGuoForce。addGuoForce这个方法有好几种存在形式，通用的形式存在于模板externalForceTemplatesImpl中。其代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
static void ‪addGuoForce( ‪Array& f, T* externalScalars,
                          ‪Array const& u,
                          T omega, T amplitude )
 {
     static const int forceBeginsAt = Descriptor::ExternalField::forceBeginsAt;
     T* force = externalScalars + forceBeginsAt;
     for (‪plint iPop=0; iPop < Descriptor::q; ++iPop) {
         T c_u = T();
         for (int iD=0; iD
             c_u += Descriptor::c[iPop][iD]*u[iD];
         }
         c_u *= Descriptor::invCs2 * Descriptor::invCs2;
         T forceTerm = T();
         for (int iD=0; iD < Descriptor::d; ++iD) {
             forceTerm +=
                 (   ((T)Descriptor::c[iPop][iD]-u[iD]) * Descriptor::invCs2
                      + c_u * (T)Descriptor::c[iPop][iD]
                 )
                 * force[iD];
         }
         forceTerm *= Descriptor::t[iPop];
         forceTerm *= 1-omega/(T)2;
         forceTerm *= amplitude;
         f[iPop] += forceTerm;
     }
 }

这个计算来自于externalForceTemplatesImpl模板。可以看到，这里面的计算和公式1.4是一致的。输出的，或者说计算得到的是输入的f。输入包含一个外部双精度的指针变量：externalScalars，以及速度矢量u，松弛因子omega和幅值amplitude。

这个forceBeginsAt默认是0。属于描述子的一个静态的public属性。目前不知道其用途。

T* force = externalScalars + forceBeginsAt;表明了作用力的计算来自于外部标量externalScalars和这个forceBeginsAt

至于后面的循环，其实就是为了计算Fi。

当基于D2Q9的时候，externalForceTemplatesImpl模板延伸出了一个externalForceTemplates2D的模板，其中重新定义了一个addGuoForce:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
static void ‪addGuoForce (
                 ‪Array::q>& f,
                 T* externalScalars,
                 ‪Array::d> const& u, T omega, T amplitude )
 {
     static const int forceBeginsAt
         = ‪descriptors::ForcedD2Q9Descriptor::ExternalField::forceBeginsAt;
     T* force = externalScalars + forceBeginsAt;
     T mu = amplitude*((T)1-omega/(T)2);
     
     static const T oneOver3 = (T)1/(T)3;
     static const T oneOver12 = (T)1/(T)12;
     static const T fourOver3 = (T)4/(T)3;
     
     const T twoUx   = (T)2*u[0];
     const T threeUx = (T)3*u[0];
     
     const T twoUy   = (T)2*u[1];
     const T threeUy = (T)3*u[1];
 
     f[0] += -fourOver3*mu*(force[0]*u[0]+force[1]*u[1]);
     
     f[1] += oneOver12*mu*(force[0]*(-(T)1+twoUx-threeUy)+force[1]*((T)1+twoUy-threeUx));
     
     f[2] += oneOver3*mu*(force[0]*(-(T)1+twoUx)-force[1]*u[1]);
     
     f[3] += oneOver12*mu*(force[0]*(-(T)1+twoUx+threeUy)+force[1]*(-(T)1+twoUy+threeUx));
     
     f[4] += -oneOver3*mu*(force[0]*u[0]+force[1]*((T)1-twoUy));
     
     f[5] += oneOver12*mu*(force[0]*((T)1+twoUx-threeUy)+force[1]*(-(T)1+twoUy-threeUx));
     
     f[6] += oneOver3*mu*(force[0]*((T)1+twoUx)-force[1]*u[1]);
     
     f[7] += oneOver12*mu*(force[0]*((T)1+twoUx+threeUy)+force[1]*((T)1+twoUy+threeUx));
     
     f[8] += -oneOver3*mu*(force[0]*u[0]+force[1]*(-(T)1-twoUy));
 }

很明显，这里直接对公式1.4在D2Q9的格子描述子下面做了进一步的实现。这样可以提高计算的效率。本质上，和普通版的addGuoForce是一样的。

4.3 包含体积力的碰撞实现

在搞定了宏观速度的计算和体积力的计算之后，就需要考虑包含体积力的碰撞实现了。根据式2.2可以看到，有两点值得注意的地方：

  • 碰撞项包含了外力项：Fi*dt
  • 碰撞项的平衡态分布函数feq的计算的宏观速度要采用基于公式1.4得到的数值

我们知道，碰撞是一种局部的行为，一般在dynamic中实现。对于Guo格式的外力模型，Palabos里面已经写好了。从包含外力项的动态的继承图里可以看到，这个动态就是：GuoExternalForceBGKdynamics。实际上，还有GuoExternalForceCompleteRegularizedBGKdynamics。这个叫正则化的Guo外力动态。

我们来看看GuoExternalForceBGKdynamics里面的碰撞的实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
template<typename T, template<typename U> class Descriptor>
void ‪GuoExternalForceBGKdynamics::collide (
        ‪Cell& cell,
        ‪BlockStatistics& statistics )
{
    T rhoBar = this->‪computeRhoBar(cell);
    ‪Array::d> u, j;
    this->‪computeVelocity(cell, u);
    T rho = Descriptor::fullRho(rhoBar);
    for (‪plint iD = 0; iD < Descriptor::d; ++iD)
    {
        j[iD] = rho * u[iD];
    }
    
    T uSqr = ‪dynamicsTemplates::bgk_ma2_collision(cell, rhoBar, j, this->getOmega());
    ‪externalForceTemplates::addGuoForce(cell, u, this->getOmega(), (T)1);
    
    if (cell.‪takesStatistics()) {
        ‪gatherStatistics(statistics, rhoBar, uSqr);
    }
}

输入参数：cell，即单元，statistics是一个BlockStatistics类型的变量。

第一句T rhoBar = this->‪computeRhoBar(cell);计算宏观速度u。然后接下来是计算j，即rho*u。把这个j带入碰撞计算式dynamicsTemplates::bgk_ma2_collision(cell, rhoBar, j, this->getOmega());可以对该单元执行碰撞处理。

值得注意的是，这里面在执行完标准的碰撞之后，还执行了一项：addGuoForce。这个命令相当于把体积力加到的分布函数上。

看式2.2，前面的标准碰撞相当于把 $ f_{i}^{\star}(\mathbf{x}, t)=f_{i}(\mathbf{x}, t)-\frac{\Delta t}{\tau}\left[f_{i}(\mathbf{x}, t)-f_{i}^{e q}(\mathbf{x}, t)\right] $实现了，后面的addGuoForce相当于把计算得到的体积力的影响加到了分布函数上，实际上是Fi*dt这部分。

来看看addGuoForce的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
static void ‪addGuoForce( ‪Array& f, T* externalScalars,
                          ‪Array const& u,
                          T omega, T amplitude )
 {
     static const int forceBeginsAt = Descriptor::ExternalField::forceBeginsAt;
     T* force = externalScalars + forceBeginsAt;
     for (‪plint iPop=0; iPop < Descriptor::q; ++iPop) {
         T c_u = T();
         for (int iD=0; iD
             c_u += Descriptor::c[iPop][iD]*u[iD];
         }
         c_u *= Descriptor::invCs2 * Descriptor::invCs2;
         T forceTerm = T();
         for (int iD=0; iD < Descriptor::d; ++iD) {
             forceTerm +=
                 (   ((T)Descriptor::c[iPop][iD]-u[iD]) * Descriptor::invCs2
                      + c_u * (T)Descriptor::c[iPop][iD]
                 )
                 * force[iD];
         }
         forceTerm *= Descriptor::t[iPop];
         forceTerm *= 1-omega/(T)2;
         forceTerm *= amplitude;
         f[iPop] += forceTerm;
     }
 }

这里的核心代码在于f[iPop] += forceTerm;不难看出，分布函数包含了此前的计算结果，是在这个基础上再加上了forceTerm。而forceTerm又

因此可以看到，在GuoExternalForceBGKdynamics中的collide是实现了体积力部分的计算的。

4.4 如何添加体积力

现在，只剩下一个问题，那就是如何添加体积力。体积力是一个矢量，对应的是cell里面的External。

看看这个External是什么呢：

1
typedef ‪ExternalFieldArraytypename Descriptor::ExternalField> ‪External;

可以看到其本质是一个外部场向量。类型就是描述子里面的ExternalField。如果描述子是Force2dDescriptor，那么这个ExternalField就是Force2dDescriptor。

在Force2dDescriptor中这样定义ExternalField的。

1
2
3
struct ‪Force2dDescriptorBase {
     typedef ‪Force2dDescriptor ‪ExternalField;
 };

那么Force2dDescriptor里面有些什么呢？

1
2
3
4
5
6
struct ‪Force2dDescriptor {
    static const int ‪numScalars = 2;
    static const int ‪numSpecies = 1;
    static const int ‪forceBeginsAt = 0;
    static const int ‪sizeOfForce   = 2;
};

numScalars是标量的个数，对于力来说，当然有2个分量。‪numSpecies种类只有1个，forceBeginsAt为0，代表不发生偏置。sizeOfForce为2，是代表Force的大小为2。

在Cell里面有一个方法，setExternalField，用这个方法来对外部场进行设置。

1
2
3
4
void ‪setExternalField(‪plint pos, ‪plint size, const T* ext) {
         ‪PLB_PRECONDITION( ‪dynamics );
         ‪dynamics->setExternalField(*this, pos, size, ext);
     }

cell里面的setExternalField调用的是dynamic里面的setExternalField。

Dynamic里面的setExternalField如下：

1
2
3
4
5
6
7
8
9
10
template<typename T, template<typename U> class Descriptor>
 void ‪Dynamics::setExternalField (
         ‪Cell& cell, ‪plint pos, ‪plint size, const T* ext)
 {
     ‪PLB_PRECONDITION(pos+size <= descriptor::ExternalField::numScalars);
     T* externalData = cell.‪getExternal(pos);
     for (‪plint iExt=0; iExt
         externalData[iExt] = ext[iExt];
     }
 }

T* externalData = cell.‪getExternal(pos);获得cell的外部场的指针。pos代表的是偏置，一般就是0，即无偏置。对externalData的赋值是采用的址传递，即cell那个指向外部场的指针被使用了，结果就是改变了指向的外部场变量。

因此，在设置单个cell的外部场时，我们可以调用cell的这个setExternalField方法。

1
2

cell.setExternalField(‪plint pos, ‪plint size, const T* ext);

其中，ext是一个二维数组的指针。因为力都是有两个分量的。

至此，对于包含外力的LBM方法从原理上和Palabos的代码实现上应该有一个全局的观念了。利用这个包含外力的LBM模型，我们可以做的有3个方面：把外力赋予到单元的外部场中，根据这个外力计算宏观速度，执行碰撞。至于根据外力计算外力项的离散项，则其实在collide的内部就已经完成了。