流体力学中各参数的量纲

基本单位

质量:[latex]M[/latex]

长度:[latex]L[/latex]

时间:[latex]T[/latex]

温度:[latex]K[/latex]

衍生单位

密度:质量除以体积,所以就是[latex]ML^{-3}[/latex]

速度:单位是[latex]m/s[/latex],量纲就是[latex]LT^{-1}[/latex]

加速度:[latex]LT^{-2}[/latex]

力:众所周知F=ma,其中m是质量,a是加速度,所以就是[latex]MLT^{-2}[/latex]

动量:质量*速度。所以就是[latex]MLT^{-1}[/latex].

能量:做功是F*S,也就是力乘以距离,这个就是能量。所以能量的单位就是力乘以距离:[latex]ML^2T^{-2}[/latex].

功率:做功的速率,就是能量除以时间,所以就是[latex]ML^2T^{-3}[/latex].

比热容:让单位质量物体升温1K需要的能量。量纲就是能量/(质量*温度),也就是[latex]L^2T^{-2}K^{-1}[/latex]

压强:单位面积上受到的力。就是力除以面积。前面说了力的量纲是[latex]MLT^{-2}[/latex],面积则是[latex]L^{2}[/latex],两个一除就是[latex]ML^{-1}T^{-2}[/latex].

热体胀系数:温度增加1℃时密度的相对减小率(也是体积的相对膨胀率).所以就是[latex]-\frac{1}{\rho}\frac{\partial \rho}{\partial p}[/latex],其中[latex]\rho[/latex]的量纲约去,热膨胀系数的量纲就是[latex]K^{-1}[/latex].

等温压缩系数:压强增加一个单位时,流体密度的相对增加率。类似地,“相对增加率”无量纲,压强在分母上,所以就是[latex]M^{-1}LT^2[/latex]

杨氏模量:单位面积应力(压强)=模量*相对应变(无量纲),可以想象模量越大的物品越“硬”,越能把较小的应变换算成较大的力。所以杨氏模量的量纲就是压强的量纲[latex]ML^{-1}T^{-2}[/latex]

体积弹性模量:(其实就是等温压缩系数的倒数),这也是一种模量:压强变化=模量*相对变形率,量纲就是压强的量纲[latex]ML^{-1}T^{-2}[/latex].

动力粘度[latex]\mu[/latex]:在牛顿流体的剪切流动中,单位面积所受的粘性切应力和速度梯度成正比。速度梯度就是“速度差”,不是速度。所以它的单位就是力/(面积*速度梯度),速度梯度的量纲是[latex]LT^{-1}/L=T^{-1}[/latex],所以就变成了力/面积*时间,即压强*时间,所以它的单位是帕斯卡*秒。也就是[latex]ML^{-1}T^{-1}[/latex].

运动粘度[latex]\nu[/latex]:它的定义是[latex]\nu=\mu/\rho[/latex]也就是动力粘度/密度,没什么意义,存在目的主要是为了方便,例如测运动粘度不需要测量力。动力粘度的量纲是[latex]ML^{-1}T^{-1}[/latex],密度则是[latex]ML^{-3}[/latex],所以运动粘度的量纲就是[latex]L^2T^{-1}[/latex]

导热系数[latex]k[/latex]:根据傅里叶热传导定律:[latex]\textbf{q}=-k\nabla T[/latex],其中左边是热流密度。两边都应该是强度量,而显然时间越长,接触面积越大则导热越多,所以左边的热流密度就是导热能量/(时间*面积),其中能量是[latex]ML^2T^{-2}[/latex],时间*面积是[latex]TL^2[/latex],所以热流密度就是[latex]MT^{-3}[/latex].右边温度梯度是[latex]KL^{-1}[/latex],所以导热系数的量纲就是[latex]MLT^{-3}K^{-1}[/latex].或者用:左边热流密度就是瓦特/平方米,右边温度梯度就是开尔文/米,左除以右就是瓦特/(开尔文*米).

牛顿粘性定律、热传导都是一种输运过程,都具有输运流量=常数*势梯度的特点。

流体模拟若干名词解释

Semi Lagrangian Advection

“全拉格朗日”是:就去算一个个粒子,网格啥的一概没有。

“半拉格朗日”就“半”在,它每次用粒子方式Advect之后,就又回到了欧拉网格上。

具体地,第n+1步的A(i,j)等于第n步的A(x,y),其中(x,y)是(i,j)这个点在第n步时应该在的位置。

MacCormack Scheme

考虑这个方程

[latex]\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0[/latex]

MacCormack方法分成两步:Predictor Step和Corrector Step.

Predictor Step:

把方程离散化,就变成了

[latex]\frac{u_i^{\overline{n+1}}-u_i^n}{\Delta t}+a\frac{u_{i+1}^n-u_i^n}{\Delta x}=0[/latex]

这个[latex]u_i^{\overline{n+1}}[/latex]表示的是“预估的第n+1步的值”的意思。把它写成explicit形式:

[latex]u_i^{\overline{n+1}}=u_i^n-a\frac{\Delta t}{\Delta x}\left(u_{i+1}^n-u_i^n\right)[/latex]

Corrector Step:

假设我们知道时间n+0.5的信息,那么就应该有

[latex]\frac{u_i^{n+1}-u_i^{n+\frac{1}{2}}}{1/2\Delta t}+a\frac{u_i^{\overline{n+1}}-u_{i-1}^{\overline{n+1}}}{\Delta x}=0[/latex]

注意,Predictor Step里面,关于x偏导那一项是i+1减i(前向差分),但Corrector Step里面就是后向差分了。这两个也可以交换,先后向差分再前向差分。但是必须不一样。这是达到二阶精度所必须做的。

Corrector Step关于x偏导是用“预估值”[latex]u^{\overline{n+1}}[/latex]做的,需要注意。

我们不知道[latex]u_i^{n+\frac{1}{2}}[/latex]是什么,就用“预估值”和[latex]u^n[/latex]进行插值:

[latex]u_i^{n+\frac{1}{2}}=\frac{1}{2}\left( u_i^n+u_i^{\overline{n+1}} \right)[/latex]

代入进去就变成了

[latex]u_i^{n+1}=\frac{u_i^n+u_i^{\overline{n+1}}}{2}-a\frac{\Delta t}{2\Delta x}\left( u_i^{\overline{n+1}}-u_{i-1}^{\overline{n+1}}\right)[/latex]

BFECC Scheme

BFECC是缩写,全拼是Back and Forth Error Compensation and Correction.

假设有场[latex]\phi[/latex].

先用“正向”Advect的算子[latex]A[/latex]算出一个预估值[latex]\hat{\phi}^{n+1}=A(\phi^n)[/latex]

然后再用“反向”的算子[latex]A^R[/latex]算出第n步的“估计”值[latex]\hat{\phi}^n=A^R(\hat{\phi}^{n+1})[/latex]

然后我们估计得到一个误差[latex]e=(\hat{\phi}^n-\phi^n)/2.[/latex]

然后用这个误差修正初始值:[latex]\bar{\phi}^n=\phi^n-e[/latex]

并重新Advect:[latex]\phi^{n+1}=A(\bar{\phi}^n)[/latex]

Projection

Advect并不保不可压性。假设Advect得到流场[latex]\textbf{v}[/latex]。

根据亥姆霍兹分解,我们可以设想,求一个“压力”场[latex]q[/latex],使得[latex]\nabla^2q=\nabla\cdot\textbf{v}[/latex]

然后更新[latex]\textbf{v}_1=\textbf{v}-\nabla q[/latex]

我们设[latex]\textbf{v}=(u,v)[/latex].右边对[latex]\textbf{v}[/latex]求散度,由于拉普拉斯算子对称,所以这里也可以用中心差分:

[latex](\nabla\cdot\textbf{v})_{i,j}=(u_{i+1,j}-u_{i,j})/2\Delta x+(v_{i,j+1}-v_{i,j-1})/2\Delta y[/latex]

前面的拉普拉斯算子写出来就是

[latex](\nabla^2q)_{i,j}=(q_{i+1,j}-2q_{i,j}+q_{i-1,j})/\Delta x^2+(q_{i,j+1}-2q_{i,j}+q_{i,j-1})/\Delta y^2[/latex]

我们要解的是[latex]q[/latex].如果把它写成一个线性方程组的形式,那么对角线就是[latex]-2(1/\Delta x^2+1/\Delta y^2)[/latex],(i,j)那一行里面,(i+1,j)和(i-1,j)对应项的元素是[latex]1/\Delta x^2[/latex],(i,j-1)和(i,j+1)对应项都是[latex]1/\Delta y^2[/latex].方程右边的常数向量则是展开的[latex](\nabla\cdot\textbf{v})[/latex].

这里还要注意一个Boundary Condition的情况。假设我们的网格范围是(1,1)~(n-1,m-1),外面一圈是“假”的边界值。如果采取“复制型”的边界条件,就是例如q(0,j)=q(1,j)这样。这时候(1,j)这一行就少了一项。

Jos Stam Stable Fluids中一些问题解读

Method of Characterstics特征方法

Stable Fluids文章中使用特征方法解决Advection(对流)。

维基百科对特征方法的解释:

https://en.wikipedia.org/wiki/Method_of_characteristics

就是找一条让PDE成为ODE的线,沿着这条线解ODE,这样就搞定了原来的PDE问题。

看论文中的解释:

特征方法用来解这样的方程:

这个[latex]\frac{\partial a(\textbf{x},t)}{\partial t}+\textbf{v}(\textbf{x})\cdot\nabla a(\textbf{x},t)[/latex]其实就是随体导数。

然后论文定义[latex]\textbf{p}(\textbf{x}_0,t)[/latex]为速度场\textbf{v}的characterstics(特征),满足

这个[latex]\textbf{p}(\textbf{x}_0,t)[/latex]的意思就是在时刻0时位于[latex]\textbf{x}_0[/latex]的这个流体点在[latex]t[/latex]跑到了什么地方。

定义

为场a的“随体场”,也就是采取拉格朗日观点,定义它为0时刻在[latex]\textbf{x}_0[/latex]的流体点在t时刻的值。

上面说了,我们要解的方程相当于随体导数,那么[latex]\bar{a}[/latex]就是“随体场”,所以求解的直接就是它的时间导数

具体来说,在做Advection的时候,沿着流线(characterstic line)“上溯”一步即可。

边界条件

Jos Stam代码中的边界条件如下:

void set_bnd ( int N, int b, float * x )
{
	int i;

	for ( i=1 ; i<=N ; i++ ) {
		x[IX(0  ,i)] = b==1 ? -x[IX(1,i)] : x[IX(1,i)];
		x[IX(N+1,i)] = b==1 ? -x[IX(N,i)] : x[IX(N,i)];
		x[IX(i,0  )] = b==2 ? -x[IX(i,1)] : x[IX(i,1)];
		x[IX(i,N+1)] = b==2 ? -x[IX(i,N)] : x[IX(i,N)];
	}
	x[IX(0  ,0  )] = 0.5f*(x[IX(1,0  )]+x[IX(0  ,1)]);
	x[IX(0  ,N+1)] = 0.5f*(x[IX(1,N+1)]+x[IX(0  ,N)]);
	x[IX(N+1,0  )] = 0.5f*(x[IX(N,0  )]+x[IX(N+1,1)]);
	x[IX(N+1,N+1)] = 0.5f*(x[IX(N,N+1)]+x[IX(N+1,N)]);
}

Stable Fluids采用的边界条件是Dirichlet边界条件,也就是:在边界处,速度场的法向分量为零。即,没有任何通量“穿过”边界。

在文章中,速度场的离散点分布如下:

注意这里的坐标是1~N.想要实现“法向为零”的边界条件怎么办呢?我们取边界上的一个点,算该点处的速度值时,会取“外面”的速度点:

就像这样,红点的速度值是格子左上角黑点和格子外面蓝点的平均值。

那么,如果想让红点处速度的y分量为零,只需要把蓝点速度的y分量置为黑点速度的y分量的相反数即可。

蓝点的x分量呢?红点速度的x分量应当和黑点的x分量相同,那么蓝点的x分量也应该等于黑点的x分量。

最后,set_bnd函数的结尾四句是设置四个角的边界条件:

红点的速度是两个蓝点速度的平均值。原因是,例如我们考察y分量,令红点为y00:

y01=-y11

y10=y11

考虑到角上的速度应当为0,也就是说(y00+y01+y10+y11)/4=0,代入得到y00+y11=0,x的情况类似,也就是说红点的速度应该和黑点的速度恰好相反。

但是Jos Stam代码中,这个点的速度应该是0,原因有待实验。

Poisson Solver

在Stable Fluids里面,求解泊松方程的代码长这样:

void lin_solve ( int N, int b, float * x, float * x0, float a, float c )
{
	int i, j, k;

	for ( k=0 ; k<20 ; k++ ) {
		FOR_EACH_CELL
			x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]))/c;
		END_FOR
		set_bnd ( N, b, x );
	}
}

这是啥玩意呢?查了一圈发现,这玩意就是线性方程组的Jacobi Solver.看一眼维基百科:

https://en.wikipedia.org/wiki/Jacobi_method

注意,每次是“基于”对角元去更新解向量,这可以解释为什么更新的是x[i][j].

Jacobi Solver的条件是:对角元最大。这个怎么满足呢?把泊松方程离散化以后,就变成了类似拉普拉斯算子的玩意:

比方说一个大致长这样的算子,中间被加了好多次,可不就是最大了嘛。每个格子一个方程,自然对角线上的元素最大。

三种坐标下的流体连续方程

笛卡尔坐标形式的连续方程

《流体力学》(上)(周光埛等,第二版)的第3章给出了这个方程的推导,这里列一下最终结果:

[latex]\begin{align*}
\boxed{
\frac{\partial\rho}{\partial t}+\frac{\partial(\rho u)}{\partial x}\frac{\partial(\rho v)}{\partial y}+\frac{\partial(\rho w)}{\partial z}=0
}
\end{align*}[/latex]

柱坐标形式的连续方程

在柱坐标系下:

[latex]\begin{align*}
x&=r\cos(\theta)\\
y&=r\sin(\theta)\\
z&=z\\
\end{align*}[/latex]

取控制体[latex]\mathrm{d}r,\mathrm{d}\theta,\mathrm{d}z[/latex],

则体积的“比例系数”就是[latex]\left|\frac{\partial(x,y,z)}{\partial(r,\theta,z)}\right|=r[/latex](参考伍胜健《数学分析》(第三册)第166页)。

整体变化率

在控制体内质量的变化率是:

[latex]\frac{\partial\rho}{\partial t}r\mathrm{d}r\mathrm{d}\theta\mathrm{d}z[/latex]

然后我们分别计算质量从三对相对的面流出控制体的速率。

沿r方向

沿这个方向,两个相对面的面积是[latex]r\mathrm{d}\theta\cdot\mathrm{d}z[/latex],那么(沿r方向做泰勒一阶近似),质量流出的速率就是

[latex]\begin{align*}\frac{\partial(\rho v_r\cdot r\mathrm{d}\theta\mathrm{d}z)}{\partial r}\mathrm{d}r\end{align*}[/latex]

注意到[latex]\mathrm{d}r,\mathrm{d}\theta,\mathrm{d}z[/latex]对r取偏导的结果均为0,所以这个值等于

[latex]\begin{align*}
\frac{\partial(\rho rv_r)}{\partial r}\mathrm{d}r\mathrm{d}\theta\mathrm{d}z
\end{align*}[/latex]

沿[latex]\theta[/latex]方向

沿这个方向,两个相对面的面积是[latex]\mathrm{d}r\cdot\mathrm{d}z[/latex],所以质量流出的速率是:

[latex]\begin{align*}
\frac{\partial(\rho v_{\theta}\cdot\mathrm{d}r\mathrm{d}z)}{\partial\theta}\cdot\mathrm{d}\theta
=\frac{\partial(\rho v_{\theta})}{\partial\theta}\mathrm{d}r\mathrm{d}\theta\mathrm{d}z
\end{align*}[/latex]

沿z方向

沿这个方向,两个相对面的面积是[latex]r\mathrm{d}\theta\cdot\mathrm{d} r[/latex]

所以质量流出的速率是:

[latex]\begin{align*}
\frac{\partial(\rho v_z\cdot r\mathrm{d}\theta\mathrm{d}r)}{\partial z}\mathrm{d}z
\end{align*}[/latex]

由于是对z求偏导,所以r也应该拿到外面:

[latex]\begin{align*}
\frac{\partial(\rho v_z)}{\partial z}r\mathrm{d}r\mathrm{d}\theta\mathrm{d}z
\end{align*}[/latex]

把这四个式子联立起来:

[latex]\begin{align*}
-\frac{\partial\rho}{\partial t}r\mathrm{d}r\mathrm{d}\theta\mathrm{d}z=
&\frac{\partial(\rho rv_r)}{\partial r}\mathrm{d}r\mathrm{d}\theta\mathrm{d}z\\
+&\frac{\partial(\rho v_{\theta})}{\partial\theta}\mathrm{d}r\mathrm{d}\theta\mathrm{d}z\\
+&\frac{\partial(\rho v_z)}{\partial z}r\mathrm{d}r\mathrm{d}\theta\mathrm{d}z
\end{align*}[/latex]

把微分项和r除掉,就得到了柱坐标形式的连续方程:

[latex]\begin{align*}\boxed{
\frac{\partial\rho}{\partial t}+\frac{1}{r}\frac{\partial(\rho rv_r)}{\partial r}+\frac{1}{r}\frac{\partial(\rho v_{\theta})}{\partial\theta}+\frac{\partial(\rho v_z)}{\partial z}=0
}\end{align*}[/latex]

球坐标形式的连续方程

在球坐标系下:

[latex]\begin{align*}
x=&r\sin\phi\cos\theta\\
r=&r\sin\phi\sin\theta\\
z=&r\cos\phi
\end{align*}[/latex]

取控制体[latex]\mathrm{d}r,\mathrm{d}\phi,\mathrm{d}\theta[/latex],体积的“比例系数”为[latex]\left|\frac{\partial(x,y,z)}{\partial(r,\phi,\theta)}\right|=r^2\sin\phi[/latex](伍胜健《数学分析》(第三册)第167页).

整体变化率

整体变化率为:

[latex]\begin{align*}
\frac{\partial\rho}{\partial t}r^2\sin\phi\mathrm{d}r\mathrm{d}\phi\mathrm{d}\theta
\end{align*}[/latex]

沿r方向

对于r固定的面,x,y,z由[latex]\phi,\theta[/latex]确定。我们计算出:

[latex]\begin{align*}
E=&(x_\phi^{\prime})^2+(y_\phi^\prime)^2+(z_\phi^\prime)^2\\
F=&x_\phi^\prime x_\theta^\prime+y_\phi^\prime y_\theta^\prime+z_\phi^\prime z_\theta^\prime\\
G=&(x_\theta^\prime)^2+(y_\theta^\prime)^2+(z_\theta^\prime)^2
\end{align*}[/latex]

那么这个面的面积就是

[latex]\begin{align*}
\sqrt{EG-F^2}\mathrm{d}\phi\mathrm{d}\theta
\end{align*}[/latex]

(伍胜健《数学分析》(第三册)第204页)

经计算得出

[latex]\begin{align*}
E=&r^2\\
F=&0\\
G=&r^2\sin^2\phi\\
\sqrt{EG-F^2}=&r^2\sin\phi
\end{align*}[/latex]

[latex]\phi[/latex]对r求导得0,所以流出速率就是

[latex]\begin{align*}
\frac{\partial(\rho r^2v_r)}{\partial r}\sin\phi\mathrm{d}r\mathrm{d}\phi\mathrm{d}\theta
\end{align*}[/latex]

沿[latex]\phi[/latex]方向

此时,x,y,z由r,[latex]\theta[/latex]决定。类似地:

[latex]\begin{align*}
E=&1\\
F=&0\\
G=&r^2\sin^2\phi\\
\sqrt{EG-F^2}=&r\sin\phi
\end{align*}[/latex]

流出速率为

[latex]\begin{align*}
\frac{\partial(\rho\sin\phi v_\phi)}{\partial\phi}r\mathrm{d}r\mathrm{d}\phi\mathrm{d}\theta
\end{align*}[/latex]

沿[latex]\theta[/latex]方向

此时,x,y,z由r,[latex]\phi[/latex]决定:

[latex]\begin{align*}
E=&1\\
F=&0\\
G=&r^2\\
\sqrt{EG-F^2}=&r
\end{align*}[/latex]

流出速率为

[latex]\begin{align*}
\frac{\partial(\rho v_\theta)}{\partial\theta}r\mathrm{d}r\mathrm{d}\phi\mathrm{d}\theta
\end{align*}[/latex]

合起来,把微分项除掉,就是:

[latex]\begin{align*}\boxed{
\frac{\partial\rho}{\partial t}+\frac{1}{r^2}\frac{\partial(\rho r^2 v_r)}{\partial r}+\frac{1}{r\sin\phi}\frac{\partial(\rho\sin\phi v_\phi)}{\partial\phi}+\frac{1}{r\sin\phi}\frac{\partial(\rho v_\theta)}{\partial\theta}=0
}\end{align*}[/latex]

基于MacCormack方法的Quasi-1D Isentropic Nozzle Flow数值模拟

物理模型

这篇文章是我学习《COMPUTATIONAL FLUID DYNAMICS: The Basics with Applications》(作者John D. Anderson, Jr.)的心得体会。

这里的Quasi-1D Nozzle Flow是书上第七章的内容:用MacCormack算法求解拉瓦尔喷管(convergent-divergent nozzle)中的流场。这种喷管的示意图如下:

这种喷管的横截面是漏斗状,先变窄再变宽(这就是’convergent-divergent’的由来),它可以把亚音速气流(从左端进入)加速到超音速(从右端流出),在最窄处恰好是音速。原因是亚音速气流和超音速气流的性质不同:亚音速流截面越小流速越快,超音速流则截面越大流速越快。

在这个问题里面,我们求解二维的问题,并假定流场的性质只和x坐标有关(也就是说在上图中,一‘竖条’上的流体属性是相同的),这就是“Quasi-1D”的含义。当然,这种假设只是对现实的一个近似,目的是计算方便。

微分方程

既然是一维求解,就需要用到一维的方程。假定流场内是理想气体,故可以从方程中消去压强场,剩下三个未知场:密度、速度(x方向)、温度。书上有相关推导,最终的微分方程是(7.42)~(7.44):

在这一章节中,我们采用无量纲量进行计算,以消除物理单位。各个参数对应的无量纲量定义在书的第297页:

其中gamma是比热容(对于空气,取1.4),下标为0代表在喷管的入口处(最左边),A*是喷管最窄处的面积。

经过推导,可以得到对应无量纲量的微分方程,即(7.46),(7.48),(7.50).

MacCormack

当然,计算机程序无法直接求解微分方程,必须设法将其转化。在这里我们采用MacCormack方法。这是一种explicit方法。所谓explicit,就像(7.42)~(7.44)一样,每一个方程里恰好含有一个关于时间的导数,这样在递推时就可以显式(explicitly)计算出来每一个场关于时间的差分,然后递推。

但直接用差分乘以时间步的精度是一阶的(它对应了泰勒展式中的头两项),这就是用到MacCormack方法的原因:它可以在只处理一阶导数的情况下将explicit方程的求解精度提升至二阶。

MacCormack方法的步骤叫做predictor-corrector step:先按照上述的一阶方法求出来下一个时刻的“预测流场”(predictor),然后用“预测流场”算出来一个“预测时间差分”,再回过头,把“预测差分”和在本时间步求得的时间差分做平均,用这个平均值去求下一个时间步的流场,作为最终结果(corrector).有一个关键点:在第一次用现有流场算时间差分时需要用forward difference(就是(f[i+1]-f[i])/dx),第二次用预测流场算时间差分时则用rearward difference(就是(f[i]-f[i-1])/dx),或者第一次用rearward,第二次用forward.

MacCormack方法是《The Effect of Viscosity in Hypervelocity Impact Cratering》,Robert W. MacCormack提出的,该论文里称用傅里叶级数进行分析可以证明,它的精度是二阶。

explicit方程本身在书中有推导.用forward difference的(predictor step)是(7.51)~(7.53):

用rearward difference的(corrector step)是(7.57)~(7.59):

变量上面加横线代表这个是“预测流场”。

时间步长

对于explicit方法,第四章的分析表示,要想获得稳定的解,时间、空间步长必须满足关系:

C是Courant number(这个理论叫Courant–Friedrichs–Lewy condition),它不应该超过1,在本例中我们取0.5.注意到当地音速a和流体速度V是随着流场改变的,对于本例,我们取计算出的所有dt的最小值。

算例1:亚音速-超音速

边界条件

书中称,由于喷管最窄处左边是亚音速流场,也就是说扰动会向左传播,直至喷管的入口,所以只能固定喷管入口处的两个变量(密度和温度),但允许速度值自由“浮动”。严谨的数学分析“超越了本书的范围”。

对于入口处的速度值和喷管出口处的所有值,由于forward difference和rearward difference总有一个在该处无法计算,所以它们的值无法直接从差分方程中获得。对于这些地方,我们运用linear extrapolation(假设每个值随x坐标线性):

 

我们规定入口处的密度值和温度值是固定的常数。由于采用了无量纲量,它们就固定为1:

参数设定

在这个算例中,我们规定喷管的长度是3L,截面函数是:

它在1.5L的地方达到最窄.

x方向的步长是0.1,这意味着一共31个格点。

在t=0的时候(也就是初始值),密度、温度、速度场如下:

实验结果

代码在此:

https://github.com/wmdcstdio/DigitalWindTunnel/blob/master/NozzleFlowQuasi1D/NozzleFlowQuasi1D/main.cpp

这份代码会先输出亚音速-超音速算例的结果。

书上的Table 7.1是初始的各个参数值,Table 7.3是进行1400次迭代后的各个参数值。可以发现,进行1400次迭代后,马赫数在格点16和17(最中间的两个格点)之间跨越了1,说明在最窄处流速确实恰为音速。

我的代码计算出的数值和书上结果完全一样,除了几处printf的舍入误差。限于篇幅(懒),我就不贴具体数据了,有兴趣者请自己下载代码运行:)

算例2:纯亚音速

在这个算例中,管子的形状稍有不同:

纯亚音速的边界情况和上面讨论的亚音速-超音速略有不同。原因是解不再唯一,实际上喷管出入口的压强比值和流场方程的解一一对应。因此我们需要规定这个比值,根据标准气体方程,实际上就是规定出口处[latex]\rho’T'[/latex]的值(为一个常数,在这个算例中等于0.93)。我们可以extrapolate出其中的一个(比如[latex]\rho[/latex]),然后用[latex]\frac{常数}{\rho}[/latex]求出另一个。

这个算例的初始值也不同:

其他和亚音速-超音速的算例相同,包括入口处的extrapolation等。

代码还是这份代码:

https://github.com/wmdcstdio/DigitalWindTunnel/blob/master/NozzleFlowQuasi1D/NozzleFlowQuasi1D/main.cpp

它应该在输出第一个算例之后,输出第二个算例迭代5000步后的结果,也就是书中的table 7.7,请自行验证:)