用于滑动网格模拟的体积 LBM 边界方案的扩展

heu_shd

已于 2025-02-24 15:18:38 修改

阅读量989

点赞数 23

文章标签：算法 LBM 计算流体力学体积边界滑动网格

于 2025-02-24 11:52:33 首次发布

本文链接：https://blog.youkuaiyun.com/weixin_75184155/article/details/145823598

版权

在计算流体力学（CFD）领域，Lattice Boltzmann Method（LBM）作为一种新兴的数值模拟方法，近年来受到了广泛关注。与传统的 Navier-Stokes 方程求解方法不同，LBM 基于粒子动力学，能够更自然地处理复杂几何和多物理问题。

1.1 滑移网格算法的动机

在工业应用中，许多流体力学问题涉及旋转机械，例如无人机旋翼、风扇、泵、涡轮等。这些设备的流场具有复杂的三维非定常特性，传统的固定网格方法难以准确模拟。为了解决这一问题，滑移网格技术应运而生。滑移网格技术通过将计算域划分为多个子域，每个子域使用独立的网格系统，从而能够灵活地处理旋转和静止部件之间的相对运动。

1.2 滑移网格的基本概念

1.2.1 域划分

滑移网格算法的核心思想是将计算域划分为两个或多个子域，每个子域可以独立地进行网格划分和求解。例如，一个包含旋转部件和静止部件的系统可以划分为两个子域：一个子域随旋转部件旋转，另一个子域固定在静止部件上。这种划分方式使得每个子域的网格系统可以独立优化，从而提高计算效率和精度。

1.2.2 局部旋转参考系和网格

在滑移网格算法中，旋转部件所在的子域通常使用局部旋转参考系。通过引入旋转参考系，旋转部件在该参考系中表现为静止，从而简化了流场的求解。具体来说，局部旋转参考系中的流场方程需要考虑由于旋转引起的惯性力，例如离心力、科里奥利力、欧拉力。

1.2.3 接口耦合

滑移网格算法的关键在于如何处理不同子域之间的接口耦合。由于不同子域的网格系统可能不完全对齐，因此需要一种有效的算法来传递流体信息，同时确保质量和动量的守恒。LBM 滑移网格算法通过引入体积边界算法（volumetric boundary algorithm）来实现这一目标。

1.3 LBM 滑移网格算法

1.3.1 在非惯性参考系中求解 LBE

在非惯性参考系中，LBM 的基本方程需要考虑外部力场的影响。具体来说，D3Q19LBGK的LBE（Lattice Boltzmann Equation）可以写为：

$f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{1}{\tau} \left[ f_i(\mathbf{x}, t) - f_i^{eq}(\mathbf{x}, t) \right] + F_i(x, t)$

其中：

$f_i(x, t)$ 是粒子分布函数。
$c_i$ 是离散速度。
$\tau$ 是松弛时间。
$f_i^{eq}(x, t)$ 是平衡分布函数。
$F_i(x, t)$ 是群力。

1.3.2 力群体的计算

群力 $F_i(x, t)$ 与外部体力场 $A(x, t)$ 之间的关系可以表示为：

$F_i = \left(1 - \frac{1}{2\tau}\right) \rho w_i \left( \frac{(c_i - u'(x, t)) \cdot A(x, t)}{T} + \frac{(c_i - u') \cdot (c_i - A(x, t))}{T^2} \right)$

其中：

$\rho$ 是流体密度。
$w_i$ 是离散速度的权重。
$u'(x, t)$ 是非惯性参考系中的速度。
$A(x, t)$ 是外部体力场。
$T$ 是温度。

1.3.3 速度转换

在非惯性参考系中，流体的速度不仅受到惯性参考系的影响，还需要考虑参考系的加速度或旋转效应。因此，速度的转换公式需要进行修正，以考虑参考系的相对运动。根据半力矫正的思想，转换公式如下：

$u'(x, t) = u(x, t) + \frac{A(x, t)}{2}$

其中：

$u(x, t)$ 是惯性参考系中的速度。
$u'(x, t)$ 是非惯性参考系中的速度。
$A(x, t)$ 是外部体力场。

1.3.4 非惯性参考系中的惯性力

非惯性参考系中的惯性力可以表示为：

$A(x, t) = -\Omega(x, t) \times (\Omega(x, t) \times r(x, t)) - 2\Omega(x, t) \times u'(x, t) - \frac{d\Omega(x, t)}{dt} \times r(x, t)$

其中：

$\Omega(x, t)$ 是旋转角速度。
$r(x, t)$ 是位置向量。
$u'(x, t)$ 是非惯性参考系中的速度。

其中公式的第一项表示离心力，第二项是科里奥利力，第三项对应于欧拉力。对于旋转恒定的情况，欧拉力将变为零。

1.3.5 速度转换公式

在惯性参考系和非惯性参考系之间，速度的转换公式为：

$u_B(x, t) = u_A(x, t) - \Omega(x, t) \times r(x, t)$

其中：

$u_A(x, t)$ 是惯性参考系中的速度。
$u_B(x, t)$ 是非惯性参考系中的速度。
$\Omega(x, t)$ 是旋转角速度。
$r(x, t)$ 是位置向量。

1.3.6 界面动力学

在滑移网格界面处，流体信息的传递需要确保质量和动量的守恒。具体来说，需要在两个参考系之间传递密度和速度信息，并确保这些信息在界面处的连续性和守恒性。

为了构建这样的界面动力学，我们可以直接使用镜面反射算法，因为它具有良好的特性，可以保持局部质量和切向动量守恒。可以证明，通过简单地将构建的表面平衡中的切向速度替换为全速度矢量，所得算法能够守恒滑动网格界面上的质量和总动量通量。这种修改后的镜面反射算法构成了当前滑动网格接口算法的基线。

1.3.6.1 构造界面处的平衡状态

$f_{\alpha, eq}^{i, A}(t) = \bar{\rho}_\alpha w_i \left[ 1 + \frac{c_i \cdot \bar{u}_\alpha^A}{T} - \frac{(c_i \cdot \bar{u}_\alpha^A)^2}{2T^2} + \frac{(c_i \cdot \bar{u}_\alpha^A)^3}{6T^3} - \frac{(c_i \cdot \bar{u}_\alpha^A)^2}{2T^2} (\bar{u}_\alpha^A)^2 \right]$

如前所述，此处使用全速度矢量对于滑动网格界面上的质量和总动量通量守恒至关重要。其中基于两侧采样值的平均密度为：

$\bar{\rho}_\alpha = \frac{1}{2} \left( \rho_{\alpha, A} + \rho_{\alpha, B} \right)$

基于两侧采样值的质量平均速度：

$\bar{u}_\alpha A = \frac{1}{2 \bar{\rho}_\alpha} \left( \rho_{\alpha, A} u_{\alpha, A} A + \rho_{\alpha, B} u_{\alpha, B} A \right)$

（质量平均速度（mass-averaged velocity） 是指在某一特定区域内，不同组分的流体速度按照各自的质量权重加权计算得到的平均速度，在处理流体和物体表面之间的相互作用时，质量平均速度用于描述不同流体或组分之间的相对速度，通过这种方式，质量平均速度能够综合考虑多种流体组分的影响，提供一个更准确的速度描述。）

其中：

$f_{\alpha, eq}^{i, A}(t)$ 是在参考系 A 中，界面处的平衡分布函数。
$\bar{\rho}_\alpha$ 是界面处的平均密度。
$w_i$ 是离散速度的权重。
$c_i$ 是离散速度。
$u_{\alpha}^{A}$ 是在参考系 A 中，界面处的平均速度。
$T$ 是温度。

1.3.6.2 界面处的匹配算法

在 T = 0 时，表面网格的内层和外层彼此重合，当 T > 0 时，内网格旋转一个 θ 的角度，这导致表面网格的内层和外层不匹配：

根据几何信息：

$\omega(t) = \frac{A_{\beta 1 B \rightarrow \alpha A}}{A_\alpha}$

$A_{\alpha} = A_{\beta_1 B \rightarrow \alpha A} + A_{\beta_2 B \rightarrow \alpha A}$

$1 - \omega(t) = \frac{A_{\beta_2 B \rightarrow \alpha A}}{A_{\alpha}}$

其中：

$\omega(t)$ 是匹配比例。
$A_{\beta 1 B \rightarrow \alpha A}$ 是在参考系 B 中的表面元素 $\beta_1$ 投影到参考系 A 中的面积。
$A_\alpha$ 是在参考系 A 中的表面元素 $\alpha$ 的面积。

1.3.6.3 界面处的流体信息交换

有了上述定义，曲面 $\alpha$ （位于域 A）从域 B 侧对密度和速度的采样可以直接链接到其匹配的对曲面 $\beta_1$ 和 $\beta_2$ 的密度和速度采样：

$\rho_{\alpha, B} = \omega(t) \rho_{\beta_1, B} + (1 - \omega(t)) \rho_{\beta_2, B}$

$u_{\alpha, B}^A = \omega(t) \rho_{\beta_1, B} u_{\beta_1, B}^A + (1 - \omega(t)) \rho_{\beta_2, B} u_{\beta_2, B}^A$

其中：

$\rho_{\alpha, B}$ 是在参考系 B 中，表面元素 $\alpha$ 的密度。
$\rho_{\beta_1, B}$ 和 $\rho_{\beta_2, B}$ 是在参考系 B 中，表面元素 $\beta_1$ 和 $\beta_2$ 的密度。
$u_{\alpha, B}^A$ 是在参考系 A 中，表面元素 $\alpha$ 的速度。
$u_{\beta_1, B}^A$ 和 $u_{\beta_2, B}^A$ 是在参考系 A 中，表面元素 $\beta_1$ 和 $\beta_2$ 的速度。

1.4 操作流程

标准的 surfel（表面元素）采集、碰撞和散射过程可以直接应用于滑网界面两侧的 surfel 元素。整体的解决过程可以按照以下步骤进行：

收集入射状态：从滑网界面两侧的相邻单元格中采样质量和动量，并收集入射状态。
匹配滑网界面元素：在滑网的外部和内部表面之间进行滑网界面元素的匹配。如前所述，在任何给定时间，某一侧的每个表面元素与另一侧的两个表面元素进行通信。两个表面元素的匹配比例依赖于滑网的旋转位置。
转换表面速度：将表面速度从局部旋转参考系转换为滑网内的静止实验室参考系。
计算质量和动量的平均值：在固定的参考框架下，计算每对匹配滑网界面元素的质量和动量的平均值。然后，也计算滑网内部的动量值。
构造平衡状态：在两侧领域的界面 surfel 元素上构造平衡状态。
应用 surfel 碰撞：在所有的 surfel 元素上执行碰撞过程。
确保局部质量守恒：在整个过程中确保质量守恒。
散射出射粒子：将出射粒子散射回其相邻的单元格。