基于CT圖像的混凝土三維微觀結(jié)構(gòu)在ANSYS中的實(shí)現(xiàn)

摘要：基于混凝土的CT平面圖像信息，用APDL語言生成的程序?qū)炷恋娜S微觀結(jié)構(gòu)在大型商業(yè)有限元飲件ANSYS中進(jìn)行三維重建。重建過程沒有生成點(diǎn)、線、面、體等幾何元素，而是利用CT平面圖像信息直接生成節(jié)點(diǎn)，連接節(jié)點(diǎn)生成單元，避免了網(wǎng)格無法剖分現(xiàn)象的出現(xiàn)。生成的單元形狀規(guī)則，長寬比合理。通過建立的更接近實(shí)際?昆凝土微觀結(jié)構(gòu)的數(shù)值模型，可利用有限元方法進(jìn)一步研究基于微觀結(jié)構(gòu)的混凝土各項(xiàng)力學(xué)特性。

關(guān)鍵詞：三維重建；混凝土；CT；ANsYs；微觀結(jié)構(gòu)

中圖分類號：TU528.01    文獻(xiàn)標(biāo)志碼：  A    文章編號： 1002—3550(2009)03一0013—03

0引言

　　我國很多高拱壩位于強(qiáng)震區(qū)，對大體積混凝土的動(dòng)態(tài)特性研究尤為重要。由于全級配大壩混凝土強(qiáng)度試驗(yàn)，既費(fèi)錢費(fèi)時(shí)又受到設(shè)備條件的限制，且少量試驗(yàn)往往還難以得出規(guī)律性成果，為此近年來國內(nèi)外開始探索基于微觀力學(xué)的有限元數(shù)值模型對混凝土內(nèi)部強(qiáng)度破壞機(jī)理的研究?；诨炷廖⒂^結(jié)構(gòu)的有限元建模得到了迅速的發(fā)展：Bazant等人提出了隨機(jī)粒子模型^[1](random particle model)，考慮了粒子分布的隨機(jī)性以模擬混凝土的骨料。大連理工大學(xué)的王寶庭、宋玉普提出了剛體一彈簧模型(rigid body spring：model)^[2]；邢紀(jì)波等提出了梁一顆粒模型^[3]；清華大學(xué)的劉光廷等將混凝土看做是水泥砂漿、骨料和兩者間的結(jié)合帶構(gòu)成的三相復(fù)合材料^[4]。唐春安利用基于統(tǒng)計(jì)數(shù)學(xué)描述材料非均勻性方法建立數(shù)值仿真模型^[5]。以上方法大多是基于隨機(jī)抽樣方法和統(tǒng)計(jì)學(xué)的知識，所建立的模型皆為隨機(jī)模型。雖然在骨料的分布和形狀上已越來越接近混凝土的真實(shí)結(jié)構(gòu)，但不可避免的與混凝土的真實(shí)微觀結(jié)構(gòu)存在著差異。CT技術(shù)的出現(xiàn)，實(shí)現(xiàn)了在無損狀態(tài)下將混凝土內(nèi)部微觀結(jié)構(gòu)用數(shù)字化方法進(jìn)行呈現(xiàn)。如果能利用CT平面圖像的信息，建立真實(shí)或近似真實(shí)的微觀結(jié)構(gòu)有限元模型，將對混凝土數(shù)值模擬研究產(chǎn)生重要的作用。李曉軍通過CT圖像得到了瀝青混凝土的二維真實(shí)微觀結(jié)構(gòu)^[6]；岳中琦等利用CT圖像建立了有限元模型并進(jìn)行了平面受力分析^[7]。但以上研究大多基于二維模型，不能反應(yīng)材料真實(shí)的空間結(jié)構(gòu)。本文運(yùn)用CT圖像提供的平面數(shù)字信息，用大型商業(yè)有限元軟件ANSYS自帶的APDL語言編制程序，建立了混凝土三維微觀結(jié)構(gòu)的有限元模型，生成的模型更接近混凝土的真實(shí)結(jié)構(gòu)，有利于將數(shù)直計(jì)算結(jié)果與物理試驗(yàn)結(jié)果對比，進(jìn)一步修改計(jì)算模型和計(jì)算參數(shù)，研究大體積混凝土的各項(xiàng)力學(xué)特性。

1三維重建的方法

　　現(xiàn)有的三維重建方法可分為兩類㈣：基于表面的方法和基于體數(shù)據(jù)的方法。

1.1基于表面的方法

　　基于表面的方法是確定物體形狀的最基本方法，其具體形式有邊界輪廓線表示和表面曲面表示。邊界輪廓線表示是根據(jù)各斷層的平面圖像確定輪廓線，并將各層輪廓線按序排列來表示物體的表面形狀。表面曲面表示是在已知各斷層圖像輪廓線的基礎(chǔ)上，采用這些輪廓線來重建物體表面?；诒砻娴姆椒▋?yōu)點(diǎn)是可以采用比較成熟的計(jì)算機(jī)圖形學(xué)方法進(jìn)行顯示，計(jì)算量相對較小，重建速度陜。缺點(diǎn)是重建出三維圖形的保真性不好。

1.2基于體數(shù)據(jù)的方法

　　基于體數(shù)據(jù)的方法不需要確立物體的表面幾何形狀，無需構(gòu)造幾何曲面，而是直接基于體素進(jìn)行重建?；隗w數(shù)據(jù)的方法重建出三維圖形的保真性好，不丟失物體形狀結(jié)構(gòu)的細(xì)節(jié)，能準(zhǔn)確的反應(yīng)映出所包含的形體結(jié)構(gòu)，避免了重建過程中所造成的偽像痕跡。缺點(diǎn)是數(shù)據(jù)量相對較大，數(shù)據(jù)的處理時(shí)間長，對計(jì)算機(jī)硬件要求較高。

　　以上的三維重建方法，目的大多是建立起物體的三維視覺效果。要進(jìn)行有限元計(jì)算，若采用以上方法進(jìn)行三維重建，生成的三維模型的單元數(shù)量將非常巨大，現(xiàn)有的計(jì)算機(jī)硬件條件無法滿足其計(jì)算要求，甚至由于實(shí)際結(jié)構(gòu)形式復(fù)雜，無法進(jìn)行網(wǎng)格剖分。因此，本文采用三維重建中基于表面方法的思想，提取平面圖像信息，并作合理的簡化，不生成點(diǎn)、線、面、體等幾何元素而是直接生成節(jié)點(diǎn)，連接節(jié)點(diǎn)生成單元，在ANSYS建立三維有限元模型。

2 CT圖像的數(shù)字化處理

　　本文重建中所用混凝土試什尺寸為Φ60 mm×120 mm的圓柱體，用標(biāo)準(zhǔn)混凝土成形后鉆取的方式獲得。CT機(jī)采用宜昌仁和醫(yī)院的菲利普16排螺旋CT，分辨率為512×512。

2.1  CT圖像格式的轉(zhuǎn)換與讀寫

　　原始的CT圖像是采用DICOM標(biāo)準(zhǔn)進(jìn)行存儲的，以DICOM標(biāo)準(zhǔn)存儲的圖像讀人MATLAB后，不易進(jìn)行處理，因此本文進(jìn)行圖像的格式轉(zhuǎn)換，也就是將原始的CT圖像轉(zhuǎn)換為256色的BMP。圖像。本文中選用了連續(xù)的20層圖像進(jìn)行重建(以下各種格式的圖片皆為頂端，中間及底端的3張圖片)。原始DIc0M格式的CT圖像經(jīng)轉(zhuǎn)換后的BMP恪式圖像如圖l所示。

[Page]

2.2圖像增強(qiáng)

　　圖像增強(qiáng)就是根據(jù)某種應(yīng)用的需要，人為的突出圖像中的某些信息，從而抑制或消除另一些信息的處理過程。使輸入的圖像具有更好的圖像質(zhì)量，有利于分析與識別。在本義中刻意突出圖像中骨料與砂漿的信息，利用形態(tài)學(xué)濾波器(交替開閉運(yùn)算)，濾除圖像中的高灰度或低灰度的小干擾區(qū)域。用閥值分割的方法將骨料與砂漿區(qū)別開米。具體為取像素閥值d，像素值大于d的認(rèn)為足骨料區(qū)，相應(yīng)付置處的像素屬性值賦為1；小于d則認(rèn)為是砂漿及背景區(qū)，相應(yīng)位置處的像素屬性值賦為0。經(jīng)閥值分割后效果如圖2所示。

3 ANSYS中的有限元建模

3.1  將CT圖像讀入MATLAB

　　在MATLAB中，用imread(‘二維圖像地址’，‘bmp’)命令依次將連續(xù)BMP格式的CT圖像讀入，每層圖像以二維矩陣(345×345)的形式保存，數(shù)組元素值為CT圖像對應(yīng)位置處的像素屬性值，即經(jīng)閥值分割處理后骨料區(qū)域數(shù)組元素值為l，砂漿及背景區(qū)域?yàn)?(本文中忽略了空洞區(qū)域)。

3.2生成單元節(jié)點(diǎn)

　　用APDL語言編制程序，將MATLAB中的像素屬性值矩陣讀入ANSYS的數(shù)組中，數(shù)組名為MSG。若每個(gè)像素點(diǎn)處生成1個(gè)節(jié)點(diǎn)，去掉邊緣背景部分則每層大約有l(wèi)00000節(jié)點(diǎn)，三維模型建成后單元數(shù)目很大，單機(jī)無法滿足其計(jì)算要求，需要對生成的節(jié)點(diǎn)數(shù)目進(jìn)行壓縮。本試驗(yàn)中試件的直徑為30mm，每行大約有345個(gè)像素格子，故每像素格的寬度約為0.1 mm。在每相鄰的9個(gè)像素點(diǎn)所占的區(qū)域中心處生成一個(gè)節(jié)點(diǎn)，則節(jié)點(diǎn)信息矩陣MSG被壓縮為(115×115)，節(jié)點(diǎn)與節(jié)點(diǎn)的間隔為0.3 mm，共生成約10000個(gè)節(jié)點(diǎn)。如圖3所示。

　　依據(jù)數(shù)組MSG的信息編制程序用APDL中的*if命令，對9個(gè)像素格所占的區(qū)域進(jìn)行判斷，如圖3所示，相鄰的9個(gè)區(qū)域中，若有多于5個(gè)像素區(qū)域?qū)儆诠橇蠀^(qū)域則生成的節(jié)點(diǎn)為骨料單元的節(jié)點(diǎn)；若多于5個(gè)點(diǎn)屬于砂漿區(qū)域則生產(chǎn)的節(jié)點(diǎn)為砂漿單元的節(jié)點(diǎn)。將生成節(jié)點(diǎn)的坐標(biāo)及屬性生信息保存在數(shù)組RMSG中，每行前三列元素值對應(yīng)節(jié)點(diǎn)的X、Y、Z坐標(biāo)值，第四個(gè)元素為節(jié)點(diǎn)的屬性直。通過壓縮節(jié)點(diǎn)后每層圖像大約生成l0000個(gè)節(jié)點(diǎn)。[Page]

3.3 ANSYS按層生成節(jié)點(diǎn)

　　將RMSG中坐標(biāo)信息讀入，運(yùn)用APDL命令：N，INUM，X，Y，Z生成第一層節(jié)點(diǎn)。然后運(yùn)用APDL命令：NGEN，ITIME，INE：，NODE1，NODE2，NINC，DX，DY，DZ，SPACE，復(fù)制其余各層節(jié)點(diǎn)，層問問隔為CT掃掃描像的間隔，本例中為0.84 mm，層問對應(yīng)位置的節(jié)點(diǎn)編號間隔為14 000。程序代碼為：

　　*DO，i，l，115    !生成第一層節(jié)點(diǎn)

　　*D0，j，1，115

　?。疨REP7

　　n，，RMSG1(X)，RMSG1(Y)，

　　*ENDDO

　　*ENDD)

　　NGEN，20，14 000，1，13 225…，-1.68 1復(fù)制其余19層節(jié)點(diǎn)

3.4連接相應(yīng)節(jié)點(diǎn)生成單元

　　選用ANSYS自帶的SOLID45單元。SOLID45為8節(jié)點(diǎn)單元，用APDL命令，E，i，j，k，l，m，n，o，p，由每層的4個(gè)節(jié)點(diǎn)與其相鄰的一層的對應(yīng)4個(gè)節(jié)點(diǎn)共8個(gè)節(jié)點(diǎn)生成一個(gè)SOLID45單元。程序代碼為：

　　*DO，f，0，18    !層控制

　　*DO，m，0，1 13   !行控制

　　*DO，n，1，114   !列控制

　　i=14 000*f+m*l 15+n  !層問節(jié)點(diǎn)編號增量為14 000

　　j=14 000*f+(m+1)*115+n

　　k=14 000*f.+m* l15+n十1

　　l=14 000*f+(m+1)*115+n+1

　　ii=14 000+i

　　jj=14 000+j

　　kk=14 000+k

　　11=14 000+1

　?。痯rep7

　　E，ii，jj，11，kk，i，j，l，k

　　*ENDDO

　　*ENDDO

　　*ENDDO

　　程序運(yùn)行后生成的單元為規(guī)則的長療體，長寬高之比約為3：3：8。規(guī)則的單元形狀和相對理想的長寬比，為后續(xù)數(shù)值計(jì)算的精度提供了保障。

3.5判別每個(gè)單元的8個(gè)節(jié)點(diǎn)屬性，賦單元力學(xué)參數(shù)

　　用APDL中的*get命令提取每個(gè)單元的節(jié)點(diǎn)編號，從數(shù)組RMSG中提取對應(yīng)的節(jié)點(diǎn)屬性值，編制程序節(jié)點(diǎn)屬性進(jìn)行判斷，若同一單元的8個(gè)節(jié)點(diǎn)都為骨料節(jié)點(diǎn)，則該單元賦骨料的力學(xué)參數(shù)；若8個(gè)節(jié)點(diǎn)都為砂漿的節(jié)點(diǎn)，則諺單元賦砂漿材料的力學(xué)參數(shù)；其他情況賦界面的力學(xué)參數(shù)。以下是用連續(xù)的20幅CT圖像生成的混凝土三維有限元模型圖。

4結(jié)論

　　(1)不同于現(xiàn)有基于隨機(jī)抽樣原理的混凝土微觀建模方法，提出了基于CT圖像信息三維有限元建模方法。骨料單元的位置依據(jù)真實(shí)試件的CT平面圖像信息產(chǎn)生，從而使有限元模型更接近于真實(shí)的}昆凝土結(jié)構(gòu)，數(shù)值模擬的精度進(jìn)一步提高。

　　(2)三維重建中沒有生成點(diǎn)、線、面、體等幾何元素，而是由圖像信息直接轉(zhuǎn)化為節(jié)點(diǎn)信息，連接節(jié)點(diǎn)生成單元。生成的單元形狀規(guī)則，長寬高比合理，避免了單元剖分中出現(xiàn)過多的尖角單元及單元無法剖分現(xiàn)象。 

　　(3)將界面單元引入模型，為混凝土的損傷斷裂分析提供了方便。雖然目前生成的界面單元尺寸仍然很大，但通過改進(jìn)程序可達(dá)到滿意的尺寸。

　　(4)5安照本文現(xiàn)有方法進(jìn)行有限元建模，若試件尺寸較大則會產(chǎn)生單機(jī)硬件條件無法滿足的計(jì)算量，但通過改進(jìn)程序，可將單元尺寸放大或在一定的尺寸要求范圍內(nèi)通過把同種材料的相鄰單元進(jìn)行合并，來減少單元數(shù)量。也可基于本文的建模方法進(jìn)行并行計(jì)算，提高計(jì)算能力。