王冬梅，張 偉，劉寅立

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

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

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

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

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

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

1.1 微分求積法的基本原理

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

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

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

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

即

1.2 權系數(shù)的確定

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

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

表示為

若選用拉格朗日插值多項式基

其中

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

則有

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

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

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

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

2 有限差分法的基本思想及差分格式的推導

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

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

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

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

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

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

第二個方程也是等價的，即

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

利用二項展開式

令1==ab，可以得到

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

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

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

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

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

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

綜上所述，微分求積法實質上是最高精度的有限差分法．

4 微分求積法應用于動力學偏微分方程

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

4.1 求解原理

設某一動力學問題的控制方程及邊界條件為

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

用矩陣形式表示為

4.2 具體步驟

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

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

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

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

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

(1)數(shù)值離散

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

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

(2)處理邊界條件

在計算四階權系數(shù)矩陣時，應用矩陣乘積得

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

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

(3)求解常微分方程組

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

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

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

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

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

5 結 語

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

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

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

［3］ 邱吉寶. 加權殘值法的理論與應用[M]. 北京：宇航出版社，1991.

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

［5］ 張文生. 科學計算中的偏微分方程有限差分法[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］ 劉金堂，楊曉東，聞邦椿. 基于微分求積法的軸向運動板的橫向振動分析[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］ 鄭兆昌. 機械振動[M]. 北京：機械工業(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.