基于修正形函數(shù)的Euler-Bernoulli開口裂紋梁單元剛度矩陣

徐 訓，朱亞杉，吳 浩

(武漢理工大學 土木工程與建筑學院，武漢 430070)

裂紋是土木工程結構最常見的損傷形式，裂紋會造成結構局部剛度降低，影響結構的動力特性。帶裂紋構件的受力特性可采用有限元軟件進行研究，如ABAQUS常用無厚度的Cohesive單元對裂紋進行模擬[1-2]，也可用指派裂紋或VCCT(virtual crack closure technique)等方法進行裂紋模擬[3-4]，上述方法在裂紋處需要進行網(wǎng)格加密以適應裂紋尖端應力場的奇異性，這將大幅增加計算自由度。工程中以結構固有頻率或模型修正等為代表的損傷識別方法要求裂紋模型能準確的描述損傷特性[5-7]，并具有模型簡單、參數(shù)少的特點。因此，以結構損傷識別為導向的裂紋梁有限元建模問題仍然是充滿挑戰(zhàn)的課題。

裂紋單元剛度矩陣的建立是有限元分析的關鍵。目前裂紋單元剛度建立方法主要包括基于應變能釋放的方法和傳統(tǒng)方法。在基于應變能釋放的方法中，裂紋產(chǎn)生的彈性應變能被加到非裂紋單元的彈性應變能中[8]，單元剛度矩陣可采用柔度計算法和直接剛度計算法兩種方法。柔度計算法是根據(jù)裂紋單元的應變能，應用卡式第二定理先推導柔度矩陣，再通過傳遞矩陣或直接求逆等方法得到裂紋單元的剛度矩陣，Papadopoulos等[9-11]對一維裂紋梁的剛度進行了研究。此外，文獻[12-14]對其在二維有限元中的應用也進行了討論。直接剛度計算法是直接利用卡式第一定理求得裂紋單元的剛度矩陣，該方法由Potirniche等[15]推導出二維裂紋單元，然后由Hall等[16]推導出帶有單邊裂紋的三維裂紋單元。Yazdi等[17]充分比較了基于應變能釋放的兩種方法，他們的研究表明用卡式第一定理生成裂紋單元剛度矩陣的不足之處，也證明了由Potirniche等[16]提出的方法在裂紋寬度較大時可能會導致解出現(xiàn)較大的誤差。

傳統(tǒng)方法是直接利用剛度折減[18]、單元截面參數(shù)折減[19]或彈性模量折減[20]等方法得到裂紋單元剛度矩陣。前者在非裂紋單元剛度矩陣中直接通過折減系數(shù)與剛度矩陣的表征裂紋對剛度矩陣的影響，無法表征裂紋深度和位置對剛度矩陣影響。后兩者的實質(zhì)是通過截面慣性矩、彈性模量的變化反映裂紋單元的抗彎剛度的變化，再通過抗彎剛度與形函數(shù)的二次導乘積的定積分得到裂紋單元剛度矩陣。如Shen等[21-23]通過指數(shù)型函數(shù)、階躍函數(shù)來描述抗彎剛度變化，這類方法認為形函數(shù)是始終不變的。

事實上，裂紋會影響裂紋附近一定距離的應力場變化，但這一變化不會直觀對應構件截面慣性矩的改變。且構件的彈性模量應只與材料相關，折減彈性模量物理意義不明確，則裂紋對單元剛度的影響可認為主要是改變了形函數(shù)。本文以Euler-Bernoulli梁為研究對象，在Euler-Bernoulli梁單元中，形函數(shù)的物理意義為兩端固定梁只產(chǎn)生一個單位位移(或轉角)時梁彎曲成的形狀[24]。有限元方法中，位移場用形函數(shù)和結點位移的乘積表示，裂紋對單元位移場的影響，必然導致裂紋單元形函數(shù)的改變。

有限元分析中通常利用位移場的改變?nèi)ッ枋鲂魏瘮?shù)的變化。Babuska等[25]利用無裂紋單元形函數(shù)的單位分解性質(zhì)提出了單位分解有限元法，通過附加位移自由度向量描述裂紋單元位移場，裂紋對位移場的影響通過增強函數(shù)和無裂紋單元形函數(shù)的乘積產(chǎn)生作用[26-27]。值得注意的是，這種方法丟失了常規(guī)有限元的插值屬性，會把有限元中的邊界條件應用復雜化，但增強函數(shù)這一概念極大促進了擴展有限元和廣義有限元方法的發(fā)展[28]。通過增加額外自由度，廣義有限元法有效描述了裂紋單元位移場函數(shù)變化[29]，并給出了廣義形函數(shù)的具體形式，但廣義有限元法需要設置與局部加強函數(shù)數(shù)目相等的額外自由度，這種依賴額外自由度的方法會造成線性相關性問題[30]，部分基于增加額外自由度的有限元方法還需考慮單元交界面的連續(xù)性問題[31]。針對該問題，Tian等[32]提出了無額外自由度廣義有限元法，該方法中結點自由度與標準有限元相同，增強函數(shù)直接附加在形函數(shù)上。

本文以單邊貫穿非擴展Euler-Bernoulli開口裂紋梁為研究對象，本文研究的最終目標是通過修正形函數(shù)，再利用抗彎剛度與修正后形函數(shù)的二次導乘積最終得到裂紋單元剛度矩陣。首先根據(jù)無額外自由度和增強函數(shù)的思想得到Euler-Bernoulli裂紋單元位移場，再通過位移場得到裂紋單元的形函數(shù)，最后利用虛位移原理導出Euler-Bernoulli裂紋的單元剛度矩陣。本文將討論裂紋單元形函數(shù)的性質(zhì)和所提單元剛度矩陣的有效性，研究結果將為基于靜力分析和振動分析的損傷識別方法提供新的研究思路。

1 Euler-Bernoulli裂紋梁單元位移場函數(shù)

單位分解有限元法通過在有限元插值函數(shù)中引入不連續(xù)函數(shù)來考慮裂紋的影響，其位移場函數(shù)uh(x)是傳統(tǒng)有限元位移場函數(shù)uFE(x)和增強位移函數(shù)uenh(x)之和[25]，如式(1)所示：

(1)

式中：ui是傳統(tǒng)有限元i結點位移自由度向量，Ni(x)是與i結點對應的形函數(shù)；ak是附加位移自由度向量，m為增強結點數(shù)，ψ(x)是增強k結點形函數(shù)的支撐域內(nèi)的增強函數(shù)。Nk(x)是與k結點對應的形函數(shù)。利用形函數(shù)單位分解的性質(zhì)，ψ(x)與形函數(shù)Nk(x)相乘可以全局或局部表征裂紋對uh(x)的影響。但增強函數(shù)ψ(x)依賴于附加自由度ak，附加自由度會導致求解規(guī)模增大等問題。針對此問題，Tian等提出了無額外自由度的有限元方法，其位移場函數(shù)如式(2)所示。

(2)

上述無額外自由度的有限元方法形函數(shù)為NiG(x)=∑Ni(x)∑ψ(x)，而增強函數(shù)ψ(x)需要由定義在該結點附近的單元結點進行構造，計算比較復雜。結合單位分解有限元單元增強函數(shù)和無額外自由度的想法，本文Euler-Bernoulli梁裂紋單元計算模型如圖1所示。

圖1中，le代表裂紋單元的長度，x0代表裂紋位置與單元i結點的距離(x0∈[0,le])，a代表裂紋的深度，h代表計算單元高度，b代表計算單元寬度，ui，θi，uj，θj分別代表i，j結點位移和轉角，Vi，Mi，Vj，Mj分別代表i，j結點的豎向力和彎矩，假定i，j結點處所示的位移與轉角均為正方向，增強函數(shù)ψ(x)只“附加”在本單元的i，j結點上，結合馮新等[33-34]的研究，梁單元位移場uh(x)用式(3)表示。

圖1 Euler-Bernoulli梁裂紋單元分析模型Fig.1 Euler-Bernoulli beam crack element analysis model

(3)

式中，F(xiàn)0=EI0表征未損傷構件的抗彎剛度，E為彈性模量，I0為完好梁截面的截面慣性矩。式中H(·)為Heaviside函數(shù)，C1，C2，C3，C4與梁的邊界條件相關。V是表征裂紋深度a的函數(shù)，根據(jù)文獻[33-35]，其具體形式如式(4)所示。

35.84s3+13.125s4)

(4)

式中，s=a/h，代表裂縫深度的歸一化參數(shù)。

將式(3)轉化為Euler-Bernoulli裂紋的撓度表達式如下

(F0C1Vx0-C2V)xH(x-x0)+

(5)

相較式(3)，式(5)清晰的表明位移場函數(shù)uh(x)是傳統(tǒng)非裂紋單元位移場Hermitian三次多項式F0(C1x3/6+C2x2/2+C3x+C4)與增強位移函數(shù)g(x)= (F0C1Vx0-C2V)xH(x-x0)-(F0C1Vx02-C2Vx0)H(x-x0)之和。可通過修改g(x)以適用多裂紋及更符合裂紋構件工程實際的情形。由式(5)知，如圖1所示的單元裂紋左端的位移場用Hermitian三次多項式擬合，而在右端位移場則由Hermitian三次多項式和g(x)共同擬合。g(x)與裂紋位置x0和裂紋深度a均有關系。值得注意的是，裂紋影響Hermitian三次多項式系數(shù)值，以適應增強函數(shù)在裂紋處出現(xiàn)的跳躍情形。

2 Euler-Bernoulli裂紋梁的形函數(shù)和單元剛度矩陣

2.1 Euler-Bernoulli裂紋梁的形函數(shù)

不引入額外自由度，根據(jù)圖1所示的計算模型，邊界條件可以寫為

x=0,uh(0)=ui,uh′(0)=θi

x=le,uh(le)=uj,uh′(le)=θj

(6)

假定i，j結點的位移ui，uj與轉角θi，θj為已知，則由上述邊界條件可以求出式(3)中的常數(shù)C1，C2，C3，C4如下

(7)

其中

(8)

(9)

Euler-Bernoulli梁非裂紋單元的形函數(shù)如式(10)所示。

(10)

式(10)和式(9)對比知，裂紋單元形函數(shù)曲線Nie(x)(i=1,2,3,4)是三次多項式曲線fi(x3,x2,x)與函數(shù)ψi(x)的疊加。fi(x3,x2,x)與Ni(x)具有相同的函數(shù)階次，但函數(shù)fi(x3,x2,x)的系數(shù)與裂紋位置x0和表征裂紋深度影響的參數(shù)V有關，其主要原因是適應函數(shù)ψi(x)在裂紋處出現(xiàn)的跳躍情形。函數(shù)ψi(x)主要表現(xiàn)裂紋對形函數(shù)曲線Nie(x)的附加影響，這一特征通過線性函數(shù)V(x-x0)和階躍函數(shù)H(x-x0)的乘積實現(xiàn)。

當x分別為0和le時，代入式(9)得到的計算結果和式(10)計算的結果一致，表明在單元邊界上裂紋引起的位移和轉角為0，式(9)的形函數(shù)并不改變式(10)的邊界值，自然滿足單元截面連續(xù)性，無需任何手段來解決單元交界面不連續(xù)問題。

值得注意的是，當裂紋深度a=0時，表征裂紋深度的參數(shù)V=0，式(9)和式(10)等價，式(10)只是式(9)的一個特例。即當裂紋深度a=0時，Nie(x)=Ni(x)(i=1,2,3,4;x∈[0,le])。

式(3)可進一步寫為

(11)

2.2 單元剛度矩陣

如圖1所示的Euler-Bernuolli裂紋梁單元，假定該單元在單元結點力Vi，Mi，Vj，Mj作用下處于平衡狀態(tài)。此時梁單元的撓曲線如式(11)所示，在此平衡狀態(tài)下，給該單元任一虛位移δu(x)，則單元i，j結點處產(chǎn)生虛位移δui，δθi，δuj，δθj。根據(jù)虛位移原理有

(12)

式中，δu(x)=δδeTNeT，I為截面慣性矩，將上式寫為矩陣形式如式(13)所示。

(13)

式中，F(xiàn)e=(Vi，Mi，Vj，Mj)，結合式(11)～(13)得到單元剛度矩陣的計算公式如式(14)所示。

(14)

本文認為裂紋對彈性模量E和截面慣性矩I沒有影響，式(14)中的抗彎剛度EI(x)為常數(shù)，再將式(9)帶入式(14)就可得到裂紋單元剛度方程。值得注意的是，根據(jù)Heaviside函數(shù)的數(shù)學性質(zhì)，式(14)計算過程中需要計算Dirac函數(shù)平方的積分。但Dirac函數(shù)平方的積分并沒有解析解，根據(jù)Dirac函數(shù)的定義，得到Dirac函數(shù)的乘積可表示如式(15)所示。

(15)

同樣，將式(15)與x2、x在單元內(nèi)的積分以及式(15)自身在單元內(nèi)的積分計算結果列為式(16)，并令三者的值分別為α1、α3和α2。

(16)

再令

(17)

最終得到裂紋梁單元各剛度項如式(18)所示。

(18)

Euler-Bernoulli梁裂紋單元剛度矩陣如式(19)所示

(19)

3 裂紋參數(shù)對形函數(shù)的影響

3.1 裂紋位置的影響

假定構件長L為43 mm，構件高h為10 mm。單元長度le=L。圖2給出了裂紋深度a=5 mm時，裂紋位置x0從i結點移動到j結點Nie(x)曲線變化情況。圖中x/L表征歸一化的位置(下同)，x0/L表示歸一化的裂紋位置(下同)。

3.2 裂紋深度的影響

同等構件長度和截面高度條件下，裂紋位置x0取10 mm，圖3給出了裂紋深度a從0變化到9 mm時Nie(x)曲線的變化情況，圖中a/h代表歸一化的裂紋深度。為方便查看Nie(x)曲線在不同裂紋深度a下的幅值變化范圍，在圖3中一并給出了其在平面xoNie上的投影情況。

4 動靜力分析驗證

為說明本文修改形函數(shù)的正確性和裂紋單元剛度矩陣的有效性，分別討論兩者在實際構件靜力分析和振動分析中的應用。

4.1 靜力撓度分析

結構撓度是橋梁檢測中常用的損傷識別指標，其主要原理是通過損傷前后結構撓度的變化對損傷進行定量和定性分析。本文將以二維純彎簡支梁和懸臂端承受集中荷載的懸臂梁為例來說明本文所提方法對結構撓度的擬合情況，簡支梁和懸臂梁計算模型如圖4所示，圖中xc為裂紋在構件中的位置。

圖4 二維簡支裂紋梁和懸臂裂紋梁模型Fig.4 The 2-D simply supported and 2-D cantilevered beams with crack

為說明本文的有效性，同樣定義撓度相對誤差Re表征擬合效果的優(yōu)劣。以ureal(x)表示通過有限元軟件ABAQUS仿真得到的構件撓度曲線，uiden(x)表示利用式(11)的得到的構件撓度曲線，其具體表達式如下。

(20)

如圖5所示，利用商業(yè)軟件ABAQUS建立了Euler-Bernuolli梁的二維有限元模型，并通過數(shù)值分析方法進行了計算。簡支梁和懸臂梁模型均采用CPS4R單元。其中，簡支梁模型整體尺寸為0.5 mm，懸臂梁網(wǎng)格的橫向尺寸為5 mm，豎向尺寸為1 mm。裂紋通過ABAQUS直接指派，考慮到裂紋尖端處應力場的奇異性，中間結點參數(shù)設為0.27，裂紋尖端退化單元設置為重復節(jié)點。裂紋部分的設置如圖5(c)所示。

(a) 簡支梁網(wǎng)格劃分模型

4.1.1 純彎簡支梁撓度分析

如圖4所示簡支梁，簡支長度L為43 mm，構件高度h為10 mm。彈性模量E=203.8 GPa，泊松比υ=0.3。兩端施加大小為1 075 N·m的彎矩M。裂紋位置xc分別取5 mm，10 mm，15 mm，21.5 mm。裂紋深度a分別取2 mm，3.5 mm，5 mm。

整個簡支梁作為一個單元，根據(jù)裂紋位置和深度由式(9)得到單元形函數(shù)Nie(x)，單元結點的轉角和位移由有限元軟件ABAQUS給出，根據(jù)式(11)得到本文撓度計算值uiden(x)。本文撓度計算值與有限元仿真撓度計算值比較如圖6所示。

(a) xc=5 mm

圖6表明不同裂縫位置處，隨著裂紋深度a的增大，有限元仿真值與本文計算值整體吻合越好，裂縫附近擬合效果一般，且整體來看本文計算值在裂紋附近范圍小于有限元仿真值。隨著裂紋位置xc的增加，有限元仿真值與本文計算值的擬合效果的變化情況不明顯。為準確描述擬合效果隨裂紋位置xc的變化，將上述工況的相對誤差計算結果列于表1中。

表1表明，同等裂紋深度下，越靠近支座處，有限元仿真撓度值和本文計算值的相對誤差越大，越靠近跨中，有限元仿真撓度值和本文計算值的相對誤差越小。當xc=5 mm時，相對誤差分別為1.443%，1.714%以及1.496%。而當裂紋位于跨中(xc=21.5 mm)時，相對誤差分別為0.862%，0.973%和0.731%。但表1相對誤差最大值為1.714%，可以認為本文值與仿真值吻合很好。

表1 各工況條件下簡支梁的撓度相對誤差計算Tab.1 RE of simply supported beam with different conditions

4.1.2 懸臂梁撓度分析

如圖4所示懸臂梁，該懸臂梁高h為20 mm。長L為300 mm，彈性模量E=203.8 GPa，泊松比v=0.3，裂縫距固定端距離xc分別為30 mm，80 mm，100 mm，150 mm，裂縫深度a分別為4 mm，7 mm，10 mm，在自由端承受大小為100 N的集中力，方向豎直向下。

整個懸臂梁作為一個單元，根據(jù)裂紋位置和深度由式(9)得到單元形函數(shù)Nie(x)，單元結點的轉角和位移由有限元軟件ABAQUS給出，根據(jù)式(11)得到本文撓度計算值uiden(x)，將本文撓度計算值與有限元仿真撓度計算值比較如圖7所示。

(a) xc=30 mm

由圖7知，隨著裂紋向懸臂端移動，裂紋對懸臂梁撓度的影響越小。整體來看，有限元仿真值與本文計算值整體吻合很好。同樣，將不同工況下的撓度值代入式(20)中量化，計算結果列于表2中。

表2 各工況條件下懸臂梁的相對誤差計算Tab.2 RE of cantilever beam with different conditions

整體來看，裂紋越深，本文方法擬合效果越好。越靠近支座處，擬合效果越差。與簡支梁不同的是，裂紋深度a為4 mm時，懸臂梁撓度曲線在裂紋處不再出現(xiàn)突變的情況，而是隨著裂紋位置的不同呈現(xiàn)出形狀相似、各點位移值差距小的特征。

a=4 mm(a/h=0.2)條件下，不同裂紋位置xc時，懸臂梁自由端位移最大差值為0.496 mm，而a=10 mm(a/h=0.5)條件下，其值為3.666 mm。這表明淺裂紋(a/h≤0.2)條件下，基于靜力撓度的損傷識別很難對結構損傷進行定位和分析，基于測點撓度變化分析的傳統(tǒng)損傷識別方法更適合裂紋深度更大的工程實際。

4.2 裂紋梁振動分析

為了驗證單元剛度矩陣的有效性，需要對Euler-Bernoulli梁單元振動分析。假定裂紋不引起質(zhì)量矩陣的變化，截面為等截面時，沿梁長方向梁線密度ρA為常量，單元長度為le，Euler-Bernoulli梁單元的質(zhì)量矩陣為：

(21)

根據(jù)裂紋幾何參數(shù)，由式(19)和式(21)分別求得單元剛度矩陣和單元質(zhì)量矩陣，將所有單元的單元剛度矩陣和質(zhì)量矩陣分別進行疊加，以形成總體剛度矩陣K和質(zhì)量矩陣M，結構的固有頻率可以通過下式求得

|K-ω2M|=0

(22)

式中，ω為結構的圓頻率，由圓頻率ω可得結構的自振頻率f如下

f=ω/2π

(23)

構件選取Sinha等用于試驗的懸臂鋁制梁，彈性模量E為69.79 GN/m2，質(zhì)量密度ρ為2 600 kg/m3，泊松比v為0.33，梁長L為1 m，寬b為0.05 m，高h為0.025 m。根據(jù)試驗過程，將懸臂邊界條件轉化為豎向彈簧和旋轉彈簧，其中豎向支撐剛度kt為26.6 MN/m，轉動剛度kθ為150 kNm/rad，計算模型如下圖所示。

圖8中不帶圈的數(shù)字代表結點號，帶圈數(shù)字代表單元號。計算工況根據(jù)Sinha等的工況選取以方便比較計算結果。同樣，懸臂梁分為16個單元，共34個自由度，得到不同裂縫深度條件下的鋁制懸臂梁的固有頻率如表3所示(表中試驗值來自文獻[23])。

圖8 懸臂裂紋梁有限元計算模型Fig.8 The FE model for the cantilever cracked beam

表3 不同裂紋深度下懸臂梁的固有振動頻率Tab.3 The natural frequencies of the cantilever beam with different depth

表4可知，a/h<0.4時，本文計算值和Sinha等的計算值小于結構的實際頻率；而a/h=0.48時，本文計算值和Sinha等的計算值大于結構的實際頻率；a/h<0.32時，本文求解一階固有頻率的誤差小于Sinha等利用抗彎剛度EI(x)折減得到一階固有頻率的誤差。當裂紋深度為a/h=0.48時，本文求解一階固有頻率的誤差為0.936%，大于Sinha等的計算結果0.863%，但本文一階固有頻率值的變化接近試驗測得的一階固有頻率值的變化。因此可以認為本文求解結構固有頻率在a/h<0.5時計算結果準確，計算結果較文獻[23]的方法更好。

表4 懸臂梁的第一階固有頻率兩種求解方法對比Tab.4 Comparison of two methods for solving 1st natural frequency of cantilever beam

下面具體討論裂紋位置和深度對裂紋構件固有頻率的基本影響以分析懸臂梁的自由振動特性和剛度矩陣的正確性，以fc表示裂紋構件的固有頻率，fib表示完整構件的固有頻率，設置參數(shù)fr表征裂縫位置和深度對構件頻率的影響，表達如下。

fr=1-fc/fib

(24)

圖9給出了上述懸臂構件的前三階固有頻率隨著裂紋位置、深度的變化情況。

(a) 懸臂裂紋梁的一階固有頻率變化率

圖9表明在裂紋位置相同條件下裂紋越深，其fr值越大，固有頻率降低越明顯。結構一階固有頻率呈現(xiàn)單調(diào)變化，裂紋距離固定端越近，裂紋越深fr值越大，反之越小。而結構二階、三階固有頻率的fr值并不隨著裂紋位置的改變呈單調(diào)性變化，二階固有頻率fr值出現(xiàn)兩處局部最大值，三階固有頻率fr值出現(xiàn)三處局部最大值。上述計算結果與已有文獻[35]中對懸臂構件前三階段的振動分析結論保持一致，證明了本文裂紋單元剛度矩陣的正確性和有效性。

5 結 論

本文根據(jù)單位分解有限元思想，把增強函數(shù)附加在單元結點上，在不增加自由度的基礎的給出了Euler-Bernoulli裂紋梁單元形函數(shù)和剛度矩陣的具體形式，并進行了靜力分析驗證和振動分析驗證，主要結論如下：

(1) 該形函數(shù)能直觀表現(xiàn)裂紋位置和深度的影響。裂紋在該形函數(shù)中表現(xiàn)為裂紋處幅值的突變，裂紋位置、深度與該形函數(shù)呈非線性關系。且該形函數(shù)自動滿足單元協(xié)調(diào)條件，當裂紋深度為0時，該形函數(shù)退化為傳統(tǒng)形函數(shù)。

(2) 通過簡支開口裂紋梁和懸臂開口裂紋梁數(shù)值仿真，該形函數(shù)能準確擬合單元位移場。對于裂紋構件，裂紋深度越大，裂紋位置越遠離結點，擬合效果就越好。

(3) 通過懸臂開口裂紋梁自由振動分析，本文所提的剛度矩陣模型能準確的描述開口裂紋對結構頻率的影響。