王俊清 袁艷平 曹曉玲 秦 萍 孫亮亮

?

地埋管井群全年非穩(wěn)態(tài)換熱特性研究

王俊清 袁艷平 曹曉玲 秦 萍 孫亮亮

（西南交通大學(xué)機械工程學(xué)院 成都 610031）

采取數(shù)值與解析計算相結(jié)合的方法，在構(gòu)建井群地埋管換熱器與地上熱泵機組耦合模型的基礎(chǔ)上，研究了動態(tài)負(fù)荷下土壤源熱泵系統(tǒng)全年運行的非穩(wěn)態(tài)特性。研究結(jié)果表明：各井換熱量隨空調(diào) 冷/熱負(fù)荷變化而上下波動，空調(diào)季各井平均壁溫總體逐漸升高，采暖季各井平均壁溫總體逐漸降低，但各井間存在一定差值；在過渡季各井平均壁溫得到一定程度恢復(fù)，秋季至少恢復(fù)了8.5℃；機組COP隨運行時間的持續(xù)，先是逐漸減小后趨于平穩(wěn)；經(jīng)過一年運行各井平均壁溫最大溫升為1.76℃，井群平均地溫升高了0.93℃。

土壤源熱泵；井群；地埋管傳熱模型；換熱量；井壁溫度；COP

0 前言

土壤源熱泵是一種綠色、節(jié)能的空調(diào)技術(shù)，對節(jié)能減排、構(gòu)建綠色建筑具有重要意義。在土壤源熱泵系統(tǒng)中，埋管換熱器的研究一直是土壤源熱泵技術(shù)的難點，同時也是該項技術(shù)研究的核心和應(yīng)用的基礎(chǔ)，實際上地源熱泵系統(tǒng)的性能不僅與地?zé)釗Q熱器的配置有關(guān)，而且也與匹配的熱泵機組的性能及建筑負(fù)荷特性有關(guān)，因此，只有對埋管換熱器、熱泵機組以及負(fù)荷特性三者進行耦合分析，才能更準(zhǔn)確合理地評價土壤源熱泵系統(tǒng)的運行性能。其中王景剛等[1]依據(jù)圓柱源理論，建立了耦合地上熱泵機組和地下埋管換熱器的模擬模型，該模型可用于長期運行的地源熱泵系統(tǒng)的短時間步長運行特性模擬；曲云霞等[2,3]采用理論解析方法建立了豎直U型埋管地?zé)釗Q熱器的準(zhǔn)三維模型，采用確定性模型法建立了熱泵機組的模型，然后通過能量和質(zhì)量守衡方程式建立了系統(tǒng)動態(tài)耦合模型，該模型可預(yù)測在不同的地?zé)釗Q熱器配置情況下的地源熱泵系統(tǒng)的各項性能；王華軍等[4]采用針對一個地源熱泵空調(diào)系統(tǒng)進行了系統(tǒng)運行特性的長期測試實驗，對冬夏季機組COP，埋管進出水溫度及地溫分布進行了研究；楊衛(wèi)波等[5]建立了土壤源熱泵系統(tǒng)及復(fù)合土壤源熱泵系統(tǒng)動態(tài)仿真模型，對土壤源熱泵系統(tǒng)各項性能參數(shù)進行了分析；劉逸等[6]建立了土壤源熱泵系統(tǒng)各部分的數(shù)學(xué)模型，結(jié)合寒冷地區(qū)的氣候特點和供暖要求，對供暖系統(tǒng)運行特性進行模擬研究，其中對地埋管換熱器性能的模擬混合應(yīng)用變熱流線源理論和圓柱源理論模型，管內(nèi)流體和鉆孔壁面之間的熱阻通過邊界元法進行求解，并考慮長期運行時地埋管換熱器之間的相互影響，從而使模擬結(jié)果更加接近實際工況；花莉等[7]利用模擬軟件TRNSYS搭建土壤源熱泵系統(tǒng)模型，對土壤源熱泵系統(tǒng)長期運行性能進行模擬分析，主要研究了土壤源熱泵長期運行的熱平衡問題；張長興等[8,9]分別通過采用有限長線熱源模型和鉆孔內(nèi)準(zhǔn)三維傳熱模型對U型管換熱器年逐時負(fù)荷工況下的換熱進行了重點分析，為采用線熱源模型進行地埋管換熱器的精確計算提供必要的參考，以及在分析地埋管換熱器傳熱基本原理的基礎(chǔ)上，在鉆孔外以“柱熱源”為基礎(chǔ)，引入格林函數(shù)法，實現(xiàn)了變負(fù)荷工況作用下，鉆孔壁溫的逐時動態(tài)卷積計算，鉆孔內(nèi)傳熱采用穩(wěn)態(tài)傳熱的形狀因子理論，通過地埋管換熱器內(nèi)的水將熱泵機組和地埋管換熱器耦合在一起，提出了一種熱泵系統(tǒng)運行特性的快速預(yù)測方法；Fordsmand M等[10]提出了垂直埋管系統(tǒng)在埋管傳熱和熱泵機組運行工況都進入穩(wěn)定情況下的土壤熱源熱泵系統(tǒng)數(shù)學(xué)模型，可計算每月的熱泵機組性能參數(shù)；Geeraert B等[11]提出了二維的線熱源、半無限大傳熱的水平盤管系統(tǒng)，在埋管傳熱和熱泵機組運行都進入穩(wěn)態(tài)工況情況下的土壤熱源熱泵系統(tǒng)數(shù)學(xué)模型，該模型與實測數(shù)據(jù)相比，計算盤管長度偏大，該穩(wěn)態(tài)模型一般作為估算用。

袁艷平等人[12-15]建立了快速求解的地埋管傳熱模型，以此為基礎(chǔ)對單井在連續(xù)運行和間歇運行下的換熱特性進行了研究，并耦合熱泵機組及負(fù)荷建立系統(tǒng)模型，對系統(tǒng)動態(tài)運行特性進行了分析。

從文獻綜述情況來看，地埋管傳熱模型大都針對單井，但在實際工程中埋管均是以群井形式出現(xiàn)，而目前關(guān)于土壤源熱泵系統(tǒng)的研究又大致分為兩種思路：一是針對某一具體實驗系統(tǒng)或工程進行分析。這種方法經(jīng)濟投入較大，理論性較差，提供的基礎(chǔ)數(shù)據(jù)又較少，難于指導(dǎo)實際的工程設(shè)計。二是通過建立系統(tǒng)模型進行仿真計算。這種方法可以任意改變工況參數(shù)或結(jié)構(gòu)參數(shù)，降低了研究成本，縮短了研究周期。

本文考慮了井群中各井間熱干擾影響，結(jié)合地埋管換熱器傳熱的基本原理，在各鉆孔外以“有限體積法”為基礎(chǔ)，實現(xiàn)了變負(fù)荷工況下，各鉆孔壁溫的逐時計算，各鉆孔內(nèi)傳熱采用準(zhǔn)三維傳熱模型計算，通過地埋管換熱器內(nèi)的水將熱泵機組和地埋管換熱器耦合在一起，建筑負(fù)荷直接加載于熱泵機組，作為系統(tǒng)的輸入，實現(xiàn)了系統(tǒng)長期運行特性的準(zhǔn)確快速預(yù)測，為地源熱泵系統(tǒng)的設(shè)計及優(yōu)化提供參考。

1 井群傳熱模型及實驗驗證

1.1 井群數(shù)學(xué)模型

井群地埋管傳熱模型采取數(shù)值計算與解析計算相結(jié)合的方法求解，以鉆孔壁為界將井群換熱空間區(qū)域分為鉆孔內(nèi)（包含多個鉆孔）和鉆孔外兩個區(qū)域；各鉆孔內(nèi)傳熱采用穩(wěn)態(tài)解析解計算，鉆孔外土壤傳熱區(qū)域采用非穩(wěn)態(tài)有限體積法進行計算，兩者通過各鉆孔壁溫度耦合連接計算。

對鉆孔內(nèi)傳熱進行以下簡化：忽略埋管與回填材料及回填材料與孔洞壁間的接觸熱阻；忽略埋管內(nèi)介質(zhì)軸向?qū)岷蚒型地埋管底部彎管的影響；管內(nèi)流體流速均勻一致，任意截面內(nèi)流體溫度均勻恒定，只沿井深方向變化；回填土、管內(nèi)流體的熱特性參數(shù)恒定；忽略熱濕遷移的影響，認(rèn)為回填土中的傳熱為純導(dǎo)熱問題。

井群中每個鉆孔內(nèi)傳熱情況完全一樣，在此僅對其中一個鉆孔為對象分析其傳熱情況，鉆孔內(nèi)微元體傳熱示意如圖1所示。

圖1 鉆孔內(nèi)微元體傳熱示意圖

基于以上簡化，考慮地埋管內(nèi)流體溫度沿程變化及兩支管“熱短路”現(xiàn)象，參照圖1對于埋管深度處的微元體，可根據(jù)鉆孔內(nèi)能量平衡得控制方程組[16]：

對方程組（1）進行簡化，通過Laplace變換可求解得出地埋管出水溫度及單位延米管長換熱量。

鉆孔外土壤區(qū)域簡化為均一純導(dǎo)熱介質(zhì)，忽略垂直方向傳熱及其他因素影響，則土壤區(qū)域傳熱控制方程為[17]：

式中，為土壤密度，；為土壤比熱，；為土壤導(dǎo)熱系數(shù)，；為土壤溫度，℃；是源項，求解中進行線性化處理，分解為常數(shù)項及隨時間和溫度變化項。

邊界條件：各鉆孔壁為變熱流邊界，遠(yuǎn)邊界為絕熱邊界。

對鉆孔外土壤計算區(qū)域進行二維網(wǎng)格離散，在控制容積內(nèi)對控制方程（2）進行空間和時間積分，組建差分方程組，結(jié)合邊界條件和初始條件對方程組進行求解，即可得到鉆孔外土壤區(qū)域溫度分布。

1.2 實驗驗證

為判定井群傳熱模型預(yù)測結(jié)果是否符合實際情況，在此通過實驗進行驗證。由井群傳熱模型可知，井群傳熱模型的數(shù)學(xué)描述與單井傳熱模型僅在邊界條件方面有差別，因此只要確保單井傳熱模型的正確，即可證明井群傳熱模型預(yù)測結(jié)果的準(zhǔn)確性。實驗系統(tǒng)實景圖如圖2所示。

圖2 實驗系統(tǒng)圖

實驗測得土壤熱物性參數(shù)及實驗系統(tǒng)運行參數(shù)如下：

土壤熱物性：=1322kg/m3，=0.49W/(m·K)，=1016J/(kg·℃)；土壤初溫：0=18.9℃；流體參數(shù)：=105L/h，fin=48.6℃。

實驗系統(tǒng)連續(xù)運行7h，每分鐘采集一次U型地埋管進出口溫度，通過計算獲得地埋管單位管長換熱量，實驗值與模型預(yù)測值對比結(jié)果如圖3，模型預(yù)測結(jié)果的誤差分析如圖4所示。

圖3 實驗值與模擬計算值對比

圖4 絕對誤差與相對誤差分析

由圖3可看出，預(yù)測結(jié)果與實驗結(jié)果變化趨勢完全一致；通過圖4進一步分析可知，運行初期（40min內(nèi)）相對誤差為5%～9%，運行穩(wěn)定后相對誤差在5%以內(nèi)，運行時間內(nèi)各運行點的絕對誤差平均值為0.9W/m，相對誤差平均值為2.4%，這進一步表明模型預(yù)測結(jié)果的可靠性。產(chǎn)生上述誤差的主要原因為：（1）實驗系統(tǒng)本身存在誤差，（2）鉆孔內(nèi)傳熱近似處理為穩(wěn)態(tài)過程，實際情況是剛啟動時鉆孔內(nèi)材料會蓄熱，傳熱為非穩(wěn)態(tài)。

2 建筑動態(tài)負(fù)荷、熱泵機組模型及計算參數(shù)

2.1 建筑動態(tài)負(fù)荷

本節(jié)以成都地區(qū)一棟普通賓館建筑為研究對象，采用Dest-c軟件模擬計算該建筑的全年逐時冷熱負(fù)荷。建筑全年逐時冷熱負(fù)荷如圖5。

圖5 建筑全年逐時空調(diào)負(fù)荷

從計算結(jié)果可知，成都地區(qū)某普通旅館峰值熱負(fù)荷為24.6kW，峰值冷負(fù)荷為52.7kW，兩者之比為0.47，全年空調(diào)熱負(fù)荷與冷負(fù)荷之比為0.43，其冷熱負(fù)荷不平衡率較大，這將對地源熱泵系統(tǒng)的長期高效運行產(chǎn)生影響。

2.2 熱泵機組模型

以建筑冷負(fù)荷為依據(jù)，選擇GSHP53型地源熱泵機組,根據(jù)樣本提供的數(shù)據(jù)，夏/冬兩季與機組冷凝器進口水溫的關(guān)系可擬合為：

制冷工況：COP=5.784+0.056t-0.002（3）

供熱工況：COP=3.257+0.133t-0.001（4）

熱泵機組制冷量與散釋到土壤中的熱量有以下關(guān)系：

同理，熱泵機組供熱量與吸收土壤中的熱量有以下關(guān)系：

式（5）和（6）中：Q為熱泵機組制冷量，kW；Q為熱泵機組供熱量，kW；t為熱泵機組進水溫度，℃；t為熱泵機組出水溫度，℃；C為循環(huán)流體比熱，kJ/(kg·℃)；為機組流量，kg/h。

2.3 計算參數(shù)

參考工程實際應(yīng)用，利用管材“換熱能力”計算所需地埋管長度，換熱能力取35W/m（管長），通過計算所需地埋管總長度為1806m，采用順排方形布置的16口井群，埋管間距為3m，井深為56m。因方形井群存在對稱性和影響共性，故相似位置的井其換熱規(guī)律是一致的，各井可依據(jù)在井群中的不同位置分為三類[18]：“中井”、“邊井”、“角井”，圖6為三類井位置平面圖。

圖6 井群中三類井位置

全年時間一般可分為采暖季、空調(diào)季與過度季，成都地區(qū)三個時間區(qū)段的劃分結(jié)果如下[19]：空調(diào)季：6月10日－8月28日；采暖季：12月1日—次年3月5日；過渡季：全年其余時間。該土壤源熱泵空調(diào)系統(tǒng)從空調(diào)季開始運行，運行一年后停止。

3 結(jié)果分析

3.1 空調(diào)季結(jié)果分析

（a）三類井單位井深換熱量變化情況

（b）特定時間三類井單位井深換熱量變化情況

圖7 空調(diào)季三類井換熱量隨運行時間變化情況

Fig.7 The heat exchange of three kinds of wells changes along with the running time in air conditioning season

圖8 空調(diào)季三類井平均壁溫隨運行時間變化情況

圖7和圖8分別為空調(diào)季三類井單位井深換熱量和井壁平均溫度隨運行時間的變化情況。從圖中可看出，三類井換熱量受建筑冷負(fù)荷影響在0～74W/m范圍內(nèi)波動，隨著土壤蓄熱的進行，三類井平均壁溫總體趨勢逐漸升高，當(dāng)空調(diào)冷負(fù)荷小到一定程度，地埋管傳熱速度小于土壤熱擴散速度，則井壁平均溫度有所下降；隨著換熱進行一段時間后，三類井換熱量及壁溫均出現(xiàn)差值，換熱量從大到小的順序為角井、邊井、中井，而井壁溫度從高到低的順序為中井、邊井、角井，其原因是隨著換熱進行各井間產(chǎn)生熱干擾，中位于井群中部其所受影響最大，而邊井與角井則相對較小。

圖9 空調(diào)季熱泵機組COP逐時變化情況

圖9為空調(diào)季熱泵機組COP逐時變化情況。從圖中可看出，機組COP隨運行時間增加，總體趨勢先逐漸減小后逐漸趨于平穩(wěn)，穩(wěn)定之后其值在4.9～5.4左右。

3.2 采暖季結(jié)果分析

（a）三類井單位井深換熱量變化情況

（b）特定時間三類井單位井深換熱量變化情況

圖10 采暖季三類井換熱量隨運行時間變化情況

Fig.10 The heat exchange of three kinds of wells changes along with the running time in heating season

圖11 采暖季三類井平均壁溫變化情況

圖10和圖11分別為采暖季三類井單位井深換熱量和平均井壁溫度隨運行時間的變化情況。從圖中可看出，三類井單位井深換熱量受建筑熱負(fù)荷影響在0～25W/m范圍內(nèi)波動，隨著土壤放熱過程的進行，三類井平均井壁溫度總體趨勢逐漸降低，整個采暖季內(nèi)三類井換熱量從大到小的順序為中井、邊井、角井，而三類井平均井壁溫從高到低的順序與其換熱量相同，但三類井換熱量和井壁平均溫度的差值隨著取熱的進行逐漸減小，這是因為經(jīng)過空調(diào)季后，井群中部地溫較高，有利地埋管換熱，但隨著取熱的進行，井群地溫趨于均衡，中部與邊沿處地溫差減小。

圖12 采暖季熱泵機組COP逐時變化情況

圖12為采暖季熱泵機組COP逐時變化情況。從圖中可知，機組COP隨運行時間增加，總體趨勢先是逐漸減小后逐漸趨于平穩(wěn)，穩(wěn)定后其值在4.85～5.15左右。

3.3 全年地溫分析

圖13 全年三類井平均壁溫變化情況

在此以井壁平均溫度為指標(biāo)探討全年地溫變化規(guī)律。圖13為全年地溫逐時變化情況。由圖可知，井群初始地溫為17.1℃，經(jīng)過空調(diào)季后，中井壁溫升高至31.15℃，邊井壁溫升高至30.74℃，角井壁溫升高至30.38℃；經(jīng)過秋季后，中井壁溫降至22.65℃，邊井壁溫降至21.37℃，角井壁溫降至20.35℃；經(jīng)過采暖季后，中井壁溫降至18.37℃，邊井壁溫降至18.30℃，角井壁溫降至18.23℃；經(jīng)過春季后，中井壁溫回升至18.86℃，邊井壁溫回升至18.73℃，角井壁溫回升至18.56℃；則經(jīng)過一年運行后，中井平均壁溫升高1.76℃，邊井平均壁溫升高1.63℃，角井平均壁溫升高1.47℃；若以平均地溫為指標(biāo)，則經(jīng)過一年運行，井群平均地溫升高了0.93℃；其原因是建筑冷熱負(fù)荷不平衡率較大，導(dǎo)致土壤蓄、放熱不平衡，對此可采用增加輔助散熱設(shè)備的復(fù)合式土壤源熱泵空調(diào)系統(tǒng)。

4 結(jié)論

（1）提出了一種可快速計算的混合解井群傳熱模型，并通過實驗驗證了其可靠性。

（2）在井群傳熱模型的基礎(chǔ)上建立了熱泵空調(diào)系統(tǒng)耦合仿真模型，從本文的計算過程來看,該系統(tǒng)耦合仿真模型具有準(zhǔn)確快速求解的特點。

（3）各井換熱量隨空調(diào)冷/熱負(fù)荷變化而上下波動，各井壁平均溫度隨蓄熱過程進行總體趨勢逐漸升高，隨取熱過程進行總體趨勢逐漸降低，因熱干擾存在各井間有些許差值；在過渡季井壁平均溫度獲得一定程度恢復(fù)，秋季各井壁平均溫度至少恢復(fù)了8.5℃；機組COP隨運行時間的持續(xù),先是逐漸減小后趨于平穩(wěn)。

（4）經(jīng)過一年運行中井平均壁溫升高1.76℃，邊井平均壁溫升高1.63℃，角井平均壁溫升高1.47℃，井群平均地溫升高0.93℃，這對土壤源熱泵空調(diào)系統(tǒng)長期高效運行產(chǎn)生不利影響，對此可采用增加輔助散熱設(shè)備的復(fù)合式土壤源熱泵空調(diào)系統(tǒng)（如增加輔助冷卻塔）。

[1] 王景剛,馬一太,張子平,等.地源熱泵的運行特性模擬研究[J].工程熱物理學(xué)報,2003,24(3):361-366.

[2] 曲云霞,方肇洪,李安桂,等.閉環(huán)地源熱泵系統(tǒng)建模[J].西安建筑科技大學(xué)學(xué)報,2007,39(2):235-239.

[3] 曲云霞.地源熱泵系統(tǒng)模擬與仿真[D].西安:西安建筑科技大學(xué),2004.

[4] 王華軍,趙軍,沈亮,等.地源熱泵系統(tǒng)長期運行特性的實驗研究[J].華北電力大學(xué)學(xué)報,2007,34(2):52-54.

[5] 楊衛(wèi)波,施明恒.混合地源熱泵系統(tǒng)(HGSHPS)的研究[J].建筑熱能通風(fēng)空調(diào),2006,25(3):20-6.

[6] 劉逸,李炳熙,付忠斌,等.寒冷地區(qū)土壤源熱泵系統(tǒng)供暖運行特性研究[J].太陽能學(xué)報,2012,33(10):1789-1795.

[7] 花莉,潘毅群,范蕊.基于TRNSYS的土壤源熱泵熱平衡問題的影響因素分析[J].建筑節(jié)能,2012,40(3):23-29.

[8] 張長興,胡松濤,李緒泉.基于線熱源模型的地源熱泵系統(tǒng)運行模擬[J].建筑科學(xué),2008,24(12):71-76.

[9] 張長興,郭占軍,劉玉峰,等.土壤源熱泵系統(tǒng)運行特性的快速預(yù)測方法[J].農(nóng)業(yè)工程學(xué)報,2012,28(24):173-178.

[10] Fordsmand M, Eggers-Lura A. Feasibility Study on the improvement of the cop of a heat pump[J]. Newways to saw energy,1980:174-180.

[11] Geeraert B, Steffens J C. The use of soil with andwithout artificial heat regeneration as a heat source for the heat pump[J]. Newways to saw energy,1980:94-101.

[12] 袁艷平,雷波,余南陽,等.地源熱泵地埋管換熱器傳熱研究(2):傳熱過程的完全數(shù)學(xué)描述[J].暖通空調(diào),2009,39(7):9-12.

[13] 袁艷平,雷波,曹曉玲,等.地源熱泵地埋管換熱器傳熱研究(3):變熱流條件下單U型地埋管換熱器的非穩(wěn)態(tài)傳熱特性[J].暖通空調(diào),2009,39(12):10-15.

[14] Yuan Yanping, Cao Xiaoling, Sun Liangliang, et al.Ground source heat pump system: A review of simulation inChina[J]. Renewable and sustainable energyreviews, 2012,16(9):6814-6822.

[15] Cao Xiaoling, Yuan Yanping, Sun Liangliang, et al.Restoration performance of vertical ground heat exchanger with various intermittent ratios[J]. Geothermics,2015,54:115-121.

[16] 楊衛(wèi)波,施明恒.二區(qū)域U型地埋管傳熱模型及其實驗驗證[J].工程熱物理學(xué)報,2008,29(5):857-860.

[17]李人憲.有限體積法基礎(chǔ)[M].北京:國防工業(yè)出版社,2008.

[18] 張丹,王發(fā)輝.基于數(shù)值計算的U型地埋管井群換熱器計算模型研究[J].可再生能源,2011,29(5):120-123.

[19] 付祥釗,黃光德,楊李寧,等.夏熱冬冷地區(qū)采暖期/空調(diào)期劃分對居住建筑能耗限值的影響分析[J].建筑科學(xué),2007,23(8):27-31.

Research on Unsteady Operating Characteristics of the Ground Source Heat Pump System

Wang Junqing Yuan Yanping Cao Xiaoling Qin Ping Sun Liangliang

( The College of Mechanical Engineering, Southwest Jiaotong University, Chengdu, 610031 )

Combining numerical and analytical calculation method, and on the basis of the establishment of wells buried tube heat exchanger and the heat pump coupling model to study the unsteady characteristics of ground source heat pump system under dynamic load throughout the year to run. Research indicates, heat exchange of all wells along with the air conditioning cold/heat load change and fluctuate, the average wall temperature of all wells increased in the air-conditioning season, the average wall temperature of each well of heating season reduce gradually, but there is a difference between the wells; the average wall temperature of each well in the transition season was recovered to some extent, and fall back to at least 8.5 ℃; the COP with the running time of the last, the first is gradually reduced, and then stabilized; after a year of operation, the maximum temperature rise of average wall temperature of each well is 1.76 ℃, the average temperature of well group increased by 0.93℃.

soil source heat pump; wells group; the heat transfer model of buried pipe; heat transfer; wall temperature; COP

1671-6612（2017）03-317-07

TK523

A

建筑環(huán)境與能源高效利用四川省青年科技創(chuàng)新研究團隊項目（2015TD0015）

王俊清（1988-），男，在讀碩士研究生，E-mail：yourongxinan@163.com

袁艷平（1973-），男，博士后，教授，博士生導(dǎo)師，E-mail：ypyuan@home.swjtu.edu.cn

2015-10-20