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

汪基偉，付 宇，王亞獎

?

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

汪基偉1，付 宇1，王亞獎2

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

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

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

現(xiàn)有鋼筋混凝土結(jié)構(gòu)有限元模型[1]中，分離式單元模型采用塊體單元模擬混凝土，桿單元模擬鋼筋[2-3]，桿單元必須在混凝土單元邊界上，結(jié)構(gòu)配筋復(fù)雜時網(wǎng)格剖分困難，有時甚至難以實現(xiàn)。整體式單元模型[4]將鋼筋彌散于整個單元中，利用單元參數(shù)(如配筋率)求得考慮鋼筋作用后的單元折算彈性模量，再以折算彈性模量求單元剛度矩陣，以考慮鋼筋對單元剛度矩陣的貢獻(xiàn)，該模型無法考慮鋼筋的具體位置，也不能求出鋼筋應(yīng)力分布。埋置式單元模型[5-6]將鋼筋埋置于混凝土單元內(nèi)部，網(wǎng)格剖分方便且能求得鋼筋應(yīng)力，單元剛度為混凝土單元剛度與鋼筋單元剛度之和，其中混凝土單元剛度矩陣仍采用等參單元的剛度矩陣，而鋼筋單元剛度矩陣則利用混凝土與鋼筋應(yīng)變相同的原則得到。文獻(xiàn)[7]和[8]分別推導(dǎo)出了二維和三維埋置式單元模型，鋼筋單元剛度矩陣s和應(yīng)力矩陣s表達(dá)式為

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

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

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

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

1.1 牛頓迭代法

1.1.1 迭代法方程的建立

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

圖1 結(jié)點編號圖

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

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

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

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

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

同時對中間結(jié)點相鄰的兩個角結(jié)點的形函數(shù)進行修改，各減該中間結(jié)點形函數(shù)的一半。

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

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

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

當(dāng)有中間結(jié)點5時，因中間結(jié)點5對角結(jié)點1和2有影響，式(8)中的系數(shù)11、13、21、23不再為0，同時式(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é)點時，按照上述思路就可寫出相應(yīng)的表達(dá)式，即

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

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

1.1.2 迭代法方程的求解

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

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

圖2 迭代初值點分布圖

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

圖3 鋼筋分布圖

表1 迭代次數(shù)

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

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

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

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

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

圖4 象限劃分圖

表2 迭代總次數(shù)

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

1.2 分塊解析法

1.2.1 等參單元表面方程建立

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

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

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

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

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

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

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

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

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

圖5 分塊結(jié)點圖

圖6 逐次分塊示意圖

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

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

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

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

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

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

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

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

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

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

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

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

3 算例及處理方法的比較

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

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

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

圖8 鋼筋布置圖

4 結(jié)束語

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

(2) 兩種求交方法均能得到正確解。牛頓迭代法需分別取表面4個角點作為迭代初值進行迭代，當(dāng)僅有一個交點時多迭代了3次；分塊解析法利用方程直接求出表面交點，計算簡便。從適用性而言，分塊解析法要優(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)格自動生成[J]. 計算機輔助設(shè)計與圖形學(xué)學(xué)報, 2000, 12(10): 761-764.

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

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

[11] 王明強, 朱永梅, 劉文欣. 有限元網(wǎng)格劃分方法應(yīng)用研究[J]. 機械設(shè)計與制造, 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

國家重點研發(fā)計劃項目(2016YFC0402002)

汪基偉(1962–)，男，浙江安吉人，教授，博士，博士生導(dǎo)師。主要研究方向為鋼筋混凝土。E-mail：wjw2903918@126.com