自由表面马蹄涡作为界面流动中一类典型的涡旋结构，广泛存在于自然界（如水黾水面滑行^[1]）与工程场景（如自由表面射流^[2]、仿生推进装置）中. 该涡旋结构随时间演化，气相中的涡丝会附着在界面上，其水下部分形态与马蹄铁相似，因此被形象地称为马蹄涡^[3]. 它们通常由固体边界脱落的涡层卷绕形成，在表面张力、黏性效应与涡动力学耦合作用下，呈现出涡量输运、拓扑重构、二次涡生成等复杂演化特性. 自由表面马蹄涡不仅直接影响动量传递效率与流动稳定性，还关联着界面波动、射流生成、能量耗散等关键物理现象，在船舶推进、深海资源开发、流体混合等工程领域具有重要应用价值.

早期研究中，Denny^[4]提出“Denny悖论”，指出传统界面波推进理论无法解释水黾之类生物的高效运动，Suter等^[5]的实验研究也支持这一结论. Hu等^[6]通过高速摄影与粒子追踪技术证实，自由表面马蹄涡主导的动量传递是解决该悖论的关键，为后续涡旋与界面相互作用研究奠定了重要基础. 此后，相关研究从实验观测、数值模拟到理论分析逐步深入，形成了多维度的探索体系.

在数值模拟方面，已有大量研究围绕自由表面条件下涡结构的形成、演化及其与界面动力学的耦合机制展开. Sarpkaya^[7]从基础流体力学角度探讨了表面张力与涡量之间的相互作用机理，指出界面曲率变化可显著影响局部涡量分布与耗散特性. 张霄月等^[8]采用数值方法研究了自由表面涡流动的形成与发展过程，揭示了界面变形对涡结构稳定性和流动非定常性的调制作用. Walker等^[2]对自由表面射流中的湍流结构进行了系统分析，阐明了界面附近湍动能分布与各向异性特征. Herlina等^[9]则模拟了强湍流驱动下的气液界面传质过程，强调了自由表面附近涡结构对传质增强的主导作用. Xiong等^[10-11]通过数值分析系统研究了涡管的螺旋特性及其演化规律，指出螺旋度在刻画涡与界面耦合过程中具有重要意义. 陈红勋等^[12] 基于格子Boltzmann方法模拟了科氏力作用下吸水口前自由表面旋涡的生成与演化过程. 这些数值工作为深入理解自由表面涡及其参数依赖性奠定了重要基础.

在实验研究方面，自由表面与涡结构相互作用的非线性与拓扑变化特征得到了大量直观验证. Gharib等^[13]通过流动可视化实验定性描述涡管在自由表面处发生断开与重联的现象，首次揭示了界面在涡拓扑演化中的主动作用. Ji 等^[14]通过实验研究涡环空化过程，系统分析了环状空泡的振荡周期、壁面冲击特性以及环量计算模型，为涡环空化动力学提供了实验依据. 在工程应用背景下，Gu 等^[15]通过实验研究核电站余热排出系统中自由表面涡诱导的气体卷吸行为，揭示了界面变形与旋涡强度之间的耦合关系. Mondal 等^[16]在多出口排液实验中发现，自由表面涡的相互干扰会显著改变吸气特征与界面拓扑形态. 针对马蹄涡结构本身的形成与演化特征，Polat等^[17]利用 PIV 技术实验研究往复旋转圆柱周围马蹄涡的生成过程. Zhang等^[18]通过泵吸水池实验定量分析自由表面涡的演化特征及其对水力稳定性的影响. 这些实验工作为数值模拟结果提供了重要验证，并为揭示涡与自由表面相互作用的本质机理奠定了基础.

尽管上述研究涵盖了丰富的流场特征与工程现象，但随着对微观机理认识的深入，表面张力（即韦伯数）被证实在界面动力学中扮演着至关重要的角色. Suerich-Gulick等^[19]在自由表面吸气涡的研究中发现，表面张力显著决定了涡核的生成尺度；Bandyopadhyay^[20]探讨界面附近的涡结构不稳定性，指出表面张力在其中起到了关键的抑制与调控作用. 特别是Zhang等^[21]针对气泡涡环的研究进一步表明，韦伯数直接决定了涡环的流态转捩路径，证实了其在多相流拓扑演化中的核心控制地位. 在此基础上，对经典实验的再审视也支撑了这一观点：Zhang等^[22]的实验虽侧重于揭示涡重联机制，但其观察到的界面涡量重组过程，本质上反映了表面张力对涡量输运的约束机制. 在更复杂的耦合流动现象方面，Terrington等^[23]的研究揭示了不同韦伯数主导下，界面拓扑连接特性的时序差异；Zhang等^[24]关于空化涡环驱动自由表面射流的研究则表明，当惯性力突破表面张力束缚时，即高韦伯数条件下，涡结构能够反向诱导自由表面响应.

这些研究共同指向一个结论：在涉及剧烈界面变形的流场中，韦伯数是调控界面破碎阈值、拓扑重构及能量耗散过程的核心物理参数. 然而，当前关于不同韦伯数条件下，自由表面马蹄涡这一特定拓扑结构的动能演化规律、耗散机制特征、螺旋度守恒特性及表面能转换效率等关键物理量的量化研究仍有待进一步拓展.

为了研究不同$ We $数下自由表面马蹄涡的演化特性及调控机制，本研究以界面流动涡动力学理论为基础，采用分段线性界面重构（piecewise linear interface construction, PLIC）的体积分数（volume of fluid, VOF）^[25]数值方法，系统研究$ We=0.24\sim 0.84 $范围内自由表面马蹄涡的演化差异. 通过可视化分析捕捉相界面波动特征、涡核迁移轨迹及涡量场时空演化规律，结合总动能衰减趋势、耗散率分布、半域总螺旋度守恒偏差及表面能累积量的定量统计，建立$ We $数与各关键物理量的量化关联，阐明$ We $数对涡环与界面作用机制的调控规律.

两相气泡环系统的动力学行为由不可压缩 Navier-Stokes 方程描述，辅以界面捕捉方法和表面张力模型.

流体运动遵循质量守恒和动量守恒定律. 对于不可压缩流动，连续性方程和动量方程分别表示为

(1)$ \left.\begin{split} & \nabla \cdot \boldsymbol{u}=0,\\&\rho \left(\dfrac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\cdot \nabla \boldsymbol{u}\right)=-\nabla p+\nabla \cdot \left[\mu \left(\nabla \boldsymbol{u}+\nabla {\boldsymbol{u}}^{{\mathrm{T}}}\right)\right]+{\boldsymbol{f}}_{\sigma }+\rho \boldsymbol{g}.\end{split}\right\} $

式中：$ \boldsymbol{u} $为速度矢量，$ p $为压力，$ \rho $和$ \mu $分别为密度和动力黏度，$ \boldsymbol{g} $为重力加速度，$ {\boldsymbol{f}}_{\sigma } $为表面张力项.

气液两相通过体积分数函数$ \varphi (x,t) $区分，其中$ \varphi =1.0 $表示纯液相，$ \varphi =0 $表示纯气相，$ 0 \lt \varphi \lt 1.0 $表示界面过渡区域. 物性参数采用线性混合模型：

(2)$ \left.\begin{split} \rho (\varphi )=\varphi {\rho }_{\text{g}}+(1-\varphi ){\rho }_{\text{l}},\\\mu (\varphi )=\varphi {\mu }_{\text{g}}+(1-\varphi ){\mu }_{\text{l}}.\end{split}\right\} $

式中：下标$ \text{g} $和$ \text{l} $分别代表气相和液相.

界面演化通过体积分数输运方程描述：

(3)$ \frac{\partial \varphi }{\partial t}+\nabla \cdot (\varphi \boldsymbol{u})=0. $

该方程描述体积分数随流场对流，维持界面运动的几何一致性.

表面张力效应采用连续表面力模型（continuum surface force , CSF）^[26]处理：

(4)$ {\boldsymbol{f}}_{\textit{σ}}=\sigma \kappa {\delta }_{\text{s}}\boldsymbol{n}. $

式中：$ \sigma $为表面张力系数，$ \kappa $为界面曲率，$ \boldsymbol{n} $为界面单位法向量，$ {\delta }_{\text{s}} $为界面狄拉克函数. 界面法向量和曲率由体积分数梯度计算：

(5)$ \boldsymbol{n}\mathbf=\frac{\mathbf{\nabla }\alpha }{|\mathbf{\nabla }\alpha |},\;\;\kappa \mathbf=-\mathbf{\nabla \cdot }\boldsymbol{n}. $

狄拉克函数$ {\delta }_{\text{s}} $通过$ |\mathbf{\nabla }\alpha | $近似，将表面张力体积力集中于界面附近区域.

控制方程的数值求解基于有限差分法，在交错网格上离散，并采用投影法^[27]处理压力与速度的耦合问题. 时间推进使用显式格式，空间离散兼顾精度与稳定性.

计算域离散为均匀三维网格，变量采用交错存储：速度分量存储于网格面心，压力、密度、黏度和体积分数存储于网格体中心. 对流项离散采用三阶迎风格式，以有效抑制数值振荡；扩散项采用二阶中心差分，以保障格式对称性.

时间积分采用二阶Adams-Bashforth格式^[28]进行时间推进. 首先，计算中间预测速度$ {u}^{*} $，其中忽略压力梯度项（标准投影法）. 动量方程离散形式如下：

(6)$ \begin{split} \rho^n \frac{\boldsymbol{u}^*-\boldsymbol{u}^n}{\Delta t}= & \frac{3}{2}\left[-\nabla \cdot\left(\rho^n \boldsymbol{u}^n \boldsymbol{u}^n\right)+\nabla \cdot\left(\mu ^ { n } \left(\nabla \boldsymbol{u}^n+\right.\right.\right. \\& \left.\left.\left.\left(\nabla \boldsymbol{u}^n\right)^{\mathrm{T}}\right)\right)\right]-\frac{1}{2}\left[-\nabla \cdot\left(\rho^{n-1} \boldsymbol{u}^{n-1} \boldsymbol{u}^{n-1}\right)+\right. \\& \left.\nabla \cdot\left(\boldsymbol{u}^{n-1}\left(\nabla \boldsymbol{u}^{n-1}+\left(\nabla \boldsymbol{u}^{n-1}\right)^{\mathrm{T}}\right)\right)\right]+\boldsymbol{f}_\sigma^n .\end{split}$

式中：上标n表示时间层，$ {\boldsymbol{u}}^{*} $为中间速度. 随后，压力-速度耦合通过投影法实现. 首先求解压力泊松方程以获得$ n+1 $时刻的压力场：

(7)$ \mathbf{\nabla \cdot }\left(\frac{1}{{\rho }^{n}}\mathbf{\nabla }{p}^{n\mathbf+1}\right)\mathbf=\frac{1}{\Delta t}\mathbf{\nabla \cdot }{{\boldsymbol{u}}}^{*}. $

接着利用压力梯度校正速度场，得到满足无散条件的速度$ {\boldsymbol{u}}^{n+1} $：

(8)$ {\boldsymbol{u}}^{n\mathbf+1}\mathbf={\boldsymbol{u}}^{*}-\frac{\Delta t}{{\rho }^{n}}\mathbf{\nabla }{p}^{n\mathbf+1}. $

泊松方程采用多重网格法求解，确保计算效率.

此外，体积分数$ \varphi $的输运通过施主−受主方法进行时间推进，其离散形式可写为

(9)$ \left.\begin{split} \varphi^{n+1} =&\varphi^n-\frac{\Delta t}{\Delta V} \sum_f F_f^n, \\F_f^n =&\int_f^n \varphi^u \boldsymbol{u}^{n+1} \cdot \boldsymbol{n} {\mathrm{d}} S .\end{split}\right\} $

式中：$F_f^n $表示通过单元面$ f $的体积分数通量，$ \Delta \textit{V} $为单元体积，$ \Delta t$为时间步长，S为单元面的面积，n为法向量.

界面重构和通量计算采用PLIC方法^[25,29-31]. 该方法的优势在于良好的质量守恒、良好的界面几何恢复能力以及在工程尺度算例上的鲁棒性. 同时，下个时间步的混合密度和混合动力黏度依据局部体积分数进行线性插值更新.

最后，数值模拟的时间步长受库朗-弗里德里希斯-刘维（Courant-Friedrichs-Lewy CFL）条件和毛细波稳定性限制：

(10)$ \Delta t\leq \min \left(\dfrac{{C}_{\text{CFL}}\Delta x}{|{\boldsymbol{u}}{|}_{\max }},\;\;{C}_{\sigma }\sqrt{\frac{\rho \Delta {x}^{3}}{2{\text{π}} \sigma }}\right). $

式中：$ {C}_{\text{CFL}} $和$ {C}_{\sigma } $为稳定性常数，通常取 0.3~0.5，本研究采用0.4；${|{\boldsymbol{u}}{|}_{\max }} $为最大速度模长.

为了系统研究韦伯数对气泡环界面重联的影响，设计一系列对照算例，覆盖表面张力主导至惯性主导的不同流态. 算例参数基于无量纲分析确定，确保物理机制的普适性.

计算域是边长为$ L $的立方盒子，网格采用$ {L}_{x}\times {L}_{y}\times {L}_{z}=512\times 512\times 512 $（无量纲单位）的均匀离散. 算例基础设置参数见表1.

如图1所示，对于初始时刻的相态分布，以$ y=0.5L $为分界面，上半部分（$ y=0.5L\sim 1.0L $）为气相（$ \varphi =0 $），下半部分（$ y=0 \sim 0.5L $）为液相（$\varphi = 1 .0$）. 垂直于分界面的涡环是由强加旋度的速度场构成，相应地也被分界面分成上半部分和下半部分. 该环半径$ {R}=0.1L $，涡核半径$ a=0.01L $，以获得涡环结构的几何合理性. 涡核半径与网格尺度的比值为$ a/\Delta x=5.12 $，该比值满足网格尺度小于涡核尺度，有利于能量集中同时避免因网格过密导致的计算量冗余. 环半径与网格尺度的比值为$ R/\Delta x=51.2 $，涡环周向覆盖网格数$ 2{\text{π}} R/\Delta x\approx 321 $，能够较好分辨涡环与气液界面相互作用的形态演变. 因此，该网格尺度能保障分辨流场初始条件的可靠性. 此外，环量$ \textit{Γ} $根据雷诺数确定，固定黏性效应以专注于韦伯数的影响. 通过调节表面张力系数来实现不同韦伯数的工况设置，以探究惯性力与表面张力竞争关系对界面动力学行为的影响.

气泡涡环的初始速度场采用Rosenhead-Moore涡丝模型生成. 由于本研究采用欧拉网格计算，须将连续的积分形式转化为离散求和形式. 在拉格朗日框架下将涡环中心线离散为$ M=2\;000 $个离散涡元，基于毕奥-萨伐尔定律的离散化形式，对欧拉交错网格上的每个面心速度分量进行积分求和：

(11)$ \boldsymbol{u}(\boldsymbol{x})=\sum\limits_{{\mathrm{l}}=1}^{M}-\frac{\varGamma }{4{\text{π}} }\frac{(\boldsymbol{x}-{\boldsymbol{x}}_{{\mathrm{l}}})\times ({\boldsymbol{t}}_{{\mathrm{l}}}\Delta s)}{{\left(|\boldsymbol{x}-{\boldsymbol{x}}_{{\mathrm{l}}}{|}^{2}+{\delta }^{2}\right)}^{3/2}}. $

式中：$ \boldsymbol{x} $为欧拉网格面心位置，$ {\boldsymbol{x}}_{\rm{l}} $为涡元位置，$ {\boldsymbol{t}}_{\rm{l}}\Delta s $ 近似于涡元段的切向向量，$ \delta =a{{\mathrm{e}}}^{-0.75} $为去奇异化的涡核平滑尺度. 这一离散化的数值计算方法有效避免了涡丝中心处的数值奇点. 具体而言，先遍历交错网格中的每一个面心，先确定该面心的空间位置，并初始化诱导速度为零向量；接着逐个遍历$ M $个涡元段，计算面心到当前涡元段的位置向量及其模长，通过式(11)更新诱导速度；最后将计算得到的诱导速度投影到该面的法向轴上，以此更新面心处的速度.

体积分数的初始化基于超采样方法. 首先遍历每一个网格单元$ (i,j,k) $，初始化该单元内位于气相气泡内部的子采样点计数$ {N}_{{\mathrm{inside}}} $=0；随后遍历该单元内$ N_{{\mathrm{sub}}}^{3} $个超采样子网格中心点，计算每个子采样点$ {\boldsymbol{x}}_{{\mathrm{sub}}} $相对于气泡涡环（圆环面）的符号距离函数值：

(12)$ \phi (\boldsymbol{x})=\left[{{\left[\left[{{(x-{{x}_{\rm{c}}})}^{2}+{(y-{{y}_{\rm{c}}})}^{2}}\right]^{1/2}-R\right]}^{2}+{(z-{{z}_{\rm{c}}})}^{2}}\right]^{1/2}-a. $

式中：$\boldsymbol{x}=[x, y, z] \in {{\bf{R}}}^3 $为三维空间中任意点的位置矢量，对应计算网格单元中心的坐标； $\boldsymbol{x}_{\mathrm{c}}=[x_{\mathrm{c}}, y_{\mathrm{c}}, z_{\mathrm{c}}]\in {{\bf{R}}}^3 $为环心的位置矢量.

若$ \phi \leq 0 $则说明该点位于气相气泡内部，此时将$ {N}_{{\mathrm{inside}}} $加1. 待所有子采样点遍历完成后，通过如下公式计算该网格单元的液相体积分数：

(13)$ {\varphi }_{i,j,k}=1-{{N}_{{\mathrm{inside}}}}/{N_{{\mathrm{sub}}}^{3}}. $

最后，基于式(2)，根据体积分数$ {\varphi }_{i,j,k} $更新该单元的密度$ {\rho }_{i,j,k} $和黏度$ {\mu }_{i,j,k} $.

最后，为了确保初始场物理一致性，在初始诱导速度场构建后额外进行无散投影处理，强制满足无散条件$ \nabla \cdot \boldsymbol{u}=0 $. 由于初始条件不是稳定状态，为此，本研究延长了数值模拟的时间，将预计算时间延长至无量纲时间$ {t}^{*}=0.1 $. 该时长基于初始动能、耗散率演化曲线确定. 此时，物理量的演化达到大致稳态.

网格分辨率是影响数值模拟结果精度与计算效率的关键因素，为了验证本研究所采用数值方法的计算准确性、算法精度及网格收敛性，确保模拟结果的可靠性与严谨性，通过多分辨率网格对比及经典基准算例对比的方式进行验证分析.

定义广义特征尺度如下：选取环半径$ R $为特征长度$ {L}_{0} $，则特征速度和特征时间分别为

(14)$ {U}_{0}=\varGamma /\left(2{\text{π}} {L}_{0}\right),\;\;{T}_{0}={L}_{0}/{U}_{0}. $

基于上述特征尺度，定义2个主导流体动力学特性的核心无量纲数.

1）表征惯性力与黏性力平衡关系的雷诺数：

(15)$ Re={{\rho }_{0}{U}_{0}{L}_{0}}/{{\mu }_{\rm{l}}}. $

2）量化惯性力与表面张力效应竞争关系的韦伯数：

(16)$ We={{\rho }_{0}U_{0}^{2}{L}_{0}}/{\sigma }. $

为了确定兼顾计算精度与效率的最优网格分辨率，开展多尺度网格收敛性验证. 选取所模拟的典型算例为验证对象，设置3种不同的三维网格分辨率，分别为N³=256³，512³，768³，保持其余数值计算参数、初始条件与边界条件完全一致，以归一化总动能的时间衰减曲线为评价指标，对比不同网格分辨率下的计算结果差异. 总动能表示为

(17)$ {E}_{{\mathrm{k}}}=\frac{1}{2}{\int}_{V}\rho (\boldsymbol{x},t)|\boldsymbol{u}(\boldsymbol{x},t){|}^{2}\,\text{d}V. $

式中：$ \rho $表示流体密度，$ \boldsymbol{u}(\boldsymbol{x},t) $为速度矢量，$ V $为整个计算域.

网格无关性验证结果如图2所示. 其中，纵坐标中的$ {E}_{{\mathrm{k}},0} $表示初始时刻的总动能，采用$ {E}_{{\mathrm{k}}}/{E}_{{\mathrm{k}},0} $对总动能进行归一化处理，有利于消除初始条件差异对结果对比的影响，从而更直观地评估不同网格分辨率下能量衰减的收敛特性；横坐标采用无量纲时间$ {t}^{*}=t/{T}_{0} $来衡量时间变化. 可以看出，512³与768³网格分辨率下的数值差异均小于 1%，表明当网格分辨率提升至512³后，继续加密网格对计算结果的影响可忽略不计. 该分辨率可在保证计算效率的前提下，满足数值模拟的精度要求，由此确定本研究后续所有模拟计算均采用512³的网格分辨率.

为了验证本研究数值方法对表面张力效应的捕捉能力，选取经典基准算例低黏度液体球滴聚并震荡^[32]为参考，将本研究模拟结果与该基准算例的经典结果进行定量与定性对比.

验证算例的关键物理参数设置与基准算例保持一致. 在雷诺数$ Re=1.9\times {10}^{4} $条件下，模拟2个等径自由液滴在表面张力驱动下的聚并震荡的过程. 验证结果如图3所示. 其中，左半边黑色实线为Sprittles的结果，右半边青色界面为本研究的结果. 从液滴聚并的形态演化特征来看，本研究模拟结果复现了基准算例中，液桥形成、自由表面曲率演化、液滴融合及后期振荡的动力学特征；从定量指标来看，自由表面关键位置的曲率变化及聚并特征时间等参数与基准算例的相对误差均小于3%. 上述结果显示，本研究数值方法可合理描述表面张力驱动下的液滴流动与聚并动力学过程，对表面张力效应的模拟效果良好，进一步佐证了方法的可靠性.

通过网格无关性验证与表面张力效应验证，确定了本研究数值模拟的较为合适的网格分辨率为512³，同时验证了所采用的界面捕捉方法、离散格式和表面张力模型能够较好地捕捉表面张力主导的液滴聚并震荡的动力学特征. 数值方法的计算精度、算法稳定性及结果可靠性均满足研究要求，为本研究后续的数值模拟研究提供了可靠的基础.

如图4 所示，展示了2种典型韦伯数（$ We=0.28 $和$ We=0.42 $）下相界面在不同无量纲时间$ {t}^{*} $的演化形态. 其中，粉红色界面对应 ${We=}{0.28} $，蓝色界面对应${We=}{0.42} $. 视角为与$ \textit{y=}0.5L $平面平行的水平中线投影，其中粉色界面对应低$ We $数（高表面张力），蓝色界面对应高$ We $数（低表面张力）. 在早期阶段（$ {\textit{t}}^{*}=0.3 $），2种工况的界面均在 $ \textit{y=}0.5\textit{h} $ 平面下方形成2个对称的角状突起. 低$ We $情况下突起较为圆钝，界面保持紧致；高 $ We $情况下尖角更为细长，显示出更强的惯性驱动拉伸特征. 至$ {\textit{t}}^{*}=0.6 $，高$ We $界面两侧尖角发生接触并产生拓扑重联，形成闭合的环状结构，而低$ We $界面在表面张力的约束下未发生角尖接触，仍保持开口构型. 在中期演化过程中（$ {\textit{t}}^{*}=0.6\sim 2.1 $），高$ We $工况下形成的环状界面能够稳定维持，并在$ {\textit{t}}^{\textit{*}}\approx 2.1 $之后开始显现局部变薄及分裂趋势. 相比之下，低$ We $界面的角状结构自$ {\textit{t}}^{*}=1.2 $起逐渐退化并被表面张力平滑，界面趋向于扁平化演化，未出现闭合环形成. 在晚期阶段（$ {\textit{t}}^{*}\geq 2.4 $），高$ We $界面表现出更复杂的拓扑变化，包括局部断裂及子结构生成，而低$ We $界面虽存在起伏与褶皱，但整体连通性保持良好. 上述结果表明，韦伯数通过调控惯性力与表面张力之间的竞争关系，决定了相界面是否经历尖角接触、拓扑重联及后续裂解过程：低 $ We $时表面张力主导界面稳定性，高$ We $时惯性效应主导界面拉伸与拓扑重构.

如图5所示，展示了特定韦伯数$ We=0.42 $下，相界面与涡量场的同步演化过程，包含相界面与涡量等值面的三维特征及切平面内的涡量分布与涡线结构. 为理解多相流中界面与涡量耦合机制提供了直观的物理图像.

图5（a）中白色透明结构为相界面，彩色等值面对应涡量场的空间分布. 在早期阶段（$ {\textit{t}}^{*}=0.5 $），相界面呈现非对称的角状突起，涡量等值面紧密附着于界面附近的强剪切区域，体现界面变形与涡量生成的耦合关系；随时间推进至$ {\textit{t}}^{*}=1.0 $，相界面的角尖因惯性拉伸进一步细长化，涡量等值面的分布范围随界面扩展同步增大；在中期（$ {\textit{t}}^{*}=1.5 $），相界面发生拓扑重联后形成的环状结构，对应涡量等值面的局部聚集；至$ {\textit{t}}^{*}=2.0 $，相界面趋向扁平化，涡量等值面的聚集程度显著降低，反映界面动力学对涡量输运的调控作用.

图5（b）展示了$ \textit{y=}0.5L $切平面上的涡量模$ |\boldsymbol{\omega }| $分布及关键涡线（黑色线圈）的演化过程. $ {\textit{t}}^{*}=0.5 $时，涡量高值区集中于相界面在该平面的投影区域，涡线呈紧凑螺旋构型，这一形态对应界面局部的强剪切效应. $ {\textit{t}}^{*}=1.0 $时，相较于$ {\textit{t}}^{*}=0.5 $时刻涡量高值区集中在相界面投影区，此时涡量高值区从涡核中心向四周扩散；同时，代表涡线的黑色线圈底部缠绕加剧，这一特征对应角尖碰撞后的涡线缠绕过程，呈现出局部紊乱的形态. $ {\textit{t}}^{*}=1.5 $时，气液界面上的涡量高值区进一步向外扩展，这一变化与图5（a）中相界面在气液界面处形成的2个卷吸洞口的扩张过程直接对应. $ {\textit{t}}^{*}=2.0 $时，涡量高值区趋于弥散，涡线底部的缠绕结构基本消失，体现了黏性耗散对涡结构的衰减作用.

总动能作为表征流场机械能累积效应的核心物理量，其演化规律直接反映涡环与界面相互作用中的能量传递与耗散机制. 总动能的计算沿用2.1节的式(17).

如图6所示为不同We数下总动能$E_{\mathrm{k}} $随无量纲时间的演化. 初始时刻t^*=0，源于相同的初始涡环配置与诱导速度场，所有算例的总动能完全一致. 随着无量纲时间的推进，各工况动能均呈衰减趋势，但其衰减速率与We数呈现出相关性. 在低We工况$ (W e=0.24,0.28,0.34)$中，动能衰减过程较为平缓，至演化晚期t^*=2.7仍保留初始动能的60%以上. 这一特征与表面张力在该参数范围内的主导作用密切相关. 低We下涡量场稳定且无明显破碎，流体运动的剪切效应较弱，动能向黏性耗散的转化效率低，因此动能得以维持在较高水平. 相比之下，高We工况$(W e=0.42,0.56,0.84) $下的动能衰减速率加快，且We数越大，衰减越剧烈 . 在We=0.84时，t^*=2.7处剩余的初始动能已在25%以下. 进一步观察发现，高We动能曲线在$t^* \approx 0.6 $处出现了明显的衰减加剧现象. 这与可视化结果中“高We界面在$t^* \approx 0.6 $发生角尖接触与拓扑重联”的特征时刻相对应. 在界面重联过程中，涡结构发生剧烈破碎，强剪切效应引发了黏性耗散的激增. 同时，界面变形做功消耗了大量机械能，导致动能进入快速衰减阶段.

如图7所示，展示了2种不同韦伯数（$ We=0.54 $和$ We=0.24 $）工况中，局部动能密度分布随无量纲时间$ {t}^{*} $的演化过程，其中${k_{\mathrm{e}}} = 0.5\rho {\mkern 1mu} |{{\boldsymbol{u}}}{|^2}$表征单位体积流体所具有的动能. 在初始时刻$ {t}^{*}=0 $，2种工况下的动能密度均集中于初始液滴的核心区域，表现为一个紧凑的高值区. 随着时间推进，高$ We $工况下的液滴因惯性力主导，发生剧烈变形与破碎，局部动能密度高值区迅速向界面及周围流场扩散，并在$ {t}^{*}=3.0 $时伴随涡环结构的形成呈现出明显的空间分离；而在低$ We $工况下，表面张力效应占优，液滴的变形受到抑制，局部动能密度始终保持在液滴内部，演化过程更为平缓.

耗散率表征机械能向热能转化的强度，其演化特征与涡结构破碎程度、界面剪切作用直接相关. 流动中的黏性耗散率定义为

(18)$\varepsilon(t)=\int_V 2 \mu(\boldsymbol{x}) \boldsymbol{S}: \boldsymbol{S} \mathrm{d} V.$

式中：$: $表示张量的双点积收缩运算.

(19)$\boldsymbol{S}=\left(\nabla \boldsymbol{u}+(\nabla \boldsymbol{u})^{\mathrm{T}}\right) / 2 .$

如图8所示为不同${We} $数下平均耗散率随无量纲时间的演化. 可以看出，黏性耗散整体随时间呈持续递减的趋势. 初始时刻涡环具有较强的剪切运动，对应较高的耗散速率；随着涡环向界面迁移、拉伸并逐渐破碎，涡结构的尺度逐步细化，但剪切强度随动能衰减而持续降低，因此黏性耗散速率也随着呈整体递减趋势.

$ We $数对耗散速率的影响显著：高$ We $数下，惯性力主导的涡拉伸与破碎过程更剧烈，初始耗散速率显著高于低$ We $数工况下的；但由于高$ We $数下动能衰减更快（见图6），耗散速率的递减幅度也更明显；而低$ We $数下，表面张力约束了涡与界面的相互作用强度，涡结构演化更平缓，耗散速率维持在较低水平并缓慢衰减.

值得注意的是，尽管在固定雷诺数$ Re=1\;000 $条件下，即流体黏性系数保持不变，但$ We $数对耗散率演化仍表现出调控作用. 这是因为$ We $数决定了界面的微观几何特征：高$ We $数下界面剧烈变形产生的细丝与液滴结构极大地增强了局部的速度梯度$ \mathbf{\nabla }\boldsymbol{u} $. 根据牛顿流体耗散机制$ \epsilon \;\mathbf{\propto }\;\mu {(\mathbf{\nabla }\boldsymbol{u})}^{2} $，在$ \mu $恒定的情况下，速度梯度的增大直接提升了机械能向热能转化的效率. 这解释了为何高$ We $工况下动能衰减速率（图6）高于低$ We $工况下的.

半域总螺旋度$ {H}_{{\mathrm{l}}} $（左半域总螺旋度）和$ {H}_{{\mathrm{r}}} $（右半域总螺旋度）作为表征半空间涡结构拓扑完整性与螺旋特性的关键参数^[10]，其演化规律直接反映$ We $数对涡动力学特性的调控. 螺旋度的守恒性与涡结构的连续性密切相关，低扰动下紧凑涡结构可维持螺旋度稳定，而涡丝断裂、拓扑重联会导致螺旋度显著损失^[33].

总螺旋度可表示为

(20)$ h = {{\boldsymbol{u}}}\cdot {{\boldsymbol{\omega }}},\; H(t) = \int_{V} {h} {\mkern 1mu} {\text{d}}V. $

式中：涡量$ \boldsymbol{\omega }=\nabla \times \boldsymbol{u} $，$ h $为螺旋度密度. 由于初始设置与流场演化均关于$ \textit{x=}\text{0.5}\textit{L} $平面对称，为了研究流场在空间上的局部螺旋特征，将计算域沿$ \textit{x=}\text{0.5}\textit{L} $平面划分为左右2个子区域$ {V}_{\rm{l}} $与$ {V}_{{\mathrm{r}}} $，并分别定义半域螺旋度：

(21)$ {H_{\rm{l}}}(t) = \int_{{V_{\rm{l}}}} h {\mkern 1mu} {\text{d}}V, \;\;{H_{\mathrm{r}}}(t) = \int_{{{V}_{{\mathrm{r}}}}} h {\mkern 1mu} {\text{d}}V. $

如图9所示呈现了不同$ We $数下计算域的左半域与右半域总螺旋度的无量纲演化特征. 整体流场以 $ \textit{x=}0.5\textit{L} $为明确对称轴，左右两侧的螺旋度呈现正负对称分布，即一侧为正值时另一侧必为绝对值近似相等的负值，这种对称抵消效应使得整个计算域的总螺旋度维持在 0 附近，符合无外部力矩输入时流体系统总螺旋度守恒的基本流体力学规律. 同时，不同$ We $数对半域螺旋度的演化过程产生显著影响. 在高$ We $数工况下，半域螺旋度呈现先小幅增长后持续衰减的特征，初始阶段涡环的定向迁移与拉伸作用会短暂强化局部涡结构的旋拧程度，推动半域螺旋度出现小幅上升，而随着涡结构逐渐破碎，旋拧形态趋于分散，螺旋度随之逐步衰减并趋近于 0. 在低$ We $数工况下，表面张力的约束作用显著抑制了涡结构的旋拧与拉伸效应，半域螺旋度未出现明显的增长阶段，仅伴随涡结构的黏性耗散呈现平缓递减的演化趋势，始终维持在较低水平.

如图10所示，展示了2种不同韦伯数$ We=0.54 $和$ We=0.24 $）工况中，局部螺旋度密度分布随无量纲时间$ {t}^{*} $的演化过程，在初始演化阶段($ {t}^{*}=0.5 $)，2种工况下的螺旋度密度均集中于液滴与周围流体的剪切界面，呈现出正负交替的双极结构. 在演化的中后期，低$ We $下强大的表面张力会试图抚平界面. 这种“抚平”并不是瞬间弹回，而是将垂直方向的深坑转化为水平方向的浅坑，即扁平化 . 当一个坑变浅变宽时，它在$ y=0.5L $这个固定高度水平面上的截交线，自然会向外扩张，看起来像是在沿水平方向上被“拉伸”得较长.

归一化表面能与界面面积、表面张力系数直接相关，其演化特征反映了$ We $数对界面变形程度的调控. 表面能的表达式为

(22)$ {E}_{\sigma }=\sigma A . $

式中：A为界面总表面积. A按如下体积分数形式估算（常用于 VOF 方法中）：

(23)$ A = \int_V | {\nabla }\varphi |{\mkern 1mu} {\text{d}}V. $

以初始表面能为基准归一化处理后，有利于消除初始条件差异对结果对比的影响，聚焦界面面积变化的主导作用.

如图11所示为不同$ {We} $数下归一化表面能随无量纲时间的演化. 可以看出，在低$ We $工况（$ We\leq 0.34 $）中，归一化表面能始终维持在 1.0 附近小幅波动（波动幅度 < 5%），表明表面张力主导下，相界面保持紧凑环状结构，界面面积变化极小，表面能处于动态平衡状态. 尽管低$ We $对应较高的$ \sigma $，但界面约束作用抑制了面积扩张，因此表面能未出现显著升高，体现了$ We $数对表面能演化的主导调控作用. 高$ We $工况（$ We\geq 0.42 $）下，归一化表面能呈现“先升后降”的特征，与界面演化阶段较为同步：第1峰值出现在$ {t}^{*}\approx 0.6 $，对应界面角尖接触与拓扑重联过程，界面拉伸导致面积显著增大，表面能升高；第2峰值出现在$ {t}^{*}\approx 1.5\sim 1.8 $，对应界面局部断裂与子结构生成（图4中期演化），新界面的形成使面积再次增大；$ {t}^{*}\geq 2.4 $后表面能逐渐下降，源于断裂后的子界面相互融合或趋于稳定. 值得注意的是，表面张力系数$ \sigma $的增大（对应$ We $降低）并未导致表面能升高，反而因界面约束作用抑制了面积扩张. 此外，在高$ We $工况下，表面能的峰值高度随$ We $增大而升高，进一步印证了$ We $数越大，界面变形越剧烈，界面面积变化越显著.

基于PLIC-VOF界面捕捉方法与 Rosenhead-Moore 涡丝初场构造，通过直接数值模拟求解不可压Navier-Stokes方程（覆盖 $ We=0.24\sim 0.84 $），研究自由表面马蹄涡与自由界面耦合演化中由韦伯数主导的动力学差异. 研究表明，韦伯数作为惯性力与表面张力竞争关系的控制量，在所研究条件下对界面形态演化的路径与涡场响应起主导调控作用.

在低韦伯数的工况下，表面张力对界面起到约束与平滑作用，界面保持相对紧致，涡结构以整体迁移和缓慢衰减为主，能量输运与耗散呈现温和演化；当韦伯数增大时，惯性力主导的驱动作用使界面发生强烈拉伸、断裂与拓扑重联，同时伴随涡结构的显著碎裂与尺度级转移，导致剪切增强、耗散激增与螺旋度的持续丧失. 上述规律从界面可视化、总动能与耗散率、半域螺旋度以及归一化表面能等多种物理量的演变曲线中得到一致印证.

本研究也存在局限性，研究范围主要集中在 $ We=0.24\sim 0.84 $ 这一过渡区间. 对于更宽的参数范围，可基于本研究结果进行合理推测：当$ We \lt 0.24 $时，极强的表面张力将使界面近似于弹性边界，涡环可能仅引起界面波动而无法发生角尖卷吸重联；当$ We \gt 0.84 $时，流场将进入飞溅与雾化区域，界面破碎将生成大量微尺度液滴，此时流动将由液滴动力学主导. 未来工作将探索更宽的韦伯数范围与不同初始条件或边界条件，结合更高阶数的数值模拟格式和更精细的界面捕捉方法，以及引入实验观测数据，来确定结论的普适性，并发展适用于工程设计的简化判据，为跨尺度工程问题的调控提供支撑.