杜軍虎， 韓 丁

（合肥工業(yè)大學交通運輸工程學院，安徽合肥 230009）

對于瀝青混合料疲勞性能的評價，傳統(tǒng)方法是通過大量的室內(nèi)試驗獲取相關(guān)數(shù)據(jù)，并與相應(yīng)的規(guī)范指標進行比對來分析瀝青混合料在不同條件下的疲勞性能［1－3］。由于在試件成型過程中存在一定的誤差，試驗結(jié)果的離散性較大，加之這類試驗耗時較長，使得工作量變大且成本提高。相比之下，用有限元軟件仿真來替代真實試驗有更大的優(yōu)越性。文獻［4］利用擴展有限元法研究了有限元模擬裂紋尖端網(wǎng)格劃分的方法；文獻［5］用有限元法對瀝青路面在初始裂紋形成后的裂紋擴展規(guī)律進行了研究；文獻［6］利用非線性有限元方法分析了瀝青混合料小梁疲勞試驗中彎拉應(yīng)力、撓度、損傷變量、裂紋擴展速率等的變化規(guī)律，并預(yù)測了瀝青混合料試件的疲勞壽命及失穩(wěn)斷裂時的裂紋長度；文獻［7］基于黏彈性連續(xù)損傷理論用有限元分析法預(yù)測了瀝青混凝土路面的疲勞性能；文獻［8］通過ABAQUS有限元模擬，分析了級配碎石夾層對基層反射裂縫的影響。在有限元軟件仿真的過程中，結(jié)果的精確性是主要考慮的問題之一。在確定工況下，有限元模型網(wǎng)格密度的大小決定了計算結(jié)果的精確性［9］。一般來說，網(wǎng)格越細化，得到的計算結(jié)果越精確，但過密的網(wǎng)格劃分會以較長的計算時間為代價。確定出合理的網(wǎng)格密度，可使有限元的仿真結(jié)果既能滿足精確性的要求，又具有合理的計算成本。

本文對瀝青混合料小梁彎曲破壞試驗進行有限元仿真；應(yīng)用線彈性斷裂力學理論，計算瀝青混合料小梁試件梁底張開型裂縫對應(yīng)的應(yīng)力強度因子KⅠ值［10－12］和J積分值［13－14］；獲取一定試驗工況時3種幾何尺寸模型在不同網(wǎng)格密度下的KⅠ值及J積分值，分析瀝青混合料小梁試件有限元模型的網(wǎng)格密度對計算結(jié)果的影響，確定對應(yīng)不同尺寸試件的合理網(wǎng)格密度，最后利用Matlab軟件對所得數(shù)據(jù)進行交互界面數(shù)據(jù)庫的建立，為使用者提供方便直觀的查詢。

1 小梁彎曲試驗的有限元模擬

1.1 模擬工況

利用有限元軟件對50 mm×50 mm× 250 mm、50 mm×50 mm×300 mm、100 mm× 100 mm×300 mm 3種幾何尺寸小梁試件的彎曲破壞進行模擬，小梁試樣重量分別為0.014 4、0.018、0.072 k N，采用三維有限元模型，模擬三分點加載模式，如圖1所示。

圖1 有限元模型及試件尺寸

為了分析網(wǎng)格劃分密度對仿真結(jié)果的影響，固定加載值大小為8.4 k N，底部兩支撐位置距小梁側(cè)面的距離為整個小梁長度的3／16，通過約束支撐位置x、y方向位移自由度來模擬對小梁的支撐。加載位置為支撐間三分點處。

1.2 網(wǎng)格劃分及單元類型確定

對于小梁模型的網(wǎng)格劃分策略，采用結(jié)構(gòu)化網(wǎng)格劃分方法。在通常情況下，小梁產(chǎn)生開裂的位置為底部中間位置，因此該部位是著重關(guān)心的區(qū)域。如果采用常規(guī)單元來劃分裂縫尖端附近區(qū)域，不論尖端附近的網(wǎng)格劃分有多細，計算精度都不會高。一些用于模擬裂縫尖端附近特殊應(yīng)力、應(yīng)變場的特殊單元隨著有限元的發(fā)展相繼出現(xiàn)。目前，簡潔、常用的辦法是直接采用1／4節(jié)點單元（奇異單元），即在裂縫尖端處將中間結(jié)點向裂縫尖端靠攏，距裂縫尖端1／4邊長處，如圖2所示，這樣的單元可以較好地反映裂縫尖端附近的應(yīng)力場。因此，在底中部取10 mm×20 mm的區(qū)域進行網(wǎng)格細化，并采用奇異單元（圖1標注位置）模擬虛擬裂縫尖端的奇異性［15］，以虛擬裂縫尖端為中心向外共設(shè)置3個圍道。對于其他區(qū)域，采用普通的六面體單元，并設(shè)置相對較大的網(wǎng)格尺寸，在試件寬度方向劃分10層網(wǎng)格單元。這樣劃分網(wǎng)格既考慮到了小梁底部中間位置計算結(jié)果的精確性，同時又能避免對不重要區(qū)域的大量無用計算。

圖2 常規(guī)單元與奇異單元的比較

除了網(wǎng)格的劃分策略，在有限元的模擬中，選用的單元類型對于計算結(jié)果同樣有很大影響。通常二次完全積分單元對應(yīng)力的計算結(jié)果較準確，適用于模型應(yīng)力集中問題，且一般情況下不會出現(xiàn)一次積分單元常見的剪力自鎖及沙漏等問題［16］。采用奇異單元對瀝青混合料小梁的斷裂力學指標進行計算時，虛擬裂縫位置會產(chǎn)生應(yīng)力集中現(xiàn)象［17］，因此適合選用二次完全積分單元。

1.3 計算結(jié)果

在一定的模擬試驗工況以及網(wǎng)格劃分方法下，網(wǎng)格密度是影響輸出結(jié)果的主要因素。為了尋求合理的網(wǎng)格密度，對3種不同尺寸小梁模型的底中部細化區(qū)域分別取12、24、48、96、192個／cm35種不同的網(wǎng)格密度，進行有限元仿真，其中3種網(wǎng)格密度下得到的計算云圖經(jīng)局部放大后如圖3所示。

圖3 不同網(wǎng)格密度下的計算云圖

在有限元模型的奇異單元區(qū)域沿試件寬度方向由前向后設(shè)置了21層圍道，每層圍道由內(nèi)向外又分別設(shè)置了3個圍道。為了使計算值更具代表性，取中間層圍道第2圍道的中間點作為數(shù)據(jù)輸出節(jié)點。通過計算，得到了3種尺寸（長×寬×高）數(shù)值試件在不同網(wǎng)格密度下該節(jié)點的應(yīng)力強度因子KⅠ值及J積分值，結(jié)果見表1所列。

表1 不同網(wǎng)格密度下的KⅠ值及J積分值計算結(jié)果

2 網(wǎng)格密度對仿真結(jié)果的影響分析

通過表1中數(shù)據(jù)可得出如下規(guī)律：應(yīng)力強度因子KⅠ值及J積分值均與網(wǎng)格密度值成正相關(guān)變化，但隨網(wǎng)格密度的增加逐漸趨于平穩(wěn)；計算時間與網(wǎng)格密度值成負相關(guān)變化，且無收斂趨勢。隨著網(wǎng)格的不斷細化，計算結(jié)果逐漸趨向于準確解，但計算時間卻明顯增長。以尺寸為250 mm×50 mm×50 mm的模型為例，當細化部分奇異單元密度達到96個／cm3時，應(yīng)力強度因子KⅠ值以及J積分值變化都已接近平穩(wěn)，如圖4所示。

圖4 250 mm×50 mm×50 mm尺寸模型計算結(jié)果隨網(wǎng)格數(shù)的變化趨勢

如果將細化區(qū)網(wǎng)格密度再增加1倍，即奇異單元密度為192個／cm3，與前者相比，KⅠ值變化了0.89%，J積分變化了1.76%，而計算成本卻增加了89.5%。

從圖4可以看出，當奇異單元密度較小時，計算值變化較大，難以滿足精度要求，隨著單元網(wǎng)格密度的增大，計算數(shù)據(jù)逐漸穩(wěn)定。以奇異單元網(wǎng)格密度為自變量x，分別以KⅠ值及J積分值為因變量y，對圖4中的各曲線進行擬合。為了使擬合系數(shù)R2更接近1，結(jié)合圖4中的曲線變化趨勢，最終確定采用y＝axb＋c形式的擬合函數(shù)，其中，y為計算值；x為網(wǎng)格密度；a、b、c均為擬合系數(shù)。

3種尺寸（尺寸1，250 mm×50 mm× 50 mm；尺寸2，300 mm×50 mm×50 mm；尺寸3，300 mm×100 mm×100 mm）下的KⅠ擬合方程分別如下：

尺寸1，y＝－19.04x－0.5463＋36；

尺寸2，y＝－18.21x－0.5726＋41.37；

尺寸3，y＝－8.315x－0.9668＋7.515。

3種尺寸下的J積分擬合方程分別如下：

尺寸1，y＝－2.048x－0.5274＋2.169；

尺寸2，y＝－2.342x－0.5687＋2.879；

尺寸3，y＝－0.185 6x－0.9521＋0.089 13。

取x趨向正無窮大時的極限值為穩(wěn)定值，結(jié)果見表2和表3所列。

表2 KⅠ值曲線擬合結(jié)果與理論解對比

表3 J積分值曲線擬合結(jié)果

3組小梁試樣的數(shù)據(jù)輸出點距小梁底面距離均為0.45 mm，由經(jīng)驗公式［18］計算得到在8.4 k N荷載作用下，各試件在數(shù)據(jù)輸出點處的應(yīng)力強度因子KⅠ值的理論值見表2所列，從表2中可以看出，由經(jīng)驗公式得出的理論解與有限元軟件得出的數(shù)值解相一致，由此可以看出利用有限元軟件得到的計算結(jié)果是可靠的。

以數(shù)值解的穩(wěn)定值為準確解，取5%為可接受的最大誤差，通過曲線方程，可以反算出各計算結(jié)果為穩(wěn)定值的95%時所對應(yīng)的奇異單元密度，即為合理網(wǎng)格密度，見表4所列。

由于模型寬度方向每層網(wǎng)格單元尺寸為0.5 cm，奇異單元區(qū)截面積為2 cm2，故每層奇異單元的體積為1 cm3，另外奇異單元自虛擬裂縫向外共有3個圍道，每一層圍道的第1圍道內(nèi)的奇異單元數(shù)為虛擬裂縫周圍的奇異單元數(shù)，故虛擬裂縫周圍的合理單元數(shù)即為合理奇異單元密度的1／3，由此得出不同尺寸模型所對應(yīng)的虛擬裂縫周圍合理單元數(shù)見表4所列。

表4 不同尺寸模型的合理網(wǎng)格密度

3 基于Matlab GUI的數(shù)據(jù)庫開發(fā)

將得到的不同尺寸模型的合理單元尺寸基于Matlab GUI開發(fā)一套交互界面的數(shù)據(jù)查詢系統(tǒng)。系統(tǒng)的主界面如圖5所示。

圖5 數(shù)據(jù)庫主界面

系統(tǒng)由2塊主面板組成，一個為條件選擇面板，另一個為查詢及結(jié)果輸出面板。以模型尺寸以及所要模擬計算的數(shù)據(jù)類型為輸入條件，以奇異單元的合理密度以及虛擬裂縫周圍的合理單元數(shù)為輸出結(jié)果，單擊“OK”按鈕，即可抽取數(shù)據(jù)庫中滿足查詢條件的結(jié)果值。例如在模型尺寸下拉菜單中選擇“300 mm×50 mm×50 mm”，在計算值類型下拉菜單中選擇“J積分值”，單擊“OK”按鈕，得到的查詢結(jié)果為：“奇異單元合理網(wǎng)格密度為135個／cm3”、“虛擬裂縫周圍合理奇異單元個數(shù)為45個”，如圖5b所示。

4 結(jié) 論

（1）對網(wǎng)格密度值和計算結(jié)果值進行函數(shù)擬合，以網(wǎng)格密度趨于無窮時該函數(shù)的極限值為精確值。分析表明，達到該精確值95%時的網(wǎng)格密度可平衡計算結(jié)果精度和計算成本。

（2）對于瀝青混合料小梁彎曲試驗，不同尺寸模型依據(jù)KⅠ值和J積分值確定的合理網(wǎng)格密度也不相同。

（3）通過Matlab GUI對計算結(jié)果建立交互界面查詢系統(tǒng)，能直接給出不同尺寸模型依據(jù)KⅠ值和J積分值確定的合理網(wǎng)格密度。

［1］ 秦 昊，梁乃興，陸兆峰.水－溫作用下瀝青混合料疲勞性能分析［J］.中南大學學報：自然科學版，2011，42（4）：1126－1132.

［2］ Xu Ouming，Han Sen，Gao Shijun，et al.Effect of asphalt content on fatigue limit for asphalt mixture［J］.Journal of Traffic and Transportation Engineering，2009，9（6）：1－6.

［3］ 瞿 鑫，田小革，欒利強，等.瀝青混合料疲勞壽命及其影響因素研究［J］.中外公路，2010，30（5）：266－269.

［4］ 茹忠亮，朱傳銳，張友良，等.斷裂問題的擴展有限元法研究［J］.巖土力學，2011，32（7）：2171－2176.

［5］ 羅 輝，朱宏平，陳傳堯.預(yù)切口瀝青混合料小梁疲勞試驗與數(shù)值研究［J］.土木工程學報，2009，42（6）：126－132.

［6］ 周志剛，張清平，袁秀湘.瀝青混凝土彎曲疲勞試驗疲勞損傷分析［J］.中南大學學報：自然科學版，2011，42（6）：1743－1750.

［7］ Kim Y R，Baek C，Underwood B S，et al.Application of viscoelastic continuum damage model based finite element analysis to predict the fatigue performance of asphalt pavements［J］.KSCE Journal of Civil Engineering，2008，12（2）：109－120.

［8］ 蔣育紅，黃曉明，廖公云.級配碎石夾層路面結(jié)構(gòu)的斷裂力學分析［J］.合肥工業(yè)大學學報：自然科學版，2009，32（4）：511－514.

［9］ 張少雄，王利永，孔 泉.網(wǎng)格粗細對于有限元模態(tài)分析計算的影響［J］.武漢理工大學學報，2006，28（5）：92－94.

［10］ 詹永祥，姚海林，盧正.表面含裂縫瀝青路面低溫收縮斷裂分析［J］.交通科學與工程，2010，26（3）：20－24.

［11］ 羅 輝，朱宏平，陳傳堯，等.基于FE／EFG耦合方法的瀝青路面裂紋擴展模擬［J］.固體力學學報，2009，30（3）：292－300.

［12 李文成，鄭 斌，張青軍，等.瀝青路面多裂紋溫度應(yīng)力的數(shù)值模擬［J］.華中科技大學學報：城市科學版，2010，27（3）：16－20.

［13］ Kuai H，Lee H J，Zi G，et al.Application of generalized J－integral to crack propagation modeling of asphalt concrete under repeated loading［J］.Transportation Research Record，2009，2127：72－81.

［14］ 劉敬輝，王端宜，劉 宇.采用J積分對瀝青混合料抗裂性能進行評價［J］.固體力學學報，2010，31（1）：16－21.

［15］ Gross D，Seelig T.Fracture mechanics with an introduction to micromechanics［M］.New York：Springer Press，2006：59－134.

［16］ 石亦平，周玉蓉.ABAQUS有限元分析實例詳解［M］.北京：機械工業(yè)出版社，2006：9－63.

［17］ Seweryn A.Modeling of singular stress fields using finite element method［J］.International Journal of Solids and Structures，2002，39（18）：4787－4804.

［18］ 王小召.水泥穩(wěn)定碎石基層斷裂性能的試驗研究［J］.鐵道建筑，2008，14（7）：117－120.