考慮鍵槽咬合作用的高拱壩非線性橫縫模型研究

彭 剛， 洪海豐， 陳燈紅， 楊紫輝， 趙藝園， 曹文昱

(1. 三峽大學 防災減災湖北省重點實驗室，湖北 宜昌 443002；2. 三峽大學 土木與建筑學院，湖北 宜昌 443002；3. 三峽大學 水利與環境學院，湖北 宜昌 443002)

我國在西南地區建設了一批300 m級的大庫容超高拱壩，如錦屏一級(305 m)、小灣(294.5 m)、白鶴灘(289 m)等重大水利水電工程[1]。西南地區地震活動頻繁，且該地區擁有我國約80%的水能資源[2-3]，在該區修建超高拱壩其抗震安全問題日益突出。橫縫工作性態是超高拱壩抗震安全評價的重要組成部分[4]，考慮橫縫動接觸對大壩的抗震安全[5]和抗震設計具有十分重要的意義。

在拱壩建設過程中通常用帶鍵槽的橫縫把壩體分割成柱狀澆筑塊，待壩體冷卻到一定的溫度后，進行接縫灌漿處理，使壩體形成整體，提高超載能力[6]。Niwa等[7]于1982年首次指出拱壩在強震作用下橫縫可能會發生張開、閉合及滑動的現象，這種現象將對大壩的位移和應力分布產生影響。橫縫的變形和受力過程看作是摩擦接觸現象，使得壩體靜動力響應呈現非線性特性。此后，眾多學者根據橫縫的幾何、變形特性進行了更深入的研究。在實際工程分析中，拱壩橫縫的模擬力學模型主要有聯結單元模型[8]、裂縫涂抹式模型和動接觸力模型[9]。陳厚群等[10]對Fenves提出的界面單元模型中進行了改進，對小灣拱壩進行了動力分析。Kuo[11]于1982年提出了裂縫涂抹模型，引進接縫無滑動假定，認為接觸面有足夠的摩擦以阻止界面滑動。牛智勇等[12-13]在B-可微方程組形式的橫縫接觸模型基礎上，針對法向、切向約束條件進行了修正，提出了考慮橫縫初始抗拉強度、鍵槽切向抗剪力與張開度相關的修正接觸模型。Alembagheri等[14]通過非線性增量動力分析，研究了有橫縫拱壩的抗震性能和安全性。田碩等[15-16]采用擴大基礎等措施會降低壩體中部的橫縫開度；采用考慮鍵槽咬合作用、橫縫設置等的非線性動力分析模型，并對白鶴灘拱壩進行地震易損性分析。Liang等[17-18]采用拉格朗日乘子法建立了動態接觸模型以考慮橫縫接觸的非線性。李海楓等[19]采用帶開-閉迭代功能的接觸單元模擬橫縫的力學特性。上述學者基本是采用單一的橫縫力學模型。聯結單元模型可以模擬橫縫的非線性關系，但不能保證接觸體之間不發生相互嵌入，李靜等[20]根據經驗和判斷選取單元的法向和切向剛度系數。動接觸模型把界面影響看作是接觸問題，通過在接觸面上施加約束條件，來確定滿足接觸約束條件和平衡條件的接觸力，常忽略鍵槽的咬合作用。因此，本文將探究動接觸-聯結混合模型模擬考慮鍵槽咬合作用的橫縫接觸非線性。

目前，工程實際常用的鍵槽形狀有矩形鍵槽[21]、梯形鍵槽、弧形鍵槽、球形鍵槽。考慮梯形鍵槽咬合作用的橫縫模型有精細鍵槽模型、鍵槽深度和坡度不變的簡化鍵槽模型、鍵槽坡度不變的簡化鍵槽模型[22-23]、平縫模型(附加切向彈簧剛度)[24]四類。本文在文獻[25]的基礎上，將從動接觸-聯結混合模型的力學模型出發，以精細鍵槽模型在位移荷載作用下參考點的相對位移為基準，對比分析了簡化模型一、簡化模型二和平縫模型的鍵槽咬合作用，并對平縫模型附加切向彈簧剛度系數的取值及其穩定性展開研究，以錦屏一級拱壩為研究對象，建立三維有限元分析模型，考慮鍵槽咬合作用的接觸橫縫模型研究其對拱壩的整體性作用和拱壩橫縫的切向位移及橫縫的法向開度影響。

1 動接觸-聯結混合模型

橫縫接觸狀態的改變在力學上屬于邊界非線性。如圖1所示的動接觸邊界，由主面和從面構成，從面上的結點S′在主面上有且僅有一個確定的結點M′(錨點)與之對應，錨點與該結點的連線方向即為接觸的法線方向；從面結點與錨點之間沿法線和切線方向的相對距離，分別表征動接觸的橫縫開度與切向位移。通過從面結點與主面錨點之間的相對位移可以唯一確定接觸的狀態。

對于任意從面結點S′在主面區間(M-1,M)的對應點M′，參數βS可以通過式(1)來確定。

(1)

從面結點S′到主面區間(M-1,M)的法向和切向距離可由式(2)確定。當h>0時，主從面接觸；當h<0時，主從面分離；當h=0時，主從面接觸且相互之間沒有力的作用。

(2)

由變分法推導出橫縫開度h的變分為

δh=-n(δuM′-δxM-1-ξδv)

(3)

為改善迭代計算的收斂性，并且保證收斂條件滿足精度要求，引入指數型應力-位移模型來描述接觸面的法向接觸本構，其數學表達式為

(4)

(5)

式中：c為初始間隙；h為接觸面之間的相對位移(以嵌入為正)；p為接觸點對上的接觸壓力；p0為特征接觸力。式(5)反映了由法向開度的微小攝動引起的接觸力變化，即法向的剛度貢獻。若選取c為一合理的微量，并保證p0值足夠大，可以將接觸面間的嵌入深度控制為可接受的微量(近似滿足接觸面間法向不可嵌入)，能夠滿足動力計算縫兩側壩體在張開狀態下無相互作用，在縫面接觸時保持接觸面之間的法向特性光滑過渡，從而使得采用隱式迭代算法時的收斂性相對較好。

應用變分原理，得到接觸邊界上法向勢能的變分方程

δ∏=pδh+hδp=[p+(?p/?h)h]δh

(6)

結合式(4)、式(5)，代入式(6)中，并運用鏈式求導法則，可推導出

(7)

將式(6)中的節點位移變分向量的系數集成到總體剛度矩陣中，形成包含接觸邊界的總體剛度矩陣，并由此推導出基于牛頓迭代的有限元隱式求解算法。在平縫模型中，通過附加切向彈簧約束相鄰壩段間的順河向位移，來考慮橫縫鍵槽的咬合作用。采用聯結單元模型中的點對彈簧單元(如圖2所示)，沿縫面布置附加切向彈簧剛度為KS的線性彈簧單元，通過調整KS可達到控制切向位移的目的。本文第3章將對KS的取值進行了討論。

2 梯形鍵槽咬合作用

2.1 梯形鍵槽的傳力特點

梯形鍵槽傳遞剪力主要是通過嚙合來進行的，單向梯形鍵槽(如圖3(a)所示)傳遞剪力時，鍵槽兩邊斜面相互接觸。壩體橫縫灌漿前，接觸面法向存在初始開度d0(如圖3(b)所示)，梯形鍵槽的傾斜角為β，則鍵槽平行于縫面的初始間隙d′可由式(8)求得。

m=d′=d0tanβ

(8)

梯形鍵槽的三種接觸狀態：張開、接觸、閉合。橫縫接觸面切向發生相對錯動時：當|Δu|-(Δwn+d0)/m<0時，鍵槽面之間不接觸，縫面自由(如圖4(a)所示)；當|Δu|-(Δwn+d0)/m≥0時，鍵槽面之間會發生接觸(如圖4(b)所示)；當Δu=0且Δwn+d0=0時，縫面閉合不產生切向咬合力(如圖4(c)所示)。

如圖3(b)所示，橫縫面上的荷載由法向荷載和切向荷載。其中，法向荷載引起的應力利用第1章的指數型應力-位移模型來描述；將切向荷載引起的應力分為兩部分：切向摩擦引起的應力τf和附加切向約束(即鍵槽咬合)引起的切向應力τtt；將τf和τtt在橫縫面上均勻化后用切向應力表達，即

τt=τf+τtt
τf=μσn

(9)

式中：τf為橫縫面切向摩擦力引起的應力，可用庫侖摩擦來表征；σn為橫縫面法向接觸力引起的應力；μ為摩擦因數。

由式(9)可知，當橫縫張開時，σn=0,τf=0，因此，橫縫面切向應力有鍵槽承擔，體現單鍵槽的的傳剪作用的附加切向應力τtt為

(10)

式中：Ktt反應鍵槽咬合作用的附件切向剛度系數；Δu為橫縫面切向相對位移；Δwn為橫縫面法向相對位移；d0為橫縫初始開度；m為鍵槽坡度。當Δu>1，sign(Δu)=1，當Δu<1，sign(Δu)=-1。

2.2 簡化鍵槽模型與精細鍵槽模型咬合作用對比

為探究簡化鍵槽模型與精細鍵槽模型咬合作用，計算了鍵槽深度和坡度不變的簡化模型(以下簡稱：簡化模型一，如圖5(a)所示)、鍵槽坡度不變的簡化模型(以下簡稱：簡化模型二，如圖5(b)所示)、反映鍵槽實際尺寸的精細模型(以下簡稱：精細模型，如圖5(c)所示)的參考點相對位移。

(11)

(12)

鍵槽為單向梯形鍵槽，鍵槽兩邊混凝土塊體的彈性模量均為25.5 GPa，泊松比為0.189。模型的右側面水平和豎直方向為固定邊界；左側為自由邊界，對左側面施加如圖6的豎向位移時程。

計算模型網格尺寸為0.122 m×0.200 m，間距d0=0.005 m，n=2時參考點BB’、CC’、DD’的相對位移如圖7(a)所示。無論采用簡化模型一還是簡化模型二，不考慮附加切向剛度影響，就無法使簡化模型達到與精細模型同等的咬合作用。當n=2時，簡化模型一、簡化模型二、精細模型豎向位移的負向峰值分別為-0.028 8 m，-0.025 0 m，-0.018 0 m，簡化模型一、簡化模型二與精細模型豎向位移的相對誤差分別為60.0%，38.9%。說明簡化模型二的咬合作用優于簡化模型一，但兩種簡化模型的相對誤差較大。參考點BB’，CC’，DD’的法向開度如圖7(b)所示，簡化模型一、簡化模型二、精細模型法向開度的峰值分別為0.016 1 m，0.016 5 m，0.016 8 m，簡化模型一、簡化模型二與精細模型法向開度的相對誤差分別為4.2%，1.8%，對比三個模型法向開度變化很小。

3 平縫模型附加切向彈簧剛度取值研究

通過含有兩個鍵槽的精細模型與同尺寸帶附加彈簧平縫模型的數值試驗來確定附加切向彈簧剛度KS的取值。帶鍵槽橫縫接觸面切向相互作用力由兩部分組成：庫倫摩擦力[27]和附加切向力(鍵槽咬合力)。當橫縫張開時，切向摩擦力為0，只有附加切向力發生作用。因此，橫縫張開情況下，根據兩側的變形來確定平縫模型的附加切向彈簧剛度KS的值。

3.1 附加切向彈簧剛度的取值

為確定合理的剛度系數，建立了精細鍵槽模型(示意圖如圖5(c))和平縫模型(示意圖如圖8(a)，計算模型網格如圖8(b))。計算模型分為左右兩部分，中間設橫縫，橫縫接觸面的摩擦因數取為0.8，橫縫間隙d0=0.005 m，混凝土彈性模量為25.5 GPa，泊松比為0.189。施加的邊界條件和位移加載方式同2.2節算例，豎向位移加載時程如圖6所示。精細鍵槽模型豎直面不發生接觸，左右塊體間的作用力通過鍵槽傾斜面來傳遞。以精細鍵槽模型中橫縫兩側的豎向變形為依據，為平縫模型選取合適的附加切向彈簧剛度，使兩模型的豎向變形趨于一致。

當n=2時，平縫模型和精細模型參考點的相對豎向位移曲線，如圖9所示，沿縫面共布置17個點對彈簧單元，其中附加切向彈簧剛度KS=0.15 GPa/m，單位與文獻[28-29]中橫縫屬性的單位一致。此時參考點A，A′豎向相對位移正負向峰值為0.026 5 m,-0.019 4 m，與精細模型參考點豎向相對位移正負向峰值+0.024 2 m，-0.018 0 m的相對誤差分別為8.6%，7.8%。

相較于簡化模型一、簡化模型二60.0%，38.9%的相對誤差，附加切向彈簧平縫模型的相對誤差明顯減小，即在10%以內。對平縫模型附加切向彈簧剛度的試算中，參考點豎向相對位移正負峰值的相對誤差很難控制在5%以內。因此，本文在這里將正負峰值相對誤差調整至10%。綜上所述，考慮到拱壩有限元分析時建模過程中的便利性，本文采用便于建模的平縫模型來模擬橫縫，采用附加切向彈簧的方式來模擬鍵槽的咬合作用。

3.2 附加切向彈簧剛度取值穩定性的驗證

為探究平縫模型附加切向彈簧剛度與單元特征尺寸(接觸面單元豎向尺寸)的關系，改變圖8(b)所示計算模型網格尺寸，由0.122 m×0.200 m改為0.244 m×0.200 m，0.367 m×0.200 m，橫縫面豎向尺寸分別為0.122 m，0.244 m，0.367 m。通過調整附加切向彈簧剛度值，使得修改單元尺寸模型參考點相對位移與精細模型參考點相對位移的相對誤差在10%以內，得出的附加切向彈簧剛度值分別為0.15 GPa/m(沿橫縫面布置17個彈簧)、0.32 GPa/m(布置8個彈簧)、0.50 GPa/m(布置5個彈簧)。三種網格尺寸模型橫縫面布置附加彈簧總剛度值分別為2.55 GPa/m，2.56 GPa/m，2.50 GPa/m，取平均值2.54 GPa/m作為模擬兩個橫縫鍵槽平縫模型的附加切向彈簧剛度值，則單個鍵槽咬合作用可由切向剛度值為1.27 GPa/m的附加彈簧替代，不同網格尺寸平縫模型與精細模型豎向相對位移對比，如圖10所示。

3.3 三維平縫模型附加切向彈簧剛度系數取值

附加切向彈簧剛度系數是對實際鍵槽切向咬合作用在橫縫接觸面上的均勻化反映，與連接在橫縫接觸面對應節點上的集中彈簧相等效，其變形方向沿鍵槽橫向，垂直于鍵槽法向。集中彈簧的附加彈簧剛度與共用節點單元面積有關。

在考慮橫縫接觸面法向指數模型和切向庫倫摩擦模型的基礎上，鑒于實際應用中橫縫兩側單元網格剖分較一致的特點，利用兩側節點構造8節點單元，如圖11，節點1～4和節點5～8分別分布在橫縫的兩個接觸面上。局部坐標下該單元剛度矩陣中的元素可由式(13)求得

Kij=KttiAiδi(j-4)(i,j=1,…,4)

(13)

式中：Ktti為單位面積鍵槽模型咬合作用的等效附加切向彈簧剛度系數Ktti=1.27 GPa/m；Ai為節點1～4所圍成的面積對節點i的貢獻面積值；δij為克羅內克函數。

4 工程應用

考慮到空庫工況為高拱壩橫縫張開-閉合的控制工況，以錦屏一級高拱壩為研究對象，建立三維有限元分析模型。其中，采用考慮鍵槽咬合作用的動接觸-聯結混合模型模擬地震作用下橫縫的非線性接觸，研究鍵槽咬合作用對于拱壩地震作用下橫縫張開-閉合-滑移的影響。

4.1 工程概況及計算模型

錦屏一級為雙曲拱壩，壩頂高程1 885 m，建基面高程1 580 m，最大壩高305 m。電站正常蓄水位1 880 m，死水位1 800 m，拱冠梁頂厚16 m，拱冠梁底厚63 m，最大中心角93.12°，頂拱中心線弧長552.23 m，厚高比0.207，弧高比1.811。三維有限元分析模型沿壩體高度方向分為17層，沿壩體厚度方向分位5層，共離散6 105個單元，其中壩體單元2 205個，地基單元3 900個，單元均為八節點六面體單元。為突出研究重點，假定該拱壩設置3條橫縫，將壩體分給4個壩段，橫縫間距為136～142 m，橫縫采用平縫模型，通過施加附加切向彈簧模擬鍵槽的咬合作用。有限元模型如圖12所示，分縫位置如圖13所示。坐標系采用笛卡爾坐標系，其中x軸正方向為橫河向左岸指向右岸(上游面)，y軸正方向為順河向，z軸正方向為豎直向。

計算模型材料參數如下：壩體不考慮材料分區，混凝土密度為2 400 kg/m3、靜彈性模量為28.8 GPa、泊松比為0.167；為突出研究重點，地基按經典、實用的無質量地基模型考慮，靜彈性模量為30.0 GPa、泊松比為0.2。動力計算時，拱壩結構阻尼采用瑞利阻尼，阻尼系數α=2ξω1ω2/(ω1+ω2)，β=2ξ/(ω1+ω2)，其中ξ為阻尼比，ω1和ω2為選取的模態頻率上下限。取ξ=5%，結合前兩階自振頻率得到空庫拱壩的阻尼系數α=0.778 8 s-1，β=0.002 9 s。空庫工況是橫縫開度的控制工況，本模型不考慮庫水的影響，計算荷載為壩體自質量和地震荷載。地震荷載采用邊界輸入的方式，輸入地震荷載加速度時程曲線如圖14所示，地震動峰值加速度為0.20g，時間步長為0.02 s，地震持時20 s。

橫縫為單向鍵槽設置，模型計算考慮兩種情況：①不考慮鍵槽咬合作用時，取KS=0；②考慮鍵槽咬合作用時，附加切向彈簧剛度系數取1.27 GPa/m，結合有限元網格每層壩體厚度和節點對橫縫面豎向尺寸的貢獻來取值。橫縫接觸面的非線性相互作用，其中，法向采用指數接觸模型；切向采用庫倫摩擦模型及考慮鍵槽咬合作用的附加剛度彈簧的點對聯結單元模型。通過試算確定指數模型的c=5×10-7m，p0=5 GPa，可較好模擬橫縫法向張開-閉合的力學行為；切向動摩擦因數f取0.6[30]、附加切向彈簧剛度系數為1.27 GPa/m，可較好模擬鍵槽咬合作用。取兩種計算模型左、中、右三條橫縫壩頂部位的法向開度與切向位移進行對比分析。

4.2 鍵槽咬合作用對橫縫開度和切向位移的影響

地震作用下，上述兩種工況下三條橫縫法向開度和切向位移峰值，如表1所示。以中縫為例，上述兩種工況下其上游壩頂參考點的法向開度和切向位移時程，如圖15所示。

由表1可見，在地震動荷載作用下，考慮鍵槽咬合作用能增強分縫拱壩的整體性，減小橫縫的法向開度，且能效限制橫縫的切向位移。考慮鍵槽咬合作用后，左縫的法向開度峰值由2.609 cm減小至2.184 cm，降低幅度約16.3%；右縫的法向開度峰值由3.868 cm減小至3.003 cm，降低幅度約22.4%。如圖15(a)所示，考慮鍵槽咬合作用，即以動接觸-聯結混合模型模擬橫縫的非線性接觸，中縫法向開度峰值降低約20.1%。如圖15(b)、圖15(c)所示，考慮鍵槽咬合作用后，中縫的切向位移大幅降低，中縫的切向位移限制在毫米級，該結論與江守燕等、Jiang等、苑舉衛研究中的結論相一致。同樣，考慮鍵槽咬合作用，對左、右縫有相同的作用效果。綜上，地震荷載作用下，考慮鍵槽的咬合作用對有縫拱壩的橫縫開度和切向位移響應均有較大影響。

表1 地震動作用下橫縫法向開度和切向位移極值

5 結 論

本文根據鍵槽的傳剪特征，建立了帶鍵槽橫縫的小試件有限元精細模型。分析結果表明：簡化模型一和簡化模型二與精細模型相比豎向位移的相對誤差較大，而通過調節平縫模型的附加切向彈簧剛度系數可以實現與精細化模型精度相當的咬合作用。因此，推薦采用帶附加切向彈簧剛度的平縫模型來模擬橫縫鍵槽的咬合作用。

在錦屏一級拱壩的地震響應分析中，研究了鍵槽咬合作用對橫縫的法向開度和切向位移的影響。研究表明，聯結單元模型的引入使得橫縫模型具備了模擬鍵槽咬合作用的能力。同時，聯結單元的引入改善了僅以動接觸力模型模擬動力作用下橫縫的非線性接觸中出現的兩側壩體互相嵌入。采用考慮鍵槽咬合作用的非線性接觸橫縫模型能夠較好地模擬橫縫拱壩的力學性能，并有效地限制橫縫的切向位移，使得橫縫的法向開度小幅降低。