1

我正在尝试以数值方式求解两个相互倾斜的紧密椭圆轨道之间转移的 delta-V 成本。我使用的方法本质上是在一个节点处计算初始轨道的速度矢量，在相对节点处计算最终轨道，然后从初始飞行路径角度、初始半径和最终半径计算转移轨道。

一个关键步骤是计算转移轨道的比角动量矢量和偏心率矢量，从而计算转移轨道的周向惯性方向余弦矩阵。然而，当我从惯性系中位置和速度矢量的叉积计算惯性系中转移轨道的角动量矢量h时，我发现幅度之间存在显着误差（相对误差为-3.9521e-8）这个向量和前面在代码中计算的标量特定角动量。

这对我来说很奇怪，因为标量角动量用于计算速度矢量。我对精度损失发生在哪里感到困惑。

我尝试提供更精确的输入，特别是我一直在使用的 mu 值，但这根本没有改变相对误差。当我使用相同的叉积方法计算轨道1和2的特定角动量时，误差在机器精度的数量级上。


        mu = 3.98600437823e+14;

        thetaNT = -55.1582940061466; % deg
        eT = 0.022905923178296;
        aT = 7.243582592195826e+06; % m

        r1A = 7.146263097977215e+06; % m
        v1RA = -1.390985544431790e+02; % m/s
        v1ThetaA = 7.494958913236144e+03; % m/s

        eR1 = [0.355828643065080;-0.934551216774375;0];
        eTheta1 = [0.934551216774375;0.355828643065080;0];

        nCpf1 = [0.263190394679355,-0.840751409136755,0.473146789255815;
            0.880932410956014,0.00949753358184791,-0.473146789255815;
            0.393305102275257,0.541338032000730,0.743144825477394];
        nCpf2 = [0.107314578042381,-0.875080710676727,0.471929370924401;
            0.879361618777851,-0.137938482815824,-0.455736896003458;
            0.463903788257849,0.463903788257849,0.754709580222772];


        v1A = sqrt(v1RA^2 + v1ThetaA^2); % Total speed of orbit 1 at A

        hT = sqrt(aT*mu*(1-eT^2)); % Specific angular momentum of transfer orbit

        eRTB = [-cosd(thetaNT);sind(thetaNT+180);0];
        eThetaTB = [-sind(thetaNT+180);-cosd(thetaNT);0];

        % Calculation of radial speed and tangential speed
        vTRA = mu/hT*eT*sind(thetaNT);
        vTThetaA = mu/hT*(1+eT*cosd(thetaNT));

        vTA = sqrt(vTRA^2+vTThetaA^2);

        vTRB = mu/hT*eT*sind(thetaNT+180);
        vTThetaB = mu/hT*(1-eT*cosd(thetaNT));

        % Conversion of radius and speeds into radius and velocity vectors
        % in perifocal frames
        r1APF1 = r1A.*eR1;
        v1APF1 = v1RA.*eR1 + v1ThetaA.*eTheta1;

        vTBPFT = vTRB.*eRTB + vTThetaB.*eThetaTB;

        v2BPF2 = v2RB.*eR2 + v2ThetaB.*eTheta2;

        % Conversion to inertial reference frame
        r1AN = nCpf1*r1APF1;
        v1AN = nCpf1*v1APF1;

        v2BN = nCpf2*v2BPF2;

        rTAN = r1AN;
        vTAN = v1AN.*(vTA/v1A);

        % Calculation of angular momentum and eccentricity vectors in
        % inertial frame
        hTN = cross(rTAN, vTAN);
        eTN = cross(vTAN, hTN)./mu - rTAN./norm(rTAN);
        diffh = (norm(hTN)-hT)/hT
        diffe = (norm(eTN)-eT)/eT

我希望 diffh 和 diffe 大约是机器精度，大约为 2.2e-16，但它们要大得多。具体来说，diffh = -3.9689e-08，diffe = 7.5474e-05。

更新：该错误似乎出现在我计算径向和速度矢量的某处，如果这有助于集中搜索。

4

2 回答 2

0

我从惯性系中位置和速度矢量的叉积计算惯性系中转移轨道的角动量矢量h

这是最有可能出现错误的地方。我不遵循您使用的数学，因为我以不同的方式处理这些问题。我会找到你的转移轨道的轨道元素，然后你的卫星在你的转移轨道中的位置作为真正的异常，然后计算 delta V 应该是相当微不足道的，因为你知道你在每个时间步的速度向量在燃烧前后。如果您需要我对此进行扩展，请告诉我在哪里可以澄清，我会的。

于 2019-07-25T20:52:14.120 回答
0

机器精度是相对的，它由给定指数处的浮点数之间的距离组成。例如，在 Matlab 中，您可以检查给定数字的机器精度：

eps(1)
  ans = 
        2.22044604049250313e-16
eps(12345)
  ans = 
        1.818989403545856e-12.

通常，您应该期望大约e-16是您正在使用的价值的倍数。1e-8位置（ 的阶1e6）和速度（ 的阶）乘积的机器精度的数量级也是如此1e2。即产品应该是这样的，1e8因此它的机器精度是8-16=-8的数量级。

另请注意，如果您关心准确性，您可能会考虑使用弧度而不是度数。

于 2019-07-01T20:32:35.870 回答