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