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

李立云, 吳浩瀚, 李雷

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

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

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

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

s為雙絞絲長度

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

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

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

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

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

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

格賓網(wǎng)片模型通過復制裝配相互交織纏繞且做圓柱螺旋運動的單根網(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,法向硬接觸以嚴格控制模型穿透。為模擬物理試驗中夾具作用于網(wǎng)片的拉伸荷載,對數(shù)值模型邊界條件做如下處理:首先,設(shè)定圖4中標識為RP的參考點,耦合上部端絲到參考點RP-1,并使該點1 000 s內(nèi)沿y軸正方向拉伸170 mm,加載方式采用平滑分析步以減少初始荷載施加時對網(wǎng)片造成的影響;耦合下部端絲到參考點RP-2,該點完全固定。其次,約束黑色圓圈內(nèi)雙絞絲x方向的位移。

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

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

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

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

(1)

(2)

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

(3)

(4)

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

(5)

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

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

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

3 結(jié)果與討論

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

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

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

(6)

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

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

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

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

圖7 偽應變能與應變能比率Fig.7 Ratio of artificial to total strain energy

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

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

(7)

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

網(wǎng)片拉伸過程中,荷載引起部分單元極度變形,最小特征尺寸快速減小,從而消耗極多的運算資源。通過質(zhì)量縮放增大單元質(zhì)量,可以提高模型的收斂性與計算效率。但質(zhì)量增加導致系統(tǒng)動能增加,進而可能會改變網(wǎng)片的受力狀態(tài),影響模型計算結(jié)果的可靠性。因此,確保試驗的準靜態(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%時波動下降,該種動能歷程曲線特征由平滑分析步的加載方式導致。當拉伸率達到17%時,網(wǎng)絲斷裂失效,動能曲線呈現(xiàn)出急劇漲跌姿態(tài)。對比圖6(b)和圖8(a)可知,數(shù)值試驗中動能與應變能相比幾乎可以忽略不計,表明本研究中采用的質(zhì)量縮放沒有改變網(wǎng)片系統(tǒng)的準靜態(tài)條件。

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

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

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

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

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

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

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

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

鑒于數(shù)值試驗中各工況格賓網(wǎng)片的力學響應相似,故選擇工況1的結(jié)果進行分析。圖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可知邊界約束處的鋼絲應力集中現(xiàn)象最為嚴重,最大應力約為644 MPa;此外,最靠近端絲的六邊形網(wǎng)孔開口處有較為嚴重的應力集中,約為600 MPa;頁面網(wǎng)絲的應力平均為300 MPa。

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

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

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

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

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

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