高毅超，梅 真

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

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

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

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

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

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

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

考慮如圖1所示的典型大壩-庫水系統(tǒng)，其中向上游無限延伸的層狀庫水被豎直截?cái)噙吔绶指畛山鼒?chǎng)和遠(yuǎn)場(chǎng)兩個(gè)部分。具有不規(guī)則幾何形狀的近場(chǎng)庫水和壩體構(gòu)成近場(chǎng)域采用有限單元離散；遠(yuǎn)場(chǎng)規(guī)則庫水簡(jiǎn)化成半無限等截面層狀介質(zhì)，遠(yuǎn)場(chǎng)庫水對(duì)近場(chǎng)域的輻射阻尼效應(yīng)通過截?cái)噙吔缟系南嗷プ饔昧?{r}表現(xiàn)。

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

豎直截?cái)噙吔绮捎帽壤吔缬邢迒卧x散，一個(gè)典型的單元如圖2所示，其比例相似中心取在下游無窮遠(yuǎn)處。

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

考慮豎直邊界截?cái)嗟陌霟o限等截面層狀庫水，將庫水視為理想聲學(xué)流體介質(zhì)，并且不考慮水庫庫底的吸收作用，動(dòng)水壓力表示的比例邊界有限單元方程為

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

(1)

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

(2)

(3)

(4)

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

[M0]=[E0]/c2

(5)

在頻域內(nèi)，截?cái)噙吔缟系刃Ч?jié)點(diǎn)荷載{R(ω)}和節(jié)點(diǎn)動(dòng)水壓力{P(ω)}之間的關(guān)系可以用動(dòng)力剛度[S∞(ω)]表示(其中ω表示頻率)，即

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

(6)

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

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

(7)

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

2 雙漸近透射邊界單元

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

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

(8)

(9)

(10)

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

(11)

(12)

(13)

(14)

(15)

(16)

(17)

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

(18)

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

(19)

單元?jiǎng)偠染仃嘯Kh]和阻尼矩陣[Ch]分別為

(20)

(21)

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

(22)

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

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

步驟1 初始化變量；

步驟2 從二進(jìn)制文件中讀取單元?jiǎng)偠染仃嚭妥枘峋仃嚕?/p>

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

3 算 例

3.1 重力壩算例

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

(a)幾何模型

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

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

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

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

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

3.2 拱壩算例

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

計(jì)算中，壩體混凝土和庫水材料參數(shù)、雙漸近透射邊界單元階數(shù)同重力壩算例，并與自由邊界(不考慮截?cái)噙吔绲奈者吔鐥l件)的計(jì)算結(jié)果進(jìn)行對(duì)比。

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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

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