運用塊體單元法的疊層梁應力分析及破壞模擬

陸曉敏，袁 濤

(河海大學工程力學系，南京 210098)

運用塊體單元法的疊層梁應力分析及破壞模擬

陸曉敏，袁 濤

(河海大學工程力學系，南京 210098)

應用平面塊體單元法建立由剛體位移、常應變、曲率組成的變形模式，推導了系統的平衡方程，開發了相應的數值模擬程序，計算了簡支梁在均布荷載作用下的位移和應力。結果表明:在較少單元情況下位移和應力非常接近于理論解;相對于有限單元法，塊體單元法的位移精度更高，一般是理論值的上限;此外，接觸面中的應力精度很高，確保了模擬接觸面破壞的可靠性。為模擬疊層結構變形的非線性，在塊體之間的接觸面內引入材料的非線性，分別考慮剪切屈服和開裂2種破壞形式。屈服判據采用M-C準則，開裂采用最大拉應力強度準則，并編制了計算機程序。最后，計算模擬了一疊層式、受均布荷載懸臂梁的應力、變形及層間剪切的破壞。結果表明:此方法能較好地模擬結構的不連續變形及破壞。

疊層梁;塊體單元法;數值模擬

目前，疊層結構的復合材料在生活和工程中逐漸得到廣泛應用，如家具材料中的復合板、木工板、工程結構中的復合梁等。這些材料具有層狀結構，其破壞的特點常呈現層間錯動和開裂、層內彎曲斷裂等，與均勻連續各向同性的材料有明顯的不同。所以，建立起能反映層狀結構變形及破壞特點的數值模擬方法具有一定的工程應用價值。

數值模擬是力學中研究結構變形和安全度的一種有效方法。其中有限單元法［1］是最常用的一種數值方法，但其在處理層狀結構時常需進行特殊處理，異常麻煩，尤其是其常用的位移模式不能很好地反映層體材料彎曲變形的特點。在塊體單元法中，由于其以塊體的變形為求解未知量，塊體與塊體之間有單獨處理的接觸面單元，因而在模擬非連續變形時具有一定的優勢。本文據此采用塊體單元法來研究疊層結構復合材料的變形、應力及破壞形態。

塊體單元法的基本思想:將固體分割成塊體系統，假設塊體的變形模式，在塊體之間人為引入接觸面單元;接觸面單元的變形由相鄰塊體的相對變形而決定;系統的平衡方程由最小勢能原理導出。早期的塊體單元法采用剛體位移模式，塊體的變形常作簡化處理。剛體位移模式的塊體單元法增加了系統的剛度，求得的變形量往往偏小，因而也影響了計算的精度。有文獻提出了一些簡化處理的方法［2-3］，但過于簡單，尤其是沒有考慮材料變形的泊松比效應，給實際工程結構分析帶來不少困難。

本文借鑒DDA的思想［4-5］，把塊體的應變作為獨立未知量，并在塊體的變形模式中進一步引入曲率作為獨立變量，建立起具有彎曲變形特點的塊體單元法理論，并在接觸面內引入材料剪切屈服和開裂的非線性模型，使之能很好地模擬層狀結構的彎曲變形。此法的優點一方面能提高計算精度;另一方面，由于考慮了塊體的變形，從而可以模擬材料非線性［6］的本構關系，尤其在塊體單元的接觸面中可以較方便地引入各種材料性質，如彈塑性、黏性流變等。本文算例表明:塊體單元法理論具有較高的精度，能很方便地模擬層狀結構的層間屈服、開裂的破壞過程。

1 平面塊體單元法基本理論

1.1 塊體的變形模式

對于平面的塊體元，塊體的位移和變形模式可以用3個剛體位移、3個平均應變和2個曲率來表示，用矩陣表示為

式中:ui，vi為i號塊體在整體坐標系下x，y軸向的位移;θi為繞形心C的轉角，在小變形的前提下這些位移是微小量;εxi，εyi，γxyi為i號塊體在整體坐標下的平均線應變和切應變;κx'i，κy'i為i號塊體在單元局部坐標系Cx'y'下的平均彎曲曲率，為使用方便，曲率定義在局部坐標系中。約定局部坐標系中的x'軸與單元某邊平行，此方向常定義為層狀結構的接觸面方位。曲率與曲率半徑的關系為

式中ρx'i，ρy'i分別為i號塊體在x'和y'軸向的平均曲率半徑。

如圖1所示:i號塊體的形心坐標用(xCi，yCi)表示;塊體上某一點A的整體坐標為(xi，yi)，局部坐標為(x'i，y'i);變形前局部坐標軸x'與x軸的夾角用φi表示。在小變形假設的前提下，A點的位移可表示為

式中Ni稱為塊體的形函數矩陣，其各元素為:

1．2 塊體的應變及應力

根據式(3)可以將塊體內任一點A的位移用

式中Bi稱為i號塊體單元應變矩陣，各元素為:

假設塊體是各向同性的材料，彈性模量為E，泊松比為υ，則考慮平面應力時的彈性矩陣為

任一點A的應力為

圖1 塊體單元與局部坐標示意圖

1．3 接觸面單元的應力和應變

在塊體單元法中，人為在塊體接觸面之間引入很薄的接觸面單元，見圖2。設在塊體i和j之間有一薄層單元(簡稱為接觸面單元，編號為k)，厚度為h(很小)，C為形心，β為接觸面單元的局部坐標x'軸的方位角。接觸面上的初始接觸點B和B'的位移可以用式(3)分別進行計算。

將位移轉換至接觸面的局部坐標系中，則

式中Ti為局部坐標與整體坐標間的轉換矩陣。假設沿接觸面法向(y'軸)的切應變和線應變均勻變化，x'向的線應變為塊體i，j在此線應變的平均值，即

圖2 接觸面單元與局部坐標示意圖

據此可以推導出接觸面單元內的應變矩陣ε'k為

式中:(x'B'，y'B')為B'點在塊體i局部坐標系中的坐標;(x'B，y'B)為B點在塊體j局部坐標系中的坐標。式(14)中B'k稱為k號縫面單元的應變矩陣。

假設接觸面單元材料的彈性模量為E'，泊松比為υ'，引入彈性矩陣

則接觸面內的應力為

1．4 系統的平衡方程

系統的平衡方程可利用最小勢能原理導出。系統的內力勢能由塊體的變形勢能和接觸面單元的變形勢能兩部分組成。塊體i單元的變形勢能為:

式中kBi稱為塊體i的單元勁度矩陣，具體計算式為

接觸面單元k的變形勢能為

式中kCk為接觸面單元k的勁度矩陣，其計算式為

系統總的彈性變形勢能為:

式中:Nc，NB分別為系統塊體單元和接觸縫面單元的個數;d為系統的未知量(位移和應變)列陣;K稱為系統的勁度矩陣，由所有塊體單元和接觸面單元的勁度矩陣組裝而成。

系統的外力勢能為

式中:Fi，pi，qi分別為作用在i塊體上的集中力、面力和體力矢量(矩陣);為相應的等效荷載列陣

式中:s，v分別為面力作用的面積和體力作用的體積。式(21)中R為系統的等效荷載列陣，由各塊體單元的等效荷載列陣組裝而成

1．5 彈塑性塊體單元法

在塊體元模型中，既可以考慮塊體材料的塑性，又可以考慮接觸縫面的塑性。考慮到對于疊層梁的屈服和開裂破壞主要發生在結構面內，故本文暫只考慮縫面材料的非線性。如接觸面單元屈服，則在計算其勁度矩陣和縫面應力時，需用彈塑性矩陣Dep代替彈性矩陣D'。設:縫面材料的抗拉強度為σt;抗壓強度為σc;粘結強度為c;摩擦因數為f。縫面單元的破壞形式主要是剪切屈服和法向開裂，剪切屈服采用M-C準則判別，開裂由法向拉應力判別。

如果接觸面單元法向應力σy'＞σt，接觸面單元開裂，則認為其剛度喪失。如果

則縫面發生剪切屈服。根據相關聯的流動法則可導出其塑性矩陣［7］:

式中:λ=E+f2G，f=tanφ，δ=sign(τx'y')，E，G，υ分別為接觸面材料的彈性模量、剪切模量和泊松比。

當接觸縫面單元進入屈服或開裂狀態后，支配方程將成為非線性方程。此時，勁度矩陣K將是位移δ的函數。為了模擬加載過程，并加快收斂速度，本文對非線性方程的求解采用子增量變剛度迭代法，在每級荷載增量中剛度矩陣保持不變，并進行迭代計算。

2 數值模擬結果

2.1 算例1:受均布荷載的簡支梁的線彈性計算分析

如圖3所示，取一單位厚度的矩形截面簡支梁，頂部受均勻分布荷載作用。分布荷載的集度為10 kN/m2;梁的跨度為l=8 m;截面高度為h= 0.8 m。跨中截面為y軸(向下為正);中性層為x軸。梁材料的彈性模量E=2.0×104MPa，泊松比υ=0.20。計算時接觸面單元厚度人為取1.0× 10-3m，網格水平方向為40等分，豎直為8層等分，塊體單元的邊長為0.2 m×0.1 m。

利用材料力學［8-14］的方法計算梁跨中撓度的值，本文計算的最大撓度發生在跨中，其值為6.411 1×10-3m，與材料力學法結果的誤差為0.549%。塊體元變形計算結果略偏大的原因主要是人為引入接觸面單元降低了梁的剛度。

圖3 簡支梁網格與荷載約束圖

表1給出了靠1/4跨截面處(x=-1.1 m)塊體單元形心的x，y向正應力和切應力的計算值及對應的理論解［15］。比較塊體元結果與理論解可以發現兩者非常接近，誤差不超過1.0%，說明塊體元計算結果正確且精度足夠。

表2給出了靠1/4跨截面處(x=-1.0 m)接觸面單元形心的x，y向正應力和切應力的計算值及對應的理論解［10］。比較兩者可以發現:接觸面單元內的切應力及法向應力與理論解非常接近，誤差不超過1.0%;但接觸面方向(y向)的正應力與理論解有一定的誤差，變化規律相同，最大誤差為14.7%。這說明:塊體元法中接觸面內的切應力和法向應力計算結果可靠且精度足夠，接觸面內切向正應力誤差略大。

表1 截面x=-1.1 m塊體形心應力分量

表2 截面x=-1.0 m接觸縫面形心應力分量

2.2 算例2:受均布荷載的懸臂疊層梁屈服非線性分析

如圖4所示，單位厚度的矩形截面懸臂梁長6.0 m，截面高2.0 m，由4層等厚的板通過粘結膠合疊加而成。網格圖見圖10，板材料的彈性模量均取8.0 MPa，泊松比取0.22，接觸面材料彈性模量和泊松比的取值與板材一致。模擬接觸面2種特性:①粘結力取0 kPa，抗拉強度取2.0 MPa，摩擦角取25°;②粘結力取300 kPa，抗拉強度取2.0 MPa，摩擦角取25°。計算時僅考慮接觸面的非線性。

圖4 懸臂梁網格荷載約束圖

①屬于無粘結的疊層梁問題，接觸面上有法向約束力及摩擦產生的切應力。由于相對滑動，粘結面兩側的軸向正應力及位移不再連續。計算了5種荷載情況:P=30，80，120，150，200 kN/m2。圖5給出了懸臂端最大撓度與荷載的關系(撓度向下定義為負值)其中:“”為整體梁計算結果(不考慮接觸面的非線性);“”為疊層梁的計算結果。由變形結果可見:接觸面由于摩擦屈服，梁軸向線應變沿截面高度不再連續，梁的剛度明顯降低。當取均布荷載集度為150 kN/m2時，懸臂端的最大撓度分別為 -0.737 8 cm(疊梁)和-0.675 45 cm(整體梁)，疊梁的撓度增加了10.0%。圖6給出梁變形后的網格圖。從圖6可以明顯看出:粘結層處變形不連續，層與層之間有相互錯動。

圖5 懸臂梁端部結點位移分量與荷載圖

圖6 懸臂梁網格變形圖(P=150 kN/m2)

圖7、8給出了截面A(截面位置x=4.20 m)處結點在5種不同荷載情況下的軸向正應力和切應力，數據見表3。軸向正應力圖表明:由于粘結層的滑移，軸向應力在粘結層兩側不再連續;隨著荷載增加，突變量也隨之增加，在層內正應力基本呈線性變化。切應力圖表明:在粘結層摩擦屈服后，切應力降低，不再呈拋物線變化。

圖7 截面A結點正應力與荷載圖

圖8 截面A結點切應力與荷載圖

表3 截面A(x=4.02 m)結點應力分量

②屬于有粘結的疊合梁問題，接觸面上有法向約束力及摩擦和粘結產生的切應力。同樣計算了5種荷載情況:P=150，200，300，500，800 kN/m2。圖9給出了P=150 kN/m2時疊層梁的變形及粘結層的屈服區范圍。由圖9可見:粘結屈服區范圍及疊層梁變形比無粘結疊層梁小。由此得出:有粘結的疊層梁結構相對更穩定。

圖10、11給出了截面A(截面位置x=4.20 m)處單元形心在5種荷載情況下的軸向正應力和切應力。軸向正應力圖表明:粘結層一旦屈服，軸向應力在粘結層兩側不再連續，切應力也有所降低。

圖9 懸臂梁網格變形圖(P=150 kN/m2)

圖10 截面A塊體單元形心正應力圖

圖11 截面A塊體單元形心切應力圖

2.3 算例3:受集中荷載兩端固定疊層梁開裂非線性分析

將上述梁的兩端改為固定端約束，在第2層(從底到高編號)作用2個向下對稱的集中荷載F(見圖12)，網格示意圖見圖12。粘結面的抗拉強度取2.0 MPa，其他參數不變。隨著荷載的增加，簡支疊層梁中間的粘結層的最大法向應力接近材料抗拉強度，接觸面層開裂。圖12～14給出了3種不同荷載集度網格變形及粘結層的開裂情況。由圖可知:隨著荷載的逐漸增大，粘結層的開裂范圍相繼擴大。圖15給出梁跨中底部的豎向位移分量與荷載的關系。結果表明:粘結層一旦發生開裂，梁的擾度就會發生突變。

圖12 荷載作用位置與網格變形圖(P=2 000 kN)

圖13 荷載作用位置與網格變形圖(P=3 000 kN)

圖14 荷載作用位置與網格變形圖(P=4 000 kN)

圖15 懸臂梁跨中底部豎向位移分量與荷載關系圖

3 結束語

本文研究表明:在塊體的變形模式中引入應變和曲率后能顯著提高計算位移和應力的精度，且接觸面的應力也有足夠的精度，確保了模擬接觸面材料非線性力學行為的可靠性。

考慮接觸面材料的屈服和開裂后，非線性塊體單元法能較好地模擬有結構面材料的非線性變形及破壞規律，如疊層復合板、層狀巖體等。本文算例驗證了非線性塊體單元法的可靠性和合理性。如在塊體內引入材料的非線性，使塊體單元法模擬的結果更加接近實際。

［1］陳國榮.有限單元法原理及應用［M］.北京:科學出版社，2009.

［2］任青文，余天堂.塊體單元法的理論和計算模型［J］.工程力學，1999，16(1):67-77.

［3］任青文，余天堂.邊坡穩定的塊體單元法分析［J］.巖石力學和工程學報，2001，20(1):20-24.

［4］石根華.數值流形方法與非連續變形分析［M］.北京:清華大學出版社，1997.

［5］陸曉敏.邊坡的彈粘塑變形及穩定性研究［J］.巖石力學與工程學報，2002，21(4):493-496.

［6］殷有泉.非線性有限元基礎［M］.北京:北京大學出版社，2007.

［7］李詠偕，施澤華.塑性力學［M］.北京:水利電力出版社，1987.

［8］孫訓芳.材料力學［M］.北京:高等教育出版社，2003.

［9］TIMOSHENKO S P.高等材料力學［M］.北京:科學出版社，1979.

［10］陳杰.考慮層間接觸變形時疊層梁層間接觸壓力分析［J］.力學與實踐，2001(4):45-47.

［11］黃傳躍，諸德超.基于剖面翹曲修正理論的復合材料疊層梁動態特性分析［J］.航空學報，1996(5):96 -101.

［12］周一勤.夾層疊層梁的內力計算及其應用［J］.華東公路，1991(5):36-39.

［13］鄧梁波.地基上復合材料疊層梁的穩定性［J］.華南理工大學學報:自然科學版，1995(4):77-80.

［14］VIJAY K G.非守恒荷載作用下疊層梁的動力穩定性［J］.鋼結構，2009(2):74.

［15］徐芝倫.彈性力學［M］.北京:高等教育出版社，1988.

(責任編輯陳 艷)

Analysis of Stress and Deformation of Layered Beam Based on Block Element Method

LU Xiao-min，YUAN Tao
(Department of Engineering Mechanics，Hohai University，Nanjing 210098，China)

Based on block element method，deformation mode included rigid displacement，constant strains and curvatures was constructed.The equilibrium equations were conducted by principle of minimum potential energy.The displacement and stress under blanced load of simple beam were calculated.The results of block element method show that stress and displacement are very close to theoretical values only with fewer elements，and the stresses in interface elements have more accuracy yet alway it will be the upper limit of theoretical value.Besides，the stress accuraly is very high，which can ensule the reliability in simulating the damage of the contact surface.To analyze nonlinear material mechanics，yielding and cracking in interface elements were simulated.Yielded criterion usedM-Cformulation principle and cracked condition used the principle that the max-tensile stress in interfaceelement exceeds maximum tensile strength.Using the computer program of upper method，the stress and deformation of a layered cantilever beam were studied under uniform pressure.The results are correct and show that the method of the paper is feasible.

layered beam;block element method;numerical analysis

TB301

A

1674-8425(2015)11-0051-09

10.3969/j.issn.1674-8425(z).2015.11.009

2015-06-22

國家自然科學基金重點資助項目(11132003)

陸曉敏(1963—)，男，江蘇無錫人，博士，教授，主要從事水工結構、巖土工程等的力學分析及穩定性研究;袁濤(1989—)，男，江蘇鎮江人，碩士研究生，主要從事計算力學與工程仿真研究。

陸曉敏，袁濤.運用塊體單元法的疊層梁應力分析及破壞模擬［J］.重慶理工大學學報:自然科學版，2015 (11):51-59.

format:LU Xiao-min，YUAN Tao.Analysis of Stress and Deformation of Layered Beam Based on Block Element Method［J］.Journal of Chongqing University of Technology:Natural Science，2015(11):51-59.