張耀鋒，尹 遠(yuǎn)，張春雷，蔣文麗，蘇 俊

(北京師范大學(xué) 核科學(xué)與技術(shù)學(xué)院，北京 100875)

時間投影電離室(TPC)是典型的氣體探測器，依據(jù)電離信號收集與放大部件的不同，TPC主要可分為絲狀電極TPC探測器和氣體電子倍增膜(GEM)類型的TPC探測器[1-3]。在TPC探測器設(shè)計及后續(xù)實驗階段，對探測器進行模擬研究，在加快設(shè)計進度、獲得優(yōu)化實驗參數(shù)、提高實驗效率等方面具有重要意義。GARFIELD程序是由歐洲核子研究中心(CERN)開發(fā)的適用于氣體探測器的模擬程序，為當(dāng)前國內(nèi)外主流氣體探測器模擬程序之一[4]。該程序主要采用有限邊界法，實現(xiàn)對具有細(xì)絲電極的氣體探測器進行電場計算以及后續(xù)的模擬優(yōu)化。對于細(xì)絲電極的電場計算，GARFIELD程序僅適用于具有無限長細(xì)絲結(jié)構(gòu)的二維探測器模型。此外，GARFIELD程序也可讀入外部有限元程序計算得到的電場數(shù)據(jù)，以進行后續(xù)的探測器模擬計算。讀入外部程序計算數(shù)據(jù)時往往受到許多限制，如對于ANSYS有限元電場計算數(shù)據(jù)，二維模型只能采用plane 121元素單元，三維模型只能采用solid 123元素單元。在GARFIELD程序讀入有限元電場數(shù)據(jù)后，還需進行插值計算，以便進行數(shù)據(jù)顯示與后續(xù)粒子漂移模擬?，F(xiàn)有GARFIELD程序中插值部分存在數(shù)據(jù)完備性檢驗過于嚴(yán)格致使插值計算失敗、電場數(shù)據(jù)經(jīng)插值后連續(xù)性不好、部分位置無法輸出有效電場數(shù)據(jù)等問題，導(dǎo)致程序?qū)μ綔y器的模擬無法進行。

有限元程序ANSYS常用于包含細(xì)絲電極結(jié)構(gòu)的探測器的電場模擬。采用solid 122元素單元以及掃掠方式進行探測器主體結(jié)構(gòu)模型構(gòu)建，可有效避免細(xì)絲結(jié)構(gòu)由于大的長寬比(長度與直徑比)帶來的因元素數(shù)目過多導(dǎo)致的計算失敗的問題。鑒于GARFIELD模擬程序不能直接讀入由solid 122元素單元計算獲得的電場數(shù)據(jù)，本文擬對GARFIELD程序進行修改編譯，并引入Quadratic Shepard插值法，以獲得連續(xù)性很好的電場分布數(shù)據(jù)，用于后續(xù)的探測器模擬。

1 細(xì)絲結(jié)構(gòu)探測器的3D電場計算

常見的TPC探測器采用數(shù)百根甚至上千根不同直徑的細(xì)絲電極作為探測器的信號收集裝置，可實現(xiàn)電離信號的雪崩放大以及收集測量[5]。細(xì)絲電極的主要特點是長度大、直徑小，具有大的長寬比。典型的電極細(xì)絲長約1 m、直徑為25 μm，其長寬比大于104。大的長寬比給細(xì)絲電極的3D電場計算帶來了巨大挑戰(zhàn)。以有限元程序ANSYS為例，用于靜電場計算的3D元素單元主要有solid 122和solid 123兩類[6]，如圖1所示?？梢妔olid 123元素單元為基本的三菱錐體類型，solid 122元素單元為類四棱柱類型。

圖1 ANSYS用于3D電場計算的元素單元結(jié)構(gòu)Fig.1 Structure of element type for 3D electric field calculation in ANSYS

圖2 用于絲狀電極的兩種元素單元網(wǎng)格劃分方法Fig.2 Two meshing ways of element type for wire electrode

這兩種元素單元常通過不同的網(wǎng)格劃分方法得到。如圖2a所示，通過自由網(wǎng)格劃分方法可得到solid 123元素單元，而通過掃掠式網(wǎng)格劃分方法則得到solid 122元素單元。在細(xì)絲結(jié)構(gòu)的網(wǎng)格劃分過程中，由于三棱錐體表面方向的不規(guī)則性，致使細(xì)絲結(jié)構(gòu)整體劃分網(wǎng)格的元素基本尺寸較小，約為細(xì)絲直徑的1/10。因此，若采用solid 123元素單元進行細(xì)絲結(jié)構(gòu)的網(wǎng)格模型構(gòu)建，會造成整體元素數(shù)目過多。對于包含數(shù)千根細(xì)絲的探測器模型，總的元素數(shù)目更會急劇增加，導(dǎo)致網(wǎng)格劃分失敗或電場計算失敗。而采用掃掠方式進行細(xì)絲結(jié)構(gòu)劃分，細(xì)絲縱向solid 122元素單元的基本尺寸與截面位置的元素基本尺寸無關(guān)，可大幅減小細(xì)絲縱向上的元素網(wǎng)格密度，從而將探測器整體模型的元素數(shù)目控制在可接受的范圍，以便實現(xiàn)探測器的整體3D電場計算。

由于探測器包含數(shù)目較多的細(xì)絲電極，常規(guī)的有限元電場計算方法在模型構(gòu)建、網(wǎng)格劃分、邊界條件施加等過程中需要大量重復(fù)性的微區(qū)操作，極大地降低了整體電場的計算效率。利用有限元程序的腳本代碼如ANSYS程序的APDL語言，可實現(xiàn)模擬程序的自動運行，從而提高模擬計算效率。但腳本語言的運行需預(yù)先對各電極絲進行準(zhǔn)確的幾何邊界定位并進行特定單元的模擬運行核驗，并在此基礎(chǔ)上對細(xì)絲電極進行逐一建模以及網(wǎng)格劃分。因此腳本語言的引入并不能提高太多的計算效率，只是方便了具有相似配置結(jié)構(gòu)的探測器模型的重復(fù)計算。為進一步提高TPC多絲探測器的3D電場計算效率，采用節(jié)點組合法進行電場計算更為有效。該方法首先構(gòu)建不包含任何細(xì)絲單元的宏觀探測器結(jié)構(gòu)模型，并進行網(wǎng)格劃分。然后確定包含所有細(xì)絲電極的微小區(qū)域，對區(qū)域內(nèi)的網(wǎng)格進行多次加密，直至元素節(jié)點分布能完整表征細(xì)絲截面形狀。最后依據(jù)細(xì)絲電極的位置和半徑信息，采用APDL腳本語言自動選取邊界對應(yīng)位置的節(jié)點，進行邊界條件施加。對單絲模型采用節(jié)點組合法處理并進行邊界條件施加后模型截面節(jié)點分布結(jié)果和采用常規(guī)掃掠式網(wǎng)格劃分及邊界條件施加后截面節(jié)點分布結(jié)果示于圖3。由圖3可知，在適當(dāng)進行節(jié)點加密后，可得到與常規(guī)方法幾乎完全一致的細(xì)絲單元的截面節(jié)點分布。這種逐步加密的方法，使得節(jié)點分布沿徑向向外存在一定規(guī)律的節(jié)點密度衰減，這種不同的節(jié)點密度分布與細(xì)絲電極徑向電場分布的指數(shù)衰減規(guī)律基本保持一致，從而有助于使用有限節(jié)點數(shù)獲得較準(zhǔn)確的計算結(jié)果。

圖3 節(jié)點組合法(a)和常規(guī)掃掠式網(wǎng)格劃分法(b)在邊界施加后的細(xì)絲截面節(jié)點分布Fig.3 Distribution of node for wire section after boundary condition applied by node-construction method (a) and regular sweeping method (b)

2 Quadratic Shepard插值法優(yōu)化應(yīng)用

2.1 有限元數(shù)據(jù)讀入

在采用solid 122元素單元對包含細(xì)絲結(jié)構(gòu)的TPC探測器進行3D電場計算后，可獲得節(jié)點位置信息以及電勢分布數(shù)據(jù)。為使GARFIELD程序能讀入電場數(shù)據(jù)，對GARFIELD程序進行了修改編譯，以便能使用該程序通過電場計算函數(shù)得到任意位置處電場矢量及電勢的值，具體過程如下：1) 在程序即將進行電場計算時，將程序代碼中的ICTYPE變量直接賦值為0，使得程序能轉(zhuǎn)向處理有限元數(shù)據(jù)的代碼；2) 針對處理有限元數(shù)據(jù)的代碼，增加能讀入有限元電場計算數(shù)據(jù)的相應(yīng)代碼，并進行插值計算，最終得到任意位置的電場和電勢數(shù)據(jù)，返回指定函數(shù)。所有新增加代碼采用Fortran語言編寫，并與GARFIELD程序進行交叉編譯。編譯完成后，程序可正常運行。

2.2 Quadratic Shepard插值法應(yīng)用

GARFIELD程序需給出任意位置的電場數(shù)據(jù)，以進行探測器模擬。由于ANSYS有限元程序計算得到的結(jié)果為所有節(jié)點處的數(shù)據(jù)，因此在數(shù)據(jù)讀入GARFIELD程序后需采用插值法進行數(shù)據(jù)處理。GARFIELD程序原始的插值方法為雙線性插值法，選取當(dāng)前位置所處的體元素的所有節(jié)點為插值數(shù)據(jù)來源。如對于solid 123元素單元，則插值數(shù)據(jù)最多選取10個節(jié)點。如此小數(shù)目的插值數(shù)據(jù)源會給插值結(jié)果帶來影響。如數(shù)據(jù)源點空間分布較均勻，則插值結(jié)果較為準(zhǔn)確。反之，則計算結(jié)果誤差較大。而ANSYS對于細(xì)絲電極電場的計算采用節(jié)點加密方法實現(xiàn)，因此節(jié)點數(shù)目密度分布非常不均勻。如，在接近細(xì)絲表面處，節(jié)點數(shù)目密度較大；而遠(yuǎn)離細(xì)絲表面的區(qū)域，節(jié)點數(shù)目密度較小。這兩處區(qū)域內(nèi)的節(jié)點密度相差數(shù)十倍到上百倍。相應(yīng)的電場計算數(shù)據(jù)如采用GARFIELD程序原始的插值代碼會導(dǎo)致完備性檢驗失敗而無法輸出結(jié)果或輸出的結(jié)果連續(xù)性差。因此，針對細(xì)絲電極的電場分布特點，需引入新的插值方法，以獲得連續(xù)性較好的電場分布結(jié)果。

圖4 solid 122元素單元對應(yīng)的插值處理方法Fig.4 Interpolation method for data from solid 122 element type

Quadratic Shepard插值法采用距離倒數(shù)加權(quán)的最小二乘法實現(xiàn)插值計算[7]，該方法一般不限定插值所需的數(shù)據(jù)源點的數(shù)目，可使用多個參數(shù)進行插值的圓滑性調(diào)節(jié)，最終獲得連續(xù)性分布較好的插值結(jié)果。在實際插值應(yīng)用過程中，常選取插值點臨近區(qū)域的數(shù)十至上百個節(jié)點進行插值。依據(jù)電場計算時采用的元素單元類型的不同，分別采用不同的插值數(shù)據(jù)點選取方式。對于solid 122元素單元構(gòu)建的模型，由于在縱向上采用掃掠方式劃分網(wǎng)格，在縱向上分布的節(jié)點數(shù)密度相對于截面方向小得多?？刹捎脠D4所示的處理方法，插值點位于(x0,y0,z0)處，先在與插值節(jié)點相鄰的兩個橫向平面z=z1和z=z2上進行二維插值，得到相應(yīng)的兩組插值結(jié)果f(x0,y0,z1)和f(x0,y0,z2)，然后在縱向上進行一維線性插值，則最終結(jié)果為：

f(x0,y0,z0)=f(x0,y0,z1)+

(1)

對于solid 123元素單元構(gòu)建的模型，可直接采用Quadratic Shepard的3D插值法，選取所插值節(jié)點附近適當(dāng)數(shù)目的原始數(shù)據(jù)點進行插值計算，最終得到準(zhǔn)確的插值結(jié)果。

2.3 插值結(jié)果分析

圖5 電勢沿徑向分布的插值計算結(jié)果與理論計算結(jié)果Fig.5 Interpolation and theory result of potential distribution along radial direction

為檢驗Quadratic Shepard插值法的準(zhǔn)確性，選取單絲電極結(jié)構(gòu)模型進行電場有限元計算及插值結(jié)果分析。細(xì)絲電極直徑為20 μm，長度為60 cm，施加電勢為1 000 V。在距中心10 cm處同心圓柱外表面設(shè)為0 V電勢區(qū)域，該區(qū)域縱向(z方向)長度為100 cm。對該模型采用節(jié)點組合方法進行3D電場有限元計算，之后將數(shù)據(jù)讀入GARFIELD程序?qū)崿F(xiàn)Quadratic Shepard插值計算以及電場數(shù)據(jù)輸出。在模型中心平面處，電勢沿徑向分布的插值計算結(jié)果與理論計算結(jié)果的對比示于圖5，可見兩者符合很好，最大相對偏差在0.5%以內(nèi)。此外，還對該平面處電場強度沿徑向分布進行分析比較，以便對電勢分布的連續(xù)性進行評價，結(jié)果示于圖6。由圖6可見，由插值計算得到的電場結(jié)果與理論計算結(jié)果符合較好，兩者之間的相對偏差較圖5有所增大，最大相對偏差在3%以內(nèi)。其中偏差較大的位置均對應(yīng)連續(xù)加密的分界區(qū)域，在該區(qū)域內(nèi)由于節(jié)點加密導(dǎo)致節(jié)點空間分布呈現(xiàn)嚴(yán)重偏心表象。由于Quadratic Shepard插值法在計算過程中需選取計算目標(biāo)點周圍一定數(shù)目的節(jié)點數(shù)據(jù)，而這些選取的節(jié)點數(shù)據(jù)分布的空間不均勻性是影響插值結(jié)果準(zhǔn)確程度的重要因素。上述插值計算過程中，實際選取其周圍距離最近的120個數(shù)據(jù)點進行插值。對于某一數(shù)據(jù)點，采用最小距離法所選取的原始節(jié)點空間分布如圖7a所示，數(shù)據(jù)分布的偏心程度異常嚴(yán)重。為改善插值結(jié)果，對最小距離法進行修正，加入修正因子，如式(2)所示：

圖6 電場強度沿徑向分布的插值計算結(jié)果與理論計算結(jié)果(a)及比較(b)Fig.6 Interpolation and theory calculation result of electric field strength radial distribution (a) and comparison (b)

十字標(biāo)線中心為所需計算的數(shù)據(jù)位置圖7 最小距離法(a)和最小距離+修正因子法(b)所得節(jié)點分布Fig.7 Node distribution for minimum distance method (a) and modified minimum distance method (b)

(2)

式中：(x0,y0)為所需要插值計算的數(shù)據(jù)點坐標(biāo)；ri為每個節(jié)點與數(shù)據(jù)點(x0,y0)之間的距離；mi為每個節(jié)點的修正因子。mi的計算方法為：選取距離該節(jié)點最近的周圍n個節(jié)點，取其中的最大距離。將最終獲得的fi最小的120個數(shù)據(jù)點作為插值點陣列。其中n是優(yōu)化的參數(shù)。如果n=30，則圖7a的節(jié)點選取結(jié)果改變?yōu)閳D7b。此時插值計算結(jié)果的相對偏差由2.9%減小為0.17%。在實際插值計算中，可針對不同位置的節(jié)點選取不同的優(yōu)化參數(shù)n，來獲得較準(zhǔn)確的插值結(jié)果。

采用上述電場數(shù)據(jù)插值計算方法，對包含多絲電極結(jié)構(gòu)的TPC探測器進行3D電場的有限元計算。在中心平面處分別進行GARFIELD程序原始插值方法的電場計算輸出和Quadratic Shepard插值法的電場計算輸出，結(jié)果示于圖8。TPC探測器門電極細(xì)絲位于y=-1.4 cm處，相鄰細(xì)絲x方向間距0.1 cm。細(xì)絲電極電勢施加按照計數(shù)分為奇、偶兩組，分別施加Vg+V0和Vg-V0電勢。其中Vg=-115 V、V0=30 V。Vg為門電極細(xì)絲處于測量打開狀態(tài)下施加的基礎(chǔ)電勢，而±V0為探測器向測量關(guān)閉狀態(tài)下轉(zhuǎn)換時在電極細(xì)絲上按照奇偶兩組分別施加的偏壓電勢。兩種插值方法所得結(jié)果的差別主要表現(xiàn)在細(xì)絲電極臨近區(qū)域。GARFIELD程序原始插值方法計算數(shù)據(jù)分布連續(xù)性較差，部分區(qū)域發(fā)生錯誤輸出，會導(dǎo)致后續(xù)模擬程序中斷；而Quadratic Shepard插值法對于電極臨近區(qū)域能獲得連續(xù)性較好的數(shù)據(jù)分布，有助于獲得準(zhǔn)確的探測器模擬結(jié)果。圖8結(jié)果表明，Quadratic Shepard插值法是一種適用于細(xì)絲電極結(jié)構(gòu)電場有限元程序計算的有效的數(shù)據(jù)后處理方法。

圖8 TPC探測器門電極細(xì)絲附近區(qū)域等電勢分布輸出結(jié)果Fig.8 Output of equal potential distribution near gating wire of TPC detector

3 結(jié)論

利用有限元程序ANSYS并結(jié)合節(jié)點組合方法實現(xiàn)了包含細(xì)絲電極結(jié)構(gòu)的TPC探測器的3D電場計算。在GARFIELD模擬程序讀入電場數(shù)據(jù)后，采用Quadratic Shepard插值法對電場數(shù)據(jù)進行插值計算，在周圍插值點分布極端不均勻情況下獲得了較為準(zhǔn)確的結(jié)果。Quadratic Shepard插值法有助于實現(xiàn)GARFIELD程序?qū)Π?xì)絲結(jié)構(gòu)的TPC探測器的3維尺度的準(zhǔn)確模擬。