工程流体力学的计算方法CFD基础

第13章

工程流体力学的计算方法

（CFD基础）第13章

工程流体力学的计算方法

（CFD基础）§13.1代数方程的牛顿迭代法

牛顿迭代法用于求解超越方程

的根，在曲线

上取一点求显然是方程的一个比更精确的解，重复以上计算可以得到任意精确的解。§13.1代数方程的牛顿迭代法牛顿迭代法用于求解超越例：水从池中经管道流出，已知管长

沿程阻力损失系数局部阻力损失系数,水径设计流量试求管径d解：列水面和管道出口截面的伯努利方程：代入数据化简得：令

则上式化为：

选作为初值

经3次迭代后得

误差小于

因此取

退出例：水从池中经管道流出，已知管长沿程阻力损失系数局部阻力§13.2差分法解析函数

可以在点

领域展开成泰勒级数

设有

三个差分节点，

其坐标为

设函数在这三个节点的值为：

设节点间距为

则有泰勒展开式①

②

退出§13.2差分法解析函数可以在点领域展开成泰勒级数设有则有

③

一阶导数向后差分式④

一阶导数向前差分式可见

具有

的一阶精度①②上述两式相减则有：

⑤

①②上述两式相加则有：

⑥

对于形如

的微分方程也可以求出y的泰勒展开式退出则有③一阶导数向后差分式④一阶导数向前差分式可见两式相减得：可见：

具有三阶精度。退出两式相减得：可见：具有三阶精度。退出在平面势流中，流函数和速度势函数均满足拉普拉斯方程：现将计算区域分成若干网格，每个网格的边长都是

，节点

简记为其二阶导数可以用式⑥近似表示，则拉普拉斯的差分式为：

令

则：

⑦

在平面势流中，流函数和速度势现将计算区域分成若干网格，每个网

对每个网格节点都建立形如上式的差分方程，就得到各节点的流函数的代数方程组，给出边界条件，用迭代法可求出其数值解。例：水在两平板间流动，上板壁的渗透速度v0=1m/s下壁不可渗透，入口和出口速度均匀分布，分别为u1=3m/s和u2=1m/s和，设板长h=3m宽h=1.5m将长和宽分成3等份：

退出对每个网格节点都建立形如上式的差分方程，就得①下壁面是一条流线，取其流函数为零，即

②左边入口处：

因而

同理

③右边出口处：

④在上壁面，

退出①下壁面是一条流线，取其流函数为零，即②左边入口处：因而各节点的代数方程，

由⑦式得：

点（2，2）

点（2，3）

点（3，2）

点（3，3）

退出上面4方程可用矩阵表示：

-104104-100110-104014-10各节点的代数方程，由⑦式得：点（2，2）点（2，3）利用高斯法解此线性方程组得：

于是各节点的流函数的数值为：

4.503.502.501.503.002.331.671.001.501.170.830.500.000.000.000.00※流函数的物理意义：平面流动中，流过两条流线间任一曲线的体积流量（单位厚度）等于两条流函数之差。退出利用高斯法解此线性方程组得：于是各节点的流函数的数值为：§13.3特征线法

特征线法用于求解一维非定常可压缩流动问题的数值解，水击压力波在管道内的传播，高速列车进入隧道时所产生的压力波的传播，都属于这种流动。

圆截面中可压缩粘性的非定常流动的运动微分方程：①

管壁粘性摩擦应力

可以用沿程阻力损失系数

表示：

②

②代入①得：

③

退出§13.3特征线法特征线法用于求解一维非定常可压连续性方程：

④

压力波的传播速度：

音波（微压波）

水击波：

D：管道直径E：流体体积弹性系数E固：管壁材料的弹性模量δ：管壁厚度ρ：流体密度水击波的传播速度C=1200～1400m/s

退出连续性方程：④压力波的传播速度：音波（微压波）水击这样连续性方程可改写成:

③，⑤两式相加，减得：

如果：

则：

如果：

则：

退出这样连续性方程可改写成:③，⑤两式相加，减得：如果：表示一条x–t的关系曲线，记作C+，称为特征曲线C+。同样

称为特征线C－。这样⑥⑦式可以写成：

沿特征线C+：⑧

沿特征线C－：⑨

退出表示一条x–t的关系曲线，记作C+，称为特征曲线C+。同样下面用特征线法研究水击波和隧道压力波问题

1.水击波

水击波传播速度C=1200～1400m/s，远大于水速度u，因此，特征线方程可以写成：

沿C+：⑽

沿C－：⑾

设管长为L，管壁为D，截面

管内流速为u，阀门过流截面A随时间而变，过流速为V，V>u

退出下面用特征线法研究水击波和隧道压力波问题1.水击波在初始时刻（t=0），管道定常出流，根据伯努利方程：

此时，管流速度为定常出流速度，记作u0

距管道入口x的截面上压强水头

的伯努利方程：

即：压强水头h沿x的分布为：

退出在初始时刻（t=0），管道定常出流，根据伯努利方程：此时，当t>0时，管道入口的压强水头恒为h0即在阀门处，V和h的关系：

此时管内流速u：

取L为特征长度，h0为特征水头，u0取为特征速度，L0/C特征时间将式⑧⑨⑩⑾⑿写成无量纲的形式：又

沿C+：⒀

沿C－：⒁

式中：沿（边界条件：）

退出当t>0时，管道入口的压强水头恒为h0即在阀门处，V和h的关初始条件：

边界条件：

⒂

（初始条件：）

目的是要计算出管道每一截面处的流动参数从时刻t=0到t=T任一时间里的变化值。

将管长L分为N等分，

即无量纲的空间步长

时间步长与空间步长相等

退出初始条件：边界条件：⒂（初始条件：）计算x–t图上的网络如上图所示。

如果

时刻管轴上每一个节点

的流动参数u，h已经算出，利用特征线方程就可以计算下一时刻

各个节点

的流动参数

将时刻

的节点（i+1，j）记为P，过P作特征线C+和C－必通过节点和（i，j-1）（i-1，j）和（i，j-1）分别记为W和E退出计算x–t图上的网络如上图所示。如果时刻管轴上每一个节点根据特征线方程⒀⒁有：

或

这样就可以求出P点的h和u，即：

其中：

退出根据特征线方程⒀⒁有：或这样就可以求出P点的h和u，即：在x=0处，只有一条特征线C－，但由于压强水头恒为1，因此：

在x=1处，即j=N处，只有一条特征线C+但u可以由式⒂给出：

由以上两式联立可解出

和

逐个时刻进行计算，就可以得到各时刻管道上个节点的流动参数，u和h的值。

退出在x=0处，只有一条特征线C－，但由于压强水头恒为1，因此：2.隧道压力波

波速C是小扰动波的传播，

气流速度u和C相比不一定是微量，不能忽略，另外C是变化的，小扰动波的传播可视为等熵过程。

由

得到：

退出2.隧道压力波波速C是小扰动波的传播，气流速度u和C相比于是特征线方程可以写成：

沿C+：⒃

沿C－：⒄式中：

特征线不经过节点，因为波动方程的差分要满足：

因空间步长△x要大于时间步长△t

所以w点落在节点（i,j-1）的右侧

E点落在节点（i,j+1）的左侧

退出于是特征线方程可以写成：沿C+：⒃沿C－：⒄特征方程⒃⒄的差分形式为：

的位置未知，要反复迭代才能确定其位置，当

确定后，点的参数用节点

的参数内插得到，符号表示两点间的平均值。当各节点的u和c求出后，由等熵关系求出节点的ρ，R，T

其初始条件和边界条件的确定比较复杂。

退出特征方程⒃⒄的差分形式为：的位置未知，要反复迭代才能确定其§13.4有限元的插值函数一．线性插值如果已知曲线上的几个点坐标

求曲线方程，可以近似用折线表示：如果将坐标原点放在x1处，即x1=0且令x2–x1=L则有：

令

退出§13.4有限元的插值函数一．线性插值如果已知曲线上的几∴显然：

称为样条函数或插值函数

已知n个点曲线可以表示为：

一.加权余量法，设有流动问题的方程和边界条件是：

退出∴显然：称为样条函数或插值函数已知n个点曲线可以表一般化为：

设u近似为：

定义余量

为：

要使

最小