【文章內(nèi)容簡介】 自由水面，本質(zhì)上是一種有限差分方法。為了描述波浪自由表面，Nichols（1975，1980）引入了流體體積函數(shù)的概念。定義函數(shù)F=F（x，y，z，t），表示單元體內(nèi)流體的體積占該單元體體積的比值。對于某個計算單元，F(xiàn)=1表示該單元體為流體所充滿；F=0表示它是一個空單元；F值介于0與1之間，表示該單元體被流體部分充滿，這又有兩種可能情況，它或者與自由表面相交，或者含有比單元體尺度小的空氣泡。F的梯度用來決定自由邊界的法線方向。得出每個計算單元的F數(shù)值及其梯度之后，就可確定自由邊界的近似位置。F函數(shù)的控制方程在形式上可表示為： （221） 該方法的最大優(yōu)點在于只需要一個函數(shù)就可描述自由表面的各種復(fù)雜變化。VOF方法基于最基本的流體力學方程，求解含有自由表面的流動問題具有較大的方便性和靈活性。該方法克服了高度函數(shù)法無法處理自由表面是坐標的多值函數(shù)情況的缺點，而且克服了所用計算機內(nèi)存和計算時間較多的缺點。另外，由于又經(jīng)許多學者的完善和改進，目前成為研究波浪問題中應(yīng)用廣泛、很有發(fā)展前途的數(shù)值模擬方法之一。本文采用VOF法來追蹤自由表面。 建筑物的構(gòu)建研究波浪與建筑物的相互作用，在模型中正確構(gòu)建建筑物也是很重要的。本模型中采用曲線來表示二維建筑物與水體交界面。設(shè)曲線為： （222） 數(shù)值造波的方法數(shù)值造波是波浪數(shù)學模型中一個重要的方面。過去人們常采用邊界造波的方法，但是邊界造波同物理模型試驗一樣存在造波邊界的二次反射問題。在邊界上造波的同時吸收來自域內(nèi)的反射波，無論對物理模型試驗還是對數(shù)值模型試驗都是具有挑戰(zhàn)性的難題。在基于NS方程和雷諾平均應(yīng)力NS方程的數(shù)值波浪模型中處理吸收造波邊界甚至更為困難。因為這種模型在造波邊界上需要自由水面位移的同時也需要速度分布，因此，模型對來自邊界的誤差比深度平均的波浪模型更為敏感。人們雖然作了很大的努力來改進造波吸收邊界條件，然而結(jié)果還不是很令人滿意，如Petit采用所謂的弱反射邊界條件的同時吸收弱反射波。由于在推導(dǎo)中引入了線性疊加的假定，該方法僅適用于小振幅波。而且Petit在長時間的計算中也經(jīng)歷了較大的數(shù)值誤差。Lin提出在計算域內(nèi)部通過控制來自邊界的質(zhì)量流保證質(zhì)量守恒，這樣可長時間計算而沒有導(dǎo)致數(shù)值不穩(wěn)定，但這種方法也僅適用于弱反射或反射波能夠事先估計的情況。另一種改進的造波方法是計算域內(nèi)源函數(shù)造波。設(shè)計的源函數(shù)不與反射波相互作用，此時在邊界上可以采取輻射條件或完全吸收條件。Wei 和Kirby將這種思想應(yīng)用于平面二維的波浪模型。Lin導(dǎo)出了在垂向二維波浪模型中各種永形波的源函數(shù)形式，并研究了這種方法的造波特性。該方法適用于長時間波浪與海岸或建筑物作用的研究，尤其是在反射波較大時的情況。由于可以根據(jù)波浪條件選擇合適的永形波波浪理論，該方法適用的造波范圍很寬。在定常水深情況下，對各種波浪理論的波形和速度的對比，得出源函數(shù)造波精度很高。下面先介紹這種方法的原理，再介紹各種永形波的造波源函數(shù)。描述波浪水流運動的基本控制方程見式（21）和（22）。為了用質(zhì)量源函數(shù)產(chǎn)生波浪，將方程（21）式改寫成： （223）式中為在域Ω上為非零的質(zhì)量源函數(shù)。在采用有限差分的計算中，僅使用一個有限的質(zhì)量源域，最小的源域為一個計算網(wǎng)格單元。通常使用一個長方形的源域，其實際長度為、高度為。設(shè)為目標波的波長，當時，源域可以看作是點源。在試驗中，通常假定。按原理，只要指定合理的源函數(shù)，任何波都能夠生成。得到的是，當源域很窄時，目標波可以看作是從源域內(nèi)產(chǎn)生的。假定所有質(zhì)量的增加或減小完全由質(zhì)量源產(chǎn)生，則下列關(guān)系必須得到滿足 （224）式中C為目標波的波速，系數(shù)2表示波浪在源域的兩邊同時產(chǎn)生，為水面位移。 邊界條件的種類造波板（即入射）邊界條件在入射邊界，令邊界虛擬網(wǎng)格中與邊界平行的流速分量等于邊界第一排物理網(wǎng)格中相應(yīng)的流速分量，虛擬網(wǎng)格中與邊界垂直的流速分量等于造波板的速度。可滑移邊界條件這種邊界條件是令邊界虛擬網(wǎng)格中與邊界平行的流速分量等于邊界第一排物理網(wǎng)格相應(yīng)的流速分量，虛擬網(wǎng)格中與邊界垂直的流速分量等于零。開放邊界條件 這種邊界條件主要用于數(shù)值波浪水槽右端。開放邊界是一種基于波浪輻射條件的邊界條件，它能夠?qū)⒃觳C端傳播過來的行進波透射出去。吸收波浪邊界條件在實際計算中對于吸收波浪邊界，必須有效的吸收波浪能量，使波浪能量不能或較少的被反射傳回計算區(qū)域內(nèi)，從而保證模擬波浪的精度。一般的波浪吸收措施是在右端邊界前設(shè)置海綿層。實際計算表明，海綿層可以很好的吸收波浪的能量，即使有少量的反射，反射波也需再次通過海綿層而損耗貽盡，從而在有效的計算區(qū)域內(nèi)反射波浪的能量可以忽略。本文中的計算程序是在右端邊界面前設(shè)置一段衰減函數(shù)，并且只對垂向速度進行衰減，而水平速度不做衰減，從而保證流動的連續(xù)性。 本章小結(jié)本章從垂向二維波浪數(shù)學模型的各個方面綜合敘述了該模型的概況，其中包括基本控制方程，數(shù)值模型的實現(xiàn)，建筑物的構(gòu)建，數(shù)值造波方法以及邊界條件的處理等，總體上反映了該模型的特點和理論背景。第3章 NEWFLUM軟件說明NEWFLUM軟件需要在Visual Fortran ，由于該程序較早開發(fā)，共有13725行，且使用Fixed Format（固定格式）編寫的，相對于Free Format（自由格式）需要遵循較多的格式規(guī)范，并規(guī)定了程序代碼每一行中每一個字段的意義。用戶若有看不懂之處可以查閱《Fortran 95程序設(shè)計》。本章將運用第二章所介紹的垂向二維波浪數(shù)學模型原理，熟悉和理解數(shù)值模擬的程序，編制直立堤的INPUT文件，運行該程序來加深理解數(shù)值模擬的程序及各個變量參數(shù)的數(shù)值和實際含義，為最終分析弧面防波堤的反射透射特性特性奠定基礎(chǔ)。 數(shù)值模型程序說明 直立堤的數(shù)模輸入程序根據(jù)不同的應(yīng)用需要，用戶既可將這一模型作為一個研究工具，也可將其作為一個黑箱。但無論那種方式，用戶首先需要了解INPUT文件的結(jié)構(gòu)，這也是正確使用模型的關(guān)鍵。本節(jié)將詳細解釋INPUT文件，作為連接用戶和模型的橋梁。實例為破碎波越過一個不可滲透的沉箱，由有孔的護面層保護著的沉箱將被用來說明怎樣使用INPUT文件來解決指定的的問題。問題的素描見圖31和圖32。圖31 波浪越堤的試驗?zāi)Ｐ褪疽鈭D在這實例中，規(guī)則波的波浪條件為：T= s，H= m，由左邊界造波，造波處水深h= m，為一常數(shù)。水槽右邊界為不可滲海灘，其坡度為1/25。右邊界水深hr = m， m（見圖31），沉箱的周圍是由三種不同孔隙率材料構(gòu)成的保護層（詳細的保護層構(gòu)成見圖32）。圖32 試驗的沉箱及有孔護面層的構(gòu)建為了模擬這個問題，相應(yīng)的INPPUT文件如下：Wave breaking and overtopping on porous structure and caisson ZLD $numparam , ==delt,twfin,autot 3,3,1,1 ==kl,kr,kt,kb$end$fldparam ,, ==xnu,gx,gy,ui,vi,rhof$end$mesh 1,400,600 ==nkx,xl(1~nkx),xl(nkx+1),xc(1~nkx),nxl(1~nkx),nxr(1~nkx) ==dxmn(1~nkx) 1,30,44 ==nky,yl(1~nky),yl(nky+1),yc(1~nky),nyl(1~nky),nyr(1~nky) ==dymn(1~nky)$end$obstcl 2 ==nobstype 2 ==nobs , , , , ==oa2(),oa1(),ob2(),ob1() , , 1,0 ==oc2(),oc1(),ioh() , , , , ==od2(),od1(),oe2(),oe1() 0,0, 0,0, 0,0, 0,0 ==nxo(),mxo(),nyo(),myo() 4 ==nobs 0.,0.,0.,0., 0.,1.,0.,1., 0.,0.,0.,0., 1.,0.,1.,0. ==oa2(),oa1(),ob2(),ob1() 0.,0.,0.,0., , 1,0,0,0 ==oc2(),oc1(),ioh() 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0. ==od2(),od1(),oe2(),oe1() 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0 ==nxo(),mxo(),nyo(),myo()$end$porous material information$ 1 ==npor 9 ==nportype ==d50(1),xporosity(1) [armor layer1] 4 ==nporous(1) 0.,0.,0.,0., 0.,1.,0., 0.,0.,0.,0., 1.,1.,0.,1. ==pa2,pa1,pb2,pb1(1~nporous) 0.,0.,0.,0., , 1,0,0,0 ==pc2,pc1,ipr(1~nporous) 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0. ==pd2,pd1,pe2,pe1(1~nporous)0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0 ==nxo,mxo,nyo,myo(1~nporous) ==d50(1),xporosity(2) [armor layer2] 4 ==nporous(2) 0.,0.,0.,0., 0.,0., 0.,0.,0.,0., 1.,1.,1.,1. ==pa2,pa1,pb2,pb1(1~nporous) 0.,0.,0.,0., , 1,0,0,0 ==pc2,pc1,ipr(1~nporous) 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0. ==pd2,pd1,pe2,pe1(1~nporous) 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0 ==nxo,mxo,nyo,myo(1~nporous) ==d50(3),xporosity(3) [armor layer3] 4 ==nporous(3) 0.,0.,0.,0., .04,.6667,.0, 0.,0.,0.,0., 1.,1.,1.,1. ==pa2,pa1,pb2,pb1(1~nporous) 0.,0.,0.,0., , 1,0,0,0 ==pc2,pc1,ipr(1~nporous) 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0. ==pd2,pd1,pe2,pe1(1~nporous) 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0 ==nxo,mxo,nyo,myo(1~nporous) ==d50(4),xporosity(4) [front filter layer1] 4 ==nporous(4) 0.,0.,0.,0., 0.,.7034,1., 0.,0.,0.,0., 1.,1.,1.,0. ==pa2,pa1,pb2,pb1(1~nporous) 0.,0.,0.,0., , 1,0,0,0 ==pc2,pc1,ipr(1~nporous) 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0., 0.,0.,0.,0. ==pd2,pd1,pe2,pe1(1~nporous) 0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0 ==nxo,mxo,nyo,myo(1~nporous) ==d50(5),xporosity(5) [front filter layer2] 4 ==nporous(5) 0.,0.,0.,0., .04,.7034,.0, 0.,0.,0.,0., 1.,1.,1.,1.==pa2,pa1,pb2,pb1(1~nporous) 0.,0.,0.,0., ,