拉壓剛度不同桁架的動力參變量變分原理和保辛算法

高 強，張洪武，張 亮，鐘萬勰

(大連理工大學 工程力學系，工業裝備結構分析國家重點實驗室，大連 116024)

經典彈性理論描述的是拉壓彈性模量相同材料的力學問題，但實際工程中由于不同的功能需求，從而對結構形式或構成這些結構的材料有不同要求。譬如，工程中常用的混凝土材料，表現出拉、壓性質不同的特點;而對于一些展開結構，為達到其設計性能，必須采用特殊的索、膜結構，這些索、膜部件同樣表現出不同的拉壓性質。具有拉、壓不同性質的材料或結構的力學分析，體現出較強的非線性特征，需要針對這類問題發展有效的求解算法。

拉、壓模量不同問題的理論研究始于上個世紀40年代。Timoshenko［1］提出了雙模量的概念。Ambartsumyan［2］研究了拉、壓模量不同圓柱殼的軸對稱性問題，并依據主應力的正負，給出了此類問題在主方向上的本構關系。Ambartsumyan［3］總結分析了大量新型材料的試驗數據，將拉、壓不同模量材料的本構關系總結為雙直線模型，并且詳細論述了彈性系數的選取。拉壓不同模量問題的研究有著廣泛的應用背景，張拉整體結構是典型的含有不同拉壓剛度的結構。文獻［4］推導了張拉整體結構平衡構型附近的線性化動力學模型。文獻［5］研究了張拉整體結構的動力學行為和控制策略。文獻［6］對張拉整體結構的動力學特點、求解方法和存在的問題作了綜述。

在數值求解方面，國內學者進行了大量的相關研究。張允真等［7］構造了雙模量問題求解的有限元格式，進行迭代求解。楊海天等［8］提出了用初應力法迭代求解雙模量彈性問題。劉相斌等［9］進一步討論了雙模量問題的剪切模量，提出了加速收斂因子的思想。He等［10］詳細討論了傳統迭代求解方法的收斂問題。楊海天等［11］討論了拉壓雙模量問題的動力分析問題。葉志明等［12］簡述了不同模量彈性問題理論及其有限元法的研究與發展。楊海天等［13］利用光滑函數技術，提出光滑化的拉壓不同彈性模量問題的應力應變關系，與有限元方法相結合，建立了拉壓不同模量一維連續體與桁架結構的數值求解模型。對于這類結構的特殊非線性問題，如何保證算法的穩定性是數值模擬的主要問題之一。

根據結構力學與最優控制的模擬關系，鐘萬勰等［14-15］建立和發展了參變量變分原理和基于此變分原理的一系列數值算法，已成功應用于彈塑性［14-15］、接觸［16-18］、摩擦［17-18］和納米管范德華力模擬［19］等非線性問題分析。本文在參變量變分原理基礎上，建立了由拉壓剛度不同桿單元組成的桁架結構的動力學參變量變分原理。將拉壓剛度不同桁架問題的非線性動力分析轉換為線性互補問題求解。結合時間有限元方法構造了求解此問題的保辛數值積分方法。此方法不需要剛度矩陣更新和迭代，計算過程高效、穩定性好。

1 拉壓不同剛度桿單元的動力參變量變分原理

圖1 拉壓剛度不同桿的本構關系Fig.1 The constitutive relation for rod with different module in tension and compression

考慮具有不同拉壓剛度的桿，設拉伸和壓縮剛度分別為K(+)和K(-)，K(+)≠K(-)，可分為兩種情況，如圖1。在圖1(b)中，如果壓縮模量等于零，則成為繩索。桿的本構關系為:

其中:

首先假設K(+)＜K(-)，即圖1(a)所示情況。將本構關系寫為如下的統一形式:

則要求:

上式與如下的互補關系等價，即:

則桿的勢能為:

因此，當K(-)＞K(+)時，桿的參變量變分原理為:

其中:下標u表示只對Δu進行變分，而λ作為參變量不變分。同理，如果K(-)＜K(+)，可類似寫出參變量變分原理為:

當然，也可將本構關系寫為如下的形式，即:

則當K(-)＞K(+)時，參變量變分原理為:

而當K(-)＜K(+)時，參變量變分原理為:

方程(7)，(8)，(10)和(11)描述的四種參變量變分原理可統一寫為:

符號sign表示取符號，其定義為:

容易證明方程(12)和(13)給出的參變量變分原理與方程(1)和(2)給出的平衡方程等價，具體的證明過程見文獻［9-10］。

以上給出的是單根桿的參變量變分原理，將桿組合成桁架時，涉及到局部坐標和整體坐標之間的轉換，局部坐標和整體坐標之間的關系如圖2所示。

圖2 坐標轉換Fig.2 The coordinate transformation

在局部坐標下有:

局部坐標與整體坐標下的位移之間的關系為:

其中:

因此有:

其中:

對于空間結構，轉換關系依然是方程(19)，只是:

其中:

給出位移在局部坐標和整體坐標之間的關系后，既可在總體坐標下，將桁架所有的桿單元按照有限元的方式進行集合，則得到整個桁架結構的參變量變分原理為:

其中:

其中Ne表示桿單元的數量。

下面考慮動力學問題，本文考慮的是無阻尼系統，因此動力學方程可通過Euler-Lagrange方程給出。Lagrange函數L的定義為:

其中T和U分別表示動能和勢能。上文已經給出了不同拉壓剛度桿的勢能U，而動能T為:

其中M是質量矩陣。因此有:

則動力學方程可由Hamilton變分原理給出，即:

變分后可給出動力學運動方程為:

當然，動力學方程還受互補條件約束，即方程(24)。

2 保辛方法

動力系統積分時，選取一個時間步長η，于是得到一系列等步長的時刻

在一個典型的積分步長t∈［tk-1，tk］內，將位移q(t)、參變量λ(t)和外力f(t)用線性函數近似，即:

將它們代入方程(29)中的作用量S，并積分得到近似作用量為:

其中:

根據離散Hamilton正則方程:

并通過方程(35)可得到:

其中:

將方程(36)帶入互補條件(24)得到:

求解線性互補問題(39)，可求得λ1，然后根據方程(36)和(37)可計算出q1和p1，從而完成一個時間步的積分。

3 數值算例

圖3 含有繩單元的單質點結構Fig.3 One DOF system with string elements

其中:

采用時間步長η=0.1 s積分到100 s，得到的位移和動量如圖4，其中實線表示本文方法計算結果，圓圈表示解析解，可以看到本文方法積分得到的結果與解析解非常吻合，證明了本文方法的正確性。若在質點上施加外載荷f(t)=sin(t)N，其它參數同上，則積分得到的位移和動量如圖5。

圖4 自由振動質點的位移和動量Fig.4 The displacement and momentum responses for free vibration

圖5 強迫振動質點的位移和動量Fig.5 The displacement and momentum responses for forced vibration

算例2:考慮如圖6所示的桁架結構，圖中實線表示拉壓模量相同的桿，它們具有相同的楊式模量和密度，分別為E=105N/m2和ρ=3×103kg/m3。虛線表示繩，即具有拉伸模量，而受壓時模量為0，左右兩根斜的繩具有相同拉伸模量E=2×105N/m2，中間的豎繩拉伸模量E=105N/m2，所有繩的質量忽略，所有桿和繩的面積為A=10-3m2。

圖6 具有繩索單元的桁架結構Fig.6 The truss structure with string elements

考慮結構的自由振動，取所有節點初始位移為0，而初始動量為:節點1水平動量為1，豎直動量為0，節點2水平動量為0，豎直動量為1，節點3水平動量為-1，豎直動量為0。采用步長η=0.2(s)積分到400 s，得到的各節點位移和動量響應如圖7，其中圖7(a)和(b)分別給出了3個節點x和y方向的位移，圖7(c)和(d)分別給出了3個節點x和y方向的動量。由于結構的對稱性和初始條件的對稱性，其響應也具有對稱性，而計算結果很好的體現了這種對稱性。3條繩索的參變量變化如圖8，參變量不等于零表示繩索處于松弛狀態。

圖7 平面桁架的響應Fig.7 The responses of the truss structure

采用步長η=0.2(s)積分到1 000 s，得到的系統能量的相對誤差分別如圖9。圖9表明本文方法給出的積分方法不會引入人工阻尼，能量始終在一定范圍內變化，這是保辛方法的典型特征。

圖8 參變量隨時間變化Fig.8 The parametric variables

圖9 能量的相對誤差Fig.9 The relative error of the Hamiltonian function

4 結論

本文建立了由拉、壓剛度不同桿單元組成的桁架結構的動力學參變量變分原理，將拉、壓剛度不同桁架問題的非線性動力分析轉換為線性互補問題求解。結合時間有限元方法構造的保辛數值積分方法，在時間步內不需要迭代和剛度矩陣更新，可精確、高效的求解此類問題，且計算過程穩定。該方法也可自然推廣到3-D桁架結構。

［1］ Timoshenko S.Strength of Materials，part II:Advanced theoryand problems ［M］. 2nded. Princeton:Van Nostrand，1941.

［2］Ambartsumyan S A.The Axisymmetric problem of circular cylindrical shell made of materials with different stiffness in tension and compression［J］. AkademiiNauk SSSR，Mekhanika，1965(4):77-85.

［3］Ambartsumyan S A著，張允真譯.不同模量彈性理論［M］.北京:中國鐵道出版社，1986.

［4］ Cornel S，Martin C，Robert E S.Linear dynamics of tensegrity structures［J］.Engineering Structures，2002，24:671-685.

［5］ Bel H A N，Smith I F C.Dynamic behavior and vibration control of a tensegrity structure［J］.International Journal of Solids and Structures，2010，47:1285-1296.

［6］ Mirats T，Hernandez J.Tensegrity Frameworks:Dynamic analysis review and open problems［J］.Mechanism and Machine Theory，2009，44:1-18.

［7］張允真，王志峰.不同拉壓模量彈性力學問題的有限元法［J］.計算結構力學及其應用，1989，6(1):236-254.

ZHANG Yun-zhen， WANG Zhi-feng.The finite element method for elasticity with different moduli in tension and compression［J］.Computational Structural Mechanics and Applications，1989，6(1):236-254.

［8］楊海天，鄔瑞鋒，楊克儉，等.初應力法解拉壓雙彈性模量問題［J］.大連理工大學學報，1992，32(1):35-39.

YANG Hai-tian， WU Rui-feng， YANG Ke-jian， et al.Solution to problem of dual extension-compression elastic modulus with initial stress method［J］.Journal of Dalian University of Technology，1992，32(1):35-39.

［9］劉相斌，張允真.拉壓不同模量有限元法剪切彈性模量及加速收斂［J］.大連理工大學學報，1994，34(6):641-645.

LIU Xiang-bin，ZHANG Yun-zhen.Modulus of elasticity in shear and accelerate convergence ofdifferentextension compression elastic modulus finite element method ［J］.Journal of Dalian University of Technology，1994，34(6):641-645.

［10］ He X T，Zheng Z L，Sun J Y.Convergence analysis of a finite element method based on different moduli in tension and compression ［J］. InternationalJournalofSolids and Structures，2009，46:3734-3740.

［11］楊海天，楊克儉，張錫成，等.拉壓雙模量問題的動力分析［J］.計算結構力學及其應用，1993，10(4):438-448.

YANG Hai-tian，YANG Ke-jian，ZHANG Xi-cheng，et al.Dynamic analysis forthe problem ofdualextensioncompression elastic modulus［J］.Computational Structural Mechanics and Applications，1993，10(4):438-448.

［12］葉志明，陳 彤，姚文娟.不同模量彈性問題理論及有限元法研究進展［J］.力學與實踐，2004，26(2):9-14.

YE Zhi-ming，CHEN Tong，YAO Wen-juan.Progress in elasticity theory with differentmoduliin tension and compression and related FEM［J］.Mechanics and Practice，2004，26(2):9-14.

［13］楊海天，張曉月，何宜謙.基于敏度分析的拉壓不同模量桁架問題的數值分析［J］.計算力學學報，2011，28(2):237-242.

YANG Hai-tian，ZHANG Xiao-yue，HE Yi-qian.Sensitivity analysis based numerical solution for truss structures with bimodulus［J］.Journal of Computational Mechanics，2011，28(2):237-242.

［14］鐘萬勰，張洪武，吳承偉.參變量變分原理及其在工程中的應用［M］.北京:科學出版社，1997.

［15］張洪武.參變量變分原理與材料和結構力學分析［M］.北京:科學出版社，2009.

［16］ Zhang H W，Xu W L，Di S L，et al.Quadratic programming method in numerical simulation of metal forming process［J］.Computer Methods in Applied Mechanics and Engineering，2002，191:5555-5578.

［17］ Zhan H W，He S Y，Li X S，Wriggers P.A new algorithm for numerical solution of 3D elastoplastic contact problems with orthotropic friction law［J］.Computational Mechanics，2004，34:1-14.

［18］ Zhang H W，He S Y，Li X S.Two aggregate-function-based algorithms for analysis of 3D frictional contact by linear complementarity problem formulation［J］.Computer Methods in Applied Mechanics and Engineering， 2005， 194:5139-5158.

［19］ Zhang H W，Wang J B，Ye H F，et al.Parametric variational principle and quadratic programming method for van der Waals force simulation of parallel and cross nanotubes［J］.International Journal of Solids and Structures，2006，44:2783-2801.