楊文廣,蔣東翔

(清華大學電力系統(tǒng)國家重點實驗室 北京，100084)

行星齒輪典型斷齒故障的動力學仿真*

楊文廣,蔣東翔

(清華大學電力系統(tǒng)國家重點實驗室 北京，100084)

基于一種改進的行星齒輪箱集總參數(shù)模型，將斷齒故障等效到時變嚙合剛度中，建立了直齒行星齒輪的太陽輪斷齒故障、行星輪斷齒故障和內(nèi)齒圈斷齒故障的動力學模型。以某單級行星齒輪為研究對象，考慮扭轉(zhuǎn)方向外部激勵，對其正常和各斷齒故障狀態(tài)下的動力學模型進行求解，對內(nèi)齒圈垂直振動加速度信號進行分析。研究結(jié)果表明：只考慮扭轉(zhuǎn)方向的外部激勵時，正常狀態(tài)下內(nèi)齒圈最高點處的垂直振動加速度信號頻譜中幅值最高的共振峰與模態(tài)能量分布相關(guān)；在各斷齒故障狀態(tài)下，內(nèi)齒圈垂直振動加速度的有效值均大于正常值，時域波形中可看到不同形式的沖擊，包絡(luò)譜中可看到對應(yīng)的故障特征頻率及其倍頻成分。

風力機；行星齒輪箱；斷齒故障；動力學模型；故障機理

引 言

齒輪箱是雙饋風力發(fā)電機組的關(guān)鍵部件之一。統(tǒng)計表明，風力機故障停機時間20%由齒輪箱故障引起，而且齒輪箱故障維修費用高，嚴重影響風力機的運行安全和經(jīng)濟效益[1]。由于行星齒輪箱結(jié)構(gòu)復(fù)雜，同時承受時變風載和電氣負載，其振動信號為非線性時變信號，因此研究齒輪箱在各種故障情況下的振動特性對于齒輪箱的振動故障診斷具有重要意義。

根據(jù)數(shù)據(jù)來源的不同，對齒輪箱振動故障特性的研究可大致分為動力學模型驅(qū)動和數(shù)據(jù)驅(qū)動。本研究內(nèi)容屬于前者。文獻[2]將行星齒輪箱的動力學模型分為3類，即純扭轉(zhuǎn)模型、剛性多體動力學模型(也稱集總參數(shù)模型)和柔性多體動力學模型。龍泉[3]建立了一個一級行星傳動兩級平行傳動的齒輪箱的純扭轉(zhuǎn)模型，考慮剛度激勵和誤差激勵，研究了齒根裂紋和齒面點蝕故障下的齒輪箱動態(tài)響應(yīng)。Chaari等[4]建立了齒面點蝕和齒根裂紋故障狀態(tài)下的行星齒輪集總參數(shù)模型，研究了齒輪在這兩種故障激勵下的振動響應(yīng)。文獻[5]建立了嚙合剛度隨齒根裂紋的裂紋深度變化模型，采用集總參數(shù)模型研究不同裂紋深度下平行齒輪副的振動響應(yīng)。陳裴[6]利用ADAMS軟件進行船用行星齒輪太陽輪斷齒故障的仿真，指出太陽輪斷齒故障時嚙合頻率周圍出現(xiàn)調(diào)制現(xiàn)象，調(diào)制頻率為其旋轉(zhuǎn)頻率乘以行星輪數(shù)。韓振南等[7]研究了平行傳動齒輪系統(tǒng)中存在齒輪剝離故障時的動力學特性。文獻[8]基于行星齒輪的振動響應(yīng)特性和嚙合沖擊的傳遞路徑，綜合考慮各嚙合成分之間的相位差，建立了行星齒輪箱的振動信號仿真模型。與純扭轉(zhuǎn)模型和柔性多體動力學模型相比較，集總參數(shù)模型可以以相對較小的計算復(fù)雜度、獲得較準確的振動信號細節(jié)[9]，非常適合于齒輪箱故障機理研究。傳統(tǒng)的集總參數(shù)模型建立在以行星架轉(zhuǎn)速旋轉(zhuǎn)的旋轉(zhuǎn)坐標系中，存在殼體的旋轉(zhuǎn)自由度難以約束等問題。

筆者提出了一種改進的剛性多體動力學模型[10]。在此基礎(chǔ)上建立直齒行星齒輪的太陽輪斷齒故障、行星輪斷齒故障和內(nèi)齒圈斷齒故障的動力學模型，考慮扭轉(zhuǎn)方向的外部激勵和時變剛度等內(nèi)部激勵，采用隱式Newmark算法對模型求解。通過對仿真信號進行分析，從動力學角度揭示了行星齒輪箱在各種斷齒故障狀態(tài)下的動力學特性。

1 改進的行星齒輪箱集總參數(shù)模型

改進的集總參數(shù)模型直接選擇慣性坐標系作為除行星輪外的其他部件的運動參考坐標系，避免了殼體旋轉(zhuǎn)自由度難以約束的問題。對各行星輪，選擇其運動參考坐標系作為繞太陽輪公轉(zhuǎn)、同時以相同角速度自轉(zhuǎn)的動坐標系。與傳統(tǒng)模型相比，每個行星齒輪新增一個公轉(zhuǎn)方向的角位移自由度。在該模型中，各行星輪的公轉(zhuǎn)角位移與行星架的自轉(zhuǎn)角位移相同。

筆者以直齒的行星輪和太陽輪嚙合的嚙合剛度矩陣為例，簡單介紹該模型的推導(dǎo)過程，詳細推導(dǎo)過程可參考文獻[10]。記太陽輪為第1個齒輪，行星輪為第2個齒輪，則這兩個齒輪間的嚙合變形協(xié)調(diào)方程為

(1)

其中：θrevp為行星輪的公轉(zhuǎn)角位移；φi為第i個齒輪對中心線到該齒輪運動參考坐標系x軸的夾角；xi,yi,θzi分別為第i個齒輪的水平振動、垂直振動和繞軸向(z軸)的扭轉(zhuǎn)角位移；ri為第i個齒輪的節(jié)圓半徑；rc為行星架的半徑；廣義自由度q=[x1,y1,θz1,x2,y2,θz2,θrevp]。

嚙合力在第1個齒輪的各運動方向上的投影Fx1，F(xiàn)y1，Tz1的計算公式為

(2)

嚙合力在行星輪的各運動方向上的投影Fx1，F(xiàn)y1，Tz1，Trevp的計算公式為

(3)

其中：Trevp為行星輪受到的嚙合力在其公轉(zhuǎn)自由度方向上的投影。

定義廣義載荷為

F=[Fx1，F(xiàn)y1，Tz1，F(xiàn)x2，F(xiàn)y2，Tz2，Trevp]

(4)

利用式(5)，得到行星輪與太陽輪嚙合的嚙合剛度矩陣Ke，該矩陣的大小為7×7

(5)

其他的剛度矩陣和阻尼矩陣可采用類似方法推導(dǎo)得到。根據(jù)各個齒輪間的相互作用關(guān)系，進行單元裝配和約束處理，得到行星齒輪級動力學數(shù)學模型為

(6)

其中：M，C，K分別為總體質(zhì)量、總體阻尼和總體剛度矩陣；F(t)為外部激勵。

當阻尼項C中各個阻尼取值比較困難時，采用瑞利阻尼近似

C=αM+βK

(7)

其中：α，β為瑞利阻尼系數(shù)，根據(jù)具體情況選擇。

為提高數(shù)值求解的穩(wěn)定性，將式(6)進行變換

(8)

其中：Feq(t)為等效激勵；F(t)為外部激勵；K(t)q為內(nèi)部激勵；Kman為時不變的綜合剛度矩陣。

2 斷齒故障時變嚙合剛度模型

齒輪副在一個時刻可能有多對齒參與嚙合，一對嚙合的直齒的嚙合剛度隨嚙合位置在不斷變化。對于行星齒輪系統(tǒng)，還需要考慮各對嚙合齒輪之間的相對相位關(guān)系。因此，采用一個簡化的行星齒輪時變嚙合剛度模型[10]。對于斷齒故障，當斷齒參與嚙合時，其提供的嚙合力為0，可通過在時變嚙合剛度模型中加入這一過程來模擬斷齒故障。

以使用的行星齒輪為仿真對象，假設(shè)行星架轉(zhuǎn)速恒為120 r/min，基于上述時變剛度模型仿真得到正常狀態(tài)、太陽輪單齒斷全齒故障、內(nèi)齒圈單齒斷全齒故障以及1號行星輪單齒斷全齒故障的部分齒輪副的時變嚙合剛度，如圖1～4所示。將故障狀態(tài)下的時變嚙合剛度與正常狀態(tài)下的相比，可以看到斷齒會導(dǎo)致時變剛度中出現(xiàn)周期性的負沖擊。從圖2看到，當太陽輪發(fā)生斷齒故障時，太陽輪與不同行星輪的時變嚙合剛度中由于斷齒導(dǎo)致的時變剛度負沖擊之間存在120°相位差。從圖4發(fā)現(xiàn)，斷齒行星輪分別與太陽輪和行星輪嚙合時的時變嚙合剛度的負沖擊之間存在180°相位差。

圖1 正常行星齒輪中部分齒輪對的時變嚙合剛度Fig.1 Mesh stiffness of some gear pairs in the healthy planetary gear

圖2 太陽輪斷齒故障時部分齒輪對的時變嚙合剛度Fig.2 Mesh stiffness of some gear pairs in the planetary gear with sun tooth breakage fault

圖3 內(nèi)齒圈斷齒故障時部分齒輪對的時變嚙合剛度Fig.3 Mesh stiffness of some gear pairs in the planetary gear with ring tooth breakage fault

圖4 行星輪斷齒故障時部分齒輪對的時變嚙合剛度Fig.4 Mesh stiffness of some gear pairs in the planetary gear with planet tooth breakage fault

3 行星齒輪參數(shù)及固有頻率分析

選擇某單級行星齒輪進行研究[11]，該行星齒輪系統(tǒng)有3個行星輪，等間隔布置，詳細參數(shù)如表1所示。其中：下標s，p，c，r分別為太陽輪、行星輪、行星架和內(nèi)齒圈；Ksn為太陽輪與行星輪的綜合嚙合剛度；Kpw為行星輪軸承的扭轉(zhuǎn)剛度，忽略軸承的扭轉(zhuǎn)剛度。通過計算可知，該級行星輪的傳動比約為4.56，輸入側(cè)等效轉(zhuǎn)動慣量。根據(jù)經(jīng)驗選擇瑞利阻尼系數(shù)α=0．012 56，β=2．12×10-7。

表1 行星齒輪的參數(shù)Tab.1 Parameters of the planetary gear

將正常及斷齒故障狀態(tài)下的時變剛度模型分別帶入到改進的行星齒輪集總參數(shù)模型中，得到該行星齒輪在正常和各斷齒故障狀態(tài)下的動力學模型。

對該行星齒輪進行模態(tài)分析，得到的固有頻率如表2所示。該行星齒輪系統(tǒng)共6個旋轉(zhuǎn)固有頻率和6個平動固有頻率。由于模型中各軸承的水平和垂直方向的剛度相同，因此12個平動固有頻率兩兩相同。

表2 固有頻率計算結(jié)果Tab.2 Results of the modal frequency analysis

4 正常狀態(tài)仿真與信號分析

只考慮扭轉(zhuǎn)方向上的外部激勵，取行星架上驅(qū)動扭矩F1(t)=10 Nm,太陽輪上負載扭矩F2為

(9)

其中：Tr為齒輪箱的傳動比；v(t)為太陽輪的當前轉(zhuǎn)速；C=120 r/min；t0為太陽輪轉(zhuǎn)速首次達到轉(zhuǎn)速C時的時間。

定性分析可知，t0時刻以前傳動鏈空載運行，t0時刻負載接入，之后齒輪箱上的總扭矩達到動態(tài)平衡，轉(zhuǎn)速將在轉(zhuǎn)速C上下波動?？梢怨浪鉻0≈CJ/F1=0．156s。選擇定步長隱式Newmark算法[12]對模型進行求解，仿真步長取1×10-6s，時長為5 s，采樣頻率為25 kHz。

圖5 正常狀態(tài)下各部件轉(zhuǎn)速Fig.5 The rotational speed of some components in the healthy planetary gear

進行正常狀態(tài)下行星齒輪的瞬態(tài)動力學分析，得到部分組件轉(zhuǎn)速的前2 s的時域波形如圖5所示。仿真得到t0≈0．158 s，與理論估計值之間誤差約為1.3%。取各轉(zhuǎn)速最后1 s的數(shù)據(jù)進行平均，作為各部件的平衡轉(zhuǎn)速，如表3所示。其中，各個理論值是以行星架轉(zhuǎn)速120 r/min為基礎(chǔ)，通過傳動比關(guān)系計算得到?？梢钥吹剑鞣抡嬷蹬c理論值之間的誤差在1%以內(nèi)。

表3 轉(zhuǎn)速的計算值與理論值對比Tab.3 Comparisons between the calculation rotational speeds and theoretical rotational speeds

實際振動傳遞器采集到的振動信號應(yīng)對應(yīng)仿真的垂直振動加速度和旋轉(zhuǎn)振動加速度信號的合成，為簡便起見，筆者取內(nèi)齒圈最高點的垂直振動加速度1～5s的信號進行信號分析，該點處旋轉(zhuǎn)運動在垂直方向上投影為0。其時域波形、頻譜和包絡(luò)譜分別如圖6所示。由于仿真模型只考慮了扭轉(zhuǎn)方向的外部激勵以及內(nèi)齒圈柔性等因素，因此仿真得到的垂直振動加速度很小?？梢钥吹剑畲蟮姆逯翟?103 Hz處，其次為743，1894和2277 Hz處，高頻段也有兩個共振區(qū)。與表2中的固有頻率分析結(jié)果比較發(fā)現(xiàn)，這些頻率值與橫振固有頻率相吻合。

圖7為包絡(luò)譜幅值的局部放大圖,可以看到微弱的192 Hz嚙合頻率成分，低頻區(qū)中的成分比較雜亂。這說明該頻譜中的主要能量集中在固有頻率附近，驗證了行星齒輪振動信號中的共振調(diào)制現(xiàn)象。

圖6 正常狀態(tài)下內(nèi)齒圈垂直振動的仿真結(jié)果Fig.6 The simulated signal of the ring vertical vibration acceleration in normal state

圖7 正常狀態(tài)下內(nèi)齒圈垂直加速度頻譜和包絡(luò)譜局部放大圖Fig.7 Partial enlarged view of spectrum and envelop amplitude spectrum of the ring vertical vibration acceleration in normal state

對頻譜中的共振調(diào)制現(xiàn)象進行進一步分析。圖8為系統(tǒng)各個平動模態(tài)的能量隨組件的分布。模態(tài)能量的定義可參考文獻[13]，圖(a)～(f)分別對應(yīng)第1個到最后一個平動固有頻率。每一副子圖均有4列，第1～4列分別對應(yīng)該固有頻率的模態(tài)能量在行星架、內(nèi)齒圈、太陽輪和所有行星輪上的分布?？梢园l(fā)現(xiàn)，1103 Hz(對應(yīng)著子圖(b))的模態(tài)能量在內(nèi)齒圈上分布的能量最大，其次為743 Hz。這說明內(nèi)齒圈振動信號上的各固有頻率附近的能量值與該固有頻率對應(yīng)的模態(tài)分布在該部件上的能量大小具有一定的關(guān)系。

圖8 各平動模態(tài)的能量隨組件的分布Fig.8 Energy distribution over components of each modal

5 斷齒故障仿真與信號分析

采用相同的外部激勵，對各斷齒狀態(tài)下的行星齒輪進行仿真，仿真時長為5 s，得到的各個轉(zhuǎn)速時域波形圖與圖5類似。行星架轉(zhuǎn)速為120 r/min時，該行星齒輪太陽輪斷齒故障特征頻率fs約為7.11Hz，內(nèi)齒圈斷齒故障特征頻率fr約為2 Hz，該頻率同時也是行星架的旋轉(zhuǎn)頻率fc，行星輪斷齒故障特征頻率fp約為5.48 Hz。仍采用1～5 s的內(nèi)齒圈垂直振動加速度信號進行信號分析，各狀態(tài)下加速度信號的時域波形如圖9所示。對比發(fā)現(xiàn)：各斷齒故障產(chǎn)生的沖擊導(dǎo)致該加速度的峰值明顯增大；不同的齒輪斷齒故障沖擊造成的波形不同。

圖9 各個狀態(tài)下該加速度信號的時域波形Fig.9 Comparisons among the waveforms of in different states

以太陽輪斷齒故障為例進行信號分析。該狀態(tài)下內(nèi)齒圈垂直振動加速度的時域波形、頻譜和包絡(luò)譜如圖10所示。通過與正常狀態(tài)下該信號的對比發(fā)現(xiàn)，頻譜上1 103 Hz頻率處的振動幅值增長最多，導(dǎo)致高頻部分的成分被淹沒。圖11為故障狀態(tài)下該信號的頻譜和包絡(luò)譜的局部放大圖。與正常狀態(tài)下的信號對比發(fā)現(xiàn)，故障狀態(tài)下信號的頻譜中存在明顯的共振調(diào)制現(xiàn)象。從包絡(luò)譜的局部放大圖可以看到，除行星架的轉(zhuǎn)頻外，被調(diào)制的頻率在低頻段主要為7.11Hz以及其倍頻，而且其3倍頻較為突出，與行星輪的個數(shù)一致。

內(nèi)齒圈斷齒故障和行星輪斷齒故障的內(nèi)齒圈垂直振動加速度信號中也存在類似的共振調(diào)制現(xiàn)象。圖12為不同狀態(tài)下包絡(luò)譜局部放大圖的對比。可以發(fā)現(xiàn)，不同故障狀態(tài)下被調(diào)制的頻率成分不同，太陽輪和行星輪斷齒故障狀態(tài)下，被調(diào)制的主要頻率是對應(yīng)的故障特征頻率，在其倍頻成分中，3倍頻成分較為突出，與行星輪的個數(shù)一致。行星輪斷齒故障狀態(tài)下，故障特征頻率主要為其故障特征頻率的偶數(shù)倍頻，其中2倍頻成分最突出，這是由于斷齒的行星輪在一個旋轉(zhuǎn)周期中先后與太陽輪和內(nèi)齒圈嚙合的原因。

圖10 太陽輪斷齒故障狀態(tài)下內(nèi)齒圈垂直振動加速度仿真結(jié)果Fig.10 The simulated signal of ring vertical vibration acceleration in the sun tooth break fault state

圖11 太陽輪斷齒故障狀態(tài)下內(nèi)齒圈垂直振動頻譜和包絡(luò)譜的局部放大圖Fig.11 Partial enlarged view of spectrum and envelop spectrum of ring vertical vibration acceleration in the sun tooth break fault state

圖12 內(nèi)圈圈垂直振動加速度的包絡(luò)譜局部放大圖Fig.12 Partial enlarge view of the envelop spectra of the ring vertical vibration acceleration

共振解調(diào)分析(demodulated resonance analysis，簡稱DRA)是齒輪箱故障信號分析常用的方法，筆者選擇高頻區(qū)5 000～10 000 Hz的頻帶進行共振解調(diào)分析，得到各狀態(tài)下包絡(luò)譜信號的局部放大圖如圖13所示。與包絡(luò)譜對比可以發(fā)現(xiàn)，共振解調(diào)得到的頻譜更加干凈，故障特征頻率更加突出；對行星輪斷齒故障，共振解調(diào)分析得到的故障特征頻率的偶數(shù)倍頻成分的幅值相對較大。這些現(xiàn)象與文獻[14]的定性分析和實驗結(jié)論一致。

圖13 內(nèi)圈圈垂直振動加速度的共振解調(diào)分析局部放大圖Fig.13 Partial enlarge view of the DRA of the ring vertical vibration acceleration

6 結(jié) 論

1) 只考慮扭轉(zhuǎn)方向的外部激勵時，齒輪箱最高點處垂直振動加速度信號中共振頻率處的能量與其對應(yīng)的固有頻率的模態(tài)能量分布存在一定的關(guān)系。

2) 正常狀態(tài)下，該振動加速度的調(diào)制頻率中存在嚙合頻率，但較為微弱，低頻區(qū)中的信號無明顯規(guī)律。

3) 當出現(xiàn)斷齒故障時，該振動加速度的有效值增大，頻譜中出現(xiàn)了明顯的共振調(diào)制現(xiàn)象。包絡(luò)譜分析顯示，被調(diào)制的頻率成分為對應(yīng)的故障特征頻率及其倍頻。

4) 對太陽輪和內(nèi)齒圈斷齒故障，除1倍頻外，其3倍頻成分也較為突出，與行星輪的個數(shù)一致；對行星輪斷齒故障，故障特征頻率的2倍頻成分最為突出，這些結(jié)論與相關(guān)文獻中的分析一致。

5) 建立的行星齒輪典型斷齒故障的動力學模型可推廣用于復(fù)雜外部激勵作用下的故障機理研究以及多級行星齒輪箱和風力機傳動鏈的故障機理研究。

[1] 陳雪峰，李繼猛，程航，等.風力發(fā)電機狀態(tài)監(jiān)測和故障診斷技術(shù)的研究與進展[J].機械工程學報，2011,47(9):45-52.

Chen Xuefeng,Li Jimeng,Cheng Hang,et al.Research and application of condition monitoring and fault diagnosis technology in wind turbines[J].Journal of Mechanical Engineering,2011,47(9):45-52.(in Chinese)

[2] Peeters J.Simulation of dynamic drive train loads in a wind turbine[D].Belgium:Katholieke Universiteit Leuven,2006.

[3] 龍泉.風電機組齒輪傳動系統(tǒng)動態(tài)特性及故障診斷方法研究[D].保定：華北電力大學,2012.

[4] Chaari F,Fakhfakh T,Haddar M.Dynamic analysis of a planetary gear failure caused by tooth pitting and cracking[J].Journal of Failure Analysis and Prevention，2006,6(2):73-78.

[5] Wu Siyan,Zuo Mingjian,Parey A.Simulation of spur gear dynamics and estimation of fault growth[J].Journal of Sound and Vibration,2008,317(3-5):608-624.

[6] 陳裴.船用行星齒輪斷齒故障機理研究[D].上海：上海交通大學,2014.

[7] 韓振南，孫文婷，高建新.含輪齒剝落的齒輪系統(tǒng)動力學故障模擬[J].振動、測試與診斷，2012,32(1):101-104.

Han Zhennan,Sun Wenting,Gao Jianxin.Dynamics fault simulation of gear transmission system including spalling[J].Journal of Vibration,Measurement &Diagnosis,2012,32(1):101-104.(in Chinese)

[8] Lei Yaguo,Tang Wei,Kong Detong,et al.Vibration signal simulation and fault diagnosis of planetary gearboxes based on transmission mechanism analysis[J].Journal of Mechanical Engineering,2014,50(17):61-68.

[9] Ambarisha V K,Parker R G.Nonlinear dynamics of planetary gears using analytical and finite element models[J].Journal of Sound and Vibration,2007,302(3):577-595.

[10] Yang Wenguang,Jiang Dongxiang.An improved rigid multibody model for the dynamic analysis of the planetary gearbox in a wind turbine[J].Shock and Vibration,2016,2016(8):1-18.

[11] Lin Jian,Parker R G.Analytical characterization of the unique properties of planetary gear free vibration[J].Journal of Vibration &Acoustics，1999,121:316-321.

[12] Gasmi A,Sprague M A,Jonkman J M,et al.Numerical stability and accuracy of temporally coupled multi-physics modules in wind-turbine CAE tools[R].[S.l.]:Office of Scientific &Technical Information Technical Reports,2013.

[13] Lin Jian,Parker R G.Sensitivity of planetary gear natural frequencies and vibration modes to model parameters[J].Journal of Sound and Vibration,1999,228(1):109-128.

[14] 馮志鵬，趙鐳鐳，褚福磊.行星齒輪箱齒輪局部故障振動頻譜特征[J].中國電機工程學報,2013,33(5):119-127.

Feng Zhipeng,Zhao Leilei,Chu Fulei.Vibration spectral characteristics of localized gear fault of planetary gearboxes [J].Proceedings of the CSEE,2013,33(5):119-127.(in Chinese)

10.16450/j.cnki.issn.1004-6801.2017.04.019

* 國家自然科學基金資助項目(51174273)

2015-12-16；

2016-02-23

TK83;TH17

楊文廣，男，1988年9月生，博士生。主要研究方向為風力機故障診斷。曾發(fā)表《An improved rigid multibody model for the dynamic analysis of the planetary gearbox in a wind turbine》(《Shock and Vibration》2016,No.8)等論文。 E-mail:ywg16@mail.tsinghua.edu.cn

蔣東翔，男，1963年2月生，教授、博士生導(dǎo)師。主要研究方向為動力系統(tǒng)設(shè)備故障診斷技術(shù)研究與應(yīng)用。 E-mail:jiangdx@mail.tsinghua.edu.cn