土層地震反應(yīng)顯式計算中阻尼矩陣系數(shù)的選取

樓夢麟，邵新剛

（同濟大學(xué)土木工程防災(zāi)國家重點實驗室，上海200092）

在進行土層地震反應(yīng)的確定性分析時，通常有兩種分析方法，即時域分析和頻域分析.在進行時域分析時，一般采用粘滯阻尼假定，即阻尼力與質(zhì)點的運動速度成正比，通過求解實系數(shù)微分方程組，可求得各自由度反應(yīng)的時程.而在頻域分析時，通常采用Fourier變換，通過求解復(fù)代數(shù)方程組，可求得土層地震反應(yīng)的穩(wěn)態(tài)解［1］，此時一般采用滯后阻尼假定.這兩種分析方法各有不同的適用范圍，采用頻域法時，可以方便地對線性體系的地震反應(yīng)進行計算，但對于非線性地震反應(yīng)的分析只能采用等效線性化方法［2］，無法進行真非線性分析.而時域分析方法能夠適應(yīng)非線性體系地震反應(yīng)分析的不同要求，既可采用等效線性化的分析方法，也可采用真非線性分析方法［3］.除了簡單結(jié)構(gòu)如梁、柱等外，一般工程動力體系地震反應(yīng)的時域計算多采用直接積分法［4］，直接積分法一般又分為隱式算法和顯式算法［5－6］.例如：對于無阻尼體系或阻尼矩陣為一對角矩陣的情況，中心差分方法所給出的格式為顯式.但是對于一般非對角線的阻尼矩陣，中心差分方法所給出的格式為隱式.在LS－DYNA［7］程序中對土層非線性動力分析時具有良好的非線性本構(gòu)模型，且采用時域顯式算法，具有較高計算效率，得到較為廣泛的應(yīng)用.然而，程序中土層計算體系的阻尼矩陣采用了僅與質(zhì)量矩陣成線性比例的假定.由于集中質(zhì)量矩陣為對角陣，阻尼矩陣也為對角陣，滿足了顯式算法的要求.

已有研究工作［8］表明，對于深覆蓋土層地震反應(yīng)計算采用土動力學(xué)中常用的與剛度矩陣成線性比例的阻尼矩陣將會產(chǎn)生低估土層地震反應(yīng)的問題.本文通過線性條件下的算例分析，討論LS－DYNA計算程序中基于質(zhì)量比例阻尼矩陣顯式算法所存在的問題，并進一步以頻域反應(yīng)的解為準(zhǔn)確解，通過對比不同工況的時域解，討論確定質(zhì)量比例阻尼系數(shù)的取值問題.

1 質(zhì)量比例阻尼矩陣存在的問題

對于實際結(jié)構(gòu)，在進行動力反應(yīng)分析時，采用有限元方法對計算土層體系進行離散化處理后，可以得到如下動力方程式［9］：

式中：M，C 和K 分別為土層的質(zhì)量矩陣、阻尼矩陣和剛度矩陣分別為對應(yīng)于土層下臥基巖面的相對位移向量，速度向量和加速度向量；p（t）為作用于土層體系上的地震等效動力荷載向量，其 中：ej為位移變換向量，aj（t）是作用在基巖上的加速度歷程.

在有限元分析中，一般取質(zhì)量矩陣M 為對角矩陣.當(dāng)C＝a0M 時（a0為比例系數(shù)），阻尼矩陣也為對角矩陣，因此，式（1）中可采用顯式算法求解.DYNA程序中建議：

式中：ω1為土層的基頻；ξ為土層的振型阻尼比；一般假定土層各階振型阻尼比相同.以ω1來確定比例系數(shù)a0符合人們的常規(guī)的認(rèn)識，即在結(jié)構(gòu)體系地震反應(yīng)中，第1階振型頻率貢獻最大.

設(shè)：

式中：Φ 為計算土層的振型模態(tài)矩陣；q（t）為廣義坐標(biāo)向量.應(yīng)用式（3）對式（1）進行坐標(biāo)變換，并利用振型正交性可得：

令

式中：ηix，ηiy，ηiz分別稱為i振型在土層x，y和z方向地震反應(yīng)中的參與系數(shù)（或貢獻系數(shù)），簡稱i振型j方向地震參與系數(shù)，φi為振型向量.一般來說低階振型的地震參與系數(shù)的絕對值大于高階振型的地震參與系數(shù).根據(jù)式（2），可得第i階振型的計算阻尼比為

很顯然，當(dāng)假定體系各階振型阻尼比都為ξ時，按質(zhì)量比例阻尼矩陣的假定，則高階振型阻尼比的計算值小于假定值ξ，折減比例系數(shù)為ω1／ωi，振型階序越高，振型阻尼比計算值越小，使得該振型地震反應(yīng)的計算值就越大.

當(dāng)土層基頻接近于或高于輸入地震的主要頻率分量時，顯然由于土層第1階振型的振型參與系數(shù)大，其對體系地震反應(yīng)的貢獻遠大于其他高階振型，因此當(dāng)?shù)?階振型阻尼比與實際阻尼比相等時，在第1階振型地震反應(yīng)計算正確得到保證的前提下，其他高階振型阻尼比偏小時使得土層地震反應(yīng)的計算誤差非常有限.而當(dāng)計算體系的基頻遠低于輸入地震波的主要頻率分量時，很顯然與地震波中主要頻率分量相近的土層高階振型將存在共振效應(yīng).眾所周知，共振放大因子為Di＝1／（2ξi），低估振型阻尼比，將高估共振放大因子，從而使得該振型的地震反應(yīng)計算值增大.如當(dāng)ω1／ωi＝1／2時，則Di增大1倍，依此類推，此時高階振型的共振效應(yīng)不可忽略.下面通過具體算例，討論質(zhì)量比例阻尼矩陣給土層地震反應(yīng)計算帶來的影響.由于適應(yīng)LS－DYNA 的程序時域顯式算法中C＝a0M的規(guī)定，本文建議質(zhì)量比例阻尼系數(shù)a0表示如下：

式中：ωs為控制圓頻率；fs為控制頻率.一般來說，當(dāng)土層基頻接近于或高于輸入地震波的主要分量的頻率時，按常規(guī)取ωs＝ω1，而本文的重點討論當(dāng)土層基頻遠低于輸入地震波主要分量的迫振頻率時ωs（或fs）的合理取值.

2 數(shù)值分析

2.1 土層計算模型

文中土層模型采用港珠澳隧道場地模型的縱剖面為分析對象，水平范圍為1 016m，底部為基巖面，左端厚度為85m，右端厚度為112m，假定為均勻土層，土介質(zhì)的計算參數(shù)選取為：彈性模量391 MPa、密度2 000kg·m－3、泊松比0.35，這一土層模型稱為土層模型1；為了考察土層厚度對ωs值取值的影響，在這一土層模型的基礎(chǔ)上，采用比例按同一比例減小土層左右兩端厚度另外形成兩個土層模型，分別稱為土層模型2和土層模型3.此外，另建立一個4層土介質(zhì)組成的水平均勻分層土層模型（簡稱土層模型4）與均質(zhì)斜坡土層相對應(yīng)，自表層向下各層土介質(zhì)的計算參數(shù)為：彈性模量185.1，391，561.3和2 492.4 MPa，質(zhì)量密度1 870，2 000，2 220和2 100 kg·m－3，泊松比0.4，0.35，0.35和0.35.土層模型的有限元網(wǎng)格如圖1所示.

圖1 土層有限元網(wǎng)格圖Fig.1 Finite element mesh of soil layer

4個土層計算模型的幾何尺寸和土層前5階自振頻率見表1，從表中可以看出，土層模型1和土層模型4的自振頻率較為接近.

表1 土層計算模型特征長度與自振頻率Tab.1 Characteristic length and natural frequency of the soil layer calculation model

2.2 基巖面輸入地震波

為了比較全面地考察基巖面輸入地震波的影響，在土層地震反應(yīng)計算中共選用了5條地震波，包括該工程場地地震安全性評估報告給出的基巖人工波和El Centro、Kobe、Taft、汶川基巖波等4條地震實測波，它們的時程、Fourier幅值譜和反應(yīng)譜分別如圖2—6所示.

圖2 基巖人工波Fig.2 Bedrock artificial wave

為了考察地震波頻譜分布對質(zhì)量比例阻尼系數(shù)取值的影響，選取了一下幾個頻率參數(shù)作為地震波的特征頻率參數(shù)，分別為：①反應(yīng)譜峰值頻率fR；②反應(yīng)譜譜曲線重心頻率fRg；③Fourier幅值譜峰值頻率fF；此外還引入考慮結(jié)構(gòu)基頻和地震波特征頻率綜合影響的頻率參數(shù)fa和fb；④fa＝（f1＋fR）／2，f1＝w1／2π；⑤fb＝（f1＋fF）／2；⑥文獻［5］中滯后阻尼轉(zhuǎn)換頻率fc.計算中輸入地震波的加速度峰值都調(diào)整為0.1g（g為重力加速度），土層各階振型頻率阻尼比取為0.08，對應(yīng)于5條輸入地震波的不同特征頻率值見表2—4.

分析表中5條地震波的特征頻率參數(shù)fR，fF和fRg可以看出所選的汶川基巖地震波由于頻譜成分分布很廣，所以fR、fF值很高，是比較特殊的一條地震波，表明高頻成分很豐富，同樣基巖人工波中高頻成分也較為豐富.

表2 對應(yīng)于土層模型的特征頻率值Tab.2 Characteristic frequency value corresponding to the model of the soil layer Hz

表3 土層模型1中地表中點加速度與位移反應(yīng)幅值Tab.3 Acceleration and displacement response amplitude of the surface midpoint on Soil Model 1

表4 土層模型1中地表加速度與位移計算誤差Tab.4 Acceleration and displacement errors of calculation of the surface midpoint on Soil Model 1 %

2.3 土層地震反應(yīng)分析

對于同一土層有限元模型，分別采用時域和頻域分析方法，頻域計算時采用了滯后阻尼假定，因此以頻域計算所得土層反應(yīng)為準(zhǔn)確值，作為不同時域分析結(jié)果的比較基準(zhǔn).在時域分析中，不同工況計算所得的反應(yīng)結(jié)果與頻域的計算結(jié)果相比較，時域近似分析結(jié)果的相對誤差可表示為

式中：a和a＊為土層同一節(jié)點處地震位移或加速度反應(yīng)峰值，其中a為時域分析結(jié)果，a＊為頻域分析結(jié)果.

表3和表4中列出了土層模型1在不同地震波激勵下的土層地表中點加速度和位移反應(yīng)幅值及相應(yīng)的時域計算結(jié)果的誤差.

從表4中可以看出：由于土層模型1 的基頻遠低于表2中5條地震波的特征頻率fR、fF以及fRg，所以當(dāng)取ωs＝ω1形成土層地震反應(yīng)計算中的阻尼矩陣時，顯然高估了土層地表的加速度峰值以及位移峰值.相比之下取其他特征頻率時，土層地震反應(yīng)峰值計算結(jié)果要好于ωs＝ω1時的結(jié)果.進一步對比表2和表4中的數(shù)據(jù)可以看出，總體上ωs＝ωa＝2πfa時，計算精度要好一些，且基本能保證地震加速度峰值為計算誤差為正，即計算值略大于實際值.這里主要考慮到：一般來說土層地震反應(yīng)計算的目的是為地面結(jié)構(gòu)地震反應(yīng)分析提供結(jié)構(gòu)基底地震動加速度時程.

圖7和圖8中比較了在El Centro波和汶川波激勵下的土層模型1 地表中點加速度反應(yīng)時程曲線，從圖中可以看出，采用ωs＝ω1形成阻尼矩陣后，計算所得的土層地表中點加速度反應(yīng)時程與頻域計算所得的時程相差較大，而采用ωs＝ωa形成阻尼矩陣后（ωa＝2πfa），計算所得的土層地表中點加速度反應(yīng)時程得到了較好改善.其他地震波激勵下的土層地表中點加速度反應(yīng)時程曲線與El Centro波和汶川波激勵下的情況類似，限于篇幅，不一一列出.

表5中列出了土層模型2在不同地震波激勵下土層地表中點加速度和位移反應(yīng)的時域計算結(jié)果的相對誤差.

表5 土層模型2中地表加速度與位移計算誤差Tab.5 Acceleration and displacement errors of calculation of the surface midpoint on Soil Model 2 %

從表5可以看出：由于土層模型2的基頻總體 上略低于表2中5條地震波的特征頻率，因此誤差變化不及表4 中顯著，但取ωs＝ω1時，依然高估了土層的地震反應(yīng)值，從加速度峰值誤差分布來看，取ωs＝ωa時誤差水平最低.采用ωs＝ω1及ωs＝ωa形成阻尼矩陣后，El Centro波和汶川波激勵下的土層地表中點加速度反應(yīng)時程在相位與波形上非常接近，區(qū)別體現(xiàn)在峰值上，且在土層模型2中的誤差不如模型1顯著了.限于篇幅時程圖不再一一列出.

表6中列出了土層模型3在不同地震波激勵下的土層地表中點加速度和位移反應(yīng)的時域計算結(jié)果的誤差.

表6 土層模型3中地表加速度與位移計算誤差Tab.6 Acceleration and displacement errors of calculation of the surface midpoint on Soil Model 3 %

從表6可以看出，由于土層基頻總體上高于或遠高于（汶川地震波除外）表2中5條地震波的特征頻率，因此當(dāng)取ωs＝ω1時總體上看計算誤差最小.這與目前LS－DYNA 程序中約定的方法相一致.采用ωs＝ω1及ωs＝ωa形成阻尼矩陣后，El Centro波和汶川波激勵下的土層地表中點加速度時程在波形上同樣接近，區(qū)別依然體現(xiàn)在峰值上.且在土層模型3中取基頻的計算結(jié)果的誤差水平已經(jīng)和ωs＝ωa時相當(dāng).同樣，限于篇幅時程圖不再一一列出.

表7中列出了土層模型4在不同地震波激勵下的土層地表中點加速度和位移反應(yīng)的時域計算結(jié)果的誤差.

表7 土層模型4中地表加速度與位移計算誤差Tab.7 Acceleration and displacement errors of calculation of the surface midpoint on Soil Model 4 %

從表7可以看出：均勻分層土層模型4所體現(xiàn)出來的規(guī)律與模型1大體相同，并沒有因為土層的分層及坡度的改變發(fā)生太大的變化，總的來說加速度反應(yīng)誤差在10%以內(nèi)略好于均勻土層1.體現(xiàn)了該方法的實際應(yīng)用價值.

采用ωs＝ω1及ωs＝ωa形成阻尼矩陣后，El Centro波和汶川波激勵下的土層地表中點加速度時程反應(yīng)的規(guī)律與模型1相類似.

此外，從表4—7可以看出，式（8）中ωs取值對加速度反應(yīng)的影響遠大于位移反應(yīng)的影響，對于3種不同類型的土層，ωs取為土層基頻時，土層地表位移峰值的計算誤差最小，且誤差都為正.其原因是：對應(yīng)于加速度簡諧波a（t）＝A0sinωt（A0為振幅）的位移波為u（t）＝－（A0／ω2）sinωt＝v0sinωt.很顯然，隨著激勵頻率ω的升高，位移簡諧波幅值按激勵頻率平方減小.因此，土層地震位移反應(yīng)更敏感于低頻激勵，因此準(zhǔn)確計算土層第1階振型位移反應(yīng)是最為關(guān)鍵，當(dāng)ωs＝ω1時可保土層位移反應(yīng)計算的合理性.

3 結(jié)論

從本文所選三種不同土層計算模型在不同輸入地震波的地震反應(yīng)數(shù)值結(jié)果比較可看出：

（1）在土層地震反應(yīng)時域顯示算法中，所形成質(zhì)量比例阻尼矩陣對土層地震反應(yīng)有很大影響，當(dāng)土層基頻等于或高于地震波的主要特征頻率時，常規(guī)做法是取土層基頻來形成阻尼矩陣是可行的.

（2）當(dāng)土層基頻略低于地震波的主要特征頻率值時，用土層基頻來形成比例質(zhì)量矩陣時，會高估土層的地震加速度反應(yīng).本文算例中，地表加速度峰值有時高估達30%和80%不等.

（3）當(dāng)土層基頻遠低于地震波的主要特征頻率值時，采用常規(guī)方法形成質(zhì)量比例阻尼矩陣式，將更加高估土層地震加速度反應(yīng)值，本文算例中地震加速度峰值誤差有時高達60%和90%不等.

（4）由于地震波為不規(guī)則波，頻率分量豐富，很難用一個頻率來形成質(zhì)量比例阻尼矩陣以合理反映各階振型阻尼比的數(shù)值，因此應(yīng)重視在形成質(zhì)量比例阻尼矩陣中頻率的選擇.本文建議按式（8）來形成質(zhì)量比例阻尼系數(shù)a0，其中頻率ωs將針對土層和輸入地震波的特征頻率來選擇.一般來說，①當(dāng)土層基頻高于輸入地震波的特征頻率時，可取ωs為土層基頻，與目前常規(guī)認(rèn)識一致；②當(dāng)土層基頻接近于特別是遠低于地震波的特征頻率時，建議ωs取為土層基頻和輸入地震波反應(yīng)譜峰值頻率的算術(shù)平均值ωa，基本可保證計算所得土層加速度峰值略高于真實解.

（5）由于土層地震反應(yīng)計算的主要目的之一是為地面結(jié)構(gòu)提供結(jié)構(gòu)基底輸入地震加速度，為此本文上述結(jié)論主要是針對土層地震加速度反應(yīng)計算，而對于地震位移反應(yīng)計算時，ωs仍取為土層基頻為好，因為土層第1階振型的位移反應(yīng)在土層地震位移反應(yīng)中起控制作用.

［1］ 廖振鵬.近場波動問題的有限元解法［J］.地震工程與工程振動，1984，4（2）：1.LIAO Zhenpeng.A finite element method for near－field wave motion in heterogeneous materials［J］.Earthquake Engineering and Engineering Vibration，1984，4（2）：1.

［2］ Martin P P，Seed H B.One－dimensional dynamic ground response analyses ［J］.Journal of Geotechnical Engineering Division，ASCE，1982，108（7）：935.

［3］ 李小軍，廖振鵬，關(guān)慧敏.粘彈性場地地形對地震動影響分析的顯式有限元－有限差分方法［J］.地震學(xué)報，1995，17（3）：362.LI Xiaojun，LIAO Zhenpeng，GUAN Huimin.An explicit finite element－finite difference method for analyzing the effect of visco－elastic Local topograph on the earthquake motion［J］.Acta Seismologica Sinica，1995，17（3）：362.

［4］ 張汝清，殷學(xué)綱，董明.計算結(jié)構(gòu)動力學(xué)［M］.重慶：重慶大學(xué)出版社，1987.ZHANG Ruqing，YIN Xuegang，DONG Ming.Computational structural dynamics［M］.Chongqing：Chongqing University Press，1987.

［5］ 李小軍，廖振鵬，杜修力.有阻尼體系動力問題的一種顯式差分解法［J］.地震工程與工程振動，1992，12（4）：75.LI Xiaojun，LIAO Zhenpeng，DU Xiuli.An explicit finite difference method for viscoelastic dynamic problem ［J］.Earthquake Engineering and Engineering Vibration，1992，12（4）：75.

［6］ 楊柏坡，陳慶彬，王琦萍.地震波動研究中的顯式有限元法［C］／／第三屆全國地震工程會議論文集.大連：大連理工大學(xué)出版社，1990：211－216.YANG Baipo，CHEN Qingbin，WANG Qiping.An explicit finite element method of the seismic wave study［C］／／Proceedings of the Third National Conference on Earthquake Engineering.Dalian：Dalian University of Technology Press，1990：211－216.

［7］ 尚曉江，蘇建宇，王化峰.ANSYS／LS－DYNA 動力分析方法與工程實例［M］.北京：中國水利水電出版社，2008.SHANG Xiaojiang，SU Jianyu，WANG Huafeng.ANSYS／LSDYNA dynamic analysis method and projects［M］.Beijing：China Water ＆Power Press，2008.

［8］ 樓夢麟，潘旦光.滯后阻尼在土層時域分析中的應(yīng)用［J］.同濟大學(xué)學(xué)報：自然科學(xué)版，2004，32（3）：281.LOU Menglin，PAN Danguang.Hysteretic damping application in time domain analysis of soil layer［J］.Journal of Tongji University：Natural Science，2004，32（3）：281.

［9］ Clough R W，Penzien J.Dynamics of structures［M］.New York：Mc－Graw Hill Inc，1993.