楊 揚(yáng),徐 緋,張岳青,湯忠斌

(1.西北工業(yè)大學(xué)航空學(xué)院結(jié)構(gòu)工程系，陜西 西安 710072; 2.北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京 100081)

?

平紋編織C/SiC復(fù)合材料低速沖擊數(shù)值模擬

楊 揚(yáng)1,2,徐 緋1,張岳青1,湯忠斌1

(1.西北工業(yè)大學(xué)航空學(xué)院結(jié)構(gòu)工程系，陜西 西安 710072; 2.北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京 100081)

首先，基于空氣炮裝置進(jìn)行了2D-C/SiC薄板在沖擊速度為79～219 m/s范圍內(nèi)的低速沖擊實(shí)驗(yàn)，對碎片云團(tuán)發(fā)展過程進(jìn)行高速攝影記錄；其次，基于Autodyn軟件正交各向異性復(fù)合材料模型，推導(dǎo)2D-C/SiC材料相關(guān)參數(shù)；選取SPH求解器建立二維計算模型，對實(shí)驗(yàn)工況進(jìn)行數(shù)值模擬，并基于碎片云結(jié)構(gòu)、B掃描檢測結(jié)果和碎片云軸向發(fā)展速度驗(yàn)證了該模型可以很好地描述C/SiC材料在沖擊載荷作用下的脆性特征和軟化行為。最后，基于數(shù)值模擬結(jié)果推導(dǎo)得出了鋼彈丸沖擊C/SiC材料的極限侵徹深度預(yù)測公式。

固體力學(xué)；C/SiC材料沖擊特性；SPH；熱防護(hù)結(jié)構(gòu)；碎片云；高速攝影

高超聲速飛行器熱防護(hù)系統(tǒng)(thermal protection system, TPS)的結(jié)構(gòu)材料必須具有輕質(zhì)、耐高溫、環(huán)境穩(wěn)定等特點(diǎn)，以滿足長航程安全服役的要求。近年來，陶瓷基復(fù)合材料在飛行器大面積熱防護(hù)領(lǐng)域應(yīng)用廣泛，這類材料可以比傳統(tǒng)的金屬TPS減重50%，同時可以延長使用壽命，降低制造成本[1]。

基于化學(xué)氣相沉積技術(shù)(chemical vapor infiltration, CVI)[2]制造的碳纖維增韌碳化硅陶瓷基復(fù)合材料(C/SiCS)以其低密度、抗氧化[3]、耐高溫[4]等優(yōu)異特性逐漸得到重視，目前已經(jīng)成功應(yīng)用于飛行器的鼻錐和機(jī)翼前緣的防熱結(jié)構(gòu)中。然而，與金屬材料不同，陶瓷基材料不具備產(chǎn)生塑性變形的能力，材料的內(nèi)部分層、基體裂紋、纖維斷裂等破壞是其吸收沖擊能量的主要方式，這些損傷會降低材料的力學(xué)性能，從而嚴(yán)重威脅材料的服役安全。因此，為了飛行器內(nèi)部結(jié)構(gòu)及設(shè)備在服役過程中的安全性考慮，這類結(jié)構(gòu)材料在沖擊載荷作用下的破壞行為及其數(shù)值預(yù)測一直是眾多學(xué)者關(guān)注的熱點(diǎn)。

目前，針對平紋編織C/SiC(2D-C/SiC)復(fù)合材料的研究并不多見，已有的相關(guān)文獻(xiàn)多數(shù)僅限于對材料簡單力學(xué)性能的實(shí)驗(yàn)研究[5-6]。楊揚(yáng)等[7]基于電炮加載裝置驅(qū)動Mylar飛片完成了沖擊速度在3 400～9 300 m/s范圍內(nèi)2D-C/SiC材料的超高速撞擊實(shí)驗(yàn)，并對材料損傷演變和破壞特征進(jìn)行了分析研究。然而，關(guān)于材料在低速沖擊載荷作用下的實(shí)驗(yàn)及其數(shù)值模擬目前尚未發(fā)現(xiàn)。

材料的各向異性特征使得在對其進(jìn)行模擬計算時面臨著巨大的挑戰(zhàn)，R.A.Clegg等[8]提出的正交各向異性模型可以較好地描述復(fù)合材料在沖擊載荷作用下的力學(xué)響應(yīng)，已經(jīng)在Kevlar-Epoxy復(fù)合材料的沖擊模擬中得到了成功的應(yīng)用[9]。然而，該模型重點(diǎn)關(guān)注了材料的非線性應(yīng)變硬化行為，在模擬C/SiC這種典型的脆性材料時無法得到理想的數(shù)值結(jié)果。本文中在已完成的球形彈丸低速沖擊2D-C/SiC薄板的實(shí)驗(yàn)基礎(chǔ)上，引入材料的各向異性本構(gòu)模型，基于Autodyn商用軟件對2D-C/SiC復(fù)合材料在鋼球彈丸沖擊作用下的數(shù)值模擬問題進(jìn)行研究。

1 2D-C/SiC材料低速沖擊實(shí)驗(yàn)

基于空氣炮加載裝置對2D-C/SiC復(fù)合材料薄板在鋼球撞擊作用下的力學(xué)響應(yīng)問題進(jìn)行了實(shí)驗(yàn)研究。在實(shí)驗(yàn)過程中，使用高速攝影技術(shù)記錄了2D-C/SiC薄板受撞擊后碎片云團(tuán)的發(fā)展過程。

1.1 實(shí)驗(yàn)材料與計劃

實(shí)驗(yàn)中選用的2D-C/SiC復(fù)合材料由西北工業(yè)大學(xué)超高溫復(fù)合材料重點(diǎn)實(shí)驗(yàn)室研制。首先，T300碳纖維正交編織成碳布預(yù)制體，在低壓爐內(nèi)基于CVI技術(shù)沉積6次SiC基體，得到孔隙率約15%、體積分?jǐn)?shù)約45%的2D-C/SiC材料，而后基于CVI沉積2次SiC涂層。最終加工成型后，切割為大小115 mm×115 mm、厚度3 mm、密度約2.11 g/cm3的正方形薄板作為目標(biāo)靶板。實(shí)驗(yàn)中選用5、6 mm兩種直徑的普通鋼球作為沖擊彈丸。具體實(shí)驗(yàn)計劃及主要結(jié)果列于表1中，表中D表示彈丸直徑，vp表示彈丸的沖擊速度，E表示彈丸的動能。

表1 低速沖擊實(shí)驗(yàn)計劃

1.2 碎片云團(tuán)結(jié)構(gòu)

限于實(shí)驗(yàn)條件，高速攝像機(jī)置于靶板斜后方位置，圖1給出了試樣8的高速攝影視頻截圖。由圖1可見，2D-C/SiC碎片云團(tuán)的發(fā)展過程大致分為2個階段：(1) 簇狀飛散(圖1(a))。彈丸接觸靶板時，在作用點(diǎn)局部區(qū)域迅速發(fā)生纖維斷裂、基體粉碎等多種形式破壞，大量固體顆粒成簇狀形式堆積飛出；(2) 兩區(qū)域形成(圖1(b)～(d))。簇狀顆粒群運(yùn)動發(fā)展，逐漸在碎片云團(tuán)中形成了界限清晰的兩區(qū)域結(jié)構(gòu)，即柱狀區(qū)和分散區(qū)。在柱狀區(qū)內(nèi)，顆粒分布密度較高、顆粒尺寸較小、運(yùn)動速度較快，并沿軸線向前運(yùn)動；在分散區(qū)內(nèi)，顆粒分布較為稀疏、顆粒尺寸較大、運(yùn)動速度較慢，且沿一定的飛散角輻射飛出。2D-C/SiC碎片云結(jié)構(gòu)與以往的金屬材料及樹脂基復(fù)合材料截然不同，是其脆性特征的典型表現(xiàn)。在數(shù)值模擬中，建立合理的本構(gòu)模型對材料特有的力學(xué)響應(yīng)特征進(jìn)行模擬計算是本文研究的重點(diǎn)。

圖1 碎片云團(tuán)發(fā)展過程Fig.1 Expansion of the debris cloud

2 耦合的各向異性材料模型

2.1 正交各向異性彈性模型

復(fù)合材料在沖擊載荷作用下的力學(xué)響應(yīng)過程十分復(fù)雜，通常是諸如體積響應(yīng)、材料各向異性等多種因素的耦合。在數(shù)值計算中，不僅球應(yīng)變會影響偏應(yīng)力的計算，偏應(yīng)變同樣會在球應(yīng)力的計算中起作用，通常需要面臨求解狀態(tài)方程和本構(gòu)關(guān)系這兩個高度耦合的子模型[10]。傳統(tǒng)的計算方法是將這兩個子模型獨(dú)立計算，然而這樣的過程較為繁瑣，C.E.Anderson等[11]首先提出了一種適用于求解各向異性材料中多重響應(yīng)耦合問題的計算方法。本文中也以此方法為基礎(chǔ)，建立如下的計算模型。

由于各向同性壓力P是3個正應(yīng)力的平均值的相反數(shù)，因此由上式可以推導(dǎo)得出：

(1)

從(1)式中可以清晰地看出，與各向同性材料中球應(yīng)力僅依賴于體積應(yīng)變不同，各向異性材料中偏應(yīng)變也會在一定程度上影響球應(yīng)力的計算。其中ε0一項的系數(shù)可稱為材料的等效體積模量，記為A。

2.2 正交各向異性失效模型

失效模型主要包括失效初始化和后失效響應(yīng)兩個階段。在失效初始化階段，采用最大應(yīng)力判斷準(zhǔn)則，即預(yù)先設(shè)置材料在3個拉伸方向和3個剪切方向的失效應(yīng)力，當(dāng)任一方向的應(yīng)力值達(dá)到失效應(yīng)力時，材料即在該方向失效，進(jìn)入后失效響應(yīng)階段。

對于材料的后失效響應(yīng)，需要在計算開始前按照如下準(zhǔn)則進(jìn)行相關(guān)設(shè)置：當(dāng)材料在任一方向進(jìn)入后失效階段時，材料在該方向的承載能力瞬間喪失，但在其他方向的承載能力保持不變；同時，材料在3個剪切方向的承載能力下降[8]。例如，假設(shè)1方向?yàn)楸“搴穸确较颍?、3方向?yàn)槊鎯?nèi)的2個正交方向。那么，當(dāng)1方向拉伸應(yīng)力大于拉伸失效應(yīng)力或12(或13)方向剪切應(yīng)力大于剪切失效應(yīng)力時，材料會發(fā)生分層破壞。此時，材料在1方向應(yīng)力立即置為0，不再具有承載能力，但在面內(nèi)的2個方向承載能力不變；同時由于分層破壞的發(fā)生，材料在剪切方向的承載能力也將有所下降，需要對剛度矩陣中的剪切剛度系數(shù)進(jìn)行退化，可對材料的剛度矩陣作如下修正(α為剛度退化因子)：

具體的材料參數(shù)及失效設(shè)置分別見表2和表3。

表2 正交彈性模型材料參數(shù)

表3 正交失效模型材料參數(shù)

3 數(shù)值模擬及分析

圖2 Autodyn計算模型Fig.2 Computational model in Autodyn

3.1 計算模型

由于正撞擊的對稱性，基于Autodyn軟件建立如圖2所示的二維半模型計算。2D-C/SiC靶板選用SPH求解器；鋼彈丸steel 4340基于FEM求解器計算；彈丸與靶板之間設(shè)置接觸，接觸間距0.01 mm；在靶板邊緣3 mm×4 mm范圍內(nèi)施加全局固定邊界，即限定此區(qū)域內(nèi)的粒子在x和y方向速度恒為0。

下面通過對碎片云結(jié)構(gòu)和典型破壞形式等實(shí)驗(yàn)結(jié)果的數(shù)值驗(yàn)證，證明該模型對描述2D-C/SiC材料在沖擊載荷作用下力學(xué)響應(yīng)的合理性和可行性。

3.2 數(shù)值合理性的實(shí)驗(yàn)對比

3.2.1 碎片云團(tuán)結(jié)構(gòu)的數(shù)值驗(yàn)證

圖3 (b)、(c)分別給出了試樣8的數(shù)值模擬結(jié)果與對應(yīng)的速度矢量圖。通過對比可以看出，實(shí)驗(yàn)結(jié)果和數(shù)值模擬結(jié)果吻合較好，具有以下相似特征：(1)碎片云團(tuán)均呈現(xiàn)明顯的雙區(qū)域結(jié)構(gòu)，即均存在沿沖擊軸線方向的柱狀區(qū)和軸線兩側(cè)的飛散區(qū)，這是2D-C/SiC在沖擊載荷作用下的特殊響應(yīng)特征；(2)從速度矢量圖可以看出，柱狀區(qū)粒子運(yùn)動方向均為沿軸線方向；而飛散區(qū)內(nèi)的粒子速度方向與軸線有大小不同的夾角。這與實(shí)驗(yàn)結(jié)果是完全一致的，表明該模型可以較好地表征2D-C/SiC在彈丸撞擊下碎片云團(tuán)的結(jié)構(gòu)特征。

圖3 試樣8計算結(jié)果Fig.3 Simulation results of sample 8

3.2.2 B掃描結(jié)果的數(shù)值驗(yàn)證

B掃描結(jié)果顯示的實(shí)際上是試樣內(nèi)部損傷的二維截面圖。圖4(a)給出試樣8的B掃結(jié)果，其中虛線區(qū)域的黃色部分表示試樣；圖4(b)～(c)給出工況8靶板擴(kuò)展穩(wěn)定后內(nèi)部損傷的模擬計算結(jié)果。

圖4 B掃描結(jié)果與模擬結(jié)果對比圖Fig.4 Comparison between B-scan and simulation results

從圖4看出:(1)數(shù)值模擬與B掃描顯示的穿孔形狀十分相似，穿孔直徑近似為彈丸直徑(6 mm)；(2)從B掃結(jié)果看，靶材背面顯示黃色較正面偏淡，這是由于材料在靠近背面的區(qū)域損傷嚴(yán)重，導(dǎo)致超聲波信號無法完全被探頭接收。相似的現(xiàn)象可以在數(shù)值模擬中顯示得更加清晰，圖4(c)給出了靶材穿孔區(qū)域的局部放大圖。由圖可見，與正面相比，靶材背面存在明顯的材料分層與剝落損傷，這與實(shí)驗(yàn)結(jié)果比較吻合。上述相似性也表明該模型可以較好地表征2D-C/SiC在彈丸撞擊下的破壞特征。

3.2.3 碎片云軸向速度的數(shù)值驗(yàn)證

碎片云的軸向速度是表征碎片云團(tuán)發(fā)展?fàn)顟B(tài)的重要指標(biāo)。實(shí)驗(yàn)中，在靶板表面標(biāo)記了刻度，據(jù)此可以根據(jù)高速攝影記錄計算出碎片云團(tuán)柱狀區(qū)前端的運(yùn)動速度，而將此速度視為碎片云軸向運(yùn)動速度vaxis。同樣，在數(shù)值模擬中，取靶板自由面與沖擊軸線的交點(diǎn)為特征點(diǎn)近似柱狀區(qū)端點(diǎn)。選取靶板穿透工況中沖擊速度相差較大的4種情況進(jìn)行驗(yàn)證，圖5給出了不同沖擊速度vaxis下隨時間的變化歷程。

圖5 不同沖擊速度下特征點(diǎn)軸向速度歷程Fig.5 Axial velocity histories under different impact velocities

圖6 特征點(diǎn)軸向速度對比結(jié)果Fig.6 Comparison of the axial velocities of simulation and experiment

從圖5可以看出，在彈丸接觸靶板后的短時間內(nèi)，特征點(diǎn)的軸向速度迅速增加并逐漸趨于穩(wěn)定。將特征點(diǎn)最終的穩(wěn)定速度視為對vaxis的近似，圖6給出了實(shí)驗(yàn)數(shù)據(jù)與數(shù)值計算的對比結(jié)果。可以看出，二者的吻合度較好，并且隨著沖擊速度的增高，相對誤差逐漸減小，其中沖擊速度219 m/s時二者的相對誤差僅為2.25%，表明該模型可以較好地表征2D-C/SiC在沖擊載荷作用下碎片云團(tuán)的速度特性。

4 極限侵徹深度預(yù)測模型

盡管碎片云團(tuán)中柱狀區(qū)內(nèi)的粒子運(yùn)動速度較高，但由于顆粒質(zhì)量較小，對機(jī)體的損傷程度有限。碎片云團(tuán)中對機(jī)體威脅較大的仍是穿透靶板的彈丸，因此可以通過對彈丸剩余速度的分析研究，建立合理有效的預(yù)測公式評估C/SiC靶板的抗沖擊性能。以上通過3方面對比驗(yàn)證了文中使用的材料模型及相關(guān)參數(shù)的正確性和合理性，下面將基于此模型對彈丸侵徹下C/SiC靶板的抗沖擊特性進(jìn)行預(yù)估。

考慮同一彈丸(直徑6 mm)以不同沖擊速度對不同厚度靶板進(jìn)行撞擊。設(shè)置靶板厚度δ為1、2、3、4、5 mm，彈丸沖擊速度vi為50、100、150、200、250、300、350、400 m/s，計算得到彈丸剩余速度變化規(guī)律如圖7所示?？梢钥闯觯?1)對同一厚度的靶板，剩余速度隨著沖擊速度的減小不斷減??；而對相同沖擊速度，剩余速度會隨著靶板厚度的減小而不斷增加；(2)剩余速度在速度較低靶板較厚時出現(xiàn)了負(fù)值，這是由于彈丸在這些工況下未能穿透靶板，出現(xiàn)回彈現(xiàn)象造成的。

評估材料抗沖擊性能比較關(guān)心的重要參數(shù)是靶板的極限穿透速度，即上述曲面與剩余速度為0的平面的交線。然而，如果固定靶板厚度分別觀察剩余速度隨沖擊速度的變化趨勢(圖8)，可以看出在臨界穿透狀態(tài)附近(圖中虛線)，剩余速度變化曲線斜率陡然增高，說明在臨界穿透狀態(tài)附近剩余速度變化較劇烈；而當(dāng)沖擊速度高于極限穿透速度之后，彈丸剩余速度變化基本穩(wěn)定，呈準(zhǔn)線性增長趨勢。因此，如果直接對剩余速度曲面進(jìn)行擬合來求解極限穿透速度，會導(dǎo)致在臨界穿透狀態(tài)附近的計算誤差較大。

圖7 剩余速度隨沖擊速度與靶板厚度的變化規(guī)律Fig.7 Variation of the residual velocity with impact velocity and plate thickness

圖8 不同厚度靶板剩余速度變化圖Fig.8 Variation of residual velocity with impact velocity for different plate thicknesses

為了解決上述問題，采用分步擬合方法：(1)假設(shè)靶板無限厚，首先固定沖擊速度，分別擬合得到彈丸剩余速度隨靶板厚度的關(guān)系，進(jìn)而求解在此沖擊速度下彈丸的極限侵徹深度；(2)其次基于Poncelet型侵徹深度預(yù)測模型[12]，利用得到的極限侵徹深度T與彈丸沖擊速度vi進(jìn)行擬合，得到如下關(guān)系：

此時，便可以在已知靶板厚度的情況下，通過反解上式預(yù)測彈丸的極限穿透速度。同時，在工程設(shè)計中，如果工程人員已經(jīng)根據(jù)飛行器的服役環(huán)境預(yù)估出外來沖擊物的撞擊速度，那么同樣可以利用上式預(yù)估彈丸的極限侵徹深度，也就可以為防護(hù)板厚度參數(shù)的確定提供一定的參考。

5 結(jié) 論

基于Autodyn軟件選取正交各向異性材料模型，推導(dǎo)了相關(guān)的材料參數(shù)，并從碎片云團(tuán)結(jié)構(gòu)、B掃描內(nèi)部損傷形貌和碎片云軸向速度三個方面對比了實(shí)驗(yàn)與模擬結(jié)果，表明參數(shù)的合理性和正確性，較好地描述了2D-C/SiC在鋼球彈丸撞擊下的主要損傷特征和力學(xué)特性。在此基礎(chǔ)上，基于數(shù)值模擬計算結(jié)果推導(dǎo)得到了鋼球彈丸撞擊2D-C/SiC材料的極限侵徹深度預(yù)測公式，這可為TPS結(jié)構(gòu)設(shè)計提供一定的工程參考價值。

[1] 張立同,成來飛,徐永東.新型碳化硅陶瓷基復(fù)合材料的研究進(jìn)展[J].航空制造技術(shù),2003(1):24-32. Zhang Li-tong, Cheng Lai-fei, Xu Yong-dong. Progress in research work of new CMC-SiC[J]. Aeronautical Manufacturing Technology, 2003(1):24-32.

[2] Chawla K K. Ceramic matrix composite[M]. 2nd ed. USA: Kluwer Academic Publisher, 2003.

[3] Lamouroux F, Bertrand S, Pailer R, et al. Oxidation-resistant carbon-fiber-reinforced ceramic-matrix composites[J]. Composites Science and Technology, 1999,59(7):1073-1085． [4] 殷曉光.C/SiC陶瓷基復(fù)合材料的力學(xué)及高溫性能研究[D].北京:清華大學(xué),2011.

[5] Camus G, Guillaumat L, Baste S. Development of damage in a 2D woven C/SiC composite under mechanical loading: Ⅰ. Mechanical characterization[J]. Composites Science and Technology, 1996,56(12):1363-1372.

[6] Bouazzaoui R E, Baste S, Camus G. Development of damage in a 2D woven C/SiC composite under mechanical loading: Ⅱ. Ultrasonic characterization[J]. Composites Science and Technology, 1996,56(12):1373-1382.

[7] 楊揚(yáng)，徐緋，張岳青,等.二維平紋編織C/SiC復(fù)合材料的超高速碰撞實(shí)驗(yàn)[J].爆炸與沖擊，2013,33(2):156-162. Yang Yang, Xu Fei, Zhang Yue-qing, et al. Hypervelocity impact experiment on two-dimensional plain-woven C/SiC composites[J]. Explosion and Shock Waves, 2013,33(2):156-162.

[8] Clegg R A, Hayhurst C J, Leahy J G, et al. Application of a coupled anisotropic material model to high velocity impact of composite textile armour[C]∥18th International Symposium and Exhibition on Ballistics, USA: San Antonio, 1999.

[9] Hayhurst C J, Hiermaier S J, Clegg R A, et al. Development of material models for Nextel and Kevlar-Epoxy for high pressures and strain rates[J]. International Journal of Impact Engineering, 1999,23(1):365-376.

[10] Hiermaier S, Riedel W, Clegg R A, et al. Advanced material models for hypervelocity impact simulations[R]. EMI-Report, E43/99, ESA CR (P) 4305, 1999.

[11] Anderson C E, Cox P A, Johnson G R, et al. A constitutive formulation for anisotropic materials suitable for wave propagation computer program-Ⅱ[J]. Computational Mechanics, 1994,15(3):201-223.

[12] 錢偉長.穿甲力學(xué)[M].北京:國防工業(yè)出版社,1984.

(責(zé)任編輯 曾月蓉)

Numerical simulation on low-speed impact response of 2D plain-woven C/SiC composite

Yang Yang1,2, Xu Fei1, Zhang Yue-qing1, Tang Zhong-bin1

(1.SchoolofAeronautics,NorthwesternPolytechnicalUniversity,Xi’an,Shaanxi710072,China; 2.StateKeyLaboratoryofExplosionScienceandTechnology,BeijingInstituteofTechnology,Beijing100081,China)

First, the experiments that the steel balls impact to 2D-C/SiC composite under velocity of 79～ 219 m/s are investigated by using the air gun. Second, material parameters of 2D-C/SiC are obtained based on an orthotropic constitutive material model in Autodyn,and numerical simulation corresponding to experimental conditions are conducted based on smooth particle hydrodynamics(SPH) solver. The comparisons between the calculation results and experimental data of the debris cloud structure, the B scan results and the axis velocity of debris cloudvalidate the ability of this model for describing the brittle characteristics and the softening behaviour of 2D-C/SiC under impact load. Finally, the limit penetration depth of 2D-C/SiC under the impact of steel ball is predicted.

solid mechanics; impact resistance of C/SiC; smooth particle hydrodynamics; thermal protection structure; debris cloud; high-speed photography

10.11883/1001-1455(2015)01-0022-07

2014-01-16；

2014-05-19

國家自然科學(xué)基金重點(diǎn)項目(90916027); 北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室開放基金項目(KFJJ12-14M)

楊 揚(yáng)(1986— )，男，博士研究生,nwpu_yang@mail.nwpu.edu.cn。

O347.3 國標(biāo)學(xué)科代碼： 13015

A