水平圓形與方形微小通道內(nèi)R134a冷凝數(shù)值模擬

張井志，李　蔚

（1浙江大學能源工程學院，浙江 杭州 310027；2先進航空發(fā)動機協(xié)同創(chuàng)新中心，浙江 杭州 310027）

?

水平圓形與方形微小通道內(nèi)R134a冷凝數(shù)值模擬

張井志1，2，李蔚1

（1浙江大學能源工程學院，浙江 杭州 310027；2先進航空發(fā)動機協(xié)同創(chuàng)新中心，浙江 杭州 310027）

利用數(shù)值模擬研究了水平圓形與方形微小通道內(nèi)R134a的冷凝換熱阻力特性，制冷劑飽和溫度為320 K。結(jié)果表明：傳熱系數(shù)與摩擦壓降梯度隨著質(zhì)量流量、干度的升高而升高，而干度大于0.85時，摩擦壓降梯度隨著干度的升高而降低。方形通道的換熱與阻力均高于圓形通道，數(shù)值結(jié)果與文獻冷凝換熱、阻力公式吻合較好。圓形通道內(nèi)冷凝液膜集聚在通道下部，而方形通道內(nèi)液膜集中在角落區(qū)域。薄液膜區(qū)域所占的比例隨著干度的增大而增大，方形通道內(nèi)的液膜厚度要小于圓形通道，換熱效果優(yōu)于圓形通道。

凝結(jié)；微小通道；數(shù)值模擬；傳熱；壓降

楊英英等［1］觀測了內(nèi)徑2 mm水平圓管內(nèi)R32冷凝流型，環(huán)形流區(qū)域隨制冷劑流量的增大而增大。Illán-Gómez等［2］研究了微小通道（dh=1.16 mm）內(nèi)R1234yf和R134a的冷凝換熱過程，實驗結(jié)果表明熱導率、汽液密度比與汽液黏度比在冷凝換熱過程中起著重要作用。Sakamatapan等［3］測定了微小通道內(nèi)R134a的冷凝傳熱系數(shù)，指出傳熱系數(shù)隨著干度、質(zhì)量流量、壁面熱流的增大而增大。Liu等［4］分析了圓形與方形通道（dh=1.152～0.952 mm）內(nèi)R152a冷凝換熱阻力特性，結(jié)果表明傳熱系數(shù)與壓降隨著飽和溫度的升高而降低，方形通道的傳熱系數(shù)要優(yōu)于圓形通道。

相比實驗方法，數(shù)值計算可以獲得內(nèi)部流動情況，進而更詳細地分析冷凝換熱特性。Wang等［5-7］研究了微小通道內(nèi)的冷凝換熱，指出在方形和三角形通道內(nèi)，表面張力的作用要強于重力及剪切力作用。Da Riva等［8-9］研究了水平光管內(nèi)的冷凝流動，指出考慮液膜區(qū)域的湍流擾動可以更好吻合實驗結(jié)果。劉納等［10］研究了R32在水平微細管（dh=1 mm）內(nèi)的冷凝，結(jié)果表明管頂部液膜厚度基本不變而底部液膜厚度逐漸變厚。El Mghari等［11］分析了微小方形通道內(nèi)的冷凝換熱情況，結(jié)果表明傳熱系數(shù)隨水力直徑的降低、壁面接觸角及通道縱橫比的增大而增大。

相對實驗研究，微小通道內(nèi)的冷凝換熱數(shù)值模擬工作仍然處于起步階段。前人針對冷凝的研究主要針對于換熱特性，而對于阻力特性的研究仍有很大的不足。本研究主要對微小圓管與方形管內(nèi)冷凝換熱現(xiàn)象進行數(shù)值研究，討論截面形狀、質(zhì)量流量、干度對其換熱阻力特性的影響，并分析了圓管與方形管內(nèi)冷凝換熱時截面液膜厚度分布。

1　數(shù)值計算模型

汽液界面采用Ansys Fluent 14.5的VOF模型進行追蹤。VOF模型質(zhì)量守恒特性較好，并具有較高的界面精度，廣泛應用于兩相流數(shù)值模擬。采用SST k-ω湍流模型考慮汽、液兩相的湍流運動。

1.1控制方程

計算區(qū)域流體控制方程如下：

連續(xù)性方程

動量方程

能量方程

體積分數(shù)方程

其中

式中，αl代表計算單元格內(nèi)液體的體積分數(shù)，αl=1，代表單元格全為液體，αl=0代表單元格全為氣體。

表面張力的作用采用Continuum Surface Force（CSF）模型考慮，其中動量方程中的體積力項F定義為

1.2汽液界面?zhèn)鳠醾髻|(zhì)模型

對于相變換熱的數(shù)值研究，關(guān)鍵點在于汽液界面?zhèn)鳠醾髻|(zhì)問題。本研究采用Lee［12］模型考慮界面的傳熱傳質(zhì)情況，此模型廣泛應用于冷凝［9-10］和沸騰［13-14］模擬。Lee［12］模型假設汽液界面溫度為飽和溫度，如式（11）所示，當汽相溫度低于飽和溫度或液相溫度高于飽和溫度時，界面產(chǎn)生質(zhì)量與熱量傳遞。由質(zhì)量傳遞所引起的能量變化，通過能量方程中汽化潛熱與質(zhì)量傳遞之積考慮。

Da Riva等［9］指出過小的r系數(shù)會導致汽液界面的溫度偏離飽和溫度，而過高的r系數(shù)則會引起計算的收斂問題。相對而言，文獻中對于蒸發(fā)、沸騰［13-14］的模擬，采用較小的r，而對于冷凝模擬則采用較大的r系數(shù)。本文模擬過程中r取值為1.5×106，其中汽液界面的溫度與飽和溫度之差小于0.5 K，傳熱傳質(zhì)模型利用Ansys Fluent中的用戶自定義函數(shù)（UDF）實現(xiàn)。

1.3幾何模型、邊界條件及離散格式

采用三維模型分析圓形與方形通道內(nèi)的冷凝情況，計算域如圖1所示，計算域長度L設定為方形通道200 mm，圓形通道400 mm。方形通道（dh=1.57 mm）與圓形通道（dh=2 mm）的單位長度換熱面積保持一致?？紤]流動換熱的對稱性，采用1/2模型，YZ截面設定為對稱邊界條件。計算域入口設定為速度入口邊界條件，出口設定為壓力出口邊界條件。管壁設定為無滑移的恒溫壁面（T=310 K），流體飽和溫度設定為320 K，制冷劑工質(zhì)的物性采用REFPROP 8.0確定。

圖1　計算模型Fig.1　Computational model

計算域采用結(jié)構(gòu)化網(wǎng)格進行劃分，壁面附近采用加強壁面函數(shù)考慮近壁區(qū)域液膜黏性底層的湍流，近壁面附近采用邊界層網(wǎng)格，以滿足壁面y+≈1.0。對網(wǎng)格進行無關(guān)性驗證，對于圓形通道，網(wǎng)格數(shù)目為191.24萬個及273.52萬個時，壁面熱流量的誤差小于0.8%?？紤]計算時間，最終圓形和方形通道的網(wǎng)格數(shù)目為191.24萬個和173.33萬個。

速度入口的速度分布、湍動能及湍流耗散利用長度為20 mm，截面網(wǎng)格分布與最終冷凝模擬模型一致的單相模擬獲得，其中單相模擬采用周期性邊界條件，以獲得充分發(fā)展狀態(tài)下入口值。入口溫度設定為飽和溫度320 K。

表1　R計算方法及離散格式Table 1　Solution method and discretization scheme

考慮模擬中重力與表面張力的計算，激活VOF模型中的隱式體積力選項，操作密度設定為汽相密度。由于采用相變模型，欠松弛因子設置為0.1～0.2，以提高計算穩(wěn)定性，汽-液界面構(gòu)造及其他參數(shù)的離散方法詳見表1。

2　結(jié)果分析與討論

2.1換熱特性

圖2為質(zhì)量流量G為317、531、715 kg·m-2·s-1時，圓形與方形通道內(nèi)冷凝傳熱系數(shù)隨干度x的變化規(guī)律。由圖2可得，傳熱系數(shù)隨著質(zhì)量流量、干度的增大而增大。在高干度（通道入口段附近）時，傳熱系數(shù)隨干度的增大趨勢要高于低干度時。在相同的質(zhì)量流量與干度下，方形通道的傳熱系數(shù)高于圓形通道。相對圓形通道，方形通道內(nèi)表面張力作用更強，促使液膜聚集在通道角落位置，導致在壁面其他位置液膜厚度大幅降低，換熱熱阻減小，換熱效果提升。同時，在相同的換熱面積下，方形通道的當量直徑要低于圓管，較小的當量直徑也是方形通道傳熱系數(shù)較高的一個原因。

圖2　傳熱系數(shù)隨干度與質(zhì)量流量的變化Fig.2　Heat transfer coefficients vs vapor quality and mass flux

2.2阻力特性

對于相變換熱，總的壓降由摩擦壓降、重力壓降及加速壓降構(gòu)成。加速壓降是由于進出口汽液含量變化而導致的壓降損失，如式（13）所示。對于橫管而言，重力壓降可以忽略，則阻力壓降梯度為總壓降梯度與加速壓降梯度之差，如式（14）所示。

圖3為不同質(zhì)量流量下，兩通道內(nèi)摩擦壓降梯度隨干度x的變化規(guī)律。由圖3可得，摩擦壓降隨著質(zhì)量流量的升高而升高，方形通道的壓降損失要高于圓形通道，且二者之差隨著質(zhì)量流量的提高而增大。當x＜0.85時，摩擦壓降隨著x上升而上升，當x＞0.85時，摩擦壓降隨著x上升而逐漸降低。由于考慮液膜區(qū)域的湍流度，在入口段附近由于液膜Reynolds數(shù)相對較小，湍流度不強，導致進口段的壓降稍低，Del Col等［15］的實驗結(jié)果也得到類似現(xiàn)象。

圖3　壓降梯度隨干度與質(zhì)量流量的變化Fig.3　Pressure drop gradients vs vapor quality and mass flux

2.3模擬結(jié)果與經(jīng)驗公式對比

為驗證數(shù)值模擬、相變傳熱傳質(zhì)模型及系數(shù)r選擇的可靠性，將數(shù)值模擬所得到的傳熱系數(shù)及摩擦壓降梯度與文獻經(jīng)驗公式進行對比。采用平均絕對誤差（MAD）與平均相對誤差（MRD）量化模擬結(jié)果與實驗結(jié)果的偏差，MAD與MRD的定義如下

式中，U代表傳熱系數(shù)與摩擦壓降梯度。

圖4為模擬傳熱系數(shù)與Cavallini等［16］和Shah［17］的經(jīng)驗公式對比情況。二者公式均對冷凝流動利用Froude數(shù)進行分區(qū)，對不同的區(qū)域采用不同的計算模型，而方法基本一致，即在管內(nèi)單相湍流強制對流的換熱關(guān)聯(lián)式（Dittus-Boelter公式）的基礎(chǔ)上添加兩相修正因子。由圖4可得，經(jīng)驗公式可以較好預測模擬結(jié)果，圓形與方形通道模擬結(jié)果基本在預測值的±30%以內(nèi)。由表2可得，Shah［17］公式可以較好預測模擬結(jié)果，圓形與方形通道的MAD分別為8.54%和10.32%。

圖4　傳熱系數(shù)模擬結(jié)果與公式的對比Fig.4　Comparison of numerical heat transfer coefficients with empirical correlations

表2　R換熱與壓降經(jīng)驗公式的MAD和MRDTable 2　MAD and MRD for heat transfer and pressure drop correlations

圖5為摩擦壓降梯度的模擬值與Zhang等［18］、Kim等［19］、Li等［20］、Zhang等［21］的經(jīng)驗公式的對比情況。4組公式均采用兩相分離模型，其中Zhang等［18］公式為實驗值擬合得到，而其他3組均為提取大量文獻實驗結(jié)果，構(gòu)建足夠豐富的壓降數(shù)據(jù)庫，分析修正Chisholm因子獲得。由圖5可得4組壓降公式可以較好預測模擬結(jié)果，其中圓形與方形通道模擬結(jié)果基本在預測值的±40%與±20%以內(nèi)。由表2可得，Li等［20］公式（MAD=14.88%，MRD=-8.76%）可以較好預測圓形通道模擬結(jié)果，而Zhang等［21］公式（MAD=10.90%，MRD=3.85%）可以較好預測方形通道模擬結(jié)果。綜合傳熱系數(shù)與壓降梯度與經(jīng)驗公式的對比，可以得到模擬傳熱傳質(zhì)模型的選擇是可靠的。

2.4汽液界面輪廓圖

本文模擬中假設壁面溫度與汽液界面溫度恒定，傳熱系數(shù)主要由液膜厚度以及液膜區(qū)域的湍流黏度確定。圖6為圓形與方形通道內(nèi)，G為317、715 kg·m-2·s-1時，汽液界面的分布情況，其中對方形通道進行了比例縮放。由圖6可得，對于圓形通道，在高干度情況下液膜厚度基本均勻。隨著干度的降低，液膜在重力的作用下逐漸向通道底部集聚，而通道頂部的液膜厚度基本不變。隨著質(zhì)量流量的增大，頂部的液膜厚度增大，而底部的液膜厚度降低。對于方形通道，液膜主要集聚在通道的角落位置，而四周管壁中間附近的液膜厚度較薄。隨著質(zhì)量流量的增大，在相同干度下，角落附近的液膜厚度逐漸降低，同時薄液膜所占通道周長的比例逐漸增大，平均換熱效果提升。

圖5　壓降梯度模擬結(jié)果與公式的對比Fig.5　Comparison of numerical pressure drop gradients with empirical correlations

2.5液膜厚度

圖7為圓形與方形通道內(nèi)，G為531 kg·m-2·s-1時，液膜厚度沿圓心角的分布情況。由圖7（a）可得，液膜厚度在圓管上部基本穩(wěn)定，約為25 μm。由于重力的作用，液膜逐漸積聚在圓管下部，厚度隨著θ的增大逐漸提高。隨著干度的降低，蒸汽冷凝量增加，重力對液膜的積聚作用增強，底部的液膜厚度逐漸增加。同時，液膜厚度驟增的位置對應的θ隨著干度的增大而逐漸增高，表明在高干度時，薄液膜區(qū)域所占比例更高。

圖6　圓形與方形通道內(nèi)汽液界面形狀Fig.6　Vapor-liquid interfaces in circular and square channels

由圖7（b）可得，方形通道內(nèi)液膜主要積聚在角落位置（θ為45°、135°），而管壁中間位置區(qū)域（θ為0°、90°、180°）液膜較薄。相對圓形通道，方形通道內(nèi)表面張力的作用更強，而重力的作用相對較小。隨著干度的降低，角落位置的液膜厚度逐漸增大，而管壁中間區(qū)域的液膜厚度以及薄液膜所占的區(qū)域逐漸降低，導致傳熱系數(shù)隨著干度的降低而降低。

對于冷凝流動，傳熱系數(shù)與液膜厚度呈反比。對于圓管，通道上側(cè)區(qū)域傳熱系數(shù)較大，蒸汽更容易在此處凝結(jié)為液體，而對于方形通道，整個換熱部分主要集中在四周管壁的中間位置。由圖7可得，方形通道內(nèi)液膜厚度小于圓形通道，換熱效果優(yōu)于圓形通道。

圖7　圓形與方形通道內(nèi)液膜厚度Fig.7　Liquid film thicknesses in circular and square channels

3　結(jié)　論

通過數(shù)值模擬研究了水平微細圓形與方形通道內(nèi)的冷凝流動換熱阻力特性，得到如下結(jié)論。

（1）傳熱系數(shù)隨著質(zhì)量流量、干度的增大而增大，方形通道的傳熱系數(shù)高于圓形通道。

（2）摩擦壓降梯度隨著質(zhì)量流量增大而提高，當x＜0.85時，隨著干度的增大而增高。同樣，方形通道的壓降要高于圓形通道。

（3）模擬所獲得的傳熱系數(shù)與摩擦壓降梯度可以較好地吻合文獻經(jīng)驗公式。

（4）圓形通道內(nèi)的液膜集聚于通道底部，而方形通道的液膜主要集中于通道角落區(qū)域。相對圓形通道，方形通道內(nèi)表面張力的作用更強，而重力的作用相對較小。

（5）薄液膜區(qū)域所占的比例隨著干度的增大而增大，方形通道內(nèi)的液膜厚度要小于圓形通道，換熱性能優(yōu)于圓形通道。

符號說明

cp——比定壓熱容，J·kg-1·K-1

dh——當量直徑，mm

E——比熱力學能，J·kg-1

F——動量方程源項，N·m-3

G——質(zhì)量流量，kg·m-2·s-1

g——重力加速度，m·s-2

h——傳熱系數(shù)，W·m-2·K-1

hlv——汽化潛熱，J·kg-1

k——熱導率，W·m-1·K-1

m——相變質(zhì)量源項，kg·m-3·s-1

MAD——平均絕對誤差

MRD——平均相對誤差

r——質(zhì)量源項系數(shù)，s-1

T——溫度，K

v——速度矢量，m·s-1

x——干度

α——單元格體積分數(shù)

δ——液膜厚度，μm

θ——角度，（°）

μ——動力黏度，kg·m-1·s-1

ρ——密度，kg·m-3

σ——表面張力，N·m-1

下角標

a——加速壓降

cal——計算值

f——摩擦壓降

l——液相

sat——飽和狀態(tài)

sim——模擬值

tp——兩相

v——汽相

References

［1］ 楊英英， 李敏霞， 馬一太. 水平光滑細管內(nèi)R32冷凝換熱的流型特性［J］. 化工學報， 2014， 65（2）： 445-452. DOI： 10.3969/j.issn. 0438-1157.2014.02.011. YANG Y Y， LI M X， MA Y T. Characteristics of flow pattern for condensation heat transfer of R32 in horizontal small tube ［J］. CIESC Journal， 2014， 65（2）： 445-452. DOI： 10.3969/j.issn.0438-1157. 2014.02.011.

［2］ ILLáN-GóMEZ F， LóPEZ-BELCHí A， GARCíA- CASCALES J R，et al. Experimental two-phase heat transfer coefficient and frictional pressure drop inside mini-channels during condensation with R1234yf and R134a ［J］. Int. J. Refrig.， 2015， 51： 12-23.

［3］ SAKAMATAPAN K， KAEW-ON J， DALKILIC A S， et al. Condensation heat transfer characteristics of R-134a flowing inside the multiport minichannels ［J］. Int. J. Heat Mass Transfer， 2013，64（0）： 976-985

［4］ LIU N， LI J M， SUN J， et al. Heat transfer and pressure drop duringcondensation of R152a in circular and square microchannels ［J］. Exp. Therm. Fluid. Sci.， 2013， 47： 60-67.

［5］ WANG H S， ROSE J W. A theory of film condensation in horizontal noncircular section microchannels ［J］. J. Heat Transfer， 2005，127（10）： 1096-1105.

［6］ WANG H S， ROSE J W. Film condensation in horizontal microchannels： effect of channel shape ［J］. Int. J. Therm. Sci.， 2006，45（12）： 1205-1212.

［7］ WANG H S， ROSE J W. Theory of heat transfer during condensation in microchannels ［J］. Int. J. Heat Mass Transfer， 2011， 54（11/12）：2525-2534.

［8］ DA RIVA E， DEL COL D. Numerical simulation of laminar liquid film condensation in a horizontal circular minichannel ［J］. J. Heat Transfer， 2012， 134（5）： 051019-051019.

［9］ DA RIVA E， DEL COL D， GARIMELLA S V， et al. The importance of turbulence during condensation in a horizontal circular minichannel ［J］. Int. J. Heat Mass Transfer， 2012， 55（13/14）：3470-3481.

［10］ 劉納， 李俊明. R32 在水平微細圓管內(nèi)凝結(jié)換熱的數(shù)值模擬［J］.化工學報， 2014， 65（11）： 4246-4253. DOI： 10.1080/10407782. 2015.1081029. LIU N， LI J M. Numerical simulation of R32 condensation heat transfer in horizontal circular microchannel ［J］. CIESC Journal， 2014，65（11）： 4246-4253. DOI： 10.1080/10407782.2015.1081029.

［11］ EL MGHARI H， ASBIK M， LOUAHLIA-GUALOUS H， et al. Condensation heat transfer enhancement in a horizontal non-circular microchannel ［J］. Appl. Therm. Eng.， 2014， 64（1/2）： 358-370.

［12］ LEE W H. A pressure iteration scheme for two-phase flow modeling［J］. Multiphase Transport Fundamentals， Reactor Safety， Applications，1980， 1： 407-431.

［13］ WEI J， PAN L， CHEN D， et al. Numerical simulation of bubble behaviors in subcooled flow boiling under swing motion ［J］. Nucl. Eng. Des.， 2011， 241（8）： 2898-2908.

［14］ YANG Z， PENG X， YE P. Numerical and experimental investigation of two phase flow during boiling in a coiled tube ［J］. Int. J. Heat Mass Transfer， 2008， 51（5）： 1003-1016.

［15］ DEL COL D， BORTOLATO M， AZZOLIN M， et al. Condensation heat transfer and two-phase frictional pressure drop in a single minichannel with R1234ze（E） and other refrigerants ［J］. Int. J. Refrig.，2015， 50： 87-103.

［16］ CAVALLINI A， DEL COL D， DORETTI L， et al. Condensation in horizontal smooth tubes： a new heat transfer model for heat exchanger design ［J］. Heat Transf. Eng.， 2006， 27（8）： 31-38.

［17］ SHAH M M. An improved and extended general correlation for heat transfer during condensation in plain tubes ［J］. HVAC&R Res.， 2009，15（5）： 889-913.

［18］ ZHANG M， WEBB R L. Correlation of two-phase friction for refrigerants in small-diameter tubes ［J］. Exp. Therm. Fluid. Sci.， 2001，25（3）： 131-139.

［19］ KIM S， MUDAWAR I. Universal approach to predicting two-phase frictional pressure drop for adiabatic and condensing mini/micro-channel flows ［J］. Int. J. Heat Mass Transfer， 2012，55（11/12）： 3246-3261.

［20］ LI W， WU Z. A general correlation for adiabatic two-phase pressure drop in micro/mini-channels ［J］. Int. J. Heat Mass Transfer， 2010，53（13/14）： 2732-2739.

［21］ ZHANG W， HIBIKI T， MISHIMA K. Correlations of two-phase frictional pressure drop and void fraction in mini-channel ［J］. Int. J. Heat Mass Transfer， 2010， 53（1/2/3）： 453-465.

Numerical simulation of condensation in horizontal circular and square minichannels using R134a

ZHANG Jingzhi1，2， LI Wei1
（1College of Energy Engineering， Zhejiang University， Hangzhou 310027， Zhejiang， China；2Collaborative Innovation Center of Advanced Aero-Engine， Hangzhou 310027， Zhejiang， China）

Heat transfer and pressure drop characteristics of condensation for R134a at saturation temperature of 320 K inside horizontal circular and square minichannels were investigated numerically. The results indicated that the heat transfer coefficients and frictional pressure drop gradients increased with mass flux and vapor quality，while the pressure drops decreased with increasing vapor quality at x＞0.85. Compared with the circular tube， the square channel with the same perimeter as the corresponding circular tubes can enhance heat transfer coefficients and increase the pressure losses. The numerical results agreed well with the well-known empirical correlations for condensation. A thicker liquid film was obtained at the bottom of the round tube， while the liquid film for the square channel accumulated at the bottom and the middle top of the channels at lower vapor quality. The proportion of the thinner liquid film region， which corresponded to a higher heat transfer coefficient， increased with increasing vapor quality， leading to a higher heat transfer coefficient at higher vapor quality. The liquid film thicknesses in square channels were lower than that in circular channels， thus enhancing the heat transfer.

condensation； minichannels； numerical simulation； heat transfer； pressure drop

引　言

冷凝換熱廣泛應用于空調(diào)、化工及電子芯片散熱等領(lǐng)域，相對單相流動，相變流動可以利用工質(zhì)的潛熱，大幅提高換熱效果。冷凝流動可以滿足電子器件微小化所帶來的散熱需求，受到相關(guān)行業(yè)的廣泛關(guān)注。相對常規(guī)通道，微小通道內(nèi)的冷凝換熱規(guī)律仍有一定不足，對其進行研究可以滿足工業(yè)以及科研的需求。

date： 2015-08-17.

Prof. LI Wei， weili96@zju.edu.cn

supported by the Natural Science Foundation of Zhejiang Province （LZ13E060001） and the National Natural Science Foundation of China （51210011）.

10.11949/j.issn.0438-1157.20151306

TK 124

A

0438—1157（2016）05—1748—07

2015-08-17收到初稿，2015-11-04收到修改稿。

聯(lián)系人：李蔚。第一作者：張井志（1988—），男，博士研究生。

浙江省自然科學基金項目（LZ13E060001）；國家自然科學基金國際合作項目（51210011）。