一、 問(wèn)題描述方腔頂蓋驅(qū)動(dòng)流動(dòng)如圖1所示的一個(gè)簡(jiǎn)化兩維方腔（高，寬都等于L)，內(nèi)部充滿水分。上表面為移動(dòng)墻，非維化速度為0 =1。其他三面為固定墻。試求方腔內(nèi)水分流動(dòng)狀態(tài)。 u=1, v=0u=0, v=0 u=0,v=0u=0, v=0圖1常微分方程理論只能求解極少一類常微分方程；實(shí)際中給定的問(wèn)題不一定是解析表達(dá)式，而是函數(shù)表，無(wú)法用解析解法.二、離散格式數(shù)值解法:求解所有的常微分方程計(jì)算解函數(shù) y(x) 在一系列節(jié)點(diǎn) a = x0 x1 xn= b處的近似值節(jié)點(diǎn)間距 為步長(zhǎng)，通常采用等距節(jié)點(diǎn)，即取 hi = h (常數(shù))。步進(jìn)式：根據(jù)已知的或已求出的節(jié)點(diǎn)上的函數(shù)值計(jì)算當(dāng)前節(jié)點(diǎn)上的函數(shù)值，一步一步向前推進(jìn)。因此只需建立由已知的或已求出的節(jié)點(diǎn)上的函數(shù)值求當(dāng)前節(jié)點(diǎn)函數(shù)值的遞推公式即可。歐拉方法幾何意義截?cái)嗾`差: 實(shí)際上，y(xn) yn， yn 也有誤差，它對(duì)yn+1的誤差也有影響，見(jiàn)下圖。但這里不考慮此誤差的影響，僅考慮方法或公式本身帶來(lái)的誤差，因此稱為方法誤差或截?cái)嗾`差。局部截?cái)嗾`差的分析：由于假設(shè)yn = y(xn) ，即yn準(zhǔn)確，因此分析局部截?cái)嗾`差時(shí)將y(xn+1) 和 yn+1都用點(diǎn)xn上的信息來(lái)表示，工具：Taylor展開(kāi)。在假設(shè) yn = y(xn)，即第 n 步計(jì)算是精確的前提下，考慮公式或方法本身帶來(lái)的誤差: Rn = y(xn+1) - yn+1 , 稱為局部截?cái)嗾`差.顯式歐拉公式一階向前差商近似一階導(dǎo)數(shù)推導(dǎo)如下：隱式歐拉公式xn+1點(diǎn)向后差商近似導(dǎo)數(shù)推導(dǎo)如下：幾何意義設(shè)已知曲線上一點(diǎn) Pn (xn , yn ),過(guò)該點(diǎn)作弦線，斜率為(xn+1 , yn +1 ) 點(diǎn)的方向場(chǎng)f(x,y)方向,若步長(zhǎng)h充分小，可用弦線和垂線x=xn+1的交點(diǎn)近似曲線與垂線的交點(diǎn)。xnxn+1Pn比較 顯式公式和隱式公式及其局部截?cái)嗾`差顯式公式隱式公式Pn+1xyy(x)中點(diǎn)歐拉公式中心差商近似導(dǎo)數(shù)控制方程交錯(cuò)網(wǎng)格因?yàn)榉角豁斏w驅(qū)動(dòng)流動(dòng)的流動(dòng)不均與性，u、v及壓力p的變化存在交錯(cuò)的現(xiàn)象。P 點(diǎn)位置u 點(diǎn)位置v 點(diǎn)位置P 點(diǎn)控制微元體u點(diǎn)控制微元體v點(diǎn)控制微元體守恒形式N-S方程動(dòng)量方程離散（x-方向）方程（1）推導(dǎo)過(guò)程：綜上有：動(dòng)量方程離散（y-方向）方程（2）三、壓力修正的基本思想壓力修正方程推導(dǎo)方程（1）歐拉顯式方程（3）方程（2）歐拉顯式方程（4）壓力修正方程U*,V*,P* 中間值方程（5）U,V,P修正值 （N+1時(shí)間步）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程序語(yǔ)言: 已知條件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é)點(diǎn)在Y方向序號(hào)變量為 i。 1 =i= MX方向網(wǎng)格數(shù)為： N。 節(jié)點(diǎn)在X方向序號(hào)變量為 j。 1= j= N節(jié)點(diǎn)矩陣：OT=zeros(M,N); OT 舊值 某個(gè)節(jié)點(diǎn) OT(i, j)NT=zeros(M,N); NT 新值間距dx=L/(N-1);