黃河口黏性非均勻泥沙沉降速度試驗(yàn)研究

陳柏文，拾 兵，王 霞，馬俊峰，王俊杰，巴彥斌

（1.中國(guó)海洋大學(xué) 工程學(xué)院，山東 青島 266100；2.黃河河口管理局 利津黃河河務(wù)局，山東 東營(yíng) 257091）

0 引言

黃河口屬弱潮型河口，黃河流域來(lái)沙量比潮流、波浪的泥沙搬運(yùn)量大很多［1］，約有80%的黃河來(lái)沙絮凝沉降在河口區(qū)域［2］，極易形成淤積型三角洲。黃河口泥沙屬于黏性細(xì)顆粒，具有復(fù)雜特殊的電化學(xué)性質(zhì)，受電化學(xué)作用、化合物作用以及網(wǎng)捕作用的影響，黏性細(xì)顆粒泥沙易在水中形成絮團(tuán)，會(huì)對(duì)河口海岸的變遷、河道沖淤等產(chǎn)生影響［3－4］。除上述作用影響外，張金鳳等［5］還通過(guò)碰撞模擬，認(rèn)為絮凝體的形成在于顆粒間的碰撞，并提出幾種主要顆粒碰撞機(jī)制：布朗運(yùn)動(dòng)、懸浮顆粒不等速沉降、速度梯度或流速剪切。國(guó)內(nèi)外研究已發(fā)現(xiàn)諸如泥沙濃度、泥沙粒徑、鹽度、溫度、分形維數(shù)以及紊動(dòng)剪切強(qiáng)度等影響泥沙沉降的因素［6－11］。

當(dāng)前測(cè)量泥沙沉降速度的方法主要有Mclaughlin法、中值沉速法、Rouse 公式擬合法、基于圖像灰度的沉速試驗(yàn)法、LISST 和OBS 結(jié)合法等［12－13］。隨著學(xué)科不斷交叉，研究人員從不同的角度研究影響泥沙沉速的因素：通過(guò)分形理論與分形維數(shù)在泥沙沉降研究中的應(yīng)用，周銀軍等［14］指出分形維數(shù)能很好地描述泥沙顆粒物理特性、群體泥沙組成、沉積物排列方式以及細(xì)顆粒泥沙絮凝體結(jié)構(gòu)；張宇卓［15］基于介觀尺度的耗散粒子動(dòng)力學(xué)方法，建立了包含力場(chǎng)勢(shì)函數(shù)的均勻各向同性紊流模型，分析了泥沙濃度、電荷量、紊動(dòng)剪切對(duì)泥沙絮凝的影響。但已有研究成果大都不能揭示絮凝后的絮團(tuán)顆粒粒徑的級(jí)配情況，且在絮凝體沉降速度的計(jì)算中采用絮凝體中值粒徑來(lái)代替實(shí)際的絮團(tuán)粒徑，實(shí)際上不同粒徑的絮凝體所占的比例并不相同，計(jì)算絮凝體沉降速度時(shí)采用中值粒徑會(huì)導(dǎo)致結(jié)果發(fā)生偏差。本文應(yīng)用分形理論以及激光粒度分析儀對(duì)斯托克斯公式進(jìn)行完善，采用加權(quán)粒徑代表絮凝體的真實(shí)粒徑，并將完善后的斯托克斯公式應(yīng)用到黃河口的泥沙沉降速度計(jì)算中，以期為更好地認(rèn)識(shí)黏性細(xì)泥沙的絮團(tuán)發(fā)育和沉降過(guò)程提供可靠的依據(jù)。

1 原理與方法

1.1 公式推導(dǎo)

Kranenburg 指出，在一個(gè)粒徑為D的絮凝體中有k個(gè)初始粒徑為Dp的黏性細(xì)泥沙顆粒［16］，則：

式中：nf為絮凝體的分形維數(shù)。

由式（1）可以得到單顆絮凝體的質(zhì)量：

式中：m1為單顆絮凝體的質(zhì)量，α為初始泥沙顆粒球形度系數(shù)，ρs為初始泥沙顆粒密度。

假設(shè)某層水體中含沙量為C，則該層水體中絮凝體質(zhì)量為

式中：l為沉降桶的長(zhǎng)度，d為沉降桶的寬度，h為該層水體的厚度。

由式（1）～式（3）可以得到某層水體中絮凝體個(gè)數(shù)為

式中：β為絮凝體顆粒球形度系數(shù)，ρf為絮凝體密度。

由式（4）～式（5）可得絮凝體密度：

式中：ρw為水的密度，Δρf為絮凝體有效密度。

隨著激光粒度分析儀的應(yīng)用普及，王國(guó)梁等［17］推導(dǎo)出體積分形維數(shù)計(jì)算公式，該公式不需要作不同粒級(jí)土壤具有相同密度這一假設(shè)，所以比質(zhì)量分維模型更具有合理性，其計(jì)算公式為

斯托克斯公式：

式中：ωs為泥沙顆粒的沉降速度，μ為沉降桶中水的動(dòng)力黏滯系數(shù)，g為重力加速度。

應(yīng)用斯托克斯公式時(shí)把球形度系數(shù)α和β作為1處理［18－19］。激光粒度分析儀可以測(cè)出48 組絮凝體粒徑以及它們所占的體積百分含量，且每組絮凝體粒徑所占比例各異［20］，為更好地體現(xiàn)絮凝體的實(shí)際粒徑，嘗試采用加權(quán)粒徑代替中值粒徑來(lái)完善斯托克斯公式。完善后的斯托克斯公式如下：

由式（6）～式（10）可推導(dǎo)出泥沙沉速公式為

式中：Pi、Pj分別為絮凝后和絮凝前第i、j組顆粒粒徑所占的體積百分含量；Di、Dpj分別為絮凝后和絮凝前第i、j組顆粒粒徑。

而加權(quán)粒徑Dw的計(jì)算公式為

令Dw為絮凝體的加權(quán)粒徑，Dpw為初始泥沙顆粒加權(quán)粒徑，則最終的沉速計(jì)算公式為

1.2 公式驗(yàn)證

采用實(shí)測(cè)泥沙沉速資料［21］對(duì)式（13）的合理性以及計(jì)算結(jié)果的精度進(jìn)行驗(yàn)證，采用均方根誤差（RMSE）、平均絕對(duì)誤差（MAE）和平均絕對(duì)百分比誤差（MAPE）進(jìn)行有效性評(píng)定，各項(xiàng)誤差的值越小，表明有效性越高。

由式（13）計(jì)算的沉速與實(shí)測(cè)資料對(duì)比如圖1 所示，計(jì)算值與實(shí)測(cè)值關(guān)系的決定系數(shù)R2為0.943，RMSE、MAE 和MAPE 均在10%以內(nèi)，即計(jì)算值與實(shí)測(cè)值吻合良好，有效性高，說(shuō)明式（13）可用于細(xì)泥沙顆粒沉降速度的計(jì)算。孟若霖等［22］運(yùn)用斯托克斯公式計(jì)算得到上述實(shí)測(cè)資料的泥沙沉降速度（采用中值粒徑計(jì)算），計(jì)算值與實(shí)測(cè)值的決定系數(shù)R2僅為0.782。主要原因是加權(quán)粒徑代替中值粒徑后，絮凝體粒徑體積占比不均的問(wèn)題得以解決，泥沙沉降速度不會(huì)因中值粒徑體積百分含量過(guò)小而出現(xiàn)偏差，加權(quán)粒徑使得結(jié)果代表性更好。因此，采用加權(quán)粒徑來(lái)推求泥沙沉降速度比采用中值粒徑來(lái)推求具有更高的準(zhǔn)確性。

圖1 沉速計(jì)算值與實(shí)測(cè)值對(duì)比

2 采樣地點(diǎn)及測(cè)量手段

本試驗(yàn)采用黃河口的黏性細(xì)顆粒泥沙進(jìn)行以含沙量和鹽度為影響因素的靜水沉降試驗(yàn)研究。試驗(yàn)沙樣取自黃河口三角洲水域。利用激光粒度分析儀測(cè)得試驗(yàn)時(shí)懸浮泥沙中值粒徑為13.28 μm。

沉降試驗(yàn)采用自主設(shè)計(jì)開(kāi)發(fā)的大型溫控自動(dòng)攪拌沉降試驗(yàn)矩形桶以及山東耐克特激光粒度分析儀。沉降桶尺寸為1.1 m×0.3 m×1.5 m（長(zhǎng)×寬×高），內(nèi)置溫度計(jì)、水位計(jì)、OBS3＋濁度儀（量程0～4 000 NTU）、可控型升溫保溫裝置以及攪拌裝置曝氣泵等。試驗(yàn)將整個(gè)水深均勻劃分為6 層，并采用OBS3＋濁度儀和激光粒度分析儀對(duì)不同懸沙濃度（2、5 kg／m3和7 kg／m3）、不同鹽度的水樣進(jìn)行濁度和泥沙顆粒粒徑分布的測(cè)量，分別記錄底層初始濁度、初始泥沙顆粒粒徑、底層濁度為初始濁度50%時(shí)的時(shí)間及對(duì)應(yīng)的絮凝體顆粒粒徑。

3 試驗(yàn)結(jié)果與分析

3.1 OBS3＋的參數(shù)率定與使用

本次試驗(yàn)絮凝泥沙顆粒取樣采用陳曦［23］的中值沉速法，即底層含沙量為初始含沙量的50%時(shí)進(jìn)行取樣。試驗(yàn)前對(duì)OBS3＋濁度儀進(jìn)行參數(shù)率定，濁度與含沙量的關(guān)系如圖2 所示，擬合曲線決定系數(shù)為0.994。

圖2 OBS3＋測(cè)量濁度值與含沙量的關(guān)系

根據(jù)濁度與含沙量的關(guān)系以及OBS3＋記錄的濁度，可得沉降桶內(nèi)隨時(shí)間變化的泥沙相對(duì)濃度，為使泥沙相對(duì)濃度隨時(shí)間的變化過(guò)程更加直觀，以含沙量為7 kg／m3的泥沙沉降過(guò)程為例，利用數(shù)學(xué)模型進(jìn)行曲線擬合，如圖3 所示。

圖3 泥沙相對(duì)濃度及曲線模型擬合

1）鹽度為0～11‰時(shí)，泥沙相對(duì)濃度隨時(shí)間變化可用反比例函數(shù)來(lái)描述，但低鹽度條件下泥沙相對(duì)濃度變化較為緩慢，且隨著鹽度增大泥沙相對(duì)濃度變化加快，因此低鹽度下泥沙相對(duì)濃度表達(dá)式為

式中：SRC為泥沙相對(duì)濃度，S為鹽度，t為泥沙沉降時(shí)間，tmax、tmin分別為泥沙絮凝沉降的最大時(shí)間、最小時(shí)間。

擬合曲線決定系數(shù)R2＝0.996。

2）鹽度為11‰～25‰時(shí)，中高鹽度下泥沙相對(duì)濃度隨時(shí)間變化相對(duì)較快，隨著鹽度增大泥沙相對(duì)濃度變化會(huì)減緩，但其過(guò)程同樣可用反比例函數(shù)來(lái)進(jìn)行描述：

擬合曲線決定系數(shù)R2＝0.984。

利用數(shù)學(xué)模型，可得出在確定溫度和含沙量條件下不同鹽度和時(shí)間對(duì)應(yīng)的準(zhǔn)確SRC值，并預(yù)測(cè)泥沙相對(duì)濃度隨時(shí)間的變化趨勢(shì)。圖3 表明，同一沉降桶中不同鹽度的水樣底層泥沙相對(duì)濃度隨時(shí)間逐漸減小，且隨著鹽度的增加底層泥沙相對(duì)濃度達(dá)到50%的時(shí)間呈現(xiàn)先減小后增大的趨勢(shì)，達(dá)到泥沙相對(duì)濃度為50%的時(shí)間基本為20～60 min，在鹽度11‰下時(shí)間達(dá)到最小值。

根據(jù)數(shù)學(xué)模型表達(dá)式可得到不同鹽度下底層泥沙相對(duì)濃度達(dá)到50%的時(shí)間，運(yùn)用前文所提出的方法，可以得到相對(duì)應(yīng)時(shí)間下絮凝體顆粒粒徑分布、加權(quán)粒徑以及分形維數(shù)，進(jìn)而利用式（13），可以求得不同含沙量和鹽度下絮凝體的沉降速度。

3.2 鹽度與含沙量對(duì)絮凝體沉降速度和粒徑的影響

1985 年，Azzalini［24］首次提出比正態(tài)分布適用范圍更廣的偏正態(tài)分布，其隨機(jī)變量x的密度函數(shù)如下：

式中：φ、Φ分別為標(biāo)準(zhǔn)正態(tài)分布的密度函數(shù)和分布函數(shù)，η、σ分別為正態(tài)分布的位置參數(shù)和形狀參數(shù)，λ為偏態(tài)參數(shù)，當(dāng)λ＞0 時(shí)右偏。

為系統(tǒng)描述以及預(yù)測(cè)泥沙絮凝沉降速度隨含沙量和鹽度的變化趨勢(shì)，采用考慮鹽度和含沙量影響因子的偏正態(tài)模型對(duì)低含沙量的絮凝體沉降速度與鹽度關(guān)系進(jìn)行參數(shù)以及曲線擬合（ωs為絮凝體沉降速度，取其數(shù)值部分并記為ω），則帶參數(shù)的數(shù)學(xué)模型表達(dá)式為

式中：ω為絮凝體沉降速度數(shù)值，C為含沙量，S為鹽度，k、b、d為待定系數(shù)，t為積分變量。

基于式（17）的擬合結(jié)果如圖4 所示。

圖4 低含沙量下鹽度與沉降速度關(guān)系

在確定初始粒徑的條件下，考慮鹽度和含沙量影響因子的絮凝體沉降速度數(shù)學(xué)模型表達(dá)式為

擬合曲線決定系數(shù)R2＝0.905，λ值取6。

而由試驗(yàn)數(shù)據(jù)以及以往研究表明的鹽度與絮凝體粒徑關(guān)系可知［6，8］，偏正態(tài)分布模型仍然可用于描述鹽度對(duì)絮凝體粒徑的影響過(guò)程。前文指出Dw為泥沙絮凝體加權(quán)粒徑，現(xiàn)取其數(shù)值部分并記為Df，因此在確定初始粒徑的條件下由式（17）可得其擬合曲線圖（見(jiàn)圖5）。

圖5 鹽度與絮凝體粒徑關(guān)系

其數(shù)學(xué)模型表達(dá)式如下：

式中：Df為泥沙粒徑數(shù)值。

擬合曲線決定系數(shù)R2＝0.947，λ值取3.5。

依據(jù)偏態(tài)分布模型的擬合曲線，可以看到部分計(jì)算值高于擬合點(diǎn)，原因是本試驗(yàn)尚未完全考慮影響泥沙絮凝沉降的全部因素，且斯托克斯公式中把球形度系數(shù)近似取為1 不盡合理。為驗(yàn)證模型的適用性，現(xiàn)將偏正態(tài)模型用陳曦試驗(yàn)數(shù)據(jù)進(jìn)行曲線擬合，結(jié)果如圖6 所示。

圖6 陳曦試驗(yàn)數(shù)據(jù)對(duì)偏正態(tài)模型的驗(yàn)證結(jié)果

由圖4～圖6 可知，本試驗(yàn)得到的絮凝體沉降速度為0.11～0.45 mm／s，加權(quán)粒徑為21.1～46.5 μm。在2、5、7 kg／m3含沙量的沉降桶中，隨著鹽度的增加，絮凝體沉降速度和粒徑均呈現(xiàn)先增大后減小的趨勢(shì)，從鹽度0 到鹽度11‰，沉降速度和粒徑逐漸增大，從鹽度11‰到鹽度25‰，沉降速度和粒徑逐漸減小。在含沙量7 kg／m3、鹽度11‰處取得最大沉降速度0.45 mm／s，對(duì)應(yīng)最大加權(quán)粒徑46.5 μm，在含沙量2 kg／m3、鹽度0 處達(dá)到最小沉速0.11 mm／s，對(duì)應(yīng)最小加權(quán)粒徑21.1 μm，達(dá)到泥沙濃度相對(duì)含量50%的時(shí)間基本為20～60 min。這些試驗(yàn)結(jié)果與陳曦試驗(yàn)的結(jié)果基本一致，且偏態(tài)分布模型與陳曦試驗(yàn)數(shù)值擬合良好，擬合曲線決定系數(shù)R2＝0.944，說(shuō)明模型準(zhǔn)確度較高，適用于其他河口或試驗(yàn)的絮凝沉降速度研究。

泥沙顆粒表面帶負(fù)電荷，在靜電作用下呈現(xiàn)獨(dú)特的雙電子層結(jié)構(gòu)，而鹽度的增大使得沉降桶中陽(yáng)離子的濃度增加，雙電子層厚度減小，從而使得泥沙顆粒間的靜電斥力作用減弱，在同等范德華力作用下，泥沙顆粒相互吸引碰撞形成粒徑更大的絮凝體，沉降速度隨之增大。而鹽度的繼續(xù)增大則會(huì)導(dǎo)致電勢(shì)逆轉(zhuǎn)，顆粒間斥力不斷增大，抑制泥沙絮凝體的發(fā)育，沉降速度也隨之減小。而在低含沙量條件下，含沙量的增加提高了泥沙顆粒的碰撞概率，有利于發(fā)育形成更大的絮凝體，沉降速度也因此增大。李學(xué)凱等［25］指出隨著硫酸鉀濃度的不斷增大，沉降速度也不斷增大，在硫酸鉀濃度增加到60 mmol／L（換算鹽度為10.45‰）時(shí)，達(dá)到最大沉速，期間沉降速度從0.23 mm／s 增加到1.42 mm／s。與本文結(jié)果相比可知，兩者沉降速度變化趨勢(shì)基本相同，只是在整體上絮凝體的沉降速度本文結(jié)果略小一些。王尚毅［26］則提出靜水條件下絮凝體沉降速度隨著鹽度的增大呈現(xiàn)出先增大、到達(dá)一定鹽度后再減小的趨勢(shì)，沉降速度為0.13～0.35 mm／s。王家生等［27］研究陽(yáng)離子濃度對(duì)黃河泥沙沉降速度的影響，提出絮凝體沉降速度在0.25 mm／s 附近。這些結(jié)論均印證了本試驗(yàn)結(jié)果的正確性。

3.3 黃河口泥沙最佳絮凝鹽度

黃河三角洲河口泥沙的沉降、淤積、輸移與擴(kuò)散對(duì)黃河河口的變遷與下游河道的安全具有重要意義，而黃河口泥沙最佳絮凝鹽度則是河口整治的重要一環(huán)。由圖7 可知，隨著鹽度的增大，絮凝體粒徑和絮凝體沉降速度變化趨勢(shì)保持一致，即鹽度在0～11‰之間時(shí)，絮凝體粒徑和沉降速度都隨著鹽度的增大而增大，鹽度在11‰～25‰之間時(shí)，絮凝體粒徑以及沉降速度都呈現(xiàn)減小趨勢(shì)，且不同含沙量條件下泥沙絮凝體的粒徑、沉降速度均在鹽度為11‰處達(dá)到最大值，故可認(rèn)為泥沙絮凝的最佳鹽度為11‰。馮曦等［28］在對(duì)黃河口泥沙特性的研究中提出在不同鹽度下黃河泥沙靜水沉降的最佳絮凝鹽度為15‰左右。李芳等［29］在對(duì)高嶺土的靜水沉降試驗(yàn)中指出其最佳絮凝鹽度為13‰左右。李學(xué)凱等［25］指出黃河泥沙最佳絮凝鹽度為10.45‰左右。

圖7 絮凝體粒徑以及沉降速度與鹽度的關(guān)系

4 結(jié)論

1）采用更具有整體代表性的絮凝體加權(quán)粒徑代替中值粒徑進(jìn)行沉速的推求，進(jìn)一步擴(kuò)展了斯托克斯沉速計(jì)算公式的顆粒適用范圍。加之所推導(dǎo)沉速公式只涉及初始及絮凝后粒徑分布，故不僅適用于黃河口，而且適用于其他河口的泥沙沉速計(jì)算。

2）偏正態(tài)分布模型在描述泥沙絮凝沉降過(guò)程、分析預(yù)測(cè)沉降速度、分析絮凝體粒徑隨鹽度和含沙量的變化趨勢(shì)等方面具有良好的適用性，可以考慮引入更多數(shù)學(xué)模型來(lái)研究泥沙沉降規(guī)律。

3）在含沙量為2、5 kg／m3和7 kg／m3的條件下，黃河泥沙絮凝沉降速度隨鹽度增大呈現(xiàn)先增大，達(dá)到最佳絮凝鹽度后減小，且隨著鹽度繼續(xù)增大呈現(xiàn)出整體變化不大的趨勢(shì)；絮凝的最佳鹽度為11‰；沉降速度為0.11～0.45 mm／s；加權(quán)粒徑為21.1～46.5 μm。