混凝土拱坝坝体优化方法

权利要求书： 1.一种混凝土拱坝坝体优化方法，其特征在于：包括以下步骤：步骤1：选择用以确定拱坝体型的基本参数以及确定坝体厚度的基本参数，作为优化设计变量；所述拱坝体型的基本参数包括拱冠梁处拱圈中心线的曲率半径和拱端中心角；所述坝体厚度的基本参数包括拱冠梁厚度、拱端厚度；

步骤2：选择适当的目标函数，以安全性为关键性指标，即以应力为目标函数；

步骤3：为保证工程安全和施工方便，确定坝体几何尺寸和工作性态两方面的约束，包括几何约束条件、输入应力约束条件和稳定约束条件；

步骤4：将拱坝的体型抽象成几何模型并用一组参数来完整地描述，基于快速多级边界单元进行坝体网格自动化划分以及荷载的自动施加；

步骤5：基于快速多级边界元法自动化计算应力，利用等效应力法求出拱坝上下游面的三个应力分量，然后由该处微分体的平衡条件求得另外三个应力分量，进一步决定拱坝上下游面的主应力，进而确定出最优拱坝坝形。

2.根据权利要求1所述的混凝土拱坝坝体优化方法，其特征在于：所述步骤3中，输入应力约束条件包括坝轴线位置移动范围的限制、坝顶最小厚度的限制、倒悬度约束和保凸性要求。

3.根据权利要求1或2所述的混凝土拱坝坝体优化方法，其特征在于：所述步骤3中，输入几何约束条件为坝体内部的最大主拉应力及最大主压应力均小于坝体混凝土的允许强度，即

|σt|max≤[σt],|σc|max≤[σc](1)；

式中：[σt]为容许拉应力，[σc]为容许压应力。

4.根据权利要求1或2所述的混凝土拱坝坝体优化方法，其特征在于：所述步骤4中，输入拱坝几何模型参数包括拱冠梁和水平拱圈的参数化描述；所述拱冠梁参数包括冠梁中心线曲线方程Yc(z)和拱冠梁厚度Tc(z)；将Yc(z)和Tc(z)都设为z的3次方程，故只需有4个控制高程(z＝z1、z2、z3、z4)处对应的Yc和Tc便可按插值方法确定Yc(z)和Tc(z)方程中的系数；拱冠梁剖面用以下12个参数来描述：上述Yci和Tci分别为控制高程zi处的拱冠梁中心线y坐标及拱冠梁厚度；水平拱圈只要确定了拱圈中心线的方程及拱冠至拱端的厚度变化规律便能确定水平拱圈的形状。

5.根据权利要求3所述的混凝土拱坝坝体优化方法，其特征在于：所述步骤4中，输入拱坝几何模型参数包括拱冠梁和水平拱圈的参数化描述；所述拱冠梁参数包括冠梁中心线曲线方程Yc(z)和拱冠梁厚度Tc(z)；将Yc(z)和Tc(z)都设为z的3次方程，故只需有4个控制高程(z＝z1、z2、z3、z4)处对应的Yc和Tc便可按插值方法确定Yc(z)和Tc(z)方程中的系数；拱冠梁剖面用以下12个参数来描述：上述Yci和Tci分别为控制高程zi处的拱冠梁中心线y坐标及拱冠梁厚度；水平拱圈只要确定了拱圈中心线的方程及拱冠至拱端的厚度变化规律便能确定水平拱圈的形状。

6.根据权利要求1或2或5所述的混凝土拱坝坝体优化方法，其特征在于：所述步骤5中，利用快速多级边界元法自动化计算应力，具体公式为：cij(x)ui(x)＝∫Γ[Uij(x,y)tj(y)?Tij(x,y)uj(y)]dΓ(y)(3)；

其中，x,y为研究域内的点，i,j＝1,2,3，x, 2π·cij(x)为点x在边界内部的角，当点x所在的边界处平滑时，有cij(x)＝0.5；Γ为研究域的边界，ui(x)和ti(y)为位移和面力，以及Uij(x,y)和Tij(x,y)为开尔文基本解，具体表达式为：其中i,j＝1,2,3，r为点x和y的距离， 代表r关于xi的偏导数，以及

7.根据权利要求6所述的混凝土拱坝坝体优化方法，其特征在于：所述步骤5中，利用的快速多级加速边界元法，快速多级算法的具体计算步骤为：(5.1)对求解域进行计算域边界的离散；

(5.2)采用最小的立方体包含求解域，将其命名为根立方体，也就是第0层立方体，设定每个叶子立方体中最多包含的节点个数，构造八叉树结构；

(5.3)指定多级展开系数，对于每个层数大于1的立方体，定义存储多极矩的矩阵；对于每个层数大于2的立方体，定义存储局部展开系数的矩阵；对于每一个叶子立方体，定义矩阵存储系数；定义矩阵存储每个叶子近点直接计算的系数；

(5.4)对于每个叶子立方体分别计算其相关系；

(5.5)计算每个叶子立方体的多极矩；

(5.6)对于处于第L层每一个非叶子立方体计算其多极矩，这个多极矩由其子立方体传递而来，这个过程称为向上传递；

(5.7)对于第L层的每一个立方体，其中L从2到最终层，通过其相互作用列表中的立方体传递而来；

(5.8)分别计算近点和远点的贡献系数；

(5.9)更新迭代向量；

(5.10)不断重读以上过程直到达到迭代算法的控制残差，进而求得应力和位移的未知量。

8.根据权利要求1或2或5或7所述的混凝土拱坝坝体优化方法，其特征在于：所述步骤5中，等效应力法为：

τtz＝τzxcosα+τzysinα(7)；

式中：α为x轴到t轴的转角，以逆时针为正，为拉应力，σx,σy,σz,τxy,τyz,τzx代表直角坐标系—xyz中的应力分量；σt和τtz为空间坐标系—rtz中的应力分量；rtz中r轴指向拱轴线的法线，t轴指向拱轴线的切线，z轴铅直向上；

若求出拱坝某水平截面上沿坝厚三个单元的形心处的梁应力σz，记为σ1,σ2,σ3；将x轴设在此3个单元的连线上，已知x＝?1、0、1三点处的正应力σz，上下游处的等效应力：式中，轴力 和弯矩 为：

同理，求出σz,τzt对应的另两组等效应力；

利用求得的上下游处的等效应力，再根据坝面微分体的平衡条件，求出该处其它三个等效应力：

式中：p为坝面所受的水压力强度，以压为正；θ为坝面与r轴的交角； 为不同方向上的等效应力；

有了6个等效应力分量之后，便能计算得出等效应力的主应力。

9.根据权利要求8所述的混凝土拱坝坝体优化方法，其特征在于：所述步骤5中，根据下式，得3个实数解便为等效应力的3个主应力：

3 2

σ?I1σ+I2σ?I3＝0(14)；

其中，

式中，I1，I2和I3为应力不变分量；σx,σy,σz,τxy,τyz,τzx代表应力分量。

说明书： 混凝土拱坝坝体优化方法技术领域[0001] 本发明属于坝体优化技术领域，具体涉及一种混凝土拱坝坝体优化方法。背景技术[0002] 所谓拱坝体型优化设计，就是在已知的给定参数条件下，寻找满足工程的施工、运用、强度、稳定性等约束条件，并使目标函数取最小值的体型方案。因此，优化设计的首要一

步是将工程实际问题表示成数学表达式，即建立数学模型，包括选定设计变量、选择目标函

数、建立约束方程，然后采用适宜的优化方法，求出问题的解。

[0003] 关于拱坝体型优化的研究，国外是从20世纪60年代开始的。1969年，R.Sharp提出了用数学规划法对拱坝进行体型优化设计，从而使拱坝体型优化设计作为一个专门课题，

引起学术界和工程界的注意；1973年，加拿大学者W.Stensch提出了自由型拱坝设计的模

型；1974年G.A.Harley和G.M.Mcneice形成了有限元分析用的网格自由划分系统，进一步改

进了W.Stcnsch的工作；1975年R.E.Rhichetts和D.C.Liewicz发表了《混凝土坝的形状优化

设计》，其中也研究了拱坝；1982年，在上海举行的国际有限元会议上，K.Wassermann发表了

一篇拱坝体型优化的论文等。

[0004] 我国拱坝体型优化设计的研究是在20世纪70年代末期由朱伯芳院士提出来的，尽管起步较晚，但发展很快。我国拱坝优化设计研究已经进行了30多年，在优化理论和实践上

均取得了显著的成绩，目前我国在拱坝优化领域已处于国际领先地位。从基于材料力学优

化发展到有限元优化，从静力优化发展到动力优化，从线性优化发展到非线性全过程分析，

拱圈的线性从单心圆发展到二次曲线、混合线型等，逐步发展和完善了我国提出的高混凝

土坝优化理论与方法。

[0005] 目前，有限元法作为一种很受欢迎的数值分析方法，逐步成为坝型优化应力计算的主要方法，但一般而言，拱坝为超大规模水利工程，有限元分析需要划分大量的三维体网

格，计算量巨大，耗时长，需要高额的计算成本。

发明内容[0006] 针对现有技术存在的不足，本发明提出了一种混凝土拱坝坝体优化方法，基于快速多级边界元法，模型重建过程更为简单，同等自由度下，精度更高。

[0007] 为了解决上述技术问题，本发明提供的混凝土拱坝坝体优化方法，其特征在于：包括以下步骤：

[0008] 步骤1：选择用以确定拱坝体型的基本参数，如拱冠梁处拱圈中心线的曲率半径、拱端中心角，以及确定坝体厚度的基本参数如拱冠梁厚度、拱端厚度等作为优化设计变量；

[0009] 步骤2：选择适当的目标函数，一般以安全性为关键性指标，即以应力为目标函数；[0010] 步骤3：为保证工程安全和施工方便，确定坝体几何尺寸和工作性态两方面的约束，主要包括几何约束条件、应力约束条件和稳定约束条件；

[0011] 步骤4：将拱坝的体型抽象成几何模型并用一组参数来完整地描述，基于快速多级边界单元进行坝体网格自动化划分以及荷载的自动施加。

[0012] 步骤5：基于快速多级边界元法自动化计算应力，利用等效应力法求出拱坝上下游面的三个应力分量，然后由该处微分体的平衡条件求得另外三个应力分量，进一步决定拱

坝上下游面的主应力，进而确定出最优拱坝坝形。

[0013] 作为优选，所述步骤3中输入应力约束条件为坝轴线位置移动范围的限制，坝顶最小厚度的限制，倒悬度约束，保凸性要求。

[0014] 进一步地，所述步骤3中输入几何约束条件为坝体内部的最大主拉应力及最大主压应力均小于坝体混凝土的允许强度，即

[0015] |σt|max≤[σt],|σc|max≤[σc](1)[0016] 式中：[σt]为容许拉应力，[σc]为容许压应力。[0017] 进一步地，所述步骤4中输入拱坝几何模型参数包括拱冠梁和水平拱圈的参数化描述，拱冠梁参数主要包括冠梁中心线曲线方程Yc(z)和拱冠梁厚度Tc(z)，一般将Yc(z)和

Tc(z)都设为z的3次方程，故只需有4个控制高程(z＝z1、z2、z3、z4)处对应的Yc和Tc便可按

插值方法确定Yc(z)和Tc(z)方程中的系数。这样拱冠梁剖面便可用以下12个参数来描述：

[0018][0019] 上述Yci和Tci分别为控制高程zi处的拱冠梁中心线y坐标及拱冠梁厚度。水平拱圈只要确定了拱圈中心线的方程及拱冠至拱端的厚度变化规律便可确定水平拱圈的形状。

[0020] 进一步地，所述步骤5中的利用快速多级边界元法自动化计算应力，具体基本计算公式为：

[0021] cij(x)ui(x)＝∫Γ[Uij(x,y)tj(y)?Tij(x,y)uj(y)]dΓ(y)(3)[0022] 其中 2π·cij(x)为点x在边界内部的角，当点x所在的边界处平滑时，有cij(x)＝0.5。ui(x)和ti(y)为位移和面力，以及

[0023][0024] 其中i,j＝1,2,3，r为点x和y的距离， 代表r关于xi的偏导数，v为泊松比，μ为剪切模量，以及

[0025][0026] 进一步地，所述步骤5中的利用的快速多级加速边界元法，快速多级算法的具体计算步骤为：

[0027] (5.1)对求解域进行计算域边界的离散；(5.2)尽可能采用最小的立方体包含求解域，将其命名为根立方体，也就是第0层立方体，设定每个叶子立方体中最多包含的节点个

数，构造八叉树结构；(5.3)指定多级展开系数，对于每个层数大于1的立方体，定义存储多

极矩的矩阵。对于每个层数大于2的立方体，定义存储局部展开系数的矩阵。对于每一个叶

子立方体，定义矩阵存储系数。定义矩阵存储每个近点直接计算的系数；(5.4)对于每个叶

子立方体分别计算其相关系；(5.5)计算每个叶子立方体的多极矩；(5.6)对于处于第L层每

一个非叶子立方体计算其多极矩，这个多极矩由其子立方体传递而来，这个过程可以称为

向上传递。(5.7)对于第L层的每一个立方体，其中L从2到最终层，通过其相互作用列表中的

立方体传递而来；(5.8)分别计算近点和远点的贡献系数；(5.9)更新迭代向量；(5.10)不断

重读以上过程直到达到迭代算法的控制残差，进而求得应力和位移的未知量。

[0028] 进一步地，所述步骤5中的等效应力法为：[0029][0030] τtz＝τzxcosα+τzysinα(7)[0031] 式中：α为x轴到t轴的转角，以逆时针为正。[0032] 若求出拱坝某水平截面上沿坝厚三个单元的形心处的梁应力σz，记为σ1,σ2,σ3(如图1和图2所示)。将x轴设在此3个单元的连线上，已知x＝?1、0、1三点处的正应力σz，上下游

处的等效应力：

[0033][0034] 式中，轴力 和弯矩 为：[0035][0036][0037] 同理，可以求出σz,τzt对应的另两组等效应力。[0038] 利用求得的上下游处的等效应力，再根据坝面微分体的平衡条件(如图3和图4所示)，可以求出该处其它三个等效应力：

[0039][0040][0041][0042] 式中：p为坝面所受的水压力强度，以压为正，θ为坝面与r轴的交角。[0043] 有了6个等效应力分量之后，根据下式，得3个实数解便为等效应力的3个主应力。[0044] σ3?I1σ2+I2σ?I3＝0(14)[0045] 其中，[0046][0047] 式中，σx,σy,σz,τxy,τyz,τzx代表应力分量。[0048] 本发明的优点及有益效果如下：[0049] 本发明采用的快速多级边界元法，一方面基于边界元法，只需要对坝体表面进行二维网格划分，而有限元法却需要建三维网格，建模过程工作量更大，因此，只需要建二维

网格的边界元法模型重建过程更为简单快捷，计算量更小；另一方面将快速多极算法和边

界元法相耦合，进一步提高了计算效率。总而言之，本发明建立的以应力(安全性)为目标函

数的坝形优化求解模型，相较于当前的有限元坝体优化方法，在同等自由度的情况下，本发

明快速多级边界元法计算效率更快且精度更高，具有不可忽视的优势。

附图说明[0050] 图1是本发明等效应力位置俯视图；[0051] 图2是本发明等效应力计算示意图[0052] 图3是本发明坝面微分体的上游平衡视图[0053] 图4是本发明坝面微分体的下游平衡视图具体实施方式[0054] 为了便于本领域技术人员理解和实施本发明，下面将结合附图及实施例对本发明作进一步的详细描述，此外，此处所描述的实施示例仅用于说明和解释本发明，并不用于限

定本发明。

[0055] 实施例1：【以国内某碾压混凝土双曲拱坝为例】[0056] 采用本发明的混凝土拱坝坝体优化方法，具体步骤如下：[0057] 步骤1：选择用以确定拱坝体型的基本参数，如拱冠梁处拱圈中心线的曲率半径、拱端中心角，以及确定坝体厚度的基本参数如拱冠梁厚度、拱端厚度等作为优化设计变量。

在实施例中，选取拱冠梁坝顶、坝底厚度，上、下游面上部倒悬度，坝高，坝体方量，控制高程

最大中心角为基本参数；

[0058] 步骤2：选择适当的目标函数，本实施例以安全性为关键性指标，即以应力为目标函数；

[0059] 步骤3：为保证工程安全和施工方便，确定坝体几何尺寸和工作性态两方面的约束，主要包括几何约束条件、应力约束条件和稳定约束条件，应力约束条件为坝顶最小厚度

的限制，倒悬度约束，保凸性要求。具体来讲，坝顶最小厚度的限制根据防护、交通、施工、抗

屈折等方面的具体要求而定，在本实施例中坝顶最小厚度为9m；对于混凝土拱坝，其倒悬度

一般控制在0.3以内；由微分几何知，可微曲面上某一点临近的凸凹性，可由高斯曲率K来判

定，坝体曲面的保凸性即为要求曲面上每一点均有K>0。几何约束条件为坝体内部的最大主

拉应力|σt|max及最大主压应力|σc|max均小于坝体混凝土的允许强度，即

[0060] |σt|max≤[σt]，|σc|max≤[σc](1)[0061] 式中：[σt]为容许拉应力，[σc]为容许压应力。[0062] 在本实施例中，稳定约束条件为要求拱坝坝体在正常蓄水位+坝体自重+泥沙压力+正常温降的组合条件下保持稳定，具体来说，保持稳定的具体要求为大坝发生的滑移距离

在设计要求容许范围之内，并且不发生倾倒屈折。应力约束条件为上、下游面上部倒悬度，

拱冠梁坝顶、坝底厚度和控制高程最大中心角的最大容许值分别为：0.3，0.3，9.0，15.0米，

21.5米，110度。此外，根据《混凝土拱坝设计规范》，对于压应力，规范要求拱坝的设计最大

压应力不得大于容许压应力[σc]，在本实施例中，[σc]＝5.0MPa；对于压应力，规范要求坝体

主拉应力不得大于1.5MPa。

[0063] 步骤4：将拱坝的体型抽象成几何模型并用一组参数来完整地描述，基于快速多级边界单元进行坝体网格自动化划分以及荷载的自动施加。拱坝几何模型参数包括拱冠梁和

水平拱圈的参数化描述，拱冠梁参数主要包括冠梁中心线曲线方程Yc(z)和拱冠梁厚度Tc

(z)，一般将Yc(z)和Tc(z)都设为z的3次方程，故只需有4个控制高程(z＝z1、z2、z3、z4)处

对应的Yc和Tc便可按插值方法确定Yc(z)和Tc(z)方程中的系数。这样拱冠梁剖面便可用以

下12个参数来描述：

[0064][0065] 上述Yci和Tci分别为控制高程zi处的拱冠梁中心线y坐标及拱冠梁厚度。水平拱圈只要确定了拱圈中心线的方程及拱冠至拱端的厚度变化规律便可确定水平拱圈的形状。

[0066] 本实施例中，以上12个参数的值为：[0067][0068] 步骤5：基于快速多级边界元法自动化计算应力，利用等效应力法求出拱坝上下游面的三个应力分量，然后由该处微分体的平衡条件求得另外三个应力分量，进一步决定拱

坝上下游面的主应力，进而确定出最优拱坝坝形。利用快速多级边界元法自动化计算应力，

具体公式为：

[0069] cij(x)ui(x)＝∫Γ[Uij(x,y)tj(y)?Tij(x,y)uj(y)]dΓ(y)(3)[0070] 其中，x,y为研究域内的点，i,j＝1,2,3，2π·cij(x)为点x在边界内部的角，当点x所在的边界处平滑时，有cij(x)＝0.5。Γ为研究域的边界，ui(x)和ti(y)为位移和面力，以

及Uij(x,y)和Tij(x,y)为开尔文基本解，具体表达式为

[0071][0072] 其中，r为点x和y的距离， 代表r关于xi的偏导数，v为泊松比，μ为剪切模量，以及

[0073][0074] 快速多级算法的具体计算步骤为：[0075] (5.1)对求解域进行计算域边界的离散；(5.2)尽可能采用最小的立方体包含求解域，将其命名为根立方体，也就是第0层立方体，设定每个叶子立方体中最多包含的节点个

数，构造八叉树结构；(5.3)指定多级展开系数，对于每个层数大于1的立方体，定义存储多

极矩的矩阵。对于每个层数大于2的立方体，定义存储局部展开系数的矩阵。对于每一个叶

子立方体，定义矩阵存储系数。定义矩阵存储每个叶子近点直接计算的系数；(5.4)对于每

个叶子立方体分别计算其相关系；(5.5)计算每个叶子立方体的多极矩；(5.6)对于处于第L

层每一个非叶子立方体计算其多极矩，这个多极矩由其子立方体传递而来，这个过程可以

称为向上传递。(5.7)对于第L层的每一个立方体，其中L从2到最终层，通过其相互作用列表

中的立方体传递而来；(5.8)分别计算近点和远点的贡献系数；(5.9)更新迭代向量；(5.10)

不断重读以上过程直到达到迭代算法的控制残差，进而求得应力和位移的未知量。

[0076] 等效应力法为：[0077] 首先，从快速多级边界元计算模块中提取出各单元的形心点处在整体坐标系(x－y－z)下的6个应力分量，之后按转轴公式计算局部坐标系(t－r－z)下的的3个应力分量

σt,σz,τzt，其中σz不变。

[0078][0079] τtz＝τzxcosα+τzysinα(7)[0080] 式中：α为x轴到t轴的转角，以逆时针为正，为拉应力，σx,σy,σz,τxy,τyz,τzx代表直角坐标系—xyz中的应力分量。σt和τtz为空间坐标系—rtz(r轴指向拱轴线的法线，t轴指向

拱轴线的切线，z轴铅直向上。)中的应力分量。

[0081] 若求出拱坝某水平截面上沿坝厚三个单元的形心处的梁应力σz，记为σ1,σ2,σ3(如图1和图2所示)。将x轴设在此3个单元的连线上，已知x＝?1、0、1三点处的正应力σz，上下游

处的等效应力：

[0082][0083] 式中，轴力 和弯矩 为：[0084][0085][0086] 同理，可以求出σz,τzt对应的另两组等效应力。[0087] 利用求得的上下游处的等效应力，再根据坝面微分体的平衡条件(如图3和图4所示)，可以求出该处其它三个等效应力：

[0088][0089][0090][0091] 式中：p为坝面所受的水压力强度，以压为正，θ为坝面与r轴的交角，为不同方向上的等效应力。

[0092] 有了6个等效应力分量之后，根据下式，得3个实数解便为等效应力的3个主应力。[0093] σ3?I1σ2+I2σ?I3＝0(14)[0094] 其中，[0095][0096] 式中，I1，I2和I3为应力不变分量；σx,σy,σz,τxy,τyz,τzx代表应力分量。[0097] 结合以上公式，本实施例可以得出优化坝体的上、下游最大主拉、压应力分别，并保证满足规范对坝体最大拉压应力的要求。

[0098] 应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权

利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发

明的请求保护范围应以所附权利要求为准。