拱壩伏格特地基變位系數公式推導修訂與驗證

劉國華，黃芯怡，王 昱，苗智超

(浙江大學建筑工程學院，浙江 杭州 310058)

0 引 言

拱壩結構應力分析是拱壩設計和評價中的一個關鍵環節[1]，拱壩應力分析的常用方法有圓筒法、純拱法、拱梁分載法(含拱梁法和拱冠梁法)、殼體理論計算方法、有限單元法[2]和結構模型試驗法等[3]，其中拱梁分載法[4-5]仍然是現階段我國拱壩應力分析的主要方法[6-7]。拱梁分載法一般都采用挪威學者伏格特(F.Vogt)[8]于100年前提出的伏格特變位計算假定及相應的伏格特地基變位系數計算公式來計算基礎變位，其中δ′的計算公式是由美國墾務局補充的[9]。一般拱壩文獻只列出伏格特地基計算假定和計算公式，缺少具體的推導過程。

美國墾務局在《壩論》(1948年)中給出了基于伏格特假定的伏格特地基變位系數的計算公式，在《拱壩設計》(1976年)中對其中變位系數的計算方法做了修正(以下簡稱“墾務局文”，相應公式簡稱“墾務局公式”)[10]。中國水利水電科學研究院汪景琦先生在1965年出版的《拱壩的設計和計算》中，較完整地列出了伏格特變位系數的簡要推導過程和最終計算公式(以下簡稱“汪景琦文”，相應公式簡稱“汪景琦公式”)[9]。東北勘測設計研究院陳玉夫先生1996年6月在《東北水利水電》“伏格特理論的研究”一文(以下簡稱“陳玉夫文”，相應公式簡稱“陳玉夫公式”)中，列出了墾務局的伏格特系數計算公式及日本編寫的《在電子計算機上計算拱壩的解析分析》中的計算公式(以下簡稱“日本公式”)，并指出了其中存在的一些差異；同時，陳玉夫文也對汪景琦文中的推導過程提出了一些質疑[11]。中國水利水電科學研究院朱伯芳院士在2002年出版的《拱壩設計與研究》中，簡要列出了伏格特變位系數的主要推導過程和最終計算公式[12](以下簡稱“朱伯芳文”，相應公式簡稱“朱伯芳公式”)[13]。

上述多個出處的伏格特變形系數計算公式不完全一致。由于伏格特變形系數為三重積分式，被積函數比較復雜，推求顯式解析表達式的求積過程十分繁復，存在很大難度。各種出處均缺少詳細的中間推導過程，不清楚是否存在推導錯誤或是否采用了某些近似假定才得以導出顯式的解析表達式。為此，本文借助于Matlab工具，結合一些推導技巧，從伏格特基本假定和波辛涅斯克Boussinesq[14]理論解出發，重新對伏格特變位系數的顯式解析表達式進行推導，并采用高精度的數值積分手段，對顯式解析表達式的正確性進行驗證。

1 伏格特系數公式推導

1.1 伏格特系數公式的理論基礎

伏格特公式是在波辛涅斯克Boussinesq關于平面無限地基的變形公式的基礎上求得的基礎表面受矩形荷裁的平均變形方程式。半空間表面受力示意如圖1所示。根據波辛涅斯克Boussinesq理論解，在半無限空間表面受法向集中力P(見圖1a)，在y=0表面點(x，0，z)上的位移為[15]

圖1 半空間表面受力示意(波辛涅斯克Boussinesq理論解作用力)

(1)

式中，ux、uz、uy分別為x、z、y方向上的位移；μ為泊松比；E為彈性模量。

半無限空間表面受z方向的切向集中力S(見圖1b)，根據波辛涅斯克Boussinesq理論解，在y=0表面點(x，0，z)上的位移為

(2)

1.2 α′計算公式的推導

拱壩基礎變形系數(伏格特系數)共有7個，其中α′為由基礎表面上單位彎矩(單寬Mx=1)產生的在x=0、y=0處沿壩厚方向(z=-0.5a～0.5a)的平均角變位(繞x軸)。推導α′的作用力示意如圖2所示。

圖2 推導α′的作用力示意

沿x方向均勻分布且單寬彎矩Mx=1作用下，(x=0，y=0，z=z′)處的鉛錘向y方向的位移為v(z′)，則

(3)[1]

位移v(z′)的點位為(x=0，y=0，z=z′)，該點位相對于(x，y=0，z)(微面積dx×dz上的荷載作用點)的相對坐標為(Δx=-x，Δz=z′-z)，固有

(4)

代入得

(5)

(6)

應用Matlab工具，直接對上式進行積分，無法得出解析表達式。先就積分式的最里層，利用Matlab工具對ξ積分，得到

(7)

將上式中ηζ分為兩部分，即ηζ=(η-ζ)ζ+ζζ。第一部分(η-ζ)ζlog(...)和第二部分ζ2log(...)，利用Matlab工具分別對η積分，求和得到

(8)

對上述被積函數中的五項分別拆分后，利用Matlab工具分別進行積分，求和并化簡，再輔以人工進一步化簡后，得到

(9)

其余6個伏格特系數公式的推導過程本文不再詳述。

1.3 伏格特系數公式推導結果及比較

以下列出目前主流的伏格特系數公式及本文伏格特系數公式的推導結果。

(1)α′為基礎表面上單位彎矩(Mx=1)產生的平均角變位。墾務局、汪景琦、日本和朱伯芳公式

(10)

本文和陳玉夫公式

(11)

(2)β′為基礎表面上單位法向力(Ny=1)產生的平均法向變位。墾務局、汪景琦和朱伯芳公式

(12)

日本公式

(13)

本文和陳玉夫公式

(14)

(3)γ′為基礎表面上單位徑向剪力(Qz=1)產生的平均徑向變位。墾務局、汪景琦、日本和朱伯芳公式

(15)

本文和陳玉夫公式

(16)

(17)

本文和陳玉夫公式

γ′=(1+μ)Eπ2msinh-12m()+2sinh-1m2(){-mm2+42+m22+μmm2+4-2msinh-12m()-m2ù?úúé?êê}

(18)

(5)δ′為基礎表面上單位扭矩(My=1)產生的平均扭轉角變位。墾務局、汪景琦和朱伯芳公式

(19)

日本公式

(20)

本文和陳玉夫公式

(21)

(6)α″為基礎表面上單位徑向剪力(Qz=1)產生的平均轉角變位。墾務局、汪景琦、日本和朱伯芳公式

(22)

本文和陳玉夫公式

(23)

(7)γ″為基礎表面上單位彎矩(Mx=1)產生的平均徑向變位。墾務局、汪景琦、日本和朱伯芳文中γ″與α″公式相同。陳玉夫公式

γ″=0

(24)

本文γ″與α″公式相同。

2 伏格特系數公式驗證

伏格特系數積分表達式在多個來源中是一致的，但在不同的來源中，基于同一積分式推求得到的顯式解析表達式有所不同。

為了驗算數值積分中積分步長對積分結果的影響，分別取1/10 000、1/5 000和1/200等3種積分步長進行計算，結果見圖4。從圖4可以看出，前2種的積分步長很小，計算結果很接近。而1/200的積分步長雖然已屬相當小的小步長，但仍存在一定的數值誤差。主要原因是被積函數存在數值趨于無窮大的奇異點，容易產生數值誤差。以積分步長1/10 000為例，圖4中一個驗算點的三重數值積分的被積函數計算次數就高達1萬億次(考慮被積函數對稱性后計算次數可有所減少)。好在現在計算機的性能比較好，高性能微機在幾十分鐘時間內即可完成圖3所示的全部伏格特系數公式的對比計算。

3 不同來源伏格特系數差異對拱壩應力分析結果的影響

為比較不同的伏格特系數對拱壩應力分析的影響，本文使用ADAO軟件，對兩座壩高200～300 m級特高拱壩進行應力分析，其中伏格特系數分別采用傳統公式(主要來源于美國墾務局)和本文推導的公式進行計算；并對徑向位移、切向位移、豎向位移、豎向扭轉位移及切向轉角位移5種壩基位移分量和壩踵主拉應力、壩踵主壓應力、壩趾主拉應力及壩趾主壓應力4種壩基應力進行差異對比，分析兩種伏格特系數對拱壩應力分析的影響。

計算得出，兩座特高拱壩算例因兩種伏格特系數差異導致的壩基線位移最大誤差大致為0.01～0.22 cm，相對誤差大致為-5.5%～6.5%；壩基角位移最大誤差大致為0.1×10-5～1.2×10-5rad，相對誤差大致為-1.3%～2.5%；壩體線位移最大誤差大致為0.01～0.31 cm，相對誤差大致為-2.3%～3.2%。

兩座特高拱壩算例因兩種伏格特系數差異導致的壩基結點壩體應力最大誤差大致為0.01～0.22 MPa，相對誤差大致為-0.2%～-1.4%；拱冠梁應力最大誤差大致為0.01～1.63 MPa，相對誤差大致為-0.73%～0.35%。

算例結果表明，如果以本文修訂后的計算結果為基準，來源于墾務局傳統伏格特系數公式所得的計算結果，壩基變位的最大相對誤差約為6.5%，壩基結點壩體應力的相對誤差約為0.5%。

4 結 論

本文借助于Matlab工具，從伏格特基本假定和波辛涅斯克Boussinesq理論解出發，對伏格特變位系數重新進行推導。將本文導出的顯式解析表達式與陳玉夫公式，墾務局、汪景琦、日本、朱伯芳之公式進行對比，并采用高精度的數值積分手段，對解析表達式正確性進行驗證。對比得出：全部7個伏格特系數，數值積分結果與本文顯式解析表達式一致；數值積分結果與陳玉夫的6個伏格特系數解析表達式一致，1個伏格特系數不一致；墾務局、汪景琦和朱伯芳顯式解析公式基本一致，日本拱壩分析的β′、、δ′等3個基礎變形系數的顯式解析公式與墾務局、汪景琦和朱伯芳顯式解析公式不同，這些顯式解析表達式均與數值積分結果存在一定差異。通過兩座特高拱壩作為算例進行應力分析得出，修訂前后的應力、位移計算結果差異不大。修訂后的伏格特公式的正確性得到了可靠的驗證，使伏格特方法更加完善。