孫強(qiáng)張智霖伍悅濱徐瑩

(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 引言

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

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

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

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

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

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

一維黏彈性管道瞬變流控制方程由連續(xù)性方程和動(dòng)量方程組成[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為時(shí)間,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個(gè)彈性單元和N個(gè)黏彈性單元串聯(lián)構(gòu)成,其中彈性單元由彈性元件彈簧表示,遵循虎克定律,黏彈性單元由N個(gè)彈簧和黏壺并聯(lián)組成的K-V元件構(gòu)成,用于描述黏彈性管道的蠕變和松弛特性。

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

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

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

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

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

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

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

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

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

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

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

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

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

其波速計(jì)算公式采用如下公式

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

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

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

2 試驗(yàn)與數(shù)據(jù)

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

圖3 實(shí)驗(yàn)裝置示意圖[24]Fig.3 Schematic diagram of the experimental setup

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

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

3 結(jié)果與分析

3.1 波速的確定

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

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

3.2 本構(gòu)曲線

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

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

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

3.3 壓力波動(dòng)曲線

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

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

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

4 結(jié)論

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

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

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

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