身管精鍛過程跨尺度多晶體塑性有限元模擬與織構預測

徐笑，樊黎霞，王亞平，董曉彬

（南京理工大學機械工程學院，江蘇南京210094）

身管精鍛過程跨尺度多晶體塑性有限元模擬與織構預測

徐笑，樊黎霞，王亞平，董曉彬

（南京理工大學機械工程學院，江蘇南京210094）

徑向精鍛是一種先進的身管制造方法。身管的多晶體材料經鍛打后，晶粒會發生擇優取向形成織構，在宏觀上產生各向異性現象，這是普通有限元方法所無法模擬的。利用數學軟件Matlab和程序設計語言Python聯合編程并結合有限元軟件ABAQUS，建立身管精鍛宏觀模型和細觀多晶體模型；應用跨尺度方法，開發程序實現了宏觀力學至微觀力學邊界條件的繼承，使研究的晶粒尺寸達到了真實尺寸。分析微觀多晶體模型的晶粒取向變化，使用Matlab編程畫出鍛后多晶體材料的極圖和取向分布函數圖，并用X射線衍射試驗驗證了晶粒取向模擬結果的正確性。建立了將身管徑向精鍛從宏觀鍛造到細觀織構聯系起來的一套研究方法，為預測鍛后身管性能變化提供了思路。

兵器科學與技術；身管徑向鍛造；有限元分析；晶體塑性模型；織構

0 引言

身管是自動武器最重要的部件之一，在射擊中承受著高溫、高壓的火藥氣體的沖刷和彈丸的劇烈摩擦。與傳統身管制造工藝相比，徑向精鍛技術自動化程度高，提高了線膛和彈膛的同軸度，降低了成本并且提高了壽命，很大程度提高了自動武器的性能。有限元分析方法是身管徑向精鍛研究的一種常用方法。Tszeng等［1］最早建立了徑向鍛造的軸對稱模型。Domblesky等［2］使用ABAQUS有限元軟件對徑向鍛造過程進行了數值模擬，并且研究了鍛件進給速度和鍛件不同點的等效塑性應變對其溫度的影響。在國內，樊黎霞等［3］使用ABAQSU有限元軟件模擬了身管的徑向鍛造過程，對徑向鍛壓力及應力進行了數值分析，并用試驗驗證了模擬結果。劉力力等［4］建立了身管徑向鍛造及回彈過程的軸對稱模型，并分析了不同工藝參數對鍛后殘余應力分布的影響。

大量研究表明，在金屬塑性變形中晶體取向對變形和斷裂等行為是敏感的。變形后多晶體的擇優取向稱為織構。Bunge［5］提出織構影響材料性能的比例達到了20%～50%之高。這就為研究身管鍛后性能提供了新思路，以微觀織構的研究預測宏觀精鍛產品的鍛后性能。在板材沖壓領域已有類似先進的研究，馬全倉等［6］根據取向、織構和簡單的滑移理論預測了3104鋁合金板的初始塑性應變比，并且計算結果與實測值變化趨勢一致。陳攀宇等［7］利用Deform-3D軟件模擬和試驗研究AZ31鎂合金管軋制過程，并使用金相顯微鏡對軋制前后織構進行分析，表明多道次軋制使得晶粒再結晶和細化。晶體塑性有限元能夠模擬材料的各向異性行為，此方法更接近材料的本質特性。Balasubrmanian等［8］研究了多晶體塑性有限元在優化深沖工藝的可行性。Zhou等［9］用晶體塑性有限元研究了面心立方材料深沖制耳與織構的關系。在國內，李大永等［10］以軋制鋁板為研究對象，利用塑性有限元方法對方盒沖壓成形和織構演化進行了數值模擬。

跨尺度分析方法作為解決宏觀和微觀耦合問題的新方法在多晶體塑性領域已經得到廣泛關注。郭一娜等［11］提出一種宏觀和微觀尺度結合的建模方法來研究鑄態42CrMo環件的鑄輾復合成形工藝。在國外，Bronkhorst等［12］在研究材料受大變形和不同熱力的微觀結構的變化時，建立了鉭的基于狀態變量的宏觀內部連續模型，并提出了一種多尺度的多晶體塑性模型，對比兩種模型在進行剪切實驗時剪切區的VonMises應力，等效塑性應變等。在國內，杜鳳山等［13］針對Cr5鍛件鍛壓過程中晶體間的滑移進行了多尺度有限元仿真，得到了宏觀滑移趨勢和細觀滑移趨勢的對比。王國春［14］利用有限元重合網格方法開發程序，并對帶圓孔的金屬板等進行多尺度應力分析。

本文將以材料的本質屬性入手研究身管徑向精鍛，結合多晶體塑性有限元方法和跨尺度方法，建立一種將宏觀身管精鍛過程和微觀織構聯系起來的建模分析方法，為研究身管徑向精鍛的織構與宏觀性能的定量關系打下基礎。結合前人的理論框架和經驗，采用非線性分析功能強大的ABAQUS有限元軟件和Matlab程序語言、Python腳本編程聯合建模，分別建立身管精鍛過程的軸對稱宏觀模型和微觀多晶體有限元模型，利用線性插值理論實現宏觀至微觀邊界條件的繼承。最后，編譯程序提取所需晶粒取向繪制極圖和取向分布函數（ODF）圖并與X射線衍射（XRD）試驗所測的實際圖形作對比，驗證本方法的可行性。

1 精鍛過程宏觀仿真模擬

身管的徑向精鍛加工是一個較為復雜的金屬材料塑性流動過程。毛坯的底部由精鍛機的夾頭加持，并給予其一定的夾持力以實現軸向進給。毛坯的另一端與擋塊配合。毛坯四周對稱分布了4個錘頭，加工時4個錘頭對毛坯進行高頻率的鍛打。毛坯內部置有芯棒，芯棒具有身管線膛和彈膛的特征。為了保證圓形的橫截面和較高的表面光滑度，毛坯與芯棒在兩次鍛打間隙一起做旋轉進給。

錘頭的速度對材料應變率及再結晶有較大影響，在精鍛加工中，此速度為一諧函數，如（1）式所示:

式中:v為錘頭速度（mm/s）；n為鍛打頻率（次/min）；h為距上止點距離（mm）；s為錘頭全行程（mm）。

1.1 軸對稱有限元模型建立

運用ABAQUS對身管的徑向鍛造過程建立有限元模型。身管徑向鍛造在忽略內部膛線、身管旋轉軸向進給和錘頭的閉合間隙的情況下，是一個軸對稱的加工過程建立身管徑向鍛造的軸對稱模型。重要部件幾何尺寸如表1所示，身管材料力學性能參數如表2所示。

表1 身管毛坯幾何尺寸Tab.1 Geometric sizes of barrel workpiece

表2 材料力學性能Tab.2 Mechanical properties of material

身管的整個鍛造過程為冷精鍛，擋塊、錘頭和芯棒視為剛體。身管毛坯材料為30SiMn2MoVA，對其作網格劃分，采用4節點軸對稱縮減積分單元CAX4R，力學性能使用應變硬化模型σ= 1 225ε0.065，σ為流變應力，ε為應變，所建立宏觀模型圖1所示。

圖1 宏觀軸對稱有限元模型Fig.1 Macroscopic axisymmetric FEA model

1.2 仿真結果

在ABAQUS/CAE中設定上述參數后，于Job模塊中創建分析并求解，得到宏觀模型的Mises等效應力分布云圖（見圖2）。鍛壓力是精鍛實際生產中非常重要的一個工藝參數，加工時不能超過機器的最大鍛壓力值，所以在圖2中選取內表面應力較大的一點作為研究對象。

圖2 Mises應力布云圖Fig.2 Mises stress nephogram

2 微觀尺度上多晶體材料建模仿真

微觀尺度多晶體建模采用ABAQUS/CAE建模與INP文件建模兩種方法結合的方式。INP文件包含了對整個模型的完整描述，在前處理器和求解器之間建立了一座數據橋梁。具體建模方法分為4步:1）建立與宏觀模型所選單元尺寸相同的模型，劃分網格，設置分析步等重要信息；2）使用Voronoi方法進行晶粒劃分；3）將宏觀模型中得到的節點位移信息通過程序進行插值計算，得到新的微觀模型上4條邊界上所有節點的位移作為新的邊界條件；4）修改UMAT子程序信息，于求解器求解時提交。

圖3 微觀有限元模型建立流程Fig.3 Microscopic FEA modeling process

2.1 微觀基礎模型的建立

微觀有限元模型的建立流程如圖3所示。選取一個在圖3（a）中編號為3 366的單元進行跨尺度分析，即為圖3（b）的標識出的單元。使用ABAQUS/ CAE建立與宏觀模型中與此單元尺寸完全一致的模型。所建模型單元尺寸為250 μm×249 μm，共劃分為15 625個單元，即為圖3（c）。為了達到完全模擬在宏觀模型中的受力情況，單元采用與宏觀模型一樣的CAX4R單元。分析步等建模信息也與宏觀模型設置同樣的參數。

假設在鍛造前，模型內的晶粒大小相似、取向隨機并且沒有取向擇優的織構。這與Voronoi圖形生成方法隨機、均勻的思想相符合。所以采用該方法劃分晶粒。使用Matlab軟件編譯程序，讀取微觀模型的基本信息，如網格劃分和節點坐標等。Voronoi圖形的結構受種子點的數量及坐標的影響非常大。為了使種子點的信息受均勻性的控制，在程序中編寫控制函數，并調用Voronoi函數生成晶粒的拓撲信息，再使用inpolygon函數循環判斷出每個單元中心點屬于哪個晶粒，最后生成晶粒分組信息，并修改入模型的INPUT文件中。模型中的晶粒取向用米勒指數法來表示。

圖4為微觀多晶體模型建模的過程示意圖，右圖已經成功劃分晶粒，其含有258個晶粒，平均晶粒尺寸為15.5 μm，符合金屬晶粒的真實大小。

圖4 微觀多晶體有限元模型Fig.4 Microscopic polycrystalline FEA model

圖5 單元頂點位移Fig.5 Displacement of element's vertex

圖6 矩形單元廣義坐標Fig.6 Rectangular element generalized coordinates

2.2 跨尺度邊界條件的繼承

為使宏觀模型計算過程中該單元的邊界條件繼承至微觀模型中，使用有限元線性插值理論，并結合材料內部的位移連續理論，進行宏觀至微觀跨尺度邊界條件的繼承。

從宏觀模型計算結果中可以得到該單元的4個頂點在整個精鍛過程中隨時間變化的位移曲線。圖5為單元第一頂點的x、y方向位移圖。x軸的正方向為軸向進給方向，y軸的正方向為毛坯的徑向。由于冷鍛工藝成形是利用金屬的塑性流動來實現的，當所取單元位置受到鍛打后，金屬會向進給方向的反方向流動，開始階段流動速度大于進給速度，所以在開始階段y軸會出現負值，隨著所取單元通過軸向進給離鍛打區域越來越遠，其位移主要受軸向進給控制，所以位移值斜率會趨近于進給速度。

微觀建模型共1 000個邊界節點。使用二維4節點線性單元形函數和位移函數聯合計算得出任意節點的位移曲線。首先提取出1 000個邊界點的坐標值，如圖6所示，先將其化成廣義坐標，再代入（2）式～（6）式中，利用得到的廣義坐標值計算每個點的4個原單插值函數Ni，再將4個頂點的位移值ui和4個單元插值函數Ni代入（7）式位移函數即可得到所需邊界點的相應位移值｛r｝，具體公式為

式中:i（ξi，ηi）分別為1（1，1）、2（-1，1）、3（-1，-1）、4（1，-1）；x0和y0為自然坐標系的原點坐標；a、b為圖6中矩形單元的長和寬。

采用Python腳本語言計算插值結果，程序生成INPUT格式的文本文件，將其修改入模型文件中，最后達到添加邊界條件的目的。首先，讀取微觀模型的OBD文件，于根裝配中提取模型邊界的節點信息，包括節點編號和坐標，并將每個節點定義為一個單獨的集合，最后以INPUT格式輸出。具體格式為

讀取宏觀模型的OBD結果文件，使用（1）式的計算方法得到1 000個節點的隨時間變化的位移曲線。在ABAQUS中添加約束，為上一步中定義的單個節點集合設置位移載荷約束，位移大小設置為1，再將求得的任意節點的時間位移數據以Amplitude方式加于約束中。具體格式為

如上例所示，約束類型為位移/旋轉型，施加目標為set-1集合，即前面所創建的邊界的單節點集合，后面的1代表約束方向為x方向，最后的一個1代表位移值為1.此約束受Amplitude控制，上例所示的Amplitude模塊中左邊一列為時間，單位為s，右邊一列為位移，單位為mm.將上述生成的集合信息、約束信息和振幅信息添加入微觀模型的INPUT文件中，完整的微觀材料模型就建成了。

當然，需要指出的是，目前沒有任何跡象表明有火星微生物在那里“游泳”。事實上，這樣環境也不適宜生命的存在。火星南極冰層下的水溫預計在華氏零下90度，遠低于水的冰點。

2.3 運用UMAT實現多晶體自定義材料屬性

微觀多晶體模型的塑性本構關系使用ABAQUS中的子程序UMAT來實現。本文采用Huang［15］的UMAT程序框架，硬化模型采用Asaro模型。模型中的某些參數如初始屈服時的硬化模量h0、應變率敏感系數n、參考剪切應變率γ·0等材料系數很難準確測量到，即使個別方法如潛硬化實驗等也只能得到單晶的個別材料參數。現存的文獻中也沒有人做過本文所用材料30SiMn2MoVA的這些材料屬性的研究。所以本文采用擬合實驗所得應力-應變曲線的方法來獲得這些材料系數。

3 數值模擬結果與試驗驗證

3.1 微觀模型模擬結果

進行劃分晶粒、材料本構關系賦予和邊界條件的繼承后就能得到微觀多晶體塑性有限元模型。提交解算器計算后得到微觀模型的真實應變分布如圖7所示。

圖7 應變分布云圖Fig.7 Strain nephogram

從圖7中可以看出，模型內部真實應變分布嚴重的不均勻，在應變帶之間交替分布著大小不同的小應變區域和應變集中區域，這是由于不同晶粒的晶粒取向和滑移系的差別在塑性變形過程中引發應變的累積，從而使得相鄰的晶粒的應變也可能會有很大的差別，最后導致應變及其不均勻的分布。有些晶粒內部應變大小差別在5倍以上，這是傳統的有限元模擬達不到的效果。應變的集中不僅僅出現在晶界附近，也能出現于晶粒的內部。由于晶粒本身的各向異性、幾何形狀和晶粒間的變形協調導致晶粒的形狀變得更加的不規則。

3.2 極圖繪制與試驗驗證

晶體學一般采用極圖和ODF圖來描述織構。身管的材料為30SiMn2MoVA，屬于體心立方（BCC）晶體，具有48個滑移系，研究發現在冷塑性變形中只有一部分滑移系有大量的開動，主要開動的滑移系分別為12個｛110｝〈111〉滑移系和12個｛112｝〈111〉滑移系。編譯Python腳本程序提取出結果文件中上述24個滑移系的密勒指數信息，并將其轉換成歐拉角（φ1，φ，φ2）格式。使用Matlab的插件MTEX編譯程序，讀取上述歐拉角信息，繪制可以較好表達變形中晶粒取向變化的｛100｝方向極圖和以φ1為視角的ODF圖。

為了驗證模擬結果的正確性，在5.8 mm口徑鍛后身管上切割了表面積為10 mm×10 mm正方形的試樣進行了XRD試驗。該試驗從材料、身管尺寸及精鍛加工參數等采用與本文一致的工藝參數。XRD分析是利用晶體形成的XRD，對物質進行內部原子在空間分布狀況的結構分析方法。本試驗采用德國Bruker公司產D8 Advance XRD儀。試驗樣品采用在鍛后身管中采用線切割方法得到邊長10 mm的塊狀材料，并在干燥環境下進行多次打磨，如圖8所示。試驗完成后得到XRD圖譜和3個方向上的不完全極圖兩種數據文件。

圖8 XRD試驗試樣Fig.8 XRD experimental samples

使用Jade5.0分析衍射圖譜得到物相信息，再配合不完全極圖3個不完全極圖導入Multex軟件處理得到相同方向的極圖和ODF圖。

圖9 BCC金屬取向線Fig.9 Orientations of BCC metal

對比觀察模擬結果和試驗結果所得到的極圖（見圖10）和ODF圖（見圖11），其強度較大的織構分布位置和對稱性都較高的一致性。BCC金屬變形纖維織構主要的織構類型有｛111｝〈UVW〉和｛HKL｝〈110〉兩類，而這些織構大都處在α取向線和γ取向線。從DOF圖中可以看出，在BCC中重要取向線α取向線上（0°，54.7°，45°）附近以及γ取向線上（90°，54.7°，45°）處有織構集中，如黑色十字標記所示。所以判定這兩個位置有織構生成，將這兩個位置的織構轉化為密勒指數法表示，分別為｛111｝〈110〉和｛111｝〈112〉。

圖10 極圖比較Fig.10 Comparison of pole figures

圖11 ODF圖比較Fig.11 Comparison of ODF figures

將模擬結果通過Horta加權法計算可以得到各｛HKL｝的軸密度，同樣XRD試驗測得各｛HKL｝的衍射強度也可計算得到軸密度。將軸密度使用地質等高繪圖軟件Winsurf繪制反極圖后可以得到模擬和試驗中較強的織構強度如表3和表4所示。

表3 模擬所得較強織構的織構強度Tab.3 Simulated texture intensity of strong texture

表4 模擬所得較強織構的織構強度Tab.4 Experimental texture intensity of strong texture

可以得到30SiMn2MoVA材料身管徑向精鍛鍛后織構成分為｛111｝〈112〉和｛111｝〈110〉，其中前者為主要織構成分。

4 結論

利用多晶體塑性有限元方法建立了身管精鍛的宏觀至微觀多尺度有限元模型，使用線性插值理論編譯程序實現了宏觀和微觀兩個模型邊界條件的統一，使研究的晶粒尺寸達到了實際大小。繪制極圖和ODF圖并與試驗結果對比，將宏觀鍛造和微觀織構聯系起來。為以后研究宏觀性能和微觀織構的關系提供了一套研究方法。

1）基于多晶體塑性有限元方法和跨尺度方法，開發程序，建立從宏觀到微觀的跨尺度身管多晶體塑性有限元模型，可以更好地反映身管精鍛中的各向異性現象，使身管精鍛的數值模擬更加符合實際。

2）基于所建立的跨尺度模型分析微觀織構，提取分析數據繪制極圖和ODF圖并與XRD試驗所得圖形對比具有較好的一致性。

3）通過多晶體塑性有限元法和跨尺度方法的結合，不僅可以更加準確地模擬金屬宏觀變形，而且能夠預測塑性變形中的織構變化。

（References）

［1］ Tszeng T C，Kobayashi S.Determination of residual stresses in radial forging，manufacturing processes simulation，PED-20［R］. Berkeley，US:Department of Mechanical Engineering，University of California，1986:31-45.

［2］ Domblesky J P，Shivpuri R，Painter B.Application of the finite-element method to the radial forging of large diameter tubes［J］. Journal of Materials Processing Technology，1995，49（1）:57-74.

［3］ 樊黎霞，劉力力，劉慶東，等.身管線膛精鍛加工過程的數值分析［J］.兵工學報，2009，30（8）:1098-1102. FAN Li-xia，LIU Li-1i，LIU Qing-dong，et al.The numerical simulation of the precision radial forging processing of the gun barrel［J］.Acta Armamentarii，2009，30（8）:1098-1102.（in Chinese）

［4］ 劉力力，樊黎霞，董雪花.盲孔法測定徑向鍛造身管表面殘余應力方法研究［J］.兵工學報，2012，33（6）:712-717. LIU Li-li，FAN Li-xia，DONG Xue-hua.Study on measurement of residual stresses in forged barrel surface［J］.Acta Armamentarii，2012，33（6）:712-717.（in Chinese）

［5］ Bunge H J.Texture analysis in materials science mathematical methods［M］.London，UK:Butterworth Heinemann Ltd，1982.

［6］ 馬全倉，毛衛民，馮惠平.3104深沖鋁板織構對初始R值的影響［J］.北京科技大學學報，2004，26（1）:79-81. MA Quan-cang，MAO Wei-min，FENG Hui-ping.Influence of texture on initial R-value of deep drawing 3104 aluminum sheet［J］.Journal of University of Science and Technology Beijing，2004，26（1）:79-81.（in Chinese）

［7］ 陳攀宇，李彩霞，李立州，等.基于Deform-3D的鎂合金管軋制成形數值模擬與試驗研究［J］.輕合金加工技術，2015，43（7）:33-37. CHEN Pan-yu，LI Cai-xia，LI Li-zhou，et al.Numerical simulation and experimental study of roll forming of magnesium alloy tubes based on Deform-3D［J］.Light Alloy Fabrication Technology，2015，43（7）:33-37.（in Chinese）

［8］ Balasubramanian S，Anand L.Single crystal and polycrystal elasto-viscoplasticity:application to earing in cup drawing of FCC materials［J］.Computational Mechanics，1996，17（4）:209-225.［9］ Zhou Y，Jonas J J，Savoie J，et al.Effect of texture on earing in FCC metals:finite element simulations［J］.International Journal of Plasticity，1998，14（1/2/3）:113-138.

［10］ 李大永，張少睿，彭穎紅，等.板材沖壓成形的晶體塑性有限元模擬［J］.機械工程學報，2008，44（1）:190-194. LI Da-yong，ZHANG Shao-rui，PENG Ying-hong，et al.Finite element simulation of sheet metal stamping with polycrystalline plasticity［J］.Journal of Mechanical Engineering，2008，44（1）: 190-194.（in Chinese）

［11］ 郭一娜，丁雙鳳，李永堂，等.42CrMo環件鑄輾復合成形的跨尺度建模［J］.機械工程學報，2014，50（14）:81-88. GUO Yi-na，DING Shuang-feng，LI Yong-tang，et al.Multiscale modeling for 42CrMo ring during blank-casting and rolling compound forming process［J］.Journal of Mechanical Engineering，2014，50（14）:81-88.（in Chinese）

［12］ Bronkhorst C A，Maudlin P J，Gray G T，et al.Accounting for microstructure in large deformation models of polycrystalline metallic materials［M］∥Ghosh S，Dimiduk D.Computational Methods for Microstructure Property Relationships.NY，US: Springer，2010:239-276.

［13］ 杜鳳山，李源，王珉，等.鍛壓過程宏-細觀跨尺度仿真研究［J］.計算力學學報，2014，31（6）:799-810. DU Feng-shan，LI Yuan，WANG Min，et al.Macro-meso crossscale simulation of forging process［J］.Chinese Journal of Computational Mechanics.2014，31（6）:799-810.（in Chinese）

［14］ 王國春.基于有限元重合網格法的結構多尺度分析研究［D］.新疆:新疆大學，2006. WANG Guo-chun.Study on multi-scale analysis of structure based on the finite element mesh superposition method［D］. Xinjiang:Xinjiang University，2006.（in Chinese）

［15］ Huang Y.A user-material subroutine incorporating single crystal plasticity in the ABAQUS finite element program［R］.Cambridge，MA，US:Division of Engineering and Applied Sciences，Harvard University，1991:2-9.

Trans-scale Polycrystalline Finite Element Simulation of Radial Forging Process for Barrel and Prediction of Texture

XU Xiao，FAN Li-xia，WANG Ya-ping，DONG Xiao-bin
（School of Mechanical Engineering，Nanjing University of Science and Technology，Nanjing 210094，Jiangsu，China）

Radial forging is an advanced manufacturing method for barrel.After the polycrystalline material of barrel is forged，the preferred orientation of the materials grain forms texture which presents the anisotropy phenomenon macroscopically.This phenomenon can not be simulated using the common finite element analysis（FEA）method.A macroscopic barrel radial forging FEA model and a microscopic polycrystalline plasticity FEA model are established based on ABAQUS code，the mathematical software Matlab and the programming language Python.The multi-scale method is used to program the boundary conditions inheriting from macro-mechanics state to meso-mechanics state which make the size of grain as reality.The pole figures of forged polycrystalline material and the orientation distribution function（ODF）figures are drawn by Matlab，which are used to analyze the orientation change of grain in the microscopic polycrystalline plasticity FEA model.The XRD experiment is used to validate the accuracy of the simulated results.A research method to connect macroscopic forging process and microscopic texture is presented，which provides a new way to predict the mechanical property of barrel after forging.

ordnance science and technology；radial forging of barrel；finite element analysis；crystalline plasticity model；texture

TJ05

A

1000-1093（2016）07-1180-07

10.3969/j.issn.1000-1093.2016.07.004

2015-08-06

國家自然科學基金項目（51575279）

徐笑（1989—），男，碩士研究生。E-mail:xuxiaoam@sina.com；

樊黎霞（1965—），女，教授，碩士生導師。E-mail:fanlixia151@njust.edu.cn