模拟落在水面上的球体

把比水轻的球丢入装有水的小烧杯中后，当球在水面上漂浮时，如何模拟水面的形状和球体的运动？本文我们演示了如何在 COMSOL Multiphysics® 软件中对该系统进行建模。尽管我们考虑的是小球体的情况，但此处所用讨论的技术同样适用于其他较大的结构。

浮球的振荡运动：一个 计算流体力学（CFD）问题

首先，将一个浮球放入一个盛满水的烧杯中。球体的半径为 0.15cm，密度为 800kg/m³，最开始悬浮在空气中。 圆柱形烧杯的半径为 0.65cm，高度为 1.7cm，装有部分水。由于该系统的尺寸足够小，因此假设流动模式为层流是合理的，这也使我们能够在普通计算机上实现快速计算。我们可以用完全相同的方式对较大的系统建模，但需要更多的计算资源。

为了模拟球体穿过水面的运动过程，我们需要对空气相以及水相中的流体流动进行建模。如果将球体视为刚体，则可以通过求解常微分方程(ODE)来计算其运动。常微分方程描述了球体从空气中掉落，撞击水面，并最终与空气和水相互作用过程中的受力情况。流体在固体边界上施加的总流体应力为：

(1)

式中，u_fluid表示流体速度场，p表示流体压力，μ 表示流体的动态黏度，n表示边界的向外法线，I表示单位矩阵。

在一个充满水的烧杯中，浮球悬浮在水面上方的空气中。隐藏了烧杯壁的一部分以显示球体。

假设球体仅沿垂直方向（z）运动，我们可以建立如下方程：

(2)

式中，m表示质量，v_s表示在z轴方向上的速度，u_s表示球体在z方向上的位移。

在球体表面上，对方程1中的流体应力的z方向分量进行积分可得到力F_z。将上述运动方程（常微分方程）与有限元模型一起求解，可以计算球体的位移和速度。

由于球体会在空气以及水中运动，为了同时模拟两种流体，我们在随时间变化的空间坐标系中建立了两相流模型。空间坐标系是一种常见的、固定的全局欧式系统，其坐标系为(x，y)。这些与流体域的网格节点相对应的空间坐标系是时间的函数，表示空间某一点当前的位置。这可以通过在动网格接口中求解流体流动问题来实现。

球与流体接触面处的网格位移等于球体的位移。在流体域内部，基于动网格接口中合适的平滑算法来计算网格变形：

(3)

在球体与流体边界上，流体速度由下式计算：

(4)

在 COMSOL Multiphysics® 中建立模型

对于上述问题，我们可以建立二维轴对称模型。为了对此问题进行建模，我们将方程1–4作为表达式，将两相流接口与动网格接口以及全局方程式进行耦合。

二维轴对称模型的几何形状，用水平线表示水面的初始位置。在常微分方程中，我们将球体视为刚体，因此可以定义球体为模型域的边界。

在此模型中，空气相以及水相中的流体流动均为层流，因此我们可以使用层流两相流，相场接口（阅读此博客，了解有关如何选择多相流接口的更多信息。）该多物理场接口自动定义了层流和相场之间的耦合，可以用来追踪水面形状变化的细节。

通过层流两相流，相场接口，我们可以将球体的边界条件看作移动的润湿壁进行模拟。通过相场接口的润湿壁边界条件，我们可以设定水面与球体壁之间的接触角；而在层流接口中，通过使用壁边界条件中的壁移动设置，我们可以设定球体壁的移动速度，在该烧杯的顶部使用出口边界条件。

下图显示了模型设置的结构和各种接口（即层流，相场，和动网格）：

在相场接口中，对于初始设置，我们将相场变量的值设定为 Φ=1（流体1为空气）和 Φ=-1（流体 2 为水），水平线代表初始时刻的水平面，即初始时刻的两相界面。水平线上方为空气，下方为水。初始设置如上图所示。

我们利用润湿壁边界条件来设定两相界面在烧杯壁以及球体壁上的接触角。为简单起见，我们假定烧杯壁以及球体壁上的接触角相等。

流体中的速度场会自动耦合至相场中来解释对流。在层流与相场接口耦合的多物理场接口中，设定流体的属性和表面张力系数。
动网格接口计算了流体域中由于球体运动而导致的网格变形。我们还使用了边界条件来描述流体边界处的网格应如何变形。这里重要的边界条件是网格位移，它等于球体边界上的位移(方程3)。

球体的运动完全归因于重力载荷。然而，流体施加了方程1中所定义的曳力。该曳力z方向上的分量可以通过在整个球体表面上对内置变量进行积分得到。内置变量 spf.T_stressz 用于计算周围流体施加在球体表面上的总应力的z方向分量。

了解更过信息，请阅读博客如何计算升力和曳力。

使用全局方程式接口求解 常微分方程（方程2）。下图显示了在全局方程式接口中使用的方程式之一，其中vs是球的速度，如方程2所述。vst是速度（或加速度）的时间导数，而 intop1 是在实体边界上定义的积分算子。

设置全局方程式接口求解球的速度。

由于球体的移动较大，因此网格变形将很大。为了保持较高的单元质量并避免单元反转，我们可以使用“自动重新划分网格” 功能。该功能可以基于网格质量暂停瞬态计算，并在继续计算之前重新划分当前网格。

仿真结果分析

下图显示了落在水面上的球体，水-空气界面（Φ= 0.5）以及流体速度的横截面图。

充满水的烧杯中的浮球的运动，水面和流体速度大小的横截面动画。为了在仿真的后期看到速度幅度的变化，图例中颜色的最大限值为 0.1[m/s]。

下图为球体位移与时间的关系图。与预期一样，球体的运动受周围的流体阻碍作用，位移幅度随时间呈指数下降。

结语

在这篇博客中，我们演示了如何对浮球的振荡运动进行建模。文中，球体及其运动被看作是轴对称的，没有旋转运动。通常，如果需要，可以将此处讨论的建模策略应用于旋转和平移模型。对这种情况进行建模，需要使用常微分方程求解刚性物体的其他旋转运动。相反，使用固体力学接口的内置功能来自动解决刚体运动会容易得多。该接口也允许模拟一个更大的系统，在该系统中，气流会变为湍流，并且可以在所有方向上产生运动。

单击下面的按钮，访问此博客中介绍的模型的文件：