FLUENT中Courant Number的含义?

  • A+
所属分类:学霸天地

FLUENT中,经常会遇到Courant数(库朗数,Courant Number)这样一个概念。例如,在基于密度的求解器的隐式算法中(图1),在基于密度的求解器的显式算法中(图2),在基于压力的求解器的耦合算法中(图3),都有Courant数的设置。

 

FLUENT中Courant Number的含义?

1 基于密度的求解器的隐式算法中的Courant

FLUENT中Courant Number的含义?

2 基于密度的求解器的显式算法中的Courant

FLUENT中Courant Number的含义?

3 基于压力的求解器的耦合算法中的Courant

Courant数的含义是什么呢?在计算流体力学中,Courant数是由波传播速度a、网格尺寸Δx以及时间步长Δt定义的一个无量纲数

FLUENT中Courant Number的含义?

它代表的是在每个时间步内,流体中传播的波通过了多少个网格尺寸

Courant数对于计算过程的数值稳定性有着至关重要的影响。描述流体运动的Navier-Stokes方程组比较复杂,这里我们可以通过一个简化的例子来理解Courant数。考虑输运方程

FLUENT中Courant Number的含义?

其中a是波传播速度。这个方程的解非常简单,即波形沿着x轴以速度a移动。例如,在初始条件

FLUENT中Courant Number的含义?

下,并假定波速a=1,方程(1)的解如视频1所示。

视频1 输运方程(1)的解

现在我们采用有限差分法来求解方程(1)。在空间和时间建立离散网格(图4),

 

FLUENT中Courant Number的含义?

 

4 在空间和时间上离散

我们依次采用两种不同的离散格式来求解。首先,我们采用最简单的,即时间方向为一阶显式,空间方向是一阶向后差分

FLUENT中Courant Number的含义?

之所以叫做显式,是因为在下一个时间层(带有上标n+1的)上,不同空间位置(不同的下标)的u的值是可以各自独立地求出来的。

使用MATLAB编写的程序代码见附录1。计算结果如视频2所示。这里我们采用Δx=0.1Courant数则采用了三种,分别为0.50.8以及1.2,对应于从小到大的三种不同的时间步长。

视频2 使用格式(2)来求解方程(1)的结果。

可以看出,Courant数取0.50.8的时候,算法是数值稳定的(尽管两者和精确解的差别不同,即精度不同),而Courant数取1.2的时候,计算结果就发散了。实际上可以证明,对于离散格式(2),其临界稳定的Courant数等于1

接着我们换一种离散格式。我们将时间方向换成一阶隐式的

FLUENT中Courant Number的含义?

格式(3)与格式(2)的不同在于,在下一个时间层(带有上标n+1的)上,不同空间位置(不同的下标)的u的值是相互联系的,这就是它被称作隐式格式的原因。使用MATLAB编写的程序代码见附录2,由于下一个时间层上不同空间位置的u的值相互联系,所以采用fsolve来联立求解。计算结果如视频3所示。这里我们仍采用Δx=0.1Courant数则分别使用0.51.5以及3

视频3 使用格式(3)来求解方程(1)的结果。

可以看出,三种Courant数下算法都是稳定的。实际上可以证明,对于离散格式(3),在任何Courant数下都是稳定的,也就是说临界Courant数是无穷大。

实际上,每种计算格式都有自己的Courant数稳定范围。像上面的第一种格式,其数值稳定所对应的Courant数范围是[0,1],而第二种格式则为[0,+),如果再换一种格式,其数值稳定所对应的Courant数范围又可能会不一样。换句话说,要使算法稳定Courant数必须满足一定的条件,即Courant数必须在特定的范围之内,具体范围取决于离散格式。在计算流体力学中,这个条件被称为CFLCourant – Friedrichs – Lewy)条件,因为它最早是由Richard CourantKurt Friedrichs以及Hans Lewy1928年合作发表的一篇论文([1])中提出的。

为什么Courant数与数值稳定性有密切的联系呢?要理解这个问题,必须理解解析域和数值域的概念。一方面,从物理上来说,由于波是以有限速度传播的,所以对于n+1时间层上某个节点(例如节点i)来说,它的u值(uin+1)只依赖于某个特定区域内的u的数值,这个区域叫做解析域analytical domain)。例如,图5绿色粗实线就是解析域,显然其斜率为1/a,因为波的传播速度是a另一方面,从数值计算来说,uin+1值依赖于周围某些特定节点的u值,至于具体依赖哪些节点,就取决于离散格式了,这些被依赖的节点所构成的区域叫做数值域numerical domain)。例如,图5的离散格式是格式(2),即一阶显式格式,对于uin+1来说,它依赖于uinui-1n,所以数值域如图中黑色粗实线围成的区域所示。

有了解析域和数值域的概念,就容易理解Courant数与数值稳定性的密切联系了。第一Courant数的大小决定了解析域和数值域之间的位置关系。例如图5、6,其对应的离散格式是显式格式(2),这种离散格式的数值域在水平方向的宽度是一个网格尺寸,又由于Courant数反映了每个时间步流体中传播的波通过了多少个网格尺寸,所以,Courant数小于1意味着数值域包含了解析域,Courant数大于1则意味着数值域没有包含解析域。第二从物理上分析,要使数值计算稳定,必须要让数值域包含解析域这是因为,只有数值域包含解析域的格式才符合物理实际;如果数值域没有包含解析域,例如像图6那样,从物理上来说uin+1值应该受到解析域(绿色实线上)u值的影响,而数值计算格式却违背了这一点,所以是不符合物理实际的

FLUENT中Courant Number的含义?

5 格式(2),Courant<1,稳定

FLUENT中Courant Number的含义?

6 格式(2),Courant>1,不稳定

7是格式(3)(一阶隐式格式)的分析。由于是隐式格式,uin+1依赖于ui-1n+1以及uin。但是这还没完,根据离散格式,ui-1n+1又依赖于ui-2n+1以及ui-1n,这样一直下去,可以看出对于这种格式来说,其数值域是向左无限延伸的。所以,对于格式(3)来说,无论Courant数多大,数值域都始终包含解析域。因此,对于这种格式,无论Courant数多大都是稳定的,或者说数值稳定所对应的Courant数范围是[0,+)

 

FLUENT中Courant Number的含义?

7 格式(3),稳定

 

总而言之,Courant数决定了解析域和数值域之间的位置关系,从而影响计算的数值稳定性。每种离散格式有各自的Courant数稳定范围,通常显式格式的稳定范围小而隐式格式的稳定范围很大。当然,上面只是从稳定性进行讨论,其实Courant数还影响着计算效率。从计算效率上来说,Courant数取太小的话,会导致时间步长太小,计算所需的时间步数量太多,降低了计算效率。所以,在保证计算稳定的前提下,又应该取尽量大的Courant数

FLUENT中的Courant数的含义也是类似的,只不过FLUENT中求解的是流体运动的控制方程组(Navier-Stokes方程组),流体中存在着几种不同的波速(流速本身,以及流速分别与向正、负方向传播的声速的叠加),所以具体处理的时候复杂一些。对于显式格式(图2),正如上面所分析的那样,Courant数应当取比较小的值(默认值是1,必要的时候还应再减小一些);对于隐式格式(图1、图3),则Courant数可以取较大的值(理论上,如果方程是线性的,例如上面例子中的对流方程,可以取任意大的Courant数,但是由于实际流体的控制方程有非线性效应,所以实际上还是有上限的)。

顺便说一下一点经验,对于图1中基于密度的隐式算法的Courant数,其默认值是5,实际计算时,刚开始迭代可以调小一些,而快要收敛的时候可以增大(10以上,甚至几十)。

北航航空科学与工程学院的研究生李亮阅读了初稿并提出了很好的建议,在此表示感谢。



附录1:格式(2)的MATLAB代码

clear m=101; % 空间方向的网格节点数量 L=10; % 计算域的长度 a=1; % 波速 C=0.5; % Courant数 x=linspace(-L/2,L/2,m); % 网格节点的空间坐标 dx=x(2)-x(1); % 网格尺寸 u=zeros(1,m); u_new=u; % 设置初始条件 for i=1:m    if x(i)<0     u(i)=2;   else     u(i)=1;       end end t=0; dt=dx/a*C; % 时间步长 while t<4   % 求下一个时间层的解   for j=2:m     u_new(j)=-a*(u(j)-u(j-1))/dx*dt+u(j);   end   u_new(1)=2;   u=u_new;   t=t+dt;     % 绘图   plot(x,u);    axis([-L/2 L/2 0 3]);     drawnow end

附录2:格式(3)的MATLAB代码

function implicit m=101; % 空间方向的网格节点数量 L=10; % 计算域的长度 a=1; % 波速 C=2; % Courant数 x=linspace(-L/2,L/2,m); % 网格节点的空间坐标 dx=x(2)-x(1); % 网格尺寸 u=zeros(1,m); % 设置初始条件 for i=1:m    if x(i)<0     u(i)=2;   else     u(i)=1;       end end t=0; dt=dx/a*C; % 时间步长 while t<4   % 求下一个时间层的解   u_new=fsolve(@(ustar)f1(ustar,u,dt,dx,a,m),u);   u=u_new;   t=t+dt;     % 绘图   plot(x,u);    axis([-L/2 L/2 0 3]);     drawnow end function residual=f1(u_new,u,dt,dx,a,m) residual=zeros(1,m); residual(1)=u_new(1)-2; % 离散方程 for j=2:m   residual(j)=(u_new(j)-u(j))/dt+a*(u_new(j)-u_new(j-1))/dx; end

参考文献

[1] Courant, R.; Friedrichs, K.; Lewy, H. (1928), "Über diepartiellen Differenzengleichungen der mathematischen Physik",Mathematische Annalen (in German), 100 (1): 32–74

[2] J. D. Anderson Jr. Computational fluid dynamics: the basics with applications.McGraw-Hill, 1995

始发于微信公众号:流体那些事儿

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: