權利要求書： 1.一種混凝土拱壩壩體優(yōu)化方法，其特征在于：包括以下步驟：步驟1：選擇用以確定拱壩體型的基本參數以及確定壩體厚度的基本參數，作為優(yōu)化設計變量；所述拱壩體型的基本參數包括拱冠梁處拱圈中心線的曲率半徑和拱端中心角；所述壩體厚度的基本參數包括拱冠梁厚度、拱端厚度；

步驟2：選擇適當的目標函數，以安全性為關鍵性指標，即以應力為目標函數；

步驟3：為保證工程安全和施工方便，確定壩體幾何尺寸和工作性態(tài)兩方面的約束，包括幾何約束條件、輸入應力約束條件和穩(wěn)定約束條件；

步驟4：將拱壩的體型抽象成幾何模型并用一組參數來完整地描述，基于快速多級邊界單元進行壩體網格自動化劃分以及荷載的自動施加；

步驟5：基于快速多級邊界元法自動化計算應力，利用等效應力法求出拱壩上下游面的三個應力分量，然后由該處微分體的平衡條件求得另外三個應力分量，進一步決定拱壩上下游面的主應力，進而確定出最優(yōu)拱壩壩形。

2.根據權利要求1所述的混凝土拱壩壩體優(yōu)化方法，其特征在于：所述步驟3中，輸入應力約束條件包括壩軸線位置移動范圍的限制、壩頂最小厚度的限制、倒懸度約束和保凸性要求。

3.根據權利要求1或2所述的混凝土拱壩壩體優(yōu)化方法，其特征在于：所述步驟3中，輸入幾何約束條件為壩體內部的最大主拉應力及最大主壓應力均小于壩體混凝土的允許強度，即

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

式中：[σt]為容許拉應力，[σc]為容許壓應力。

4.根據權利要求1或2所述的混凝土拱壩壩體優(yōu)化方法，其特征在于：所述步驟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坐標及拱冠梁厚度；水平拱圈只要確定了拱圈中心線的方程及拱冠至拱端的厚度變化規(guī)律便能確定水平拱圈的形狀。

5.根據權利要求3所述的混凝土拱壩壩體優(yōu)化方法，其特征在于：所述步驟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坐標及拱冠梁厚度；水平拱圈只要確定了拱圈中心線的方程及拱冠至拱端的厚度變化規(guī)律便能確定水平拱圈的形狀。

6.根據權利要求1或2或5所述的混凝土拱壩壩體優(yōu)化方法，其特征在于：所述步驟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所述的混凝土拱壩壩體優(yōu)化方法，其特征在于：所述步驟5中，利用的快速多級加速邊界元法，快速多級算法的具體計算步驟為：(5.1)對求解域進行計算域邊界的離散；

(5.2)采用最小的立方體包含求解域，將其命名為根立方體，也就是第0層立方體，設定每個葉子立方體中最多包含的節(jié)點個數，構造八叉樹結構；

(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所述的混凝土拱壩壩體優(yōu)化方法，其特征在于：所述步驟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所述的混凝土拱壩壩體優(yōu)化方法，其特征在于：所述步驟5中，根據下式，得3個實數解便為等效應力的3個主應力：

3 2

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

其中，

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

說明書： 混凝土拱壩壩體優(yōu)化方法技術領域[0001] 本發(fā)明屬于壩體優(yōu)化技術領域，具體涉及一種混凝土拱壩壩體優(yōu)化方法。背景技術[0002] 所謂拱壩體型優(yōu)化設計，就是在已知的給定參數條件下，尋找滿足工程的施工、運用、強度、穩(wěn)定性等約束條件，并使目標函數取最小值的體型方案。因此，優(yōu)化設計的首要一

步是將工程實際問題表示成數學表達式，即建立數學模型，包括選定設計變量、選擇目標函

數、建立約束方程，然后采用適宜的優(yōu)化方法，求出問題的解。

[0003] 關于拱壩體型優(yōu)化的研究，國外是從20世紀60年代開始的。1969年，R.Sharp提出了用數學規(guī)劃法對拱壩進行體型優(yōu)化設計，從而使拱壩體型優(yōu)化設計作為一個專門課題，

引起學術界和工程界的注意；1973年，加拿大學者W.Stensch提出了自由型拱壩設計的模

型；1974年G.A.Harley和G.M.Mcneice形成了有限元分析用的網格自由劃分系統(tǒng)，進一步改

進了W.Stcnsch的工作；1975年R.E.Rhichetts和D.C.Liewicz發(fā)表了《混凝土壩的形狀優(yōu)化

設計》，其中也研究了拱壩；1982年，在上海舉行的國際有限元會議上，K.Wassermann發(fā)表了

一篇拱壩體型優(yōu)化的論文等。

[0004] 我國拱壩體型優(yōu)化設計的研究是在20世紀70年代末期由朱伯芳院士提出來的，盡管起步較晚，但發(fā)展很快。我國拱壩優(yōu)化設計研究已經進行了30多年，在優(yōu)化理論和實踐上

均取得了顯著的成績，目前我國在拱壩優(yōu)化領域已處于國際領先地位。從基于材料力學優(yōu)

化發(fā)展到有限元優(yōu)化，從靜力優(yōu)化發(fā)展到動力優(yōu)化，從線性優(yōu)化發(fā)展到非線性全過程分析，

拱圈的線性從單心圓發(fā)展到二次曲線、混合線型等，逐步發(fā)展和完善了我國提出的高混凝

土壩優(yōu)化理論與方法。

[0005] 目前，有限元法作為一種很受歡迎的數值分析方法，逐步成為壩型優(yōu)化應力計算的主要方法，但一般而言，拱壩為超大規(guī)模水利工程，有限元分析需要劃分大量的三維體網

格，計算量巨大，耗時長，需要高額的計算成本。

發(fā)明內容[0006] 針對現(xiàn)有技術存在的不足，本發(fā)明提出了一種混凝土拱壩壩體優(yōu)化方法，基于快速多級邊界元法，模型重建過程更為簡單，同等自由度下，精度更高。

[0007] 為了解決上述技術問題，本發(fā)明提供的混凝土拱壩壩體優(yōu)化方法，其特征在于：包括以下步驟：

[0008] 步驟1：選擇用以確定拱壩體型的基本參數，如拱冠梁處拱圈中心線的曲率半徑、拱端中心角，以及確定壩體厚度的基本參數如拱冠梁厚度、拱端厚度等作為優(yōu)化設計變量；

[0009] 步驟2：選擇適當的目標函數，一般以安全性為關鍵性指標，即以應力為目標函數；[0010] 步驟3：為保證工程安全和施工方便，確定壩體幾何尺寸和工作性態(tài)兩方面的約束，主要包括幾何約束條件、應力約束條件和穩(wěn)定約束條件；

[0011] 步驟4：將拱壩的體型抽象成幾何模型并用一組參數來完整地描述，基于快速多級邊界單元進行壩體網格自動化劃分以及荷載的自動施加。

[0012] 步驟5：基于快速多級邊界元法自動化計算應力，利用等效應力法求出拱壩上下游面的三個應力分量，然后由該處微分體的平衡條件求得另外三個應力分量，進一步決定拱

壩上下游面的主應力，進而確定出最優(yōu)拱壩壩形。

[0013] 作為優(yōu)選，所述步驟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坐標及拱冠梁厚度。水平拱圈只要確定了拱圈中心線的方程及拱冠至拱端的厚度變化規(guī)律便可確定水平拱圈的形狀。

[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層立方體，設定每個葉子立方體中最多包含的節(jié)點個

數，構造八叉樹結構；(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] 本發(fā)明的優(yōu)點及有益效果如下：[0049] 本發(fā)明采用的快速多級邊界元法，一方面基于邊界元法，只需要對壩體表面進行二維網格劃分，而有限元法卻需要建三維網格，建模過程工作量更大，因此，只需要建二維

網格的邊界元法模型重建過程更為簡單快捷，計算量更??；另一方面將快速多極算法和邊

界元法相耦合，進一步提高了計算效率。總而言之，本發(fā)明建立的以應力(安全性)為目標函

數的壩形優(yōu)化求解模型，相較于當前的有限元壩體優(yōu)化方法，在同等自由度的情況下，本發(fā)

明快速多級邊界元法計算效率更快且精度更高，具有不可忽視的優(yōu)勢。

附圖說明[0050] 圖1是本發(fā)明等效應力位置俯視圖；[0051] 圖2是本發(fā)明等效應力計算示意圖[0052] 圖3是本發(fā)明壩面微分體的上游平衡視圖[0053] 圖4是本發(fā)明壩面微分體的下游平衡視圖具體實施方式[0054] 為了便于本領域技術人員理解和實施本發(fā)明，下面將結合附圖及實施例對本發(fā)明作進一步的詳細描述，此外，此處所描述的實施示例僅用于說明和解釋本發(fā)明，并不用于限

定本發(fā)明。

[0055] 實施例1：【以國內某碾壓混凝土雙曲拱壩為例】[0056] 采用本發(fā)明的混凝土拱壩壩體優(yōu)化方法，具體步驟如下：[0057] 步驟1：選擇用以確定拱壩體型的基本參數，如拱冠梁處拱圈中心線的曲率半徑、拱端中心角，以及確定壩體厚度的基本參數如拱冠梁厚度、拱端厚度等作為優(yōu)化設計變量。

在實施例中，選取拱冠梁壩頂、壩底厚度，上、下游面上部倒懸度，壩高，壩體方量，控制高程

最大中心角為基本參數；

[0058] 步驟2：選擇適當的目標函數，本實施例以安全性為關鍵性指標，即以應力為目標函數；

[0059] 步驟3：為保證工程安全和施工方便，確定壩體幾何尺寸和工作性態(tài)兩方面的約束，主要包括幾何約束條件、應力約束條件和穩(wěn)定約束條件，應力約束條件為壩頂最小厚度

的限制，倒懸度約束，保凸性要求。具體來講，壩頂最小厚度的限制根據防護、交通、施工、抗

屈折等方面的具體要求而定，在本實施例中壩頂最小厚度為9m；對于混凝土拱壩，其倒懸度

一般控制在0.3以內；由微分幾何知，可微曲面上某一點臨近的凸凹性，可由高斯曲率K來判

定，壩體曲面的保凸性即為要求曲面上每一點均有K>0。幾何約束條件為壩體內部的最大主

拉應力|σt|max及最大主壓應力|σc|max均小于壩體混凝土的允許強度，即

[0060] |σt|max≤[σt]，|σc|max≤[σc](1)[0061] 式中：[σt]為容許拉應力，[σc]為容許壓應力。[0062] 在本實施例中，穩(wěn)定約束條件為要求拱壩壩體在正常蓄水位+壩體自重+泥沙壓力+正常溫降的組合條件下保持穩(wěn)定，具體來說，保持穩(wěn)定的具體要求為大壩發(fā)生的滑移距離

在設計要求容許范圍之內，并且不發(fā)生傾倒屈折。應力約束條件為上、下游面上部倒懸度，

拱冠梁壩頂、壩底厚度和控制高程最大中心角的最大容許值分別為：0.3，0.3，9.0，15.0米，

21.5米，110度。此外，根據《混凝土拱壩設計規(guī)范》，對于壓應力，規(guī)范要求拱壩的設計最大

壓應力不得大于容許壓應力[σc]，在本實施例中，[σc]＝5.0MPa；對于壓應力，規(guī)范要求壩體

主拉應力不得大于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坐標及拱冠梁厚度。水平拱圈只要確定了拱圈中心線的方程及拱冠至拱端的厚度變化規(guī)律便可確定水平拱圈的形狀。

[0066] 本實施例中，以上12個參數的值為：[0067][0068] 步驟5：基于快速多級邊界元法自動化計算應力，利用等效應力法求出拱壩上下游面的三個應力分量，然后由該處微分體的平衡條件求得另外三個應力分量，進一步決定拱

壩上下游面的主應力，進而確定出最優(yōu)拱壩壩形。利用快速多級邊界元法自動化計算應力，

具體公式為：

[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層立方體，設定每個葉子立方體中最多包含的節(jié)點個

數，構造八叉樹結構；(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] 結合以上公式，本實施例可以得出優(yōu)化壩體的上、下游最大主拉、壓應力分別，并保證滿足規(guī)范對壩體最大拉壓應力的要求。

[0098] 應當理解的是，上述針對較佳實施例的描述較為詳細，并不能因此而認為是對本發(fā)明專利保護范圍的限制，本領域的普通技術人員在本發(fā)明的啟示下，在不脫離本發(fā)明權

利要求所保護的范圍情況下，還可以做出替換或變形，均落入本發(fā)明的保護范圍之內，本發(fā)

明的請求保護范圍應以所附權利要求為準。