線性阻尼桿振動問題基于變量變換的多辛算法

鮑四元，鄧子辰

（1.蘇州科技大學 土木工程學院，江蘇 蘇州 215011；2.西北工業(yè)大學 工程力學系，西安 710072）

線性阻尼桿振動問題基于變量變換的多辛算法

鮑四元1，鄧子辰2

（1.蘇州科技大學 土木工程學院，江蘇 蘇州 215011；2.西北工業(yè)大學 工程力學系，西安 710072）

提出線性阻尼桿振動問題基于變量變換的多辛離散算法。首先通過變量變換把耗散系統(tǒng)化為保守系統(tǒng)。其次以變換變量組成狀態(tài)向量并采用中點離散方法構造中點Box多辛離散格式，然后，基于空間層的狀態(tài)變量建立矩陣形式的遞推關系，最后結合空間邊界條件和初始條件建立線性方程組求解。研究結果表明，構造的多辛算法不僅能夠保持系統(tǒng)守恒型幾何性質，通過狀態(tài)變量合理表示了邊界條件，而且能夠較準確地體現(xiàn)系統(tǒng)的耗散效應。

振動與波；變量變換；多辛；保結構；耗散

作為保結構算法的奠基人，馮康先生敏銳地覺察到數(shù)值算法構造過程中算法在保持原系統(tǒng)固有幾何性質方面具有重要的性能，并于1984年在雙微國際會議上提出了基于辛幾何原理的辛算法［1］，用于解決有限維Hamilton體系數(shù)值求解過程中的保結構問題。保結構算法的基本出發(fā)點是使得數(shù)值離散盡可能與原連續(xù)系統(tǒng)保持在統(tǒng)一框架下，即盡可能多地保持系統(tǒng)的固有幾何性質。自此，馮康先生開辟了保結構算法研究的新領域，至今國內外學者在辛算法方面取得了重大而豐富的研究成果[1–5]。

在多年的研究過程中，辛算法的一個不足之處在于空間依賴的連續(xù)系統(tǒng)并不能采用辛算法進行數(shù)值求解。因此，Bridges教授在辛算法的基礎上，針對無窮維Hamilton系統(tǒng)提出了多辛算法［6］，用以在數(shù)值求解過程中保持系統(tǒng)的局部幾何性質。最近，多辛算法在動力學方面得到不少應用和發(fā)展[2，7–9]。

實際的力學系統(tǒng)往往要考慮耗散效應，耗散效應的存在破壞了系統(tǒng)的多辛對稱性。近年來，胡偉鵬教授和鄧子辰教授針對耗散Hamilton系統(tǒng)，提出了廣義多辛算法［10］，將保結構算法理論體系拓寬至耗散系統(tǒng)，從而突破了多辛算法的一個“瓶頸”，但廣義多辛算法存在一定的近似和誤差。

變量變換是力學中重要理論和方法。例如，在振動問題的研究中，變量變換能夠將耦合振動解耦，分離成無耦合的簡正振動[11]。文獻[12]引入了三類含時的變量變換，將阻尼諧振子(DHO)中的阻尼從形式上消除，變換成簡諧振子(SHO)。

本文基于變量變換的方法，把阻尼項從桿振動微分方程中形式上消去，使得以新變量表示的振動微分系統(tǒng)是保守的，從而在多辛理論框架下，采用多辛算法研究了大阻尼桿的振動特性，即首先通過正則變換，構造桿振動控制方程的廣義多辛形式，并研究多辛形式中的數(shù)學表述。其次，采用中點離散方法離散多辛形式，構造其中點Box廣義多辛離散格式；隨后利用該格式求解，模擬阻尼桿的振動，最后研究振動過程中的阻尼耗散效應。本文的研究相對于廣義多辛算法而言，采用變量變換把耗散系統(tǒng)轉化為保守系統(tǒng)，因此精度更高，將進一步完善保結構算法的理論體系，同時為耗散動力學系統(tǒng)的數(shù)值求解提供新的途徑。

1 阻尼桿振動方程基于變換變量的的多辛形式

由于材料本構關系中含有阻尼項，桿件的振動需要考慮阻尼。單位長度懸臂阻尼桿的一維自由振動控制方程可以表述如下

其中u是縱向位移，α是縱向應力波的傳遞速度，c為不可忽略的阻尼系數(shù)，由于阻尼項的存在，桿振動振幅將會逐漸減小，常規(guī)的數(shù)值算法難以得到精確的衰減速度和振動時間。為了把系統(tǒng)變換為保守的，消去微分方程中的阻尼項，把式（1）改寫為

式（3）代入式（2），化簡得

式（4）中不再含U的1階偏導數(shù)，故易于寫成多辛形式。引入新變量和正交變換：V=Ut，W=Ux，得到與式（4）等價的多辛方程組形式為

式（5）可寫為Hamilton體系下的1階形式

其中Hamilton函數(shù)為

多辛算法的特色在于其滿足三種基本的守恒律，即：多辛守恒律、局部能量守恒律和局部動量守恒律。這也是多辛形式的優(yōu)勢所在。

2 基于變換新變量的多辛離散

中心差分法是最簡單的差分離散方法，在構造多辛格式時也最簡單，計算量最小。本文選取空間方向的步長為h，空間方向節(jié)點數(shù)為m，對應于i=1，2，…，m；時間方向的步長為τ，時間方向節(jié)點數(shù)為n，對應于j=1，2，……，n。對計算區(qū)域進行均勻劃分，狀態(tài)變量z在(xi,tj)網格點上的值記為對式（5）在空間和時間方向上分別用中點差分離散法近似，得到方程（5）的中點Box格式為

其中部分符號的含義如下

式（9）中z表示新變量U、V和W。故式（8）可寫為

其中i=1，2，……，m-1。

空間上桿左右兩端處均有邊界條件，宜對建立的中點Box廣義多辛格式在空間層上用矩陣表示。格式（10）可寫成如下的向量形式

式（12）矩陣下方的標注表示該矩陣的維數(shù)。而式（11）中i取不同值匯總后可化為如下矩陣形式

式（13）反映了不同空間層上狀態(tài)向量的迭代關系，圖1給出輔助直觀的理解，可由圖1中相鄰豎直線上的狀態(tài)向量的迭代體現(xiàn)。

圖1 空間層上的狀態(tài)變量迭代法示意圖（i對應于空間，j對應于時間變量）

由式（13），當l取不同值所得式子匯總后寫成矩陣 形 式 ，令記則得矩陣形式的線性方程組

其中R是塊三對角矩陣，以塊形式表示時維數(shù)是(m-1)×m，代入矩陣P、Q展開后的實際維數(shù)是3(m-1) (n-1)×3 mn。

方程（15）的Y向量含3 mn個待定元素，而矩陣R的行數(shù)為3(m-1)(n-1)小于3 mn。因此需結合初始條件和邊界條件求解出全部未知量。實際問題中可調節(jié)參數(shù)m和n的值，使得方程（15）可解。

3 初始條件和邊界條件的變換

利用方程（15）求解桿振動問題時，需考慮邊界條件和初始條件。但式（3）變換了位移變量，故邊界條件和初始條件需采用變換后的新變量表示。

對于零初始條件，一般采用新變量仍表示為零初始條件。

桿的空間約束包括自由端和固定端兩種。若約束是自由端，其邊界條件為把式（3）代入化簡得

若約束是固定端，其邊界條件為u(x,t)=0，代入式（3）得

對于給定的位移初始或力初始條件，也可根據(jù)式（3）得到新變量表示的形式。

4 數(shù)值實驗

以一端固定一端自由的懸臂均勻桿為例，其他約束可類似研究。桿振動控制方程仍為式（1），受初始位移（或外力）邊界條件。設桿長為1，參數(shù)α=1，阻尼系數(shù)c=0,2和4。桿的空間約束邊界條件如下

同時給定自由端單位縱向位移作為初值條件

本例屬于給定初值的自由振動問題，但不易得到精確解。

按式（16），邊界條件式（18）變換為

圖2 不同m，n值所得自由端位移隨時間變化圖

圖2中，(m，n)分別取(11，13)、(17，19)和(23，25)時，自由端的位移均逐漸衰減，而三者中取(11，13)時誤差相對大一些，尤其后半段時間位移出現(xiàn)波動，實際上后部分的位移可由平滑處理或擬合趨近為零。一般認為，圖1的計算網格越細，結果應更準確一些，即圖2（c）的結果相對更準確。

阻尼系數(shù)c對桿自由端位移的影響也反映在圖2中。阻尼系數(shù)大時，自由端位移衰減得更快。另外，圖3和圖4給出不同阻尼系數(shù)時桿中點的位移和應變變化，其中參數(shù)m取16，n取18。

圖3 桿中點的位移變化（m=16，n=18）

圖4 桿中點的應變變化（m=16，n=18）

5 結語

本文給出了線性阻尼桿振動基于變量變換的多辛算法。所提算法滿足多辛算法的多辛守恒律、局部能量守恒律和局部動量守恒律。計算最終轉換為針對線性方程組求解出所有狀態(tài)變量，通過變量變換描述線性阻尼引起的耗散，即表述阻尼時未采用近似，從而更好地模擬并解決了實際問題。

所提方法對于某些類似的非保守系統(tǒng)的多辛算法具有重要的借鑒意義。

[1]馮康，秦孟兆.哈密爾頓系統(tǒng)的辛幾何算法[M].杭州：浙江科學技術出版社，2003

[2]秦孟兆，王雨順.偏微分方程中的保結構算法[M].杭州：浙江科技出版社，2011.

[3]鐘萬勰.力、功、能量與辛數(shù)學[M].大連：大連理工大學出版社，2007.

[4]鐘萬勰，辛破繭.辛拓展新層次[M].大連：大連理工大學出版社，2011.

[5]YAO WEIAN,ZHONG WANXIE.Symplectic elasticity [M].Singerpore:world scientific press,2009.

[6]BRIDGES T J.Multi-symplectic structures and wave propagation[J].MathematicalProceedingsofthe Cambridge Philosophical Society,1997,121(1):147-190.

[7]胡偉鵬，鄧子辰.無限維動力學系統(tǒng)的保結構分析方法[M].西安：西北工業(yè)大學出版社，2015.

[8]HU W P,HAN S M,DENG Z C,et al.Analyzing dynamic response of non-homogeneous string fixed at both ends[J].International Journal of Non-Linear Mechanics,2012, 47(10):1111-1115.

[9]王俊杰.彈性波方程的多辛Preissmann格式計算[J].地球物理學進展，2014，29(4)：1758-1765.

[10]胡偉鵬，鄧子辰.大阻尼桿振動的廣義多辛算法[J].動力學與控制學報，2013，11(1)：1-4.

[11]丁光濤.阻尼諧振子的拉格朗日函數(shù)和哈密頓函數(shù)[J].大學物理，2009，28(5)：13-14.

[12]丁光濤.阻尼諧振子到簡諧振子的變換[J].大學物理，2013，32(1)：1，17.

[13]鐘萬勰.應用力學的辛數(shù)學方法[M].北京：高等教育出版社，2006.

Multi-symplecticAlgorithm for VibrationAnalysis of the Rods with Linear Damping

BAO Si-yuan1,DENG Zi-chen2
(1.School of Civil Engineering,University of Science and Technology of Suzhou, Suzhou 215011,Jiangsu China; 2.Department of Engineering Mechanics,Northwestern Polytechnic University,Xi’an 710072,China)

The discrete multi-symplectic alogrithim for vibration problem analysis of the rods with linear damping is proposed.First of all,the dissipation system is transformed into a conservative system by variable transformation.Then,the state vector is formed by the transformed variables,and the midpoint box multi-symplectic form is generated based on the midpoint discrete method.Finally,the matrix-type recurrence formula is established based on the state vector and the linear equation set including the boundary condition and the initial condition is formed and solved.The numerical results show that the proposed multi-symplectic algorithm can reserve the conservative geometrical properties of the system,reasonably express the boundary condition by the state variables,and accurately reflect the dissipation effect of the system.

vibration and wave;variable transformation;multi-symplectic;structure-preserving;dissipation

O326；O422.6

：A

：10.3969/j.issn.1006-1335.2017.01.004

1006-1355(2017)01-0016-04

2016-08-26

國家自然科學基金資助項目（11202146）

鮑四元（1980－），男，安徽省六安市人，博士，副教授，研究方向為非線性動力學及計算。

鄧子辰，男，博士生導師，長江學者。E-mail:dweifan@nwpu.edu.cn