方腔顶盖驱动流动

一、
二、问题描述
方腔顶盖驱动流动
如图1所示的一个简化两维方腔（高，宽都等于L)，内部充满水分。

上表面为移动墙，非维化速度为ｕ／ｕ0 =1。

其他三面为固定墙。

试求方腔内水分流动状态。

u=1, v=0
u=0, v=0
u=0,v=0
u=0, v=0
图1
常微分方程理论
只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法.
二、离散格式
数值解法:求解所有的常微分方程 计算解函数 y(x) 在一系列节点
a = x 0< x 1<…<
x n = b 处的近似值
)
,...,1()
(n i x y y i i =≈
节点间距
为步长，通常采
用等距节点，即取 hi = h (常数)。

步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。

因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。

欧拉方法
1(,) 0,1,...
n n n n y y h f x y n +=+=
几何意义
在假设 y n = y (x n )，即第 n 步计算是精确的前提下，考虑公式或方法本身带来的误差: R n = y (x n +1)
y n +1 , 称为局部截
断误差.
显式欧拉公式
一阶向前差商近似一阶导数
22
3
1112
3
2
()[()()()()]
[ (,)] ()()
h n n n n n n n n n h n R y x y y x hy x y x O h y hf x y y x O h +++'''=-=+++-+''=
+
推导如下：
隐式欧拉公式
x n +1点向后差商近似导数 推导如下：
1()()
()n n n y x y x y x h
+-'≈
111()()() ()()(,)
n n n n n
n n n n n y x y x hy x y x y y x y y h f x y +++'≈+↑≈≈=+11()()
()n n n y x y x y x h
++-'≈
11()()()
()n n n n n
y x y x hy x y x y ++'≈+↑≈
几何意义
设已知曲线上一点 P n (x n , y n ),过该点作弦线，斜率为(x n +1 , y n +1 ) 点的方向场f (x ,y )方向,若步长h 充分小，可用弦线和垂线x =x n +1的交点近似曲线与垂线的交点。

比较 显式公式和隐式公式及其局部截断误差
显式公式
1(,) 0,1,...
n n n n y y h f x y n +=+=23
1112
()()()
h n n n n R y x y y x O h +++''=-=
+
隐式公式
中点欧拉公式
中心差商近似导数
h x y x y x y 2)
()()(021-≈
'
1
,...,1),(211-=+=-+n i y x f h y y i i i i
111(,)
n n n n y y h f x y +++=+2
3
1112
()()()
h n n n n R y x y y x O h +++''=-=-+))
(,(2)()(1102x y x f h x y x y +≈P n+1
x y
y(
控制方程
交错网格
因为方腔顶盖驱动流动的流动不均与性，u、v及压力p的变化存在交错的现象。

P 点位置
u 点位置
v 点位置
P 点控制微元体u点控制微元体v点控制微元体
守恒形式N-S方程
动量方程离散（x-方向）
方程（1）
推导过程：
⎪⎪
⎭
⎫
⎝
⎛
+
=
+
-j
i
j
i
j
i
u
u
u
,
2
1
,
2
1
,2
1
综上有：
动量方程离散（y-方向）
方程（2）
三、压力修正的基本思想
压力修正方程推导
方程（1）
欧拉显式
方程（3）
方程（2）欧拉显式
方程（4）
压力修正方程
U*,V*,P* 中间值
方程（5）U’,V’,P’修正值（N+1时间步）U N+1，V N+1，P N+1，满足连续方程的值
方程（6）
方程（7）
方程（8）
综合方程3,5,6，可得
方程（9）
综合方程4,5,7，可得
方程（10）将方程9,10,5 代入方程8，可得
方程（11）
四、SIMPLE算法流程图
M程序语言:
已知条件
L=H; 方腔的宽，高
u=1, v=0；方腔上表面的速度
u=0, v=0；方腔左表面的速度
u=0, v=0；方腔右表面的速度
u=0, v=0；方腔下表面的速度
网格划分
Y方向网格数为：M。

节点在Y方向序号变量为i。

1< =i<= M
X方向网格数为：N。

节点在X方向序号变量为j。

1<= j<= N
节点矩阵：
OT=zeros(M,N); OT 旧值某个节点OT(i, j)
NT=zeros(M,N); NT 新值间距
dx=L/(N-1); X方向间距dy=H/(M-1); Y方向间距控制方程离散
方程（11）。