Причиной несоответствия между методами расчета удельного момента импульса? - PullRequest
Новая
       41

Причиной несоответствия между методами расчета удельного момента импульса?

0 голосов
/ 01 июля 2019

Я пытаюсь численно определить для дельта-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 будет отличаться и будет порядкаТочность машины, примерно 2,2е-16, но они намного больше.В частности, diffh = -3.9689e-08 и dif = 7.5474e-05.

ОБНОВЛЕНИЕ: Кажется, ошибка появляется где-то в моих вычислениях радиальных и скоростных векторов, если это помогает сконцентрировать ваш поиск.

1 Ответ

0 голосов
/ 01 июля 2019

Точность станка равна относительно , она состоит из расстояния между числами с плавающей запятой при данном показателе степени. Например, в Matlab вы можете проверить точность станка для заданного числа: 

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

Обычно вы должны ожидать чего-то примерно в e-16 раз превышающее значение, с которым вы работаете. Таким образом, 1e-8 находится в порядке величины точности станка для произведения положения (порядка 1e6) и скорости (порядка 1e2). То есть продукт должен быть порядка 1e8, следовательно, его точность машины порядка 8-16 = -8.

Также обратите внимание, что если вас беспокоит точность, вы можете рассмотреть возможность использования углов в радианах, а не в градусах.

...