声空化条件下传动液体积弹性模量的时空演变

陈益宏，许沧粟，李孝禄，李运堂，陈源，金杰，庞文

(1.中国计量大学机电工程学院，浙江杭州 310018；2.浙江省智能制造质量大数据溯源与应用重点实验室，浙江杭州 310018；3.浙江大学能源工程系，浙江杭州 310027；4.中广核太阳能德令哈有限公司，青海德令哈 817000)

0 前言

体积弹性模量是指液体所受压力与体积变化量的比值，表征了传动系统中传动液的抗压缩能力[1]。体积弹性模量的值越大，传动液的抗压缩性能力越强。当系统处于高速运作的启停环节时，体积弹性模量会直接影响传动系统的动态特性，甚至会造成机器零部件的损坏[2]。因此，需引入有效体积弹性模量来准确反映系统实际工况。而有效体积弹性模量主要受传动液工作压力、含气率和温度的影响[3-4]。由于空化效应、低压和高温会导致传动管中的气体含量升高，有效体积弹性模量变小，从而对压力精准传递产生影响[5]。

传动液压力和含气率的时空演变过程较为复杂。CATANIA等[6]使用均相流模型，避免了计算气相和液相之间的滑移速度，比较了预测数据与实验数据(利用Moehwald-Bosch MEP2000-CA4000高性能试验台)，结果表明：即使在空化引起的不连续情况下，也可以高精度地模拟压力波传播。CATANIA等[7]还利用流体控制方程，考虑压缩效应引起的温度变化，提出了一种新的正压空化模型来预测压力波动、含气率和温度的变化，并与实验结果进行了比较，结果证实了所研究的模型在模拟声空化开始和结束时的有效性和鲁棒性。在此基础上，FERRARI[8]提出一种将Zielke层流频率相关摩擦模型扩展到湍流的原始方法，并将其应用于动量平衡方程，结果表明：所提出的模型能正确地预测压力波振幅相对于时间的衰减。

传动液压力和含气率的预测对研究体积弹性模量的时空演变起到关键作用。魏超等人[9]为提升油液压缩与膨胀过程动态特性分析的准确性，将含气量变化纳入有效体积弹性模量的动态模型中，并对比了4种稳态模型(Wylie、Nykanen、Ruan和AMESim模型)在高、低压时的适用性和变化特性，结果表明：在Henry定律下，动态模型预测的精度更高。唐东林等[10]推导了含气油液有效体积弹性模量理论模型，并将理论模型的数值模拟结果与实验数据对比，结果表明：在油液压力小于空气分离压时，有效弹性模量随初始含气量升高而降低；在大于空气分离压的范围内，有效弹性模量随升压时间的增加而降低。李孝禄等[11-12]基于有效体积弹性模量理论模型，利用高速摄像机获得不同类别流体(泡状流、弹状流、塞状流和环状流)的状态图，通过图像处理分析了含气率和有效弹性模量之间的关系，结果表明：制动液的有效弹性模量随含气率的升高而降低，且在含气率大于1%时，为满足制动效果需提高其制动压力。袁晓明等[13]基于完全空化模型，以改进的Henry定律提出了更高精度的动态体积弹性模量预测模型(Model 1)，与3种稳态模型和Sakama动态模型进行对比，并和实验数据进行比较，结果表明：Model 1动态模型与实验结果更接近，提高了预测的准确性。

上述学者仅在时间域上对体积弹性模量进行了预测，并未在空间上进行预测。因此，本文作者考虑空气空化、蒸汽空化和伪空化3种空化效应和由压缩效应引起的温度变化，分别建立了传动管中体积弹性模量均相流模型和动态模型，利用气液混合物的流体状态方程来连接压力项，并结合Roe格式分解和Steger Warming通量分裂方法，提出一种新的数值求解方法来预测不同空化区中的压力波传播和瞬态气体含量的变化，并预测体积弹性模量的时空演变。讨论压力、含气率和温度对动态体积弹性模量的影响，并比较两种模型的区别。

1 数学模型

1.1 均相流模型

在传动管中，考虑传动管内压缩前后传动液的质量守恒以及流动对管壁剪切应力的影响，以密度ρ和动量ρu作为守恒变量，采用质量守恒方程和动量守恒方程对可压缩气液两相流动进行数学建模，得到守恒型流体动力学方程组，并用矩阵的形式表达：

(1)

式中：ρ为混合物流体密度，kg/m3；u为传动液流速，m/s；p为传动液压力，Pa；t为时间变量，s；f为传动液流动阻力系数，当雷诺数Re≤2 000时，流动为层流，f=64/Re；当雷诺数Re>2 000时，流动为紊流，f=0.316/Re0.25。

为研究动态体积弹性模量与传动液在压缩膨胀过程中气体含量的关系，引入含气率的表达式：

(2)

式中：χ是指空气g或蒸汽v；α为含气率，%；ρl为传动液密度，kg/m3。

传动液的密度按照流型的不同可分别表示为均匀泡状流和分层流的函数：

(3)

理想气体的状态方程由气体实验定律和阿伏伽德罗定律给出：

P=ρχrχT

(4)

式中：r为气体常数，J/(mol·K)；T为温度，K。

对于正压流流动，通常在空化区，温度T是压力p的函数，T=T(p)。在式(3)中传动液的密度是压力p和温度T的函数，ρl=ρl(T,p)。

对于可压缩流体的一维流动，在忽略边界交换热量的前提下，温度场用能量守恒方程来表示：

(5)

式中：h为系统单位质量的焓，J/kg；Q为系统内单位质量的黏性功率耗散，J/(kg·s)。

在等熵条件下，考虑纯液体流动时液压管内的压缩带来温度的变化，将任何过程中液体的压力与温度联系起来，得到以下关系式[14]：

(6)

式中：Cp为定压比热容，J/(kg·K)；γT为等温声速，m/s；φ为热膨胀系数，1/K；ET为等温体积弹性模量，Pa。

传动液声速γ与体积弹性模量E有着密切关系，因此给出纯液体流动、空气空化和蒸汽空化时的声速表达式：

(7)

在气液两相条件下，给出体积弹性模量E与声速γ和含气率α之间的关系：

(8)

均相流模型需要增加气液混合物的状态方程，才能使得方程组封闭。此研究中，将流体看作是纯液体和一定量的空气和蒸汽的均质混合物，为求解气液两相流中的压力项，需要用气液混合物流体状态方程连接。给定一个气液质量流量比μ，那么混合物流体密度ρ、气体密度ρχ和含气率α均为压力p和温度T的函数。式(6)给出了压力与温度的关系，从而可以得出混合物流体密度和压力p、温度T之间的关系：

(9)

式中：μ为气液质量流量比，%。

通过混合物流体状态方程将压力p与混合物流体密度ρ联系起来，由于管内流动会产生压降，为此需给出压力梯度的表达式：

(10)

1.2 体积弹性模量动态模型

在实际工作过程中，有效体积弹性模量表征传动液的抗压缩性，其定义为

E=-V(dp)/(dV)

(11)

抗压缩特性主要受流体的压力、含气率和温度影响，且随着流体被压缩，其温度会增高，流体趋于膨胀，进而产生附加压力。因此，出于对上述3种因素的考虑，将纯传动液(不含空气与蒸汽)的有效体积弹性模量[15]简化为

Ec=E0+κΔp+ηΔT

(12)

式中：Ec为纯传动液的有效体积弹性模量，Pa；E0为标准状态下的有效体积弹性模量，Pa；κ为压力变化系数；η为温度变化系数。

在传动液传动过程中，由于传动液本身溶解了部分空气和传动管中存在着部分空气，在空化过程中产生的空气与蒸汽会使含气率升高，导致由压力变化引起的传动液体积改变，影响有效体积弹性模量。

考虑到纯传动液、空气和蒸汽3种成分在过程中的变化，且压力、含气率和温度随时间变化，引入动态体积弹性模量的表达式，将式(11)改写为

(13)

在标准大气压p0下的气体压缩方程为

(14)

式中：Vg0和Vv0分别为标准大气压p0下的空气体积和蒸汽体积，m3。

纯传动液的体积由以下表达式[16]给出：

Vl(t)=Vl0e(p-p0)/Ec

(15)

将式(14)、(15)关于压力p求导，代入式(13)得：

(16)

2 求解方法

(17)

采用Roe格式分解，对特征向量进行处理，得到流速u和密度ρ的平均表达式：

(18)

采用Steger-Warming通量分裂法在空间上进行流通量的离散，对时间项采用欧拉向前差分。差分格式如下：

(19)

3 数值结果与分析

3.1 传动液的物性参数

CATANIA等[7]以高压燃油喷射系统为研究对象，包括燃油泵、压力管和单弹簧喷油器，并利用ISO 4113试验油模拟柴油。在燃油泵出口和喷油器进口附近设置了一个压阻传感器，通过实验收集两端的压力数据。在此研究中，ISO 4113试验油[7]用于模拟纯传动液的密度、声速和热膨胀系数，这些参数是温度T和压力p的函数。

式(1)的初值条件为

(20)

此研究使用了图1所示的压力边界条件[7]。

图1 压力边界条件[7]

模拟中使用的其他参数如表1所示。初始气体含量设置为1.425%，包括溶解在传动液中的气体和气泡形式存在的气体。模拟中使用的网格大小为0.01 m×2.5 μs(480万个网格)。通过入口和出口压力波的压力差来模拟传动液的压缩和膨胀过程。

表1 模拟计算所需的参数

3.2 体积弹性模量的时空演变

为了清晰地表达传动液体积弹性模量的时空演变，在图2中绘制了沿时间和管道轴向方向的体积弹性模量变化的三维表示和映射投影。可以看出：体积弹性模量最大达到1 570 MPa。体积弹性模量主要受压力、含气率和温度的影响，在时空中发生变化，且在空化区域的变化尤为剧烈。在x=0～0.51 m、t=0.1～0.14 s以及t=0～0.06 s、t=0.16～0.24 s的部分区域，体积弹性模量为较低值。

图2 体积弹性模量变化的三维表示和映射投影

3.3 压力对动态体积弹性模量的影响

图3为传动液压力变化的三维表示和映射投影。可以看出：压力波在传动管中传播，在传动液初始压力为0.6 MPa的条件下，最大压力可达到约50 MPa。当t=0.1～0.14 s、x=0～0.51 m时，传动液压力显著下降，其压力接近甚至低于饱和蒸汽压，蒸汽的局部空化效应发生在传动管入口附近。随后，沿着管道轴向方向，传动液压力开始增加。从压力波动的角度来看，传动液压力的变化趋势与CATANIA等[7]的结果相同。

图3 压力变化的三维表示和映射投影

图4为均相流模型与动态模型下的动态体积弹性模量在管内不同位置随压力的变化。可以看出：当压力小于1 MPa，管内任意位置的动态体积弹性模量随着压力的变化相同；当压力在1～10 MPa内时，动态体积弹性模量随压力增大而快速增大，其变化率在不断减小；当压力大于10 MPa时，在均相流模型中，动态体积弹性模量随压力缓慢增加，而动态模型下的动态体积弹性模量最终稳定在1 290 MPa。

图4 管内不同位置的体积弹性模量随压力的变化

图5为均相流模型与动态模型下的动态体积弹性模量比较结果。可以看出：在空气空化后的高压区，均相流模型得到的动态体积弹性模量要高于动态模型。这是由于均相流模型直接以能量守恒方程得到温度变化，在这个阶段，压力从低压区快速传播至高压区，温度迅速升高；而动态模型以升温系数与温差的比例关系得到温度变化，在这个阶段温度趋于平稳。这使得纯液体的有效体积弹性模量在低压区向高压区转变的过程中会有较大差异，最终影响混合物的动态体积弹性模量。在低压区，两种模型的预测结果吻合，而在蒸汽空化后的高压区，两种模型的预测结果存在较小的差别。

3.4 含气率对动态体积弹性模量的影响

图6为含气率变化的三维表示和映射投影，其被分为3个区域。深黄色区域代表空气空化区，主要发生在整个管道的0～0.04 s和0.2～0.24 s时间间隔内。压力下降至空气分离压以下，气泡迅速产生，最大含气率约为1.425%。橙红色区域表示蒸汽空化区。这一过程主要发生在x=0～0.51 m、t=0.1～0.14 s。随着时间的推移，空化区移动到管的中部并不断收缩，蒸汽空化程度在x= 0.45 m左右达到最大，含气率最大达到1.925%，这与文献[7]中的最大含气率相近。紫色区域代表伪空化区，占据了大部分时空区域，含气率随压力增大而减小。

图6 含气率变化的三维表示和映射投影

图7为均相流模型与动态模型下的动态体积弹性模量在不同初始含气率下随压力的变化。可以看出，给定初始含气率的条件下，动态体积弹性模量随着压力的增大而增大。当压力小于1 MPa时，动态体积弹性模量变化率随初始含气率增大而减小。当压力在1～10 MPa内时，动态体积弹性模量随压力增大而快速增大，其变化率逐渐减小。当初始含气率小于5%，动态模型下的动态体积弹性模量的变化曲线要比均相流模型所得曲线更加平滑。而当初始含气率为10%，动态模型下的动态体积弹性模量变化率在5 MPa和7.5 MPa的压力下突增。当压力大于10 MPa，均相流模型中，动态体积弹性模量随压力缓慢增加。而在动态模型中，随着初始含气率的降低，动态体积弹性模量更快趋于平稳，这与文献[9]中的趋势一致。

图7 不同初始含气率下动态体积弹性模量随压力的变化

图8为均相流模型与动态模型下管内不同位置的动态体积弹性模量随含气率的变化。可知：管道入口与出口附近的动态体积弹性模量随含气率的增加而逐渐减小。当α<0.2%时，管道中间部分的动态体积弹性模量迅速降低，这是因为这些位置的压降过快。当α=0.2%～1%时，动态体积弹性模量随含气率的增加呈线性下降的趋势。当α>1%时，管内任意位置的动态体积弹性模量变化一致。由此可以看出：在伪空化区，含气率会影响动态体积弹性模量的变化率，而在空气和蒸汽空化区几乎没有影响。但在空气和蒸汽空化区的动态体积弹性模量更低，更不利于压力的传递。

图8 管内不同位置的动态体积弹性模量随含气率的变化

3.5 温度对动态体积弹性模量的影响

图9为传动液温度变化的三维表示和映射投影。温度随时间变化的最大幅值接近6 K，而CATANIA等[7]获得了约7 K的温度变化。在此研究中，未考虑空化区气液相变和外部交换热产生的热量，但模拟温度变化的最大幅值接近CATANIA等的结果。映射投影表明，正向流动中的温度增速比反向流动更快。此外，蒸汽空化区中的温度变化小于伪空化区中温度变化。

图9 温度变化的三维表示和映射投影

图10为均相流模型与动态模型下的动态体积弹性模量在不同初始温度下随压力的变化。可知，当压力小于1 MPa时，均相流模型中，不同初始温度下的动态体积弹性模量变化相差不大。在动态模型中，初始温度为293.15、353.15 K时，动态体积弹性模量呈线性增长，而初始温度为308.15、323.15、338.15 K时，动态体积弹性模量仅在0～ 0.4 MPa内呈线性增长。当压力在1～7.5 MPa内时，均相流模型中，动态体积弹性模量随初始温度增大而减小；而在动态模型中，当压力在0.4～1.0 MPa内时，初始温度在338.15～353.15 K内，动态体积弹性模量随初始温度增大而增大。当压力大于7.5 MPa时，均相流模型中，动态体积弹性模量呈缓慢线性增长；而在动态模型中，动态体积弹性模量趋于稳定。

图10 不同初始温度下体积弹性模量随压力的变化

4 结论

传动液中常含有空气和蒸汽，在压力波动的情况下传动液含气率随之发生变化。含气率的变化影响传动液动态体积弹性模量，进而影响传动系统的精准调控。为此，建立了均相流模型预测体积弹性模量的时空演化以及动态模型预测动态体积弹性模量，分析压力、含气率和温度对动态体积弹性模量的影响，并将均相流模型和动态模型的预测结果进行对比。结果表明，均相流模型适用于低压条件，而动态模型在高压区也同样适用。

(1) 当压力小于1 MPa时，动态体积弹性模量随压力增大而增大。当压力在1～10 MPa内时，动态体积弹性模量随压力增大而快速增大，其变化率不断减小。当压力大于10 MPa时，在均相流模型中，动态体积弹性模量随压力增大而呈线性缓慢增长趋势；而在动态模型中，气体将不断溶解，动态体积弹性模量随压力增大而稳定在1 290 MPa。

(2) 动态体积弹性模量随含气率增大而减小，其变化率不断减小。当初始含气率小于5%时，动态模型预测动态体积弹性模量-压力曲线相比均相流模型更平滑，而初始含气率在10%的情况下，均相流模型的预测曲线更平滑。

(3) 动态体积弹性模量随温度增大而减小，其变化率不断减小。在低压区，当初始温度为293.15 、353.15 K时，两种模型预测结果相同。当压力大于7.5 MPa，均相流模型中，动态体积弹性模量呈缓慢线性增长，而在动态模型中，动态体积弹性模量趋于稳定。