Здравствуйте.
Подскажите, в чем может быть проблема? Всю голову уже сломал, пытаясь найти ошибку.
Код:
In[1]:= Param0 = {n -> 16, N -> 11, m1t -> 70, m2t -> 80, O1C1t -> 0.1, at -> 0.15,
bt -> 0.20, ct -> 0.15, M0t -> 8, k1t -> 1.3, k2t -> 0.8}
Out[1]= {n -> 16, N -> 11, m1t -> 70, m2t -> 80, O1C1t -> 0.1, at -> 0.15, bt -> 0.2,
ct -> 0.15, M0t -> 8, k1t -> 1.3, k2t -> 0.8}
In[2]:= Param = {k1 -> k1t + 0.01 N, k2 -> k2t (1 + 0.01 N) 10^-2,
M0 -> M0t (1 + 0.01 n) 10^2, m1 -> m1t (1 + 0.01 N),
m2 -> m2t (1 + 0.01 N), O1C1 -> O1C1t (1 + 0.001 N) 10^-3,
a -> at + 0.001 n, b -> bt + 0.001 N, c -> ct + 0.001 n, R1 -> 0.12,
R2 -> 0.18, \[Alpha] -> 0.0024, \[Beta] -> 0.54, \[Omega]z0 ->
590, \[Tau] -> 1.20, \[CapitalDelta]\[Tau] -> 0.2, zA -> 0,
zB -> (bt + 0.001 N) + (ct + 0.001 n)} /. Param0
Out[2]= {k1 -> 1.41, k2 -> 0.00888, M0 -> 928., m1 -> 77.7, m2 -> 88.8,
O1C1 -> 0.0001011, a -> 0.166, b -> 0.211, c -> 0.166, R1 -> 0.12,
R2 -> 0.18, \[Alpha] -> 0.0024, \[Beta] -> 0.54, \[Omega]z0 -> 590, \[Tau] ->
1.2, \[CapitalDelta]\[Tau] -> 0.2, zA -> 0, zB -> 0.377}
In[3]:= Param1 = {MDz -> M0 - k1*\[Phi]'[t], MCz -> -k2*\[CurlyPhi]''[t]^2} /. Param
Tfin = \[Tau] /. Param
Out[3]= {MDz -> 928. - 1.41 Derivative[1][\[Phi]][t],
MCz -> -0.00888 (\[CurlyPhi]^\[Prime]\[Prime])[t]^2}
Out[4]= 1.2
Нахождение массы ротора, координат центров масс дисков и всего ротора относительно системы координат Axyz
In[5]:= M = m1 + m2 /. Param
rС1 = {{O1C1*Sin[\[Beta]]}, {-O1C1*Cos[\[Beta]]}, {b}} /. Param
rС2 = {{0}, {0}, {-a}} /. Param
xC = (m1 rС1[[1, 1]] + m2 rС2[[1, 1]])/M /. Param
yC = (m1 rС1[[2, 1]] + m2 rС2[[2, 1]])/M /. Param
Out[5]= 166.5
Out[6]= {{0.0000519791}, {-0.0000867143}, {0.211}}
Out[7]= {{0}, {0}, {-0.166}}
Out[8]= 0.0000242569
Out[9]= -0.0000404667
Нахождение тензоров инерции дисков
In[10]:= I1C1x1y1z1 = m1 R1^2/4 {{1, 0, 0}, {0, 1, 0}, {0, 0, 2}} /. Param
I2C2\[Xi]\[Eta]\[Zeta] = m2 R2^2/4 {{1, 0, 0}, {0, 1, 0}, {0, 0, 2}} /. Param
Out[10]= {{0.27972, 0, 0}, {0, 0.27972, 0}, {0, 0, 0.55944}}
Out[11]= {{0.71928, 0, 0}, {0, 0.71928, 0}, {0, 0, 1.43856}}
Нахождение матрицы перехода S=\[Gamma]^T от системы координат С2\[Xi]\[Eta]\[Zeta] к
системе координат C2x2y2z2
In[12]:= S = {{1, 0, 0}, {0, Cos[\[Alpha]], -Sin[\[Alpha]]}, {0, Sin[\[Alpha]],
Cos[\[Alpha]]}} /. Param
Out[12]= {{1, 0, 0}, {0, 0.999997, -0.0024}, {0, 0.0024, 0.999997}}
Нахождение тензора инерции диска 2 в системе координат C2x2y2z2 в
виде I2 C2 x2 y2 z2
=STIC2 \[Xi]\[Eta]\[Zeta] S
In[13]:= I2C2x2y2z2 = Transpose[S].I2C2\[Xi]\[Eta]\[Zeta].S
Out[13]= {{0.71928, 0., 0.}, {0., 0.719284, 0.00172627}, {0., 0.00172627, 1.43856}}
Нахождение тензоров инерции дисков в системе координат Axyz
In[14]:= I1Axyz = I1C1x1y1z1 +
m1 ((Transpose[rС1].rС1)[[1, 1]] IdentityMatrix[3] -
rС1.Transpose[rС1]) /. Param
Out[14]= {{3.739, 3.5022*10^-7, -0.000852183}, {3.5022*10^-7, 3.739,
0.00142166}, {-0.000852183, 0.00142166, 0.559441}}
In[15]:= I2Axyz = I2C2x2y2z2 +
m2 ((Transpose[rС2].rС2)[[1, 1]] IdentityMatrix[3] -
rС2.Transpose[rС2]) /. Param
Out[15]= {{3.16625, 0., 0.}, {0., 3.16626, 0.00172627}, {0., 0.00172627, 1.43856}}
Нахождение тензора инерции ротора в системе координат Axyz
In[16]:= IAxyz = I1Axyz + I2Axyz
Out[16]= {{6.90526, 3.5022*10^-7, -0.000852183}, {3.5022*10^-7, 6.90526,
0.00314792}, {-0.000852183, 0.00314792, 1.998}}
Составление уравнений
In[17]:= eqMz = IAxyz[[3, 3]] *\[Phi]''[t] == MDz + MCz /. Param1 /. Param
eqQx = -M yC *\[Phi]''[t] - M xC *\[Phi]'[t]^2 == XA + XB
eqQy = M xC *\[Phi]''[t] - M yC *\[Phi]'[t]^2 == YA + YB
eqMx = IAxyz[[1, 3]]* \[Phi]''[t] - IAxyz[[2, 3]]* \[Phi]'[t]^2 == -zA *YA -
zB *YB /. Param
eqMy = IAxyz[[2, 3]] *\[Phi]''[t] + IAxyz[[1, 3]]*\[Phi]'[t]^2 ==
zA *XA + zB *XB /. Param
Out[17]= 1.998 (\[Phi]^\[Prime]\[Prime])[t] ==
928. - 1.41 Derivative[1][\[Phi]][t] -
0.00888 (\[CurlyPhi]^\[Prime]\[Prime])[t]^2
Out[18]= -0.00403878 Derivative[1][\[Phi]][t]^2 +
0.0067377 (\[Phi]^\[Prime]\[Prime])[t] == XA + XB
Out[19]= 0.0067377 Derivative[1][\[Phi]][t]^2 +
0.00403878 (\[Phi]^\[Prime]\[Prime])[t] == YA + YB
Out[20]= -0.00314792 Derivative[1][\[Phi]][t]^2 -
0.000852183 (\[Phi]^\[Prime]\[Prime])[t] == -0.377 YB
Out[21]= -0.000852183 Derivative[1][\[Phi]][t]^2 +
0.00314792 (\[Phi]^\[Prime]\[Prime])[t] == 0.377 XB
Решение уравнений.
Выражение реакций
In[22]:= solR = Solve[{eqQx, eqQy, eqMx, eqMy}, {XA, XB, YA, YB}]
Out[22]= {{XA -> -0.00177835 Derivative[1][\[Phi]][t]^2 -
0.00161222 (\[Phi]^\[Prime]\[Prime])[t],
XB -> -0.00226043 Derivative[1][\[Phi]][t]^2 +
0.00834992 (\[Phi]^\[Prime]\[Prime])[t],
YA -> -0.00161222 Derivative[1][\[Phi]][t]^2 +
0.00177835 (\[Phi]^\[Prime]\[Prime])[t],
YB -> 0.00834992 Derivative[1][\[Phi]][t]^2 +
0.00226043 (\[Phi]^\[Prime]\[Prime])[t]}}
In[23]:= XA[t_] := XA /. solR
XB[t_] := XB /. solR
YA[t_] := YA /. solR
YB[t_] := YB /. solR
RA[t_] := Sqrt[XA[t]^2 + YA[t]^2]
RB[t_] := Sqrt[XB[t]^2 + YB[t]^2]
Численное решение дифференциального уравнения вращения ротора
In[33]:= sol = NDSolve[{eqMz, \[Phi][0] == 0, \[Phi]'[0] == \[Omega]z0} /.
Param, {\[Phi], \[Phi]', \[Phi]''}, {t, 0, Tfin}]
During evaluation of In[33]:= NDSolve::underdet: There are more dependent variables, {\[Phi][t],\[CurlyPhi][t]}, than equations, so the system is underdetermined.
Out[33]= NDSolve[{1.998 (\[Phi]^\[Prime]\[Prime])[t] ==
928. - 1.41 Derivative[1][\[Phi]][t] -
0.00888 (\[CurlyPhi]^\[Prime]\[Prime])[t]^2, \[Phi][0] == 0,
Derivative[1][\[Phi]][0] == 590}, {\[Phi], Derivative[
1][\[Phi]], \[Phi]^\[Prime]\[Prime]}, {t, 0, 1.2}]
Ошибка
Цитата:
NDSolve::underdet: There are more dependent variables, {\[Phi][t],\[CurlyPhi][t]}, than equations, so the system is underdetermined.
Очень надеюсь на вашу помощь