孤立波作用下邊界層的模擬與分析

張 健，劉昭偉，包洪福

（1.河海大學(xué) 水文水資源與水利工程科學(xué)國家重點(diǎn)實驗室，江蘇 南京 210098；

2.清華大學(xué) 水利水電工程系，北京 100084；3.中國長江三峽集團(tuán)公司，北京 100038）

孤立波作用下邊界層的模擬與分析

張 健1，2，劉昭偉1，2，包洪福2，3

（1.河海大學(xué) 水文水資源與水利工程科學(xué)國家重點(diǎn)實驗室，江蘇 南京 210098；

2.清華大學(xué) 水利水電工程系，北京 100084；3.中國長江三峽集團(tuán)公司，北京 100038）

孤立波作用下的邊界層內(nèi)的剪切力以及渦量變化對海嘯傳播和海底地形塑造十分重要。本文基于多區(qū)域譜方法，利用直接模擬（DNS）數(shù)值模型，對在具有矩形斷面的U形水洞內(nèi)的孤立波下的邊界層流動進(jìn)行了模擬。將數(shù)值模擬結(jié)果與解析解以及試驗結(jié)果進(jìn)行了對比，發(fā)現(xiàn)數(shù)值結(jié)果與后兩者吻合得較好。模擬結(jié)果顯示，在低雷諾數(shù)下，擾動不會改變流態(tài)，而隨著雷諾數(shù)的增大，流態(tài)會變得十分復(fù)雜。中等雷諾數(shù)情況下，邊界層內(nèi)會產(chǎn)生正向的渦，并進(jìn)行穩(wěn)定的傳播。在較高雷諾數(shù)情況下，流動進(jìn)入層流向紊流發(fā)展的過渡期，此時邊界層內(nèi)會產(chǎn)生正渦以及負(fù)渦，并會在水深方向進(jìn)行不規(guī)則運(yùn)動。

孤立波；邊界層；剪切力；渦管

1 研究背景

海嘯波進(jìn)入大陸架淺水區(qū)時，會受到水深急劇變化的影響，從而能量瞬時集中，波高驟然增大。在這一過程中非線性作用突顯，傳統(tǒng)線性Stokes波浪理論已不再適用，需要采用孤立波模型對其進(jìn)行研究。在高雷諾數(shù)假設(shè)下，遠(yuǎn)離水底的水體內(nèi)黏性作用很小，可忽略不計，因而可以利用有勢流動的解析方法進(jìn)行分析。但海嘯波的長時間、長距離傳輸過程則會對其進(jìn)行精確預(yù)測產(chǎn)生非常明顯的負(fù)面作用，如傳播過程中的阻力影響難以估算、邊界層分離影響程度等等。因而需要利用孤立波模型對這些問題進(jìn)行進(jìn)一步的探究。

通過假設(shè)孤立波作用下邊界層內(nèi)流動為層流狀態(tài)，可以對這一問題進(jìn)行理論分析。例如，Liu and Orfila利用勢流理論推導(dǎo)了瞬時長波傳播時，邊界層內(nèi)的流速分布以及剪切力分布的解析解。［1］Liu等通過迭代的方法對邊界層內(nèi)非線性的流動方程進(jìn)行求解，發(fā)現(xiàn)非線性方程結(jié)果與線性方程結(jié)果相差很小，得到了在一定情況下可以只考慮線性項在邊界層內(nèi)的影響的結(jié)論。［2］

除理論分析外，數(shù)值模擬也是對邊界層進(jìn)行探究的有效手段。Vittoriand Blondeaux利用RANS方法對在自由表面的孤立波作用下的底部邊界層進(jìn)行模擬。結(jié)果顯示在減速期、波高大于某一臨界值情況下會出現(xiàn)紊動現(xiàn)象，該臨界值與邊界層厚度與當(dāng)?shù)厮钪扔嘘P(guān)，而且邊界層內(nèi)剪切系數(shù)與多種參數(shù)有關(guān)［3］。Suntoyo and Tanaka同樣利用RANS方法對孤立波邊界層問題進(jìn)行探究，采用BSL k-ω紊流模型計算線性化后的邊界層方程。模擬結(jié)果顯示隨著雷諾數(shù)的增加，紊動能以及底面剪切力會隨之增加，但反向剪切力的峰值卻會減小。在完全發(fā)展的紊流狀態(tài)下，邊界層內(nèi)剪切力與邊界層外紊動擴(kuò)散所產(chǎn)生的速度同時產(chǎn)生，而且反向剪切力峰值同樣減小，這一結(jié)論對近岸泥沙輸移影響巨大［4］。

利用試驗可以對理論分析和數(shù)值計算進(jìn)行驗證。Liu等在波浪水槽中進(jìn)行了一系列孤立波下邊界層速度場測量試驗，試驗中孤立波由活塞式造波機(jī)產(chǎn)生，并利用PIV對速度場進(jìn)行測量。試驗條件下邊界層為層流狀態(tài)。試驗結(jié)果與Liu and Orfila和Liu etal對層流流場的分析結(jié)論相一致。［2］

僅利用活塞制造的波動很難在底部形成紊動，為對邊界層內(nèi)的相干結(jié)構(gòu)進(jìn)行探究，Sumer et al在具有矩形斷面的U形水洞（water tunnel）中進(jìn)行了一系列探究性的試驗。試驗利用氣壓系統(tǒng)模擬制造孤立波運(yùn)動形態(tài)，利用水洞是U形的這一特殊構(gòu)造，模擬孤立波中流體微團(tuán)的速度變化，從而使底部出現(xiàn)不穩(wěn)定狀態(tài)，邊界層內(nèi)不再保持層流［5］。試驗得到了雷諾數(shù)介于9.4×104和2.0×106之間的工況下的邊界層內(nèi)速度和底面剪切力的數(shù)據(jù)。此處雷諾數(shù)的定義為：

式中：U0m為邊界層外自由流動的最大速度，m/s；ν為運(yùn)動黏滯系數(shù)，m2/s；2a為孤立波作用下自由流動區(qū)域中一個周期相位內(nèi)流體質(zhì)點(diǎn)的最大位移的1/2，m。

Sumer et al對試驗現(xiàn)象進(jìn)行分析，發(fā)現(xiàn)當(dāng) Re≤2×105時，流場為層流狀態(tài)；當(dāng) 2×105＜Re≤5×105時，流動基本保持層流狀態(tài)，但當(dāng)壓力梯度變向時會在邊界層內(nèi)出現(xiàn)一系列渦管，渦管的出現(xiàn)，使得底面剪切力分布也變得不再光滑；當(dāng)雷諾數(shù)達(dá)到5×105時，流場由層流開始向紊流轉(zhuǎn)變，進(jìn)入過渡期。當(dāng)雷諾數(shù)達(dá)到試驗條件下的最大值2.0×106情況時，約有半個周期進(jìn)入了紊流狀態(tài)，而開始的半個周期則仍保持層流狀態(tài)。根據(jù)雷諾數(shù)增長引起的邊界層內(nèi)剪切力變化規(guī)律來看，雷諾數(shù)繼續(xù)增加會導(dǎo)致紊動相位持續(xù)的時間持續(xù)增加，甚至遍布幾乎整個周期，即流動進(jìn)入完整的紊動流動狀態(tài)，但現(xiàn)在試驗還無法確定這一位置的具體雷諾數(shù)情況。［5］

有別于Vittori和 Blondeaux以及Suntoyo和 Tanaka的大渦模擬方法，本文采用DNS法，利用Di?amessis etal提出的多區(qū)域譜方法模型［6］，對Liu和Orfila、Liu et al的理論分析以及Sumer et al的試驗分析進(jìn)行驗證，并對不同時期的邊界層內(nèi)流場分布的特性進(jìn)行探究。

2 模型建立

2.1 控制方程二維U型水槽中的流動可以用Navier-Stokes方程進(jìn)行描述：

其中，u=（u，w）為速度矢量；p為壓力項，Pa；ρ為水的密度，kg/m3，假設(shè)水體不可壓縮，則ρ可以取常數(shù)1.0×103kg/m3；ν為運(yùn)動黏滯系數(shù)，取值為1.0×10-6m2/s。

在邊界層內(nèi)由于黏性作用而無法保持有勢運(yùn)動，根據(jù)這一特性，我們可以將邊界層內(nèi)的壓強(qiáng)場分解為兩部分進(jìn)行求解，即有勢流動引起的平均項和有旋運(yùn)動引起的擾動項，這樣壓強(qiáng)場可以表示為：

其中下標(biāo)為w的表示的是平均項；帶有上標(biāo)波浪線的為擾動項。壓強(qiáng)平均項利用Liu和Orfila的解析方法進(jìn)行推導(dǎo)

2.2 失穩(wěn)觸發(fā)通過對非恒定流的振蕩邊界層進(jìn)行測量以及數(shù)值模擬顯示紊動產(chǎn)生的動力學(xué)條件十分復(fù)雜。Costamagna，Vittor和Blondeaux認(rèn)為在振蕩流動的水體中產(chǎn)生紊動現(xiàn)象與穩(wěn)定水體中的紊動現(xiàn)象有相同的機(jī)理［7］。在加速期末段的近壁區(qū)域中會出現(xiàn)低速帶，而在減速期這些低速帶就會扭轉(zhuǎn)、振蕩，最終破碎，產(chǎn)生小型渦。由于黏性作用，這些渦又會在加速期的初段耗盡。因而一般的紊流模型只能顯示不穩(wěn)定邊界層的整體特性，卻不能顯示這一過程的變化特征。

在Sumer etal的U型水洞試驗中，我們可以看到邊界層內(nèi)的渦流管會沿一定距離分布，并同時產(chǎn)生和消失，換言之渦流管的生成和發(fā)展具有一定周期性。在數(shù)值模擬中，如果沒有引入任何擾動的話，流場幾乎不會進(jìn)入失穩(wěn)狀態(tài)。為了使數(shù)值模擬中的流場迅速進(jìn)入失穩(wěn)狀態(tài)，并可以對不穩(wěn)定的變化過程有更清晰的探究，可以在最易失穩(wěn)的相位下引入一個微小的速度擾動分量，此相位下流動速度達(dá)到最大值，壓力梯度開始變向。速度擾動分量的縱向和垂向表達(dá)形式為：

其中，m表示擾動速度的峰值，m值要大到足夠引起流場失穩(wěn)，也要小到不會繼續(xù)對失穩(wěn)后的流場產(chǎn)生影響，因而在本文的模擬中m=0.001；b表示初階不穩(wěn)定的相對強(qiáng)度，Smyth和 Moum認(rèn)為b=0.4177［8］；h0表示邊界層的厚度；z0表示在0°相位位置時邊界層的垂向坐標(biāo)；u0=U0m；k0為沿縱向方向的波數(shù)。

本文所采用的數(shù)值模型在縱向方向上采用切比雪夫多項式進(jìn)行擬合，因而擾動速度最終表達(dá)形式也是以三角級數(shù)的形式展現(xiàn)出來，而且可以式（6）、式（7）滿足連續(xù)方程以及下文的邊界條件。

為了計算h0、z0以及k0三個參數(shù)，需要引入剪切層渦厚度這一參數(shù)，其定義為：

δω雖然會隨相位改變而有所不同，但差異并不大，可以用ωt≈32°的位置計算渦厚度。于是上述參數(shù)的計算公式為：

式中λ表示波長。

2.3 邊界條件與初始條件根據(jù)Sumer et al的試驗，我們將計算域定義在一個Lx×H的矩形區(qū)域內(nèi)，其中H為U型水槽從槽底到中心線的高度；而Lx為沿波流動方向的可變長度，對于不同試驗有所不同，主要由渦流管的個數(shù)和尺寸決定，以便滿足速度和壓力的周期性特性：

計算域的頂部為無滑移、無變形的邊界：

底部為固定邊界：

2.4 數(shù)值方法控制方程利用Diamessis et al提出的多區(qū)域譜方法模型進(jìn)行求解［6］。模型中時間導(dǎo)數(shù)項利用向后分步差分進(jìn)行離散，共分3步交替進(jìn)行，以保證時間計算的準(zhǔn)確性。二維模型中兩個方向具有不同的性質(zhì)，因而也采用不同的離散方法對空間導(dǎo)數(shù)項進(jìn)行離散。在周期性方向（波運(yùn)動方向）采用傅立葉譜法進(jìn)行離散，基函數(shù)選用切比雪夫多項式，模態(tài)階數(shù)根據(jù)不同問題進(jìn)行選?。ū疚挠嬎惆咐胁捎们?階）。在非周期性方向（水深方向），將計算域劃分為多個子區(qū)域，每個區(qū)域的垂向長度不同，為更好的對邊界層進(jìn)行模擬，由下至上區(qū)域高度逐漸提高。計算時，每個子區(qū)域內(nèi)采用勒讓得譜法進(jìn)行離散，模態(tài)階數(shù)同樣由問題精度要求決定（本文計算案例中采用了前4階）。較高雷諾數(shù)以及較薄邊界層的情況下，子區(qū)域間數(shù)值傳遞往往會出現(xiàn)Gibbs振蕩，為避免誤差的繼續(xù)發(fā)展，模型依次采用了罰函數(shù)、譜過濾以及界面平均三種方法來保證計算的穩(wěn)定以及誤差最小化。

本文的模擬都從靜止?fàn)顟B(tài)出發(fā)：

3 算例分析

本文共對4種工況進(jìn)行模擬，分別對應(yīng)三種不同的流場狀態(tài)：無渦層流、有渦層流和過渡期。利用前文所建立的數(shù)值模型，對不同流態(tài)的底面剪切力、渦量分布等進(jìn)行分析，并與理論分析和試驗結(jié)果進(jìn)行比較。表1為4種工況的參數(shù)設(shè)置情況，各個參數(shù)的物理含義如前文所述。

表1 算例參數(shù)

對于具體的算例來說，縱向計算長度為Lx，垂向長度為0.145m。計算時為了能更好地反映邊界層以及水體底部的信息，在垂向上將計算域劃分為11個子區(qū)域，由下至上高度分別為0.01、0.02、0.04、0.08、0.16、、0.32、0.64、1.28、2.56、5.12以及4.27 cm。經(jīng)過網(wǎng)格無關(guān)性分析之后發(fā)現(xiàn)縱向96個節(jié)點(diǎn)，垂向每個子區(qū)域內(nèi)25個節(jié)點(diǎn)就可以對問題進(jìn)行較為精確地模擬，這種劃分方式可以使得效率和精度得到較好平衡。

3.1 無渦層流式（6）、式（7）表示的擾動速度對渦管的迅速產(chǎn)生有重要的影響。但在低雷諾數(shù)情況下，引入擾動是否也會產(chǎn)生渦管，則需要進(jìn)一步的分析。本節(jié)將對表1的工況1進(jìn)行模擬來說明無渦層流的流場分布情況。數(shù)值模擬時網(wǎng)格參數(shù)如上文所述。

低雷諾數(shù)情況下底面剪切力分布如圖1所示。圖中實線部分為數(shù)值計算結(jié)果，虛線部分為Liu和 Orfila、Liu等提出的解析解，兩者對比可以發(fā)現(xiàn)具有較好的吻合度。而且在任何相位位置都未發(fā)生震蕩現(xiàn)象，流動保持穩(wěn)定，直至速度降到0。

圖1 Re=2.56×104時底面剪切力數(shù)值結(jié)果與解析結(jié)果相位對比

圖2為數(shù)值計算得到的3種不同相位下渦量分布等勢圖。從圖中可以看出流場未出現(xiàn)渦管，維持層流狀態(tài)，與剪切力分布相互驗證。而且，隨著相位的發(fā)展，渦量也逐漸減小，直至變?yōu)?，滿足層流特性。

3.2 有渦層流Sumer等根據(jù)試驗結(jié)果將孤立波作用下邊界層內(nèi)的流動分為4個階段，并給出了每個階段具體的雷諾數(shù)分布情況［5］，但實際上兩個階段之間很難有較為明顯的界限，尤其是有渦管層流與過渡區(qū)之間的界限更加不明顯。利用工況2、工況3進(jìn)行對比來分析中等雷諾數(shù)情況下邊界層內(nèi)流動情況。與3.1節(jié)不同的是，進(jìn)入不穩(wěn)定階段后很難進(jìn)行理論求解，因而將數(shù)值計算結(jié)果與試驗進(jìn)行對比。

圖3分別為Re=4.20×105，Re=6.06×105時底面剪切力分布數(shù)值解與試驗結(jié)果對比情況。從圖中可以看出，中等雷諾數(shù)下在加速期過后流動進(jìn)入不穩(wěn)定狀態(tài)，剪切力出現(xiàn)輕微震蕩。由于這種震蕩具有隨機(jī)性，因而試驗結(jié)果與數(shù)值解果略有差距，但二者的趨勢卻是一致的。另一方面，有渦層流和過渡期的界限并不十分明顯，但雷諾數(shù)越大，震蕩幅度越大，延續(xù)的相位也越大。

圖2 Re=2.56×104時3種相位時邊界層內(nèi)渦量分布等勢圖

圖3 中等雷諾數(shù)時底面剪切力數(shù)值解與試驗點(diǎn)相位對比

圖4 Re=4.20×105時3種相位時邊界層內(nèi)渦管分布等勢圖

圖5 Re=6.06×105時3種相位時邊界層內(nèi)渦管分布等勢圖

圖4、圖5分別為Re=4.20×105，Re=6.06×105時邊界層內(nèi)渦管發(fā)展等勢圖。從圖中可以看出，中等雷諾數(shù)情況下，渦管產(chǎn)生后會與波進(jìn)行反向傳播，但渦管中心高度不會改變，并穩(wěn)定的在邊界層內(nèi)發(fā)展。而有渦層流與過渡前期的渦管分布區(qū)別并不大。

3.3 過渡期隨著雷諾數(shù)進(jìn)一步的增加，邊界層內(nèi)流動會發(fā)生巨大的變化。接下來將對層流向紊流轉(zhuǎn)變的過渡期進(jìn)行模擬分析，所選工況為工況4。

圖6為大雷諾數(shù)下底面剪切力分布情況。與中等雷諾數(shù)情況對比可以看到，此時剪切力的震蕩情況更加劇烈，且振幅已經(jīng)超過了穩(wěn)定狀態(tài)下由于壓力梯度引起的最大剪切力。試驗結(jié)果與數(shù)值模擬結(jié)果有相同的趨勢，但相位略有差距。隨著雷諾數(shù)的增加，試驗測量的難度也急劇增大，因而這種情況下的測量結(jié)果也很難與實際情況完全相同，因而可以認(rèn)為數(shù)值結(jié)果在一定程度上是合理的。

圖6 Re=1.22×106時底面剪切力數(shù)值解與試驗點(diǎn)相位對比

圖7為工況的渦管分布情況。從圖中可以看出，大雷諾數(shù)下渦管的運(yùn)動趨勢與低雷諾數(shù)相似，但復(fù)雜程度卻加大了，渦管的分布有逐漸上移的趨勢，這使得邊界層的界限也跟著上移。

圖7 Re=1.22×106時3種相位時邊界層內(nèi)渦管分布等勢圖

4 結(jié)論與討論

本文利用Diamessis etal.提出的多區(qū)域譜方法模型，對孤立波作用下邊界層內(nèi)流動進(jìn)行了分析，并將計算結(jié)果與Liu和Orfila、Liu et al的解析解以及Sumer etal的試驗分析進(jìn)行了對比。通過數(shù)值模擬得到：（1）利用二維模型可以對U型水槽試驗進(jìn)行模擬，模擬結(jié)果也較為理想；（2）低雷諾數(shù)（Re≤2×105）情況下流動不會受到擾動的影響，可以完全保持層流狀態(tài)，邊界層內(nèi)并不產(chǎn)生渦管；（3）中等雷諾數(shù)（2×105＜Re≤5×105）情況下，邊界層內(nèi)會產(chǎn)生渦管，但流動仍可以維持在層流狀態(tài)，底面剪切力會小幅振蕩，但會隨相位逐漸減小，渦產(chǎn)生后會與波進(jìn)行反向運(yùn)動，但渦管中心一般不會改變，流動較為穩(wěn)定；（4）高雷諾數(shù)（Re＞5×105）情況下，流動進(jìn)入過渡狀態(tài)，底面剪切力發(fā)生劇烈振蕩，振幅較大，渦管的產(chǎn)生與發(fā)展出現(xiàn)復(fù)雜情況。

自然狀態(tài)下，流動多處于高雷諾數(shù)情況，而低、中雷諾數(shù)則只在波產(chǎn)生的初期出現(xiàn)，因而關(guān)于高雷諾數(shù)情況下渦管產(chǎn)生、發(fā)展以及相互作用仍需要做進(jìn)一步的探究。

參 考 文 獻(xiàn)：

［1］ Liu P L-F，Orfila A.Viscous effects on transient long-wave propagation［J］.J.Fluid Mech.，2004，520：83-92.

［2］ Liu P L-F，Park Y S，Cowen E A.Boundary layer flow and bed shear stress under a solitary wave［J］.J.Fluid Mech.，2007，574：449-463.

［3］ Blonderaux P，VittoriG.RANSmodelling of the turbulent boundary layer under a solitary wave［J］.Coastal Engi?neering，2012，60：1-10.

［4］ Suntoyo，Tanaka H.Numericalmodeling of boundary layer flow for a solitarywave［J］.J.Hydro-environmentRe?search，2009，3：129-137.

［5］ Sumer B M，Jensen PM，Sorensen L B，et al.Coherent structures in wave boundary layers.Part2.solitarymo?tion［J］.J.Fluid Mech.，2010，646：207-231.

［6］ Diamessis P J，Domaradzki JA，Hesthaven JS.A spectralmultidomain penaltymethod model for the simulation of high Reynolds number localized incompressible stratified turbulence［J］.J.Compu.Phys.，2004，202（1）：298-322.

［7］ Costamagna P，VittoriG，Blondeaux P.Coherent structures in oscillatory boundary layers［J］，J.Fluid Mech.，2003，474：1-33.

［8］ Smyth W D，Moum JN.Length scales of turbulence in stably stratified mixing layers［J］.J.Phys.Fluids.，2000，12（6）：1327-1342.

Simulation and analysis of boundary layer under the solitary wave

ZHANG Jian1，2，LIU Zhaowei1，2，BAO Hongfu2，3
（1.State Key Laboratory of Hydrology-Water Resourcesand Hydraulic Engineering，HohaiUniversity，Nanjing 210098，China；2.Departmentof Hydraulic Engineering，Tsinghua University，Beijing 100084，China；3.China Three GorgesCorporation，Beijing 100038，China）

The shear stress and vortex tube in the boundary layer are very important to tsunami propaga?tion and underwater terrain under the solitary wave.This paper takes advantage of direct numerical simula?tion method（DNS） to simulate the boundary layer under the solitary wave in a U shape water tunnel，based on the spectral multi-domain model.By comparing numerical results with analytical solutions and ex?perimental results，it is found that they agree with each other well.The numerical results show that with a low Reynolds number，the disturbance does not change the flow pattern，and with the increase of Reynolds number，flow pattern will become very complicated.Within a medium Reynolds number，positive vortex tube will appear in the boundary layer，and transport stably.Under the condition of a high Reynolds num?ber，layer flow changes into a transition period of turbulence flow and positive vorticity and negative vortici?ty appear in the boundary layer，and have oscillating development.

solitary wave；boundary layer；shear stress；vortex tube

TV139.2

：Adoi：10.13244/j.cnki.jiwhr.2015.04.003

1672-3031（2015）04-0254-07

（責(zé)任編輯：李福田）

2015-04-30

國家自然科學(xué)基金資助項目（51279079）；“十二五”國家科技支撐計劃（2013BAB05B05，2010BAE00739）

張?。?993-），男，吉林人，碩士生，主要從事水力學(xué)研究。zhang0249@126.com

劉昭偉（1973-），男，河北人，副教授。liuzhw@tsinghua.edu.cn