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

高 強(qiáng)，張洪武，張 亮，鐘萬勰

(大連理工大學(xué) 工程力學(xué)系，工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024)

經(jīng)典彈性理論描述的是拉壓彈性模量相同材料的力學(xué)問題，但實(shí)際工程中由于不同的功能需求，從而對(duì)結(jié)構(gòu)形式或構(gòu)成這些結(jié)構(gòu)的材料有不同要求。譬如，工程中常用的混凝土材料，表現(xiàn)出拉、壓性質(zhì)不同的特點(diǎn);而對(duì)于一些展開結(jié)構(gòu)，為達(dá)到其設(shè)計(jì)性能，必須采用特殊的索、膜結(jié)構(gòu)，這些索、膜部件同樣表現(xiàn)出不同的拉壓性質(zhì)。具有拉、壓不同性質(zhì)的材料或結(jié)構(gòu)的力學(xué)分析，體現(xiàn)出較強(qiáng)的非線性特征，需要針對(duì)這類問題發(fā)展有效的求解算法。

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

在數(shù)值求解方面，國(guó)內(nèi)學(xué)者進(jìn)行了大量的相關(guān)研究。張?jiān)收娴龋?］構(gòu)造了雙模量問題求解的有限元格式，進(jìn)行迭代求解。楊海天等［8］提出了用初應(yīng)力法迭代求解雙模量彈性問題。劉相斌等［9］進(jìn)一步討論了雙模量問題的剪切模量，提出了加速收斂因子的思想。He等［10］詳細(xì)討論了傳統(tǒng)迭代求解方法的收斂問題。楊海天等［11］討論了拉壓雙模量問題的動(dòng)力分析問題。葉志明等［12］簡(jiǎn)述了不同模量彈性問題理論及其有限元法的研究與發(fā)展。楊海天等［13］利用光滑函數(shù)技術(shù)，提出光滑化的拉壓不同彈性模量問題的應(yīng)力應(yīng)變關(guān)系，與有限元方法相結(jié)合，建立了拉壓不同模量一維連續(xù)體與桁架結(jié)構(gòu)的數(shù)值求解模型。對(duì)于這類結(jié)構(gòu)的特殊非線性問題，如何保證算法的穩(wěn)定性是數(shù)值模擬的主要問題之一。

根據(jù)結(jié)構(gòu)力學(xué)與最優(yōu)控制的模擬關(guān)系，鐘萬勰等［14-15］建立和發(fā)展了參變量變分原理和基于此變分原理的一系列數(shù)值算法，已成功應(yīng)用于彈塑性［14-15］、接觸［16-18］、摩擦［17-18］和納米管范德華力模擬［19］等非線性問題分析。本文在參變量變分原理基礎(chǔ)上，建立了由拉壓剛度不同桿單元組成的桁架結(jié)構(gòu)的動(dòng)力學(xué)參變量變分原理。將拉壓剛度不同桁架問題的非線性動(dòng)力分析轉(zhuǎn)換為線性互補(bǔ)問題求解。結(jié)合時(shí)間有限元方法構(gòu)造了求解此問題的保辛數(shù)值積分方法。此方法不需要?jiǎng)偠染仃嚫潞偷?jì)算過程高效、穩(wěn)定性好。

1 拉壓不同剛度桿單元的動(dòng)力參變量變分原理

圖1 拉壓剛度不同桿的本構(gòu)關(guān)系Fig.1 The constitutive relation for rod with different module in tension and compression

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

其中:

首先假設(shè)K(+)＜K(-)，即圖1(a)所示情況。將本構(gòu)關(guān)系寫為如下的統(tǒng)一形式:

則要求:

上式與如下的互補(bǔ)關(guān)系等價(jià)，即:

則桿的勢(shì)能為:

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

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

當(dāng)然，也可將本構(gòu)關(guān)系寫為如下的形式，即:

則當(dāng)K(-)＞K(+)時(shí)，參變量變分原理為:

而當(dāng)K(-)＜K(+)時(shí)，參變量變分原理為:

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

符號(hào)sign表示取符號(hào)，其定義為:

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

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

圖2 坐標(biāo)轉(zhuǎn)換Fig.2 The coordinate transformation

在局部坐標(biāo)下有:

局部坐標(biāo)與整體坐標(biāo)下的位移之間的關(guān)系為:

其中:

因此有:

其中:

對(duì)于空間結(jié)構(gòu)，轉(zhuǎn)換關(guān)系依然是方程(19)，只是:

其中:

給出位移在局部坐標(biāo)和整體坐標(biāo)之間的關(guān)系后，既可在總體坐標(biāo)下，將桁架所有的桿單元按照有限元的方式進(jìn)行集合，則得到整個(gè)桁架結(jié)構(gòu)的參變量變分原理為:

其中:

其中Ne表示桿單元的數(shù)量。

下面考慮動(dòng)力學(xué)問題，本文考慮的是無阻尼系統(tǒng)，因此動(dòng)力學(xué)方程可通過Euler-Lagrange方程給出。Lagrange函數(shù)L的定義為:

其中T和U分別表示動(dòng)能和勢(shì)能。上文已經(jīng)給出了不同拉壓剛度桿的勢(shì)能U，而動(dòng)能T為:

其中M是質(zhì)量矩陣。因此有:

則動(dòng)力學(xué)方程可由Hamilton變分原理給出，即:

變分后可給出動(dòng)力學(xué)運(yùn)動(dòng)方程為:

當(dāng)然，動(dòng)力學(xué)方程還受互補(bǔ)條件約束，即方程(24)。

2 保辛方法

動(dòng)力系統(tǒng)積分時(shí)，選取一個(gè)時(shí)間步長(zhǎng)η，于是得到一系列等步長(zhǎng)的時(shí)刻

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

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

其中:

根據(jù)離散Hamilton正則方程:

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

其中:

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

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

3 數(shù)值算例

圖3 含有繩單元的單質(zhì)點(diǎn)結(jié)構(gòu)Fig.3 One DOF system with string elements

其中:

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

圖4 自由振動(dòng)質(zhì)點(diǎn)的位移和動(dòng)量Fig.4 The displacement and momentum responses for free vibration

圖5 強(qiáng)迫振動(dòng)質(zhì)點(diǎn)的位移和動(dòng)量Fig.5 The displacement and momentum responses for forced vibration

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

圖6 具有繩索單元的桁架結(jié)構(gòu)Fig.6 The truss structure with string elements

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

圖7 平面桁架的響應(yīng)Fig.7 The responses of the truss structure

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

圖8 參變量隨時(shí)間變化Fig.8 The parametric variables

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

4 結(jié)論

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

［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著，張?jiān)收孀g.不同模量彈性理論［M］.北京:中國(guó)鐵道出版社，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］張?jiān)收妫踔痉?不同拉壓模量彈性力學(xué)問題的有限元法［J］.計(jì)算結(jié)構(gòu)力學(xué)及其應(yīng)用，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］楊海天，鄔瑞鋒，楊克儉，等.初應(yīng)力法解拉壓雙彈性模量問題［J］.大連理工大學(xué)學(xué)報(bào)，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］劉相斌，張?jiān)收?拉壓不同模量有限元法剪切彈性模量及加速收斂［J］.大連理工大學(xué)學(xué)報(bào)，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］楊海天，楊克儉，張錫成，等.拉壓雙模量問題的動(dòng)力分析［J］.計(jì)算結(jié)構(gòu)力學(xué)及其應(yīng)用，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ìn)展［J］.力學(xué)與實(shí)踐，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］楊海天，張曉月，何宜謙.基于敏度分析的拉壓不同模量桁架問題的數(shù)值分析［J］.計(jì)算力學(xué)學(xué)報(bào)，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］鐘萬勰，張洪武，吳承偉.參變量變分原理及其在工程中的應(yīng)用［M］.北京:科學(xué)出版社，1997.

［15］張洪武.參變量變分原理與材料和結(jié)構(gòu)力學(xué)分析［M］.北京:科學(xué)出版社，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.