有限單元法和ANSYS程序(FiniteElementMethodandANSYSProgram)§

5.1變分法(VariationMethods)簡介

變分法的定義:主要是研究泛函(functional)及其極值的求解方法.所謂泛函,就是以函數作為自變量的函數,即泛函就是函數的函數.彈性力學變分法所研究的泛函,主要是指能量(energy),能量是應力應變的函數,而應力應變又是坐標、力、材料性質等的函數，所以能量是泛函。因此，彈性力學中的變分法又稱為能量法(energymethod).一、形變勢能

按照材料力學(themechanicsofmaterials),當彈性體僅受一個方向力作用,例如x方向，其stresscomponent為

x,相應的

straincomponent為

x,那末單位體積的應變能（

strainenergy）

is

x

x/2.即U1又稱為變形勢能密度或單位比能（specificenergy）．同樣，當單元體只受剪stresscomponents

xy

作用時，相應的剪straincomponent

xy時,單位體積的應變能is

xy

xy/2.

設彈性受全部六個應力分量的作用，對于平面問題,,而如果是平面應力問題,σz=0,如果是平面應變問題，εz=0,所以無論對任何平面問題，單位體積的應變能都為：(a)(b)因此，物體總的應變能為：利用物理方程可得到用應變表示的單位體積的應變能為(c)(e)(e)式分別對三個應變分量求導,得到上式說明，彈性體每單位體積的應變能，對于任一形變分量的微分就等于相應的應力分量，將幾何方程代入方程(e),得用位移分量表示的單位比能的形式：(f)(5.1)又由方程(c),得在以上式子中，將Ｅ變為，將就得平面應變問題的相應方程．二、外力勢能（potentialenergyofexternalforce)

若彈性體受外力作用，平面區域Ａ內的體力分量為Ｘ，Ｙ，邊界sc上的面力分量為，則外力所做的功稱為外力功(workofexternalforce)：(5.2)由于外力做了功，消耗了外力勢能，因此，在發生實際位移時，彈性體的外力勢能是三、變分方程(variation

equation)

１、虛位移(virtualdisplacements)的概念設有平面問題中的任一單位厚度的彈性體，在一定外力作用下處于平衡狀態．命u,v為實際的位移分量，它們滿足平衡方程和邊界條件．現假想這些位移分量發生了邊界條件所容許的微小位移，即所謂的虛位移或位移變分，則有

2、位移變分方程(equationofdisplacements

variation)

由于位移的變分，引起外力功的變分δＷ和外力勢能的變分δＶ為由于位移的變分，引起應變的變分為從而引起形變勢能的變分為根據能量守恒定理，形變勢能的增加應該等于外力勢能的減少，也就等于外力所做的功，故有因而得這就是位移變分方程,又稱為位移變分Lagrange’s方程３、虛功方程(virtualworkequation)

(5.3)將(5.1)代進上式,得將上式代入位移變分方程(5.3)中,得(5.4)(5.4)式就叫做虛功方程(virtualworkequation).它表示：如果在虛位移發生之前，彈性體處于平衡狀態，那么，在虛位移過程中，外力在虛位移上所做的功就等于應力在虛應變上所做的虛功．§5.2基本量和基本方程的矩陣表示(MatrixRepresentationsofBasicQuantitiesandEquations)

在有限單元法中,為了簡潔清晰地表示各個基本量以及它們之間的關系,也為了編寫程序以便利用電子計算機進行計算,廣泛地采用矩陣表示和矩陣計算.1體力(bodyforces)的矩陣表示:設物體受體力P的作用,它的兩個分量為X和Y,則有:

２

面力(faceforces)的矩陣表示:設物體受面力的作用,它的兩個分量為

則有:(5.5)(5.6)３

應力(stress)的矩陣表示:４

應變(strian)的矩陣表示:５

位移(displacement)的矩陣表示:(5.7)(5.8)(5.9)６幾何方程(geometricalequations)的矩陣表示:７

物理方程(physicalequations)的矩陣表示:平面應力情況下簡寫的形式：(5.10)(5.11)(5.12)這里［Ｄ］稱為彈性矩陣(elasticmatrix)對于平面應變問題，可將上式中的Ｅ變為，即得平面應變問題的彈性矩陣（８）虛功方程(virtualworkequation)的矩陣表示虛功方程(5.4)可以表示為：(5.13)(5.1４)這里，在finiteelementmethod中,作用于彈性體的分布力（distributedforces）

一般用作用在某些點等價的集中力（equivalentconcentratedforces）來代替，如

Fig5.1所示．假定作用在

ipoint的集中力是

FxiandFyi

，假定作用在

jpoint的集中力

是

FxjandFyj

respectively,etc.這些分量和相應的虛位移分別用以下列陣表示：Fig.5.1(5.17)

(5.1８)

于是各外力在虛位移的虛功為：

于是由（５.14),得這就是在集中力作用下的虛功方程的矩陣表示形式(5.19)§5.3有限元導言

(IntroductionofFiniteElementMethod)

有限單元法出現于40年代，被應用于飛機結構分析，有限元這個術語是1956年Turner首先使用的。一般來說，任何能用微分方程描述的物理現象，都能夠通過變分原理建立的有限元方法來模擬，所以有限元的應用非常廣泛。不但可以解決像荷載—位移問題，還可以解決滲流、熱傳導等問題。有限元法的基本思想是用分段逼近(piecewiseapproximations)的方法,如利用多項式(polinomials)來代替連續函數(continuousfunctions),即是將連續體無限個節點變成有限個足夠多的節點，通過建立節點力與節點位移之間的關系，求解節點位移，又通過給出一定的假設建立單元內任一點位移與節點位移之間的關系，計算出應變，運用應力—應變關系，從應變就可以計算出應力。Fig.5.2

在Fig.5.1(a)中，表示一個承受初始應力px,py,pxy

的一個地下洞室(undergroundopening)；在Fig.5.1(b)中，表示將地下洞室周圍區域劃分為很多三角形單元(triangularelements)；在Fig.5.1(c)中,表示節點(nodes)為I,j,m

的一個三角形單元．有限單元法的基本解題步驟為：1.劃分單元；2.建立位移模型，即建立單元內任一點位移與節點位移之間的關系

設三角形單元三個節點的位移分別為：(ui,vi),(uj,vj),(um,vm),則三角形單元任何一點的位移u與v為：這樣，這里［Ｎ］稱為插值函數矩陣(amatrixofinterpolationfunctions,{Ue}稱為節點位移列陣(acolumnvectorlistingthenodalDisplacements)或(5.20)3.推導單元的剛度矩陣方程。既建立節點力與節點位移之間的關系．如果知道任一點的位移，則其應變為：或(5.21)將(5.20)式代入(5.21)式，得：(5.22)這里，(5.23)[B]稱為應變與節點位移關系矩陣(therelationshipmatrixbetweenstrainsandnodaldisplacements.而應力由下式確定：(5.2４)[Ｄ]稱為彈性矩陣（theelasticitymatrix）.在有限單元法中，作用在單元所有的力都用等價的節點力矢量(equivalentnodalforcevector)來表示:根據虛功原理:在虛位移過程中，外力在虛位移上所做的功就等于應力在虛應變上所做的虛功．而外力的虛功為:(a)應力在虛應變上所做的虛功為：(b)根據虛功原理，這兩個功應相等，故有：

稱為單元的剛度矩陣（thestiffnessmatrixofaelement）.

我們就得到單元的剛度方程(thestiffnessequationofaelement):

如果記：

4.建立有限單元法的基本方程：將所有單元組合起來，就得到有限單元法的總體剛度方程

其中：[K]—總體剛度矩陣，

—位移列陣，各個節點的位移，

—荷載列陣，將荷載化為節點力列陣5.求解節點位移(考慮邊界條件)。6.通過位移求解應力和應變。§5.4有限元網格(FiniteElementNetwork)1.常用單元(1)線性單元(linearelements)：三角形、矩形或其他四邊形。形函數是線性函數，即單元內任一點的坐標可用單元節點坐標的線性函數來表示。(2)拋物形單元(parabolaelements)：除了角上有節點（主外節點），邊緣上也有節點（副外節點）的單元。有時在內部也有節點（內節點），直邊或曲邊均可。單元內任一點的坐標，可用一個拋物線內插法來求得。121231234圖5.3線形單元及節點位置2.單元劃分要注意的幾個問題

(1)相鄰兩個單元的節點要與節點重合（外節點與外節點、內節點與內節點），不能與無節點邊重合。

(2)單元不必是相同尺寸，應力有突變的地方，單元劃分應較小。

(3)任何一個單元必須只能在一種材料區，即它不能跨越兩種材料區。

(4)同一單元的各個邊長，一般不要相差太大。5216436312456783124578125463圖5.4拋物線形單元及節點位置§5.5位移插值函數DisplacementInterpolationFunction(形函數,插值函數）(ShapeFunction)1.三節點三角單元的位移插值函數

所謂位移插值函數是指表是節點位移與任一點位移之間的關系的函數,即:從上一節我們知道,三節點三角單元的位移插值函數矩陣為:現在討論如何確定[N]的各個分量,我們知道,組成[N]各分量的插值函數,當它的變量等于某個節點的值時,它應該等于該節點的位移值,這就要求:這里,等這里[I]和[0]表示單位矩陣和零矩陣.同時,一個點位移的兩個分量可通過同樣的方式插值而得,因此很明顯有:最簡單的位移插值函數為線性插值函數,即假定為(5.28)

這些系數等是待定系數,通過解下列聯立方程(simultaneousequations)確定:

從以上六個方程解出6個待定系數這里一定是的線性函數,即(5.29)因而可將(5.28)寫成節點位移函數的形式：而一定是的線性函數,即這里，

(i,j,m)(5.30)而(5.3１a)(i,j,m)(5.31c)(5.3１b)A=areaofthetriangularelement

2.四節點四邊形單元插值函數對于四節點四邊形單元，單元內任意一點的位移與節點位移的關系，也是通過形函數建立的，以下以一個四節點的正方形單元為例，說明形函數的含義。取局部坐標如下圖所示，四個節點的坐標分別是(1,1),

(1,-1),

(-1,-1),

(-1,1).（1，1）（-1，1）（-1，-1）（1，-1）1234（0，0）ξη圖5.5四邊形單元映射為面積2×2個單位的正方形單元假定四個節點分別有位移，則單元內任一點的位移可用下式表示：即：（1）則是分別與節點1，2，3，4相關的形函數，形函數Ni（i=1,2,3,4）是坐標點的函數，即

形函數有如下的性質，當把第i個節點坐標代入時：如將節點1的坐標（1，1）代入，那么就有形函數怎么得來的：（1）從觀察和直覺中得到例如：統一的式子可以寫為：

（2）用數學過程來導出形函數的表達式對于矩形單元，可假設位移是坐標的雙線性函數。

式中，是待定系數。將4個節點位移值代入上述方程得:

這里8個方程8個未知數，可以求得，并將其代入上述假定式（2），與（1）式相比較，就可以求得值。

3.坐標變換(CoordinateTransformation)

在進行每個單元的計算時往往采用局部坐標，而每個單元的局部坐標與整體坐標往往是不一致的，故需要進行坐標變換。為什么要對單元采用局部坐標，這里因為在局部坐標下很多積分計算比較容易，例如函數在正方形區域上積分，可寫為：比在一個不規則的原四方形域內的積分容易的多。

從整體坐標到局部坐標的變換，可以看成是一種映射，例如從整體坐標下任意形狀的四節點單元映射成局部坐標下的兩個單位邊長的正方形單元。我們可以通過形函數來映射，即可假定單元內點的坐標滿足：式中，

分別為整體坐標下的四節點單元四個節點坐標，

代表形函數，它是局部坐標的函數，這就建立了整體坐標與的變換式，上式可寫為：

或：以下，我們介紹等參數單元的概念：對于位移，我們已有

如果對于單元內任一點的坐標也存在如下關系：即位移的插值公式與坐標的插值公式具有相同的形狀（位移插值函數=坐標插值函數），則這樣的單元就叫做等參數單元。我們一般都采用等參數單元.由數學我們知道，對于無限小區域dA有，

雅可比矩陣的元素很容易通過下式求得：式中

叫做雅可比矩陣的行列式。這樣整體坐標下的積分就可以變成局部坐標下的積分：在求剛度矩陣中，我們還需要求形函數對笛卡兒系坐標的導數，如求：寫成矩陣形式：⒋應變位移關系(Therelationshipbetweenthestrainanddisplacement)①.平面應力或平面應變：②.軸對稱情況：r為徑向，y為對稱軸§

5.6等效節點力（EquivalentNodalForces)

在有限單元法中，所有的力，無論是載荷重力、溫度應力、初應力還是由初應變、孔隙壓力、離心力、慣性力等引起的外力都要轉換到作用于節點上的等效力，以下我們來具體分析。5.6.1集中力引起的等效節點力設單元邊上作用集中力R,則利用力矩平衡條件可求得:

圖5.6集中力的等效節點力(5.32)5.6.2面力的等效節點力(Theequivalentnodalforcescausingbydistributedforces)

由作用在單元邊界或表面上的面力引起的等效節點力(這是由虛功方程推得的,證明從略)

如果知道表面正應力σ和剪應力τ,則其等效節點力為:θ為單元的荷載邊界與總體坐標x軸的夾角。式中,T是轉換矩陣，(5.33)5.6.3單元體力引起的等效節點力（Theequivalentnodalforcescausingbybodyforces

）

根據虛功原理,可求出單元體力引起的等效節點力.設體力為則等效節點力為:(5.34)

5.6.4初應變引起的等效節點力（Theequivalentnodalforce仿

causingbyinitialstrains

）初應變通常當作零來處理，但當我們試圖模擬洞室開挖、巖石加固、灌漿等過程，在這些過程中的各個階段，其前一階段產生的應變就被當作初應變，初應變也可能是外力以外的其他因素引起的，如蠕變、滲流。

根據虛功原理可以推導出初應變引起的等效節點力

若某種原因，初應變產生變化，則采用增量公式來計算由產生初應變的變化引起的節點力增量。(5.35)5.6.5初應力（或初應力增量）引起的等效節點力（節點力增量）（Theequivalentnodalforcescausingbyinitialstresses

）根據虛功原理可以推導出初應力引起的等效節點力或寫成增量的形式：

(5.36)5.6.6孔隙壓力引起的等效節點力5.6.7溫度應力在受到約束的巖體中，溫度的改變會引起應力和應變。如果巖土的熱膨脹系數為，溫度的改變產生的線應變為單元由于溫度改變而引起的等效節點力以類似于初應變方式處理。

與初應力的處理方法一樣。

這里

為孔隙壓力的增量，，，

§5.7單元剛度矩陣(StiffnessMatrixofanElement)

假如設插值函數為線性,將

Eq.(5.30)代入

Eq.(5.23),我們就得到三節點三角形單元應變與節點位移的關系矩陣[B](relationshipmatrixbetweenstrainsandnodaldisplacementsforthethreenodestriangularelement),(5.37)這里A為三角形的面積,bi,bj

等由(5.3１)確定．將(5.13)和(5.37)代如方程(5.25)，就得(5.38)這里

[krs]isamatrix(r=i,j,k;s=i,j,k),對于

planestress問題:

(5.39)對于planestrainproblems,彈性矩陣(elasticitymatrix)變成:(5.40)所以在planestrain情況下,三角單元的剛度矩陣(stiffnessmatrix)變為:

(5.41)

§5.8有限單元法的總體剛度矩陣

(

GlobalStiffnessEquationofFiniteElementMethod)有限單元法包含著對每個單元都要計算它的剛度矩陣(stiffnessmatrix[ke])和等效節點力vector{Fe}.對于單元剛度矩陣，簡單的單元可以直接積分解出；對于高階單元，還需要數值積分。將這些剛度“疊加”在一起，就可以得到總體結構的剛度矩陣。相加時是把對應點的對應項相加。設整個結構的節點數為