基于四面体质元划分的纯引力轨道万有引力干扰计算方法

摘要

本发明公开了一种涉及航天动力学和科学计算技术领域的基于四面体质元划分的纯引力轨道万有引力干扰计算方法，包括以下步骤：首先，选择四面体作为质元形状对航天器进行划分，并定义如下参数：尺度参数SR、长宽比参数AR、力几何因子FGP、梯度几何因子GGP，然后计算不同SR、AR下的FGP、GGP，得到四面体质元尺度、形状对计算精度的影响曲线，再根据所要求的计算精度及影响曲线确定四面体质元尺度，并根据该质元尺度对航天器模型采用四面体质元进行划分，最后对万有引力进行计算。本发明的方法能够按照要求的精度计算航天器对验证质量的万有引力。在给定精度的情况下，可以得到对应的四面体质元尺度要求，并依此要求划分质元，完成万有引力作用的计算，满足所要求的精度。

说明书

技术领域

本发明属于航天动力学和科学计算技术领域，涉及一种纯引力轨道万有引力干扰的数值计算方法。

背景技术

在一些基础物理空间任务中，需要验证质量沿着纯引力轨道飞行提供测量基准以探测引力波和检验广义相对论，包括LISA和ASTROD任务(参见期刊《经典与量子引力》(Classicaland Quantum Gravity)2003年第20卷的文章“LISA的集成模型(The LISA integrated model)”和期刊《原子核物理B》(Nuclear Physics B)2007年第166卷153-158页文章“ASTROD(激光天文动力学)and ASTROD I”)。同样，利用验证质量沿着近地纯引力轨道飞行，并获取验证质量的纯引力轨道，能够用于精确测量地球重力场(参见2010年在台湾举行的第六届卫星星座和编队飞行国际研讨会会刊中的文章“采用精密编队飞行技术获取纯引力轨道”)。外航天器作用在验证质量上的万有引力是一个主要的干扰力，影响纯引力轨道的性能水平(参见《经典与量子引力》(Classical and Quantum Gravity)2004年第21卷第5期S653-S660页的文章“当前的LISA残余加速度误差估计(Current error estimates for LISA spuriousaccelerations)”)。为此，需要对纯引力轨道的万有引力干扰进行分析，为提高纯引力轨道性能打下基础。

一般而言，万有引力干扰无法直接测量，也不能通过解析方法进行整体计算。为此，LISA模型团队建立了数值计算方法。该方法采用航天器有限单元模型提供的结点质量和位置，并将每个单元近似为质点计算其对验证质量的引力、力矩和梯度作用，然后对所有单元求和得到整体量(参见《经典与量子引力》(Classical and Quantum Gravity)2005年第22卷第10期S395-S402页的文章“LISA自引力分析模型(Self-gravity modeling for LISA)”)。由于对单元采用了质点近似，必然会引起误差。并且不同的单元划分方式对于计算结果的影响也不同。

发明内容

本发明的目的是在纯引力轨道万有引力干扰数值计算中，能够实现满足精度要求的准确和有效计算。

本发明基于四面体质元划分的纯引力轨道万有引力干扰计算方法包括：选择四面体作为质元形状对航天器进行划分、分析四面体质元对计算精度的影响、提出满足精度要求的质元划分准则、外航天器的万有引力计算。所述的选择四面体作为质元形状，是基于四面体能够实现自动对复杂结构进行单元划分的优势的。分析四面体质元对计算精度的影响，是在给定的直角坐标系下，推导了正四面体对于外部一点的万有引力及其梯度表达式；然后定义了表示四面体尺度与其距离外部一点的相对大小的尺度参数SR，定义了表征四面体形状变化的长宽比参数AR；将四面体近似为质量集中于质心的质点，计算其对外部一点的万有引力及万有引力梯度；万有引力计算质点近似与精确表达式之间的相对误差，定义为力几何因子FGP；万有引力梯度计算质点近似与精确表达式之间的相对误差，定义为梯度几何因子GGP；分别计算在不同SR、AR值下的FGP和GGP，得到了质元尺度、形状对计算精度的影响曲线。所述的提出满足精度要求的质元划分准则，是在得到质元尺度、形状对计算精度的影响曲线后，根据所要求的万有引力干扰计算精度，参照曲线得到相应的质元划分要求。所述的外航天器的万有引力计算，是根据质元划分要求对航天器进行质元划分，将各质元近似为质点计算其对验证质量的万有引力及梯度，然后求和得到航天器对验证质量的万有引力及其梯度。

本发明的基于四面体质元划分的纯引力轨道万有引力干扰计算方法，步骤如下：

第一步：选择四面体作为质元形状对航天器进行划分；

第二步：定义如下参数：

(1)表示四面体尺度与其距离外部一点的相对大小的尺度参数SR；

(2)表示四面体形状变化的长宽比参数AR；

(3)表示万有引力计算质点近似与精确表达式之间的相对误差的力几何因子FGP；

(4)表示万有引力梯度计算质点近似与精确表达式之间的相对误差的梯度几何因子GGP；

第三步：计算不同尺度参数SR、长宽比参数AR下的力几何因子FGP、梯度几何因子GGP，得到四面体质元尺度、形状对计算精度的影响曲线；

第四步：根据所要求的万有引力干扰计算精度及所述的影响曲线确定相应的四面体质元尺度；

第五步：根据第四步所确定的质元尺度对航天器模型采用四面体质元进行划分；

第六步：将每一个四面体质元近似为点质量计算其对验证质量的万有引力及梯度；

第七步：将所有四面体质元的计算结果求和，得到航天器对验证质量的万有引力干扰。

优选的，所述表示四面体尺度与其距离外部一点的相对大小的尺度参数SR定义如下：

$>\mathrm{SR}=\frac{r_{P/G}}{L_1}$>

式中，r_P/G是点P到四面体质心的距离，L₁为四面体最长边的长度。

优选的，所述表示四面体形状变化的长宽比参数AR定义如下：

$>\mathrm{AR}=\frac{L_1}{L_2}$>

式中，L₁为四面体最长边的长度，L₂为四面体最短边的长度。

优选的，表示万有引力计算质点近似与精确表达式之间的相对误差的力几何因子FGP定义如下：

>

式中，F^P/C为四面体C对位于其外部一点P的万有引力，为将四面体C近似为位于质心的质点时对点P的万有引力。

优选的，所述的根据如下公式计算：

>

式中，G是万有引力常数、m^P为质点P的质量、ρ^C是四面体C的密度，V是四面体C的体积，点P到四面体质心的距离。

优选的，所述的F^P/C通过如下方法计算：

将四面体置于一个坐标系中，然后根据如下公式计算：

$>\left|F^{P/C}\right|=\sqrt{{\left(F_x^{P/C}\right)}^2+{\left(F_y^{P/C}\right)}^2+{\left(F_z^{P/C}\right)}^2}$>

四面体C对点P的万有引力在坐标轴x、y、z方向的向分量。

优选的，所述第三步中，计算不同尺度参数SR、长宽比参数AR下的力几何因子FGP通过如下方法实现：

将正四面体按照如下方法置于一个坐标系中：正四面体的一个顶点位于原点，该顶点所对的底面平行于x-y平面且位于上方，并且底面内的一条边平行于x轴；将所述正四面体做如下变形获得一系列不同长宽比参数AR的四面体：将正四面体拉伸，且与x-y平面平行的底面边长保持不变作为所述L₂，而其它三个边被拉长同等长度作为所述L₁，拉伸不同的长度即可得不同的长宽比参数AR，计算不同长宽比参数AR、尺度参数SR下的力几何因子FGP，所述的计算公式如下：

$>\frac{F_x^{P/C}}{{\mathrm{Gm}}^P}={\rho }^C{\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\xi -x}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\rho }^C{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\xi -x}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{F_y^{P/C}}{{\mathrm{Gm}}^P}={\rho }^C{\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\eta -y}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\rho }^C{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\eta -y}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{F_z^{P/C}}{{\mathrm{Gm}}^P}={\rho }^C{\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\zeta -z}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\rho }^C{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\zeta -z}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }$>

其中G是万有引力常数、m^P为质点P的质量、ρ^C是四面体C的密度，h为四面体高度，参数a为四面体C中与x-y平面平行的底面边长的一半，参数b为四面体C在x-y平面的投影三角形中平行于x轴的底边与y轴的交点到原点的距离，参数c为四面体C在x-y平面的投影三角形中平行于x轴的底边对应的顶点到原点的距离，(x、y、z)为质点P的坐标，(ξ、η、ζ)为四面体C内部一点的直角坐标。

优选的，表示万有引力梯度计算质点近似与精确表达式之间的相对误差的梯度几何因子GGP定义如下：

$>\mathrm{GGP}=\frac{\left|{\lambda }_{\mathrm{point},p}\right|-\left|{\lambda }_p\right|}{\left|{\lambda }_p\right|}$>

式中，λ_p为四面体C对位于其外部一点P的引力梯度张量矩阵T^P/C的主特征值，λ_point，p为将四面体C近似为位于质心的质点时对点P的引力梯度张量矩阵的主特征值。

优选的，所述第三步中，计算不同尺度参数SR、长宽比参数AR下的梯度几何因子GGP通过如下方法实现：

将正四面体按照如下方法置于一个坐标系中：正四面体的一个顶点位于原点，该顶点所对的底面平行于x-y平面且位于上方，并且底面内的一条边平行于x轴；将所述正四面体做如下变形获得一系列不同长宽比参数AR的四面体：将正四面体拉伸，且与x-y平面平行的底面边长保持不变作为所述L₂，而其它三个边被拉长同等长度作为所述L₁，拉伸不同的长度即可得不同的长宽比参数AR，计算不同长宽比参数AR、尺度参数SR下的梯度几何因子GGP，所述的引力梯度张量矩阵T^P/C通过如下公式计算：

$>T^{P/C}=\left(\begin{array}{ccc}\frac{{\partial F}_x^{P/C}}{\partial x}& \frac{{\partial F}_x^{P/C}}{\partial y}& \frac{{\partial F}_x^{P/C}}{\partial z}\\ \frac{{\partial F}_y^{P/C}}{\partial x}& \frac{{\partial F}_y^{P/C}}{\partial y}& \frac{{\partial F}_y^{P/C}}{\partial z}\\ \frac{{\partial F}_z^{P/C}}{\partial x}& \frac{{\partial F}_z^{P/C}}{\partial y}& \frac{{\partial F}_z^{P/C}}{\partial z}\end{array}\right)$>

式中，四面体C对点P的万有引力在坐标轴x、y、z方向的向分量；所述引力梯度张量矩阵通过如下公式计算：

>

式中，

>>

>>

>>

其中G是万有引力常数、m^P为质点P的质量、ρ^C是四面体C的密度，V是四面体C的体积，点P到四面体质心的距离，(x、y、z)为质点P的坐标，z_Q为四面体C在的质心在z轴上的坐标。

优选的，所述引力梯度张量矩阵T^P/C的计算公式中的各参量通过分别如下公式计算：

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_x^{P/C}}{\partial x}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{2{(\xi -x)}^2-{(\eta -y)}^2-{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{2{(\xi -x)}^2-{(\eta -y)}^2-{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_y^{P/C}}{\partial y}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-{(\xi -x)}^2+2{(\eta -y)}^2-{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-{(\xi -x)}^2+2{(\eta -y)}^2-{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_z^{P/C}}{\partial z}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-{(\xi -x)}^2-{(\eta -y)}^2+2{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-{(\xi -x)}^2-{(\eta -y)}^2+2{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_x^{P/C}}{\partial y}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\xi -x)(\eta -y)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\xi -x)(\eta -y)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_x^{P/C}}{\partial z}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\xi -x)(\zeta -z)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\xi -x)(\zeta -z)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_y^{P/C}}{\partial z}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\eta -y)(\zeta -z)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\eta -y)(\zeta -z)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

其中G是万有引力常数、m^P为质点P的质量、ρ^C是四面体C的密度，h为四面体高度，参数a为四面体C中与x-y平面平行的底面边长的一半，参数b为四面体C在x-y平面的投影三角形中平行于x轴的底边与y轴的交点到原点的距离，参数c为四面体C在x-y平面的投影三角形中平行于x轴的底边对应的顶点到原点的距离，(x、y、z)为质点P的坐标，(ξ、η、ζ)为四面体C内部一点的直角坐标。

本发明的基于四面体质元划分的纯引力轨道万有引力干扰计算方法的优点在于：

本发明的方法能够按照所要求的精度计算航天器对验证质量的万有引力作用。在给定精度的情况下，可以得到对应的四面体质元尺度要求，并依此要求划分质元，完成万有引力作用的计算，满足所要求的精度。

附图说明

图1为本发明的方法流程图；

图2为正四面体对外部一点的引力示意；

正四面体的一个顶点位于原点，该顶点所对的底面平行于x-y平面且位于上方，并且底面内的一条边平行于x轴；

图3为图2中的正四面体在x-y平面的影图；

图4为在不同的参数AR下，FGP绝对值的最大值随着参数SR的变化情况；

图5为在不同的参数AR下，GGP绝对值的最大值随着参数SR的变化情况。

具体实施方式

下面结合实施例和附图对本发明的基于四面体质元划分的纯引力轨道万有引力干扰计算方法做进一步说明。

本发明的方法在对航天器模型结构进行质元划分时，采用四面体作为质元形状，并定义了如下参数：

(1)定义力几何因子(force geometry parameter，FGP)以表征万有引力计算质点近似与精确表达式之间的相对误差

>

式中，F^P/C为四面体C对位于其外部一点P的万有引力，为将四面体C近似为位于质心的质点时对点P的万有引力，FGP的值越小，意味着对四面体质元进行质点近似所引起的误差越小。

(2)采用引力梯度张量的主特征值的绝对值定义梯度几何因子(gradient geometryparameter，GGP)为

>

其中，λ_p、λ_point，p分别是梯度张量T^P/C和的主特征值。

(3)定义四面体尺度与其距离外部一点的相对大小的尺度参数SR：

$>\mathrm{SR}=\frac{r_{P/G}}{L_1}---\left(4\right)$>

式中，r_P/G是点P到四面体质心的距离，L₁为四面体最长边的长度。

(4)定义表示四面体形状变化的长宽比参数AR：

$>\mathrm{AR}=\frac{L_1}{L_2}---\left(5\right)$>

式中，L₁为四面体最长边的长度，L₂为四面体最短边的长度。

参见附图1，通过计算不同尺度参数SR、长宽比参数AR下的力几何因子FGP、梯度几何因子GGP，能够得到四面体质元尺度、形状对计算精度的影响曲线；在确定了所要求的万有引力干扰计算精度后，可以对照影响曲线确定相应的四面体质元尺度；再根据该尺度对航天器进行划分，进而计算出航天器对验证质量的万有引力干扰。

为了便于计算，本发明的方法以正四面体为基础，并将正四面体进行一定的变形得到不同AR值的四面体，以便对不同SR、AR下的力几何因子FGP、梯度几何因子GGP进行计算。

计算时，如图2、图3所示，首先将正四面体按照如下方法置于一个坐标系中：正四面体的一个顶点位于原点，该顶点所对的底面平行于x-y平面且位于上方，并且底面内的一条边平行于x轴；将所述正四面体做如下变形获得一系列不同长宽比参数AR的四面体：将正四面体拉伸，且与x-y平面平行的底面边长保持不变作为所述L₂，而其它三个边被拉长同等长度作为所述L₁，拉伸不同的长度即可得不同长宽比参数AR的四面体，此时正四面体的质心在z轴上，为点Q＝(0，0，z_Q)。

1、计算力几何因子FGP，根据公式(1)，首先需要分别计算F^P/C和具体如下：

(1)计算四面体C对位于其外部一点P的万有引力F^P/C

首先，可以将四面体C置于一个直角坐标系中，并分别计算四面体C对位于其外部一点P的万有引力F^P/C在坐标轴x、y、z方向的分量然后采用如下公式计算万有引力F^P/C：

$>\left|F^{P/C}\right|=\sqrt{{\left(F_x^{P/C}\right)}^2+{\left(F_y^{P/C}\right)}^2+{\left(F_z^{P/C}\right)}^2}---\left(6\right)$>

对于均质四面体，其对位于其外部一点P的万有引力F^P/C在坐标轴x、y、z方向的分量可以分别通过如下公式计算：

$>\frac{F_x^{P/C}}{{\mathrm{Gm}}^P}={\rho }^C{\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\xi -x}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\rho }^C{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\xi -x}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }---\left(7\right)$>

$>\frac{F_y^{P/C}}{{\mathrm{Gm}}^P}={\rho }^C{\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\eta -y}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }$>

$>+{\rho }^C{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\eta -y}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{F_z^{P/C}}{{\mathrm{Gm}}^P}={\rho }^C{\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\zeta -z}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }---\left(9\right)$>

$>+{\rho }^C{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{\zeta -z}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{3/2}}\mathrm{d\eta d\xi d\zeta }$>

其中G是万有引力常数、m^P为质点P的质量、ρ^C是四面体C的密度，h为四面体高度，参数a为四面体C中与x-y平面平行的底面边长的一半，参数b为四面体C在x-y平面的投影三角形中平行于x轴的底边与y轴的交点到原点的距离，参数c为四面体C在x-y平面的投影三角形中平行于x轴的底边对应的顶点到原点的距离，(x、y、z)为质点P的坐标，(ξ、η、ζ)为四面体C内部一点的直角坐标。

(2)计算

将正四面体近似为位于质心的质点，则其对点P的引力变为

>

式中，G是万有引力常数、m^P为质点P的质量、ρ^C是四面体C的密度，V是四面体C的体积，点P到四面体质心的距离。

对于正四面体，其边长已知，根据参数a的定义可知其为边长的一半，因此可将其他参数用参数a进行换算，换算公式如下：

$>b=-\frac{\sqrt{3}}{3}a,$>$>c=\frac{2\sqrt{3}}{3}a,$>$>h=\frac{2\sqrt{6}}{3}a$>$>z_Q=\frac{\sqrt{6}}{2}a,$>$>V=\frac{2\sqrt{2}}{3}{a}^3---\left(11\right)$>

对于如前所述拉伸后的四面体，公式(10)中各参数和参数a之间的公式会发生变化，具体如下：

数b和c保持不变，参数h变化为

$>h=2a\sqrt{{\mathrm{AR}}^2-\frac{1}{3}}---\left(12\right)$>

以及

$>z_Q=\frac{3}{4}h=\frac{3}{2}a\sqrt{{\mathrm{AR}}^2-\frac{1}{3}},$>$>V=\frac{1}{3}h\left(\sqrt{3}{a}^2\right)=\frac{2\sqrt{3}}{3}{a}^3\sqrt{{\mathrm{AR}}^2-\frac{1}{3}}---\left(13\right)$>

2、计算梯度几何因子GGP，根据公式(2)，λ_p、λ_point，p分别是梯度张量T^P/C和的主特征值，因此计算GGP首先需要计算梯度张量T^P/C和

(1)计算精确表达的引力梯度张量T^P/C，公式如下：

$>T^{P/C}=\left(\begin{array}{ccc}\frac{{\partial F}_x^{P/C}}{\partial x}& \frac{{\partial F}_x^{P/C}}{\partial y}& \frac{{\partial F}_x^{P/C}}{\partial z}\\ \frac{{\partial F}_y^{P/C}}{\partial x}& \frac{{\partial F}_y^{P/C}}{\partial y}& \frac{{\partial F}_y^{P/C}}{\partial z}\\ \frac{{\partial F}_z^{P/C}}{\partial x}& \frac{{\partial F}_z^{P/C}}{\partial y}& \frac{{\partial F}_z^{P/C}}{\partial z}\end{array}\right)---\left(14\right)$>

此矩阵具有对称性且迹为零，因此仅有5个独立分量，5个独立分量及全部对角量的计算公式如下：

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_x^{P/C}}{\partial x}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{2{(\xi -x)}^2-{(\eta -y)}^2-{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }---\left(15\right)$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{2{(\xi -x)}^2-{(\eta -y)}^2-{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_y^{P/C}}{\partial y}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-{(\xi -x)}^2+2{(\eta -y)}^2-{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }---\left(16\right)$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-{(\xi -x)}^2+2{(\eta -y)}^2-{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_z^{P/C}}{\partial z}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-{(\xi -x)}^2-{(\eta -y)}^2+2{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }---\left(17\right)$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-{(\xi -x)}^2-{(\eta -y)}^2+2{(\zeta -z)}^2}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_x^{P/C}}{\partial y}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\xi -x)(\eta -y)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }---\left(18\right)$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\xi -x)(\eta -y)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_x^{P/C}}{\partial z}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\xi -x)(\zeta -z)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }---\left(19\right)$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\xi -x)(\zeta -z)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

$>\frac{1}{{\mathrm{Gm}}^P{\rho }^C}\frac{\partial F_y^{P/C}}{\partial z}={\int }_0^h{\int }_0^{\frac{\zeta }{h}a}{\int }_{\frac{\zeta }{h}b}^{\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\eta -y)(\zeta -z)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }---\left(20\right)$>

$>+{\int }_0^h{\int }_{-\frac{\zeta }{h}a}^0{\int }_{\frac{\zeta }{h}b}^{-\frac{b-c}{a}\xi +\frac{\zeta }{h}c}\frac{-3(\eta -y)(\zeta -z)}{{[{(\xi -x)}^2+{(\eta -y)}^2+{(\zeta -z)}^2]}^{5/2}}\mathrm{d\eta d\xi d\zeta }$>

其中G是万有引力常数、m^P为质点P的质量、ρ^C是四面体C的密度，h为四面体高度，参数a为四面体C中与x-y平面平行的底面边长的一半，参数b为四面体C在x-y平面的投影三角形中平行于x轴的底边与y轴的交点到原点的距离，参数c为四面体C在x-y平面的投影三角形中平行于x轴的底边对应的顶点到原点的距离，(x、y、z)为质点P的坐标(ξ、η、ζ)为四面体C内部一点的直角坐标。

(2)计算将正四面体近似为位于质心的质点时的则引力梯度张量公式如下：

>

其中

>>

>>

>>

同样，在进行和T^P/C的计算时，对于正四面体，可以将各参数按照公式(11)进行换算，对于拉伸后的四面体可以公式(12)、公式(13)进行换算。

根据以上各表达式，可以计算出均质四面体质元在体外任意一点的FGP和GGP值。在计算中取正四面体的边长为单位长度，P点的点质量设为单位质量。

图4是在不同的参数AR值下，FGP绝对值的最大值随着参数SR增加的变化情况，表现了四面体质元尺度和形状对FGP值的综合影响。随着质元尺度的减小，FGP绝对值的最大值减小的更快；随着参数AR的增加，FGP绝对值的最大值逐渐增加并在参数AR＞10之后趋于稳定。这表明四面体质元的尺度要比形状对FGP的影响严重。并且，如果参数SR不小于4，也就是限制四面体质元尺度不大于其质心到验证质量距离的1/4，则对于各种形状的四面体质元FGP的绝对值总是小于10^-2。

图5是在不同AR值下，GGP绝对值的最大值随SR的变化情况，表现了四面体质元尺度和形状的综合影响。同样，四面体质元的尺度要比形状对GGP的影响严重，且若SR不小于5，即限制质元尺度不大于其质心到验证质量距离的1/5，则对于各种形状的四面体质元GGP的绝对值总小于10^-2。

在限制四面体质元尺度不大于其质心到验证质量距离的1/5的情况下，对于各种形状的四面体质元计算误差总是低于1％。例如，对于具有半径25cm腔体的内编队系统，要求距离验证质量最近的部位划分质元的尺度要小于5cm；或者，对于距离验证质量1m的外卫星部件，如果其长度不超过20cm，可以直接将其近似为质点来计算对验证质量的万有引力，而不需要再进一步划分质元。

本发明的方法计算过程中采用正面体均匀拉伸进行计算，其结论可以适用于四面体的一般变形，因为正四面体的计算结果说明只要质心距达到边长的一定倍数之后，计算精度能够满足；且拉伸后的四面体出现了不同长度的边，计算结果表明只要质心距达到长边的一定倍数之后，计算精度能够满足。因此对于一般的四面体，只要质心距达到最长的边的一定倍数之后，计算精度就能够满足。