流体力学第10章

⒃ ⒄
退出
特征方程⒃⒄的差分形式为：
C ：x j
xw
(u c)W ,P t,uP
uW
2
(C 1
P
CW ) GW ,P t
0
C ：x j
xE
(u
c)E,P t,uP
uE
2
1(CP
CE ) GE,Pt
0
xW 和xE 的位置未知，要反复迭代才能确定其位置，当 xW , xE
确定后，点的参数用节点 x x , x j1, j j1 的参数内插得到，符号表
第十章 工程流体力学的计算方法（CFD基础）
§6.1代数方程的牛顿迭代法
牛顿迭代法用于求解超越方程 f (x) 0 的根，在曲线
y f (x) 上取一点 (x0 , y0 )
求
x1
x0
f (x0 ) f (x0 )
显然 x1 是方程 f (x) 0 的一个比 x0
更精确的解，
重复以上计算可以得到任意精确的解。
目 的 是 要 计 算 出 管 道 每 一 截 面 处 的 流 动 参 数 从 时 刻 t=0 到
t=T任一时间里的变化值。
将管长L分为N等分，x
L N
即无量纲的空间步长
x
1 N
时间步长与空间步长相等 t x 1
N
退出
计算x–t图上的网络如上图所示。
如果 ti 时刻管轴上每一个节点 xj 的流动参数u，h已经算出， 利用特征线方程就可以计算下一时刻 ti1 各个节点 x j 的流动参数
hi1,0 1
u i 1, 0
1
B Fr
(hi
1,
0
B F ur i1,0 )
在x=1处，即j=N处，只有一条特征线C+但u 可以由式⒂给出：
hi1,N A F ur i1,N
ui1,N (t ) hi1,N
由以上两式联立可解出 ui1,N 和 hi1,N
逐个时刻进行计算，就可以得到各时刻管道上个节点的流动 参数，u和h的值。
t
t
x t dt x dt
如果：u c x
u 则：
(u
c)
u
u
dx
u
du
t
t
x t dt x dt
退出
u c x t
表示一条x–t的关系曲线，记作C+，称为特征曲线C+。
同样
u c x t
称为特征线C－ 。
这样⑥⑦式可以写成：
沿特征线C+：
dx dt du
u dP
C
u u dt 0
在阀门处，V和h的关系：V 2gh
A(t) A(t)
此时管内流速u： t 0, x L,u A0 V A0
2gh (t)
2 gh
取L为特征长度，h0为特征水头，u0取为特征速度，L0/C特征时间
将式⑧⑨⑩⑾⑿写成无量纲的形式：又 P gh
沿C+：dx 1, du dh ku u dt 0 ⒀
t
C 2D
⑾
设管长为L，管壁为D，截面
A0
D2 4
管内流速为u，阀门过流截面A随时间而变
A A(t) ，过流速为V，V>u
退出
在初始时刻（t=0），管道定常出流，根据伯努利方程：
h0
V2 2g
L D
u2 2g
u 2 [
2g
L D
( A0 )2 ] A
此时，管流速度为定常出流速度，记作u0
u0
退出
例：水从池中经管道流出，已知管长l 50m 沿程阻力损失系数
0.03 局部阻力损失系数 6 ,水径 H 2m
设计流量 Q 0.06 m3 h 试求管径d
解：列水面和管道出口截面的伯努利方程：
H (a l ) 1 ( 4Q )2 d 2g d 2
代入数据化简得：
6721d 5 7d 1.5 0 令 x 10d 则上式化为：f (x) 0.06721x5 0.7x 1.5 0 选 x0 2 作为初值 x1 x0 f (x0 ) f (x0 )
i1, j )
1 y 2
( i, j1
2 i, j
) i, j1
0
令 x y
则： i1, j
2 i, j
2(1 2 ) i1, j
2 i, j 1
i1, j
0
⑦
退出
对每个网格节点都建立形如上式的差分方程，就得到各 节点的流函数的代数方程组，给出边界条件，用迭代法可求 出其数值解。 例：水在两平板间流动，上板壁的渗透速度v0=1m/s下壁不可 渗透，入口和出口速度均匀分布，分别为u1=3m/s和u2=1m/s 和，设板长h=3m 宽h=1.5m
dt
Fr
沿C－：dx 1, du dh ku u dt 0 ⒁
dt
Fr
式中：沿 k
L0
2CDu0
,
Fr
Ch0 gu0
（边界条件：）
退出
初始条件：t
0：u
1,
h
1
L
D
D 1
2
(t)
X
t 0, x 0：h 1
边界条件： x 1：u
L
D
1
2 (0)
(t
)
h (t)
h
⒂
（初始条件：）
⑧
C 2D
dx 沿特征线C－： dt
du
u
dP
C
C
2D
u
u
dt
0
⑨
退出
下面用特征线法研究水击波和隧道压力波问题
1.水击波 水击波传播速度C=1200～1400m/s，远大于水速度u，因此
，特征线方程可以写成：
沿C+： x C, du dP u u dt 0
t
C 2D
⑽
沿C－：x C, du dP u u dt 0
y2
如果将坐标原点放在x1处，即x1=0且令x2–x1=L
则有：
y
(1
x L
)
y1
x L
y2 ,
L
xi
x1
令
1
(
x)
1
x L
2
(x)
x L
退出
∴ y y11(x) y22 (x)
显然：
i
(x)
x L
1
L x 0 x 0 x L
i i
(xi ) 1 (x j, j i)
压力波的传播，都属于这种流动。
圆截面中可压缩粘性的非定常流动的运动微分方程：
u u u 1 P 4 0
①
t x x D
管壁粘性摩擦应力 τ0 可以用沿程阻力损失系数 表示：
0
1 u 2
8
1 u u
8
②
②代入①得：u u u 1 P u u 0 t x x 2D
③
退出
连续性方程： t
这样就可以求出P点的h和u，即：
hP
hi1, j
A 2
B ,uP
ui1, j
AB 2Fr
其中： A hi, j1 Fr ui, j1 (1 k ui, j1 t ) B hi, j1 Fr ui, j1 (1 k ui, j1 t )
退出
在x=0处，只有一条特征线C－，但由于压强水头恒为1，因此：
将时刻 ti1 的节点（i+1，j）记为P，过P作特征线C+和C－必
通过节点和（i，j-1）（i-1，j）和（i，j-1）分别记为W和E
退出
根据特征线方程⒀⒁有：
C ：uP
uW
1 Fr
(hP
hW ) kuW
uW
t
0
C ：uP
uE
1 Fr
(hP
hE ) kuE
uE
t
0
或 C ：hP A Fr uP C ：hP B Fr uP
设函数在这三个节点的值为：
fi1 f (xi1 ), fi f (xi ), fi1 f (xi1 ) 设节点间距为 x 则有泰勒展开式
f i 1
fi
f
i
x
1 2!
fi x2
o(x3 )
①
f i 1
fi
fix
1 2!
fi x
2
o(x3 )
②
退出
则有
fi
fi
fi1 x
o(x)
③ 一阶导数向后差分式
示两点间的平均值。 当各节点的u和c求出后，由等熵关系求出节点的ρ ，R，T 其初始条件和边界条件的确定比较复杂。
退出
§6.4有限元的插值函数
一．线性插值
如果已知曲线上的几个点坐标
(x1, y1 ),, (xn , yn ) 求曲线
方程，可以近似用折线表示：
y
(1
x x2
x1 x1
)
y1
x x1 x2 x1
2dc
1
退出
于是特征线方程可以写成：
沿C+： dx u c, du 2dC Gdt 0
dt
1
沿C－： dx u c, du 2dC Gdt 0
dt
1
式中：G u u
2D
特征线不经过节点，因为波动 程的差分要满足：x u c
t 因空间步长△x要大于时间步长△t 所以w点落在节点（i,j-1）的右侧 E点落在节点（i,j+1）的左侧
点（3，2） 2,2 10 3,2 0.5 4 3,3 0
点（3，3） 2,3 4 3,2 104 3,3 110 0
上面4方程可用矩阵表示：
退出
利用高斯法解此线性方程组得： 2,2
7 6
,
2,3
7 3
,
3,
2
5 6
,
3,
3
5 3
于是各节点的流速函数的数值为：
4.50 3.00 1.50 0.00
退出
2.隧道压力波
波速C是小扰动波的传播，C 2 RT 气流速度u和C相比不一定
是微量，不能忽略，另外C是变化的，小扰动波的传播可视为等熵 过程。
由
C 2 RT和P RT和 P
常数
得到：
2dC dT CT
dP d dT P T
dP d P
dP P
2 1
dC C
dP
C
C
dP P
经3次迭代后得 x3 2.31707 误差小于 106
因此取 d 0.214m
退出
§6.2差分法
解析函数 f (x) 可以在点 x0 领域展开成泰勒级数
f (x)
f (x0 )
f
( x0
)( x
x0
)
1 2!
f
( x
x0
)2
设有 i 1,i,i 1 三个差分节点， 其坐标为 xi1 , xi , xi1
2gh0
L ( A0 )2
DA
P
距管道入口x的截面上压强水头 h (g) 的伯努利方程：
h0
h
x D
u2 2g
u2 , 2g
L
h0 ( A0
)2
DA
x
即：压强水头h沿x的分布为：h0
h[1
L
D] ( A0 )2
DA
退出
当t>0时，管道入口的压强水头恒为h0即 t 0, x 0, h h0
(u)
x
0或
t
u
x
u x
0
④
压力波的传播速度：C 2 dP
d
音波（微压波） C 2 RT
E
水击波：
C2
1
ED
E因
D：管道直径
E：流体体积弹性系数
E因：管壁材料的弹性模量 ρ：流体密度
δ：管壁厚度
水击波的传播速度C=1200～1400m/s
退出
P 1 P
t P t C 2 t
P 1 P
具有三阶精度。
退出
在平面势流中，流函数和速度势函数均满足拉普拉斯方程：
2 2 0
x2 y 2
现将计算区域分成若干网格，每个
网格的边长都是 (x, y) ，节点
(xi , yi ) 简记为 (i, j)
其二阶导数可以用式⑥近似表示，则拉普拉斯的差分式为：
1 x 2
( i1, j
2 i, j
3.50 2.33 1.17 0.00
2.50 1.67 0.83 0.00
1.50 1.00 0.50 0.00
※ 流函数的物理意义：平面流动中，流通两条流线间任一 曲线的体积流量（单位厚度）等于两条函数之差。
来自百度文库
退出
§6.3特征线法 特征线法用于求解一无非定常可压缩流动问题的数值解
，水击压力波在管道内的传播，高速到与进入隧道时所产生的
fi
fi1 x
fi
o(x)
④ 一阶导数向前差分式
可见 f 具有 x 的一阶精度
①②上述两式相减则有：
f
i
fi1 x
fi
o(x2 )
⑤
①②上述两式相加则有：
fi
fi1 2 fi x 2
fi1 o(x2 )
⑥
对于形如 y f (x, y) 的微分方程也可以求出y的泰勒展开式
yi1
yi
yix
4,
j
y
1
4,2 0.5. 4,3 1.0, 4,4 1.5
④在上壁面，y h 1.5
v0
x
1
2,4 3.5, 3,4 2.5
退出
各节点的代数方程，(
x y
2)
由⑦式得：
点（2，2）1.5 10 2,2 3,2 4 2,3 0
点（2，3）3 4 2,2 10 2,3 3,3 14 0
1 2
yix 2
o(x3 )
yi
yi1
yi1x
1 2
yi1x 2
o(x3 )
退出
两式相减得：
yi
yi1
1 2
(
yi1
yi)x
1 2
(
yi1
yi)x 2
o(x)
yi1 yi yi x o(x4 )
yi
yi1
1 2
(
yi1
yi)x
1 4
yi x
3
o(x4
可见：yi
yi1
1 2
(
yi1
y)i x
将长和宽分成3等份：x 1m, y 0.5m
退出
①下壁面是一条流线，取其流函数为零，即 i,1 0,i 1,2,3,4
②左边入口处：
u1
u y
3
因而
1,2
1,1
1,1y
1
2
1,1y
2
o(y3 )
0 3 0.5 0
1
同理 1,3 3, 1,4 4.5
③右边出口处：u2
x P x C 2 x
这样连续性方程可改写成:
C U X
1 (P
C t
u P) x
0
③，⑤两式相加，减得：
u (u c) u 1 [P (u c) p] u u 0
t
x C t
x 2D
u (u c) u 1 [P (u c) p] u u 0
t
x C t
x 2D
如果：u c x 则：u (u c) u u dx u du