CN110287506A - 一种液压泵柱塞副温度场的流固热耦合分析方法 - Google Patents

Abstract

本发明公开了一种液压泵柱塞副温度场的流固热耦合分析方法，提出一种合理的数值模拟法求解能量方程，最终模拟出柱塞副稳态时温度场的分布情况；首先分析柱塞副的油膜特性，然后进行温度场的分析与计算，包括整体流程、边界条件的确定、流体和固体边界热传导的计算方式、理论模型的确定和模型的求解等；最终基于MATLAB实现该求解过程所得的结果，过程中考虑了边界向外传热的情况，因此更贴近实际工况。

Description

一种液压泵柱塞副温度场的流固热耦合分析方法

技术领域

本发明属于流体传动与控制技术领域，具体涉及一种液压泵柱塞副温度场的流固热耦合 分析方法。

背景技术

液压泵在工业领域具有广泛的应用，其性能、可靠性和寿命直接影响系统相关特性。影 响液压泵可靠性的主要失效模式有老化、疲劳和磨损三种，油液的温度对柱塞泵部件的失效 有着直接影响。柱塞-转子副(柱塞副)是液压柱塞泵主要的运动副之一，柱塞副油膜特性分 析问题涉及流场、结构运动和受力(结构)、和温度场分析，是典型的流固热耦合问题，国内 外针对这个问题有实验测量法、经验公式法、解析法和数值模拟法四种分析方法，现阶段最 有研究意义的为数值模拟法。温度场分析是油膜特性分析的主要环节，对于掌握液压泵运动 副失效机理，进而针对典型故障模式(如过热、磨损、密封件老化导致的泄漏等)采取有效 的控制措施，保证其长期工作的稳定性和可靠性具有至关重要的意义。

温度场的分析涉及两个方面：(1)理论分析模型的建立。国内外研究者针对温度场分析 的理论模型大都基于能量方程，加入了针对具体问题的分析，确定了相应的理论模型，由于 侧重点不同，理论分析模型也是多种多样。(2)模型的解算方法。温度场的计算模型通常利 用数值模拟法计算其近似解，常用方法有限差分法、有限元法和边界元法。

确定理论模型时，适当地将问题简化，但尽可能考虑多种因素，逼近实际工况。对于柱 塞副这种简单几何形状中的换热问题的模型解算，有限差分法最易实施，也有足够的准确性。 但能量方程本身是一个非常复杂的能量方程，数值法解算结果很难收敛并合理，因此国内外 还没有很好地解决这个问题的办法。

发明内容

有鉴于此，本发明的目的是提供一种液压泵柱塞副温度场的流固热耦合分析方法，能模 拟出柱塞副稳态时温度场的分布情况。

一种液压泵柱塞副温度场的流固热耦合分析方法，包括如下步骤：

步骤1、确定柱塞温度场的理论模型，具体为：

S11、确定柱塞副的油膜半径方向两个边界的边界温度T_solid；油膜轴向两个边界的边界温 度T_oil，以及柱塞腔出口处边界温度T_ref；

S12、柱塞副间隙油膜的能量方程表示成公式(2)所示的形式：

式中：ρ₀表示油液的密度；c₀表示油液的比热容；k为油液的热导率；u,v,w分别为油液 在θ,r,z方向的速度；T为油膜温度分布；Φ_D为热源项，又称耗散函数；R_p表示柱塞底部半径； q为泵的流量；△T为泵出入口油液的温差；V为泵内油液总体积；

S13、柱塞副油膜边界的热传导模型如公式(4)所示：

T＝T₀-λ(T₀-T_a)△l (4)

式中：T₀为初始油膜温度；T_a为与油膜发生热传导的外界物体的温度；λ为热传导系数， 其大小根据发生热传导的材料而定；△l为热传导方向每个网格的长度；

步骤2、对步骤1确定的能量方程进行无量纲化，得到无量纲的能量方程，并写成标准 形式；然后用近似形式代替偏导形式，确定离散方程中的各项系数，通过迭代得到仅考虑边 界绝热的无量纲温度分布；之后加入边界传热的计算，得到最终的温度场分布，具体步骤如 下：

S21、各物理量与其无量纲量之间的对应关系如公式组(5)所示：

式中：T^*,H,Θ,R,Z,U,W分别为T,h,θ,r,z,u,v的无量纲形式；w_p为轴向运动速度的无量纲化的常 量，R_c表示柱塞腔半径；ω表示柱塞泵主轴转动角速度；L_t为t时刻柱塞在柱塞腔内部分的长 度；

将无量纲形式代入能量方程(2)中，得到无量纲的能量方程，如公式(6)所示：

S22、求解温度分别对θ,r,z偏导的近似形式：

温度在θ方向偏导的近似形式，如公式组(7)所示：

温度对r方向偏导的近似形式如公式组(8)所示：

温度对z方向是偏导的近似形式，如公式组(9)所示：

其中，将θ方向等分为n_θ份，r方向等分为n_r份，z方向等分为n_z份，则计算域被划分为 n_θ×n_r×n_z个小格，i,j,k分别表示θ,r,z方向小格的位置；T^* _i,j,k表示i,j,k位置上的温度；

S23、将求解域分为9种情况，分别得到每种情况下能量方程的数值形式，具体为：

1.当2≤i≤n_θ+1,1<j<n_r,1<k<n_z时，公式(6)表示的能量方程最终整理得如公式(10)所示 的形式：

式中各迭代系数分别表示如公式(11)-(17)所示：

2.当2≤i≤n_θ+1,j＝1,1<k<n_z时，公式(6)表示的能量方程最终整理得如公式(18)所示的 形式：

其中，迭代系数I₁,I₂,J₂,K₁,K₂,F分别如公式(12),(13),(14),(15),(16),(17)所示，迭代 系数A,J₁分别如公式(19),(20)所示：

3.当2≤i≤n_θ+1,j＝n_r,1<k<n_z时，公式(6)表示的能量方程最终整理得如公式(21)所示的 形式。

其中，迭代系数A,I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(19)、(12)、(13)、(20)、(14)、(15)、 (16)和(17)所示：

4.当2≤i≤n_θ+1,1<j<n_r,k＝1时,公式(6)表示的能量方程最终整理得如公式(22)所示的形 式：

其中，迭代系数I₁,I₂,J₁,J₂,K₁,F分别如公式(12),(13),(14),(14),(15),(17)所示，迭代 系数A,K₂分别如公式(23),(24)所示：

5.当2≤i≤n_θ+1,1<j<n_r,k＝n_z时，公式(6)表示的能量方程最终整理得如公式(25)所示的 形式：

其中，迭代系数I₁,I₂,J₁,J₂,K₂,F分别如公式(12),(13),(14),(14),(24),(17)所示，迭代 系数A,K₁分别如公式(26),(27)所示：

6.当2≤i≤n_θ+1,j＝1,k＝1时，公式(6)表示的能量方程最终整理得如公式(28)所示的形式。

其中，迭代系数I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(12),(13),(20),(14),(15),(24),(17)所 示，迭代系数A如公式(29)所示：

7.当2≤i≤n_θ+1,j＝1,k＝n_z时，公式(6)表示的能量方程最终整理得如公式(30)所示的形 式：

其中，迭代系数I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(12),(13),(20),(14),(27),(24),(17)所 示，迭代系数A如公式(31)所示：

8.当2≤i≤n_θ+1,j＝n_r,k＝1时，公式(6)表示的能量方程最终整理得如公式(32)所示的形 式：

其中，迭代系数A,I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式 (29),(12),(13),(20),(14),(15),(24),(17)所示：

9.当2≤i≤n_θ-1,j＝n_r,k＝n_z时，最终整理得如公式(33)所示的形式：

其中，迭代系数A,I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式 (31),(12),(13),(20),(14),(27),(24),(17)所示：

采用数值模拟法分别对于以上九种情况下的能量方程的数值形式进行迭代计算，每次迭 代均遍历求解域内每个格点，求出绝热条件下的温度场分布；

S24、判断求解域内各格点属于径向边界还是轴向边界上的格点，如果属于径向边界，将 该边界上格点求解得到的温度按照公式(34)进行更新：

式中：表示径向边界上各格点的温度初值；λ_solid为油膜与固体之间的传热系数；

如果属于轴向边界，将该边界上格点求解得到的温度按照公式(35)进行更新：

式中：表示轴向边界上各格点的温度初值，λ_oil为油膜与外界油液之间的传热系数；

S25、根据设定的迭代结束条件，结束S23和S24的迭代计算，得到最终的温度场分布。

较佳的，所述迭代结束条件为：

其中，ε为所设定的相对精度标准，取10^-2～10^-6；T^* _i,j,k ⁿ⁺¹表示本次迭代计算得到的温度 的滤波结果，计算方法为：对于数值模拟法每次迭代计算出的温度，进行滤波处理，即：将 本次迭代结果T^* _i,j,k与上一次的滤波结果T^* _i,j,k ⁿ进行加权相加，作为第n+1次的滤波结果，如公 式(36)所示。

T^* _i,j,k ⁿ⁺¹＝γT^* _i,j,k+(1-γ)T^* _i,j,k ⁿ (36)

式中：γ为加权系数，取值范围为0～1。

本发明具有如下有益效果：

本发明的一种液压泵柱塞副温度场的流固热耦合分析方法，提出一种合理的数值模拟法 求解能量方程，最终模拟出柱塞副稳态时温度场的分布情况；首先分析柱塞副的油膜特性， 然后进行温度场的分析与计算，包括整体流程、边界条件的确定、流体和固体边界热传导的 计算方式、理论模型的确定和模型的求解等；最终基于MATLAB实现该求解过程所得的结果， 过程中考虑了边界向外传热的情况，因此更贴近实际工况。

附图说明

图1为柱塞分布示意图；

图2为柱塞结构示意图；

图3为本发明的模型求解方案图；

图4为本发明的模型求解流程图；

图5为柱塞底部中心相对偏心距变化曲线；

图6为柱塞相对倾覆角度变化曲线；

图7为柱塞副油膜厚度分布；

图8为柱塞副油膜压力场分布；

图9为柱塞副油膜温度场分布。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

柱塞泵一般有多个柱塞，各个柱塞的球头通过滑靴固定在斜盘上，且均匀分布，如图1 所示。每个柱塞做相同的周期性运动，柱塞-转子副是指柱塞与缸体之间所形成的运动副。柱 塞与缸体共同随主轴转动而转动，同时，由于斜盘有一定倾角，柱塞还在柱塞腔内做往复运 动，即柱塞与缸体之间在共同转动的同时还存在相对运动。因此柱塞与柱塞腔不可能完全压 紧，存在很小的缝隙，充满油膜，这部分油膜的温度特性直接影响到柱塞泵的性能、磨损状 况及寿命的预测，因此柱塞副的温度场分析至关重要。

图1中，R_s＝34.55mm为柱塞的分布圆半径；θ_s为柱塞在分布圆上的位置，取值范围为 (0,2π)其中分布圆最下面的点定为参考点，即θ_s0＝0；L₀＝45.5mm为柱塞长度；L_s为柱塞球头 和滑靴的总长度，大小约为柱塞长度的一半L_s＝L₀/2；R_p＝9.4mm为柱塞底部的半径；R_c＝9.41mm为柱塞腔的半径；β为斜盘倾角，一般为12°或15°，本发明取15°。

不考虑时间差，每个柱塞的运动是完全相同的，因此选取其中一个柱塞进行温度场分析 即可。由于受力不均，柱塞会在柱塞腔内倾覆，倾覆角度非常小，但由于油膜本身很薄，这 个角度不能忽略。图2显示了柱塞和柱塞腔之间的结构和位置关系。

图2中，e为柱塞底部中心的偏心距；α为柱塞轴线与柱塞腔轴线之间的夹角，即柱塞倾 覆角度。

以柱塞的底部中心为原点，柱塞底部圆周逆时针方向为θ方向，柱塞半径方向为r方向， 柱塞轴线方向为z方向，建立柱坐标系。图2中阴影部分所示为柱塞副油膜分析的计算域， 包括柱塞与柱塞腔之间的间隙油膜以及两个固体与油膜接触的表面。

油膜温度场分析考虑流固热耦合影响，涉及到流体、固体和热场，其中，流体部分为柱 塞与缸体之间的油膜，固体部分为柱塞在柱塞腔内部分的外壁和缸体的内壁，热场的计算除 油膜部分外，还需要考虑油膜边界处的热交换；具体的分析过程如下：

步骤1、确定柱塞温度场的理论模型，具体为：

S11、边界条件

油膜半径方向两个边界分别与柱塞和缸体接触，因此两个固体端面的温度为T_solid＝25℃； 轴向的两个边界分别与柱塞腔内油液和泵内油液接触，柱塞底部边界温度设置为T_oil＝40℃， 柱塞腔出口处边界温度设置为T_ref＝25℃。

S12、边界绝热的流固热耦合模型

能量方程(1)是分析计算热量传递过程的基本方程之一，通常表述为：流体微元的内能增 量等于通过热传导进入微元体的热量、微元体中产生的热量及周围流体对微元体所做功之和。 该方程假设流体与固体的边界绝热，只考虑了流体的热传导、热对流及热源的情况。

式中：ρ₀＝775kg/m³表示油液的密度；c₀＝2100J/(kg·K)表示油液的比热容； k＝0.145W/(m·K)为油液的热导率；u,v,w分别为油液在θ,r,z方向的速度；T为油膜温度分布；Φ_D为热源项，又称耗散函数。

针对柱塞副间隙油膜稳态情况下的温度场分布分析与计算，可做假设：(1)油膜处于热 稳定状态，即各物理状态不随时间变化而变化；(2)比热容、密度和热导率均为常数；(3) 由于膜厚方向的尺寸相比于其余两个方向小得多，因此可以忽略膜厚方向的热对流。将柱塞 副间隙油膜的能量方程表示成公式(2)所示的形式。

式(2)中，R_p表示柱塞底部半径；左侧整项表示热对流，右侧第一项为热传导项，第二项 为热源项，又称耗散函数。热源项一般根据油液的流动速度求得，考虑摩擦生热的情况，但 实际情况中，热源不只来源于摩擦生热，还有很多其它形式的热量耗散。因此做如下假设： 假设柱塞泵出入口油液内能之差为柱塞泵整体的热源，而该热源主要由柱塞副、滑靴副和配 流盘副三个运动副分担，其中柱塞副占其中的1/3。在该假设的条件下，热源的计算如公式(3) 所示。

式中：q＝0.145W/(m·K)为泵的流量；△T＝20℃为泵出入口油液的温差；V＝1L为泵内油 液总体积。

S13、边界传热模型

能量方程(1)假设油膜边界均绝热，但实际工况中，柱塞副油膜边界不可能绝热，因此需 考虑油膜边界处与外界的热传导。针对边界的热传导所建立的模型如公式(4)所示。

T＝T₀-λ(T₀-T_a)△l (4)

式中：T₀为初始油膜温度；T_a为与油膜发生热传导的外界物体的温度；λ为热传导系数， 其大小根据发生热传导的材料而定；△l为热传导方向每个网格的长度。

步骤2、对温度场分析理论模型进行求解是油膜特性分析中非常重要的步骤，同时也是 油膜特性仿真的难点。针对所建立的理论模型，利用有限差分法提出一种基于能量方程(2)的 求解方法，如图3和图4所示：首先，对能量方程(1)进行无量纲化，得到无量纲的能量方程， 并写成标准形式；然后用近似形式代替偏导形式，确定离散方程中的各项系数，通过迭代得 到仅考虑边界绝热的无量纲温度分布；之后加入边界传热的计算，并进行加权处理，得到一 次迭代后的无量纲温度分布；比较前后两次的迭代结果，若满足迭代精度要求，则该温度场 即为最终结果，否则继续迭代；最后将所得结果进行量纲化，得到最终的温度场分布，具体 步骤如下：

S21、无量纲化是数值模拟法求解微分方程很重要的一步，其目的是将所有物理量转化为 无量纲的量，且均归于同一数量级，有效地减小了计算误差。无量纲化的形式不一，理论上 不会影响最终计算结果。本发明中各物理量与其无量纲量之间的对应关系如公式组(5)所示。

式中：T^*,H,Θ,R,Z,U,W分别为T,h,θ,r,z,u,v的无量纲形式；H₀为初始油膜厚度；L_t为t时刻 柱塞在柱塞腔内部分的长度。w_p为轴向运动速度的无量纲化的常量，取w_p＝1m/s。

将无量纲形式代入能量方程(2)中，得到无量纲的能量方程，如公式(6)所示。

上式中，k₀表示油液的热导率,为常数；

S22、求解温度分别对θ,r,z偏导的近似形式：

数值模拟法求解微分方程有多种方法，例如有限差分法、有限元法、体积元法等，其中 最常用的就是有限差分法，既简单也能基本满足计算精度要求。有限差分法即将计算域划分 为有限个网格，利用差分法将微分形式表达成差商形式，中心差分形式是最常用的差分格式。 对于柱塞副来讲，以中心差分格式为基础，再结合各方向边界特点，可总结出其近似形式。

数值模拟法的核心在于计算域的离散化，将θ方向等分为n_θ份，r方向等分为n_r份，z方 向等分为n_z份，这样计算域被划分为了n_θ×n_r×n_z个小格，用i,j,k分别表示θ,r,z方向小格的位 置，△θ,△r,△z分别表示θ,r,z方向小格的长度，则将各偏导项的近似形式总结如下。

θ方向是圆周方向，因此i＝1和i＝n_θ可看作一个格点，则i＝n_θ+1就是i＝2所对应的物理状 态的值，i＝0就是i＝n_θ-1所对应的物理状态的值，据此得出温度在θ方向偏导的近似形式， 如公式组(7)所示。

注：为了减少计算情况划分，便于计算，在MATLAB中编程时，在温度场的圆周方向两 个边界外又人为加上一层网格，即温度场求解域用2≤i≤n_θ-1表示，i＝1为i＝n_θ的值，i＝n_θ+2 为i＝3的值，i＝1和i＝n_θ+2所代表的两个面没有实际意义，只是为了方便计算。

r方向是半径，对温度偏导的近似形式如公式组(8)所示。

z方向是柱塞轴线方向，对温度偏导的近似形式，如公式组(9)所示。

S23、迭代计算过程

数值法求解微分方程的过程就是一个迭代计算的过程，不断迭代直至达到收敛条件，得 到温度场分布。求解温度场分析模型首先要解算能量方程，最主要的就是确定能量方程的迭 代系数。将公式(7)至(9)中的差商形式代入公式(6)，代替偏导形式，然后整理为可以求解的数 值形式。由于数值求解不能超出所规定的格点范围，因此，将求解域分为以下9种情况进行处 理。

1.当2≤i≤n_θ+1,1<j<n_r,1<k<n_z时，代表油膜内部(即不处于边界上)各格点，公式(6) 表示的能量方程最终整理得如公式(10)所示的形式。

式中各迭代系数分别表示如公式(11)-(17)所示。

2.当2≤i≤n_θ+1,j＝1,1<k<n_z时，油膜与柱塞接触的边界上的各点，公式(6)表示的能量方 程最终整理得如公式(18)所示的形式。

其中，迭代系数I₁,I₂,J₂,K₁,K₂,F分别如公式(12),(13),(14),(15),(16),(17)所示，迭代系数A,J₁分 别如公式(19),(20)所示。

其中，△θ,△R,△Z分别表示θ,r,z方向小格的长度；

3.当2≤i≤n_θ+1,j＝n_r,1<k<n_z时，油膜与柱塞腔内壁接触的边界上的各点，公式(6)表示 的能量方程最终整理得如公式(21)所示的形式。

其中，迭代系数A,I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(19),(12),(13),(20),(14),(15),(16),(17)所示。

4.当2≤i≤n_θ+1,1<j<n_r,k＝1时，代表油膜与柱塞底部油液的边界上各点，公式(6)表示的 能量方程最终整理得如公式(22)所示的形式。

其中，迭代系数I₁,I₂,J₁,J₂,K₁,F分别如公式(12),(13),(14),(14),(15),(17)所示，迭代系数A,K₂分别 如公式(23),(24)所示。

5.当2≤i≤n_θ+1,1<j<n_r,k＝n_z时，代表油膜与刚出柱塞腔时油液的边界上各点，公式(6) 表示的能量方程最终整理得如公式(25)所示的形式。

其中，迭代系数I₁,I₂,J₁,J₂,K₂,F分别如公式(12),(13),(14),(14),(24),(17)所示，迭代系数A,K₁分别如公式(26),(27)所示。

6.当2≤i≤n_θ+1,j＝1,k＝1时，公式(6)表示的能量方程最终整理得如公式(28)所示的形式。

其中，迭代系数I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(12),(13),(20),(14),(15),(24),(17)所示，迭代 系数A如公式(29)所示。

7.当2≤i≤n_θ+1,j＝1,k＝n_z时，公式(6)表示的能量方程最终整理得如公式(30)所示的形式。

其中，迭代系数I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(12),(13),(20),(14),(27),(24),(17)所示，迭代 系数A如公式(31)所示。

8.当2≤i≤n_θ+1,j＝n_r,k＝1时，公式(6)表示的能量方程最终整理得如公式(32)所示的形式。

其中，迭代系数A,I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(29),(12),(13),(20),(14),(15),(24),(17)所示。

9.当2≤i≤n_θ-1,j＝n_r,k＝n_z时，最终整理得如公式(33)所示的形式。

其中，迭代系数A,I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(31),(12),(13),(20),(14),(27),(24),(17)所示。

采用数值模拟法对于以上九种情况下的能量方程(公式(10)、(18)、(21)、(22)、(25)、 (28)、(30)、(32)、(33))进行迭代计算，每次迭代可以遍历计算域内每个格点，求出绝热 条件下的温度场分布。

S24、求解能量方程时确定的迭代系数仅考虑了边界绝热时，为了更贴近实际工况，还要 考虑边界向外传热的情况。

油膜在圆周方向形成一个闭环，不存在与外界的热交换。

油膜在径向与固体接触，向固体进行热传导，因此径向边界进行传热后的无量纲温度计 算方法如公式(34)所示。

式中：表示径向边界上各格点的温度初值；λ_solid为油膜与固体之间的传热系数。

油膜在轴向与外界油液接触，向外界油液进行热传导，轴向边界进行传热后的温度计算 方法以公式(35)表示。

式中：表示轴向边界上各格点的温度初值，λ_oil为油膜与外界油液之间的传热系数。

对公式(10)、(18)、(21)、(22)、(25)、(28)、(30)、(32)和(33)中求解得到在径 向边界和轴向边界上格点的温度T^* _i,j,k分别利用式(34)和(35)进行修正，由此将两个方向 的边界传热均考虑在内，求出考虑边界传热的温度场分布。

S25、确定迭代结束条件：

由于数值模拟法是近似求解，需要不断迭代以逼近最终结果，因此得到整个求解域的温 度分布之后，为了保证计算结果的收敛性，通常还需采用松弛方法。对于数值模拟法每次迭 代计算出的温度，进行滤波处理，即：将本次迭代结果T^* _i,j,k与上一次的滤波结果T^* _i,j,k ⁿ进行 加权相加，作为第n+1次的滤波结果，如公式(36)所示。

式中：γ为加权系数，取值范围为0～1，通常取0.3。

每次迭代计算，油膜的温度场分布都会更新，并向最终解收敛。理论上，通过不断的迭 代计算，最终可以得到稳定的温度场分布，即每次得到的温度场分布不再随迭代次数变化而 变化，则认为此时的温度场即为能量方程(1)的解。但实际计算中，由于有一定误差，一般会 设置适当的迭代精度，允许一定的容错率，当第n+1次的迭代结果与第n次的迭代结果的差满 足该精度，则得到最终解。具体判断标准如公式(37)所示。

式中：ε为所设定的相对精度标准，通常取10^-2～10^-6，这里取10^-3。

当先后两次迭代所得温度场分布满足公式(37)，此时的温度场分布即为温度分布计算模型 解得的无量纲的温度场分布。再根据无量纲化时所用关系，将其量纲化，求得实际的温度场 分布。

实施例：

针对所建立的柱塞副的倾覆模型，基于MATLAB对柱塞副油膜特性的计算过程进行实 现，完成了计算机对油膜状态的仿真模拟，得出一系列仿真数据结果，包括动态平衡的过程 以及稳态时的油膜状态。

1、动态平衡

图5显示了柱塞副中柱塞底部中心偏离初始位置的距离随时间的变化曲线，纵轴为柱塞 底部中心相对偏心距，横轴为时间。其中，正值表示柱塞底部远离分布圆圆心，负值表示柱 塞底部靠近分布圆圆心，则可看出，柱塞底部几乎始终靠向分布圆圆心。而且柱塞的偏心运 动可视作周期性运动，始终在负向最大偏心值附近波动，总体来看，可达到动态平衡。该平 衡并不严格，这可能是由于控制方程为非线性复杂微分方程，求解过程中运用了很多近似方 法，会造成一定误差，但总体趋势符合柱塞的运动情况。

图6展示了柱塞的倾覆角度随时间的变化曲线，纵轴为柱塞相对倾覆角度，横轴为时间。 由图可知，柱塞的倾覆角度变化情况也基本可以看作周期性变化，达到一种动态平衡柱塞泵 开始运行后，柱塞会迅速达到最大倾覆角度，然后在最大倾覆角度附近不断波动。

综合图5和图6来看，柱塞底部偏心运动及柱塞产生倾覆的动态平衡过程有一定联系， 变化周期及趋势均基本一致。另外，由于结构的限制，柱塞底部会在柱塞腔内最大偏心位置 不断波动，且柱塞的倾覆角度也会在最大角度附近不断波动，两者互相影响、相辅相成。

2、油膜状态

在此动态平衡过程中，可认为柱塞达到了相对稳定状态。图7展示了稳态时某时刻柱塞 与柱塞腔之间油膜厚度的分布，z轴是油膜厚度，x轴为油膜θ方向的位置，y轴为油膜r方 向的位置。图7中油膜厚度分布沿θ轴对称，这是由于圆周方向对称位置结构与工况均相同， 则油膜状态相同；对于z轴不同的点的油膜厚度与底部偏心距和倾覆角度直接相关，在有倾 覆的情况下，一端油膜厚度较小，另一端油膜厚度较大，中间各格点油膜厚度必由小到大分 布，与实际情况相符。

图8是柱塞泵运行达到近似稳定状态后某时刻的压力场分布，轴为压力值，轴为油膜方 向的位置，轴为油膜方向的位置。图8中显示的压力分布与图7中的油膜厚度分布是相对应 的：油膜厚度越薄，油膜压力值越大。

油膜厚度最薄的点对应于压力分布三个峰值。当油膜厚度较厚时，油膜压力变化并不十 分明显，这也是由于挤压效应降低所致，符合实际情况。

图9是柱塞泵运行达到近似稳定状态后的柱塞副油膜在柱塞某一径向圆柱截面上的温度 场分布，z轴为温度值，x轴为油膜θ方向的位置，y轴为油膜z方向的位置。图9中显示油 膜温度场分布的变化趋势与油膜厚度分布相关油膜厚度较薄的位置温度也会较高，可能由于 摩擦产生的热量更多所致，与实际相符。

3、结论

本发明着重研究柱塞副内部油膜的特性，寻求一套利用计算机实现的油膜仿真计算方法， 指导柱塞泵的设计与改进。通过仿真分析得到以下结论：

1)油膜的形成均是由于两个固体之间存在相对运动，而油膜的作用除润滑之外，还会防 止干摩擦导致的固体磨损严重。

2)固体的运动一般分为轴向运动及倾覆。这两者看似没有关联，实际上都会达到相同的 动态平衡，即两者的周期性变化、变化的趋势及产生突变的时间节点均相同。

3)油膜的压力场、温度场分布均与油膜厚度分布息息相关，一般来讲，油膜厚度越薄的 位置，压力和温度值越大，会产生局部高压高温。

4)油膜的温度场分布一般在边界位置会较低，而油膜内部出现局部高温。这是由于油膜 内部由于摩擦及其他原因产生的热量不能及时传导出去导致的，而边界与外界进行热量交换 更加方便，温度不会过高。

针对柱塞副进行油膜特性仿真的结果从理论上分析是有参考价值的，可基本得到油膜特 性的变化趋势。但是其准确性及是否可完全代替实验方法指导柱塞泵的设计与改进仍须经实 验方法进一步验证。

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

Claims (2)

1.一种液压泵柱塞副温度场的流固热耦合分析方法，其特征在于，包括如下步骤：

步骤1、确定柱塞温度场的理论模型，具体为：

S11、确定柱塞副的油膜半径方向两个边界的边界温度T_solid；油膜轴向两个边界的边界温度T_oil，以及柱塞腔出口处边界温度T_ref；

S12、柱塞副间隙油膜的能量方程表示成公式(2)所示的形式：

式中：ρ⁰表示油液的密度；c⁰表示油液的比热容；k为油液的热导率；u,v,w分别为油液在θ,r,z方向的速度；T为油膜温度分布；Φ_D为热源项，又称耗散函数；R_p表示柱塞底部半径；q为泵的流量；△T为泵出入口油液的温差；V为泵内油液总体积；

S13、柱塞副油膜边界的热传导模型如公式(4)所示：

T＝T₀-λ(T₀-T_a)△l (4)

式中：T₀为初始油膜温度；T_a为与油膜发生热传导的外界物体的温度；λ为热传导系数，其大小根据发生热传导的材料而定；△l为热传导方向每个网格的长度；

步骤2、对步骤1确定的能量方程进行无量纲化，得到无量纲的能量方程，并写成标准形式；然后用近似形式代替偏导形式，确定离散方程中的各项系数，通过迭代得到仅考虑边界绝热的无量纲温度分布；之后加入边界传热的计算，得到最终的温度场分布，具体步骤如下：

S21、各物理量与其无量纲量之间的对应关系如公式组(5)所示：

式中：T^*,H,Θ,R,Z,U,W分别为T,h,θ,r,z,u,v的无量纲形式；w_p为轴向运动速度的无量纲化的常量，R_c表示柱塞腔半径；ω表示柱塞泵主轴转动角速度；L_t为t时刻柱塞在柱塞腔内部分的长度；

将无量纲形式代入能量方程(2)中，得到无量纲的能量方程，如公式(6)所示：

S22、求解温度分别对θ,r,z偏导的近似形式：

温度在θ方向偏导的近似形式，如公式组(7)所示：

温度对r方向偏导的近似形式如公式组(8)所示：

温度对z方向是偏导的近似形式，如公式组(9)所示：

其中，将θ方向等分为n_θ份，r方向等分为n_r份，z方向等分为n_z份，则计算域被划分为n_θ×n_r×n_z个小格，i,j,k分别表示θ,r,z方向小格的位置；T^* _i,j,k表示i,j,k位置上的温度；

S23、将求解域分为9种情况，分别得到每种情况下能量方程的数值形式，具体为：

1）当2≤i≤n_θ+1,1<j<n_r,1<k<n_z时，公式(6)表示的能量方程最终整理得如公式(10)所示的形式：

式中各迭代系数分别表示如公式(11)-(17)所示：

2）当2≤i≤n_θ+1,j＝1,1<k<n_z时，公式(6)表示的能量方程最终整理得如公式(18)所示的形式：

其中，迭代系数I₁,I₂,J₂,K₁,K₂,F分别如公式(12),(13),(14),(15),(16),(17)所示，迭代系数A,J₁分别如公式(19),(20)所示：

3）当2≤i≤n_θ+1,j＝n_r,1<k<n_z时，公式(6)表示的能量方程最终整理得如公式(21)所示的形式。

其中，迭代系数A,I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(19)、(12)、(13)、(20)、(14)、(15)、(16)和(17)所示：

4）当2≤i≤n_θ+1,1<j<n_r,k＝1时,公式(6)表示的能量方程最终整理得如公式(22)所示的形式：

其中，迭代系数I₁,I₂,J₁,J₂,K₁,F分别如公式(12),(13),(14),(14),(15),(17)所示，迭代系数A,K₂分别如公式(23),(24)所示：

5）当2≤i≤n_θ+1,1<j<n_r,k＝n_z时，公式(6)表示的能量方程最终整理得如公式(25)所示的形式：

其中，迭代系数I₁,I₂,J₁,J₂,K₂,F分别如公式(12),(13),(14),(14),(24),(17)所示，迭代系数A,K₁分别如公式(26),(27)所示：

6）当2≤i≤n_θ+1,j＝1,k＝1时，公式(6)表示的能量方程最终整理得如公式(28)所示的形式。

其中，迭代系数I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(12),(13),(20),(14),(15),(24),(17)所示，迭代系数A如公式(29)所示：

7）当2≤i≤n_θ+1,j＝1,k＝n_z时，公式(6)表示的能量方程最终整理得如公式(30)所示的形式：

其中，迭代系数I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(12),(13),(20),(14),(27),(24),(17)所示，迭代系数A如公式(31)所示：

8）当2≤i≤n_θ+1,j＝n_r,k＝1时，公式(6)表示的能量方程最终整理得如公式(32)所示的形式：

其中，迭代系数A,I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(29),(12),(13),(20),(14),(15),(24),(17)所示：

9）当2≤i≤n_θ-1,j＝n_r,k＝n_z时，最终整理得如公式(33)所示的形式：

其中，迭代系数A,I₁,I₂,J₁,J₂,K₁,K₂,F分别如公式(31),(12),(13),(20),(14),(27),(24),(17)所示：

采用数值模拟法分别对于以上九种情况下的能量方程的数值形式进行迭代计算，每次迭代均遍历求解域内每个格点，求出绝热条件下的温度场分布；

S24、判断求解域内各格点属于径向边界还是轴向边界上的格点，如果属于径向边界，将该边界上格点求解得到的温度按照公式(34)进行更新：

式中：表示径向边界上各格点的温度初值；λ_solid为油膜与固体之间的传热系数；

如果属于轴向边界，将该边界上格点求解得到的温度按照公式(35)进行更新：

式中：表示轴向边界上各格点的温度初值，λ_oil为油膜与外界油液之间的传热系数；

S25、根据设定的迭代结束条件，结束S23和S24的迭代计算，得到最终的温度场分布。

2.如权利要求1所述的一种液压泵柱塞副温度场的流固热耦合分析方法，其特征在于，所述迭代结束条件为：

其中，ε为所设定的相对精度标准，取10^-2～10^-6；T^* _i,j,k ⁿ⁺¹表示本次迭代计算得到的温度的滤波结果，计算方法为：对于数值模拟法每次迭代计算出的温度，进行滤波处理，即：将本次迭代结果T^* _i,j,k与上一次的滤波结果T^* _i,j,k ⁿ进行加权相加，作为第n+1次的滤波结果，如公式(36)所示。

T^* _i,j,k ⁿ⁺¹＝γT^* _i,j,k+(1-γ)T^* _i,j,k ⁿ (36)

式中：γ为加权系数，取值范围为0～1。