计算流体力学.docxVIP

数值计算一维非稳态导热，长度一米，两端初始温度分别为0度和300度的不锈钢棒，采用显示差分方法和Crank-Nicolson方法，分别计算不锈钢棒内的温度分布变化情况，显示差分方法的时间步长选取符合稳定性条件和非稳定性条件两种，观察计算结果并分析。首先，我们做出假设：令不锈钢棒的初始温度为0度，两端的热源为恒温热源，热量不断地从两端进行热传导过程；棒的导热系数设定为1，棒为均匀的导热系数保持不变的不锈钢棒。我们可以由题意和假设条件，给出棒的初始和边界条件： （1）空间一维的热传导方程： （2）显示差分方法为： （3）克兰克-尼克尔森的隐式差分方法为： （4）从两个种差分格式的表达上，我们明显看出差分格式的不同，在显示方法中，每一个差分方程只包含一个未知数，因此这个未知数可以用直接计算的方式显示的求解，在隐式方法中，我们必须对排列在同一时间层上的所有网格点上的未知参量，必须将他们联立起来同时求解，才能求出这些未知量，这种方法定义为隐式方法。本文的主要工作是通过Matlab编程，将计算得到的矩阵数据进行整理分析，并在Matlab中画出特定的图像。对显示差分方法和隐式差分方法进行比较，给出显示差分方法分稳定性和不稳定性两种结果的对比，画出不稳定性条件下的震荡曲线，证明出隐式差分格式有着更好的时间上的稳定性，对时间步长要求不严格，有很好的计算效率。图1通过MATLAB将数据进行三维的展示，我只将隐式的最后分布图表示出来，主要是因为，显示和隐式的最终分布基本一致。如上图1所示，红色代表高温区，蓝色代表低温区，最后的温度分布趋近于一条斜线，这与我们所设想的是一致的，高温区向低温热传导，并且对于在高温区附近的温度变化来看，在开始的时候呈现快速而又复杂的变化情况。 （a） （b）图2 显示（a）和隐式（b）中点温度变化情况不锈钢棒的中点温度都是由初始条件的273.15k，慢慢升到最后的平衡温度423.15度，隐式差分格式计算时，数值略小于423.15k，但是误差极小。如图2所示，证明了两种差分格式的一致正确性。图3不稳定条件下的震荡在之前用显示计算的时候都是保证了稳定性条件，即 ，我们改变时间的步长，将时间步长从0.001增加到0.01，这就导致稳定系数超过稳定条件，即的时候，会出现计算不稳定的情况，如图3所示，这就证明了，对于显示有限差分格式来说，稳定性要求对推进变量有所限制。要保证时间步长很小，才能保持稳定。 （a） （b） （c） （d） 图 4 对于隐式计算时间步长0,001（a），时间步长0,01（b），时间步长0.05（c）和时间步长0.1（d）中点的温度变化。我们可以发现，隐式的差分格式，用大的多的时间步长也能够保持稳定性，要将时间变量计算到给定的时间值，只需要较少的时间步，这样，计算机的运行程序时间更短。但是，对于隐式的差分格式，它的计算程序略有复杂，对于复杂的运算情况，我们应该自己做出最佳的判断。本文从实际问题出发，对热传导方程的显示差分格式和隐式差分格式进行编程计算，提高了对有限差分计算的认识。并且比较了两种差分格式的优缺点，为之后的学习，提供了很好的铺垫。也能够更好的熟悉商业化软件的计算特点，给之后的fluent模拟打下基础。附录：显示差分格式分MATLAB计算程序：dx=0.1;%设定空间步长dt=0.001;%设定时间步长m=10;n=2000;t=zeros(n+1,m+1);t(1,:)=273.15;%设定边界条件t(:,1)=273.15;t(:,m+1)=573.15; %设定初始条件arf=1;stable=arf*dt/dx^2; %稳定性条件 for j= 1:n+1 %进行循环计算 for i=2:m t(j+1,i)=stable*(t(j,i-1)+t(j,i+1))+(-2*stable+1)*t(j,i); endend %最后得到矩阵T，可以再通过Matlab进行画图隐式差分格式的MATLAB计算程序：dx=0.1;%设定空间步长dt=0.001;%设定时间步长arfa=1;N=1/dx+1;%设定空间计算点位n=2000;%设定时间计算点位T=zeros(n,N);K=zeros(n,N);D=zeros(n,N);T(1,:)=0;%设定边界条件T(:,1)=273.15;%设定初始条件T(:,N)=573.15; A=(arfa*dt)/(2*dx*dx);B=1+(arfa*dt)/(dx*dx);for j=1:n%进行循环计算（托马斯算法） for i=2:(N-1) K(j,i)=-T(j,i)-((arfa*dt)/(2*dx*dx))*(T(j,i+1)-2*T(j,i)+T(j,i-1)); end K(j,2)=K(j,2)-A*T