高毅超，梅 真

(1.華僑大學(xué) 土木工程學(xué)院，福建 廈門 361021；2.福建省結(jié)構(gòu)工程與防災(zāi)重點(diǎn)實(shí)驗(yàn)室，福建 廈門 361021)

在強(qiáng)震作用下，大壩與庫水之間發(fā)生的顯著動力相互作用將給大壩的動力響應(yīng)帶來重大影響。其中，如何準(zhǔn)確模擬動水壓力波在半無限庫水中傳播引起的輻射阻尼效應(yīng)是該研究領(lǐng)域的熱點(diǎn)和難點(diǎn)。在大壩-庫水動力相互作用分析中，有限元法是有效的數(shù)值方法，但是需要在計算域的截斷邊界施加人工邊界以避免波動的虛假反射，或者與其他無限域數(shù)值方法耦合以考慮無限介質(zhì)輻射阻尼的影響。Kü?ükarslan等[1]結(jié)合有限元法和Sommerfeld輻射邊界提出了結(jié)構(gòu)-流體動力相互作用的數(shù)值分析方法，但是Sommerfeld輻射邊界在理論上僅僅具有零階精度，需要增加近場域離散范圍以獲得良好的計算精度。邊界元法[2-3]只需離散邊界使得問題維數(shù)降低一維，并且能夠自動滿足無窮遠(yuǎn)處的輻射邊界條件，因而也被應(yīng)用于結(jié)構(gòu)-流體動力相互作用分析。然而邊界元求解所需的基本解可能十分復(fù)雜或者不存在，大大限制了邊界元在實(shí)際工程中的應(yīng)用。

比例邊界有限單元法(Scaled Boundary Finite Element Method,SBFEM)是Song等[4]提出的一種半解析數(shù)值方法，它只需離散邊界而在徑向保持解析，兼具有限元法和邊界元法的優(yōu)點(diǎn)。它不需要基本解即可自動滿足無窮遠(yuǎn)處的輻射邊界條件，因此非常適用于無限域動力問題的模擬。杜建國等[5]基于SBFEM提出了一種求解壩面動水壓力的半解析方法；Lin等[6-8]將SBFEM應(yīng)用于大壩-庫水動力相互作用分析，然而其時域的求解均是基于耗時的卷積積分運(yùn)算完成的。

高階局部透射邊界[9]具有高階精度和較高的計算效率，然而絕大部分高階透射邊界屬于高頻單向漸近邊界[10]，無法同時模擬層狀介質(zhì)中行波和快衰波的傳播。Prempramote等[11]基于動力剛度連分式雙漸近解提出了一類高階雙漸近透射邊界，隨著漸近階數(shù)的增加，該透射邊界在全頻范圍內(nèi)迅速收斂到準(zhǔn)確解，具有很高的計算精度和計算效率。Wang等[12-14]結(jié)合SBFEM將該透射邊界推廣應(yīng)用到層狀庫水的模擬，構(gòu)建了動水壓力波高階雙漸近透射邊界，并建立了有限元-動水壓力波高階雙漸近透射邊界的時域耦合分析模型。數(shù)值算例分析結(jié)果表明，該耦合分析模型具有很高的計算精度和較高的計算效率。

ABAQUS是一套功能強(qiáng)大的通用有限元分析軟件，它包括了豐富的單元庫和材料模型庫，并且具有良好的前后處理程序和強(qiáng)大的非線性方程求解器，在許多工程領(lǐng)域都得到了廣泛應(yīng)用。此外，ABAQUS提供了豐富的用戶子程序接口，用戶可以通過子程序接口自定義單元、材料、荷載和邊界條件等，具有良好的可擴(kuò)展性。

本文基于ABAQUS提供的用戶子程序UEL接口，編寫了動水壓力波高階雙漸近透射邊界單元，以超單元的形式直接將雙漸近透射邊界嵌入到近場有限元方程，實(shí)現(xiàn)了有限元-雙漸近透射邊界的時域直接耦合分析模型。通過數(shù)值算例驗(yàn)證雙漸近透射邊界單元程序編寫的正確性，并將其應(yīng)用到大壩-庫水動力相互作用分析。

1 庫水的比例邊界有限單元方程

考慮如圖1所示的典型大壩-庫水系統(tǒng)，其中向上游無限延伸的層狀庫水被豎直截斷邊界分割成近場和遠(yuǎn)場兩個部分。具有不規(guī)則幾何形狀的近場庫水和壩體構(gòu)成近場域采用有限單元離散；遠(yuǎn)場規(guī)則庫水簡化成半無限等截面層狀介質(zhì)，遠(yuǎn)場庫水對近場域的輻射阻尼效應(yīng)通過截斷邊界上的相互作用力-{r}表現(xiàn)。

圖1 大壩-庫水耦合系統(tǒng)Fig.1 A coupled dam-reservoir system

豎直截斷邊界采用比例邊界有限單元離散，一個典型的單元如圖2所示，其比例相似中心取在下游無窮遠(yuǎn)處。

圖2 典型的比例邊界有限單元Fig.2 A typical scaled boundary finite element

考慮豎直邊界截斷的半無限等截面層狀庫水，將庫水視為理想聲學(xué)流體介質(zhì)，并且不考慮水庫庫底的吸收作用，動水壓力表示的比例邊界有限單元方程為

[E0]{p},ξξ-[E2]{p}-[M0]{p},tt=0

(1)

式中:{p}為動水壓力;{p},ξξ為動水壓力對比例邊界坐標(biāo)ξ的二次偏導(dǎo);{p},tt為動水壓力對時間的二次導(dǎo)數(shù)；系數(shù)矩陣[E0],[E2]和[M0]的定義為

(2)

(3)

(4)

式中:ρ和K分別為庫水的密度和體積模量；[N]為單元插值函數(shù)，[B1]，[B2]和|J|只與邊界的幾何信息有關(guān)，其定義詳見Gao等的研究。對于等截面層狀庫水，[M0]和[E0]存在以下關(guān)系

[M0]=[E0]/c2

(5)

在頻域內(nèi)，截斷邊界上等效節(jié)點(diǎn)荷載{R(ω)}和節(jié)點(diǎn)動水壓力{P(ω)}之間的關(guān)系可以用動力剛度[S∞(ω)]表示(其中ω表示頻率)，即

{R(ω)}=[S∞(ω)]{P(ω)}

(6)

在豎直邊界上應(yīng)用虛功原理，并考慮動力剛度的定義式(6)，可以將式(1)改寫成動力剛度表示的比例邊界有限單元方程

[S∞(ω)][E0]-1[S∞(ω)]-[E2]+ω2[M0]=0

(7)

基于式(7)，結(jié)合動力剛度的連分式漸近解和輔助變量技術(shù)，即可構(gòu)建動水壓力波高階雙漸近透射邊界。

2 雙漸近透射邊界單元

對于等截面層狀庫水，動力剛度表示的比例邊界有限單元方程可以通過模態(tài)變換的方式解耦，考慮如下廣義特征值分解

[E2][Φ]=[E0][Φ][Λ2]/h2
[Φ]T[E0][Φ]=[I]
[Φ]T[E2][Φ]=[Λ2]/h2

(8)

(9)

(10)

根據(jù)Gao等的推導(dǎo)結(jié)果，動水壓力波高階雙漸近透射邊界在時域內(nèi)的表達(dá)式為

(11)

(12)

(13)

(14)

(15)

(16)

(17)

式(17)即為高階雙漸近透射邊界單元的控制方程，其中{ph}和{fh}分別為高階雙漸近透射邊界單元的動水壓力向量及荷載向量，定義分別為

(18)

{fh}={{r}T;{0}T;…;{0}T}T

(19)

單元剛度矩陣[Kh]和阻尼矩陣[Ch]分別為

(20)

(21)

式中:0為零矩陣，黑體字母表示分塊矩陣，其定義分別為

(22)

雙漸近透射邊界單元的剛度矩陣、阻尼矩陣為常數(shù)矩陣，在分析過程中只需要計算一次，因此可以預(yù)先求解，再根據(jù)需要讀入ABAQUS，提高計算效率。本文在程序?qū)崿F(xiàn)上，將預(yù)先計算得到的雙漸近透射邊界單元的剛度矩陣和阻尼矩陣分別保存到兩個二進(jìn)制文件中，在ABAQUS調(diào)用UEL的過程中根據(jù)分析需要讀取。

ABAQUS提供的用戶子程序接口UEL采用FORTRAN語言編寫，并在每一迭代步中被調(diào)用。用戶在UEL中需要定義的變量主要包括AMATRX(剛度矩陣、質(zhì)量矩陣等)和RHS(右端荷載向量)。本文在UEL的框架下，編寫了雙漸近透射邊界單元子程序UEL_DAOB.for，其程序流程如下：

步驟1 初始化變量；

步驟2 從二進(jìn)制文件中讀取單元剛度矩陣和阻尼矩陣；

步驟3 根據(jù)標(biāo)志變量IFLAGS(1)的取值，組裝矩陣AMATRX和右端荷載向量RHS，并返回ABAQUS計算主程序。

3 算 例

3.1 重力壩算例

考慮一個典型的重力壩-庫水系統(tǒng)，重力壩壩高72 m，近場庫離散范圍取18 m(約為1/4庫水深度)，其幾何模型參數(shù)和有限元網(wǎng)格如圖3所示。重力壩和近場庫水均采用4節(jié)點(diǎn)四邊形單元離散，其中壩體單元84個，庫水聲學(xué)流體單元33個；流-固耦合界面采用11個兩節(jié)點(diǎn)界面單元離散；為了同近場庫水網(wǎng)格協(xié)調(diào)，豎直截斷邊界采用11個兩結(jié)點(diǎn)線單元離散。材料參數(shù)取值如下：壩體混凝土彈性模量Ed=30 GPa，泊松比vd=0.2，質(zhì)量密度為2 400 kg/m3；庫水密度為1 000 kg/m3，動水壓力波波速c=1 438.7 m/s。庫水表面的壓力邊界條件為動水壓力p=0，庫底為剛性庫底。

(a)幾何模型

(b)有限元網(wǎng)格圖3 重力壩-庫水系統(tǒng)Fig.3 A gravity dam-reservoir system

為了驗(yàn)證雙漸近透射邊界單元的正確性及計算精度，采用有限元擴(kuò)展網(wǎng)格解作為參考解，并與ABAQUS自帶的黏性邊界模型進(jìn)行對比。在擴(kuò)展網(wǎng)格模型中，遠(yuǎn)場庫水有限單元離散長度為3 600 m，動水壓力波在該范圍內(nèi)往返傳播所需的時間約為5 s(對應(yīng)的無量綱時間為tc/h=100)。在重力壩基底輸入如圖4所示的水平向三角脈沖荷載，計算無量綱時間步長取0.1；時步積分算法采用Newmark-β法，積分常數(shù)取γ=0.5，β=0.25，即平均常加速度法。高階雙漸近透射邊界的漸近階數(shù)取MH=ML=4。

圖4 三角脈沖荷載時程曲線Fig.4 Time history of triangular impulse

在水平向三角脈沖荷載作用下，重力壩壩底動水壓力時程如圖5所示。從圖5中可以看出，高階雙漸近透射邊界單元(MH=ML=4)的計算結(jié)果與有限元擴(kuò)展網(wǎng)格解幾乎完全吻合，而ABAQUS自帶黏性邊界的計算結(jié)果與擴(kuò)展網(wǎng)格解偏差較大，并且很快衰減到零。因此，雙漸近透射邊界單元在程序編寫上是正確的，并且具有很高的計算精度，在實(shí)際應(yīng)用中可以減小近場庫水的離散范圍。

圖5 三角脈沖荷載作用下重力壩壩底動水壓力時程Fig.5 Hydrodynamic pressure responses at heel of the gravity dam subject to triangular impulse

3.2 拱壩算例

考慮圖6所示的拱壩-庫水系統(tǒng)，混凝土雙曲拱壩壩高265 m，水庫蓄水深度為252 m，近場庫水離散范圍約為0.5倍壩高。拱壩和近場庫水統(tǒng)一采用8節(jié)點(diǎn)六面體單元離散，其中壩體單元6 400個，庫水聲學(xué)流體單元9 792個，流-固耦合界面采用816個4節(jié)點(diǎn)四邊形單元離散，單元數(shù)量為816個；豎直截斷邊界采用816個4節(jié)點(diǎn)四邊形單元離散。

計算中，壩體混凝土和庫水材料參數(shù)、雙漸近透射邊界單元階數(shù)同重力壩算例，并與自由邊界(不考慮截斷邊界的吸收邊界條件)的計算結(jié)果進(jìn)行對比。

圖6 拱壩-庫水系統(tǒng)Fig.6 An arch dam-reservoir system

考慮拱壩基底受到順河向(Y向)El Centro地震波作用，其中El Centro地震波的加速度時程曲線參見Wang等的研究。地震波作用下，A點(diǎn)(見圖6)的動水壓力時程曲線如圖7所示。從圖中看出，與不考慮輻射阻尼效應(yīng)的自由邊界計算結(jié)果相比，考慮了雙漸近透射邊界后，拱壩壩踵的動水壓力響應(yīng)大大減小了；其中，采用自由邊界和雙漸近透射邊界模型計算，動水壓力的峰值分別為3.2 MPa和0.8 MPa。

圖7 El Centro地震波作用下A點(diǎn)動水壓力時程Fig.7 Hydrodynamic pressure responses at point A subject to El Centro earthquake

壩踵動水壓力達(dá)到峰值時，拱壩上游面的動水壓力等值線圖如圖8所示。從中可以看出，動水壓力數(shù)值隨著庫水深度的增加而增大；自由邊界模型高估了壩面的動水壓力，從而導(dǎo)致拱壩壩體應(yīng)力的增加，增大了壩體的動力響應(yīng)。

圖8 拱壩上游面最大動水壓力等值線圖Fig.8 Contour map of the maximum hydrodynamic pressure distribution on the upstream surface of the arch dam

由上述結(jié)果可知，采用自由邊界模型，即不考慮無限庫水的輻射阻尼，在一定程度上夸大了壩面動水壓力的數(shù)值，由此夸大了壩體的動力響應(yīng)。無限庫水的輻射阻尼效應(yīng)對壩體動力響應(yīng)有重要影響，是大壩地震響應(yīng)分析中應(yīng)當(dāng)考慮的重要部分。

4 結(jié) 論

高階雙漸近透射邊界具有良好的計算精度和計算效率。本文基于ABAQUS提供的UEL接口，采用FORTRAN語言開發(fā)了雙漸近透射邊界單元UEL_DAOB.for，將高階雙漸近透射邊界嵌入到ABAQUS，通過算例分析得到以下結(jié)論：

(1)重力壩算例驗(yàn)證了雙漸近透射邊界單元程序的正確性，并且表明高階雙漸近透射邊界具有很高的計算精度。

(2)拱壩算例結(jié)果表明，與不考慮無限庫水輻射阻尼的自由邊界計算模型相比，考慮高階雙漸近透射邊界后，壩面動水壓力大大減小了，其峰值從3.2 MPa減低到了0.8 MPa。

(3)不考慮無限庫水輻射阻尼，在一定程度上夸大了壩面動水壓力和壩體動力響應(yīng)，因此在進(jìn)行大壩動力響應(yīng)分析時，考慮無限庫水的輻射阻尼效應(yīng)是十分必要的。

本文編寫的雙漸近透射邊界單元用戶子程序具有良好的計算精度和適用性，可用于實(shí)際大壩的地震響應(yīng)分析。