田志昌,趙波,趙根田

內(nèi)蒙古科技大學(xué)建筑與土木工程學(xué)院,內(nèi)蒙古包頭014010

基于六面體單元程序計(jì)算地基沉降

田志昌,趙波*,趙根田

內(nèi)蒙古科技大學(xué)建筑與土木工程學(xué)院,內(nèi)蒙古包頭014010

本文根據(jù)線彈性力學(xué)的有限單元法，選用六面體單元，借助Fortran 90具有模塊化、封裝機(jī)制、自定義等編程特性，開(kāi)發(fā)線彈性階段的六面體單元計(jì)算程序，為計(jì)算地基沉降提出一種可行、有效的方法。并利用該方法能自行修訂所需參數(shù)，以達(dá)到計(jì)算精度要求。最后，通過(guò)與ANSYS軟件計(jì)算模擬地基沉降的結(jié)果比較與分析，驗(yàn)證其可靠性和適用性，并對(duì)比不同節(jié)點(diǎn)六面體單元的精確度。

Fortran 90;六面體單元;地基沉降

地基沉降量的控制是防止產(chǎn)生不均勻沉降而出現(xiàn)裂縫、局部?jī)A斜等工程安全隱患的關(guān)鍵[1]，因此地基沉降的計(jì)算和預(yù)測(cè)成為建筑地基設(shè)計(jì)的核心問(wèn)題。目前其計(jì)算方法主要分為兩種：一種是把分層總和法做為基礎(chǔ)的實(shí)用計(jì)算方法；另一種是以建立材料本構(gòu)關(guān)系為主的有限元等數(shù)值方法[2]。而我國(guó)在《建筑地基基礎(chǔ)設(shè)計(jì)規(guī)范》[3]中，建議使用對(duì)分層總和法修訂后的規(guī)范法。雖然該方法的適用性較好，但即便是同一規(guī)范中，由于Ψs，壓縮層厚度確定方法的不同[4,5]，都會(huì)導(dǎo)致不同的計(jì)算值。且規(guī)范法以一維變形的假定也影響了計(jì)算精度。因?yàn)橥翆拥膸缀螛?gòu)造和受力情況以三維性為主，其一維壓縮只是一種特殊情況，如在內(nèi)蒙古包頭地區(qū)高層建筑厚筏板基礎(chǔ)沉降計(jì)算中，規(guī)范法計(jì)算結(jié)果與實(shí)測(cè)沉降值存在較大誤差，個(gè)別工程竟高達(dá)允許誤差的5倍，故該方法并不具有參考價(jià)值[6]。

數(shù)值理論方法較為完善，可以考慮非線性、彈塑性等，伴隨著計(jì)算機(jī)技術(shù)的發(fā)展，產(chǎn)生了許多具備計(jì)算地基沉降的有限元軟件，如AbAQUS，ANSYS但土性參數(shù)的確定比較復(fù)雜和困難，且選取參數(shù)的種類(lèi)也不同[7]。因此，我國(guó)許多的軟件，如YJK，GS等，都還沒(méi)有嵌入三維實(shí)體單元計(jì)算模塊。所以，本文試圖先在線彈性階段結(jié)合六面體等參數(shù)單元理論和Fortran 90[8]編程語(yǔ)言，開(kāi)發(fā)應(yīng)用程序來(lái)計(jì)算地基沉降的方法。為進(jìn)一步探索非線性階段的編程尊定基礎(chǔ)。首先，三維六面體單元具有塊體形狀的特性，適合模擬巖土這種固體，故用它來(lái)計(jì)算地基沉降和承載力分析具有一定的優(yōu)越性。再次，這種小型計(jì)算程序可以使具備Fortran 90編程知識(shí)的使用者，針對(duì)地基實(shí)際情況，自行在程序中添加或減少影響地基沉降的控制參數(shù)和條件，來(lái)達(dá)到所需的準(zhǔn)確度。

1 六面體單元程序的理論基礎(chǔ)

1.1 等參數(shù)單元

母單元和子單元是所有等參數(shù)單元的兩種單元形式。在母單元內(nèi)主要進(jìn)行計(jì)算工作，這是由于它具備計(jì)算精度高和適用性好的特點(diǎn)。而子單元充分反應(yīng)了來(lái)自結(jié)構(gòu)的實(shí)際特征。這兩種單元的相互映射則是通過(guò)形函數(shù)進(jìn)行坐標(biāo)變換來(lái)達(dá)到的，因此就存在著存在著兩種不同的坐標(biāo)系，即全局坐標(biāo)（X，Y，Z）和局部坐標(biāo)(ξ，η，ζ)。子單元的直角坐標(biāo)系始終是全局坐標(biāo)；而局部坐標(biāo)不僅扮演著母單元的直角坐標(biāo)系，也扮演著子單元的曲線坐標(biāo)系兩種角色。并且曲線坐標(biāo)系只適用于單個(gè)獨(dú)立的子單元，而直角坐標(biāo)被整個(gè)結(jié)構(gòu)中的所有子單元中共同使用。

所以等參數(shù)單元對(duì)于曲線邊界和曲面邊界適應(yīng)性更強(qiáng)，模擬結(jié)構(gòu)形狀更加準(zhǔn)確，也因其位移模式階數(shù)較高而更好地反映結(jié)構(gòu)中復(fù)雜的應(yīng)力分布，并且還可以靈活的增減結(jié)點(diǎn)來(lái)構(gòu)造各種過(guò)度單元。

1.2 單元形態(tài)的布置

由于等參數(shù)單元進(jìn)行了坐標(biāo)變換，所以其單元內(nèi)應(yīng)變的變化規(guī)律，不僅取決于形函數(shù)的階次，并且還與單元形狀及各邊節(jié)點(diǎn)位置相關(guān)，因此需要確定單元形態(tài)。同時(shí)單元節(jié)點(diǎn)編號(hào)，相應(yīng)的形函數(shù)位置和順序須與其規(guī)定的單元形態(tài)（母單元）相符合，否則會(huì)導(dǎo)致雅克比矩陣行列式值為負(fù)（單元形態(tài)出現(xiàn)內(nèi)凹現(xiàn)象），進(jìn)而使計(jì)算結(jié)果錯(cuò)誤或程序不能運(yùn)行。

六面體等參數(shù)單元包括線性和二次兩種單元，其單元形態(tài)布置分別采用圖1和圖2的單元（母單元）形態(tài)，全局坐標(biāo)系采用右手坐標(biāo)系，Z軸向上。

圖1 8節(jié)點(diǎn)空間等參元Fig.1 The isoparametric elements in 8 nodes

圖2 20節(jié)點(diǎn)空間等參元Fig.2 The isoparametric elements in 20 nodes

1.3 本構(gòu)矩陣

彈性D矩陣也稱(chēng)本構(gòu)矩陣，對(duì)于不同的材料性質(zhì)，如極端各向異性體、正交各向異性體等，都可化為6*6彈性常數(shù)矩陣，實(shí)質(zhì)是材料不同其中包含的獨(dú)立彈性常數(shù)不同。由于只考慮線彈性階段，因此它的本構(gòu)關(guān)系符合廣義的胡克定律，即彈性體材料的彈性模量E和μ泊松比決定了它的獨(dú)立彈性常數(shù)，如式(1)。

1.4 應(yīng)變矩陣和雅克比矩陣

通過(guò)把單元位移函數(shù)式代入到空間應(yīng)變問(wèn)題的公式中，即可得到：

式中B就為應(yīng)變矩陣。根據(jù)微分法則并對(duì)其進(jìn)行集合，得到：

式中i為節(jié)點(diǎn)數(shù)，(x，y，z)為全局坐標(biāo)，形函數(shù)Ni(ζ，η，ξ)則是用局部坐標(biāo)表達(dá)的，

通過(guò)公式(4)和(5)的分解與回代，可得到應(yīng)變B矩陣。從積分學(xué)可知，雅克比矩陣｜J｜＞0是兩個(gè)坐標(biāo)之間一一變換的充要條件，所以等參變換也必須服從此條件。如果｜J｜=0，則[J]-1不存在，產(chǎn)生導(dǎo)數(shù)和微元轉(zhuǎn)換也就不存在，變換不成立。其中雅克比率的概念就是是給定單元偏離理想單元（正常單元）的一個(gè)度量，其取值范圍是-1.0到1.0。1.0代表理想單元形狀。測(cè)量方法本質(zhì)就是通過(guò)映射的方式把局部坐標(biāo)下的理想單元形狀轉(zhuǎn)化為全局坐標(biāo)下的實(shí)際形狀。

1.5 單元?jiǎng)偠扰c數(shù)值積分

整個(gè)結(jié)構(gòu)矩陣計(jì)算中的一個(gè)核心環(huán)節(jié)就是單元?jiǎng)偠染仃嚨男纬伞卧獎(jiǎng)偠染仃囋诤艽蟪潭壬蠜Q定著計(jì)算結(jié)果的精度，它隨著單元幾何形狀、節(jié)點(diǎn)數(shù)目及節(jié)點(diǎn)自由度選用的不同而發(fā)生變化。單元?jiǎng)偠染仃囍械膹椥跃仃噭t取決于材料的物性參數(shù)。由此可見(jiàn)，只要單元的幾何形狀及其材料物理特性確定，單元的剛度矩陣就可確定。

式中ξi，ηj，ζk為積分點(diǎn)，Wi，Wj，Wk為權(quán)系數(shù)，nip為積分次數(shù)。

數(shù)值積分有兩類(lèi)積分方法，包括所取積分點(diǎn)是等間距和不等間距的。此處采用積分點(diǎn)不等間距的高斯積分法，因?yàn)樗梢栽诜e分點(diǎn)數(shù)目確定的情況下通過(guò)選擇位置更加合理的積分點(diǎn)來(lái)減少積分點(diǎn)的索取數(shù)量并達(dá)到更高的精度，從而節(jié)省機(jī)器時(shí)間和便于程序的實(shí)現(xiàn)。雖然各維數(shù)上的積分點(diǎn)數(shù)目由各個(gè)自變量在被積函數(shù)中可能出現(xiàn)的最高次數(shù)分別決定但為使用方便，常常統(tǒng)一為最高值。因此，對(duì)于線性六面體單元采用2*2*2規(guī)則計(jì)算單元?jiǎng)偠染仃?，二次單元?jiǎng)t采用3*3*3規(guī)則來(lái)計(jì)算。

1.6 結(jié)構(gòu)剛度矩陣組裝、存貯和分解

結(jié)構(gòu)剛度矩陣是由單元?jiǎng)偠染仃嚱M裝而成，包含了每個(gè)單元對(duì)結(jié)構(gòu)產(chǎn)生的貢獻(xiàn)作用，這種貢獻(xiàn)作用的疊加則是通過(guò)單元定位向量來(lái)完成的。單元定位向量不僅使每個(gè)單元的各個(gè)未知量在結(jié)構(gòu)未知量中確定相應(yīng)的位置還考慮了每個(gè)單元的邊界條件，很方便地解決了程序編寫(xiě)時(shí)計(jì)算機(jī)自動(dòng)化中所需的信息。對(duì)于結(jié)構(gòu)剛度矩陣的存貯，采用節(jié)省存貯單元的一維變帶寬存貯，這種存貯方法不僅可以提高計(jì)算效率還對(duì)減小舍入誤差有著明顯的效果。其存貯方式有按行存貯和按列存貯兩種[9]，存貯序號(hào)分別按公式（8）和（9）計(jì)算。

上式中字母I，J分別代表元素的行列號(hào)，數(shù)組KD存放結(jié)構(gòu)剛度矩陣中的主對(duì)角線上的元素。

本程序中結(jié)構(gòu)剛度矩陣主對(duì)角元素選用按行存貯。因?yàn)榉纸饨Y(jié)構(gòu)剛度矩陣時(shí)采用消元分解法，而消元分解法求解時(shí)是按行進(jìn)行肖元的，比按列存貯更有效。最后通過(guò)結(jié)構(gòu)靜力學(xué)有限元方程（10）即可解出節(jié)點(diǎn)位移量。

其中{F}為節(jié)點(diǎn)載荷向量；[K]為結(jié)構(gòu)剛度矩陣；syggg00為節(jié)點(diǎn)位移量。

圖3為主程序的結(jié)構(gòu)圖，揭示了程序編寫(xiě)的方案和各個(gè)子程的作用,并且給出了計(jì)算線性六面體單元?jiǎng)偠染仃嚨氖纠绦颉?/p>

圖3 主程序流程圖Fig.3 The process of main program

計(jì)算線性六面體單元?jiǎng)偠染仃嚦绦?

2 計(jì)算模擬地基沉降

2.1 六面體單元程序計(jì)算地基沉降

本計(jì)算程序所適用的單元材料屬性由其本構(gòu)矩陣D來(lái)確定，采用有限元法中的小片檢驗(yàn)的原理[10]來(lái)檢查、驗(yàn)證。模擬地基的模型建立和有限單元的劃分借助ANSYS軟件。為此采用簡(jiǎn)化實(shí)用的地基模型進(jìn)行沉降計(jì)算，其模型分析步驟如下：

(1)地基表面受6種不同的工況荷載，分別為：120 KN/m2、130 KN/m2、150 KN/m2、170 KN/m2、180 KN/m2、200 KN/m2；(2)不計(jì)地基重力作用；(3)邊界條件：①所有垂直邊界水平位移為零；②底部網(wǎng)格所有節(jié)點(diǎn)垂直、水平位移為零；③所有垂直邊界不排水。(4)雜填土下為基巖；(5)表1提供模擬所需計(jì)算參數(shù)；(6)圖4為地基沉降計(jì)算簡(jiǎn)圖。

圖4 地基沉降計(jì)算簡(jiǎn)圖Fig.4 the diagram calculating foundation settlement

表1 地基模擬的計(jì)算參數(shù)Table 1 Parameters of foundation simulation

依據(jù)表1和模擬分析步驟的內(nèi)容，表2和表3分別表示模擬地基在6種不同工況荷載下，利用六面體單元程序中不同節(jié)點(diǎn)單元計(jì)算的最大沉降值。

2.2 ANSYS計(jì)算地基沉降

ANSYS軟件因具有多場(chǎng)耦合分析功能，不僅能使巖土力學(xué)性能的模擬更加準(zhǔn)確，還可以考慮非線性的本構(gòu)關(guān)系及分期施工過(guò)程，使實(shí)際情況在模擬中得到較好地反映。其中Drucker-Prager屈服準(zhǔn)則可為巖土介質(zhì)的模擬計(jì)算提供精確的結(jié)果。采用DP本構(gòu)模型,采用SOLID45(線性六面體單元)和SOLID95(二次六面體單元)兩種單元形式進(jìn)行劃分，計(jì)算2.1中模擬地基在6種不同荷載下的沉降值，計(jì)算結(jié)果為表4。

表4 ANSYS計(jì)算地基沉降值Table 4 The value of foundation settlement calculating byANSYS

3 結(jié)論

(1)本文用ANSYS、六面體單元程序?qū)δM地基的沉降進(jìn)行了計(jì)算。從圖5和圖6中得出，相對(duì)于ANSYS計(jì)算模擬地基的沉降值，使用六面體單元程序的線性六面體單元計(jì)算模擬地基的沉降值與ANSYS計(jì)算值之間的誤差約為3%。而使用其中的二次六面體單元計(jì)算模擬地基的沉降值與ANSYS計(jì)算值之間的誤差約2%。說(shuō)明應(yīng)用程序單元形態(tài)的布置與節(jié)點(diǎn)形函數(shù)的布置一致，程序編制正確，調(diào)試、運(yùn)行流暢，以及使用的理論公式和數(shù)學(xué)方法基本正確。同時(shí)也說(shuō)明六面體單元程序中使用二次單元計(jì)算精度要高于使用線性單元，其單元適應(yīng)能力更強(qiáng)，實(shí)際運(yùn)用中能用其較少或適量的單元?jiǎng)澐值玫骄_的計(jì)算結(jié)果。

圖5 ANSYS和線性六面體單元程序計(jì)算結(jié)果對(duì)比圖Fig.5 The comparison betweenANSYS and the program of linear hexahedral element

(2)通過(guò)與ANSYS計(jì)算結(jié)果的比較，六面體單元應(yīng)用程序的整體性、融合性、各個(gè)子程序之間的協(xié)調(diào)性、計(jì)算結(jié)果的可靠性和適用性良好，不僅為程序非線性階段的編制奠定了基礎(chǔ)，也可視為計(jì)算地基沉降的一種方法，并與其它方法相結(jié)合來(lái)有效地防范和應(yīng)對(duì)施工沉降中出現(xiàn)的一些特殊情況，也可為設(shè)計(jì)、施工人員或科研工作者解決或預(yù)測(cè)地基沉降提供有效的參考和對(duì)比。

[1]周德泉,楊志華,鄧超.復(fù)合地基沉降分析與計(jì)算方法研究[J].公路與汽運(yùn),2015,1(170):85-88

[2]楊光華.地基沉降計(jì)算的新方法[J].巖石力學(xué)與工程學(xué)報(bào),2008,27(4):679-686

[3]中華人民共和國(guó)建設(shè)部.GB50007-2002建筑地基基礎(chǔ)設(shè)計(jì)規(guī)范[S].北京:中國(guó)建筑工業(yè)出版社,2002

[4]張欽喜,劉鴻哲,樊紹峰.地基沉降計(jì)算方法的研究與改進(jìn)[J].北京工業(yè)大學(xué)學(xué)報(bào),2009(1):78-83

[5]劉全林,魏煥衛(wèi).地基沉降計(jì)算中壓縮層厚度確定方法的比較[J].巖土工程技術(shù),2001(4):208-209

[6]張寧.類(lèi)比法預(yù)測(cè)包頭市民營(yíng)經(jīng)濟(jì)發(fā)展中心公寓基礎(chǔ)沉降[D].包頭:內(nèi)蒙古科技大學(xué),2013

[7]曹強(qiáng).基于施工沉降觀測(cè)值預(yù)測(cè)深層回填土地基最終沉降[D].包頭:內(nèi)蒙古科技大學(xué)報(bào),2013

[8]顏慧軍,王友海.程序設(shè)計(jì)語(yǔ)言FORTRAN90研究與應(yīng)用[J].哈爾濱建筑大學(xué)學(xué)報(bào),2000,33(2):105-109

[9]趙超燮.結(jié)構(gòu)矩陣分析原理[M].北京:人民教育出版社,1982:148-149

[10]朱伯芳.有限單元法原理與應(yīng)用[M].第三版.北京:中國(guó)水利水電出版社,2009:154-157

Foundation Settlement Calculated by the Program of Hexahedral Element

TIAN Zhi-chang,ZHAO Bo*,ZHAO Gen-tian
College of Architecture and Civil Engineering/Inner Mongolia University of Science and Technology,Baotou 014010,China

According to the linear elastic mechanics for finite element,this paper chose the most common hexahedral element to resort to Fortran 90 with the characteristics of modularization,encapsulation mechanism and custom to develop the program of hexahedral element in the linear elastic stage so as to propose a feasible and effective method to calculate the foundation settlement.This method was used to revise the required parameters in order to meet the precision of the calculation.Finally,the results between this method and ANSYS simulation were compared and analyzed to verify the reliability and applicability and contrast the accuracy of hexahedral element in different nodes.

Fortran 90;hexahedral element;foundation settlement

U446.1

A

1000-2324(2016)06-0841-05

2016-06-28

2016-08-30

國(guó)家自然科學(xué)基金(51268042)

田志昌(1962-),男,內(nèi)蒙古包頭人,博士,教授,主要從事有限元軟件開(kāi)發(fā)和防災(zāi)減災(zāi).E-mail:1435227030@qq.com

*通訊作者:Author for correspondence.E-mail:308012774@qq.com