梁 超

(昆明理工大學(xué) 國土資源工程學(xué)院，云南 昆明 650093)

?

采空區(qū)頂板動力響應(yīng)數(shù)值模擬研究

梁 超

(昆明理工大學(xué) 國土資源工程學(xué)院，云南 昆明 650093)

介于巖體材料及爆破地震波的復(fù)雜性，采用理論分析方法難以對采空區(qū)頂板動力響應(yīng)做出準(zhǔn)確計算，而且求解出的計算公式通常比較復(fù)雜，工程上難以直接應(yīng)用。FLAC3D數(shù)值模擬軟件中的動力分析功能在解決巖土工程方面的非線性動力分析問題方面有巨大優(yōu)勢，能更深入的研究采空區(qū)頂板在爆破震動作用下的動應(yīng)力場、動位移場和塑性區(qū)的變化規(guī)律及動靜應(yīng)力場的疊加作用機(jī)理。在理論研究的基礎(chǔ)上，借助FLAC3D提供的動力分析功能，建立數(shù)值計算模型，對靜力荷載和爆破動載共同作用下采空區(qū)頂板失穩(wěn)過程和機(jī)理進(jìn)行數(shù)值模擬分析。

爆破地震波；動力響應(yīng)；采空區(qū)頂板；數(shù)值模擬

0 前 言

爆破作業(yè)是金屬類礦山在開采過程中非常重要且必不可少的工序，不論是生產(chǎn)還是掘進(jìn)都需要進(jìn)行大量爆破作業(yè)，而爆破必然產(chǎn)生震動，對周圍巖體形成動力加載環(huán)境，如果動荷載強(qiáng)度達(dá)到某一程度會直接造成圍巖的動力失穩(wěn)破壞。另一方面，采空區(qū)頂板在爆破地震波的反復(fù)動力加載作用下，其宏觀力學(xué)性能將逐漸被弱化乃至最終失效或形成連續(xù)非線性的累積損傷。所以，爆破震動效應(yīng)對采空區(qū)圍巖造成的危害不能被忽略[1-3]。

介于巖體材料及爆破地震波的復(fù)雜性，采用理論分析方法難以對采空區(qū)頂板動力響應(yīng)做出準(zhǔn)確計算，而且求解出的計算公式通常比較復(fù)雜，工程上難以直接應(yīng)用。而隨著計算機(jī)軟件技術(shù)的廣泛應(yīng)用，為研究巖體結(jié)構(gòu)動力響應(yīng)提供了一種新方法。鑒于數(shù)值模擬的特點(diǎn)和求解非線性動力問題的優(yōu)越性，用數(shù)值計算的方法來研究爆破地震作用下采空區(qū)頂板結(jié)構(gòu)動力響應(yīng)的規(guī)律已成為爆破領(lǐng)域研究方法的一個必然趨勢[57]。FLAC3D數(shù)值模擬軟件中的動力分析功能在解決巖土工程方面的非線性動力分析問題具有巨大的優(yōu)勢。為更深入的研究采空區(qū)頂板在爆破震動作用下的動應(yīng)力場、動位移場和塑性區(qū)的變化規(guī)律及動靜應(yīng)力場的疊加作用機(jī)理，本章在前面理論研究的基礎(chǔ)上，借助FLAC3D提供的動力分析功能，建立數(shù)值計算模型，對靜力荷載和爆破動載共同作用下采空區(qū)頂板失穩(wěn)過程和機(jī)理進(jìn)行數(shù)值模擬分析。

1 數(shù)值模擬軟件概述

FLAC是連續(xù)介質(zhì)快速拉格朗日差分分析方法的英文縮寫，是專門為巖土工程穩(wěn)定性分析設(shè)計和服務(wù)的程序，屬于有限差分軟件。他能夠準(zhǔn)確模擬材料的塑性破壞和流動，能夠?qū)ν临|(zhì)巖石材料進(jìn)行力學(xué)特性模擬和塑性流動分析，具有與巖土工程實(shí)際擬合好，計算穩(wěn)定性好、收斂快的特點(diǎn)，因此目前是包括水電、礦業(yè)、公路、鐵路等巖土工程領(lǐng)域求解三維問題時普遍采用的理想分析工具[4]。

FLAC3D的基本計算原理是將計算區(qū)域劃分為若干四面體單元，用節(jié)點(diǎn)將各個單元相互連接，當(dāng)某個節(jié)點(diǎn)受到荷載作用后，其運(yùn)動方程可以用時間步長的有限差分形式表示，并且荷載只對該節(jié)點(diǎn)周圍的若干節(jié)點(diǎn)有影響。程序根據(jù)單元節(jié)點(diǎn)的速度變化來計算單元之間的相對位移，進(jìn)而可以求出單元應(yīng)變，然后根據(jù)單元材料的本構(gòu)關(guān)系即可求出單元應(yīng)力。

FLAC3D中為巖土體工程問題的求解提供了11種力學(xué)模型和5種計算模式，其中巖土工程常用的力學(xué)模型有Null model開挖模型、Mohr-Coulomb塑性模型和Drucker-Prager塑性模型；而常用的計算模式有靜力模式、動力模式、滲流模式及蠕變模式。其中，靜力模式是是FLAC3D默認(rèn)的計算模式。動態(tài)分析時，要先進(jìn)行靜力計算，并在此基礎(chǔ)上打開動力計算功能，設(shè)置動態(tài)邊界然后直接施加速度或應(yīng)力荷載進(jìn)行動力擾動分析。

2 數(shù)值模型的建立及計算方案

2.1計算模型區(qū)域及幾何尺寸

本文根據(jù)某礦豎向排列采空區(qū)建立數(shù)值計算模型。根據(jù)礦體開挖尺寸和上部圍巖特點(diǎn)，設(shè)置模型尺寸的長度為290 m，寬度為220 m，高度為200 m。模型由上下兩個橫截面為矩形的采空區(qū)組成，并且不考慮上方采空區(qū)跨度的變化，兩個采空區(qū)的長寬高尺寸均為42 m×36 m×18 m。以豎直方向?yàn)槟Ｐ偷腪軸，沿礦體走向方向?yàn)閅軸，水平方向內(nèi)垂直礦體走向設(shè)為X軸。模型的Z軸的正方向范圍即0 m≤Z≤200 m為數(shù)值計算區(qū)域，X方向的計算范圍是 ，Y方向的范圍為0 m≤Z≤290 m。網(wǎng)格單元的劃分采用不等劃分，根據(jù)需要對模型開挖部分的網(wǎng)格劃分的細(xì)一些，對模型的外部較遠(yuǎn)的區(qū)域劃分相對稀疏。模型內(nèi)共劃分了79 296個單元，85 140個節(jié)點(diǎn)，三維計算模型及網(wǎng)格劃分如圖1所示。

圖1 計算模型及網(wǎng)格分布

本次模擬分析中，靜力計算和動力計算均采用Mohr-Coulomb強(qiáng)度準(zhǔn)則，其屈服函數(shù)為：

(1)

ft=σ3-σt

(2)

式中：Nφ=(1+sinφ)/(1-sinφ)，其中φ為內(nèi)摩擦角；σ1，σ3分別為最大和最小主應(yīng)力；c為黏聚力；σt為巖石抗拉強(qiáng)度。當(dāng)巖體內(nèi)某一點(diǎn)應(yīng)力滿足fs<0時發(fā)生剪切破壞；當(dāng)巖體內(nèi)某一點(diǎn)應(yīng)力滿足ft>0時發(fā)生剪切破壞。計算所需的巖體力學(xué)參數(shù)見表1。

表1 巖體力學(xué)參數(shù)取值

2.2初始條件與邊界設(shè)置

2.2.1 初始地應(yīng)力場

原巖應(yīng)力是巖土工程數(shù)值模擬計算中的一個基本參數(shù)。其數(shù)值的大小，尤其是水平應(yīng)力與垂直應(yīng)力的比值，對計算結(jié)果的影響往往大于計算模型等因素的影響。因此，確定原巖應(yīng)力是計算分析中的一個基本問題。一般來說，原巖應(yīng)力大小和方向的確定，有效的方法是進(jìn)行現(xiàn)場實(shí)測。由于沒有進(jìn)行現(xiàn)場地應(yīng)力測試，初始地應(yīng)力場僅按自重應(yīng)力場考慮。根據(jù)彈性原理，豎向應(yīng)力σv=γH，水平應(yīng)力σh=kσv。其中：γ為巖體容重，H為埋深；k為側(cè)壓力系數(shù)，k=v/(1-v)；v為泊松比。

2.2.2 靜力計算邊界條件

靜力分析過程中采用位移邊界條件，由于采動的影響范圍有限，離采場較遠(yuǎn)處巖體位移值很小，將模型邊界處位移視為零，即在計算模型的左右(X方向)邊界、前后(Y方向)邊界和底邊界(Z=0)均施加位移約束條件，上邊界為自由邊界。模型頂部距離地面的平均距離為200 m，將模型上方巖體的作用轉(zhuǎn)化成面力垂直向下施加在模型的上部(Z=200)。

2.2.3 動力計算邊界條件

在進(jìn)行動力分析時，合理設(shè)置模型周圍的邊界條件是非常重要的，因?yàn)樵谶吔缟蠒a(chǎn)生波的反射影響，在一定程度上會對動力分析模擬的結(jié)果產(chǎn)生較大影響。如果把計算模型設(shè)置的很大，雖然可以消除反射波的影響，分析結(jié)果也會得到很大的改善，但是較大的計算模型會導(dǎo)致單元數(shù)與節(jié)點(diǎn)數(shù)增多，從而使計算時間大大增加。FLAC3D提供了靜態(tài)(黏性)和自由場兩種邊界條件設(shè)置來減小波在模型邊界上反射[61]。自由場邊界一般應(yīng)用于地面結(jié)構(gòu)的動力分析，考慮到采空區(qū)位于地表以下位置，因此本次計算模型動力分析在底邊界和側(cè)向邊界(包括X方向和Y方向)設(shè)置黏性邊界條件來吸收邊界上的入射波。上方采空區(qū)的底板設(shè)置為動荷載邊界，其他采空區(qū)邊界均設(shè)為自由邊界。

2.3力學(xué)阻尼與中心頻率的選取

采用FLAC3D求解動力方程的方法可解決準(zhǔn)靜力問題和動力問題兩類力學(xué)問題。在這兩類問題中都要使用阻尼，但準(zhǔn)靜力問題需要更多的阻尼使得動力方程能夠更快的收斂平衡。對于動力問題中的阻尼，需要在數(shù)值模擬中重現(xiàn)自然系統(tǒng)在動荷載作用下的阻尼大小。阻尼的存在使得能量產(chǎn)生了消耗，并且引起質(zhì)點(diǎn)位移和速度的衰減。在動力分析中，F(xiàn)LAC3D提供了瑞利阻尼、局部阻尼和人工粘滯阻尼3種形式[61]。常用臨界阻尼值范圍一般在2%～10%。對于巖土材料而言，在進(jìn)行動力計算時多數(shù)能量耗散在塑性流動階段，因此在進(jìn)行大應(yīng)變的動力分析時只需設(shè)置一個較小的阻尼比即可。本文采用的是瑞利阻尼。對于比較簡單的模型可用自振頻率作為瑞利阻尼的中心頻率，通過重力作用下無阻尼的自振計算求得。本文計算模型中心頻率取為30 Hz，阻尼比取5%時系統(tǒng)振幅方能得到較好的減弱。經(jīng)測試，地震波在0.3 s后完成基本衰減，故系統(tǒng)模型受動態(tài)擾動時間選取0.3 s。

2.4爆破震動載荷的施加

在FLAC3D動力計算中，動載荷輸入可以采用加速度時程、速度時程、位移時程和應(yīng)力時程4種不同方式，若采用粘滯邊界條件，則必須輸入速度時程進(jìn)行分析。數(shù)值計算中可以將施加的爆破地震荷載簡化為諧波的形式或者簡化為三角形荷載的形式。根據(jù)第四章的理論分析，本章采用諧波形式的應(yīng)力激勵動荷載模型，并將簡化后的諧波動荷載施加于上方采空區(qū)的底板，沿Z軸負(fù)方向傳播，進(jìn)行爆破動荷載作用下采空區(qū)頂板動力響應(yīng)狀況計算，包括應(yīng)力、位移及塑性區(qū)的變化與分布情況。

2.5計算方案與監(jiān)測點(diǎn)

為研究豎上下交疊采空區(qū)隔層頂板在開挖形成的靜力荷載以及爆破震動荷載作用下的穩(wěn)定性分析采空區(qū)頂板在不同振幅、不同頻率的爆破震動荷載作用下的動力響應(yīng)變化規(guī)律設(shè)計了以下4種計算方案：a 初始應(yīng)力條件下對礦體模型上下兩個空區(qū)同時進(jìn)行開挖，計算至靜力平衡，得到采空區(qū)隔板在自重靜荷載作用下的應(yīng)力、位移及塑性區(qū)并監(jiān)測隔層頂板下表面中間部位應(yīng)力和位移變化；b 在爆破地震波荷載峰值為1 MPa時，頻率分別為20，50，100 Hz的不同地震波頻率下隔層頂板的動力響應(yīng)；c 在頻率為50 Hz，地震波荷載峰值分別為0.5，1，1.5，2 MPa時的動荷載作用下頂板的動力響應(yīng)規(guī)律；d 隔層頂板在頻率為100 Hz，動荷載峰值為0.5 MPa的爆破動荷載作用下的響應(yīng)。

由固體力學(xué)理論可知，在自重和動荷載作用下采空區(qū)頂板的最大位移發(fā)生在空區(qū)間隔板的底部中心位置。所以為研究采空區(qū)隔板在爆破荷載作用下的動力響應(yīng)變化規(guī)律，選取隔板中間豎直方向截面上的節(jié)點(diǎn)作為位移和應(yīng)力的跟蹤監(jiān)測點(diǎn)，在隔板厚度為10 m時該點(diǎn)在模型上的坐標(biāo)為(145,1,73)。

整個計算過程分二步進(jìn)行。先進(jìn)行靜力分析，然后在此基礎(chǔ)上再施加動力荷載進(jìn)行動力問題的分析。在第1步靜力分析前，先關(guān)閉動態(tài)分析模式，在模型邊界施加水平和豎直方向應(yīng)力，形成初始地應(yīng)力場，然后將采空區(qū)一次性開挖并計算至平衡狀態(tài)，分析采空區(qū)頂板在靜荷載作用下的變形，并將其圍巖靜態(tài)應(yīng)力場結(jié)果作為第2步動態(tài)擾動的初始應(yīng)力場。第2步輸入動力荷載簡化為諧波形式的應(yīng)力激勵施加于模型上部空區(qū)底板沿Z軸負(fù)方向傳播，分析其在爆破荷載作用下的動態(tài)響應(yīng)特征。

3 模擬結(jié)果分析與討論

FLAC在分析輸出的應(yīng)力結(jié)果時，他給出的最大主應(yīng)力，實(shí)際上是最小主應(yīng)力，最小主應(yīng)力實(shí)際上是最大主應(yīng)力，用這兩個應(yīng)力指標(biāo)，可表征地壓活動在介質(zhì)應(yīng)力狀態(tài)變化方面產(chǎn)生應(yīng)力集中或應(yīng)力松馳的程度。且負(fù)應(yīng)力為壓應(yīng)力，正應(yīng)力為拉應(yīng)力。根據(jù)以上設(shè)置建立的模型首先進(jìn)行原巖應(yīng)力計算，然后對礦體開挖進(jìn)行靜力計算，開挖后頂板在靜荷載作用下的最大主應(yīng)力云圖、最小主應(yīng)力云圖、豎直方向位移云圖以及拉伸剪切破壞區(qū)域如圖2所示。監(jiān)測點(diǎn)的應(yīng)力和位移隨開挖過程變化的關(guān)系見圖3。

從圖2和圖3可以看出：采空區(qū)開挖后圍巖原始應(yīng)力平衡被打破，引起應(yīng)力重分布和局部應(yīng)力集中，并產(chǎn)生了向采空區(qū)方向的位移。在上部采空區(qū)的頂板處和下部采空區(qū)的底板處均出現(xiàn)較大的拉應(yīng)力，在空區(qū)邊界和拐角處出現(xiàn)拉伸剪切破壞區(qū)，隔板比較完好中間未出現(xiàn)拉伸破壞區(qū)域。當(dāng)頂板巖體達(dá)到再次平衡時，在底部中心處產(chǎn)生的位移和拉應(yīng)力分別為2 mm、0.85 MPa。

圖2 開挖后隔層頂板的位移和應(yīng)力云圖及塑性區(qū)

圖3 監(jiān)測點(diǎn)的應(yīng)力和位移隨開挖過程變化曲線

圖4分別為在荷載峰值為1 MPa，頻率分別為20，50，100 Hz的爆破動荷載作用下隔層頂板的最大主應(yīng)力和豎向位移局部放大云圖。

圖4 頂板在不同頻率動荷載作用下的最大主應(yīng)力和位移云圖

圖5分別為在荷載峰值為1 MPa，頻率分別為20，50，100 Hz的爆破動荷載作用下隔層頂板底部監(jiān)測點(diǎn)的位移和應(yīng)力隨動力計算時間的變化曲線。

圖5 頂板在不同頻率動荷載作用下的動力響應(yīng)曲線

從圖4和圖5可以看出：在爆破荷載峰值相同，頻率不同的爆破地震波作用下，隔層頂板的動力響應(yīng)也不同。頂板的動力響應(yīng)隨著頻率的增加先增大后減小，爆破地震荷載頻率越大，頂板產(chǎn)生的動拉應(yīng)力和動位移越小。當(dāng)頻率為50 Hz時頂板的動力反應(yīng)最大，說明50 Hz頻率的爆破地震波對采空區(qū)頂板的穩(wěn)定在影響最大，隔層頂板結(jié)構(gòu)的固有頻率與爆破地震波的頻率接近，產(chǎn)生共振效應(yīng)，可能使頂板發(fā)生突然失穩(wěn)破壞。

圖6分別為在荷載峰值為0.5 MPa，頻率50 Hz的動荷載作用下隔層頂板最大主應(yīng)力和豎向位移局部放大云圖以及隔層頂板底部監(jiān)測點(diǎn)的位移和應(yīng)力隨動力計算時間的變化曲線。

圖6 頻率為50 Hz，荷載峰值為0.5 MPa時頂板的動力響應(yīng)規(guī)律

圖7分別為在荷載峰值為0.5 MPa，頻率100 Hz的動荷載作用下隔層頂板最大主應(yīng)力和豎向位移局部放大云圖以及隔層頂板底部監(jiān)測點(diǎn)位移和應(yīng)力隨動力計算時間的變化曲線。

圖7 頻率為100 Hz，荷載峰值為0.5 MPa時頂板動力響應(yīng)規(guī)律

圖8是在頻率為50 Hz，地震荷載峰值分別為0.5，1，1.5，2 MPa時的隔層頂板塑性區(qū)變化。

從圖8可以看出：在不同的爆破地震荷載強(qiáng)度作用下頂板的動力破壞程度不同，爆破荷載峰值越大，頂板的動力響應(yīng)就越大，產(chǎn)生的拉應(yīng)力也越大。當(dāng)動荷載峰值為0.5 MPa時，只有隔層頂板上下空區(qū)的邊界出現(xiàn)塑性區(qū)，而當(dāng)荷載峰值大于1 MPa時隔層頂板受到的拉應(yīng)力大于其抗拉強(qiáng)度，底部中間區(qū)域首先發(fā)生拉伸破壞，開始出現(xiàn)小范圍塑性區(qū)，并隨著爆破荷載峰值的增大，頂板發(fā)生拉伸破壞的范圍也越大，并且破壞區(qū)域的形狀為拱形。

圖8 不同動荷載峰值作用下頂板塑性區(qū)變化規(guī)律

為了更清楚的觀察爆破地震波的頻率和振動幅值對頂板動力響應(yīng)的規(guī)律，提取監(jiān)測點(diǎn)在不同頻率不同荷載幅值下的最大拉應(yīng)力變化值，并將其導(dǎo)入到Matlab中繪制出頂板的動力響應(yīng)與爆破荷載及頻率的關(guān)系曲線如圖9所示。

圖9頂板的動力響應(yīng)與爆破荷載及頻率的關(guān)系

從圖9可以清晰的看出：采空區(qū)隔層頂板在不同頻率、不同荷載峰值的地震波作用下的應(yīng)力變化過程。從以上數(shù)值模擬分析可以得出以下結(jié)論：爆破地震荷載峰值大小對采空區(qū)隔板穩(wěn)定性的影響最明顯，隔板最大位移和最大拉應(yīng)力隨峰值增大而增大，當(dāng)動荷載峰值超過1 MPa時隔板底部中心部位開始出現(xiàn)拉伸破壞區(qū)，并且動荷載峰值越大破壞區(qū)域越大；隨頻率的增大，隔板最大位移和最大應(yīng)力先增大后減小，當(dāng)爆破荷載頻率接近頂板的固有頻率時，頂板產(chǎn)生的動力反應(yīng)達(dá)到最大值，即50 Hz頻率的動荷載對采空區(qū)隔板的穩(wěn)定性最不利。

4 結(jié) 論

本文以上下交疊地下礦山采空區(qū)隔層頂板為研究對象，采用數(shù)值模擬方法對采空區(qū)頂板在爆破震動作用下的動力響應(yīng)進(jìn)行了數(shù)值模擬研究。應(yīng)用FLAC3D計算軟件，建立豎向排列的采空區(qū)模型。首先模擬出空區(qū)開挖以后在靜力載荷下頂板的位移和應(yīng)力變化情況，并在靜荷載作用的基礎(chǔ)上對模型施加不同峰值的爆破動荷載，對頂板在動靜荷載共同作用下的位移場、應(yīng)力場和塑性區(qū)的響應(yīng)進(jìn)行了研究，通過對比不同爆破荷載峰值情況下采空區(qū)頂板的動力響應(yīng)來分析對穩(wěn)定性的影響程度，揭示了動靜荷載共同作用下豎向排列采空區(qū)隔層頂板失穩(wěn)破壞機(jī)理。

[1] 劉亞群, 李海波. 爆破荷載作用下黃麥嶺磷礦巖質(zhì)邊坡動態(tài)響應(yīng)的UDEC模擬研究[J]. 巖石力學(xué)與工程學(xué)報, 2004, 23(11): 42-45.

[2] 陳占軍, 朱傳云, 周小恒. 爆破荷載作用下巖石邊坡動態(tài)響應(yīng)的FLAC3D模擬研究[J]. 爆破, 2005, 22(4): 8-15.

[3] 逄煥東, 陳士海. 彈性介質(zhì)中爆破地震波傳播的分區(qū)變化規(guī)律研究[J]. 振動與沖擊, 2009(3).

[4] 閆長斌, 徐國元, 李夕兵. 爆破震動對采空區(qū)穩(wěn)定性影響的FLAC3D分析[J]. 巖石力學(xué)與工程學(xué)報, 2005, 8(16): 2894-2899.

[5] Khandelwal M, Singh T N. Prediction of blast induced ground vibrations and frequency in opencast mine: A neural network approach[J]. Journal of sound and vibration, 2006, 289(4): 711-725.

[6] 傅洪賢, 趙勇, 謝晉水, 等. 隧道爆破近區(qū)爆破振動測試研究[J]. 巖石力學(xué)與工程學(xué)報, 2011, 30(2): 335-340.

[7] 吳姍, 陳何, 孫宏生, 等. 大紅山銅礦束狀孔采礦及爆破振動監(jiān)測技術(shù)[J]. 有色金屬工程, 2015, 5(5): 126-129.

[8] 劉波, 韓彥輝. FLAC原理,實(shí)例與應(yīng)用指南[M]. 北京: 人民交通出版社, 2005.

AStudyonNumericalSimulationofDynamicResponseofRoof

LIANG Chao

(FacultyofLandResourceEngineering,KunmingUniversityofScienceandTechnology,Kunming,Yunnan650093,China)

The roof of mined area is relatively weak, and the goaf roof is one of the most serious accidents in underground metal mines. The underground workers will not only pose a threat to the safety, but also have a significant impact on the production and management of underground pressure. There are many factors to affect the goaf roof stability. The stress distribution and the blasting vibration are the two most important factors. This paper mainly uses the method of theoretical analysis to be combined with numerical simulation, from static and dynamic two aspects of goaf roof stability of the stress in the weight and load of blasting under dynamic loads. The main contents of the thesis are the following: using the elastic vibration theory of thin plate down the top plate under uniform harmonic pulse displacement under load to calculate formulas for stress analysis. The damping, including the dynamic load strength, influence of dynamic frequency response of roof strata, and the dynamic response in the empty mining district under dynamic and static loads of roof are all analyzed.

Blasting earthquake wave; Roof of mined void; Dynamic response; Numerical simulation

2016-07-24

梁超(1991-)，男，寧夏石嘴山市人，在讀碩士研究生，研究方向：安全工程，手機(jī)：18695200970，E-mail：376677756@qq.com.

TD327.2

：Adoi：10.14101/j.cnki.issn.1002-4336.2016.04.019