徐建平， 朱耀庭

(1. 重慶交通大學(xué) 土木工程學(xué)院，重慶 400074； 2. 江西省高速公路投資集團(tuán)有限責(zé)任公司，江西 南昌 330025；3. 江西省交通科學(xué)研究院，江西 南昌 330038)

0 引 言

瀝青混合料的力學(xué)性能研究在道路工程中處于基礎(chǔ)性地位，對(duì)全面掌握瀝青路面結(jié)構(gòu)行為具有重要作用。瀝青混合料是典型的粘、彈、塑性綜合體，其力學(xué)行為與加載速率、荷載時(shí)間以及溫度等因素密切相關(guān)[1]。

為能對(duì)瀝青路面結(jié)構(gòu)進(jìn)行更精確地應(yīng)力應(yīng)變分析，有必要對(duì)串聯(lián)型粘彈-粘塑性本構(gòu)模型進(jìn)行有限元實(shí)現(xiàn)。目前學(xué)界對(duì)瀝青混合料與時(shí)間的相關(guān)模型進(jìn)行了研究與探討[2-7]，將與時(shí)間率相關(guān)模型植入有限元分析程序中，實(shí)現(xiàn)了瀝青路面結(jié)構(gòu)有限元計(jì)算與分析。然而，對(duì)于瀝青混合料多階段串聯(lián)型粘彈-粘塑性本構(gòu)模型的有限元實(shí)現(xiàn)還少見(jiàn)報(bào)道。

筆者針對(duì)目前瀝青混合料串聯(lián)型粘彈-粘塑性本構(gòu)模型有限元實(shí)現(xiàn)方面的不足，在分析串聯(lián)型模型力學(xué)特性的基礎(chǔ)上，分別對(duì)粘彈、粘塑本構(gòu)模型的增量離散化，采用數(shù)值穩(wěn)定的徑向回退與向后歐拉積分方法[8-9]，并結(jié)合連續(xù)迭代R-M理論[10]，將編寫的材料本構(gòu)模型子程序移植至有限元軟件ABAQUS中，從而實(shí)現(xiàn)了對(duì)瀝青混合料粘彈-粘塑性力學(xué)行為的數(shù)值模擬和分析。筆者通過(guò)建立隱式應(yīng)力積分方程，并推導(dǎo)迭代所需的一致切線剛度矩陣，最后利用算例驗(yàn)證了有限元實(shí)現(xiàn)的有效性。

1 粘彈性本構(gòu)方程有限元

針對(duì)瀝青混合料的粘彈塑性力學(xué)特性，筆者在熱力學(xué)理論基礎(chǔ)上建立了串聯(lián)型粘彈-粘塑性本構(gòu)模型，一維模型如圖1，本構(gòu)方程如式(1)。利用室內(nèi)試驗(yàn)方法和數(shù)學(xué)計(jì)算對(duì)本構(gòu)模型進(jìn)行了參數(shù)獲取和敏感性評(píng)價(jià)[11]。

圖1 串聯(lián)型粘彈-粘塑性本構(gòu)模型Fig. 1 Serial viscoelastic-viscoplastic constitutive model

(1)

式中：pr、lr、h、m、c、γ、q、b、β、s分別為粘彈-粘塑性材料常數(shù)。

根據(jù)Boltzmann線性疊加理論，可通過(guò)一個(gè)遺傳積分方程將各向同性粘彈性材料應(yīng)力應(yīng)變關(guān)系表示如式(2)[12]。

(2)

式中：σij(t)為應(yīng)力張量；Sij(t)為偏應(yīng)力張量；σkk(t)/3為應(yīng)力張量第一不變量；Gijkl(t)為剪切松弛張量；Kijkl(t)為體積松弛張量；ekl為偏應(yīng)變張量；εv為體積應(yīng)變；上標(biāo)符號(hào)表示為對(duì)時(shí)間的一次求導(dǎo)。

假定tr時(shí)刻應(yīng)力狀態(tài)已知，將荷載作用時(shí)間劃為離散時(shí)間間隔為tr+1=tr+Δt，通過(guò)Prony級(jí)數(shù)性質(zhì)表示廣義Maxwell模型力學(xué)特性，以此推導(dǎo)出遺傳積分應(yīng)力增量在Prony級(jí)數(shù)條件下的表現(xiàn)形式如式(3)。

(3)

式(3)包含了彈性及粘性部分偏應(yīng)力。對(duì)粘性部分積分進(jìn)行單獨(dú)分析時(shí)，每個(gè)子模型粘性應(yīng)變化可以通過(guò)式(4)來(lái)表示。

(4)

在時(shí)間增量為Δt條件下，式(4)變化可表示為式(5)：

(5)

由式(4)、(5)可得式(6)：

(6)

假定在時(shí)刻tr≤t≤tr+Δt內(nèi)，粘彈性行為服從等應(yīng)變率為Rε變化，如式(7)：

(7)

式(7)等號(hào)兩邊各減去tr時(shí)刻粘性應(yīng)變，可得到粘性應(yīng)變?cè)隽咳缡?8)：

(8)

則有式(9)：

(9)

同理可推導(dǎo)出體積應(yīng)力增量，其表達(dá)式如式(10)：

(10)

上述遞推關(guān)系式可轉(zhuǎn)化為式(11)：

(11)

定義偏應(yīng)力對(duì)偏應(yīng)變的變化率為剪切切線模量GT，體積應(yīng)力對(duì)體積應(yīng)變的變化率為體積切線模量KT，法向應(yīng)力增量可通過(guò)應(yīng)力的體積增量和偏增量進(jìn)行表示，得出GT和KT表達(dá)式如式(12)：

(12)

在有限元計(jì)算中，還需推導(dǎo)一致切線剛度矩陣，根據(jù)Jacobian矩陣的定義[13]可得出如式(13)：

(13)

2 粘塑性本構(gòu)方程有限元

由于粘塑性本構(gòu)模型方程是用時(shí)間導(dǎo)數(shù)微分方程表示，如何利用精確和有效的數(shù)學(xué)算法對(duì)模型方程進(jìn)行積分運(yùn)算，是有限元實(shí)現(xiàn)過(guò)程中遇到的最大困難。有限元計(jì)算采用等參單元，對(duì)本構(gòu)模型方程積分過(guò)程是通過(guò)高斯點(diǎn)計(jì)算來(lái)實(shí)現(xiàn)，而更新的材料參數(shù)為應(yīng)力和粘塑性兩個(gè)變量。所采用數(shù)值算法必須滿足以下4個(gè)基本條件[14]：① 被積分后的表達(dá)方程與粘塑性本構(gòu)方程相一致；② 數(shù)值方法是穩(wěn)定的；③ 粘塑性增量是一致的；④ 具備二階精確度。因此，隱式積分算法可保證應(yīng)力結(jié)果準(zhǔn)確性，能滿足上述基本條件，筆者將在向后歐拉積分算法基礎(chǔ)上，建立粘塑性本構(gòu)方程的一致切線算子。

2.1 本構(gòu)方程的離散化和返回映射

粘塑性本構(gòu)模型采用一組微分方程和其初值條件進(jìn)行表述，微分方程組由應(yīng)變?cè)隽客ㄟ^(guò)積分離散計(jì)算進(jìn)行實(shí)現(xiàn)，由于向后歐拉方法具有數(shù)值有效性和絕對(duì)收斂性，因而被廣泛應(yīng)用于(粘)塑性有限元數(shù)值計(jì)算中[15-16]。采用后退歐拉積分的方法將粘塑性本構(gòu)方程離散，并對(duì)硬化方程進(jìn)行表述，方程離散后表示如式(14)：

(14)

圖2 徑向回退法幾何闡述Fig. 2 Geometrical view of the radial regression method

若tr時(shí)刻σij(tr)、εij(tr)、εvpij(tr)、Xij(tr)等力學(xué)狀態(tài)變量均為已知，通過(guò)給定Δεij(tr+1)和Δt，則可采用彈性預(yù)估-塑性校正方法求出tr+1時(shí)刻的各個(gè)變量。假設(shè)給定的應(yīng)變?cè)隽喀う舏j(tr+1)全部為彈性應(yīng)變，則引入彈性嘗試應(yīng)力如式(15)：

(15)

式中：Cijkl為彈性模量4階張量。

此時(shí)，屈服準(zhǔn)則由彈性嘗試應(yīng)力改為式(16)：

(16)

(17)

式中：CijklΔεvpkl(tr+1)稱為粘塑性校正項(xiàng)。

由此可見(jiàn)，只需求得εvpkl(tr+1)，則σij(tr+1)等變量將容易獲得。

假定粘塑性材料各向同性，且塑性不可壓縮，則粘塑性校正項(xiàng)偏量可通過(guò)材料剪切模型進(jìn)行表示，如式(18)：

(18)

則由式(18)可得式(19)：

(19)

則得出有效應(yīng)力第2不變量如式(20)：

R[ΔP(tr+1)]

(20)

式(20)可表示為累計(jì)粘塑性應(yīng)變?chǔ)(tr+1)的方程，對(duì)有效應(yīng)力第2不變量進(jìn)行表示如式(21)：

(21)

式(20)和式(21)得到一個(gè)關(guān)于ΔP(tr+1)的非線性方程，如式(22)：

(22)

當(dāng)ΔP(tr+1)求解得出后，則可求解出時(shí)間子步條件下各變量增量，從而得出各變量。值得指出的是，式(22)為一非線性方程組，可采用N-R迭代方法對(duì)未知變量進(jìn)行求解，計(jì)算流程如圖3。

圖3 N-R迭代算法流程Fig. 3 N-R iterative algorithmic process

2.2 一致切線算子

有限元程序在計(jì)算積分過(guò)程中，在時(shí)間增量結(jié)束時(shí)，一致切線剛度矩陣d△σij/d△εij用于表示應(yīng)變?cè)隽康臒o(wú)窮小變化情況下的應(yīng)力分量變化。對(duì)式(14)進(jìn)行微分推導(dǎo)，可得Jacobican行列式如式(23)：

(23)

3 串聯(lián)型粘彈-粘塑性模型有限元

串聯(lián)型粘彈-粘塑性本構(gòu)模型有限元計(jì)算與彈塑性材料相類似，計(jì)算之初需要判斷材料是否已經(jīng)進(jìn)入屈服階段，在未進(jìn)入屈服之前只表現(xiàn)為粘彈性行為，而進(jìn)入屈服階段以后則表現(xiàn)為粘彈性和粘塑性的共同響應(yīng)，應(yīng)變則表現(xiàn)為粘彈性與粘塑性應(yīng)變之和。因此，串聯(lián)型粘彈-粘塑性模型的有限元計(jì)算本質(zhì)上和粘彈性或粘塑性分析本質(zhì)上時(shí)相同的，需根據(jù)屈服法則對(duì)各個(gè)有限元單元做出屈服判斷，然后按不同應(yīng)力應(yīng)變狀態(tài)做出相應(yīng)的判定和計(jì)算處理。對(duì)于已進(jìn)入屈服的有限元單元采用粘塑性計(jì)算方法，對(duì)未進(jìn)入屈服的有限元單元?jiǎng)t采用粘彈性方法進(jìn)行分析。其有限元計(jì)算流程如圖4。

圖4 串聯(lián)型粘彈-粘塑性有限元計(jì)算流程Fig. 4 The calculation process of serial viscoelastic-viscoplastic finite element model

4 數(shù)值驗(yàn)證

筆者利用ABAQUS有限元平臺(tái)進(jìn)行有限元用戶材料(UMAT)子程序的二次開(kāi)發(fā)和編程[18]。通過(guò)對(duì)室內(nèi)實(shí)驗(yàn)進(jìn)行模擬，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對(duì)比，以驗(yàn)證數(shù)值實(shí)現(xiàn)的準(zhǔn)確性。由于文中所涉及本構(gòu)模型為串聯(lián)型，要準(zhǔn)確獲取材料模型參數(shù)要將粘彈性和粘塑性進(jìn)行解耦和分離，對(duì)材料行為時(shí)間相關(guān)性進(jìn)行分配，基于此提出的策略為：通過(guò)小應(yīng)力幅值動(dòng)態(tài)試驗(yàn)獲取粘彈性材料參數(shù)，利用靜態(tài)蠕變?cè)囼?yàn)對(duì)粘彈性和粘塑性解耦獲取粘塑性力學(xué)行為，結(jié)合本構(gòu)模型特征編制數(shù)據(jù)擬合程序，得出粘塑性模型參數(shù)，具體模型參數(shù)和實(shí)驗(yàn)數(shù)據(jù)請(qǐng)見(jiàn)文獻(xiàn)[2]。

瀝青混合料動(dòng)態(tài)模量和靜態(tài)蠕變室內(nèi)試驗(yàn)采用的是圓柱形試件，試件高度為150 mm，直徑為100 mm。試件幾何尺寸、荷載和邊界條件均軸對(duì)稱，可采用軸對(duì)稱單元模型以簡(jiǎn)化為二維模型進(jìn)行計(jì)算，但筆者編制的UMAT子程序仍考慮到了三維狀態(tài)條件。

建模中，有限元模型尺寸和試件相同，剛體模擬加載壓頭和支撐底座。有限元模型網(wǎng)格劃分通過(guò)剖分和掃掠方法得到圓心向圓周發(fā)散的三維網(wǎng)格，單元類型為20結(jié)點(diǎn)減縮積分單元C3D20R。在實(shí)際的試驗(yàn)過(guò)程中，對(duì)試件底面并沒(méi)有進(jìn)行約束，同時(shí)由于實(shí)驗(yàn)過(guò)程采用雙層乳膠膜以保障試件端面的自由變形，只進(jìn)行了Z方向約束，即U3=0。試件頂面為在施加均布應(yīng)力的加載方式。

4.1 粘彈性結(jié)果

粘彈性材料響應(yīng)的表征可體現(xiàn)為蠕變和應(yīng)力松弛，由于瀝青混合料材料的特殊性，一般室內(nèi)對(duì)粘彈性采用蠕變實(shí)驗(yàn)獲取蠕變應(yīng)變或柔度，再通過(guò)蠕變和松弛之間的粘彈性關(guān)系，可計(jì)算出松弛應(yīng)力或模量[19]。計(jì)算結(jié)果如圖5。

圖5 應(yīng)力松弛解析解和蠕變數(shù)值解與UMAT計(jì)算對(duì)比Fig. 5 Comparison between analytical solution of stress relaxation and numerical solution of creep and UMAT calculation

由松弛和蠕變計(jì)算結(jié)果可看出：粘彈性UMAT應(yīng)力松弛計(jì)算解和解析解基本上相同，結(jié)果幾乎無(wú)差別，而蠕變有限元解和數(shù)值解也基本上相同，表明有限元子程序有非常好的準(zhǔn)確性。

4.2 粘塑性結(jié)果

粘塑性本構(gòu)方程由一組隱形方程組構(gòu)成，無(wú)法直接得出解析解，通過(guò)文獻(xiàn)[2]中改進(jìn)的有限差分格式可獲取本構(gòu)方程的數(shù)值曲線，本節(jié)采用蠕變模擬進(jìn)行差分解和有限元解的對(duì)比分析。同時(shí)，由于粘塑性需要采用達(dá)到屈服條件才可驗(yàn)證，模擬加載應(yīng)力需大于材料的粘塑性屈服應(yīng)力。不同條件蠕變計(jì)算結(jié)果如圖6。

圖6 有限元計(jì)算和改進(jìn)差分法計(jì)算結(jié)果對(duì)比Fig. 6 Comparison between the calculation results of finite element calculation and the improved difference method

由蠕變實(shí)驗(yàn)?zāi)M和計(jì)算結(jié)果對(duì)比可知，粘塑性子程序是準(zhǔn)確和可靠的，能正確地反映出粘塑性材料的力學(xué)響應(yīng)和內(nèi)變量的發(fā)展趨勢(shì)。

4.3 串聯(lián)型粘彈-粘塑性結(jié)果

串聯(lián)型粘彈-粘塑性材料的有限元計(jì)算之初要判定材料是否已進(jìn)入屈服階段，在未進(jìn)入屈服以前只表現(xiàn)為粘彈性行為，進(jìn)入屈服以后則表現(xiàn)為粘彈性及粘塑性共同響應(yīng)。因此，串聯(lián)型粘彈-粘塑性計(jì)算結(jié)果驗(yàn)證分為兩個(gè)階段，一個(gè)為材料未達(dá)到屈服狀態(tài)；另一個(gè)為材料進(jìn)入屈服狀態(tài)而發(fā)生了粘塑性變形。未進(jìn)入屈服和進(jìn)入屈服狀態(tài)的模擬結(jié)果如圖7、8。

圖7 未達(dá)屈服條件下計(jì)算結(jié)果對(duì)比Fig. 7 Comparison of calculation results under non-yield conditions

圖8 屈服條件下計(jì)算結(jié)果對(duì)比Fig. 8 Comparison of calculation results under yield condition

由圖7、8可見(jiàn)：無(wú)論屈服狀態(tài)還是未達(dá)屈服狀態(tài)，計(jì)算結(jié)果和理論分析都具有很好一致性，串聯(lián)型粘彈-粘塑性本構(gòu)模型有限元實(shí)現(xiàn)是有效的，能準(zhǔn)確反映出瀝青混合料與時(shí)間相關(guān)力學(xué)行為發(fā)展和演化規(guī)律。

5 結(jié) 論

筆者以瀝青混合料與時(shí)間相關(guān)本構(gòu)模型為研究對(duì)象，通過(guò)編制用戶子程序，對(duì)串聯(lián)型粘彈-粘塑性模型進(jìn)行了有限元分析，通過(guò)瀝青混合料試件在典型室內(nèi)試驗(yàn)條件力學(xué)行為進(jìn)行了有限元模擬和室內(nèi)試驗(yàn)數(shù)據(jù)的對(duì)比，驗(yàn)證了筆者提出方法的正確性和可靠性，并得出以下結(jié)論：

1)串聯(lián)型本構(gòu)模型的有限元可依據(jù)模型特點(diǎn)表示為粘彈性和粘塑性兩部分的應(yīng)力相同而應(yīng)變相加；

2)通過(guò)遞推關(guān)系得出粘彈性、粘塑性本構(gòu)模型的離散化方程，數(shù)值解和解析解進(jìn)行對(duì)比，其結(jié)果基本一致，驗(yàn)證了文中模型離散化法有效性；

3)粘塑性模型采用徑向回退與向后歐拉積分方法，并結(jié)合連續(xù)迭代R-M法可有效地解決模型離散化非線性方程組計(jì)算過(guò)程的穩(wěn)定性和收斂問(wèn)題；

4)當(dāng)將串聯(lián)型模型用于瀝青路面結(jié)構(gòu)計(jì)算中，只需將瀝青結(jié)構(gòu)層材料定義為模型屬性即可快速實(shí)現(xiàn)路面結(jié)構(gòu)永久變形計(jì)算和分析。