微分求積法在工程結(jié)構(gòu)動力學(xué)中的應(yīng)用研究

王冬梅，張 偉，劉寅立

(1. 天津科技大學(xué)理學(xué)院，天津 300457；2. 北京工業(yè)大學(xué)機(jī)電學(xué)院，北京 100124)

隨著科學(xué)技術(shù)的發(fā)展，現(xiàn)代工程結(jié)構(gòu)，無論是航空航天還是動力工程都向大型、高速、大功率、高性能、高精度和輕結(jié)構(gòu)的方向發(fā)展，使得結(jié)構(gòu)動力學(xué)問題尤為突出和重要．在結(jié)構(gòu)系統(tǒng)中存在的材料彈塑性、構(gòu)建大變形等非線性因素，在產(chǎn)品的設(shè)計(jì)、結(jié)構(gòu)的控制中都需要考慮．對這些問題進(jìn)行分析的第一步就是建立一個(gè)與實(shí)際結(jié)構(gòu)基本相似的理想分析模型．任何一個(gè)實(shí)際結(jié)構(gòu)是一個(gè)連續(xù)系統(tǒng)，用偏微分方程才能夠精確地描述動力學(xué)問題．如果能求得偏微分方程的解析解，也就得到了問題的精確解；但是只有極少數(shù)的偏微分方程能夠得到精確解．客觀事物是多種多樣的，尤其工程技術(shù)中的問題更是如此．由于實(shí)際需要，必須尋求各種數(shù)值方法以獲得問題的近似解．動力學(xué)偏微分方程(組)含有空間變量和時(shí)間變量．?dāng)?shù)值求解含有空間變量和時(shí)間變量的偏微分方程(組)的近似解一般有空間離散和時(shí)間離散兩個(gè)步驟．對方程空間離散后，產(chǎn)生一組關(guān)于時(shí)間變量的常微分方程，然后利用時(shí)間離散方法求解常微分方程組．常用的空間離散方法有有限差分法、加權(quán)殘值法、伽遼金(Galerkin)法、有限元法、譜方法等[1–4]．

有限差分法是發(fā)展最早、比較成熟且應(yīng)用最為廣泛的求解偏微分方程的一種數(shù)值計(jì)算方法．有限差分法的數(shù)學(xué)基本原理是泰勒展開．任何差分格式都可以利用泰勒展開得到[5]．必須注意的是，對于內(nèi)部點(diǎn)和邊界點(diǎn)、均勻網(wǎng)格和不均勻網(wǎng)格，差分格式的表達(dá)是不一樣的．比如，對于內(nèi)部點(diǎn)，可以用中心差分格式，而邊界點(diǎn)只能用向前差分或向后差分格式；對于非均勻網(wǎng)格點(diǎn)，差分格式的表達(dá)式非常復(fù)雜，尤其是高階差分格式．由于以上特點(diǎn)，差分格式在實(shí)際應(yīng)用中僅限于低階的和均勻網(wǎng)格點(diǎn)．然而，對實(shí)際應(yīng)用中許多問題的模擬常常需要非均勻網(wǎng)格和高階差分格式．另外，低階差分格式在求解問題時(shí)為達(dá)到所需精度，需要大量的網(wǎng)格點(diǎn)，這無疑會增大計(jì)算量．

20世紀(jì) 70年代初發(fā)展起來的微分求積法(differential quadure method，簡稱 DQ 法)則是一種高精度、易于實(shí)施的求解偏微分方程(組)的數(shù)值方法．微分求積法是于 1971年由美國學(xué)者 Bellman和Casti[6]作為對積分求積思想的推廣而提出來的，其本質(zhì)是將函數(shù)在求解區(qū)域內(nèi)的每個(gè)網(wǎng)格處的導(dǎo)數(shù)值用域內(nèi)全部網(wǎng)格點(diǎn)上的函數(shù)值的加權(quán)線性和近似表示．由于微分求積法相對于傳統(tǒng)的數(shù)值計(jì)算方法具有原理簡單(不依賴于變分原理)、計(jì)算量小、精度高、易于編程等優(yōu)點(diǎn)，Bert等[7]于1987年首先將該法用于結(jié)構(gòu)元件的分析，從而掀起了對微分求積法及其在工程應(yīng)用中的研究熱潮．2008年，Malekzadeh等[8]用微分求積法研究了疊層厚圓拱的面內(nèi)自由振動．2009年，Matbuly等[9]用微分求積法對功能梯度裂縫梁的自由振動進(jìn)行分析，計(jì)算了其固有頻率．同年，劉金堂等[10]利用微分求積法對受面內(nèi)載荷軸向運(yùn)動板的橫向振動進(jìn)行分析，研究了軸向運(yùn)動速度、板材料剛度及長寬比對板橫向振動固有頻率及臨界速度的影響．2010年，Alibeigloo等[11]利用微分求積法對功能梯度圓柱殼進(jìn)行了靜力分析；2011年，Setoodeh等[12]利用微分求積法對功能梯度軸對稱圓柱殼的瞬態(tài)動力學(xué)性質(zhì)進(jìn)行了分析．2012年，Jam等[13]對中等厚度的功能梯度板進(jìn)行了非線性彎曲分析．2013年，Wang等[14]利用微分求積法研究了自由邊界條件下的菱形板的自由振動．2015年，Tornabene等[15]利用微分求積法研究了復(fù)合橢圓及橢圓柱的自由振動．2017年，Arani等[16]利用微分求積法對碳納米管電流變的夾層板的自由振動進(jìn)行了分析．

截至目前，國內(nèi)外學(xué)者將微分求積法用于動力學(xué)問題的研究還很少．以上微分求積法的發(fā)展和應(yīng)用都是針對各種工程結(jié)構(gòu)的靜力、穩(wěn)定性和自由振動分析，求解的方程是只含有空間變量的偏微分方程．本文將從理論上證明微分求積法是具有最高精度的高階差分法，是低階差分法的推廣，精度高，計(jì)算量小，并且其權(quán)系數(shù)可以求出顯示表達(dá)式，容易編制程序．從理論上把微分求積法推廣到動力學(xué)偏微分方程的求解，為動力學(xué)問題的研究提供一種簡單有效的分析方法．

1 微分求積法基本原理及其權(quán)系數(shù)的確定

1.1 微分求積法的基本原理

微分求積法的基本原理是將函數(shù)在求解區(qū)域內(nèi)的每個(gè)網(wǎng)格點(diǎn)處的導(dǎo)數(shù)值用域內(nèi)全部網(wǎng)格點(diǎn)上的函數(shù)值的加權(quán)線性和近似表示．簡單起見，以一維函數(shù)為例說明DQ法的基本思想．

設(shè)函數(shù) ()f x在區(qū)間[,]a b上有m階導(dǎo)數(shù)．將區(qū)間[,]a b用N個(gè)點(diǎn)分成 1N?個(gè)區(qū)間，如圖1所示．

圖1 一維問題節(jié)點(diǎn)分布Fig. 1 One dimensional problem node distribution

則函數(shù) ()f x在各節(jié)點(diǎn)處的各階導(dǎo)數(shù)值表示為

即

1.2 權(quán)系數(shù)的確定

在微分求積法發(fā)展的初期，其權(quán)系數(shù)矩陣均采用范德蒙矩陣求逆得到[6–7]．實(shí)踐表明，在節(jié)點(diǎn)數(shù)大于15時(shí)對應(yīng)的范德蒙矩陣是病態(tài)的，難以精確求出其逆矩陣，所以在運(yùn)用中一般節(jié)點(diǎn)數(shù)小于15．為了克服DQ法在計(jì)算權(quán)系數(shù)時(shí)的困難，Shu等[17]提出廣義微分求積法(generalized differential quadrature)，即通過對多項(xiàng)式線性向量空間的分析，利用拉格朗日插值多項(xiàng)式求導(dǎo)得出微分求積法權(quán)系數(shù)的顯式表達(dá)式．微分求積法權(quán)系數(shù)的確定方法如下：

首先假定一個(gè)偏微分方程的解函數(shù)可以用高次多項(xiàng)式來近似．可以選用不同的多項(xiàng)式的基來表達(dá)函數(shù) f( x)，比如選用基可以近似

表示為

若選用拉格朗日插值多項(xiàng)式基

其中

f( x)還可以用牛頓插值多項(xiàng)式、勒讓德多項(xiàng)式近似．

則有

將式(9)代入式(8)得

下面采用拉格朗日插值多項(xiàng)式基或式(7)確定權(quán)系數(shù)．為了敘述的簡明，令

從表達(dá)式(13)中可以看出，如果給定 xi，由式(8)就很容易計(jì)算出，進(jìn)而計(jì)算出，i≠j．但是，計(jì)算需要先求出二階導(dǎo)數(shù)．這是一項(xiàng)非常復(fù)雜的工作．根據(jù)線性向量空間理論：多項(xiàng)式的一組基可以唯一地由另外一組基來表達(dá)．因此，如果一組基滿足一個(gè)線性約束關(guān)系，比如式(1)，則另外一組基也滿足此關(guān)系．因此，由拉格朗日插值多項(xiàng)式確定的與多項(xiàng)式基確定的是相等的．因此滿足由多項(xiàng)式基時(shí)得到的關(guān)系式

式(7)兩端分別對 x求二階及高階導(dǎo)數(shù)，并利用式(13)就得到高階導(dǎo)數(shù)的權(quán)系數(shù)顯式的表達(dá)式

2 有限差分法的基本思想及差分格式的推導(dǎo)

一般使用的差分格式都是低階的[5]，即給定網(wǎng)格點(diǎn)的各階導(dǎo)數(shù)用一些鄰近點(diǎn)的函數(shù)值的加權(quán)和來表示，其系數(shù)由泰勒展式得到．以一維函數(shù)為例，對于m階導(dǎo)數(shù)，N個(gè)網(wǎng)格點(diǎn)可以產(chǎn)生?Nm階精度的差分近似，如果在整個(gè)求解區(qū)域上只有N個(gè)網(wǎng)格點(diǎn)，那么對于m階導(dǎo)數(shù)，?Nm階差分格式應(yīng)是最高精度的差分近似．1?N階差分格式是最高精度的差分近似，比如，設(shè)

將式(18)代入式(17)得到

將式(18)代入式(1)可得到高階導(dǎo)數(shù)的最高精度的高階差分格式權(quán)系數(shù)滿足的代數(shù)方程組

要得出高精度的高階差分格式，需要求解式(22)代數(shù)方程組，當(dāng)階數(shù)較高時(shí)，其求解是非常困難的．

3 微分求積法和最高精度有限差分法的等價(jià)性

為說明微分求積法和最高精度有限差分法的等價(jià)性．首先需滿足式(10)和式(21)方程組是等價(jià)的．顯然它們的第一個(gè)方程是等價(jià)的，即

第二個(gè)方程也是等價(jià)的，即

假設(shè)這兩個(gè)系統(tǒng)的前1p+個(gè)方程都是等價(jià)的，即

利用二項(xiàng)展開式

令1==ab，可以得到

利用式(26)，方程組(21)的前2p+ 個(gè)方程可以寫為

將式(25)、式(27)代入式(28)得

式(29)說明兩個(gè)系統(tǒng)的第2+p個(gè)方程也是等價(jià)的．因此，式(10)系統(tǒng)和式(21)系統(tǒng)是等價(jià)的．

下面證明式(16b)和式(16c)與式(22)方程組是等價(jià)的．顯然，式(16b)與式(22)的第一個(gè)方程是完全一樣的．下面證明滿足式(22)方程組．假設(shè)滿足式(22)，即

將式(31)代入式(22)得

由此證明了式(16b)和式(16c)與式(22)方程組的等價(jià)性．

綜上所述，微分求積法實(shí)質(zhì)上是最高精度的有限差分法．

4 微分求積法應(yīng)用于動力學(xué)偏微分方程

由前面討論可知微分求積法(DQ)與最高精度的高階差分法是等價(jià)的．高階差分格式具有高階截?cái)嗾`差，用高階差分格式求解問題時(shí)只需要較少的網(wǎng)格點(diǎn)就能達(dá)到所需精度，而且無論是求解區(qū)域的內(nèi)部點(diǎn)還是邊界點(diǎn)，其差分格式的表達(dá)形式是一樣的，因而高階差分法的實(shí)施要比低階差分法的實(shí)施容易．然而，用泰勒展開確定高階差分格式的權(quán)系數(shù)需要求解代數(shù)方程組，這一工作非常復(fù)雜，而微分求積法的權(quán)系數(shù)已給出顯式表達(dá)式，因此，解決實(shí)際問題時(shí)微分求積法則更易于實(shí)施．本節(jié)將微分求積法推廣到動力學(xué)偏微分方程的求解．

4.1 求解原理

設(shè)某一動力學(xué)問題的控制方程及邊界條件為

式中：N是一線性或非線性微分算子；u( x, t)是未知函數(shù)；B是一邊界算子；Γ是區(qū)域Ω的邊界；x和t分別是獨(dú)立的空間變量和時(shí)間變量．考慮用微分求積法對其進(jìn)行空間離散，以一維函數(shù)為例，設(shè)式(33)方程的解函數(shù)的空間變量x的變化范圍是[a, b]，取N個(gè)互異節(jié)點(diǎn)

用矩陣形式表示為

4.2 具體步驟

工程結(jié)構(gòu)中的動力學(xué)方程大都是四階及四階以上的高階偏微分方程．一個(gè)端點(diǎn)不止一個(gè)邊界條件，因此邊界條件的處理不能簡單進(jìn)行．下面用示例進(jìn)一步說明微分求積法求解動力學(xué)方程的步驟、邊界條件的處理方法以及微分求積法的計(jì)算量小、易于實(shí)施、高精度等特點(diǎn)．

例1 一受簡諧作用力的簡支梁的控制方程為

式中：EI是抗彎剛度；ρ是梁的密度；A為截面積；l為梁的長度；sinωPt為簡諧激勵．邊界條件為

例 2 由 Wickert等[18]給出的一線性軸向移動梁的控制方程為

式中：u( x, t)為梁振動的橫向位移；v為軸向移動速度；P為軸向緊力．邊界條件同例1．

(1)數(shù)值離散

假設(shè)梁的長度為 1,m，則上述問題的求解區(qū)間為[0,1]．在區(qū)間[0,1]取N個(gè)互異節(jié)點(diǎn)實(shí)例中N=11，微分求積法常取切比雪夫多項(xiàng)式的根作為節(jié)點(diǎn)坐標(biāo)

利用微分求積法則，方程中的未知函數(shù)在各節(jié)點(diǎn)處的各階導(dǎo)數(shù)值表示為

(2)處理邊界條件

在計(jì)算四階權(quán)系數(shù)矩陣時(shí)，應(yīng)用矩陣乘積得

將式(48)分別代入控制方程(40)和(42)得

式(49)和式(50)都是二階常微分方程組．

(3)求解常微分方程組

用時(shí)間離散方法如差分法、龍格庫塔法求解式(49)和式(50)，即得梁橫向振動的響應(yīng)．如果是非線性的動力學(xué)方程，可以結(jié)合非線性動力學(xué)理論對工程結(jié)構(gòu)進(jìn)行周期、混沌、分叉等非線性動力學(xué)性質(zhì)進(jìn)行分析．

為了驗(yàn)證上述微分求積法的正確性和精確性，分別將上述兩個(gè)示例的振動響應(yīng)曲線的微分求積解和解析解作了比較，第一個(gè)例子的解析解可以用分離變量法[19]得到，第二個(gè)例子的解析解見文獻(xiàn)[20]．圖 2和 3分別描繪了系統(tǒng)(40)和(42)的梁橫向振動的時(shí)間歷程曲線．圖 2中的參數(shù)取值為：圖 3中的參數(shù)取值為：由圖2和圖3可以看出，微分求積解曲線和解析解曲線幾乎是重合的，因而表明微分求積法是有效的，且只用 11個(gè)節(jié)點(diǎn)就達(dá)到較高的精度．

這兩個(gè)示例中的邊界條件是簡支邊界條件，用權(quán)系數(shù)矩陣修正法處理的邊界條件，即在形成高階權(quán)系數(shù)矩陣時(shí)就引入了邊界條件．當(dāng)然也可以用文獻(xiàn)所述的δ法來處理，相比較就不如權(quán)系數(shù)矩陣修正法方便．只要邊界條件中有低階導(dǎo)數(shù)為 0，那么就可用權(quán)系數(shù)矩陣修正法處理，讓此邊界條件精確滿足．

圖2 式(40)系統(tǒng)的時(shí)間歷程曲線Fig. 2 Time history curve of system(40)

圖3 式(42)系統(tǒng)的時(shí)間歷程曲線Fig. 3 Time history curve of system(42)

5 結(jié) 語

本文從理論上證明了微分求積法和最高精度的高階有限差分法的等價(jià)性．高階差分格式具有高階截?cái)嗾`差，用高階差分格式求解問題時(shí)只需要較少的網(wǎng)格點(diǎn)就能達(dá)到所需精度，而且無論是求解區(qū)域的內(nèi)部點(diǎn)還是邊界點(diǎn)，其差分格式的表達(dá)形式是一樣的，因而高階差分法的實(shí)施要比低階差分法的實(shí)施容易．然而，用泰勒展開確定高階差分格式的權(quán)系數(shù)，需要求解代數(shù)方程組，非常復(fù)雜，而微分求積法的權(quán)系數(shù)已給出顯式表達(dá)式，因此，實(shí)際應(yīng)用時(shí)利用微分求積法更易于實(shí)施．進(jìn)一步地將微分求積法推廣到動力學(xué)偏微分方程的求解，用數(shù)值算例闡明了微分求積法求解工程結(jié)構(gòu)動力學(xué)偏微分方程的方法步驟，并將微分求積解與解析解進(jìn)行比較，證明了微分求積法的有效性．以上研究結(jié)果表明了微分求積法具有高精度、易于實(shí)施、計(jì)算量小等優(yōu)點(diǎn)，為工程結(jié)構(gòu)的動力學(xué)分析提供了一種簡單有效的方法．

［1］ 李榮華，馮國忱. 微分方程數(shù)值解法[M]. 北京：高等教育出版社，2002.

［2］ Morton K W，Mayer D F. 偏微分方程數(shù)值解法[M].李治平，門大力，許現(xiàn)民，等譯. 北京：人民郵電出版社，2006.

［3］ 邱吉寶. 加權(quán)殘值法的理論與應(yīng)用[M]. 北京：宇航出版社，1991.

［4］ 向新民. 譜方法的數(shù)值分析[M]. 北京：科學(xué)出版社，2000.

［5］ 張文生. 科學(xué)計(jì)算中的偏微分方程有限差分法[M].北京：高等教育出版社，2006：40–55.

［6］ Bellman R E，Casti J. Differential quadrature and longterm integration[J]. Journal of Mathematical Analysis and Applications，1971，34(2)：235–238.

［7］ Bert C W，Jang S K，Striz A J. New methods for analyzing vibration of structural components[C]//Proceedings of the 28th Structures，Structural Dynamics and Materials Conference. Reston，VA：AIAA，1987：936–943.

［8］ Malekzadeh P，Setoodeh A R，Barmshouri E. A hybrid layerwise and differential quadrature method for in-plane free vibration of laminated thick circular arches[J]. Journal of Sound and Vibration，2008，315(1)：212–225.

［9］ Matbuly M S，Ragb O，Nassar M. Natural frequencies of a functionally graded cracked beam using the differential quadrature method[J]. Applied Mathematics and Computation，2009，215(6)：2307–2316.

［10］ 劉金堂，楊曉東，聞邦椿. 基于微分求積法的軸向運(yùn)動板的橫向振動分析[J]. 振動與沖擊，2009，28(3)：178–181.

［11］ Alibeigloo A，Nouri V. Static analysis of functionally graded cylindrical shell with piezoelectric layers using differential quadrature method[J]. Composite Structures，2010，92(8)：1775–1785.

［12］ Setoodeh A R，Tahani M，Selahi E. Hybrid layerwisedifferential quadrature transient dynamic analysis of functionally graded axisymmetric cylindrical shells subjected to dynamic pressure[J]. Composite Structures，2011，93(11)：2663–2670.

［13］ Jam J E，Maleki S，Andakhshideh A. Non-linear bending analysis of moderately thick functionally graded plates using generalized differential quadrature method[J]. International Journal of Aerospace Sciences，2012，1(3)：49–56.

［14］ Wang X W，Wu Z. Differential quadrature analysis of free vibration of rhombic plates with free edges[J]. Applied Mathematics and Computation，2013，225：171–183.

［15］ Tornabene F，F(xiàn)antuzzi N，Bacciocchi M，et al. Free vibrations of composite oval and elliptic cylinders by the generalized differential quadrature method[J]. Thin-Walled Structures，2015，97：114–129.

［16］ Arani A G，Jamali S A，Zarei H B A. Differential quadrature method for vibration analysis of electro-rheological sandwich plate with CNT reinforced nanocomposite facesheets subjected to electric field[J]. Composite Structures，2017，180：211–220.

［17］ Shu C，Richards B E. Application of generalized differential quadrature to solve two-dimensional incompressible Navier-Stokes equations[J]. International Journal for Numerical Methods in Fluids，1992，15(7)：791–798.

［18］ Wickert J A，Mote C D. Response and discretization methods for axially moving materials[J]. Applied Mechanics Reviews，1991，44(11S)：279–284.

［19］ 鄭兆昌. 機(jī)械振動[M]. 北京：機(jī)械工業(yè)出版社，1980.

［20］ Kong L，Parker R G. Approximate eigensolutions of axially moving beams with small flexural stiffness[J]. Journal of Sound and Vibration，2004，276(1)：459–469.