黃文強(qiáng)
(蘭州市城市建設(shè)設(shè)計院,甘肅 蘭州 730050)
某縣城工業(yè)園區(qū)供水工程的供水管道項(xiàng)目中,水源地至凈水廠的地形高差大、輸送距離遠(yuǎn),屬于長距離輸水管道工程。該項(xiàng)目設(shè)計供水規(guī)模為29000 m3/d,水源地位于縣城外的水庫,該水庫每天可提供30000 m3的水量,滿足產(chǎn)業(yè)園區(qū)每天29000 m3的水量要求。水源地高程為2170.300 m,凈水廠高程為2017.600 m,兩地高差152.7 m,輸送距離為13.64 km。由于水源地與凈水廠間地形高差較大,因此優(yōu)先選擇有壓重力輸水方式。
給水管道管材的選擇,應(yīng)根據(jù)管徑、內(nèi)壓、外部荷載和管道敷設(shè)場地的地形、地質(zhì),按照運(yùn)行安全、耐久、減少漏損、施工和維修方便、經(jīng)濟(jì)、安全等因素綜合分析來確定。根據(jù)該工程的特點(diǎn)和選擇管材的原則,該工程由于施工期短,球墨鑄鐵管等管材較重運(yùn)輸困難,而且在當(dāng)?shù)毓I(yè)園區(qū)中就有生產(chǎn)高密度聚乙烯管的企業(yè),因此該工程采用PE100級高密度聚乙烯管(HDPE),管道采用熱熔連接。
輸水管線最高日量29000 m3/d,最高日平均時供水量為1209 m3/h,輸水管道應(yīng)根據(jù)技術(shù)經(jīng)濟(jì)比較確定管徑,應(yīng)考慮到充分利用地形高差,使輸送設(shè)計流量所采用的管徑最小,以減少工程投資,但是當(dāng)流速越大,運(yùn)行時引起的水錘升壓就越高,一般按經(jīng)濟(jì)流速0.9~1.5 m/計算管徑,且不宜大于3 m/s。經(jīng)計算敷設(shè)一根輸水管,管徑為DN560時,流速為1.711 m/s。
有壓重力輸水管道要求水源與凈水廠的水頭差大于輸水管道的水頭損失與凈水廠所需水頭之和,因此輸水管的末端存在靜水壓。輸水管道的公稱壓力應(yīng)根據(jù)最大使用壓力確定,在最大使用壓力上附加0.2~0.4 MPa的安全余量,提高管道公稱壓力等級會相應(yīng)提高工程投資,經(jīng)過技術(shù)經(jīng)濟(jì)比較,該工程采用公稱壓力為Pn=1.0 MPa的HDPE管道。
3.3.1 柯爾勃洛克-懷特
管道水力計算包括沿程水頭損失計算和局部水頭損失計算,沿程水頭損失 hf=λ(L/d)×(v2/2g)是通用公式,而λ系數(shù)的取值與管道斷面形狀、管材、水流狀態(tài)、水溫等因素有關(guān)。供水管道水力計算公式常用的有滿寧公式、海曾-威廉公式、達(dá)西公式、柯爾勃洛克-懷特公式、舍列維夫公式等,不同的公式有不同的適用范圍,計算結(jié)果也略有不同。資料顯示,柯爾勃洛-懷特公式最適合于PE、PP管的水力計算,現(xiàn)行《埋地塑料給水管道工程技術(shù)規(guī)程》(CJJ 101—2016)所采用的就是柯爾勃洛-懷特公式,管道沿程水頭損失按下式計算:
其中水力摩阻系數(shù)λ系數(shù)按下式計算:
式中:di為管道內(nèi)徑,m;g 為重力加速度,9.81 m/s2;L為管段長度,m;Re為雷諾數(shù);t為水溫,℃;v為平均流速,m/s;λ為水力摩阻系數(shù);γ為水的運(yùn)動粘滯度,m2/s;Δ為管道當(dāng)量粗糙度,一般取0.010~0.013×0.001;
由于給水管道采用熱熔連接,接口處在管內(nèi)壁形成了凸起的內(nèi)環(huán),由于缺少計算資料,因此管道局部水頭損失按沿程水頭損失的10%計算。
3.3.2 數(shù)值解法的選擇
求解水力摩阻系數(shù)計算公式中的λ,可將公式變?yōu)榍蠼夥匠蹋╢λ)=0的方程,λ為變量,其余均為已知量,公式變?yōu)槌椒匠糖蠼?,對于此類一元非線性實(shí)函數(shù),一般沒有解析解,只能采用求近似解的數(shù)值解法。非線性方程數(shù)值解法有二分法、簡單迭代法、牛頓迭代法及割線法等。
(1)二分法是將含根區(qū)間[a,b]逐次分半,檢查小區(qū)間端點(diǎn)函數(shù)值符號變化,以確定更小的含根區(qū)間。具體作法是先取x0=(a+b)/2,如果 (fx0)=0,則x*=x0;否則 (fa)f(x0)<0或 (fx0)f(b)<0,且僅有一個不等式成立,若 (fa)(fx0)<0,記a1=a,b1=x0,則 x*∈(a1,b1)。再取 x1=( a1+b1)/2,若 f(x1)=0,則x*=x1;否則,可得 x*∈(a2,b2)。重復(fù)上述過程,經(jīng) k 次等分后有根區(qū)間為(ak,bk),且(bk-ak)/2≤ε。取xk=(ak+bk)/2作為所求根x*的近似值。二分法的優(yōu)點(diǎn)是程序簡單,對函數(shù)(fx)光滑性要求不高,收斂速度與公比為1/2的等比級數(shù)相同,收斂速度較慢,適于求初始近似值,然后再用其它方法進(jìn)一步精確化。
(2)簡單迭代法是用極限過程來逐步逼近所給問題精確解的計算方法,是在計算過程中常用的一種方法。用迭代法求非線性方程(fx)=0的近似根,首先需要將方程(fx)=0改寫成等價方程x=g(x),g(x)為連續(xù)函數(shù),當(dāng)|g(x)′|<1,g(x)為迭代函數(shù),對于方程 x=g(x),選取初選值 x0,構(gòu)造序列 xk+1=g(xk)(k=0,1,2...),當(dāng) k→∞ 時,如果迭代過程收斂,xk→x*,則x*為方程的解。其收斂速度為線性收斂。
迭代過程的幾何意義就是在xy平面上確定曲線 y=g(x)與直線 y=x 的交點(diǎn),首先在曲線 y=g(x)上可確定一點(diǎn)P0,過P0引平行x軸的直線,設(shè)此直線與直線y=x交于點(diǎn)Q1,然后過Q1再作平行于y軸的直線,它與曲線y=g(x)的交點(diǎn)記作P1,過P0引平行x軸的直線與直線y=x交于Q2,如此繼續(xù)做下去可得到點(diǎn)Pk趨近于P*點(diǎn),其橫坐標(biāo)xk→x*。
(3)牛頓迭代法也稱切線法,對于任意x∈(x*,δ),若 f′(x)≠0,初選值 x0∈(x*,δ)做迭代,xk+1=xk-f(xk)/f′(xk),(k=0,1,2...)
牛頓迭代法的幾何意義:在曲線y=f(x)上從初值所對應(yīng)的點(diǎn)P0(x0,f(x0))做切線,以切線與x軸的交點(diǎn)x1=x0-f(x0)/f′(x0)作為f(x)=0的根x*的近似值,反復(fù)重復(fù)這個過程,得到點(diǎn)列{xk},xk趨近x*,其收斂速度不低于二階。當(dāng)初值x0離x*較近時,收斂速度很快,但當(dāng)x0離x*較遠(yuǎn)時,有時會出現(xiàn)發(fā)散的情況,此時可以采用牛頓下山法,xk+1=xk-βf(xk)/f′(xk),(k=0,1,2...),β 稱為下山因子,為了保 證滿足|f(xk+1)|<|f(xk)|,依次取β=1,1/2,1/22,...,試算直到上式成立,再轉(zhuǎn)向下一步迭代,牛頓下山法的優(yōu)點(diǎn)是對初值x0的選取無限制,即使使用牛頓法不收斂的x0,運(yùn)用下山法幾次迭代后,其迭代值進(jìn)入收斂區(qū)域,但下山法不再具有平方收斂的速度。
(4)對于一個復(fù)雜的 f(x)函數(shù),當(dāng)求導(dǎo)困難時,可以采用割線法,割線法分為單點(diǎn)割線法和兩點(diǎn)割線法。將牛頓迭代公式中的f(x)′用差商來代替,令f(x)′≈(f(xk)-f(x0))/(xk-x0),得到單點(diǎn)割線 法 ,選 初 值 x0,x1,xk+1=xk-f(xk)/(f(xk)-f(x0)),(k=1,2...);若令 f(x)′≈(f(xk)-f(xk-1))/(xkxk-1),則得到兩點(diǎn)割線法,選初值x0,x1,xk+1=xk-f(xk)/(f(xk)-f(xk-1))。
單點(diǎn)割線法是用過點(diǎn)(x0,f(x0))和(xk,f(xk))兩點(diǎn)的割線與x軸交點(diǎn)的橫坐標(biāo)xk+1作為x*的近似值。選初值x0,函數(shù)在f(x1)f(x0)<0,則方程在區(qū)間(x0,x1)內(nèi)至少有一個實(shí)根,單點(diǎn)割線法要求f′(x)和 f"(x)在區(qū)間內(nèi)不變號,而對于兩點(diǎn)割線法則無此要求。兩點(diǎn)割線法的幾何意義是用過點(diǎn)(x0,f(x0))和(x1,f(x1))兩點(diǎn)的割線與x軸交點(diǎn)的橫坐標(biāo)x2作為x*的新近似值,重復(fù)此過程,用過點(diǎn)(xk,f(xk))和(xk-1,f(xk-1))的割線與x軸交點(diǎn)的橫坐標(biāo)xk+1作為x*的下一近似值。
單點(diǎn)割線法的迭代函數(shù),滿足 0<|g′(x*)|<1,是局部收斂且為線性收斂;雙點(diǎn)割線法的收斂階p=1.618...,具有較快的收斂速度,但隨著xk→x*時,(fxk)-(fxk-1)不斷減小,有效位數(shù)大大減少,用它做分母將導(dǎo)致很大的誤差,從而影響計算結(jié)果。
3.3.3 用牛頓下山法計算柯爾勃洛-懷特公式
數(shù)值計算中選用的算法不同,其計算的效率也不相同,選用效率高的算法不僅可以提高解題速度,而且可以減少誤差產(chǎn)生的機(jī)會,減少誤差的積累,同時考慮到初值選取對計算的影響,本文采用牛頓下山法計算柯爾勃洛-懷特公式。
將公式 1/λ1/2=-2log[2.51/Re λ1/2+Δ/3.72di]寫為函數(shù)f(λ)=0的形式,得到方程
f(λ)=1/λ1/2+2log[2.51/(Reλ1/2)+Δ/3.72di]=0(5)式中:f(λ)中λ為自變量,其余均為常量,方程f(λ)=0有實(shí)根時λ的取值區(qū)間為(0,∞),函數(shù)f(λ)在定義域(0,∞)內(nèi)連續(xù),令 x=1/λ1/2,a=2.51/Re,b=Δ/3.72di,則函數(shù) f(λ)變?yōu)?f(x)=x+2log(ax+b),對函數(shù)f(x)求導(dǎo)得到
由于變量x及常量a、b均大于0,則在定義域內(nèi)f(x)′>0,f(x)"<0,函數(shù)在(0,∞)上的圖形是凸的。由于函數(shù)f(x)單調(diào)遞增,則函數(shù)與橫軸x交點(diǎn)唯一,方程f(x)=0僅有單根。迭代公式為
依次取β=1,1/2,1/22,...試算,直到|f(xk)|<|f(xk-1)|,之后繼續(xù)由迭代公式進(jìn)行計算,直到|xk+1-xk|滿足絕對誤差限ε的要求。計算出方程的近似解 x*=xk+1后,再由 x=1/λ1/2算出 λ。
用自然語言描述牛頓下山法算法如下:
第一步,取值 x0,ε,N(ε 為絕對誤差限,N 為設(shè)定的最大迭代次數(shù));
第二步,令0≥k;
第三步,令1≥β;
第四步,當(dāng)k≤N時,計算f′(x0);
第五步,判斷 f′(x0)=0 是否成立,若是,則 x0為解,輸出 x*=x0,停止;若不是,輸出 f′(x0);
第六步,計算f(x0);
第七步,計算x=x0-βf(x0)/f′(x0);
第八步,計算f(x);
第九步,判斷|f(x)|<|f(x0)|是否成立,若是,則輸出x;
若不是,則β/2≥β,返回第七步;
第十步,判斷|x-x0|<ε 是否成立,若是,則 x 為解,輸出 x*=x,停止;若不是,則 k+1≥k,x≥ x0,返回第三步;
第十一步,已迭代N次仍未達(dá)到精度要求,停止;
上述算法計算時需要先人工計算出函數(shù)的導(dǎo)數(shù)f(x)′,之后再由C語言、FORTRAN語言等進(jìn)行編程計算,如果使用Matlab軟件編程計算,可由程序語句直接求得函數(shù)導(dǎo)數(shù)f(x)′,使得復(fù)雜函數(shù)的編程計算更加便捷,在缺少以上工具時也可使用Excel軟件進(jìn)行計算。
根據(jù)該工程已知條件,v=1.711 m/s,di=0.5 m,t=12℃,ε=0.00005
計算可得 γ=1.236×106 m2/s,Δ=0.00001,Re=692150.107;
使用Excel軟件進(jìn)行計算,當(dāng)x0取初值x0=5時,經(jīng)過3次迭代就得到了滿足絕對誤差限ε要求的|x3-x2|,3.64×10-5<0.00005,此時 λ=0.01276,管道沿程水頭損失為每公里3.8 m,局部損失為沿程損失的10%,總水頭損失為每公里4.18 m。
根據(jù)該工程水力計算結(jié)果,輸水管道起點(diǎn)k0+000(高程為 2170.300 m),終點(diǎn) k13+640(高程為2017.600 m),兩地高差150.7 m,需要在k3+910(高程為 2104.300 m),k11+310(高程為 2034.800 m)處分別設(shè)置減壓水池一座,共計兩座減壓水池,k3+910處末端靜水壓為49.656 m,k11+310處末端靜水壓為38.568 m,k13+640,k13+640處末端靜水壓為7.461 m,可滿足管道內(nèi)最大使用壓力不超過設(shè)計值,同時滿足輸水所需壓頭。為了消除水錘的危害,在管道沿線高點(diǎn)以及直線管段每隔約1 km設(shè)置彌合性水錘預(yù)防閥一個。
輸水管道在條件允許時采用有壓重力流輸水,通過計算合理選擇管材及管徑,并且根據(jù)管道公稱壓力在適當(dāng)位置設(shè)置減壓池,具有投資省、運(yùn)行成本低、管理方便等諸多優(yōu)點(diǎn),是一種理想的輸水方式。