方腔顶盖驱动流动

一、问题描述

方腔顶盖驱动流动

如图1所示的一个简化两维方腔（高，宽都等于L)，内部充满水分。上表面为移动墙，非维化速度为ｕ／ｕ0 =1。其他三面为固定墙。试求方腔内水分流动状态。

u=1, v=0

u=0,v=0

u=0, v=0

图1

常微分方程理论

只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法.

二、离散格式

数值解法:求解所有的常微分方程

计算解函数 y(x) 在一系列节点a = x0< x1<…<

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的交点近似曲线与垂线的交点。

比较 显式公式和隐式公式及其局部截断误差

显式公式

x n x n+1

P n

P n+1

x

y y(x )

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 +++=+23

1112

()()()

h n n n n R y x y y x O h +++''=-=-

+))

(,(2)()(1102x y x f h x y x y +≈

x 0 x 2

x 1

控制方程

交错网格

因为方腔顶盖驱动流动的流动不均与性，u、v及压力p的变化存在交错的现象。

P 点位置

u 点位置

v 点位置

P 点控制微元体u点控制微元体v点控制微元体

守恒形式N-S方程

动量方程离散（x-方向）

方程（1）

推导过程：

⎪⎪⎭

⎫

⎝⎛

+=+-j

i j i j

i u u u ,21,21,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；方腔上表面的速度