Моделирование полёта тела с учетом скорости ветра

[Leandoer]

Добрый день! Решаю задачу о моделировании полёта тела с учетом сопротивления воздуха и скорости ветра.
Составил систему дифференциальных уравнений.

2-й Закон Ньютона: $m\overrightarrow{a} = m\overrightarrow{g} - \frac{1}{2} c_{f} \rho S (V-V_{w})^2 \overrightarrow{e}$
$e_{x} = \frac{\left\lvert{V_{x}-V_{w_{x}}}\right\rvert}{\left\lvert V-V_{w}\right\rvert}, e_{y} = \frac{\left\lvert{V_{y}-V_{w_{y}}}\right\rvert}{\left\lvert V-V_{w}\right\rvert}, e_{z} = \frac{\left\lvert{V_{z}-V_{w_{z}}}\right\rvert}{\left\lvert V-V_{w}\right\rvert}$

Где $(V-V_{w})^2 = (V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + (V_{z}-V_{w_{z}})^2$
$V $ - скорость тела, $V_{w}$ - скорость ветра

Проекции скоростей:
$V_{w_{x}} = V_{w}\cos(\alpha), V_{w_{y}} = V_{w}\sin(\alpha), V_{w_{z}} = 0 $
$V_{x}=V\cos(\varphi)\sin(\theta), V_{y}=V\sin{\varphi}\sin(\theta), V_{z}=V\cos(\theta)$

Таким образом,
$ \frac{dV_{x}}{dt} = -\frac{c_{f} \rho S}{m} \left| V_{x}-V_{w_{x}} \right| \sqrt{(V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + V_{z}^2} \\ \frac{dV_{y}}{dt} = -\frac{c_{f} \rho S}{m} \left| V_{y}-V_{w_{y}} \right| \sqrt{(V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + V_{z}^2} \\ \frac{dV_{z}}{dt} = -g -\frac{c_{f} \rho S}{m} \left| V_{z}-V_{w_{z}} \right| \sqrt{(V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + V_{z}^2}$

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

[svv]

Ошибки:
1) Потерян коэффициент $\frac 1 2$ в итоговых формулах.
2) В процитированном фрагменте выделенные величины надо писать как векторы:

Leandoer в сообщении #1608210 писал(а):

2-й Закон Ньютона: $m\overrightarrow{a} = m\overrightarrow{g} - \frac{1}{2} c_{f} \rho S ({\color{magenta}V}-{\color{magenta}V_w})^2 \overrightarrow{e}$
$e_{x} = \frac{\left\lvert{V_{x}-V_{w_{x}}}\right\rvert}{\left\lvert {\color{magenta}V}-{\color{magenta}V_w}\right\rvert}, e_{y} = \frac{\left\lvert{V_{y}-V_{w_{y}}}\right\rvert}{\left\lvert {\color{magenta}V}-{\color{magenta}V_w}\right\rvert}, e_{z} = \frac{\left\lvert{V_{z}-V_{w_{z}}}\right\rvert}{\left\lvert {\color{magenta}V}-{\color{magenta}V_w}\right\rvert}$

Где $({\color{magenta}V}-{\color{magenta}V_w})^2 = (V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + (V_{z}-V_{w_{z}})^2$
${\color{magenta}V}$ - скорость тела, ${\color{magenta}V_w}$ - скорость ветра

Математически осмысленны оба выражения, $|\overrightarrow{V}-\overrightarrow{V}_w|$ и $|V-V_w|$, где $V=|\overrightarrow{V}|, V_w=|\overrightarrow{V}_w|$. Но они не совпадают, и в Ваши уравнения должно входить первое.

Мелкие придирки и рекомендации:
Вместо $V, V_w$ я бы писал $\mathbf v, \mathbf w$.
Уравнение движения:
$\dfrac{d\mathbf v}{dt}=\mathbf g-k|\mathbf v-\mathbf w|(\mathbf v-\mathbf w)$ ,
где $k=\dfrac{c_{f} \rho S}{2m}$ (константа, если все величины в дроби постоянны).
Если скорость ветра $\mathbf w$ постоянна в пространстве и времени, уравнение движения упростится, если записать его через относительную скорость $\mathbf v_r=\mathbf v-\mathbf w$ (что эквивалентно переходу в систему отсчёта, в которой скорость ветра равна нулю):
$\dfrac{d\mathbf v_r}{dt}=\mathbf g-k|\mathbf v_r|\mathbf v_r$
Вектор $\mathbf e$ не нужен совсем.

Leandoer в сообщении #1608210 писал(а):

Проекции скоростей:
$V_{w_{x}} = V_{w}\cos(\alpha), V_{w_{y}} = V_{w}\sin(\alpha), V_{w_{z}} = 0 $
$V_{x}=V\cos(\varphi)\sin(\theta), V_{y}=V\sin{\varphi}\sin(\theta), V_{z}=V\cos(\theta)$

Не уверен, что будет удобно работать с таким представлением скоростей.

[svv]

Ошибки 3)

Leandoer в сообщении #1608210 писал(а):

$e_{x} = \frac{{\color{magenta}|}{V_{x}-V_{w_{x}}}{\color{magenta}|}}{\left\lvert V-V_{w}\right\rvert}, e_{y} = \frac{{\color{magenta}|}{V_{y}-V_{w_{y}}}{\color{magenta}|}}{\left\lvert V-V_{w}\right\rvert}, e_{z} = \frac{{\color{magenta}|}{V_{z}-V_{w_{z}}}{\color{magenta}|}}{\left\lvert V-V_{w}\right\rvert}$
...
$\frac{dV_{x}}{dt} = -\frac{c_{f} \rho S}{m} {\color{magenta}|} V_{x}-V_{w_{x}} {\color{magenta}|} \sqrt{(V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + V_{z}^2}\\\frac{dV_{y}}{dt} = -\frac{c_{f} \rho S}{m} {\color{magenta}|} V_{y}-V_{w_{y}} {\color{magenta}|} \sqrt{(V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + V_{z}^2}\\\frac{dV_{z}}{dt} = -g -\frac{c_{f} \rho S}{m} {\color{magenta}|} V_{z}-V_{w_{z}} {\color{magenta}|} \sqrt{(V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + V_{z}^2}$

Модули, выделенные цветом, надо убрать (заменив, где надо, на скобки). Они не просто излишни, а дают ошибочные выражения. В Вашем варианте все компоненты ускорения всегда будут отрицательными (точнее, неположительными), хотя это вовсе не обязательно так.

[Leandoer]

svv в сообщении #1608225 писал(а):

Ошибки 3)

Leandoer в сообщении #1608210 писал(а):

$e_{x} = \frac{{\color{magenta}|}{V_{x}-V_{w_{x}}}{\color{magenta}|}}{\left\lvert V-V_{w}\right\rvert}, e_{y} = \frac{{\color{magenta}|}{V_{y}-V_{w_{y}}}{\color{magenta}|}}{\left\lvert V-V_{w}\right\rvert}, e_{z} = \frac{{\color{magenta}|}{V_{z}-V_{w_{z}}}{\color{magenta}|}}{\left\lvert V-V_{w}\right\rvert}$
...
$\frac{dV_{x}}{dt} = -\frac{c_{f} \rho S}{m} {\color{magenta}|} V_{x}-V_{w_{x}} {\color{magenta}|} \sqrt{(V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + V_{z}^2}\\\frac{dV_{y}}{dt} = -\frac{c_{f} \rho S}{m} {\color{magenta}|} V_{y}-V_{w_{y}} {\color{magenta}|} \sqrt{(V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + V_{z}^2}\\\frac{dV_{z}}{dt} = -g -\frac{c_{f} \rho S}{m} {\color{magenta}|} V_{z}-V_{w_{z}} {\color{magenta}|} \sqrt{(V_{x}-V_{w_{x}})^2 + (V_{y}-V_{w_{y}})^2 + V_{z}^2}$

Модули, выделенные цветом, надо убрать (заменив, где надо, на скобки). Они не просто излишни, а дают ошибочные выражения. В Вашем варианте все компоненты ускорения всегда будут отрицательными (точнее, неположительными), хотя это вовсе не обязательно так.

Большое спасибо! В конечной системе исправил модули. Однако, численно исследуя случай, когда скорость ветра направлена строго по оси X против скорости тела, я получаю смещение по оси Y, которого быть не должно. Очень странно...

[svv]

Leandoer в сообщении #1608318 писал(а):

В конечной системе исправил модули.

А коэффициент $\frac 1 2$ тоже добавили?

Leandoer в сообщении #1608318 писал(а):

Однако, численно исследуя случай, когда скорость ветра направлена строго по оси X против скорости тела, я получаю смещение по оси Y, которого быть не должно. Очень странно...

Ну, если Вы покажете программу, шансы найти ошибку увеличатся.

[sergey zhukov]

Leandoer в сообщении #1608318 писал(а):

смещение по оси Y, которого быть не должно.

Видимо, эффект Магнуса...

[Leandoer]

svv в сообщении #1608339 писал(а):

Leandoer в сообщении #1608318 писал(а):

В конечной системе исправил модули.

А коэффициент $\frac 1 2$ тоже добавили?

Leandoer в сообщении #1608318 писал(а):

Однако, численно исследуя случай, когда скорость ветра направлена строго по оси X против скорости тела, я получаю смещение по оси Y, которого быть не должно. Очень странно...

Ну, если Вы покажете программу, шансы найти ошибку увеличатся.

Да, 1/2 тоже добавлена.

Перед тем как скинуть программу уточню, что отклонение по оси Y очень-очень маленькое. Оно имеет порядок $10^{-14}$. Может, это какая-то нестрашная ошибочка, которую можно проигнорировать и дальше продолжать исследование с текущей реализацией. Что скажете?

[svv]

Пока ничего, недостаточно информации.

[Leandoer]

svv в сообщении #1608356 писал(а):

Пока ничего, недостаточно информации.

Вот основной фрагмент кода, без вывода графиков:

Используется синтаксис C++

#include <cmath>
#include <string>
#include <map>

using namespace std;

class Frem{
public:
    Frem ( map<string, double> data ) : t(0.0), m(data.at("m")), r(data.at("r")), x(0.0), y(0.0), z(data.at("h")),
        vx(data.at("v")*cos(data.at("phi")*(M_PI/180))*sin(data.at("theta")*(M_PI/180))), vy(data.at("v") * sin(data.at("phi")*(M_PI/180))*sin(data.at("theta")*(M_PI/180))), vz(data.at("v") * cos(data.at("theta")*(M_PI/180))),
        kappa(data.at("kappa")), g(data.at("g")),
        wx(data.at("w") * cos(data.at("alpha")*(M_PI/180))), wy(data.at("w")*sin(data.at("alpha")*(M_PI/180))), wz(0.0){};

    double get_t() const {return t;}
    double get_x() const {return x;};
    double get_y() const {return y;};
    double get_z() const {return z;};

    void make_step(double dt);


private:
    double g;
    double t;
    double m, r;
    double x, y, z;
    double vx, vy, vz;
    double wx, wy ,wz;
    double kappa;

    double ax(double vx, double vy, double vz, double wx, double wy, double wz) const;


    double ay(double vx, double vy, double vz, double wx, double wy, double wz) const;


    double az(double vx, double vy, double vz, double wx, double wy, double wz) const;

    void Frem::make_step(double dt)
{
    double kx1=vx*dt;
    double ky1=vy*dt;
    double kz1=vz*dt;
    double kvx1=ax(vx, vy, vz, wx, wy, wz)*dt;
    double kvy1=ay(vx, vy, vz, wx, wy, wz)*dt;
    double kvz1=az(vx, vy, vz, wx, wy, wz)*dt;

    double kx2=(vx+kvx1/2)*dt;
    double ky2=(vy+kvy1/2)*dt;
    double kz2=(vz+kvz1/2)*dt;
    double kvx2=ax(vx+kvx1/2, vy+kvy1/2, vz+kvz1/2, wx, wy,wz)*dt;
    double kvy2=ay(vx+kvx1/2, vy+kvy1/2, vz+kvz1/2, wx, wy,wz)*dt;
    double kvz2=az(vx+kvx1/2, vy+kvy1/2, vz+kvz1/2, wx ,wy,wz)*dt;

    double kx3=(vx+kvx2/2)*dt;
    double ky3=(vy+kvy2/2)*dt;
    double kz3=(vz+kvz2/2)*dt;
    double kvx3=ax(vx+kvx2/2, vy+kvy2/2,vz+kvz2/2 ,wx ,wy,wz)*dt;
    double kvy3=ay(vx+kvx2/2 ,vy+kvy2/2,vz+kvz2/2 ,wx ,wy,wz)*dt;
    double kvz3=az(vx+kvx2/2 ,vy+kvy2/2,vz+kvz2/2 ,wx ,wy,wz)*dt;

    double kx4=(vx+kvx3)*dt;
    double ky4=(vy+kvy3)*dt;
    double kz4=(vz+kvz3)*dt;
    double kvx4=ax(vx+kvx3 ,vy+kvy3,vz+kvz3 ,wx ,wy,wz)*dt;
    double kvy4=ay(vx+kvx3 ,vy+kvy3,vz+kvz3 ,wx ,wy,wz)*dt;
    double kvz4=az(vx+kvx3 ,vy+kvy3,vz+kvz3 ,wx ,wy,wz)*dt;

    t+=dt;
    x+=(kx1 + 2*kx2 + 2*kx3 + kx4)/6;
    y+=(ky1 + 2*ky2 + 2*ky3 + ky4)/6;
    z+=(kz1 + 2*kz2 + 2*kz3 + kz4)/6;
    vx+=(kvx1 + 2*kvx2 + 2*kvx3 + kvx4)/6;
    vy+=(kvy1 + 2*kvy2 + 2*kvy3 + kvy4)/6;
    vz+=(kvz1 + 2*kvz2 + 2*kvz3 + kvz4)/6;

}


double Frem::ax(double vx, double vy, double vz, double wx, double wy, double wz) const
{
    return -kappa * 1.29 * M_PI * r * r /(2*m) * (vx - wx) * sqrt(pow(vx - wx, 2) + pow(vy - wy, 2) + pow(vz, 2));
}

double Frem::ay(double vx, double vy, double vz, double wx, double wy, double wz) const
{
    return  - kappa * 1.29 * M_PI * r * r /(2*m) * (vy - wy) * sqrt(pow(vx - wx, 2) + pow(vy - wy, 2) + pow(vz, 2));
}

double Frem::az(double vx, double vy, double vz, double wx, double wy,double wz) const
{
    return -g -kappa * 1.29 * M_PI * r * r /(2*m) * (vz - wz) * sqrt(pow(vx - wx, 2) + pow(vy - wy, 2) + pow(vz, 2));
}

};
 

[wrest]

Leandoer в сообщении #1608355 писал(а):

что отклонение по оси Y очень-очень маленькое. Оно имеет порядок $10^{-14}$.

Увеличьте точность расчетов. double это около 14-16 значащих десятичных цифр и ошибка как раз в последней? Возьмите точность расчетов например 100 значащих цифр и посмотрите какое будет отклонение.

[svv]

Leandoer
Я посмотрел Вашу программу. Она написана понятно. Вы решаете систему ДУ первого порядка методом Рунге-Кутта. Метод реализован правильно. Я очень старался найти хоть какую-то ошибку, но не смог. Остаётся вопрос — какие Вы брали параметры и начальные условия (m,r,h,v,theta,phi,kappa,g,w,alpha). А также конечное значение времени. Хотелось бы взять точно те значения, при которых у Вас получилось отклонение по $y$ порядка $10^{-14}$, чтобы воспроизвести результат.

Leandoer в сообщении #1608318 писал(а):

скорость ветра направлена строго по оси X против скорости тела

Это значит, что один из углов $\varphi, \alpha$ равен $0$, а второй $180$ (в градусах). И для того угла, который $180$, с низкой точностью вычисляется выражение sin(180*(M_PI/180)), входящее в формулу для $v_y(0)$ или $w_y$. Математически это нуль, но у меня в программе, например, получается 1.2246063538e-16.

[Leandoer]

Например, $m=1, R=1, V=90, W=40, \alpha=180, \varphi=0, \theta=45, \kappa=0.0005, x_{0}=y_{0}=z_{0}=0, g=10, timestep=0.001$

[svv]

Да, механизм такой, как я описал. Проделайте такие опыты:
1) Сейчас экземпляр класса Frem инициализируется с помощью конструктора
Frem(map<string, double>)
в котором есть список инициализации, но тело конструктора пустое. Напишите в теле { wy=0; }
В результате смещение по $y$ исчезнет.

2) Теперь напишите в теле { wy=w*sin(M_PI); }
Смещение вернётся к прежнему значению.

3) Просто вычислите sin(M_PI) и посмотрите на точность.

Используя тип long double, можно уменьшить эффект на три порядка, но средство, с помощью которого Вы строите 3d графики, вытянет и эту погрешность. Если есть возможность установить равный масштаб по всем осям, сделайте это.

[Leandoer]

svv в сообщении #1608489 писал(а):

Да, механизм такой, как я описал. Проделайте такие опыты:
1) Сейчас экземпляр класса Frem инициализируется с помощью конструктора
Frem(map<string, double>)
в котором есть список инициализации, но тело конструктора пустое. Напишите в теле { wy=0; }
В результате смещение по $y$ исчезнет.

2) Теперь напишите в теле { wy=w*sin(M_PI); }
Смещение вернётся к прежнему значению.

3) Просто вычислите sin(M_PI) и посмотрите на точность.

Используя тип long double, можно уменьшить эффект на три порядка, но средство, с помощью которого Вы строите 3d графики, вытянет и эту погрешность. Если есть возможность установить равный масштаб по всем осям, сделайте это.

Большое спасибо!