Quadratic Shepard插值法在TPC探測(cè)器電場(chǎng)計(jì)算中的應(yīng)用

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

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

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

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

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

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

圖1 ANSYS用于3D電場(chǎng)計(jì)算的元素單元結(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

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

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

圖3 節(jié)點(diǎn)組合法(a)和常規(guī)掃掠式網(wǎng)格劃分法(b)在邊界施加后的細(xì)絲截面節(jié)點(diǎn)分布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元素單元對(duì)包含細(xì)絲結(jié)構(gòu)的TPC探測(cè)器進(jìn)行3D電場(chǎng)計(jì)算后，可獲得節(jié)點(diǎn)位置信息以及電勢(shì)分布數(shù)據(jù)。為使GARFIELD程序能讀入電場(chǎng)數(shù)據(jù)，對(duì)GARFIELD程序進(jìn)行了修改編譯，以便能使用該程序通過(guò)電場(chǎng)計(jì)算函數(shù)得到任意位置處電場(chǎng)矢量及電勢(shì)的值，具體過(guò)程如下：1) 在程序即將進(jìn)行電場(chǎng)計(jì)算時(shí)，將程序代碼中的ICTYPE變量直接賦值為0，使得程序能轉(zhuǎn)向處理有限元數(shù)據(jù)的代碼；2) 針對(duì)處理有限元數(shù)據(jù)的代碼，增加能讀入有限元電場(chǎng)計(jì)算數(shù)據(jù)的相應(yīng)代碼，并進(jìn)行插值計(jì)算，最終得到任意位置的電場(chǎng)和電勢(shì)數(shù)據(jù)，返回指定函數(shù)。所有新增加代碼采用Fortran語(yǔ)言編寫(xiě)，并與GARFIELD程序進(jìn)行交叉編譯。編譯完成后，程序可正常運(yùn)行。

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

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

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

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

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

(1)

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

2.3 插值結(jié)果分析

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

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

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

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

(2)

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

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

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

3 結(jié)論

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