田志昌,趙波,趙根田

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

基于六面體單元程序計算地基沉降

田志昌,趙波*,趙根田

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

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

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

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

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

1 六面體單元程序的理論基礎

1.1 等參數(shù)單元

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

所以等參數(shù)單元對于曲線邊界和曲面邊界適應性更強，模擬結構形狀更加準確，也因其位移模式階數(shù)較高而更好地反映結構中復雜的應力分布，并且還可以靈活的增減結點來構造各種過度單元。

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

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

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

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

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

1.3 本構矩陣

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

1.4 應變矩陣和雅克比矩陣

通過把單元位移函數(shù)式代入到空間應變問題的公式中，即可得到：

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

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

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

1.5 單元剛度與數(shù)值積分

整個結構矩陣計算中的一個核心環(huán)節(jié)就是單元剛度矩陣的形成。單元剛度矩陣在很大程度上決定著計算結果的精度，它隨著單元幾何形狀、節(jié)點數(shù)目及節(jié)點自由度選用的不同而發(fā)生變化。單元剛度矩陣中的彈性矩陣則取決于材料的物性參數(shù)。由此可見，只要單元的幾何形狀及其材料物理特性確定，單元的剛度矩陣就可確定。

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

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

1.6 結構剛度矩陣組裝、存貯和分解

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

上式中字母I，J分別代表元素的行列號，數(shù)組KD存放結構剛度矩陣中的主對角線上的元素。

本程序中結構剛度矩陣主對角元素選用按行存貯。因為分解結構剛度矩陣時采用消元分解法，而消元分解法求解時是按行進行肖元的，比按列存貯更有效。最后通過結構靜力學有限元方程（10）即可解出節(jié)點位移量。

其中{F}為節(jié)點載荷向量；[K]為結構剛度矩陣；s8co0ac為節(jié)點位移量。

圖3為主程序的結構圖，揭示了程序編寫的方案和各個子程的作用,并且給出了計算線性六面體單元剛度矩陣的示例程序。

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

計算線性六面體單元剛度矩陣程序:

2 計算模擬地基沉降

2.1 六面體單元程序計算地基沉降

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

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

圖4 地基沉降計算簡圖Fig.4 the diagram calculating foundation settlement

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

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

2.2 ANSYS計算地基沉降

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

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

3 結論

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

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

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

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

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

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

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

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

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

[7]曹強.基于施工沉降觀測值預測深層回填土地基最終沉降[D].包頭:內(nèi)蒙古科技大學報,2013

[8]顏慧軍,王友海.程序設計語言FORTRAN90研究與應用[J].哈爾濱建筑大學學報,2000,33(2):105-109

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

[10]朱伯芳.有限單元法原理與應用[M].第三版.北京:中國水利水電出版社,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

國家自然科學基金(51268042)

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

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