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

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

0 引言

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

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

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

1 預備知識

定義1變分數(shù)階Coimbra微分為[18]

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

定義2變分數(shù)階Caputo微分為[19]

當f(t)=xmtn時，有

當xm=1時，有

(4)

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

(5)

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

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

用一系列移位Bernstein多項式表示矩陣

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

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

2 控制方程的建立

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

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

(11)

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

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

其中σ(x,t)是應力(Pa).粘彈性材料的變分數(shù)階模型為

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

初始條件為

(18)

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

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

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

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

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

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

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

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

其中

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

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

(23)

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

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

(24)

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

3.3 變分數(shù)階微分算子矩陣

顯然，有

其中

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

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

邊界條件可以重新寫為

初始條件可以重新寫為

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

(33)

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

4 數(shù)值分析

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

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

邊界條件

初始條件

其中

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

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

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

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

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

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

下面利用本文算法對混凝土梁進行模擬運算，混凝土材料的相關(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,變分數(shù)階函數(shù)α(t)=-0.1t+1,t∈(0,1).

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

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

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

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

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

5 結(jié)束語

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