0 引 言
自引力是一种弱相互作用力，其地面实测可采用高精度扭秤系统实现[10]，但对测量环境以及装置要求较高，目前多通过仿真建模来计算。Nico Brandt等[11]采用端到端的建模方法，模拟检验质量周围的结构质量扰动，构建LISA Pathfinder自引力模型，计算检验质量受到的自引力扰动噪声。Merkowitz等[12]提出STOP-G模型，分析了航天器的自引力效应以满足LISA引力波探测计划的设计要求。Swank[13]给出了均匀密度直角平行六面体在一点上的自引力和自引力刚度的解析表达式，并以LISA检验质量为测试案例展开研究。国内对引力波的研究起步较晚，高志勇等[14]基于有限元方法计算了引力参考传感器中引力源对检验质量的引力效应，并提出“类自适应”网格划分方法以减小计算时间。林志勇[15]提出基于距离控制网格加密的自引力仿真方法，并使用立方体质量块和电极笼模型进行测试。
考虑到国内关于航天器整星建模来计算分析检验质量自引力的相关研究较少，本文依据“太极计划”空间引力波设计方案，构建航天器整星有限元模型，从而计算检验质量受到的自引力大小，分析检验质量近似为单质点、多质点、立方体、类立方体情况下的自引力计算偏差，并比较积分函数与理论解析公式的计算效果。
1 问题描述
TM1和TM2受自引力作用产生加速度噪声扰动，限制了空间引力波探测系统性能。结合空间引力波探测的基本原理，本文重点关注检验质量敏感轴即TM1 x1、TM2 x2方向，计算航天器对检验质量的引力作用，以服务于后续的自引力分析、配平、科学数据处理等实际需求。
2 自引力计算
2.1 整体流程
自引力计算流程如图2 所示，从航天器的实际构型出发，搭建几何模型，划分网格将结构离散成大小适当的有限元单元，实现航天器有限元模型的构建。从有限元模型中可提取各节点质量和位置参数信息，并根据引力计算公式计算每个节点对检验质量的引力，最后将所有节点上的自引力矢量求和，即可得出检验质量受到的总自引力。
2.2 检验质量自引力计算
以TM1为例，在其坐标系O1x1y1z1中，存在质量为M的点质量Q(x0, y0, z0)，计算点质量Q对检验质量的引力作用。考虑到自引力效应主要贡献检验质量加速度噪声，故直接计算该点质量对检验质量的引力加速度，即该点质量每单位质量上受到的自引力大小。根据对检验质量的不同近似，计算方式如下。
（1）若将TM1近似为单一质点，则点Q对检验质量的引力加速度为
式中，G为引力常数。
（2）若将TM1近似为多个等质量的质点，则点Q对检验质量的引力加速度为
式中，N为质点的个数；(xi, yi, zi)为第i个质点的坐标。
（3）若将TM1近似为长、宽、高分别为2c、2b、2a且密度均匀的长方体（多考虑为立方体），P(x, y, z)为检验质量内部任意一点，则点Q对检验质量的引力加速度为
式中，m为TM1的质量；ρ为TM1的均匀密度；V为TM1所占据的三维空间。
（4）根据引力波探测实际方案，检验质量并非规则的立方体，而是8个顶角缺失，z方向两个面上带有凹槽的类立方体，计算点Q对检验质量的引力加速度时只需要从立方体计算结果中扣除缺失的顶角和凹槽的引力贡献，点Q对检验质量的引力加速度为
式中，Fc为相应立方体受到的引力；Fi、Fj分别为第i个缺失顶角和第j个凹槽部分受到的引力。
若将检验质量视为单一或多个质点，根据上述公式直接计算即可；若将检验质量视为立方体或者类立方体，则计算过程中涉及多重积分求解。可以直接使用多重积分函数或者推导的理论解析公式计算，理论解析公式如下：
式中，ax为检验质量x方向受到的引力加速度。
y、z方向计算公式与x方向类似，这里不重复说明。
基于上述方法可以计算出空间内任一点质量对检验质量的引力加速度，结合有限元模型中提取的点质量信息，叠加计算即可得到检验质量受到的总自引力。
3 仿真实验
3.1 实验参数设置
本实验的有限元模型构建基于Hyperworks软件，计算过程基于MatlabR2020b软件实现，硬件仿真条件为i7-13700 16核并行运算，设置检验质量为边长46 mm，质量1.95 kg的均匀密度立方体，引力常数G取值6.67×10−11 （N∙m2）/kg2。
3.2 航天器有限元模型构建
本实验依据“太极计划”空间引力波探测器的设计方案，构建航天器整星有限元模型，如图3 所示。该模型包含了检验质量、太阳电池阵、气瓶、MOSA、固定机构等组件，从有限元模型中共提取出3161325 个点质量信息，包括点质量的质量和空间位置信息。
3.3 自引力计算结果
将检验质量近似为立方体，通过理论解析公式求解得出两检验质量的自引力计算结果，见表1 。
表 1 自引力计算结果Table 1 Self-gravity calculation results |
| 编号 | x方向 | y方向 | z方向 |
| TM1 | 6.14×10−9 | 6.86×10−9 | −7.14×10−9 |
| TM2 | 4.95×10−9 | −5.79×10−9 | −7.63×10−9 |
3.4 检验质量近似不同质点数目的计算偏差
如图4 所示，近似检验质量为等质量质点的个数分别为13、23、33、43、53，我们分别计算了两个检验质量敏感轴即TM1 x1、TM2 x2方向上受到的自引力加速度，计算偏差如图5 所示。
将检验质量近似为单个质点时，两检验质量敏感轴方向自引力计算偏差均大于5×10−14 m/s2，随着质点数目的增加，自引力计算偏差逐渐减小，当质点数N取53时，两检验质量敏感轴方向自引力计算偏差均小于1×10−15 m/s2。此外，质点数目的增加也会导致总的计算时间的增加，如图6 所示。
3.5 顶角缺失及凹槽对自引力计算的偏差
假设检验质量缺失的8个顶角形状均为边长为4 mm的立方体，凹槽为长8 mm、宽8 mm、高4 mm的立方体。计算该情况下检验质量受到的自引力并与近似为立方体的情况进行对比，其中类立方体与立方体质量相等，计算偏差见表2 。
表 2 立方体和类立方体的自引力计算偏差Table 2 Calculation deviation of self-gravity for cube and quasi-cube |
| 编号 | 类别 | x方向 | y方向 | z方向 |
| TM1 | 偏差/m∙s−2 | 7.94×10−13 | 1.48×10−14 | 9.77×10−14 |
| 百分比/% | 0.0129 | 0.0002 | 0.0014 | |
| TM2 | 偏差/m∙s−2 | 7.83×10−13 | 5.72×10−15 | 8.89×10−14 |
| 百分比/% | 0.0158 | 0.0001 | 0.0012 |
由于实际缺失的顶角和凹槽相较于检验质量的整体尺寸较小，且存在一定的对称关系，所以将检验质量近似为规则的立方体所计算的自引力偏差较小，两检验质量各轴方向自引力加速度大小计算偏差均小于8×10−13 m/s2，相对计算偏差均小于0.02%。
3.6 积分函数和理论解析公式求解的计算偏差
将检验质量视为立方体或者类立方体的计算过程涉及多重积分运算，可直接调用MATLAB软件中的积分函数，将计算结果与理论解析公式的计算结果进行对比，两者计算偏差见表3 。
表 3 积分函数与理论解析公式的自引力计算偏差Table 3 Calculation deviation of self-gravity for integral function and theoretical analytic formula |
| 编号 | 类别 | x方向 | y方向 | z方向 |
| TM1 | 偏差/m∙s−2 | 2.78×10−21 | 1.89×10−19 | 1.45×10−20 |
| 百分比/% | 4.53×10−11 | 2.76×10−9 | 2.03×10−10 | |
| TM2 | 偏差/m∙s−2 | 1.51×10−19 | 6.97×10−20 | 5.40×10−20 |
| 百分比/% | 3.06×10−9 | 1.20×10−9 | 7.08×10−10 |
从表中结果可知，积分函数与理论解析公式计算的自引力各轴计算相差均不超过10−18 m/s2，计算偏差满足实际使用需求，表明理论解析公式计算结果有效。同时，我们比较了积分函数和理论解析公式计算自引力所需的时间，如图7 所示。与积分函数相比，理论解析公式计算速度显著提升，极大地节省了计算资源。
综上所述，基于有限元方法计算自引力时，若将检验质量视为单一或多个质点，则计算精度受质点数量限制，对于自引力分析、量级估计等精度要求不高的任务，可以采用该方法。然而，对于自引力精确配平、科学数据处理等高精度需求而言，则需要考虑检验质量的实际构型，通过使用积分函数或理论解析公式进行计算，这两种方法精度较高，计算偏差小且易于实现。但是与积分函数计算相比，理论解析公式计算速度更快。此外，考虑到检验质量实际上是顶角缺失且带有凹槽的类立方体，若将其视为规则立方体进行计算则会存在偏差，但偏差相对数值较小。因此在计算检验质量受到的自引力时，需要权衡精度要求、计算效率以及实现难度等因素，根据具体情况选择合适的方法。
4 结语
本文针对空间引力波探测自引力需求，依据“太极计划”空间引力波设计方案，构建了航天器整星有限元模型，从而计算检验质量受到的自引力大小，分析检验质量不同近似状态下的自引力计算偏差，为工程实践提供参考。此外，比较了积分函数与理论解析公式的计算效果，两种方法计算结果相差不超过10−18 m/s2，理论解析公式计算精度可靠且速度更快，可有效应用于后续自引力计算。
后续研究将进一步完善有限元模型，提取充足的点质量信息，在更强大的仿真算力支持下实现更高精度、更加全面的自引力计算及分析。