计算流体力学算例

有限差分法的一个算例

——计算流体力学大作业

作者：郝柏函2010011545

指导：李嵩

1. 题目

编程计算热传导方程

22u u

t x

α??=??, 0,[0,2],=2t x α≥∈ 边界条件：0

2

0,0x x u

u x ==?==?

初始条件：

=0=sin

4

t x

u π

（1） 用FTCS 格式

()n

j n j n j n j su u s su

u 11

121-+++-+=

()n

j n j n

j n

j

n j u u u

x

t

u u 1

12

12-+++-?=

?-α

分别在满足和不满足稳定性条件两种情况下计算，给出结果比较和分析。 （2）自选一种其他格式编程计算，并给出结果和分析 注：原题中给出的初始条件=0=sin 2

t x

u π与边界条件

2

0x u

x

=?=?是矛盾的，所以将其改

为=0=sin

4

t x

u π

2. FTCS 格式

2.1. 计算方法

2.1.1. 差分格式及其相容性

对于方程22u u t x

α??=??，采用FTCS 差分格式()11122n n j j n n n

j j j u u u u u t x α++--=-+??，

即

()11112n n n n

j j j j u su s u su ++-=+-+

(1.1)

其中，2

t

s x α?=

?。以下讨论这一格式的相容性。

()()()12

2

11

22

12n n j

j n n n j j j u u u

O t t

t

u u u u O x x x ++--?=

+????-+=+???

所以

()()()1

2211222n n j j n n n

j j j u u u u u u u O t O x t x t

x ααα++-??-????--+--=?+? ? ? ????????? 因此，该格式与原微分方程是相容的，而且对于时间精度是一阶的，对于空

间，精度是二阶的。 2.1.2. 稳定性与收敛性

对于适定的线性微分方程，格式如果差分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了。

设()

i j n j n u r e

θ= 则(1.1)式可写为

(

)

(

)()(

)

(

)11112i j i j i j i j n n n n r e sr e s r e sr e θθ

θθ

+-+=+-+

即

()()112i i n n r r se s se θθ-+=+-+

放大因子

()()1

2

12122cos 121cos 14sin 2

n i i n

se s se s r s r s s G θθθθθ

-++-+=-+=--=-=

=

所以

141s G -≤≤

为保证1G ≤,应有

2

12

t

x s α?=

?≤

(1.2)

只要满足(1.2)式，差分格式就是稳定的。

2.1.

3. 初始、边界条件处理以及全部计算过程

初始条件=0=sin

2

t x

u π，差分格式为

1=sin

2

j j x

u π? (1.3)

边界条件02

0,

0x x u u x

==?==?，为了保证空间的二阶精度，采用二次多项式来构造差分

格式，结果为12

1

40,3

n n n n m m m

u u u u ---==，其中m 是x 方向上位置的格点数。

给出了初始、边界条件，以及之前的差分格式，就可以给出完整地算法： （1）首先用(1.3)式计算出第一个时层的温度；

（2）然后使用(1.1)式就算出下一时层的温度值，但是，此时11

1,n n m u u ++还没有就算出，

然后利用边界条件求出111

112

1

40,3

n n n n m m m

u u u

u

++++---==

（3）不断使用第（2）步，直至计算出所要求时层所对应的温度值 注：为了保证计算效率，2

s t

x

α?=

?不应过小。如果要求计算结果是稳定的，应满足(1.2)

式，如果要求不稳定，应不满足(1.2)式。

2.2. 计算结果与分析

本文采用matlab 编程，程序见于第4小节。

在不稳定的差分格式下，计算结果是不可采信的，如图 1所示。

图1 不稳定格式计算得到温度分布，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 稳定格式计算得到温度分布，s=0.0005,其中，时间采用5,000,000步，空间采用100步

这三个计算结果相对于之前的不稳定计算结果，只是增大改变了时间步数。使得(1.2)式得以满足。但仅仅是这一条件的改变，使得之前所描述的定性结果是正确的。

但是，这三个计算结果也是有微弱的差别的。仅仅看0.5t s =时，2x m =处，三个解算结果温度值是不同的，分别为0.5226，0.5380，0.5395。

尽管，这三个数值之间有微弱的差别，但是整体上来说趋近于0.54这个数值。而且温度分布的整体趋势、数值之间的差距也几乎为0.这说明计算确实是稳定的、收敛的。而且，时间步数越多，结果会越精确。

3. FTCS 隐格式

3.1. 计算方法

3.1.1. 差分格式及其相容性

采用FTCS 隐格式()11

11

1

12

2n n j

j n n n j j j u u u

u u t

x α

+++++--=

-+??，即

()1

11

11112n n n n j

j j j u su u su s

++++-=+++ (1.4)

其中，2

t

s x α?=

?。以下讨论这一格式的相容性。

()()()

12

111

2

1122

12n n j

j n n n j j j u u u

O t t

t

u u u u O x x x

+++++--?=

+????-+=+???

所以

()()()1

2111

211222n n j j n n n j j j u u u u u u u O t O x t x t

x ααα+++++-??-????--+--=?+?

? ? ????????? 其中，22,u u t x

????均在1

n j u +处取值

因此，该格式与原微分方程是相容的，而且对于时间精度是一阶的，对于空间，精度是二阶的。 3.1.2. 稳定性与收敛性

对于适定的线性微分方程，格式如果差分格式，那么稳定和收敛是等价的。所以只需要讨论稳定性就可以了。

设()

i j n j n u r e

θ= 则(1.1)式可写为

(

)

()()()()

11111112i j i j i j i j n n n n r e sr e r e sr e s

θθθθ

+-+++=

+++ 即

()1122cos n n s s r r θ++-=

放大因子

()

11

122cos 1

121cos n n s s s r G r θ

θ++-=

+-=

=

所以

1

0121G s

<

+≤≤ 该差分格式是无条件稳定的。

3.1.3. 初始、边界条件处理以及全部计算过程

初始条件=0=sin

2

t x

u π，差分格式为

1=sin

2

j j x

u π? (1.5)

边界条件02

0,

0x x u u x

==?==?，为了保证空间的二阶精度，采用二次多项式来构造差分

格式，结果为12

1

40,3

n n n n m m m

u u u u ---==，其中m 是x 方向上位置的格点数。

给出了初始、边界条件，以及之前的差分格式，就可以给出完整地算法： （1）首先用(1.3)式计算出第一个时层的温度； （2）然后使用(1.4)式以及111112

1

40,3

n n n n m m m

u u u

u

++++---==

，组成一个m 元一次方程组， 对于这个方程组，以第n 时层的温度值为初值，采用迭代法计算，直至误差足够小。

即

()1,01,11,11,1,1

11

111,112

,1,2,,0

1,1,2,,11243

n n j j n k n k n k n n k j

j j j n n n k m m m u u j m

u u su u su j m s u u u ++++++++-++++--====++=-+-=

其中，k 为迭代次数。

（3）不断使用第（2）步，直至计算出所要求时层所对应的温度值。

3.2. 计算结果与分析

为了与之前的显格式进行比较，这里取空间步数为100，时间步数则分别取50,000， 500,000， 5,000,000，得到如所示的结果。每一个时间步长的计算精度都是16110-?。这三个结果都满足：正温度；越靠近原点处，温度越低；原点处温度为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时刻的温度分布，时间步长5,000,000相较于显格式，隐格式很快收敛到了0.5396，然而，显格式即使步数为5,000,000，也只能收敛到0.5397。

由此可见，隐格式具有如下优越性

（1）无条件稳定，不需要考虑时间步长和空间步长的关系；

（2）收敛步数短。

另外，我们从以上数据还可以看出，对本算例而言，有如下特点：

显格式的计算结果相对于精确值偏小，而隐格式的计算结果偏大。

4.本文程序

在压缩包中，应该有.m文件的程序，fluid_conpution.m为显格式的计算程序，fluid_conpution2.m为隐格式的计算程序。

4.1.FTCS格式程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% FTCS格式计算热传导方程

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc

clear 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)^2

end

% 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

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'),ylabel('u/^{o}C')

legend('t=0时刻的温度分布','t=0.5s时刻的温度分布')

4.2.FTCS隐格式程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% FTCS隐格式计算热传导方程

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc

clear 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));

end

u=w;

end

% 4.绘制计算结果

plot(x,u,'k-')

title('隐格式下t=0.5s时刻的温度分布，时间步长500,000')

xlabel('x/m'),ylabel('u/^{o}C')

legend('t=0时刻的温度分布','t=0.5s时刻的温度分布')