格賓網(wǎng)片拉伸試驗的有限元法數(shù)值實現(xiàn)

李立云, 吳浩瀚, 李雷

(北京工業(yè)大學(xué)城市建設(shè)學(xué)部, 北京 100124)

格賓網(wǎng)片是鍍鋅低碳鋼絲織成的金屬網(wǎng),廣泛應(yīng)用于土木水利工程領(lǐng)域內(nèi),如防護(hù)巖崩落石[1]、邊坡防護(hù)[2]、路堤穩(wěn)定性增強(qiáng)[3]和堤岸防護(hù)[4]中。目前,格賓網(wǎng)片的力學(xué)性能研究多采用物理試驗展開,以獲得材料特性或校準(zhǔn)數(shù)值模型[5-6]。基于物理試驗,人們明晰了格賓網(wǎng)片宏觀力學(xué)性能和變形破壞特征,發(fā)現(xiàn)影響網(wǎng)片拉伸強(qiáng)度的關(guān)鍵因素有網(wǎng)孔尺寸、網(wǎng)絲直徑和絞邊方式[7-8]。然而,通過物理試驗所獲的力學(xué)指標(biāo)有限,若進(jìn)行更細(xì)致研究則需布設(shè)大量傳感器,導(dǎo)致物理試驗實施存在諸多困難和不確定性。數(shù)值模擬方法可以消除物理試驗中的施作困難和不確定性,為研究網(wǎng)片各組分對其整體力學(xué)性能的影響帶來便利性,進(jìn)而優(yōu)化網(wǎng)片設(shè)計或改進(jìn)物理試驗方案[9]。

網(wǎng)片力學(xué)性能數(shù)值研究中最常采用的為離散元法[10-12]。但是,格賓網(wǎng)片的離散元模型過度依賴物理試驗進(jìn)行校準(zhǔn)且不能反映網(wǎng)片的三維變形特征。 眾所周知,有限元法適用于具有非線性幾何形狀、復(fù)雜力學(xué)行為和各種接觸條件的連續(xù)介質(zhì)問題,為降低計算成本,一些使用桁架單元或膜單元的簡化網(wǎng)片有限元模型被提出[13-14]。但是,簡化網(wǎng)片模型不能準(zhǔn)確描述荷載作用下網(wǎng)絲的力學(xué)響應(yīng),實際風(fēng)險點(diǎn)位可能被忽略。現(xiàn)存少數(shù)可用的三維數(shù)值模型雖已較為成熟,但都十分復(fù)雜。因此,構(gòu)建一個正確可信、方便快捷的三維格賓網(wǎng)片有限元模型,是對網(wǎng)片開展更深入力學(xué)性能研究中亟待解決的問題,還具有推廣格賓結(jié)構(gòu)實際工程應(yīng)用的價值。

鑒于此,聚焦于格賓網(wǎng)片拉伸試驗的數(shù)值實現(xiàn),探討格賓網(wǎng)片精細(xì)化三維模型[15]在有限元數(shù)值模擬中的適用性和擇取關(guān)鍵參數(shù)的可靠性,首先基于三維精細(xì)化建模方法建立格賓網(wǎng)片數(shù)值模型,并將其移植于有限元軟件Abaqus中模擬恒定速率下的準(zhǔn)靜態(tài)拉伸試驗。之后,分析網(wǎng)片系統(tǒng)的能量曲線,探討保證模擬結(jié)果可靠的情況下控制計算成本的最優(yōu)方案。最后,通過對比分析數(shù)值模擬與物理試驗結(jié)果,驗證模擬方法的合理性及模擬結(jié)果的正確性,揭示拉伸荷載作用下格賓網(wǎng)片的變形特征以及應(yīng)力分布情況。為進(jìn)行不同規(guī)格網(wǎng)片的拉伸試驗創(chuàng)造有利條件,為后續(xù)格賓結(jié)構(gòu)力學(xué)性能的深入研究奠定基礎(chǔ)。

1 網(wǎng)片特征及夾持方式

格賓網(wǎng)片的幾何特征按空間結(jié)構(gòu)分為六邊形網(wǎng)孔、邊絲和端絲,如圖1(a)所示。幾何尺寸為長幅L、寬幅W、網(wǎng)絲直徑d和網(wǎng)孔尺寸M×N,網(wǎng)孔尺寸由雙絞合網(wǎng)絲軸線之間的距離M與上下兩結(jié)點(diǎn)之間的距離N控制,如圖1(b)所示。

s為雙絞絲長度

以尺寸L×W=1 000 mm×1 000 mm的網(wǎng)片物理拉伸試驗[15]為基準(zhǔn)。物理試驗中網(wǎng)片被夾具夾持的尺寸范圍為930 mm×660 mm,夾持方式如圖2所示,根據(jù)《工程用機(jī)編鋼絲網(wǎng)及組合體》(YB/T 4221—2016)夾持[16]。

黑色圓點(diǎn)表示夾具和網(wǎng)片之間的連接點(diǎn)

2 格賓網(wǎng)片數(shù)值模型

2.1 格賓網(wǎng)片三維有限元建模

格賓網(wǎng)片的數(shù)值模型采用筆者提出的三維有限元精細(xì)化建模方法[15]實現(xiàn)。該方法基于SolidWorks軟件和ABAQUS軟件的結(jié)合,充分考慮了網(wǎng)片的空間復(fù)雜性,通過捕捉網(wǎng)絲的運(yùn)動軌跡實現(xiàn)復(fù)雜網(wǎng)片建模。網(wǎng)片建模過程如圖3所示。

圖3 格賓網(wǎng)片建模流程圖Fig.3 The modelling procedure of gabion mesh

格賓網(wǎng)片模型通過復(fù)制裝配相互交織纏繞且做圓柱螺旋運(yùn)動的單根網(wǎng)絲得到,構(gòu)建單根網(wǎng)絲模型所需的參數(shù)為:雙絞絲長度s=30 mm、網(wǎng)絲直徑d=2.5 mm、網(wǎng)孔尺寸M×N=80 mm×100 mm以及雙絞絲纏繞數(shù)n=3次。根據(jù)上述參數(shù),利用SolidWorks軟件建立單根網(wǎng)絲模型,并在Abaqus中完成裝配,得到格賓網(wǎng)片數(shù)值模型如圖4所示。

網(wǎng)絲間接觸定義為通用接觸,切向摩擦系數(shù)0.3,法向硬接觸以嚴(yán)格控制模型穿透。為模擬物理試驗中夾具作用于網(wǎng)片的拉伸荷載,對數(shù)值模型邊界條件做如下處理:首先,設(shè)定圖4中標(biāo)識為RP的參考點(diǎn),耦合上部端絲到參考點(diǎn)RP-1,并使該點(diǎn)1 000 s內(nèi)沿y軸正方向拉伸170 mm,加載方式采用平滑分析步以減少初始荷載施加時對網(wǎng)片造成的影響;耦合下部端絲到參考點(diǎn)RP-2,該點(diǎn)完全固定。其次,約束黑色圓圈內(nèi)雙絞絲x方向的位移。

2.2 網(wǎng)絲材料參數(shù)標(biāo)定

低碳鋼絲是典型具有長塑性階段的材料,為在數(shù)值模擬中準(zhǔn)確地描述其變形過程,定義材料屬性時要對在物理試驗中得到的本構(gòu)模型數(shù)據(jù)進(jìn)行標(biāo)定。參考朱凱[17]的2.7 mm網(wǎng)絲物理拉伸試驗結(jié)果,如圖5所示。

圖5 網(wǎng)絲的應(yīng)力-應(yīng)變曲線Fig.5 Stress-strain curves of the wire

數(shù)值模擬中的材料屬性參數(shù)包括網(wǎng)絲彈性模量E、真實應(yīng)力σt和塑性應(yīng)變εp。根據(jù)圖5初始直線段,可以計算得出彈性模量E=20 GPa。但是物理試驗得到的應(yīng)力-應(yīng)變數(shù)據(jù)為名義應(yīng)力σn和名義應(yīng)變εn,由式(1)、式(2)獲得。數(shù)值計算中所需的真實應(yīng)力σt和塑性應(yīng)變εp需要通過名義應(yīng)變εn和名義應(yīng)力σn轉(zhuǎn)換獲得:由式(3)、式(4)計算得出真實應(yīng)變εt和真實應(yīng)力σt,并由式(5)得到塑性應(yīng)變εp。

(1)

(2)

式中:Δl為網(wǎng)片鋼絲長度的變化量;l0為網(wǎng)片鋼絲的初始長度;F為荷載;A0為網(wǎng)片鋼絲的初始截面面積。

(3)

(4)

式中:l為網(wǎng)絲的當(dāng)前長度;A為網(wǎng)絲的當(dāng)前截面面積。

(5)

式(5)中:εe為彈性應(yīng)變。

根據(jù)式(1)～式(5)對網(wǎng)絲拉伸物理試驗數(shù)據(jù)進(jìn)行處理后得到材料塑性參數(shù),如表1所示。采用表1中的真實應(yīng)力σt和塑性應(yīng)變εp賦予材料屬性,進(jìn)行單根網(wǎng)絲的拉伸數(shù)值模擬后得到應(yīng)力-應(yīng)變曲線,如圖5所示。圖5中數(shù)值模擬與物理試驗結(jié)果吻合程度非常高,說明材料的彈塑性參數(shù)設(shè)置正確。

表1 網(wǎng)絲材料塑性參數(shù)Table 1 Plasticity parameters of wire material

3 結(jié)果與討論

3.1 數(shù)值試驗可靠性分析

所依托的網(wǎng)片拉伸試驗是一個典型低加載速率下的準(zhǔn)靜態(tài)問題。準(zhǔn)靜態(tài)求解分析一般用ABAQUS/Standard算法解決,但Standard算法解決非線性問題時,往往需要很長的計算時間,甚至無法收斂到最終解。因此,采用廣泛用于瞬態(tài)響應(yīng)分析的ABAQUS/Explicit算法,該算法很好地改善了非線性問題的收斂性。同時,研究中采用三維8節(jié)點(diǎn)減縮積分六面體實體單元(C3D8R)和質(zhì)量縮放功能,前者可以克服剪切自鎖與模型扭曲問題,后者降低了Explicit算法的計算成本。

減縮積分單元在各個方向上的積分點(diǎn)數(shù)目較少,能節(jié)約計算成本。但較少的積分點(diǎn)數(shù)目會導(dǎo)致沙漏現(xiàn)象發(fā)生,采用偽應(yīng)變能Ea進(jìn)行沙漏控制。偽應(yīng)變能、塑性變形耗散能、可恢復(fù)的應(yīng)變能和單元扭曲控制耗散能等能量的總和為應(yīng)變能Ei,偽應(yīng)變能占應(yīng)變能的多少是判斷計算結(jié)果是否可靠的重要指標(biāo)。偽應(yīng)變能越大,計算結(jié)果可靠度越低。為了便于討論,定義偽應(yīng)變能與總應(yīng)變能的比率為μ,如式(6)所示。根據(jù)文獻(xiàn)[18-19],μ需控制在10%以下且越小越好。

(6)

偽應(yīng)變能的大小與模型單元數(shù)量密切相關(guān),單位體積內(nèi)單元數(shù)量越多,偽應(yīng)變能越小,計算結(jié)果越可靠,但計算成本也會急劇上升。表2列出了5種計算工況下的單元數(shù)量和計算耗時情況,可以看出,隨著單位體積內(nèi)單元數(shù)量的增加,單位耗時基本呈線性增加。

表2 數(shù)值試驗工況Table 2 Numerical test conditions

網(wǎng)片拉伸過程中的偽應(yīng)變能與應(yīng)變能歷程曲線如圖6所示。偽應(yīng)變能與應(yīng)變能曲線隨著拉伸率的增大呈現(xiàn)出相似的上升趨勢。圖6(a)表明,偽應(yīng)變能在拉伸率到達(dá)12%前增長緩慢,之后迅速增大,至試驗結(jié)束時依然保持較大的增長速率。造成這種現(xiàn)象的原因是:試驗開始階段的單元變形較小,沙漏不明顯;但隨著網(wǎng)片被逐漸拉長,單元扭曲變形增大,導(dǎo)致偽應(yīng)變能急劇增加。此外,工況5～工況1隨單元細(xì)分偽應(yīng)變能逐漸減小,且增長幅度亦逐漸降低,工況5與工況1的偽應(yīng)變能在試驗結(jié)束時最大相差了115 J。同時,由圖6(b)可知,各工況下應(yīng)變能變化趨勢大致相同。應(yīng)變能在拉伸率到達(dá)10%前增長較慢,之后快速增長。但在試驗結(jié)束前,由于網(wǎng)絲斷裂后應(yīng)力重分布,應(yīng)變能呈現(xiàn)出減小再增加的波動狀態(tài)。

能量比率μ變化如圖7所示。試驗伊始,荷載的施加使μ曲線呈現(xiàn)出激增又下跌的起伏現(xiàn)象;之后,μ稍微增加并保持穩(wěn)定,直至拉伸率為11%時曲線發(fā)生波動;拉伸率達(dá)到17%時網(wǎng)絲失效,曲線增幅變大。雖然各工況下曲線發(fā)展趨勢相近,但能量比率的平均值差別較大。隨著單元的細(xì)分,能量比率的平均值從工況5的9%降低至工況1的2.5%,整體呈遞減趨勢;雖然工況3與工況4的能量比率曲線基本重合,但這屬于特殊情況。分析整體試驗情況,細(xì)分單元是減小偽應(yīng)變能的一種可靠手段。此外,工況2與工況1平均值相差僅為0.5%,但隨著單元總數(shù)的增加,計算成本會急劇上升,故在保持較低計算成本的情況下工況2的單元劃分方式有效控制了沙漏,較為合理。

圖7 偽應(yīng)變能與應(yīng)變能比率Fig.7 Ratio of artificial to total strain energy

對質(zhì)量縮放功能的必要性和可靠性進(jìn)行討論。

在Explicit算法中,算法穩(wěn)定所需的時間增量可表示為

(7)

式(7)中:Le和ρ分別為元素最小特征尺寸和密度;E為元素的彈性模量。

網(wǎng)片拉伸過程中,荷載引起部分單元極度變形,最小特征尺寸快速減小,從而消耗極多的運(yùn)算資源。通過質(zhì)量縮放增大單元質(zhì)量,可以提高模型的收斂性與計算效率。但質(zhì)量增加導(dǎo)致系統(tǒng)動能增加,進(jìn)而可能會改變網(wǎng)片的受力狀態(tài),影響模型計算結(jié)果的可靠性。因此,確保試驗的準(zhǔn)靜態(tài)條件尤為重要。采用半自動質(zhì)量縮放將最小增量步放大為0.003 s。圖8展示了不同工況下網(wǎng)片系統(tǒng)的動能以及能量平衡的歷程曲線。

圖8 動能與能量平衡歷程曲線Fig.8 Kinetic energy and energy balance history

由圖8(a)可知,網(wǎng)片系統(tǒng)動能在荷載施加后至10%拉伸率前先平穩(wěn)增加,然后至拉伸率為17%時波動下降,該種動能歷程曲線特征由平滑分析步的加載方式導(dǎo)致。當(dāng)拉伸率達(dá)到17%時,網(wǎng)絲斷裂失效,動能曲線呈現(xiàn)出急劇漲跌姿態(tài)。對比圖6(b)和圖8(a)可知,數(shù)值試驗中動能與應(yīng)變能相比幾乎可以忽略不計,表明本研究中采用的質(zhì)量縮放沒有改變網(wǎng)片系統(tǒng)的準(zhǔn)靜態(tài)條件。

能量平衡Et為網(wǎng)片系統(tǒng)產(chǎn)生與耗散的能量總合,即Et=Ek+Ei-外力做的功-因質(zhì)量縮放所產(chǎn)生的能量(Ek為動能),其歷程曲線如圖8(b)所示。由其定義可知,能量平衡越接近于0,則模擬過程的收斂性越好,數(shù)值模擬結(jié)果的可靠性越強(qiáng)。數(shù)值模擬中的能量平衡曲線在拉伸率為13%前幾乎為0;當(dāng)拉伸率超過13%以后,能量平衡快速增大,但試驗結(jié)束時能量平衡最大數(shù)值為工況5中56 J,最小數(shù)值僅為工況1中13 J,所有工況下能量平衡與應(yīng)變能的比值都可以忽略不計。因此,數(shù)值試驗中的能量收斂情況良好。

綜上所述,數(shù)值試驗獲得的結(jié)果是可靠的。

3.2 網(wǎng)片的軸力響應(yīng)分析

數(shù)值試驗得到了單軸拉伸作用下格賓網(wǎng)片的力學(xué)響應(yīng)特征,將其與在物理試驗[15]中得到的結(jié)果進(jìn)行對比分析。圖9給出了網(wǎng)片的軸力-伸長率曲線,可以看出,所有工況下軸力-伸長率曲線發(fā)展趨勢基本相同,都能較好地重現(xiàn)格賓網(wǎng)片拉伸破壞過程的兩個階段,即六邊形網(wǎng)孔機(jī)械變形階段和網(wǎng)絲彈塑性變形階段。第一階段的伸長率小于13%,軸力與伸長率近似為線性關(guān)系且曲線斜率較平緩,網(wǎng)片的伸長量主要為六邊形網(wǎng)孔機(jī)械變形與網(wǎng)絲彈性變形。第二階段的伸長率為13%～17%,網(wǎng)孔機(jī)械變形停止,軸力-伸長率曲線仍近似為線性關(guān)系但斜率較第一階段增大良多,網(wǎng)片伸長量以網(wǎng)絲塑性變形為主。網(wǎng)絲達(dá)到極限承載力后失效斷裂,表現(xiàn)為脆性破壞和軸力呈快速回落。

圖9 網(wǎng)片拉伸試驗軸力曲線Fig.9 Axial force curves of mesh tensile test

從峰值軸力(圖9)上來看,不同工況有明顯的差異。工況5軸力峰值與物理試驗相差6 kN,相對誤差達(dá)到了23%,且曲線波動情況明顯呈鋸齒狀,尤其在網(wǎng)絲斷裂前這一現(xiàn)象最為嚴(yán)重。工況5～工況1曲線發(fā)展趨勢逐漸趨于平滑,軸力峰值不斷增加,工況1軸力峰值與物理試驗僅相差2 kN,相對誤差為7%。數(shù)值試驗結(jié)果得到優(yōu)化的原因是單元劃分逐漸精細(xì),網(wǎng)片系統(tǒng)的偽應(yīng)變能控制得到改善。數(shù)值結(jié)果表明,當(dāng)采取減縮積分單元進(jìn)行模擬時,需要嚴(yán)格控制單元大小以控制偽應(yīng)變能。一旦偽應(yīng)變能過大,即能量比率μ>5%時,會使模擬結(jié)果曲線波動異常和軸力誤差增大。當(dāng)能量比率μ<5%時,模擬結(jié)果曲線平滑且與物理試驗結(jié)果較為符合。3.1節(jié)結(jié)果已表明,單元劃分過細(xì)會導(dǎo)致計算成本急劇上升,同時模擬結(jié)果并無較大改善,因此,此類數(shù)值試驗的能量比率μ可控在約5%。

由圖9可知,數(shù)值試驗與物理試驗的最終網(wǎng)片拉伸率相差2%,該誤差由兩點(diǎn)導(dǎo)致。首先,邊界條件的不同是造成這一誤差的主要原因。施加拉伸荷載時,數(shù)值試驗耦合端絲結(jié)點(diǎn)進(jìn)行加載,而物理試驗采用夾具夾持網(wǎng)片加載。夾具拉伸網(wǎng)片時會造成與其接觸的雙絞絲壓縮變形,該段變形并未被記錄處理,反而直接將夾具的位移距離作為網(wǎng)片被拉伸的絕對距離,導(dǎo)致計算獲得的網(wǎng)片拉伸率大于其實際拉伸率。數(shù)值模擬試驗通過耦合端絲加載,使網(wǎng)片整體受拉,很好地發(fā)現(xiàn)了這一問題,為后續(xù)物理試驗的改進(jìn)提供了經(jīng)驗:物理試驗應(yīng)記錄網(wǎng)片變形前后的絕對位置以減少誤差。其次,采用Explicit算法計算時會不可避免地累積誤差,但在節(jié)約計算成本的條件下,這種累積誤差可以接受。相較于物理試驗與離散元數(shù)值模擬,有限元數(shù)值模擬更為經(jīng)濟(jì)快捷,且很好地揭示了網(wǎng)片的力學(xué)性能。

3.3 網(wǎng)片變形特征與應(yīng)力分布

鑒于數(shù)值試驗中各工況格賓網(wǎng)片的力學(xué)響應(yīng)相似,故選擇工況1的結(jié)果進(jìn)行分析。圖10顯示了工況1中網(wǎng)絲失效前一時刻的狀態(tài)。網(wǎng)片在拉伸過程中X軸方向上出現(xiàn)了類似“頸縮”的現(xiàn)象,表現(xiàn)為六邊形網(wǎng)孔在寬幅方向上縮短的同時在長度方向上伸長。同時,寬幅邊緣處的六邊形網(wǎng)孔產(chǎn)生極大形變,網(wǎng)絲幾乎緊貼在一起。物理試驗[15]中網(wǎng)片的拉伸破壞情況如圖11所示。對比圖10與圖11發(fā)現(xiàn)格賓網(wǎng)片的宏觀變形特征與物理試驗的宏觀變形特征基本吻合。由圖10可知邊界約束處的鋼絲應(yīng)力集中現(xiàn)象最為嚴(yán)重,最大應(yīng)力約為644 MPa;此外,最靠近端絲的六邊形網(wǎng)孔開口處有較為嚴(yán)重的應(yīng)力集中,約為600 MPa;頁面網(wǎng)絲的應(yīng)力平均為300 MPa。

圖10 網(wǎng)片拉伸數(shù)值試驗應(yīng)力云圖Fig.10 Numerical tensile test stress diagram

圖11 網(wǎng)片拉伸物理試驗的變形破壞Fig.11 Tensile deformation failure of gabion mesh

數(shù)值試驗中網(wǎng)片的應(yīng)力分布情況與網(wǎng)絲斷裂位置相對應(yīng)。圖12給出了網(wǎng)絲斷裂時刻的變形圖以及各模擬工況下網(wǎng)絲的斷裂情況。網(wǎng)絲斷裂位置都位于應(yīng)力集中現(xiàn)象最嚴(yán)重的邊界約束處。其中,工況4和5僅有一個斷裂點(diǎn)位于網(wǎng)片右側(cè),工況1和2則同時有兩個斷裂點(diǎn)位于網(wǎng)片左側(cè),在工況3中有4個點(diǎn)處的網(wǎng)絲同時斷裂。這種網(wǎng)絲失效模式也體現(xiàn)在物理試驗中,但不同的是物理試驗中與夾具相連的六邊形網(wǎng)孔開口處也時常斷裂。數(shù)值試驗和物理試驗結(jié)果存在上述差異的主要原因歸納如下:一方面,數(shù)值模擬的各項條件過于理想,故斷裂位置集中在應(yīng)力集中最嚴(yán)重處;另一方面,物理試驗中網(wǎng)絲材料性能、網(wǎng)片制作過程等存在的差異使得有較嚴(yán)重應(yīng)力集中的六邊形開口處也較易斷裂。

圖12 網(wǎng)絲斷裂示意圖Fig.12 Schematic diagram of wire failure

4 結(jié)論

通過建立格賓網(wǎng)片三維數(shù)值模型并將其移植于ABAQUS有限元軟件,數(shù)值實現(xiàn)了格賓網(wǎng)片拉伸試驗,得到如下結(jié)論

(1)格賓網(wǎng)片力學(xué)性能數(shù)值研究中的材料參數(shù)應(yīng)該采用彈性模量、真實應(yīng)力和塑性應(yīng)變。

(2)采用減縮積分單元與質(zhì)量縮放功能節(jié)約了計算成本。為保證數(shù)值試驗結(jié)果的可靠性,需要控制沙漏情況,沙漏嚴(yán)重會使偽應(yīng)變能過大,導(dǎo)致結(jié)果不夠準(zhǔn)確。細(xì)化單元能夠有效控制沙漏情況。

(3)數(shù)值試驗結(jié)果與物理試驗吻合良好,軸力-拉伸率曲線的兩階段特點(diǎn)明顯,網(wǎng)片的應(yīng)力分布很好地反映了物理試驗的宏觀變形特征與破壞特點(diǎn)。

(4)造成數(shù)值試驗與物理試驗差異的可能原因是約束條件不同以及Explicit算法的誤差累積。