洪耀武，王鐵軍，韓大平，許慶彥，柳百成

(1. 清華大學 機械工程系，北京 100084；2. 清華大學 先進成形制造教育部重點實驗室，北京 100084;3. 沈陽黎明航空發(fā)動機(集團)有限責任公司，沈陽 110043)

調(diào)節(jié)片熔模鑄造過程的應力數(shù)值模擬

洪耀武1,2，王鐵軍3，韓大平3，許慶彥1,2，柳百成1,2

(1. 清華大學 機械工程系，北京 100084；2. 清華大學 先進成形制造教育部重點實驗室，北京 100084;3. 沈陽黎明航空發(fā)動機(集團)有限責任公司，沈陽 110043)

采用商用軟件 ProCAST對調(diào)節(jié)片熔模鑄造過程的應力進行數(shù)值模擬，分析調(diào)節(jié)片凝固過程中的應力變化情況，結(jié)果表明：裂紋缺陷處應力較大，凝固較慢。針對調(diào)節(jié)片裂紋的成因，提出在裂紋處增加內(nèi)澆道的改進工藝方案，并對其在改進工藝下的鑄造應力進行數(shù)值模擬。模擬結(jié)果表明：改進方案可明顯減小原方案裂紋處的應力，消除該處裂紋缺陷，提高產(chǎn)品合格率。實際生產(chǎn)表明：改進方案明顯消除了原方案的裂紋缺陷，合格率提高了30%。

鑄造應力；數(shù)值模擬；調(diào)節(jié)片；熔模鑄造

Abstract:The commercial software ProCAST was used to simulate the stress of adjustment sheets during the solidification, and then the stress change tendency of the adjustment sheets during the solidification was analyzed. The results indicate that the stress is bigger and the solidification is slower in the places where the crack defects may easily happen. The casting stress of an improved scheme for adjustment sheets was also simulated. The simulated results show that the improved scheme can reduce the casting stress and get higher yield. The real production indicates that the improved scheme can eliminate defect and increase the yield by 30%.

Key words:casting stress; numerical simulation; adjustment sheet; investment casting

調(diào)節(jié)片作為航空發(fā)動機矢量噴口的重要部件，其在熔模鑄造過程中容易產(chǎn)生裂紋缺陷，本文作者利用ProCAST軟件對鑄造過程產(chǎn)生的應力進行了數(shù)值模擬和分析，其研究對解決實際問題具有指導意義。計算機的發(fā)展，為鑄造過程的模擬計算提供了硬件基礎。通過數(shù)值模擬技術人們不僅可以“看到”鑄件的充型及凝固過程，還能通過各種判據(jù)如 Niyama判據(jù)來進行鑄件的縮孔、縮松缺陷預測，從而使改進工藝所需的試驗次數(shù)大大減少，縮短試制周期，降低生產(chǎn)成本[1]。目前，鑄造凝固過程應力場的模擬技術發(fā)展迅速，并逐漸應用于生產(chǎn)實踐。李輝等[2]采用商用軟件ProCAST對HT250應力框試件的熱應力進行了數(shù)值模擬研究。朱慧等[3]對HT200發(fā)動機剛體在兩種不同約束條件下的應力進行了模擬研究。康進武等[4]等基于FDM/FEM集成鑄件熱應力分析系統(tǒng)對HT200應力框試件的熱應力進行了數(shù)值模擬研究。SABAU等[5?7]對A365鋁合金階梯試塊熔模鑄件最終尺寸進行了模擬預測。在鑄造應力模擬中，鑄件與鑄型間的相互作用具有重要影響，基于 Ansys接觸單元法，劉小剛等[8?9]等將面面接觸單元法用于處理鑄件與鑄型間界面的相互作用。針對熔模鑄造過程模擬人們也進行了相應的研究，如 SCHEPPE和 SAHM[10]利用MAGMASOFT計算分析了不同工藝條件下NiAl基金屬間化合物熔模鑄件的溫度場，STEMMLER等[11]在熔模鑄件和型腔之間耦合了接觸模型和傳熱計算，考慮了間隙的形成對傳熱過程的影響，利用 CASTS對型殼的預熱和鑄件的凝固進行了計算。此外，GRIFFITHS[12]、HALLAM[13]、SCHENEFELT[14]和WOODBURRY等[15]對界面?zhèn)鳠岬挠嬎惴椒ê徒缑鎿Q熱系數(shù)的確定方法進行了研究。由于高溫合金的高溫力學性能尤其是糊狀區(qū)的力學性能很難得到，目前對于高溫合金鑄件的鑄造應力模擬還較少，故對調(diào)節(jié)片熔模鑄造過程中應力進行數(shù)值模擬研究具有重要意義。

本文作者采用商業(yè)軟件 ProCAST對調(diào)節(jié)片進行熔模鑄造應力數(shù)值模擬，分析調(diào)節(jié)片在凝固過程中的溫度場及應力場變化，從而得到其產(chǎn)生裂紋缺陷的原因；同時對調(diào)節(jié)片工藝改進方案進行鑄造過程數(shù)值模擬。

1 力學模型與熱學邊界條件

1.1 力學模型

目前熱應力數(shù)值模擬主要采用熱彈塑性模型，該模型不直接計入黏性效應，它認為材料屈服前為彈性，屈服后為塑性，彈性模量與屈服應力是溫度的函數(shù)，且當材料接近熔點時，彈性模量與屈服應力均變?yōu)榱?。本研究采?ProCAST軟件中提供的熱彈塑性模型。

對于彈性模型，應力和應變符合如下關系：

式中：{dσe}、[De]和{dεe}分別為彈性應力增量、彈性模量矩陣、彈性應變增量，下標e以及后面的p和t分別表示彈性、塑性及傳熱。

在熱彈塑性模型中，應變增量可表示為

且應力與總應變符合如下關系：

式中：[D]ep為彈塑性模量矩陣。

熱應變增量包括溫度變化產(chǎn)生的收縮以及彈性模量和線膨脹系數(shù)隨溫度變化而引起的增量，可以表示為

式中：T0是初始溫度；T是瞬時溫度；α是線膨脹系數(shù)。

塑性應變增量的大小和方向由流動準則確定：

式中：H′為材料的塑性硬化模量，可由簡單的拉伸曲線得到。由以上各式可得出熱彈塑性模型的本構(gòu)方程：

1.2 熱學邊界條件

熔模型殼的存在對鑄件的散熱和凝固過程具有重要影響，建立鑄件、型殼及環(huán)境相互之間的傳熱模型是對熔模鑄件凝固過程進行數(shù)值模擬的前提。澆注后，鑄件、鑄型和環(huán)境之間存在動態(tài)復雜的換熱過程：在鑄件、型殼及其二者之間的界面處，保溫氈、砂箱及其二者之間界面處，熱量(ΔQ)依靠熱傳導形式進行傳遞：

式中：h為界面換熱系數(shù)；T1和T2分別為界面兩側(cè)的溫度。

在型殼表面，保溫氈表面及砂箱表面與環(huán)境之間存在輻射和對流換熱，在鑄件冒口與環(huán)境之間也存在輻射和對流換熱過程：

式中：h1表示模殼或者砂箱等表面與環(huán)境的對流換熱系數(shù)；T表示模殼或者砂箱等表面的溫度；Ta表示環(huán)境溫度；ε表示輻射換熱系數(shù)。

關于具體的熱學邊界賦值、鑄件/鑄型力學模型研究的具體內(nèi)容及力學邊界條件設置可見文獻[16]。

2 調(diào)節(jié)片結(jié)構(gòu)分析

所研究的調(diào)節(jié)片鑄件為薄壁結(jié)構(gòu)件，采用 K424合金熔模精鑄而成，其鑄件結(jié)構(gòu)如圖1所示。該鑄件長約為360 mm、寬約為110 mm、最小壁厚0.9 mm，調(diào)節(jié)片生產(chǎn)時容易在圖 1(b)中所指位置出現(xiàn)裂紋缺陷，合格率較低，因此，有必要對調(diào)節(jié)片進行模擬研究，分析裂紋缺陷的原因，為工藝改進提供理論依據(jù)。

圖1 調(diào)節(jié)片結(jié)構(gòu)示意圖Fig.1 Schematic diagram of adjustment sheet casting: (a) Side view; (b) Top view

調(diào)節(jié)片澆注時兩件一組，對稱排布，如圖 2(a)所示。生產(chǎn)時型殼放在置有保溫氈的砂箱中。型殼、砂箱及保溫氈預熱溫度為1 070 ℃，澆注溫度為1 480℃。模擬時根據(jù)對稱性，取整個澆注系統(tǒng)的1/4，如圖2(b)所示。采用四面體網(wǎng)格，最長網(wǎng)格邊長為1 mm，網(wǎng)格總體單元數(shù)量為 7.8×105。模殼的力學模型采用線彈性模型，鑄件的力學模型采用熱彈塑性模型。

圖2 調(diào)節(jié)片鑄造工藝簡圖及有限元模型Fig.2 Casting process diagram of adjustment sheet (a) and FEM model (b)

調(diào)節(jié)片鑄件的材質(zhì)為K424，表 1～6所列分別為K424的主要合金成分、K424及型殼的部分熱物性參數(shù)及其高溫力學性能[17]。

表1 K424合金的化學成分[17]Table 1 Chemical composition of K424 alloy (mass fraction, %)[17]

表2 K424合金的熱物性參數(shù)[17]Table 2 Thermo-physical parameters of K424 alloy[17]

表3 型殼材料的熱物性參數(shù)[17]Table 3 Thermo-physical parameters of mold material[17]

表4 K424合金的彈性模量隨溫度的變化[17]Table 4 Elastic modulus of K424 alloy at different temperatures[17]

表5 K424合金的屈服強度隨溫度的變化[17]Table 5 Yield strength of K424 alloy at different temperatures[17]

表6 K424合金的硬化系數(shù)隨溫度的變化[17]Table 6 Hardening coefficient of K424 alloy at different temperatures[17]

3 模擬結(jié)果及分析

3.1 溫度場模擬結(jié)果與分析

溫度場的模擬是應力場模擬的基礎，因此，首先分析葉片溫度場分布情況。圖3所示為調(diào)節(jié)片在不同溫度區(qū)間的冷卻速率。從圖3中可以明顯看出，凝固前期調(diào)節(jié)片的底板中無肋板的中心部位冷卻速率較快，此時調(diào)節(jié)片從該處往外冷卻，后期調(diào)節(jié)片地板邊緣冷卻速率較快，調(diào)節(jié)片從底板邊緣向內(nèi)冷卻。圖 4所示為調(diào)節(jié)片冷卻不同時刻的溫度場分布。從圖 4(a)和(b)可以看出，冷卻前期調(diào)節(jié)片由底板中心向外冷卻，從圖4(c)～(e)可以看出冷卻后期調(diào)節(jié)片由底板邊緣向內(nèi)冷卻。同時從圖3(a)可以看出，1、2兩處冷速慢于周圍的，故其凝固慢于周圍的的凝固。圖5所示為調(diào)節(jié)片凝固不同時刻的固相分數(shù)分布。

圖3 調(diào)節(jié)片不同溫度區(qū)間的冷卻速率Fig.3 Cooling rate of adjustment sheet between different temperature zones: (a) 1257?1480 ℃; (b) 800?1 000 ℃

從圖 5(a)和(b)可以看出，1、2兩處(結(jié)構(gòu)上處于調(diào)節(jié)片中心兩側(cè)的懸空梁處)在周圍已經(jīng)凝固時還未完全凝固，尚處于糊狀區(qū)，而此時受周圍收縮產(chǎn)生的拉應力，當拉應力達到此兩處在該溫度下的強度極限時，就會產(chǎn)生微裂紋，在后續(xù)凝固中由于拉應力的進

一步累積，微裂紋擴展，因而1與2兩處為調(diào)節(jié)片在生產(chǎn)中容易產(chǎn)生裂紋缺陷的位置，與實際生產(chǎn)中的裂紋位置一致。對于調(diào)節(jié)片在圖5(a)中的1處截取截面YZ，從結(jié)構(gòu)上對調(diào)節(jié)片1處的應力進行分析，結(jié)果如圖6所示。

圖4 調(diào)節(jié)片的溫度場分布Fig.4 Evolution of temperature field of adjustment sheet:(a)t=152 s; (b)t=182 s; (c)t=234 s; (d)t=330 s; (e)t=1 300 s

圖5 不同時刻調(diào)節(jié)片固相分數(shù)分布Fig.5 Solid fraction of adjustment sheet at different times:(a)t=152 s; (b)t=182 s

圖6 調(diào)節(jié)片的YZ截面Fig.6 YZ cross-section of adjustment sheet casting

圖6中兩個矩形框A和B，為橫澆道、內(nèi)澆道與調(diào)節(jié)片底板之間組成的框架結(jié)構(gòu)，由框架的應力變化規(guī)律可知，橫澆道(相當于框架的粗桿)的收縮會通過內(nèi)澆道(相當于橫梁)對調(diào)節(jié)片平板(相當于細桿)產(chǎn)生拉伸作用，且由于型殼的阻礙作用，調(diào)節(jié)片底板始終處于受拉狀態(tài)。因此，圖3(a)的1處于此框架的受拉位置，該處的應力比較大，當超過該溫度下的屈服極限時就會產(chǎn)生裂紋。

3.2 應力場模擬結(jié)果與分析

基于上述溫度場的模擬分析結(jié)果，進一步對調(diào)節(jié)片凝固過程中鑄造應力開展了模擬研究。通過數(shù)值模擬得到了調(diào)節(jié)片在凝固過程中的應力場及其變化，定義Y方向正應力為σY，等效應力為σeff。

圖7所示為調(diào)節(jié)片在不同時刻的σY分布云圖。結(jié)合圖5不同時刻調(diào)節(jié)片固相分數(shù)分布圖對圖7的Y方向應力σY進行分析，由于1、2兩處周圍金屬液先凝固收縮會受到周圍阻礙其收縮的拉應力，而力的作用是相互的，所以1、2兩處會受到周圍先凝固部位的拉應力，圖7(b)顯示1、2兩處在Y方向上受拉應力，而此時圖5(b)顯示1、2兩處仍處于糊狀區(qū)，此時這兩處受的拉應力如果超過材料的屈服極限，則這兩處可能出現(xiàn)微裂紋。在后續(xù)的凝固過程中，調(diào)節(jié)片底板進一步冷卻收縮，由于型殼的阻礙作用，1、2兩處的拉應力會逐漸積累，可能使微裂紋擴展。從結(jié)構(gòu)上看，1、2兩處均處在調(diào)節(jié)片底板與橫澆道和內(nèi)澆道間組成的框架結(jié)構(gòu)中，因此，此兩處最終在Y方向上的應力由于框架使得其受拉而顯得比較大，屬于應力集中部位。

圖8所示為調(diào)節(jié)片冷卻末期的Y向應力與等效應力分布云圖。從圖8中可看出，1、2兩處屬于應力集中區(qū)，根據(jù)上述分析，1、2兩處可能會有微裂紋存在，所以，1、2兩處應力的較高水平可能使得微裂紋擴展，故此兩處容易在生產(chǎn)中產(chǎn)生裂紋，與實際生產(chǎn)裂紋位置一致。

圖7 不同凝固時刻σY分布圖Fig.7 Evolution of σYat different solidification times: (a) t=152 s; (b) t=182 s; (c) t=134 s

圖8 調(diào)節(jié)片冷卻末期應力分布Fig.8 Evolution of stress field for adjustment sheet casting at final solidification stage: (a) σY; (b) σeff

4 改進工藝方案及其模擬結(jié)果

在調(diào)節(jié)片的實際生產(chǎn)中，容易在圖3(a)所示的1、2兩處產(chǎn)生裂紋。從上述調(diào)節(jié)片凝固過程中溫度場與應力場分析看，此兩處產(chǎn)生裂紋的主要原因是這兩處相對周圍凝固較慢，受來自周圍的拉應力，且均處于調(diào)節(jié)片底板、橫澆道和內(nèi)澆道組成的框架中，凝固過程中會受到較大的拉應力。因此，對調(diào)節(jié)片提出改進的工藝方案，在圖9的A處增加內(nèi)澆道(圖9為調(diào)節(jié)片局部圖)，可去除1、2兩處的框架結(jié)構(gòu)，又可進行補縮，減小這兩處微裂紋形成的可能性，從而減少這兩處裂紋的產(chǎn)生。

圖9 調(diào)節(jié)片工藝改進方案Fig.9 Improved scheme of adjustment sheet casting

圖10所示為調(diào)節(jié)片冷卻末期原方案與工藝改進方案的等效應力分布云圖。從圖10可以看出，由于去除了上述框架結(jié)構(gòu)，1、2兩處在凝固末期的等效應力明顯減小，因此，這兩處在實際生產(chǎn)中產(chǎn)生的裂紋基本消除。

進一步從熱裂指數(shù)(具體計算方程可參考ProCAST軟件介紹)對兩種方案進行分析，圖11所示

圖10 兩種方案下調(diào)節(jié)片冷卻末期的應力分布Fig.10 Evolution of stress field for adjustment sheet casting under different schemes at final solidification stage: (a) Initial scheme; (b) Improved scheme

圖11 兩種方案下調(diào)節(jié)片熱裂指數(shù)圖Fig.11 Hot tearing indicator of adjustment sheet casting under different scheme: (a) Initial scheme; (b) Improved scheme

為兩種方案下調(diào)節(jié)片各處熱裂指數(shù)。從圖 11可以看出，改進方案可明顯的消除2處的熱裂趨勢；而1處的熱裂趨勢并未消除，但是由于1處加了內(nèi)澆道，可對1處進行補縮，且加了內(nèi)澆道之后，該處的應力明顯減小，所以，實際生產(chǎn)中，1處產(chǎn)生的裂紋減少。

經(jīng)過上述模擬分析可知，改進的工藝方案可有效減少調(diào)節(jié)片裂紋產(chǎn)生，提高調(diào)節(jié)片的合格率。沈陽黎明航空發(fā)動機(集團)有限責任公司的實際生產(chǎn)表明，改進的工藝方案可提高調(diào)節(jié)片的產(chǎn)品合格率30%。由此可見，鑄造應力數(shù)值模擬技術可以對鑄造工藝的改進提供科學依據(jù)。

5 結(jié)論

1) 對調(diào)節(jié)片凝固過程溫度場進行了模擬，冷卻前期調(diào)節(jié)片由底板心部向外冷卻，冷卻后期調(diào)節(jié)片由底板邊緣向內(nèi)冷卻。

2) 調(diào)節(jié)片1、2兩處受拉應力，且屬于較晚凝固部位，冷卻末期等效應力也較大，該處可能產(chǎn)生裂紋。

3) 改進方案可以明顯較小1、2兩處的應力和熱裂指數(shù)，減少這兩處的裂紋產(chǎn)生。實際生產(chǎn)表明，改進方案可以消除 1、2兩處裂紋，提高產(chǎn)品合格率30%。

REFERENCES

[1] 孫長波, 唐 寧, 史鳳嶺, 費 晶, 許慶彥, 柳百成. 機匣件真空熔模鑄造的數(shù)值模擬[J]. 鑄造, 2010, 59(2): 169?173.SUN Chang-bo, TANG Ning, SHI Feng-ling, FEI Jing, XU Qing-yan, LIU Bai-cheng. Numerical simulation of vacuum investment casting cartridge receiver [J]. Foundry, 2010, 59(2):169?173.

[2] 李 輝, 時建松, 張愛琴. 應力框熱應力數(shù)值模擬及變形分析[J]. 鑄造, 2010, 59(1): 38?41.LI Hui, SHI Jian-song, ZHANG Ai-qin. Numerical simulation and deformation analysis of thermal stresses in stress frame [J].Foundry, 2010, 59(1): 38?41.

[3] 朱 慧, 黃天佑, 康進武. 兩種不同約束條件下發(fā)動機缸體鑄件熱應力場的數(shù)值模擬[J]. 鑄造, 2006, 55(6): 604?607.ZHU Hui, HUANG Tian-you, KANG Jin-wu. Numerical simulation of thermal stress of diesel engine cylinder block castings under two different constraints [J]. Foundry, 2006, 55(6):604?607.

[4] 康進武, 朱 慧, 黃天佑. 基于有限差分網(wǎng)格的鑄件熱應力分析[J]. 鑄造, 2006, 55(3): 259?261.KANG Jin-wu, ZHU Hui, HUANG Tian-you. Thermal stress analysis of castings based on finite difference meshes [J].Foundry, 2006, 55(3): 259?261.

[5] SABAU A S. Alloy shrinkage factors for the investment casting process [M]. Metallurgical and Materials Transactions B, 2006,37(1): 131?140.

[6] SABAU A S, VISWANATHAN S. Material properties for predicting wax pattern dimensions in investment casting [J].Materials Science and Engineering A, 2003, 362(1): 125?134.

[7] CANNELL N, SABAU A S. Predicting pattern tooling and casting dimensions for investment casting, phase II, ORNL/TM-2005/228 [R]. Oak Ridge, Tennessee: Oak Ridge National Laboratory, 2005.

[8] 劉小剛, 康進武, 黃天佑, 邊敦亭, 張良紅. 大型軋鋼機機架凝固過程溫度場應力場模擬分析[J]. 鑄造, 2006, 55(9): 922?926.LIU Xiao-gang, KANG Jin-wu, HUANG Tian-you, BIAN Dun-ting, ZHANG Liang-hong. Numerical simulation of the solidification process of a heavy roll stand for temperature field and thermal stress field [J]. Foundry, 2006, 55(9): 922?926.

[9] 劉小剛, 康進武, 羅通國, 楊小兵, 王長建, 程建國, 黃天佑,柳百成. 水輪機下環(huán)鑄件凝固過程熱應力模擬分析[J]. 鑄造技術, 2006, 27(9): 984?987.LIU Xiao-gang, KANG Jin-wu, LUO Tong-guo, YANG Xiao-bing, WANG Chang-jian, CHENG Jian-guo, HUANG Tian-you, LIU Bai-cheng. Thermal stress analysis of the solidification process of francis turbine band casting [J].Foundry Technology, 2006, 27(9): 984?987.

[10] SCHEPPE F, SAHM P R. Comparison of the numerical simulation and the cast process of nickel aluminides [C]//Modeling of Casting, Welding and Advanced Solidification Process. Warrendale: TMS, 2000: 207?214.

[11] STEMMLER M, LASCHET G, HAAS L. Simulation of investment casting: Coupling of contact and heat transfer calculation between casting and mold during solidification [C]//Modeling of Casting, Welding and Advanced Solidification Process. Warrendale: TMS, 2000: 222?229.

[12] GRIFFITHS W D. The determination of the heat transfer coefficient in experiments involving unidirectional solidification[C]// Modeling of Casting, Welding and Advanced Solidification Process. Warrendale: TMS, 2000: 143? 150.

[13] HALLAM P C. Modeling of the interfacial heat transfer between an Al-Si alloy casting and a coated die steel [C]// Modeling of Casting, Welding and Advanced Solidification Process.Warrendale: TMS, 2000: 159?166.

[14] SCHENEFELT J R, LUCK R, TAYLOR R P. Reducing uncertainty in the use of the inverse method of determining heat transfer coefficient data [C]// Modeling of Casting, Welding and Advanced Solidification Process. Warrendale: TMS, 2000:215?221.

[15] WOODBURRY K A, PIWONKA T S, KE Q. Interfacial heat transfer coefficients during solidification of aluminum in resin-bonded and green sand molds [C]// Modeling of Casting,Welding and Advanced Solidification Process. Warrendale: TMS,2000: 270?277.

[16] 洪耀武. 鎳基高溫合金熔模鑄件鑄造應力數(shù)值模擬[D]. 北京:

清華大學, 2011: 13?22.

HONG Yao-wu. Numerical simulation of casting stress of Ni-based superalloy casting during investment casting process[D]. Beijing: Tsinghua University, 2011: 13?22.

[17] 《工程材料實用手冊》編輯委員會. 工程材料實用手冊[M].

第二卷. 北京: 中國標準出版社, 2002: 645?651.

《Engineering materials practical manual》Editorial Committee.Engineering materials practical manual [M]. 2nd ed. Beijing:China Standards Press, 2002: 645?651.

(編輯 陳衛(wèi)萍)

Numerical simulation of stress of adjustment sheet during investment casting process

HONG Yao-wu1,2, WANG Tie-jun3, HAN Da-ping3, XU Qing-yan1,2, LIU Bai-cheng1,2
(1. Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China;2. Key Laboratory for Advanced Materials Processing Technology, Ministry of Education,Tsinghua University, Beijing 100084, China;3. Shenyang Liming Aero-Engine (Group) Corporation, Shenyang 110043, China)

TG21+3

A

1004-0609(2012)07-1897-07

國家高技術研究發(fā)展計劃資助項目(2007AA04Z141)

2011-06-28；

2012-03-09

許慶彥，教授，博士；電話：13701087143；E-mail: scjxqy@tsinghua.edu.cn