Вывод уравнений движения ЛА типа "самолет"
Допущения:
1) Конструкция ЛА абсолютно жесткая (ЛА - твердое тело);
2) Атмосфера спокойная (W=0 - скорость ветра);
3) Масса ЛА постоянная т.к. заданное время переходных процессов меньше 60 с;
4) Земля плоская и не вращается;
5) Угол тангажа не равен +/-90 градусов;
6) Спарведлива гипотеза стационарности (согласно этой гипотезе, аэродинамическая сила и аэродиномический момент в нестационарном режиме полета полностью определяются величинами параметров в данный момент времени).
Кинематические уравнения поступательного движения

R={L[t],H[t],dZ[t]}(*Радиус вектор центра масс ЛА, L[t] - дальность полета, H[t] - высота полета, dZ[t] - боковое отклонение*)

vg={vgx[t],vdy[t],vgz[t]} (*Вектор земной сокрости центра масс ЛА в нормальной системе координат*)

MatrixForm[D[R,t]]==MatrixForm[vg] (*Кинематическое уравнение движения центра масс самолета в векторной форме*)

MatrixForm[M1={{Cos[Bteta[t]],-Sin[Bteta[t]],0},{Sin[Bteta[t]],Cos[Bteta[t]],0},{0,0,1}}] (*Матрица поворота на угол наклона траектории*)

MatrixForm[M2={{Cos[BPsi[t]],0,Sin[BPsi[t]]},{0,1,0},{-Sin[BPsi[t]],0,Cos[BPsi[t]]}}] (*Матрица поворота на угол пути*)

MatrixForm[M=M2.M1] (*Матрица перехода от траекторной системы координат к нормальной системе координат *)

MatrixForm[ReplaceAll[M,{BPsi[t]->0,Bteta[t]->0}]] (*Проверка подстановкой нулевых углов поворотов. Должна получиться единичная матрица*)

MatrixForm[vzem={vk,0,0}] (*Вектор земной скорости в траекторной системе координат*)

MatrixForm[vg=M.vzem] (*Проецируем вектор земной скорости на нормальные оси координат*)

MatrixForm[D[R,t]]==MatrixForm[vg] (*Кинематическое уравнение движения центра масс самолета в векторной форме*)

u1=D[R[[1]],t]==vg[[1]]

u2=D[R[[2]],t]==vg[[2]]

u3=D[R[[3]],t]==vg[[3]]

MatrixForm[MM1={{Cos[LPsi[t]],0,Sin[LPsi[t]]},{0,1,0},{-Sin[LPsi[t]],0,Cos[LPsi[t]]}}] (*Матрица поворота на угол рыскания*)

MatrixForm[MM2={{Cos[LTeta[t]],-Sin[LTeta[t]],0},{Sin[LTeta[t]],Cos[LTeta[t]],0},{0,0,1}}](*Матрица поворота на угол тангажа*)

MatrixForm[MM3={{1,0,0},{0,Cos[LGamma[t]],-Sin[LGamma[t]]},{0,Sin[LGamma[t]],Cos[LGamma[t]]}}] (*Матрица поворота на угол крена*)

MatrixForm[MM=MM1.MM2.MM3] (*Матрица перехода от связанной системы координат к нормальной системе координат*)

MatrixForm[SM1={{Cos[Alpha[t]],-Sin[Alpha[t]],0},{Sin[Alpha[t]],Cos[Alpha[t]],0},{0,0,1}}] (*Матрица поворота на угол атаки*)

MatrixForm[SM2={{Cos[Betta[t]],0,Sin[Betta[t]]},{0,1,0},{-Sin[Betta[t]],0,Cos[Betta[t]]}}] (*Матрица поворота на угол скольжения*)

MatrixForm[SM=Transpose[SM2.SM1]] (*Матрица перехода от скоростной системы координат к связанной системе координат*)

vsk={vvozd[t],0,0} (*Воздушая скорость ЛА в скоростной системе координат*)

vsv=SM.vsk  (*Проецируем воздушную скорость ЛА на оси связанной системы координат*)

MatrixForm[D[R,t]]==MatrixForm[MM.vsv] (*Кинематическое уравнение движения центра масс самолета в векторной форме*)

uu1=D[R[[1]],t]==(MM.vsv)[[1]]

uu2=D[R[[2]],t]==(MM.vsv)[[2]]

uu3=D[R[[3]],t]==(MM.vsv)[[3]]

Динамические уравнения поступательного движения

v={vx[t],vy[t],vz[t]} (*Вектор воздушной скорости ЛА в связанной системе координат*)

f={fx,fy,fz} (*Вектор результирующей сил, действующих на ЛА в полете, в связанной системе координат*)

w={wx[t],wy[t],wz[t]} (*Вектор угловой скорости ЛА в связанной системе координат*)

dupd=massa*(D[v,t]+Cross[w,v])==f (*Динамическое уравнение движения центра масс самолета в векторной форме*)

dupdfk=Solve[dupd,{D[vx[t],t],D[vy[t],t],D[vz[t],t]}] (*Приводим систему дифференциальных уравнений к стандартному виду (форме Коши),разрешая её относительно производных*)

u4=dupdfk[[1,1]]

u5=dupdfk[[1,2]]

u6=dupdfk[[1,3]]

dupd2=massa*(D[vsv,t]+Cross[w,vsv])==f (*Подставим в уравнение проекции воздушной скорости на оси связанной системы координат*)

uu4=dupd2[[1,1]]==dupd2[[2,1]]

uu5=dupd2[[1,2]]==dupd2[[2,2]]

uu6=dupd2[[1,3]]==dupd2[[2,3]]

uuu=Simplify[Solve[{uu4,uu5,uu6},{D[vvozd[t],t],D[Alpha[t],t],D[Betta[t],t]}]] (*Приводим систему дифференциальных уравнений к стандартному виду (форме Коши),разрешая её относительно производных*)

uuu4=uuu[[1,1]]

uuu5=uuu[[1,2]]

uuu6=uuu[[1,3]]

Кинематические уравнения вращательного движения

proekc={D[LGamma[t],t],0,0}+Transpose[MM2.MM3].{0,D[LPsi[t],t],0}+Transpose[MM3].{0,0,D[LTeta[t],t]} (*Проецируем угловые скорости рыскания, тангажа и крена на оси связанной системы координат*)

u7=w[[1]]==proekc[[1]]

u8=w[[2]]==proekc[[2]]

u9=w[[3]]==proekc[[3]]

ufk=Simplify[Solve[{u7,u8,u9},{D[LGamma[t],t],D[LPsi[t],t],D[LTeta[t],t]}]] (*Приводим систему дифференциальных уравнений к стандартному виду (форме Коши),разрешая её относительно производных*)

uu7=ufk[[1,1]]

uu8=ufk[[1,2]]

uu9=ufk[[1,3]]

Динамические уравнения вращательного движения

MatrixForm[TI={{Ixx,-Ixy,-Ixz},{-Iyx,Iyy,-Iyz},{-Izx,-Izy,Izz}}] (*Тензор инерции*)

Ixz=Izx=Iyz=Izy=0 (*Ось OZ перпендикулярна плоскости симметрии XOY*)

Iyx=Ixy

MatrixForm[TI]

MatrixForm[K=TI.w] (*Главный момент количеств движения в связанной системе координат*)

M={Mx,My,Mz} (*Результирующий момент аэродинамических сил и тяги в связанной системе координат(Момент силы тяжести, центробежных и кориолисовых сил инерции вокруг центра масс ЛА равен нулю)*)

lchduvd=D[K,t]+Cross[w,K] (*Полная производная по времени главного вектора момента количеств движения*)

Полная производная главного вектора момента количеств движения равна главному моменту внежних сил, действующих на тело.

u10=lchduvd[[1]]==M[[1]]

u11=lchduvd[[2]]==M[[2]]

u12=lchduvd[[3]]==M[[3]]

udvdfk=Simplify[Solve[{u10,u11,u12},{D[wx[t],t],D[wy[t],t],D[wz[t],t]}]] (*Приводим систему дифференциальных уравнений к стандартному виду (форме Коши),разрешая её относительно производных*)

uu10=udvdfk[[1,2]]

uu11=udvdfk[[1,3]]

uu12=udvdfk[[1,1]]