杜軍虎， 韓 丁

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

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

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

1 小梁彎曲試驗(yàn)的有限元模擬

1.1 模擬工況

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

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

為了分析網(wǎng)格劃分密度對(duì)仿真結(jié)果的影響，固定加載值大小為8.4 k N，底部?jī)芍挝恢镁嘈×簜?cè)面的距離為整個(gè)小梁長(zhǎng)度的3／16，通過約束支撐位置x、y方向位移自由度來模擬對(duì)小梁的支撐。加載位置為支撐間三分點(diǎn)處。

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

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

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

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

1.3 計(jì)算結(jié)果

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

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

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

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

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

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

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

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

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

尺寸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趨向正無窮大時(shí)的極限值為穩(wěn)定值，結(jié)果見表2和表3所列。

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

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

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

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

由于模型寬度方向每層網(wǎng)格單元尺寸為0.5 cm，奇異單元區(qū)截面積為2 cm2，故每層奇異單元的體積為1 cm3，另外奇異單元自虛擬裂縫向外共有3個(gè)圍道，每一層圍道的第1圍道內(nèi)的奇異單元數(shù)為虛擬裂縫周圍的奇異單元數(shù)，故虛擬裂縫周圍的合理單元數(shù)即為合理奇異單元密度的1／3，由此得出不同尺寸模型所對(duì)應(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塊主面板組成，一個(gè)為條件選擇面板，另一個(gè)為查詢及結(jié)果輸出面板。以模型尺寸以及所要模擬計(jì)算的數(shù)據(jù)類型為輸入條件，以奇異單元的合理密度以及虛擬裂縫周圍的合理單元數(shù)為輸出結(jié)果，單擊“OK”按鈕，即可抽取數(shù)據(jù)庫中滿足查詢條件的結(jié)果值。例如在模型尺寸下拉菜單中選擇“300 mm×50 mm×50 mm”，在計(jì)算值類型下拉菜單中選擇“J積分值”，單擊“OK”按鈕，得到的查詢結(jié)果為：“奇異單元合理網(wǎng)格密度為135個(gè)／cm3”、“虛擬裂縫周圍合理奇異單元個(gè)數(shù)為45個(gè)”，如圖5b所示。

4 結(jié) 論

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

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

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

［1］ 秦 昊，梁乃興，陸兆峰.水－溫作用下瀝青混合料疲勞性能分析［J］.中南大學(xué)學(xué)報(bào)：自然科學(xué)版，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］ 瞿 鑫，田小革，欒利強(qiáng)，等.瀝青混合料疲勞壽命及其影響因素研究［J］.中外公路，2010，30（5）：266－269.

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

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

［6］ 周志剛，張清平，袁秀湘.瀝青混凝土彎曲疲勞試驗(yàn)疲勞損傷分析［J］.中南大學(xué)學(xué)報(bào)：自然科學(xué)版，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］ 蔣育紅，黃曉明，廖公云.級(jí)配碎石夾層路面結(jié)構(gòu)的斷裂力學(xué)分析［J］.合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2009，32（4）：511－514.

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

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

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

［12 李文成，鄭 斌，張青軍，等.瀝青路面多裂紋溫度應(yīng)力的數(shù)值模擬［J］.華中科技大學(xué)學(xué)報(bào)：城市科學(xué)版，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積分對(duì)瀝青混合料抗裂性能進(jìn)行評(píng)價(jià)［J］.固體力學(xué)學(xué)報(bào)，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有限元分析實(shí)例詳解［M］.北京：機(jī)械工業(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)定碎石基層斷裂性能的試驗(yàn)研究［J］.鐵道建筑，2008，14（7）：117－120.