计算流体力学算例

有限差分法的一个算例——计算流体力学大作业作者：郝柏函2010011545指导：李嵩题目编程计算热传导方程,边界条件：初始条件：用FTCS格式分别在满足和不满足稳定性条件两种情况下计算，给出结果比较和分析。（2）自选一种其他格式编程计算，并给出结果和分析注：原题中给出的初始条件与边界条件是矛盾的，所以将其改为MACROBUTTONMTEditEquationSection2SEQMTEqn\r\hSEQMTSec\hFTCS格式计算方法差分格式及其相容性对于方程，采用FTCS差分格式，即 其中，。以下讨论这一格式的相容性。所以 因此，该格式与原微分方程是相容的，而且对于时间精度是一阶的，对于空间，精度是二阶的。稳定性与收敛性对于适定的线性微分方程，格式如果差分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了。设则GOTOBUTTONZEqnNum641444REFZEqnNum641444\*Charformat\!(2.1)式可写为 即 放大因子 所以 为保证,应有 只要满足GOTOBUTTONZEqnNum292638REFZEqnNum292638\*Charformat\!(2.2)式，差分格式就是稳定的。初始、边界条件处理以及全部计算过程初始条件，差分格式为 边界条件，为了保证空间的二阶精度，采用二次多项式来构造差分格式，结果为，其中是方向上位置的格点数。给出了初始、边界条件，以及之前的差分格式，就可以给出完整地算法：（1）首先用GOTOBUTTONZEqnNum224877REFZEqnNum224877\*Charformat\!(2.3)式计算出第一个时层的温度；（2）然后使用GOTOBUTTONZEqnNum641444REFZEqnNum641444\*Charformat\!(2.1)式就算出下一时层的温度值，但是，此时还没有就算出，然后利用边界条件求出（3）不断使用第（2）步，直至计算出所要求时层所对应的温度值注：为了保证计算效率，不应过小。如果要求计算结果是稳定的，应满足GOTOBUTTONZEqnNum292638REFZEqnNum292638\*Charformat\!(2.2)式，如果要求不稳定，应不满足GOTOBUTTONZEqnNum292638REFZEqnNum292638\*Charformat\!(2.2)式。计算结果与分析本文采用matlab编程，程序见于第REF_Ref377047466\r\h4小节。在不稳定的差分格式下，计算结果是不可采信的，如REF_Ref377052026\h图1所示。图SEQ图\*ARABIC1不稳定格式计算得到温度分布，s=1.04,其中，时间采用2400步，空间采用100步采用不稳定格式虽然也能得到比较光滑的温度分布图，但是，根据本算例的物理意义，左端为恒定温度0，右端为绝热壁面，所以计算结果应该是，温度始终大于0，别且距左端越近，温度越低。可见，非稳定格式的计算结果是定性错误的。而稳定格式的计算结果是可以采信的。如REF_Ref377054156\h图2，REF_Ref377054159\h图3，REF_Ref377054161\h图4所示。图SEQ图\*ARABIC2稳定格式计算得到温度分布，s=0.05,其中，时间采用50,000步，空间采用100步图SEQ图\*ARABIC3稳定格式计算得到温度分布，s=0.005,其中，时间采用500,000步，空间采用100步图SEQ图\*ARABIC4稳定格式计算得到温度分布，s=0.0005,其中，时间采用5,000,000步，空间采用100步这三个计算结果相对于之前的不稳定计算结果，只是增大改变了时间步数。使得GOTOBUTTONZEqnNum292638REFZEqnNum292638\*Charformat\!(2.2)式得以满足。但仅仅是这一条件的改变，使得之前所描述的定性结果是正确的。但是，这三个计算结果也是有微弱的差别的。仅仅看时，处，三个解算结果温度值是不同的，分别为0.5226，0.5380，0.5395。尽管，这三个数值之间有微弱的差别，但是整体上来说趋近于0.54这个数值。而且温度分布的整体趋势、数值之间的差距也几乎为0.这说明计算确实是稳定的、收敛的。而且，时间步数越多，结果会越精确。FTCS隐格式计算方法差分格式及其相容性采用FTCS隐格式，即 其中，。以下讨论这一格式的相容性。所以 其中，均在处取值因此，该格式与原微分方程是相容的，而且对于时间精度是一阶的，对于空间，精度是二阶的。稳定性与收敛性对于适定的线性微分方程，格式如果差分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了。设则GOTOBUTTONZEqnNum641444REFZEqnNum641444\*Charformat\!(2.1)式可写为 即 放大因子 所以 该差分格式是无条件稳定的。初始、边界条件处理以及全部计算过程初始条件，差分格式为 边界条件，为了保证空间的二阶精度，采用二次多项式来构造差分格式，结果为，其中是方向上位置的格点数。给出了初始、边界条件，以及之前的差分格式，就可以给出完整地算法：（1）首先用GOTOBUTTONZEqnNum224877REFZEqnNum224877\*Charformat\!(2.3)式计算出第一个时层的温度；（2）然后使用GOTOBUTTONZEqnNum771453REFZEqnNum771453\*Charformat\!(2.4)式以及，组成一个m元一次方程组，对于这个方程组，以第n时层的温度值为初值，采用迭代法计算，直至误差足够小。即 其中，k为迭代次数。（3）不断使用第（2）步，直至计算出所要求时层所对应的温度值。计算结果与分析为了与之前的显格式进行比较，这里取空间步数为100，时间步数则分别取50,000，500,000，5,000,000，得到如所示的结果。每一个时间步长的计算精度都是。这三个结果都满足：正温度；越靠近原点处，温度越低；原点处温度为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。图SEQ图\*ARABIC5隐格式下t=0.5s时刻的温度分布，时间步长50,000图SEQ图\*ARABIC6隐格式下t=0.5s时刻的温度分布，时间步长500,000图SEQ图\*ARABIC7隐格式下t=0.5s时刻的温度分布，时间步长5,000,000相较于显格式，隐格式很快收敛到了0.5396，然而，显格式即使步数为5,000,000，也只能收敛到0.5397。由此可见，隐格式具有如下优越性（1）无条件稳定，不需要考虑时间步长和空间步长的关系；（2）收敛步数短。另外，我们从以上数据还可以看出，对本算例而言，有如下特点：显格式的计算结果相对于精确值偏小，而隐格式的计算结果偏大。本文程序在压缩包中，应该有.m文件的程序，fluid_conpution.m为显格式的计算程序，fluid_conpution2.m为隐格式的计算程序。FTCS格式程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%FTCS格式计算热传导方程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clcclearall%1.定义参数alpha=2;%热传导微分方程的系数T=0.5;%计算时长，可修改D=2;%计算区域宽度L1=500000;%时间网格数L2=100;%空间网格数s=1;%差分格式的系数whiles>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,'k--'),holdon%3.迭代计算下一时层的温度forii=1:L1forjj=2:L2u(jj)=s*((u(jj+1)-u(jj))-(u(jj)-u(jj-1)))+u(jj);%下一时层几乎所有格点处的温度值endu(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'),ylabel('u/^{o}C')legend('t=0时刻的温度分布','t=0.5s时刻的温度分布')FTCS隐格式程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%FTCS隐格式计算热传导方程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clcclearall%1.定义参数alpha=2;%热传导微分方程的系数T=0.5;%计算时长，可修改D=2;%计算区域宽度L1=5e2;%时间网格数L2=100;%空间网格数eps=1e-16;%每一时层的精度要求dt=T/