方腔顶盖驱动流动

一、问题描述方腔顶盖驱动流动如图1所示的一个简化两维方腔（高，宽都等于L)，内部充满水分。上表面为移动墙，非维化速度为ｕ／ｕ0=1。其他三面为固定墙。试求方腔内水分流动状态。u=1,v=0u=0,v=0u=0,v=0u=0,v=0图1常微分方程理论只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法.二、离散格式数值解法:求解所有的常微分方程计算解函数y(x)在一系列节点a=x0x1…xn=b处的近似值节点间距为步长，通常采用等距节点，即取hi=h(常数)。步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。),...,1()(nixyyii欧拉方法几何意义1(,)0,1,...nnnnyyhfxyn在假设yn=y(xn)，即第n步计算是精确的前提下，考虑公式或方法本身带来的误差:Rn=y(xn+1)yn+1,称为局部截断误差.截断误差:实际上，y(xn)yn，yn也有误差，它对yn+1的误差也有影响，见下图。但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此称为方法误差或截断误差。局部截断误差的分析：由于假设yn=y(xn)，即yn准确，因此分析局部截断误差时将y(xn+1)和yn+1都用点xn上的信息来表示，工具：Taylor展开。显式欧拉公式一阶向前差商近似一阶导数推导如下：223111232()[()()()()][(,)]()()hnnnnnnnnnhnRyxyyxhyxyxOhyhfxyyxOh1()()()nnnyxyxyxh111()()()()()(,)nnnnnnnnnnyxyxhyxyxyyxyyhfxy隐式欧拉公式xn+1点向后差商近似导数推导如下：几何意义设已知曲线上一点Pn(xn,yn),过该点作弦线，斜率为(xn+1,yn+1)点的方向场f(x,y)方向,若步长h充分小，可用弦线和垂线x=xn+1的交点近似曲线与垂线的交点。11()()()nnnyxyxyxh111111()()()()()(,)nnnnnnnnnnyxyxhyxyxyyxyyhfxy比较显式公式和隐式公式及其局部截断误差显式公式隐式公式xnxn+1Pn1(,)0,1,...nnnnyyhfxyn231112()()()hnnnnRyxyyxOh111(,)nnnnyyhfxy231112()()()hnnnnRyxyyxOh中点欧拉公式中心差商近似导数控制方程hxyxyxy2)()()(0211,...,1),(211niyxfhyyiiii))(,(2)()(1102xyxfhxyxyPn+1xyy(交错网格因为方腔顶盖驱动流动的流动不均与性，u、v及压力p的变化存在交错的现象。P点位置u点位置v点位置P点控制微元体u点控制微元体v点控制微元体守恒形式N-S方程动量方程离散（x-方向）方程（1）推导过程：jijijiuuu,21,21,21综上有：动量方程离散（y-方向）方程（2）三、压力修正的基本思想压力修正方程推导方程（1）欧拉显式方程（3）方程（2）欧拉显式方程（4）压力修正方程U*,V*,P*中间值方程（5）U’,V’,P’修正值（N+1时间步）UN+1，VN+1，PN+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=MX方向网格数为：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）