本发明涉及一种模拟空泡溃灭诱发水垢空蚀的数值方法，采用Fortran软件，对构建的流场中的流体质量守恒方程与动量方程进行编程，构建空泡及流场区域模型，对构建的流场中的流体质量守恒方程与动量方程进行求解，采用Tecplot处理软件对求解的文件进行数据分析，提取气泡溃灭时的压力和形态数据，采用计算流体力学软件Ansys/workbench软件平台构建水垢的结构模型并设置水垢的结构参数，通过软件的post后处理模块中得到水垢的形变量以及应力分布云图，以水垢的形变量代替水垢的空蚀程度。与现有技术相比，本发明的优点在于：通过耦合流体控制方程与水垢控制方程，使得该数值方法可以模拟空泡在近壁面处的运动行为，所模拟的空泡形态变化与真实空泡形态变化之间误差较小。
1.一种模拟空泡溃灭诱发水垢空蚀的数值方法，其特征在于：包括如下步骤：
步骤(1)、构建流体控制方程，该流体控制方程包括流场中质量守恒方程与动量方程，分别如方程(11)和(12)所示：上述方程中ρ＝ρl+(ρl-ρg)Hε(φ)；μ＝μl+(μl-μg)Hε(φ)；
其中，为矢量微分算符，u是流体质点的速度矢量，*代表无量纲参数，为求偏导数的符号；ρ是液体和气的混合密度；φ为相函数，τ为无量纲时间，p是气泡溃灭时的液体质点的压力，g是重力加速度，Re是表征液体粘性力的雷诺数，μ是液体和气的混合动力学粘度，T表示计算的时刻；Bo为表征气泡表面张力的邦德数；κ是气泡表面曲率，F为流体体积，n为气泡法向向量，ρl是液体密度，ρg是气体密度，μl是液体动力粘度，μg是气体动力粘度，Hε(φ)为带有气液交界面厚度ε的Heaviside函数，定义为：D0为气泡的直径，t为有量纲 时间；
步骤(2)、气液交界面的捕捉方式采用CLSVOF方法；
步骤(3)、流场求解：采用Fortran软件，利用步骤(1)中构建的流场中的流体质量守恒方程与动量方程进行编程，构建空泡及流场区域模型，对方程(12)中的 进行显式二阶迎风离散， 进行中心差分格式离散；输入计算初始值：Re＝
5000；Bo＝300，ρl＝1000kg/m3，μl＝1.308*10-3Pa*s；D0＝1mm；然后对步骤(1)中构建的流场中的流体质量守恒方程与动量方程进行求解，输出后缀名为dat的文件；
步骤(4)、液-固交界面处理：将步骤(3)求解完成后输出的后缀名为dat的文件采用Tecplot处理软件进行数据分析，提取气泡溃灭时的压力和形态数据，将气泡溃灭时的压力和形态数据作为计算固定水垢变形的初始参数；
步骤(5)、构建水垢控制方程：水垢控制方程用方程(13)表示：
其中，δ，及 分别是拉格朗日水垢节点的位移、速度和加速度；σ是气泡的表面张力系数，on表示“在”；Ω表示水垢的计算域，ρs是水垢的密度；T表示计算的时刻，Γg和Γh分别是Dirichlet和Neumann边界；j和h分别是规定的位移和牵引边界条件；δ0和 分别是水垢初始状态的位移和速度向量；
在方程(13)中，通过方程(14)获得δ，及 的矩阵数据
其中M＝∫ΩNTρsNdV，K＝∫ΩBTCBdV， Γ是水垢受力面周线，N是水垢的质量矩阵，B是水垢的位移差值矩阵，p是步骤(4)中提取出的气泡溃灭时的液体质点压力数据，V表示水垢的体积，C表示水垢的刚度矩阵；
采用计算流体力学软件Ansys/workbench软件平台对方程(14)进行编程，构建水垢的结构模型并设置水垢的结构参数；
步骤(6)、设置水垢材料为碳酸钙，碳酸钙的弹性模量为70GPa，泊松比为0.3，相对密度为2.710；将步骤(4)中得到气泡溃灭时的压力和形态作为水垢空蚀计算的初始条件，求解方程(13)，求解完成后，在计算流体力学软件Ansys/workbench软件的post后处理模块中得到水垢的形变量以及应力分布云图，以水垢的形变量代替水垢的空蚀程度。
一种模拟空泡溃灭诱发水垢空蚀的数值方法
技术领域
[0001] 本发明涉及一种模拟空泡溃灭诱发水垢空蚀的数值方法。
背景技术
[0002] 工业及家庭净液体及加热设备在长时间使用后，会在液体容器表面形成致密的金属氧化膜。该金属氧化膜(即水垢，主要包含碳酸钙和碳酸镁)的形成不仅对净液体设备造成永久性损伤，而且极大地损害人体健康。利用空泡溃灭清除水垢是一种高效、环保、且无二次污染的技术手段。空泡溃灭是指液体内部局部压力降低到液体饱和蒸汽压时，液体内部或液固交界面上出现的蒸汽或气体空泡形成、生长、收缩、溃灭，以及由此产生的一系列物理和化学变化的过程。空蚀是指空泡溃灭时形成的微射流与微激波，攻击壁面形成损伤的过程。微射流除垢技术即为利用空泡溃灭时形成的微激波与微射流，及其造成的极端高温高压环境，对材料表面的杂质及污垢进行清除的过程。利用空泡溃灭并冲击壁面实现表面除垢的过程是典型的流固耦合现象。了解并掌握空泡溃灭诱发材料空蚀的流固耦合作用机理能够有效地解决微射流除垢技术的关键问题。
[0003] 空化与空蚀问题一直是微射流除垢技术领域的核心问题，由于空化与空蚀现象的复杂性，单一空泡溃灭诱发材料形变的机理研究一直是流体力学的重要研究方面之一。研究空泡溃灭诱发材料形变过程不仅是研究多泡空蚀的起点，而且是研究整个气泡-壁面流固耦合问题的基础。目前，现有的对空泡运动模型做过理论研究的技术方案中，大多都是基于假设气泡在运动发展过程中始终保持为球形的基础上进行的，然而实际上，气泡在临近液体边界时，其运动受边界条件影响，在生长、收缩、溃灭阶段不再保持球形，会出现环形溃灭，甚至出现高速射流，而高速射流是微射流除垢的重要因素；因此建立一种模拟空泡溃灭诱发水垢空蚀的数值方法，以实现捕捉空泡溃灭的瞬态演化规律，以及固体材料的应力应变过程，在液体处理行业、微射流除垢等领域具有重要的工程应用价值。
发明内容
[0004] 本发明所要解决的技术问题是针对上述现有技术提供一种模拟空泡溃灭诱发水垢空蚀的数值方法，该方法能够捕捉空泡溃灭的瞬态演化规律，同时还能捕捉固体材料的应力应变过程。
[0005] 本发明解决上述技术问题所采用的技术方案为：一种模拟空泡溃灭诱发水垢空蚀的数值方法，其特征在于：包括如下步骤：
[0006] 步骤(1)、构建流场中的无量纲控制方程，该无量纲控制方程包括流场中质量守恒方程与动量方程，分别如方程(11)和(12)所示：
[0007]
[0008]
[0009] 上述方程中ρ＝ρl+(ρl-ρg)Hε(φ)；μ＝μl+(μl-μg)Hε(φ)；
[0010] 其中，为矢量微分算符，u是流体质点的速度矢量，*代表无量纲参数，为求偏导数的符号；ρ是液体和气的混合密度；φ为相函数，τ为无量纲时间，p是气泡溃灭时的液体质点的压力，g是重力加速度，Re是表征液体粘性力的雷诺数，μ是液体和气的混合动力学粘度，T表示计算的时刻；Bo为表征气泡表面张力的邦德数；κ是气泡表面曲率，F为流体体积，n为气泡法向向量，ρl是液体密度，ρg是气体密度，μl是液体动力粘度，μg是气体动力粘度，Hε(φ)为带有气液交界面厚度ε的Heaviside函数，定义为：
[0011]
[0012]
[0013] D0为气泡的直径，t为有量纲 时间；
[0014] 步骤(2)、气液交界面的捕捉方式采用CLSVOF方法；
[0015] 步骤(3)、流场求解：采用Fortran软件，利用步骤1中构建的流场中的流体质量守恒方程与动量方程进行编程，构建空泡及流场区域模型，对公式(12)中的 进行显式二阶迎风离散， 进行中心差分格式离散；输入计算初始值：Re＝5000；Bo＝300，ρl＝1000kg/m3，μl＝1.308*10-3Pa*s；D0＝1mm；然后对步骤1中构建的流场中的流体质量守恒方程与动量方程进行求解，输出后缀名为dat的文件；
[0016] 步骤(4)、液-固交界面处理：将步骤(3)求解完成后输出的后缀名为dat的文件采用Tecplot处理软件进行数据分析，提取气泡溃灭时的压力和形态数据，将气泡溃灭时的压力和形态数据作为计算固定水垢变形的初始参数；
[0017] 步骤(5)、水垢空蚀的计算：构建水垢的结构模型，水垢的结构模型由方程(13) 表示：
[0018]
[0019] 其中，δ，及 分别是拉格朗日水垢节点的位移、速度和加速度；σ是气泡的表面张力系数，on表示“在”；Ω表示水垢的计算域，ρs是水垢的密度；T表示计算的时刻，Γg和Γh分别是Dirichlet和Neumann边界；j和h分别是规定的位移和牵引边界条件；δ0和 分别是水垢初始状态的位移和速度向量；
[0020] 在方程(13)中，通过方程(14)获得δ，及 的矩阵数据
[0021]
[0022] 其中M＝∫ΩNTρsNdV，K＝∫ΩBTCBdV， Γ是水垢受力面周线， N是水垢的质量矩阵，B是水垢的位移差值，p是步骤(4)中提取出的气泡溃灭时的液体质点压力数据，V表示水垢的体积，C表示水垢的刚度矩阵；
[0023] 采用计算流体力学软件Ansys/workbench软件平台构建水垢的结构模型并设置水垢的结构参数，设置水垢材料为碳酸钙，碳酸钙的弹性模量为70GPa，泊松比为0.3，相对密度为2.710；将步骤(4)中得到气泡溃灭时的压力和形态作为水垢空蚀计算的初始条件，求解方程(13)，求解完成后，在计算流体力学软件Ansys/workbench软件的 post后处理模块中得到水垢的形变量以及应力分布云图，以水垢的形变量代替水垢的空蚀程度。
[0024] 与现有技术相比，本发明的优点在于：
[0025] (1)本发明耦合流体控制方程与水垢控制方程，使得该数值方法可以模拟空泡在近壁面处的运动行为，所模拟的空泡形态变化与真实空泡形态变化之间误差较小，可捕捉空泡瞬态演化规律、空泡的溃灭特性以及空泡溃灭处的壁面材料的空蚀量以及应力分布；
[0026] (2)本发明能够有效地控制流体的物性参数(例如雷诺数、邦德数)以及水垢的物性参数(例如杨氏模量、刚度、密度等)，对全面认识和分析空泡溃灭诱发材料空蚀具有重要的现实意义；
[0027] (3)本发明利用数值模拟结果指导并优化水处理设备的结构参数，以实现有效地实现微射流除垢，确保工业及家庭水处理设备的安全，以及极大地提高饮用水品质。
附图说明
[0028] 图1为本发明实施例中模拟空泡溃灭诱发水垢空蚀的数值方法的流程图；
[0029] 图2为本发明实施例中气泡及流场区域模型示意图；
[0030] 图3为本发明实施例中气泡瞬态演化规律和空泡溃灭特性；
[0031] 图4为本发明实施例中水垢空蚀模型示意图。
具体实施方式
[0032] 以下结合附图实施例对本发明作进一步详细描述。
[0033] 本发明提供的模拟空泡溃灭诱发水垢空蚀的数值方法，参见图1所示，包括如下步骤：
[0034] 步骤(1)、构建流场中的无量纲控制方程，该无量纲控制方程包括流场中质量守恒方程与动量方程，分别如方程(11)和(12)所示：
[0035]
[0036]
[0037] 上述方程中ρ＝ρl+(ρl-ρg)Hε(φ)；μ＝μl+(μl-μg)Hε(φ)；
[0038] 其中，为矢量微分算符，u是流体质点的速度矢量，*代表无量纲参数，为求偏导数的符号；ρ是液体和气的混合密度；φ为相函数，τ为无量纲时间，p是气泡溃灭时的液体质点的压力，g是重力加速度，Re是表征液体粘性力的雷诺数，μ是液体和气的混合动力学粘度，T表示计算的时刻；Bo为表征气泡表面张力的邦德数；κ是气泡表面曲率，F为流体体积，n为气泡法向向量，ρl是液体密度，ρg是气体密度，μl是液体动力粘度，μg是气体动力粘度，Hε(φ)为带有气液交界面厚度ε的Heaviside函数，定义为：
[0039]
[0040]
[0041] D0为气泡的直径，t为有量纲 时间；
[0042] 步骤(2)、气液交界面的捕捉方式采用CLSVOF方法；在该界面捕捉方法中，对于不可压缩流动，相函数φ保证了对流输运方程的质量守恒，相函数φ的零等值面为相界面，同时保证了气-液交界面的顺滑；
[0043] 步骤(3)、流场求解：采用Fortran软件，利用步骤1中构建的流场中的流体质量守恒方程与动量方程进行编程，构建空泡及流场区域模型，如图2所示，对公式(12) 中的进行显式二阶迎风离散， 进行中心差分格式离散；输入计算初始值：Re＝5000；Bo＝300，ρl＝1000kg/m3，μl＝1.308*10-3Pa*s；D0＝1mm；然后对步骤1中构建的流场中的流体质量守恒方程与动量方程进行求解，输出后缀名为 dat的文件；
[0044] 步骤(4)、液-固交界面处理：将步骤(3)求解完成后输出的后缀名为dat的文件采用Tecplot处理软件进行数据分析，提取气泡溃灭时的压力和形态数据，将气泡溃灭时的压力和形态数据作为计算固定水垢变形的初始参数；气泡瞬态演化规律和气泡溃灭特性如图3所示；
[0045] 步骤(5)、水垢空蚀的计算：构建水垢的结构模型，水垢的结构模型由方程(13) 表示：
[0046]
[0047] 其中，δ， 及 分别是拉格朗日水垢节点的位移、速度和加速度；σ是气泡的表面张力系数，on表示“在”；Ω表示水垢的计算域，ρs是水垢的密度；T表示计算的时刻，Γg和Γh分别是Dirichlet和Neumann边界；j和h分别是规定的位移和牵引边界条件；δ0和 分别是水垢初始状态的位移和速度向量；
[0048] 在方程(13)中，通过方程(14)获得δ，及 的矩阵数据
[0049]
[0050] 其中M＝∫ΩNTρsNdV，， Γ是水垢受力面周线，N是水垢的质量矩阵，，B是水垢的位移差值矩阵，p是步骤(4)中提取出的气泡溃灭时的液体质点压力数据，V表示水垢的体积，C表示水垢的刚度矩阵；
[0051] 采用计算流体力学软件Ansys/workbench软件平台构建水垢的结构模型并设置水垢的结构参数，参见图4所示；设置水垢材料为碳酸钙，碳酸钙的弹性模量为70GPa，泊松比为0.3，相对密度为2.710；将步骤(4)中得到气泡溃灭时的压力和形态作为水垢空蚀计算的初始条件，求解方程(13)，求解完成后，在计算流体力学软件 Ansys/workbench软件的post后处理模块中得到水垢的形变量以及应力分布云图，以水垢的形变量代替水垢的空蚀程度。
