模拟落在水面上的球体

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

浮球的振荡运动：一个 CFD 问题

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

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

(1)

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

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

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

(2)

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

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

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

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

(3)

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

(4)

在COMSOL Multiphysics® 中建立模型

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

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

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

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

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

在“ 相场” 接口中，对于初始设置，我们将相场变量的值设定为Φ= 1（流体1为空气）和 Φ= -1（流体 2 为水），水平线代表初始时刻的水平面，即初始时刻两相界面。水平线上方为空气，水平线下方为水。初始设置如上图所示。
我们利用润湿壁 边界条件来设定两相界面再烧杯壁以及球体壁上的接触角。为简单起见，我们假定烧杯壁以及球体壁上的接触角相等。

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

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

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

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

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

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

仿真结果分析

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

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

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

结语

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