含偏心热源的矩形散热板双面对流散热性能的解析方法

技术领域

本发明属于电子器件的封装散热技术领域，尤其涉及一种器件热源偏心封装在矩形 散热板的上表面且上下两个表面均处于对流冷却状态下的矩形散热板散热性能的解析方 法。

背景技术

电子器件在按小型化趋势发展的同时，其电子功能也在日益提升，越来越高的芯片 集成度导致单位面积的功耗增加，芯片结温逐步攀升。为避免结温过高而缩短器件的使 用寿命和可靠度，将器件封装在肋片式散热器的表面是最常见的散热方式，并且散热器 通常采用大面积矩形散热板以承载更多数量的肋片。这种器件热源明显小于矩形散热板 的情况，也普遍存在于器件与电路板封装散热的应用中。在这种小面积热源向大面积矩 形散热板进行传导散热时，矩形散热板的热阻将包含一维传导热阻和扩散热阻，且热源 与矩形散热板的面积差距过大时，扩散热阻很可能远大于一维传导热阻，导致矩形散热 板贴装有热源一侧的上表面的温度出现显著的不均匀分布。因此，在同类电子器件封装 结构的散热设计中并不能单纯考虑容易计算的一维传导热阻，而必须考察扩散热阻的影 响，以避免低估局部最高温度而造成器件的过热失效。

在扩散热阻研究的早期，人们多以半无限尺寸的散热体为解析计算的对象，或者针 对有限厚度的矩形散热板，在其无热源的下表面施加等温边界条件后进行解析计算。Song 和Lee等人对这些早期工作进行了很好地总结，并提出以对流边界条件代替等温边界条 件来研究矩形散热板的散热性能则更符合实际应用，进而发展出一种简单的解析计算方 法。然而，上述解析计算的对象均属于热源与矩形散热板中心重合的情况，为使解析模 型更适于工程应用，Muzychka等人提出了一种含偏心热源的矩形散热板温度分布和扩散 热阻的解析计算方法，能够准确求解出矩形散热板的表面温度分布和扩散热阻的数值。

虽然Muzychka等人的方法能计算偏心热源情况，但前提是假设热源所在矩形散热板 的上表面对散热的贡献可以忽略而能施加绝热边界条件，从而回避设定为混合边界条件 时所带来的解析计算困难。然而针对矩形散热板的上表面同样参与对流散热的更为普遍 化的问题时，这种绝热边界的处理则会引入较大的计算误差。

为解决矩形散热板的上表面处于混合边界条件时给解析计算带来的困难，Luo等人 将分离变量法与等效热路分析法相结合，提出了含偏心热源的矩形散热板双面对流冷却 的解析方法，并针对同一种随机案例完成了解析解和COMSOL软件数值计算结果的对 比，验证了该解析计算方法的有效性。然而，上述方法由于运用了等效热路分析法，只 能计算出热源的平均温度和平均扩散热阻值，并不能准确得到任意指定坐标上的温度数 值，因此仍然存在低估局部最高温度的可能。

发明内容

本发明的目的在于解决上述的技术问题而提供一种适用于偏心热源矩形散热板的上 下表面均参与对流散热的更为普遍化的解析方法，能有效处理矩形散热板的上表面施加 的混合边界条件，准确计算出任意指定坐标上的温度数值，获得详细的温度分布和局部 最高温度，并能计算出扩散热阻和总热阻的最大值。

为实现上述目的，本发明采用如下技术方案：

一种含偏心热源的矩形散热板双面对流散热性能的解析方法，包括以下步骤：

(1)确定散热性能解析计算用物理参数，包括：二维热源的面积A_s＝c×d及其均匀 的发热功率Q；矩形散热板的上、下表面积A_b＝a×b及其厚度t；矩形散热板的导热系 数k；热源中心坐标(x_c,y_c,0)；矩形散热板的上表面在热源区域之外的对流换热系数h₀， 以及矩形散热板的整个下表面的对流换热系数h_t；环境温度T_f；

(2)对热源区域划分网格，将每条网格线交点作为计算节点并获得计算节点坐标；

在热源区域x和y方向上划分出网格数分别为nx和ny的均匀网格，nx和ny均为偶数；将 x和y方向上的网格步长分别记为Δx和Δy，则每个网格面积A_source-mesh＝Δx·Δy＝c/nx·d/ ny，沿x和y方向上的计算节点数分别为nx+1和ny+1；任意一计算节点的坐标为(x,y,z)， 其中(x_c–c/2)≤x≤(x_c+c/2)，(y_c–d/2)≤y≤(y_c+d/2)，z＝0。

(3)计算出热源区域中所有节点的过余温度θ(x,y,z)的初始值；节点的过余温度θ(x, y,z)初始值计算公式如下：

$\left(\begin{array}{c}\theta (x,y,z)=A_0+B_{0}z+\underset{m=1}{\overset{\infty }{\Sigma }}\cos \left({\lambda }_{m}x\right)[A_1\cosh \left({\lambda }_{m}z\right)+B_1\sinh \left({\lambda }_{m}z\right)]\\ +\underset{n=1}{\overset{\infty }{\Sigma }}\cos \left({\delta }_{n}y\right)[A_2\cosh \left({\delta }_{n}z\right)+B_2\sinh \left({\delta }_{n}z\right)]\\ +\underset{m=1}{\overset{\infty }{\Sigma }}\underset{n=1}{\overset{\infty }{\Sigma }}\cos \left({\lambda }_{m}x\right)\cos \left({\delta }_{n}y\right)[A_3\cosh \left({\beta }_{mn}z\right)+B_3\sinh \left({\beta }_{mn}z\right)]\end{array}\right)$

其中，λ_m＝mπ/a，δ_n＝nπ/b，β_mn＝(λ_m²+δ_n²)^0.5，其中m,n＝1,2,3…为无限累加变量； 式中A_i与B_i，i＝0,1,2,3，为傅立叶系数；

(4)统计出热源区域所有节点的过余温度初始值的平均值θ_mean,0，作为在后续迭代计 算中判断收敛与否的初始解；

(5)建立热源区域的热量损失矩阵Q'

式中，Q'_IJ表示热源区域中每四个相邻节点所包围的一个网格内的平均热量损失，I ＝1,…,nx；J＝1,…,ny；

(6)利用热量损失矩阵Q'，建立热源区域的非均匀热功率矩阵Q”，

(7)在步骤(6)的基础上，进行多次迭代运算，且每次迭代计算时，将热源区域的节 点坐标，以及利用步骤(6)所确定的热源区域非均匀热功率矩阵Q”计算获得的傅立叶系数 C_i、D_i分别替代傅立叶系数A_i与B_i，全部代入步骤(3)中的过余温度初始值计算公式，计 算出热源区域中所有节点的过余温度θ(x,y,z)的迭代值；

(8)统计出热源区域所有节点的过余温度迭代值的平均值θ_mean,L，L为迭代次数；

(9)根据过余温度迭代值的平均值θ_mean,L，判断迭代结果是否收敛；若收敛，则进入 步骤(10)；否则，返回步骤(5)继续增加一次迭代计算；

(10)在最后一次迭代计算热源区域过余温度的结果中，统计出热源区域内过余温度 的局部最高值θ_max(x_s,y_s,0)；

(11)按步骤(2)的网格划分方法，对矩形散热板不含热源区域的上表面、下表面划分 均匀网格，获得对应的计算节点坐标；将该计算节点坐标，以及利用最后一次迭代所确 定的热源区域非均匀热功率矩阵Q”计算获得的傅立叶系数C_i、D_i分别替代傅立叶系数 A_i与B_i，全部代入步骤(3)中的过余温度初始值计算公式，计算得到不含热源区域的上表 面、下表面上所有节点的过余温度；

(12)统计出矩形散热板不含热源区域的上表面、下表面所有计算节点的过余温度的 平均值θ_mean(x₀,y₀,0)、θ_mean(x_t,y_t,t)；

(13)利用步骤(11)中获得的节点的过余温度结果，计算矩形散热板不含热源区域的 上表面、下表面的对流散热总功率Q_up、Q_low，具体是先分别计算矩形散热板不含热源区 域的上表面、下表面的每四个相邻节点所包围的一个网格内的对流散热功率Q_up-mesh、 Q_low-mesh，再将所有计算出的网格内的对流散热功率分别求和，即为矩形散热板不含热源 区域的上表面、下表面的对流散热总功率Q_up、Q_low；

(14)利用热源区域内过余温度的局部最高值θ_max(x_s,y_s,0)，不含热源区域的上表面、 下表面所有计算节点的过余温度的平均值θ_mean(x₀,y₀,0)、θ_mean(x_t,y_t,t)，以及不含热源 区域的上表面、下表面的对流散热总功率Q_up、Q_low，计算热源向矩形散热板的上、下表 面散热的最大扩散热阻R_{s_0,max}和R_{s_t,max}；并利用热源区域内过余温度的局部最高值θ_max(x_s, y_s,0)以及二维热源的均匀的发热功率Q计算热源到环境间的最大总热阻R_th,max。

步骤(3)中的过余温度初始值计算公式中，傅立叶系数A_i，i＝0,1,2,3，分别由下式计算 得到；

$A_0=\frac{Q}{abk(\frac{h_t}{{\mathrm{th}}_t+k}+\frac{h_0}{k})}$

$A_1=\frac{2Q}{{\mathrm{abck\lambda }}_m[\phi \left({\lambda }_m\right){\lambda }_m+\frac{h_0}{k}]}[sin\left({\lambda }_m\frac{2x_c+c}{2}\right)-sin\left({\lambda }_m\frac{2x_c-c}{2}\right)]$

$A_2=\frac{2Q}{{\mathrm{abdk\delta }}_n[\phi \left({\delta }_n\right){\delta }_n+\frac{h_0}{k}]}[sin\left({\delta }_n\frac{2y_c+d}{2}\right)-sin\left({\delta }_n\frac{2y_c-d}{2}\right)]$

$\left(\begin{array}{c}{A}_3=\frac{4Q}{{\mathrm{abcdk\lambda }}_m{\delta }_n[\phi \left({\beta }_{mn}\right){\beta }_{mn}+\frac{h_0}{k}]}\\ ·[sin\left({\lambda }_m\frac{2x_c+c}{2}\right)-\sin \left({\lambda }_m\frac{2x_c-c}{2}\right)][\sin \left({\delta }_n\frac{2y_c+d}{2}\right)-sin\left({\delta }_n\frac{2y_c-d}{2}\right)]\end{array}\right)$

傅立叶系数B_i由其与傅立叶系数A_i，i＝0,1,2,3的下列对应关系式计算获得：

$B_0=-\left(\frac{h_t}{{\mathrm{th}}_t+k}\right)A_0$

B_i＝-φ(ζ)A_i,i＝1,2,3

$\phi \left(\zeta \right)=\frac{\zeta \sinh \left(\zeta t\right)+\frac{h_t}{k}\cosh \left(\zeta t\right)}{\zeta \cosh \left(\zeta t\right)+\frac{h_t}{k}\sinh \left(\zeta t\right)}$

当i取1、2或3时，ζ分别替换成λ_m、δ_n或β_mn。

步骤(5)中每四个相邻节点所包围的一个网格内的平均热量损失Q'_IJ，的计算公式为：

$\left(\begin{array}{c}{Q^{\prime }}_{LJ}=h_0{A}_{cource-mesh}\frac{1}{4}[\theta (x_0+(I-1)\Delta x,y_0+(J-1)\Delta y,0)+\theta (x_0+(I-1)\Delta x,y_0+J\Delta y,0)\\ +\theta (x_0+I\Delta x,y_0+(J-1)\Delta y,0)+\theta (x_0+I\Delta x,y_0+J\Delta y,0)]\end{array}\right)$

式中，x₀和y₀为热源区域最靠近坐标原点的计算节点在x-y平面上的坐标值，分别 为(x_c–c/2)和(y_c–d/2)；I＝1,…,nx，J＝1,…,ny；节点的过余温度θ(x,y,z)的取值分 为两种情况：若刚完成热源区域过余温度的初始值计算，则取初始值；若已实施后续的 迭代计算，则取最近一次迭代计算获得的迭代值。

步骤(7)中，在采用步骤(3)中过余温度初始值计算公式计算热源区域中所有节点的过 余温度θ(x,y,z)的迭代值时，利用步骤(6)所确定的热源区域非均匀热功率矩阵Q”计算获 得的分别替代余温度初始值计算公式中所述傅立叶系数A_i、B_i的傅立叶系数C_i、D_i， i＝0,1,2,3，分别由下式计算得到：

$C_0=\frac{{Q^{\prime \prime }}_{11}+...+{Q^{\prime \prime }}_{nxny}}{abk(\frac{h_t}{{\mathrm{th}}_t+k}+\frac{h_0}{k})}$

$\left(\begin{array}{c}{C}_1=\frac{2}{{\mathrm{abk\lambda }}_m\Delta x[\phi \left({\lambda }_m\right){\lambda }_m+\frac{h_0}{k}]}\\ ·\left\{\right[\sin \left({\lambda }_m\frac{2x_c-c+2\Delta x}{2}\right)-\sin \left({\lambda }_m\frac{2x_c-c}{2}\right)]({Q^{\prime \prime }}_{11}+\mathrm{...}+{Q^{\prime \prime }}_{1ny})+\mathrm{...}\\ +[\sin \left({\lambda }_m\frac{2x_c-c+2nx·\Delta x}{2}\right)-\sin \left({\lambda }_m\frac{2x_c-c+2(nx-1)·\Delta x}{2}\right)]({Q^{\prime \prime }}_{nx1}+\mathrm{...}+{Q^{\prime \prime }}_{nxny})\}\end{array}\right)$

$\left(\begin{array}{c}{C}_2=\frac{2}{{\mathrm{abk\delta }}_n\Delta y[\phi \left({\delta }_n\right){\delta }_n+\frac{h_0}{k}]}\\ ·\left\{\right[\sin \left({\delta }_n\frac{2y_c-d+2\Delta y}{2}\right)-\sin \left({\delta }_n\frac{2y_c-d}{2}\right)]({Q^{\prime \prime }}_{11}+\mathrm{...}+{Q^{\prime \prime }}_{nx1})+\mathrm{...}\\ +[\sin \left({\delta }_n\frac{2y_c-d+2ny·\Delta y}{2}\right)-\sin \left({\delta }_n\frac{2y_c-d+2(ny-1)·\Delta y}{2}\right)]({Q^{\prime \prime }}_{1ny}+\mathrm{...}+{Q^{\prime \prime }}_{nxny})\}\end{array}\right)$

$\left(\begin{array}{c}{C}_3=\frac{4}{{\mathrm{abk\lambda }}_m{\delta }_n\Delta x\Delta y[\phi \left({\beta }_{mn}\right){\beta }_{mn}+\frac{h_0}{k}]}\\ ·\left\{{Q^{\prime \prime }}_{11}\right[\sin \left({\lambda }_m\frac{2x_c-c+2\Delta x}{2}\right)-\sin \left({\lambda }_m\frac{2x_c-c}{2}\right)]\\ ·\sin [\left({\delta }_n\frac{2y_c-d+2\Delta y}{2}\right)-\sin \left({\delta }_n\frac{2y_c-d}{2}\right)+\mathrm{...}\\ +{Q^{\prime \prime }}_{nxny}[\sin \left({\lambda }_m\frac{2x_c-c+2nx·\Delta x}{2}\right)-\sin \left({\lambda }_m\frac{2x_c-c+2(nx-1)·\Delta x}{2}\right)]\\ ·[\sin \left({\delta }_n\frac{2y_c-d+2ny·\Delta y}{2}\right)-\sin \left({\delta }_n\frac{2y_c-d+2(ny-1)·\Delta y}{2}\right)]\}\end{array}\right)$

$D_0=-\left(\frac{h_t}{{\mathrm{th}}_t+k}\right)C_0$

D_i＝-φ(ζ)C_i,i＝1,2,3

$\phi \left(\zeta \right)=\frac{\zeta \sinh \left(\zeta t\right)+\frac{h_t}{k}\cosh \left(\zeta t\right)}{\zeta \cosh \left(\zeta t\right)+\frac{h_t}{k}\sinh \left(\zeta t\right)}$

当i取1、2或3时，ζ分别替换成λ_m、δ_n或β_mn。

步骤(9)中，迭代结果是否收敛的判据式为：

|θ_mean,L-θ_mean,L-1|<0.001，

若计算结果满足判据式的要求，则迭代结果收敛；否则，迭代结果不收敛。

步骤(11)中，在采用步骤(3)中过余温度初始值计算公式计算不含热源区域的上表面、 下表面上所有节点的过余温度时，利用最后一次迭代所确定的热源区域非均匀热功率矩 阵Q”计算获得的分别替代余温度初始值计算公式中所述傅立叶系数A_i、B_i的傅立叶系数 C_i、D_i，i＝0,1,2,3，的计算公式与步骤(7)相同。

步骤(13)中，计算矩形散热板不含热源区域的上表面上每四个相邻节点所包围的一个 网格内的对应散热功率Q_up-mesh的公式如下：

$Q_{up-mesh}=h_0{A}_{up-mesh}\frac{1}{4}({\theta }_1+{\theta }_2+{\theta }_3+{\theta }_4)$

式中，A_up-mesh为步骤(11)中划分的不含热源区域的上表面每个均匀网格的面积，θ₁、 θ₂、θ₃和θ₄为不含热源区域的上表面每四个相邻节点的过余温度；

计算矩形散热板的下表面上每四个相邻节点所包围的一个网格内的对应散热功率 Q_low-mesh的公式如下：

$Q_{low-mesh}=h_t{A}_{low-mesh}\frac{1}{4}({\theta }_{11}+{\theta }_{22}+{\theta }_{33}+{\theta }_{44})$

式中，A_low-mesh为步骤(11)中划分的下表面每个均匀网格的面积，θ₁₁、θ₂₂、θ₃₃和θ₄₄为下表面每四个相邻节点的过余温度。

步骤(14)中，所述热源向矩形散热板的上、下表面散热的最大扩散热阻R_{s_0,max}和 R_{s_t,max}，及热源到环境间的最大总热阻R_th,max的计算公式如下：

$R_{s\_0,max}=\frac{{\theta }_{max}(x_s,y_s,0)-{\theta }_{mean}(x_0,y_0,0)}{Q_{up}},$

$R_{s\_t,max}=\frac{{\theta }_{max}(x_s,y_s,0)-{\theta }_{mean}(x_t,y_t,t)}{Q_{low}},$

$R_{th,max}=\frac{{\theta }_{\max }(x_s,y_s,0)}{Q}.$

在解析过程中，步骤(3)中所述无限累加变量m和n分别预先设定一最小累加次数， 当累加次数大于设定的最小值后，判断最后一次累加结果与上一次结果的差值是否小于 预定判据值，若满足则完成本次累加计算，否则继续累加计算。

在执行步骤(2)和步骤(11)中，通过描绘热源区域过余温度的最高值θ_max(x_s,y_s,0)随网 格数量增加的变化数据图，找出数据变化趋于饱和的起点，将该起点对应的网格数定为 网格划分的最少数量。

与现有技术相比，本发明的有益效果是：

(1)本发明解析计算方法，能够针对上表面封装偏心热源的矩形散热板结构，准确地 定量计算矩形散热板的上、下表面同时参与对流散热的这一电子封装散热中更为普遍化 的问题。

(2)本发明解析计算方法所用到的计算公式，是将分离变量法与网格划分法相结合而 推导得到的，并通过多次迭代计算来保证矩形散热板上表面的对流边界条件更接近于真 实情况，从而巧妙地实现了混合边界条件的有效处理，解决了混合边界条件的解析困难 问题，且具有较高的计算精度。

(3)本发明解析计算方法能准确得到任意指定坐标上的温度数值，进而描绘出详细的 温度分布图，方便获得局部最高温度，并且能计算出封装散热系统的扩散热阻和总热阻 的最大值。因此，本发明的解析计算方法对散热问题的描述更直观，获得的数据更全面， 能为电子器件的封装散热设计工作提供更充分的分析数据。

附图说明

图1A～1B分别为俯视以及正视状态下的含偏心热源的矩形散热板双面对流散热的示 意图；

图2为含偏心热源的矩形散热板双面对流散热性能的解析方法流程图。

具体实施方式

下面，结合实例对本发明的实质性特点和优势作进一步的说明，但本发明并不局限 于所列的实施例。

参见图1～2所示，含偏心热源的矩形散热板双面对流散热性能的解析方法，包括下 述步骤：

(1)确定散热问题或性能解析用的物理参数

在解析计算含偏心热源的矩形散热板双面对流散热的性能之前，需要首先确定的物 理参数为：二维热源1的面积A_s＝c×d及其均匀的发热功率Q；矩形散热板2的上、下 表面积A_b＝a×b及其厚度t；矩形散热板的导热系数k；热源的中心坐标(x_c,y_c,0)；矩形 散热板的上表面在热源区域之外的对流换热系数h₀，以及矩形散热板的整个下表面的对 流换热系数h_t；环境温度T_f。在此，矩形散热板2的四个侧面和热源区域1向环境的散 热量均较小，可忽略不计。

(2)对热源区域划分网格并获得计算节点坐标

在热源区域的x和y方向上划分出网格数分别为nx和ny的均匀网格，且nx和ny均为偶 数；将x和y方向上的网格步长分别记为Δx和Δy，则每个网格的面积为A_source-mesh＝Δx·Δy＝ c/nx·d/ny。热源区域中每条网格线的交点即为计算节点，则沿x和y方向上的计算节点 数分别为nx+1和ny+1；任意一个计算节点的坐标记为(x,y,z)，其中(x_c–c/2)≤x≤(x_c+c /2)，(y_c–d/2)≤y≤(y_c+d/2)，z＝0。

(3)计算出热源区域中所有节点的过余温度的初始值

假设热源区域向环境存在对流换热系数为h₀的额外的对流散热量，则热源区域中任 意一个计算节点的过余温度θ(x,y,z)的初始值，均可以将其节点坐标代入到式(1)中计算 得到。

$\begin{array}{c}\theta (x,y,z)=A_0+B_{0}z+\underset{m=1}{\overset{\infty }{\Sigma }}\cos \left({\lambda }_{m}x\right)[A_1\cosh \left({\lambda }_{m}z\right)+B_1\sinh \left({\lambda }_{m}z\right)]\\ +\underset{n=1}{\overset{\infty }{\Sigma }}\cos \left({\delta }_{n}y\right)[A_2\cosh \left({\delta }_{n}z\right)+B_2\sinh \left({\delta }_{n}z\right)]\\ +\underset{m=1}{\overset{\infty }{\Sigma }}\underset{n=1}{\overset{\infty }{\Sigma }}\cos \left({\lambda }_{m}x\right)\cos \left({\delta }_{n}y\right)[A_3\cosh \left({\beta }_{mn}z\right)+B_3\sinh \left({\beta }_{mn}z\right)]\end{array}---\left(1\right)$

在此，λ_m＝mπ/a，δ_n＝nπ/b，β_mn＝(λ_m²+δ_n²)^0.5，其中m,n＝1,2,3…。

式(1)中A_i与B_i(i＝0,1,2,3)为傅立叶系数，其中A_i可分别由式(2)～式(5)计算得到。

$A_0=\frac{Q}{abk(\frac{h_t}{{\mathrm{th}}_t+k}+\frac{h_0}{k})}---\left(2\right)$

$A_1=\frac{2Q}{{\mathrm{abck\lambda }}_m[\phi \left({\lambda }_m\right){\lambda }_m+\frac{h_0}{k}]}[sin\left({\lambda }_m\frac{2x_c+c}{2}\right)-sin\left({\lambda }_m\frac{2x_c-c}{2}\right)]---\left(3\right)$

$A_2=\frac{2Q}{{\mathrm{abdk\delta }}_n[\phi \left({\delta }_n\right){\delta }_n+\frac{h_0}{k}]}[sin\left({\delta }_n\frac{2y_c+d}{2}\right)-sin\left({\delta }_n\frac{2y_c-d}{2}\right)]---\left(4\right)$

$\begin{array}{c}{A}_3=\frac{4Q}{{\mathrm{abcdk\lambda }}_m{\delta }_n[\phi \left({\beta }_{mn}\right){\beta }_{mn}+\frac{h_0}{k}]}\\ ·[sin\left({\lambda }_m\frac{2x_c+c}{2}\right)-\sin \left({\lambda }_m\frac{2x_c-c}{2}\right)][\sin \left({\delta }_n\frac{2y_c+d}{2}\right)-sin\left({\delta }_n\frac{2y_c-d}{2}\right)]\end{array}---\left(5\right)$

B_i可由其与A_i(i＝0,1,2,3)的对应关系式(6)～式(8)计算获得。

$B_0=-\left(\frac{h_t}{{\mathrm{th}}_t+k}\right)A_0---\left(6\right)$

B_i＝-φ(ζ)A_i,i＝1,2,3(7)

$\phi \left(\zeta \right)=\frac{\zeta \sinh \left(\zeta t\right)+\frac{h_t}{k}\cosh \left(\zeta t\right)}{\zeta \cosh \left(\zeta t\right)+\frac{h_t}{k}\sinh \left(\zeta t\right)}---\left(8\right)$

当i取1、2或3时，ζ分别替换成λ_m、δ_n或β_mn。

(4)统计出热源区域所有节点的过余温度初始结果的平均值

将步骤(3)得到的热源区域所有节点过余温度的初始结果进行求平均，记为θ_mean,0， 作为在后续迭代计算中判断收敛与否的初始解。

(5)建立热源区域的热量损失矩阵

由于步骤(3)中假设了热源区域向环境的额外对流散热，将导致计算节点的过余温度 初始值比实际值偏低。此外，扩散热阻会使热源区域的过余温度呈现中心高而边缘低的 非均匀分布，进而使整个热源上的热量损失同样为非均匀分布。因此在本步骤中，需要 通过式(9)建立热源区域参与对流散热的热量损失矩阵Q'。

式(9)中，任意一个元素Q'_IJ均表示为每四个相邻节点所包围的一个网格内的平均热量 损失，其计算式为：

$\begin{array}{c}{Q^{\prime }}_{LJ}=h_0{A}_{cource-mesh}\frac{1}{4}[\theta (x_0+(I-1)\Delta x,y_0+(J-1)\Delta y,0)+\theta (x_0+(I-1)\Delta x,y_0+J\Delta y,0)\\ +\theta (x_0+I\Delta x,y_0+(J-1)\Delta y,0)+\theta (x_0+I\Delta x,y_0+J\Delta y,0)]\end{array}---\left(10\right)$

式(10)中，x₀和y₀为热源区域最靠近坐标原点的计算节点在x-y平面上的坐标值，分 别为(x_c–c/2)和(y_c–d/2)；I＝1,…,nx；J＝1,…,ny。每一个节点的过余温度θ(x,y,z) 的取值分为两种情况：若刚刚完成热源区域过余温度的初始值计算，则取初始值；若已 实施后续的迭代计算，则取最近一次迭代计算获得的迭代值。

(6)建立热源区域的非均匀热功率矩阵

为了弥补热源区域被假设为参与对流散热的热量损失，将步骤(5)中计算得到的热量 损失矩阵叠加在原本均匀发热的热源区域，从而得到热源区域的非均匀热功率矩阵Q”， 见式(11)所示。

(7)计算出热源区域中所有节点的过余温度的迭代值

在步骤(6)的基础上，需要进行多次迭代运算才能最终准确地得到热源区域中所有节 点的过余温度。在每次迭代计算时，将热源区域的节点坐标，以及利用步骤(6)所确定的 热源区域非均匀热功率矩阵Q”计算获得的傅立叶系数C_i、D_i分别替代傅立叶系数A_i与B_i， 全部代入步骤(3)中的过余温度初始值计算公式(1)，计算出热源区域中所有节点的过余温 度θ(x,y,z)的迭代值。其中，傅立叶系数C_i与D_i(i＝0,1,2,3)分别由式(12)～式(15)和式(16)～ 式(18)计算得到。

$C_0=\frac{{Q^{\prime \prime }}_{11}+...+{Q^{\prime \prime }}_{nxny}}{abk(\frac{h_t}{{\mathrm{th}}_t+k}+\frac{h_0}{k})}---\left(12\right)$

$\begin{array}{c}{C}_1=\frac{2}{{\mathrm{abk\lambda }}_m\Delta x[\phi \left({\lambda }_m\right){\lambda }_m+\frac{h_0}{k}]}\\ ·\left\{\right[\sin \left({\lambda }_m\frac{2x_c-c+2\Delta x}{2}\right)-\sin \left({\lambda }_m\frac{2x_c-c}{2}\right)]({Q^{\prime \prime }}_{11}+\mathrm{...}+{Q^{\prime \prime }}_{1ny})+\mathrm{...}\\ +[\sin \left({\lambda }_m\frac{2x_c-c+2nx·\Delta x}{2}\right)-\sin \left({\lambda }_m\frac{2x_c-c+2(nx-1)·\Delta x}{2}\right)]({Q^{\prime \prime }}_{nx1}+\mathrm{...}+{Q^{\prime \prime }}_{nxny})\}\end{array}---\left(13\right)$

$\begin{array}{c}{C}_2=\frac{2}{{\mathrm{abk\delta }}_n\Delta y[\phi \left({\delta }_n\right){\delta }_n+\frac{h_0}{k}]}\\ ·\left\{\right[\sin \left({\delta }_n\frac{2y_c-d+2\Delta y}{2}\right)-\sin \left({\delta }_n\frac{2y_c-d}{2}\right)]({Q^{\prime \prime }}_{11}+\mathrm{...}+{Q^{\prime \prime }}_{nx1})+\mathrm{...}\\ +[\sin \left({\delta }_n\frac{2y_c-d+2ny·\Delta y}{2}\right)-\sin \left({\delta }_n\frac{2y_c-d+2(ny-1)·\Delta y}{2}\right)]({Q^{\prime \prime }}_{1ny}+\mathrm{...}+{Q^{\prime \prime }}_{nxny})\}\end{array}---\left(14\right)$

$\begin{array}{c}{C}_3=\frac{4}{{\mathrm{abk\lambda }}_m{\delta }_n\Delta x\Delta y[\phi \left({\beta }_{mn}\right){\beta }_{mn}+\frac{h_0}{k}]}\\ ·\left\{{Q^{\prime \prime }}_{11}\right[\sin \left({\lambda }_m\frac{2x_c-c+2\Delta x}{2}\right)-\sin \left({\lambda }_m\frac{2x_c-c}{2}\right)]\\ ·\sin [\left({\delta }_n\frac{2y_c-d+2\Delta y}{2}\right)-\sin \left({\delta }_n\frac{2y_c-d}{2}\right)+\mathrm{...}\\ +{Q^{\prime \prime }}_{nxny}[\sin \left({\lambda }_m\frac{2x_c-c+2nx·\Delta x}{2}\right)-\sin \left({\lambda }_m\frac{2x_c-c+2(nx-1)·\Delta x}{2}\right)]\\ ·[\sin \left({\delta }_n\frac{2y_c-d+2ny·\Delta y}{2}\right)-\sin \left({\delta }_n\frac{2y_c-d+2(ny-1)·\Delta y}{2}\right)]\}\end{array}---\left(15\right)$

$D_0=-\left(\frac{h_t}{{\mathrm{th}}_t+k}\right)C_0---\left(16\right)$

D_i＝-φ(ζ)C_i,i＝1,2,3(17)

$\phi \left(\zeta \right)=\frac{\zeta \sinh \left(\zeta t\right)+\frac{h_t}{k}\cosh \left(\zeta t\right)}{\zeta \cosh \left(\zeta t\right)+\frac{h_t}{k}\sinh \left(\zeta t\right)}---\left(18\right)$

当i取1、2或3时，ζ分别替换成λ_m、δ_n或β_mn。

(8)统计出热源区域所有节点的过余温度迭代结果的平均值

将步骤(7)得到的热源区域所有节点过余温度的迭代结果进行求平均，记为θ_mean,L。 其中迭代次数L＝1,2,3…。

(9)判断迭代结果是否收敛

每一次迭代计算完成后，均需要判断迭代结果是否收敛，其判据式为：

|θ_mean,L-θ_mean,L-1|<0.001,L＝1,2,3,…(19)

若满足判据式(19)的要求，则迭代结果收敛，进入步骤(10)；否则，迭代结果不收敛， 返回步骤(5)继续增加一次迭代计算。

(10)统计出热源区域内过余温度的局部最高值

针对本散热问题，过余温度的局部最高值必定出现在热源区域内。因此，可在最后 一次迭代计算热源区域过余温度的结果中统计出最高值，记为θ_max(x_s,y_s,0)。

(11)对矩形散热板的上表面(不含热源区域)和下表面划分网格并计算所有节点的过 余温度

按步骤(2)的网格划分方法，对矩形散热板的上表面(不含热源区域)和下表面分别划 分均匀网格，并获得对应的计算节点坐标。将该计算节点坐标，以及利用最后一次迭代 所确定的热源区域非均匀热功率矩阵Q”计算获得的傅立叶系数C_i、D_i分别替代傅立叶系 数A_i与B_i，全部代入步骤(3)中的过余温度初始值计算公式(1)，计算得到不含热源区域的 上表面、下表面上所有节点的过余温度。其中，傅立叶系数C_i与D_i(i＝0,1,2,3)分别由式(12)～ 式(15)和式(16)～式(18)计算得到。

(12)统计出矩形散热板的上表面(不含热源区域)和下表面所有节点的过余温度的平 均值

将步骤(11)中获得的矩形散热板的上表面(不含热源区域)所有节点的过余温度进行 求平均，记为θ_mean(x₀,y₀,0)；同理，统计出矩形散热板的下表面所有节点的过余温度的 平均值，记为θ_mean(x_t,y_t,t)。

(13)计算矩形散热板的上表面(不含热源区域)和下表面的对流散热总功率

在步骤(11)中获得的矩形散热板的上表面(不含热源区域)所有节点的过余温度结果 中，将每四个相邻节点的过余温度记为θ₁、θ₂、θ₃和θ₄，并按照式(20)计算出每四个相邻 节点所包围的一个网格内的对流散热功率Q_up-mesh。

$Q_{up-mesh}=h_0{A}_{up-mesh}\frac{1}{4}({\theta }_1+{\theta }_2+{\theta }_3+{\theta }_4)---\left(20\right)$

式(20)中A_up-mesh为步骤(11)中划分的不含热源区域的上表面每个均匀网格的面积。将矩形 散热板的上表面(不含热源区域)所有网格的对流散热功率求和，得到其对流散热总功率 Q_up。

同理，可计算出矩形散热板下表面的对流散热总功率Q_low。计算矩形散热板的下表 面上每四个相邻节点所包围的一个网格内的对应散热功率Q_low-mesh的公式见式(21)所示。

$Q_{low-mesh}=h_t{A}_{low-mesh}\frac{1}{4}({\theta }_{11}+{\theta }_{22}+{\theta }_{33}+{\theta }_{44})---\left(21\right)$

式(21)中，A_low-mesh为步骤(11)中划分的下表面每个均匀网格的面积，θ₁₁、θ₂₂、θ₃₃和θ₄₄为下表面每四个相邻节点的过余温度。

(14)计算最大扩散热阻和最大总热阻

通过式(22)和式(23)可分别计算出热源向矩形散热板的上、下表面散热的最大扩散热 阻R_{s_0,max}和R_{s_t,max}，并且通过式(24)可计算出热源到环境间的最大总热阻R_th,max。

$R_{s\_0,max}=\frac{{\theta }_{max}(x_s,y_s,0)-{\theta }_{mean}(x_0,y_0,0)}{Q_{up}}---\left(22\right)$

$R_{s\_t,max}=\frac{{\theta }_{max}(x_s,y_s,0)-{\theta }_{mean}(x_t,y_t,t)}{Q_{low}}---\left(23\right)$

$R_{th,max}=\frac{{\theta }_{\max }(x_s,y_s,0)}{Q}---\left(24\right)$

需要说明的是，在执行步骤(3)、步骤(7)和步骤(11)中，均须对式(1)进行计算。为保 证计算结果收敛且不耗费过多的计算资源，可对无限累加变量m和n分别预先设定一个 最小的累加次数。当累加次数大于设定的最小值后，通过判断此次累加结果与上一次结 果的差值是否小于0.001，若满足则完成本次累加计算。

需要说明的是，在执行步骤(2)和步骤(11)中，均须对所计算的区域划分一定数量的 均匀网格。理论上，划分的网格数量越多，则计算结果越趋近于真实值，但计算时长却 急剧增加。因此可通过描绘热源区域过余温度的最高值θ_max(x_s,y_s,0)随网格数量而增加的 变化数据图，找出数据变化趋于饱和的起点，将该起点对应的网格数量定为网格划分的 最少数量。

需要说明的是，所述含偏心热源的矩形散热板双面对流散热性能的解析计算方法， 可在具有计算功能的软件中对所有步骤进行编程计算。

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