含多裂紋損傷圓弧曲梁自由振動擾動的有限元網(wǎng)格自適應(yīng)分析

王永亮，王建輝，張 磊

(1. 中國礦業(yè)大學(xué)(北京)，力學(xué)與建筑工程學(xué)院，北京 100083；2. 中國礦業(yè)大學(xué)(北京)，煤炭資源與安全開采國家重點實驗室，北京 100083)

圓弧型曲梁作為基本構(gòu)件，廣泛應(yīng)用于土木工程、機械工程、航空航天工程等領(lǐng)域中[1-3]。曲梁在工程實際中常帶裂紋損傷，準確評估帶裂紋損傷曲梁的動力性能是結(jié)構(gòu)設(shè)計的重要考慮因素；裂紋損傷深度、數(shù)目和分布均會改變曲梁基本特性，擾動梁的頻率和振型[4-6]，明確裂紋損傷對動力性能的影響，可以有效確保結(jié)構(gòu)的安全使用和有針對性的加固改造。同時，利用含損傷曲梁的實際自振頻率和振型可以進行裂紋識別和定位[7-8]，裂紋損傷深度、數(shù)目、位置的精準識別依賴于高精度的頻率和振型解答[9-10]。為了獲得梁構(gòu)件的高精度自由振動解答，一些解析方法和理論模型得到發(fā)展，但仍難有效應(yīng)用于變曲線線型、多裂紋損傷、各類邊界條件等復(fù)雜工況[1,11-12]。

數(shù)值計算成為分析復(fù)雜結(jié)構(gòu)動力性能的合理選擇和重要技術(shù)，有限元法被發(fā)展和應(yīng)用于求解含裂紋損傷曲梁的自振頻率和振型[13-15]，但解答精度依賴于網(wǎng)格劃分質(zhì)量，解答因網(wǎng)格劃分難免引入誤差[16]。特別是裂紋損傷形成各階振型的擾動影響，對于非均勻分布網(wǎng)格的有效性提出較高要求。有限元網(wǎng)格自適應(yīng)分析方法可有效地優(yōu)化網(wǎng)格分布，在直線梁彈性屈曲[17]、曲梁振動[18-19]、板殼振動[20]、含損傷梁振動[21]和屈曲[22]、巖體變形和斷裂[23]等問題求解中展示出很好的求解效力。本文將建立圓弧形曲梁裂紋的截面損傷缺陷比擬方案，進行裂紋深度、位置、數(shù)目的模擬，引入變截面Timoshenko 梁的h型有限元網(wǎng)格自適應(yīng)分析方法[18-19]，求解含裂紋損傷圓弧曲梁自由振動問題，得到優(yōu)化的網(wǎng)格和滿足預(yù)設(shè)誤差限Tol的高精度連續(xù)階自振頻率和振型。文中給出求解含多裂紋損傷圓弧曲梁自由振動數(shù)值算例，利用得出的有無損傷頻率差、振型差分析了多裂紋損傷深度、數(shù)目、分布對圓弧曲梁自振頻率和振型的擾動影響。

1 圓弧曲梁裂紋損傷表征方法

考慮圖1 所示含裂紋損傷平面曲梁，曲梁中性軸坐標為s，坐標系為xyz，其中x、y為曲梁平面內(nèi)坐標，x沿軸線切向，y沿軸線法向，z垂直于軸線所在平面。面內(nèi)振動的位移為：沿x軸位移振幅u、沿y軸位移振幅v和繞z軸的轉(zhuǎn)角振幅ψz。記曲梁曲率半徑為R(s)，截面剪切剛度修正系數(shù)為κ，截面面積為A(s)，對z軸慣性矩為I(s)，長度為l，梁高度為h，梁厚度為b。記材料彈性模量為E，剪切模量為G，泊松比為ν，密度為 ρ。

圖1 含裂紋損傷曲梁坐標系和符號Fig. 1 Coordinate systems and symbols of cracked curved beam

本文研究曲梁中的微裂紋損傷，使得梁截面產(chǎn)生弱化、梁的截面屬性衰減。本研究采用裂紋截面損傷缺陷比擬方法[21-22]，裂紋處截面損傷定義為：

式中，Tol為自由振動解答的預(yù)設(shè)誤差限。

2 圓弧曲梁自由振動

本文研究的平面曲梁面內(nèi)自由振動，該特征值問題的微分控制方程為[18-19,24]：

3 網(wǎng)格自適應(yīng)細分加密

有限元計算存在相比當前網(wǎng)格解答具有更高收斂階的超收斂點[25]，利用超收斂點結(jié)合單元拼片、高階形函數(shù)插值技術(shù)，可以提高當前有限元解的精度，得到全域的超收斂解[21,26-27]。本文對于圓弧曲梁的自由振動問題，求得當前網(wǎng)格下振型(位移)的有限元解后，利用有限元后處理超收斂拼片恢復(fù)方法，得到振型的超收斂解：

式中：P為給定函數(shù)向量；a為待定系數(shù)向量。隨后，利用振型解答并通過Rayleigh 商計算可以獲得自振頻率值[28]。引入振型超收斂解，可對當前網(wǎng)格下振型有限元解進行能量模形式下的誤差估計[21, 25]：

利用振型誤差估計，網(wǎng)格可以進行優(yōu)化處理來降低和控制振型的誤差，達到預(yù)設(shè)的解答精度。本文方法對每個有限元單元e上的振型誤差進行判斷，如果誤差控制式(10)不滿足，則表明該單元上振型解答的誤差過大，需要通過進行網(wǎng)格優(yōu)化處理，本文采用單元均勻細分加密的h型網(wǎng)格自適應(yīng)方式來增加模型自由度、降低單元上解答的誤差[21]。當前單元細分生成的新單元長度與目前誤差和單元階次相關(guān)，即利用當前誤差可以估計新單元的長度：

綜合以上各方法可形成如下整體計算分析方案，獲得含裂紋損傷圓弧曲梁的各階頻率和振型高精度解答：

1)含裂紋損傷圓弧曲梁模型。利用裂紋損傷表征方法(式(1)～式(5))，模擬多裂紋在曲梁中的深度、數(shù)目、分布，形成含裂紋損傷圓弧曲梁模型。

2)當前網(wǎng)格下頻率和振型有限元解。在當前有限元網(wǎng)格下，利用曲梁自由振動問題的有限元逆冪迭代分析方法(式(6)～式(8))，求解含裂紋損傷圓弧曲梁模型，得到頻率和振型的有限元解答。

3)誤差估計并加密更新有限元網(wǎng)格。利用網(wǎng)格自適應(yīng)細分加密方法，對當前振型解答進行誤差估計，不滿足預(yù)設(shè)誤差限Tol時，在裂紋損傷擾動振型區(qū)域進行網(wǎng)格細分加密，獲得更新的加密網(wǎng)格(式(9)～式(11))。在更新的有限元網(wǎng)格下，返回步驟2)、步驟3)進行循環(huán)計算和誤差估計，直到獲得一套充分優(yōu)化的網(wǎng)格和滿足誤差限的解答。

4 數(shù)值算例

本文方法已經(jīng)編制相應(yīng)的Fortran 90 語言程序代碼，程序開發(fā)實施基于Microsoft Visual Studio和Intel Visual Fortran 編程軟件平臺。本節(jié)給出求解具有代表性的多種含裂紋圓弧曲梁自由振動數(shù)值算例，對網(wǎng)格自適應(yīng)劃分以及頻率、振動解答的精確性進行討論；對裂紋損傷深度、數(shù)目、分布等因素影響自振頻率和振型擾動進行了分析，檢驗了本文算法的可靠性和實用性。本節(jié)所有算例均采用3 次元，初始網(wǎng)格采用2 個單元，給定的初始誤差限為Tol=10-4。

例1. 無裂紋損傷圓弧曲梁

為檢驗本文方法求解無損傷圓弧曲梁的精確性和有效性，本研究對圖2 所示兩端固定的常截面1/4 圓弧曲梁進行求解。為便于檢驗計算結(jié)果的數(shù)值精度，該算例采用無量綱的純數(shù)值計算，曲梁的基本幾何與物理參數(shù)如下：

圖2 無裂紋損傷1/4 圓弧曲梁模型Fig. 2 Model of a quarter of uncracked circularly curved beam

表1 所示為使用本文方法求解得到的前5 階自振頻率解答，將文獻[29]中結(jié)構(gòu)力學(xué)求解器采用9000 個常截面直線型單元進行求解得到的高精度解答進行對比分析，同時給出了各階求解使用的最終單元數(shù)目和頻率誤差，可見本文方法在自適應(yīng)網(wǎng)格下得到的各階解答遠小于預(yù)設(shè)誤差限要求。需要指出的是，本研究對振型進行誤差控制，使用振型解答并通過Rayleigh 商計算得出具有更高收斂階的頻率[28]，確保頻率值亦能嚴格滿足誤差限。

表1 無裂紋損傷1/4 圓弧曲梁自振頻率值Table 1 Natural frequencies of a quarter of uncracked circularly curved beam

圖3 所示為使用本文方法求解得到的第1 階、第5 階振型解答，并在橫坐標軸上標記出自適應(yīng)網(wǎng)格的最終分布情況。為方便直觀顯示和對比分析，圖中振型結(jié)果均進行歸一化處理(令最大振型值為1)?？梢钥闯?，振型在兩固定端的位移均為0 值；本文方法求解各階振型均劃分出非均勻網(wǎng)格，且在振型變化平緩區(qū)域使用稀疏網(wǎng)格、在振型變化劇烈處采用了相對細密的網(wǎng)格，避免了全域使用一致細密網(wǎng)格的冗余性。隨著階次的增加，振型復(fù)雜程度增強，第5 階振型比第1 階使用了更多的單元。

圖3 裂紋損傷1/4 圓弧曲梁振型Fig. 3 Vibration modes of a quarter of uncracked circularly curved beam

例2. 單一裂紋損傷曲梁不同裂紋深度

為檢驗本文方法分析含不同裂紋深度圓弧曲梁自由振動問題的有效性，本研究對圖4 所示兩端簡支的單裂紋損傷曲梁進行求解。該曲梁裂紋深度分別取為α=0.0、0.16、0.5，裂紋位置取為β=0.6125。該圓弧曲梁夾角為120°，裂紋位置角度為73.5°，其余的基本幾何與物理參數(shù)如下：

圖4 單裂紋損傷曲梁不同裂紋深度(α=0.0、0.16、0.5，β=0. 6125)模型Fig. 4 Model of curved beam with single crack in different depth cases (α=0.0, 0.16, 0.5, β =0.6125)

使用本文方法分別計算了該曲線梁在三種裂紋損傷深度工況下面內(nèi)自由振動的連續(xù)前5 階特征對，計算頻率值列于表2。文獻[14]結(jié)合能量方法和有限元模型、文獻[6]采用物理模型對上述問題進行分析，得到頻率值如表2 所示。通過對比本文方法和能量方法、物理模型求解結(jié)果，可以看出隨著問題復(fù)雜程度增加(如階次增加、裂紋深度增大)，二者因基本分析方法不同導(dǎo)致個別階次的本文方法與能量方法解答誤差(下劃線標出)略有增大，在各類裂紋深度工況的其余階次下均展示出良好的一致性。

表2 單裂紋損傷曲梁不同裂紋深度自振頻率值Table 2 Natural frequencies of curved beam with single crack in different depth cases

圖5 給出了本文方法求解裂紋深度對自振頻率擾動影響結(jié)果。為方便直觀顯示和對比分析，圖中給出頻率值(左側(cè)縱軸標識刻度)和頻率差(右側(cè)縱軸標識刻度，含裂紋損傷情況頻率值與無裂紋損傷情況頻率值的差值)。可以看出，頻率值隨階次增加，沒有出現(xiàn)顯著性差異；頻率差均為負值，可知裂紋損傷的出現(xiàn)降低了各階頻率值；α=0.5 時的各階頻率差降低幅度均比α=0.16 大，可見裂紋損傷深度越大，則梁截面產(chǎn)生弱化、梁截面屬性衰減程度越大，表現(xiàn)為頻率值的顯著減低。

圖5 裂紋深度對自振頻率擾動影響Fig. 5 Disturbance influence of crack depth on natural frequencies

圖6 所示為使用本文方法求解得到單裂紋損傷深度α=0.5 下第1 階、第5 階振型解答?？梢钥闯觯裥驮诹鸭y損傷附近區(qū)域出現(xiàn)擾動，裂紋損傷對轉(zhuǎn)動位移 ψz擾動最為明顯，本文方法求解自適應(yīng)劃分出最終非均勻的網(wǎng)格，在裂紋附近區(qū)域使用了相對密集的網(wǎng)格來適應(yīng)裂紋損傷引起振型的變化，體現(xiàn)了本文方法自適應(yīng)劃分網(wǎng)格對各階振型變化的適應(yīng)性。

圖6 單裂紋損傷曲梁裂紋深度α=0.5 下振型Fig. 6 Vibration modes of curved beam with single crack in depth case α=0.5

為分析裂紋損傷程度對振型的擾動行為，圖7所示為裂紋損傷α=0.16 時各振型與相應(yīng)階次無損傷振型的差值曲線，可以看出裂紋損傷所在局部區(qū)域?qū)φ裥妥兓酗@著影響，裂紋損傷是影響各振型分量擾動變化的主要因素。在本例工況條件下，裂紋損傷對各振型分量均出現(xiàn)擾動，其中轉(zhuǎn)動位移 ψz擾動最大。

圖7 單裂紋損傷曲梁裂紋深度α=0.16 下振型擾動Fig. 7 Vibration modes disturbance of curved beam with single crack in depth case α=0.16

圖8 所示為裂紋損傷α=0.5 時各振型與相應(yīng)階次無損傷振型的差值曲線，振型差值相比α=0.16時各振型更大，裂紋損傷程度越大則振型擾動愈加劇烈。通過上述結(jié)果可以看出，振型差幅值與損傷程度相關(guān)，通過定量控制裂紋損傷量，可有效控制振型擾動。本算例檢驗了本文方法求解含裂紋損傷圓弧曲梁解答的精確性，以及對各類裂紋損傷深度問題的適用性。

圖8 單裂紋損傷曲梁裂紋深度α=0.5 下振型擾動Fig. 8 Vibration modes disturbance of curved beam with single crack in depth case α=0.5

例3. 多裂紋損傷曲梁不同裂紋數(shù)目

含裂紋損傷曲梁的裂紋數(shù)目n是影響裂紋振動特性的又一重要因素，本研究設(shè)置如表3 所示的典型多裂紋數(shù)目(n=2、3、4)和裂紋位置工況。本例采用圖4 所示的兩端簡支圓弧曲梁，該曲梁的幾何模型和基本物理參數(shù)同式(13)。

表3 多裂紋損傷數(shù)目和位置工況Table 3 Number and location of multiple cracks damage

各工況中多裂紋為圖9 所示的均勻分布形式，各裂紋間夾角分別為40°(工況Ⅰ，n=2)、30°(工況Ⅱ，n=3)、24°(工況Ⅲ，n=4)。

圖9 多裂紋損傷曲梁不同裂紋數(shù)目(n=2、3、4)模型Fig. 9 Model of curved beam with multiple cracks in different number cases (n=2, 3, 4)

使用本文方法分別計算了該曲線梁在三種裂紋損傷數(shù)目工況下面內(nèi)自由振動的連續(xù)前50 階特征對，遴選典型的計算頻率值列于表4?？梢钥闯?，頻率值隨裂紋數(shù)目的增加，除個別階次(下劃線標出)頻率值略有增加，整體上呈現(xiàn)逐漸降低的趨勢。

表4 多裂紋損傷曲梁不同裂紋數(shù)目自振頻率值Table 4 Natural frequencies of curved beam with multiple cracks in different number cases

圖10 給出了本文方法求解裂紋數(shù)目對前5 階自振頻率擾動影響結(jié)果，可以看出頻率值隨階次增加，各裂紋數(shù)目工況沒有出現(xiàn)顯著性差異；裂紋損傷數(shù)目的增加降低了各階頻率值，一般情況下，裂紋數(shù)目越多則降低程度越大；在第3 階、第4 階時，出現(xiàn)頻率差為正值的情況，即裂紋損傷數(shù)目增多反而提高頻率值；在第4 階時，出現(xiàn)裂紋數(shù)目為4 工況的頻率差比裂紋數(shù)目為3 工況的頻率差更小的現(xiàn)象，即裂紋增多并沒有顯著降低頻率值。綜合以上結(jié)果，可知多裂紋數(shù)目與位置同時影響頻率值，增加裂紋數(shù)目整體上有降低各階頻率的趨勢，但因為裂紋位置的改變在某些階次上會出現(xiàn)頻率值增加的現(xiàn)象。因此，在原多裂紋損傷位置基礎(chǔ)上，繼續(xù)增加新的裂紋或增大原有裂紋深度(如本文例2)，才能出現(xiàn)各階均降低的頻率值。

圖10 裂紋數(shù)目對自振頻率擾動影響Fig. 10 Disturbance influence of crack number on natural frequencies

圖11 所示為使用本文方法求解得到多裂紋損傷曲梁不同裂紋數(shù)目下首階振型擾動解答?？梢钥闯?，振型在各多裂紋損傷附近區(qū)域出現(xiàn)擾動，轉(zhuǎn)動位移 ψz擾動最為顯著，本文方法求解自適應(yīng)劃分出最終非均勻的網(wǎng)格，在裂紋附近區(qū)域使用了相對密集的網(wǎng)格來適應(yīng)裂紋損傷引起振型的變化，體現(xiàn)了本文方法自適應(yīng)劃分網(wǎng)格對多裂紋損傷曲梁各階振型變化的適應(yīng)性。

圖11 多裂紋損傷曲梁不同裂紋數(shù)目振型擾動Fig. 11 Vibration modes disturbance of curved beam with multiple cracks in different number cases

例4. 多裂紋損傷曲梁不同裂紋分布

為進一步分析多裂紋損傷分布對曲梁自由振動的影響，本例仍采用圖4 所示的兩端簡支圓弧曲梁，該曲梁的幾何模型和基本物理參數(shù)同式(13)。該曲梁設(shè)置5 條裂紋損傷，考慮裂紋損傷沿曲梁均勻分布(各裂紋位置β 為0.1667、0.3333、0.5000、0.6667、0.8333)、裂紋損傷集中于曲梁左側(cè)集中分布(各裂紋位置β 為0.0833、0.1667、0.2500、0.3333、0.4167)2 種工況，各工況的多裂紋如圖12 所示，各裂紋間夾角分別為20°(裂紋均勻分布)、10°(裂紋左側(cè)集中分布)。

圖12 多裂紋損傷曲梁不同裂紋分布模型Fig. 12 Model of curved beam with multiple cracks in different distribution cases

使用本文方法分別計算了該曲線梁在不同裂紋分布工況下面內(nèi)自由振動的連續(xù)前50 階特征對，遴選典型的計算頻率值列于表5，裂紋分布對自振頻率擾動影響如圖13 所示?？梢钥闯?，裂紋左側(cè)集中分布相比均勻分布，頻率值在低階(如第1 階～第4 階)時具有更高的數(shù)值，而在高階時具有更低的數(shù)值，相同數(shù)目多裂紋的不同分布形式成為影響振動特性的重要因素。因此，需要同時精準檢測出裂紋損傷數(shù)目和各裂紋位置，才能準確估計含裂紋損傷時的頻率值。

圖13 裂紋分布對自振頻率擾動影響Fig. 13 Disturbance influence of crack distribution on natural frequencies

表5 多裂紋損傷曲梁不同裂紋分布自振頻率值Table 5 Natural frequencies of curved beam with multiple cracks in different distribution cases

圖14 所示為使用本文方法求解得到多裂紋損傷曲梁不同裂紋數(shù)目下首階振型擾動解答?？梢钥闯?，振型在各均勻分布和集中分布的裂紋損傷附近區(qū)域出現(xiàn)擾動，本文方法在裂紋附近區(qū)域使用了相對密集的網(wǎng)格來適應(yīng)裂紋損傷引起振型的變化，體現(xiàn)了本文方法自適應(yīng)劃分網(wǎng)格均能很好適應(yīng)不同裂紋密集程度誘發(fā)的振型擾動。需要指出的是，在圖14(b)所示的裂紋左側(cè)集中分布工況下，振型在右側(cè)區(qū)域同時出現(xiàn)較大幅值的振型擾動，體現(xiàn)了多裂紋損傷對整體振型的強擾動行為。

圖14 多裂紋損傷曲梁不同裂紋分布振型擾動Fig. 14 Vibration modes disturbance of curved beam with multiple cracks in different distribution cases

5 結(jié)論

本文建立圓弧形曲梁裂紋的截面損傷缺陷比擬方案和h型有限元網(wǎng)格自適應(yīng)分析方法，求解含裂紋損傷圓弧曲梁自由振動問題，得到優(yōu)化的網(wǎng)格和滿足預(yù)設(shè)誤差限的高精度自振頻率和振型解答，定量研究多裂紋損傷深度、數(shù)目、分布形式等對圓弧曲梁自振頻率和振型的擾動行為。本文的主要結(jié)論如下：

(1)裂紋區(qū)域自適應(yīng)網(wǎng)格。自適應(yīng)網(wǎng)格算法對無損傷、含損傷曲梁分析具有良好適用性，振型在裂紋損傷附近區(qū)域出現(xiàn)擾動，本文自適應(yīng)優(yōu)化出非均勻網(wǎng)格，在裂紋附近區(qū)域使用了相對密集的網(wǎng)格來適應(yīng)裂紋損傷引起振型的變化。

(2)裂紋損傷深度。裂紋損傷的出現(xiàn)降低了各階頻率值，損傷深度越大則降低程度越大；裂紋損傷對轉(zhuǎn)動位移 ψz擾動最大，損傷程度越大越加劇擾動幅值。

(3)多裂紋損傷數(shù)目。裂紋數(shù)目與位置同時影響頻率值，增加裂紋數(shù)目整體上有提高頻率的趨勢，但因為裂紋位置的改變在某些階次頻率上亦會降低頻率值，各裂紋損傷附近區(qū)域的振型均出現(xiàn)擾動。

(4)多裂紋損傷分布。裂紋一側(cè)集中分布相比均勻分布，頻率值在低階時具有更高的數(shù)值，而在高階時具有更低的數(shù)值；振型在各均勻分布和集中分布的裂紋損傷附近區(qū)域均出現(xiàn)擾動，相同數(shù)目多裂紋的不同分布形式成為影響振動特性的重要因素。