玻璃準靜態破壞現象隱式和顯式組合仿真研究

王 程，臧孟炎

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

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

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

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

1 內聚力模型

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

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

圖1 雙線性型內聚力模型的法向及切向本構曲線

雙線性型固有內聚力模型的本構方程為

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

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

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

2 組合仿真分析方法

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

基于ABAQUS軟件的玻璃準靜態破壞現象的組合仿真分析方法如下：

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

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

（3）以玻璃破壞開始時間為計算結束時間，定義新的玻璃準靜態破壞有限元模型M2，重新進行隱式計算。

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

3 數值仿真與評價

3.1 仿真模型

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

圖2 玻璃網格劃分

圖3 厚度8 mm玻璃同軸雙環加載模型

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

表1 玻璃、剛體環的材料物性

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

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

3.2 組合仿真計算

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

圖4 破壞發生前的等效應力云圖

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

3.3 仿真結果分析及比較

3.3.1 破壞模式分析

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

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

圖5 兩種仿真方法得到的破壞現象與試驗結果

圖6 厚度6 mm玻璃雙環試驗仿真與試驗結果

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

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

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

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

3.3.3 計算效率分析

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

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

4 結論

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

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