有限元分析教案

第一章有限元法概述在機械設計中,人們常常運用材料力學、結(jié)構力學等理論知識分析機械零構件的強度、剛度和穩(wěn)定性問題.但對一些復雜的零構件,這種分析常常就必須對其受力狀態(tài)和邊界條件進行簡化.否那么力學分析將無法進行.但這種簡化的處理常常導致計算結(jié)果與實際相差甚遠,有時甚至失去了分析的意義.所以過去設計經(jīng)驗和類比占有較大比重.由于這個原因,人們也常常在設計中選擇較大的平安系數(shù).如此也就造成所設計的機械結(jié)構整體尺寸和重量偏大,而局部薄弱環(huán)節(jié)強度和剛度又缺乏的設計缺陷.近年來,數(shù)值計算機在工程分析上的成功運用,產(chǎn)生了一門全新、高效的工程計算分析學科一一有限元分析方法.該方法徹底改變了傳統(tǒng)工程分析中的做法.使計算精度和計算領域大大改善.§1.1有限元方法的開展歷史、現(xiàn)狀和將來一,歷史有限元法的起源應追溯到上世紀40年代〔20世紀40年代〕.1943年R.Courant從數(shù)學的角度提出了有限元法的根本觀點.50年代中期在對飛機結(jié)構的分析中,誕生了結(jié)構分析的矩陣方法.1960年R.W.Clough在分析彈性力學平面問題時引入了"FiniteElementMethod"這一術語,從而標志著有限元法的思想在力學分析中的廣泛推廣.60、70年代計算機技術的開展,極大地促進了有限元法的開展.具體表現(xiàn)在：1〕由彈性力學的平面問題擴展到空間、板殼問題.2〕由靜力平衡問題一一穩(wěn)定性和動力學分析問題.3〕由彈性問題一一彈塑性、粘彈性等問題.二,現(xiàn)狀現(xiàn)在有限元分析法的應用領域已經(jīng)由開始時的固體力學,擴展到流體力學、傳熱學和電磁力學等多個傳統(tǒng)的領域.已經(jīng)形成了一種非常成熟的數(shù)值分析計算方法.大型的商業(yè)化有限元分析軟件也是層出不窮,如：SAP系列的代表SAP2000〔StructureAnalysisProgram〕美國安世軟件公司的ANSYS大型綜合有限元分析軟件美國航天航空局的NASTRAN系列軟件除此以外,還有MASTER、ALGO、ABIQUES、ADINA、COSMOS等.三,將來有限元的開展方向最終將和CAD的開展相結(jié)合.運用“四個化〞可以概括其今后的開展趨勢.那就是：可視化、集成化、自動化和網(wǎng)絡化.限元法的特點機械零構件的受力分析方法總體說來分為解析法和數(shù)值法兩大類.如大家學過的材料力學、結(jié)構力學等就是經(jīng)典的解析力學分析方法.在這些解析力學方法中,彈性力學的分析方法在數(shù)學理論上是最為嚴謹?shù)囊环N分析方法.其解題思路是：從靜力、幾何和物理三個方面綜合考慮,建立描述彈性體的平衡、應力、應變和位移三者之間的微分方程,然后考慮邊界條件,從而求出微分方程的解析解.其最大的有點就是,嚴密精確.缺點就是微分方程的求解困難,很多情況下,無法求解.數(shù)值方法是一種近似的計算方法.具體又分為“有限差分法〞和“有限元法〞.“有限差分法〞是將得到的微分方程離散成近似的差分方程.通過對一系列離散的差分方程求解,得到最終的力學問題近似解.其優(yōu)點就是：計算簡單收斂性好.缺點是：計算程序無法標準化,在不能獲得整個問題的微分方程時,該方法不能運用.由于其是將微分方程轉(zhuǎn)為差分方程,所以它是一種數(shù)學近似.“有限元法〞的根本思想就是“先分后合〞或者“化整為零,又積零為整〞.與有限差分不同,它是在力學模型上進行近似處理,也就是〔分塊近似〕.具體做法：把連續(xù)體模型轉(zhuǎn)為由有限個單元組成的離散體模型,離散體模型之間通過一些節(jié)點聯(lián)系.對于每一個離散體個體選擇簡單的函數(shù)近似表示其中的物理變化規(guī)律〔如位移等〕,運用力學方程推導單元的平衡方程組,然后集合所有的方程組形成表征整體結(jié)構的方程組,引入邊界條件,求取最后問題的解.優(yōu)點：概念清楚、易于學習理解,適用性強,便于電算化.缺點：計算精度受單元劃分的影響較大.限元分析的一般過程為了能夠了解有限元分析的全貌,我們就一個簡單的例子,來分析一下有限元分析的三個過程：結(jié)構離散化、單元分析、整體分析.一,結(jié)構離散化在該階段中,要完成把連續(xù)結(jié)構的力學模型轉(zhuǎn)變?yōu)殡x散的力學模型.處理的好壞,直接影響到最后分析結(jié)果的正確與否、計算的精度和計算的效率.根據(jù)模型的傳力特性和分析的目標,正確選擇單元類型.通常單元分為：一維單元、二維單元和三維單元.所謂一維單元就是指所求物理量僅隨一個坐標變量而變化的單元.如桁架、平面剛架和空間剛架單元.一維單元：桿單元、梁單元.二維單元：三角形單元、四邊形單元〔平面類問題〕三維單元：四面體單元、六面體單元等〔空間問題〕計算精度和計算效率：取決于單元劃分的形狀、大小和分布狀況.通常單元愈多、愈密集,計算精度愈高,但計算效率愈低.有限元分析工作就是要在精度和效率兩者之間做到有機的統(tǒng)一.二,單元分析進行單元分析的目的是為了到處表征單元力學特性的“單元剛度矩陣〞.一般說來該過程有三種方法：1,直接法.2,虛功原理法〔變分法〕.3,加權余數(shù)法.直接法概念淺顯,易于理解物理含義.變分法需要泛函的數(shù)學知識,其推導過程具有嚴謹?shù)臄?shù)學概念.加權余數(shù)法適用于泛函不存在的應用范圍.本教材將運用虛功原理方程結(jié)合彈性力學和材料力學中的知識來推求幾種常見單元的單剛計算公式.現(xiàn)在先看一個簡單的階梯軸的軸向拉伸問題例：如下圖的變截面直桿,受拉力P,運用有限元方法分析其變形.取任意單元,長度為1,面積為A,單元內(nèi)任意一點的軸向位移和其位置坐標成正比,即

u=a0+aix其中a.、a1為待定系數(shù).由于桿的兩個端點節(jié)點1、2是單元上的點,所以它們應該滿足上述方程.節(jié)點1,x1=0,u1=ao+a1X0=a.節(jié)點2,x2=i,u2=ao+a1xla1=〔u2-u1〕/l將求出的結(jié)果帶入方程并整理,就得：U=Uix=U=Uix=1--u1+-u2

<l)l-NiN2PUiU2=NU式中：m、也是形函數(shù)[N]形函數(shù)矩陣{8}e節(jié)點位移向量由位移與應變的關系知道：udu-ududxdxdx將上面推出的位移表達式代入,可得:udu-udu-u_du_￡N；?一dxdxdx=B?：、/e—x——1--dxILl上式中的[B]稱為應變矩陣或幾何矩陣.運用材料力學中的虎克定律,可以將應變和應力聯(lián)系起來.單向應力狀態(tài)的虎克定律為EVEVu1i」J2二SL.：,e[ke=B[ke=BTDIBdvv[D]矩陣是彈性矩陣.對于一維單元來說,就是所以我們這兒討論的例題：-11E.1mEA11Adx「111其中[S]稱為應力矩陣.利用虛功方程可以建立力與位移之間的關系,也就是單剛方程.在后面我們將會推導出它的一般形式如下：式中：{F}e為單元節(jié)點力向量,對我們這個例子應為[UiU2]T.[K]e為單元剛度矩陣.后面將推導出它的計算公式為求得單剛矩陣,也就完成了單元分析.總結(jié)單剛矩陣推導的步驟,應該分為四步：1〕假定單元內(nèi)位移變化的近似規(guī)律,即選擇位移模式.2〕運用幾何關系,推求位移與應變的關系.

3〕應用物理規(guī)律,把應變與應力聯(lián)系起來.4〕運用虛功方程的力與位移關系,求出單剛矩陣.單元分析是整個有限元分析的核心.不同的單元由于其力學特性不同,而具有不同的單元剛度矩陣,我們這本教材就是要學習幾種常用單元矩陣的推導和計算.了解各種單元的力學特性,為以后選擇單元類型打好根底.三,整體分析1,由各單元剛度矩陣組集成整個結(jié)構的總剛度矩陣.ki=aki=a：1tli1-11一1-1il2[-11一整個結(jié)構有三個節(jié)點,首先將單元剛度矩陣擴充為3X3的矩陣,移動各元素使之與單剛矩陣中的元素位置相對應,如下:11-10K1=EA1-110I1000_然后直接相加.000k1=受01—1I2-0-11_EAEAI1—生EAEAI1—生I1I1EA.EA2I1l2EA2120EA2I2ea2I22,把各單元的節(jié)點力向量組集成總的節(jié)點載荷向量.riR=<0,03,根據(jù)邊界條件,修改總剛度矩陣,獲得總剛方程組.邊界條件修改之前的總剛方程：EA1EA2I1EA1EA2I1I2EA2I20EA2I2EA2U1"2?

lu3,I2修改以后〔采用置“0,1〞法〕EAiiEAiiiEAi110EA尸EAi於7T0|Ui0<U211"J4,求解方程組,得出總的節(jié)點位移向量.解得的解是：Ui—Pll*Pl2eA1eAT.fUi—ea2J3,0有了節(jié)點位移,再回代到前面單元推導過程中的公式X和XX,就可以求得每個單元的應變和應力了.從這個簡單的例子,我們了解了有限元法求解力學問題三大步驟中的內(nèi)容,想必很多同學會說,這樣復雜,如果運用材料力學的知識,我還來得快些.但是大家不要忘記,有限元的計算很多都是編程完成,而且現(xiàn)在很多的商業(yè)軟件都已經(jīng)完成了很多的工作.我們學習有限元主要是了解它的原理,并對常見單元的力學特性有所了解,這樣對于以后運用有限元起到幫助作用.所以下面章節(jié)的內(nèi)容,就是圍繞這個主題展開.要到達這個目的,我們還必須學習必要的彈性力學知識.對彈性力學知識的學習,也對我們以后把握問題的本質(zhì)有幫助.第二章平面問題平面問題在力學研究的課題中屬彈性力學的范疇.該類問題不僅本身具有典型性,而且在機械零構件的分析中,也是應用得非常廣泛.所以這類問題也稱之為經(jīng)典的力學問題.我們知道,實際的機械零構件都是具有三維空間尺寸的物體,理應作為三維對象處理,但是當物體的幾何形狀和受力狀態(tài)處于某些特定的情況下,近似地簡化為平面問題,不僅可以大大簡化計算的工作量,而且其精度也完全能夠滿足所要求.如：直齒圓柱齒輪可在垂直與孔軸線的截平面內(nèi)作平面應力分析就足以了解整個齒輪的受力狀態(tài)；大壩的橫斷面可作平面應變分析來了解整個大壩受力情況等.本章是全書的重點,在這里不僅介紹彈性力學的根本知識.還將系統(tǒng)地講解有限元的基本概念、原理和方法.是學習以后各章節(jié)的根底.§2.1外力、應力、應變和位移外力、應力、應變和位移的概念在材料力學中已經(jīng)學習過,由于這些概念在彈性力學、有限元法中具有和在材料力學中不同的規(guī)定,彈性力學中的規(guī)定和有限元法是完全相同的,所以在這里我們將根據(jù)彈性力學的習慣表達方法把他們集中的加以闡述.一,外力外界作用在物體上的作用力,可以分為兩大類：1〕體積力分布在物體體積內(nèi)的力.如：重力、慣性力和磁性力等.單位體積的體力在坐標軸上的分量X、Y、Z,稱為體力分量.符號規(guī)定為沿坐標軸正向的為正,沿負向的為負.2〕面力作用在物體外表的力.如輪壓、水壓等.面力在坐標軸上的投影,表示為X、Y面力在坐標軸上的投影,表示為X、Y、Z.符號沿正彈性體受到外力作用后,內(nèi)部產(chǎn)生的反抗變形的內(nèi)軸為正,負軸為負.力.二,應力力.以彈性體中P點為定點的微單元體來考察.所謂微單元體,就是圖中PA、PB、PC的邊長分別為dx、dy和dzo以下簡稱這樣的微單元體為微元體.微元體每個面上的應力都可以分解為三個應力分量.以圖中紅面為例,分別是八、.xy、Txz0應力命名的規(guī)那么：以應力所在面垂直的坐標軸為第一個下標,應力指向為第二下標.如果下標相同就用一個下標表示.符號規(guī)定：正面上的應力與坐標軸同向為正,反之為負.負面上的應力與坐標軸反向為正.反之為負.所謂正面就是面的外法向與坐標軸同向為正.反之為負面.作用在兩個相互垂直的面上,并且垂直與該兩面交線的剪應力互等.即：6個,它們是:Txy=Pyx；Pyz=Tzy；P6個,它們是:如此以來,代表P點應力狀態(tài)的應力分量應有J'-'二y二zxyxz-yzT三,位移任一點的位置移動用u、V、w表示它在坐標軸上的三個投影分量.符號規(guī)定：沿坐標軸正向為正,反之為負.四,應變彈性體內(nèi)各點的位移在受力后一般是不相同的.各點之間距離的改變,從而使物體形狀發(fā)生變化,即所謂的變形.而物體的形狀總可以用它各局部的長度和角度表示.長度的改變稱為正應變￡,角度的改變稱為剪應變丫.以微元體三個棱邊的線伸長和角度的變化,就分別有和6個應力分量相對應的6個應變分量,即：L」；x；y；zxyxzyzT為與前面符號規(guī)定一致,這里對剪應變的符號規(guī)定如下：正應變伸長為正,縮短為負；剪應變使直角變小為正,變大為負.§2.2兩類平面問題前面我們講過,實際受力物體都是三維的空間物體,作用在其上的外力,通常也是一個空間力系,其應力分量、應變分量和位移也都是x、v、z三個變量的函數(shù).但是當所考察的物體具有某種特殊的形狀和特殊的受力狀態(tài)時,就可以簡化為平面問題處理.彈性力學中的平面問題有兩類.,平面應力問題當物體的長度與寬度尺寸,遠大于其厚度〔高度〕尺寸,并且僅受有沿厚度方向均勻分布的、在長度和寬度平面內(nèi)的力作用時,該物體就可以簡化為彈性力學中的平面應力問題.我們分析以下其應力特征.當2=±1/2時,有〔Tz=0、Pzx=0、Tzy=0o由于板較薄〔相對于長度和寬度尺寸〕,外力沿板厚又是均勻分布的,根據(jù)應力應連續(xù)的假定〔彈性力學中的根本假定〕,所以可以認為,整個板的各點均有bz=0、°zx=0、Tzy=0o如此以來,描述空間問題的6個應力分量也就變?yōu)榱?個,即七:'-Lx二yxyT而且這些應力分量僅是x、y兩個變量的函數(shù).,平面應變問題和平面應力相比擬,平面應變是當物體是一個很長、很長的柱形體,其橫截面沿長度方向保持不變,物體承受平行于橫截面且沿長度方向均勻分布的力時,該問題就可以簡和平面應力相比擬,平面應變是分析其應力特征.假定其長度方向為無限長,

那么任一橫截面都可以看作是物體的對稱面,如此那么有該面上的點都有w=0,也就是橫截面上的所有點都不會發(fā)生Z方向的位移.由這一點可以推出也就有sz=0、Tzx=0、Tzy=0o￡z=0,那么是否也就有bz=0呢可能有同學想d=E￡,當然也就有bz=0,這是錯誤的.平面應變狀態(tài)下bzW0的.雖然不等于零,但它也不是一個獨立的變量了,它由

by的大小而決定.如此以來,獨立的應力分量同平面應力問題一樣也是3個:三,兩類問題的比擬1,幾何特征平面應力厚度〈V長度、寬度平面應變厚度＞＞長度、寬度為便于說明可講上述長度看作為厚度2,受力特征外力都必須在其面內(nèi)且不沿厚度方向變化3,應力特征平面應力6r0、Tzx=0、Tzy=0.&產(chǎn)0自由變形〔無約束〕平面應變,3zW0但不是自變量、Tzx=0、7rzy=0o￡z=0總以上比擬可以看出,平面應力是真正的2維〔平面〕應力狀態(tài),而平面應變卻不是,而是3維應力狀態(tài),只不過bz不是獨立變量而是隨橫截面平面應力分量而定.獨立變化的應力分量只有3個,類似于平面應力狀態(tài).衡微分方程彈性力學求解問題是從靜力學、幾何學和物理學三方面綜合考慮的.所以我們首先微元體應該滿足的平衡條件一一平衡微分方程.我們以平面問題為例推導,看看它應該具有什么形式.首先對平面問題的微元體進行受力分析圖,如左所示.物體靜力平衡的條件是：匯Fx〔y〕=0;匯M=0.先看EFx=0「x,±F1dxdy1-「xdy1-%,二Tyxdydx1-yxdx1,Xdydx1=0展開化簡得同理可求得匯Fy=0滿足得條件,色〕+色〕我=0N＞x由匯M=0,列出方程如下：卜熊町1喈十山咚-1管尸1巖-…巖=.化簡后得:xy2fxdxyxy2fxdxyx2：:ydy略去微量項,可得：“y=3這就是前面所將的剪應力互等.對于平面應變問題,微元體的前后面還有正應力bz,不過它們是互等的.對于推導出來的結(jié)果,沒有任何影響.所以平面問題的平衡微分方程就是：三.工.丫二0X=0寫成矩陣形式何方程考察平衡微分方程,其中具有三個未知變量bx、by、Pxy,,而只有兩個方程,方程具有無數(shù)個解.說明僅從靜力學關系無法求解該方程.我們必須從其它方面尋求幫助.彈性體在受到外力后,會發(fā)生位移和形變,從幾何上描述彈性體各點位移于應變之間的關系,就是彈性力學中的又一個重要方程一一幾何方程.1坐標軸x的變化率,1坐標軸x的變化率,偏導數(shù)—：:y移為u、..,v,那么B點的位移就是:二uuR'=u——dxBex..'同理的D點的位移分量vB'二v=v—dx-:xuD'=uFu—dy

:y：:vvd''v—dy二y由于a角在位移和形變很微小的情況下非常小,所以AB'^AB線段AB位移后的總伸長量為仍然取截面的微元體ABCD,AB、CD邊長為dx、dy,厚度為“1〞.位移u、v都是x、y的函數(shù),即u(x,y)、v(x,y),偏導數(shù)出、0表示位移分量u、v沿x:x￡v表示位移分量u、v沿坐標軸y的變化率,設A點的位-y

一,一一一〞一',cu,cu,AB-AB=AB-AB=ub-ua=u+—dx-u=—dx:x=—dy/dy=—：:y；:y；x=-udx/dx=-u=—dy/dy=—：:y；:y剪應變由a、3兩個角度組成：：t：：tgABV-xdx-vdxu—dx-u

一：x由于里Y1,所以口=生,同理可得日=包一x二x二y:v：:u…xv=■=.xy--二x二y0Li.Cy綜合以上幾何方程,并將它們寫成0Li.Cyyyxy由以上方程可以看出,當彈性體的位移分量確定以后,由幾何方程可以完全確定應變,反過來,應變卻不能完全確定彈性體的位移.這是由于物體產(chǎn)生位移的原因有兩點：1〕變形產(chǎn)生的位移.2〕因運動產(chǎn)生的位移.因此彈性體有位移不一定有應變,有應變就一定有位移.理方程描述彈性體內(nèi)應力與應變關系的方程,我們稱之為物理方程,也叫材料的本構方程.彈性力學通常研究的是各向同性材料,在三維應力狀態(tài)下的應力應變關系.當彈性體處于小變形條件下,正應力只會引起微元體各棱邊的伸長或縮短,而不會影響棱邊之間角度的變化,剪應力只會引起角度的變化而不會引起各棱邊的伸長或縮短.因此運用力的疊加原理、單向虎克定律和材料的橫向效應〔泊松效應〕,我們就可以很容易的推導出材料在三向應力狀態(tài)下的虎克定律,也就是通常所說的廣義虎克定律.式中,式中,EG材料線彈性模量-材料剪切彈性模量-材料橫向收縮系數(shù),即泊松系數(shù).三者不是獨立的.具有以下關系:

這些參數(shù)都是材料的固有屬性系數(shù),可以通過查材料手冊獲得.如：鋼材的彈性模量=196?206GPa之間,通常取2.1x105MPa,科=0.24?0.28之間,也取為0.3進行計算,=79GPa.將以上空間問題的物理方程運用到平面問題,其形式如下：Zz=0、Tzx=0、Tzy=0.所以:1〕平面應力Zz=0、Tzx=0、Tzy=0.所以:1二-二x一二01=之x1二-七；x.」二yE121匕；zyTOC\o"1-5"\h\z=—T=TzycxyxyGE從以上物理方程,也論證了我們前面說的￡zw0的結(jié)論,但由于它是由x和y方向應力產(chǎn)生的附加無約束變形,所以通常不予以考慮.在有限元分析中更多的是運用應變表示的應力關系,所以我們將上式變形一下：xyxy二E21」xy以上方程的矩陣表達形式為:ax,ayax,ayyxyE1-」201-口2yxy簡記為：七〕-D,'.；〕式中:｛b｝、｛￡｝為該問題的應力、應變向式中:[D]為彈性矩陣.它是一個對稱矩陣,且只與材料的彈性常數(shù)有關.2〕平面應變問題的物理方程由于&z=0所以由空間物理方程的第三式得:〕,代入〔1〕、〔2〕式得一2CT1-」y11-」2P二一x1-」zy1TxyG21日一■xyE-xy同理變型為應變表示應力的形式：十—十——1-J-yE1-」22xyxy+z1」x1-UE1-j22E1-」2xy21+——1-U矩陣形式:□xVxy；x無\y轉(zhuǎn)轉(zhuǎn)換到平面應幾何方程和物的.也可簡記為七D「二比擬可以發(fā)現(xiàn)只需將平面應平面應變問題的彈性矩陣不同于平面應力問題的彈性矩陣,力問題彈性矩陣［D］中的材料常數(shù)E換為E/〔1-^2〕,科換為科/〔1-科〕比擬可以發(fā)現(xiàn)只需將平面應的彈性矩陣.其實彈性矩陣的這種轉(zhuǎn)換方法,是彈性力學中將平面應力結(jié)果,變問題結(jié)論的一般方法.由于在兩種平面問題的描述方程中〔平衡微分方程、理方程〕,只有物理方程是不同界條件求解彈性力學問題實際就c是在確定邊界面問題而言〕,學的角度看,就融T條件下,求解8個根本c是在確定邊界面問題而言〕,學的角度看,就融T是求解偏微分方程的邊值問題.邊界條件的給出通常是各式各樣的.大體可以分為三類:1,第一類邊值問題給定物體的體力和面力條件,確定彈性體的應力場和位移場.此類問題邊界以力的形式給出,所以也稱為應力邊界條件.我們可以來考察一下應力邊界的一般形式：ajivj=Ti五是在S,面上給出的力的分量.平面問題如左圖所示,設陰影局部的微元體弧長為ds,厚度為單元厚度“1〞,其法線與X軸的夾角為由陰影局部微元體的平衡條件可以推出：Xds父1-crx?dscos6x1—e,dssine父1=0"人xyYdsx1-try,dssin6x1-zxy,dscos日x1=0'化簡后得：◎xcos9+ixysin日=XJcy、_>此即為平面問題應力邊界方程.仃ysine+%ycos@=Y2,第二類邊值問題給出彈性體的體力和物體外表各點得位移條件,確定彈性體得應力場和位移場.由于以位移給出得邊界條件,所以也稱為位移邊界問題.一般得位移邊界條件為：5=也在Su面上.3,第三類邊值問題給定彈性體得體力和一定邊界上得面力,其余邊界上的位移,確定其應力場和位移場.由于邊界以力和位移兩種形式給出,所以也稱為混合邊界問題.針對不同的邊界條件,彈性力學求解的方法也有所不同.彈性力學的解題方法〔解析法〕1,應力法由于第一類問題的邊界條件以應力形式給出,所以以應力作為根本的未知量的求解過程,就是人們通常所說的應力法.由于平衡方程中有三個未知量,而只有兩個平衡微分方程,必須找出另外一個包含應力分量的方程,才能求得方程解.考慮到彈性體變形前是一個連續(xù)體,變形后也應是連續(xù)體的根本假設,所以要求微元體的變形一定要協(xié)調(diào),才能使變形前、后,不會發(fā)生裂縫、重疊等現(xiàn)象.要使變形協(xié)調(diào),就要研究幾何方程.前面介紹的平面問題幾何方程如下：分別對ex、ey、求y、x的二階偏導,然后相加:fV二2fV二2xy十—I=.x:x:y上式說明三個應變分量之間應滿足的連續(xù)性條件,我們稱之為形變協(xié)調(diào)方程〔相容方程〕通過物理方程,使上述的形變協(xié)調(diào)方程換成應力表示的形式,使之與平衡微分方程就構成了應力法中需要求解的方程組.具體我們來看看：1〕利用物理方程消去相容方程中的形變分量〔以平面應力為例〕1f2.1.E守〞…y妥、…x2(1+N)3、xyEj1f2.1.E守〞…y妥、…x2(1+N)3、xyEjxjy_2_2「*-『『「yx=21?口：：2xy.:x；:y2〕利用平衡微分方程,消去上述公式中的剪應力xy-:y二二x-X_:xxxy:x/y—y-Y-:y二22"y；xyX求偏導,C2式對y求偏導,然后兩者相加9.2_產(chǎn)二x:X二二y三Y———21f~~—x:x二y二y代入相容方程,化簡二2——(er-kkr.2xy-y-2-

exd2OxdX32OyY.2一.2.x二x二y■:y-2_二:-x-2x.2C+—y

2~■：2c：__y-2x(仃y-N仃X)=I1+N)-1二十——ySo*ySo*y?二y11?日Y+——對于平對于平面應變而言,就可以了.運用前面講過的物理方程的轉(zhuǎn)換方法,只需將上式中的代以(1/(1-)1+CT1+CT)=y1-」.XFY-f-:x：y3〕最終求解的方程組平面應力問題三-三X二0一一Y二0,：y二x二x.二y=-1.」平平面應變問題三個微分方程,三個未知變量,再考慮邊界條件,即可求得.如果是單連體〔只具有唯一的封閉邊界〕的對象,滿足了以上方程組后就是實際的解.但對于多連體〔具有多個封閉邊界〕的對象,中還包含有待定系數(shù),這些待定系數(shù)會導致位移的解出現(xiàn)多值性.所以對于多連體的問題,還應考慮位移的單值條件,才能最終確定.該局部的內(nèi)容可以參見徐芝綸編的?簡明彈性力學教程?中圓環(huán)受均布壓應力的情況〔P87〕2,位移法位移法主要針對第二類邊界條件問題求解.解題步驟：1〕改寫物理方程使之成為應變表示應力的形式2〕應用幾何方程表示以上得到公式中的應變3〕將它們代入平衡微分方程經(jīng)整理最后得到的位移法求解平面應力問題方程為：E1-口E1-口222.uu十c2改1-J-2u1-2v+2cy22dxdy)E(d2v上1—□d2vJ+□d2u"1-N2(W22加2fixcy；兩個未知量,兩個方程,再加以邊界條件即可求得問題的解.以上介紹的解析法中,應力法和位移法是求解彈性力學問題的根本方法.但都需要解聯(lián)立的偏微分方程組.求解過程中的數(shù)學難度,常常導致這種求解是無法進行的.由于應力法在體力為常量的情況下可以進一步簡化為求解一個單獨的微分方程的問題,所以應力法在解析法中相對應用較多.但即使這樣,在應力法中,也常常采用逆解法或半逆解法.體力情況下應力法的簡化、應力函數(shù)及實例分析我們前面講述了彈性力學的三大方程,及應用這三大方程的應力法和位移法解題步驟.但是也說了要求解這些聯(lián)立的偏微分方程在數(shù)學上是存在很大難度的.很多情況下,根本無法進行.那么彈性力學如何在實際中進行應用,它們和我們前面學過的材料力學區(qū)別究竟在哪里我們將通過這一節(jié)的學習,一方面了解如何應用這些彈性力學的方程求解問題,另一方面加深對力學概念的理解,建立力學分析問題的直觀感覺,為建立有限元模型打好根底.我們知道在大多數(shù)情況下,我們分析的對象,體力是常數(shù),它不隨x、y坐標變化.如此以來,前面講解的第三個方程〔應力表示的相容方程〕,就可以簡化為了：二2二2J+￡y[&X+Oy〕=0簡記為：▽2〔仃、+仃〕=0-2-2yxy：xn以上方程稱為拉普拉斯微分方程,數(shù)學上也稱之為調(diào)和方程,滿足調(diào)和方程的函數(shù)稱之為調(diào)和函數(shù),及這里的〔仃*+byV2是拉普拉斯算子.

這樣以來常體力情況下的應力法方程就是:二.吐X=0*：y三.三.丫一.:yY、'一二x?二y=0以上方程都不含有材料常數(shù)E、W,所以平面應力和平面應變兩類問題具有相同的方程,這說明：在單連體問題中,只要邊界相同、受同樣的分布外力,應力分布與材料無關；也與是平面應力還是平面應變的狀態(tài)無關.以上結(jié)論的意義：⑷彈性力學平面解答的應用范圍加寬.②為實驗應力分析提供了理論依據(jù)〔光彈實驗〕下面我們考察平衡方程：其解由隹〕+匡=0其次微分方程的通解,加上任意一組特解組成.三■士工0特解特解我們可以很容易找到.如：:-x=-Xx二y=-Yyxy=0所以所以現(xiàn)在關鍵是找其次方程的通解.由第一個方程,可得：尚、2G已/1T"萬=ly…分的充要條件.所謂全微分就是有一個函數(shù)由數(shù)學微分理論,該式是一個函數(shù)全微dA=：xdy+JdA=：xdy+J〞xydx且:二：:A-7一■xy：:Aex同理由第二式可得：;〕::Bexxy:By由剪由剪應力公式又知存在一個函數(shù)可以使d=BdyAdxB=—故::x由于應力與函數(shù).存在這樣的關系,因此函數(shù).即是應力函數(shù).我們用應力函數(shù)來表示相容方程：'J';=\4=0上式說明.是重調(diào)和函數(shù).前面講過在彈性力學中,常用逆解法和半逆解法.所謂逆解法就是設定各種滿足相容方程的應力函數(shù),運用bx、by與.的關系,求得應力分量,再考察其滿足何種邊界條件,從而知曉這樣的應力函數(shù)可以解決什么問題.所謂半逆解法就是根據(jù)彈性體的邊界形狀和受力關系,設定局部應力分量為何種形式的函數(shù),從而確定應力函數(shù).,在運用應力函數(shù)求出所有的應力分量,根據(jù)邊界條件確定應力分量應具有的最終形式.下面我們來看一個半逆解法的例子.運用逆解法求簡支梁受均布載荷的應力分布.II由材料力學知,彎曲應力主要由彎矩M引起,剪應力由剪力引起,而擠壓應力由分布載荷q引起.現(xiàn)在q為不隨x變化的常量.因此我們設6y不隨x坐標變化,即仃y=f(y),c2(P,、因此%=fyex我們對x積分:J=xf(y)+f1(y)；:x2=Yfyxfiyf2y上式中fi(y)、f2(y)是待定函數(shù).由于應力函數(shù)必須滿足相容方程,所以:.2..2.2二二4.x：14：d2fy一2一2一,2x二ydy代入到式O中獨二/d4f(y)+d4fi(y)+d4f2(y)：:y42dy4dy4d代入到式O中4_.442_dfyx2xdfiydf2y2dfy=0dy4dy4dy4dy2考察上式可以看出它是一個x的二次方程,所以一般情況下只有兩個根.也就是說只有兩個位置能夠滿足上式.但我們對相容方程的要求是絕對滿足.也就是要求在整個梁的范圍內(nèi)都滿足.所以只有該方程的系數(shù)項和自由項全部為零.即：d4fy_nd4fiyd4f2y,d2f2y_nTOC\o"1-5"\h\z404042.20dydydydy?1-f(y)=Ay3+By2+Cy+D②fiy=Ey3Fy2Gy@d4f2d4f2y

dy4d2fydy2--12Ay-4BA5B432f2y=)yyHyKy⑷106公式中的A.、B、C…K都是待定系數(shù).公式.3、O中分別省掉了常數(shù)項和一次項、常數(shù)項.這是由于fi(y)和f2(y)分別是應力函數(shù)中x的一次項和常數(shù)項的原因,這樣處理不會對應力分量產(chǎn)生影響.最后求出的應力函數(shù)為：2=—Ay3By2CyDxEy3Fy2Gy--y3--y4Hy3Ky2TOC\o"1-5"\h\z2106由應力與應力函數(shù)的關系,可以求出各個應力分量：?2?:x232二x=r二一6Ay2Bx6Ey2F-2Ay3-2By26Hy2Ky232二y-Ay3By2CyDxy=-x2Ay22ByC-3Ey22FyG由于以上求得的應力分量滿足了平衡方程和相容方程,所以只需根據(jù)邊界確定A…K的系數(shù),就求得了該問題的解.根據(jù)對稱性,知道為偶函數(shù),為奇函數(shù),所以有E=F=G=0通常梁的跨度遠大于梁的深度,因此上下邊界是主要邊界,它們必須滿足.二yy4=0二yy#=-q,xyy：；h=0將它們代入的表達式,并且考慮E=F=G=0—A—BhCD-0842

,3,2hAhA—32-hAhBC=043h2A-hBC=04以上四個方程解四個未知數(shù),求得:9將他們代回到9將他們代回到應力分量的表達式中,也就2A=-―3B=0C=-h32h有:■xy6qh32q26q34_■xy6qh32q26q34_xyq3-y6Hy2Kh36q2

二hyh3qq,—y--2h2四x2h左右兩個邊界,由于前面已經(jīng)考慮了對稱性,沒有水平力.要x=l時,bx=0,考察bx所以這個僅考慮優(yōu)邊界.的表達式,除非q=0.而這和條件相違背.所以在這個邊界上我們只能要求局部滿足.運用圣維南原理運用等效力系代替它O的誤差只在力作用點附近較大〕.運用的等效力系就是合成力系為平衡力系:h2〔這樣產(chǎn)生.二xx'dy=0h~2h2.二xx^ydy=0h合力等于0合力矩等于0由第一個條件得K=0,〔奇函數(shù)在對稱區(qū)間上的積分為零〕；由第二個條件得ql2

h3q

10h可以證實剪應力的合力為-ql.即h2,xyx'dy=「qlh最終求得的結(jié)果,加以整理：-x2yq-I

--xyTOC\o"1-5"\h\zhh\/由于厚度為“1〞,此時其慣性矩I=一,由于厚度為“1〞,此時其慣性矩1282任意一點的彎矩M=qll-x--l-x2=-l2-x222剪力Q--qlql-x--qx所以上式中的應力分量可以改寫為:2y2yQSxy各項應力的分布bx第一項為主應力項,與材料力學中的結(jié)果完全一致.第二項為應力修正項.當L/h>4時,僅占主項的1/60;當L/h=2時,僅占主項的1/15.所以對于深梁的工程構件不容無視第二項的影響.§2.8虛功方程從前一節(jié)深梁的例子,我們可以看到,彈性力學解析求解的過程是非常復雜的.這樣的求解對實際工程情況說來,很多情況根本是不可能的.所以長期以來,技術人員就一直探求數(shù)值求解的方法.有限元法是其中最成功的方法.為分析單元特性和簡化分析過程,我們還需了解單元的能量關系.由于在力學上,很多時候從能量的角度分析,可以大大簡化分析的過程.一,應變能的概念有材料力學知,彈性體在受到外力作用發(fā)生變形的過程中,彈性體內(nèi)部會存儲變形勢能——應變能.在單向應力場中,單位體積的應變能的計算可以表示為：…1dU=一；二；2對于平面問題,有三個應力分量和與之對應的應變分量.由于在小變形情況下,正交力系互不影響,由力的疊加原理,所以該種情況的應變能為：TOC\o"1-5"\h\z1：,T：■1：,T：■dU=-x-x1y-y'-xyxy=-?22其中："‘」x二yxyT'：=x；yxyT整個彈性體的應變能：U=11二:=二卅二1J)Ti；,dv2v2v上式也表示應力在應變上所做的總功.二,虛功原理和虛功方程在理論力學中,我們曾經(jīng)學到一個虛功原理,也稱虛位移原理.其根本思想就是：假定加在物體上一個可能的、任意的、微小的位移,在平衡條件下,物體的約束反力所做虛功應等于外力所做虛功.由于能量必須守恒.在這里所說的可能的虛位移是指位移必須滿足的約束條件.任意的是指位移的方向和類型是隨意的.把這一原理運用到現(xiàn)在的彈性體中,衡量彈性體應滿足的平衡能量關系就是：假定加在彈性體上一個可能的、任意的、微小的位移,在平衡條件下,彈性體內(nèi)的應變能應等于外力所做虛功.同樣是由于能量必須守恒.運用這一原理,我們可以推到有限元中廣泛用到的虛功方程.

假定彈性體發(fā)生u、v的虛位移,那么由幾何方程得:-**二-**二ux二Fx:y*cu二v,,,,,xy-Fy；:x現(xiàn)考察彈性體微元體和邊界處微元體上的力所做虛功:⑷內(nèi)部微元體上的力所做虛功左面的應力虛功-Jxdy11u右面的應力虛功左、右兩面上的虛功之和〔略去高右面的應力虛功左、右兩面上的虛功之和〔略去高階微量,并考慮*x二xxudv1；x同理得剪應力的虛功之和二u二,xy*xy——udvi二y二y,rr一一,,1.r*體力X的虛功Xudv1同樣的考慮Y方向的仃y以及Eyx以及體力Y的虛功,然后疊加成內(nèi)部微元體上的虛功如下：dwi:-y；式y(tǒng)x***YvdviT.,.,x；x二y；y二y二x+Txy?xydvi⑸邊界上的微元體**設斜邊中點處的虛位移u、v,形心處的應力為仃xby、Txy那么在直角邊上的應力和位移均有一個負的增量,如下圖,虛功計算為::二xdx」,—xdy1：x2,*v**dx定一戶xu-u+%%I-dyMi〔略去了局階微量〕IL「：x2同理同理的dy直角邊上的剪應力虛功為:r*一."xyui；*：xy*：:u*dyu?xydx1.n：：y2代入已有的幾何關系:dx=dssin3,dy代入已有的幾何關系:斜面上面力的虛功體積力的虛功同樣的方法求得另一面上正應力與剪應力的虛功,全部相加即得斜邊微元體上的虛功之和為：'x二.xydw2=；x7二:.'x二.xydw2=；x7,,YVdV27,1xJxxyxydV2Nfx+iX-：xcos、：-xysin-u*Y-0ysin：-.xycos-v*1支反力處的虛位移為零,所以支反力不做功,將dWi+dW2,并對整個體積積分,可以得到整個彈性體內(nèi)的總虛功：**一?X-：xcos；1〞xysin「uY—'sin；一：；xycos「vdss+1V根據(jù)平衡微分方程和靜力邊界條件,上式的第一、第二項都是零,所以彈性體的總虛功為:***Wlz—x-xy-y-xyxyv根據(jù)能量守恒,它應與外力的虛功相等Wz：jHx；x二y；y,xy；ydv=Xu*Yv*dv-Ii：Xu*Yv*dsvvs由于該等式引入了平衡方程和邊界方程,所以上式虛功方程等價于靜力平衡條件〔內(nèi)部和邊界微元體〕.不同之處在于它是一種能量的表示形式.為了便于有限元中方便運用,引入廣義力和物理方程,虛功方程變形為：；、,.*：*：F〕-〔；*書?：；加v綜合以上推到過程,虛功方程表達的物理概念就是：“假設彈性體處于平衡狀態(tài),那么外力在虛位移上作的虛功應等于應力在應變上作的虛應變功,或者說等于虛應變能〞.面問題單元劃分有限元法在有限元法在平面問題進行分析時,角形單元的優(yōu)點是簡單且對結(jié)構的不規(guī)那么邊界逼近好,部的應力應變變化.這兩點我們會逐漸向大家說明.或者矩形單元,三而矩形單元卻更能反映實際彈性體內(nèi)所以一般說來,有限元分析,單元劃分的密度和單元種類選取,對計算結(jié)果起重要作用.一般單元劃分越密集,結(jié)果越精確.單元多也導致求解的線性方程組階數(shù)增高,要求計算機的內(nèi)存也更大,計算的時間也越長,分析的效率就越低.解決才采用三角形單元和四邊形單元、這一矛盾的方法就是在應力集中區(qū)域單元劃分密集一些,應力變化梯度小的位置,劃稀疏些,這樣就能兼顧精度與效率的關系.一般的原那么是：根據(jù)結(jié)構的受力和支承特點,按對稱和反對稱的性質(zhì),簡化分析模型,以減少計算分析的規(guī)模.合理布局單元的密集程度,以使計算結(jié)果精度高而計算量小.在同一單元內(nèi),單元的特性數(shù)據(jù)和材質(zhì)數(shù)據(jù)應保持一致.集中載荷的作用點和載荷密度突變處應有節(jié)點.在欲知道應力狀態(tài)、內(nèi)力情況和位移值的位置應有節(jié)點.單元的選取欲分析的目標密切相關.模型的單元劃分好后,把所有的單元和節(jié)點按一定的規(guī)律和順序進行編號,選擇適當?shù)淖鴺讼?直角、柱面和球面),以方便確定各節(jié)點的坐標值.點位移、節(jié)點力和節(jié)點載荷彈性體在承受外力作用后,其內(nèi)力的傳遞實際是通過單元之間的邊界來實現(xiàn)的.但我們把結(jié)構離散化后,如果單元劃分得足夠小時,可以看成為其內(nèi)力的傳遞通過單元與單元之間的節(jié)點進行傳遞.對于平面問題而言,每個節(jié)點都有位移和力兩個未知量,這兩個量又都是x、y的函數(shù),注意平面問題的節(jié)點是不能傳遞力矩的,為什么一,節(jié)點位移對三節(jié)點三角形單元而言,因有三個節(jié)點,每個節(jié)點的位移都有x,y兩個分量,所以一共有6個自由度.單元節(jié)點位移向量可表示為：Je=LViUjVjUmVmT二,節(jié)點力所謂節(jié)點力,就是單元對節(jié)點或節(jié)點對單元作用的力,它是彈性體內(nèi)部的作用力,也就是我們常說的內(nèi)力.和上面相同的道理,節(jié)點力向量為：T?=UiViUjVjUmVmT■三,節(jié)點載荷■節(jié)點載荷是指作用在節(jié)點上的外力,包括：J⑷直接作用在節(jié)點上的外力②經(jīng)等效處理后,移植到節(jié)點上的等效力.可用Xi、Yi.表示.由力平衡條件知,節(jié)點要保持平衡,那么作用在節(jié)點上的載荷應等于節(jié)點—n/.內(nèi)力的合力.即:Xi=nUi,Yi=nVi,ee所有的節(jié)點載荷向量表示為：R；Fje§2.11三節(jié)點三角形單元的位移模式和形函數(shù)彈性結(jié)構受外載荷后,內(nèi)部各點的位移變化規(guī)律,一般都是很復雜的.很難找到一個函數(shù)來描述整個結(jié)構內(nèi)部各點的位移變化規(guī)律.但當把整個結(jié)構離散以后,在一個相當小的單元內(nèi)部,卻可以用簡單的函數(shù)來近似描述單元內(nèi)部位移的這種變化規(guī)律.而不致造成很大的誤差.這就好比一條復雜的曲線,用一個函數(shù)很難描述,但在把這條曲線分段以后,對于每一條分段曲線,卻可以用直線或拋物線來描述.一,三角形單元的位移模式設單元內(nèi)任意一點的位移分量u、v是其位置坐標x、y的線性函數(shù),那么：ai…a6是待定系數(shù).u=ai…a6是待定系數(shù).v=a4+a5x+a6y改寫方程組成為矩陣的形式:單元的三個頂點i、j、m的坐標,分別為〔xiyj〔Xjyj懷口〔xmym〕,由于它們也是單元上的點,所以應該滿足以上假定的位移變化規(guī)律.代入上式:UiUi=a〔a?xia3yiVi=a4a5xa6yiUjUj=a〔a?xja3yjVj=a'a§xja6yjUmUm:a-a?xia3yiVm=a4a5xma6ym解以上6個方程,求得6個待定系數(shù).aiUi為yaiUi為yiujxjyjUmxmymxiyii?xjym-xmyj2AUi'"-.'xmyi_xiymUjxiyj_xjyiUm'1xjyj1xmym12AajUjamUmUiV\UjXiXjXmV」ViVj=2A-Vj-VmUiVm-ViUjVi-VUm-biUibjUjbmUm12Aa2同理得XiUiXjUjXmXiXjViV」2Vxm-XjUiXi-XmUjXi-XiUmXm_工一2AiuiCjUjCmUma4a5a62A

12A

12AiUiiUiiUiaiajUjbjUj,CjUj=XjVm-XmYjb=Vj-VmCi=—Xj-Xm'amUmbmUmCmUm—fijm順序輪換,a是三角形的面積公在單元節(jié)點的順序號i,j,m必須是按逆時針排列,否那么系數(shù)行列式是負值,而三角形的面積為負值,是不合理的.求得的6個系數(shù)可以用以下矩陣表示:■ai0aj0am01‘Uir1bi0bj0bm0Vi忸11Ci0Cj0Cm0Uj2A0ai0aj0amVja6UJ0bi0bj0bmUm0Ci0Cj0CmLVm二,形函數(shù)

將所求得的6個待定系數(shù)代入位移模式表達式中:■LJa0aj0am0"[|Ui'b0bj0bn0Vixy000];0cj0cm0uUjI001xy0ai0aj0am1Vj0b0>0bmUm0ci0cj0cmVmJ11x012Ay000也01xya6UiV1

UiV1

2Aabxcy

00bxcy0abxciyat^xciy0——fi,j,m順序輪換.人一1,令Ni=——aibixciy2A0NiNm01：K:0Nm」lVm,ui卜人＞上式就是假定位移模式下導出的單元內(nèi)任一點位移表達式.其中的插值該式的數(shù)學意義就是單元內(nèi)任一點的位移可以由單元節(jié)點的某種形式插其中的插值基函數(shù)就是Ni、Nj、Nm.對于我們目前假定的位移模式是線性函數(shù),所以得出的插值基函數(shù)也數(shù)也是類似的線性函數(shù).由此可以看出,以也稱之為位移形態(tài)函數(shù),簡稱形函數(shù).三,形函數(shù)的性質(zhì)插值基函數(shù)具有反映單元位移變化形態(tài)的特征,[N]就是形函數(shù)矩陣.③單③單元內(nèi)任一點的三個形函數(shù)之和恒等于1.NiNjNm=1CD在單元的三個頂點處,有i節(jié)點處i節(jié)點處Ni=1Nj=0Nm=0j節(jié)點處j節(jié)點處Nj=1Ni=0Nm=0m節(jié)點處m節(jié)點處Nm=1Ni=0Nj=0以上這些,可以通過簡單的數(shù)學運算進行證實.四,位移模式收斂性的分析由于位移模式的選取是人為假定的,這種假定只能近似模擬單元內(nèi)位移的變化規(guī)律,由于單元剛度矩陣的推導是以假定的位移模式展開的,那么這種假定的位移模式能否使有限元數(shù)值解收斂與精確解,在很大程度上就取決于所選的位移模式,通過數(shù)學證實,可以找出位移模式滿足收斂性的幾個條件：A〕完備性1〕位移模式必須包含單元的常應變狀態(tài).

將u、v的表達式代入幾何方程,得：；:u；x=——=—a1a2xa3y=a2Fx；x.:v-y.:v-y=—a4a5xa6y=a6.x.:vFu一：：xy=———=—a4a5xa6y—a〔a?xa3y=a3axy:xfy：:xfy由于系數(shù)a2-a6都是常數(shù),所以上面的應變分量也是常量.這也說明所選的位移模式中包含有彈性體的常應變狀態(tài).在上面的表達式中不含x、y的變量,說明單元的應變是常量,這也說明這種單元是一種常應力單元.2〕位移模式必須包含單元的剛體位移.彈性體位移一般包含兩個局部,即變形位移和沒有彈性變形的剛體位移.那么作為模擬單元位移狀態(tài)的位移模式,也就應該能同時反映這兩局部的位移.在③中已經(jīng)證實了彈性應變的變形,下名說明他也包含有剛體位移的特征.改寫位移表達式如下：u=au=a1a2xa3y=a1v=a4a5xa6y=a4a5-a3a5a3a2xy22a5-a32當%當%=%=/y=0時,由O知,a2=a6=a3+a§=0,所以上式:u二a1v—a4a5-u二a1v—a4a5-a32a5-a3x2我們再來看看一個點作剛體位移的運動方程.點M先轉(zhuǎn)到Mi,再由Mi移到M2,如左圖所示：u=u0-r切cosa=u0-wzy

v=v0+r6sina=v0+wzx**比擬上面的*和**式,可以看出a5=u0,a4=v0,a5-a32=Wz由此得出在三角形線性位移模式中,也包含了單元的剛體位移.B〕協(xié)調(diào)性3〕位移模式必須能夠反映位移的連續(xù)性.在彈性力學求解問題時,曾經(jīng)講到過變形協(xié)調(diào)方程,也說過它是彈性體變形后仍保持連續(xù)和不發(fā)生撕裂、侵入缺陷的條件.那么位移模式的選取,也應該保證單元之間不出現(xiàn)撕裂和侵入的缺陷.由于上面假定的位移模式是線性函數(shù),而兩點就能決定一條直線.由于相鄰單元的公共邊界上的位移由與之相連的兩個節(jié)點插值獲得,而相鄰的單元具有兩個公共的節(jié)點,所以通過這兩個節(jié)點所得的插值值,不可能出現(xiàn)不同.也就是不可能出現(xiàn)下列圖~a和圖b的情況.以上三個條件是選取位移模式必須考慮的.完備形條件是收斂的必要條件；協(xié)調(diào)條件是充分條件.在有限元中,滿足完備條件的的單元是完備單元,滿足協(xié)調(diào)條件的是協(xié)調(diào)單元.〔三節(jié)點三角形〕單元剛度矩陣的推導上一節(jié)我們已經(jīng)建立了三角形單元的位移插值模式,并求得了形函數(shù)的方程,這樣就完成了單元內(nèi)任一點的位移由單元節(jié)點位移表示〔插值〕的工作,接下來運用我們已學過的一系列知識,我們就可以完成單元剛度矩陣的推導了.一,推導過程1〕由位移插值函數(shù)導出單元應變的單元節(jié)點位移表達式「￡

「￡

除0的0Ni曳|ex0Ni0NNxN曳|ex0Ni0NNxNjx0-Nj0Nj：:y.Njx:NNm：:y：Nm二y：NmUi、vmj,1,、??一將Nj=—(ai+bix+ciyM弋入上式,可得:2AUibi0Uibi0bj0bmtie10zf0——0c0cj02A!jCibCjbjCm0]cm'bmViujvjum1Vm01Cibii=i,j,m如01Cibii=i,j,m-bi{/=BiBjBm^}e其中fei0j2A矩陣[B]稱為應變矩陣,或稱為幾何矩陣.由以上計算公式知,它與單元的節(jié)點坐標有關,但不隨點的坐標變化,就是說在這一單元內(nèi)所有點的應變是相同的.2〕求得單元應力的單元節(jié)點位移的表達式將幾何矩陣代入單元的物理方程,就有：J)e二D；)e=D1B:e彈性彈性矩陣[D]是由材料常數(shù)組成的矩陣.令[S]=[D][B],代入平面應力的物理方程,就有七產(chǎn)"St1e=—E71-」21sL2152A■七產(chǎn)"St1e=—E71-」21sL2152A■bi1七

1-NrC寫成分塊矩陣的形式rbi%1一」Ci001-」2」CCi

1-<

-Vbi」CiC

1-\

Tbi12Abj」bj1—J-bi10:G—Cj2jbj0CjCj1—Jbji=i,j,m0Cjbjbm

0CmbmJbm1-JC2Ui1%Cm1-\bm2一也可以3〕用虛功方程導出單元剛度矩陣〔單剛矩陣〕虛功方程%*,.｝=1｛屋yd及｝dvV假定單元的厚度為t,上式改寫為單元的虛功方程形式,

也*?T{f}e=J[Q*TIdLltdxdy虛應變也可以用幾何方程表示qf=ib七eq*tT=〔fefc*ri=〔fe*y>fei代入上式隹*FT{FF=H(%*}eTlB『bIfek}etdAA由于虛位移元素是常量,所以可以提到積分號以外,并與左邊的消去〔為什么〕.于是上式變?yōu)椋骸寒a(chǎn)=BTbIBLVtdxdy,A\令ke=以BTb】Btdxdy,虛功方程就成為了單剛方程由于［B］、［D］都是不含有x,y的常數(shù)矩陣,所以雙重積分實際就是對面積積分了.ke=BTbIBtAa——是三角形面積將前將前面求得的應變矩陣和彈性矩陣代入,然后作矩陣乘法.就得到我們要求得的矩陣計算公式.我們在這里采用分塊矩陣的方法記憶.2kj^^41-J22kj^^41-J2AHbsCr1-」—2"CrCs1-brCs2,,1-<JbrCs——bsCr21-」--2-brbsr,s=i,j,m注意以上是平面應力狀態(tài)的單剛矩陣,如果是平面應變問題呢二,單元剛度矩陣的性質(zhì)單元剛度矩陣［K］e表示單元反抗變形的水平.它與通常的彈簧剛度系數(shù)k的物理意義本質(zhì)相同,只不過［K］e是一個6X6階的矩陣.共有36個元素.這是由于三角形的節(jié)點力向量和節(jié)點位移向量都為6的緣故.1〕物理意義分塊矩陣［kj］表示的物理含義是：節(jié)點j處產(chǎn)生單位位移,而節(jié)點i、m被約束,此時在節(jié)點i處產(chǎn)生的節(jié)點力.我們寫出它的4個分量元素：12ij22在上面的矩陣中,元素的第一個下標表示產(chǎn)生節(jié)點力的節(jié)點,第二個下標表示產(chǎn)生節(jié)點位移的節(jié)點.上標“1〞表示水平分量,“2〞表示豎直分量.而且上標和下標的關系是對應的.也就是說第一個下標對應第一個上標；第二個下標對應第

二個上標.如此就有:12ki；2就表示節(jié)點j處產(chǎn)生豎直單位位移,在節(jié)點i處產(chǎn)生的水平方向的節(jié)點力.2〕單元剛度矩陣是對稱矩陣一〔keT二k?這一點可以通過簡單的數(shù)學證實如下：〔kH=〔BTb】B[aT=BTb1伯TT二bTDIB[a二k7單元矩陣的對稱性,從物理學角度反映出的道理就是,“功的互等〞.也就是在節(jié)點j處產(chǎn)生某一位移引起節(jié)點i處的節(jié)點力,應等于在節(jié)點i處產(chǎn)生相同位移引起節(jié)點j處的節(jié)點力.3〕單元剛度矩陣中的元素只與單元的材料性質(zhì)、幾何形狀、尺寸大小有關,而與單元的位置無關.單元剛度矩陣中不含有araj、am,上節(jié)對位移模式收斂性分析中,曾經(jīng)說明了ana4分別表示單元的平移分量,而1a；=--aiUi-ajUjamum2A1a4二一aiViajVjamVm2A由上式知單元的平移運動與ah小am有關,而[K]e又與ar電、am無關,所以說它不隨坐標軸的平動而變.4〕單元剛度矩陣是奇異矩陣,即KT=0從力學的角度理解單兀的剛度方程-K升療,當給定位移時,可以求得力；當給定力時,卻不能求得位移.由于[K]e不存在逆矩陣,在單元沒有給出任何約束的情況下,除有應變的可能性外,還同時有剛體位移的可能性.所以方程無解.荷的節(jié)點移置前面對于有限元模型的分析時,曾經(jīng)說過,單元之間的力傳遞是通過節(jié)點進行的.所以不在節(jié)點上的力,必須按靜力等效原那么,把它們移置到節(jié)點上.靜力等效原那么：原載荷在任何虛位移上所做的虛功,與移置到節(jié)點上的節(jié)點載荷所做的虛功相等.這種處理方法,和我們前面講到的圣維南原理相同.它們只會影響模型局部的應力分布,而不會影響整個結(jié)構的應力.下面根據(jù)力的類型,分類說明處理的方法.1〕集中力的移置如下圖的受力示意圖,M處的力為MbhPyT,移置后的節(jié)點載荷為{r*=XiYiXjYjXmym】T,虛功相等就是：{?=TM,=「：*P〕???忖〕=N'ide,所以有〔和*]￡r}=〔NKsMe"p}=〔fe*}eTInI{p}?.：R」NHP〕由于我們現(xiàn)在一直局限于單元的載荷移置,所以上式中的{R}應記為{R}eo如果將它寫成分量的形式：{R}e=XYXjYjXmYm『=NiRNiPyNjBNjPyNmPxNmPy1從上面公式可以清楚的看到,載荷移置結(jié)果與單元位移模式密切相關.2〕體力移置體積力密度為{p}=反Y1,將單元中微元體的體力{pXdxdy看作集中力,那么=N■ptdxdy一,A3〕面力的移置設在單元的一個邊界上作用有分布力,面力的密度為{p}=!XY『,將微面積tds上的面力也看作是集中力,那么有R=sN1ptds以上面力,體力的積分運算較為復雜,在線性模式下,可根據(jù)理論力學中靜力學平行力分解的原理,直接求得等效的節(jié)點載荷：③均質(zhì)等厚的三角形單元,受W的體力,那么三個節(jié)點分別受到W的節(jié)點載荷.3⑶i、j邊受到集度為q的均布面力載荷,那么i、j邊節(jié)點受到鄴的等效節(jié)點載荷.21③i、j邊受到線性分布的面力,i處集度為零,j處集度為q其合力為R=-qlt,那么i2一一一1一一2處的節(jié)點載荷為1R,j處的載荷為2R.33這些簡單的規(guī)那么,對于今后實際中的應用,可以提升效率.§2.14整體分析該過程是將離散別離的各個單元組集成離散的結(jié)構物,從而建立模型的總剛方程.一,總剛方程和總剛矩陣的組集1,總剛度矩陣的組集原那么A,整個離散結(jié)構變形后,各個單元在節(jié)點處仍然協(xié)調(diào)地相互連接.即環(huán)繞某個節(jié)點的n個單元,在節(jié)點i處具有相同的位移.數(shù)學公式的描述：W[=如￥="??=也F=也}B,各個節(jié)點應滿足靜力平衡條件.即每個節(jié)點上的節(jié)點力合力應等于該節(jié)點的節(jié)點載荷.

數(shù)學公式描述:此處的￡代表環(huán)繞節(jié)點i的所有單元節(jié)點力求和.e2,在該原那么指導下,實例的組集過程如圖a所示的平面問題,采用圖b的方法劃分網(wǎng)格,單元分析完成后,現(xiàn)將它們組集成原問題的有限元離散網(wǎng)格,求該問題的總剛矩陣.3為i,j,m=1,23為i,j,m=1,2,111單元G〕節(jié)點號〔i節(jié)點從最小號開始,然后逆時針排列節(jié)點號〕,單剛方程：kkk將它們展開311汨：1二憶/⑥)1?K2k2』?k3k3,下2、卜211Y?卜22匚2』IN"」22112222335j31aJ?K32k2?-K331C.3?完全相同的道理,得其他三個單元展開的單剛方程單元Q節(jié)點號為i,j,m=1,3,4KJKdYEJKlJJkHJkF、〕2單元單元o節(jié)點號為i,j,m=3,5,4g3=K33k3：K35k513K34濘/『5：'3=K53HG3-K55Hf3?K54k4?下d二K43HdK45HJK44N/單元Q節(jié)點號為i,j,m=2,3,5X二K22"、"K23S"K25k.5)5M=K32k2：4K33I〞/K35跳J(F514=K52"、LK"、/K"/實際上,剛度方程可以寫成如下的一般形式：{Fj=K工卜［瓦}k表示單元的節(jié)點組成k4,j,m首先我們運用規(guī)那么.2,即各個節(jié)點力的合力等于節(jié)點載荷“,括:e=Me所以也就有：；FJ,岸二口){f2}2+k2}4={r2}匕明岸『3~g3二代)下4)2?尸4節(jié)=R:'仁手?Is)4=g:'我們將上面得到的方程按這個規(guī)律相加：匕汜：2=4：'=匕1工：屋k2)1413憶3)1—2."憶3匕3〞憶4憶4)2對這個等式運用規(guī)那么.i,即%了=也F=…=電}n=%}并整理有：叱Kii1KI?〔KM%".2"：,代)=憶1〞「靖工仁：'"小14.)％"一代：'=^311-K312''ll卜321,K32I'J)k331?K3312k333k334",K342K3433-^35-K354X代)二k1憶JK4312K43M3"K4412K4414；憶585)

應｝=k52M&｝+收537+卜531近｝+-543J+fKJ+1^55/砧5｝將上述的這些方程寫為矩陣的形式,就得到了該問題的總剛方程為:RiR2」RRiR2」R3,=R4R5「ki產(chǎn)Ik2ilIkJ-0ki21

卜22F

卜32「0

k524kJ2kJk23140k33l2-23匕3卜3-2324k54『卜35/IS'1〞與花3?占40J我們對這個總剛方程進行分析,可以得出以下的規(guī)律3,總剛矩陣的規(guī)律得到結(jié)構的節(jié)點載荷,節(jié)點載總剛方程就是運用所有節(jié)點的位移向量和總剛矩陣相乘,荷向量〔外界對彈性體作用的載荷〕很容易求得,節(jié)點位移向量得到結(jié)構的節(jié)點載荷,節(jié)點載陣如何得出.如都根據(jù)上述的步驟,成千上萬個節(jié)點,該如何作呢所以有必要總結(jié)總剛度矩陣的規(guī)律..當r=s時,即子塊矩陣［Krs］是總剛矩陣主對角線上的子塊矩陣.由環(huán)繞節(jié)點r〔s〕各單元剛度矩陣的相應對角線子塊相加.②當rws時,但r、s是相鄰單元的公共邊上的節(jié)點時,子塊矩陣［Krs］等于兩相鄰單元剛矩陣的相應子塊矩陣相加.如上例中的［K13］、［K31］等.③當rws時,且r、s只是一個單元的邊上節(jié)點時,子塊矩陣［Krs］就等于該單元剛矩陣的相應子塊矩陣.如上例中的［K12］、［K21］等.②當rws時,且r、s不屬于一個單元的邊上的節(jié)點時,子塊矩陣［Krs］等于零.如上例中的［K15］、［K51］等.利用上述規(guī)律不僅可以檢驗總剛組集的正確與否,而且可以直接組集總剛矩陣〔手工組集〕.但這種方法不利于計算機編程.下面介紹計算機組集總剛矩陣的方法4,總剛矩陣組集的步驟③擴大各單元剛度矩陣,使之成為與總剛矩陣相同的階數(shù).②根據(jù)總體節(jié)點的編號,將各單元剛度矩陣的各個子塊移到相應的位置上.其余位置充零.⑶把各個改造過的矩陣直接相加,就得到總剛矩陣.我們以單元.4為例講解步驟.1、Q單元4的i、j、m=2、5、3,所以單元剛度方程是K25K55K3K25K55K35h's"偉2川K4k52K4K33電｝4值｝43-4.1,.4?222532kkkTOC\o"1-5"\h\z000000k224k2340K2J0K32TK3310K35700000-0K524K5340叱557_二,總剛度矩陣的性質(zhì)總剛度矩陣由單元剛度矩陣組集而成,所以也具有單元剛度矩陣的一些性質(zhì),如相同的物理意義,位置無關性、對稱性和奇異性等,還具有以下性質(zhì)：1〕稀疏性由規(guī)律4知,當rs,且r、s不屬于同一單元白兩個節(jié)點是[K「s]=0,說明互不相關的節(jié)點數(shù)愈多,零子塊矩陣就愈多.一般說來相關節(jié)點數(shù)不超過9,而整個分析對象常常成百上千、上萬或幾十萬.如果整體有100個節(jié)點,那么百分之九是非零子塊,而百分之九十多都是零子塊.所以在總剛度矩陣中非零的子塊是很稀少的.2〕帶狀性所謂帶狀性,就是指總剛矩陣中的非零子塊,集中分布在主對角線的兩側(cè),呈帶狀分布.半帶寬值就是計算這帶狀寬度的數(shù)值.它是以排列元素最長的一行,從第一個非零元素起至主對角線元素止,所有元素的個數(shù).其數(shù)值可以由節(jié)點的總體編號算出.B=〔相鄰節(jié)點號的最大差值+1〕X〔單個節(jié)點的自由度個數(shù)〕對于平面問題說來,就是：B=〔相鄰節(jié)點號的最大差值+1〕X2下列圖中的兩種編號方式,可以的分別計算其半帶寬值:圖a,B=〔2+1〕x2=6圖b,B=〔6+1〕X2=141士―彳國身彳國身兇1廠71工半帶寬值影響計算機存取總剛矩陣所需的內(nèi)存大小.愈小愈好.由于半帶寬值是直接受節(jié)點總體編號的影響.所以我們在建立有限元模型時,應慎重考慮,采用優(yōu)化的方法.〔現(xiàn)在通常在程序中都配有這樣的優(yōu)化程序〕.三,總剛矩陣的壓縮存取技術由于總剛矩陣具有對稱性,所以我們只需存入主對角線上半帶或下半帶的元素,就可以完成解方程的運算.此即為總剛矩陣的半帶寬存儲方法.假定取主對角線上半帶元素存儲,具體做法就是每行元素以主對角線上的元素開始,存儲每行半帶寬數(shù)值的元素個數(shù).如此各元素的行號不變,改變的只是列號.新列號和原列號的關系式如下：新列號=原列號-行號+1前面圖示的例子,如果采用圖示表示總剛矩陣的存儲,就是：如下列圖所示的情況.可以看出原總剛矩陣存儲需要24X24的2維數(shù)組,改為半帶寬存儲后,就成為了24X6的數(shù)組了.

另有一維的壓縮存儲方法,比這種二維的半帶寬存儲方法壓縮更多.在這兒我們就不做介紹了.有興趣的同學可以看相關的參考書籍.講到這兒,實際我們已經(jīng)知道,每個元素在總剛矩陣中的位置,在節(jié)點編號完成后,就完全確定下來了.所以計算機對總剛矩陣的計算,是一部完成的.即“邊對號移置、邊改列號、邊累加〞.四,總體邊界條件的處理前面介紹總剛矩陣的性質(zhì)時,說明了總剛矩陣是奇異矩陣.即|k】二0.就是說總剛矩陣不存在逆矩陣.要求得節(jié)點位移的位移解,還必須引入邊界條件.1,邊界條件的類型①節(jié)點固定,即Ui=Vj=00給定節(jié)點位移值.即Ui=Ui,Vj=Vj2,處理方法1置“0、1〞法即首先在總剛矩陣［K］中與位移分量相對應的行和列元素改為0,但主對角線上的元素改為1,然后在節(jié)點載荷向量的列陣中,與對應元素的位移用代替,其余元素減去位移分別乘［K］中相應元素.例：某結(jié)構的總剛方程為卜］0乳0缶｝10刈=010M,3=5、v2=c4,根據(jù)上面的方法修改如下：如果展開上式,馬如果展開上式,馬上就有U1=c1、V2=C4第二行k22V1k23U2k25U3k210V5二丫1一k21c1一k24c4一1000…0]lrU1L0k22k230k210V1丫1—k21c1—k24c40k32k330k310U2>=<X2-k31c1-k34c4>0001…0V2c4iaaasi-a——■■■■■■■■■■■■■■J0k102k1030k10101V5,丫5—k101c1—k104c4.移項k21.X22V1k23U2k24c4k25U3k21V5=丫1從中可以看出,這樣處理后不僅可以直接得到U1=c1、V2=c4,而且它們產(chǎn)生的效果也計入到了其他方程中.如果是固定位移情況,那么載荷向量中的變化是什么呢〔好好想想〕2乘大數(shù)法在總剛矩陣［K］中,把與位移位移相對應的行與列主對角線上的元素乘一個很大的數(shù)101°,然后把載荷向量中的對應元素代以給定位移乘以相應主對角線上的元素,再同樣乘以一個很大的數(shù).同上例./-1/-10,,,,「父10kI2k〔3k〔4…ku.C[k[1M10k21k22k23k24…卜210V1Yk31k32k33k34卜310U2X2,,,,“10,<>=<k41k42k43k44M10…k410V2C4k44父10amm*mm■.ak101k102k1030k1010V5丫51010、、考察第一個方程:k1k12V1k110V5=c1k111010⑤求解⑤求解應力上乒=D.Lje-DI〔b,’e兩邊同除以"黑1010,由于k11x1010>>k1j,所以U1=c1同理可得V2=C4,〔思考該方法不能用于固定位移的情況,為什么〕必須說明的是,當總剛矩陣采用半帶寬存儲的時候,以上邊界條件的引入也應當根據(jù)半帶寬存儲的格式修改,否那么就是錯誤的.實際應用中,當是第一種情況時,采用方法.1；當?shù)诙N情況時,用方法O2.五,應力計算總剛方程在引入足夠的邊界條件后,就可以進行求解了.所謂足夠的邊界條件就是要使系統(tǒng)是一個幾何不變的系統(tǒng).數(shù)值解方程的方法大多采用高斯消元法.當求解出位移向量以后,就可以通過單元的幾何方程求解應變,通過物理方程求解應力了.這個過程稱之為回代過程.@求解應變行=B㈠由于三角形單元是常應力與常應變單元,所以由上述方法求得的應力或應變看作是形心處的應力或應變.而且很明顯,這樣求得的應力或應變是不連續(xù)的.為了推算彈性體內(nèi)某一點的接近實際應力,我們通常采用以下兩種方法,來平滑應力的突變.①繞節(jié)點平均法一一就是將繞某一節(jié)點的各單元形心應力加以平均,來表示該節(jié)點的應力.如左圖所示.二x1J二Ax.二Bx12二--1i0AxBx,■.C,■.D,■.E,1.Fx22xxxx該方法在各單元面積相差不大,單元的內(nèi)部節(jié)點時較為準確.而在外部邊界節(jié)點上,那么不理想.所以對于邊界節(jié)點多采用三點插值的方法求得.也就是用內(nèi)部的節(jié)點來推算.,兩單元平均法一一把相鄰兩單元的應力相加平均,表示兩單元公共邊中點的應力值.

Ox2二；二晨,二Bx邊界邊上的應力值也采用插值方法得到.卜面通過一個例子來看看,具體至此我們對于有限元求解平面卜面通過一個例子來看看,具體如何計算的.題計算,有限元法求解平面問題的步驟:1〕根據(jù)分析對象和給定的條件,根據(jù)一定的比例繪出計算簡圖,該圖應有各部位的尺寸、外力和支承情況.2〕選擇適宜的坐標系,劃分單元準備數(shù)據(jù).數(shù)據(jù)包括節(jié)點數(shù)、編號和坐標,單元數(shù)、編號和節(jié)點組成；材料特性常數(shù)；載荷大小及位置；邊界條件.3〕計算單元的面積、應變矩陣、彈性矩陣和單元的剛度矩陣.4〕組集總剛矩陣,處理載荷移置,形成載荷向量.5〕處理邊界條件,修改總剛矩陣,得到最后的剛度方程.6〕求解線性方程組,得到位移向量.7〕求解應力值.三,例題設有如下圖的正方形薄板,在對角線上作用有沿厚度均勻分布的載荷,其合力為2kN,板厚為1個單元厚度,為使計算簡便,假定材料的科=0,求該板的變形.分析：該薄板對稱于對角線,受力也是如此,所以可取其1/4來計算.由于是取1/4,所以該模型的兩個直角邊,都是對稱線.X的對稱線不能有x方向的位移,y軸的對稱線不能有y方向的位移.直角的頂點既是x軸對稱線上的點,又是y軸對稱線上的點,所以該點在x、y方向的位移都應該約束.最后得出的計算模型如下.解：1準備的數(shù)據(jù)A單元節(jié)點數(shù)據(jù)坐標123456X001012Y211000

B單元信息數(shù)據(jù)①鳴⑷⑷I1223J2455M3536所有單元的板厚是“1〞C材料常數(shù)彈性模量E,泊送比科=0D載荷數(shù)據(jù)00000000TE位移邊界1、-10Vl0v2u3V3E位移邊界1、-10Vl0v2u3V300u50u601T2計算單元剛度矩陣1xiyi一.i.單兀1面積A=-1x2y221X3y3b1=y2-y3=1-1=0C1=-X2X3=01=1形成的幾何矩陣b1B1=1020C1b1C1彈性矩陣10單元1的剛度矩陣111210201=0.513b2=y3-y1=1-2=一1b3=y1-y2=2-1=1c2--x3-x1--1-0--1C3=一x1x2=00=0b20b3000-1cc1Lcc0c20c3=0000..2屋一八,c2b2c3b310-1010-100-1010.0.50-0.5010k1=B「b】B1A=E-0.501.52|-0.5-10.500-10.500.5-1000.5-1-0.51.50-0.5010「0「0.50-0.5-0.500.5同理可求得k2、k產(chǎn)和k4.需要說明的是可以利用位置無關性直接得出kI2和k了,1000.5k3=E.0.5200由于它們是i1000.5k3=E.0.5200TOC\o"1-5"\