張志秋,陳振華,聶旭濤,姚程煒
(中國空氣動力研究與發展中心,四川綿陽 621000)
基于流固熱耦合低溫風洞擴散段熱力學特性分析
張志秋*,陳振華,聶旭濤,姚程煒
(中國空氣動力研究與發展中心,四川綿陽 621000)
低溫風洞降溫過程中,溫度變化范圍大,容易在結構內部引發較大熱應力,影響設備運行安全。以中國空氣動力研究與發展中心0.3m跨聲速低溫風洞擴散段為研究對象,基于流固熱耦合方法,采用多物理場數據交換接口MpCCI,聯合結構有限元軟件Abaqus和計算流體動力學軟件Fluent,建立擴散段的流固熱耦合仿真模型,分析低溫內流場的換熱特性,計算低溫風洞結構的溫度及應力分布。通過低溫風洞試驗發現,流固熱耦合仿真結果接近于實際的測量結果,能夠準確反映低溫風洞結構的熱力學特性,可為低溫風洞的結構安全性能優化提供可靠的仿真分析方法。
流固熱耦合;低溫風洞;擴散段;熱力學特性分析
低溫風洞通過降低試驗氣體溫度可以提高試驗雷諾數,能夠滿足飛行器氣動性能準確模擬對雷諾數的需求,是十分重要的飛行器氣動性能地面模擬設備。試驗氣體降溫過程中,低溫風洞結構與氣流間存在著不均勻熱交換,結構局部易出現較大溫度梯度,引發熱變形,且受各類邊界約束以及氣動載荷的作用,導致結構應力超出許用范圍,從而威脅設備運行安全。另外,風洞結構變形會引起風洞回路軸線偏離、氣動型面變形和內部氣流表面出現階差等諸多問題,影響風洞的氣動性能。因此,分析和控制低溫環境下風洞結構的熱變形是低溫風洞研制的關鍵技術問題。
復雜結構的熱變形分析通常采用試驗測量和數值仿真等方法。相對而言,數值仿真能夠準確、全面地反映結構溫度場和應力場分布情況,為結構熱力學特性優化提供可靠、充分的參考依據,且成本較低,分析周期較短。所以,基于數值仿真方法,國內外學者對低溫結構的熱力學特性開展了廣泛研究。朱鴻梅[1-2]采用有限元方法,對冷卻過程中非磁性低溫杜瓦的螺紋粘接接頭進行瞬態熱應力仿真,分析了等效應力與冷卻時間的變化關系。Kim、Kang[3-4]等對復合材料/鋁環狀試件的熱力學特性進行了研究,并將低溫熱應力與有限元仿真結果作出對比分析。朱立偉[5]應用有限元分析軟件分析了LNG船用超低溫球閥在超低溫條件下的應力特性,并提出改進應力集中的措施。上述研究均采用有限元分析方法,關注的是結構的熱力學特性,而低溫環境下考慮流固熱耦合作用的熱力學特性分析尚不多見。考慮流固熱耦合作用的研究,也大多關注流體傳熱[6-7],或溫度變化較小的結構傳熱[8-9]。
本文以風洞典型結構部段——擴散段為研究對象,基于流固熱耦合方法,建立低溫風洞擴散段的流固熱耦合仿真模型,重點分析低溫下擴散段的換熱特性,計算擴散段結構的溫度及應力分布,并在風洞低溫運行過程中對擴散段溫度和應力進行測試研究。
低溫試驗在中國空氣動力研究與發展中心的0.3m跨聲速低溫風洞中進行,試驗主要研究擴散段的熱力學特性。低溫風洞正常情況外壁面有保溫層,風洞調試期間,壁面暫未加保溫層,故本文計算按未加保溫層考慮。風洞運行時,擴散段將氣流動能恢復為壓力能,以減少氣流在擴散段下游各段的能量損失[10]。擴散段結構關于垂直中心面左右對稱,由內流道、駐室、駐室進氣管道和支座組成。4個支座均為滑動支座,通過雙頭螺柱壓緊于絕熱塊上,從而限制結構在豎直方向的位移,同時保證低溫工況下在水平方向可自由收縮。內流道中通入快速流動的低溫氮氣,使流場和結構得以快速降溫。為降低駐室與內流道之間的溫度梯度,需從風洞回路向駐室引入冷卻氣體,其溫度與內流道中的氮氣接近,但流速相對較小。結構外壁面直接與空氣接觸。擴散段模型如圖1所示。
降溫過程中,擴散段換熱不均勻,導致局部存在溫度梯度,引發熱應力。本文采用流固熱耦合方法,求解結構的溫度分布;采用順序熱力耦合方法,求解結構的應力分布。

圖1 低溫風洞擴散段模型Fig.1 Structure model of cryogenic wind tunnel diffuser section
2.1 傳熱理論
傳熱有3種方式:導熱、對流和熱輻射。本文研究對象處于低溫工況,可忽略輻射效應。
流體內部的傳熱涉及3個重要物理定律:質量守恒定律、動量守恒定律和能量守恒定律,此處不再贅述。
流體的管內流動狀態取決于當地雷諾數Re,Re>104時為完全湍流[11]。內流道入口雷諾數為2.96 ×106,采用標準k-ε兩方程模型進行模擬。
本文中流體通過耦合界面與固體發生熱交換,壁溫是時間的函數,為第一類邊界條件。
結構內部傳熱方式為導熱,直角坐標系中的導熱方程通常寫作:

式中:q·為單位體積產生的熱能。無內熱源,則
邊界條件為第三類邊界條件:

式中:Tw為固體壁面溫度,Tf為膜層溫度,即對流邊界層中流體的溫度。考慮熱邊界層對換熱的影響,用T∞表示遠場溫度,Tf可用下式計算:

結構外壁面與大氣傳熱方式為自然對流;內流道結構與流體之間的傳熱方式為強制對流。駐室腔體內氮氣流速很低,與駐室殼體間的傳熱方式可視為自然對流。
等溫水平長圓柱的大空間自然對流換熱,可采用以下關聯式[11]估算:

式中:Ra為瑞利數,Pr為普朗特數,要求Ra≤1012。
2.2 流固熱耦合計算方法
多物理場耦合計算通常分為直接耦合法和間接耦合法。直接耦合是對問題涉及的所有數學物理方程進行聯合求解,一般適用于較簡單的物理問題。本文采用間接迭代耦合法,分別對流場和結構進行分析,再將耦合界面的物理量通過映射與差值[12-13]傳遞給外場進行計算,如此循環迭代。
采用多物理場代碼耦合軟件Mp CCI,聯合Abaqus與Fluent進行流固熱耦合仿真。在結構與流場耦合界面上,交換的物理量為膜層溫度Tf、換熱系數h和壁溫Tw。壁溫由Abaqus計算,經Mp CCI傳入Fluent。耦合壁面的換熱系數與膜層溫度由Fluent計算,經Mp CCI傳入Abaqus。計算采用串行算法,即Abaqus收到Fluent傳入的初始數據后開始迭代,經過固定耦合時間步長后,將計算得到的耦合量(節點溫度)傳入Fluent;Fluent收到數據后經相同時間步長的迭代,將計算得到的耦合量(換熱系數與膜層溫度)傳入Abaqus,依此循環,典型串行算法的數據傳遞過程如圖2所示。

圖2 串行耦合算法Fig.2 Serial coupling algorithms
2.3 熱應力計算
熱應力是由于部件受熱不均勻而存在溫度差異,導致各處膨脹變形或收縮變形不一致,且相互制約而產生內應力。對各項同性材料,自由膨脹的變分量為:

式中:εx、εy、εz為沿X、Y、Z方向應變分量;α為線膨脹系數。
如果物體邊界有約束條件,自由膨脹受到約束限制,微元體會產生熱應力。根據線性應力原理和胡克定律,描述應力和應變關系的物理方程為:


式中:E為材料楊氏模量;μ為泊松比;σx、σy、σz為沿X、Y、Z方向應力分量;G為材料剪切模量,G=2E(1+μ);ΔT為從無應力狀態開始的溫差。
擴散段中溫差較大區域為駐室與內流道的過渡區域,這也是應力分析重點關心的區域。
本節分別建立流體與結構傳熱模型,在Mp CCI中設置相關耦合參數,聯合Abaqus與Fluent進行流固熱耦合分析。將結構溫度場作為載荷導入Abaqus中,進行順序熱力耦合分析。
3.1 流體傳熱模型
擴散段內的流場比較復雜,對其做出以下幾點簡化。首先,駐室內氣流壓差小,且大部分位置氣流速度低于1m/s,對換熱能力的促進作用較小,因此以自然對流作為后期結構熱固耦合的邊界條件。流體傳熱模型中僅考慮內流道的流場,內流道壁面為耦合界面。其次,擴散段內流道實際入口氣流較紊亂,速度分布復雜,此處簡化為簡單的均勻來流。最后,不考慮介質的摩擦損失以及介質的物性參數隨溫度的變化。
流體傳熱模型入口馬赫數為0.2;入口溫度隨時間變化,通過Matlab對實測入口溫度進行多項式擬合,并寫入UDF,如圖3所示。出口邊界設置為outflow;壁面邊界條件設置為不滑移,溫度由Abaqus傳入。

圖3 降溫曲線與擬合降溫曲線Fig.3 Cooling curve and fitting cooling curve
湍流模型采用k-ε兩方程模型,壁面函數采用標準壁面函數,要求第一層網格厚度位于對數律區域,且Y*>11.225。經Fluent計算后,調整第一層網格厚度,直到滿足湍流模型對Y*的要求。
網格為六面體結構網格,數量為293 436,如圖4所示。為驗證網格無關性,對以下工況——入口馬赫數為0.2,入口流體溫度為270K,壁面溫度為300K——進行驗證,對比壁面平均對流換熱系數,結果如圖5所示。

圖4 流體網格Fig.4 Fluid gridding

圖5 網格無關性驗證Fig.5 Grid independence verification
3.2 結構傳熱模型
對結構作出以下幾點簡化。風洞擴散段是薄壁殼體,有限元采用實體單元對實際模型簡化少,能較好反映實際結構,但厚度方向節點不易控制,網格質量不好,計算量大。因此本文采用殼體單元建模,可提高計算效率。此外,降溫過程中,空氣中的水蒸氣會在擴散段外壁面凝華結霜,對傳熱會產生一定影響。建模中通過設置較小的對流換熱系數來近似模擬,不考慮結霜過程。最后,不考慮支座處摩擦阻力對位移的約束。
擴散段結構關于垂直中心面左右對稱,對于有限元計算,只需建立半模即可,如圖6所示。采用Abaqus前處理器對結構劃分網格,并細化駐室與內流道交界處網格,保證殼單元厚度方向至少有5個積分點[14-15],總網格數為19 714。傳熱分析采用傳熱單元DS4。結構傳熱分析中只需考慮模型的換熱問題。初始溫度設置為300K;外壁面自然對流換熱系數可按公式(4)進行估算,駐室外表面不同截面處的直徑、溫度均變化較大,計算時采用算數平均值。降溫末期,結構外表面結冰結霜會減弱傳熱效果,因此換熱系數的取值應適當減小。結構外壁面換熱系數取為5 W/(m2·K);駐室內冷卻氣體引自風洞回路,溫度與內流道氮氣溫度接近,因此將內流道入口測量溫度作為駐室氣體溫度。駐室內氣流速度小,在大部分位置速度低于1m/s,因此簡化為自然對流邊界處理。自然對流換熱系數通常在10W/(m2·K)以內,參照封閉腔內自然對流換熱工況[16]下的換熱水平,將駐室內換熱系數取為3.5W/(m2·K)。

圖6 結構網格Fig.6 Structure gridding
3.3 流固熱耦合模型
在流固熱耦合模型中,耦合界面為內流道壁面。總的降溫時間為6400s,耦合物理量為流體溫度、換熱系數和固壁溫度。
溫度從300K降至110K,每秒鐘溫度下降0.03K。軟件中,流固熱耦合時間步長設置為1s。在該時間步長內,耦合量的變化都足夠小,從而可準確反映瞬態過程中的換熱特性。
3.4 結構應力模型
熱應力分析中采用S4R殼體單元,支座處限制Y方向自由度,出口端面限制Z方向自由度,將傳熱分析得到的溫度場作為溫度載荷加載到應力分析步。
測試系統包括低溫電阻溫度傳感器、低溫應變傳感器、數據采集器和顯示系統。
溫度和應變傳感器布置在駐室與內流道外壁面的水平位置,實時監測結構的溫度與應力變化。粘貼傳感器之前,先要對粘貼表面進行處理,保證傳感器與表面的緊密貼合。然后在傳感器與鋼結構之間涂上剛性較好的低溫膠,在傳感器背面涂上彈性較好的防護膠。通過在測點附近粘貼補償塊的方式,消除溫度變化造成的傳感器誤差。
數據采集儀型號為吉時利2701,將其與工控機箱串聯,通過電腦顯示器即可設置相關參數,實時監測結構的溫度與應力變化。
試驗中測點分布如圖7所示,其中T1~T13為溫度測點。測點T1~T9位于駐室椎形體外表面水平母線上,測點T10位于冷氣入口管道與駐室連接處附近,測點T11~T13位于冷氣入口管道外壁;示意圖中三角形代表應力測點S1,位于駐室與內流道連接處。在內流道、駐室靠近內表面處布置溫度探針,監測與結構發生換熱的氣體溫度。

圖7 測點分布圖Fig.7 Diagram of measuring points distribution
降溫時間為6400s,氣流溫度由300K降至130K。此間,結構表面會結霜。隨后,風洞恢復常溫,降溫過程如圖8所示。

圖8 降溫過程Fig.8 Cooling process
本節首先將低溫風洞試驗溫度測量結果與流固熱耦合仿真結果進行對比,分析擴散段結構的傳熱特性。然后將低溫風洞試驗應力測量結果與順序熱力耦合仿真結果進行對比,分析結構應力狀態。
5.1 耦合傳熱分析
由降溫曲線的斜率可判斷降溫速率:0~2500s,溫度均勻下降;2500~3400s,降溫速率逐漸減慢;3400~4500s溫度均勻下降;4500~4800s,降溫速率突然變大,隨后恢復之前的降溫速率。
流固熱耦合仿真計算得到結構在1000、2500、5000和6400s的溫度場,如圖9所示。由溫度分布可知,駐室與內流道過渡區域的溫度梯度較大。駐室內氣流速度低,換熱效果差,且熱量傳導至溫度較低的內流道需經過結構導熱,以上原因導致駐室外壁溫度明顯高于內流道外壁溫度。

圖9 降溫過程結構溫度分布Fig.9 Distribution of structural temperature when cooling down
測點T2位于內流道外壁,T4位于駐室與內流道過渡區域。通過對比發現,仿真曲線趨勢與測量曲線保持一致,最大偏差為20K左右。2500~3400s,降溫速率減慢,溫度測量與仿真曲線的斜率也隨時間減小,符合實際情況。4500~4800s,降溫速率突然變大,在仿真曲線中可明顯看到斜率增大,而實際測量曲線并未反映出該降溫速率的變化,如圖10和11所示。

圖10 測點T2溫度變化曲線Fig.10 Temperature change curve of T2

圖11 測點T4溫度變化曲線Fig.11 Temperature change curve of T4
測點T7、T9位于駐室外壁面,由本節前面的分析可知該處降溫速率較慢。實際測量發現,溫度在250K以上,計算結果與試驗數據吻合較好,如圖12和13所示。

圖12 測點T7溫度變化曲線Fig.12 Temperature change curve of T7

圖13 測點T9溫度變化曲線Fig.13 Temperature change curve of T9
測點T10位于駐室上,靠近冷氣入口管道底部。測點T12位于冷氣入口管道上。由于管道內氣流速度較低,換熱不充分,所以該處結構溫度較高。試驗測量與仿真計算偏差較小,在10K以內,如圖14和15所示。

圖14 測點T10溫度變化曲線Fig.14 Temperature change curve of T10

圖15 測點T12溫度變化曲線Fig.15 Temperature change curve of T12
總的來說,溫度測量與仿真曲線吻合較好。最大偏差在20K以內,發生在降溫末期。
5.2 結構應力分析
結構應力的大小與溫度梯度成正比。由5.1節傳熱分析可知,駐室與內流道過渡區域應力較大。圖16為結構在1000、2500、5000和6400s時刻的應力分布。由應力分布可知,支座、法蘭以及駐室與內流道過渡區域應力較大。由于采用殼體單元建模,支座、法蘭處過于簡化,可能導致應力計算不準確。
對比應力仿真曲線與結構降溫曲線可知,應力變化與溫度變化基本對應。在0~2500s,均勻降溫,應力隨時間增大;2500~3400s,降溫速率逐漸減小,由于結構內部導熱,局部溫度梯度減小,熱應力減小;3400~4500s,均勻降溫,應力隨時間增大;4500~4800s,降溫速率突然變大后恢復,S1處應力曲線也出現一個凸起。
由圖17可知,應力仿真值與實際測量值吻合較好。最大偏差為10MPa左右,發生在降溫末期。溫度應力最大值121MPa出現在降溫末期,在駐室與支座之間。使仿真與計算產生偏差的原因,可能有以下3點:(1)風洞回路中,低溫氮氣經過試驗段后均勻性較差,導致結構傳熱不均勻;(2)降溫末期,擴散段外結霜結冰,影響應力測量;(3)防護膠經過低溫后粘合度下降,甚至出現脆裂。
通過仿真分析和試驗測試,擴散段溫度應力低于材料的許用應力(115MPa),在低溫和給定的試驗溫度變化條件下,擴散段結構是安全的。

圖16 降溫過程結構應力分布Fig.16 Distribution of structural stress when cooling down

圖17 降溫過程測量與仿真應力變化曲線Fig.17 Simulation and measuring curve of stress when cooling down

圖18 最大應力位置Fig.18 The position of maximum stress
通過對低溫風洞擴散段的數值計算與試驗得出如下2點結論:
(1)采用殼體單元對風洞擴散段的結構進行有限元建模是可行的,可以減小計算規模。盡管局部位置可能過于簡化,但應力計算值與仿真值基本在一個量級。
(2)利用第三方軟件MpCCI聯合Abaqus與Fluent進行流固熱耦合仿真,得到結構的溫度場分布;再將溫度場作為載荷,利用Abaqus進行順序熱力耦合分析,求解結構的應力場。通過風洞低溫試驗測試結果表明,該方法能夠準確反映低溫風洞結構的熱力學特性,可為低溫風洞的結構安全性能優化提供可靠的仿真分析方法。
[1]朱鴻梅,徐烈,孫恒,等.低溫粘接技術的應用及熱應力分析[J].低溫與超導,2004,32(2):29-30.Zhu H M,Xu L,Sun H,et al.Applications and thermal stress analysis of cryogenic adhesive bonding technology[J].Cryogenics and superconductivity,2004,32(2):29-30.
[2]Zhu H M,Sun H,Xu L,et al.Transient thermal stress and analysis of threaded-adhesive joints applied in non-magnetic Dewar[J].International Journal of Adhesion &Adhesives,2007,27(8):621-628.
[3]Kim M G,Kang S G,Kim C G,et al.Thermally induced stress analysis of composite/aluminum ring specimens at cryogenic temperature[J].Composites Science and Technology,2008,68(3):1080-1087.
[4]Allen R F,Bowen P.Thermo elastic analysis of a type 3 cryogenic tank considering curing temperature and autofrettage pressure[J].Journal of Reinforced Plastics and Composites,2008,27(5):459-471.
[5]朱立偉,柳建華,張良,等.LNG船用超低溫球閥的低溫應力分析及數值模擬[J].低溫與超導,2010,(5):11-14.Zhu L W,Liu J H,Zhang L,et al.Numerical simulation of stress and tightness of cryogenic valve used in LNG carrier[J].Cryogenics and Superconductivity,2010,(5):11-14.
[6]陳尊敬,王雷雷,孟華.考慮發動機冷卻通道固壁內耦合導熱影響的低溫甲烷超臨界壓力傳熱研究[J].航空學報,2013,(1):8-17.Chen Z J,Wang L L,Meng H.Study of heat transfer of cryogenic methane undersupercritical pressure with consideration of thermal conduction[J].Acta Aeronautica et Astronautica Sinica,2013,(1):8-17.
[7]王雷雷,徐可可,孟華.航空煤油超臨界壓力流/固/熱耦合數值模擬[J].工程熱物理學報,2013,(7):1357-1360.Wang L L,Xu K K,Meng H.Numerical study of conjugate heattransfer of acviation kerosene at supercritical pressures[J].Journal of Engineering Thermophysics,2013,(7):1357-1360.
[8]朱康,厲彥忠,王磊,等.飽和氫氣加注過程中低溫貯箱降溫特性及熱應力分布的數值研究[J].西安交通大學學報,2014,(5):1-7.Zhu K,Li Y Z,Wang L,et al.Investigation on cool-down behavior andthermal stress of cryogenic tank during saturated hydrogen gas filling process[J].Journey of Xi’an Jiaotong University,2014,(5):1-7.
[9]李陽,魏蔚,王彩莉,等.低溫絕熱氣瓶頸管傳熱的數值模擬與試驗[J].低溫與超導,2008,(1):9-12. Li Y,Wei W,Wang C L,et al.Numerical simulation and experimental of heat transfer in cryogenic insulated cylinders neck tube[J].Cryogenics and Superconductivity,2008,(1):9-12.
[10]劉政崇,廖達雄,董誼信,等.高低速風洞氣動與結構設計[M].北京:國防工業出版社,2003.Liu Z C,Liao D X,Dong Y X,et al.Aerodynamics and structural design of high and low speed wind tunnel[M].Beijing:Defense Industry Press,2003.
[11]Incropear F P,Dewitt D P.Fundamentals of heat and mass transfer[M].5th ed.New York:John Wiley &Sons,2002.
[12]Felippa C A,Park K C,Farhat C.Partitioned analysis of coupled mechanical systems[J].Computer Methods in Applied Mechanics and Engineering,2011,0190(24):3247-3270.
[13]Matthies H G,Steindorf J.Fully coupled fluid-structure interaction using weak coupling[J].PAMM,2002,1(1):37-38.
[14]趙騰倫.ABAQUS6.6在機械工程中的應用[M].北京:中國水利水電出版社,2007.Zhao T L.Application of ABAQUS6.6 on mechanic engineering[M].Beijing:China Water Resources and Hydropower Publishing House,2007.
[15]莊茁.基于ABAQUS的有限元分析及應用[M].北京:清華大學出版社,2009.Zhuang Z.Finite element analysis and application based on ABAQUS[M].Beijing:Tsinghua University Press,2009.
[16]楊小川.復雜熱環境中大型薄殼體內的自然對流數值模擬[D].哈爾濱工業大學,2008.Yang X C.Numerical simulation of natural convectionm inside huge thin-skin enclosers exposed to complicated thermal environment[D].Harbin:Harbin Institute of Technology,2008.
Thermodynamic characteristic analysis of the cryogenic wind tunnel diffuser section based on fluid-thermal-structural coupling
Zhang Zhiqiu*,Chen Zhenhua,Nie Xutao,Yao Chengwei
(China Aerodynamics Research and Development Center,Mianyang Sichuan 621000,China)
In cooling down process,the cryogenic wind tunnel has a wide range of temperature variation,which may lead to strong structural thermal stress and affect equipment smooth and safe operation.The cryogenic wind tunnel diffuser section is investigated in this paper based on fluid-thermal-structural coupling.The data exchange interface for multiple physical fields,MpCCI,is adopted to combine the structural finite element software,Abaqus,and Computational Fluid Dynamics(CFD)software,Fluent.The fluid-thermal-structural coupling model of the diffuser section is established to investigate heat transfer characteristics of the cryogenic internal flow field and calculate the distribution of temperature and stress in the cryogenic wind tunnel structure.The cryogenic wind tunnel experiment indicates that the results of fluid-thermal-structural cosimulation agree well with the test results,which reflect correctly the thermodynamic characteristic of cryogenic wind tunnel structure and provide reliable analysis method for safety performance optimization of wind tunnel structure.
fluid-thermal-structural interaction;cryogenic wind tunnel;diffuser section;thermodynamic characteristic analysis
V211.754
A

(編輯:楊 娟)
1672-9897(2016)06-0018-08
10.11729/syltlx20160100
2016-06-15;
2016-08-18
*通信作者E-mail:zhiqiuz@gmail.com
Zhang Z Q,Chen Z H,Nie X T,et al.Thermodynamic characteristic analysis of the cryogenic wind tunnel diffuser section based on fluidthermal-structural coupling.Journal of Experiments in Fluid Mechanics,2016,30(6):18-25.張志秋,陳振華,聶旭濤,等.基于流固熱耦合低溫風洞擴散段熱力學特性分析.實驗流體力學,2016,30(6):18-25.
張志秋(1992-),男,湖北黃岡人,碩士研究生。研究方向:風洞設計。通信地址:四川省綿陽市二環路南段6號14信箱402分箱(621000)。E-mail:zhiqiuz@gmail.com