基于真实粗糙表面的交叉金属线接触热阻有限元分析法的制作方法

本发明涉及接触热阻分析领域，尤其是基于真实粗糙表面的交叉金属线接触热阻有限元分析法

背景技术：

随着科学技术的发展，接触热阻在微电子器件领域中的影响越来越大，是影响微电子器件性能的关键因素之一，因此对接触热阻的研究具有很重要的意义。目前对于接触热阻的研究，国内外学者基于实验测量总结出一系列经验公式，这些经验公式的应用范围较小。近年来，有学者提出接触的有限元分析，但是大多研究是基于光滑表面进行相关模拟，忽略了金属线表面形貌因素，获得的数据相对于显示情况准确性较差。

技术实现要素：

发明目的：为解决上述技术问题，考虑金属线表面形貌因素对接触热阻分析的影响，令实验结果更贴近实际值，本发明提供一种基于真实粗糙表面的交叉金属线接触热阻有限元分析法。

技术方案：本发明提供的技术方案为：基于真实粗糙表面的交叉金属线接触热阻有限元分析法，该方法包括步骤：

(1)用原子力显微镜扫描真实金属线的表面获得金属线表面形貌图像；将金属线表面形貌图像导入Gwyddion图像处理软件中，通过Gwyddion图像处理软件中的Roughness工具分析其表面形貌，并且得到金属线表面均方根粗糙度Ra；

(2)以金属微米线中间点为受力点，以金属微米线的一端为端点，长度方向为x轴构建挠度与作用力的关系公式：

$w=\underset{n=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{w}_{n}sin\left({\mathrm{\α}}_{n}l\right)$

其中，

α_n＝nπ/L

$w_n=\sqrt[3]{-\frac{q_n}{2}+\sqrt{\frac{{q_n}^2}{4}+\frac{p^3}{27}}}+\sqrt[3]{-\frac{q_n}{2}-\sqrt{\frac{{q_n}^2}{4}+\frac{p^3}{27}}}$

$p=\frac{4I}{A_c}$

$q_n=-\frac{8L^3}{{\mathrm{EA}}_c{n}^4{\mathrm{\π}}^4}Fsin\left({\mathrm{\α}}_{n}l\right)$

式中，l为力的作用点距离固定端的距离，w为挠度，n为正整数，L为金属微米线的长度，F为作用力，E为金属微米线的弹性模量，I为金属微米线横截面惯性矩，A_c为金属微米线的横截面积；

(3)基于步骤(1)得到的金属线表面均方根粗糙度Ra，利用MATLAB软件模拟符合高斯分布的随机粗糙表面，获取随机粗糙表面的高度坐标数据，包括步骤：

(3-1)利用MATLAB中的白噪声产生器函数randn()生成白噪声序列a(x，y)，对a(x，y)进行傅里叶变换，得到函数A(ω_x，ω_y)；

(3-2)构建指数型自相关函数：

$R(x,y)={\left(R_a\right)}^2\·\exp (-2.3\×\sqrt{{\left(\frac{x}{{\mathrm{\β}}_x}\right)}^2+{\left(\frac{y}{{\mathrm{\β}}_y}\right)}^2})$

其中，β_x、β_y分别为x方向和y方向上的自相关长度；

(3-3)计算Z(ω_x，ω_y)＝H(ω_x，ω_y)·A(ω_x，ω_y)；其中，

P(ω_x，ω_y)＝H²(ω_x，ω_y)·C

式中，P(ω_x，ω_y)为R(x，y)的功率频谱密度函数，C为预设的白噪声序列a(x，y)的功率频谱密度值；

(3-4)对Z(ω_x，ω_y)进行反傅里叶变换，得到随机粗糙表面的高度坐标沿x方向、y方向的变化关系函数z(x，y)；引入步骤(2)中确定的挠度与作用力关系公式，将函数z(x，y)修改为z(x-(w₁-w)，y)；w₁为距离端点d₁mm处的挠度；

(3-5)取x方向上的采样点K个，设为x_k，k＝[1，2，…，K]，y方向上的采样点M个，设为y_m，m＝[1，2，…，M]，根据函数z(x-(w₁-w)，y)获取每个点坐标(x_k，y_m)的高度值；

(4)根据步骤(2)中得到的挠度与作用力的关系公式，通过Ansys软件构建接触模型，该接触模型由两个长度方向相互垂直的半圆柱构成；

(5)将步骤(3)中得到的各个点坐标(x_k，y_m)及每个点坐标对应的高度值导入步骤(4)构建的触模型中，得到带有粗糙表面的接触模型；

(6)对带有粗糙表面的接触模型进行有限元网格划分；分别在带有粗糙表面的接触模型中两个半圆柱上加载温度，根据公式：计算接触热阻；式中，R_C圆为接触热阻，T₁为加载在一个半圆柱上的平均温度，T₂为加载在另一个半圆柱上的平均温度，q为热流密度。

进一步的，所述步骤(5)中粗糙表面的创建步骤为：

(5-1)基于各坐标(x_k，y_m)和每个坐标(x_k，y_m)对应的高度值构建矩阵A＝[z_k，m]，其中z_k，m表示坐标(x_k，y_m)的点的高度值；

(5-2)将矩阵A＝[z_k，m]中的每一行元素创建为一条曲线、每一列元素创建为一条曲线；

(5-3)以曲线之间的交点为切割点对每条曲线进行分段，将每四条形成一个封闭圈的曲线创建为一个曲面，创建的所有曲面即为所述粗糙表面。

进一步的，所述步骤(2)中挠度与作用力关系公式的参数w_n确定的方法为：

1)计算金属线的弯曲应变能：

金属线的弯矩M为：

金属线的弯曲应变能U_B为：

2)计算金属线的轴向应变能：

计算轴向拉伸长度Δl为：根据Δl计算金属线的轴向应变能为：

3)计算施加在金属线上的作用力所做的功U_F为：

$U_F=-\frac{8L^3}{{\mathrm{EA}}_c{n}^4{\mathrm{\π}}^4}Fsin\left({\mathrm{\α}}_{n}l\right)$

4)根据最小作用势原理：δ(U_S+U_B+U_F)＝0得到

$F\sin \left({\mathrm{\α}}_{n}l\right){\mathrm{dw}}_n-\frac{{\mathrm{EA}}_c}{8}\frac{{\mathrm{\π}}^4}{L^3}{w_n}^3{n}^4{\mathrm{dw}}_n=\frac{{\mathrm{EI\π}}^4}{2L^3}{n}^4{w}_n{\mathrm{dw}}_n$

求解可得：

$w_n=\sqrt[3]{-\frac{q_n}{2}+\sqrt{\frac{{q_n}^2}{4}+\frac{p^3}{27}}}+\sqrt[3]{-\frac{q_n}{2}-\sqrt{\frac{{q_n}^2}{4}+\frac{p^3}{27}}}$

$p=\frac{4I}{A_c},q_n=-\frac{8L^3}{{\mathrm{EA}}_c{n}^4{\mathrm{\π}}^4}Fsin\left({\mathrm{\α}}_{n}l\right).$

有益效果：与现有技术相比，本发明具有以下优势：

本发明基于真实粗糙表面进行相应的计算，相比于基于理想光滑表面所得到的结果更加准确，更加贴近真实情况；与利用二维粗糙表面分析并生成三维粗糙表面相比，本发明所提供的方法失真性更小。

附图说明

图1为实施例中的力学模型图；

图2为实施例力学模型中两端固定梁受力变形示意图；

图3为实施例中正弦模式1下的挠度曲线示意图；

图4为实施例中正弦模式2下的挠度曲线示意图；

图5为实施例中正弦模式3下的挠度曲线示意图；

图6为实施例中使用AFM扫描获得的真实形貌图；

图7为实施例中计算机模拟的粗糙表面示意图；

图8为实施例中的三维模型示意图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

(1)模型的力学分析

如图1所示为本实施例的力学模型图，该模型所采用的结构是两根两端固定的梁结构交叉接触，两根梁均为Pt线，长度为20mm，直径分别为20μm、30μm，当梁受较大的集中力作用时，其挠度的推导考虑弯曲和拉伸的耦合变形，此时，可以运用能量守恒原则来推导。

在推导挠度公式时，将两段固定梁的挠度分解成多个正弦曲线叠加而成，本实施例中图2所示的挠度曲线由图3、图4、图5所示的正弦模式1、2、3挠度曲线叠加而成。

如图2所示，假设固定梁的中点处受到集中载荷作用，以固定梁的一端为固定点，沿着固定梁长度方向的挠度可以表示w为

$w=\underset{n=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{w}_{n}sin\left({\mathrm{\α}}_{n}l\right),$

${\mathrm{\α}}_n=\frac{n\mathrm{\π}}{L}$

其中：n为正整数，w_n为正弦波的波峰值，l为力的作用点距离固定点的距离，L为金属线的长度。

在能量守恒项中，主要包括弯曲应变能、拉伸应变能和节点作用力所做的功。

1)对于弯曲应变能的分析如下：

梁的弯曲转角θ为：

弯矩M为：其中，E为材料的弹性模量，I为截面惯性矩

梁的弯曲应变能U_B为：

对挠度公式求二阶导数可得：

$\frac{d^{2}w}{{\mathrm{dl}}^2}=-w_1\frac{{\mathrm{\π}}^2}{L^2}sin\frac{\mathrm{\π}l}{L}-w_2\frac{2^2{\mathrm{\π}}^2}{L^2}sin\frac{\mathrm{\π}l}{L}-\mathrm{...}=-{\mathrm{\Σw}}_n{{\mathrm{\α}}_n}^2\sin \left({\mathrm{\α}}_{n}l\right)$

将上述的挠度公式带入弯曲应变能公式可得

$U_B=\frac{{\mathrm{EI\π}}^4}{4L^3}\underset{n=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{n}^4{w_n}^2$

2)对于轴向应变能的分析如下：

轴向拉伸长度Δl为：由挠度公式得代入上式得到

$\mathrm{\Δ}l=\frac{{\mathrm{\π}}^2}{4L}\underset{n=1}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}{w_n}^2{n}^2$

在假设图3、图4、图5中三种模式之间不耦合的情况下可得梁的拉伸应变能U_S为：

其中，A_c为截面面积

3)节点处的作用力所作的功U_F为：

$q_n=-\frac{8L^3}{{\mathrm{EA}}_c{n}^4{\mathrm{\π}}^4}Fsin\left({\mathrm{\α}}_{n}l\right)$

根据最小作用势原理：δ(U_S+U_B+U_F)＝0，对于每种模式都可以解得：

$F\sin \left({\mathrm{\α}}_{n}l\right){\mathrm{dw}}_n-\frac{{\mathrm{EA}}_c}{8}\frac{{\mathrm{\π}}^4}{L^3}{w_n}^3{n}^4{\mathrm{dw}}_n=\frac{{\mathrm{EI\π}}^4}{2L^3}{n}^4{w}_n{\mathrm{dw}}_n$

其中：“-”表示是由于拉应力做功

可以解得：

$w_n=\sqrt[3]{-\frac{q_n}{2}+\sqrt{\frac{{q_n}^2}{4}+\frac{p^3}{27}}}+\sqrt[3]{-\frac{q_n}{2}-\sqrt{\frac{{q_n}^2}{4}+\frac{p^3}{27}}}$

其中，

通过以上计算，得到了挠度与作用力的关系，在构建三维模型时，加入挠度的影响。

(2)真实粗糙表面的获取

1)利用AFM扫描金属线的表面，得到表面的真实形貌图像如图6所示。

2)将图像文件导入到Gwyddion图像处理软件中，分析其表面形貌，并且得到表面轮廓的算术平均偏差Ra

3)计算机利用MATLAB软件基于分析所得到的平均算术粗糙度计算出符合高斯分布的随机粗糙表面，并且输出表面高度坐标到txt文件，计算粗糙表面包括如下步骤：

利用MATLAB中的白噪声产生器函数randn()生成一个白噪声序列a(x，y)，并且对其进行傅里叶变换A(ω_x，ω_y)；

a)对指数型自相关函数进行离散化和傅里变换，得到信号的功率谱密度P(ω_x，ω_y)，其中：δ为表面均方根粗糙度Ra，β_x、β_y分别为x，y方向上的自相关长度；

b)确定输入矩阵a(x，y)的功率谱密度C，此处取C＝1；

c)根据公式P(ω_x，ω_y)＝H^２(ω_x，ω_y)·C求得滤波器的传递函数H(ω_x，ω_y)

d)由Z(ω_x，ω_y)＝H(ω_x，ω_y)·A(ω_x，ω_y)获得输出矩阵的傅里叶变换Z(ω_x，ω_y)

e)通过对Z(ω_x，ω_y)反傅里叶变换即可获得输出高度矩阵z(x，y)

交叉梁由于力的作用产生一定的挠度，所以导入的点阵需要根据挠度公式对坐标进行相应的修改。修改后的坐标点应为z(x-(w₁-w)，y)，w₁为距离固定端10mm处的挠度，w为距离固定端lmm处的挠度：w＝w_n|_l＝10000-y，此处n取20。

利用上述关系式，取x，y方向上的自相关长度均为10μm，Ra由Gwyddion测得为10.2μm，生成一个100*100的随机粗糙表面，并且对生成的10000个点进行坐标系转换得到柱状表面，如图7所示。

将点云分三列输出到txt文件中，第一列为x坐标，第二列为y坐标，第三列为高度矩阵所储存的值z(x-(w₁-w)，y)。

4)通过notepad文本编辑器对点云进行处理，使其变成能够被ansys读取的命令流文件。

(3)三维模型的建立

本发明中交叉金属线十字相交是对称的，故各自取关于接触处对称的半圆柱进行建模，如图8所示。

曲面的建立如下：

S1：曲面是100x100的点阵，也就是100x100的矩阵，ansys通过循环命令，将每一行的100个点创建曲线；

S2：点阵的第一列和最后一列同样创建曲线，利用LCSL命令，将上述所有曲线在交点处互相切割，对区线进行分段；

S3：通过AL命令，将每四条封闭的曲线创建为一个曲面，这样，共创建99块曲面组成整体的曲面形貌

创建半圆柱体的其他侧面以后就可以通过面创建立体模型。

(4)网格的划分

考虑到计算成本与精确度等因素，对模型进行网格划分的时候，采用自由划分与比例划分相结合的方式，。在接触处的网格较为密集，其他地方较为稀疏，这样做即可减少计算成本又可以保证精确度。网格划分如图8所示

(5)仿真计算

1)建立接触对，直径为30μm的圆柱侧面作为接触面，直径为20μm的圆柱侧面作为目标面。在初始条件下，利用修改接触处的x方向坐标使得两接触面发生初始接触的效果。

2)施加载荷和约束

对垂直于y、z方向上的面均施加自由度约束，仅允许x方向上的位移；同时，对图8中的直线A1A2进行全约束，即x、y、z三个方向上均不可发生位移。载荷F作用在B2点，利用CP命令耦合直线B1B2在x方向上的自由度。

3)求解计算结果。接触热阻R_C的计算式为：其中，T₁、T₂分别为上下圆柱体的平均温度，q为热流密度。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。