CN105975700A - 一种模拟超声空泡动力学行为的数值方法 - Google Patents

Abstract

本发明提供一种模拟超声空泡动力学行为的数值方法，具体步骤包括：一，构建控制方程，所述控制方程包括超声场中质量守恒方程与动量方程；二，构建超声辐射力方程；三，将超声辐射力方程带入动量方程中，并对控制方程进行无量纲化；四，求解无量纲化后的控制方程，提取所需的流场数据，实现对超声空泡动力学行为的模拟。该方法所模拟的超声空泡形态变化与真实超声空泡形态变化之间误差较小，可捕捉超声空泡瞬态演化规律和空泡的溃灭特性。

Description

一种模拟超声空泡动力学行为的数值方法

技术领域

本发明涉及一种模拟超声空泡动力学行为的数值方法，属于流体机械工程、多相流和计算流体力学技术领域。

背景技术

超声空化是声与液体媒介质相互作用而产生的一种非线性效应。液体内部由于声波作用导致局部压力降低到液体饱和蒸汽压时，会产生汽化现象，同时溶解于液体中的气体也会析出，形成气泡(又称空泡、空穴)。在声波的膨胀和压缩的过程中，这些气泡表现出的振荡、生长、收缩以及溃灭等动力学行为，以及由此产生的一系列物理和化学变化过程称为超声空化。超声空化引起的热效应以及机械效应等已被应用于众多领域，如医学治疗、工业污水处理、纳米材料制备等。另一方面，超声空泡溃灭产生的冲击力，造成化工设备不同程度的损伤，严重影响化工业生产的可靠性与安全性。

超声空化问题一直是流体机械工程和多相流领域的关键核心问题，由于超声空化现象的复杂性，单一空化泡的动力学机理研究一直是超声空化研究的重要方面之一。研究单一空化气泡的动力学过程不仅是研究多泡空化的起点，而且是研究整个超声空化现象的基础。目前，关于超声空泡动力学行为及其影响因素的数值研究，多基于Rayleigh-Plesset(R-P)方程开展。然而研究表明：单个气泡受超声作用会呈现多种振荡形态，即形态振荡(Shape Oscillation)、体积振荡(Volume Oscillation)、分裂振荡(SplittingOscillation)以及混乱振荡(Chaotic Oscillation)。因此，以球型溃灭假设为理论基础的R-P方程并不能有效地捕捉超声空泡溃灭的瞬态演化规律和空泡溃灭特征。所以建立一种模拟超声空泡的数值方法，以实现捕捉超声空泡瞬态演化规律和空泡溃灭特性，具有重要的现实意义。

发明内容

针对现有技术中存在的问题，本发明的目的是提供一种模拟超声空泡动力学行为的数值方法，该方法能够瞬态演化规律和空泡溃灭特性。

为了达到上述目的，本发明采用以下技术方案：

一种模拟超声空泡动力学行为的数值方法，具体步骤包括：

一，构建控制方程

所述控制方程包括超声场中质量守恒方程与动量方程，分别如式(1)和(2)所示：

$\▿\·u=0---\left(1\right)$

$\mathrm{\ρ}\left(\mathrm{\φ}\right)(\frac{\∂u}{\∂t}+u\·\▿u)=-\▿p+\mathrm{\ρ}\left(\mathrm{\φ}\right)g+\▿\·\[\mathrm{\μ}(\▿u+{\▿}^{T}u)\]+\mathrm{\σ}\mathrm{\κ}\left(F\right)\▿H_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)+F_a---\left(2\right)$

其中，u是流体质点的速度矢量，p是流体质点的压力，g是重力加速度，σ是气泡的表面张力系数，F_a为超声辐射力，κ是表面曲率，H_ε(φ)为Heaviside函数，F为流体体积，φ为相函数，ρ是混合密度，μ是混合动力学粘度，其表达式如下：

二，构建超声辐射力方程

将超声辐射力表示为F_a≈(0，F_y)：

$F_y=\frac{p_a^2}{{\mathrm{\ρ}}_l{c}_0^2}{e}^{-2\mathrm{\α}\left(\mathrm{\ω}\right)y}\[\mathrm{\α}cos2(\mathrm{\ω}t-ky)-ksin2(\mathrm{\ω}t-ky)+\mathrm{\α}\]---\left(11\right)$

$\mathrm{\α}\left(\mathrm{\ω}\right)=2{\mathrm{\ω}}^2\mathrm{\η}/\left(3{\mathrm{\ρc}}_0^3\right)---\left(10\right)$

其中，c₀为声波在液体中的传播速度，y为距离声源的距离，ω＝2πf表征圆频率，f为声频率，p_a为声压幅值，k＝ω/c₀为声波数，ρ_l为液体相密度，η＝1.002*10^-3Pa·s为粘性系数；

三，将所述F_a带入动量方程式(2)中，并对控制方程进行无量纲化；

四，求解无量纲化后的控制方程，提取所需的流场数据，实现对超声空泡动力学行为的模拟。

进一步地，本发明在步骤(4)中，针对无量纲化后的控制方程，对流项进行显示二阶迎风离散，粘性项进行中心差分格式离散，根据实际应用情况对计算流域的边界条件和初始化条件进行设置；然后提取所需的流场数据。

进一步地，本发明根据公式(6)和公式(7)分别计算流体体积F和相函数φ；

$\frac{\∂F}{\∂t}+u\·\▿F=0---\left(6\right)$

$\frac{\∂\mathrm{\φ}}{\∂t}+u\·\▿\mathrm{\φ}=0---\left(7\right)$

进一步地，根据所述公式(4)和公式(5)分别计算混合密度ρ和动力学粘度μ，

ρ＝ρ_g+(ρ_l-ρ_g)H_ε(φ) (4)

μ＝μ_g+(μ_l-μ_g)H_ε(φ) (5)

其中，ρ_l和ρ_g分别表示液体相和气体相密度，μ_l和μ_g分别表示液体相和气体相动力学粘度μ。

进一步地，本发明以超声波参数为基准，如超声频率、超声幅值等，对控制方程进行无量纲化。

采用上述的技术方案之后，本发明具有以下有益效果：

(1)本发明一种模拟超声空泡动力学行为的数值方法，在现有控制方程基础上，增加超声辐射力F_a，使得该控制方程可以在超声环境下适用，所模拟的超声空泡形态变化与真实超声空泡形态变化之间误差较小，可捕捉超声空泡瞬态演化规律和空泡的溃灭特性。

(2)本发明的一种模拟超声空泡动力学行为的数值方法，能够实现超声空泡在多影响因素(雷诺数、邦德数、欧拉数和声波数)下的数值模拟，对全面认识和分析超声空泡动力学行为具有重要的现实意义。

(3)本发明的一种模拟超声空泡动力学行为的数值方法，利用数值模拟结果指导并优化超声空泡在医学治疗、工业污水处理、纳米材料制备等应用中的工况设计，以及降低或避免造成化工设备振动及损伤，确保工业生产的可靠性与安全性。

附图说明

图1为模拟超声空泡动力学行为的数值方法的流程图；

图2为CLSVOF方法的计算流程图；

图3为超声空泡及流场区域模型示意图；

图4为超声空泡参数设置与模型验证；

图5为超声空泡瞬态演化规律和空泡溃灭特性。

具体实施方式

下面根据说明书附图和具体实施例对本发明做进一步详细描述。

如图1所示，本发明一种模拟超声空泡动力学行为的数值方法，具体步骤包括：

(1)构建控制方程：

所述控制方程包括超声场中质量守恒方程与动量方程，分别如公式(1)和(2)所示：

$\▿\·u=0---\left(1\right)$

$\mathrm{\ρ}\left(\mathrm{\φ}\right)(\frac{\∂u}{\∂t}+u\·\▿u)=-\▿p+\mathrm{\ρ}\left(\mathrm{\φ}\right)g+\▿\·\[\mathrm{\μ}(\▿u+{\▿}^{T}u)\]+\mathrm{\σ}\mathrm{\κ}\left(F\right)\▿H_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)+F_a---\left(2\right)$

其中，u是流体质点的速度矢量，p是流体质点的压力，g是重力加速度，σ是气泡的表面张力系数，F_a为超声辐射力，κ是表面曲率，H_ε(φ)为Heaviside函数，作用是将单相和混合相流动控制方程表达成统一形式：

$H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)=\left\{\begin{array}{cc}0,& \mathrm{\&phi;}<-\mathrm{\&epsiv;}\\ \frac{1}{2}\&lsqb;1+\frac{\mathrm{\&phi;}}{\mathrm{\&epsiv;}}+\frac{1}{\mathrm{\&pi;}}sin\left(\mathrm{\&pi;}\frac{\mathrm{\&phi;}}{\mathrm{\&epsiv;}}\right)\&rsqb;,& \left|\mathrm{\&phi;}\right|\&le;\mathrm{\&epsiv;}\\ 1& \mathrm{\&phi;}>\mathrm{\&epsiv;}\end{array}\right.---\left(3\right)$

其中,ε为界面厚度。

方程(2)中，ρ是混合密度，μ是混合动力学粘度，其表达式如下：

ρ＝ρ_g+(ρ_l-ρ_g)H_ε(φ) (4)

μ＝μ_g+(μ_l-μ_g)H_ε(φ) (5)

其中，l和g分别表示液体相和气体相。

方程(2)中，F为流体体积，φ为相函数。

VOF方法规定两相流中的一种相态(通常为液态)为“目标相态”，定义一个相函数F(VOF函数)，每个单元上的F函数定义为目标相态所占体积与单元总体积的比值。对于不可压缩流动，根据质量守恒定律推得VOF函数对流输运方程为：

$\frac{\&part;F}{\&part;t}+u\&CenterDot;\&dtri;F=0---\left(6\right)$

Level Set方程中相函数φ(x,0)为x到界面的符号距离函数，为保证相函数φ的零等值面为相界面，在任意时刻，对于相界面上的点，均要求相函数φ值为零，所以LevelSet函数对流输运方程可以写为：

$\frac{\&part;\mathrm{\&phi;}}{\&part;t}+u\&CenterDot;\&dtri;\mathrm{\&phi;}=0---\left(7\right)$

CLSVOF方法计算两相流的主要过程包括相函数初始化、流动控制方程求解、相函数对流运输方程求解、相界面重构以及φ函数再次重新初始化等。通过相界面重构和φ函数再次重新初始化等过程，使数值求解过程中的φ函数满足质量守恒特性，从而更加精确地计算表面张力。CLSVOF方法的具体实施流程如图2所示。

(2)构建超声辐射力方程。

方程(2)中的F_a表征由于声波作用产生的，作用在单位体积元上的平均声辐射力。当超声辐射力近似认为由平面波造成时，其沿声束轴线的表达式为F_a≈(0，F_y)：

$F_y=-\frac{1}{c_0}\frac{\&part;I}{\&part;y}---\left(8\right)$

其中，I＝ρc₀v²(y,t)为瞬时声强，c₀为声波在液体中的传播速度，y为距离声源的距离。质点的竖直速度分量为：

$v(y,t)=\frac{p_a}{{\mathrm{\&rho;}}_l{c}_0}{e}^{-\mathrm{\&alpha;}\left(\mathrm{\&omega;}\right)y}cos(\mathrm{\&omega;}t-ky)---\left(9\right)$

其中，ω＝2πf表征圆频率，f为声频率；p_a为声压幅值；k＝ω/c₀为声波数。由斯托克斯定律得衰减项表达式为：

$\mathrm{\&alpha;}\left(\mathrm{\&omega;}\right)=2{\mathrm{\&omega;}}^2\mathrm{\&eta;}/\left(3{\mathrm{\&rho;c}}_0^3\right)---\left(10\right)$

其中，η＝1.002*10^-3Pa·s为粘性系数；

联立式(8)-(10)最终得有耗散效应的超声辐射力的瞬态表达式为：

$F_y=\frac{p_a^2}{{\mathrm{\&rho;}}_l{c}_0^2}{e}^{-2\mathrm{\&alpha;}\left(\mathrm{\&omega;}\right)y}\&lsqb;\mathrm{\&alpha;}cos2(\mathrm{\&omega;}t-ky)-ksin2(\mathrm{\&omega;}t-ky)+\mathrm{\&alpha;}\&rsqb;---\left(11\right)$

(3)对控制方程(2)进行无量纲化：

$x^{*}=\frac{x}{D};u^{*}=\frac{u}{Df};{\mathrm{\&tau;}}^{*}=\frac{t}{f^{-1}};{\mathrm{\&rho;}}^{*}=\frac{\mathrm{\&rho;}}{{\mathrm{\&rho;}}_1};P*=\frac{P}{{\mathrm{\&rho;}}_1{D}^2{f}^2};{\mathrm{\&mu;}}^{*}=\frac{\mathrm{\&mu;}}{{\mathrm{\&mu;}}_1};k^{*}=\frac{k}{D^{-1}};g^{*}=\frac{g}{{\mathrm{Df}}^2};{\mathrm{\&alpha;}}^{*}=\frac{\mathrm{\&alpha;}}{D^{-1}}$

其中D表示有效直径的气泡，定义为D＝(6V_b/π)^1/3，V_b是气泡体积。

无量纲动量方程方程重新表示为:

$\frac{\&part;{\mathrm{\&rho;}}^{*}{u}^{*}}{\&part;{\mathrm{\&tau;}}^{*}}+\&dtri;{\mathrm{\&rho;}}^{*}{u}^{*}{u}^{*}=-\&dtri;p^{*}+{\mathrm{\&rho;}}^{*}{g}^{*}+\frac{1}{\mathrm{Re}}\&dtri;\&CenterDot;\&lsqb;{\mathrm{\&mu;}}^{*}(\&dtri;u^{*}+{\&dtri;}^T{u}^{*})\&rsqb;+\frac{1}{Bo}{\mathrm{\&kappa;}}^{*}n\&dtri;H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)+{F_a}^{*}---\left(13\right)$

$F_y={\mathrm{Eu}}^2{e}^{-2{\mathrm{\&alpha;}}^{*}{y}^{*}}\&lsqb;{\mathrm{\&alpha;}}^{*}cos2(2{\mathrm{\&pi;\&tau;}}^{*}-k^{*}{y}^{*})-k^{*}sin2(2{\mathrm{\&pi;\&tau;}}^{*}-k^{*}{y}^{*})+{\mathrm{\&alpha;}}^{*}\&rsqb;---\left(14\right)$

其中，影响超声单泡动力学行为的因素：雷诺数、邦德数、欧拉数以及声波数可分别表征为：

$\mathrm{Re}=\frac{D^2{\mathrm{\&rho;}}_{1}f}{{\mathrm{\&mu;}}_1}Bo=\frac{{\mathrm{\&rho;}}_1{D}^3{f}^2}{\mathrm{\&sigma;}};Eu=\frac{p_a}{{\mathrm{\&rho;}}_1{c}_{0}Df};k^{*}=\frac{k}{D^{-1}}---\left(15\right)$

数据如表1所示：

表1 参数设置

(4)对控制方程中的对流项进行显示二阶迎风离散，粘性项进行中心差分格式离散。同时，根据实际应用情况对计算流域的边界条件和初始化条件进行设置。为保证数值计算的稳定性、收敛性、同时节约计算时间，需要对时间步长进行设置。数值计算时的最大时间步长通常由粘性力、表面张力等确定。

(5)构建超声空泡及流场区域模型，如图3所示。并采用Fortran软件对上述方程编程并求解，求解完成后自动输出后缀名为dat的文件，所述后缀名为dat文件可采用Tecplot流体后处理软件进行数据分析，提取所需的流场数据，如气泡溃灭时的压力和速度。超声空泡模型验证如图4所示；超声空泡瞬态演化规律和空泡溃灭特性如图5所示。

(6)利用步骤一至五所述的一种模拟超声空泡动力学行为的数值方法的模拟结果指导并优化超声空泡在医学治疗、工业污水处理、纳米材料制备等应用中的工况设计：改变超声空泡的影响因素，确保超声空泡在实际应用中能在流体介质中产生工业所需的瞬态压力和速度。所述的改变超声空泡的实际影响因素主要指改变雷诺数、邦德数、欧拉数以及声波数。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种模拟超声空泡动力学行为的数值方法，其特征在于，具体步骤包括：

一，构建控制方程

所述控制方程包括超声场中质量守恒方程与动量方程，分别如式(1)和(2)所示：

▽·u＝0 (1)

$\mathrm{\&rho;}\left(\mathrm{\&phi;}\right)(\frac{\&part;u}{\&part;t}+u\&CenterDot;\&dtri;u)=-\&dtri;p+\mathrm{\&rho;}\left(\mathrm{\&phi;}\right)g+\&dtri;\&CenterDot;\&lsqb;\mathrm{\&mu;}(\&dtri;u+{\&dtri;}^{T}u)\&rsqb;+\mathrm{\&sigma;}\mathrm{\&kappa;}\left(F\right)\&dtri;H_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)+F_a---\left(2\right)$

其中，u是流体质点的速度矢量，p是流体质点的压力，g是重力加速度，σ是气泡的表面张力系数，F_a为超声辐射力，κ是表面曲率，H_ε(φ)为Heaviside函数，F为流体体积，φ为相函数，ρ是混合密度，μ是混合动力学粘度；

二，构建超声辐射力方程

将超声辐射力表示为F_a≈(0，F_y)：

$F_y=\frac{p_a^2}{{\mathrm{\&rho;}}_l{c}_0^2}{e}^{-2\mathrm{\&alpha;}\left(\mathrm{\&omega;}\right)y}\&lsqb;\mathrm{\&alpha;}cos2(\mathrm{\&omega;}t-ky)-ksin2(\mathrm{\&omega;}t-ky)+\mathrm{\&alpha;}\&rsqb;---\left(11\right)$

$\mathrm{\&alpha;}\left(\mathrm{\&omega;}\right)=2{\mathrm{\&omega;}}^2\mathrm{\&eta;}/\left(3{\mathrm{\&rho;c}}_0^3\right)---\left(10\right)$

其中，c₀为声波在液体中的传播速度，y为距离声源的距离，ω＝2πf表征圆频率，f为声频率，p_a为声压幅值，k＝ω/c₀为声波数，ρ_l为液体相密度，η＝1.002*10^-3Pa·s为粘性系数；

三，将所述F_a带入动量方程式(2)中，并对控制方程进行无量纲化；

四，求解无量纲化后的控制方程，提取所需的流场数据，实现对超声空泡动力学行为的模拟。

2.根据权利要求1所述模拟超声空泡动力学行为的数值方法，其特征在于，在步骤(4)中，针对无量纲化后的动量方程(2)，对流项[u·▽u]进行显示二阶迎风离散，粘性项(▽·[μ(▽u+▽^Tu)])进行中心差分格式离散，根据实际应用情况对计算流域的边界条件和初始化条件进行设置；然后提取所需的流场数据。

3.根据权利要求1所述模拟超声空泡动力学行为的数值方法，其特征在于，根据公式(6)和公式(7)分别计算流体体积F和相函数φ；

$\frac{\&part;F}{\&part;t}+u\&CenterDot;\&dtri;F=0---\left(6\right)$

$\frac{\&part;\mathrm{\&phi;}}{\&part;t}+u\&CenterDot;\&dtri;\mathrm{\&phi;}=0---\left(7\right).$

4.根据权利要求1所述模拟超声空泡动力学行为的数值方法，其特征在于，所述公式(4)和公式(5)分别计算混合密度ρ和动力学粘度μ，

ρ＝ρ_g+(ρ_l-ρ_g)H_ε(φ) (4)

μ＝μ_g+(μ_l-μ_g)H_ε(φ) (5)

其中，ρ_l和ρ_g分别表示液体相和气体相密度，μ_l和μ_g分别表示液体相和气体相动力学粘度μ。

5.根据权利要求1所述模拟超声空泡动力学行为的数值方法，其特征在于，以超声波参数为基准，对控制方程进行无量纲化。