粘性流体力学—平面驻点流动(西门茨流动)资料
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
ode.m文件程序如下： function dfdx=ode(x, f) dfdx=[f(2);f(3);-f(1)*f(3)+f(2)^2-1];
lbc.m 文件程序如下： function res=lbc(f0,finf) res=[f0(1);f0(2);finf(2)-1];
程序说明： ode.m文件描述一阶常微分方程组，由于方程是三阶微分方程， 所以需要三个一阶微分方程来描述。 lbc.m文件描述边界条件，f0表示初始值，finf表示末端值，本问 题告诉我们的是f(1)和f(2)的初始值与f(2)的末端值，其中finf(2)-1表 示finf(2)-1=0，其它以此类推。
第二节 方程推导和求解
对于平面驻点流动，为了满足粘性流动的无滑移的条件，西门茨 给出其精确解.
现假定：
式中p0表示驻点O（x=0，y=0）处的压强；p为任意点（x，y） 处的压强；a为常数。
可证明，该假定将自动满足连续性方程。
由平面运动的N-S方程可以确定f与F两个函数，将上式带入恒定 的N-S方程
第一步：将方程化为一阶常微分方程组。
f (1)
df (1) f (2) df (2) f (3) df (3) f (1) f (3) f (2)2 1
边界条件：
0，f (1) 0 0，f (2) 0 ，f (2) 1
第二步：建立ode.m和lbc.m两个M文件。
可得：
u
u x
v
u y
1
p x
2u x2
2u y 2
)
u
v x
v
v y
1
p y
2v x2
2v y 2
)
f 2 ff a2 f ff 1 a2F f
2
（2-1） （2-2）
边界条件: 由以上三个边界条件可得:
首先解 f 2 ff a2 f ，作变量置换，令
y, f ( y) A()
/*绘图命令*/
xlabel('轴\it \eta');ylabel('轴\it \phi')
legend('平面驻点流动\phi曲线','平面驻点流动d\phi/d\eta曲线','平面驻点
流动d^2\phi/d\eta^2曲线')
title('平面驻点流动的数值解')
程序说明：
(1)solinit是被指定为x和f域的范围。x是初始网格点，f表示在节点solinit.x(i)处 f(x(i))的初始值猜测解solinit.f(:,1)，一般用bvpinit实现。
ay 0
0
0.4
0.8
1.2
1.6
2
2.4
0.0881 0.3124 0.622 0.9798 1.362 1.7553
d u d U
d 2 d 2
0 1.2326
0.4145 0.8463
0.6859 0.5251
0.8467 0.2938
0.9323 0.1474
0.9732 0.0658
(2)bvp4c是MATLAB 7.0软件求解一阶常微分方程组的库函数，调用格式如程序 所述(bvp4c的调用格式有三种，程序中的只是其中一种)。但由于sol不能直接输出 数值解，所以要用bvp4c的配置函数deval。
3、平面驻点的数值解和相应的曲线
从结果中可以看出在 2.8
之后就趋近于1，在 2.4
随着渗透系数随深度衰减系数的增大，盆地内部驻点位置越来越 深，中间和区域流动系统的面积减小，其补给区和排泄区所占面积也相 应减小，而各个局部流动系统的面积越来越大。
谢谢观赏~
2.4
a
a
（2-10）
解式 ff 1 a2F f ，可得压强 p ，对其积分可得
2
F ( y)
1 a2
(2f
f
2)
（2-11）
联立式2-2、2-7可得：
p y
a
a ( ) a
a
（2-12）
在边界层内 ,, 都只是1的数量级，因而沿壁面法线的压强梯度
p ~ a a
y 当 很小时，压强梯度也很小。
此外 p p0 a a y 表明流动过程中压力逐渐增大至 p0
第三节 利用MATLAB求数值解
1、问题描述：
平面驻点流动方程： 2 1 0
边界条件：
0， 0 0， 0 ， 1
其中:
ay
d u d U
求解 在[0,4]的范围内 ,, 的数值解。
2、求解过程
以上3个式子的流速与压强均满足势流方程，并且是不可压缩粘 性流动运动方程的精确解。但是不可压缩粘性流动的运动方程中多一
粘性项 2u 。对于势流 u ， 为流速的势函数。则
2u 2 () (2) 0
则N-S方程中的粘性项对于势流而言恒等于零。但势流解却不能 满足“无滑移”这个粘性流动的边界条件.
第三步：求解方程。
在MATLAB 7.0工作窗口输入程序：
infinity=4;
solinit=bvpinit(0:0.4:infinity,[0 0 0]);
sol=bvp4c(@ode,@lbc,solinit);
x=0:0.4:infinity
f=deval(sol,x)
plot(x,f(1,:),'ob',x,f(2,:),'rp',x,f(3,:),'b*')
0.9905 0.0265
2.8 2.153 0.99ቤተ መጻሕፍቲ ባይዱ 0.009
3.2 2.5523 0.9992 0.0028
3.6 2.954 0.999
0
弗勒塞林求解的平面驻点流动和轴对称驻点流动解：
u U A U a U() a y a
则
u ()
（2-9）
U
在η=2.4左右，ϕ’ = 0.99，即此时 粘性流动的流速已接近势流流速，只差 百分之一。可以此点距固体壁面的距离 作为边界层的厚度δ，则
第四节工程实际应用盆地尺度的地下水流动过程中多个水流系统交汇可以形成滞流区是油气聚集沉积矿产形成的重要部位驻点可以分为盆地内部驻点和盆地底界驻点两大类盆地内部驻点sp1sp2位于逆向局部水流系统的下方是四个水流系统同时发生汇聚和发散的部位
平面驻点流动（西门茨流动）
指导教师：
报告人：
目录
01
02
问题概述
（2-6）
由此方程可改写为 ayy,, f ((yy)) aA())
（2-7）
式中，a为由势流解得到的常数，υ为流体的运动粘度，为已知量。
则方程可简化为
边界条件为
2 1 0 0 0 0
（2-8）
0
方程 2 1 0 仍然是非线性的，难以求得解析解。
可采用数值解法，得下表：
时 =0.99，说明流动的流 速已经达到来流速度的99%。
第四节 工程实际应用
盆地尺度的地下水流动过程中，多个水流系统交汇可以形成滞流 区，是油气聚集、沉积矿产形成的重要部位 。
驻点可以分为盆地内部驻点和盆地底界驻点两大类，盆地内部驻 点（SP 1、SP 2）位于逆向局部水流系统的下方，是四个水流系统同 时发生汇聚和发散的部位。
方程推导与求解
03 利用Matlab求数值解
04
工程实际应用
第一节 问题概述
平面驻点流动，在平面有势流动中x、y两个方向的流速分别为：
U ax
（1-1）
V ay
（1-2）
于是由伯努利方程，可得压强分布为
p0
p
1 2
a2
x2
F ( y)
式中，p0为驻点O的压强，F ( y) y 2 .
（1-3）
可得：
f f A A y y
f
2 f y 2
(A ) y
A 2
f
3 f y3
( A 2) y
A 3
带入式可得：
2 A2 (2 ) a2 A 3
（2-3） （2-4） （2-5）
如果式2-5中 2 A2 a2 A 3 ，方程将大大简化，则需要满足
A a/ a / A a
lbc.m 文件程序如下： function res=lbc(f0,finf) res=[f0(1);f0(2);finf(2)-1];
程序说明： ode.m文件描述一阶常微分方程组，由于方程是三阶微分方程， 所以需要三个一阶微分方程来描述。 lbc.m文件描述边界条件，f0表示初始值，finf表示末端值，本问 题告诉我们的是f(1)和f(2)的初始值与f(2)的末端值，其中finf(2)-1表 示finf(2)-1=0，其它以此类推。
第二节 方程推导和求解
对于平面驻点流动，为了满足粘性流动的无滑移的条件，西门茨 给出其精确解.
现假定：
式中p0表示驻点O（x=0，y=0）处的压强；p为任意点（x，y） 处的压强；a为常数。
可证明，该假定将自动满足连续性方程。
由平面运动的N-S方程可以确定f与F两个函数，将上式带入恒定 的N-S方程
第一步：将方程化为一阶常微分方程组。
f (1)
df (1) f (2) df (2) f (3) df (3) f (1) f (3) f (2)2 1
边界条件：
0，f (1) 0 0，f (2) 0 ，f (2) 1
第二步：建立ode.m和lbc.m两个M文件。
可得：
u
u x
v
u y
1
p x
2u x2
2u y 2
)
u
v x
v
v y
1
p y
2v x2
2v y 2
)
f 2 ff a2 f ff 1 a2F f
2
（2-1） （2-2）
边界条件: 由以上三个边界条件可得:
首先解 f 2 ff a2 f ，作变量置换，令
y, f ( y) A()
/*绘图命令*/
xlabel('轴\it \eta');ylabel('轴\it \phi')
legend('平面驻点流动\phi曲线','平面驻点流动d\phi/d\eta曲线','平面驻点
流动d^2\phi/d\eta^2曲线')
title('平面驻点流动的数值解')
程序说明：
(1)solinit是被指定为x和f域的范围。x是初始网格点，f表示在节点solinit.x(i)处 f(x(i))的初始值猜测解solinit.f(:,1)，一般用bvpinit实现。
ay 0
0
0.4
0.8
1.2
1.6
2
2.4
0.0881 0.3124 0.622 0.9798 1.362 1.7553
d u d U
d 2 d 2
0 1.2326
0.4145 0.8463
0.6859 0.5251
0.8467 0.2938
0.9323 0.1474
0.9732 0.0658
(2)bvp4c是MATLAB 7.0软件求解一阶常微分方程组的库函数，调用格式如程序 所述(bvp4c的调用格式有三种，程序中的只是其中一种)。但由于sol不能直接输出 数值解，所以要用bvp4c的配置函数deval。
3、平面驻点的数值解和相应的曲线
从结果中可以看出在 2.8
之后就趋近于1，在 2.4
随着渗透系数随深度衰减系数的增大，盆地内部驻点位置越来越 深，中间和区域流动系统的面积减小，其补给区和排泄区所占面积也相 应减小，而各个局部流动系统的面积越来越大。
谢谢观赏~
2.4
a
a
（2-10）
解式 ff 1 a2F f ，可得压强 p ，对其积分可得
2
F ( y)
1 a2
(2f
f
2)
（2-11）
联立式2-2、2-7可得：
p y
a
a ( ) a
a
（2-12）
在边界层内 ,, 都只是1的数量级，因而沿壁面法线的压强梯度
p ~ a a
y 当 很小时，压强梯度也很小。
此外 p p0 a a y 表明流动过程中压力逐渐增大至 p0
第三节 利用MATLAB求数值解
1、问题描述：
平面驻点流动方程： 2 1 0
边界条件：
0， 0 0， 0 ， 1
其中:
ay
d u d U
求解 在[0,4]的范围内 ,, 的数值解。
2、求解过程
以上3个式子的流速与压强均满足势流方程，并且是不可压缩粘 性流动运动方程的精确解。但是不可压缩粘性流动的运动方程中多一
粘性项 2u 。对于势流 u ， 为流速的势函数。则
2u 2 () (2) 0
则N-S方程中的粘性项对于势流而言恒等于零。但势流解却不能 满足“无滑移”这个粘性流动的边界条件.
第三步：求解方程。
在MATLAB 7.0工作窗口输入程序：
infinity=4;
solinit=bvpinit(0:0.4:infinity,[0 0 0]);
sol=bvp4c(@ode,@lbc,solinit);
x=0:0.4:infinity
f=deval(sol,x)
plot(x,f(1,:),'ob',x,f(2,:),'rp',x,f(3,:),'b*')
0.9905 0.0265
2.8 2.153 0.99ቤተ መጻሕፍቲ ባይዱ 0.009
3.2 2.5523 0.9992 0.0028
3.6 2.954 0.999
0
弗勒塞林求解的平面驻点流动和轴对称驻点流动解：
u U A U a U() a y a
则
u ()
（2-9）
U
在η=2.4左右，ϕ’ = 0.99，即此时 粘性流动的流速已接近势流流速，只差 百分之一。可以此点距固体壁面的距离 作为边界层的厚度δ，则
第四节工程实际应用盆地尺度的地下水流动过程中多个水流系统交汇可以形成滞流区是油气聚集沉积矿产形成的重要部位驻点可以分为盆地内部驻点和盆地底界驻点两大类盆地内部驻点sp1sp2位于逆向局部水流系统的下方是四个水流系统同时发生汇聚和发散的部位
平面驻点流动（西门茨流动）
指导教师：
报告人：
目录
01
02
问题概述
（2-6）
由此方程可改写为 ayy,, f ((yy)) aA())
（2-7）
式中，a为由势流解得到的常数，υ为流体的运动粘度，为已知量。
则方程可简化为
边界条件为
2 1 0 0 0 0
（2-8）
0
方程 2 1 0 仍然是非线性的，难以求得解析解。
可采用数值解法，得下表：
时 =0.99，说明流动的流 速已经达到来流速度的99%。
第四节 工程实际应用
盆地尺度的地下水流动过程中，多个水流系统交汇可以形成滞流 区，是油气聚集、沉积矿产形成的重要部位 。
驻点可以分为盆地内部驻点和盆地底界驻点两大类，盆地内部驻 点（SP 1、SP 2）位于逆向局部水流系统的下方，是四个水流系统同 时发生汇聚和发散的部位。
方程推导与求解
03 利用Matlab求数值解
04
工程实际应用
第一节 问题概述
平面驻点流动，在平面有势流动中x、y两个方向的流速分别为：
U ax
（1-1）
V ay
（1-2）
于是由伯努利方程，可得压强分布为
p0
p
1 2
a2
x2
F ( y)
式中，p0为驻点O的压强，F ( y) y 2 .
（1-3）
可得：
f f A A y y
f
2 f y 2
(A ) y
A 2
f
3 f y3
( A 2) y
A 3
带入式可得：
2 A2 (2 ) a2 A 3
（2-3） （2-4） （2-5）
如果式2-5中 2 A2 a2 A 3 ，方程将大大简化，则需要满足
A a/ a / A a