1、有限差分法的一个算例计算流体力学大作业作者：郝柏函 2010011545指导：李嵩1. 题目编程计算热传导方程, 2utx0,2=tx边界条件： 02,0xxu初始条件： =0sin4tu（1） 用 FTCS 格式 njnjnjnj suss 112 njjnjjj uxtu1121 分别在满足和不满足稳定性条件两种情况下计算，给出结果比较和分析。（2）自选一种其他格式编程计算，并给出结果和分析注：原题中给出的初始条件 与边界条件 是矛盾的，所以将其=0sin2tx20xu改为 =0sin4txu2. FTCS 格式2.1. 计算方法2.1.1.差分格式及其相容性对于方程 ，采用 FTCS 差
2、分格式 ，2utx 1112njjnnjjjuutx即* MERGEFORMAT 1 12nnjj jjuss (2.1)其中， 。以下讨论这一格式的相容性。tsx122112njjnnjjjuOtttuxx所以 1 22112njjnnjjjuuuOtxtxtx 因此，该格式与原微分方程是相容的，而且对于时间精度是一阶的，对于空间，精度是二阶的。2.1.2.稳定性与收敛性对于适定的线性微分方程，格式如果差分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了。设 ijnjure则* MERGEFORMAT (2.1) 式可写为1 112ijij ijijnnnnresrsres 即1i
3、 ine放大因子122cos14inni iserG所以 1sG为保证 ,应有1G* MERGEFORMAT (2.2)21txs只要满足* MERGEFORMAT (2.2) 式，差分格式就是稳定的。2.1.3.初始、边界条件处理以及全部计算过程初始条件 ，差分格式为=0sin2txu* MERGEFORMAT (2.3)1=sin2jjxu边界条件 ，为了保证空间的二阶精度，采用二次多项式来构造差02,0xxu分格式，结果为 ，其中 是 方向上位置的格点数。1214,3nnmux给出了初始、边界条件，以及之前的差分格式，就可以给出完整地算法：（1）首先用* MERGEFORMAT (2.3
4、) 式计算出第一个时层的温度；（2）然后使用* MERGEFORMAT (2.1)式就算出下一时层的温度值，但是，此时还没有就算出，然后利用边界条件求出1,nmu111240,3nnnmuu（3）不断使用第（2）步，直至计算出所要求时层所对应的温度值注：为了保证计算效率， 不应过小。如果要求计算结果是稳定的，应满足* 2stxMERGEFORMAT (2.2)式，如果要求不稳定，应不满足* MERGEFORMAT (2.2) 式。2.2. 计算结果与分析本文采用 matlab 编程，程序见于第 4 小节。在不稳定的差分格式下，计算结果是不可采信的，如 图 1 所示。图 1 不稳定格式计算得到温
5、度分布，s=1.04,其中，时间采用 2400 步，空间采用 100 步采用不稳定格式虽然也能得到比较光滑的温度分布图，但是，根据本算例的物理意义，左端为恒定温度 0，右端为绝热壁面，所以计算结果应该是，温度始终大于 0，别且距左端越近，温度越低。可见，非稳定格式的计算结果是定性错误的。而稳定格式的计算结果是可以采信的。如 图 2， 图 3， 图 4 所示。图 2 稳定格式计算得到温度分布， s=0.05,其中，时间采用 50,000 步，空间采用 100 步图 3 稳定格式计算得到温度分布， s=0.005,其中，时间采用 500,000 步，空间采用 100 步图 4 稳定格式计算得到温度
6、分布， s=0.0005,其中，时间采用 5,000,000 步，空间采用 100 步这三个计算结果相对于之前的不稳定计算结果，只是增大改变了时间步数。使得* MERGEFORMAT (2.2)式得以满足。但仅仅是这一条件的改变，使得之前所描述的定性结果是正确的。但是，这三个计算结果也是有微弱的差别的。仅仅看 时， 处，三个解0.5ts2xm算结果温度值是不同的，分别为 0.5226，0.5380 ，0.5395。尽管，这三个数值之间有微弱的差别，但是整体上来说趋近于 0.54 这个数值。而且温度分布的整体趋势、数值之间的差距也几乎为 0.这说明计算确实是稳定的、收敛的。而且，时间步数越多，结
7、果会越精确。3. FTCS 隐格式3.1. 计算方法3.1.1.差分格式及其相容性采用 FTCS 隐格式 ，即1112njjnnjjjuutx* MERGEFORMAT 111nnnj jjjssu(2.4)其中， 。以下讨论这一格式的相容性。2tsx1 21122njjnnjjjuOtttuxx所以 1 211 22njjnnjjjuuuOtxtxtx 其中， 均在 处取值2,t1nj因此，该格式与原微分方程是相容的，而且对于时间精度是一阶的，对于空间，精度是二阶的。3.1.2.稳定性与收敛性对于适定的线性微分方程，格式如果差分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了。设
8、ijnjure则* MERGEFORMAT (2.1) 式可写为1112ij ijijijnnnnsrerse 即1con放大因子12csonsrG所以 102Gs该差分格式是无条件稳定的。3.1.3.初始、边界条件处理以及全部计算过程初始条件 ，差分格式为=0sin2txu* MERGEFORMAT (2.5)1=sin2jjxu边界条件 ，为了保证空间的二阶精度，采用二次多项式来构造差02,0xxu分格式，结果为 ，其中 是 方向上位置的格点数。1214,3nnmx给出了初始、边界条件，以及之前的差分格式，就可以给出完整地算法：（1）首先用* MERGEFORMAT (2.3) 式计算出第
9、一个时层的温度；（2）然后使用* MERGEFORMAT (2.4)式以及 ，组成一111240,3nnnmuu个 m 元一次方程组，对于这个方程组，以第 n 时层的温度值为初值，采用迭代法计算，直至误差足够小。即 1,0,1, 1, 1, 2,2,2143jjnknknkj jjjnkmususjmu 其中，k 为迭代次数。（3）不断使用第（2）步，直至计算出所要求时层所对应的温度值。3.2. 计算结果与分析为了与之前的显格式进行比较，这里取空间步数为 100，时间步数则分别取 50,000， 500,000， 5,000,000，得到如所示的结果。每一个时间步长的计算精度都是 。这三个结果
10、都满足：正温度；越靠近原点处，温度越低；160原点处温度为 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,0
11、00图 7 隐格式下 t=0.5s 时刻的温度分布，时间步长 5,000,000相较于显格式，隐格式很快收敛到了 0.5396，然而，显格式即使步数为5,000,000，也只能收敛到 0.5397。由此可见，隐格式具有如下优越性（1）无条件稳定，不需要考虑时间步长和空间步长的关系；（2）收敛步数短。另外，我们从以上数据还可以看出，对本算例而言，有如下特点：显格式的计算结果相对于精确值偏小，而隐格式的计算结果偏大。4. 本文程序在压缩包中，应该有.m 文件的程序， fluid_conpution.m 为显格式的计算程序，fluid_conpution2.m 为隐格式的计算程序。4.1. FTCS
12、 格式程序% FTCS格式计算热传导方程%clcclear all% 1.定义参数alpha=2; % 热传导微分方程的系数T=0.5; % 计算时长，可修改D=2; % 计算区域宽度L1=500000; % 时间网格数L2=100; % 空间网格数s=1; % 差分格式的系数while s1/2 % 减小时间步长，以保证sepsv=w;for jj=2:L2w(jj)=(s*(w(jj+1)+w(jj-1)+u(jj)/(1+2*s); % 下一时层几乎所有格点处的温度值endw(L2+1)=(4*w(L2)-w(L2-1)/3; % 下一时层x=D处的温度值epsilon=max(abs(w-v);endu=w;end% 4.绘制计算结果plot(x,u,k-)title(隐格式下t=0.5s时刻的温度分布，时间步长500,000)xlabel(x/m),ylabel(u/oC)legend(t=0时刻的温度分布,t=0.5s时刻的温度分布)
