李世森，張 睿（天津大學(xué)，天津 300072）

?

泥沙顆粒沿水深方向的受力研究

李世森，張 睿
（天津大學(xué)，天津 300072）

摘要：對(duì)于研究泥沙群體性運(yùn)動(dòng)和實(shí)際工程的人來(lái)說(shuō)，單顆粒泥沙的運(yùn)動(dòng)一直是研究問(wèn)題的基礎(chǔ)和關(guān)鍵。天然泥沙顆粒幾乎都是非均勻的，其運(yùn)動(dòng)形式也具有隨機(jī)性，這就造成顆粒附近流場(chǎng)的隨機(jī)性。國(guó)內(nèi)外學(xué)者對(duì)單顆粒泥沙受力情況的研究由來(lái)已久，文章采用Flow3D軟件對(duì)單顆粒泥沙在定常流下的受力進(jìn)行了數(shù)值模擬，得出了拖曳力系數(shù)和上舉力系數(shù)關(guān)于泥沙顆粒直徑、流速等因素的一定規(guī)律。

關(guān)鍵詞：?jiǎn)晤w粒泥沙；托曳力系數(shù)；上舉力系數(shù)

引 言

河床表面的沉淀物顆粒受水流力的影響，水流力可以分為上舉力和拖曳力。Shields（1936）研究床面泥沙顆粒力的平衡狀態(tài)，進(jìn)而推導(dǎo)出了無(wú)粘性均勻顆粒的起動(dòng)拖曳力，并根據(jù)實(shí)驗(yàn)數(shù)據(jù)畫(huà)出了起動(dòng)拖曳力與雷諾數(shù)的關(guān)系圖，及希爾茲曲線［1］。而后，Einstein和El-Samni（1949）［2］實(shí)驗(yàn)測(cè)量了顆粒由上下壓力差產(chǎn)生的上舉力。經(jīng)過(guò)多年的研究，發(fā)現(xiàn)在一定范圍內(nèi)，拖曳力及上舉力之間存在一定比例，結(jié)果受雷諾數(shù)等的影響，但仍未得出統(tǒng)一公式。

泥沙顆粒在水流中受拖曳力與上舉力，拖曳力FD與上舉力FL的一般表達(dá)式為：

式中：CD和 CL為拖曳力系數(shù)和上舉力系數(shù)；u0為作用在床面泥沙顆粒上的流速。拖曳力系數(shù)與上舉力系數(shù)由實(shí)驗(yàn)得出，具體不同流速下的實(shí)驗(yàn)已很多。

拖曳力系數(shù)是用來(lái)表示在流體中物體受牽引或受阻礙的程度的無(wú)量綱指數(shù)，其應(yīng)用于拖曳力方程中，數(shù)值相對(duì)較小表示物體受拖曳力小，反之表示受拖曳力大。從物理模型實(shí)驗(yàn)到數(shù)值模擬，拖曳力系數(shù)已經(jīng)得到了廣泛大量的研究，上舉力系數(shù)較少，一般認(rèn)為兩系數(shù)與顆粒形狀、表面粗糙度、面積和水流條件等有關(guān)。然而，因?yàn)樾螤?、水流等條件的變化對(duì)其影響過(guò)于復(fù)雜，仍不能有一個(gè)適用于任意雷諾數(shù)的統(tǒng)一的解析公式，當(dāng)雷諾數(shù)較大甚至趨近于無(wú)窮時(shí)，各家結(jié)果差異較大，大量的實(shí)驗(yàn)數(shù)據(jù)分別形成了球形和非球形顆粒的拖曳力、上舉力系數(shù)與雷諾數(shù)的近似曲線。Einstein和El-Sammi （1949）［2］測(cè)量了直接作用在塑料球上的壓力差，以此作為上舉力，并發(fā)現(xiàn)上舉力系數(shù)CL對(duì)不同水流流速（測(cè)量點(diǎn)為球心所對(duì)上部0.15D處）保持一個(gè)定值。Davies和Samad（1978）［5］發(fā)現(xiàn)在邊界處顆粒下層縫隙流動(dòng)強(qiáng)烈時(shí)，當(dāng)Re＜5顆粒所受上舉力為負(fù)值，當(dāng)R6＜5上舉力為正值。Brayshaw（1983）［3］在Re*=52 000時(shí)，測(cè)得直徑為16 mm的顆粒所受上舉力與拖曳力之比為1.8。不同量級(jí)的Re*條件下，顆粒拖曳力和上舉力比值和變化范圍的研究還存在不足。McLauglin（1991）［6］、Mei（1992）［7］、Kurose和Konori（1999）［8］都測(cè)量了球形顆粒上舉力系數(shù)，做出了CL與Re*的關(guān)系曲線，其范圍主要在0＜Re*＜103，隨Re*的增大CL逐漸變小，從他們的統(tǒng)計(jì)數(shù)據(jù)看，CL的變化范圍在0～4之間，有可能出現(xiàn)負(fù)值。

上述國(guó)內(nèi)外學(xué)者進(jìn)行的各種數(shù)值、物理模型實(shí)驗(yàn)和理論研究，揭示了泥沙拖曳力和上舉力的一定關(guān)系，對(duì)泥沙輸移問(wèn)題的研究做出巨大貢獻(xiàn)。但拖曳力系數(shù)和上舉力系數(shù)在大雷諾數(shù)條件下的可靠數(shù)據(jù)仍不夠充足，影響因素研究有待深入。本文基于此背景，主要研究?jī)?nèi)容如下：利用Flow3D軟件，在不同流速、不同高度、不同粒徑條件下，進(jìn)行單顆粒球形泥沙在紊流定常流中的模擬，得到流速和壓力數(shù)據(jù)，進(jìn)而得到泥沙顆粒受力和受扭的結(jié)果，并分析拖曳力系數(shù)和上舉力系數(shù)，討論大雷諾數(shù)下泥沙顆粒的受力情況和流場(chǎng)流態(tài)。

1 數(shù)學(xué)模型的建立與分析

1.1 基本思路

天然河床上的泥沙相互遮蔽，單顆粒受力泥沙較少，但單顆粒泥沙受力情況簡(jiǎn)單經(jīng)典，天然泥沙形狀各異，一般都為不規(guī)則顆粒，研究起來(lái)難度較大，故本文從球形單顆粒泥沙入手。不考慮底床坡度的情況下，本文模擬水平平面底上單顆粒球形泥沙固定（不移動(dòng)、不旋轉(zhuǎn)）在定常流中的受力情況。利用Flow3D軟件建立三維模型，河床一頭入口：x=0處設(shè)一簡(jiǎn)單牛頓流體，顆粒固定于流場(chǎng)穩(wěn)定處，軟件計(jì)算流場(chǎng)各點(diǎn)流速、壓力，從而得出泥沙受力情況。模型采用嵌套式長(zhǎng)方體（或立方體）網(wǎng)格，大網(wǎng)格模擬整個(gè)流場(chǎng)，小網(wǎng)格加密球體周圍區(qū)域，精細(xì)刻畫(huà)顆粒受力。三維立體水槽尺寸為24 m×1 m ×0.7 m，顆粒固定在13 m處（顆粒的位置經(jīng)過(guò)建模計(jì)算檢驗(yàn)，流體進(jìn)入流場(chǎng)約10～11 m發(fā)展充分，流態(tài)穩(wěn)定），流場(chǎng)進(jìn)口邊界條件為流量邊界，出口為連續(xù)性邊界，紊動(dòng)模型采用RNG模型。模擬實(shí)驗(yàn)從流速、顆粒直徑等因素的對(duì)比出發(fā)，進(jìn)行建模和對(duì)照分析。

圖1 網(wǎng)格示意（局部）

計(jì)算初始水深設(shè)定為0.6 m，大網(wǎng)格上表面邊界設(shè)置為壓力邊界，取Fluid fraction為0，即液體上表面為空氣的自由表面。沿水深方向的Z軸，網(wǎng)格設(shè)定為非均勻漸變網(wǎng)格，在［0，0.1］和［0.55，0.65］這兩個(gè)區(qū)域內(nèi)手動(dòng)加密，保重沿Z向網(wǎng)格間距0.005 m，加密水流流速變化劇烈的水槽底部和浮動(dòng)變化的液面區(qū)域。在保證有效計(jì)算的前提下，第一組模擬：顆粒直徑為20 cm，進(jìn)口流速為0.2 m/s，顆粒布置在水槽底至0.4 m內(nèi)各個(gè)位置，各個(gè)位置間隔為0.05 m，模擬不同位置下在流場(chǎng)中的受力，目的為研究顆粒在流體上部的情況；第二組模擬：顆粒直徑為10 cm，其實(shí)情況與第一組相同，目的在于與第一組模擬在粒徑這一因素上對(duì)比；第三組模擬：顆粒直徑為20 cm，顆粒始終布置在水槽底部，改變流速分別為0.1 m/s、0.2 m/s、0.3 m/s、0.4 m/s、0.5 m/s、0.6 m/s、0.7 m/s、1.0 m/s，目的為研究不同流速對(duì)顆粒受力的影響。模擬實(shí)驗(yàn)中，顆粒雷諾數(shù)范圍為2×104～2×105，屬于紊流，且顆粒雷諾數(shù)較大。

1.2 模型準(zhǔn)確性和可行性研究

1.2.1 顆粒半徑校正

必須注意，軟件的網(wǎng)格劃分為長(zhǎng)方體（或立方體）網(wǎng)格，在描述球形顆粒時(shí)存在一定誤差，不能完全精細(xì)描繪球體，故需要進(jìn)行校準(zhǔn)。Ladd（1994）［9］針對(duì)這一情況提出了“水力半徑”這一概念進(jìn)行修正，其方法是將一個(gè)已知半徑ag的球體放置在蠕動(dòng)流中計(jì)算其拖曳力，然后根據(jù)已知的Sangani和Acrivos（1982）［10］修正的拖曳力公式反推出水力半徑a。本文采用Flow3D軟件刻畫(huà)球體時(shí)也遇到校準(zhǔn)問(wèn)題，且后處理工作用自己編制的Fortran程序計(jì)算，故綜合考慮，采用球形顆粒沉降實(shí)驗(yàn)來(lái)校準(zhǔn)顆粒計(jì)算半徑，具體方法如下：用Flow3D軟件模擬顆粒在豎直水槽中自由沉降的情況，模擬中以球體固定，水相對(duì)于球體運(yùn)動(dòng)來(lái)代替顆粒運(yùn)動(dòng)，模型進(jìn)口和出口邊界采用相同的速度邊界，保證整個(gè)流場(chǎng)的速度保持相同值，模擬顆粒勻速運(yùn)動(dòng)的情況，計(jì)算四組：第一組水流方向?yàn)閆軸正向，重力加速度取9.8；第二組水流方向?yàn)閆軸正向，重力加速度取0；第三組水流方向?yàn)閆軸負(fù)向，重力加速度取9.8；第四組水流方向?yàn)閆軸負(fù)向，重力加速度取0。計(jì)算結(jié)果發(fā)現(xiàn)，項(xiàng)目1的（Pz+TAOX）與項(xiàng)目2的（Pz+TAOX）相減得到顆粒所受浮力等于2.835378E+01，項(xiàng)目3與項(xiàng)目4相減得到受浮力等于2.835286E+01，取兩浮力值的平均值，反推得到直徑為1.767908E-01 m。觀察Flow3D模擬的結(jié)果文件，以迎水面流速為部分作為計(jì)算半徑，結(jié)果為0.17 m。對(duì)于球與水的接觸面的刻畫(huà)，球體越大、網(wǎng)格越細(xì)越精確，故根據(jù)網(wǎng)格尺寸，結(jié)合模擬結(jié)果文件，可以得到各組實(shí)驗(yàn)的計(jì)算直徑，以此直徑進(jìn)行后處理。

1.2.2 模型驗(yàn)證

拖曳力系數(shù)的研究不僅限于泥沙顆粒起動(dòng)受力，在顆粒自由沉降中，阻力系數(shù)相當(dāng)于拖曳力系數(shù)，統(tǒng)稱為drag coefficient。為了驗(yàn)證模型和后處理工作的正確性，同時(shí)避免由于網(wǎng)格刻畫(huà)不準(zhǔn)確帶來(lái)的數(shù)據(jù)差別，建立20 cm×20 cm×20 cm的立方體沉降模型，求解立方體受力，得出立方體顆粒阻力系數(shù)CD，與前人實(shí)驗(yàn)數(shù)據(jù)和經(jīng)驗(yàn)值進(jìn)行比較。

與前文球體顆粒沉降模型類似，用Flow3D軟件模擬立方體顆粒在水槽中的自由沉降，固定立方體，水相對(duì)于立方體運(yùn)動(dòng)來(lái)代替顆粒運(yùn)動(dòng)，模型進(jìn)口和出口邊界采用相同的速度邊界。模擬兩組實(shí)驗(yàn)：第一組水槽豎向垂直，水流方向?yàn)閆軸正向；第二組水槽水平放置，水流方向沿X軸正向。模擬流場(chǎng)的雷諾數(shù)為Re=1×105，通過(guò)計(jì)算，求得：第一組FD=4.937 326 N，CD=0.987；第二組FD=4.585 866 N，CD=0.917。

國(guó)內(nèi)外學(xué)者針對(duì)球形和非球形顆粒（如立方體等）的自由沉降（或自由降落）做了大量研究。阻力系數(shù)CD是流場(chǎng)和顆粒形狀構(gòu)成的復(fù)雜函數(shù)，Pettyjohn和Christiansen（1948）［11］以及Hartman和Coughlin（1993）［12］都將CD表示為CD=CD（Re，ψ），其中ψ為顆粒形狀系數(shù)。

B.Krueger、S.Wirtz、V.Scherer（2015）［13］使用數(shù)碼相機(jī)立體成像功能對(duì)顆粒自由沉降進(jìn)行拍攝，從而得到不同形狀顆粒的阻力系數(shù)， 4組立方體（cube）阻力系數(shù)分別為：1.033（Re≈2 000）、0.950（Re≈4 000）、1.018（Re≈4 000）、0.966 （Re∈8 000，10 000）。

龐啟秀、徐金環(huán)、辛海霞（2006）［14］進(jìn)行水槽實(shí)驗(yàn)研究了塊體的水流拖曳力，以長(zhǎng)方體為目標(biāo)，研究了不同長(zhǎng)高比下的拖曳力系數(shù)。實(shí)驗(yàn)水槽長(zhǎng)20 m，寬0.3 m，高0.4 m，實(shí)驗(yàn)塊體尺度長(zhǎng)寬高各組變化。實(shí)驗(yàn)測(cè)量了以上19種塊體尺度在平均流速0.3 m/s、0.6 m/s、0.8 m/s、1.2 m/s情況下的顆粒受力，可知雷諾數(shù)與數(shù)值實(shí)驗(yàn)在同一量級(jí)上，拖曳力系數(shù)CD和形狀參數(shù)（bs/L2）1/3的關(guān)系如圖2。從圖中散點(diǎn)可知，形狀參數(shù)（bs/L2）1/3在立方體時(shí)等于1，圖中CD值在0.818和1附近。

綜上，數(shù)值模擬數(shù)據(jù)與前人實(shí)驗(yàn)數(shù)據(jù)吻合，本文數(shù)值方法可行。

圖2 龐啟秀等水槽試驗(yàn)結(jié)果

2 數(shù)值模擬結(jié)果分析

2.1 水流流態(tài)

以0.2 m/s入口流速流體為例，此情況下，顆粒球心坐標(biāo)為（13，0.5，0.1），直徑為20 cm，顆粒雷諾數(shù)，屬于紊流。從XZ面看，球形顆粒在來(lái)流沖擊下，顆粒上部頂點(diǎn)附近存在一個(gè)流速較大的區(qū)域，顆粒近底部流速較小，水流到達(dá)球面后繞過(guò)球體運(yùn)動(dòng)，在水流方向上球后部有回流現(xiàn)象，且后部流速變化較前部劇烈。分別觀察球形顆粒沿Z向上移0.05 m、0.10 m、0.15 m距離（顆粒仍固定不動(dòng)）后XZ面顆粒局部流速等值線圖會(huì)發(fā)現(xiàn)，與顆粒放置于床底對(duì)比發(fā)現(xiàn)，顆粒上移過(guò)程中，球下部頂點(diǎn)附近發(fā)展出一流速較大的區(qū)域，此區(qū)域流速與上部流速較大區(qū)域的流速接近，與周圍流速值差別較大。隨顆粒繼續(xù)上升，此區(qū)域發(fā)展越來(lái)越大，當(dāng)顆粒不受床面影響時(shí)，上下兩個(gè)大流速區(qū)域沿顆粒對(duì)稱軸對(duì)稱，且水流方向上顆粒流速變化也對(duì)稱，顆粒后部為流速較小的區(qū)域。對(duì)球心在（13，0.5，0.1）處顆粒，從XY面、YZ面分別觀察Z=0 m、0.1 m、0.2 m和X=12.9 m、13.0 m、13.1 m位置的顆粒流速示意圖，發(fā)現(xiàn)：高速區(qū)域集中在上半球前側(cè)，低速區(qū)域在水流方向上球的前后均有分布，后部范圍較大，且都集中在下半球，尤其球后部靠近床面處水流流速較低，變化劇烈，且流速圖也充分顯示了球形顆粒周圍流速變化的對(duì)稱性。

圖3 流速等值線

2.2 上舉力和拖曳力

泥沙顆粒在流場(chǎng)中受壓力、由水的粘滯性帶來(lái)的表面剪切力和重力作用，其中壓力和剪切力均由水產(chǎn)生，可以綜合為水流力，分解為上舉力和拖曳力，也有不少學(xué)者并不分解，研究總的水流力。研究上舉力和拖曳力的一般公式可得上舉力系數(shù)CL和拖曳力系數(shù)CD，數(shù)碼實(shí)驗(yàn)?zāi)M大雷諾數(shù)情況下的顆粒受力，并針對(duì)不同因素進(jìn)行分析。

對(duì)顆粒水中受力分析時(shí)，應(yīng)同時(shí)考慮力和由力產(chǎn)生的扭矩作用。球體在水中受到壓力，分可分解為三個(gè)方向的分量：Px、Py、Pz，受到剪切力，分解為三個(gè)方向：TAOx、TAOy、TAOz，壓力全部指向球心，剪切力沿表面切線方向，會(huì)對(duì)顆粒產(chǎn)生繞坐標(biāo)軸的扭矩：MTAOx、MTAOy、MTAOz，扭矩會(huì)使顆粒產(chǎn)生轉(zhuǎn)動(dòng)，通過(guò)模擬計(jì)算發(fā)現(xiàn)，扭矩很小，幾乎可以忽略。在研究顆粒受力時(shí)，其合力作用點(diǎn)一直沒(méi)有明確結(jié)論。本文分析，以TAOx為例，TAOx的作用點(diǎn)不在球心時(shí)，會(huì)產(chǎn)生對(duì)Y、Z軸的扭矩，故每一單元面積積分所求得的TAOx直接相加為T(mén)AOx的合力，然后求出合力距離Y、Z軸的距離，然后求解其扭矩的求法卻不能直接累加，因?yàn)榇笮∠嗤?、方向相反的平行剪切力?duì)合力大小不產(chǎn)生影響，但會(huì)產(chǎn)生一個(gè)扭矩。為了達(dá)到這一目的，編寫(xiě)了Fortran程序?qū)δM數(shù)據(jù)進(jìn)行后處理。在達(dá)到計(jì)算精度的前提下，將球劃分了180份，每一部分進(jìn)行積分計(jì)算。

流速的選取對(duì)CL和CD的大小有不能忽略的影響。Einstein（1949）［2］采以距圓球頂面以下0.2D處作為理論床面，然后選取理論床面以上0.35D處的流速進(jìn)行計(jì)算。本文考慮顆粒遮蔽度的影響，采用球前流速不受影響處的某一斷面上、圓球迎水面直徑所對(duì)的沿深度方向的垂線平均流速。通過(guò)建模計(jì)算發(fā)現(xiàn)，模型中流場(chǎng)10 m之后流態(tài)穩(wěn)定，故計(jì)算中選取11 m處對(duì)應(yīng)流速。

2.2.1 顆粒上升高度變化

天然河床上的泥沙顆粒以群體形式存在，泥沙層層疊疊，一般河床上有若干層泥沙，顆粒之間相互支撐相互遮蔽，自然就有相對(duì)上層的顆粒和相對(duì)下層的顆粒，為了研究泥沙距離底床高度這一因素對(duì)受力的影響，進(jìn)行不同上升高度的球體受力模擬，本文從單顆粒上升入手，不考慮遮蔽問(wèn)題，模擬目的在于測(cè)算泥沙由于上升所處水流狀況的不同而造成的受力不同。本次模擬記為第一組，模擬中，首組泥沙顆粒固定在底床上，球體直徑為20 cm，球心坐標(biāo)為（13，0.5，0.1），接下來(lái)各組以0.05 m為間隔上升，上升過(guò)程中一直固定不動(dòng)，嵌套網(wǎng)格也隨之上升，控制入口流速為0.2 m/s，出口邊界為出流邊界。計(jì)算結(jié)果不在此詳述。

由計(jì)算結(jié)果可知，直徑為20 cm的球在入口流速為0.2 m/s流場(chǎng)中受力，壓力Px、Py、PZ，其中 Py較小，幾乎可以忽略；剪切力TAOx、TAOy、TAOz，與壓力相比，剪切力都非常小，相差若干個(gè)量級(jí)，其中TAOx較TAOy、TAOz較大，但其對(duì)拖曳力系數(shù)CD的影響已經(jīng)非常小，只能影響其小數(shù)點(diǎn)后第三位數(shù)字，TAOy、TAOz則影響太小完全可以忽略；剪切力對(duì)顆粒產(chǎn)生的扭矩MTAOx、MTAOy、MTAOz也非常小，幾乎可以忽略。隨著顆粒上升，垂線平均流速越來(lái)越大，在球體球心移動(dòng)至0.30 m時(shí)（項(xiàng)目5），觀察其垂線流速分布，流速全部相等，近似于自由落體沉降運(yùn)動(dòng)的情況，顆粒再上移至0.35 m、0.40 m情況也是如此。

拖曳力系數(shù)CD由Px和TAOx的合力求得，放置于床面時(shí)CD為0.524，顆粒上移CD迅速減小至0.2左右，且隨顆粒上移CD越來(lái)越小，當(dāng)顆粒接近自由落體情況時(shí)（即項(xiàng)目5、6、7）CD穩(wěn)定在0.190～0.191附近。

上舉力系數(shù)CL由Pz和TAOz的合力求得，床面時(shí)CL為3.685，顆粒上移CL逐漸減小，當(dāng)顆粒接近自由落體情況時(shí)（即項(xiàng)目5、6、7）CD穩(wěn)定在2.70附近。

圖4 CD、CL與Z關(guān)系

2.2.2 顆粒直徑變化

泥沙顆粒的形狀會(huì)直接影響其所受上舉力和拖曳力的大小，本文采用球體模擬泥沙顆粒，因此針對(duì)球體直徑這一因素進(jìn)行對(duì)照模擬。第二組模擬與第一組模擬相比，僅改變球體直徑，其他不變，設(shè)定球體直徑為10 cm，球心坐標(biāo)為（13，0.5，0.05），以0.05 m為間隔上升，入口流速為0.2 m/s。計(jì)算結(jié)果略。

經(jīng)計(jì)算發(fā)現(xiàn)，本組模擬中，Py、TAOy、TAOz、MAOX、MTAOY、MTAOZ依然與其他數(shù)據(jù)相差若干數(shù)量級(jí)，可以忽略。隨著顆粒上升，垂線平均流速變化與第一組情況一致，在項(xiàng)目5時(shí)垂線流速幾乎不變，進(jìn)入沉降模式。放置于床面時(shí)，拖曳力系數(shù)CD為0.602，顆粒上移CD迅速減小至0.3左右，且隨顆粒上移CD越來(lái)越小，當(dāng)顆粒接近自由落體情況時(shí)（即項(xiàng)目5、6、7）CD穩(wěn)定在0.303附近。

床面處上舉力系數(shù)CL為5.401，顆粒上移CL逐漸減小，且隨顆粒上移CD越來(lái)越小，當(dāng)顆粒接近自由落體情況時(shí)（即項(xiàng)目5、6、7）CD穩(wěn)定在在3.20～3.21附近。將第二組模擬得到的CD、CL與第一組對(duì)比如圖5。

圖5 第一組與第二組模擬CD、CL與Z關(guān)系對(duì)比

圖6 第一組與第二組模擬曲線對(duì)比

對(duì)比發(fā)現(xiàn)，第二組模擬得到的CD、CL均大于第一組值，即半徑越小CD、CL越大，且CD變化比CL大，但兩組模擬曲線的變化趨勢(shì)一致。

將第二組模擬的Px、Pz、TAOx、TAOz與第一組進(jìn)行對(duì)比。

由四組對(duì)比曲線可知，Px、Pz、TAOx三個(gè)相對(duì)較大的值變化趨勢(shì)相似，第一組模擬值遠(yuǎn)大于第二組，TAOz接近0，兩組值最終都穩(wěn)定在0附近。

2.2.3 流速變化

為了研究流場(chǎng)變化對(duì)顆粒受力產(chǎn)生的影響，改變?nèi)肟谒俣龋紫龋旱谌M模擬，設(shè)定球體直徑為20 cm，設(shè)定進(jìn)口流速為0.4 m/s，然后如第一組模擬一樣以0.0 5m為間隔將顆粒上移；第四組模擬，球體直徑為20 cm，始終將顆粒固定于床面，改變?nèi)肟诹魉?，分別為：0.1 m/s、0.3 m/s、0.5 m/s、0.6 m/s、0.7 m/s、1.0 m/s，0.2 m/s與0.4 m/s分別在第一組和第三組模擬中計(jì)算過(guò)，其他條件均不變。

由計(jì)算結(jié)果可知，直徑為20 cm的球在入口流速為0.4 m/s流場(chǎng)中受力，可以發(fā)現(xiàn)，與前兩組相同，壓力Py，剪切力TAOy、TAOz，扭矩MTAOX、MTAOY、MTAOZ相對(duì)很小，可以忽略。拖曳力系數(shù)CD和上舉力系數(shù)CL沿垂線變化規(guī)律與前兩組模擬一致，本模擬中球體球心移動(dòng)至0.35 m時(shí)（項(xiàng)目6）時(shí)達(dá)到近似于自由落體沉降運(yùn)動(dòng)的情況。拖曳力系數(shù)CD床面時(shí)值為0.517，顆粒上移減小至0.211，隨顆粒上移CD越來(lái)越小，項(xiàng)目4、5、6時(shí)CD穩(wěn)定在0.188附近，項(xiàng)目7略大，為0.191。上舉力系數(shù)CL由PZ和TAOZ的合力求得，此組床面時(shí)CL為1.279，小于第一組，且隨顆粒上移CD越來(lái)越小，最終CD穩(wěn)定在0.68附近，同樣小于第一組。如圖7。

圖7 第一組與第三組模擬CD、CL與Z關(guān)系對(duì)比

從圖7可以發(fā)現(xiàn)，拖曳力系數(shù)CD在流速改變時(shí)變化不大，兩條曲線幾乎吻合；上舉力系數(shù)CL則變化較大，流速增大后CL變小，變小將近1/4，但曲線變化趨勢(shì)相同。

第一組和第三組模擬得到的Px、Pz、TAOx、TAOz做對(duì)比，如圖8。

圖8 第一組與第三組模擬曲線對(duì)比

由圖8可知，流速較大的第三組模擬得到的Px、Pz、TAOx均大于第一組的值，其變化規(guī)律相似，其中Pz值兩條曲線最終均穩(wěn)定于35.5左右，相差較少，TAOz較小，兩條曲線變化趨勢(shì)類似，最終均穩(wěn)定于0附近。

第四組模擬調(diào)整流速，通過(guò)觀察結(jié)果可以發(fā)現(xiàn)，流速改變拖曳力系數(shù)CD幾乎穩(wěn)定不變，在0.52附近浮動(dòng)，上舉力系數(shù)CL隨著流速增大越來(lái)越小，兩值曲線如圖9、圖10。

圖9 CD、CL與V關(guān)系做Px、Pz、TAOx、TAOz曲線

圖10 Px、Pz、TAOx、TAOz關(guān)于流速V的曲線

由圖9、圖10可知，流速與Px、Pz、TAOx、TAOz的關(guān)系比較規(guī)律，其中流速與Px、Pz的曲線為平滑拋物線，與TAOx、TAOz的關(guān)系近似為一條斜線。

2.2.3 力臂、作用點(diǎn)分析

球體顆粒受到壓力和剪切力，壓力均指向球心，不會(huì)產(chǎn)生對(duì)球體的扭矩，剪切力沿表面切線方向，會(huì)產(chǎn)生對(duì)球體的扭矩。剪切力沿坐標(biāo)軸分解：TAOx、TAOy、TAOz、TAOx，會(huì)產(chǎn)生繞Y軸和Z軸的扭矩，本文稱其為MTAOXY、MTAOXZ，TAOx的力臂分別記為XZ、XY；與此相同，TAOy產(chǎn)生扭矩MTAOXY、MTAOXZ，力臂為YZ、YX；TAOz產(chǎn)生MTAOZX、MTAOZY，力臂為ZY、ZX；由此可得，剪切力關(guān)于三個(gè)坐標(biāo)軸的總扭矩：MTAOX、MTAOY、MTAOZ。以第一組模擬項(xiàng)目一為例，見(jiàn)表1。

表1 受力計(jì)算

由表1可知，球體受到的扭矩MTAOX、MTAOY、MTAOZ近似為零，剪切力分別產(chǎn)生的扭矩MTAOXY、MTAOXZ、MTAOYX、MTAOYZ、MTAOZX、MTAOZY也近似為零，相比之下MTAOXY、MTAOZY雖然也很小，但比其他量大若干數(shù)量級(jí)；觀察剪切力TAOx、TAOy、TAOz，相比TAOx較大，TAOz次之，TAOy近似為零；本項(xiàng)目球體半徑為0.10 m，觀察剪切力對(duì)應(yīng)的力臂，MTAOXY、MTAOZY對(duì)應(yīng)力臂XZ、ZX分別為0.013 m、-0.014 m，靠近球心，XY、ZY則近似為零，位于球心上。綜上可發(fā)現(xiàn)，剪切力與壓力相比較小，相比之下僅TAOx、TAOz需考慮，其產(chǎn)生的扭矩MTAOXY、MTAOZY較大，故總扭矩僅MTAOY需考慮，但仍然非常小，而且剪切力作用線幾乎都集中于球心附近，將三個(gè)分力合成，剪切力合力大小為0.002 kN，作用點(diǎn)為（-0.014，0，0.013）。其他組模擬中發(fā)現(xiàn)，扭矩MTAOX、MTAOY、MTAOZ也很小，可以忽略，剪切力作用點(diǎn)也在球心附近。

3 結(jié) 論

本文以球體代替泥沙，對(duì)河床上的單顆粒泥沙在定常流中的受力進(jìn)行了研究，通過(guò)數(shù)值模擬方法，分析其受力、受扭情況，并結(jié)合多個(gè)變量研究了泥沙顆粒拖曳力系數(shù)和上舉力系數(shù)的變化規(guī)律。

研究表明：1）泥沙顆粒在水流中受到水壓力和剪切力的影響，剪切力對(duì)顆粒產(chǎn)生扭矩，但扭矩較??；2）在Rep∈（1×104，2×105）時(shí)，流速改變對(duì)拖曳力系數(shù)CD影響非常小，顆粒半徑則影響較大，球體半徑越小拖曳力系數(shù)CD越大，隨著顆粒離開(kāi)床面上移，CD迅速減小，當(dāng)顆粒處于流速沿垂線變化不大的高度時(shí)，CD穩(wěn)定在一定范圍內(nèi)；在此范圍內(nèi)，上舉力系數(shù)CL大于拖曳力系數(shù)CD，CL變化范圍比CD大，流速越大CL越小，顆粒直徑越大CL也越小，與拖曳力系數(shù)相同的是，隨著顆粒沿垂線上移，CL不斷減小，當(dāng)球前流速沿垂線方向幾乎不變時(shí)，最終穩(wěn)定在一定范圍內(nèi)。

參考文獻(xiàn)：

［1］ Shields A.Anwendung der Aechlichkeitsmechanik und der Turbulenzforschung auf die Geschiebewegung ［J］.Berlin， Mit.Preussische Versuchsanstalt fur Wasserbau，1936.

［2］ Einstein H A， E A EO-Samni.Hydrodynamic Forces on A Rough Wall［J］.Rev Modern Phys， 1949， 21（3）：520-524.

［3］ Brayshaw A C， L E Frostick， L Reid.The hydrodynamics of particle clusters and sediment entrainment in course alluvial channels ［J］.Sedimentology，1983， 30：137-143.

［4］ Schmeeckle M W， Nelson J M， Shreve R L.Forces on stationary particles in near-bed turbulent flows ［J］.Geophys Res， 2007：112.

［5］ Davies T R H， M F A Samad.Fluid dynamic lift on a bed particle ［J］.Hydraul Div Am Soc Civ Eng， 1978，104：1171-1182.

［6］ McLaughlin J B.Internal migration of a small sphere in linear shear flows ［J］.Fluid Mech， 1991：261-274.

［7］ Mei R.An approximate expression for the shear lift force on a particle at a finite Reynolds number ［J］.Multiphase Flow， 1992， 18：145-147.

［8］ Kurse R， Komori S.Drag and lift forces on a rotating sphere in a linear shear flow ［J］.Fluid Mech，1999，384：183-206.

［9］ Ladd A J C.Numerical simulations of particle suspensions via a discretized Boltzmannequation.Part I.Theoretical foundation ［J］.Fluid Mech， 1994，271：285-309.

［10］ Sangani A S， Acrivos A.Slow flow through a periodic array of spheres ［J］.Multiphase Flow， 1982， 8：343-360.

［11］ Pettyjohn E S， Christiansen E B.Effect of Particle Shape on Free-Settling Rates of Isometric Particles ［J］.Chem Eng Prog， 1948， 44（2）：157.

［12］ Hartman M， Coughlin R W.On the Incipient Fluidized State of Solid Particles ［J］.Collect Czech Chem Commun，1993， 58：1213.

［13］ Krueger B， Wirtz S， Scherer V.Measurement of drag coefficients of non-spherical particles with a camera-based method ［J］.Powder Technology， 2015，278：157-170.

［14］ 龐啟秀，徐金環(huán)，辛海霞.塊體形狀對(duì)水流拖曳力的影響［J］.水道港口，2006，27（1）：5-9.

Analysis of Forces on Sediment Particles at Different Depth

Li Shisen1，Zhang Rui2
（Tianjin University， Tianjin 300072， China）

Abstract：It is the foundation and key point for one who studies the swarming movement of sediment particles and participates in the practical engineering to fully understand the movement of single sediment particle.The most natural sediment particles are non-uniform and they move at random， which results in randomness of the flow field near sediment particles.In a very long time， the stress of single sediment particle has been studied in China and abroad.Now， a numerical simulation is made for the stress of single sediment particle under the steady flow by using Flow3D software， which produces the variation rules of towing and uplifting coefficients under different diameters of sediment particle and velocities of flow.

Key words：single sediment particle； towing force coefficient； uplift force coefficient

中圖分類號(hào)：U652.2+1

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1004-9592（2016）02-0006-07

DOI：10.16403/j.cnki.ggjs20160202

收稿日期：2015-10-14

作者簡(jiǎn)介：李世森（1969-），男，副教授，主要從事泥沙運(yùn)動(dòng)規(guī)律及淤積研究。