IB-LBM浸入边界法

浸入边界法1IB-LBM的基本数学构造Hofler和Schwarzer在2000年提出了一种在对NS方程求解的有限差分方程中增加一项体积力的IBM来模拟充满粒子流场的固体粒子的运动。其粘性不可压流的控制方程表示为：+uuupu+ft（1）0u（2）f：作用在流体上的外力上述两个公式在IB-LBM中可有格子Boltzmann方程所代替，即：1,,,,τeqxtttfxtxtfxtFtffe（3）012pfftue（4）24112ssFccfeueue（5）F表示作用在物体表面上的力，其物理过程包括：碰撞过程：1,,,,τeqxttfxtxtfxtFftf（6）迁移过程：,,xtttfxtfe（7）由此可知，碰撞仅仅是局部的，但迁移和邻近点皆有关系。平衡态分布函数为：242()22122eqsssufccceueu（8）其中sc为声速，s/3cc式中表示各离散速度方向的加权系数，其值分别为：40911,2,3,4915,6,7,836（9）该模型中，流体的松弛时间和运动粘性系数的关系可表示为：212sct（10）1.1IB-LBM速度修正法在LBM中，同时满足NS方程的宏观量是通过Chapman-Enskog展开得到的。在速度修正法里其表达式方程有如下表示：*uef（11）12uft（12）其中，*u为中间层速度，u为修正速度。结合方程（11），（12），方程（4）可变为：*uuu（13）在传统的IB-LBM中，f由胡克定理、直接力法或者动量交换法来求得。其中胡克定理表达式为：,fluidwallFstkkVtVt（14）x,,,ftFstxXstds（15）其中,xXst称为狄拉克插值函数，其作用是将作用在物面上的拉格朗日坐标系下的外力插值到欧拉坐标系下对流体点的作用力，描述了流体和固体之间相互作用的关系。如果单从上述方程来求外力f的话，流体会穿越物面，导致质量不守恒。为了解决流体穿透的问题，该方法在求解f时做了相应处理，即假设f为非预知力，由满足无滑移边界条件的速度修正场求得。方程（13）中，设u为欧拉坐标系下的修正速度，有：*,,,ijijijuxtuxtuxt（16）其中，*u可由方程（11）求得。依此，由方程（14）、（15）类似推导出u的求法。假设lBU表示拉格朗日坐标系下的修正速度，相应地，则有：,,,BBBuxtUXtxXstds（17）其中狄拉克插值函数可以光滑地近似为：,,,,lllBijijBBBxXstDxXstxXstyYst(18)Peskin提出,BxXst表达式为：11cos,2420,2rrrr＞（19）因此可应用方程（18）将方程（17）离散：,,,,1,2,lllijBBijijBlluxtUXtDxXstslm…,（20）由此类推可得：,,,,lllBBijijijBijUXtuxtDxXstxy（21）将方程（20）代入方程（16）中，则有*,,,,lllijijBBijijBlluxtuxtUXtDxXsts（22）将方程（22）代入方程（21）可得*,,,,,llllBBijijBBijlijijijlUXtuxtDxyUXtDsDxy（23）将方程（23）改写成：*,,,,,llllBBijijBBijlijijijlUXtuxtDxyUXtDsDxy（24）其中,lijijijBDDxXst，12,,,,,lmBBBBXUUUU……，12,,,,,lmBuuuu……，*,,,1,2,,lllBBijijijuUXtuxtDxylm…于是可将方程（24）简化为：AXB宏观量压强和密度分别为：f，2spc（25）作用在流体上的外力可由方程（12），求得2/fut（26）计算步骤如下：（1）给定初始状态值，计算出矩阵A的值,并求得矩阵A的逆A-1;(2)给定F的初始值为0，即0F，应用式(3)求ntt时的密度分布函数值，通过式(11)和式(25)，计算该时刻的流场宏观量，如宏观量、p等;（3）求解方程组AXB，获得物面边界点的修正速度lBU，再通过式(20)将求得的物面边界点的修正速度lBU插值到欧拉点上，也即求得欧拉点的修正速度iju；（4）利用式（16）修正流场中欧拉点的速度，同时应用式(26)求得欧拉点上的力密度函数值f；（5）利用式（8）求平衡态分布函数；（6）重复(2)-(5)，直到收敛。