拱坝伏格特地基变位系数公式推导修订与验证

刘国华，黄芯怡，王 昱，苗智超

(浙江大学建筑工程学院，浙江 杭州 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个基础变形系数的显式解析公式与垦务局、汪景琦和朱伯芳显式解析公式不同，这些显式解析表达式均与数值积分结果存在一定差异。通过两座特高拱坝作为算例进行应力分析得出，修订前后的应力、位移计算结果差异不大。修订后的伏格特公式的正确性得到了可靠的验证，使伏格特方法更加完善。