模拟落在水面上的球体

2018年 2月 4日

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

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

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

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

(1)

\mathbf{f} = \mathbf{n} \cdot \left\{ -p\mathbf{I}+\left(\mu\left(\mathbf{\nabla} \mathbf{u}_{\text{fluid}}+\left(\mathbf{\nabla} \mathbf{u}_\text{fluid}\right)^{T}\right)-\frac{2}{3}\mu\left(\mathbf{\nabla}\cdot\mathbf{u}_{\text{fluid}}\right)\mathbf{I}\right)\right\}

式中, ufluid表示流体速度场, p 表示流体压力, μ 表示流体的动态黏度, n 表示边界的向外法线, I 表示单位矩阵。
A model of a sphere falling on a water surface.
在一个充满水的烧杯中,浮球悬浮在水面上方的空气中。烧杯壁的一部分已隐藏以显示出球体。

假设球体仅沿垂直方向(z)运动,我们可以建立如下方程:

(2)

\begin{align}
F_z &= m g + m\frac{d v_{s}}{d t}\\
v_s &= \frac{d u_s}{d t}
\end{align}

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

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

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

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

(3)

\left.\mathbf{u_{\text{mesh}}}\right|_{\text{Interfacing boundaries}} = \left(0,u_{\text{s}}\right)

在球体与流体边界上,流体速度由下式计算:

(4)

\left.{\mathbf{u}_{\text{fluid}}}\right|_{\text{Interfacing boundaries}} = \left(0,v_s\right)

在COMSOL Multiphysics® 中建立模型

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

The geometry of a 2D axisymmetric model of a sphere above the water surface.
二维轴对称模型的几何形状,用水平线表示水面的初始位置。在ODE中我们将球体视为刚体,因此可以定义球体为模型域的边界。

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

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

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

A screenshot of the Model Builder in COMSOL Multiphysics with three interfaces highlighted.

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

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

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

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

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

A screenshot of the Global Equations settings in COMSOL Multiphysics.

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

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

仿真结果分析

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

 

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

下图为球体位移与时间的关系图。正如预期,它的运动受周围的流体阻碍作用,位移幅度随时间呈指数下降。
A 1D plot of the sphere displacement versus time.

结语

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


评论 (0)

正在加载...
浏览 COMSOL 博客