市政污泥化學調理裝置CFD模擬

許琪，鄧超，時亞飛，徐新宇，虞文波，酈超，陳燁，梁莎，胡敬平，何姝，王榮，楊昌柱，楊家寬

（1華中科技大學環境科學與工程學院，湖北 武漢 430074；2宇星科技發展有限公司，廣東 深圳518057）

市政污泥化學調理裝置CFD模擬

許琪1，鄧超1，時亞飛1，徐新宇1，虞文波1，酈超1，陳燁1，梁莎1，胡敬平1，何姝2，王榮2，楊昌柱1，楊家寬1

（1華中科技大學環境科學與工程學院，湖北 武漢 430074；2宇星科技發展有限公司，廣東 深圳518057）

市政污泥機械脫水前多采用化學調理，調理裝置攪拌過程中的流場分布直接影響調理效果。采用多重參考系-流體體積模型對實驗室污泥調理罐（內徑為40 cm）的流場特性進行計算流體力學（computational fluid dynamics，CFD）模擬研究，驗證了實驗室調理罐模型氣液兩相分布情況。實驗室調理罐的模擬結果表明：單層直槳葉調理罐的攪拌槳直徑與調理罐內徑比優化值為1:1.75，擋流板寬度與調理罐內徑比優化值為1:20，可有效消除流體旋渦，增強湍動作用。在此基礎上，對中試調理罐（內徑為110 cm）的流場特性進行模擬并與實際測量值校驗。該研究采用的多重參考系-流體體積模型可為同類污泥調理罐設計與優化提供參考。

化學調理；氣液兩相流；計算流體力學；多重參考系方法；流體體積模型；攪拌容器

引 言

隨著中國對生活污水處理廠的出廠污泥含水率要求提高，污泥脫水技術得到了更多關注。污泥化學調理是污泥脫水技術中必不可少的步驟[1-2]。污泥調理罐作為常用的調理裝置，在為污泥與化學調理劑的混合提供場所的同時，其攪拌效果直接影響污泥調理效果及調理劑的用量[3-4]。由于調理罐形式多樣、內部流場復雜等原因，目前關于污泥調理過程中調理罐的實驗研究較少。

近年來隨著計算機技術的飛速發展，CFD技術也得到巨大發展與應用，其作為流體力學研究領域的一種先進手段，通過模擬機械設備內部的流場來研究結構對其性能的影響，具有經濟可靠、快速方便等特點[5-8]。Ranade等[9]對下壓式斜葉渦輪槳攪拌槽內流場進行了研究，發現網格尺寸對數值模擬的精確性有一定的影響；一些研究者采用k-ε湍流模型及滑移網格法對各類攪拌反應器內流場特性進行研究[10-13]；此外還有采用類似模型結合功率系數等改進攪拌性能[14-15]；Aubin等[16]使用顆粒追蹤方法對Intermig槳葉攪拌器進行了研究并通過改變槳葉結構對流場優化；Kelly等[17]使用多重參考系方法對簡單軸流槳葉的流場進行研究；陳佳等[18]采用歐拉-歐拉模型對氣液反應器內氣含率進行了研究，李良超等[19]也采用此模型對氣泡尺寸進行了研究。然而，目前有關于攪拌罐的研究大多忽略了氣液兩相界面的變化，沒有考慮攪拌過程中空氣對攪拌效果的影響，另一方面大多研究缺少模擬與實際情況的對比分析，對模型的準確性有待驗證。

在污泥調理罐的研究過程中，攪拌槳直徑（簡稱槳徑）與擋流板寬度是影響攪拌流場的重要因素[20]。一方面，槳徑與罐徑比過大或無擋流板時會造成徑向傳遞作用不明顯，攪拌過程容易形成凹液面，而凹液面吸入大量的空氣可能引起攪拌軸的脈動降低攪拌效果，嚴重時會引起攪拌軸的變形和失效。在滿足液相充分混合的前提下，盡量避免液面形成渦流造成能量的損失。另一方面，槳徑與罐徑比過小或擋流板寬度過寬時壁面附近速度較小，不利于流動及混合效果。因此選擇合適的槳徑與擋流板寬度有助增強調理效果，減少能量損失。

本文采用氣液兩相流模型對實驗室與中試污泥調理罐分別進行了流場特性研究，并在驗證模型的基礎上對實驗室污泥調理罐槳徑與罐徑比及擋流板寬度與罐徑比進行優化設計。

1 數值模擬方法

使用商用CFD軟件Fluent，采用三維模型模擬罐內流體的流動，使用有限體積法來離散方程，壓力速度耦合使用SIMPLE算法，即求解壓力耦合方程的半隱式法。假設攪拌罐內工作介質為空氣和水，且是均勻的。

1.1 控制方程

流體的流動滿足質量守恒定律、動量守恒定律和能量守恒定律。這3個定律在計算流體動力學中的體現就是相應的流體動力學基本控制方程組。

氣相和液相的分界面的追蹤采用VOF模型，體積分數函數αq定義為單元內第q相流體所占有體積與該單元的體積之比。對于攪拌罐內氣液兩相界面變化的問題，由于只有氣液兩相，故q=1或q=2。αq滿足方程

其中，ρ是流體的密度，t是時間，U是速度矢量，u、v和w是速度矢量U在x、y和z方向上的分量，η是流體的動力黏度，Su、Sv、Sw是3個動量方程的廣義源項，T是溫度，λ是流體的熱導率，cp是流體的比熱容，ST是黏性耗散項。

1.2 重要模型設置

（1）多重參考系法設置

多重參考系法（moving reference frame，簡稱MRF）是一種穩態近似法，是對旋轉流場模擬最常用的方法[21]，它需要將計算區域分為動靜兩部分。攪拌槳葉附近區域采用旋轉參考系，旋轉速度根據實際旋轉情況設置，其他區域使用靜止參考系，兩個區域交界面定義為Interface。

（2）多相流模型設置

根據調理罐的特點，調理罐內多相流的模擬屬于自由表面流[22]，因此多相流模型選用流體體積模型（volume of fluid，簡稱VOF）。VOF模型是一種在固定的歐拉網格下的表面跟蹤方法，通過求解單獨的動量方程和處理穿過每一流體的體積分數來模擬兩種或三種不能混合的流體。模型設置中，第一相設置為空氣相，第二相設置為水相，VOF模型采用隱式求解策略，啟用體積力方程，并設置模型初始液位高度。

2 計算模型

2.1 實驗室調理罐模型

2.1.1 幾何模型及優化方案 實驗室調理罐模型如圖1所示，為圓柱體結構，攪拌槳與擋流板均為4片成90°夾角布置。

實驗室用調理罐內徑（以下簡稱罐徑D0）為40 cm。對實驗室調理罐優化的方案共8種，先選取4種不同槳徑（Dp）與罐徑比優化，在此優化結果基礎上，再選取4種不同擋流板寬度（Wb）與罐徑比優化。具體優化方案見表1。

圖1 實驗室調理罐示意Fig.1 Diagram of bench scale conditioning tank (unit: cm)

表1 實驗室調理罐優化方案Table 1 Optimization programs of bench scale conditioning tank

網格劃分采用非結構四面體網格。設置3種不同尺寸網格進行模擬計算（槳葉區域網格做加密處理），槳葉區域與槳葉外區域尺寸分別為0.8、0.4、 0.3 cm與1.2、1.0、0.8 cm，對應網格數量為635773、1134344、3607747，其平行于縱軸的一條直線上的速度分布如圖2所示。其中，1134344個網格的劃分方案滿足模擬精度與時效要求。

2.1.2 邊界條件設置 對調理罐采取無進出口設置，對罐體內壁設置為無滑移固壁，攪拌軸和攪拌槳葉設置為移動壁面，根據實際情況設置旋轉速度為10.5 rad·s?1，模型初始液位高度為27 cm。湍流模型選用標準k-ε模型。

2.1.3 實驗驗證 對比模擬與實驗結果的凹液面位置作為驗證模型的方法。將凹液面最低點距罐頂的距離定義為凹液面高度。在實驗中測量此高度值與模擬得到兩相界面分布進行對比。

2.2 中試調理罐模型

2.2.1 幾何模型 中試調理罐模型如圖3，其為圓柱體加弧形底部結構，3片成120°夾角雙層槳葉，4片成90°夾角擋流板。

網格劃分采用非結構四面體網格。設置3種不同尺寸網格模擬計算，槳葉區域與槳葉外區域尺寸分別為1.2、0.8、0.8 cm與2.0、2.0、1.5 cm，對應網格數量為3158772、5340719、9188666，其平行于縱軸的一條直線上速度分布如圖4所示。其中，5340719個網格的劃分方案滿足模擬精度與時效要求。

2.2.2 邊界條件設置 根據實際情況設置旋轉速度為5 rad·s?1，模型初始液位高度為105 cm，其他邊界條件設置與實驗室調理罐相同。

2.2.3 實驗驗證 采用三維探頭的聲學多普勒流速儀（Flow Tracker，精度達到0.001 m·s?1）對調理罐內流場速度進行測試。如圖3（b）所示，受實際罐體上方開口限制，在罐體一側設置上中下3層每層A、B、C 3點共9個測點進行對比分析，上中下3層分別為水面下10、20、30 cm處。對比實測流速與模擬流速結果對模型結果進行驗證分析。

圖2 實驗室調理罐模型不同網格模擬結果Fig.2 Different grids simulation results in bench scale conditioning tank

圖3 中試調理罐示意Fig.3 Sketch of pilot-scale conditioning tank (unit: cm)

圖4 中試調理罐模型不同網格模擬結果Fig.4 Different grids simulation results in pilot-scale scale conditioning tank

3 模擬結果分析

3.1 實驗室調理罐槳徑與罐徑比模擬與優化

3.1.1 兩相界面分布 模擬兩相界面分布如圖5左側圖所示，顏色由藍色漸變為紅色代表兩相中氣相體積分數從低至高，氣相體積分數是指在調理罐中氣相體積與兩相總體積之比。根據經驗情況，認為氣相體積分數為50 %左右處的界面為模擬得到的兩相分界面，即圖中虛線處。

圖5 不同Dp/D0參數下模擬氣液兩相界面分布與實驗結果照片Fig.5 Two-phase interface maps and experimental results under differentDp/D0parameters

因此，Paddle-0在攪拌時形成了最大的凹液面，空氣相的分布范圍較大，攪拌效果不佳。隨著槳葉半徑的減小，凹液面的高度也隨著降低，對比圖5（c）、（d），槳葉的大幅度的減短并沒有導致凹液面的大幅度降低，其原因是在Paddle-2的條件下，凹液面已經降低到了一個較低水平。

圖5右側圖為與優化方案相同條件下實際調理罐攪拌情況，測量凹液面最低點距罐頂的高度與模擬結果進行對比，結果如表2所示。可以看到模擬結果與實測結果規律一致，且相對誤差均在10 %左右，驗證了模擬結果是比較準確的。

3.1.2 軸向速度 在垂直于調理罐中心位置截取縱

截面，在縱截面上軸向速度云圖對比如圖6。隨著槳葉半徑的減小，軸向速度逐漸增加，Paddle-0軸向速度最小，Paddle-2、Paddle-3的軸向速度明顯高于前兩個方案。在它們的槳葉中心區域與外側區域形成軸向速度方向相反且速度差異較大的情況，這樣能形成上下層流動循環，有助于增強湍動作用。

表2 模擬與實驗條件下凹液面高度Table 2 Concave surface height under simulation and experimental conditions

圖6 不同Dp/D0參數下軸向速度對比Fig.6 Axial velocity comparison chart under differentDp/D0parameters

圖7 不同Dp/D0參數下整體速度對比Fig.7 Overall velocity comparison chart under differentDp/D0parameters

圖8 不同Wb/D0參數下氣液兩相界面模擬分布Fig.8 Two-phase interface maps under differentWb/D0parameters

3.1.3 整體速度 在調理罐槳葉中部水平截取截面，不同方案下槳葉中部橫截面速度云圖如圖7。整體速度隨著槳葉長度的減短而逐漸減小，Paddle-0整體速度最大，在槳葉端部達到最大值1.8m·s?1左右，Paddle-1與Paddle-2整體速度差別不大，Paddle-3整體速度偏小。Paddle-1、Paddle-2、Paddle-3相比Paddle-0從槳葉端部到罐壁的速度過渡更為緩和，這是導致凹液面減小的主要原因之一。

綜上所述，隨著槳葉長度的減小，攪拌流場的整體速度降低，但軸向速度有增加的趨勢，說明減小槳葉長度有利于流體湍動和提高攪拌效果。Paddle-0和Paddle-1軸向速度過小，Paddle-2和Paddle-3的軸向速度較大且兩者相差不大，但Paddle-3的整體速度要比Paddle-2小而不利于混合傳質過程。因此，選擇1:1.75（Paddle-2）作為最優槳徑與罐徑比。

3.2 實驗室調理罐擋流板寬度與罐徑比模擬與優化

3.2.1 兩相界面分布 設置擋流板后的兩相界面分布情況如圖8。增加擋流板可以減弱凹液面的產生，且隨著擋流板寬度的增加，其作用越明顯。圖8右側圖為相同條件下實際調理罐情況，兩相界面分布隨擋流板寬度增加的變化趨勢與模擬結果是吻合的。

3.2.2 軸向速度 軸向速度對比如圖9，增加擋流板后軸向速度大小與方向發生明顯變化。從速度方向來看在軸向界面上存在一個循環，靠近軸心區域的流體向上流動，靠近壁面的流體向下流動，這樣形成了較強上下層湍動作用。對比來看，Baffle-1中軸向速度較小且主要為切向方向，Baffle-2的軸向速度最大，Baffle-3的軸向速度略微減小。

3.2.3 整體速度 增加擋流板之后整體速度有所降低，但在槳葉的端部形成了一個高速轉動的旋渦，對比沒有增加擋流板的情況攪拌效果有了明顯的改善（如圖10所示）。4種方案對比發現Baffle-2的旋渦發展最完善，而Baffle-3離擋流板近的區域速度幾乎為零，形成了較大的速度死角，并不利于攪拌效果的提高。

綜上所述，隨著擋流板寬度的增加，整體速度降低但軸向速度有較大幅度增加，Baffle-2的軸向速度最大。相比整體速度的減小，軸向速度的增加更有利于攪拌效果的提高，并且從實驗結果看Baffle-2的湍動作用比另3種方案劇烈。因此，選擇1:20（Baffle-2）作為最優擋流板寬度與罐徑比。

圖9 不同Wb/D0參數下軸向速度對比Fig.9 Axial velocity comparison chart under differentWb/D0parameters

3.3 中試調理罐結果分析

3.3.1 兩相界面分布 兩相界面分布如圖11，從圖中無法看到明顯的兩相界面分界線，可能由于模型較為復雜，模擬結果精度不高造成。

3.3.2 實驗驗證結果分析 實測速度與模擬結果對比如圖12，可以發現實測的速度大小與模擬結果差別較大，實測的速度整體較小，且并沒有與模擬結果相似的規律。其原因可能是該模型較實驗室調理罐復雜，采用該方法模擬的精度不高，而且多重參考系法中流體速度是通過角速度大小直接定義，與角速度、半徑有關，但與實際中攪拌過程的流體運動有所區別。因此，結合實測與模擬結果，認為對復雜結構調理罐模擬的方法需要進一步研究。

圖10 不同Wb/D0參數下整體速度對比Fig.10 Overall velocity comparison chart under differentWb/D0parameters

圖11 中試調理罐兩相界面分布Fig.11 Two-phase interface map of pilot-scale conditioning tank

4 結 論

根據本文對實驗室調理罐進行優化研究，結合實驗數據驗證，確定單層直葉槳葉調理罐最佳槳徑與罐徑比為1:1.75，最佳擋流板寬度與罐徑比為1:20，在此條件下可以有效消除流體旋渦，增強湍動作用，使攪拌混合調理效果達到最優。結合實驗室與中試調理罐研究發現，多重參考系-多相流模型能準確模擬流場中旋渦的形成。

符 號 說 明

cp——流體比熱容，J·kg?1·K?1

Dp——調理罐內徑，cm

D0——調理罐槳徑，cm

ST——黏性耗散項，J

Su，Sv，Sw——分別為u、v、w方向動量方程廣義源項，N

T——溫度，K

t——時間，s

U——速度矢量，m·s?1

u——直角坐標系下水平方向速度，m·s?1

v——直角坐標系下縱向方向速度，m·s?1

Wb——擋流板寬度，cm

w——直角坐標系下垂直方向速度，m·s?1

η——動力黏度，Pa·s

λ——流體熱導率，W·m?1·K?1

ρ——流體密度，kg·m?3

[1] Zhang Weijun, Xiao Ping, Liu Yuanyuan, Xu Shiwei, Xiao Feng, Wang Dongsheng, Christopher W K Chow. Understanding the impact of chemical conditioning with inorganic polymer flocculants on soluble extracellular polymeric substances in relation to the sludge dewaterability [J].Separation and Purification Technology, 2014, 132: 430-437.

[2] Wu Xuefeng (吳雪峰), Li Qingqing (李青青), Li Xiaoping (李小平). Management system of municipal sludge treatment and disposal [J].EnvironmentalScience & Technology(環境科學與技術), 2010, 33 (4): 186-195.

[3] Liu Huan (劉歡), Yang Jiakuan (楊家寬), Shi Yafei (時亞飛), Li Ye (李野), He Shu (何姝). Relationship of evaluation indexes of sludge dewatering performance under different conditioning programs [J].EnvironmentalScience(環境科學), 2011, 32 (11): 3394-3399.

[4] Wang Hongna, Jia Xiaoqiang, Wang Xue, Zhou Zhengxi, Wen Jianping, Zhang Jinli. CFD modeling of hydrodynamic characteristics of a gas-liquid two-phase stirred tank [J].AppliedMathematical Modelling, 2014, 38 (1): 63-92.

[5] Tamburini A, Cipollina A, Micale G, Brucato A, Ciofalo M. CFD simulations of dense solid-liquid suspensions in baffled stirred tanks: Prediction of solid particle distribution [J].ChemicalEngineering Journal, 2013, 223: 875-890.

[6] Jahoda M, Tomá?ková L, Mo?těk M. CFD prediction of liquid homogenisation in a gas-liquid stirred tank [J].ChemicalEngineering ResearchandDesign, 2009, 87 (4): 460-467.

[7] Wang Leqin (王樂勤), Du Hongxia (杜紅霞), Wu Dazhuan (吳大轉), Dai Weiping (戴維平). Numerical simulation of mixing process in stirred tank with multiple impellers [J].Journal of Engineering Thermophysics(工程熱物理學報), 2007, 28 (3): 418-420.

[8] Zhou Guozhong (周國忠), Wang Yingchen (王英琛), Shi Litian (施力田). CFD study of mixing process in stirred tank [J].Journal ofChemical Industry and Engineering(China) (化工學報), 2003, 54 (7): 886-890.

[9] Ranade V V, Bourne J R, Joshi J B. Fluid mechanics and blending in agitated tanks [J].Chemical Engineering Science, 1991, 46 (8): 1883-1893.

[10] Ma Qingshan (馬青山), Nie Yiqiang (聶毅強), Bao Yuyun (包雨云), Wang Yingchen (王英琛), Shi Litian (施力田). Numerical simulation of hydrodynamics in stirred tank [J].JournalofChemicalIndustry andEngineering(China) (化工學報), 2003, 54 (5): 612-618.

[11] Ma Qingshan (馬青山), Wang Yingchen (王英琛), Shi Litian (施力田), Wang Jiajun (王嘉駿). Measurement and numerical simulation of flow field in stirred tank with multiple impellers [J].Journalof ChemicalIndustryandEngineering(China) (化工學報), 2003, 54 (12): 1661-1666.

[12] Hou Shuandi (侯拴弟), Zhang Zheng (張政), Wang Yingchen (王英琛), Shi Litian (施力田). Numerical simulation of turbulent flow in stirred tank agitated by axial impeller [J].Journal of Chemical Industry and Engineering(China) (化工學報), 2000, 51 (1): 70-76.

[13] Sun Hui (孫會), Pan Jiazhen (潘家禎). Numerical investigation of flow fields in stirred vessels with novel combined inner and outer agitators [J].Journal of Chemical Industry and Engineering(China) (化工學報), 2006, 57 (1): 13-20.

[14] Liu Baoqing (劉寶慶), Zhang Yikun (張義堃), Liu Jingliang (劉景亮), Qin Fulei (秦福磊), Lin Xinghua (林興華), Jin Zhijiang (金志江). Numerical simulation of power consumption and mixing characteristic in stirred vessel with novel coaxial mixer [J].CIESC Journal(化工學報), 2013, 64 (4): 1135-1144.

[15] Pan Zhongyong (潘中永), Li Shuai (李帥), Zhang Hua (張華), Shi Weidong (施衛東). Numerical simulation of pulp suspension flow in stirred pot [J].Journal of Drainage and Irrigation Machinery Engineering(排灌機械工程學報), 2011, 29 (3): 204-208.

[16] Jo?lle Aubin, Catherine Xuereb. Design of multiple impeller stirred tanks for the mixing of highly viscous fluids using CFD [J].Chemical EngineeringScience, 2006, 61 (9): 2913-2920.

[17] Kelly W, Gigas B. Using CFD to predict the behavior of power law fluids near axial-flow impellers operating in the transitional flow regime [J].ChemicalEngineeringScience, 2003, 58 (10): 2141-2152.

[18] Chen Jia (陳佳), Xiao Wende (肖文德). Gas-liquid flow dynamics simulation in side-entering stirred tank [J].CIESCJournal(化工學報), 2013, 64 (7): 2344-2352.

[19] Li Liangchao (李良超), Wang Jiajun (王嘉駿), Gu Xueping (顧雪萍), Feng Lianfang (馮連芳), Li Bogeng (李伯耿). Computational fluid dynamics simulation of bubble size and local gas holdup in stirred vessel [J].JournalofZhejiangUniversity:EngineeringScience(浙江大學學報: 工學版), 2010, 44 (12): 2396-2400.

[20] Cláudio P Fonte, Bruno S Pinho, Vania Santos-Moreau, José Carlos B Lopes. Prediction of the induced gas flow rate from a self-inducing impeller with CFD [J].ChemicalEngineeringTechnology, 2014, 37 (4): 571-579.

[21] Zhou Guozhong (周國忠), Shi Litian (施力田), Wang Yingchen (王英琛). Computational fluid dynamics progress in stirred tank reactors [J].ChemicalEngineering(China) (化學工程), 2004, 32 (3): 28-32.

[22] Lao Liyun (勞力云), Zheng Zhichu (鄭之初), Wu Yingxiang (吳應湘), Li Donghui (李東暉). On the gas-liquid two-phase flow regimes and their recognition methods [J].AdvancesinMechanics(力學進展), 2002, 32 (2): 235-249.

CFD simulation of chemical conditioning unit of municipal sludge

XU Qi1, DENG Chao1, SHI Yafei1, XU Xinyu1, YU Wenbo1, LI Chao1, CHEN Ye1, LIANG Sha1, HU Jingping1, HE Shu2, WANG Rong2, YANG Changzhu1, YANG Jiakuan1
(1School of Environmental Science & Engineering,Huazhong University of Science & Technology,Wuhan430074,Hubei,China;2Universtar Science & Technology Co.,Ltd.,Shenzhen518057,Guangdong,China)

Chemical conditioning is commonly used before the mechanically dewatering process for sewage sludge. The flow field distribution of conditioning tank has a significant effect on conditioning when paddles are operating. Moving reference frame-volume of fluid model was used for simulating flow field distribution of a sludge conditioning tank with the inner diameter of 40 cm by CFD (computational fluid dynamics) technology, and the gas-liquid two-phase distribution was verified. The results showed that in the single straight blade paddle conditioning tank, the optimium ratio of paddle diameter to tank diameter was 1:1.75 and the optimum ratio of baffle width to tank diameter was 1:20. Under these optimum conditions, the fluid swirl could be effectively eliminated and the turbulent action could be enhanced. Furthermore, the flow field distribution simulation was conducted for pilot-scale conditioning tank with the inner diameter of 110 cm and the simulation results were compared with test results. Moving reference frame-volume of fluid model used in this study can provide references for design and optimization of similar sludge conditioning tank.

chemical conditioning; gas-liquid flow; computational fluid dynamics (CFD); moving referenceframe; volume of fluid (VOF); stirred vessel

Prof. YANG Jiakuan, jkyang@mail.hust.edu.cn

10.11949/j.issn.0438-1157.20150319

X 705

：A

：0438—1157（2015）10—4145—10

2015-03-16收到初稿，2015-05-09收到修改稿。

聯系人：楊家寬。

：許琪（1991—），男，博士研究生。

教育部科學技術研究項目（113046A）；中央高校基本科研業務費專項資金項目（2013TS071）。

Received date: 2015-03-16.

Foundation item: supported by the Project of Chinese Ministry of Education (113046A) and the Fundamental Research Funds for the Central Universities of Ministry of Education of China (2013TS071).