方腔頂蓋驅(qū)動流動

1、一、 問題描述方腔頂蓋驅(qū)動流動如圖1所示的一個簡化兩維方腔（高，寬都等于L)，內(nèi)部充滿水分。上表面為移動墻，非維化速度為0 =1。其他三面為固定墻。試求方腔內(nèi)水分流動狀態(tài)。 u=1, v=0u=0, v=0 u=0,v=0u=0, v=0圖1常微分方程理論只能求解極少一類常微分方程；實際中給定的問題不一定是解析表達式，而是函數(shù)表，無法用解析解法.二、離散格式數(shù)值解法:求解所有的常微分方程計算解函數(shù) y(x) 在一系列節(jié)點 a = x0< x1<< xn= b處的近似值節(jié)點間距 為步長，通常采用等距節(jié)點，即取 hi = h (常數(shù))。步進式：根據(jù)已知的或已求出的節(jié)點上的函數(shù)值計

2、算當前節(jié)點上的函數(shù)值，一步一步向前推進。因此只需建立由已知的或已求出的節(jié)點上的函數(shù)值求當前節(jié)點函數(shù)值的遞推公式即可。歐拉方法幾何意義截斷誤差: 實際上，y(xn) » yn， yn 也有誤差，它對yn+1的誤差也有影響，見下圖。但這里不考慮此誤差的影響，僅考慮方法或公式本身帶來的誤差，因此稱為方法誤差或截斷誤差。局部截斷誤差的分析：由于假設(shè)yn = y(xn) ，即yn準確，因此分析局部截斷誤差時將y(xn+1) 和 yn+1都用點xn上的信息來表示，工具：Taylor展開。在假設(shè) yn = y(xn)，即第 n 步計算是精確的前提下，考慮公式或方法本身帶來的誤差: Rn = y(x

3、n+1) - yn+1 , 稱為局部截斷誤差.顯式歐拉公式一階向前差商近似一階導(dǎo)數(shù)推導(dǎo)如下：隱式歐拉公式xn+1點向后差商近似導(dǎo)數(shù)推導(dǎo)如下：幾何意義設(shè)已知曲線上一點 Pn (xn , yn ),過該點作弦線，斜率為(xn+1 , yn +1 ) 點的方向場f(x,y)方向,若步長h充分小，可用弦線和垂線x=xn+1的交點近似曲線與垂線的交點。xnxn+1PnPn+1xyy(x)比較 顯式公式和隱式公式及其局部截斷誤差顯式公式隱式公式x0x2x1中點歐拉公式中心差商近似導(dǎo)數(shù)控制方程交錯網(wǎng)格因為方腔頂蓋驅(qū)動流動的流動不均與性，u、v及壓力p的變化存在交錯的現(xiàn)象。P 點位置u 點位置v 點位置P 點

4、控制微元體u點控制微元體v點控制微元體守恒形式N-S方程動量方程離散（x-方向）方程（1）推導(dǎo)過程：綜上有：動量方程離散（y-方向）方程（2）三、壓力修正的基本思想壓力修正方程推導(dǎo)方程（1）歐拉顯式方程（3）方程（2）歐拉顯式方程（4）壓力修正方程U*,V*,P* 中間值方程（5）U,V,P修正值 （N+1時間步）UN+1，VN+1， PN+1，滿足連續(xù)方程的值方程（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； 方腔下表面的速度網(wǎng)格劃分Y方向網(wǎng)格數(shù)為： M。節(jié)點在Y方向序號變量為 i。 1< =i<= MX方向網(wǎng)格數(shù)為： N。 節(jié)點在X方向序號變量為 j。 1<= j<= N節(jié)點矩陣：OT=zeros(M,N); OT 舊值 某個節(jié)點 OT(i