王 程，臧孟炎

(華南理工大學(xué) 機械與汽車工程學(xué)院，廣州 510640)

據(jù)公安部網(wǎng)站統(tǒng)計，2020年上半年，我國汽車保有量已經(jīng)達到2.7億輛。風(fēng)擋玻璃是重要的安全部件，由于生產(chǎn)過程中不可避免的質(zhì)量缺陷而導(dǎo)致強度分布不均，汽車在行駛過程中，巨大載荷所引起的車架變形或者風(fēng)壓所形成的較大風(fēng)載可能會導(dǎo)致玻璃的破壞，成為威脅司乘人員生命安全的一個重要因素[1-2]。所以，可靠高效的汽車風(fēng)擋玻璃的準(zhǔn)靜態(tài)破壞性能評價對于保證汽車行駛的安全性至關(guān)重要。

關(guān)于汽車玻璃的準(zhǔn)靜態(tài)破壞問題，相關(guān)學(xué)者做了大量研究工作。玻璃的準(zhǔn)靜態(tài)破壞問題研究常用的試驗方法是三點或四點彎曲方法，但由于玻璃存在質(zhì)量缺陷，這兩種方法無法準(zhǔn)確評價玻璃的準(zhǔn)靜態(tài)破壞性能。因此，VITMAN等[3]最早提出了用于一定邊長玻璃板的同軸雙環(huán)加載試驗，可準(zhǔn)確評價玻璃的準(zhǔn)靜態(tài)破壞性能。RITTER等[4]使用同軸雙環(huán)加載試驗，評價了不同尺寸玻璃的準(zhǔn)靜態(tài)破壞性能。GULATI等[5]對新制的和老化的汽車玻璃進行同軸雙環(huán)加載試驗，研究玻璃老化對汽車玻璃破壞強度的影響。MüLLER-BRAUN等[6]對邊長為125 mm的平板玻璃進行了同軸雙環(huán)試驗，同時使用ANSYS隱式數(shù)值分析方法，研究玻璃的徑向及切向破壞強度的分布。CASTORI等[7]對邊長為400 mm的玻璃板做了不同支撐情況下的同軸雙環(huán)試驗，同時進行隱式有限元仿真分析等效應(yīng)力分布，研究一種能夠預(yù)測玻璃破壞強度的數(shù)學(xué)模型。綜上所述，對于玻璃的準(zhǔn)靜態(tài)破壞問題研究，多數(shù)采用試驗方法。即使有限元顯式分析方法在風(fēng)擋玻璃沖擊破壞方面得到了廣泛應(yīng)用[8]，但玻璃準(zhǔn)靜態(tài)破壞的仿真分析尚無有效計算方法對應(yīng)。主要原因在于準(zhǔn)靜態(tài)破壞問題現(xiàn)象時間過長，顯式有限元方法難以實現(xiàn)，隱式有限元方法在破壞發(fā)生后的收斂性難以保證且計算效率極低。所以，如何準(zhǔn)確且高效地對玻璃準(zhǔn)靜態(tài)加載破壞進行仿真分析，對提升汽車玻璃的開發(fā)設(shè)計水平具有重大意義。

基于上述考慮，提出玻璃準(zhǔn)靜態(tài)破壞仿真的隱式和顯式組合仿真方法，充分發(fā)揮兩種計算方法的優(yōu)勢。即在破壞前的加載階段采用隱式有限元方法以較大的時間步長分析玻璃的變形，在即將發(fā)生破壞時改用顯式有限元方法分析破壞現(xiàn)象。以玻璃的同軸雙環(huán)加載試驗為研究對象，使用內(nèi)聚力模型描述玻璃裂紋的發(fā)生和傳播，比較組合仿真分析方法和單純隱式分析方法的計算精度和計算效率，為汽車風(fēng)擋玻璃準(zhǔn)靜態(tài)破壞性能評價的有限元仿真應(yīng)用提供指導(dǎo)和幫助。

1 內(nèi)聚力模型

內(nèi)聚力模型法是一種基于斷裂力學(xué)，預(yù)測脆性材料的動態(tài)斷裂的有效方法，由于可以描述裂紋的擴展和考慮破壞發(fā)生的能量耗散，近年來廣泛應(yīng)用于汽車風(fēng)擋玻璃等脆性材料的沖擊破壞仿真研究工作中[8-10]。內(nèi)聚力模型首先由DUGDALE[11]和BARENBLATT[12]提出，BARENBLATT將內(nèi)聚力模型應(yīng)用于脆性材料的斷裂中。HILLERBORG等[13]首次將內(nèi)聚力模型應(yīng)用到有限元計算中，模擬了脆性材料的斷裂過程。

雙線性型固有內(nèi)聚力模型由MI等[14]提出，其法向和切向的本構(gòu)曲線如圖1所示。其牽引-分離法則假定：在外載荷的作用下，界面上作用牽引力隨著裂紋面間相對位移的增大而線性增加。當(dāng)內(nèi)聚區(qū)域的牽引力達到最大牽引力之后，內(nèi)聚區(qū)的裂紋面相對位移會持續(xù)增加，但界面上的作用牽引力開始線性下降，此時材料剛度下降，發(fā)生軟化，裂紋開始萌生及擴展。當(dāng)牽引力降為0時，裂紋面相對位移達到最大值，此時材料已經(jīng)發(fā)生破壞。

圖1 雙線性型內(nèi)聚力模型的法向及切向本構(gòu)曲線

雙線性型固有內(nèi)聚力模型的本構(gòu)方程為

式中：Tn和Tt分別表示裂紋內(nèi)聚區(qū)域的法向牽引力和切向牽引力，對應(yīng)的破壞應(yīng)力為σmax和τmax；δn和δt分別為裂紋內(nèi)聚區(qū)域兩裂面間的法向及切向分離量；和分別為裂紋開始萌生時內(nèi)聚區(qū)域兩裂面間的法向及切向分離量；和分別為裂紋內(nèi)聚區(qū)域兩裂面間的法向及切向最終分離量。斷裂能為雙線性型固有內(nèi)聚力模型的本構(gòu)曲線與分離量坐標(biāo)軸包圍的三角形面積的大小，即：

而在實際的材料破壞中，材料的裂紋不是單一類型的破壞，而是某種混合模式。因此，為描述混合模式下的材料破壞，BENZEGGAGH和KENANE等[15]提出了一種開裂準(zhǔn)則，稱為BK開裂準(zhǔn)則。此時，內(nèi)聚區(qū)域裂紋的總斷裂能為：

式中：β為破壞模式的混合度，一般為0.5；XMU為材料常數(shù)。這種材料開裂準(zhǔn)則由于簡單易用，所以廣泛地用于復(fù)合材料的仿真破壞中。目前，有限元商用分析軟件ABAQUS和LS-DYNA的隱式和顯式計算方法中，均已嵌入表征玻璃脆性破壞的內(nèi)聚力模型。

2 組合仿真分析方法

仔細(xì)研究夾層玻璃的準(zhǔn)靜態(tài)破壞問題，發(fā)現(xiàn)可以將這個過程分解為加載階段和破壞階段。其中加載階段往往需要10～102 s甚至更長的現(xiàn)象時間，而破壞階段的現(xiàn)象時間不足1 s。這是因為作為脆性材料的玻璃，發(fā)生破壞時裂紋的擴展速度極快[16]。如果單純使用隱式求解器分析夾層玻璃的準(zhǔn)靜態(tài)破壞過程，則由于玻璃破壞的現(xiàn)象時間極短，極易導(dǎo)致玻璃破壞分析過程中出現(xiàn)收斂性問題。即便使用其他方法使其收斂，計算步長也會大幅度縮短，收斂性迭代計算效率將嚴(yán)重下降。反之，若單純使用顯式求解器分析，雖然無需考慮收斂性問題，但玻璃準(zhǔn)靜態(tài)加載的過長現(xiàn)象時間對應(yīng)顯式分析方法微秒級的時間步長，導(dǎo)致需要計算的步長數(shù)過多，發(fā)生數(shù)據(jù)溢出無法計算。所以，如果使用隱式分析和顯式分析相結(jié)合的方法，即加載階段使用隱式分析方法，而破壞階段使用顯式分析方法，就能充分發(fā)揮兩個方法的優(yōu)勢：加載階段使用隱式有限元方法可以有較大的時間步長，破壞階段使用顯式有限元方法避免了復(fù)雜的收斂性計算。

基于ABAQUS軟件的玻璃準(zhǔn)靜態(tài)破壞現(xiàn)象的組合仿真分析方法如下：

（1）在ABAQUS的UVARM子程序中，定義玻璃內(nèi)聚力單元等效拉應(yīng)力破壞強度，且滿足破壞強度后立即終止計算。通過編譯得到新的ABAQUS執(zhí)行文件。

（2）對嵌入內(nèi)聚單元的玻璃準(zhǔn)靜態(tài)破壞仿真模型M1進行隱式有限元計算。當(dāng)某個內(nèi)聚單元滿足破壞強度，UVARM子程序命令計算立即終止，得到夾層玻璃破壞開始時間。

（3）以玻璃破壞開始時間為計算結(jié)束時間，定義新的玻璃準(zhǔn)靜態(tài)破壞有限元模型M2，重新進行隱式計算。

（4）將M2計算最終時刻得到的玻璃節(jié)點位移和單元應(yīng)力作為初始條件導(dǎo)入到顯式計算仿真模型M3中，分析玻璃的破壞階段。

3 數(shù)值仿真與評價

3.1 仿真模型

基于玻璃的同軸雙環(huán)加載試驗[7]，考慮結(jié)構(gòu)和載荷的對稱性，建立邊長為400 mm，厚度為8 mm的1/4同軸雙環(huán)加載試驗仿真模型。其中厚度8 mm玻璃采用4層全積分六面體單元，網(wǎng)格分布為輻射狀網(wǎng)格，玻璃單元數(shù)量為24 256個，如圖2所示。在玻璃網(wǎng)格公共界面間插入零厚度的內(nèi)聚力單元，使用內(nèi)聚力單元模擬玻璃的破壞。加載環(huán)和支撐環(huán)的半徑分別為75 mm和150 mm，它們的截面半徑為2.5 mm。支撐環(huán)定義固定約束，給加載環(huán)定義0.023 mm/s的向下恒定速度，得到玻璃準(zhǔn)靜態(tài)加載模型M1，如圖3所示。

圖2 玻璃網(wǎng)格劃分

圖3 厚度8 mm玻璃同軸雙環(huán)加載模型

玻璃選擇線彈性材料模型。加載環(huán)和支撐環(huán)定義為剛體，材料參數(shù)均參照普通鋼材。材料物性見表1[7]。

表1 玻璃、剛體環(huán)的材料物性

玻璃的內(nèi)聚力模型定義為雙線性型固有內(nèi)聚力模型，破壞強度取99 MPa[7]， I型破壞的斷裂能為10 N/m，II型、III型破壞的斷裂能為50 N/m[8]。內(nèi)聚力單元的破壞準(zhǔn)則使用BK準(zhǔn)則，XMU系數(shù)為2，罰剛度取500 GPa/mm[8]。

加載環(huán)和支撐環(huán)與玻璃間的接觸，定義為面與面接觸（Surface to Surface Contact），摩擦因數(shù)為0.1。在破壞階段顯式分析模型中，添加玻璃單元自接觸，類型為通用接觸（General Contact），摩擦因數(shù)為0.9，已將玻璃裂紋間的自接觸考慮在內(nèi)。

3.2 組合仿真計算

首先以M1模型隱式計算確定玻璃的初始破壞發(fā)生時間，然后以M2模型第2次隱式計算得到玻璃破壞發(fā)生前的等效應(yīng)力狀態(tài)，如圖4所示。

圖4 破壞發(fā)生前的等效應(yīng)力云圖

讀入隱式分析結(jié)果，建立破壞階段顯式分析仿真模型M3，進行仿真分析，最終得到玻璃的破壞模式、破壞載荷和破壞位移。

3.3 仿真結(jié)果分析及比較

3.3.1 破壞模式分析

厚度8 mm玻璃同軸雙環(huán)加載仿真模型使用組合計算方法和純隱式計算方法得到的玻璃破壞仿真結(jié)果與試驗結(jié)果，如圖5所示。由圖5可知，組合仿真計算方法得到的玻璃破壞模式圖5a與試驗結(jié)果圖5c呈現(xiàn)良好的一致性：加載環(huán)內(nèi)部區(qū)域玻璃發(fā)生全面粉碎性破壞，加載環(huán)和支撐環(huán)之間的玻璃產(chǎn)生大量徑向裂紋且加載環(huán)和支撐環(huán)附近有顯著的周向裂紋，而加載環(huán)外部區(qū)域玻璃的徑向裂紋明顯減少。由圖5b可知，純隱式仿真所得到的玻璃破壞情況僅在加載環(huán)內(nèi)部區(qū)域與試驗結(jié)果一致，其他區(qū)域的玻璃沒有觀察到破壞現(xiàn)象。

將圖3所示模型的玻璃厚度由8 mm變?yōu)? mm（刪除一層玻璃六面體單元），使用組合仿真方法得到的厚度6 mm玻璃雙環(huán)試驗仿真結(jié)果與試驗結(jié)果如圖6所示。由圖6可知，組合仿真方法仍然得到了與對應(yīng)試驗一致性良好的分析結(jié)果。

圖5 兩種仿真方法得到的破壞現(xiàn)象與試驗結(jié)果

圖6 厚度6 mm玻璃雙環(huán)試驗仿真與試驗結(jié)果

綜上所述，就破壞模式來看，純隱式計算方法顯然不能有效預(yù)測復(fù)雜的玻璃準(zhǔn)靜態(tài)破壞問題，而組合分析方法得到了一致性的仿真計算結(jié)果。

3.3.2 破壞載荷與破壞位移分析

表2為厚度8 mm與6 mm玻璃同軸雙環(huán)加載試驗和仿真所得到的破壞位移和破壞載荷對比。由表2可知，無論是破壞位移還是破壞載荷，仿真與試驗結(jié)果的誤差都在10%以內(nèi)，呈現(xiàn)良好的一致性，從定量方面說明了組合仿真分析方法的有效性。

表2 試驗與仿真得到的破壞位移和破壞載荷

3.3.3 計算效率分析

前述仿真結(jié)果使用具有Intel(R) Xeon(R) Silver 4210 2.2 GHz處理器的計算服務(wù)器，8核并行計算。表3為厚度8 mm和厚度6 mm玻璃模型分別使用組合分析方法和單純隱式分析方法的計算時間對比。盡管組合分析方法在加載階段計算了兩次。由表3可知，無論是厚度8 mm玻璃模型還是厚度6 mm玻璃模型，組合仿真方法所消耗的CPU時間要遠(yuǎn)小于單純使用隱式計算的總計算時間。這是因為玻璃同軸雙環(huán)加載試驗的破壞現(xiàn)象時間極短，而隱式算法在分析玻璃的破壞現(xiàn)象時，需要大幅降低計算步長進行迭代計算，而且一個步長的迭代計算需要多次迭代才能收斂，導(dǎo)致計算時間延長，計算效率降低。

表3 組合仿真和純隱式仿真的計算時間對比

4 結(jié)論

（1）針對玻璃準(zhǔn)靜態(tài)破壞仿真問題，提出了隱式和顯式組合仿真分析方法。通過用戶材料子程序確定夾層玻璃初始破壞發(fā)生時間，通過兩次隱式計算和一次顯式計算，利用隱式計算方法和顯式計算方法分別在玻璃加載階段和破壞階段的優(yōu)勢，實現(xiàn)提高精度和計算效率的目的。

（2）以厚度為8 mm和厚度為6 mm的玻璃同軸雙環(huán)準(zhǔn)靜態(tài)加載破壞試驗為研究對象，分別使用組合仿真分析方法和純隱式仿真方法進行了仿真分析。組合仿真分析方法的破壞模式與試驗結(jié)果基本一致，破壞位移與破壞載荷與試驗結(jié)果誤差小于10%，計算效率顯著高于純隱式計算方法，充分說明了組合仿真方法在玻璃準(zhǔn)靜態(tài)破壞現(xiàn)象仿真分析方面的有效性。