基于ALE方法的盾構對接施工鋼環(huán)頂伸數(shù)值模擬及穩(wěn)定性分析

王歡歡1,2,楊麗麗3,楊勛1,2,金先龍1,2

(1.上海交通大學機械系統(tǒng)與振動國家重點實驗室，上海200240；2. 上海交通大學機械與動力工程學院，上海200240; 3. 上海應用技術學院，上海201400)

摘要：基于任意拉格朗日-歐拉有限元方法(ALE)，建立了盾構對接施工鋼環(huán)頂伸的三維數(shù)值模型，解決了土體變形流動導致網(wǎng)格畸變的問題，實現(xiàn)了鋼環(huán)-土體動態(tài)耦合過程。分析了鋼環(huán)頂伸過程中，鋼環(huán)前側阻力變化、鋼環(huán)應力應變狀況。研究了不同的頂伸速度下，對鋼環(huán)頂伸的影響。在充分考慮接觸及重力條件下，分析了鋼環(huán)軸向及側向的穩(wěn)定性狀況。研究結果可為盾構對接施工方案設計提供參考依據(jù)。

關鍵詞：任意拉格朗日-歐拉方法；盾構對接；鋼環(huán)頂伸；穩(wěn)定性

中圖分類號:U455.4; TP391.9

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2015.19.014

Abstract:Based on arbitrary Lagrangian-Eulerian(ALE) method, a three-dimensional finite element model for steel tube extension in shield machine butting was established. The problem of mesh deformation distortion caused by soil flow deformation was solved, the steel tube-soil dynamic coupled process was realized. Resistance change and stress-strain conditions of steel tube during the extension was analyzed. Simulation of steel tube extending under different velocities was implemented. Axial stability and lateral stability of steel tube were analyzed considering contact and gravity status. The results provided a reference for the design of shield butting construction.

基金項目：國家自然科學基金重大研究計劃(91215302)；國家自然科學基金(51408207，51478366) 國家自然科學基金項目(91216106);國家重大基礎研究項目資助

收稿日期：2014-05-05修改稿收到日期：2014-11-19 2014-05-19修改稿收到日期：2014-09-25

Simulation and stability analysis for steel tube extending in shield butting construction based on the ALE method

WANGHuan-huan1,2,YANGLi-li3,YANGXun1,2,JINXian-long1,2(1. State Key Laboratory of Mechanical System and Vibration，Shanghai Jiaotong University., Shanghai 200240, China;2. College of Mechanical Engineering, Shanghai Jiaotong University., Shanghai 200240, China3. Shanghai Institute of Technology, Shanghai 201400, China)

Key words:arbitrary Lagrangian-Eulerian method; shield butting; steel tube extension; stability

近幾十年，盾構施工技術在國內外得到了廣泛應用。對于一些超長超大直徑的大型隧道，單條隧道連續(xù)施工會造成周期過長、施工故障概率高等問題。因此，盾構地中對接施工成為解決該問題的一種有效技術方法[1-2]。盾構地中對接法即兩臺盾構機從兩側相向掘進至結合地點，在地中進行對接以完成整條隧道的盾構施工，以其安全性、經(jīng)濟性、省時高效等優(yōu)點被廣泛應用。

盾構地中對接方式主要有土木式對接法和機械式對接法：土木式對接法是通過將對接地點地層作加固處理，達到止水和預防地層失穩(wěn)的效果后，完成盾構拆卸并進行隧道襯砌；機械式對接法是通過對盾構進行特殊設計，而使兩臺盾構直接進行對接的方法。在工程策劃與施工過程中，應結合工程實際情況，根據(jù)土層狀況、工期造價、精度要求選擇對接方法。相比較機械對接法，土木對接法以其經(jīng)濟性更為常用，也是本文對接施工所采用的方法。

對于直徑接近和超過15m的超大直徑隧道，盾構對接施工仍是一項具有挑戰(zhàn)性的世界級難題，國際上也缺乏足夠的經(jīng)驗積累和深入研究。為了確保超大直徑盾構對接施工的安全和質量，有必要對對接施工中關鍵技術問題需要進行深入研究。其中，鋼環(huán)頂推過程和鋼環(huán)卸載過程是盾構對接施工中的最關鍵兩步驟。因此，鋼環(huán)頂推過程和穩(wěn)定性分析是本項目的研究重點。

鋼環(huán)在土體中的頂推過程，屬于典型的結構-土體耦合問題。對于此類問題，國內外已有一些相關研究。Bauer[3]分析了輪胎在土體上的轉動過程，通過將土體的摩擦力及粘滯阻力作為外載荷加載在輪胎上，并沒有考慮土體變形等因素的影響；Kushwaha[4]分析了農業(yè)機具與土體間的相互作用，并考慮土體的小變形及材料非線性，雖然可以求得機具在土體摩擦力及阻力作用下的動態(tài)響應，但是無法模擬土體在農機作用下產生的較大變形，因而會產生較大的偏差；Nezami[5]采用DEM方法模擬了挖掘機鏟斗工作的過程，克服了土體大變形帶來的網(wǎng)格畸變，并同試驗數(shù)據(jù)進行了比對；Asaf[6]指出利用DEM方法求解結構耦合問題時，其參數(shù)的選取缺乏穩(wěn)定性，導致求解結果失真；Fredj[7]采用ALE方法，模擬了管線鋪設的過程，充分考慮管道及土體的大變形問題，有效解決了管道-土體間的耦合作用問題，并利用試驗進行了驗證。

當薄壁圓筒型構件承受外壓時，將在筒壁上產生軸向和環(huán)向壓力。如果壓縮應力超過材料的屈服點或強度極限時，和內壓圓筒一樣，將發(fā)生強度破壞。然而往往在筒壁壓應力未達到材料的屈服極限點時，筒體突然失去原有形狀發(fā)生屈曲失穩(wěn)。在本項目中，鋼環(huán)直徑較大(15m)，而厚度僅為0.1m，故有可能在未達到材料屈服極限前而發(fā)生失穩(wěn)。對于存在接觸面的結構屈曲問題進行分析時，以往研究對于接觸區(qū)域往往進行單獨處理。Herzl[8]采用選代算法首先算出接觸區(qū)域的虛擬力，然后將這種虛擬力轉化成假想彈簧的剛度系數(shù)，并對原始剛度矩陣進行修正，以處理屈曲模態(tài)中的接觸問題。在對鋼環(huán)進行穩(wěn)定性分析時，受摩擦力及土體粘滯阻力的影響，如果直接將其簡化為虛擬力，將導致計算結果存在一定的偏差。

本文基于ALE方法，提出了盾構對接過程中鋼環(huán)頂伸的數(shù)值模擬方法，研究鋼環(huán)在頂伸過程中的應力應變狀況，進而在充分考慮接觸及重力條件下，分析鋼環(huán)軸向及側向的穩(wěn)定性狀況，為盾構對接施工提供設計參考依據(jù)。

1原理與方法

1.1ALE算法

鋼環(huán)在頂推過程中，受土體的阻滯作用下，產生一定的位移與變形，同時土體在鋼環(huán)的作用下發(fā)生大變形及塑性流動，需對兩者界面間的耦合問題進行有效處理。

為實現(xiàn)對鋼環(huán)頂推過程和卸載過程的數(shù)值模擬，應建立鋼環(huán)與土體的耦合動力學有限元模型。其中，鋼環(huán)可以簡化為線彈性材料，并且變形較小，但必須考慮土體的彈塑性大變形和塑性流動。為了既能描述鋼環(huán)的彈性小變形，又能描述土體彈塑性大變形和塑性流動，本項目擬采用ALE (Arbitrary Lagrangian-Eulerian)理論來建立鋼環(huán)與土體的耦合動力學有限元模型。

ALE方法綜合了 Lagrangian 方法與 Euler 方法的優(yōu)點，在材料域與空間域外引入了參考域，通過在參考域網(wǎng)格上的求解，既解決了Lagrangian 方法中材料的嚴重扭曲，又可解決 Eulerian 方法中移動邊界引起的復雜性問題。

ALE描述中包含三個坐標系：Lagrangian坐標系、Eulerian坐標系以及參考坐標系，三者間的轉換關系為：

(1)

式中，Xt表示Lagrangian坐標，i表示參考網(wǎng)格，xt表示Eulerian坐標。vt和ut分別為材料和網(wǎng)格速度，wt為對流速度。

圖1為三種坐標系之間的關系示意圖。

圖1　Lagrangian，Eulerian和ALE域之間的映射 Fig.1 Mappings between Lagrangian, Eulerian, ALE descriptions

在不考慮溫度變化的情況下，基于ALE描述下的土體連續(xù)性方程、動量方程為：

(2)

(3)

基于Lagrangian描述的鋼環(huán)結構質量方程、動量方程為：

ρs(X,t)J(X,t)=ρ0(X,t)J0(X,t)

(4)

(5)

式中，ρs為鋼環(huán)密度；J為空間和材料坐標間的Jacobian行列式；u為固體結構位移；X為Lagrangian坐標；bi為體力；σij為應力張量。

鋼環(huán)與土體間的相互作用通過罰函數(shù)耦合方式來定義。罰函數(shù)方式耦合通過跟蹤流體與結構之間的相對位移來施加節(jié)點力，節(jié)點力的大小與位移成比例關系。

1.2顯式計算方法

采用 ALE 建模理論建立的鋼環(huán)-土體耦合動力學有限元模型是一個典型的非線性動力學問題?？梢圆捎蔑@式計算方法來對這類復雜非線性動力學問題進行數(shù)值求解。

對于某一個時刻，非線性動力學系統(tǒng)的運動方程為：

(6)

采用顯式中心差分法對上式進行數(shù)值求解，假定0，t1，t2，…，tn時刻節(jié)點位移、速度與加速度均為已知，現(xiàn)求解tn+1(t+Δt)時刻的結構動力響應。在中心差分法中，加速度和速度可以用位移表示，即

(7)

將式(7)代入式(6)中，即可得中心差分法的遞推公式

(8)

式中

分別為有效質量矩陣與有效載荷向量。求解式(8)，即可獲得t+Δt時刻節(jié)點位移向量Xt+Δt。在整個時域范圍內，可由上述遞推公式求得各個離散時間點處節(jié)點的位移、速度和加速度。

顯式計算方法由于采用集中質量矩陣，運動方程的求解是非耦合的，不需要組集成總體剛度矩陣，并且采用中心單點積分，因此大大節(jié)省存儲空間和求解機時。

2工程概況及模型建立

2.1問題概述

盾構機對接位置及周邊土體幾何簡圖見圖2。該對接方案的具體思路為：首先1號盾構掘進至指定位置，進行隧道加固、泥水倉固化，盾構部分結構拆除，此時2號盾構同時在正常推進，在距離對接位置50m處，進行相對位置探測(多次)，并進行修正掘進(多次)，并到達對接地點。

在到達指定位置后，2號盾構進行隧道加固、刀盤倉內固化、盾構部分結構拆除，準備工作就緒后進行鋼環(huán)套體頂伸，并達到指定的頂伸長度，見圖3。鋼環(huán)頂伸到位后，進行土體凍結達到一定強度。隨后進行凍結施工，清理刀盤內泥土，對接部密封并加固、盾構機解體與清理等一系列工序。在其中鋼環(huán)頂伸作為較為關鍵的步驟，將直接關系到盾構對接能否順利完成，所以需對其頂伸過程中的應力應變狀況及穩(wěn)定性進行深入研究。

圖2　盾構機及土體幾何簡圖 Fig.2 Position of shield and soil

圖3　鋼環(huán)頂伸示意圖 Fig.3 Diagram of steel tube promote

2.2盾構頂推有限元模型

據(jù)提供的工程地質勘探資料，以鋼環(huán)頂伸區(qū)域為核心區(qū)域，采用六面體實體單元，建立了土體的分層有限元模型，見圖4。

鋼環(huán)直徑為15m，而厚度僅為0.1m，制造加工存在較大的難度。在保證足夠的強度下，為降低制造難度，初步設計中鋼環(huán)前端沒有采用刃角構造，在隨后的研究中可以考慮不同結構的影響。鋼環(huán)所在位置的局部視圖見圖5，鋼環(huán)頂伸區(qū)域的網(wǎng)格進行了細化，對密封圈等結構進行了簡化處理。

圖4　有限元整體模型 Fig.4 Integrated finite element model

圖5　鋼環(huán)處有限元局部模型 Fig.5 Partial finite element model of steel tube

2.3材料參數(shù)與邊界條件

目前廣為接受的是采用彈塑性材料模型來模擬土體材料的非線性力學行為，通常情況下常采用相對簡單的材料模型，使用時其屈服面遵循德魯克-普拉格(Drucker-Prager)屈服準則或莫爾-庫侖(Mohr-Coulomb)屈服準則。本仿真數(shù)值模型采用的是軟件中提供的Drucker-Prager材料本構，各土層分布及材料參數(shù)見圖2及表1。

表1　各土層材料參數(shù)

在有限元模型中，土體底部為全約束，側面法向約束，上部為自由端，并整體施加全局重力。

2.4初始地應力平衡

在處理初始地應力時，主要考慮重力場對土體的作用，且暫不考慮盾構與土體間的作用關系。重力初始平衡首先對土體有限元模型加載重力，然后將平衡狀態(tài)時刻的土體應力場導出，作為初始應力加載。仿真模型采用此種方法可以避免初始有限元模型在重力作用下的大變形，較好的模擬了實際變形情況。土體初始沉降及土體初始應力見圖6、圖7。

圖6　土體初始沉降 Fig.6 Initial settlement of soil

圖7　土體初始應力 Fig.7 Initial stress of soil

2.5屈曲分析有限元模型

在進行屈曲分析時，無需對鋼環(huán)后側襯砌等結構進行建模，以減小建模規(guī)模、節(jié)省計算時間。根據(jù)地質資料，以鋼環(huán)為核心區(qū)域，采用六面體實體單元，建立鋼環(huán)軸向屈曲分析的有限元模型，見圖8。軸向屈曲分析有限元模型局部視圖見圖9，并充分考慮鋼環(huán)與兩側土體間的接觸關系。

圖8　屈曲分析模型 Fig.8 Finite element model of buckling analysis

圖9　軸向屈曲分析模型局部視圖 Fig.9 Partial finite element model of axial buckling analysis

在進行側向屈曲分析時，鋼環(huán)內部土體被挖空、盾構機內部主要支撐結構被拆卸，故側向屈曲分析有限元模型鋼環(huán)內部與軸向分析時略有區(qū)別，軸向屈曲分析有限元模型局部視圖見圖10。

圖10　側向屈曲分析模型局部視圖 Fig.10 Partial finite element model of lateral buckling analysis

3結果分析與討論

3.1鋼環(huán)前側阻力分析

在獲得了土體初始沉降及土體初始應力之后，將平衡狀態(tài)時刻的土體應力場導出，作為初始應力加載；再加入盾構、鋼環(huán)等內部結構，同時施加重力，計算得到盾構以及內部結構的應力形態(tài)；然后利用非線性顯式動力學軟件LS-DYNA模擬了鋼環(huán)在土體中的頂伸過程。

隨著鋼環(huán)在土體中的深入，土體在鋼環(huán)作用下會產生一定的流動，土體流動的速度矢量圖見圖11，鋼環(huán)前側土體的流動更為明顯。

圖11　土體流動速度矢量圖 Fig.11 Velocity vector of soil fluid

在鋼環(huán)頂推方向的前側設置了四個壓力測量點，壓力時程分析時，壓力值取四點的平均值，每兩個測量點之間相隔90度，見圖12。為保證鋼環(huán)的穩(wěn)定性及軸心位置的精準度，鋼環(huán)頂推速度控制在1mm/s，頂伸長度為2585mm。鋼環(huán)尾部在液壓缸推動作用下穩(wěn)速前進，為此在尾部施加相應的邊界條件模擬鋼環(huán)前進的過程，并保持速度的穩(wěn)定，同時考慮鋼環(huán)在支撐結構下的約束條件。

圖12　壓力測量點 Fig.12 Measuring point of pressure

見圖13，隨著鋼環(huán)的深入，前側壓力逐漸增加，并在一定時刻達到峰值，隨后有所下降，并穩(wěn)定在一定范圍以內。在0～700s之間，鋼環(huán)前側壓力穩(wěn)步增加，并在700～800s之間達到峰值，最大值約為7.9MPa，隨后鋼環(huán)前側壓力開始下降，并在1500s之后基本穩(wěn)定，穩(wěn)定在6MPa左右。

圖13 壓力時程曲線(v=1mm/s)Fig.13Timehistorycurveofmeasuringpressure(v=1mm/s)圖14 壓力時程曲線(v=2mm/s)Fig.14Timehistorycurveofmeasuringpressure(v=2mm/s)圖15 壓力時程曲線(v=5mm/s)Fig.15Timehistorycurveofmeasuringpressure(v=5mm/s)

當速度改變時，會對鋼環(huán)前側阻力產生一定的影響，峰值壓力存在一定的區(qū)別，但最終的穩(wěn)定壓力基本相同。當速度為額定工況的2倍及5倍時，最大阻力約為7.4及7.2MPa，且最終穩(wěn)定在6MPa左右，見圖14、圖15。在不同的速度下，鋼環(huán)前側壓力的變化規(guī)律是一致的，盾構機頂推液壓缸采用十九組千斤頂，每組的最大頂推力為15000kN，根據(jù)鋼環(huán)截面積，頂推壓力最大可達58.85MPa，能夠滿足工作需要。

3.2鋼環(huán)應力應變分析

鋼環(huán)在頂推過程中等效應力變化見圖16～圖18。隨著鋼環(huán)在土體中的深入，鋼環(huán)等效應力逐漸增大，400s時鋼環(huán)最大等效應力為38.4MPa，在875s左右時鋼環(huán)等效應力達到最大值，最大等效應力為88.31MPa。隨著鋼環(huán)進一步的深入，鋼環(huán)的最大等效應力有所下降，并基本穩(wěn)定在70MPa左右。鋼環(huán)最大等效應力出現(xiàn)在鋼環(huán)前端，即與土體接觸的區(qū)域。不同時刻時，鋼環(huán)等效應力的分布也不盡相同，和土體初始接觸時，鋼環(huán)主要受前側阻力的影響，當深入到一定深度，鋼環(huán)不僅受四周土體的擠壓力同時也受到土體摩擦力的影響。

圖16 鋼環(huán)等效應力分布(T=400s)Fig.16Von-misesstressofsteeltube(T=400s)圖17 鋼環(huán)等效應力分布(T=900s)Fig.17Von-misesstressofsteeltube(T=900s)圖18 鋼環(huán)等效應力分布(T=1500s)Fig.18Von-misesstressofsteeltube(T=1500s)

圖19 鋼環(huán)位移分布(T=400s)Fig.19Displacementofsteeltube(T=400s)圖20 鋼環(huán)位移分布(T=900s)Fig.20Displacementofsteeltube(T=900s)圖21 鋼環(huán)位移分布(T=1500s)Fig.21Displacementofsteeltube(T=1500s)

鋼環(huán)在頂推過程中的變形見圖19～圖21，放大比例因子為50倍。鋼環(huán)的變形主要集中在與土體接觸的區(qū)域，400s時鋼環(huán)最大變形量為3.52mm，在900s時鋼環(huán)最大變形量增大到6.11mm，1500s時為6.98mm。穩(wěn)定后鋼環(huán)的最大變形量在7mm左右。

3.3鋼環(huán)屈曲失穩(wěn)分析

3.3.1鋼環(huán)側向失穩(wěn)分析

當盾構內部結果拆卸后，土層上部如存在大型建筑或大型工程機械裝備，過大的土壓力可能會導致鋼環(huán)發(fā)生側向失穩(wěn)。

為此利用ANSYS軟件分析了鋼環(huán)側向屈曲前十階特征值及屈曲模態(tài)，前十階特征值見表2，當土體上部壓力超出1.57MPa時，鋼環(huán)發(fā)生側向失穩(wěn)。

鋼環(huán)前兩階側向屈曲模態(tài)見圖22～圖23，從圖中可以看出，當土體上部壓力達到1.57MPa時，鋼環(huán)頂部首先發(fā)生失穩(wěn)。故在進行施工時，應對鋼環(huán)的側向失穩(wěn)重點加以防范，避免鋼環(huán)側向承受過大的壓力。

表2　鋼環(huán)側向屈曲前十階特征值

圖22　鋼環(huán)側向失穩(wěn)第一階屈曲模態(tài) Fig.22 1st modal of steel tube lateral buckling

圖23　鋼環(huán)側向失穩(wěn)第二階屈曲模態(tài) Fig.23 2nd modal of steel tube lateral buckling

3.3.2鋼環(huán)軸向失穩(wěn)分析

在頂伸過程中，鋼環(huán)后部受液壓缸的驅動，前側及側面受土體的阻力及摩擦力，若阻力過大可能會出現(xiàn)軸向失穩(wěn)，為此分析了軸向屈曲前十階特征值及屈曲模態(tài)。鋼環(huán)軸向屈曲前十階特征值見表3，當頂推壓力達到508MPa時，鋼環(huán)發(fā)生軸向失穩(wěn)。

表3　鋼環(huán)軸向屈曲前十階特征值

鋼環(huán)前兩階軸向屈曲模態(tài)見圖24～圖25，從圖中可以看出，當鋼環(huán)前側阻力達到508MPa時，鋼環(huán)與土體接觸的前側發(fā)生軸向失穩(wěn)。此數(shù)值已超過鋼環(huán)材料的屈服極限強度，并根據(jù)對鋼環(huán)前側阻力的分析，可知發(fā)生軸向失穩(wěn)的可能性較小。

圖24　鋼環(huán)軸向失穩(wěn)第一階屈曲模態(tài) Fig.24 1st modal of steel tube axial buckling

圖25　鋼環(huán)軸向失穩(wěn)第二階屈曲模態(tài) Fig.25 2nd modal of steel tube axial buckling

4結論

針對大型盾構對接施工，本文基于ALE方法，模擬了盾構對接過程中鋼環(huán)頂伸過程，分析了頂伸過程中鋼環(huán)的應力應變狀況，并對鋼環(huán)的軸向及側向的穩(wěn)定性進行了分析，并得到了如下幾點結論：

(1)在鋼環(huán)頂伸的過程中，隨著鋼環(huán)的深入，前側阻力逐漸增加，在達到峰值后開始下降，并穩(wěn)定在一定范圍以內；

(2)當頂伸速度提高時，鋼環(huán)前側阻力的變化規(guī)律是一致的，其峰值壓力略有下降，但最終的穩(wěn)定壓力基本相同；

(3)鋼環(huán)和土體初始接觸時，鋼環(huán)等效應力分布主要受前側阻力的影響。當頂伸到一定深度，鋼環(huán)不僅受四周土體的擠壓力同時也受到土體摩擦力的影響，其不同時刻的最大等效應力均小于屈服極限強度。另外鋼環(huán)變形主要集中在鋼環(huán)前端和土體接觸的區(qū)域；

(4)對于鋼環(huán)的穩(wěn)定性狀況，鋼環(huán)發(fā)生軸向失穩(wěn)的可能性較小。當鋼環(huán)內部土體被挖空、盾構機內部主要支撐結構被拆卸后，對鋼環(huán)的側向失穩(wěn)應當重點加以關注。

參考文獻

[1]王國安. 盾構地中對接施工技術初探[J]. 隧道建設，2007, 2(8):536-541.

WANG Guo-an. Preliminary study on construction technology of connecting shied machines in the ground [J]. Tunnel Construction, 2007, 2(8):536-541.

[2]洪開榮. 水下盾構隧道硬巖處理與對接技術[J]. 隧道建設，2012, 32(3):361-365.

HONG Kai-rong. Case study on hard rock treatment technology and shield docking technology in boring of underwater tunnels[J]. Tunnel Construction, 2012, 32(3):361-365.

[3]Bauer R, Leung W, Barfoot T. Experimental and simulation results of wheel soil interaction for planetary rovers [J].Intelligent Robots and Systems,2007,44(9):2852-2862.

[4]Kushwaha R L, Shen J, Barfoot T. Finite element analysis of the dynamic interaction between soil and tillage tool[J].American Society of Agricultural Engineers, 1995,38(5):1315-1319.

[5]Nezami E G, Youssef M A, Zhao D W, et al. Simulation of front end loader bucket-soil interaction using discrete element method[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2007,31:1147-1162.

[6]Asaf Z, Ruhinslein D, Shmulevich I. Determination of discrete element model parameters required for soil tillage[J].Soil and Tillage Esearch,2007, 92(12):227-242

[7]Fredj A, Aaron D. A 3-dimentional continuum ALE model for soil-pipe interaction[C]//7th International Pipeline Conference, Albert, 2008.

[8]Herzl C. Contact buckling and postbuckling of thin rectangular plates[J].Journal of the Mechanics and Physics of Solids,2001,49:209-230.

第一作者潘月月女，博士生,1987年生

通信作者李正農男，教授，博士生導師,1962年生

第一作者李鄭發(fā)男，博士，高工，1983年8月生