付延玲,金瑋澤,陳興賢,2,談金忠
1.河海大學(xué)地球科學(xué)與工程學(xué)院,南京 210098
2.江蘇省地質(zhì)礦產(chǎn)勘查局,南京 210018
3.江蘇南京地質(zhì)工程勘察院,南京 210041
地面沉降是一種由多種因素引起的地面高程緩慢降低的地質(zhì)現(xiàn)象,其中過量開采地下水是引發(fā)地面沉降問題的主要原因。但近年來,隨著城市建設(shè)進(jìn)程的加快,高層、超高層建筑的不斷增多,建筑荷載引發(fā)的地面沉降問題也日益突出。研究者們[1-4]認(rèn)為,在上海,高群體建筑引起的地面沉降占總沉降量的30%左右,且沉降速度與建筑面積增長(zhǎng)速度之間存在線性相關(guān)關(guān)系。
高層建筑荷載一方面直接導(dǎo)致土體有效應(yīng)力增加,發(fā)生固結(jié)壓縮,引發(fā)地面沉降;另一方面在地下水埋深較淺的地區(qū),高層建筑荷載的作用也會(huì)對(duì)地下水滲流場(chǎng)產(chǎn)生影響,導(dǎo)致滲流場(chǎng)重新調(diào)整,滲流場(chǎng)的變化同樣會(huì)引起應(yīng)力場(chǎng)的變化,引發(fā)地面沉降或隆起變形。以往學(xué)者的研究更多集中在由抽水作用引起地下水滲流場(chǎng)變化而引發(fā)的地面沉降變形方面[5-7],對(duì)由高層建筑荷載引發(fā)的土體變形特征,尤其是高層建筑荷載短時(shí)間作用下引起滲流場(chǎng)調(diào)整變化,進(jìn)而形成的土體變形的研究尚不多見。
為了準(zhǔn)確模擬預(yù)測(cè)由高層建筑荷載引發(fā)的地面沉降與隆起變形問題,筆者以比奧固結(jié)理論為基礎(chǔ),綜合考慮了含水層和隔水層的流變特性,以及水力學(xué)和土力學(xué)參數(shù)隨滲流場(chǎng)和應(yīng)力場(chǎng)發(fā)生變化的動(dòng)態(tài)變化問題,建立了高層建筑荷載引發(fā)地面沉降與隆起變形的三維有限元數(shù)值模型,對(duì)高層建筑荷載引起的地面沉降與隆起變形進(jìn)行了模擬預(yù)測(cè)分析,探討了變形過程中土體水力參數(shù)與力學(xué)參數(shù)的變化特征。
飽和土體(假定其土骨架變形為線彈性、微小變形、滲流符合達(dá)西定律、水不可壓縮或微壓縮)的三維比奧固結(jié)方程[8]如下:
式中:G為剪切模量;ν為泊松比;wx、wy、wz分別為x、y、z方向上的位移分量;u為孔隙水壓力;Kx、Ky、Kz分別為x、y、z方向上的滲透系數(shù);γ為土的重度;γw為水的重度;▽為拉普拉斯算子,▽2=為時(shí)間。
土的本構(gòu)關(guān)系是土的力學(xué)特性即應(yīng)力-應(yīng)變-強(qiáng)度-時(shí)間等關(guān)系的數(shù)學(xué)表達(dá)式。對(duì)于考慮流變特性的土體來說,其變形特征主要表現(xiàn)為變形的時(shí)間與應(yīng)力水平有關(guān),所顯示的是具有彈性、塑性和黏滯性的黏彈塑性體,若將此類土體的總應(yīng)變?cè)隽縟ε分為彈塑性應(yīng)變?cè)隽縟εep、黏彈性應(yīng)變?cè)隽縟εve和黏塑性應(yīng)變?cè)隽縟εvp,則具有流變特性的土體中任意點(diǎn)在任意時(shí)刻的應(yīng)變?cè)隽浚?,9]為
式(3)中各部分的應(yīng)變?cè)隽靠梢杂上率龇椒ù_定。
1)彈塑性應(yīng)變?cè)隽?/p>
由土體的彈塑性本構(gòu)模型可得
2)黏彈性應(yīng)變?cè)隽?/p>
由kelvin流變模型E1ε+Keε=σ,可得應(yīng)力不變時(shí)復(fù)雜應(yīng)力狀態(tài)下的黏彈性應(yīng)變?cè)隽浚?/p>
式中:ηe=E1/Ke,E1為kelvin體黏彈性模量,Ke為黏滯系數(shù);A為應(yīng)力矩陣,
E為彈性模量。
3)黏塑性應(yīng)變?cè)隽?/p>
利用黏塑性法確定黏塑性應(yīng)變?cè)隽?,在該種方法中允許材料在有限“期間”內(nèi)超越破壞準(zhǔn)則(以破壞準(zhǔn)則函數(shù)F的值大于0來表示)。在討論土體的黏塑性應(yīng)變而非塑性應(yīng)變時(shí),應(yīng)變的變化率與超越量有關(guān),即有以下關(guān)系式:
式中:Qs為塑性勢(shì)函數(shù)對(duì)于摩爾-庫(kù)倫材料來說,
其中:φ為內(nèi)摩擦角;c為黏聚力;
β=為第三偏應(yīng)力不變量,
σ與τ分別為主應(yīng)力與剪應(yīng)力。
如果將黏塑性應(yīng)變率與一個(gè)偽時(shí)間步相乘,就可以得到累加到下一個(gè)時(shí)間步的黏塑性應(yīng)變?cè)隽?,于是?/p>
數(shù)值計(jì)算與絕對(duì)穩(wěn)定的時(shí)間步與假定的破壞準(zhǔn)則有關(guān)。
對(duì)于摩爾-庫(kù)倫材料有
塑性勢(shì)函數(shù)對(duì)應(yīng)力的偏導(dǎo)數(shù)可以表示為
所以,應(yīng)用摩爾-庫(kù)倫準(zhǔn)則的土體黏塑性應(yīng)變?cè)隽靠梢员硎緸?/p>
將式(4)、(5)、(10)代入式(3),可求得彈塑-黏彈-黏塑性體的應(yīng)變?cè)隽繛?/p>
式(11)即為土體的黏彈塑性本構(gòu)方程。
利用伽遼金加權(quán)余量法離散方程,考慮到土體的非線性特性,取Δt時(shí)間內(nèi)的位移增量來代替位移,將式(1)、(2)離散成增量形式[10]:
式中:Δδ為結(jié)點(diǎn)位移增量;Δu為結(jié)點(diǎn)孔隙壓力增量;ˉK為固體剛度矩陣;K為滲透流量矩陣;K′為應(yīng)力-滲流耦合項(xiàng)矩陣;ΔQ為流量增量矩陣;B為自由面的積分矩陣;R為等效節(jié)點(diǎn)荷載,當(dāng)存在高層建筑荷載時(shí),R為包括高層建筑荷載引起的附加應(yīng)力值;Rt為t時(shí)刻已經(jīng)發(fā)生的位移所平衡了的那部分荷載。
因?yàn)闈B流取決于孔隙壓力全量的分布,而不是取決于時(shí)間內(nèi)孔隙壓力增量,所以孔壓要用全量的形式表示。記時(shí)刻tn和tn+1時(shí)單元節(jié)點(diǎn)i的孔壓全量分別為ui(n)和ui(n+1),且 Δui=ui(n+1)-ui(n),則式(12)可變換為
式(13)即為三維比奧固結(jié)有限元方程。
1.4.1 孔隙度與滲透系數(shù)的非線性
流固耦合問題實(shí)際上是孔隙應(yīng)力的消散引起土體骨架的變形,導(dǎo)致滲透系數(shù)變化,從而影響土體的滲透性,宏觀上表現(xiàn)為土體的固結(jié)變形。在比奧固結(jié)的假定條件下,根據(jù)孔隙度的相關(guān)定義和滲流力學(xué)Kozeny-Carman方程推得孔隙度n和滲透系數(shù)K的動(dòng)態(tài)表達(dá)式[11-12]:
式中:n0為初始孔隙度;K0為初始滲透系數(shù);εv為體應(yīng)變
1.4.2 土體參數(shù)的非線性
采用鄧肯-張非線性模型,將土體的本構(gòu)關(guān)系推廣到非線性,則本構(gòu)關(guān)系{Δσ}=D{Δε}中矩陣D中的彈性常數(shù)E、ν不再視為常量,而是隨著應(yīng)力狀態(tài)改變而改變,其切線彈性模量Et和切線泊松比νt的表達(dá)式如下[13]:
式中:Rf為破壞比;σ1為第一主應(yīng)力;σ3為第三主應(yīng)力;ω為彈性模量與固結(jié)壓力曲線的斜率;logα,G為土體常規(guī)三軸壓縮實(shí)驗(yàn)結(jié)果所繪曲線截距;I=0.04,D=3為土體實(shí)驗(yàn)參數(shù);p為大氣壓強(qiáng)。
1.5.1 初始條件
1)地應(yīng)力初始條件
采用土體的自重應(yīng)力估算土體的初始應(yīng)力:
式中:σx、σz為土體的初始水平向和垂向應(yīng)力;z為計(jì)算點(diǎn)深度;K0為靜止側(cè)壓力系數(shù),K0=。為有效內(nèi)摩擦角
2)位移初始條件
3)孔隙水壓力初始條件
式中:u0(x,y,z)為研究區(qū)域內(nèi)已知初始孔隙水壓力。
1.5.2 邊界條件
1)孔隙水壓力邊界條件Γ1
式中:us為水頭邊界Γ1上的已知孔隙水壓力。
2)流量邊界條件Γ2
式中:qL為邊界Γ2上的已知單位面積流量。
3)自由面邊界條件Γ3
式中:μ為土體給水度;θ為自由面外法線方向與垂線的交角;q為通過自由面邊界Γ3的單位面積流量;Z為自由面所在的高程。
4)位移邊界條件Γ4
比奧固結(jié)有限元方程結(jié)合定解條件和水力學(xué)參數(shù)及土力學(xué)參數(shù)動(dòng)態(tài)變化模型即可運(yùn)用Fortran語言編制相應(yīng)的有限元程序進(jìn)行求解[5,14]。
為了研究高層建筑荷載影響下的地面沉降情況,排除其他影響因素干擾,建立一個(gè)簡(jiǎn)單的均質(zhì)含水層模型。模型長(zhǎng)×寬×高為500m×500m×100 m,建筑物按每層自重1.2t/m2、承重每層0.3t/m2計(jì)算,地基尺寸的長(zhǎng)×寬為70m×50m。用八節(jié)點(diǎn)六面體單元離散化模型,在平面上剖分為2 500個(gè)矩形網(wǎng)格單元,垂向上從上往下概化為:潛水含水層、第1承壓含水層,第2承壓含水層、第3承壓含水層及各含水層之間的黏性土弱含水層,共7層。每層土體劃分為一個(gè)參數(shù)分區(qū),共7個(gè)參數(shù)分區(qū)。其中:潛水含水層巖性以粉砂、亞砂土為主,底板埋深20~30m,靜水位埋深0.9~1.2m;第1承壓含水層以粉砂為主,底板埋深50~52m,靜水位埋深3.2~3.8m;第2承壓含水層以亞砂土、細(xì)中砂為主,底板埋深79~83m,靜水位埋深7.6~8.1m;第3承壓含水層以細(xì)中砂、中粗砂為主,底板埋深96~100m,靜水位埋深11.5~13.4m。模型四周均概化為第一類已知水頭邊界,底部概化為隔水邊界。建筑物位置及模型分層如圖1所示。模型各參數(shù)分區(qū)物理力學(xué)性質(zhì)參數(shù)如表1所示。其中:K0x、K0y為初始水平向滲透系數(shù);K0z為初始垂向滲透系數(shù);ν0為初始泊松比;E0為初始彈性模量;Sy為儲(chǔ)水率。
圖1 模型示意圖Fig.1 Model sketch
模型建筑物高度按10層計(jì)算。模型計(jì)算周期為1a,劃分為15個(gè)應(yīng)力期。其中:第1到第4應(yīng)力期分別為5,10,20,31d;第5到15應(yīng)力期每個(gè)應(yīng)力期時(shí)間分別為第2月到第12月單月時(shí)間。假定初始情況下的各參數(shù)分區(qū)孔隙度大小均相同,且初始孔隙度n0=0.36。選取建筑物中心點(diǎn)所在剖面(y=260m)模型運(yùn)行15個(gè)應(yīng)力期時(shí)的地面沉降量,如圖2所示,其中地面變形量正值表示隆起量,負(fù)值表示沉降量。
圖2 不同應(yīng)力期末的地面沉降曲線圖Fig.2 Land subsidence curve of difference stress period
土體的固結(jié)壓縮變形較為復(fù)雜,受到土體本身性質(zhì)、邊界條件、排水條件以及上部荷載方式等因素影響。由圖2可以看出,由高層建筑荷載引起的地面沉降以建筑物中心點(diǎn)為中心呈現(xiàn)漏斗狀:到第4個(gè)應(yīng)力期末,即31d時(shí),最大沉降達(dá)到-3.8mm;在第10個(gè)應(yīng)力期末,即第7個(gè)月末時(shí),最大沉降量為-8.5mm。由圖2b可以看出:模型運(yùn)行前9個(gè)應(yīng)力期時(shí),高層建筑荷載引起的地面沉降在x=200~300m和x=400~500m區(qū)間段均存在隆起現(xiàn)象;在第10個(gè)應(yīng)力期末,即第7個(gè)月末時(shí),隆起現(xiàn)象消失。而由圖2a可知:在高層建筑荷載作用5,10,20d時(shí)隆起量依次增大,最大值分別為0.9mm,1.0mm,1.5mm;31d時(shí)隆起量為1.3mm??梢钥闯?,20d累計(jì)隆起量相對(duì)于31d累計(jì)隆起量較大。這是因?yàn)橛捎诟邔咏ㄖ奢d主要影響淺部地層,在模型中,淺部土體垂向滲透系數(shù)較小,所以淺部地層的孔隙水在受到高層建筑荷載作用時(shí),更多地以橫向流動(dòng)為主;而高層建筑荷載作用時(shí)間較短,其周圍土體中的孔隙水急劇增多,導(dǎo)致超靜孔隙水壓力在初期急劇增大,有效應(yīng)力相應(yīng)減小,土體骨架膨脹變形,必然使得建筑周圍土體發(fā)生回彈,地表出現(xiàn)隆起;且隨著應(yīng)力期的增加,建筑周圍土體的超靜孔隙水壓力緩慢消散,地面隆起緩慢消失。
隨著土體的固結(jié)壓縮與回彈變形,建筑物周圍及下部土體力學(xué)及水力學(xué)參數(shù)均會(huì)伴隨著土體的變形而發(fā)生相應(yīng)變化。選取模型第1161號(hào)單元和第1167號(hào)單元進(jìn)行孔隙度、滲透系數(shù)、彈性模量及泊松比分析。其中:第1161號(hào)單元位于建筑物中心底部;第1167號(hào)單元中心點(diǎn)坐標(biāo)為(430,260),在剖面y=200m上,距離建筑物40m,且兩單元均位于模型第一層,單元中心線位于y=260m剖面上。兩單元孔隙度及滲透系數(shù)隨應(yīng)力期的變化如圖3所示。
由于1161號(hào)單元與1167號(hào)單元所處區(qū)域地面沉降變化情況不同,1161號(hào)單元位于模型第1層地面沉降中心點(diǎn)位置,而1167號(hào)單元所處區(qū)域在前9個(gè)應(yīng)力期均出現(xiàn)隆起現(xiàn)象,所以兩單元的孔隙度變化趨勢(shì)是不同的。由圖3a可以看出:1167號(hào)單元孔隙度在20d增大到最大值,隨著應(yīng)力期的增加,孔隙度緩慢減??;而1161號(hào)單元在整個(gè)應(yīng)力期時(shí)間段內(nèi),孔隙度一直處于減小的趨勢(shì)。造成這種現(xiàn)象的原因可能是:在建筑荷載作用初期,建筑物下部土體的孔隙水短時(shí)間擴(kuò)散到周圍淺部地層,使得孔隙水大量聚集,超靜孔隙水壓力急劇增大,土體有效應(yīng)力相應(yīng)減小,淺部土體發(fā)生膨脹,使得建筑物周圍土體孔隙度在20d達(dá)到最大值;但隨后超靜孔隙水壓力逐漸消散,使得孔隙度緩慢減小,在第10個(gè)應(yīng)力期,即第7個(gè)月時(shí)恢復(fù)到初始孔隙度大小。建筑物中心點(diǎn)下部土體在受到建筑荷載作用受到擠壓后,土體孔隙水向四周擴(kuò)散,使得孔隙水壓力一直減小,有效應(yīng)力一直增加,所以孔隙度隨著應(yīng)力期的變化始終減小。
滲透系數(shù)隨著孔隙度的變化而發(fā)生變化??紫抖仍龃螅馏w出水能力增強(qiáng),導(dǎo)致滲透系數(shù)變大,而孔隙度減小,使土體出水能力變?nèi)?,從而?dǎo)致滲透系數(shù)減小。由圖3也可以看出,1161號(hào)單元和1167號(hào)單元的滲透系數(shù)與孔隙度變化趨勢(shì)相同。
表1 地層參數(shù)Table 1 Stratum parameters
圖3 1161號(hào)單元和1167號(hào)單元孔隙度、水平向滲透系數(shù)和垂向滲透系數(shù)隨應(yīng)力期變化曲線圖Fig.3 Variation curve among porosity and horizontal hydraulic conductivity and vertical hydraulic conductivity and stress periods of 1161element and 1167element
伴隨著土體的回彈及壓縮,彈性模量、泊松比必然發(fā)生變化。由圖4可以看出,1167號(hào)單元彈性模量、泊松比的變化與孔隙度的變化相對(duì)應(yīng),在20d左右土體彈性模量達(dá)到最小值,泊松比達(dá)到最大值,隨著應(yīng)力期增加彈性模量緩慢增大,泊松比緩慢減小。而1161號(hào)單元彈性模量及泊松比的變化均較為平緩,彈性模量不斷增大,泊松比不斷減小。通過對(duì)土體力學(xué)及水力學(xué)參數(shù)的變化趨勢(shì)分析可以發(fā)現(xiàn),在前9個(gè)應(yīng)力期內(nèi)各參數(shù)的變化相比后幾個(gè)應(yīng)力期要快,這是因?yàn)殡S著應(yīng)力期的增長(zhǎng),土體固結(jié)變形量逐漸減小,導(dǎo)致土體力學(xué)及水力參數(shù)的變化相對(duì)于前期而言越來越不明顯。
1)以比奧固結(jié)理論為基礎(chǔ),結(jié)合土體的非線性流變理論,將土體的本構(gòu)關(guān)系推廣到黏彈塑性,同時(shí)考慮土體水力學(xué)參數(shù)和土力學(xué)參數(shù)隨有效應(yīng)力的動(dòng)態(tài)變化問題,建立起的地下水滲流與土體變形三維全耦合模型更加符合實(shí)際,進(jìn)一步提高了模型計(jì)算的精度。
2)由高層建筑荷載引起的地面沉降呈現(xiàn)漏斗狀,以建筑物中心點(diǎn)為漏斗中心。建筑荷載作用初期會(huì)引起周圍土體隆起,建筑物周圍淺部土層土體孔隙度、滲透系數(shù)及泊松比也隨著增加,到達(dá)峰值后,隨著時(shí)間的增加又緩慢減小;而彈性模量逐漸減小,到達(dá)最小值后,隨著時(shí)間的增加緩慢增大。建筑荷載下部土層,孔隙度、滲透系數(shù)及泊松比隨時(shí)間的增加均緩慢減小,而彈性模量隨時(shí)間的增加緩慢增大。
(References):
[1]龔士良.上海城市建設(shè)對(duì)地面沉降的影響[J].中國(guó)地質(zhì)災(zāi)害與防治學(xué)報(bào),1998,9(2):108-111.Gong Shiliang.Urban Construction’s Impact on Land Subsidence in Shanghai[J].The Chinese Journal of Geological Hazard and Control,1998,9(2):108-111.
[2]許燁霜,馬磊,沈水龍.上海市城市化進(jìn)程引起的地面沉降因素分析[J].巖土力學(xué),2011,32(增刊1):578-582.Xu Yeshuang,Ma Lei,Shen Shuilong.Influential Factors on Development of Land Subsidence with Process of Urbanization in Shanghai[J].Rock and Soil Mechanics,2011,32(Sup.1):578-582.
[3]嚴(yán)學(xué)新,沈國(guó)平.上海城區(qū)建筑密度與地面沉降關(guān)系分析[J].水文地質(zhì)工程地質(zhì),2002,29(6):21-25.Yan Xuexin,Shen Guoping.Relationship Between Building Density and Land Subsidence in Shanghai Urban Zone[J].Hydrogeology and Engineering Geology,2002,29(6):21-25.
[4]唐益群,宋壽鵬,陳斌,等.不同建筑容積率下密集建筑群區(qū)地面沉降規(guī)律研究[J].巖石力學(xué)與工程學(xué)報(bào),2010,29(增刊1):3425-3431.Tang Yiqun,Song Shoupeng,Chen Bin,et al.Study of Land Subsidence Rule of Dense Buildings Under Different Floor Area Ratios[J].Chinese Journal of Rock Mechanics and Engineering,2010,29(Sup.1):3425-3431.
[5]陳興賢,駱祖江,安曉宇,等.深基坑降水三維變參數(shù)非穩(wěn)定滲流與地面沉降耦合模型[J].吉林大學(xué)學(xué)報(bào):地球科學(xué)版,2013,43(5):1572-1578.Chen Xingxian,Luo Zujiang,An Xiaoyu,et al.Coupling Model of Groundwatrer Three Dimensional Variable Parametrics Non-Steady Seepage and Land-Subsidence[J].Journal of Jilin University:Earth Science Edition,2013,43(5);1572-1578.
[6]Luo Zujiang,Zeng Feng.Finite Element Numerical Simulation of Land Subsidence and Groundwater Exploitation Based on Visco-Elastic-Plastic Biot’s Consolidation Theory[J].Journal of Hydrodynamics,2011,23(5):615-624.
[7]錢家歡,殷宗澤.土工原理與計(jì)算[M].北京:水力水電出版社,1996.Qian Jiahuan,Yin Zongze.Principle and Canculation of Geotechnics[M].Beijing:China Waterpower Press,1996.
[8]陳曉平,白世偉.軟黏土地基黏彈塑性比奧固結(jié)的數(shù)值分析[J].巖土工程學(xué)報(bào),2001,23(4):481-484.Chen Xiaoping,Bai Shiwei.The Numerical Analysis of Visco-Elastic-Plastic Biot’s Consolidation for Soft Clay Foundation[J].Chinese Journal of Geotechnical Engineering,2001,23(4):481-484.
[9]李醫(yī)民,周鳳燕.一類三維偏微分方程邊值問題的解法[J].江蘇大學(xué)學(xué)報(bào):自然科學(xué)版,2004,25(4):328-331.Li Yimin,Zhou Fengyan.Solution to a Class of Boundary Value Problem of Three Dimensional Partial Differential Equation[J].Journal of Jiangsu University:Natural Science Edition,2004,25(4):328-331.
[10]冉啟全,李士倫.流固耦合油藏?cái)?shù)值模擬中物性參數(shù)動(dòng)態(tài)模型研究[J].石油勘探與開發(fā),1997,24(3):61-65.Ran Qiquan,Li Shlun.Study on Dynamic Models of Reservoir Parameters in the Coupled Simulation of Multiphase Flow and Reservoir Deformation[J].Petroleum Exploration and Development,1997,24(3):61-65.
[11]田杰,劉先貴,尚根華.基于流固耦合理論的套損力學(xué)機(jī)理分析[J].水動(dòng)力學(xué)研究與進(jìn)展:A 輯,2005,20(2):221-225.Tian Jie,Liu Xiangui,Shang Genhua.Casing Damage Mechanism Based on Theory of Fluid-Solid Coupling Flow Through Underground Rock[J].Journal of Hydrodynamics:Series A,2005,20(2):221-225.
[12]羅剛,張建民.鄧肯-張模型和沈珠江雙屈服面模型的改進(jìn)[J].巖土力學(xué),2004,25(6):887-890.Luo Gang,Zhang Jianmin.Improvement of Duncan-Chang Nonlinear Model and Shen Zhujiang’s Elastoplastic Model for Granular Soils[J].Rock and Soil Mechanics,2004,25(6):887-890.
[13]彭國(guó)倫.Fortran 95程序設(shè)計(jì)[M].北京:中國(guó)電力出版社,2005.Peng Guolun.Fortran 95[M].Beijing:China Electric Power Press,2005.
[14]Smith I M,Griffiths.D V.有限元方法編程[M].3版.王菘,周堅(jiān)鑫,王來,等,譯.北京:電子工業(yè)出版社,2003.Smith I M,Griffiths D V.Programming the Finite Element Method[M].3rd ed.Translated by Wang Song,Zhou Jianxin, Wang Lai,et al.Beijing:Electronic Industry Press,2003.