水利工程論文-非正交曲線坐標(biāo)下二維水流計(jì)算的SIMPLEC算法.doc

水利工程論文-非正交曲線坐標(biāo)下二維水流計(jì)算的SIMPLEC算法摘要：本文采用Laplace方程坐標(biāo)變換方法生成正交曲線網(wǎng)格，并對(duì)淺水流動(dòng)的控制方程進(jìn)行坐標(biāo)變換，方程離散時(shí)采用B型交錯(cuò)網(wǎng)格。利用“水位掃描法”結(jié)合壁面函數(shù)法來處理移動(dòng)邊界，用SIMPLEC算法解非正交曲線坐標(biāo)下的k-雙方程紊流模型，修正了由網(wǎng)格的非正交性引起的誤差。通過對(duì)美國Colorado洲FallRiver的資料進(jìn)行流場(chǎng)驗(yàn)證，計(jì)算結(jié)果與實(shí)測(cè)資料基本符合，顯示了本模型在不規(guī)則水域計(jì)算中的實(shí)用價(jià)值。關(guān)鍵詞：坐標(biāo)變換k-紊流模型水位掃描法壁面函數(shù)SIMPLEC算法隨著經(jīng)濟(jì)發(fā)展和社會(huì)進(jìn)步，水利工程建設(shè)的步伐也在進(jìn)一步加快，其中港航建設(shè)、大壩建設(shè)中的泥沙問題以及近來倍受世人關(guān)注的水污染問題已經(jīng)成為制約水利發(fā)展的瓶頸問題，弄清河流、湖泊、海洋中水動(dòng)力因素，是解決以上問題的重要基礎(chǔ)。近年來，數(shù)學(xué)模型已逐步取代物理模型實(shí)驗(yàn)成為研究水流的重要手段，而淺水流動(dòng)模型是處理大區(qū)域流場(chǎng)的一種非常有效的模型。它屬于非線性方程組，在目前只能用數(shù)值方法求解，因此，有必要研究一種簡單、高效的方法來求解淺水流動(dòng)問題。自Patankar和Spalding1發(fā)展了SIMPLE算法以來，該方法被廣泛應(yīng)用于不可壓縮流場(chǎng)的數(shù)值模擬，而且該方法還得到了進(jìn)一步的發(fā)展，主要有SIMPLER算法2、SIMPLEC算法3、SIMPLEX算法4和SIMPLET算法5等。這些模型均成功地應(yīng)用于速度壓力耦合的流場(chǎng)計(jì)算，深度平均的淺水流動(dòng)模型是在靜壓假定下導(dǎo)出的，一般流體模型中的速度壓力耦合也就轉(zhuǎn)換成淺水流動(dòng)模型中的速度水深耦合6。天然河流、海灣的邊界曲折、地形復(fù)雜，采用坐標(biāo)變換是解決問題的途徑之一。目前多數(shù)N-S方程的坐標(biāo)變換中，流程全部采用逆變分量，這樣就增加了方程的復(fù)雜程度。于是忽略掉方程中的非正交項(xiàng)，利用正交變換下的方程進(jìn)行數(shù)值求解7,8。對(duì)于具有復(fù)雜邊界的海灣及彎曲的河流，坐標(biāo)變換中很難保證每個(gè)點(diǎn)都正交，特別是邊界附近。水位變化是水力計(jì)算中難點(diǎn)之一，在目前的紊流數(shù)學(xué)模型中，多簡單的利用“凍結(jié)法”，這樣做將失去對(duì)邊界出流動(dòng)模擬的準(zhǔn)確性。本文研究中，采用正交曲線坐標(biāo)變換生成數(shù)值網(wǎng)格，而數(shù)值計(jì)算中采用非正交曲線坐標(biāo)下的k-雙方程紊流模型，這樣可以自動(dòng)修正網(wǎng)格生成中的非正交項(xiàng)。流速除對(duì)流項(xiàng)中采用逆變分量，在其余各項(xiàng)中均采用原始分量，這樣使得方程書寫簡單，有利于將各方程寫成通用形式，編寫的程序變得更規(guī)范。作者受JianYe同位網(wǎng)格9的啟發(fā)，對(duì)普通交錯(cuò)網(wǎng)格做了修改，即采用B型交錯(cuò)網(wǎng)格，使得u,v,k,的計(jì)算布置在一個(gè)節(jié)點(diǎn)上，有利于節(jié)省計(jì)算程序代碼，使程序書寫更加規(guī)范。引入動(dòng)邊界掃描技術(shù)，結(jié)合紊流模型的壁面函數(shù)法，使壁面隨著真實(shí)邊界而變化。數(shù)值求解時(shí)，采用控制體積法離散方程，運(yùn)用SIMPLEC算法，使計(jì)算的流場(chǎng)更符合實(shí)際流場(chǎng)。1數(shù)值網(wǎng)格本文對(duì)計(jì)算區(qū)域用Laplace方程實(shí)施坐標(biāo)變換，生成正交的貼體網(wǎng)格，控制方程：(1)方程組(1)可采用SOR方法求解，然后用正交曲線網(wǎng)格邊界正交化處理的三次樣條插值邊界滑動(dòng)法10，處理其邊界處的網(wǎng)格，以提高其正交性。與其它方法相比，此方法雖在一定程度上提高了邊界處網(wǎng)格的正交性，但難以保證完全正交，如果再把曲線坐標(biāo)系下N-S方程中的非正交項(xiàng)去掉，由斜交網(wǎng)格引起的計(jì)算誤差是不容忽視的。為了修正此項(xiàng)誤差，本文采用非正交曲線坐標(biāo)下的平面二維k-雙方程紊流模型，模型能自動(dòng)修正網(wǎng)格生成中的非正交項(xiàng)。2數(shù)學(xué)模型笛卡兒坐標(biāo)下深度平均的k-雙方程紊流模型的通用微分方程：(2)方程(2)轉(zhuǎn)換到曲線坐標(biāo)(、)下，僅在對(duì)流項(xiàng)中使用流速的逆變分量，而在其它項(xiàng)中使用原始變量，這樣既簡化了方程，又使所有方程仍可寫為曲線坐標(biāo)下的通用方程，模型的微分方程可寫為如下通用形式：(3)式中：為所求問題的因變量；U和V分別為直角坐標(biāo)下流速u和v的逆變分量，僅在對(duì)流項(xiàng)中出現(xiàn)；為擴(kuò)散系數(shù)；S為源項(xiàng)。當(dāng)表示某一特定量時(shí)，和S對(duì)此特定量有特定的意義和表達(dá)式，這時(shí)方程(3)亦賦予特定的意義，模型的控制方程組如表1所示。式中：u和v分別為笛卡兒坐標(biāo)下x和y方向的深度平均的流速分量；k和分別為深度平均的紊動(dòng)動(dòng)能和紊動(dòng)動(dòng)能耗散率；H為水深；zs為水位；zb為河床高程；為分子運(yùn)動(dòng)粘性系數(shù)；t為漩渦運(yùn)動(dòng)粘性系數(shù)；bx和by分別為x和y方向上的底摩擦力eff為有效粘性系數(shù)；C、C1、C2、k和為經(jīng)驗(yàn)常數(shù)，其中C為0.09，C1為1.44，C2為1.92，k為1.0，為1.3。3模型離散如圖1所示，本文采用B型交錯(cuò)網(wǎng)格，即將u,v,k，及各標(biāo)量布置在網(wǎng)格中心，將水深相關(guān)項(xiàng)H、zs、zb項(xiàng)布置在網(wǎng)格的左下方節(jié)點(diǎn)上。這樣布置有以下兩個(gè)優(yōu)點(diǎn)：有利于邊界條件的引入；在編寫程序u,v,k,各方程共用一套系數(shù)，可節(jié)省程序代碼，使程序簡單明了。本文用控制容積法，采用混合格式來導(dǎo)出離散方程。其邊界條件為：在進(jìn)口邊界，所有邊界條件都按本征條件給出，即：u=u0,v=0,k=k0,=0；在出口邊界，給定出口處的水位zs，令v=0，其余變量都取法向?qū)?shù)為零，即：在水陸邊界，在計(jì)算河流、河口、海灣的非恒定流場(chǎng)時(shí)，水位不是恒定不變的，而且每一時(shí)間步長的水位可能都不一樣，這就帶來了水陸邊界的移動(dòng)，如果每一時(shí)間步長都進(jìn)行坐標(biāo)變換而生成新的網(wǎng)格是非常不經(jīng)濟(jì)的。于是本文在前人“凍結(jié)法”的基礎(chǔ)上，提出了適合k-雙方程紊流模型的“移動(dòng)邊界的壁面函數(shù)法”。如圖2、圖3，水位變化后，使邊界處的網(wǎng)格干出，原來的岸邊界12、34變成了新的岸邊界PQ、ST，而壁面函數(shù)仍然布置在12、34岸邊界處；水域內(nèi)島嶼由于水位的變化而出露，其島嶼內(nèi)部按“凍結(jié)法”來處理，其邊界ABCD周圍也需用壁面函數(shù)法來處理。本文采用“水位掃描法”來判斷新的水陸邊界：從邊界12向34掃描，如水深不為零，則令其為起始點(diǎn)；從邊界34向12掃描，如水深不為零，則令其為終點(diǎn)；然后掃描水域內(nèi)部，確定島嶼的邊界ABCD。之后對(duì)水陸邊界實(shí)施壁面函數(shù)法。在壁面邊界，k-雙方程紊流模型是高Reynolds數(shù)模型，只適用于充分發(fā)展的紊流，而在近壁處，粘性效應(yīng)起主導(dǎo)作用。另外，由于近壁處物理量變化很快，需要非常密的網(wǎng)格，才能真實(shí)地模擬實(shí)際流動(dòng)。因此本文采用壁面函數(shù)法11，即用半分析的方法得到的解來近似由壁面到紊流核心區(qū)之間的流速、紊動(dòng)動(dòng)能、紊動(dòng)動(dòng)能耗散率的分布規(guī)律，將壁面的影響(如壁面應(yīng)力)附加到差分方程中(應(yīng)先將有關(guān)的邊界系數(shù)置為零)12。4SIMPLEC算法求解模型SIMPLE算法的基本步驟如下：(1)計(jì)算坐標(biāo)變換的相關(guān)系數(shù)；(2)根據(jù)進(jìn)、出口邊界的水位，給定全場(chǎng)水位，計(jì)算全場(chǎng)初始水深H*；(3)解動(dòng)量方程求un,vn；(4)解k-紊流模型求k、；(5)解水深校正方程求H，并修正水深H=H*+hH；(6)修正流速u、v及有效粘性系數(shù)eff，然后重復(fù)(2)(5)步直至收斂。模型求解中，為了利于非線性迭代的收斂，計(jì)算中采用亞松弛技術(shù)。計(jì)算中采用ADI技術(shù)和TDMA算法，以連續(xù)方程的誤差小于給定的值作為判斷收斂的依據(jù)。在恒定流計(jì)算中，以相鄰時(shí)間步之間全場(chǎng)各點(diǎn)流速誤差的最大絕對(duì)值小于給定的值作為達(dá)到定常狀態(tài)的依據(jù)。5模型應(yīng)用為了使本文數(shù)學(xué)模型具有代表意義，本文采用天然連續(xù)彎道河流美國Colo