多軟件聯(lián)合的煤礦三維數(shù)值建模方法及其應(yīng)用*

朱權(quán)潔，劉金海，王勝開，高林生，梁　娟

(1. 華北科技學(xué)院 安全工程學(xué)院，北京 東燕郊 101601；2. 馳宏鋅鍺股份有限公司，云南 曲靖 655011；3. 防災(zāi)科技學(xué)院，北京 東燕郊 101601)

0　引言

礦井地應(yīng)力場測量與反演是確定工程巖體力學(xué)屬性、進(jìn)行圍巖穩(wěn)定性分析以及地下工程開挖設(shè)計科學(xué)化的必要前提[1-2]，也是有效防治沖擊地壓等災(zāi)害、保證礦井安全有序生產(chǎn)的關(guān)鍵。目前，地應(yīng)力測試都是在硬巖層中進(jìn)行的，煤層應(yīng)力值除可根據(jù)已知點對應(yīng)巖層地應(yīng)力實測結(jié)果進(jìn)行估算外，目前尚無其他較好辦法，但這種方法也具有局限性和不穩(wěn)定性。要全盤掌握礦區(qū)內(nèi)不同深度、不同位置、同一煤層不同深度上的應(yīng)力分布，需進(jìn)行區(qū)域地應(yīng)力場的反演分析，這其中地質(zhì)模型的數(shù)值計算是關(guān)鍵步驟。

數(shù)值計算是當(dāng)前解決科研難題、實際工程問題的有效手段[3-5]。在巖土工程領(lǐng)域，有限元法、離散元法以及有限差分法等是數(shù)值分析計算的重要方法之一。Itasca公司研發(fā)的FLAC3D數(shù)值模擬軟件是其中的優(yōu)秀代表之一，該軟件對于模擬不產(chǎn)生離散、連續(xù)變形模型具有良好效果，因此，在巖土工程、采礦工程等行業(yè)和領(lǐng)域應(yīng)用廣泛。

在具有強(qiáng)大數(shù)值分析能力的同時，F(xiàn)LAC3D的建模尤其是復(fù)雜模型的構(gòu)建能力是其最大短板。對于規(guī)則模型而言，利用fish語言編寫命令流可實現(xiàn)快速建模，但不規(guī)則、復(fù)雜的地質(zhì)模型依靠命令流實現(xiàn)是非常困難的，這也是FLAC3D面臨的一大挑戰(zhàn)[6]。現(xiàn)有研究成果中多選擇“先外部構(gòu)建幾何模型、后導(dǎo)入FLAC3D”的方式開展[7-9]。如林杭等編制了Surpac到FLAC3D的接口程序，直接實現(xiàn)Surpac三維模型到FLAC3D的轉(zhuǎn)換；紀(jì)洪廣等[[10]提出一種利用VC++開發(fā)的復(fù)雜礦山地質(zhì)形體網(wǎng)格精細(xì)劃分方法，解決了FLAC3D構(gòu)建復(fù)雜地質(zhì)體費時、繁瑣等缺陷；趙文等[11]提出利用MATLAB和Surfer劃分網(wǎng)格，并生成FLAC3D命令流語句，提高了FLAC3D建模速度。這些方法中前者在網(wǎng)格劃分方面存在缺陷，后者對使用者的編程能力要求較高。

此外，袁海平等[12]提出聯(lián)合AutoCAD，ANSYS和FLAC3D聯(lián)合構(gòu)建三維地質(zhì)模型，該方法在CAD中制作面元素，仍然受限于構(gòu)建大型復(fù)雜三維模型；李翔等[13]基于鉆孔柱狀圖資料，提出利用SURPAC，CAD和ARCGIS建立復(fù)雜三維地址模型。該方法一定程度上解決上了上述問題，但在利用ANSYS的APDL語言生成模型時涉及大量命令流操作，影響模型構(gòu)建的速度，如果遇到命令錯誤會導(dǎo)致模型無法完整生成。由此可見，如何有效、快速構(gòu)建復(fù)雜三維地質(zhì)模型進(jìn)而開展數(shù)值計算是地應(yīng)力反演的前提和重要研究課題[14-15]。

鑒于上述原因，本文提出了一種復(fù)雜三維地質(zhì)模型的構(gòu)建方法，該方法聯(lián)合3DMINE，GOCAD，SURFER，RHIRO，ANSYS以及FLAC3D以及MATLAB共7種軟件共同建模，減少了編程工作量、兼顧了ANSYS強(qiáng)大的前處理功能，彌補(bǔ)了大量命令流編寫、無法合理劃分網(wǎng)格等缺陷，實現(xiàn)了復(fù)雜三維地質(zhì)體的快速構(gòu)建。

1　建模思路

以山東鄆城煤礦為例，對整個建模思路進(jìn)行介紹。為了客觀反演地應(yīng)力的分布規(guī)律，綜合考慮地形、地質(zhì)構(gòu)造以及地層介質(zhì)等因素，結(jié)合前期開展的地質(zhì)勘探工作，最后確立了目標(biāo)研究區(qū)域的范圍：以坐標(biāo)點A(7 344.8，3 958.1)，B(9 324.8，3 958.1)，C(9 324.8，9 998.1)和D(7 344.8，9 998.1)圈定的長方形區(qū)域為計算區(qū)域，如圖1中所示。

圖1　幾何模型構(gòu)建區(qū)域Fig.1　Construction area of geometric model

為了建立該的三維地質(zhì)模型，使該模型更貼進(jìn)實際的地形地貌，減少因數(shù)值模擬不精確帶來的地應(yīng)力反演誤差，開展了礦區(qū)復(fù)雜三維地質(zhì)模型精細(xì)化構(gòu)建研究_建模原理基于“點-線-面-體”逐步構(gòu)建的思路，在構(gòu)建過程中使用了GOCAD，ANSYS，F(xiàn)LAC3D等多種軟件，具體過程及思路如圖2所示。

圖2　三維地質(zhì)模型構(gòu)建過程及思路Fig.2　Process of 3D numerical modeling

2　點陣數(shù)據(jù)和面的形成

常規(guī)數(shù)值模擬一般采用規(guī)則建模方法構(gòu)建模型，這類建模方式由于軟件自身和建模方式不靈活從而無法完整表現(xiàn)地質(zhì)體的起伏、斷層構(gòu)造等信息。為了精細(xì)地描述復(fù)雜的地層結(jié)構(gòu)，提出利用鉆孔信息生成巖層界面等高線，并“以點構(gòu)線、以線構(gòu)面、以面構(gòu)體”的三級建模法構(gòu)建模型，以形成更真實的復(fù)雜的地質(zhì)體模型，提高數(shù)值模擬的精確度和可靠性。

2.1　等高線的生成

勘探鉆孔是構(gòu)建巖層界面等值線的數(shù)據(jù)來源。鉆孔含有大量信息，通過鉆孔柱狀圖可以了解礦體、巖層界線以及斷層構(gòu)造等隱藏于地下的信息。根據(jù)這些信息，可以生成相應(yīng)的巖層界面等高線。利用CAD建立三維等高線操作較為麻煩，因此，選擇借助SURPAC或3DMINE等數(shù)字礦山三維軟件直接生成等高線，其操作步驟如下。

1)準(zhǔn)備鉆孔原始數(shù)據(jù)。從柱狀圖中讀取相應(yīng)的巖性信息表、鉆孔空間坐標(biāo)，并導(dǎo)入到三維軟件之中。

2)三維鉆孔形態(tài)的生成。根據(jù)數(shù)據(jù)庫中存儲的鉆孔信息，生成三維形態(tài)的鉆孔立體圖，如圖3(a)所示。圖中以開孔位置為起點、終孔位置為終點連接成線，并按照鉆孔巖性表中的巖性特點分段標(biāo)記。

3)三維等高線的生成。根據(jù)這些不同地層的不同標(biāo)記，提取相應(yīng)地層的點坐標(biāo)，這些點坐標(biāo)經(jīng)過分步取樣、插值即可生成等高線圖，如圖3(b)所示的煤層頂板等高線。

圖3　利用鉆孔信息生成等值線Fig.3　Contour map generation using borehole data

鉆孔數(shù)量的多少影響數(shù)值建模的精細(xì)程度，因此，可根據(jù)工程實際要求設(shè)計不同的鉆孔數(shù)量需求。如存在大量外部數(shù)據(jù)(如等高線)，可直接進(jìn)行加載，導(dǎo)入GOCAD進(jìn)行離散處理。

2.2　離散點陣的生成

由于前期鉆孔信息有限，所生成的等高線上僅存有少量離散的擬合點坐標(biāo)數(shù)據(jù)，這些數(shù)據(jù)如果直接導(dǎo)入ANSYS或FLAC3D進(jìn)行建模，一方面將造成模型偏差較大(嚴(yán)重失真)，另一方面利用少量點構(gòu)建體操作性不強(qiáng)。因此，需要對等高線進(jìn)行離散化和細(xì)分處理，生成相應(yīng)的等高線擬合點數(shù)據(jù)。

經(jīng)過3DMINE處理后生成dxf文件，并導(dǎo)入到GOCAD中。GOCAD讀取dxf數(shù)據(jù)，提取出等高線中擬合點的數(shù)據(jù)。如圖4所示，為煤層頂板等高線的網(wǎng)格點陣圖，其中包括等高線和陣列點的空間坐標(biāo)信息。

圖4　利用GOCAD生成網(wǎng)格點陣坐標(biāo)(煤層頂板)Fig.4　Grid lattice coordinates using GOCAD

2.3　數(shù)據(jù)網(wǎng)格化處理

前期得到的GOCAD導(dǎo)出的等高線坐標(biāo)點數(shù)據(jù)有限且較為離散，直接導(dǎo)入ANSYS構(gòu)建實體模型較為復(fù)雜，因此，借助SURFER軟件對離散數(shù)據(jù)進(jìn)行插值處理，使數(shù)據(jù)點均勻化，即網(wǎng)格化操作。

利用SURFER軟件“數(shù)據(jù)”選項讀取上一步獲得的離散等高線擬合點數(shù)據(jù)，選擇“克里金插值算法”對離散數(shù)據(jù)進(jìn)行插值和均勻化處理，即得到等高線的點陣數(shù)據(jù)集，生成線框圖如圖5所示。

圖5　點陣數(shù)據(jù)構(gòu)建的線框圖Fig.5　Wireframe generation using grid lattice coordinates

需要注意的是，由于地質(zhì)模型在導(dǎo)入ANSYS中時會考慮模型大小問題，因此，需要對鉆孔的絕對坐標(biāo)進(jìn)行偏移，以保證模型的順利生成。坐標(biāo)系偏移為：X偏移-20479 680，Y偏移-3911 100(這是后文進(jìn)行反演時的重要參數(shù)，以此對比CAD圖紙上任一點與幾何模型坐標(biāo)的關(guān)聯(lián)性)。

2.4　不規(guī)則面的生成

SURFER提供的點坐標(biāo)數(shù)據(jù)無法直接導(dǎo)入ANSYS模型，借助命令流過程復(fù)雜。因此，為減少編程帶來的工作量，將SURFER生成的dat文件導(dǎo)入RHINO三維軟件，并利用RHINO提供的srfptgrid命令將數(shù)據(jù)點連線成面，如圖6所示，最后將該面文件轉(zhuǎn)換為ANSYS可讀的iegs格式文件(該格式可供ANSYS直接讀取)。如圖6所示，利用RHINO軟件生成的三維面模型，分別展示了該面模型的俯視圖、前視圖、右視圖和透視圖。

圖6　模型的生成與格式轉(zhuǎn)換Fig.6　Face model generation and data format conversion

3　模型的建立與數(shù)值計算

FLAC3D是采礦、巖土工程等領(lǐng)域的專業(yè)數(shù)值模擬軟件，在數(shù)值計算和分析方面有先天優(yōu)勢，但其建模功能受制于命令流式的手動建模方式，對于構(gòu)建復(fù)雜模型過程復(fù)雜。ANSYS有限元程序具有強(qiáng)大的前處理功能，因此，借助ANSYS軟件建立三維實體模型、劃分網(wǎng)格單元，然后導(dǎo)入FLAC3D進(jìn)行數(shù)值計算，將大大簡化建模過程，提高數(shù)值模擬的效率和精度。

3.1　幾何模型構(gòu)建

ANSYS提供了2種建模方式，一是創(chuàng)建或?qū)雽嶓w模型，然后進(jìn)行網(wǎng)格劃分，并生成有限元模型；另一種是利用單元和節(jié)點編號，依靠ANSYS提供的命令流進(jìn)行生成操作。本文選擇前者，其建模過程可簡述如下。

1)將前面處理好的關(guān)鍵巖層分界面iges文件導(dǎo)入WORKBENCH的Geometry之中，形成如圖7(a)線框圖，(b)面域圖所示的面文件。

2)將相鄰兩面選取，點擊Skin命令蒙皮成體，“面面構(gòu)體”生成體單元，如圖7(c)所示，完成后組成新的Part部分。

3)轉(zhuǎn)入到“Model”模式，選擇合理參數(shù)，對模型進(jìn)行網(wǎng)格劃分。

4)將Model中劃分網(wǎng)格后的模型導(dǎo)入到Mechanical APDL模式。

5)然后利用河海大學(xué)鄭文棠博士編寫的ANSYS_TO_FLAC3D接口程序?qū)NSYS模型信息導(dǎo)出，另存為節(jié)點(Node)和單元(Element)2個文件，并經(jīng)過轉(zhuǎn)換文件轉(zhuǎn)換為*.flac3d文件。FLAC3D通過Import Grid即可直接讀取該文件，得到如圖8所示的三維地質(zhì)模型。

圖7　ANSYS中由面構(gòu)體形成的三維地質(zhì)模型Fig.7　3D model generation using the face data

3.2　網(wǎng)格劃分

礦山工程巖體應(yīng)力具有三維特征，然而三維數(shù)值計算及成圖復(fù)雜，對計算軟件要求高，再加上反演區(qū)域范圍大，劃分單元多，給計算帶來難度。為此，模型大小及網(wǎng)格劃分尤為關(guān)鍵。

ANSYS提供了強(qiáng)大的網(wǎng)格劃分功能，相比FLAC3D更為智能和便利。綜合考慮有限元數(shù)值模擬的速度和構(gòu)建效果，網(wǎng)格劃分時選定的網(wǎng)格大小為80 m，最小網(wǎng)格為0.5 m，采用六面體網(wǎng)格劃分法將計算模型進(jìn)行劃分，共劃分得到114 240個單元，488 179個節(jié)點，網(wǎng)格劃分完畢后的效果如圖8所示。

圖8　網(wǎng)格劃分結(jié)果Fig.8　Mesh partition

3.3　FLAC3D中模型建立與數(shù)值計算

在ANSYS WORKBENCH中建立好三維地質(zhì)模型后，對實體模型進(jìn)行單元劃分，將模型導(dǎo)入到ANSYS APDL中，然后利用河海大學(xué)鄭文棠博士編寫的ANSYS_TO_FLAC3D接口程序?qū)NSYS模型信息導(dǎo)出，另存為節(jié)點(Node)和單元(Element)兩個文件，并經(jīng)過轉(zhuǎn)換文件轉(zhuǎn)換為*.flac3d文件。FLAC3D通過Import Grid即可直接讀取該文件，得到如圖9所示的三維地質(zhì)模型。

圖9　有限元數(shù)值模型的建立與導(dǎo)入Fig.9　Establishment and loading of finite element model

根據(jù)實際情況對不同巖層分別賦予材料屬性，并進(jìn)行分組處理；然后依據(jù)前期分析確立相應(yīng)的邊界條件和計算參數(shù)，最后導(dǎo)入模型進(jìn)行數(shù)值計算。

4　地應(yīng)力反演結(jié)果展示與分析

以山東鄆城煤礦為例，根據(jù)該礦前期測試的4個測點的地應(yīng)力數(shù)據(jù)，通過建立三維數(shù)值計算模型，利用先進(jìn)的數(shù)值模擬軟件(FLAC3D)和應(yīng)力反分析技術(shù)，反演和重構(gòu)鄆城煤礦一采區(qū)初始地應(yīng)力場。最終獲得的鄆城煤礦初始地應(yīng)力場回歸模型，回歸模型的相關(guān)性系數(shù)為88.26%，表明擬合效果較好。

圖10所示為利用回歸系數(shù)重構(gòu)的初始地應(yīng)力場SZZ的分布情況，取其中3組巖層上覆巖層、煤層和底部巖層展示。

圖10　煤巖層SZZ應(yīng)力分布規(guī)律Fig.10　SZZ stress distribution of coal seam and rock strata

圖11　-800 m最大主應(yīng)力SIG1分布云圖Fig.11　Nephogram of maximum principal stress on -800 m level

為了更清晰展示鄆城煤礦一采區(qū)區(qū)域地應(yīng)力概況，將反演重構(gòu)的煤層地應(yīng)力場投影到CAD圖上，得到如圖11所示的最大主應(yīng)力SIG1分布云圖。可以看出，最大主應(yīng)力與埋深和構(gòu)造密切相關(guān)。圖中深色區(qū)域為主應(yīng)力較大區(qū)域，這些區(qū)域正好位于褶曲構(gòu)造位置，由此可以推斷，褶曲是引起區(qū)域內(nèi)巖體應(yīng)力高度集中的主因；同時，區(qū)域內(nèi)的水平應(yīng)力大于垂直應(yīng)力。

利用上述建模方法，實現(xiàn)了對復(fù)雜礦山地質(zhì)模型的快速構(gòu)建，這為準(zhǔn)確計算和展示地應(yīng)力場分布規(guī)律提供了依據(jù)，為區(qū)域地應(yīng)力場反演奠定了基礎(chǔ)，同時也為礦井深部開采和采掘規(guī)劃提供了依據(jù)。

5　結(jié)論

1)利用3DMINE，SUFER，ANSYS以及FLAC3D等軟件聯(lián)合建立大型精細(xì)三維有限元模型的方法，通過“點-線-面-體”的建模思路，實現(xiàn)了復(fù)雜三維地質(zhì)模型的構(gòu)建。

2)以山東鄆城煤礦為例，通過構(gòu)建了該礦的三維地質(zhì)模型，結(jié)合現(xiàn)場實測測點地應(yīng)力值，反演出該礦的地應(yīng)力場分布，并以應(yīng)力云圖的方式進(jìn)行直觀展示，效果良好。

3)實際工程應(yīng)用表明，該方法可實現(xiàn)大尺度、復(fù)雜三維地質(zhì)模型的構(gòu)建，有效解決了復(fù)雜模型構(gòu)建效率低、網(wǎng)格劃分難度大等問題，通過數(shù)據(jù)的逐步對接，減少了中間過程的人為干擾和誤差，一定程度縮短了建模時間、提高了數(shù)值計算精度。

[1]康紅普，林健，張曉. 深部礦井地應(yīng)力測量方法研究與應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報，2007，26(5)：929-933.

KANG Hongpu, LIN Jian, ZHANG Xiao. Research and application of in-situ stress measurement in deep mines[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(5): 929-933.

[2]楊德林. 巖土工程問題的反演理論與工程實踐[M]. 北京：科學(xué)出版社，1996.

[3]龔曉南．對巖土工程數(shù)值分析的幾點思考[J]．巖土力學(xué)，2011，32(2)：321-325.

GONG Xiaonan. Reflections on numerical analysis of geotechnical engineering[J]. Rock and Soil Mechanics, 2011, 32(2): 321-325.

[4]朱志潔，張宏偉，韓軍，等. 不同應(yīng)力條件下巷道穩(wěn)定性研究[J]. 中國安全生產(chǎn)科學(xué)技術(shù)，2015，11(11)：11-16.

ZHU Zhijie, ZHANG Hongwei, HAN Jun, et al. Research on stability of roadway under different conditions of stress field[J]. Journal of Safety Science and Technology, 2015, 11(11): 11-16.

[5]劉寧，朱維申，辛小麗. 雙江口水電站初始地應(yīng)力場反演回歸分析[J]. 山東大學(xué)學(xué)報(工學(xué)版)，2008，38(6)：121-126.

LIU Ning, ZHU Weishen, XIN Xiaoli. Back regression analysis on the initial geostress field of the Shuangjiangkou hydropower station[J]. Journal of Shandong university: Engineering Science, 2008, 38(6): 121-126.

[6]鄭文棠，徐衛(wèi)亞，童富果，等．復(fù)雜邊坡三維地質(zhì)可視化和數(shù)值模型構(gòu)建[J]．巖石力學(xué)與工程學(xué)報，2007，26(6)：1633-1644．

ZHENG Wentang, XU Weiya, TONG Fuguo, et al. 3D geological visualization and numerical modeling of complicated slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(6): 1633-1644.

[7]裴啟濤，李海波，劉亞群，等. 復(fù)雜地質(zhì)條件下壩區(qū)初始地應(yīng)力場二次反演分析[J]. 巖石力學(xué)與工程學(xué)報，2014，33(S1)：2779-2785.

PEI Qitao, LI Haibo, LIU Yaqun, et al. Two-stage back analysis of initial geostress field of dam areas under complex geological conditions[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(S1): 2779-2785.

[8]崔芳鵬，胡瑞林，劉照連，等．基于Surfer平臺的FLAC3D復(fù)雜三維地質(zhì)建模研究[J]．工程地質(zhì)學(xué)報，2008(5)：699-702.

CUI Fangpeng, HU Ruilin, LIU Zhaolian, et al. Surfer software platform based on complex three dimensional geological digital models for pre-processing of FLAC3D[J]. Journal of Engineering Geology, 2008(5): 699-702.

[9]廖秋林，曾錢幫，劉彤，等．基于ANSYS平臺復(fù)雜地質(zhì)體FLAC3D模型的自動生成[J]．巖石力學(xué)與工程學(xué)報，2005，24(6)：1010-1013.

LIAO Qiulin, ZENG Qianbang, LIU Tong, et al. Automatic model generation of complex geologic body with FLAC3D based on ansys platform[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(6): 1010-1013.

[10]向鵬，紀(jì)洪廣，鄒靜，等. 復(fù)雜礦山地質(zhì)形體的三維數(shù)值網(wǎng)格劃分自動化平臺[J]. 巖土力學(xué)，2015，36(4)：1211-1216.

XIANG Peng, JI Hongguang, ZOU Jing, et al. An automation platform for generating the 3D numerical grids of mining fields with complex geomorphology[J]. Rock and Soil Mechanics, 2015, 36(4): 1211-1216.

[11]趙文，段魁，王彥東，等. 基于MATLAB、SURFER平臺的FLAC3D復(fù)雜三維地質(zhì)建模研究[J]. 工程地質(zhì)學(xué)報，2016，24(s1)：208-213.

ZHAO Wen, DUAN Kui, WANG Yandong, et al. Study on complex three-dimensional geological modeling of Flac3D based on Matlab and Surfer platform[J]. Journal of Engineering Geology, 2016, 24(s1): 208-213.

[12]袁海平，王金安，蔡美峰. 復(fù)雜地應(yīng)力場幾何跨尺度反演與重構(gòu)方法研究[J]. 采礦與安全工程學(xué)報，2011，28(4)：589-595.

YUAN Haiping, WANG Jinan, CAI Meifeng. Geometric cross-scale back-analysis and reconstruction menthod for complicated geo-stress field[J]. Journal of Mining & Safety Engineering, 2011, 28(4): 589-595.

[13]李翔，王金安，張少杰. 復(fù)雜地質(zhì)體三維數(shù)值建模方法研究[J]. 西安科技大學(xué)學(xué)報，2012，32(6)：676-681.

LI Xiang, WANG Jinan, ZHANG Shaojie. 3D modeling method of complicated geological body[J]. Journal of Xi’an University of Science and Technology, 2012, 32(6): 676-681.

[14]王金安，李飛. 復(fù)雜地應(yīng)力場反演優(yōu)化算法及研究新進(jìn)展[J]. 中國礦業(yè)大學(xué)學(xué)報，2015，44(2)：189-205.

WANG Jinan, LI Fei. Review of inverse optimal algorithm of in-situ stress and new achievement[J]. Journal of China University of Mining & Technology, 2015, 44(2): 189-205.

[15]趙自強(qiáng)，王學(xué)潮，隨裕紅. 南水北調(diào)西線工程深埋隧洞巖爆與地應(yīng)力研究[J]. 華北水利水電學(xué)院學(xué)報，2002，23(1)：48-50.

ZHAO Ziqiang, WANG Xuechao, SUI Yuhong. The deep tunnel rock blasting and crustal stress research of west route project of water transfer from south to north[J]. Journal of North China Institute of Water Conservancy and Hyadroelectric Power, 2002, 23(1): 48-50.