有限差分法解热传导问题

傅里叶定律：在导热现象中，单位时间内通过给定截面的热量，正比例于垂直于该界面方向上的温度变化率和截面面积，而热量传递的方向则与温度升高的方向相反。导热微分方程：定解条件：使微分方程获得适合某一问题的的解的附加条件。边界条件：NDSolve[{D[u[x,t],t]=D[u[x,t],{x,2}],u[x,0]=x*(1-x),u[0,t]=0,u[1,t]=0,u,{x,0,1},{t,0,0.3}]NDSolve[{D[u[x,t],t]=D[u[x,t],{x,2}]+x,u[x,0]=x*(1-x),u[0,t]=0,u[1,t]=0,u,{x,0,1},{t,0,0.3}]NDSolve[{D[u[x,t],t]=D[u[x,t],{x,2}],u[x,0]=Sin[x]*Sin[x],Derivative[1,0][u][0,t]=0,Derivative[1,0][u][Pi,t]=0,u,{x,0,Pi},{t,0,1}]NDSolve[{D[u[x,t],t]=D[u[x,t],{x,2}],u[x,0]=-x*(1-x),Derivative[1,0][u][0,t]=0,Derivative[1,0][u][1,t]=3-u[1,t],u,{x,0,1},{t,0,0.3}]建立控制方程及定解条件确定节点（区域离散化）建立节点物理量的代数方程设立温度场的迭代初值求解代数方程是否收敛解的分析改进初场是否xyxynm(m,n)MN泰勒级数展开法：若取上面式右边的前三项，并将两式相加移项整理即二阶导数的中心差分：同样可得：)(222,1,,1,22xoxtttxtnmnmnmnm)(2221,,1,,22yoytttytnmnmnmnm对于二维稳态导热问题，在直角坐标中，其导热微分方程为：其节点方程为：0ytt2txtt2tj,i,v21j,ij,i1j,i2j,1ij,ij,1i0ytxtv2222热平衡法基本思想：对每个有限大小的控制容积应用能量守恒，从而获得温度场的代数方程组，它从基本物理现象和基本定律出发，不必事先建立控制方程，依据能量守恒和Fourier导热定律即可。能量守恒：流入控制体的总热流量＋控制体内热源生成热＝流出控制体的总热流量＋控制体内能的增量内部节点：01,1,,1,1nmnmnmnmΦΦΦΦ0右左下上＋ΦΦΦΦ(m,n)oyx(m-1,n)(m+1,n)(m,n-1)xxyy(m,n+1)0vΦΦΦΦΦ右左下上＋0,1,,1,,,1,,1yxΦyttxyttxxttyxttynmnmnmnmnmnmnmnm042,1,1,,1,1ΦxtttttnmnmnmnmnmΦxtttttnmnmnmnmnm21,1,,1,1,4yx171.边界节点离散方程的建立：qwxyqw(1)平直边界上的节点2,1,1,,1,224xΦttqxttnmnmnmwnmnm0222,,1,,1,,,1yxΦyttxyttxyqxttynmnmnmnmnmwnmnmyx18(2)外部角点2222,1,,1,xΦqxtttnmwnmnmnm0222222,,1,,,1yxΦyttxqxqyxttynmnmnmwwnmnmyxxyqw19(3)内部角点)22322(6122,11,1,,1,wnmnmnmnmnmqxxttttt0432222,,1,,1,,,1,,1yxΦqxyttxyttxqyxttyxttynmwnmnmnmnmwnmnmnmnmyxxyqw写出所有内节点和边界节点的温度差分方程n个未知节点温度，n个代数方程式：nnnnnnnnnnnbtatatatbtatatatbtatatat..............................................................2211222221212112121111Gauss－Seidel迭代)(1)(1)(1)1(11)1(11)1(2)(1)(323)1(12122)1(21)(1)(21211)1(1nknnnknnnknknnkkkknnkkbxaxaaxbxaxaxaaxbxaxaax)(11)(11)1()1(inijkjijijkjijiikibxaxaaxGauss－Seidel迭代)()()1(1)1(bUxLxDxkkkbDUxDxLDIkk1)(1)1(1)(bDLDIUxDLDIxkk111)(111)1()()(bLDUxLDxkk1)(1)1()()(bLDfULDB11)(,)(Gauss－Seidel迭代200℃100℃100℃Tf=0℃k=1W/(m*K)h=10W/(m2*℃)(qw=h*(Tw-Tf))t1t2t3t18t17t16Ax=bb=[300,200,200,300,100,0,0,100,100,0,0,100,100,0,0,100,0,0]’x=[t1,t2,t3,...t18]’偏微分方程工具箱Step1“Draw模式”绘制平面有界区域，通过公式把Matlab系统提供的实体模型：矩形、圆、椭圆和多边形，组合起来，生成需要的平面区域.Step2“Boundary模式”定义边界，声明不同边界段的边界条件.Step3“PDE模式”定义偏微分方程，确定方程类型和方程系数c,a,f,d，根据具体情况，还可以在不同子区域声明不同系数.Step4“Mesh模式”网格化区域，可以控制自动生成网格的参数，对生成的网格进行多次细化，使网格分割更细更合理.Step5“Solve模式”解偏微分方程，对于椭圆型方程可以激活并控制非线性自适应解题器来处理非线性方程；对于抛物线型方程和双曲型方程，设置初始边界条件后可以求出给定时刻t的解；对于特征值问题，可以求出给定区间上的特征值.求解完成后，可以返回到Step4，对网格进一步细化，进行再次求解.Step6“View模式”计算结果的可视化，可以通过设置系统提供的对话框，显示所求的解的表面图、网格图、等高线图和箭头梯形图.对于抛物线型和双曲线型问题的解还可以进行动画演示.