金素花，解加全，韓存弟，孫虹霞，陳一鳴

(1.燕山大學(xué) 理學(xué)院，河北 秦皇島 066004；2.太原師范學(xué)院 數(shù)學(xué)系，山西 晉中 030619；3.太原理工大學(xué) 先進(jìn)成形與智能裝備研究院，山西 太原 030024 )

0 引言

粘彈性材料是介于理想粘性和理想彈性之間的一種材料，由于其具有良好的吸噪和減振性能[1-3]，所以在工程實(shí)踐中得到了廣泛應(yīng)用.混凝土是生活中常見的一種粘彈性材料，經(jīng)常被用于道路工程.隨著粘彈性材料的逐漸發(fā)展，變分?jǐn)?shù)階模型得到了廣泛的研究.相較于整數(shù)階和分?jǐn)?shù)階模型，變分階模型可以利用較少的參數(shù)去描述動(dòng)態(tài)粘彈性行為.在大變形條件下，變分?jǐn)?shù)階微分算子比分?jǐn)?shù)階微分算子能更準(zhǔn)確地模擬粘彈性材料的本構(gòu)關(guān)系[4].為了準(zhǔn)確描述非晶玻璃態(tài)聚合物的壓縮變形，Meng等[5]提出了變分?jǐn)?shù)階本構(gòu)模型，Li等[6]提出了一種變分?jǐn)?shù)階模型來描述形狀記憶聚合物的粘彈性特性，并與分?jǐn)?shù)階模型進(jìn)行比較，證明了變分?jǐn)?shù)階微分方程模型比分?jǐn)?shù)階微分方程模型更適合于形狀記憶聚合物的記憶行為建模.

變分?jǐn)?shù)階模型的研究不僅包含模型的建立，還包含模型的計(jì)算.變分?jǐn)?shù)階粘彈性模型的研究主要有Adomian方法[7]、同倫攝動(dòng)法[8]、同倫分析法[9]、變分迭代法[10]和小波函數(shù)算子矩陣分析法[11,12]等.由于變分?jǐn)?shù)階偏微分方程的階數(shù)會(huì)隨著自變量時(shí)刻發(fā)生變化，這種不穩(wěn)定性極大增加了求解方程的難度，目前已有的計(jì)算方法都存在計(jì)算步驟復(fù)雜和計(jì)算精度低等缺點(diǎn).Bernstein多項(xiàng)式具有良好的穩(wěn)定性和逼近性.Kadkhoda[13]利用Bernstein多項(xiàng)式求解了流體力學(xué)中線性變分?jǐn)?shù)階微分方程，Wang等[14]采用Bernstein多項(xiàng)式研究了兩類時(shí)間分?jǐn)?shù)階廣義五階微分方程的數(shù)值解，Derakhshan等[15]提出了一種利用Bernstein多項(xiàng)式求解Caputo-Prabhakar意義下一類非線性變分?jǐn)?shù)階微分方程的數(shù)值方法.

兩端固定梁作為一個(gè)經(jīng)典的力學(xué)模型在工程中有著廣泛的應(yīng)用.近年來，梁的分?jǐn)?shù)階振動(dòng)被廣泛研究[16,17]，但是對于梁的變分?jǐn)?shù)階振動(dòng)的研究幾乎沒有，文中建立了粘彈性梁的變分?jǐn)?shù)階位移控制方程，用移位Bernstein多項(xiàng)式求解變分?jǐn)?shù)階粘彈性微分方程的數(shù)值解，并結(jié)合數(shù)值模擬分析混凝土粘彈性梁的動(dòng)力學(xué)行為.

1 預(yù)備知識

定義1變分?jǐn)?shù)階Coimbra微分為[18]

其中,t≥0,0<α(x) ≤1,α(x)是變分?jǐn)?shù)階函數(shù)，f(x,t)是連續(xù)函數(shù)，Γ是Gamma函數(shù).

定義2變分?jǐn)?shù)階Caputo微分為[19]

當(dāng)f(t)=xmtn時(shí)，有

當(dāng)xm=1時(shí)，有

(4)

定義3在區(qū)間[0,1]上定義Bernstein多項(xiàng)式

(5)

擴(kuò)展x到區(qū)間[0,L]上，移位Bernstein多項(xiàng)式的通項(xiàng)公式為

用二項(xiàng)式定理展開(L-x)k-i，有

用一系列移位Bernstein多項(xiàng)式表示矩陣

其中，k為多項(xiàng)式的項(xiàng)數(shù)，且

P是可逆的，且Hk(x)=P-1φ(x).

2 控制方程的建立

考慮一個(gè)由粘彈性材料做成的梁，受到一個(gè)垂直于梁水平方向的外部載荷，使得梁發(fā)生形變，如圖1所示.由于粘彈性材料的影響，其形變會(huì)發(fā)生一定程度上的恢復(fù)行為，因此，產(chǎn)生了垂直于梁水平方向的振動(dòng).兩端固定梁的動(dòng)力學(xué)方程為

圖1 兩端橫向固定梁的模型

(11)

其中ρ是粘彈性材料的密度(kg/m3)，S是梁的橫截面積(m2)，w(x,t)是梁沿外載荷方向的位移(m)，M(x,t)是彎矩(N·m)，f(x,t)是外部載荷(N).

梁彎矩和應(yīng)力間的關(guān)系為

其中σ(x,t)是應(yīng)力(Pa).粘彈性材料的變分?jǐn)?shù)階模型為

將(12)～(14)式代入(11)式，得到兩端固定梁的控制方程為

初始條件為

(18)

3 移位Bernstein多項(xiàng)式數(shù)值算法

應(yīng)用移位Bernstein 多項(xiàng)式對梁位移控制方程在時(shí)域內(nèi)進(jìn)行數(shù)值求解，算法結(jié)構(gòu)如圖2所示.

圖2 移位Bernstein多項(xiàng)式的數(shù)值算法結(jié)構(gòu)

3.1 基函數(shù)逼近位移函數(shù)

位移函數(shù)w(x,t)∈L2([0,L]×[0,R])用移位的Bernstein多項(xiàng)式去逼近，即

其中ui,j為二維函數(shù)，U為函數(shù)逼近系數(shù).

3.2 整數(shù)階微分算子矩陣

定義4如果存在矩陣D,使得φ′(x)=Dφ(x)，則稱D為關(guān)于x的一階微分算子矩陣，即

其中

將(20)式左右兩邊再求一次導(dǎo)數(shù)，有

φ″(x)=D2φ(x).

(23)

根據(jù)(20)和(23)式，應(yīng)用數(shù)學(xué)歸納法可得

φ(k)(x)=Dkφ(x),

(24)

根據(jù)(19)和(23)式可以得到

3.3 變分?jǐn)?shù)階微分算子矩陣

顯然，有

其中

結(jié)合(19),(24)和(26)式，有

將(25)和(29)式代入(15)式，梁位移控制方程轉(zhuǎn)化成矩陣乘積的形式

邊界條件可以重新寫為

初始條件可以重新寫為

φT(x)Uφ(0)=φT(x)UDφ(0)=0.

(33)

利用配點(diǎn)法，將(30)式的變量進(jìn)行離散得到一組代數(shù)方程組.應(yīng)用Matlab和最小二乘法求得ui,j，代入原來方程即可得到兩端固定梁的位移數(shù)值解.

4 數(shù)值分析

4.1 變分?jǐn)?shù)階微分方程數(shù)值解的比較

本節(jié)求解粘彈性梁位移控制方程算例的數(shù)值解，并對精確解與數(shù)值解進(jìn)行比較.控制方程對應(yīng)的一般算例為

邊界條件

初始條件

其中

變分?jǐn)?shù)階函數(shù)α(t)=-0.1t+1,t∈(0,1).數(shù)值算例的精確解為w(x,t)=x2(1-x2)t2,x∈(0,1),t∈(0,1).

采用移位Bernstein多項(xiàng)式置配算法時(shí)，多項(xiàng)式的項(xiàng)數(shù)選擇為n=4；用該算法求解變分?jǐn)?shù)階偏微分方程的算例時(shí)，數(shù)值解表示為wn(x,t).

圖3為一般數(shù)值算例的數(shù)值解和精確解在x=0.5處的比較.由圖可知，數(shù)值解和精確解的圖像高度一致，說明文中算法對求解變分階微分方程具有可行性.圖4顯示的是數(shù)值解和精確解的絕對誤差e(x,t)=|wn(x,t)-w(x,t)|，絕對誤差數(shù)量級可以達(dá)到10-11，進(jìn)一步證明本文算法對解決變分?jǐn)?shù)階粘彈性梁問題具有足夠高的準(zhǔn)確性.

圖3 算例的數(shù)值解與精確解

圖4 數(shù)值解的絕對誤差

4.2 不同載荷下的粘彈性梁分析

下面利用本文算法對混凝土梁進(jìn)行模擬運(yùn)算，混凝土材料的相關(guān)參數(shù)如下[20,21]：E=1 000 MPa,θ=0.18 s,ρ=4 000 kg·m-3.選取梁的長度l=5 m,橫截面積S=0.01 m2,慣性矩I=0.14/12,變分?jǐn)?shù)階函數(shù)α(t)=-0.1t+1,t∈(0,1).

將數(shù)據(jù)代入梁的控制方程(15)，利用本文算法對粘彈性梁在不同載荷下進(jìn)行數(shù)值求解.在不同均布載荷作用下，混凝土梁在不同位置和時(shí)間的位移數(shù)值解如圖5所示.粘彈性梁的兩端位移始終為零，不受時(shí)間的影響，這與邊界條件是一致的.在其它位置，粘彈性梁的位移隨時(shí)間逐漸增大，在t=1 s時(shí)達(dá)到最大值.粘彈性梁的最大位移在梁中部x=2.5 m處.梁的位移關(guān)于中間位置對稱.

比較圖5(a)和(b)可知，當(dāng)粘彈性梁上施加的均布載荷不同時(shí)，其位移也不同.在f(x,t)=10的均布載荷下，位移最大值約為0.14 m，在f(x,t)=20的均布載荷下，位移最大值為0.29 m，說明粘彈性梁的位移隨所施加的均勻載荷的增加而增加，這一結(jié)論與文獻(xiàn)[22]的結(jié)論相一致.與分?jǐn)?shù)階方程相比，變分?jǐn)?shù)階方程能夠較準(zhǔn)確地模擬大應(yīng)變下粘彈性梁的動(dòng)態(tài)行為,因而具有更廣泛的應(yīng)用.

圖5 不同均布載荷下的位移數(shù)值解

粘彈性梁在線性載荷作用下梁的位移隨線性載荷斜率的增加而增大.當(dāng)斜率為2(圖6(a))時(shí)，梁的最大位移約為0.09 m；當(dāng)斜率為4(圖6(b))時(shí)，梁的最大位移約為0.19 m.

圖6 不同線性載荷下的位移數(shù)值解

5 結(jié)束語

文中基于移位Bernstein多項(xiàng)式對變分?jǐn)?shù)階粘彈性梁提出了一種有效的數(shù)值算法.基于兩端固定梁的動(dòng)力學(xué)方程、粘彈性材料變分?jǐn)?shù)階本構(gòu)關(guān)系和應(yīng)變位移關(guān)系建立了變分?jǐn)?shù)階粘彈性梁的控制方程，數(shù)值分析結(jié)果表明，該算法求解變分?jǐn)?shù)階偏微分方程具有高效性;不同載荷下，混凝土梁的位移隨時(shí)間的增大而增加,梁的位移在梁中間位置處達(dá)到最大值且關(guān)于其對稱;混凝土梁的位移會(huì)隨著均布載荷的增加而增加，同時(shí)也會(huì)隨著線性載荷斜率的增加而增大.