基于ABAQUS的動水壓力波雙漸近透射邊界單元及應用

高毅超，梅 真

(1.華僑大學 土木工程學院，福建 廈門 361021；2.福建省結構工程與防災重點實驗室，福建 廈門 361021)

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

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

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

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

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

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

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

圖1 大壩-庫水耦合系統Fig.1 A coupled dam-reservoir system

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

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

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

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

(1)

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

(2)

(3)

(4)

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

[M0]=[E0]/c2

(5)

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

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

(6)

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

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

(7)

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

2 雙漸近透射邊界單元

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

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

(8)

(9)

(10)

根據Gao等的推導結果，動水壓力波高階雙漸近透射邊界在時域內的表達式為

(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)

雙漸近透射邊界單元的剛度矩陣、阻尼矩陣為常數矩陣，在分析過程中只需要計算一次，因此可以預先求解，再根據需要讀入ABAQUS，提高計算效率。本文在程序實現上，將預先計算得到的雙漸近透射邊界單元的剛度矩陣和阻尼矩陣分別保存到兩個二進制文件中，在ABAQUS調用UEL的過程中根據分析需要讀取。

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

步驟1 初始化變量；

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

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

3 算 例

3.1 重力壩算例

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

(a)幾何模型

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

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

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

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

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

3.2 拱壩算例

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

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

圖6 拱壩-庫水系統Fig.6 An arch dam-reservoir system

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

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

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

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

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

4 結 論

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

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

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

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

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