道路綜合管廊輸水管道參數辨識研究

孫強張智霖伍悅濱徐瑩

(1.東北林業大學 土木工程學院,哈爾濱 150040; 2.東北林業大學 人工環境控制與能源應用研究所,哈爾濱 150040; 3.哈爾濱工業大學 建筑學院,哈爾濱 150090; 4.哈爾濱工業大學 寒地建筑科學與工程研究中心,哈爾濱 150090;

5.哈爾濱商業大學 能源與建筑工程學院,哈爾濱 150028)

0 引言

在城市道路綜合管廊長距離輸水管道系統中,高分子聚合物管材應用愈來愈廣泛,并逐漸取代傳統彈性管材(如鍍鋅鋼管和鑄鐵)[1-2]。這類高聚物管材,又稱其為黏彈性管道,相較于傳統管材而言,兼具黏性和彈性力學行為[3-5]。當綜合管廊管路系統中因水泵突然啟停而發生停泵水錘時,黏彈性管道動力學行為與彈性管道截然不同,管壁應變出現延遲,這導致黏彈性瞬變流過程中波速降低、壓力峰值減小、壓力波相位延遲。因此,準確描述綜合管廊長距離輸水管道水力瞬變動力學行為是十分必要的。

要準確描述綜合管廊黏彈性管道瞬變流壓力波動,需要確知本構參數,因此,黏彈性管道本構參數辨識在瞬變流數值計算中非常重要。對于黏彈性管道本構參數辨識的研究,許多學者通常是借助優化算法(如最小二乘估計[6]、遺傳算法[7-8]等)對本構參數進行辨識。Covas等[6,9]利用最小二乘法辨識本構參數,并分析了K-V元件個數對于辨識結果準確性的影響。研究表明,3個及以上的K-V元件個數對于辨識結果準確性的影響較小。Pezzinga等[7-8]利用微遺傳算法辨識本構參數,研究表明KV元件個數過多時,無法確知本構參數與壓力波動周期間的關系。Gally等[10]基于力學拉伸實驗確定管材本構曲線,數值計算得到較好結果。此外,本構參數的辨識對不同管長[11-12]、摩阻模型[13-14]以及氣液兩相流[15-16]等情況下的黏彈性管道瞬變流的數值計算研究也產生重要影響。

對于辨識問題采用的優化算法來說,最小二乘估計由于其簡便性可能存在一定誤差,而遺傳算法辨識準確,但卻過于復雜,對工程而言耗時較長。因此,本研究基于一種簡便且精確的優化算法——序列二次規劃(Sequential Quadratic Programming ,SQP),利用一維擬穩態摩阻模型結合Kelvin-Voigt(K-V)模型,辨識綜合管廊黏彈性管道本構參數(蠕變柔量和延遲時間),并就不同溫度下的本構曲線和壓力曲線辨識結果進行討論,揭示不同溫度下辨識結果的影響機制。

1 辨識模型的構建

城市道路綜合管廊黏彈性管道參數辨識模型的構建基于黏彈性管道瞬變流模型結合K-V力學模型和一維擬穩態摩阻模型,以實驗數據與數值計算結果的壓力差值平方和最小為目標,給定本構參數的約束條件,采用優化算法進行求解的過程。

1.1 黏彈性管道瞬變流控制方程

一維黏彈性管道瞬變流控制方程由連續性方程和動量方程組成[17-19]

其中,穩態管壁剪切應力用Darcy-Weisbach公式表示為

式中:H為管道內壓力水頭,m;Q為流量,m3/s;a為波速,m/s;εr為管道的延遲應變;V為平均流速,m/s;τw為剪切應力,Pa;x為軸向距離,m;t為時間,s;g為重力加速度,m/s2;ρ為流體密度,kg/m3;f為摩擦系數;D為管道直徑,m;A為管道橫截面積,m2。

1.2 黏彈性管道本構模型

本研究所采用的黏彈性本構模型為K-V力學模型,如圖1所示,該模型[20-21]由1個彈性單元和N個黏彈性單元串聯構成,其中彈性單元由彈性元件彈簧表示,遵循虎克定律,黏彈性單元由N個彈簧和黏壺并聯組成的K-V元件構成,用于描述黏彈性管道的蠕變和松弛特性。

圖1 K-V力學模型Fig.1 Kevin-Voigt model

采用K-V模型,應用拉普拉斯變換,蠕變函數可表示為

式中:J0為瞬時柔量,J0=1/E0,E0為瞬時彈性模量;Jk為第k個K-V元件的蠕變柔量,Jk=1/Ek,Ek為第k個彈簧的彈性模量;τk=Fk/Ek為第k個黏壺的松弛時間;Fk為第k個K-V元件中黏壺的黏性系數;e為壁厚。

1.3 本構參數辨識方法

基于已知測點的瞬變流壓力波動數據,采用SQP算法,通過對實驗與模擬的壓頭之間差的平方和建立目標函數,實現誤差最小化,辨識黏彈性管道本構參數。其數學關系式如下

式中:Hi為模擬壓力水頭;H*i為實驗壓力水頭;Nx為節點總數。

約束條件:每個K-V元件的蠕變柔量Jk和延遲時間τk在合理的數量級范圍內。

SQP算法是非線性規劃算法中的一類特殊數學規劃問題,利用泰勒級數將目標函數在迭代點處簡化為二次函數,將約束條件簡化為線性函數的一種求解最優解問題的算法[22-23]。通過給定初始值、收斂精度,在約束范圍內進行遍歷搜索,直到找到滿足精度的最優解,結束算法。SQP算法收斂性好、計算效率高、邊界搜索能力強,但在應用此算法進行辨識求解時,需注意初始點的選取應結合工程實踐經驗給出。

該辨識方法的具體過程,如圖2所示。根據圖2做以下分析。

(1)確定待辨識的黏彈性管材基本信息,包括管長(L)、流速(V)、管徑(D)、壁厚(e)、密度(ρ)、阻力系數(f)和初始壓力(H0)。

(2)采用特征線法求解瞬變流模型,確定壓力H。

(3)確定黏彈性管道本構模型中K-V元件個數。本文根據Covas等[6,9]給出的建議值,選取3個元件的K-V模型。

(4)確定瞬變流壓力波的波速。

其波速計算公式采用如下公式

式中:νp為泊松比;K為流體體積模量,Pa;E0為彈性管道楊氏模量,Pa。

圖2 本構參數辨識流程路線圖Fig.2 Flow diagram of calibration of constitutive parameters

(5)用SQP算法辨識本構參數,即管道的蠕變柔量Jk和延遲時間τk。

2 試驗與數據

本文采用Gally等[10]所設計的黏彈性管路水力瞬變流室內實驗,該實驗裝置是一個由定壓罐、低密度聚乙烯管和末端快關閥門組成的輸水系統,如圖3所示。實驗管道總長43.1 m,管內徑41.6 mm,壁厚4.2 mm。管路的兩端用固定支架固定。定壓罐體積為15 m3,上下游敞口水箱的體積分別為9 m3,其中上游水箱配有加熱和水溫控制裝置。另外,該實驗的蠕變曲線由Gally等[10]通過流變振動儀進行動態測試所得,在限定的頻率范圍(0.1～11 Hz)內直接測量輸水管路的應力和應變,通過測得儲能模量E1和內摩擦角δ,最終數值計算求得蠕變曲線[10,24]。

圖3 實驗裝置示意圖[24]Fig.3 Schematic diagram of the experimental setup

實驗通過加熱和水溫控制裝置對水箱中流體加熱,從而進行不同溫度下的末端關閥實驗,關閥時間12 ms,3個高精度壓力傳感器安裝在管路的上中下游,記錄不同位置的瞬變流壓力,具體實驗參數見表1。本文主要分析下游末端閥門處的壓力變化情況,為道路綜合管廊輸水管道參數辨識提供指導。

表1 不同情況下的實驗數據[24]Tab.1 Experimental data for different cases

3 結果與分析

3.1 波速的確定

從理論計算波速的模擬結果可以看出,模擬與實驗之間的誤差較大。故波速的數值仍需進行校核。本研究的計算波速參考文獻[23],基于實驗所測得的本構曲線,將計算得到的波速進行試算,考慮到一維擬穩態摩阻模型的局限性,僅考慮壓力波第一個峰值的擬合情況,最終確定13.8 ℃時壓力波波速為379.42 m/s,25 ℃時壓力波波速為351.3 m/s,31 ℃時壓力波波速為341.4 m/s,35 ℃時壓力波波速為320.5 m/s,38.5 ℃時 壓力波波速為315.4 m/s。同時,在瞬態流動數值模擬中,將管道劃分為64個等長單元(Nx=64)。圖4為不同水溫情況下計算波速的校核結果。由圖4可以看出,經過校核得到的壓力波波速在描述第一個壓力波動峰值時較準確,可以用來作為初始波速進行后續的綜合管廊黏彈性管道本構參數辨識。

圖4 不同溫度下波速的校核結果Fig.4 Calibration of wave speed at different water temperatures

3.2 本構曲線

圖5為不同水溫的本構曲線辨識結果。由圖5可以看出,在不同水溫下,綜合管廊黏彈性管道的本構特性是不相同的,其蠕變函數受水溫的影響較大,且隨著水溫的升高,管道的蠕變參數逐漸增大。蠕變曲線辨識結果與實驗數值計算的結果存在一定的差異,這與黏彈性管道本身的應力應變歷史積累有一定關系。同時管道瞬變流實驗壓力波與模擬壓力波存在時間上的不同步性,以及壓力波波速在辨識過程中的差異,都導致了兩者之間的差異。

圖5 不同溫度下蠕變參數辨識結果Fig.5 Calibration of creep parameters at different water temperatures

另外,通過對實驗的蠕變曲線觀察可知,蠕變函數在0～0.07 s內數值陡增,隨后管道蠕變曲線隨時間緩慢增加,不同水溫的管道蠕變曲線最大值出現時間不同并且數值也不同。水溫為13.8、25、31 ℃的蠕變曲線,在4 s左右達到最大值,水溫為35 ℃和38.5 ℃的蠕變曲線,在5 s左右達到最大值。

3.3 壓力波動曲線

不同溫度下壓力波動計算結果,如圖6所示。由圖6可以看出,當水溫為13.8 ℃時,辨識出的壓力波動峰值和谷值均低于實驗結果。當水溫為25 ℃時,數值計算的第一個壓力峰值略高于實驗結果,其后的壓力峰值略高于實驗結果。當水溫為31、35、38.5 ℃時,壓力波動的峰谷值及相位的計算結果與實驗結果趨于吻合。

圖6 不同溫度下壓力曲線辨識結果Fig.6 Calibration of pressure curve at different water temperatures

對比不同溫度下的結果可知,隨著水溫升高,數值計算的壓力波動結果更趨近于實驗數據,這是由于本構曲線辨識過程中,摩阻項的影響越來越小、相應的管道黏彈性項的影響越來越大所導致的[21,25]。

4 結論

本文提出了一種基于瞬變流分析的綜合管廊黏彈性管道本構參數辨識方法,并討論研究了不同溫度下辨識結果的差異。

(1)SQP算法能夠快速準確辨識本構曲線,相比于實驗的本構曲線,得到了較準確的本構參數,同時壓力曲線的計算結果與實驗壓力結果接近。

(2)基于黏彈性管道瞬變流模型,結合一維擬穩態摩阻模型,所提出的本構參數辨識模型,可較準確辨識綜合管廊管道本構參數。

(3)通過對不同水溫下的本構曲線進行分析可知,蠕變函數受溫度的影響較大,且隨著水溫的升高,綜合管廊管道的蠕變參數逐漸增大,并且城市道路綜合管廊黏彈性管道參數辨識模型模擬壓力波動結果在峰值和相位上更接近實驗結果。