空间科学与试验学报

空间科学与试验学报 >

2024 , Vol. 1 >Issue 2: 68 - 73

DOI: https://doi.org/10.19963/j.cnki.2097-4302.2024.02.008

空间科学

“太极计划”空间引力波探测器自引力计算与分析

  • 汤宁标 , 1, 2 ,
  • 杨中光 , 1, * ,
  • 余贤圣 1 ,
  • 何涛 1 ,
  • 蔡志鸣 1 ,
  • 余金培 1, 2 ,
  • 李华旺 1, 2
展开
  • 1. 中国科学院微小卫星创新研究院, 上海201304
  • 2. 中国科学院大学, 北京100049
（1989-），男，博士，副研究员，主要从事科学卫星总体设计等方面的研究。通信地址：上海市浦东新区雪洋路1号（201304）电子邮箱:

（2000-），男，硕士生，主要从事科学卫星总体设计等方面的研究。通信地址：上海市浦东新区雪洋路1号（201304）电子邮箱:

网络出版日期: 2024-11-25

基金资助

国家重点研发计划资助项目（2020YFC2200901）

版权

版权所有 © 2024 空间科学与试验学报编辑部
收起

Self-Gravity Calculation and Analysis of "Taiji" Space Gravitational Wave Detector

  • Ningbiao TANG , 1, 2 ,
  • Zhongguang YANG , 1, * ,
  • Xiansheng YU 1 ,
  • Tao HE 1 ,
  • Zhiming CAI 1 ,
  • Jinpei YU 1, 2 ,
  • Huawang LI 1, 2
Expand
  • 1. Innovation Academy for Microsatellites of Chinese Academy of Sciences, Shanghai 201304, China
  • 2. University of Chinese Academy of Sciences, Beijing 100049, China

Online published: 2024-11-25

Copyright

Copyright © 2024 Journal of Space Science and Experiment. All rights reserved.
Fold

摘要

本文针对空间引力波探测自引力需求，依据“太极计划”空间引力波探测器设计方案，构建航天器整星有限元模型，基于有限元方法计算检验质量受到的自引力大小，分析将检验质量近似为单质点、多质点、立方体和类立方体情况下的自引力计算偏差，比较积分函数与理论解析公式的计算效果。仿真结果表明：当检验质量被近似为单一或多个质点时，其计算精度受质点数量限制；当检验质量被近似为立方体和类立方体时，其计算精度较高且两者相对计算偏差小于0.02%；当使用积分函数与理论解析公式计算自引力时，两种方法计算结果相差不超过10−18 m/s2，其中理论解析公式计算速度更快，可有效应用于后续自引力计算。

关键词： 空间引力波探测; 自引力; 有限元方法; 计算偏差分析

本文引用格式

汤宁标 , 杨中光 , 余贤圣 , 何涛 , 蔡志鸣 , 余金培 , 李华旺 . “太极计划”空间引力波探测器自引力计算与分析[J]. 空间科学与试验学报, 2024 , 1(2) : 68 -73 . DOI: 10.19963/j.cnki.2097-4302.2024.02.008

Abstract

In view of the self-gravity requirements of space gravitational wave detection, we build a finite element model of the whole spacecraft based on the design scheme of "Taiji" space gravitational wave detector, and calculate the self-gravity of the test mass based on the finite element method. The self-gravity calculation deviation is analyzed when the test mass is approximated as a single particle, multiple particles, a cube and a quasi-cube, and the calculation effects of the integral function and the theoretical analytical formula is compared. The simulation results show that when the test mass is approximated as a single or multiple particles, the calculation accuracy is limited by the number of particles; when the test mass is approximated as a cube or a quasi-cube, the calculation accuracy is higher and the relative calculation deviation between the two methods is less than 0.02%; when using the integral function and the theoretical analytical formula to calculate self-gravity, the difference in the calculation results of the two methods does not exceed 10-18 m/s2. The theoretical analytical formula is faster in calculation and can be effectively used in subsequent self-gravity calculations.

Key words： space gravitational wave detection; self-gravity; finite element method; calculation deviation analysis

0 引 言

引力波是爱因斯坦广义相对论的重要预测之一，是物质和能量剧烈运动与变化所产生的一种物质波[1]。宇宙中存在大量的引力波源，覆盖了从1 aHz（1 aHz=10−18 Hz）到10 kHz的宽阔频段，不同频段的引力波源不同，对应的科学目标也不同。空间引力波探测将打开引力波毫赫兹频段窗口，在研究黑洞并合、宇宙起源等方面具有重大科学意义[2-3]。
引力波经过时会挤压或拉伸附近的时空，从而引起空间中不同位置间的光程变化。空间引力波探测利用高精度激光干涉仪，测量距离十万公里以上的两检验质量间皮米级的光程变化，从而反演引力波信息[4]。为有效检测到引力波信号，需要测量和抑制检验质量受到除引力波外的其他非保守力导致的扰动，以实现引力波探测所需的灵敏度。目前在研究的空间引力波探测计划主要有由美国航空航天局和欧洲航天局合作的激光干涉空间天线计划[5]（Laser Interferometer Space Antenna，LISA）、国内中国科学院主导的“太极计划”[6]以及中山大学主导的“天琴计划”[7]等。
空间引力波探测系统灵敏度要求检验质量的残余加速度噪声在探测频段内优于10−15 m∙s−2∙Hz−1/2量级，系统噪声指标要求极高。航天器与检验质量之间的引力即自引力[8]，自引力和航天器质量分布等因素密切相关，由自引力效应引起的加速度扰动是噪声预算的一项重要来源，要求优于10−10 m/s2量级[9]。
自引力是一种弱相互作用力，其地面实测可采用高精度扭秤系统实现[10]，但对测量环境以及装置要求较高，目前多通过仿真建模来计算。Nico Brandt等[11]采用端到端的建模方法，模拟检验质量周围的结构质量扰动，构建LISA Pathfinder自引力模型，计算检验质量受到的自引力扰动噪声。Merkowitz等[12]提出STOP-G模型，分析了航天器的自引力效应以满足LISA引力波探测计划的设计要求。Swank[13]给出了均匀密度直角平行六面体在一点上的自引力和自引力刚度的解析表达式，并以LISA检验质量为测试案例展开研究。国内对引力波的研究起步较晚，高志勇等[14]基于有限元方法计算了引力参考传感器中引力源对检验质量的引力效应，并提出“类自适应”网格划分方法以减小计算时间。林志勇[15]提出基于距离控制网格加密的自引力仿真方法，并使用立方体质量块和电极笼模型进行测试。
考虑到国内关于航天器整星建模来计算分析检验质量自引力的相关研究较少，本文依据“太极计划”空间引力波设计方案，构建航天器整星有限元模型，从而计算检验质量受到的自引力大小，分析检验质量近似为单质点、多质点、立方体、类立方体情况下的自引力计算偏差，并比较积分函数与理论解析公式的计算效果。

1 问题描述

图1为航天器结构，其中S/C表示航天器，对应坐标系为O0x0y0z0，坐标系原点O0位于航天器底板中心处；TM1和TM2表示两检验质量，对应坐标系分别为O1x1y1z1和O2x2y2z2，坐标系原点O1、O2位于两检验质量标称位置中心处；θ为移动光学组件（Moving Optical Sub-Assembly，MOSA）夹角。
图 1 航天器结构示意图

Fig.1 Spacecraft structure diagram

TM1和TM2受自引力作用产生加速度噪声扰动，限制了空间引力波探测系统性能。结合空间引力波探测的基本原理，本文重点关注检验质量敏感轴即TM1 x1、TM2 x2方向，计算航天器对检验质量的引力作用，以服务于后续的自引力分析、配平、科学数据处理等实际需求。

2 自引力计算

2.1 整体流程

自引力计算流程如图2所示，从航天器的实际构型出发，搭建几何模型，划分网格将结构离散成大小适当的有限元单元，实现航天器有限元模型的构建。从有限元模型中可提取各节点质量和位置参数信息，并根据引力计算公式计算每个节点对检验质量的引力，最后将所有节点上的自引力矢量求和，即可得出检验质量受到的总自引力。
图 2 自引力计算流程图

Fig.2 Flow chart of self-gravity calculation

2.2 检验质量自引力计算

以TM1为例，在其坐标系O1x1y1z1中，存在质量为M的点质量Q(x0, y0, z0)，计算点质量Q对检验质量的引力作用。考虑到自引力效应主要贡献检验质量加速度噪声，故直接计算该点质量对检验质量的引力加速度，即该点质量每单位质量上受到的自引力大小。根据对检验质量的不同近似，计算方式如下。
（1）若将TM1近似为单一质点，则点Q对检验质量的引力加速度为
$ \boldsymbol{a} = GM\frac{{({x_0}{\text{ }}{y_0}{\text{ }}{z_0})}}{{{{(x_0^2 + y_0^2 + z_0^2)}^{1.5}}}} $
式中，G为引力常数。
（2）若将TM1近似为多个等质量的质点，则点Q对检验质量的引力加速度为
$ \boldsymbol{a} = \sum\limits_{i = 1}^N {\frac{{GM}}{N}\frac{{({x_0} - {x_i}{\text{ }}{y_0} - {y_i}{\text{ }}{z_0} - {z_i})}}{{{{\left[ {{{({x_0} - {x_i})}^2} + {{({y_0} - {y_i})}^2} + {{({z_0} - {z_i})}^2}} \right]}^{1.5}}}}} $
式中，N为质点的个数；(xi, yi, zi)为第i个质点的坐标。
（3）若将TM1近似为长、宽、高分别为2c、2b、2a且密度均匀的长方体（多考虑为立方体），P(x, y, z)为检验质量内部任意一点，则点Q对检验质量的引力加速度为
$ {\boldsymbol{a} = \dfrac{{GM\rho }}{m}\iiint\limits_V {\dfrac{{({x_0} - x{\text{ }}{y_0} - y{\text{ }}{z_0} - z)}}{{{{\left[ {{{({x_0} - x)}^2} + {{({y_0} - y)}^2} + {{({z_0} - z)}^2}} \right]}^{1.5}}}}{\mathrm{d}}x{\mathrm{d}}y{\mathrm{d}}z} }$
式中，m为TM1的质量；ρ为TM1的均匀密度；V为TM1所占据的三维空间。
（4）根据引力波探测实际方案，检验质量并非规则的立方体，而是8个顶角缺失，z方向两个面上带有凹槽的类立方体，计算点Q对检验质量的引力加速度时只需要从立方体计算结果中扣除缺失的顶角和凹槽的引力贡献，点Q对检验质量的引力加速度为
$ \boldsymbol{a} = \frac{{{\boldsymbol{F}_{\mathrm{c}}} - \displaystyle\sum\limits_{i = 1}^8 {{\boldsymbol{F}_i} - \sum\limits_{j = 1}^2 {{\boldsymbol{F}_j}} } }}{m} $
式中，Fc为相应立方体受到的引力；Fi、Fj分别为第i个缺失顶角和第j个凹槽部分受到的引力。
若将检验质量视为单一或多个质点，根据上述公式直接计算即可；若将检验质量视为立方体或者类立方体，则计算过程中涉及多重积分求解。可以直接使用多重积分函数或者推导的理论解析公式计算，理论解析公式如下：
${\begin{aligned} &\dfrac{{ - 8abc}}{{GM}}{a_x} = \left( {{y_0} - b} \right)\ln \left( {\dfrac{{{z_0} + c + r\left( { - a, - b,c} \right)}}{{{z_0} - c + r\left( { - a, - b, - c} \right)}} \times \dfrac{{{z_0} - c + r\left( {a, - b, - c} \right)}}{{{z_0} + c + r\left( {a, - b,c} \right)}}} \right) - \\ &\ \ \left( {{y_0} + b} \right)\ln \left( {\dfrac{{{z_0} + c + r\left( { - a,b,c} \right)}}{{{z_0} - c + r\left( { - a,b, - c} \right)}} \times \dfrac{{{z_0} - c + r\left( {a,b, - c} \right)}}{{{z_0} + c + r\left( {a,b,c} \right)}}} \right) + \\ &\ \ \left( {{z_0} + c} \right)\ln \left( {\dfrac{{{y_0} + b + r\left( {a,b,c} \right)}}{{{y_0} + b + r\left( { - a,b,c} \right)}} \times \dfrac{{{y_0} - b + r\left( { - a, - b,c} \right)}}{{{y_0} - b + r\left( {a, - b,c} \right)}}} \right) - \\ &\ \ \left( {{z_0} - c} \right)\ln \left( {\dfrac{{{y_0} - b + r\left( { - a, - b, - c} \right)}}{{{y_0} - b + r\left( {a, - b, - c} \right)}} \times \dfrac{{{y_0} + b + r\left( {a,b, - c} \right)}}{{{y_0} + b + r\left( { - a,b, - c} \right)}}} \right) + \\ &\ \ \left( {{x_0} - a} \right)\left( \begin{gathered} \left[ {\arctan \left( {\dfrac{{{z_0} + c}}{{{x_0} - a}}\dfrac{{{y_0} + b}}{{r\left( { - a,b,c} \right)}}} \right) - \arctan \left( {\dfrac{{{z_0} + c}}{{{x_0} - a}}\dfrac{{{y_0} - b}}{{r\left( { - a, - b,c} \right)}}} \right)} \right] + \\\left[ {\arctan \left( {\dfrac{{{z_0} - c}}{{{x_0} - a}}\dfrac{{{y_0} - b}}{{r\left( { - a, - b, - c} \right)}}} \right) - \arctan \left( {\dfrac{{{z_0} - c}}{{{x_0} - a}}\dfrac{{{y_0} + b}}{{r\left( { - a,b, - c} \right)}}} \right)} \right] \\ \end{gathered} \right) + \\ &\ \ \left( {{x_0} + a} \right)\left( \begin{gathered} \left[ {\arctan \left( {\dfrac{{{z_0} + c}}{{{x_0} + a}}\dfrac{{{y_0} - b}}{{r\left( {a, - b,c} \right)}}} \right) - \arctan \left( {\dfrac{{{z_0} + c}}{{{x_0} + a}}\dfrac{{{y_0} + b}}{{r\left( {a,b,c} \right)}}} \right)} \right] + \\\left[ {\arctan \left( {\dfrac{{{z_0} - c}}{{{x_0} + a}}\dfrac{{{y_0} + b}}{{r\left( {a,b, - c} \right)}}} \right) - \arctan \left( {\dfrac{{{z_0} - c}}{{{x_0} + a}}\dfrac{{{y_0} - b}}{{r\left( {a, - b, - c} \right)}}} \right)} \right] \\ \end{gathered} \right) \\ \end{aligned} }$
式中，ax为检验质量x方向受到的引力加速度。
$ r\left( {a,b,c} \right) = \sqrt {{{\left( {{x_0} + a} \right)}^2} + {{\left( {{y_0} + b} \right)}^2} + {{\left( {{z_0} + c} \right)}^2}} $
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 航天器有限元模型

Fig.3 Spacecraft finite element model

3.3 自引力计算结果

将检验质量近似为立方体，通过理论解析公式求解得出两检验质量的自引力计算结果，见表1。
表 1 自引力计算结果

Table 1 Self-gravity calculation results （m/s2）

编号 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所示。
图 4 不同质点数目设置

Fig.4 Setting of different number of particles

图 5 不同质点数目的自引力计算偏差

Fig.5 Calculation deviation of self-gravity for different number of particles

将检验质量近似为单个质点时，两检验质量敏感轴方向自引力计算偏差均大于5×10−14 m/s2，随着质点数目的增加，自引力计算偏差逐渐减小，当质点数N取53时，两检验质量敏感轴方向自引力计算偏差均小于1×10−15 m/s2。此外，质点数目的增加也会导致总的计算时间的增加，如图6所示。
图 6 不同质点数目的计算时间

Fig.6 Calculation time for different number of particles

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所示。与积分函数相比，理论解析公式计算速度显著提升，极大地节省了计算资源。
图 7 积分函数和理论解析公式计算时间

Fig.7 Calculation time of the integral function and the theoretical analytic formula

综上所述，基于有限元方法计算自引力时，若将检验质量视为单一或多个质点，则计算精度受质点数量限制，对于自引力分析、量级估计等精度要求不高的任务，可以采用该方法。然而，对于自引力精确配平、科学数据处理等高精度需求而言，则需要考虑检验质量的实际构型，通过使用积分函数或理论解析公式进行计算，这两种方法精度较高，计算偏差小且易于实现。但是与积分函数计算相比，理论解析公式计算速度更快。此外，考虑到检验质量实际上是顶角缺失且带有凹槽的类立方体，若将其视为规则立方体进行计算则会存在偏差，但偏差相对数值较小。因此在计算检验质量受到的自引力时，需要权衡精度要求、计算效率以及实现难度等因素，根据具体情况选择合适的方法。

4 结语

本文针对空间引力波探测自引力需求，依据“太极计划”空间引力波设计方案，构建了航天器整星有限元模型，从而计算检验质量受到的自引力大小，分析检验质量不同近似状态下的自引力计算偏差，为工程实践提供参考。此外，比较了积分函数与理论解析公式的计算效果，两种方法计算结果相差不超过10−18 m/s2，理论解析公式计算精度可靠且速度更快，可有效应用于后续自引力计算。
后续研究将进一步完善有限元模型，提取充足的点质量信息，在更强大的仿真算力支持下实现更高精度、更加全面的自引力计算及分析。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
1
WEBER J. General relativity and gravitational waves[M]. Chicago:Courier Corporation,2004.

2
RILES K. Gravitational waves: Sources, detectors and searches[J]. Progress in Particle and Nuclear Physics, 2013, 68, 1- 54.

DOI

3
CAPRINI C, FIGUEROA D G. Cosmological backgrounds of gravitational waves[J]. Classical and Quantum Gravity, 2018, 35 (16): 163001.

DOI

4
JENNRICH O. LISA technology and instrumentation[J]. Classical and Quantum Gravity, 2009, 26 (15): 153001.

DOI

5
AMARO-SEOANE P, AUDLEY H, BABAK S, et al. Laser interferometer space antenna[J]. arXiv preprint arXiv, 2017, 0078 (6): 1702.

6
HU W R, WU Y L. The Taiji program in space for gravitational wave physics and the nature of gravity[J]. National Science Review, 2017, 4 (5): 685- 686.

DOI

7
LUO J, CHEN L S, DUAN H Z, et al. Tianqin: A space-borne gravitational wave detector[J]. Classical and Quantum Gravity, 2016, 33 (3): 035010.

DOI

8
GOPSTEIN A M, HAILE W B, MERKOWITZ S M. Self‐gravity analysis and visualization tool for LISA[J]. AIP Conference Proceedings, 2006, 873 (1): 571- 575.

9
STEBBINS R T, BENDER P L, HANSON J, et al. Current error estimates for LISA spurious accelerations[J]. Classical and Quantum Gravity, 2004, 21 (5): S653.

DOI

10
TRENKEL C, WARREN C, WEALTHY D. Ground-based self-gravity tests for LISA Pathfinder and LISA[J]. Journal of Physics: Conference Series, 2009, 154 (1): 012010.

11
BRANDT N, FICHTER W, KERSTEN M, et al. End-to-end modeling for drag-free missions with application to LISA pathfinder[J]. IFAC Proceedings Volumes, 2004, 37 (6): 241- 246.

DOI

12
MERKOWITZ S M, CONKEY S, HAILE W B, et al. Structural, thermal, optical and gravitational modelling for LISA[J]. Classical and Quantum Gravity, 2004, 21 (5): S603.

DOI

13
SWANK A J. Gravitational mass attraction: Properties of a right-angled parallelepiped for the LISA drag-free system[J]. Classical and Quantum Gravity, 2006, 23 (10): 3437.

DOI

14
高志勇, 王上, 王智. 基于FEM的引力参考传感器自引力计算与补偿[J]. 中国空间科学技术, 2024, 44 (2): 89- 97.

15
林志勇. 空间引力波探测自引力仿真技术研究[D]. 西安:长安大学,2023.

Options
文章导航

/

〈 〉