潰壩水流沖擊彈性板模型試驗(yàn)與數(shù)值模擬

關(guān)鍵詞：流固耦合；潰壩水流；振動響應(yīng)；光滑粒子流體動力學(xué)

水利工程建設(shè)中經(jīng)常出現(xiàn)自由表面流與結(jié)構(gòu)相互作用的現(xiàn)象，典型的流固耦合（FSI）問題包括水閘泄流對閘門的沖擊、波浪荷載作用下壩體產(chǎn)生的振動響應(yīng)以及漁船行駛時水流和漁船之間的相互作用[1-3]。FSI問題均涉及劇烈的自由表面流動和結(jié)構(gòu)邊界大變形的復(fù)雜特征，因此研究此類FSI問題在維護(hù)結(jié)構(gòu)的安全性及穩(wěn)定性方面具有重要的理論和工程實(shí)踐意義。

傳統(tǒng)意義上常采用有限體積法（FVM）、有限元法（FEM）、有限差分法（FDM）來處理FSI問題[4-6]。Attili等[7]基于FVM模擬波浪沖擊柔性擋板，分析了流體和固體之間的相互作用現(xiàn)象，借助物理試驗(yàn)驗(yàn)證了程序的有效性，但FVM在處理大變形問題時存在局限性且常需要采用動網(wǎng)格計(jì)算，解決FSI問題的適應(yīng)性較差；Liao等[8]基于FDM-FEM耦合方法建立二維潰壩水流沖擊彈性板數(shù)值模型，但FDM-FEM法網(wǎng)格依賴性強(qiáng)，在處理FSI問題時收斂性不穩(wěn)定且模擬水流效果不理想。上述研究方法在處理自由液面邊界和大變形問題時易產(chǎn)生網(wǎng)格畸形，導(dǎo)致流固界面跨界插值不準(zhǔn)確，影響計(jì)算精度。然而光滑粒子流體動力學(xué)（SPH）方法作為一種粒子法，具有完全的拉格朗日特性，可有效處理自由液面大變形等問題[9-10]。Zhang等[11]將改進(jìn)SPH方法與光滑有限元法（SFEM）進(jìn)行耦合，更好地模擬彈性結(jié)構(gòu)FSI問題；Tang等[12]在拉格朗日背景下提出SPH-DEM模型，該模型可以準(zhǔn)確高效地捕捉復(fù)雜流動條件下彈性固體的變形及其破壞特征；Xue等[13]將SPH法和再生粒子法（RKPM）結(jié)合建立SPH-RKPM求解器，改進(jìn)了流固邊界處的壓力場從而提高求解器在面對高速沖擊問題時的穩(wěn)定性。然而上述耦合方法較SPH法在計(jì)算的靈活性和流固邊界處理技術(shù)上有顯著差距[14]。

因此，本研究基于SPH法建立潰壩水流沖擊彈性板FSI模型，借助模型試驗(yàn)（EXP）和FEM驗(yàn)證求解器的可靠性；通過改變潰壩的高度，分別從流場性質(zhì)、流固作用力、振動響應(yīng)等角度研究自由表面流沖擊下彈性結(jié)構(gòu)的主要動力學(xué)特征。

1基本理論

SPH方法是一種純拉格朗日型的無網(wǎng)格方法，本研究對流體和固體均采用SPH法進(jìn)行離散，所有的粒子均按照控制方程的規(guī)律進(jìn)行運(yùn)動，根據(jù)Monaghan[15]在選擇適合核函數(shù)的基礎(chǔ)上進(jìn)行核近似和粒子近似。通過對質(zhì)量守恒方程和動量守恒方程的求解，不斷更新粒子的密度和速度，通過動量方程計(jì)算得到速度的導(dǎo)數(shù)再進(jìn)行積分計(jì)算粒子的運(yùn)動軌跡[16-17]。

根據(jù)質(zhì)量守恒方程（連續(xù)性方程），粒子a（流體或固體）密度的導(dǎo)數(shù)可通過SPH插值表示為

2模型試驗(yàn)

本試驗(yàn)?zāi)M潰壩水流對下游結(jié)構(gòu)物的沖擊作用。該試驗(yàn)裝置主要由水箱、照明系統(tǒng)、測量系統(tǒng)、數(shù)據(jù)采集系統(tǒng)、可視化系統(tǒng)等5部分組成，試驗(yàn)裝置如圖1所示。試驗(yàn)的水箱由聚甲基丙烯酸甲酯（PMMA）材料制成，長0.9m，高0.65m，寬0.2m；彈性板高0.12m，寬0.2m，厚0.005m，彈性板距離右側(cè)壁面0.3m；照明系統(tǒng)采用多個工業(yè)LED燈；測量系統(tǒng)采用納米壓痕儀測量彈性板楊氏模量；數(shù)據(jù)采集系統(tǒng)采用幀率為1461fps的高速攝像機(jī)記錄水流及彈性板演變趨勢，最后將試驗(yàn)結(jié)果呈現(xiàn)在屏幕上。

在彈性板上切割出0.01m×0.01m的正方塊體，通過納米壓痕儀測量彈性板上4個測點(diǎn)的楊氏模量和硬度，取其平均值分別為5.77MPa和0.912MPa，試驗(yàn)相關(guān)流程如圖2所示。在彈性板上均勻選取3個不同顏色的測點(diǎn)，使用高速攝像機(jī)追蹤每個步長內(nèi)彈性結(jié)構(gòu)的變形情況，在試驗(yàn)過程中閘門左側(cè)形成一定高度（H）的水柱，通過抽拉擋水板形成潰壩水流，水流在重力的作用下產(chǎn)生垂直加速運(yùn)動，各個測點(diǎn)布置及閘門位移時程曲線如圖3所示。

3模型試驗(yàn)與數(shù)值仿真對比驗(yàn)證

3.1水流流態(tài)對比分析

為了更加直觀地比較試驗(yàn)和數(shù)值模擬結(jié)果，將EXP、SPH、FEM3種方法典型時間步長下的潰壩水流流態(tài)和彈性結(jié)構(gòu)的變形情況進(jìn)行可視化展示，如圖4所示。在潰壩發(fā)生初期水流整體流態(tài)呈流線型，水體粒子均勻地向前運(yùn)動，水流的前端呈斜向下的鋒面狀，隨著潰壩的進(jìn)行水流沖擊彈性板，此時水體粒子的動能轉(zhuǎn)化為彈性板的彈性勢能，水流翻越彈性板后撞擊壁面，飛濺粒子的動能轉(zhuǎn)化為重力勢能，在重力的作用下粒子向下運(yùn)動撞擊下部的水體粒子，在此過程中液體出現(xiàn)飛濺、自由液面破碎、部分水流回旋等現(xiàn)象。

通過對比不同時刻的水流流態(tài)可以發(fā)現(xiàn)，基于SPH法建立的數(shù)值模型與試驗(yàn)更加接近。t=0.26s時水流第1次沖擊下彈性板達(dá)到最大負(fù)向位移；t=0.56s時水舌接觸壁面，彈性板在振動水流作用下產(chǎn)生振蕩現(xiàn)象，水流撞擊壁面導(dǎo)致水流回流第2次撞擊彈性板，在水壓力的作用下彈性板的變形再次增大，此時水箱右側(cè)水流較多且動能較大，在水流動能的驅(qū)動下彈性結(jié)構(gòu)逐漸回彈；在t=1.0s時，水流回旋倒流向左側(cè)運(yùn)動，EXP和SPH求解器下的彈性板均停留在負(fù)方向，而FEM結(jié)果與試驗(yàn)相比有一定的差異。

在1s的整個過程中，SPH方法對每個時刻自由液面模擬結(jié)果和結(jié)構(gòu)變形程度均與試驗(yàn)結(jié)果保持一致，且SPH方法展現(xiàn)出良好的線性特點(diǎn)，而FEM較SPH方法和試驗(yàn)結(jié)果相比，每個時間步內(nèi)結(jié)構(gòu)變形更大、流速更快，基于FEM模擬大變形問題存在一定的誤差，模擬效果略差于SPH方法。

3.2自由液面及彈性板位移對比分析

為進(jìn)一步探究SPH方法的可靠性和優(yōu)越性，結(jié)合FEM和EXP結(jié)果，選取3個典型時刻（t=0.2s、t=0.4s、t=0.6s），三者自由液面輪廓曲線如圖5所示：\"

為了定量比較EXP、SPH、FEM方法模擬的差異性，對自由液面曲線提取特征點(diǎn)位置坐標(biāo)，計(jì)算三者的相似度。標(biāo)記試驗(yàn)數(shù)據(jù)點(diǎn)為N1，數(shù)值模擬數(shù)據(jù)點(diǎn)為N2，自由液面相似度（）計(jì)算公式如下：

式中：y1i為試驗(yàn)數(shù)據(jù)縱坐標(biāo)值；y2i為數(shù)值模擬結(jié)果縱坐標(biāo)值。

通過計(jì)算可得SPH、FEM和EXP的平均相似度分別為90.93％和82.07％，計(jì)算結(jié)果見表1。結(jié)果表明基于SPH方法的模擬結(jié)果與EXP吻合度較高，SPH方法能更好地模擬潰壩水流沖擊彈性板這一過程。

在彈性板上選取3個測點(diǎn)繪制1s內(nèi)位移時程變化曲線，如圖6所示。2種數(shù)值模擬方法下彈性板的變形趨勢與試驗(yàn)結(jié)果均保持一致，但SPH法與試驗(yàn)結(jié)果更加接近。由圖6可知，F(xiàn)EM數(shù)值模擬的精度略低，原因是當(dāng)水流撞擊水箱右壁回流時，水流中產(chǎn)生部分空腔，空腔內(nèi)湍流的劇烈運(yùn)動對彈性結(jié)構(gòu)產(chǎn)生振動沖擊荷載，0.35～0.55s內(nèi)彈性板開始進(jìn)行小幅度的周期性振動，此時彈性結(jié)構(gòu)展現(xiàn)出高階振動模態(tài)，在高階振型中空腔內(nèi)的湍流運(yùn)動起主導(dǎo)作用。有關(guān)振動部分的研究在4.3節(jié)具體闡述。

綜上所述，SPH方法相對于FEM在水流流態(tài)、自由液面曲線相似度、彈性板變化趨勢等方面與EXP結(jié)果更為接近，從而說明SPH求解器的有效性、準(zhǔn)確性和優(yōu)越性。

4流固耦合本質(zhì)特征分析

4.1渦量的產(chǎn)生及流速分布

渦量是描述流體旋轉(zhuǎn)性質(zhì)的物理量，速度分布表示流體中各個粒子速度的大小和方向，渦量與速度的分布密切相關(guān)，它們共同決定了流場的性質(zhì)。本研究的渦量即為！，！由式（5）計(jì)算可得。

在物理試驗(yàn)中數(shù)據(jù)的空間分辨率相對粗糙，無法用有限差分法確定式（5）中的速度的導(dǎo)數(shù)；而在數(shù)值仿真中當(dāng)拉格朗日粒子之間的距離小于0.1cm時，可以用核函數(shù)來計(jì)算速度的導(dǎo)數(shù)且不會產(chǎn)生較大的誤差。

在潰壩數(shù)值模擬過程中，在不同潰壩高度下，產(chǎn)生的渦沿不同路徑在不同時刻出現(xiàn)的大小均不相同，圖7給出不同潰壩高度下典型時刻的渦量分布圖。水流初次撞擊彈性板時，彈性板上端水流及自由液面出現(xiàn)較為明顯的渦流現(xiàn)象，且產(chǎn)生順時針的渦量；隨著水流回流再次撞擊彈性板，渦流現(xiàn)象更加顯著，且產(chǎn)生逆時針的渦量。渦量的產(chǎn)生與轉(zhuǎn)變的原因如下：①渦量是由空腔內(nèi)自由表面的強(qiáng)曲率產(chǎn)生的[19]，在空腔內(nèi)自由表面的強(qiáng)曲率會導(dǎo)致流體速度場的變化，此時渦量產(chǎn)生在自由表面處；②空腔周圍粒子的速度一般不為0且速度較大，當(dāng)空腔坍塌時根據(jù)斯托克斯循環(huán)定理，一定會產(chǎn)生渦量通量；③由于水流之間速度的跳躍導(dǎo)致渦流由順時針方向轉(zhuǎn)變?yōu)槟鏁r針方向。

同樣渦量和速度矢量場之間存在一定相關(guān)性，限于篇幅本研究僅選取0.3m潰壩高度的速度場進(jìn)行展示，模擬結(jié)果如圖8所示。將典型時刻的速度矢量場與渦量進(jìn)行對比分析，發(fā)現(xiàn)渦量的大小取決于速度場的復(fù)雜性和變化程度，而不僅僅是速度的大小，在某些時刻高速運(yùn)動下的流體可能會產(chǎn)生較大的渦量，但在某些低速流體下也會產(chǎn)生較大的渦量。總體來說SPH數(shù)值模擬結(jié)果展示了流體的渦量場和速度場，很好地模擬和預(yù)測了流體中的渦流和漩渦現(xiàn)象。

4.2流固作用力對比分析

分析彈性板的流固作用力有助于更好地研究水流與彈性板之間的耦合作用，如圖9（H0為潰壩的水頭高度，h和b分別為彈性板的高度和厚度）所示，在彈性板上均勻選取6個測點(diǎn)，測點(diǎn)A的高度為0.0855m，測點(diǎn)F的高度為0.004m，在2個測點(diǎn)間均勻地選取4個測點(diǎn)分別為B、C、D、E，每個測點(diǎn)的間距為0.0163m。限于篇幅本研究僅展示A、D、F3個測點(diǎn)在不同潰壩高度下流固作用力的變化情況。

通過式（6）可以計(jì)算結(jié)構(gòu)和流體之間的相互作用力，據(jù)計(jì)算流體粒子對固體粒子的作用力等于固體粒子對流體粒子的力，驗(yàn)證了動態(tài)邊界作用-反作用原理以及動態(tài)邊界的本質(zhì)，進(jìn)一步驗(yàn)證SPH求解器的有效性。

圖10展示了A、D、F特征點(diǎn)根據(jù)式（6）計(jì)算不同潰壩高度下流固作用力的變化趨勢。在不同潰壩高度下主要分為2個典型時刻，第1個典型時刻為水流對彈性板的首次沖擊，第2個典型時刻為水流回流第2次撞擊彈性板。通過對比A、D、F3個測點(diǎn)的流固作用力可知，0.4m的潰壩高度在勢能和動能的作用下率先撞擊彈性板且產(chǎn)生較大的流固作用力；不同潰壩高度下所產(chǎn)生最大的流固作用力均在水流第2次沖擊彈性板之后。其中，0.4m高的潰壩水流所產(chǎn)生的最大力為81.1N，0.3、0.2m工況下產(chǎn)生的最大力分別為71.89、64.7N?；赟PH方法模擬潰壩水流沖擊彈性板過程可知：①在整個過程中，流固作用力呈現(xiàn)出先增大后減小，再增大再減小的變化趨勢；②彈性板所受力的峰值出現(xiàn)在水流二次沖擊彈性板的初期；③彈性板的變形程度對流固力的大小有顯著影響，彈性板變形越小所受到的力越大；④在0.35～0.55s水流翻越彈性板時，流固作用力呈現(xiàn)周期性的變化規(guī)律。

4.3彈性板頻率響應(yīng)分析

彈性板的振動頻率反映了彈性結(jié)構(gòu)的固有特性，分析其頻率有助于了解潰壩水流對彈性板沖擊過程中板的動態(tài)響應(yīng)特性。在物理試驗(yàn)部分，測試彈性板在空氣中的自由振動，分析結(jié)構(gòu)的固有頻率。通過理論式（7）和式（8）計(jì)算可知彈性結(jié)構(gòu)一階固有頻率為3.285Hz；在物理試驗(yàn)中采用0.3m高的潰壩水流沖擊彈性板，流激振動下固有頻率為2.741Hz，由于水流附加質(zhì)量的影響導(dǎo)致激發(fā)態(tài)頻率小于結(jié)構(gòu)的固有頻率，由圖11可知流激振動下的前三階頻率分別為2.741、6.051、9.360Hz。

本研究在模型試驗(yàn)中只能分析彈性結(jié)構(gòu)自由端自然振動頻率和水流沖擊下的流激振動頻率，無法準(zhǔn)確分析彈性結(jié)構(gòu)各部分的振動特性，因此借助SPH數(shù)值方法更加具體地研究彈性結(jié)構(gòu)各部分的振動情況。限于篇幅并和模型試驗(yàn)形成對比，選取潰壩水流高度為0.3m，在彈性上等間距布置多個特征點(diǎn)，特征點(diǎn)如圖9所示。

在水流沖擊彈性板過程中，0.35～0.55s時流固作用力呈現(xiàn)周期性的變化規(guī)律，且分析彈性板的位移可知，在同時段彈性結(jié)構(gòu)的振動屬于平穩(wěn)期，因此僅針對此時段采用Zhang等[20]提出的WTEMD方法濾除強(qiáng)背景噪聲，提取彈性結(jié)構(gòu)各個測點(diǎn)的振動特征信息，識別結(jié)果如圖12和圖13所示。

根據(jù)圖12和圖13可知，采用WTEMD濾波降噪后，測點(diǎn)A—F的一階振動頻率分別為2.44、2.56、2.68、2.76、2.80和2.88Hz，二階振動頻率分別為5.68、5.76、5.84、6.08、6.16和6.24Hz，三階振動頻率分別為8.96、9.12、9.32、9.44、9.52和9.57Hz。通過分析彈性結(jié)構(gòu)前三階頻率可知，SPH各特征點(diǎn)數(shù)值結(jié)果與試驗(yàn)結(jié)果自由端前三階頻率基本保持一致；彈性板自由端沿固定端前三階頻率逐漸增大，其中，特征點(diǎn)A—F的第一階頻率由2.44Hz增大至2.88Hz，第二階頻率由5.68Hz增大至6.24Hz，第三階頻率由8.96Hz增大至9.57Hz。相關(guān)研究結(jié)果表明彈性板沿自由端振動頻率逐漸降低，彈性板在受到潰壩水流的沖擊后，其自由端變形較大振動響應(yīng)逐漸減弱，該計(jì)算結(jié)果可為大變形結(jié)構(gòu)振動問題提供借鑒與參考。

5結(jié)論

本研究基于光滑粒子流體動力學(xué)方法建立潰壩水流沖擊彈性板流固耦合數(shù)值模型，同步模型試驗(yàn)及有限元法驗(yàn)證光滑粒子流體動力學(xué)求解器的有效性，同時分析流場的渦量、流固作用力、彈性結(jié)構(gòu)振動響應(yīng)進(jìn)一步剖析彈性結(jié)構(gòu)流固耦合的本質(zhì)特征，所得結(jié)論如下：

（1）潰壩水流沖擊彈性板過程中，由于自由液面的強(qiáng)曲率和水流之間速度的跳躍導(dǎo)致渦量的產(chǎn)生；當(dāng)不同流速的水流撞擊后渦量的方向發(fā)生轉(zhuǎn)變。

（2）彈性板在水流沖擊下固定端到自由端的流固作用力逐漸減小，彈性結(jié)構(gòu)變形越小受到的流固作用力相對較大且流固作用力峰值出現(xiàn)在水流二次撞擊彈性板初期。

（3）彈性板在受到潰壩水流的沖擊后，彈性板一階振動頻率分布均在2.5Hz左右，二階振動頻率和三階振動頻率分別在6.0Hz和9.0Hz左右。彈性板固定端沿自由端的振動頻率逐漸降低，自由端變形較大其振動響應(yīng)逐漸減弱。