分數階微分雙參數黏彈性地基矩形板動力響應

寇 磊， 白 云

(同濟大學 地下與建筑工程系，上海 200092)

分數階微分雙參數黏彈性地基矩形板動力響應

寇 磊， 白 云

(同濟大學 地下與建筑工程系，上海 200092)

基于彈性地基Pasternak雙參數模型，利用分數階微分得到黏彈性地基雙參數模型，并在此基礎上建立采用分數階微分Kelvin模型的雙參數黏彈性地基上彈性和黏彈性矩形板在動荷載作用下的動力方程；利用Galerkin方法和分段處理的數值計算方法求解四邊簡支的彈性和黏彈性地基板的動力方程，通過自由振動算例驗證該求解方法的正確性；并分析沖擊動荷載作用下分數階微分Kelvin模型的分數階、粘滯系數、水平剪切系數和模量參數對位移響應的影響。結果表明：分數階微分黏彈性模型可以描述不同黏彈性材料的力學行為；分數階取值0.5前后，矩形板位移響應值出現了不同的衰減發展形態；粘滯系數、水平剪切系數和模量系數取值越大，位移響應衰減速度越快。

分數階微分；黏彈性地基；雙參數模型；動力響應；參數影響

parametric influence

動荷載作用下地基上板的動力響應是剛性路面板、機場跑道、船塢底板和碼頭平臺等實際工程中常見的一類問題，因而得到了廣泛重視和大量研究[1-6]。目前，在研究地基與上部結構相互作用的理論中，常采用的地基計算模型主要有Winkler模型、彈性半空間模型以及雙參數模型，雙參數模型即克服了Winkler模型中地基變形不連續的缺陷，又在數學計算上比半空間模型易于處理[5-7]；絕大多數研究在描述土體變形隨時間變化的流變性質時采用流變經驗公式、流變元件模型和經典黏彈性模型，而實際天然地基土體的力學特性受眾多因素影響，其應力與應變呈現明顯的非線性及隨時間變化的特性；鑒于分數階微積分非常適合于刻畫具有記憶和遺傳性質的材料和過程[8]，采用分數微積分來描述軟土的流變本構關系逐漸得到廣泛應用[9-13]。

本文基于分數階微積分，建立雙參數黏彈性地基模型，在此基礎上對地基上薄板在動荷載作用下進行分析，并對影響動力響應的各個參數進行討論，為實際工程中地基與上部結構的動力相互作用分析提供理論基礎。

1 分數階微分雙參數黏彈性地基矩形板受荷動力方程

1.1 分數階微分雙參數黏彈性地基

采用分數階微積分理論中的Riemannn-Liouville 分數階算子[14]，對于函數f(t)的α階微分定義為

(1)

設f(t)至少存在一階導數連續，并令t0=0，式(1)展開為

(2)

由式(2)可知，分數階微分作用于函數f(t)等于Abel核Ir(t)與函數f(t)的廣義Stieltjes卷積。

對于黏彈性材料，其積分型本構方程為

(3)

式中，J(t)為蠕變柔量；σ(t)為應力；ξ(t)為應變。

由式(2)和式(3)可知, R-L分數階微分與描述黏彈性材料的積分型本構方程表達形式是一致的。

天然地基在荷載作用下的沉陷變形具有明顯的時間效應，地基土體介于理想固體和理想流體之間，因而考慮采用分數階微分描述地基土體的黏彈性特征更符合實際狀況。根據Pasternak 提出的雙參數彈性地基模型[15]，假定黏彈性地基土體在水平剪切層(x-y面)為均質各向同性，并同時考慮地基在傳遞上部結構變形時隨時間變化的影響，如圖1所示，得到剪切層每單位寬度上的剪切力為:

圖1 黏彈性地基雙參數模型Fig.1 The two-parameter viscoelastic foundation

(4)

式中，Nx、Ny為剪切力；τxz、τyz為剪應力；Dα為R-L分數階微分算子；μ為水平剪切系數，與剪切層變形有關；wd(x,y,t)為地基的沉陷。

根據Padovan得到的分數階微分黏彈性本構方程[16]，即

Sσ(t)=Qξ(t)

(5)

得到黏彈性地基土體豎向彈簧層的作用力為

R=S-1Qwd(x,y,t)

(6)

同時建立z方向力的平衡方程，即

(7)

將式(4)和式(6)代入方程(7)中，得到黏彈性地基與上部結構之間的相互作用力q(x,y,t)，即:

q(x,y,t)=S-1Qwd(x,y,t)-μDα2wd(x,y,t)

(8)

若令式(5)中M=N=1,b1=0則式(5)簡化為:

σ(t)=E0ξ(t)+E1Dαξ(t)=kξ(t)+ηDαξ(t)

(9)

式(9)稱為分數階微分Kelvin模型，該模型由模量系數為k的Hook彈簧和粘滯系數為η含分數階α的Newton粘壺并聯而成，如圖2所示。當α=1時，即為經典Kelvin模型。

圖2 分數階微分Kelvin模型Fig.2 Fractional derivative Kelvin model

將式(9)代入式(8)得到采用分數階微分Kelvin模型的q(x,y,t)的表達式，即

q(x,y,t)=kwd(x,y,t)+

ηDαwd(x,y,t)-μDα2wd(x,y,t)

(10)

1.2 雙參數黏彈性地基上黏彈性薄板受荷動力方程

對于黏彈性地基上的均質各向同性黏彈性薄板，設其體積變形為彈性變形，即

σkk(t)=3Kξkk(t)

設其純剪切變形符合分數階微分Kelvin模型，即

Sij(t)=2(G+η′Dβ)eij

式中，K為體積模量；G為剪切模量；η′為黏性系數；Dβ為R-L分數階微分算子。

根據薄板ξz=0得到

(11)

基于kirchhoff薄板理論，可得

(12)

當地基板上受動荷載p(x,y,t)作用，薄板的動力控制微分方程為

(13)

式中，w(x,y,t)為板的變形；q(x,y,t)為地基與板之間的相互作用力；ρ，h為板的單位密度和厚度。

將式(12)和式(10)代入方程(13)，根據變形協調條件wd=w，得到采用分數階微分Kelvin模型的雙參數黏彈性地基上黏彈薄板的動力方程，即

(14)

1.3 雙參數黏彈性地基上彈性薄板動力方程

對于黏彈性地基上的均質各向同性彈性薄板，根據kirchhoff薄板理論，得到

(15)

將式(15)和式(10)代入方程(13)，得到采用分數階微分Kelvin模型的雙參數粘彈性地基上彈性薄板的動力方程，即

D4w+kw+ηDαw-μDα2w+

(16)

2 動力方程求解

對于板四邊簡支情況，如圖3，邊界條件為

初始條件為

圖3 四邊簡支板Fig.3 Plate with four edges simply supported

由于非線性偏微分-積分方程(14)、(16)難以求解，采用Galerkin方法進行簡化，在邊界條件下，方程(14)、(16)的解可展開為三角級數形式

(17)

將式(17)代入式(14)，得到：

AT″(t)+BT(t)+CDβT(t)+EDαT(t)=Ff(t)

(18)

式中，

將式(18)變化為

(19)

先計算方程組(19)出現的分數階微分項DβT(t)、DαT(t)，現以DαT(t)為例。

根據R-L分數階微分的定義可知，設T(t)至少存在一階導數連續，則

(20)

(21)

被積函數在S1的定義域中處處有定義，由復化梯形求積公式直接求得[17]。

n≥3

(22)

(23)

由式(20)、式(22)和式(23)求出分數階微分項DαT(t)，同理求出分數階微分項DβT(t)。對于方程組(19)，將計算時間t進行離散t=nΔt，Δt為等距步長；先用Euler方法計算起始3步的值，然后用改進的四步Adams預估-校正算法進行計算。

方程(16)同樣按上述方法進行求解。

3 算例及參數討論

3.1 求解自由振動

采用R-L分數階微分算子的黏彈性振子的自由振動，其微分方程及初始條件為[19]

(24)

式中，f(t)為位移函數；a，b為常數。當α=1/2時，方程(24)的解析解為[17]

采用本文提出的分數階微分數值算法Matlab編程求解方程(24)，計算中選取a=0.2，b=1。并將計算結果與該方程的解析解進行對比，如圖4所示。

圖4 分數階自由振動的解析解與數值解對比Fig.4 The comparison between the analytical and numerical solution of fractional free vibration

從圖4可知,本文的數值計算結果與精確解相當吻合,從而驗證了本文求解分數階微分方程數值方法的正確性。

3.2 求解沖擊動載作用

放置在黏彈性地基上的四邊簡支彈性矩形薄板，板中心處受到沖擊動載作用，沖擊荷載采用半波余弦曲線表示，其荷載表達式為

p(x,y,t)=

式中,p0為加載峰值，T為荷載沖擊周期，x0，y0為板中心坐標。

利用Matlab編程計算t=0～50 s的地基板中心的位移響應及平行x軸的板中心軸線在t=25 s、t=50 s的位移響應值。計算時選取p0=100 kN,T=0.1 s，板和地基相應參數如表1。由于一階和高階Galerkin截斷得到的動力學行為定性相同[20]，為簡化計算，令m=1，n=1。計算結果如圖5、圖6所示。

表1 矩形板及地基參數

圖5 t=25s、t=50s時板中心軸線的位移響應曲線Fig．5Thedisplacementcorrespondingcurvesofcentrallineoftheplateatt=25sandt=50s圖6 t=0-50s時板中心的位移響應曲線Fig．6Thedisplacementcorrespondingcurvesoftheplatecenteratt=0-50s圖7 不同α取值時板中心位移響應時間曲線Fig．7Thedisplacementcorrespondingcurvesattheplatecenterunderdifferentαvalue

由圖5可知，地基薄板在沖擊荷載作用下沿平行x軸的中心軸線的位移響應以板中心位置兩端對稱分布，位移響應值從端部向板中心區域發展，板中心位置達到最大值；隨著時間的發展，位移響應值衰減明顯。由圖6可知，板中心的位移響應值隨著沖擊荷載作用迅速增大，當地基板不在受沖擊荷載作用即t>0.1 s后，位移值繼續增大，在t=7 s左右達到最大值，隨后位移振幅衰減明顯，但衰減周期保持不變。

3.3 參數分析

在沖擊荷載算例的基礎上，以下將逐一研究分數階α、粘滯系數η、水平剪切系數μ和模量系數k對地基板動力位移響應的影響，在討論某一個參數時，其它參數值保持不變。

將分數階α分別取為0.1、0.3、0.5、0.7，0.9，計算t=0～50 s時薄板中心處的位移響應值及t=50 s時薄板平行x軸的板中心軸線的位移響應值，計算結果如圖7、圖8所示。

由圖7可知，分數階α值越大，板受到的阻尼越小，板中心位移響應達到初始峰值的時間越長，初始峰值也越大；隨著分數階α值的不斷增大，由于板受到的阻尼不斷變小，位移振幅衰減周期越來越長；當α<0.5時位移響應振幅的衰減明顯，當0.5<α<1位移響應振幅的衰減逐漸放緩。

圖8 t=50 s時板中心軸線的位移響應曲線圖Fig.8 The displacement corresponding curves of central line of the plate of t=50 s under different α values

由圖8可知，分數階α的取值對板的位移響應發展形態有較大影響；分數階α取值的不同，板中心軸線的位移響應均以板中心位置兩端對稱分布，板中心為絕對值最大位置，但位移響應的性質卻發生了變化：當α<0.5時位移響應為正值，α值越大，同一位置處位移響應越小；當0.5<α<1時，位移響應為負值，α值越小，同一位置處位移響應越小；其中，當分數階α→1時逐漸退化為采用經典Kelvin模型的薄板位移響應情況。由以上分析可知，通過改變分數階α值的大小可以反映不同黏彈性材料的力學行為，分數階α值可以根據材料的流變實驗數據擬合得到，采用分數階黏彈性模型比經典黏彈性模型更能很好的反映黏彈性材料的力學性質。

將粘滯系數η分別取為1 MN/m3、10 MN/m3、45 MN/m3、100 MN/m3，計算t=0～50 s時薄板中心處的位移響應值，計算結果如圖9所示。

由圖9可知，粘滯系數η值越大，板受到的阻尼越大，板中心位移響應達到初始峰值的時間越短，初始峰值也越小；隨著粘滯系數η值的不斷增大，位移響應振幅的衰減越來越劇烈，振幅衰減速度也越快，衰減周期越來越小。

將水平剪切系數μ分別取為0、100 MN/m3、250 MN/m3、500 MN/m3，計算t=0～50 s時薄板中心處的位移響應值，計算結果如圖10所示。

由圖10可知，當μ=0即為Winkler模型的位移響應值；水平剪切系數μ值越大，板受到的阻尼越大，板中心位移響應達到初始峰值的時間越短，初始峰值也越小；隨著水平剪切系數μ的不斷增大，位移響應振幅的衰減越來越劇烈，振幅衰減速度也越快，衰減周期越來越小。

將模量系數k分別取為5 MN/m3、20 MN/m3、50 MN/m3、100 MN/m3，計算t=0～50 s時薄板中心處的位移響應值，計算結果如圖11所示。

圖9 不同η取值時板中心位移響應時間曲線Fig．9Thedisplacementcorrespondingcurvesattheplatecenterunderdifferentηvalue圖10 不同μ取值時板中心位移響應時間曲線Fig．10Thedisplacementcorrespondingcurvesattheplatecenterunderdifferentμvalue圖11 不同k取值時板中心位移響應時間曲線Fig．11Thedisplacementcorrespondingcurvesattheplatecenterunderdifferentkvalue

由圖11可知，模量系數k值越大，板中心位移響應達到初始峰值的時間越短，初始峰值也越小；隨著模量系數k的不斷增大，位移響應振幅的衰減越來越明顯，振幅衰減速度也變快，衰減周期越小。

4 結 論

基于分數階微分理論、黏彈性理論和上部結構與地基相互作用理論，研究了分數階微分雙參數黏彈性地基矩形板的動力響應問題， 并討論了分數階、粘滯系數、水平剪切系數和模量參數對位移響應的影響，得到的主要結論有：

(1) 提出的求解分數階微分雙參數黏彈性地基上薄板動力方程的數值方法，可以推廣到求解黏彈性地基上其它結構形式的動力響應問題；

(2) 分數階微分黏彈性模型可以較大范圍描述不同黏彈性材料的力學行為；

(3) 薄板動荷載作用下位移響應值隨著時間的發展出現衰減現象；分數階α取值在α=0.5前后，薄板位移響應值出現了不同的衰減發展形態；粘滯系數η、水平剪切系數μ和模量系數k取值越大，薄板位移響應衰減速度越快，衰減周期越小，其中，粘滯系數η、水平剪切系數μ起主導作用；

(4) 在實際應用中，需要根據不同地基土類型的力學特性，選取合理的分數階微分黏彈性模型機確定其相應參數的取值范圍。

[ 1 ] Zaman M, Taheri M R, Alvappillal A. Dynamic response of a thick plate on viscoelastic foundation to moving loads[J]．International Journal for Numerical Analysis and Method in Geo-mechanics,199l,15(9):627-647.

[ 2 ] 祝彥知,薛保亮,王廣國.粘彈性地基上粘彈性地基板的自由振動解析[J].巖石力學與工程學報,2002,21(1):112-118. ZHU Yan-zhi, XUE Bao-liang, WANG Guang-guo. Free vibration analysis of viscoelastic foundation plate on viscoelastic foundation[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(1):112-118.

[ 3 ] Sun L. Dynamic response of kirchhoff plate on a viscoelastic foundation to harmonic circular loads[J].Journal of Applied Mechanics, 2003, 70(4):595-600.

[ 4 ] Kim S M, McCullough B F. Dynamic response of plate on viscous Winkler foundation to moving loads of varying amplitude [J].Engineering Structures, 2003, 25(9):1179-1188.

[ 5 ] 張選兵，羅先啟，葛修潤，等. 雙參數彈性地基上板承受沖擊荷載的動力響應的解析解[J].巖石力學與工程學報,2001, 20(6),846-850. ZHANG Xuan-bing, LUO Xian-qi, GE Xiu-run, et al. Analytical solution on dynamic response of rectangular plate with free edges of two parameters elastic foundation under impact load [J].Chinese Journal of Rock Mechanics and Engineering,2001(11), 20(6):846-850.

[ 6 ] 顏可珍,夏唐代,黃立葵.雙參數粘彈性地基無限長板的瞬態動力響應分析[J].巖石力學與工程學報, 2005,24(24):4576-4580. YAN Ke-zhen, XIA Tang-dai, HUANG Li-kui. Dynamic response of strip on two-parameter viscoelastic foundation under impact loading [J].Chinese Journal of Rock Mechanics and Engineering, 2005, 24(24):4576-4580.

[ 7 ] 彭麗，丁虎，陳立群.黏彈性Pasternak地基梁振動的復模態分析[J].振動與沖擊,2013,32(2):143-146. PENG Li, DING Hu, CHEN Li-qun. Complex modal analysis for vibrations of a beam on a viscoelastic Pasternak foundation[J].Journal of Vibration and Shock 2013,32(2):143-146.

[ 8 ] 王東生，曹磊. 混沌、分形及其應用[M].合肥：中國科學技術大學出版社,1995．

[ 9 ] 孫海忠，張衛. 一種分析軟土黏彈性的分數導數開爾文模型[J].巖土力學,2007,28(9)：1983-1986. SUN Hai-zhong, ZHANG Wei. Analysis of soft soil with viscoelastic fractional derivative Kelvin model [J]. Rock and Soil Mechanics, 2007, 28(9):1983-1986.

[10] 殷德順，任俊娟，和成亮,等.一種新的巖土流變模型元件[J].巖石力學與工程學報,2007,26(9):1899-1903. YIN De-shun, REN Jun-juan, HE Cheng-liang, et al. A new rheological model element for geo-materials[J].Chinese Journal of Rock Mechanics and Engineering, 2007,26(9):1899-1903.

[11] 劉林超，楊驍.，豎向集中力作用下分數導數型半無限體粘彈性地基變形分析[J].工程力學,2009,26(1):13-17. LIU Lin-chao, YANG Xiao. Analysis on settlement of semi-infinite viscoelastic ground based on fractional model[J]. Engineering Mechanics, 2009, 26(1):13-17.

[12] 何利軍，孔令偉，吳文軍，等.采用分數階導數描述軟黏土蠕變的模型[J].巖土力學,2011,32(增刊2):239-243. HE Li-jun, KONG Ling-wei, WU Wen-jun, et al. A description of creep model for soft soil with fractional derivative [J]. Rock and Soil Mechanics, 2011, 32(S2):239-243.

[13] 寇磊.分數階微分型雙參數黏彈性地基矩形板受荷響應[J].力學季刊,2013,34(1):154-160. KOU Lei. Response of rectangular plate on fractional derivative two-parameter viscoelastic foundation[J].Chinese Journal of Mechanics,2013,34(1):154-160.

[14] Samko S G, Kilbas A A, Marichev O L. Fractional integrals and derivatives: theory and application[M].USA: Gordon and Breach Science Publishers, 1993.

[15] Arnold D K. Elastic and viscoelastic foundation models [J].Journal of Applied Mechanics, 1964, 31(3):491-498.

[16] Padovan J. Computational algorithm for FE formulations involving fractional operators [J].Computational Mechanics, 1987, 2(4):271-287.

[17] 張衛，徐華，清水信行.，分數算子描述的黏彈性體力學問題數值方法[J].力學學報,2004,36(5)：617-622. ZHANG Wei, XU Hua, Shimizut N. The numerical analysis formulation of the viscoelastic solid modeled by fractional operator[J]. Chinese Journal of Theoretical and Applied Mechanics, 2004, 36(5):617-622.

[18] 朱正佑，李根國，程昌鈞. 分數積分的一種數值計算方法及其應用[J].應用數學和力學,2003,24(4):331-341. ZHU Zheng-you, LI Gen-guo, CHENG Chang-jun. A numerical method for fractional integral with applications[J]. Applied Mathematics and Mechanics, 2003, 24(4):331-341.

[19] Sakakibara S. Properties of vibration with fractional derivative damping of order 1/2[C] //JSME International Journal, Series C: Dynamics, Control, Robotics, Design and Manufacturing, Japan, Iwaki-Meisei University, 1997:393-399.

[20] 陳立群，程昌鈞. 非線性粘彈性梁的動力學行為[J].應用數學和力學，2000,21(9):897-902. CHEN Li-qun, CHENG Chang-jun. Dynamical behavior of nonlinear viscoelastic beams[J].Applied Mathematics and Mechanics,2000,21(9):897-902.

Dynamic response of rectangular plates on two-parameter viscoelastic foundation with fractional derivatives

KOU Lei, BAI Yun

(Department of Geotechnical Engineering, Tongji University, Shanghai 200092, China)

Based on a two-parameter Pasternak model of elastic foundation, a two-parameter model of viscoelastic foundation was derived by using fractional derivatives. The dynamic equations of elastic and viscoelastic rectangular plates on a two-parameter viscoelastic foundation with the fractional Kelvin model under dynamic loads were established. The dynamic equations of elastic and viscoelastic rectangular plates with four-edge simply supported were solved with Galerkin method and the segmented numerical method, the correctness of the solutions was verified with examples of free vibration. The influences of fractional order, viscosity parameter, horizontal shear coefficient and modulus of the fractional Kelvin model under an impact load on the displacement responses of the plates were analyzed. The results showed that the fractional derivative viscoelastic model can describe the mechanical behavior of different viscoelastic materials; the displacement responses of the rectangular plates have different attenuation forms before and after the fractional order value of 0.5; the attenuation speed of the displacement response increases with increase in viscosity coefficient, horizontal shear coefficient and modulus.

fractional derivative; viscoelastic foundation；two-parameter model; dynamic response；

教育部長江學者和創新團隊發展計劃 (IRT1029)

2013-03-25 修改稿收到日期：2013-05-30

寇磊 男，博士生，1983年5月生

TU44

A

10.13465/j.cnki.jvs.2014.08.025