道路綜合管廊輸水管道參數(shù)辨識研究

孫強張智霖伍悅濱徐瑩

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

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

0 引言

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

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

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

1 辨識模型的構(gòu)建

城市道路綜合管廊黏彈性管道參數(shù)辨識模型的構(gòu)建基于黏彈性管道瞬變流模型結(jié)合K-V力學(xué)模型和一維擬穩(wěn)態(tài)摩阻模型,以實驗數(shù)據(jù)與數(shù)值計算結(jié)果的壓力差值平方和最小為目標(biāo),給定本構(gòu)參數(shù)的約束條件,采用優(yōu)化算法進行求解的過程。

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

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

其中,穩(wěn)態(tài)管壁剪切應(yīng)力用Darcy-Weisbach公式表示為

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

1.2 黏彈性管道本構(gòu)模型

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

圖1 K-V力學(xué)模型Fig.1 Kevin-Voigt model

采用K-V模型,應(yīng)用拉普拉斯變換,蠕變函數(shù)可表示為

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

1.3 本構(gòu)參數(shù)辨識方法

基于已知測點的瞬變流壓力波動數(shù)據(jù),采用SQP算法,通過對實驗與模擬的壓頭之間差的平方和建立目標(biāo)函數(shù),實現(xiàn)誤差最小化,辨識黏彈性管道本構(gòu)參數(shù)。其數(shù)學(xué)關(guān)系式如下

式中:Hi為模擬壓力水頭;H*i為實驗壓力水頭;Nx為節(jié)點總數(shù)。

約束條件:每個K-V元件的蠕變?nèi)崃縅k和延遲時間τk在合理的數(shù)量級范圍內(nèi)。

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

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

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

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

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

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

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

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

圖2 本構(gòu)參數(shù)辨識流程路線圖Fig.2 Flow diagram of calibration of constitutive parameters

(5)用SQP算法辨識本構(gòu)參數(shù),即管道的蠕變?nèi)崃縅k和延遲時間τk。

2 試驗與數(shù)據(jù)

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

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

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

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

3 結(jié)果與分析

3.1 波速的確定

從理論計算波速的模擬結(jié)果可以看出,模擬與實驗之間的誤差較大。故波速的數(shù)值仍需進行校核。本研究的計算波速參考文獻[23],基于實驗所測得的本構(gòu)曲線,將計算得到的波速進行試算,考慮到一維擬穩(wěn)態(tài)摩阻模型的局限性,僅考慮壓力波第一個峰值的擬合情況,最終確定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。同時,在瞬態(tài)流動數(shù)值模擬中,將管道劃分為64個等長單元(Nx=64)。圖4為不同水溫情況下計算波速的校核結(jié)果。由圖4可以看出,經(jīng)過校核得到的壓力波波速在描述第一個壓力波動峰值時較準(zhǔn)確,可以用來作為初始波速進行后續(xù)的綜合管廊黏彈性管道本構(gòu)參數(shù)辨識。

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

3.2 本構(gòu)曲線

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

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

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

3.3 壓力波動曲線

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

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

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

4 結(jié)論

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

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

(2)基于黏彈性管道瞬變流模型,結(jié)合一維擬穩(wěn)態(tài)摩阻模型,所提出的本構(gòu)參數(shù)辨識模型,可較準(zhǔn)確辨識綜合管廊管道本構(gòu)參數(shù)。

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