第八章一般殼體問題的有限元法

第八章一般殼體問題的有限元法第1頁，共48頁，2023年，2月20日，星期三

將殼體曲面劃分為有限個單元，它們都是曲面單元。但是在單元細分時，用平面單元組成的一個單向或雙向折板來近似殼體的幾何形狀將會得到良好的結(jié)果。通常對于任意形狀的殼體，采用三角形單元比較方便，如圖8-1所示。如果在殼體上容易找到同一平面上的四個點，可以采用平面四邊形單元。例如具有正交邊界的柱面殼體，如圖8-2所示。圖8-1任意殼體作為平面三角形單元的集合圖8-2圓柱殼作為平面矩形單元的集合§8-1平面殼體單元第2頁，共48頁，2023年，2月20日，星期三殼體平面單元的應(yīng)力狀態(tài)是由平面應(yīng)力和彎曲應(yīng)力的疊加而成的，因此在構(gòu)造殼體平面單元時，只要將第二章和第七章所討論的相應(yīng)單元進行簡單的組合就可以了。同樣，前述二章所導(dǎo)出的剛度矩陣可作為建立殼體平面單元剛度矩陣的基礎(chǔ)。現(xiàn)在把平面單元的計算步驟歸納如下

1.劃分單元，選定整體坐標(biāo)系oxyz，定出節(jié)點在整體坐標(biāo)系中的坐標(biāo)值。

2.對于各個單元利用節(jié)點坐標(biāo)值，建立一個局部坐標(biāo)系例如三角形單元123，可以選取節(jié)點1為局部坐標(biāo)系的原點，并且以1-2邊為軸的正方向，如圖8-3所示。于是，方向的單位e1求得是(8-1)第3頁，共48頁，2023年，2月20日，星期三式中是矢量12的長度。取單元的外法線方向作為軸的正方向，于是它的單位矢量圖8-3三角形單元局部坐標(biāo)系第4頁，共48頁，2023年，2月20日，星期三容易看出，矢量12和13的矢性積的模等于三角形面積Δ的一倍，即|12×13|=2Δ。最后，按右手定則可以決定y軸的正方向,它的單位矢量e2是e2=e3×e1(8-3)利用上述方法確定的局部坐標(biāo)系，三角形單元123是在平面內(nèi)，它的三個角點的局部坐標(biāo)值是很容易確定的。對于柱面上的矩形單元，局部坐標(biāo)的原點選在矩形的形心，通常選軸和x軸均沿柱面母線方向。如圖8-4中所示，由矢量12確定單位矢量e1,再由矢量14確定單位矢量e2，于是e3=e1×e2。第5頁，共48頁，2023年，2月20日，星期三

3.對于各個單元，確立在局部坐標(biāo)系中的結(jié)點載荷列陣。殼體載荷可以分解成二組：一組作用在平面內(nèi)，另一組垂直于平面。為此，在計算各個單元的結(jié)點載荷列陣（包括等效結(jié)點力）可以直接引用第二章和第七章中所敘述的載荷計算的相應(yīng)公式。各個單元的結(jié)點載荷列陣圖8-4矩形單元局部坐標(biāo)求得后，建立變換矩陣公式，從而把轉(zhuǎn)換到整體坐標(biāo)系中去求出在整體坐標(biāo)下的單元節(jié)點載荷列陣，然后經(jīng)各單元的簡單疊加可以求出結(jié)構(gòu)在整體坐標(biāo)下的節(jié)點載荷列陣。第6頁，共48頁，2023年，2月20日，星期三顯然，平面單元在局部坐標(biāo)系中，結(jié)點i有五個廣義位移：即

，其中前兩個對應(yīng)于平面應(yīng)力問題，后三個對應(yīng)于平板彎曲問題。類似地，所對應(yīng)的結(jié)點力列陣

顯然，在上式中實際上總是等于零的。。為了經(jīng)坐標(biāo)變換后不影響在整體坐標(biāo)系中對各特征量的計算，我們引進(8-4)第7頁，共48頁，2023年，2月20日，星期三容易看出，把以上結(jié)點位移和結(jié)點力變換到整體坐標(biāo)中后，他的結(jié)點位移和結(jié)點力列陣具有如下形式上式右端的前三項分別表示位移和力，后三項分別表示轉(zhuǎn)角和力矩，它們都是有明顯物理意義的矢量。因此，（8-4）式和（8-5）式之間的坐標(biāo)變換公式是

式中而(b)

(8-5)(8-6)(a)第8頁，共48頁，2023年，2月20日，星期三于是，殼體單元e在局部坐標(biāo)下的結(jié)點位移列陣是或而所對應(yīng)的單元節(jié)點力（包括等效節(jié)點力）列陣是或式中n=3是對應(yīng)于三角形單元；n=4對應(yīng)于四邊形單元。本節(jié)以下的n所指的意義均是如此，不再重復(fù)說明。(c)

(d)

(e)(f)第9頁，共48頁，2023年，2月20日，星期三

4.建立局部坐標(biāo)系中的單元剛度矩陣，從而求出整體坐標(biāo)系中的單元剛度矩陣。如果將單元剛度矩陣和對應(yīng)于單元節(jié)點劃分為n×n個子矩陣，每個子矩陣都是6×6的，于是的子矩陣有如下形式式中和分別是平面應(yīng)力問題和平面彎曲問題的相應(yīng)子矩陣，它們是2×2和3×3矩陣。圖8-5示出了在局部坐標(biāo)系中三角形殼體單元剛度矩陣用平面應(yīng)力和平面彎曲剛度矩陣的構(gòu)成方法。(8-7)第10頁，共48頁，2023年，2月20日，星期三圖8-5三角形殼體單元剛度矩陣用平面應(yīng)力和平板彎曲剛度矩陣的構(gòu)成方法第11頁，共48頁，2023年，2月20日，星期三單元e中任意結(jié)點i的平衡方程，在兩個坐標(biāo)系中分別為，式中是剛度矩陣的子矩陣。而對于局部坐標(biāo)和整體坐標(biāo)之間的變換公式是把(h)式代入(g)式得(g)(h)第12頁，共48頁，2023年，2月20日，星期三將公式(g)中的第一式左乘矩陣，并且同上式進行比較，可以得到由于是正交陣，容易證明也是正交陣，即。這樣就得到關(guān)于矩陣的轉(zhuǎn)換公式(8-8)第13頁，共48頁，2023年，2月20日，星期三

5.集和單元剛度矩陣及等效結(jié)點力。線作簡單求和然后將它們放入整體剛度矩陣[K]和等效結(jié)點荷載列陣的相應(yīng)位置上去。第14頁，共48頁，2023年，2月20日，星期三

6.修改整體剛度矩陣，然后求解平衡方程式中是總的結(jié)點位移列陣。特別值得注意，在局部坐標(biāo)系中單元剛度矩陣[K]對于三角形單元它的第6、12、及18行和列全是零元素，對于四邊形單元它的第6、12、18及24行和列全是零元素，其原因是轉(zhuǎn)角，并不包含在平面應(yīng)力的單元結(jié)點位移列陣中。當(dāng)所有在一個結(jié)點相連接的單元共面時，殼體結(jié)構(gòu)的剛度矩陣將是奇異的。避免這個奇異性的一個辦法是引入關(guān)于殼體法線的轉(zhuǎn)動為零的附加條件。第15頁，共48頁，2023年，2月20日，星期三

7．計算應(yīng)力。首先是按照公式求出局部坐標(biāo)系中的結(jié)點位移，再按第二章中所給出的公式計算應(yīng)力、和；通過第七章所給出的公式計算、和進而求得應(yīng)力和。、和。于是，殼體應(yīng)力可以由簡單的疊加求得；即，第16頁，共48頁，2023年，2月20日，星期三

對于一個殼體結(jié)構(gòu)，如果采用上節(jié)所述的平面單元，將會引起幾何上的離散誤差。人們希望采用曲面單元來描述殼體的真正幾何形狀，使之用不太多的單元來替代復(fù)雜形狀的殼體，并得到具有一定精度的解答。另外，在薄殼理論中都是用中面位移來表示中面轉(zhuǎn)動。正如在第七章中所述，這將要求在單元交界面上有橫向位移及其一階導(dǎo)數(shù)的連續(xù)性，于是增加了選擇位移模式的困難。如果考慮橫向剪切變形的影響就可以認(rèn)為中面轉(zhuǎn)動是獨立變量而不依賴于位移的一階導(dǎo)數(shù)。因此，只要利用單元交界面上位移函數(shù)的連續(xù)性就可以了，并不要求其一階導(dǎo)數(shù)的連續(xù)性。

現(xiàn)在我們來論述一個考慮橫向剪切影響的曲面單元，稱為八結(jié)點四十個自由度的一般殼體單元，如圖8-6所示。8-2考慮橫向剪切變形影響的殼體單元第17頁，共48頁，2023年，2月20日，星期三

在圖8-6中所示的殼單元，象空間等參數(shù)單元一樣引進一個自然坐標(biāo)系。命為殼體中面上的曲線坐標(biāo)；對應(yīng)于的表面稱為頂面（或上表面），對應(yīng)于的表面稱為底面（或下表面）。在單元的中面上選取八個點稱為結(jié)點，過各結(jié)點i(i=1,2,…,8)作中面的法線，交頂面和底面的點稱為結(jié)點i的對點。結(jié)點i相對應(yīng)的對點，它的整體坐標(biāo)值分別記作圖8-6八結(jié)點四十個自由度的一般殼體單元一.單元幾何形狀的確定第18頁，共48頁，2023年，2月20日，星期三于是，中面上的結(jié)點i的整體坐標(biāo)值是顯然，結(jié)點i處的中面法線方向可以由下列單位矢量所確定(8-10)

(8-9)第19頁，共48頁，2023年，2月20日，星期三式中l(wèi)3i、m3i和n3i是結(jié)點i處中面法線方向?qū)τ谡w坐標(biāo)軸oxyz的方向余弦，而hi是結(jié)點i處的殼體厚度，即

(a)結(jié)點i處法線上任意點的整體坐標(biāo)值，可以通過矢量相加得到（圖8-7），即(b)第20頁，共48頁，2023年，2月20日，星期三于是，單元內(nèi)任意點的坐標(biāo)值可以通過形函數(shù)的插值得到，即式中形函數(shù)由公式（5-1）表示。這樣，我們就可以通過八對點的整體坐標(biāo)值，按照（8-11）式近似地確定了單元的形狀。當(dāng)時，分別確定上下表面各點；當(dāng)時，確定中面各點；而單元的側(cè)表面是由中面法線（或近似的中面法線）所構(gòu)成。(8-11)第21頁，共48頁，2023年，2月20日，星期三

假設(shè)中面法線變形后仍為直線，但是不再是變形后的中面法線。也就是說，中面法線有繞兩不同軸的轉(zhuǎn)動，這兩不同軸分別垂直于法線。設(shè)V1i和V2i分別是這兩軸上的單位矢量，顯然它們具有一定的任意性。為了唯一地確定起見，不妨設(shè)式中(c)

(8-12)二.位移模式第22頁，共48頁，2023年，2月20日，星期三于是有而另一軸的單位矢量可以由V3i和V1i的矢量積求得注意，在（8-12）式中如果i×V3i=0，則可用j×V3i來代替，而(j)(d)

(e)

(8-13)第23頁，共48頁，2023年，2月20日，星期三設(shè)結(jié)點i處的中面法線V3i繞V1i和V2i兩軸的轉(zhuǎn)角分別為和，如圖8-8所示。于是，轉(zhuǎn)動矢量可以寫成圖8-8結(jié)點的自由度(8-14)第24頁，共48頁，2023年，2月20日，星期三現(xiàn)在我們來計算結(jié)點i處法線上任意點的位移值。這可以通過結(jié)點i在整體坐標(biāo)中的位移ui、vi、wi以及法線繞i點轉(zhuǎn)動而產(chǎn)生的位移相加而得到；關(guān)于轉(zhuǎn)動引起的位移可以根據(jù)運動學(xué)中的公式計算。于是，在結(jié)點i處法線上各點在整體坐標(biāo)系中的位移是顯然。于是，上式就可以寫成(h)(g)第25頁，共48頁，2023年，2月20日，星期三這樣，利用插值方法可以得到單元內(nèi)任意結(jié)點的位移列陣是式中(8-16)(8-15)第26頁，共48頁，2023年，2月20日，星期三可以將（8-15）式寫成標(biāo)準(zhǔn)形式式中[I]是三階單位陣。(8-17)(i)第27頁，共48頁，2023年，2月20日，星期三

在整體坐標(biāo)系中，利用幾何方程和（8-17）式可以將應(yīng)變列陣寫成標(biāo)準(zhǔn)形式是（8-18）三.應(yīng)變計算第28頁，共48頁，2023年，2月20日，星期三式中而(j)(8-19)第29頁，共48頁，2023年，2月20日，星期三因此，為了計算應(yīng)變需要算出六個偏導(dǎo)數(shù)、、和、前三個偏導(dǎo)數(shù)可以按照公式（5-54）進行計算，后三個則按定義分別是[J]-1的第三列的三個元素，即由（5-58）式得到(8-20)第30頁，共48頁，2023年，2月20日，星期三出現(xiàn)在所提到的各式中的S、T和V分別按公式（5-55）所定義。將（8-11）式代入（5-55）式得(8-21)第31頁，共48頁，2023年，2月20日，星期三

我們還假設(shè)中面法線上的線段不伸長也不縮短。如果在單元中面各點處建立局部的直角坐標(biāo)系，而軸是中面法線方向，于是有；這里表示方向的位移分量。建立局部坐標(biāo)系的關(guān)鍵是中面上各點的法線方向余弦矢量V3，一個合理的簡單做法是由單元各結(jié)點處中面法線的方向余弦矢量V3i通過內(nèi)插法構(gòu)造而成，即(8-22)四.坐標(biāo)變換第32頁，共48頁，2023年，2月20日，星期三中面法線的方向余弦矢量構(gòu)成后，可以按照如下方法作出軸及軸的方向余弦矢量V1和V2于是，在坐標(biāo)系中的矢量都可以通過下列變換矩陣變換為整體坐標(biāo)系oxyz中的矢量。(8-23)(k)第33頁，共48頁，2023年，2月20日，星期三在局部坐標(biāo)系中，它的應(yīng)變列陣是根據(jù)應(yīng)變張量的定義，立刻可以寫出應(yīng)變的坐標(biāo)形式式中(8-26)(8-24)(8-25)第34頁，共48頁，2023年，2月20日，星期三

在殼體中還假設(shè)，于是在局部坐標(biāo)中有如下的應(yīng)力應(yīng)變關(guān)系式中以及(8-28)(8-27)(1)五.彈性矩陣的變換第35頁，共48頁，2023年，2月20日，星期三在整體坐標(biāo)系oxyz中，應(yīng)力應(yīng)變關(guān)系雖然仍可寫成如下形式式中。但是彈性矩陣[D]不再像以前那樣具有簡單的形式，而需要用坐標(biāo)將變換成[D]。由于這兩組應(yīng)力表示同一個應(yīng)力狀態(tài)，因此在任意虛位移上，兩個坐標(biāo)系中單位體積內(nèi)的應(yīng)力所做的虛功必須相等，即是(8-30)(8-29)第36頁，共48頁，2023年，2月20日，星期三將（8-27）式代入上式并利用（8-25）式，由于的任意性得和（8-29）式相比較，得到這就是坐標(biāo)變換情況下彈性矩陣的變換公式?？梢宰C明這一形式具有普遍意義。(m)(8-31)第37頁，共48頁，2023年，2月20日，星期三

習(xí)慣下，我們要計算的殼體應(yīng)力是對應(yīng)于局部坐標(biāo)系的。利用（8-25）式和（8-18）式，公式（8-27）可以寫成其中(8-33)(8-32)根據(jù)（8-19）式和（8-26）式，很容易地算出矩陣的顯式。為了使公式書寫不致太長，我們給出它的分塊形式，即(8-34)六.應(yīng)力計算第38頁，共48頁，2023年，2月20日，星期三式中(n)(o)(p)(q)第39頁，共48頁，2023年，2月20日，星期三而(j=1,2,3)(r)于是(8-35)第40頁，共48頁，2023年，2月20日，星期三式中(s)第41頁，共48頁，2023年，2月20日，星期三在進行強度校核時，只要計算殼體上下表面上的應(yīng)力，即(8-36)第42頁，共48頁，2023年，2月20日，星期三容易算出(8-37)(8-38)第43頁，共48頁，2023年，2月20日，星期三

若將單元剛度矩陣寫成如（5-61）式的形式，于是