CN110069837B - 横观各向同性多层涂层体系三维温度场的求解方法 - Google Patents

Abstract

本发明公开了一种横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场的求解方法，包括以下步骤：1)通过引入二维傅里叶积分变换在频域推导横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场的频域解析解；2)采用基于二维快速傅里叶变换的转换算法由步骤1)的频域解析解转换获得横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场分布。该方法求解速度快、精度高，具有较优的鲁棒性，适用于具有任意涂层层数和涂层厚度的横观各向同性多层涂层体系半空间表面在面分布移动热源作用下的三维温度场的求解，适用范围广。

Description

横观各向同性多层涂层体系三维温度场的求解方法

技术领域

本发明涉及表面移动摩擦热源作用下温度场模拟仿真领域，尤其涉及一种横观各向同性多层涂层体系表面在面分布移动摩擦热源作用下其热源作用微区三维温度场的模拟计算方法。

背景技术

随着航空发动机等机械动力传动系统转速和传动功率的不断提升，滚动轴承与齿轮等关键基础零部件摩擦副接触微区在表面移动摩擦热源作用下温度急剧升高，接触微区材料在高温状态下服役将出现热软化、材料微观组织恶化以及机械力学性能退化等一系列问题，此外在摩擦热源作用下如果接触微区温度升高超过材料的耐温极限，摩擦副将发生热胶合，导致机械传动系统丧失工作能力。因此求解摩擦副接触微区在表面移动摩擦热源作用下的微区温度场是评估摩擦副的服役状态的重要依据和避免出现热胶合恶性失效的关键。

现有的方法针对主要是无涂层、单层涂层、双层涂层的涂层体系或者热特性为各向同性的多层涂层体系在表面热源作用下三维稳态温度场的求解。随着材料科学和表面工程技术的发展，涂层技术已由单层涂层发展为多层复合涂层、纳米超晶格多层涂层体系，并被应用于提高航空发动机的机械传动系统摩擦副的抗磨损、抗疲劳和热胶合性能，但对于横观各向同性多层涂层体系表面在面分布移动摩擦热源作用下的三维温度场的求解尚无现成的求解方法。

发明内容

为解决现有技术中存在的问题，本发明提供一种横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场的求解方法。

为此，本发明的技术方案如下：

一种横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场的求解方法，包括以下步骤：

1)通过引入二维傅里叶积分变换在频域推导横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场的频域解析解；

2)选择一个区域作为计算域，采用基于二维快速傅里叶逆变换的转换算法，由步骤1)的频域解析解转换获得横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场分布。

上述方法中，步骤1)中的频域解析解的推导步骤如下：

步骤一、对第k层横观各向同性层状材料三维温度场的微分控制方程

/>

实施二维傅里叶积分变换

获得三维温度场微分控制方程的频域形式：

其中：

x为平行于移动热源方向的坐标，单位为m；

y为垂直于移动热源方向且平行于层状材料同性平面的坐标，m；

z_k为第k层横观各向同性层状材料垂直于同性平面的坐标，m；

ω_x为二维傅里叶积分变换与变量x对应的频域变量；

ω_y为二维傅里叶积分变换与变量y对应的频域变量；

T^(k)为温度，K；

为z方向的热传导系数，W/(m·K)；

为x方向的热传导系数，W/(m·K)；

c_k为体积比热容，J/(m³·K)；

V为热源移动速度，m/s；

i为虚数单位符号，

步骤二、求第k层层状材料三维温度场在频域的控制方程的通解：

其中：

是与ω_x和ω_y有关的待定参数，/>

步骤三、确定各层材料三维温度场微分控制方程通解的待定参数

对于基体，由于z_N+1→∞时，

所以/>

对于其它待定参数，由边界条件和各界面的连续条件建立关于各层材料三维温度场频域控制方程通解的待定参数的线性方程组：

A_{(2N+1)×(2N+1)}M_(2N+1)×1＝R_(2N+1)×1 (4)

其中：

/>

线性方程组的系数矩阵A_{(2N+1)×(2N+1)}的子矩阵分别为：

式中，h_l-1为第l-1层涂层的厚度，N为涂层层数；

线性方程组的待求变量矩阵M_(2N+1)×1的子矩阵分别为：

线性方程组的右边矩阵R_(2N+1)×1的子矩阵分别为：

其中：

为作用在多层涂层体系半空间表面的面分布移动热源Q_H(x,y)的傅里叶积分变换；

通过分析方程的系数矩阵的特殊形式推导获得关于各个待定参数的解的递推公式：

其中：

上述方法中，步骤2)的具体步骤如下：

步骤一、在任意深度z处选择一个矩形区域Ω_c＝{(x,y)|x_b≤x≤x_e,y_b≤y≤y_e}作为计算域，通常x_b＝-2a_H，x_e＝2a_H，y_b＝-2a_H，y_e＝2a_H，并把深度z处的计算域Ω_c划分为(N_x-1)×(N_y-1)个均匀网格单元，a_H为赫兹点接触的接触半径，单位为m，各单元的尺寸为Δ_x×Δ_y＝[(x_e-x_b)/(N_x-1)]×[(y_e-y_b)/(N_y-1)]，第[i,j]个单元几何中心处的温度记为T[i,j]；

步骤二、把对应频域的计算域Ω_F＝{(ω_x,ω_y)|-π/2Δx≤ω_x<π/2Δx，-π/2Δy≤ω_y<π/2Δy}划分为

个均匀网格单元，/>

E_p为频域网格细化倍数，为2的非负整数次幂，频域网格单元的大小为/>

/>

步骤三、由深度z处的频域解

计算在频域网格单元各个节点处的值：

从而构造一个具有

个元素的二维矩阵/>

步骤四、通过对二维矩阵

的元素位置进行翻转操作得到二维矩阵/>

步骤五、对二维矩阵

进行二维快速傅里叶逆变换得到新的二维矩阵T′：

步骤六、深度z处各节点的温度值T[i,j]为：

T[i,j]＝T′[i-N_x/2+1,j-N_y/2+1](N_x/2≤i≤N_x-1,N_y/2≤j≤N_y-1)，

其中涂层层数N理论上可以是任意的正整数。

本发明专利的有益效果如下：

(1)推导了热特性为横观各向同性多层涂层体系半空间表面在面分布移动热源作用下确定各层涂层频域通解待定参数的解的递推公式，获得了三维温度场在频域的封闭解析解，同时应用了二维快速傅里叶逆变换算法进行加速求解，求解速度快、精度高。

(2)此外涂层的层数N可以为任意正整数，具有较优的鲁棒性，适用于具有任意涂层层数和涂层厚度的横观各向同性多层涂层体系半空间表面在面分布移动热源作用下的三维温度场的求解，适用范围广。

附图说明

图1：横观各向同性多层涂层体系半空间表面在面分布移动热源作用下示意图；

图2：温度场求解技术路线示意图；

图3：空间计算域的网格单元划分示意图；

图4：频域的网格单元加密划分示意图；

图5：二维矩阵

的元素进行翻转操作示意图；

图6：由二维矩阵T′提取空间计算域各节点温度值示意图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细描述。

如图1所示，本发明是一种关于横观各向同性多层涂层体系半空间表面在面分布移动热源作用下稳态温度场的求解方法，图中N为涂层层数，h_k为第k层涂层的厚度，κ_k为第k层涂层材料的热传导系数，γ_k为第k层涂层材料的热扩散系数，V_s为表面热源的移动速度，Q_H(x,y)为涂层体系表面的面分布移动热源。本发明的技术路线如图2所示，其具体实施步骤如下：

步骤一、对第k层横观各向同性层状材料三维温度场的微分控制方程

实施二维傅里叶积分变换

获得三维温度场微分控制方程的频域形式：

其中：

x为平行于移动热源方向的坐标，m；

y为垂直于移动热源方向且平行于层状材料同性平面的坐标，m；

z_k为第k层横观各向同性层状材料垂直于同性平面的坐标，m；

T^(k)为温度，K；

为z方向的热传导系数，W/(m·K)；

为x方向的热传导系数，W/(m·K)；

c_k为体积比热容，J/(m³·K)；

V为热源移动速度，m/s；

ω_x为二维傅里叶积分变换与变量x对应的频域变量；

ω_y为二维傅里叶积分变换与变量y对应的频域变量；

i为虚数单位符号，

步骤二、求第k层层状材料三维温度场微分控制方程在频域的通解可得：

其中：

是与ω_x和ω_y有关的待定参数，/>

步骤三、确定各层材料三维温度场微分控制方程通解的待定参数

对于基体，由于z_N+1→∞时，

所以/>

对于其它待定参数，由表面边界条件和各界面连续条件建立关于各层材料温度控制方程频域通解的未知待定参数的线性方程组：

A_{(2N+1)×(2N+1)}M_(2N+1)×1＝R_(2N+1)×1 (4)

其中：

线性方程组的系数矩阵A_{(2N+1)×(2N+1)}的子矩阵分别为：

式中，h_l-1为第l-1层涂层的厚度，N为涂层层数；

线性方程组的待求变量矩阵M_(2N+1)×1的子矩阵分别为：

线性方程组的右边矩阵R_(2N+1)×1的子矩阵分别为：

其中：

为作用在多层涂层体系表面移动面分布热源Q_H(x,y)的傅里叶积分变换。通常摩擦热源分布可以假设为：

其二维傅里叶积分变换为：

通过分析方程的系数矩阵的特殊形式可以推导获得关于各个待定参数的解的递推公式，具体结果如下：

通过分析方程的系数矩阵的特殊形式可以推导获得关于各个待定参数的解的递推公式，具体结果如下：

/>

其中：

步骤四、选择一个矩形区域Ω_c＝{(x,y)|x_b≤x≤x_e,y_b≤y≤y_e}作为计算域，通常x_b＝-2a_H，x_e＝2a_H，y_b＝-2a_H，y_e＝2a_H，其中a_H为赫兹点接触的接触半径，单位为m。采用基于二维快速傅里叶逆变换的转换算法可以由任意深度z处的温度场的频域解转换获得其空间计算域各网格单元的温度值，其具体实过程如下：

①如图3所示，把深度z处的计算域Ω_c＝{(x,y)|x_b≤x≤x_e,y_b≤y≤y_e}划分为(N_x-1)×(N_y-1)个均匀网格单元N_x和N_y取为2的正整数次幂，单元尺寸为Δ_x×Δ_y＝[(x_e-x_b)/(N_x-1)]×[(y_e-y_b)/(N_y-1)]，第[i,j]个单元节点的温度记为T[i,j]。

②如图4所示，把对应频域的计算域Ω_F＝{(ω_x,ω_y)|-π/2Δx≤ω_x<π/2Δx，-π/2Δy≤ω_y<π/2Δy}划分为

个均匀网格单元，/>

E_p为频域网格细化倍数，为2的非负整数次幂，频域网格单元的尺寸为/>

③由深度z处的频域解

计算在频域网格单元各个节点处的值：

从而构造一个具有

个元素的二维矩阵/>

④如图5所示，通过矩阵

的元素位置进行翻转操作得到矩阵/>

即：

⑤对二维矩阵

进行二维快速傅里叶逆变换得到新的二维矩阵T′，即：

⑥如图6所示，由二维矩阵T′提取深度z处空间计算域各网格单元的温度值T[i,j]为：

T[i,j]＝T′[i-N_x/2+1,j-N_y/2+1](N_x/2≤i≤N_x-1,N_y/2≤j≤N_y-1)，

/>

/>

Claims (2)

1.一种横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场的求解方法，其特征在于包括以下步骤：

1)通过引入二维傅里叶积分变换在频域推导横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场的频域解析解；

2)选择一个区域作为计算域，采用基于二维快速傅里叶逆变换的转换算法，由步骤1)的频域解析解转换获得横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场分布；

其中，步骤1)中的频域解析解的推导步骤如下：

步骤一、对第k层横观各向同性层状材料三维温度场的微分控制方程

实施二维傅里叶积分变换

获得三维温度场微分控制方程的频域形式：

其中：

x为平行于移动热源方向的坐标，单位为m；

y为垂直于移动热源方向且平行于层状材料同性平面的坐标，m；

z_k为第k层横观各向同性层状材料垂直于同性平面的坐标，m；

ω_x为二维傅里叶积分变换与变量x对应的频域变量；

ω_y为二维傅里叶积分变换与变量y对应的频域变量；

T^(k)为温度，K；

为z方向的热传导系数，W/(m·K)；

为x方向的热传导系数，W/(m·K)；

c_k为体积比热容，J/(m³·K)；

V为热源移动速度，m/s；

i为虚数单位符号，

步骤二、求第k层层状材料三维温度场在频域的控制方程的通解：

其中：

是与ω_x和ω_y有关的待定参数，/>

步骤三、确定各层材料三维温度场微分控制方程通解的待定参数

对于基体，由于z_N+1→∞时，

所以/>

对于其它待定参数，由边界条件和各界面的连续条件建立关于各层材料三维温度场频域控制方程通解的待定参数的线性方程组：/>

A_{(2N+1)×(2N+1)}M_(2N+1)×1＝R_(2N+1)×1 (4)

其中：

线性方程组的系数矩阵A_{(2N+1)×(2N+1)}的子矩阵分别为：

式中，h_l-1为第l-1层涂层的厚度，N为涂层层数；

线性方程组的待求变量矩阵M_(2N+1)×1的子矩阵分别为：

线性方程组的右边矩阵R_(2N+1)×1的子矩阵分别为：

其中：

为作用在多层涂层体系半空间表面的面分布移动热源Q_H(x,y)的傅里叶积分变换；

通过分析方程的系数矩阵的特殊形式推导获得关于各个待定参数的解的递推公式：

其中：

/>

2.如权利要求1所述的横观各向同性多层涂层体系半空间表面在面分布移动摩擦热源作用下三维温度场的求解方法，其特征在于：步骤2)的具体步骤如下：

步骤一、在任意深度z处选择一个矩形区域Ω_c＝{(x,y)|x_b≤x≤x_e,y_b≤y≤y_e}作为计算域，x_b＝-2a_H，x_e＝2a_H，y_b＝-2a_H，y_e＝2a_H，并把深度z处的计算域Ω_c划分为(N_x-1)×(N_y-1)个均匀网格单元，a_H为赫兹点接触的接触半径，单位为m，各单元的尺寸为Δ_x×Δ_y＝[(x_e-x_b)/(N_x-1)]×[(y_e-y_b)/(N_y-1)]，第[i,j]个单元几何中心处的温度记为T[i,j]；

步骤二、把对应频域的计算域Ω_F＝{(ω_x,ω_y)|-π/2Δx≤ω_x<π/2Δx，-π/2Δy≤ω_y<π/2Δy}划分为

个均匀网格单元，/>

E_p为频域网格细化倍数，为2的非负整数次幂，频域网格单元的大小为/>

步骤三、由深度z处的频域解

计算在频域网格单元各个节点处的值：

从而构造一个具有

个元素的二维矩阵/>

步骤四、通过对二维矩阵

的元素位置进行翻转操作得到二维矩阵/>

步骤五、对二维矩阵

进行二维快速傅里叶逆变换得到新的二维矩阵T′：

步骤六、深度z处各节点的温度值T[i,j]为：

T[i,j]＝T′[i-N_x/2+1,j-N_y/2+1](N_x/2≤i≤N_x-1,N_y/2≤j≤N_y-1)，

/>

Images