中硬場地下兩種土層地震反應方法與精確解的對比①

李瑞山， 袁曉銘， 李程程

(中國地震局工程力學研究所地震工程與工程振動重點實驗室，黑龍江 哈爾濱 150080)

中硬場地下兩種土層地震反應方法與精確解的對比①

李瑞山， 袁曉銘， 李程程

(中國地震局工程力學研究所地震工程與工程振動重點實驗室，黑龍江 哈爾濱 150080)

基于諧波入射下的波動理論頻域精確解，導出線性時域精確解，并在Matlab環(huán)境中編制相應的計算程序；選取8個簡化的中硬場地剖面，用LSSRLI-1、精確解和SHAKE2000三種方法計算各場地在不同輸入條件下的地震反應。結果表明：程序計算所得地表反應譜和土體剪應變分布與SHAKE2000結果一致；LSSRLI-1方法得到的地表反應譜與前二者結果一致；LSSRLI-1方法在某些情況下得到的土體剪應變分布與另外二者結果存在較大偏差，該偏差對地表反應譜有著不可忽略甚至非常顯著的影響。

地震反應分析； LSSRLI-1； 精確解； 土表反應譜； 土體剪應變

0 引言

土層地震反應分析是巖土地震工程領域的一個重要分支，主要研究場地對土體中傳播的地震波的影響。土層場地在承載上部結構的同時還要充當地震波由基巖向地表傳播的媒介。建筑結構在地震作用下的響應分析需要以地面運動作為輸入，所以地表運動預測的精度直接關系到結構地震反應分析的可靠性。土層結構對地表加速度峰值、反應譜平臺值和特征周期都有直接的影響[1-3]，通過土層地震反應分析得到的土體剪應變分布以及地表運動特征可以從一定程度上預測場地在地震作用下液化的可能性[4]。

現有的土層地震反應分析方法主要可以分成線性、等效線性和非線性三類。土體在動剪應力作用下往往表現出強烈的應力-應變非線性關系，鑒于此，線性方法在實際分析當中很少被使用。非線性方法需要在時域當中通過微小的時間步積分來得到相應的分析結果，所以在實際操作中非常耗時。不僅如此，非線性方法的分析精度在很大程度上依賴于土體動應力-動應變本構關系的正確與否，而土體動力本構本身就非常復雜。等效線性化方法利用線性方法，數學概念簡單，物理意義明確，并將土體動力非線性性質引入到了分析當中，在不背離其本質的前提下大大方便了計算。等效線性化方法在總體動力學效應大致相當的意義上用一個等效的剪切模量和阻尼比代替所有不同應變振幅下的剪切模量和阻尼比，將非線性問題轉化為線性問題[5]。自Seed和Idriss于1970年首次將等效線性化的概念引入巖土地震工程領域，并于1972年編制了相應的計算程序——SHAKE[6]，現在SHAKE版本已更新至SHAKE2000，是國際主流的土層反應分析等效線性化計算程序。目前我國地震安全評價所使用的是LSSRLI-1程序[7-8]，該程序自提出至今已有20多年，代表著當時的國際先進水平，為我國的防震減災事業(yè)做出了卓越的貢獻。但是，在實際使用過程當中該程序也表現出了一些缺點和不足[9]，尤其是當土層較軟或輸入地震運動較強時，其計算結果不合理性較為明顯。

本文擬復現一維土層地震反應分析頻域精確解，并利用快速傅里葉變換(FFT)和逆傅里葉變換(IFFT)實現幅值和頻率成分都極不規(guī)則的隨機地震荷載在時域和頻域當中的相互轉換。基于以上理論，結合傅里葉變換的線性可疊加性，給出土層地震反應分析一維線性精確解，并在Matlab環(huán)境中編制相應的計算程序。最后利用該程序進行中硬場地在不同輸入條件下的地震反應計算，并將所得結果與SHAKE2000和LSSRLI-1進行對比。為了和本文的線性精確解保持一致性和可對比性，SHAKE2000和LSSRLI-1在計算時均采用線性不迭代方式。

1 線性精確解

假設實際場地水平成層，土性為各向同性的Kelvin-Voigt黏彈性體，輸入的地震運動為由基巖垂直向上入射的剪切波。各土層由其層厚h、重度ρ、初始剪切模量G、初始阻尼比ξ來表征，如圖1所示。各層的局部坐標原點位于層頂，方向向下。土體中的波動用數學方程可以表達為：

(1)

在頻率為ω的諧波入射下，式(1)的解為：

圖1 水平成層簡化模型Fig.1 Simplified model of the horizontally layered site

(2)

具體對第n層土，以其空間坐標zn和時間坐標t所表達的位移函數為：

(3)

式中，E和F分別為土層中所傳播的上行波和下行波波幅，k為波數，其表達式為：

(4)

根據分層界面上位移和應力的連續(xù)條件，可以導出相鄰兩層波幅之間的數量關系：

(5)

用矩陣可以簡單表達為：

(6)

根據輸入，通過遞歸便可得到各層的波幅矢量。然后根據位移場表達式(3)和幾何方程，可以求得土層中剪應變表達式為：

(7)

根據先前求得的各層波幅矢量便可很容易得到土體當中剪應變的分布。另外，為了考慮土體等效阻尼，采用復剛度G*代替原來的剪切剛度G，剪切波速υS和波數k隨之相應變化，不再贅述。

(8)

基巖輸入的剪切波經快速傅里葉變換可以分解為一系列諧波的疊加，將各自分別輸入所得到的響應量對應疊加后進行逆傅里葉變換便可得到相應的時域反應。

本文根據上述理論推導，在Matlab環(huán)境當中編制相應的計算程序，用來求解基巖輸入加速度時程時的地表響應和土體剪應變分布。為了更加直觀地反應輸入參數相同時輸出量的差異性，計算采取線性運算，不進行迭代。

2 計算模型

本文采用單層土和雙層土兩種簡化計算模型，共計8個計算剖面，其中4個單層剖面，4個雙層剖面，各剖面的土層參數詳見表1，其中υS1和υS2分別為上下兩層土的剪切波速，h1和h2分別為上下兩層厚度。為分析輸入地震動頻譜成分對計算結果差異性的影響，選取三條頻譜特性各異的地震動時程作為輸入，分別為ElCentro波、AKTH19波和KSRH09波，其中AKTH19波和KSRH09波是從KiK-net數據庫中選取的井下臺站實測數據。AKTH19波的主要頻譜特點是低頻成分絕對占優(yōu)，高頻分量很少，而KSRH09波相比而言具有十分豐富的高頻成分。計算時分別將三條加速度時程調幅至0.1g、0.2g和0.4g作為輸入，以對比三種方法計算結果差異性與輸入地震動強烈程度之間的關系。調幅至0.2g時的加速度時程曲線如圖2所示。

圖2 輸入加速度時程曲線Fig.2 Input acceleration time histories

在采用等效線性化計算時，某次循環(huán)所使用的動力非線性參數是根據上一次循環(huán)所得到的土體剪應變分布而定的，等效剪應變一般可按下式計算：

(9)

其中，γmax為剪應變時程的幅值。按式(9)算得等效剪應變之后，便可從模量比和阻尼比隨剪應變變化曲線上插值得到新的G/Gmax和ξ值進行下一次迭代計算，直到相鄰兩次計算結果相對差值小于某一容許值。

土層地震反應分析結果能夠為工程實際所直接應用的是地表運動響應，包括運動幅值及其反應譜等。地表運動與土地剪應變分布有著密切的關系。為了直觀表現剪應變差異對地表加速度反應譜的影響，本文將用線性方法得到的剪應變按等效線性化的思想插值后求得新的動力非線性參數，并利用新參數重新計算地表響應。這樣剪應變的差別便通過重新計算的地表反應譜差別表現出來。文中插值所用的是孫銳等[10]通過大量共振柱試驗統(tǒng)計得出的粉質黏土非線性平均線，具體見圖3。

圖3 剪切模量和阻尼比隨剪應變幅的變化曲線Fig.3 Variation of shear modulus and damping ratio with shear strain amplitude

表1 計算模型參數表

Table1 Parameters of the calculation models

單層模型場地編號υS/(m·s-1)h/m112010220024312050420080雙層模型場地編號vS1/(m·s-1)h1/mvS2/(m·s-1)h2/m51003150661508250167120181803282002630054

3 結果對比分析

土層地震反應分析所得到的地表加速度反應譜可以為工程實際所直接應用，其計算的正確與否直接關系到對應的計算方法的可靠性。本文共選取72種計算工況，但因為采取不迭代的計算方式，實際獨立的工況只有24種。限于篇幅，圖4中僅繪制了在輸入波峰值為0.2 g時其中8個工況所計算得到的地表加速度反應譜。從圖中可以清晰看出由精確解、SHAKE2000和LSSRLI-1三種方法計算得到的地表加速度反應譜幾乎完全相同，這充分說明LSSRLI-1計算地表加速度的方法是可靠的，同時也說明本文在Matlab環(huán)境當中編制的精確解程序在輸入參數相同的情況下所計算出的地表運動同國際主流的SHAKE2000結果是一致的。

圖4(c)、(e)和(g)工況的輸入幅值均為0.2 g的KSRH09波，得到的地表加速度反應譜從峰值個數和大體形狀上來講很相似，但很明顯峰值大小和相對應的周期各不相同，這充分體現了場地條件對輸入加速度波的頻譜特性的改變，不同場地在相同輸入下的響應不盡相同。從圖4(b)、(d)和(f)三個工況對比可以得到相同的結論。另外，從圖4中還可以看出，在不同的基巖輸入下得到的地表反應譜相互之間存在很大的差異性，這說明地表輸出同時受到基巖輸入和場地條件的影響。

圖5是土體剪應變隨深度的分布，其中輸入波峰值為0.2 g，本文選取24個獨立計算工況中的8個作為代表來呈現。從圖中可以看出，本文編制的精確解程序計算得到的土體剪應變分布與SHAKE2000結果非常吻合，這種一致性與場地和輸入條件均無關，說明SHAKE2000剪應變計算結果完全可以由本文程序來代表。綜合上面三種方法計算得到的地表反應譜的一致性可以得出結論：SHAKE2000計算出的土體剪應變和地表反應譜同本文根據精確解所編制程序的計算結果均一致，本文程序與國際主流的SHAKE2000具有相同的精度。下文中SHAKE2000結果均由精確解來代表，而不再單獨進行對比分析。

從圖5中可以明顯看出LSSRLI-1程序計算出的土體剪應變分布與精確解結果存在一定差異。圖5(d)和(f)兩工況當中剪應變的差別很小，圖5(a)、(b)和(e)三個工況呈現出不可忽略的剪應變差別，剩下的圖5(c)、(g)和(h)三個工況的剪應變差異性十分顯著，并且均為LSSRLI-1結果遠大于精確解結果。下面就影響剪應變相對差別大小的不同因素進行分析。

首先，整體剪應變水平的大小直接影響著二者剪應變差異性的大小。為此將圖示8個工況分為三組，圖5(e)和(f)為低應變組，該組呈現出的剪應變之間有一定的差異性，但并不顯著，尤其對于分層模型中的上層，三種方法得到的剪應變幾乎完全一致；圖5(a)，(b)和(h)為中應變組，該組剪應變差別較低應變組明顯有所增大，甚至如圖5(h)所示非常顯著；圖5(c)，(d)和(g)為大應變組，剪應變差別也更為明顯。一般地，剪應變水平越高時，LSSRLI-1計算得到的土體剪應變分布與精確解偏差越顯著。

圖4 地表加速度反應譜(輸入波峰值=0.2 g)Fig.4 Ground acceleration response spectrum (peak acceleration=0.2 g)

圖5 土體剪應變分布(輸入波峰值=0.2 g)Fig.5 Distribution of soil shear strains (peak acceleration=0.2 g)

另外，輸入波的頻譜特性對二者計算得到的土體剪應變差別也有很大的影響。三條輸入波中，按高頻成分多少依次為AKTH19、El Centro和KSRH09。在AKTH19波輸入時(圖5(b)、(d)和(f)，LSSRLI-1所得剪應變在一定程度上較精確解結果偏大或偏小，但偏差不十分顯著。當高頻成分相對較多的El Centro輸入時，LSSRLI-1剪應變較精確解結果均偏大，有時偏差甚至非常顯著。將三條波中高頻成分最多的KSRH09波作為輸入波時，LSSRLI-1剪應變比精確解結果嚴重偏大，二者相對差值很多時候超過100%。由此可見，LSSRLI-1剪應變計算的適用性與輸入波有很大關系，低頻成分較豐富的波輸入時，計算結果偏差較小，計算方法適用性較強；高頻成分豐富的波輸入時，計算所得剪應變偏差較為顯著，計算方法適用性較差。

但是剪應變不能直接被工程實際所使用，通過其差異性無法依據工程可接受標準直觀地判斷某種計算方法應用于土層地震反應分析時是否可靠。為了用剪應變差別來影射其對地表反應譜的影響，本文利用圖5中所示的剪應變，按式(9)算得各層的等效剪應變，然后按非線性曲線(圖3)插值得到新的土層參數，重新計算得到地表加速度反應譜(圖6)。因圖5中的SHAKE2000結果和精確解結果幾乎相同，所以在圖6中沒有繪出SHAKE2000結果。

圖6 使用新參數計算出的地表加速度反應譜Fig.6 Ground acceleration response sprctrum calculated by using the updated parameters

結合圖5，從圖6中可以明顯看出，中硬場地剪應變計算結果偏差越大，對地表反應譜的影響就越明顯。圖5中(a)、(b)、(d)、(e)和(f)所示的剪應變差別對地表反應譜有一定影響，但影響不大，在工程可接受范圍之內。圖5中(e)、(g)和(h)三個工況所表現出的剪應變偏差十分顯著，對應在圖6中表現出十分顯著的地表反應譜差別。另外圖5中(a)、(b)剪應變差別明顯大于圖5(d)，而圖6并未表現出(a)、(b)工況差別明顯大于(d)工況的現象，究其原因是因為(d)工況的剪應變水平較(a)、(b)兩工況要高很多，不同的應變水平對非線性的敏感程度不同。綜合上述兩點可以得出如下結論：剪應變差別對地表反應譜的影響同時受偏差大小和應變整體水平兩個因素控制，但其中剪應變差別占主導。

4 結論

本文復現了諧波輸入下的土層反應頻域精確解，在此基礎上導出了時域線性精確解，并在Matlab環(huán)境中編制了相應的計算程序。文中共選取了8個簡化的中硬場地剖面，選取了三條頻譜特性各異的加速度時程作為輸入。經對比分析由精確解、SHAKE2000和LSSRLI-1分別計算出的土表加速度反應譜和土體剪應變分布，可以得出如下結論：

(1) 本文根據導出的精確解所編制的計算程序具有與SHAKE2000所代表的國際主流等效線性化土層地震反應分析方法同樣的精度，二者計算得到的土表加速度反應譜和土體剪應變分布完全一致。

(2) LSSRLI-1計算得到的土表加速度反應譜與精確解和SHAKE2000結果一致，LSSRLI-1計算地表運動的相應模塊是可靠的。

(3) LSSRLI-1計算所得剪應變與精確解結果之間存在一定的偏差，但大部分情況下偏差很小，由此導致的地表反應譜偏差在工程可接受范圍之內。

(4) 當輸入波高頻成分較豐富時，LSSRLI-1計算所得剪應變可能會出現較為顯著的偏差，對地表反應譜影響十分明顯，不可忽略，說明該方法對部分情況有一定的不適用性，剪應變的計算方法亟需改進。

References)

[1] 薄景山，李秀領，劉德東，等.土層結構對反應譜平臺值的影響[J].地震工程與工程振動，2003，23(4)：29-33.BO Jing-shan，LI Xiu-ling，LIU De-dong，et al.Effects of Soil Layer Construction on Platform Value of Response Spectra[J].Earthquake Engineering and Engineering Vibration，2003，23(4)：29-33.(in Chinese)

[2] 薄景山，李秀領，劉德東，等.土層結構對反應譜特征周期的影響[J].地震工程與工程振動，2003，23(5)：42-45.BO Jing-shan，LI Xiu-ling，LIU De-dong，et al.Effects of Soil Layer Construction on Characteristic Periods of Response Spectra[J].Earthquake Engineering and Engineering Vibration， 2003，23(5)：42-45.(in Chinese)[3] 薄景山，李秀領，劉紅帥.土層結構對地表加速度峰值的影響[J].地震工程與工程振動，2003，23(3)：35-40.BO Jing-shan，LI Xiu-ling，LIU Hong-shuai.Effects of Soil Layer Construction on Peak Accelerations of Ground Motions [J].Earthquake Engineering and Engineering Vibration，2003，23(3)：35-40.(in Chinese)

[4] 孫銳，袁曉銘.液化土層地震動特征分析[J].巖土工程學報，2004，26(5)：684-690.SUN Rui，YUAN Xiao-ming.Analysis on Featrue of Surface Ground Motion for Liquefied Soil Layer[J].Chinese Journal of Geotechnical Engineering，2004，26(5)：684-690.(in Chinese)

[5] 廖振鵬.工程波動理論導論[M].北京：科學出版社，2002.LIAO Zhen-peng.Introduction to Wave Motion Theory for Engineering[M].Beijing：Science Press，2002.(in Chinese)

[6] Shnabel P，Lysmer J，Seed H B.SHAKE：A Computer Program for Earthquake Response Analysis of Horizontal Layer Sites[R].1972.

[7] 廖振鵬，李小軍.地表土層地震反應的等效線性化解法[M]//廖振鵬.地震小區(qū)劃(理論與實踐).北京：地震出版社，1989：141-153.LIAO Zhen-peng，LI Xiao-jun.Linearization Method for Computing Earthquake Response of Ground Layered Soil[M]//LIAO Zhenpeng.Seismic Microzonation.Beijing：Seismological Press，1989，：141-153.(in Chinese)

[8] 李小軍.一維土層地震反應線性化計算程序[M]//廖振鵬，地震小區(qū)劃(理論與實踐)北京：地震出版社，1989：250-265.LI Xiao-jun.A Computer Program for Calculating Earthquake Response of Ground Layered Soil[M]//Liao Zhen-peng.Seismic Microzonation.Beijing：Seismological Press，1989：250-265.(in Chinese)

[9] 劉德東，齊文浩，張宇東，等.現行土層地震反應分析存在的問題[J].防災科技學院學報，2009，21(3)：34-37.LIU De-dong，QI Wen-hao，Zhang Yudong，et al.Problems Existing in Current Seismic Response Analysis for Soil Layers[J].J of Institute of Disaster-Prevention Science and Technology，2009，11(3)：34-37.(in Chinese)

[10] 孫銳，陳紅娟，袁曉銘.土的非線性動剪切模量比和阻尼比不確定性分析[J].巖土工程學報，2010，32(8)：1228-1235.SUN Rui，CHEN Hong-juan，YUAN Xiao-ming.Uncertainty of Non-linear Dynamic Shear Modular Ratio and Damping Ratio of Soils[J].Chinese Journal of Geotechnical Engineering，2010，32(8)：1228-1235.(in Chinese)

Contrasting Study of Two Seismic Response Methods and the Exact Solutions for Moderately Stiff Soil Sites

LI Rui-shan， YUAN Xiao-ming， LI Chieng-cheng

(KeyLaboratoryofEarthquakeEngineeringandEngineeringVibration，InstituteofEngineeringMechanics，ChinaEarthquakeAdministration，Harbin，Heilongjiang100080，China)

Earthquake ground motions can significantly affect buildings when they are amplified in surface soil layers.Therefore，in order to estimate the behavior of buildings during severe earthquakes，it is essential to evaluate those characteristics of amplification.Although investigations have been carried out，many problems remain unsolved with regard to those characteristics.Ground response analysis typically involves predictions of surface ground motions，developing a spectrum of design responses，determining dynamic stress and strain，and evaluating the potential for liquefaction.Methods used to calculate the ground response are commonly divided into linear，nonlinear，and equivalent linear concepts.The linear analysis assumes that the shear modulus as constant.However，soils often exhibit nonlinear behavior，even at low levels of strain.The nonlinear method is time consuming.The equivalent linear method，represented by the SHAKE model，is one of the most frequently used methods for analysis of ground responses to earthquakes.This method is used primarily for approximations because it uses effective strain to define the material property for analysis.However，the equivalent linear method remains the most popular method of seismic response analysis.It has been 20 years since the safety evaluation program LSSRLI-1 was developed and was used for seismic response analysis.At the time，LSSRLI-1 represented the most advanced method and contributed greatly to earthquake preparedness and disaster reduction for China.However，some deficiencies in the model were identified.Based on its focus in the frequency domain with a harmonic wave incident，the exact linear solution was deduced in the time domain，and compiled into a program in the Matlab environment.For the study，eight sites with medium-hard soils were chosen.Seismic responses were calculated using three different methods：the LSSRLI-1，exact solution，and SHAKE2000.The results show that the ground response spectrums and distribution of soil shear strain，calculated using the exact solution and SHAKE2000 methods，are in agreement with each other.The ground response spectrums，calculated using LSSRLI-1，were consistent with the other two methods.However，in some cases，the distribution of soil shear strain，calculated using LSSRLI-1，showed a relatively large deviation from the SHAKE2000 and exact solution methods.The extent of these differences may significantly affect the estimates of ground motion，and cannot be ignored.

seismic response analysis； LSSRLI-1； exact solution； response spectrum of soil surface； soil shear strain

2014-08-20

國家科技部地震行業(yè)專項(201308015)；黑龍江省自然科學基金資助項目留學歸國基金(LC2013C14)；國家自然科學基金(51278472)

李瑞山(1987-)，男，博士研究生，主要從事土動力學和巖土地震工程研究.E-mail：lrshan22@hotmail.com

TU43

A

1000-0844(2015)02-0565-06

10.3969/j.issn.1000-0844.2015.02.0565