计算流体力学算例

1、有限差分法的一个算例计算流体力学大作业作者：郝柏函2010011545指导：李嵩1. 题目编程计算热传导方程, 边界条件：初始条件：（1） 用FTCS格式分别在满足和不满足稳定性条件两种情况下计算，给出结果比较和分析。（2）自选一种其他格式编程计算，并给出结果和分析注：原题中给出的初始条件与边界条件是矛盾的，所以将其改为2. FTCS格式2.1. 计算方法2.1.1. 差分格式及其相容性对于方程，采用FTCS差分格式，即 其中，。以下讨论这一格式的相容性。 所以因此，该格式与原微分方程是相容的，而且对于时间精度是一阶的，对于空间，精度是二阶的。2.1.2. 稳定性与收敛性对于适定的线性微分方程

2、，格式如果差分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了。设则式可写为 即 放大因子 所以为保证,应有 只要满足式，差分格式就是稳定的。2.1.3. 初始、边界条件处理以及全部计算过程初始条件，差分格式为 边界条件，为了保证空间的二阶精度，采用二次多项式来构造差分格式，结果为，其中是方向上位置的格点数。给出了初始、边界条件，以及之前的差分格式，就可以给出完整地算法：（1）首先用式计算出第一个时层的温度；（2）然后使用式就算出下一时层的温度值，但是，此时还没有就算出，然后利用边界条件求出（3）不断使用第（2）步，直至计算出所要求时层所对应的温度值注：为了保证计算效率，不应过小。如

3、果要求计算结果是稳定的，应满足式，如果要求不稳定，应不满足式。2.2. 计算结果与分析本文采用matlab编程，程序见于第4小节。 在不稳定的差分格式下，计算结果是不可采信的，如图 1所示。图 1 不稳定格式计算得到温度分布，s=1.04,其中，时间采用2400步，空间采用100步采用不稳定格式虽然也能得到比较光滑的温度分布图，但是，根据本算例的物理意义，左端为恒定温度0，右端为绝热壁面，所以计算结果应该是，温度始终大于0，别且距左端越近，温度越低。可见，非稳定格式的计算结果是定性错误的。而稳定格式的计算结果是可以采信的。如图 2，图 3，图 4所示。图 2 稳定格式计算得到温度分布，s=0.

4、05,其中，时间采用50,000步，空间采用100步图 3 稳定格式计算得到温度分布，s=0.005,其中，时间采用500,000步，空间采用100步图 4 稳定格式计算得到温度分布，s=0.0005,其中，时间采用5,000,000步，空间采用100步这三个计算结果相对于之前的不稳定计算结果，只是增大改变了时间步数。使得式得以满足。但仅仅是这一条件的改变，使得之前所描述的定性结果是正确的。但是，这三个计算结果也是有微弱的差别的。仅仅看时，处，三个解算结果温度值是不同的，分别为0.5226，0.5380，0.5395。尽管，这三个数值之间有微弱的差别，但是整体上来说趋近于0.54这个数值。而且

5、温度分布的整体趋势、数值之间的差距也几乎为0.这说明计算确实是稳定的、收敛的。而且，时间步数越多，结果会越精确。3. FTCS隐格式3.1. 计算方法3.1.1. 差分格式及其相容性采用FTCS隐格式，即 其中，。以下讨论这一格式的相容性。 所以其中，均在处取值因此，该格式与原微分方程是相容的，而且对于时间精度是一阶的，对于空间，精度是二阶的。3.1.2. 稳定性与收敛性对于适定的线性微分方程，格式如果差分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了。设则式可写为 即 放大因子 所以该差分格式是无条件稳定的。3.1.3. 初始、边界条件处理以及全部计算过程初始条件，差分格式为 边

6、界条件，为了保证空间的二阶精度，采用二次多项式来构造差分格式，结果为，其中是方向上位置的格点数。给出了初始、边界条件，以及之前的差分格式，就可以给出完整地算法：（1）首先用式计算出第一个时层的温度；（2）然后使用式以及，组成一个m元一次方程组，对于这个方程组，以第n时层的温度值为初值，采用迭代法计算，直至误差足够小。即其中，k为迭代次数。（3）不断使用第（2）步，直至计算出所要求时层所对应的温度值。3.2. 计算结果与分析为了与之前的显格式进行比较，这里取空间步数为100，时间步数则分别取50,000， 500,000， 5,000,000，得到如所示的结果。每一个时间步长的计算精度都是。这三

7、个结果都满足：正温度；越靠近原点处，温度越低；原点处温度为0；x=2处绝热等定性的条件。而且计算结果并没有很大的区别。三种计算结果下，在0.5s时，x=2处的温度值分别为，0.5396，0.5396，0.5396，几乎没有什么区别。在其他点处，计算结果也是几乎完全相同。其实，对于隐格式，时间步长也不需要到50,000这么多。时间步数为5,000，在0.5s时，x=2处的温度值已经是0.5397，步数为500，该温度值已经是0.5399。图 5 隐格式下t=0.5s时刻的温度分布，时间步长50,000图 6 隐格式下t=0.5s时刻的温度分布，时间步长500,000图 7 隐格式下t=0.5s时

8、刻的温度分布，时间步长5,000,000相较于显格式，隐格式很快收敛到了0.5396，然而，显格式即使步数为5,000,000，也只能收敛到0.5397。由此可见，隐格式具有如下优越性 （1）无条件稳定，不需要考虑时间步长和空间步长的关系；（2）收敛步数短。另外，我们从以上数据还可以看出，对本算例而言，有如下特点：显格式的计算结果相对于精确值偏小，而隐格式的计算结果偏大。4. 本文程序在压缩包中，应该有.m文件的程序，fluid_conpution.m为显格式的计算程序，fluid_conpution2.m为隐格式的计算程序。4.1. FTCS格式程序% FTCS格式计算热传导方程%clccl

9、ear all% 1.定义参数alpha=2; % 热传导微分方程的系数T=0.5; % 计算时长，可修改D=2; % 计算区域宽度L1=500000; % 时间网格数L2=100; % 空间网格数s=1; % 差分格式的系数while s>1/2 % 减小时间步长，以保证s<=0.5，稳定情况下适用 L1=L1*10; dt=T/L1; % 时间步长 dx=D/L2; % 空间步长 s=alpha*dt/(dx)2end % dt=T/L1;% dx=D/L2;% s=alpha*dt/(dx)2 % 2.给出初始条件x=0:dx:D;u=sin(pi*x/4);plot(x,u

10、,'k-'),hold on % 3.迭代计算下一时层的温度for ii=1:L1 for jj=2:L2 u(jj)=s*(u(jj+1)-u(jj)-(u(jj)-u(jj-1)+u(jj); % 下一时层几乎所有格点处的温度值 end u(1)=0; % 下一时层x=0处的温度值 u(L2+1)=(4*u(L2)-u(L2-1)/3; % 下一时层x=D处的温度值end % 4.绘制计算结果 plot(x,u,'k-')title('稳定格式下t=0.5s时刻的温度分布，s=0.0005')xlabel('x/m'),yla

11、bel('u/oC')legend('t=0时刻的温度分布','t=0.5s时刻的温度分布')4.2. FTCS隐格式程序% FTCS隐格式计算热传导方程%clcclear all% 1.定义参数alpha=2; % 热传导微分方程的系数T=0.5; % 计算时长，可修改D=2; % 计算区域宽度L1=5e2; % 时间网格数L2=100; % 空间网格数eps=1e-16; % 每一时层的精度要求dt=T/L1; % 时间步长dx=D/L2; % 空间步长s=alpha*dt/(dx)2% 2.给出初始条件x=0:dx:D;u=sin(pi*x/4);plot(x,u,'k-'),hold on % 3.迭代法计算下一时层的温度 for ii=1:L1 epsilon=1; % 每一时层的精度 u(1)=0; % 下一时层x=0处的温度值 w=u; while epsilon>eps v=w; for jj=2:L2 w(jj)=(s*(w(jj+1)+w(jj-1)+u(jj)/(1+2*s); % 下一时层几乎所有格点处的温度值 end w(L2+1)=(4*w(L2)-w(L2-1)/3; % 下一时层x=D处的温度值 epsilon=max(abs(w-v);