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