陆 军，武清玺

（河海大学 土木工程学院，南京 210098）

拱坝是一种重要的坝型，以其结构合理和体形优美而著称，其所有的优良性能建立在拱坝稳定可靠的基础上。坝体一旦失事，将严重影响工农业生产，因此拱坝的可靠度研究和安全分析早已得到人们的重视。本文将考虑基岩材料的非线性性质，研究非线性方程的求解方法，并将材料的力学特性、拱坝受到的外来作用等作为随机变量，以有限元计算的结果为基础，利用二次响应面法计算拱坝结构重要部位单元的可靠度，通过整理和分析得出坝体可靠指标的变化规律。

1 结构可靠度计算的响应面法

目前，结构的极限状态方程一般都基于抗力－荷载效应模型。现有可靠度计算方法多是以极限状态方程具有明确的解析表达式为基础的。但是对于某些复杂结构系统，由于结构本身的复杂性，基本随机变量的输入与输出量之间的函数关系常常是高度非线性的，有时甚至得不到明确的解析表达式，在计算这类复杂结构的可靠度时，采用一次二阶矩法、JC法等就存在一定的困难。Box和Wilson于1951年首先提出响应面法，当时对于响应面法的研究仅限于如何用统计的方法得到一个显式函数，用来逼近一个复杂的隐式函数。直到1985年，文献［1］首先用一次响应面法分析研究了土坡稳定的可靠度问题，而后文献［2－4］研究才逐步拓展了这一领域。响应面法是统计学的综合试验技术，用于处理复杂系统的输入（随机变量）与输出（系统响应）的转换关系问题。该方法采用有限的试验，通过回归拟合解析表达式Z′＝g′（X）代替真实曲面Z＝g（x）。

在可靠度计算中，设极限状态方程为：

该方程将n维随机空间分为失效区和安全区两部分。结构的失效概率为：

式中 fX（x1，x2，…，xn）为n维随机向量的联合概率密度函数。在多数实际问题中，不但极限状态方程是非线性的，而且n维随机向量的联合概率密度函数也难以得到，所以精确求解式（2）非常困难。对于大型、复杂的工程结构，可靠度计算常采用数值方法，如随机有限元法等。由于基岩的本构关系是非线性的，利用随机有限元法求解要经过多次迭代，计算繁琐且结果的误差也较大，难以满足精度要求。近年来得到迅速发展的结构可靠度计算的响应面法，可以应用确定性有限元计算程序对结构进行可靠度分析，使这一方法展现出良好的应用前景。对于非线性极限状态方程，采用在设计验算点处的一次Taylor展开式求得结构的可靠指标的精度不理想。二次可靠度分析法的计算精度虽有所提高，但工作量成倍增加。结构可靠度计算的二次响应面法与一次二阶矩法相比有2个优点：①用一个待定的显式响应面函数逼近实际的隐式极限状态函数，使可靠度计算得到简化；②计算结构可靠指标时以二次响应面取代一次近似式，使计算精度有所提高。

设二次响应面函数为：

式中 a，bi，ci为待定系数，需由2n＋1个样本点确定。 这一方法的计算过程为：第1步，以均值μX点为中心，样本点在区间（μX-fσX，μX+fσX）中选取，建议f的取值为1～3；用g（X）在2n＋1个样本点处的函数值，可以确定响应面函数中的待定系数；响应面函数得到以后，即可求出原极限状态面上设计验算点的近似值XD。第2步，选取新的样本中心点，新的样本中心点XM选在由均值点mX到XD的直线上，并保证在新的中心点处满足：

这样选取新中心点的目的是为了使所选样本点包含原极限状态面更多的信息。第3步是以为中心点选取新的一组样本点，重复上述过程，即可得到新的响应面函数，并近似求出相应的可靠指标。整个过程需求出4n＋3个g（X）值，这一方法的逼近过程如图1所示［5-8］。

2 非线性有限元

2.1 材料的破坏准则与本构模型

基于剪压破坏的理想弹塑性材料模拟基岩，屈服准则采用Mohr－Coulomb准则［9］：

坝体混凝土材料采用弹性模型，采用拉裂压碎破坏准则：

式中 c，φ分别为材料的抗剪断粘聚力和内摩擦角；σ1，σ3分别为单元的最大、最小主应力；ft，fc分别为材料的抗拉、抗压强度的设计值。

2.2 非线性有限元求解方法

考虑岩体为弹塑性材料，材料的性质与应力和变形的历史有关，本构方程用增量形式表达。按载荷作用的实际情况划分为若干个子步，在小的载荷增量下逐步计算。在某一载荷增量{Δδ}n的作用下，有限元的平衡条件为：

式中 {ΔR}n为载荷增量，求解式（8）得位移增量{Δδ}n。 在载荷增量{ΔR}n施加之前，已作用有累积载荷{R}n-1，结构中的应力、应变和位移均已知道，分别为{σ}n-1，{ε}n-1和{δ}n-1。 那么在新的累积载荷{R}n={R}n-1+{ΔR}n作用下，其总的应力、应变和位移分别见（9）式［10］：

3 工程算例分析

3.1 工程概况

福建省永泰县青龙溪大洋水库大坝坝型设计为实体混凝土三心变厚双曲拱坝，坝顶高程为236.0m，最大坝高87.2m，大坝为Ⅲ等3级建筑物。水库的正常蓄水位230m，校核洪水位236.4m。工程区位于区域地质构造相对稳定地段，未见区域性的断裂构造通过坝址，河槽段无断层构造通过，局部节理裂隙发育地带采用加密固结灌浆孔工程处理。

3.2 计算模型

本次计算采用三维模型进行计算，坐标系以坝顶拱冠梁上游点为坐标原点，左岸为X轴正方向，下游为Y轴正方向，竖直向下为Z轴正方向。模拟范围为上游0.8倍坝高，下游1.5倍坝高，左右两岸1.5倍坝高，坝基深度为1倍坝高，坝顶高程以上模拟20m山体。

网格的划分依据福建省水利水电科学研究所提供的拱坝勘察设计资料进行，网格采用八节点六面体单元，节点总数152208，单元总数为141064，其中坝体单元总数为5244。离散后的整体模型如图2所示。计算中对坝基岩体材料应力应变关系采用基于增量法的弹塑性模型，采用对岩体强度参数进行折减的方法分别模拟右岸和坝体河槽段岩体的两条软弱带如图2，坝体混凝土材料采用弹性模型。

3.3 计算参数及计算结果

由地质勘查资料，坝基网格根据材料物理参数的不同分为河槽段的岩石1和两岸的岩石2单元组；右岸和河槽段的软弱带分别定义为软弱1和软弱2单元组；坝体单元为混凝土单元组。

计算拟定的荷载组合为：自重＋校核洪水位（236.4m）＋设计温升（升温温度场与坝体稳定温度场的差值）。为了模拟大坝施工和逐级蓄水的过程，采用了分级加载的方式。其中，坝体重力荷载分6级加载；水荷载分9级加载；温度荷载分6级加载。不考虑地应力，仅计入岩体自重所形成的初始应力场，采用定值考虑的材料计算参数如表1所示。

现对各随机变量取均值时的计算结果加以分析。（应力以拉为正，压为负）坝体上、下游面单元的安全性是最为关心的，所以本文首先分析了这些区域的应力计算结果，如图3所示。由图3（a）、（b）可以看出，坝体最小主应力分布规律是：上游面单元由拱冠梁处向两侧逐渐变小，由顶部向下部先逐渐变大再逐渐变小，在距坝顶约1/3坝高处，压应力达到最大值；下游面单元的最小主应力由拱冠梁处向两侧逐渐变大，并在临近坝肩处达到最大值。由于坝体形状基本上左右对称，所以得到的应力也大致为对称分布。坝体拉应力分布区域较小，主要分布在坝踵附近，拉应力最大值为1.5MPa。

表1 材料计算参数表

图3 拟定荷载组合下拱坝应力云图

3.4 拱坝可靠度计算及分析

设计人员总是力求使混凝土设计参数取值合理、可行，尽可能使混凝土具有优良的物理、力学性能。拱坝建成运行后上游面的水压力、混凝土水化热、上游面库水和下游面气温的温差都是坝体的主要荷载。其中水荷载作为主要荷载是无法避免的，温度荷载可以通过精心组织施工减小其影响，所以本文将坝体上游水位、基岩、混凝土的某些物理力学参数作为随机变量考虑，对可靠度计算结果影响甚微和缺乏统计资料的参数作为定值考虑。计算中采用的随机变量的特征值如表2所示。

表2 计算中所采用的随机变量统计特征值

坝体混凝土的抗裂或抗压可靠度极限状态方程可表示为：

式中 R为坝体混凝土的抗裂强度或抗压强度，分别对应通过非线性有限元计算出的混凝土单元的第一和第三主应力；S为荷载效应，是关于随机变量E1，E2，E3，tg准，c，Hu的隐式函数。右岸岩体的薄弱带基于Mohr－Coulomb屈服准则的可靠度极限状态方程可表示为：

式中 c，准分别为右岸岩体薄弱带的粘聚力和摩擦角；σ1，σ3分别是通过非线性有限元计算出的该薄弱带单元的第一和第三主应力，是关于随机变量E1，E2，E3，tg准，c，Hu的隐式函数。

3.5 计算结果

利用有限元计算的应力结果，再结合响应面法可以得到坝体单元的可靠指标，得到足够多单元的可靠指标后再绘制出可靠指标的等值线图，比较直观的反映出坝体结构的重要部位的可靠度情况。以48号单元为例给出某一单元的可靠指标计算过程。

（2）响应面函数确定后再利用梯度优化法求出极限状态面上设计验算点的近似值XD为（14.9943，20.0042，19.9393，0.4659，0.7001，74.5734，1.7972）。 按照上文介绍的方法选取新的中心点XM（14.9929，20.0053，19.9243，0.4658，0.7002，76.9225，1.7966）。

（3）以新的中心点XM为中心开始选取新的一组样本点，重复第一步的工作，即可求得极限状态面上设计验算点的值和相应的可靠指标 β=4.5。

本文取坝体主要区域的320个单元计算其可靠指标，绘制出抗压可靠指标的等值线图如图4。该图清晰地反映了拱坝结构的可靠指标分布情况。拱坝上游面拱冠梁自由端附近的单元、下游面坝肩附近的坝体单元可靠指标均较小。从图4a、图4b图看出坝体单元抗压可靠指标大小的变化趋势与其主应力变化趋势与图3相同。由于拱坝结构内部以压应力为主，拉应力主要分布在建基面附近。拉应力的最大值出现在坝踵处，此处单元的抗拉可靠指标的最小值为1.1，需采取措施增大材料的抗拉能力。软弱1和软弱2单元组中单元的可靠指标均在5.0以上，说明在给定的荷载下是比较安全的。

4 结语

（1）混凝土拱坝及坝基材料均具有多样的物理、力学特性，其离散性强、变异性大，本构关系是非线性的。在对其作可靠度分析时，直接得到系统极限状态方程的解析表达式往往是非常困难的，并且联合概率密度函数也难以得到。响应面法结合有限元分析方法，为解决类似工程的可靠性分析提供了可供选用的一条新途径。

（2）可靠度计算时采用响应面法解决拱坝类复杂结构问题稳定性较好，计算可靠指标的梯度优化法收敛速度快，一般仅需3～5次迭代计算。

（3）采用商用有限元计算软件结合自编的FORTRAN程序实现结构可靠指标的计算方法在PC上解决本论文的计算具有实用价值。

［1］Wong F S.Slope reliability and response surface method［J］.Journal of Geotechnical Engineering，1985，b111（1）：32-53.

［2］Liu Y W，Moses F.A sequential response surface method and its application in the reliability analysis of aircraft structural systems［J］.Structural Safety，1994，16：39-46.

［3］武清玺，卓家寿.结构可靠度分析的变f序列响应面法及其应用［J］.河海大学学报，2001，29（2）：75-78.

［4］武清玺，俞晓正.混凝土面板堆石坝可靠度计算方法研究［J］.岩土工程学报，2004，26（4）：468-472.

［5］武清玺.结构可靠性分析及随机有限元法［M］.北京：机械工业出版社，2005.

［6］结构可靠度理论［M］.北京：中国建筑工业出版社，2005.

［7］吴世伟.结构可靠度分析［M］.北京：人民交通出版社，1990.

［8］熊铁华，常晓林.响应面法在结构体系可靠度分析中的应用［J］.工程力学，2006，23（4）：58-61.

［9］李咏偕，施泽华.塑性力学［M］.北京：中国水利水电出版社，1987.

［10］朱伯芳，有限单元法原理与应用［M］.北京：中国水利水电出版社，1998.