汪基偉，付 宇，王亞獎(jiǎng)

?

基于空間組合式單元模型的表面求交算法

汪基偉1，付 宇1，王亞獎(jiǎng)2

(1. 河海大學(xué)土木與交通學(xué)院，江蘇 南京 210098；2. 泉州市審計(jì)局，福建 泉州 362000)

在空間埋置組合式單元模型中，鋼筋單元可埋置于混凝土單元任何位置，混凝土單元網(wǎng)格剖分不受鋼筋位置的限制，方便實(shí)用，但需確定鋼筋單元兩端在混凝土單元表面的位置坐標(biāo)。因此，求解鋼筋線與混凝土單元表面的交點(diǎn)坐標(biāo)是應(yīng)用該單元模型的前提，現(xiàn)有的求解方法只適用于混凝土單元表面是平面的情況。為此，提出了牛頓迭代法和分塊解析法兩種處理方法，能求解鋼筋線與混凝土單元表面為任何形狀時(shí)的交點(diǎn)坐標(biāo)，增強(qiáng)了該模型的適用性。通過算例驗(yàn)證了這兩種方法的正確性。從適用性而言，分塊解析法要優(yōu)于牛頓迭代法。

鋼筋混凝土有限元；單元模型；前處理；曲面表面

現(xiàn)有鋼筋混凝土結(jié)構(gòu)有限元模型[1]中，分離式單元模型采用塊體單元模擬混凝土，桿單元模擬鋼筋[2-3]，桿單元必須在混凝土單元邊界上，結(jié)構(gòu)配筋復(fù)雜時(shí)網(wǎng)格剖分困難，有時(shí)甚至難以實(shí)現(xiàn)。整體式單元模型[4]將鋼筋彌散于整個(gè)單元中，利用單元參數(shù)(如配筋率)求得考慮鋼筋作用后的單元折算彈性模量，再以折算彈性模量求單元?jiǎng)偠染仃?，以考慮鋼筋對(duì)單元?jiǎng)偠染仃嚨呢暙I(xiàn)，該模型無法考慮鋼筋的具體位置，也不能求出鋼筋應(yīng)力分布。埋置式單元模型[5-6]將鋼筋埋置于混凝土單元內(nèi)部，網(wǎng)格剖分方便且能求得鋼筋應(yīng)力，單元?jiǎng)偠葹榛炷羻卧獎(jiǎng)偠扰c鋼筋單元?jiǎng)偠戎?，其中混凝土單元?jiǎng)偠染仃嚾圆捎玫葏卧膭偠染仃嚕摻顔卧獎(jiǎng)偠染仃噭t利用混凝土與鋼筋應(yīng)變相同的原則得到。文獻(xiàn)[7]和[8]分別推導(dǎo)出了二維和三維埋置式單元模型，鋼筋單元?jiǎng)偠染仃噑和應(yīng)力矩陣s表達(dá)式為

其中，s為鋼筋單元應(yīng)變矩陣；s為鋼筋彈性模量；s為鋼筋截面面積；d為桿單元任意點(diǎn)所在的微分段；為單元結(jié)點(diǎn)位移列陣。

由式(1)、(2)可知，埋置式單元模型必須得到鋼筋線與混凝土單元表面的交點(diǎn)坐標(biāo)。文獻(xiàn)[8]雖給出了求解方法，但該方法假定單元表面為平面，不適用于單元表面為曲面的情況。在鋼筋混凝土有限元網(wǎng)格剖分中，希望形成表面為平面的規(guī)則單元，但結(jié)構(gòu)形狀復(fù)雜時(shí)難以做到，剖分后必然會(huì)出現(xiàn)表面為曲面的單元，現(xiàn)有方法無法應(yīng)用。曲面網(wǎng)格一般是由圓柱、圓錐、球面等規(guī)則曲面以及Bézier、NURBS等自由曲面組合而成[9]，在鋼筋混凝土結(jié)構(gòu)中單元表面較為規(guī)則，采用高次單元能夠較好地逼近結(jié)構(gòu)的曲線和曲面邊界[10-11]，提高計(jì)算精度。為此，本文提出了兩種鋼筋線與混凝土單元表面求交方法，能有效得到鋼筋結(jié)點(diǎn)坐標(biāo)，適用于任意單元表面，擴(kuò)大了埋置組合單元模型的適用范圍。

1 鋼筋線與混凝土單元表面求交方法

在本文提出的兩種方法中，第1種是依據(jù)有限元中整體坐標(biāo)與局部坐標(biāo)的對(duì)應(yīng)關(guān)系，利用鋼筋線與混凝土單元表面交點(diǎn)對(duì)應(yīng)的局部坐標(biāo)中必有一項(xiàng)絕對(duì)值為1的條件，采用修正的牛頓法，建立求解交點(diǎn)坐標(biāo)的迭代公式，稱為牛頓迭代法。第2種稱為分塊解析法，可直接采用整體坐標(biāo)建立混凝土單元表面方程，然后與鋼筋直線方程聯(lián)立求解出交點(diǎn)的整體坐標(biāo)，再由交點(diǎn)整體坐標(biāo)求解其局部坐標(biāo)。依據(jù)局部坐標(biāo)中必有一項(xiàng)絕對(duì)值為1的條件，判斷交點(diǎn)求解精度，若精度不能滿足要求則分塊建立曲面方程求解。

1.1 牛頓迭代法

1.1.1 迭代法方程的建立

混凝土單元采用如圖1(a)所示的8~20結(jié)點(diǎn)等參單元，其中9~20結(jié)點(diǎn)為中間結(jié)點(diǎn)，可以刪除。

圖1 結(jié)點(diǎn)編號(hào)圖

由于鋼筋上任意一點(diǎn)均在混凝土單元內(nèi)，對(duì)任意表面形狀的混凝土等參單元，必然滿足

其中，N為有限元形函數(shù)；x、y和z為單元結(jié)點(diǎn)坐標(biāo)。

將式(4)代入式(3)，有

=–1表面對(duì)應(yīng)的形函數(shù)N表達(dá)式為

當(dāng)該面存在中間結(jié)點(diǎn)時(shí)，中間結(jié)點(diǎn)的形函數(shù)為

同時(shí)對(duì)中間結(jié)點(diǎn)相鄰的兩個(gè)角結(jié)點(diǎn)的形函數(shù)進(jìn)行修改，各減該中間結(jié)點(diǎn)形函數(shù)的一半。

將式(5)改寫成與、有關(guān)的方程，即

其中，16個(gè)待定系數(shù)由該平面的結(jié)點(diǎn)坐標(biāo)決定，若沒有中間結(jié)點(diǎn)，可將式(6)代入，有8個(gè)系數(shù)不為0，其中4個(gè)系數(shù)14、16、17、18為

另外4個(gè)系數(shù)24、26、27、28的表達(dá)式可由14、16、17、18表達(dá)式中的替換成得到，即

當(dāng)有中間結(jié)點(diǎn)5時(shí)，因中間結(jié)點(diǎn)5對(duì)角結(jié)點(diǎn)1和2有影響，式(8)中的系數(shù)11、13、21、23不再為0，同時(shí)式(9)中的系數(shù)17、18、27、28也需修正。式(10)給出11、13、17、18表達(dá)式，同樣將11、13、17、18表達(dá)式中的替換成可得到21、23、27、28表達(dá)式。當(dāng)有其他中間結(jié)點(diǎn)時(shí)，按照上述思路就可寫出相應(yīng)的表達(dá)式，即

采用Newton法建立迭代式求解非線性方程組式(8)，迭代式為

雅可比矩陣表達(dá)式及相應(yīng)系數(shù)為

1.1.2 迭代法方程的求解

求解時(shí)，首先要判斷鋼筋線與塊體單元表面方程有無交點(diǎn)，可將該面沿對(duì)角拆分構(gòu)成兩個(gè)平面，若鋼筋線與兩個(gè)平面的夾角都小于一定的限值(可取1°)，則認(rèn)為鋼筋與該表面平行，無交點(diǎn)。另外，還有兩種情況需特殊處理：

迭代初值的選擇直接影響到迭代法的收斂速度。鋼筋線與單元表面交點(diǎn)為一個(gè)或兩個(gè)，首先討論鋼筋與表面只有一個(gè)交點(diǎn)的情況，圖2給出了3種迭代初值的選擇方案，分別為：

圖2 迭代初值點(diǎn)分布圖

圖3是尺寸為50 mm×200 mm×50 mm的混凝土塊體，塊體表面規(guī)則，單元內(nèi)表面為曲面。塊體內(nèi)布置了3根鋼筋，其起點(diǎn)和終點(diǎn)的整體坐標(biāo)分別為1號(hào)鋼筋(25,0,25)?(35,200,25)、2號(hào)鋼筋(10,0,12)?(40,200,45)和3號(hào)鋼筋(10,0,10)?(12,200,16)。表1給出上述3種迭代初值方案的計(jì)算結(jié)果。

圖3 鋼筋分布圖

表1 迭代次數(shù)

由表1可知，3種方案都能求到真實(shí)解，采用方案1時(shí)，1、3號(hào)鋼筋的迭代次數(shù)明顯少于其余兩種方案，2號(hào)鋼筋的迭代次數(shù)與兩種方案較為接近。方案1的迭代次數(shù)整體偏小，即迭代初值點(diǎn)取等參單元表面的中心點(diǎn)時(shí)收斂速度最好。

當(dāng)表面交點(diǎn)為兩個(gè)時(shí)，給出了2種迭代初值選擇方案：

方案1. 以中心點(diǎn)作為迭代初值點(diǎn)，求得第一個(gè)交點(diǎn)后，將局部坐標(biāo)取負(fù)號(hào)再次進(jìn)行迭代。

方案2. 分別以4個(gè)角點(diǎn)作為迭代初值進(jìn)行循環(huán)，求得兩個(gè)交點(diǎn)后跳出循環(huán)。

下面采用給定結(jié)點(diǎn)坐標(biāo)構(gòu)造單元表面，考慮兩個(gè)交點(diǎn)位置的所有可能情況來討論上述兩種方案的可行性。依據(jù)等參單元表面中心點(diǎn)將表面劃分成4個(gè)象限，第1個(gè)交點(diǎn)在第3象限，第2個(gè)交點(diǎn)分別位于第1、2、3及4象限，如圖4所示，即兩個(gè)交點(diǎn)分別位于3-1、3-2、3-3、3-4象限。表2給出了其計(jì)算結(jié)果。

圖4 象限劃分圖

表2 迭代總次數(shù)

對(duì)迭代解還需驗(yàn)證是否為真實(shí)的交點(diǎn)，這時(shí)先用迭代解局部坐標(biāo)值,n+1和,n+1是否均在區(qū)間[–1,1]內(nèi)來判斷，若成立，再用該交點(diǎn)至鋼筋起點(diǎn)P和終點(diǎn)P距離之和是否等于鋼筋線總長(zhǎng)度來判斷，滿足則為真實(shí)的交點(diǎn)。

1.2 分塊解析法

1.2.1 等參單元表面方程建立

先對(duì)有中間結(jié)點(diǎn)的混凝土單元利用式(4)補(bǔ)全可能缺少的中間結(jié)點(diǎn)，使所有混凝土單元為8或20結(jié)點(diǎn)等參單元。

建立等參單元表面方程時(shí)，先區(qū)分是否為平面。在單元表面任取3個(gè)結(jié)點(diǎn)建立平面方程，若其余結(jié)點(diǎn)都滿足該平面方程則該表面為平面，該平面方程就是該面的表面方程。

若等參單元表面為曲面，則分8和20結(jié)點(diǎn)等參單元兩種情況分別建立表面方程。

(1) 8結(jié)點(diǎn)等參單元。對(duì)于8結(jié)點(diǎn)等參單元，其表面的邊界線均為直線，當(dāng)該曲面不經(jīng)過坐標(biāo)原點(diǎn)時(shí)曲面方程通式為

其中，為9個(gè)待定系數(shù)。

利用已知的4個(gè)角點(diǎn)由式(4)插值得到4個(gè)中間結(jié)點(diǎn)和1個(gè)表面中點(diǎn)的結(jié)點(diǎn)坐標(biāo)，組成9個(gè)坐標(biāo)點(diǎn)后代入式(13)，就形成包含9個(gè)待定系數(shù)的非齊次線性方程組，采用高斯消去法求解就可得這9個(gè)待定系數(shù)。若該曲面經(jīng)過原點(diǎn)，可進(jìn)行平移變換求得交點(diǎn)。

(2) 20結(jié)點(diǎn)等參單元。對(duì)于20結(jié)點(diǎn)等參單元，一般情況下采用與8結(jié)點(diǎn)單元相同的方法能求得單元表面方程，但對(duì)下列兩種情況需特別處理。

①方程通式問題。當(dāng)式(13)中的某些項(xiàng)不存在時(shí)，如母線平行于軸的柱面方程與軸無關(guān)，式(13)中待定系數(shù)為0，按式(13)求解待定系數(shù)時(shí)，求解過程中會(huì)出現(xiàn)系數(shù)矩陣主元為0的情況，這時(shí)需將主元為0所在的行與列以及對(duì)應(yīng)行的常數(shù)項(xiàng)取為0，同時(shí)將該主元取1，以保證方程順利求解。

②計(jì)算精度問題。由于20結(jié)點(diǎn)單元棱邊存在曲線，當(dāng)求得的交點(diǎn)局部坐標(biāo)不能滿足計(jì)算精度的要求時(shí)，可利用中間結(jié)點(diǎn)將該面分成4個(gè)子塊(圖5)，分別建立表面方程求解，即分塊求交。當(dāng)計(jì)算精度還不滿足時(shí)，可再次分塊求交，如圖6所示。

圖5 分塊結(jié)點(diǎn)圖

圖6 逐次分塊示意圖

1.2.2 鋼筋直線方程與交點(diǎn)真?zhèn)闻袆e

鋼筋線直線方程可利用其起點(diǎn)坐標(biāo)(x,y,z)和終點(diǎn)坐標(biāo)(x,y,z)寫成以為變量的參數(shù)方程

聯(lián)立式(13)~(14)即可求得交點(diǎn)，交點(diǎn)可能是一個(gè)或兩個(gè)，本文依次采用下述3個(gè)步驟確定真實(shí)交點(diǎn)：

步驟1.先根據(jù)交點(diǎn)與單元邊界的位置關(guān)系進(jìn)行判別，交點(diǎn)坐標(biāo)值必定在單元結(jié)點(diǎn)坐標(biāo)最大值和最小值之間。

步驟2. 再依據(jù)交點(diǎn)的局部坐標(biāo)值進(jìn)行判別。求交點(diǎn)整體坐標(biāo)(,,)所對(duì)應(yīng)的局部坐標(biāo)(,,)可采用下列迭代格式

若該交點(diǎn)為實(shí)際交點(diǎn)，則局部坐標(biāo)中必有一項(xiàng)絕對(duì)值為1，同時(shí)其余兩項(xiàng)在[–1,1]區(qū)間內(nèi)。若求得的局部坐標(biāo)有兩項(xiàng)在[–1,1]區(qū)間內(nèi)，另一項(xiàng)絕對(duì)值接近1但不滿足計(jì)算精度要求，則分塊求解。

步驟3. 最后用該交點(diǎn)至鋼筋起點(diǎn)P和終點(diǎn)P距離之和是否等于鋼筋線總長(zhǎng)度來判斷，若成立則該迭代解為真實(shí)的交點(diǎn)。

2 鋼筋線與混凝土單元求交過程

鋼筋線與混凝土單元表面具體求交過程：

(1) 對(duì)所有混凝土單元循環(huán)，利用鋼筋線起點(diǎn)和終點(diǎn)與單元邊界的位置關(guān)系確定出鋼筋起點(diǎn)P和終點(diǎn)P所在的單元號(hào)(該點(diǎn)坐標(biāo)值必定在所在單元的結(jié)點(diǎn)坐標(biāo)最大值和最小值之間)，記為N和N。

(2) 利用上節(jié)所述方法對(duì)單元N所有表面求交，找到鋼筋線與單元N的交點(diǎn)，并記下交點(diǎn)所在表面S，然后對(duì)除N外所有單元循環(huán)，找到與S共面的相鄰單元N，N即鋼筋要進(jìn)入的下一單元。

(3) 令N=N，對(duì)第2步循環(huán)，但需注意這時(shí)對(duì)S面不再求交，直至N=N。

3 算例及處理方法的比較

以一個(gè)承臺(tái)算例來驗(yàn)證本文方法，并對(duì)兩種方法進(jìn)行比較。承臺(tái)截面尺寸為800 mm×800 mm，高700 mm，承臺(tái)內(nèi)布置直徑12 mm、間距200 mm的分布鋼筋。樁沿交界位置向下截取長(zhǎng)800 mm的一段，樁身直徑為400 mm。采用空間埋置組合式單元模型，剖分混凝土單元時(shí)不必考慮鋼筋位置，剖分后共有2 176個(gè)單元，如圖7所示。由于模型中承臺(tái)網(wǎng)格剖分時(shí)受限于樁體網(wǎng)格，沿樁豎直方向的承臺(tái)網(wǎng)格與樁體網(wǎng)格保持一致，此時(shí)會(huì)出現(xiàn)鋼筋穿過混凝土單元曲面表面的情況。

通過運(yùn)算發(fā)現(xiàn)，兩種處理方法得到的結(jié)果相同，但由于牛頓迭代法需分別取表面4個(gè)角點(diǎn)作為迭代初值，當(dāng)僅有一個(gè)交點(diǎn)時(shí)多迭代了3次；相比而言，分塊解析法計(jì)算簡(jiǎn)便。由于分布鋼筋較多，本文只取靠近樁體兩個(gè)方向上有代表性的4根分布鋼筋，并將模型倒置，最終結(jié)果如圖8所示。

圖7 結(jié)構(gòu)有限元網(wǎng)格模型

圖8 鋼筋布置圖

4 結(jié)束語

(1) 對(duì)表面為任意形狀的單元，提出了兩種鋼筋線與混凝土單元表面求交方法，并給出其計(jì)算過程。利用這兩種方法能有效得到直鋼筋穿過混凝土單元表面的鋼筋結(jié)點(diǎn)坐標(biāo)信息，增強(qiáng)了空間埋置組合單元模型的適用性。

(2) 兩種求交方法均能得到正確解。牛頓迭代法需分別取表面4個(gè)角點(diǎn)作為迭代初值進(jìn)行迭代，當(dāng)僅有一個(gè)交點(diǎn)時(shí)多迭代了3次；分塊解析法利用方程直接求出表面交點(diǎn)，計(jì)算簡(jiǎn)便。從適用性而言，分塊解析法要優(yōu)于牛頓迭代法。

[1] 江見鯨, 陸新征. 混凝土結(jié)構(gòu)有限元分析[M]. 2版. 北京: 清華大學(xué)出版社, 2013: 129-133.

[2] NGO D, SCORDELIS A C. Finite element analysis of reinforced concrete beams [J]. Journal of American Concrete Institute, 1967, 64(3): 152-163.

[3] GOODMAN R E, TALYOR R L, BREKKE T L. Closure on a model for the mechanics of jointed rock [J]. Journal of the Soil Mechanics and Foundations Division, 1970, 94(3): 637-660.

[4] SUIDAN M, SCHNOBRICH W C. Finite element analysis of reinforced concrete [J]. Journal of the Structural Division, 1973, 99: 2109-2122.

[5] CHANG T Y, TANIGUCHI H, CHEN W F. Nonlinear finite element analysis of reinforced concrete panels [J]. Journal of Structural Engineering, 1987, 113(1): 122-140.

[6] ZIENKIEWICZ O C, OWEN D R J, PHILLIPS D V, et al. Finite element methods in the analysis of reactor vessels [J]. Nuclear Engineering & Design, 1972, 20(2): 507-541.

[7] 汪基偉, 張雄文, 林新志. 考慮粘結(jié)滑移的平面組合式單元模型研究與應(yīng)用[J]. 工程力學(xué), 2008, 25(1): 97-102.

[8] 巫昌海, 汪基偉. 混凝土三維鋼筋埋置組合式有限單元模型及其網(wǎng)格自動(dòng)生成[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2000, 12(10): 761-764.

[9] 關(guān)振群, 宋超, 顧元憲, 等. 有限元網(wǎng)格生成方法研究的新進(jìn)展[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2003, 15(1): 1-14.

[10] 古成中, 吳新躍. 有限元網(wǎng)格劃分及發(fā)展趨勢(shì)[J]. 計(jì)算機(jī)科學(xué)與探索, 2008, 2(3): 248-259.

[11] 王明強(qiáng), 朱永梅, 劉文欣. 有限元網(wǎng)格劃分方法應(yīng)用研究[J]. 機(jī)械設(shè)計(jì)與制造, 2004(1): 22-24.

Surface Intersection Algorithm Based on Spatial Embedded Finite Element Model

WANG Jiwei1, FU Yu1, WANG Yajiang2

(1. College of Civil and Transportation Engineering, Hohai University, Nanjing Jiangsu 210098, China; 2. Quanzhou Audit Bureau, Quanzhou Fujian 362000, China)

In the spatial embedded finite element model, the bar element can be embedded at any position in the concrete element. The element model is convenient because the mesh of the concrete element is not limited by the reinforcement. In the application of the element model, it is necessary to obtain the coordinates where the both ends of the steel element lie at the surface of the concrete element. Therefore, it is the premise of application of the model to solve the intersection of the reinforcement line and the concrete element surface. The existing method is only suit for the concrete element with plane surface. In this paper, the Newton iterative method and the divided block analytic method are proposed which can solve the intersection of the reinforcement line and the surface of concrete element with any shape. The method in the paper can enhance the applicability of the spatial embedded finite element model. The correctness of the two methods is verified by an example. In terms of applicability, the divided block analytic method is superior to the Newton iterative method.

reinforced concrete finite element method; element model; pre-processing of finite element method; curved surface

TP 391；TU 375

10.11996/JG.j.2095-302X.2018020221

A

2095-302X(2018)02-0221-07

2017-08-04；

2017-09-14

國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目(2016YFC0402002)

汪基偉(1962–)，男，浙江安吉人，教授，博士，博士生導(dǎo)師。主要研究方向?yàn)殇摻罨炷痢-mail：wjw2903918@126.com