第一类边界条件非稳态传热方程计算代码.doc

第一类边界条件非稳态传热方程计算代码%用于模拟计算第一类边界条件非稳态传热方程%Uxx(x,t)=Ut(x,t)%u(0,t)=u(1,t)-0%u(x,0)=sin(pi*x)%采用Crank-Nicolson方法，%代码根据Numerical Mathematics and Computing 4th ed伪代码编写LENGTH=1; %空间长度TIME=0.1;%瞬态模拟时间del_x=0.1;del_t=0.005;n=LENGTH/del_x; m=TIME/del_t;s=del_x2/del_t;gamma=2+s;for i=1:n-1d(i)=gamma;c(i)=-1;%t=0时刻，u赋初值u(i)=sin(pi*i*del_x);endP=zeros(m+1,n+1);%存储所有时间层与空间层交叉点温度%外部循环为时间节点，内部循环为空间节点for j=1:m+1for i=1:n-1d(i)=gamma;v(i)=s*u(i);endv=tri(n-1,c(1:n-2),d,c(1:n-2),v);u=v;P(j,:)=0,v,0;t=j*del_t;%u=v;endX,Y=meshgrid(0:del_x:LENGTH,0:del_t:TIME);mesh(X,Y,P);%下面是调用求解三对角矩阵的子程序function x=tri(n,A,D,C,b)%求解三对角方程：Ax=b%A-下对角向量%D-主对角向量%C-上对角向量%b-等式右边向量for i=2:nxmult=A(i-1)/D(i-1);D(i)=D(i)-xmult*C(i-1);b(i)=b(i)-xmult*b(i-1);endx(n)=b(n)/D(n);for i=(n-1) :-1: 1x(i)=(b(i)-C(i)*x(i+1)/D(i);endx=x;一维非稳态热传导显式计算（C+）第一类边界条件 2010年10月30日 星期六 12:45/显式计算一维非定常热传导方程#include #include using namespace std;int main() double T11; /存放初始值和边界值 double result11; /存放结果 double L=1; /几何模型 double deltaX=0.1; /空间步长 double deltaT=0.001; /时间步长 int n=L/deltaX+1; /空间节点数 for(int i=1;i10;i+)/初始化 if(abs(i+1)*deltaX) - 0.5*L 1e-9) Ti=2*(i)*deltaX; else Ti=2*(1-(i)*deltaX); T0=0; Tn=0; /显示初始数据 cout time: 0*deltaT endl; for(int j=0;j11;j+) /时间迭代 coutTj ; coutendlendl; for(int i=1;i11;i+) /时间迭代 result0=T0; resultn=Tn; for(int j=1;j10;j+) /空间迭代 resultj=Tj+deltaT*(Tj+1-2*Tj+Tj-1)/deltaX/deltaX; Tj=resultj; T0=0; Tn=0; cout time: i*deltaT endl; for(int j=0;j11;j+) coutresultj ; coutendlendl; -后记-此程序非常的简单，但是有两点需要注意：（1）空间数据用一个一维数组，时间迭代采用另外一个一维数组，这样有利于内存节省。当然也可以开辟一个二维数组将所有网格包含进去，虽然在程序编写上要方便一些，但是比较耗费内存，特别是在网格比较多的情况下。（2）显式计算需要满足CFL条件，简单的表示应该是(时间步长)/(空间步长*空间步长)0.5，否则会导致解的震荡或非物理解。本例中该值为0.1/显式计算一维非定常热传导方程#include #include #include using namespace std;int main() double T11; /存放初始值和边界值 double result11; /存放结果 double L=1; /几何模型 double deltaX=0.1; /空间步长 double deltaT=0.001; /时间步长 int n=L/deltaX+1; /空间节点数 for(int i=1;i10;i+)/初始化 if(abs(i+1)*deltaX) - 0.5*L 1e-9) Ti=2*(i)*deltaX; else Ti=2*(1-(i)*deltaX); T0=0; Tn=0; /显示初始数据 cout time: 0*deltaT endl; for(int j=0;j11;j+) /时间迭代 coutTj ; coutendlendl; for(int i=1;i11;i+) /时间迭代 result0=T0+2*deltaT*(T1-(1+deltaX)*T0)/deltaX/deltaX; resultn=Tn+2*deltaT*(Tn-1-(1+deltaX)*Tn)/deltaX/deltaX; for(int j=1;j10;j+) /空间迭代 resultj=Tj+deltaT*(Tj+1-2*Tj+Tj-1)/deltaX