錢向東,李 晨,沈人杰

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

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

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

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

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

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

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

(1)

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

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

(2)

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

(3)

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

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

(4)

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

(5)

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

(6)

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

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

(7)

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

(8)

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

(9)

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

(10)

2 動(dòng)力方程的時(shí)域求解

2.1 核函數(shù)的選擇

2.2 Newmark直接積分法

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

(11)

Newmark法計(jì)算格式為

(12)

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

(13)

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

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

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

3.1 計(jì)算模型

Koyna重力壩壩高103 m,底寬70 m，壩頂寬度為14.8 m，幾何參數(shù)如圖1所示。為了討論和比較方便，計(jì)算時(shí)略去地基對(duì)壩體動(dòng)力響應(yīng)的影響，即簡(jiǎn)單地采用剛性地基模型，同時(shí)也不考慮庫水的影響，以空庫工況進(jìn)行分析，有限元網(wǎng)格如圖2所示。壩體混凝土材料的密度為2 643 kg/m3,動(dòng)彈性模量為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)所示的廣義阻尼模型，只需一個(gè)阻尼矩陣C和一個(gè)核函數(shù)g(t)就可以表示材料的非黏滯阻尼特性。在缺乏試驗(yàn)資料的情況下，一般由質(zhì)量矩陣和剛度矩陣的線性組合計(jì)算阻尼矩陣C=a1M+b1K，系數(shù)a1、b1由系統(tǒng)的兩階振型阻尼比和自振圓頻率確定。由上述計(jì)算模型的自振特性分析，得到大壩的前2階自振圓頻率為20.20 rad/s和56.01 rad/s，從而可以確定a1=1.485、b1=1.312×10-4。

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

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

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

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

3.4 計(jì)算結(jié)果

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

圖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情況下，計(jì)算得到的壩體水平位移、速度和加速度的響應(yīng)峰值。

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

3.5 結(jié)果分析

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

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

4 結(jié) 語

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