引言随着计算机技术、计算数学以及应用数学的不断发展，特别是以并行计算机和并行计算为基础的高性能计算的兴起，使得计算能力大幅度提高，从而能够精确求解各种复杂的微分方程问题，数值模拟逐渐成为学者认识自然规律的主要手段，物理实验逐渐变成辅助手段。在21世纪高性能计算已成为高科技的核心技术，基于高性能计算的数值模拟成为经济竞争力的关键。高性能计算应用水平已成为衡量各行业高科技水平的标准。传热学研究广泛应用了数值模拟方法。传热学是研究由温差引起的热量传递规律及其应用的工程技术学科[1]。基于计算机VB语言，对传热学的二维稳态导热过程进行数值模拟计算，编写程序，直观地展示传热过程中的温度分布，分析固定边界热流密度变化和固定边界变化对导热的影响。1传热学数值求解的开发应用实例准备二维受热平板，导热系数k=1 000 W/(m·K)，板厚1 cm，二维受热平板模型如图1所示。西侧边界有稳定热流流入，热流强度q=500 kW/m2，东侧和南侧边界绝热，北侧边界保持常温，TN=100 ℃[2]。10.3969/j.issn.1004-7948.2024.02.014.F001图1二维受热平板模型控制方程为∂∂xk∂T∂x+∂∂yk∂T∂y=0，方程离散的方式采用中心差分[3]的形式。（1）控制体网格划分：采用内节点法，Δx=Δy=0.1。（2）方程的离散[4-5]（适合节点6和节点7）：∫sn∫we∂∂xk∂T∂xdxdy+∫sn∫we∂∂yk∂T∂ydxdy=k∂T∂xe-k∂T∂xw·Δy+k∂T∂yn-k∂T∂ys·Δx=kTE-TPΔx-kTP-TWΔx·Δy+kTN-TPΔy-kTP-TSΔy·Δx=0 (1)整理得：aPTP=aWTW+aETE+aSTS+aNTN+b,  b=0 (2)（3）边界节点1离散：∫sn∫we∂∂xk∂T∂xdxdy+∫sn∫we∂∂yk∂T∂ydxdy=k∂T∂xe-k∂T∂xw·Δy+k∂T∂yn-k∂T∂ys·Δx=kT5-T1Δx+kq·Δy+kT2-T1Δy+0·Δx=0 (3)整理得：aPTP=aETE+aNTN+b,  b=q·k·Δy (4)2 000T1=1 000T5+1 000T2+50 000 (5)（4）采用同样的方法离散节点2、节点3、节点4、节点5、节点8、节点9、节点10、节点11和节点12。节点2和节点3：aPTP=aETE+aSTS+aNTN+b,  b=q·k·Δy (6)节点4：aPTP=aETE+aSTS+b,  b=q·k·Δy+2kΔxΔyTN (7)节点5：aPTP=aWTW+aETE+aNTN+b,  b=0 (8)节点8：aPTP=aWTW+aETE+aSTS+b,  b=2kΔxΔyTN (9)节点9：aPTP=aWTW+aNTN+b,  b=0 (10)节点10和节点11：aPTP=aWTW+aSTS+aNTN+b,  b=0 (11)节点12：aPTP=aWTW+aSTS+b,b=2kΔxΔyTN (12)其中：aW=ΔyΔxkw (13)aE=ΔyΔxke (14)aS=ΔxΔyks (15)aN=ΔxΔykn (16)aP=aW+aE+aS+aN-SP (17)固定热流时SP=0，固定温度边界时SP=-2kΔxΔy。（5）整理得到矩阵形式为五对角矩阵[6]（其中M1为横向网格数，L1为纵向网格数）。离散矩阵形式如图2所示。10.3969/j.issn.1004-7948.2024.02.014.F002图2离散矩阵形式（6）使用Visual Basic(VB)语言建立控制方程式，确定初始条件及边界条件，划分计算网格，从而生成计算节点；随后建立离散方程，确定离散初始条件和边界条件，通过给定求解需要控制的参数进行离散方程的求解；依据方程是否收敛进行迭代计算，直至计算结果收敛。计算流程如图3所示。10.3969/j.issn.1004-7948.2024.02.014.F003图3计算流程（7）VB程序的运行界面如图4所示。10.3969/j.issn.1004-7948.2024.02.014.F004图4VB程序的运行界面2固定边界参数变化对受热平板内部温度的影响参考文献[7]，不同边界温度下受热平板内部各节点的温度分布如表1所示。当固定边界的温度由100 ℃变为200 ℃时，受热平板内部各节点的温度在100 ℃的边界温度条件下增加100 ℃；固定边界的温度变为其他数据时，各节点的温度在原来结果上增加一个固定边界温度的变化值。10.3969/j.issn.1004-7948.2024.02.014.T001表1不同边界温度下受热平板内部各节点的温度分布项目t1t2t3t4t5t6t7t8t9t10t11t12tn1=100 ℃260242206146228211178130212197166124tn2=2tn1=200 ℃360342306246328311278230312297266224tn3=tn1+50310292256196278261228180262247216174tn2-tn1100100100100100100100100100100100100tn3-tn1505050505050505050505050℃不同边界热流密度下受热平板内部各节点的温度分布如表2所示。当固定边界热流密度由500 kW/m2变为1 000 kW/m2时，受热平板内部各节点的温度变化无数值型的规律；当固定边界热流密度变成其他数据时，也未得出热流密度的变化和受热平板内部温度分布的规律。10.3969/j.issn.1004-7948.2024.02.014.T002表2不同边界热流密度下受热平板内部各节点的温度分布项目t1t2t3t4t5t6t7t8t9t10t11t12q1=500 kW/m2260242206146228211178130212197166124q2=2q1=1 000 kW/m2420385311193356322256159324293232148q3=q1+50276257216151241222186133223206173126q2-q1160143105471281117829112966624q3-q1161510513118311972℃比较固定边界温度变化与固定边界热流密度变化对受热平板内部温度分布的影响[8]。固定边界温度和固定边界热流密度均变为原来的两倍时，不同边界条件变化率相同时受热平板内部各节点的温度分布如表3所示。10.3969/j.issn.1004-7948.2024.02.014.T003表3不同边界条件变化率相同时受热平板内部各节点的温度分布项目t1t2t3t4t5t6t7t8t9t10t11t12q1260242206146228211178130212197166124q2=2q1420385311193356322256159324293232148q2-q1160143105471281117829112966624tn1260242206146228211178130212197166124tn2=2tn1360342306246328311278230312297266224tn2-tn1100100100100100100100100100100100100℃虽然固定边界温度和固定边界热流密度均变为原来的两倍，但是热流密度对受热平板内部温度分布的影响比边界温度大一些。因为边界温度的变化对内部节点的影响相同，但是热流密度的变化根据其距离边界的远近决定对内部节点的影响程度。为了证明是否仅在上述示例情况下才会出现边界温度变化值也是温度分布变化值的情况，重新编写了两侧都是温度边界的程序。固定一侧边界不变而改变另一侧时，两次各节点的温度分布之间无准确的数值关系；两者同时变化时，也未发现这样的规律。这说明规律并不在所有的情况下都成立，在特定的边界条件下才会发生。经过计算数据的分析，在固定热流密度不变的情况下，改变边界温度的值即为最终温度分布的增加值。根据建立的离散方程以及计算过程可以发现，边界温度的变化主要影响常数项的值，对其他各项的值均无影响。说明出现这种规律的主要影响因素是常数项。拓展的数值计算发现，规律必须在特定的情况下才能出现。3导热相对于不同边界条件的灵敏度灵敏度Sxy[9]为：Sxy=ΔyyΔxx=xy·ΔyΔx (18)式中：Sxy——y对x的灵敏度。3.1导热对固定边界热流密度的灵敏度计算以内节点6为例，当q=500 kW/m2时，t6=211 ℃；当q变为原来的两倍时（即Δq=500 kW/m2），此时t6'=322 ℃（即Δt6=111 ℃）。经计算[10]，Sqt6=0.53。各节点对固定边界热流密度的灵敏度如表4所示。10.3969/j.issn.1004-7948.2024.02.014.T004表4各节点对固定边界热流密度的灵敏度项目t1t2t3t4t5t6t7t8t9t10t11t12q1260242206146228211178130212197166124q2=2q1420385311193356322256159324293232148Δti160143105471281117829112966624Δq500500500500500500500500500500500500Sqti0.620.590.510.320.560.530.440.220.530.490.400.19℃3.2导热对固定边界温度的灵敏度计算当tn=100 ℃时，t6=211 ℃；当tn变为原来的两倍时（即Δtn=100 ℃），此时t6'=311 ℃（即Δt6=100 ℃）。经计算，Stnt6=0.47。各节点对固定边界温度的灵敏度如表5所示。10.3969/j.issn.1004-7948.2024.02.014.T005表5各节点对固定边界温度的灵敏度项目t1t2t3t4t5t6t7t8t9t10t11t12tn1260242206146228211178130212197166124tn2=2tn1360342306246328311278230312297266224Δti100100100100100100100100100100100100Δtn100100100100100100100100100100100100Stnti0.380.410.490.680.440.470.560.770.470.510.600.81℃3.3导热对不同边界变化的灵敏度比较将各节点的灵敏度进行比较，各节点对不同边界变化的灵敏度如表6所示。10.3969/j.issn.1004-7948.2024.02.014.T006表6各节点对不同边界变化的灵敏度项目t1t2t3t4t5t6t7t8t9t10t11t12Sqti0.620.590.510.320.560.530.440.220.530.490.400.19Stnti0.380.410.490.680.440.470.560.770.470.510.600.81℃通过数据分析可以看出，边界条件相同的内部节点，导热对固定边界热流密度的灵敏度高于对固定边界温度的灵敏度；但边界节点的灵敏度与其距边界条件的距离有关，边界节点对距离较近的边界条件的灵敏度更高。4结语边界条件相同的内部节点，导热对固定边界热流密度的灵敏度高于对固定边界温度的灵敏度；但边界节点的灵敏度与其距边界条件的距离有关，边界节点对距离较近的边界条件的灵敏度更高。文中研究暂未考虑其他参数的影响，后续研究可以在此基础上综合除固定边界热流密度与固定边界温度外的次要参数进行数据分析，进一步研究导热对固定边界热流密度的灵敏度与导热对固定边界温度的灵敏度的高低。