DEEPSOIL和FLAC3D一維土層地震反應分析對比研究

梁騰龍，駱星宇，周東*，覃釗

(1.廣西大學土木建筑工程學院， 廣西南寧530004；2.廣西防災減災與工程安全重點實驗室， 廣西南寧530004)

0 引言

土層地震反應分析是工程結構抗震設計的重要環(huán)節(jié)。在場地地震安全性評價中，場地土層地震反應分析軟件將基巖地震動時程和場地地質條件勘測轉化為場地設計反應譜，為工程建筑物提供合理的抗震設計依據(jù)。地震反應分析的計算方法可以分為頻域分析和時域分析2類。工程抗震中常用的等效線性化方法就是一種頻域分析方法。等效線性化方法的基本思想就是在總體動力學效應大致相當?shù)囊饬x上，用一個等效剪切模量及等效阻尼比替代所有不同應變幅值下的剪切模量和阻尼比，將非線性問題轉化為線性問題來求解。土層地震反應分析的另一種計算方法為時域分析方法。時域分析方法是將場地簡化為剪切運動的多自由度集中質量體系，每個質量單元間通過彈簧和阻尼器連接，通過建立剛度、阻尼和質量矩陣求解各單元的動力反應。目前，基于頻域和時域計算原理開發(fā)的土層地震反應分析軟件有SHAKE、LSSRLI-1、SOILQUAKE、DEEPSOIL等[1]。這些軟件的開發(fā)和應用基本解決了一維土層非線性地震反應的計算問題，在硬土場地上的計算結果可靠且與實測地震記錄相接近。

在對復雜工程場地進行地震反應分析時，借助于數(shù)值模擬軟件的計算方法得到了廣泛地使用。陳健云等[2]通過FLAC3D對淺埋軟土場地地鐵車站地震荷載作用下的動力響應進行數(shù)值研究。邊金等[3]借助FLAC3D對地鐵車站進行了動力分析。得益于自由場邊界和多種動力荷載輸入方式，F(xiàn)LAC3D在一維黏彈性土層場地地震反應的模擬上可以得到合理準確的結果，與傳統(tǒng)地震反應分析軟件具有良好的一致性。Mejia等[4]利用SHAKE程序與FLAC3D程序進行了半無限彈性土層模型的地震反應計算，結果表明FLAC3D的計算結果與目標曲線非常接近。蘇建鋒[5]驗證了用FLAC3D軟件方法進行土層地震反應計算的可行性，與等效線性化程序RSLEIBM計算結果的對比表明，在輸入地震動強度較小時，2種方法計算得到的結果差別不大。FLAC3D軟件提供了瑞利阻尼、局部阻尼和滯后阻尼3種阻尼形式來模擬土體的滯回特征，但與前2種阻尼相比，目前利用滯回阻尼的完全非線性分析在實際工程的應用還較少，與其他一維土層地震反應分析軟件計算結果的對比特別是與實測地震記錄的對比還亟待研究。

為了分析FLAC3D在一維土層地震反應分析中與現(xiàn)行主流地震計算程序之間的異同點，理解和掌握FLAC3D在非線性分析中阻尼參數(shù)的選擇，依據(jù)日本國家地球科學與防災研究所建立的KiK-net強震動觀測網(wǎng)提供的硬土場地地震記錄和地質數(shù)據(jù)，對DEEPSOIL頻域、時域分析和 FLAC3D動力分析進行對比研究。通過對比DEEPSOIL和FLAC3D在不同類型場地和不同地震荷載輸入條件下計算的地表地震動峰值和反應譜，檢驗FLAC3D在硬土場地上進行一維土層地震反應分析的適用性，為復雜場地的地震模擬提供依據(jù)。

1 土層地震分析軟件簡介

1.1 DEEPSOIL簡介

DEEPSOIL軟件是由美國伊利諾伊大學香檳分校( UIUC)專家研發(fā)的一維土層地震響應計算程序，可以在頻域(線性和等效線性)和時域(線性和非線性)進行分析[6]。DEEPSOIL有一個圖形用戶界面，可以通過選擇適當?shù)淖枘峁胶头蔷€性土層參數(shù)，借助自動化的擬合曲線程序計算土體的非線性行為。DEEPSOIL中包含線性、等效線性和時域非線性3種分析方法，其中時域非線性分析用Newmark β方法(隱式)或Heun方法(顯式)求解時域運動方程。本文所使用的DEEPSOIL版本是DEEPSOIL v7.0.26 Released版。

1.2 FLAC3D簡介

FLAC3D是Itasca公司開發(fā)的有限差分數(shù)值模擬軟件，通過動力模塊實現(xiàn)三維土層模型的地震反應模擬，在巖土領域中得到了廣泛運用。FLAC的動力計算采用完全非線性方法，除了常用的瑞利阻尼和局部阻尼外，還可以采用滯后阻尼的形式嵌入到程序當中的來模擬土體的滯回現(xiàn)象。由于滯后阻尼中的控制參數(shù)是通過對動剪切模量比隨剪應變的關系曲線進行擬合得到的，因此在應用中可以很好地與等效線性方法形成對應關系。

2 計算土層模型及輸入

2.1 簡化土層剖面及輸入

圖1 簡化土層場地剖面Fig.1 Physical properties and shear wave velocity of subsoil profiles

為對比DEEPSOIL頻域、時域分析與FLAC3D動力分析的異同點，按非線性程度設計了3個簡化土層剖面，以內部運動和露頭運動2種方式輸入地震波。所設計的簡化土層剖面如圖 1所示，3個土層剖面分別為非線性程度依次遞增的等波速土層剖面、變波速土層剖面和非線性土層剖面。

等波速土層剖面的土層厚度H設為30 m，土層剪切波速設為200 m/s，土層重度γ設為18 kN/m3,基巖層剪切波速設為500 m/s，基巖重度γ設為22 kN/m3,阻尼比λ設為5%，土層卓越頻率為1.667 Hz。

變波速土層剖面采用與等波速剖面相同的土層厚度和土層密度。由于實際工程中土層剪切模量通常與有效平均應力呈正相關，因此變波速土層剖面中設定土層剪切波速隨埋深的增加而增大。參考文獻[7]提出的描述剪切波速與深度之間擬合關系公式，變波速土層剖面的剪切波速Us計算公式為

圖2 剪切模量與阻尼比擬合曲線Fig.2 Shear modulus and damping ratio curves

vs=a+bz，

(1)

式中：a、b為擬合參數(shù)；z為土層深度。為使變波速土層剖面與等波速土層剖面具有相同的卓越頻率，擬合參數(shù)分別取值a=135，b=4.71。

非線性土層剖面采用與變波速土層剖面相同的剪切波速設定，通過動剪切模量比和阻尼比隨剪應變衰減的關系曲線(簡稱模量衰減曲線)體現(xiàn)土體在動力作用下的非線性現(xiàn)象。為了描述土體在剪切應變增大時剪切模量衰減和阻尼比增加的現(xiàn)象，Seed利用模量衰減曲線概括了土體在動力作用下呈現(xiàn)的非線性，并給出了不同類型土壤的模量比和阻尼比隨剪應變的實驗值[8]。非線性土層剖面中采用Seed提出的砂土模量比和阻尼比實驗值見表 1。DEEPSOIL和FLAC3D分別采用了不同的方式來擬合模量衰減曲線。DEEPSOIL程序中內設的修正雙曲線模型允許通過適當定義模型參數(shù)的方式來考慮土體的非線性行為，在選擇Seed砂土模量衰減實驗值后可以通過擬合曲線程序自動填入非線性參數(shù)[9]。FLAC中采用滯后阻尼的方式模擬巖土材料的滯后特性，通過輸入擬合參數(shù)的方式將模量衰減曲線引入模型動力計算中。采用4參數(shù)模型(Sig4)輸入FLAC3D中的滯后阻尼，計算使用的阻尼參數(shù)分別為1.001，-0.609，-1.195和-0.001 2。DEEPSOIL和FLAC3D對Seed砂土實驗值的擬合效果如圖 2所示，擬合參數(shù)見表 1。從圖 2中可以看出，為了同時滿足模量比和阻尼比的擬合精度，DEEPSOIL擬合的剪切模量比在小應變時略小于目標值，而在大應變時略大于目標值；阻尼比在應變范圍內均略大于目標值。FLAC3D的擬合結果與DEEPSOIL基本一致。

表1 土層非線性參數(shù)Tab.1 Nonlinear calculation data of soil layers

DEEPSOIL中動力荷載輸入的方式分為內部運動和露頭運動2種。DEEPSOIL手冊中建議剛性基巖模型對應土層內部記錄的地震波，而柔性基巖模型對應露頭地震波。FLAC3D中動力荷載的輸入也分為這2種方式：當模型為剛性基巖時，直接在模型底部網(wǎng)格施加加速度時程；當模型為柔性基巖時，通過安靜邊界和應力時程的方式施加地震荷載。

選擇1940年Imperial Valley地震中 El Centro 記錄的南北分量作為輸入波。該地震記錄共有1 560個數(shù)據(jù)點，采樣頻率為50 Hz，在2.04 s時達到水平加速度峰值0.296 g，主導頻率為1.172 Hz，有效持續(xù)時間為31.16 s。計算時將加速度時程調幅至0.1 g作為輸入。調幅后El Centro地震波的加速度時程和傅里葉幅值譜如圖 3所示。

青浦區(qū)是上海市唯一與江蘇、浙江接壤的郊區(qū)，加強流域聯(lián)動合作顯得尤為重要。為了解決邊界水事糾紛，確保一方平安，青浦區(qū)積極開展與流域機構及兄弟省市的工作交流，初步與江蘇、浙江相鄰市縣建立了水務執(zhí)法聯(lián)席會議制度，2011年5月18日會同蘇州市水利局、嘉興市水利局、上海水務行政執(zhí)法總隊水利支隊、江蘇省太湖地區(qū)水利工程管理處等單位共同開展太湖流域“一湖兩河”聯(lián)合巡查活動，共同維護流域正常水事秩序及防汛安全。此外，在省市邊界河湖治理、供水工程管理等方面能夠積極溝通、主動協(xié)調。2012年3月，積極做好了浙江嘉善縣作為緊急備用水源的省際湖泊長白蕩相關協(xié)調工作。

2.2 臺站場地選擇及輸入

李兆焱等[10]對巨厚場地上地震記錄進行了地震反應分析，指出雖然DEEPSOIL在強震中表現(xiàn)不理想，但在較弱地震動下，DEEPSOIL計算的地表 PGA與加速度反應譜均接近實際情況。王鸞等[11]指出在中強地震動和較弱地震動情況下，DEEPSOIL在硬場地上的計算結果與實測反應譜較為相近。基于以上結論，從KiK-net地震臺網(wǎng)中選取2處硬土場地，篩選地表峰值在0.05 g左右的地震記錄用于計算和分析。

KiK-net地震臺網(wǎng)由日本國家地球科學與防災研究所建立，包含日本地區(qū)上百個地震臺站的記錄，可以根據(jù)時間、震級、地表峰值、臺站名稱等信息來篩選地震波，同時部分臺站在井下基巖與地表均設立了加速度傳感器，可以獲得包括場地土層剖面、地震波輸入與地表響應的完整數(shù)據(jù)[12]。在所選擇的硬土場地IWTH12臺站和IBRH11臺站的地震記錄中，分別選擇了地表峰值為0.056 g和0.053 g的2次地震用于計算。將地震記錄中的地表加速度記錄作為實測值，井下加速度記錄作為剛性基巖的內部運動輸入波。同時挑選距離計算土層距離較近且坐落于堅硬巖層上的IWTH02臺站和IBRH12臺站作為露頭運動參考點，取參考點所記錄同期地震的地表加速度記錄幅值的一半作為計算場地的露頭運動輸入波。對有2個方向的地震動記錄選取EW方向的記錄參加計算。由于KiK-net中未提供場地土層剖面的密度信息，參考文獻[13]中提出的壓縮波速vp與土體密度ρ的統(tǒng)計關系式進行計算:

ρ=0.31×vp0.25，

(2)

式中vp為壓縮波速。

計算臺站場地的臺站場地土層剖面和剪切波速分布圖如圖4所示，臺站場地及地震波基本信息見表 2。

(a) IWTH12臺站

根據(jù)場地剪切波速分布圖及土層剖面圖，分別選取深度-74.00、-30.00 m作為地震動輸入界面。FLAC3D土層反應模型的單元網(wǎng)格尺寸設為1.00 m，采用自由場邊界。土層模型所需土層參數(shù)如剪切模量G和體積模量K由剪切波速進行換算：

G=ρvs，

(3)

(4)

式中v為泊松比。

表2 臺站場地及地震波基本信息Tab.2 Basic information of hard site stations for contrast test

DEEPSOIL中提供了各類土的非線性參數(shù)。根據(jù)KiK-net公開資料中的土性描述來確定不同土類的非線性曲線。土層非線性參數(shù)見表3。

表3 土層非線性參數(shù)Tab.3 Nonlinear calculation data of soil layers

(5)

式中：Sa1、Sa2分別代表計算加速度反應譜和實測加速度反應譜；lgT為反應譜中周期從0.01～10 s的對數(shù)。易知反應譜計算值越接近于實測值，則譜面積比SR越接近于1。

3 計算結果與分析

3.1 簡化土層剖面對比分析

在土層地震反應分析結果的多個指標當中，加速度反應譜能夠綜合反映地震荷載的強弱和頻譜特性，是工程結構抗震中常用的指標之一。通過FLAC3D動力計算、DEEPSOIL頻域和時域分析3種方法計算了簡化土層剖面地表加速度反應譜，以此來對比分析3種方法計算結果的異同。各類型場地計算的地表加速度反應譜如圖5—7所示。

從圖中可以明顯看出，在簡化土層的等波速剖面和變波速剖面上，F(xiàn)LAC3D模型的地表加速度反應譜計算結果與DEEPSOIL頻域及時域分析結果基本一致，說明所建立的FLAC3D模型對粘彈性土層地表響應的計算結果是準確的。在簡化土層的非線性剖面中，F(xiàn)LAC3D模型的計算結果與DEEPSOIL時域分析結果較為接近，而與DEEPSOIL頻域分析結果在0.3～0.6 s周期范圍內存在一定差異，其余周期范圍變化趨勢基本一致。當輸入地震波的主導頻率與場地卓越頻率相近時，輸入地震波中的強振分量將在頻域計算中被放大，被稱為虛假共振現(xiàn)象。DEEPSOIL頻域分析的計算結果在0.3～0.6 s周期范圍內的譜峰值大于FLAC3D與DEEPSOIL時域分析結果，可能是由上述虛假共振現(xiàn)象造成的。

土層的地表反應受地震動輸入方式的影響。在同一類型的場地上，圖5—7中內部運動和露頭運動的計算結果表現(xiàn)出明顯的差異，這體現(xiàn)了基巖層與土層波阻抗比對地表地震反應的影響。另外通過不同類型土層剖面之間的對比也可以看出，土層非線性程度對地表響應有較大的影響，隨著土層非線性程度的增加，地表反應譜的形態(tài)、峰值以及峰值所在的周期都存在明顯差異。

(a) 內部運動

(a) 內部運動

(a) 內部運動

不同場地類型地表地震動幅值(peak ground acceleration, PGA)見表 4。由表 4可以看出，在相同場地類型和地震波輸入方式下，F(xiàn)LAC3D與DEEPSOIL的計算結果較為接近。隨著土層非線性程度的增加，采用露頭運動輸入的柔性基巖場地的地表PGA變化不大，而剛性基巖場地的地表PGA變化明顯。

表4 不同場地類型地表PGA計算值對比Tab.4 Surface PGA of various layers

3.2 臺站場地對比分析

IWTH12臺站和IBRH11臺站所在的場地等效剪切波速分別為312、242 m/s，覆蓋層厚度分別為20、30 m，屬于我國規(guī)范中的II類場地。在2個臺站數(shù)據(jù)中挑選2組地震波輸入來計算反應譜，得到結果如圖8(a)—11(a)所示。同時挑選距離計算土層距離較近且坐落于堅硬巖層上的臺站IWTH02和IBRH12作為柔性基巖計算的露頭運動參考點，結果如圖 8(b)—11(b)所示。為了對比線性和非線性土層地震分析方法對計算結果的影響，分別按照阻尼比為5%的黏彈性場地和引入模量衰減曲線的非線性場地進行計算，不同計算方法反應譜面積比對比見表5。

(a) 內部運動

(a) 內部運動

(a) 內部運動

(a) 內部運動

當輸入類型為內部運動時，從圖 8(a)可以看出，DEEPSOIL計算的黏彈性場地地表反應譜曲線與實測值在形狀上較為相似，但在譜峰值上存在明顯差異。對比圖 8(a)和圖 9(a)可以看出，在引入模量衰減曲線后，非線性場地的計算結果與實測值在峰值和形狀上都較黏彈性場地更為相近。IBRH11臺站也呈現(xiàn)了相似的現(xiàn)象。同時，DEEPSOIL時域和頻域計算結果在0.03～0.1 s周期范圍內一定程度上再現(xiàn)了實測地震記錄中的高頻成分。由表5可見，譜面積比SR同樣顯示了這一點，與黏彈性場地相比，當考慮土體的非線性時，IWTH12臺站計算的譜面積比從50%增加至60%左右，IBRH11臺站計算的譜面積比從50%增加至80%左右，表明了輸入地震波類型為內部運動時，考慮土體的非線性的計算結果與實測值更為接近。

當輸入類型為露頭運動時，對比圖 8(b)和圖 10(b)可以看出，IWTH12臺站的計算結果與實測值較為接近，而IBRH11臺站的計算結果與實測值有較明顯的差異。這可能是由于所選擇的露頭運動參考點與計算場地的距離不同導致的。類似地，與黏彈性場地的計算結果相比，非線性場地的計算結果在0.1～0.3 s周期范圍內與實測值更為相似。由表5可見，譜面積比SR顯示，IWTH12臺站非線性反應譜面積比從45%增加至50%左右，IBRH11臺站反應譜面積比從40%增加至60%左右，表明輸入地震波類型為露頭運動時，非線性場地的計算效果要略優(yōu)于黏彈性場地。

地表PGA實測值與計算值對比見表 6。由表 6可以看出，F(xiàn)LAC3D的計算結果與DEEPSOIL十分接近，但均明顯小于實測地表PGA。非線性場地的計算結果要比黏彈性場地更為接近實測值，但仍然有明顯差距。

表5 不同計算方法反應譜面積比對比Tab.5 Spectrum ratio calculated by different methods

表6 Kik-net臺站場地計算PGA結果對比Tab.6 Comparisons of PGA resulfs of KiK-net sites

4 結語

基于不同非線性程度的人造剖面和日本KiK-net臺網(wǎng)的實際場地及地震記錄數(shù)據(jù)，從地表反應譜和PGA兩方面對FLAC3D和DEEPSOIL程序進行了對比分析，得出以下主要結論：

①在簡化土層剖面計算中，F(xiàn)LAC3D與DEEPSOIL計算結果較為接近。當場地卓越頻率與地震波主導頻率相接近時，DEEPSOIL頻域分析在非線性場地剖面上的反應譜計算結果大于FLAC3D與DEEPSOIL時域分析，這是頻域分析方法的虛假共振現(xiàn)象所致。

②在實測場地反應譜的對比中，考慮了土體非線性的場地計算結果較黏彈性場地更為接近實測值，其中剛性基巖場地的計算結果在短周期范圍較為準確，而采用露頭運動參考點的柔性基巖場地計算結果在長周期范圍更為合理。

③與實測場地地表PGA相比，黏彈性模型計算的地表PGA偏小，考慮了土體非線性的場地計算結果較為良好，但在土層分布復雜的硬土場地上仍然與實測記錄有著較大差距。