方腔顶盖驱动流动.doc

1、一、问题描述方腔顶盖驱动流动如图 1 所示的一个简化两维方腔（高，宽都等于L)，内部充满水分。上表面为移动墙，非维化速度为 0 =1。其他三面为固定墙。试求方腔内水分流动状态。u=1, v=0u=0, v=0u=0,v=0u=0, v=0图 1常微分方程理论只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法.二、离散格式数值解法 :求解所有的常微分方程计算解函数 y(x)在一系列节点a = x0 x1xn= b 处的近似值yiy( xi ) ( i 1, . , n)节点间距为步长，通常采用等距节点，即取hi = h (常数 )。步进式：根据已知的或已求

2、出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。 因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。欧拉方法yn 1 yn h f ( xn , yn )n 0, 1,.几何意义在假设 yn = y(xn)，即第 n 步计算是精确的前提下，考虑公式或方法本身带来的误差 : Rn = y(xn+1) yn+1 , 称为局部截断误差 .截断误差 :实际上， y( x )y ， y 也有误差，它对 y的误差也有影响，见下nnnn+1图。但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此称为方法误差或截断误差。局部截断误差的分析 ：由于假设 y = y

3、 ( x ) ，即 y 准确，因此分析局部截断误nnn差时将 y( x)和 y都用点 x 上的信息来表示，工具：Taylor 展开 。n+1n+1nRn 1 y(xn 1) yn 1 y(xn )2hy(xn) h2 y (xn) O(h3)ynhf(xn, yn )h2y (x ) O(h3)2n显式欧拉公式一阶向前差商近似一阶导数y ( xn 1 )y ( xn )y ( xn )h推导如下：y( xn1 )y( xn ) hy ( xn )y( xn )yny( xn1 )yn 1ynh f ( xn , yn )隐式欧拉公式x 点向后差商近似导数n+1推导如下：y ( x n 1 )y

4、 ( x n )y ( x n 1 )hy(xn 1)y(xn )hy (xn 1)y(xn )yny(xn 1)yn 1ynh f (xn 1, yn 1)几何意义设已知曲线上一点P ( x , y ), 过该点作弦线，斜率为 ( x,y)nnnn+1n +1点的方向场 f ( x, y) 方向 , 若步长 h 充分小，可用弦线和垂线x=x的n+1交点近似曲线与垂线的交点。xn+1xnPn比较显式公式和隐式公式及其局部截断误差显式公式y n 1y nh f ( xn , y n )n0, 1,.Rn 1 y ( xn 1 ) yn 123 )h2 y ( xn ) O ( h隐式公式yn 1

5、ynh f ( xn 1 , yn 1 )Ry(x)yh2y (x ) O(h3)n 1n 1n 12n中点欧拉公式中心差商近似导数xy(yy ( x 1y ( x 2 ) y ( x 0 )2 hy( x2 )y( x0)2h f ( x1, y( x1)yi 1 yi 1 2hf(xi,yi) i 1,.,n 1控制方程交错网格因为方腔顶盖驱动流动的流动不均与性，u、v 及压力 p 的变化存在交错的现象。P 点位置u 点位置v 点位置P 点控制微元体u 点控制微元体v 点控制微元体守恒形式 N-S 方程动量方程离散（ x-方向）方程（ 1）推导过程：ui , j1u1u12i2, ji2,

6、 j综上有：动量方程离散（ y-方向）方程（ 2）三、压力修正的基本思想压力修正方程推导方程（ 1）欧拉显式方程（ 3）方程（ 2）欧拉显式方程（ 4）压力修正方程U*,V*,P*中间值方程（ 5）U ,V 修,P正值 （N+1 时间步）N+1N+1N+1U ，V ， P ，满足连续方程的值方程（ 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= MX 方向网格数为：N。节点在 X 方向序号变量为j。 1= j= N节点矩阵：OT=zeros(M,N); OT 旧值某个节点OT(i, j)NT