錢向東,李 晨,沈人杰

(河海大學(xué)力學(xué)與材料學(xué)院，江蘇 南京 210098)

現(xiàn)代計算力學(xué)可以非常精確地分析線彈性結(jié)構(gòu)的靜力學(xué)問題，也能夠高精度地計算質(zhì)量矩陣和剛度矩陣，從而獲得具有較高精度的結(jié)構(gòu)無阻尼自振特性，特別是較低階的頻率和振型。然而，對于結(jié)構(gòu)振動分析中不可忽略的阻尼特性，還處于簡單的估算狀態(tài)。常用的阻尼模型為經(jīng)典的黏滯阻尼模型[1]，該模型無論是耗能機理還是試驗觀測都不符合混凝土、巖石、橡膠以及許多復(fù)合材料的阻尼特性[2]。具體應(yīng)用時，一般近似地采用振型阻尼比或比例阻尼表示阻尼特性或建立阻尼矩陣。這種結(jié)構(gòu)阻尼的表達方法忽視了其在結(jié)構(gòu)動力學(xué)中的重要性。因此，為了獲得可靠的結(jié)構(gòu)動力響應(yīng)，應(yīng)該采用符合材料耗能機理的阻尼模型，把阻尼矩陣的計算精度提升到與質(zhì)量矩陣和剛度矩陣同等的水平。

阻尼體現(xiàn)了結(jié)構(gòu)在振動過程中的能量損耗，一般認為主要由機械能轉(zhuǎn)化為熱能后消散，也有部分可能由于引發(fā)周圍介質(zhì)(地基、水、空氣等)的振動而逸散。如果不考慮能量逸散，則結(jié)構(gòu)的能量損耗主要來自材料阻尼，可由材料的本構(gòu)關(guān)系表示。事實上黏彈性、黏彈塑性材料的本構(gòu)關(guān)系均可反映材料在變形過程中的能量損耗。由于黏彈性材料本構(gòu)關(guān)系的多樣性和復(fù)雜性，各國學(xué)者提出了諸如復(fù)常數(shù)模量模型[3]、模態(tài)應(yīng)變能模型[4]、標(biāo)準(zhǔn)流變模型[5]、分數(shù)導(dǎo)數(shù)模型[6]、分數(shù)指數(shù)模型[7]、微振子模型[8]和廣義流變模型[9]等黏彈性結(jié)構(gòu)的動力分析模型，不同的動力分析模型對應(yīng)著不同的阻尼模型。其中最有代表性的就是廣義流變模型，該模型的阻尼力在數(shù)學(xué)上表示為速度與某一核函數(shù)的卷積，是一種非黏滯阻尼模型。當(dāng)核函數(shù)為狄拉克-δ函數(shù)時，模型就退化為黏滯阻尼模型。因此，該阻尼模型又被稱為廣義阻尼模型[9]，目前常用的核函數(shù)有指數(shù)型、高斯型[10]等。

基于廣義阻尼模型的結(jié)構(gòu)運動方程是耦合的二階積分-微分方程組。為了避免求解復(fù)雜的積分-微分方程組，通常在頻率域進行求解[11-12]，在頻率域構(gòu)造核函數(shù)，以便獲得可以實現(xiàn)時域逆變換的傳遞函數(shù)矩陣，應(yīng)用范圍受到了極大的限制。隨著計算技術(shù)的發(fā)展和工程應(yīng)用的需要，人們開始研究在時間域直接求解廣義阻尼結(jié)構(gòu)運動方程的方法。對于指數(shù)型核函數(shù)，通過引入狀態(tài)變量消除卷積積分,在狀態(tài)空間中采用直接積分法和振型疊加法進行求解[13-17]，段忠東等[18]采用精細積分大大提高了狀態(tài)空間法的計算精度，錢向東等[19]將該方法應(yīng)用于混凝土重力壩的地震響應(yīng)分析。對于其他類型的核函數(shù)，狀態(tài)空間法不再適用，于是有人嘗試采用中心差分法[20]、Newmark法[21-22]和Galerkin加權(quán)余量法[23]在時間域?qū)Ψe分-微分方程進行直接數(shù)值積分，相關(guān)研究都是針對簡單的多自由度算例開展，還未見具體的工程應(yīng)用。

本文根據(jù)黏彈性材料廣義流變模型的積分型本構(gòu)關(guān)系，導(dǎo)出了具有廣義阻尼模型的結(jié)構(gòu)動力學(xué)方程。以Koyna混凝土重力壩為例，將廣義阻尼模型應(yīng)用于混凝土重力壩的地震響應(yīng)分析。采用Newmark直接積分法求解積分-微分方程組，分析了松弛參數(shù)、時間積分步長對大壩地震響應(yīng)的影響，同時比較了廣義阻尼模型和黏滯阻尼模型情況下大壩的地震響應(yīng)。

1 廣義阻尼模型的結(jié)構(gòu)動力學(xué)方程

混凝土、巖石材料均具有黏彈性性質(zhì)，在線性、等溫和初始應(yīng)變ε(0)=0條件下， 其時間域內(nèi)的本構(gòu)關(guān)系可用廣義流變模型表示為[24]

(1)

在均勻、各向同性情況下，松弛函數(shù)矩陣為

G(t)=D+ηg(t)

(2)

式中：D、η分別為彈性矩陣和黏性矩陣；g(t)為卷積核函數(shù)。則式(1)可以表示為

(3)

當(dāng)η=0時，式(3)為線彈性模型的本構(gòu)關(guān)系；當(dāng)核函數(shù)g(t)為狄拉克-δ函數(shù)時，式(3)為Voigt-Kelvin模型的本構(gòu)關(guān)系。

采用有限元對黏彈性結(jié)構(gòu)進行離散，離散后結(jié)構(gòu)的動力虛功原理可以表示為

(4)

將單元的相關(guān)變量代入式(4)并考慮黏彈性本構(gòu)關(guān)系式(3)，則有

(5)

式中：ke為單元彈性剛度矩陣；me為單元質(zhì)量矩陣；ce為單元黏性阻尼矩陣；fe為單元等效結(jié)點荷載列陣。

(6)

式中：M為整體質(zhì)量矩陣；K為整體彈性剛度矩陣；C為整體黏性阻尼矩陣；F(t)為整體等效結(jié)點荷載列陣；δx為整體虛位移列陣。

由虛位移的任意性，可得黏彈性結(jié)構(gòu)的動力學(xué)方程為

(7)

式(7)為二階積分-微分方程組。記

(8)

如果結(jié)構(gòu)由n種材料組成，則阻尼力可以表示為

(9)

相應(yīng)的結(jié)構(gòu)的動力學(xué)方程為

(10)

2 動力方程的時域求解

2.1 核函數(shù)的選擇

2.2 Newmark直接積分法

對于指數(shù)型核函數(shù)，由于其函數(shù)的特殊性，通過引入一組輔助變量(狀態(tài)變量)，可以將積分-微分方程方程組(7)或(10)轉(zhuǎn)化為一階微分方程組[13]，再求解，稱之為狀態(tài)空間法。其優(yōu)點是求解過程中無須計算卷積積分，但缺點是將未知量(或方程數(shù))擴大了一倍，且不適應(yīng)用于其他類型的核函數(shù)。

(11)

Newmark法計算格式為

(12)

其中α、β為兩個計算參數(shù)，當(dāng)α≥0.5、β≥0.25(0.5+α)2時為無條件穩(wěn)定格式。將式(12)代入式(11)，并采用梯形法計算卷積積分，則可導(dǎo)出求解t+Δt時刻位移響應(yīng)的方程：

(13)

Newmark直接積分法適用于任何形式的核函數(shù)，關(guān)于算法的實現(xiàn)、卷積積分的討論以及算例驗證可參見文獻[22]。

3 混凝土重力壩地震響應(yīng)分析

本文以Koyna大壩為例，利用ANSYS有限元建模功能，采用自編的基于MATLAB的程序[22]進行大壩地震響應(yīng)分析，討論廣義阻尼模型及其核函數(shù)對重力壩地震響應(yīng)的影響。

3.1 計算模型

Koyna重力壩壩高103 m,底寬70 m，壩頂寬度為14.8 m，幾何參數(shù)如圖1所示。為了討論和比較方便，計算時略去地基對壩體動力響應(yīng)的影響，即簡單地采用剛性地基模型，同時也不考慮庫水的影響，以空庫工況進行分析，有限元網(wǎng)格如圖2所示。壩體混凝土材料的密度為2 643 kg/m3,動彈性模量為31 027 MPa,泊松比為0.15，各振型的阻尼比均為0.05。

圖1 Koyna重力壩模型尺寸Fig.1 Model size of the Koyna gravity dam

圖2 壩體有限元網(wǎng)格Fig.2 Finite element mesh of the dam

3.2 阻尼模型與參數(shù)

設(shè)大壩由一種混凝土材料組成，采用式(8)所示的廣義阻尼模型，只需一個阻尼矩陣C和一個核函數(shù)g(t)就可以表示材料的非黏滯阻尼特性。在缺乏試驗資料的情況下，一般由質(zhì)量矩陣和剛度矩陣的線性組合計算阻尼矩陣C=a1M+b1K，系數(shù)a1、b1由系統(tǒng)的兩階振型阻尼比和自振圓頻率確定。由上述計算模型的自振特性分析，得到大壩的前2階自振圓頻率為20.20 rad/s和56.01 rad/s，從而可以確定a1=1.485、b1=1.312×10-4。

阻尼模型分別采用指數(shù)核函數(shù)和高斯核函數(shù)，分別稱為指數(shù)型非黏滯阻尼和高斯型非黏滯阻尼。目前，關(guān)于松弛參數(shù)μ的試驗資料非常少，文獻[26]根據(jù)鋼筋混凝土懸臂梁的錘擊振動測試，識別得到μ的平均值約為104.3 s-1。為了考察μ的取值對阻尼效應(yīng)的影響，本文μ分別采用70 s-1、100 s-1、130 s-1、160 s-1、190 s-1進行對比計算。為了比較，同時也采用黏滯阻尼模型進行計算。

3.3 地震波與直接積分參數(shù)

采用實測的Koyna水平向地震波(圖3)作為地震輸入，順河向峰值加速度為4.647 8 m/s2，地震加速度記錄時間長度為10 s，間隔0.01 s。取Newmark直接積分法參數(shù)α=0.5、β=0.25,數(shù)值積分的時間步長Δt分別取0.01 s、0.005 s、0.002 5 s。

圖3 Koyna水平向地震加速度Fig.3 Horizontal acceleration of Koyna earthquake

3.4 計算結(jié)果

圖4和圖5分別給出了Δt=0.01 s,松弛參數(shù)μ為70 s-1、100 s-1、130 s-1、160 s-1、190 s-1情況下，指數(shù)型非黏滯阻尼和高斯型非黏滯阻尼對應(yīng)的壩頂水平向位移響應(yīng)和加速度響應(yīng)時程曲線。圖6給出了黏滯阻尼模型對應(yīng)的壩頂水平向位移響應(yīng)和加速度響應(yīng)時程曲線。

圖4 壩頂水平向地震響應(yīng)(指數(shù)型非黏滯阻尼)Fig.4 Horizontal seismic response at dam crest with exponential non-viscous damping model

圖5 壩頂水平向地震響應(yīng)(高斯型非黏滯阻尼) Fig.5 Horizontal seismic response at dam crest with Gaussian non-viscous damping model

圖6 壩頂水平向地震響應(yīng)(黏滯阻尼模型)Fig.6 Horizontal seismic response at dam crest with viscous damping model

表1則分別給出了Δt為0.01 s、0.005 s、0.002 5 s情況下，計算得到的壩體水平位移、速度和加速度的響應(yīng)峰值。

表1 壩體地震響應(yīng)峰值

3.5 結(jié)果分析

從圖4～6可以看出，各阻尼模型情況下，壩頂響應(yīng)的時間變化規(guī)律基本一致，表明阻尼模型及松弛參數(shù)μ的變化對壩頂水平位移、加速度響應(yīng)的時間變化規(guī)律影響不大。從圖4和圖5可以看出，雖然松弛參數(shù)μ的變化不影響響應(yīng)的時間變化規(guī)律，但對響應(yīng)的峰值有影響，無論是壩頂位移還是加速度，隨著松弛參數(shù)μ的不斷增大，響應(yīng)的峰值不斷減小，松弛參數(shù)越大，響應(yīng)峰值越接近黏滯阻尼的響應(yīng)峰值，與前面關(guān)于廣義阻尼模型的理論分析相符。

由表1可以看出：①各阻尼模型情況下，隨著直接積分時間步長的減小，響應(yīng)峰值不斷增大，說明Newmark直接積分法存在算法阻尼，以保證算法的無條件穩(wěn)定性。因此，為了保證計算結(jié)果的精度，在條件允許的情況下，工程分析應(yīng)盡可能采用較小的時間步長。另外，對于非黏滯阻尼模型，從表1中響應(yīng)峰值隨松弛參數(shù)μ的變化也可以看出，μ的數(shù)值越大，響應(yīng)峰值越接近黏滯阻尼的響應(yīng)峰值。②位移、速度和加速度響應(yīng)峰值對阻尼模型、松弛參數(shù)和時間積分步長的敏感性不同，加速度響應(yīng)峰值的敏感性最大，其次是速度響應(yīng)峰值，位移響應(yīng)峰值的敏感性最小，符合函數(shù)及其導(dǎo)數(shù)的變化規(guī)律。③比較高斯型非黏滯阻尼的結(jié)果，松弛參數(shù)μ從70 s-1變化至190 s-1，3種時間步長下的響應(yīng)峰值始終大于黏滯阻尼的響應(yīng)峰值，呈現(xiàn)從右側(cè)一致逼近黏滯阻尼結(jié)果的趨勢。當(dāng)時間步長Δt≥0.005 s時，指數(shù)型非黏滯阻尼模型的響應(yīng)峰值并不隨著松弛參數(shù)μ的增大從右側(cè)一致逼近黏滯阻尼結(jié)果的趨勢，不符合核函數(shù)一致收斂于δ函數(shù)的趨勢。當(dāng)Δt=0.002 5 s時，指數(shù)型非黏滯阻尼模型的響應(yīng)峰值才表現(xiàn)出隨μ的增大從右側(cè)一致逼近黏滯阻尼結(jié)果的趨勢。由此說明，當(dāng)采用Newmark直接積分法時，指數(shù)型非黏滯阻尼模型對時間步長的要求高于高斯型非黏滯阻尼模型。

4 結(jié) 語

本文根據(jù)黏彈性材料的廣義流變模型導(dǎo)出了廣義阻尼結(jié)構(gòu)動力學(xué)方程，指出了經(jīng)典的黏滯阻尼模型可以從Voigt-Kelvin本構(gòu)模型導(dǎo)出，是廣義阻尼模型的一種特例。以Koyna混凝土重力壩為例，進行了基于廣義阻尼模型的地震響應(yīng)分析。分析時分別采用指數(shù)核函數(shù)和高斯核函數(shù)2種廣義阻尼模型，并與黏滯阻尼模型的結(jié)果進行了比較，表明地震作用下廣義阻尼模型的響應(yīng)峰值均大于黏滯阻尼模型的響應(yīng)峰值。因此，有必要采用非黏滯阻尼模型進行混凝土大壩的地震響應(yīng)分析。采用的Newmark直接積分法對卷積核函數(shù)的類型沒有限制，是一種通用的時域求解方法，與狀態(tài)空間法只適用于指數(shù)核函數(shù)相比，更便于嵌入現(xiàn)有的結(jié)構(gòu)動力分析系統(tǒng)中。從本文算例的結(jié)果來看，指數(shù)型非黏滯阻尼模型對數(shù)值積分時間步長的要求高于高斯型非黏滯阻尼模型。本文對卷積積分采用了梯形積分格式，沒有討論其他積分格式對計算精度和耗時的影響，存在進一步改進和優(yōu)化的空間。非黏滯阻尼模型中松弛參數(shù)μ的取值對計算結(jié)果有較大的影響，為了更加科學(xué)合理地分析混凝土結(jié)構(gòu)的動力響應(yīng)，應(yīng)進一步開展相關(guān)試驗和動力參數(shù)識別的研究，以確定混凝土材料的松弛參數(shù)。