本發(fā)明涉及巖石力學(xué)與工程、采礦工程領(lǐng)域，具體涉及一種脆性材料熱破裂的數(shù)值模擬方法。

背景技術(shù)：

熱應(yīng)力是引起巖石等脆性材料破壞的一個(gè)重要因素，當(dāng)熱應(yīng)力超過(guò)材料自身的強(qiáng)度極限，就會(huì)使材料產(chǎn)生破裂。熱破裂會(huì)對(duì)巖石的彈性和機(jī)械破裂性質(zhì)造成很大的影響，并引起巖石物性參數(shù)(如孔隙度、滲透率等)的變化，從而對(duì)許多工程造成重大影響。例如：在核廢料的存儲(chǔ)中，由于核廢料裂變而使圍巖溫度顯著升高，導(dǎo)致巖石的熱破裂，地下水的滲入會(huì)進(jìn)一步加劇圍巖破壞、甚至導(dǎo)致核素遷移，造成地下水污染；在石油開(kāi)采中，利用巖石的熱破裂，增加巖石的滲透性。目前，脆性材料熱破裂的數(shù)值模擬研究方法主要有不連續(xù)變形分析法、有限元法、顆粒法等，這些方法普遍存在部分力學(xué)參數(shù)很難準(zhǔn)確估算，無(wú)法對(duì)工程實(shí)際進(jìn)行準(zhǔn)確的模擬和仿真。

傳統(tǒng)有限元方法模擬裂紋擴(kuò)展時(shí)通常需要在裂尖構(gòu)建奇異單元，且隨著裂紋擴(kuò)展不斷地進(jìn)行網(wǎng)格重構(gòu)，該方法實(shí)現(xiàn)困難且低效。邊界元法雖然避免了網(wǎng)格重新劃分的問(wèn)題，但嚴(yán)重依賴(lài)于問(wèn)題的基本解，對(duì)于非線性、非均質(zhì)等問(wèn)題其優(yōu)勢(shì)大大降低。本發(fā)明的裂紋處理采用三角單元網(wǎng)格開(kāi)裂技術(shù)，該技術(shù)只對(duì)裂尖局部單元進(jìn)行網(wǎng)格調(diào)整，不涉及不連續(xù)位移場(chǎng)的處理問(wèn)題，也無(wú)須引入額外自由度，更加簡(jiǎn)便、高效，具有較好的適用性。

技術(shù)實(shí)現(xiàn)要素：

針對(duì)現(xiàn)有技術(shù)中的上述不足，本發(fā)明旨在提供一種脆性材料熱破裂的數(shù)值模擬方法，對(duì)脆性材料在熱應(yīng)力作用下產(chǎn)生的拉張破壞過(guò)程進(jìn)行模擬，為分析脆性材料的熱破壞過(guò)程和機(jī)制提供一種新的有限元方法。

為了達(dá)到上述發(fā)明目的，本發(fā)明采用的技術(shù)方案為：

提供一種脆性材料熱破裂的數(shù)值模擬方法，其包括以下步驟，

采用工程對(duì)象數(shù)字模型對(duì)脆性材料進(jìn)行計(jì)算網(wǎng)格劃分；

設(shè)置模擬時(shí)的材料參數(shù)、溫度邊界、應(yīng)力限定邊界和位移限定邊界；

計(jì)算應(yīng)力場(chǎng)：

{σ}＝[d]([b]·{δ}^e-{ε0})

其中，[d]為彈性矩陣；[b]為應(yīng)變矩陣；{ε0}為溫度變化引起的變形；{δ}^e為單元位移；

當(dāng)計(jì)算網(wǎng)格中存在一個(gè)節(jié)點(diǎn)的第一主應(yīng)力大于等于脆性材料抗拉強(qiáng)度時(shí)，對(duì)該節(jié)點(diǎn)所在處的待開(kāi)裂單元進(jìn)行開(kāi)裂處理得到新形成的計(jì)算網(wǎng)格，并返回計(jì)算應(yīng)力場(chǎng)步驟；

當(dāng)存在多個(gè)節(jié)點(diǎn)的第一主應(yīng)力大于等于脆性材料抗拉強(qiáng)度時(shí)，對(duì)承受最大第一主應(yīng)力的節(jié)點(diǎn)所在處的待開(kāi)裂單元進(jìn)行開(kāi)裂處理得到新形成的計(jì)算網(wǎng)格，并返回計(jì)算應(yīng)力場(chǎng)步驟；

當(dāng)所有節(jié)點(diǎn)的第一主應(yīng)力均小于脆性材料抗拉強(qiáng)度，且模擬時(shí)間未達(dá)到設(shè)置模擬時(shí)長(zhǎng)時(shí)，按設(shè)定值遞增溫度邊界后返回計(jì)算應(yīng)力場(chǎng)步驟，否則停止模擬操作；

所述待開(kāi)裂單元為包含第一主應(yīng)力大于等于其抗拉強(qiáng)度的節(jié)點(diǎn)的單元、且位于裂紋擴(kuò)展方向上的單元。

本發(fā)明的有益效果為：本方案能夠通過(guò)設(shè)置的材料參數(shù)和各種邊界條件模擬脆性材料由不同熱膨脹系數(shù)產(chǎn)生的熱應(yīng)力破壞模式，并對(duì)遭到破壞的單元進(jìn)行開(kāi)裂處理，以降低熱破裂對(duì)脆性材料的彈性和機(jī)械破裂性質(zhì)造成的影響；采用該數(shù)值模擬方法能夠準(zhǔn)確反映出在溫度荷載下，脆性材料拉張破壞過(guò)程，并能有效地分析裂紋開(kāi)裂規(guī)律和應(yīng)力轉(zhuǎn)移變化過(guò)程。

本發(fā)明的裂紋處理過(guò)程中，裂紋可以直接劈開(kāi)一個(gè)單元，或沿單元邊界擴(kuò)展，因此裂紋能夠不受初始網(wǎng)格的限制沿任意路徑擴(kuò)展；與現(xiàn)有的網(wǎng)格重構(gòu)算法相比，該方法只須對(duì)裂尖局部單元進(jìn)行網(wǎng)格開(kāi)裂或節(jié)點(diǎn)移動(dòng)，更加簡(jiǎn)便、高效，具有較好的適用性。

附圖說(shuō)明

圖1為脆性材料熱破裂的數(shù)值模擬方法一個(gè)實(shí)施例的流程圖。

圖2為脆性材料熱破裂的數(shù)值模擬方法的三角形單元開(kāi)裂示意圖。

圖3為脆性材料中內(nèi)嵌顆粒圓形試樣模型示意圖。

圖4a為t＝350℃時(shí)內(nèi)嵌顆粒圓形試樣第一主應(yīng)力示意圖。

圖4b為t＝370℃時(shí)內(nèi)嵌顆粒圓形試樣第一主應(yīng)力示意圖。

圖5為t＝390℃時(shí)內(nèi)嵌顆粒圓形試樣裂紋擴(kuò)展過(guò)程和應(yīng)力轉(zhuǎn)移過(guò)程第一主應(yīng)力示意圖。

圖6為t＝400℃時(shí)內(nèi)嵌顆粒圓形試樣裂紋擴(kuò)展過(guò)程和應(yīng)力轉(zhuǎn)移過(guò)程第一主應(yīng)力示意圖。

圖7為試樣模擬結(jié)果與實(shí)驗(yàn)結(jié)果效果對(duì)比示意圖。

具體實(shí)施方式

下面對(duì)本發(fā)明的具體實(shí)施方式進(jìn)行描述，以便于本技術(shù)領(lǐng)域的技術(shù)人員理解本發(fā)明，但應(yīng)該清楚，本發(fā)明不限于具體實(shí)施方式的范圍，對(duì)本技術(shù)領(lǐng)域的普通技術(shù)人員來(lái)講，只要各種變化在所附的權(quán)利要求限定和確定的本發(fā)明的精神和范圍內(nèi)，這些變化是顯而易見(jiàn)的，一切利用本發(fā)明構(gòu)思的發(fā)明創(chuàng)造均在保護(hù)之列。

參考圖1，圖1示意性地示出了脆性材料熱破裂的數(shù)值模擬方法一個(gè)實(shí)施例的流程圖；如圖1所示，該方法包括步驟101至步驟106。

在步驟101中，采用工程對(duì)象數(shù)字模型對(duì)脆性材料進(jìn)行計(jì)算網(wǎng)格劃分；其中，模型中的網(wǎng)格被劃分為三角形三節(jié)點(diǎn)單元。本方案中構(gòu)建的模型是采用數(shù)值分析通用圖形用戶(hù)界面進(jìn)行創(chuàng)建的。

在步驟102中，設(shè)置模擬時(shí)的材料參數(shù)、溫度邊界、應(yīng)力限定邊界和位移限定邊界；其中的材料參數(shù)包括材料密度、彈性模量、泊松比、熱膨脹系數(shù)、熱傳導(dǎo)系數(shù)和抗拉強(qiáng)度。

在步驟103中，計(jì)算應(yīng)力場(chǎng)：

{σ}＝[d]([b]·{δ}^e-{ε0})

其中，[d]為彈性矩陣，其為由彈性模量e和泊松比ν控制的各向同性矩陣；[b]為應(yīng)變矩陣；{ε0}為溫度變化引起的變形；{δ}^e為單元位移。

由于應(yīng)力場(chǎng)是基于溫度場(chǎng)和位移場(chǎng)進(jìn)行獲得的，本方案在獲得應(yīng)力場(chǎng)之前還需要對(duì)溫度場(chǎng)和位移場(chǎng)進(jìn)行求解。

在本發(fā)明的一個(gè)實(shí)施例中，溫度場(chǎng)的獲取方法為：

溫度場(chǎng)計(jì)算采用fourier導(dǎo)熱定律：

其中，ρ為材料密度kg/m³；k為在x、y方向的熱傳導(dǎo)系數(shù)(w/m·℃)；c為比熱值(j/m³·℃)；t為時(shí)間；t為溫度(℃)。

由于本方案主要針對(duì)穩(wěn)態(tài)熱傳導(dǎo)進(jìn)行模擬，所以溫度不隨時(shí)間變化時(shí)，則kii▽²t＝0，此時(shí)，i＝1,2，將得到的熱傳導(dǎo)系數(shù)帶入fourier導(dǎo)熱定律公式，得到溫度場(chǎng)t。

位移場(chǎng)的獲取方法為，采用常規(guī)有限單元法求解位移場(chǎng)：

[k]·{δ}＝{q}t

其中，[k]為剛度矩陣，{q}t為熱荷載(受熱膨脹而形成的相當(dāng)荷載)。

實(shí)施時(shí)，計(jì)算應(yīng)力采用的由于溫度變化引起的變形{ε0}的表達(dá)式為：

{ε0}＝αt[110]

其中，α為材料的熱膨脹系數(shù)，t為溫度的變化。

在步驟104中，當(dāng)計(jì)算網(wǎng)格中存在一個(gè)節(jié)點(diǎn)的第一主應(yīng)力大于等于脆性材料抗拉強(qiáng)度時(shí)，對(duì)該節(jié)點(diǎn)所在處的待開(kāi)裂單元進(jìn)行開(kāi)裂處理得到新形成的計(jì)算網(wǎng)格，并返回計(jì)算應(yīng)力場(chǎng)步驟。

在步驟105中，當(dāng)存在多個(gè)節(jié)點(diǎn)的第一主應(yīng)力大于等于脆性材料抗拉強(qiáng)度時(shí)，對(duì)承受最大第一主應(yīng)力的節(jié)點(diǎn)所在處的待開(kāi)裂單元進(jìn)行開(kāi)裂處理得到新形成的計(jì)算網(wǎng)格，并返回計(jì)算應(yīng)力場(chǎng)步驟。

其中，待開(kāi)裂單元為包含第一主應(yīng)力大于等于其抗拉強(qiáng)度的節(jié)點(diǎn)的單元、且位于裂紋擴(kuò)展方向上的單元。

本方案的待開(kāi)裂單元是基于拉張破壞原理而確定的，其中拉張破壞原理為：當(dāng)節(jié)點(diǎn)第一主應(yīng)力大于等于巖石的抗拉強(qiáng)度時(shí)，巖石發(fā)生拉張破壞；如果有多個(gè)節(jié)點(diǎn)的第一主應(yīng)力大于等于抗拉強(qiáng)度時(shí)，承受最大第一主應(yīng)力的點(diǎn)破裂。

在本發(fā)明的一個(gè)實(shí)施例中，第一主應(yīng)力的表達(dá)式為：

其中，σ1為第一主應(yīng)力；σx為x軸方向上的應(yīng)力；σy為y軸方向上的應(yīng)力；τxy為y軸上切應(yīng)力；

第一主應(yīng)力與x軸間夾角為：

實(shí)施時(shí)，本方案優(yōu)選對(duì)待開(kāi)裂單元進(jìn)行開(kāi)裂處理進(jìn)一步包括：獲取裂紋擴(kuò)展方向與待開(kāi)裂單元的開(kāi)裂邊形成的交點(diǎn)；增加與第一主應(yīng)力大于等于脆性材料抗拉強(qiáng)度的節(jié)點(diǎn)相重合的新增節(jié)點(diǎn)；根據(jù)新增節(jié)點(diǎn)和交點(diǎn)，將待開(kāi)裂單元沿裂紋擴(kuò)展方向裂開(kāi)。

由于單元開(kāi)裂后，可能引起部分單元幾何性質(zhì)不良，特別當(dāng)裂紋擴(kuò)展方向與單元邊之間的夾角很小時(shí)，將會(huì)形成較差的計(jì)算網(wǎng)格。

實(shí)施時(shí)，本方案優(yōu)選在新形成的計(jì)算網(wǎng)格與返回計(jì)算應(yīng)力場(chǎng)步驟之間還包括對(duì)新形成的計(jì)算網(wǎng)格進(jìn)行修正處理。采用修正處理后，能夠保證計(jì)算網(wǎng)格計(jì)算的精度。

在本發(fā)明的一個(gè)實(shí)施例中，修正處理時(shí)的具體實(shí)現(xiàn)方法為：

當(dāng)單元與裂紋擴(kuò)展方向之間的夾角小于設(shè)定值時(shí)，判斷單元短邊與長(zhǎng)邊的比值是否小于設(shè)定最小比值(小于1的實(shí)值)：

若小于，則將計(jì)算網(wǎng)格的邊移動(dòng)至裂紋擴(kuò)展方向，通過(guò)合并節(jié)點(diǎn)刪除夾角小于設(shè)定值的單元。

在步驟106中，當(dāng)所有節(jié)點(diǎn)的第一主應(yīng)力均小于脆性材料抗拉強(qiáng)度，且未達(dá)到設(shè)置模擬時(shí)長(zhǎng)時(shí)，按設(shè)定值(時(shí)間步)遞增溫度邊界后返回計(jì)算應(yīng)力場(chǎng)步驟，否則停止模擬操作。

為實(shí)現(xiàn)對(duì)開(kāi)裂處理過(guò)程的數(shù)值模擬結(jié)果的處理，本方案根據(jù)時(shí)間步來(lái)自動(dòng)保存數(shù)據(jù)文件，即每隔設(shè)定時(shí)間步長(zhǎng)自動(dòng)保存該時(shí)間步長(zhǎng)內(nèi)計(jì)算網(wǎng)格的單元數(shù)據(jù)，之后再為下一溫度載荷計(jì)算做準(zhǔn)備。

下面結(jié)合附圖2詳細(xì)描述三角形三節(jié)點(diǎn)單元的開(kāi)裂過(guò)程：

①判斷開(kāi)裂單元。假如節(jié)點(diǎn)n5滿(mǎn)足拉張破壞原理，查找包含節(jié)點(diǎn)n5的所有單元，在裂紋擴(kuò)展方向上的單元即為需要開(kāi)裂的單元，即單元e1-3-5和e7-5-9；

②增加交點(diǎn)。沿裂紋擴(kuò)展方向找到需要開(kāi)裂單元開(kāi)裂邊上的交點(diǎn)n2和n8；

③增加重合節(jié)點(diǎn)。增加與開(kāi)裂節(jié)點(diǎn)n5重合的節(jié)點(diǎn)n10；

④劈開(kāi)單元。根據(jù)新增節(jié)點(diǎn)(包括重合節(jié)點(diǎn))，將單元沿裂紋擴(kuò)展方向劈開(kāi)，即開(kāi)裂單元e1-3-5和e7-5-9會(huì)在新增加的節(jié)點(diǎn)n2，n8和n10的基礎(chǔ)上劈開(kāi)為四個(gè)獨(dú)立的單元e1-2-5，e10-2-3和e7-5-8，e8-10-9；

⑤修正單元網(wǎng)格。非開(kāi)裂單元e3-5-6，e9-5-6單元信息會(huì)在新增重合節(jié)點(diǎn)的基礎(chǔ)上修改為e3-10-6，e9-10-6，單元e1-5-4和e7-5-4則保持不變。

下面結(jié)合具體的實(shí)例對(duì)本方案的數(shù)值模擬方法進(jìn)行說(shuō)明：

本實(shí)施例為一種由熱膨脹系數(shù)不同所導(dǎo)致的材料破裂模式的數(shù)值模擬形式，其具體步驟為：

(1)建立工程對(duì)象數(shù)字模型，如圖3所示，計(jì)算實(shí)例為帶有顆粒的圓形試樣，模型是由兩個(gè)同心圓組成，內(nèi)圓直徑為12mm，外圓直徑為25mm，內(nèi)、外同心圓分別表示兩種不同的材料，分別命名為顆粒和基質(zhì)。

(2)施加材料參數(shù)，顆粒和基質(zhì)的具體材料參數(shù)見(jiàn)表1；

表1材料力學(xué)參數(shù)

(3)設(shè)置時(shí)間控制步長(zhǎng)，總時(shí)間t＝40s，時(shí)間步(每隔設(shè)定時(shí)間步長(zhǎng))δt＝1s；

(4)施加溫度邊界條件，實(shí)驗(yàn)設(shè)置模型基質(zhì)和中央的顆粒施加相同的溫度條件，初始溫度t＝0℃，以δt＝10℃的溫度增量不斷升溫。

(5)計(jì)算模型溫度場(chǎng)、位移場(chǎng)和應(yīng)力場(chǎng)，根據(jù)步驟(2)、步驟(3)所獲取的材料參數(shù)和邊界條件，通過(guò)有限元程序自動(dòng)生成系統(tǒng)fepg計(jì)算模型的各個(gè)場(chǎng)值。

(6)計(jì)算第一主應(yīng)力和方向角。將步驟(5)計(jì)算所得的應(yīng)力值帶入，進(jìn)一步計(jì)算第一主應(yīng)力及其方向角。

(7)判斷是否有單元開(kāi)裂。由強(qiáng)度準(zhǔn)則判斷是否有單元滿(mǎn)足開(kāi)裂條件，若有單元滿(mǎn)足開(kāi)裂條件，劈開(kāi)單元(開(kāi)裂流程詳見(jiàn)圖2)，重復(fù)步驟(5)～(7)；若沒(méi)有單元滿(mǎn)足開(kāi)裂條件，保存當(dāng)前時(shí)間步的開(kāi)裂結(jié)果，并判斷是否達(dá)到設(shè)置的總時(shí)間，若沒(méi)有達(dá)到設(shè)置的總時(shí)間，增加溫度增量δt，即溫度tn+1＝tn+δt，重復(fù)步驟(4)～(7)；若達(dá)到設(shè)置的總時(shí)間，結(jié)束程序計(jì)算。

(8)查看厚壁圓筒試樣在第一主應(yīng)力場(chǎng)下的開(kāi)裂結(jié)果。

圖4a為t＝350℃時(shí)內(nèi)嵌顆粒圓形試樣第一主應(yīng)力示意圖。圖4b為t＝370℃時(shí)內(nèi)嵌顆粒圓形試樣第一主應(yīng)力示意圖。圖5為t＝390℃時(shí)內(nèi)嵌顆粒圓形試樣開(kāi)裂過(guò)程第一主應(yīng)力示意圖。圖6為t＝400℃時(shí)內(nèi)嵌顆粒圓形試樣開(kāi)裂過(guò)程第一主應(yīng)力示意圖。

從這些視圖可以看出，在升溫過(guò)程中，由于基質(zhì)和顆粒具有不同的熱膨脹系數(shù)，造成了基質(zhì)和顆粒的界面上變形的不協(xié)調(diào)，繼而在顆粒周?chē)a(chǎn)生了應(yīng)力集中現(xiàn)象，升溫的初始階段，模型中應(yīng)力處于積聚過(guò)程，但不滿(mǎn)足開(kāi)裂條件，未有破壞現(xiàn)象發(fā)生；t＝390℃時(shí)，試樣基質(zhì)中靠近顆粒的區(qū)域開(kāi)始發(fā)生破壞，t＝400℃裂紋在原有的基礎(chǔ)上進(jìn)一步擴(kuò)展，最終形成一條主裂紋將基質(zhì)貫通破壞。

圖7為試樣模擬結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比圖，由圖7可以看出試樣最終破壞模式與實(shí)驗(yàn)結(jié)果相符合。