二维稳态导热问题-数值解法

二维稳态导热问题的数值解法1二维稳态导热问题的数值解法1.问题重述1.1．第一题如图所示，一个无限长矩形柱体，其横截面的边长分别为L1和L2，常物性。该问题可视为二维稳态导热问题，边界条件如图中所示，其中L1=0.6m，L2=0.4m，Tw1=60℃，Tw2=20℃，λ=200W/(m•K)。1)编写程序求解二维导热方程。2)绘制x=L1/2和y=L2/2处的温度场，并与解析解进行比较。已知矩形内的温度场的解析解为：1211w2w1shshsin,LLLyLxttyxt。1.2．第二题将第一题中2yL处的边界条件变为2wtt，其他条件不变。1)编写程序求解二维导热方程并计算从y=0处导入的热量2。2)当时，该二维导热问题可简化为一维导热问题。在一维的近似下，计算从y=0处导入的热量1，并比较不同L2/L1下21的比值。由该问题的解析解可知：L2/L10.0070.010.050.080.1210.99870.99120.9560.930.9122.第一题二维稳态导热问题的数值解法22.1第一问2.1.1问题分析对题干中所给的区域划分网格，如图2-1所示，其中，有（）（）个未知温度节点，由题意可知，该区域内无内热源，所有未知节点皆为内节点，将所有节点分为三类，如下表所示。对于虚线框内的节点定义为第一类节点，为内节点，由热平衡法可得：其中：，当即时，节点类别下标变化范围离散方程①②a②b②c③2.1.2计算二维稳态导热问题的数值解法3因为，取M=7，N=5，步长，共15个内节点，经过32次迭代后，所有数值6位有效数字已经稳定，计算结果如下表所示，计算程序见附录1n=560.00000069.99999277.32049980.00000077.32052670.00003860.000000n=460.00000065.81187670.06647271.62376970.06648965.81189960.000000n=360.00000063.18105165.50975166.36212165.50976663.18107060.000000n=260.00000061.40260262.42938562.80521962.42939661.40261560.000000n=160.00000060.00000060.00000060.00000060.00000060.00000060.000000m=1m=2m=3m=4m=5m=6m=7将上述数据用MATLAB画图，拟合三维曲面如图2-2所示，MATLAB作图代码见附录22.1.3进一步精确计算因为网格划分的步长越小离散方程所得的结果越精确，所以，重新取M=31，N=21，步长，内节点为551个。经过1000次迭代后，所有数值6位有效数字已经稳定。将所得数据用MATLAB进行三维曲面拟合，所得图形如图2-3所示。2.2第二问图2-2图2-3二维稳态导热问题的数值解法42.2.1时，m=4，，温度场分布如图2-4所示，将二维稳态传热方程所得的解与理论解进行比较，如下表所示y00.10.20.30.4计算值60.00000062.80521966.36212171.62376980.000000理论值60.00000062.74016466.24888071.51026080.000000相对误差00.1037%0.1709%0.1587%02.2.2时，n=3，，温度场分布如图2-5所示，将二维稳态传热方程所得的解与理论解进行比较，如下表所示x00.10.20.30.40.50.6计算值60.00000063.18105165.50975166.36212165.50976663.18107060.000000理论值60.00000063.12444065.41168966.24888065.41168963.12444060.000000相对误差0%0.0897%0.1499%0.1709%0.1499%0.0897%0.0000%图2-4处的温度场图2-5处的温度场二维稳态导热问题的数值解法53.第二题3.1第一问3.1.1问题分析网格依旧沿用第一题的划分方式，则节点类型总共有2种，所有未知节点均为内节点，所有边界上的节点均已知，如下表所示：节点类别下标变化范围离散方程①②a②b②c②d3.1.2计算：因为，取M=7，N=5，步长，共15个内节点，经过32次迭代后，所有数值6位有效数字已经稳定，计算结果如下表所示，计算程序见附录3n=520.00000020.00000020.00000020.00000020.00000020.00000020.000000n=460.00000041.64665235.51149934.00759035.51149941.64665260.000000n=360.00000051.07511046.39175345.00736446.39175351.07511060.000000n=260.00000056.26203753.97303753.23836053.97303756.26203760.000000n=160.00000060.00000060.00000060.00000060.00000060.00000060.000000m=1m=2m=3m=4m=5m=6m=7二维稳态导热问题的数值解法6将上述数据用MATLAB画图，拟合三维曲面如图3-1所示从y=0处导入的热量可以近似看作从y=0向y=0.1处传递的热量，y=0.1处的温度分别为x00.10.20.30.40.50.6t60.00000056.26203753.97303753.23836053.97303756.26203760.000000则整理，得由于所以()图3-1二维稳态导热问题的数值解法73.2第二问当时，对该区域重新划分网格，取每个网格都为正方形，横向分为M份，纵向分为N份，则内节点方程为：则：[∑（）]当时，,令M=1000，N=7，由此可求得，计算程序见附录4:∑63（𝐾）[∑（）]633366因为，若将该导热问题看作是一维导热问题则所以33663当时，，可令M=1000，N=10，求得当时，，可令M=1000，N=50，求得当时，，可令M=1000，N=80，求得二维稳态导热问题的数值解法8当时，，可令M=1000，N=100，求得结果如下表所示，计算程序与附录4所示的算法程序相似0.0070.010.050.080.15679.1883243964.800897764.702490464.699940364.6993511135837.6648792960.1794152940.498092939.988072939.87021142857.1429800000.0000160000.0000100000.000080000.00000.993860.991200.955880.929400.91175（解析解）0.99870.99120.9560.930.912二维稳态导热问题的数值解法9附录1第一题第一问算法程序#includestdafx.hint_tmain(intargc,_TCHAR*argv[]){return0;}#definePI3.14159#includemath.hvoidmain(){inti,m,n;doublet[7][5];for(m=0;m7;m++){t[m][0]=60;t[m][4]=60+20*sin(PI*0.1*m/0.6);}for(n=0;n5;n++){t[0][n]=60;t[6][n]=60;}for(m=1;m6;m++){for(n=1;n4;n++)二维稳态导热问题的数值解法10t[m][n]=0;}for(i=0;i90;i++){for(m=1;m6;m++){for(n=1;n4;n++)t[m][n]=0.25*(t[m+1][n]+t[m-1][n]+t[m][n+1]+t[m][n-1]);}}for(n=0;n5;n++){for(m=0;m7;m++)printf(%10f,t[m][n]);printf(\n);}}二维稳态导热问题的数值解法11附录2第一题第一问MATLAB作图代码A=[0060.00000000.160.00000000.260.00000000.360.00000000.460.00000000.560.00000000.660.0000000.1060.0000000.10.161.4026020.10.262.4293850.10.362.8052190.10.462.4293960.10.561.4026150.10.660.0000000.2060.0000000.20.163.1810510.20.265.5097510.20.366.3621210.20.465.5097660.20.563.1810700.20.660.000000二维稳态导热问题的数值解法120.3060.0000000.30.165.8118760.30.270.0664720.30.371.6237690.30.470.0664890.30.565.8118990.30.660.0000000.4060.0000000.40.169.9999920.40.277.3204990.40.380.0000000.40.477.3205260.40.570.0000380.40.660.000000];x=A(:,1);y=A(:,2);z=A(:,3);minx=min(x);maxx=max(x);miny=min(y);maxy=max(y);[X,Y,Z]=griddata(x,y,z,linspace(minx,maxx)',linspace(miny,maxy),'v4');figure,surf(X,Y,Z)二维稳态导热问题的数值解法13附录3第二题第一问算法程序#includestdafx.hint_tmain(intargc,_TCHAR*argv[]){return0;}#includemath.hvoidmain(){inti,m,n;doublet[7][5];for(m=0;m7;m++){t[m][0]=60;t[m][4]=20;}for(n=0;n4;n++){t[0][n]=60;t[6][n]=60;}for(m=1;m6;m++){for(n=1;n4;n++)t[m][n]=0;二维稳态导热问题的数值解法14}for(i=0;i300;i++){for(m=1;m6;m++){for(n=1;n4;n++)t[m][n]=0.25*(t[m+1][n]+t[m-1][n]+t[m][n+1]+t[m][n-1]);}}for(n=0;n5;n++){for(m=0;m7;m++)printf(%10f,t[m][n]);printf(\n);}}二维稳态导热问题的数值解法15附录4第二题第二问算法程序#includestdafx.hint_tmain(intargc,_TCHAR*argv[]){return0;}voidmain(){inti,m,n;doublet[1001][8];doublej=0.0;for(m=0;m1001;m++){t[m][0]=60;t[m][7]=20;}for(n=0;n7;n++){t[0][n]=60;t[1000][n]=60;}for(m=1;m1000;m++){for(n=1;n7;n++)t[m][n]=0;二维稳态导热问题的数值解法16}for(i=0;i3000;i++){for(m=1;m1000;m++){for(n=1;n7;n++)t[m][n]=0.25*(t[m+1][n]+t[m-1][n]+t[m][n+1]+t[m][n-