基于節(jié)點(diǎn)坐標(biāo)變量的桁架大位移非線性有限元法

2014-09-23 22:13:55劉樹堂

建筑科學(xué)與工程學(xué)報(bào) 2014年2期

劉樹堂

文章編號(hào)：16732049(2014)02013805

收稿日期：20131117

摘要：為求解桁架大位移問題, 提出了一種基于節(jié)點(diǎn)坐標(biāo)變量的非線性有限元法——以桿端節(jié)點(diǎn)坐標(biāo)向量為顯函變量寫出單元桿端力向量表達(dá)式,由單元桿端力向量裝配結(jié)構(gòu)非線性平衡方程。求解時(shí),首先根據(jù)矩陣微分理論求出單元桿端力向量關(guān)于桿端節(jié)點(diǎn)坐標(biāo)向量的導(dǎo)數(shù)矩陣,由該導(dǎo)數(shù)矩陣裝配結(jié)構(gòu)非線性微分平衡方程;然后按照Newton切線法原理建立等效線性逼近方程,引入邊界約束條件得到結(jié)構(gòu)節(jié)點(diǎn)坐標(biāo)的迭代公式。研究結(jié)果表明：該方法穩(wěn)定性好、精度高、收斂速度快且簡(jiǎn)單易用, 為求解桁架大位移問題提供了一種有效方法。

關(guān)鍵詞：非線性分析;Newton切線法;大位移;有限元法;平面扁桁架；節(jié)點(diǎn)

中圖分類號(hào)：TU311.4 文獻(xiàn)標(biāo)志碼：A

Nonlinear Finite Element Method Solving Large Displacement Problems of Trusses Subjected to Node Coordinate Variable

LIU Shutang

Abstract: In order to solve the large displacement problems of trusses, a nonlinear finite element method was proposed, which used the node coordinate variables as unknowns. An expression of the member end force vector was written in terms of members end coordinates, by which the global nonlinear equilibrium equation was fabricated. In solution, firstly, based on matrix differential theory, element derivative matrix was obtained with respect to the member end force vector, by which the global differential nonlinear equilibrium equation was fabricated. Secondly, based on Newton tangent method theory, the global equivalent linear matrix equation was established, then the structures boundary restraints was introduced, the iterative formulae for structure node coordinates were obtained. The research results show that the present method has good stability, high precision, quick convergence and easiness in use, and is very efficient for solving the large displacement problems of trusses.

Key words: nonlinear analysis; Newton tangent method; large displacement; finite element method; planar flat truss; node

0 引言

在桿系結(jié)構(gòu)的設(shè)計(jì)分析中,結(jié)構(gòu)通常假設(shè)為桁架(節(jié)點(diǎn)鉸接結(jié)構(gòu)),桿件軸線為直線,桿件材料為線彈性狀態(tài)。當(dāng)結(jié)構(gòu)位移較小時(shí),可忽略位移對(duì)結(jié)構(gòu)平衡方程的影響,并通常以結(jié)構(gòu)零態(tài)構(gòu)形建立結(jié)構(gòu)平衡方程。當(dāng)結(jié)構(gòu)位移較大時(shí)(如高聳、大跨及柔性結(jié)構(gòu)受荷載作用后[HJ]其位移通常已達(dá)到與結(jié)構(gòu)尺度同一量級(jí)的情況),位移對(duì)結(jié)構(gòu)平衡方程的影響不能忽略,需要基于荷載態(tài)構(gòu)形建立結(jié)構(gòu)非線性平衡方程,且需要采用非線性求解方法來求解。

桁架大位移問題非線性分析涉及以下3個(gè)方面的問題:①用何變量作為結(jié)構(gòu)非線性平衡方程的未知量;②采用何種方法建立結(jié)構(gòu)非線性平衡方程;③利用何種方法求解結(jié)構(gòu)非線性平衡方程。

結(jié)構(gòu)非線性平衡方程的未知量常用的有2種:節(jié)點(diǎn)位移和節(jié)點(diǎn)坐標(biāo)。傳統(tǒng)上,人們通常取節(jié)點(diǎn)位移為未知量建立結(jié)構(gòu)非線性平衡方程，但是近年來一些研究則選取節(jié)點(diǎn)坐標(biāo)作為未知量［13］。Liu等［4］研究也發(fā)現(xiàn),以節(jié)點(diǎn)坐標(biāo)為未知量建立的結(jié)構(gòu)非線性平衡方程列式簡(jiǎn)單,使求解方法也變得簡(jiǎn)單。

結(jié)構(gòu)非線性平衡方程的建立方法常用的有2種:直接建立方法［5］和勢(shì)能最小化方法［13］。直接建立方法是根據(jù)桿件內(nèi)力和外荷載直接建立總節(jié)點(diǎn)力向量平衡方程,該方法概念明確，應(yīng)用簡(jiǎn)單。勢(shì)能最小化方法是：先建立結(jié)構(gòu)總勢(shì)能泛函,將非線性幾何方程和應(yīng)力應(yīng)變關(guān)系代入總勢(shì)能泛函中,經(jīng)離散化處理得到關(guān)于未知量的泛函表達(dá)式,再根據(jù)總勢(shì)能最小化條件得到關(guān)于未知量的非線性平衡方程。勢(shì)能最小化方法也常忽略一些無法處理的高階項(xiàng)。

結(jié)構(gòu)非線性平衡方程的求解方法可分為2種主要類型:直接迭代法和微分線性逼近法［6］。直接迭代法的特點(diǎn)是:不涉及求導(dǎo)過程,只是利用原非線性方程建立迭代公式,它適用于不能或很難求導(dǎo)的非線性方程。該方法適用性較強(qiáng),但是收斂速度慢。微分線性逼近法的特點(diǎn)是:求出結(jié)構(gòu)非線性方程關(guān)于未知量的導(dǎo)數(shù)矩陣(Jacobi矩陣),據(jù)此建立原方程的線性逼近方程,得到關(guān)于未知量的迭代公式，進(jìn)而求解。目前建立桁架結(jié)構(gòu)非線性平衡方程時(shí)習(xí)慣于以節(jié)點(diǎn)位移作為未知量,很難實(shí)現(xiàn)關(guān)于節(jié)點(diǎn)位移未知量的求導(dǎo)過程,并且也沒有相關(guān)的研究報(bào)道。

對(duì)于目前桁架大位移問題的分析方法,孫煥純等［7］和許強(qiáng)等［8］認(rèn)為“以往對(duì)桁架結(jié)構(gòu)的大變形非線性分析,都是應(yīng)用最小勢(shì)能原理建立關(guān)于節(jié)點(diǎn)位移的非線性聯(lián)立平衡方程,求解的工作量大,尤其對(duì)多自由度的大型復(fù)雜桁架更為突出”。為了解決這些問題, 文獻(xiàn)［7］,［8］中提出了“交替迭代線性逼近法”,并有效應(yīng)用于2個(gè)典型的桁架大位移問題,但是對(duì)于一些荷載工況也不能得到問題的收斂解。

為了求解桁架大位移問題,筆者提出一種新的非線性有限元分析方法：①以桿端節(jié)點(diǎn)坐標(biāo)向量為顯函變量寫出單元桿端力向量表達(dá)式,由單元桿端力向量裝配結(jié)構(gòu)總節(jié)點(diǎn)力向量平衡方程;②根據(jù)矩陣微分理論,求出單元桿端力向量關(guān)于桿端節(jié)點(diǎn)坐標(biāo)向量的導(dǎo)數(shù)矩陣,由單元導(dǎo)數(shù)矩陣裝配結(jié)構(gòu)關(guān)于總節(jié)點(diǎn)坐標(biāo)向量的導(dǎo)數(shù)矩陣;③按照Newton切線法原理,建立結(jié)構(gòu)等效線性逼近方程;④引入邊界約束條件得到節(jié)點(diǎn)坐標(biāo)的迭代公式。本文方法的結(jié)構(gòu)非線性平衡方程列式簡(jiǎn)單且精確,求解過程也簡(jiǎn)單,收斂速度快,可有效求解桁架大位移問題。

1 結(jié)構(gòu)非線性平衡方程的建立方法

設(shè)荷載態(tài)下單元的方向向量和長(zhǎng)度分別為

νc(c)=［-I I］c

(1)

lc(c)=cTI66c

(2)

I66= I3 -I3-I3I3

(3)

式中:c為荷載態(tài)桿端節(jié)點(diǎn)坐標(biāo)向量;νc(c)，lc(c)分別為單元的方向向量和長(zhǎng)度;I為3階單位向量。

荷載態(tài)下單元軸力為

N(c)=EAlc(c)-lz lz

(4)

式中:N(c)為荷載態(tài)下單元軸力;lz為零態(tài)單元長(zhǎng)度;A為單元截面面積;E為材料彈性模量。

荷載態(tài)下單元桿端力向量f(c)為

f(c)=1 lc(c)-νc(c) νc(c)N(c)

(5)

將式(1)~(4)代入式(5),得到單元桿端力向量f(c)為

f(c)=EA(1 lz-1 cTI66c)I66c

(6)

由單元桿端力向量f(c),根據(jù)節(jié)點(diǎn)編號(hào)和對(duì)號(hào)入座原則,寫出以總節(jié)點(diǎn)坐標(biāo)向量為未知量的結(jié)構(gòu)非線性平衡方程為

F(C)=f(c)-B=0

(7)

式中:F(C)為結(jié)構(gòu)荷載態(tài)構(gòu)形總節(jié)點(diǎn)坐標(biāo)向量C的函數(shù);B為結(jié)構(gòu)荷載態(tài)總節(jié)點(diǎn)力向量;f(c)為由f(c)按節(jié)點(diǎn)編號(hào)和對(duì)號(hào)入座原則進(jìn)行裝配得到的向量。

2 結(jié)構(gòu)非線性平衡方程的求解方法

令式(7)關(guān)于結(jié)構(gòu)總節(jié)點(diǎn)坐標(biāo)向量C的導(dǎo)數(shù)為FC(C),即

FC(C)=dF(C) dC=df(c) dc

(8)[HJ1.9mm]

根據(jù)矩陣微分理論,方程式(8)中df(c) dc的求導(dǎo)結(jié)果為

df(c) dc=EA (cTI66c)3/2I66ccTI66+EA(1 lz-1 cTI66c)I66

(9)

根據(jù)Newton切線法原理,結(jié)構(gòu)非線性平衡方程的線性逼近方程式為

F(C(k))=F(C(k-1))+FC(C(k-1))? (C(k)-C(k-1))=0

(10)

根據(jù)節(jié)點(diǎn)坐標(biāo)約束條件,C可分為非約束節(jié)點(diǎn)坐標(biāo)Cf和約束節(jié)點(diǎn)坐標(biāo)Cr。設(shè)有一元素為0或1的正交變換矩陣It可完成以下變換

CfCr=ItC或

C=ITtCfCr

(11)

將式(11)代入式(10),然后方程兩邊再左乘It,有

ItF(C(k-1))+ItFC(C(k-1))ITt?（CfCr(k)-CfCr(k-1)）=0

(12)

令P=ItF(C(k-1)),K=ItFC(C(k-1))ITt,則式(12)可簡(jiǎn)寫為

P+K（CfCr(k)-

CfCr(k-1)）=0

(13)

對(duì)應(yīng)Cf和Cr的分塊,P和K也進(jìn)行相應(yīng)分塊,對(duì)于約束節(jié)點(diǎn)坐標(biāo)向量，有C(k)r-C(k-1)r=0，則式(13)可寫為

PfPr+Kff KfrKrf Krr

C(k)f-C(k-1)f 0=0

(14)

按矩陣乘法展開式(14),得到節(jié)點(diǎn)坐標(biāo)向量的迭代公式為

C(k)f=C(k-1)f-（Kff）-1Pf

(15)

3 算例分析

一平面扁桁架及其節(jié)點(diǎn)、桿件編號(hào)與節(jié)點(diǎn)荷載P見圖1。扁桁架桿件截面均為圓形,桿件截面直徑見表1，彈性模量E=200 GPa。扁桁架零態(tài)構(gòu)形節(jié)點(diǎn)x,y坐標(biāo)見表2，取右半部的結(jié)構(gòu)進(jìn)行分析。

圖1 平面扁桁架（單位:cm）

Fig.1 Planar Flat Truss (Unit:cm)

采用本文方法對(duì)1.0P荷載系（圖1）進(jìn)行分析,節(jié)點(diǎn)位移計(jì)算結(jié)果見表3,桿件軸力計(jì)算結(jié)果見表4。為便于比較,文獻(xiàn)［7］和ANSYS的計(jì)算結(jié)果也列于表3，4。通過比較可以看出,本文結(jié)果與文獻(xiàn)［7］及ANSYS計(jì)算結(jié)果是一致的。

1.0P荷載系節(jié)點(diǎn)y方向位移和桿件軸力迭代[CM)][HJ]

表1 扁桁架桿件截面直徑

Tab.1 Diameters of Member Section for Flat Truss

桿件編號(hào) 1,2,4,8,9 3,5,7,10 6

]截面直徑/cm 3.8 2.0 5.2

表2 扁桁架零態(tài)構(gòu)形節(jié)點(diǎn)x,y坐標(biāo)

Tab.2 Node Coordinates in x，y Directions for Flat Truss

節(jié)點(diǎn)編號(hào) x坐標(biāo) y坐標(biāo)

1 0.000 000 00 50.000 000 00

2 1 000.000 000 00 0.000 000 00

3 300.166 574 15 44.997 224 54

4 299.833 425 85 35.002 435 47

5 700.166 574 15 24.997 224 54

6 699.833 425 85 15.002 775 41

曲線分別如圖2，3所示。從圖2，3可以看出,迭代僅4次就滿足了節(jié)點(diǎn)坐標(biāo)容差ε=10-7的要求,文獻(xiàn)［7］中迭代20次達(dá)到相同結(jié)果,本文方法收斂速度較快。

由圖2，3還可以看出,迭代收斂前的曲線為較長(zhǎng)一段水平線,說明迭代收斂解是穩(wěn)定和可靠的。

對(duì)2.5P荷載系進(jìn)行分析,通過對(duì)比可以看出,本文結(jié)果與文獻(xiàn)［7］及ANSYS計(jì)算結(jié)果基本是一致的,只是在小數(shù)點(diǎn)后3位后有一些差別，這是由于方法不同產(chǎn)生的誤差。

2.5P荷載系節(jié)點(diǎn)y方向位移和桿件軸力迭代變化曲線分別如圖4，5所示。迭代8次就滿足了節(jié)點(diǎn)坐標(biāo)容差ε=10-7的要求,文獻(xiàn)［7］中迭代20次達(dá)到相同結(jié)果,本文方法收斂速度較快。由圖4，5可以看出,迭代收斂前的曲線為較長(zhǎng)一段水平線,說明迭代收斂解是穩(wěn)定和可靠的。

對(duì)2.6P荷載系進(jìn)行分析,對(duì)于2.6P荷載系,文獻(xiàn)［7］方法和ANSYS均不能求解。本文2.6P荷載系下節(jié)點(diǎn)位移和桿件軸力迭代曲線如圖6所示。迭代20次滿足了節(jié)點(diǎn)坐標(biāo)容差ε=10-7的要[CM(22]求[KG-*7]。[KG*3]由圖6可以看出,迭代收斂前的曲線為較長(zhǎng)一

表3 節(jié)點(diǎn)位移

Tab.3 Node Displacements

1.0P荷載系 2.5P荷載系 2.6P荷載系

節(jié)點(diǎn)位移本文方法文獻(xiàn)［7］方法 ANSYS 本文方法文獻(xiàn)［7］方法 ANSYS 本文方法

v1 -1.046 800 -1.047 400 -1.047 -1.620 90 -1.619 -1.619 -103.900 0

u3 -0.018 554 -0.018 563 -0.019 -0.149 12 -0.150 -0.150 -0.739 2

v3 -0.920 740 -0.921 225 -0.921 -4.583 30 -4.607 -4.606 -77.410 0

u4 -0.013 704 -0.013 711 -0.014 -0.180 38 -0.182 -0.182 0.255 5

v4 -0.920 910 -0.921 390 -0.921 -4.582 30 -4.605 -4.605 -77.390 0

u5 -0.020 650 -0.020 660 -0.021 -0.218 67 -0.220 -0.220 -0.837 6

v5 -0.502 240 -0.150 250 -0.503 -3.973 10 -3.998 -3.998 -35.840 0

u6 -0.004 703 -0.004 705 -0.005 -0.155 15 -0.156 -0.156 0.243 5

v6 -0.502 780 -0.503 030 -0.503 -3.975 00 -4.000 -4.000 -35.820 0

注:ui,vi分別為節(jié)點(diǎn)i在x,y方向的位移。

表4 桿件軸力

Tab.4 Axial Forces of Members

1.0P荷載系 2.5P荷載系 2.6P荷載系

桿件軸力本文方法文獻(xiàn)［7］方法 ANSYS 本文方法文獻(xiàn)［7］方法 ANSYS 本文方法

N1 -155 84.000 -15 584.300 -15 584.502 -64 293.400 -64 485.182 -64 484.939 -5 568.84

N2 -150 75.000 -15 075.100 -15 074.934 -13 266.800 -13 071.182 -13 071.001 77 937.50

N3 77.098 77.102 77.104 359.310 358.438 358.449 -1 301.68

N4 -128 98.000 -12 898.300 -12 898.483 -56 338.500 -56 521.954 -56 521.950 -9 621.14

N5 -2 704.200 -2 704.200 -2 704.233 -8 006.830 -8 015.200 -8 015.164 4 140.54

N6 -123 78.000 -12 378.300 -12 378.148 -4 942.960 -4 735.709 -4 735.696 73 808.40

N7 77.213 77.217 77.220 351.755 350.715 350.723 -1 304.25

N8 -156 28.000 -15 627.900 -15 628.101 -64 740.200 -64 935.253 -64 935.106 -5 420.73

N9 -150 80.000 -15 080.400 -15 080.215 -12 940.600 -12 741.762 -12 741.534 77 929.50

N10 -2 696.400 -2 696.400 -2 696.384 -8 321.230 -8 332.917 -8 332.953 4 105.19

注:N1~N10分別為1~10號(hào)桿件的軸力。

圖2 1.0P荷載系下節(jié)點(diǎn)y方向位移迭代曲線

Fig.2 Iteration Curves for Node ydirection Displacements Under 1.0P Load System

圖3 1.0P荷載系下桿件軸力迭代曲線

Fig.3 Iteration Curves for Axial Forces of Members Under 1.0P Load System

圖4 2.5P荷載系下節(jié)點(diǎn)y方向位移迭代曲線

Fig.4 Iteration Curves for Node ydirection Displacements Under 2.5P Load System

圖5 2.5P荷載系下桿件軸力迭代曲線

Fig.5 Iteration Curves for Axial Forces of Members Under 2.5P Load System

圖6 2.6P荷載系下節(jié)點(diǎn)y方向位移迭代曲線

Fig.6 Iteration Curves for Node ydirection Displacement Under 2.6P Load System

圖7 2.6P荷載系下結(jié)構(gòu)平衡狀態(tài)

Fig.7 Equilibrium State of Structure Under 2.6P Load System

段水平線,說明迭代收斂解是穩(wěn)定和可靠的。

2.6P荷載系下結(jié)構(gòu)平衡狀態(tài)如圖7所示。對(duì)大于2.6P荷載系也進(jìn)行了分析,均能夠獲得穩(wěn)定收斂解(分析結(jié)果略)。

4 結(jié) 語

(1)以節(jié)點(diǎn)坐標(biāo)為未知量列出桁架大結(jié)構(gòu)的非線性平衡方程,建立方程過程沒有任何近似,所得到的非線性平衡方程是精確的列式。

(2)基于本文的桁架結(jié)構(gòu)非線性平衡方程列式,可有效利用微分線性逼近這種精確方法進(jìn)行求解,所形成的求解方程很簡(jiǎn)單,易于應(yīng)用。

(3)大位移平面扁桁架算例研究表明,本文方法對(duì)于結(jié)構(gòu)非常接近失穩(wěn)的穩(wěn)定平衡狀態(tài)能夠獲得精[LL][LL]確解,且具有很高的穩(wěn)定性和很快的收斂速度。

(4)本文方法穩(wěn)定性好、精度高、收斂速度快且簡(jiǎn)單易用,為求解桁架大位移問題提供了一種有效方法,具有重要理論意義和工程應(yīng)用價(jià)值。

參考文獻(xiàn)：

References:

［1］OHKUBO S,WATADA Y,FUJIWAKI T.Nonlinear Analysis of Truss by Energy Minimization［J］.Computers & Structures,1987,27(1):129145.

［2］ KAZBERUK A,MIEDZIALOWSKI C Z,TRIBILLO R.Finite Element Discretization by Minimization of Elastic Strain Energy Method［J］.Finite Elements in Analysis and Design,1999,32(2):6370.

［3］ GRECOA M,GESUALDOA F A R,VENTURINIB W S.Nonlinear Positional Formulation for Space Truss Analysis［J］.Finite Elements in Analysis and Design,2007,42(12):10791086.

［4］ LIU S T,LONG Q L.A New Method Tracing Loaddeflection Equilibrium Path of a Doubly Nonlinear Truss［J］.Applied Mechanics and Materials,2012,166169:6872.

［5］沈世釗,徐崇寶,趙臣,等.懸索結(jié)構(gòu)設(shè)計(jì)［M］.2版.北京:中國(guó)建筑工業(yè)出版社,2006.

 SHEN Shizhao,XU Chongbao,ZHAO Chen,et al.Design of Cable Structures［M］.2nd ed.Beijing:China Architecture & Building Press,2006.

［6］熊洪允,曾紹標(biāo),毛云英.應(yīng)用數(shù)學(xué)基礎(chǔ)［M］.天津:天津大學(xué)出版社,1994.

 XIONG Hongyun,ZENG Shaobiao,MAO Yunying.Fundamentals of Applied Mathematics［M］.Tianjin:Tianjin University Press,1994.

［7］孫煥純,許強(qiáng),龍武智.桁架結(jié)構(gòu)幾何大變形分析的精確方法［J］.應(yīng)用力學(xué)學(xué)報(bào),2009,26(1):4550.

 SUN Huanchun,XU Qiang,LONG Wuzhi.Accurate Algorithm for Geometrically Large Deflection Analysis of Truss Structures［J］.Chinese Journal of Applied Mechanics,2009,26(1):4550.

［8］許強(qiáng),陳慶,孫煥純.大跨度桁架幾何大變形結(jié)構(gòu)分析的一種數(shù)值方法［J］.土木工程學(xué)報(bào),2009,42(1):1622.

 XU Qiang,CHEN Qing,SUN Huanchun.An Algorithm for Analysis of Large Span Truss Structures with Large Deformation［J］.China Civil Engineering Journal,2009,42(1):1622.