数值传热学第一次大作业

目录1题目..............................................................22问题初始化........................................................23FTCS格式.........................................................33.1离散差分方程....................................................33.2计算程序........................................................33.3计算结果........................................................54BTCS格式.........................................................94.1离散差分方程................................................94.2计算程序...................................................104.3计算结果...................................................115BTCS格式........................................................125.1离散差分方程...............................................125.2计算程序...................................................135.3计算结果...................................................15非定常热传导问题1题目求方程:𝜕𝑢𝜕𝑡=𝜈𝜕2𝑢𝜕𝑥2初边值问题的有限差分解。求解域为𝑥,𝑡∈0,1×[0,∞]，初始条件和边界条件为𝑢𝑥,0=𝑓(𝑥)𝑢0,𝑡=𝑎(𝑡)𝑢1,𝑡=𝑏(𝑡)其中𝜈=1，初值条件的具体取法为:𝑓x=01−10/3𝑥+(10/3)0𝑥0.30.3≤𝑥≤0.70.7𝑥1.0取网格点数M=100，要求计算t=0.01，t=0.1，t=1，t=10时的数值解，计算的时间步长取法为：σ=𝜈Δ𝑡Δ𝑥2=0.1,0.5,1.0差分格式：FTCS格式，BTCS格式，CNCS格式。2问题初始化问题为一维，要求网格点数M=100，故x方向按空间步长：Δ𝑥=1M=0.01而根据σ规定时间步长Δ𝑡取为1×10−5，5×10−5，1×10−4；且要求计算t=0.01，t=0.1，t=1，t=10时的数值解，故时间网格数𝑁=𝑡𝛥𝑡，表1为不同时间t与时间步长Δ𝑡所对应的网格数𝑁分布。𝜎0.10.51.0𝑡Δ𝑡1×10−55×10−51×10−40.011×1032001000.11×1042×1031×10311×1052×1041×104101×1062×1051×1053FTCS格式3.1离散差分方程内部节点：𝑈𝑘𝑛+1−𝑈𝑘𝑛Δ𝑡=𝜈𝑈𝑘+1𝑛−2𝑈𝑘𝑛+𝑈𝑘−1𝑛Δ𝑥2，𝑘=1,2,3,…𝑀−1;𝑛=0,1,2,…,𝑁−1定解条件：初始值：𝑢𝑥,0=𝑓𝑥，即：𝑈𝑘0=𝑓𝑥=𝑓𝑘∙Δ𝑥，𝑘=0,1,2,…𝑀由初始条件的取值函数𝑓𝑥，将边界条件设为定值0，即：𝑈0𝑛=𝑈𝑀𝑛=0，𝑛=0,1,…,𝑁3.2计算程序#includeiostream//输入输出流头文件#includecmath//包含了许多数学函数的库文件usingnamespacestd;//使用std这个名字空间intmain()//返回整形变量，若voidmain不用返回{doublec=0.1;//参数设定doubledeltaX=0.01;//空间步长表1不同的计算时间t与时间步长Δ𝑡所对应的时间网格数𝑁分布doublet=0.01;//选定时刻doubledeltaT=c*deltaX*deltaX;//时间步长intconstM=100;//空间网格数intN=t/deltaT;//10s内的时间网格数doubleT[M+1];//存放初始值和边界值doubleresult[M+1];//存放结果inti,j;//初始化for(i=1;iM;i++){if(i*deltaX0.3){T[i]=0;}elseif(i*deltaX0.7){T[i]=1;}else{T[i]=10/3*(1-i*deltaX);}T[0]=0;T[M]=0;}//节点差分for(i=1;i=N+1;i++)//时间迭代{T[0]=0;T[M]=0;for(j=0;j=M;j++){result[j]=T[j];}for(j=1;j100;j++)//空间迭代{T[j]=result[j]+deltaT*(result[j+1]-2*result[j]+result[j-1])/deltaX/deltaX;}if(i*deltaT==t){coutFTCSendl;coutTime=tendl;coutσ=cendlendl;for(j=0;j=M;j++){coutresult[j]endl;//输出结果}coutendlendl;}}}3.3计算结果本非定常热传导问题设置不同取值时刻与不同时间步长取法来探究𝑡与𝜎对数值解的影响。图1、2、3、4显示不同时刻t=0.01，t=0.1，t=1，t=10时的温度分布，取𝜎为定值𝜎=0.1。图1图2图3图1、2、3、4分别是一维非稳态导热随时间变化的过程。随着时间从初始时刻到第10s时刻，温度分布是逐渐趋向于零的，且越远离边界，温度越高，这些符合实际分度分布情况。下图5是相同计算时间t=0.01时，时间步长取法σ=0.1、0.5结果比较，由图知，两个计算结果基本相同。说明模拟计算具有一定的稳定性。下图6是计算时间t=0.01时，时间步长取法σ=1的计算结果。由图可知方程的解出现了振荡，是因为差分方程能得到稳定解的条件是方程中各项系数大于零。当σ=1时，𝑈𝑘𝑛的系数1−2σ=−10，所以会出现振荡。图4图5图64BTCS格式4.1离散差分方程内部节点：𝑈𝑘𝑛−𝑈𝑘𝑛−1Δ𝑡=𝜈𝑈𝑘+1𝑛−2𝑈𝑘𝑛+𝑈𝑘−1𝑛Δ𝑥2，𝑈𝑘𝑛+1−𝑈𝑘𝑛Δ𝑡=𝜈𝑈𝑘+1𝑛+1−2𝑈𝑘𝑛+1+𝑈𝑘−1𝑛+1Δ𝑥2，𝑘=1,2,3,…𝑀−1;𝑛=0,1,2,…,𝑁−11+2𝜎𝑈𝑘𝑛+1=𝜎𝑈𝑘+1𝑛+1+𝑈𝑘−1𝑛+1+𝑈𝑘𝑛；其中σ=𝜈Δ𝑡Δ𝑥2。定解条件：初始值：𝑢𝑥,0=𝑓𝑥，即：𝑈𝑘0=𝑓𝑥=𝑓𝑘∙Δ𝑥，𝑘=0,1,2,…𝑀边界条件同FTCS格式设为定值0。差分方程的求解：构造n+1时刻的线性方程组：𝐴𝑖𝑈𝑖=𝐵𝑖𝑈𝑖+1+𝐶𝑖𝑈𝑖−1+𝐷𝑖，𝑖=1,2,…,𝑀−1;端点条件：i=1，𝐶𝑖=0i=𝑀−1，𝐵𝑖=0其中，各项系数：𝐴𝑖=1+2𝜎，𝑖=1,2,…,𝑀−1𝐵𝑖=𝜎，𝑖=1,2,…,𝑀−2𝐶𝑖=𝜎，𝑖=1,2,…,𝑀−1𝐷𝑖=𝑈𝑖𝑛，𝑖=1,2,…,𝑀−2再由THOMAS算法将三对角阵消为容易解的二对角阵。4.2计算程序#includeiostream//输入输出流头文件#includecmath//包含了许多数学函数的库文件usingnamespacestd;//使用std这个名字空间intmain(){doublec=0.1;//参数设定doubledeltaX=0.01;//空间步长doublet=0.01;//选定时刻doubledeltaT=c*deltaX*deltaX;//时间步长intconstM=100;//空间网格数intN=t/deltaT;//时间网格数doubleT[M+1],P[M],Q[M],A[M],B[M],C[M],D[M];intj,i;//初始化for(i=1;iM;i++){if(i*deltaX0.3){T[i]=0;}elseif(i*deltaX0.7){T[i]=1;}else{T[i]=10/3*(1-i*deltaX);}T[0]=0;T[M]=0;}//节点差分for(j=0;j=N+1;j++){T[0]=0;//赋边界值T[M]=0;//计算内点值for(i=1;i=M-1;i++){A[i]=1+2*c;}for(i=1;i=M-2;i++){B[i]=c;}B[M-1]=0;for(i=2;i=M-1;i++)4.3计算结果下图显示不同时刻t=0.01，t=0.1，t=1，t=10时的温度分布，取𝜎为定值𝜎=0.1。下图是计算时间t=0.01时，𝜎取0.1、0.5、1.0时的结果，可以看出由于时间的差分采用隐式格式，具有较好的稳定性，没有出现振荡，只是有部分数值偏离。5BTCS格式5.1离散差分方程显示格式：𝑈𝑘𝑛+1−𝑈𝑘𝑛Δ𝑡=𝜈𝑈𝑘+1𝑛−2𝑈𝑘𝑛+𝑈𝑘−1𝑛Δ𝑥2，𝑘=1,2,3,…𝑀−1;𝑛=0,1,2,…,𝑁−1隐式格式：𝑈𝑘𝑛+1−𝑈𝑘𝑛Δ𝑡=𝜈𝑈𝑘+1𝑛+1−2𝑈𝑘𝑛+1+𝑈𝑘−1𝑛+1Δ𝑥2，𝑘=1,2,3,…𝑀−1;𝑛=0,1,2,…,𝑁−11+2𝜎𝑈𝑘𝑛+1=𝜎𝑈𝑘+1𝑛+1+𝑈𝑘−1𝑛+1+𝑈𝑘𝑛；其中σ=𝜈Δ𝑡Δ𝑥2。显示、隐式相加，得到：2+2𝜎𝑈𝑘𝑛+1=𝑈𝑘𝑛2𝜎−2+𝑈𝑘+1𝑛+𝑈𝑘+1𝑛+1+𝑈𝑘−1𝑛+1+𝑈𝑘−1𝑛定解条件：初始值：𝑢𝑥,0=𝑓𝑥，即：𝑈𝑘0=𝑓𝑥=𝑓𝑘∙Δ𝑥，𝑘=0,1,2,…𝑀边界条件同FTCS格式设为定值0。差分方程的求解：构造n+1时刻的线性方程组：𝐴𝑖𝑈𝑖=𝐵𝑖𝑈𝑖+1+𝐶𝑖𝑈𝑖−1+𝐷𝑖，𝑖=1,2,…,𝑀−1;端点条件：i=1，𝐶𝑖=0，D𝑖=𝑈2𝑛+𝑈0𝑛+𝑈1𝑛2𝜎−2i=𝑀−1，𝐵𝑖=0，DM−1=𝑈M𝑛+𝑈𝑀−2𝑛+𝑈𝑀−1𝑛2𝜎−2其中，各项系数：𝐴𝑖=2+2𝜎，𝑖=1,2,…,𝑀−1𝐵𝑖=1，𝑖=1,2,…,𝑀−2𝐶𝑖=1，𝑖=1,2,…,𝑀−1𝐷𝑖=𝑈i+1𝑛+𝑈𝑖−1𝑛+𝑈𝑖𝑛2𝜎−2，𝑖=1,2,…,𝑀−2同样再由THOMAS算法将三对角阵消为容易解的二对角阵。5.2计算程序#includeiostream//输入输出流头文件#includecmath//包含了许多数学函数的库文件usingnamespacestd;//使用std这个名字空间intmain(){doublec=0.1;//参数设定doubledeltaX=0.01;//空间步长doublet=0.01;//选定时刻doubledeltaT=c*deltaX*deltaX;//时间步长constintM=100;//空间网格数intN=t/deltaT;//时间网格数doubleT[M+1],P[M],Q[M],A[M],B[M],C[M],D[M]