张耀锋，尹 远，张春雷，蒋文丽，苏 俊

(北京师范大学 核科学与技术学院，北京 100875)

时间投影电离室(TPC)是典型的气体探测器，依据电离信号收集与放大部件的不同，TPC主要可分为丝状电极TPC探测器和气体电子倍增膜(GEM)类型的TPC探测器[1-3]。在TPC探测器设计及后续实验阶段，对探测器进行模拟研究，在加快设计进度、获得优化实验参数、提高实验效率等方面具有重要意义。GARFIELD程序是由欧洲核子研究中心(CERN)开发的适用于气体探测器的模拟程序，为当前国内外主流气体探测器模拟程序之一[4]。该程序主要采用有限边界法，实现对具有细丝电极的气体探测器进行电场计算以及后续的模拟优化。对于细丝电极的电场计算，GARFIELD程序仅适用于具有无限长细丝结构的二维探测器模型。此外，GARFIELD程序也可读入外部有限元程序计算得到的电场数据，以进行后续的探测器模拟计算。读入外部程序计算数据时往往受到许多限制，如对于ANSYS有限元电场计算数据，二维模型只能采用plane 121元素单元，三维模型只能采用solid 123元素单元。在GARFIELD程序读入有限元电场数据后，还需进行插值计算，以便进行数据显示与后续粒子漂移模拟。现有GARFIELD程序中插值部分存在数据完备性检验过于严格致使插值计算失败、电场数据经插值后连续性不好、部分位置无法输出有效电场数据等问题，导致程序对探测器的模拟无法进行。

有限元程序ANSYS常用于包含细丝电极结构的探测器的电场模拟。采用solid 122元素单元以及扫掠方式进行探测器主体结构模型构建，可有效避免细丝结构由于大的长宽比(长度与直径比)带来的因元素数目过多导致的计算失败的问题。鉴于GARFIELD模拟程序不能直接读入由solid 122元素单元计算获得的电场数据，本文拟对GARFIELD程序进行修改编译，并引入Quadratic Shepard插值法，以获得连续性很好的电场分布数据，用于后续的探测器模拟。

1 细丝结构探测器的3D电场计算

常见的TPC探测器采用数百根甚至上千根不同直径的细丝电极作为探测器的信号收集装置，可实现电离信号的雪崩放大以及收集测量[5]。细丝电极的主要特点是长度大、直径小，具有大的长宽比。典型的电极细丝长约1 m、直径为25 μm，其长宽比大于104。大的长宽比给细丝电极的3D电场计算带来了巨大挑战。以有限元程序ANSYS为例，用于静电场计算的3D元素单元主要有solid 122和solid 123两类[6]，如图1所示。可见solid 123元素单元为基本的三菱锥体类型，solid 122元素单元为类四棱柱类型。

图1 ANSYS用于3D电场计算的元素单元结构Fig.1 Structure of element type for 3D electric field calculation in ANSYS

图2 用于丝状电极的两种元素单元网格划分方法Fig.2 Two meshing ways of element type for wire electrode

这两种元素单元常通过不同的网格划分方法得到。如图2a所示，通过自由网格划分方法可得到solid 123元素单元，而通过扫掠式网格划分方法则得到solid 122元素单元。在细丝结构的网格划分过程中，由于三棱锥体表面方向的不规则性，致使细丝结构整体划分网格的元素基本尺寸较小，约为细丝直径的1/10。因此，若采用solid 123元素单元进行细丝结构的网格模型构建，会造成整体元素数目过多。对于包含数千根细丝的探测器模型，总的元素数目更会急剧增加，导致网格划分失败或电场计算失败。而采用扫掠方式进行细丝结构划分，细丝纵向solid 122元素单元的基本尺寸与截面位置的元素基本尺寸无关，可大幅减小细丝纵向上的元素网格密度，从而将探测器整体模型的元素数目控制在可接受的范围，以便实现探测器的整体3D电场计算。

由于探测器包含数目较多的细丝电极，常规的有限元电场计算方法在模型构建、网格划分、边界条件施加等过程中需要大量重复性的微区操作，极大地降低了整体电场的计算效率。利用有限元程序的脚本代码如ANSYS程序的APDL语言，可实现模拟程序的自动运行，从而提高模拟计算效率。但脚本语言的运行需预先对各电极丝进行准确的几何边界定位并进行特定单元的模拟运行核验，并在此基础上对细丝电极进行逐一建模以及网格划分。因此脚本语言的引入并不能提高太多的计算效率，只是方便了具有相似配置结构的探测器模型的重复计算。为进一步提高TPC多丝探测器的3D电场计算效率，采用节点组合法进行电场计算更为有效。该方法首先构建不包含任何细丝单元的宏观探测器结构模型，并进行网格划分。然后确定包含所有细丝电极的微小区域，对区域内的网格进行多次加密，直至元素节点分布能完整表征细丝截面形状。最后依据细丝电极的位置和半径信息，采用APDL脚本语言自动选取边界对应位置的节点，进行边界条件施加。对单丝模型采用节点组合法处理并进行边界条件施加后模型截面节点分布结果和采用常规扫掠式网格划分及边界条件施加后截面节点分布结果示于图3。由图3可知，在适当进行节点加密后，可得到与常规方法几乎完全一致的细丝单元的截面节点分布。这种逐步加密的方法，使得节点分布沿径向向外存在一定规律的节点密度衰减，这种不同的节点密度分布与细丝电极径向电场分布的指数衰减规律基本保持一致，从而有助于使用有限节点数获得较准确的计算结果。

图3 节点组合法(a)和常规扫掠式网格划分法(b)在边界施加后的细丝截面节点分布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插值法优化应用

2.1 有限元数据读入

在采用solid 122元素单元对包含细丝结构的TPC探测器进行3D电场计算后，可获得节点位置信息以及电势分布数据。为使GARFIELD程序能读入电场数据，对GARFIELD程序进行了修改编译，以便能使用该程序通过电场计算函数得到任意位置处电场矢量及电势的值，具体过程如下：1) 在程序即将进行电场计算时，将程序代码中的ICTYPE变量直接赋值为0，使得程序能转向处理有限元数据的代码；2) 针对处理有限元数据的代码，增加能读入有限元电场计算数据的相应代码，并进行插值计算，最终得到任意位置的电场和电势数据，返回指定函数。所有新增加代码采用Fortran语言编写，并与GARFIELD程序进行交叉编译。编译完成后，程序可正常运行。

2.2 Quadratic Shepard插值法应用

GARFIELD程序需给出任意位置的电场数据，以进行探测器模拟。由于ANSYS有限元程序计算得到的结果为所有节点处的数据，因此在数据读入GARFIELD程序后需采用插值法进行数据处理。GARFIELD程序原始的插值方法为双线性插值法，选取当前位置所处的体元素的所有节点为插值数据来源。如对于solid 123元素单元，则插值数据最多选取10个节点。如此小数目的插值数据源会给插值结果带来影响。如数据源点空间分布较均匀，则插值结果较为准确。反之，则计算结果误差较大。而ANSYS对于细丝电极电场的计算采用节点加密方法实现，因此节点数目密度分布非常不均匀。如，在接近细丝表面处，节点数目密度较大；而远离细丝表面的区域，节点数目密度较小。这两处区域内的节点密度相差数十倍到上百倍。相应的电场计算数据如采用GARFIELD程序原始的插值代码会导致完备性检验失败而无法输出结果或输出的结果连续性差。因此，针对细丝电极的电场分布特点，需引入新的插值方法，以获得连续性较好的电场分布结果。

图4 solid 122元素单元对应的插值处理方法Fig.4 Interpolation method for data from solid 122 element type

Quadratic Shepard插值法采用距离倒数加权的最小二乘法实现插值计算[7]，该方法一般不限定插值所需的数据源点的数目，可使用多个参数进行插值的圆滑性调节，最终获得连续性分布较好的插值结果。在实际插值应用过程中，常选取插值点临近区域的数十至上百个节点进行插值。依据电场计算时采用的元素单元类型的不同，分别采用不同的插值数据点选取方式。对于solid 122元素单元构建的模型，由于在纵向上采用扫掠方式划分网格，在纵向上分布的节点数密度相对于截面方向小得多。可采用图4所示的处理方法，插值点位于(x0,y0,z0)处，先在与插值节点相邻的两个横向平面z=z1和z=z2上进行二维插值，得到相应的两组插值结果f(x0,y0,z1)和f(x0,y0,z2)，然后在纵向上进行一维线性插值，则最终结果为：

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

(1)

对于solid 123元素单元构建的模型，可直接采用Quadratic Shepard的3D插值法，选取所插值节点附近适当数目的原始数据点进行插值计算，最终得到准确的插值结果。

2.3 插值结果分析

图5 电势沿径向分布的插值计算结果与理论计算结果Fig.5 Interpolation and theory result of potential distribution along radial direction

为检验Quadratic Shepard插值法的准确性，选取单丝电极结构模型进行电场有限元计算及插值结果分析。细丝电极直径为20 μm，长度为60 cm，施加电势为1 000 V。在距中心10 cm处同心圆柱外表面设为0 V电势区域，该区域纵向(z方向)长度为100 cm。对该模型采用节点组合方法进行3D电场有限元计算，之后将数据读入GARFIELD程序实现Quadratic Shepard插值计算以及电场数据输出。在模型中心平面处，电势沿径向分布的插值计算结果与理论计算结果的对比示于图5，可见两者符合很好，最大相对偏差在0.5%以内。此外，还对该平面处电场强度沿径向分布进行分析比较，以便对电势分布的连续性进行评价，结果示于图6。由图6可见，由插值计算得到的电场结果与理论计算结果符合较好，两者之间的相对偏差较图5有所增大，最大相对偏差在3%以内。其中偏差较大的位置均对应连续加密的分界区域，在该区域内由于节点加密导致节点空间分布呈现严重偏心表象。由于Quadratic Shepard插值法在计算过程中需选取计算目标点周围一定数目的节点数据，而这些选取的节点数据分布的空间不均匀性是影响插值结果准确程度的重要因素。上述插值计算过程中，实际选取其周围距离最近的120个数据点进行插值。对于某一数据点，采用最小距离法所选取的原始节点空间分布如图7a所示，数据分布的偏心程度异常严重。为改善插值结果，对最小距离法进行修正，加入修正因子，如式(2)所示：

图6 电场强度沿径向分布的插值计算结果与理论计算结果(a)及比较(b)Fig.6 Interpolation and theory calculation result of electric field strength radial distribution (a) and comparison (b)

十字标线中心为所需计算的数据位置图7 最小距离法(a)和最小距离+修正因子法(b)所得节点分布Fig.7 Node distribution for minimum distance method (a) and modified minimum distance method (b)

(2)

式中：(x0,y0)为所需要插值计算的数据点坐标；ri为每个节点与数据点(x0,y0)之间的距离；mi为每个节点的修正因子。mi的计算方法为：选取距离该节点最近的周围n个节点，取其中的最大距离。将最终获得的fi最小的120个数据点作为插值点阵列。其中n是优化的参数。如果n=30，则图7a的节点选取结果改变为图7b。此时插值计算结果的相对偏差由2.9%减小为0.17%。在实际插值计算中，可针对不同位置的节点选取不同的优化参数n，来获得较准确的插值结果。

采用上述电场数据插值计算方法，对包含多丝电极结构的TPC探测器进行3D电场的有限元计算。在中心平面处分别进行GARFIELD程序原始插值方法的电场计算输出和Quadratic Shepard插值法的电场计算输出，结果示于图8。TPC探测器门电极细丝位于y=-1.4 cm处，相邻细丝x方向间距0.1 cm。细丝电极电势施加按照计数分为奇、偶两组，分别施加Vg+V0和Vg-V0电势。其中Vg=-115 V、V0=30 V。Vg为门电极细丝处于测量打开状态下施加的基础电势，而±V0为探测器向测量关闭状态下转换时在电极细丝上按照奇偶两组分别施加的偏压电势。两种插值方法所得结果的差别主要表现在细丝电极临近区域。GARFIELD程序原始插值方法计算数据分布连续性较差，部分区域发生错误输出，会导致后续模拟程序中断；而Quadratic Shepard插值法对于电极临近区域能获得连续性较好的数据分布，有助于获得准确的探测器模拟结果。图8结果表明，Quadratic Shepard插值法是一种适用于细丝电极结构电场有限元程序计算的有效的数据后处理方法。

图8 TPC探测器门电极细丝附近区域等电势分布输出结果Fig.8 Output of equal potential distribution near gating wire of TPC detector

3 结论

利用有限元程序ANSYS并结合节点组合方法实现了包含细丝电极结构的TPC探测器的3D电场计算。在GARFIELD模拟程序读入电场数据后，采用Quadratic Shepard插值法对电场数据进行插值计算，在周围插值点分布极端不均匀情况下获得了较为准确的结果。Quadratic Shepard插值法有助于实现GARFIELD程序对包含细丝结构的TPC探测器的3维尺度的准确模拟。