Динамика ротора с многошариковым автобалансировочным устройством Бакалаврская работа

Динамика ротора с многошариковым автобалансировочным устройством Бакалаврская работа

Санкт-Петербургский государственный университет
Механика и математическое моделирование
Теоретическая механика

Динамика ротора с многошариковым автобалансировочным устройством

Бакалаврская работа

Научный руководитель:

Ассистент, к. ф.м.-н.

Рецензент:

Доцент, к. ф.м.-н.

Санкт-Петербург

2017

SAINT-PETERSBURG STATE UNIVERSITY
Mechanics and Mathematical modeling
Theoretical Mechanics

Maksimov Anton Sergeevich

Dynamics of rotor with multiball autobalancer

Bachelor’s Thesis

Scientific supervisor:

Assistant, PhD A. S. Kovachev

Reviewer:

Аssociate professor, PhD V. kov

Saint-Petersburg

2017

Оглавление

Введение. Механическая модель ротора с многополостным многошариковым АБУ. Вывод уравнений движения ротора с АБУ. Модель соударения шариков АБУ. Алгоритм численного интегрирования с учетом соударений шариков в АБУ Примеры расчетов движения ротора с многополостным АБУ

1). Автобалансировка ротора с трехполостным АБУ с 9 шариками при отсутствии соударений

2) Автобалансировка однополостного АБУ с 2 шариками при наличии соударения

3) Автобалансировка двуполостного АБУ с 8 шариками при наличии соударения

Вывод. Список литературы. Приложение 1. Код программы.

Введение

Современные технические средства позволят довести конструкционную и технологическую составляющие дисбаланса ротора до допустимого минимума. Однако изменение величины дисбаланса может происходить во время эксплуатации роторной машины, и для того, чтобы заняться балансировкой придется остановить механизм, зачастую это невозможно или нежелательно. В таких случаях используются автоматические балансировочные устройства (АБУ). Они находят применение во многих механизмах с быстровращающимися элементами, например CD, DVD приводы, HDD накопители, электродвигатели, различные центрифуги. В промышленной, транспортной, бытовой и прецизионной технике все шире применяются шаровые автобалансирующие устройства, предназначенные для полной балансировки высокооборотных роторов с переменным дисбалансом.

Исследованиям данных машин посвящены многие работы, как отечественных [1-3], так и зарубежных ученых[4-6]. В большинстве работ авторы пренебрегают возможными столкновениями балансировочных шариков АБУ. И если в исследовании стационарных режимов [8] это справедливо, то при исследовании влияния динамического дисбаланса на амплитуду прецессионных колебаний [6] возникает вопрос, сохранятся ли полученные результаты, если учесть взаимодействие шаров АБУ.

Целью данной работы является разработка алгоритма построения движения простейшей модели ротора (модель Джеффкотта [4]) с многополостным, многошариковым автобалансировочным устройством принимая во внимание взаимодействие шаров АБУ между собой. Сделать вывод о влиянии столкновения шариков на прецессионные колебания ротора.

Механическая и математическая модель ротора с АБУ

Рис. 1

Рассматривается модель статически неуравновешенного ротора, представляющего собой жесткий диск радиуса R со смещенным центром тяжести (точка G) , симметрично закрепленный посередине невесомого изотропного упругого вала. Вал вращается вертикально в шарнирных опорах O1 ,O2 . Диск оснащен многополостным, многошариковым автобалансировочным устройством (АБУ), представляющим обойму, содержащую несколько кольцеобразных полостей, заполненных вязкой жидкостью, в которых могут свободно передвигаться балансировочные шарики. Геометрический центр АБУ совпадает с геометрическим центром диска (точка С). Полагаем, что движения диска происходить в горизонтальной плоскости, перпендикулярной оси вала.

Введем неподвижную систему координат OXYZ следующим образом: начало координат O находится в середине отрезка O1O2; ось OZ направленна вертикально вверх и проходит через точки O1 ,O2 ; оси OX и OY взаимоперпендикулярны и лежат в горизонтальной плоскости. (Рис. 1)

Рис. 2

Обойма АБУ содержит N не пересекающихся полостей (Рис. 2), равномерно распределенных между минимальным (0 < Rmin) и максимальным (Rmax < R) радиусом, i-ая полость имеет радиус ri и в ней находятся ni шариков плотности с, массы m и радиуса rb. Для балансировки ротора с эксцентриситетом s = ||CG|| необходимо, чтобы балансировочный коэффициент у, введенный в работе [работа 2006 года], удовлетворял следующему условию:

Отсюда можно рассчитать геометрические параметры АБУ:

При этом должны выполняться следующий условия:

Первое условие означает что полости в обойме АБУ не накладываются друг на друга. Второе условие гарантирует наличие зазора между шариками в каждой полости.

Вывод уравнений движения ротора с АБУ.

Введем следующие обозначения:
X, Y положение центра ротора (точка C);
и – угол поворота вала, отсчитывается от оси OX против часовой стрелки;
шi, j – угол отклонения j-го шарика в i-ой полости АБУ относительно оси, проходящей через точки C и G;
M – масса диска;
m – масса одного шарика АБУ;
k – коэффициент упругости вала;
c, d, f – коэффициенты диссипации, учитывающие потерю энергии, для поперечного смещения диска, для поворота вала, для движения шарика в АБУ;
Mвр – внешний вращающий момент, приложенный к валу;
Jc – момент инерции диска;
– радиус инерции диска.

Для получения уравнения движения системы воспользуемся метод Лагранжа. Выражения для кинетической, потенциальной энергии и диссипативной функции Релея имеют вид:

В которых:

Запишем Уравнения Лагранжа второго рода, относительно обобщенных координат X, Y, и, шi, j :
(1)                

Введем безразмерные координаты и время: . Тогда система (1) перепишется в виде:
(2)                        

Где введены безразмерные параметры:

Модель соударения шариков АБУ

a)
б)
в)
Рис. 3

В качестве модели взаимодействия шариков АБУ между собой возьмем модель абсолютно упругого удара. Рассмотрим ее на примере одновременного соударения трех шаров одинаковой массы и радиуса, движущихся по кольцеобразной полости. Примем допущение, что радиусы шариков много меньше радиуса окружности, по которой они движутся. В момент времени, предшествующий удару, скорости шариков ш1, ш3 направлены как показано на Рис. 3,a, скорость шарика ш2 равна 0. Непосредственно в момент удара скорости всех шариков равны 0 (Рис.3 б). После удара, согласно закону сохранения количества движения и принятому допущению о размерах шаров, скорость шарика ш1 станет равна скорости шарика ш3 до удара, для и наоборот, а скорость шарика ш2 не изменится (Рис.3 в). Аналогично рассматриваем соударение для произвольного числа шаров.

Алгоритм построения движения ротора с учетом взаимодействия шариков АБУ

е – точность вычислений;
– начальные значения;

Для того, чтобы воспользоваться методом интегрирования Рунге-Кутта приведем систему (2) к системе дифференциальных уравнений первого порядка. Для этого сделаем замену:

С учетом введенных обозначений получаем систему N+6 дифференциальных уравнений.

Эта система линейна относительно производных и автономна, значит ее можно привести к виду:

Где z – вектор состояния:

Обозначим за K(z, h) функцию приращения вектора состояния за один шаг h по времени:        

Схема алгоритма построения движения ротора с АБУ

Вычисления m-го вектора состояния

Проверка точности вычислений.

В случае если, ||zm-1 – zm|| > е, то уменьшаем шаг в два раза и повторяем пункт I Если же было найдено 10 значений zm-10…zm без изменения шага h, то увеличиваем его в два раза.

Расстояние между шариками АБУ.

Введем функцию минимального расстояния между шарами АБУ на m-м шаге интегрирования:

Если выполнено условие , то произошел удар на промежутке между zm-1 и zm состоянием системы и надо перейти к пункту IV, иначе возвращаемся к пункту I.

Поиск точки столкновения

Если д(zm-1 , h) отрицательно, то ищем he такое, чтобы: |д(zm-1 , he)| < е. Как далее будет показано, всегда положительна.
Так как , и функция непрерывна по , то зафиксируем и введем функцию . Теперь можно воспользоваться численными методами для нахождения корня he уравнения:

(3)                                        

Например можно воспользоваться методом хорд, описанным в [].
При нахождения корня he по этому методу попутно получаем h1 , h2:

Обработка данных после нахождения точки столкновения

Имеем: h1 , he, h2 из предыдущего пункта, zm-1 – предыдущее состояние системы.

Вычисляем по методу Рунге-Кутта следующие три вектора состояния, соответственно: до удара, во время удара, после удара:

Здесь полностью соответствует вектору zm+1 , только у него изменены скорости шариков, участвующих в столкновении согласно введённой модели соударения шаров для движения после удара. После вычисления меняем скорости соответствующих шариков на ноль в векторе zm+1 .

На данном этапе возможны случаи, при которых , тогда координаты шаров, участвующих в соударении и соседних к ним, вернем к значениям, соответствующим состоянию zm так, чтобы расстояния между всеми шарами были положительны.

Теперь меняем шаг итерации h на (he - h1) , это значительно ускоряет вычисления, в случае когда шарики в АБУ начинают тереться друг о друга (когда скорости сталкивающихся шаров близки). После чего возвращаемся к пункту I

Примеры расчетов движения ротора с многополостным АБУ

1). Автобалансировка ротора с трехполостным АБУ с 9 шариками при отсутствии соударений.

Рассмотрим ротор с трехполостным АБУ при следующих параметрах:

Параметры системы: M=1; Jc=0.01; с = 7800; c=20; d=0.5; f=0.02; Mвр = 500; k=0.00002; Щ = 447; Rmax=0.1; Rmin=0.01; N=3; n1=2; n2=3; n3=4; у=2.

Все начальные скорости равны нулю, начальное положение как на Рис. 4

Рис.4 Начальное положение. Красная точка – центр масс шариков, зеленая – диска.
a).
б).
в).
г).
Рис. 5 a). График движения шариков в первой полости; б). График движения шариков во второй полости; в). График движения шариков в третьей полости; г). График амплитуды прецессионных колебаний.

Из графиков видно, что результаты расчетов, полученные по модели учитывающей взаимодействия шариков, при отсутствии соударений, полностью согласуется с результатами расчетов по модели не учитывающей взаимодействия.

2) Автобалансировка однополостного АБУ с 2 шариками при наличии соударения

Рассмотрим ротор с двухшариковым АБУ при следующих параметрах:

M=1; Jc=0.01; с = 7800; c=20; d=0.5; f=0.02; Mвр = 500; k=0.00002; Щ = 447; Rmax=0.1; Rmin=0.01; N=1; n1=2; у=2.

Все начальные скорости равны нулю, начальное положение как на Рис. 6

Рис.6 Начальное положение

а)
б)
Рис. 7 Графики движения шариков a). с учетом взаимодействия шаров; б). без учета взаимодействия
Рис.8 Амплитуда прецессионного движения ротора.

По графикам (Рис. 7, 8) видно, что взаимодействие шариков между собой влияет на траектории движения шариков в АБУ и на амплитуду прецессионного колебания ротора.

3) Автобалансировка двуполостного АБУ с 8 шариками при наличии соударения

Рассмотрим ротор с двуполостным АБУ при следующих параметрах:

M=1; Jc=0.01; с = 7800; c=20; d=0.5; f=0.02; Mвр = 500; k=0.00002; Щ = 447; Rmax=0.1; Rmin=0.01; N=2; n1=6; n2=2; у=2.

Все начальные скорости равны нулю, начальное положение как на Рис. 9

Рис. 9 Начальное положение

Рис.10 Движение шариков в первой полости.

Рис. 11 Амплитуда прецессионных колебаний.

Табл.1 Положение шариков в АБУ после стабилизации ротора
По модели при учете взаимодействия шариков	По модели без учета взаимодействия шариков

Как видно взаимодействие шариков между собой влияет не только на траектории движения шариков и амплитуду колебаний, но и на конечное положение шаров внутри АБУ после стабилизации ротора

Вывод

В ходе работы был построен алгоритм численного интегрирования системы, описывающей динамику ротора, оснащённого многополостным, многошариковым автобалансировочным устройством с учетом взаимодействия шариков между собой. Получены численные и графические результаты, анализ которых дает возможность утверждать, что взаимодействие шариков между собой внутри АБУ слабо влияет на прецессионные колебания ротора, и при этом не увеличивает время автобалансировки.

Список литературы

1). Bykov V. G. Kovachev A. S. Dynamics of a rotor with an eccentric ball auto-balancing device. // Vestnik St. Petersburg University: Mathematics 2014.Vol 47(4) P. 173–180

2). Dimentberg F. M. Shatalov K. T. Gusarov A. A. Oscillation of machines. Moscow: Mechanical engineering. 1964 year.[In Russian]

3). Dimentberg F. M. Bending vibrations of rotating shafts. Moscow. Acad. 1959.[In Russian]

4). Genta G. Dynamics of Rotating Systems. Springer Science 2005

5).J Chung, I Jang Dynamic response and stability analysis of an automatic ball balancer for a flexible rotor. // Journal of Sound and Vibration, 2003. Vol. 1(259) P. 31-43

6). L. Sperling, B. Ryzhik, H. Duckstein Single-Plane Auto-Balancing of Rigid Rotors// Techniche mechanic, band 24, Heft 1,(2004), 1-24

7) D. J. Rodrigues a, A. R. Champneys a,, M. I. Friswell b, R. E. Wilson, Experimental investigation of a single-plane automatic balancing mechanism for a rigid rotor// Journal of Sound and Vibration 330 (2011) 385–403

8). CТАЦИОНАРНЫЕ РЕЖИМЫ ДВИЖЕНИЯ НЕУРАВНОВЕШЕННОГО РОТОРА C АВТОБАЛАНСИРОВОЧНЫМ МЕХАНИЗМОМ // Вестник СПбГУ. Сер. 1, 2006, вып. 2

Приложение 1. Код программы

Rotor parameters

(*Disk parameters*)

M=1;  (*Disk weight*)

h=0.002;  (*Disk thicness*)

с0=7800;  (*Disk density*)

Rd=Sqrt[M/(Pi с0 h)];  (*Disk radius*)

Jc = M Rd^2/2;  (*Disk moment of inertia*)

radI=Sqrt[Jc/M]; (*Disk radius of inertia*)

(*MBA parameters*)

n2 = 1;  (*Amount of MBA radius*)

n = {8};  (*Amount of ball on each radius*)

N0=;

Rmax=0.1;  (*MBA maximum radius*)

Rmin=0.05;  (*MBA minimum radius*)

r=Table[Rmin+(Rmax-Rmin)i/n2,{i,1,n2}]; (*List of MBA radius*)

s=0.005;

у=2;

с0b=7800;  (*Ball density*)

m=у M s/; (*Ball weight*)

radball = (3 m/(4 Pi с0b))^(1/3);  (*Ball radius*)

For[j=1,j<=n2,j++,

  If[2Pi/ArcSin[radball/r[[j]]]<n[[j]],

  Print["---------------------------------- ALARM--------------------------------------------"];

  Print["With ball radius = ",radball," you can't place ",n[j]," balls on radius №",j," = ",r[[j]]];

  ]

  ];

If[2radball>(Rmax-Rmin)/n2,

  Print["---------------------------------- ALARM--------------------------------------------"];

  Print["With ball radius = ",radball," you can't place ",n2," bag in autobalancer"];

  ];

c=20;  (*внешнее сопротивление поперечному движению ротора*)

b=0.01;  (* внутреннее трение в вале*)

d=0.5;  (*внешнее сопротивление вращательному движению *)

f=0.02;  (*сопротивление движению шариков*)

Mv=500;  (*внешний вращающий момент, приложен к диску*)

k=2 10^5;  (*коэф упругости вала*)

Щ=Sqrt[k/M]; 

(*Безразмерные параметры*)

ϵ=s/radI; м=m/M; с=Table[r[[i]]/radI,{i,1,n2}]; д=c/(M Щ); в=b/(M Щ); г=d/(M radI^2 Щ); ч=f/(m radI^2 Щ); к=Mv/(M radI^2 Щ^2); сball=radball/radI;

Equations of the system

(*Уравнения ротора*)

urn1 := (1+м )x''[t]+(д+в)x'[t]+в и'[t]y[t]+x[t]-ϵ(и'[t]^2 Cos[и[t]]+и''[t]Sin[и[t]])-м ;

urn2 := (1+м )y''[t]+(д+в)y'[t]-в и'[t]x[t]+y[t]-ϵ(и'[t]^2 Sin[и[t]]-и''[t]Cos[и[t]])-м ;

urn3 := и''[t]+г и'[t]+в(x'[t]y[t]-x[t]y'[t]+и'[t](x[t]^2+y[t]^2))-к-ϵ(x''[t] Sin[и[t]]-y''[t]Cos[и'[t]])-м ч ;

(*Уравнения шаров АБУ*)

urnbi_,j_:=ч шi, j'[t]+с[[i]](шi, j''[t]+и''[t])-с[[i]](x''[t]Sin[и[t]+шi, j[t]]-y''[t]Cos[и[t]+шi, j[t]]);

(*List of system coordinates*)

variables={x[t],y[t],и[t]};

variables2={};

For[j=1,j<=n2,j++,

  variables2 = Join[variables2,Table[шj, i[t],{i, n[[j]]}]];

  ];

variables=Join[variables2,D[variables2,t],variables, D[variables, t]];

(*List of system equations*)

urn={urn1==0,urn2==0,urn3==0};

For[j=1,j<=n2,j++,

  urn = Join[urn, Table[urnbj, i==0,{i, n[[j]]}]];

  ];

A={};

G1={};

G2={};

G3={};

G4={};

Clear[i, j];

For[i=1,i<=n2,i++,

  G1=Join[G1,Table[Coefficient[urn1,шi, j''[t]],{j, n[[i]]}]];

  G2=Join[G2,Table[Coefficient[urn2,шi, j''[t]],{j, n[[i]]}]];

  G3=Join[G3,Table[Coefficient[urn3,шi, j''[t]],{j, n[[i]]}]];

  ];

G1=Join[G1,{Coefficient[urn1,x''[t]],Coefficient[urn1,y''[t]],Coefficient[urn1,и''[t]]}];

G2=Join[G2,{Coefficient[urn2,x''[t]],Coefficient[urn2,y''[t]],Coefficient[urn2,и''[t]]}];

G3=Join[G3,{Coefficient[urn3,x''[t]],Coefficient[urn3,y''[t]],Coefficient[urn3,и''[t]]}];

A=Join[{G1},{G2},{G3}];

For[i=1,i<=n2,i++,

  For[j=1,j<=n[[i]],j++,

  G4={};

  For[num=1,num<=n2,num++,

  G4=Join[G4,Table[Coefficient[urnbi, j,шnum, num2''[t]],{num2,n[[num]]}]];

  ];

  G4=Join[G4,{Coefficient[urnbi, j,x''[t]],Coefficient[urnbi, j,y''[t]],Coefficient[urnbi, j,и''[t]]}];

  A=Join[A,{G4}];

  ];

  ];

B={{urn1},{urn2},{urn3}};

For[i=1,i<=n2,i++, B=Join[B, Table[{urnbi, j},{j, n[[i]]}]]];

variab=Join[D[variables[[1;;Sum[n[[i]],{i, n2}]]],{t,2}],D[variables[[-3;;-1]],t]];

B=-Simplify[Expand[B-A. variab]];

nv=Length[variables];

function[cord_]:=Module[{urnl=A, urnr=B, sol, out, p,i, t1},

  t1=SessionTime[];

  For[i=1,i<=nv, i++,

  urnl=urnl/.variables[[i]]->cord[[i]];

  urnr=urnr/.variables[[i]]->cord[[i]];

  ];

  sol=Transpose[LinearSolve[urnl, urnr]][[1]];

  out=Join[Table[cord[[i+N0]],{i, N0}],Delete[sol,{{-1},{-2},{-3}}],{cord[[-3]],cord[[-2]],cord[[-1]]},{sol[[-3]],sol[[-2]],sol[[-1]]}];

  (*Print["out = ",out];*)

  tfunc=Join[tfunc,{SessionTime[]-t1}];

  out

  ];

Initial Conditions

primarydisk={0,0,0,0,0,0};

primaryball={};

(*Ball initial position*)

For[i=1,i<=n2,i++,

  For[j=1,j<=n[[i]],j++,

  primaryball=Join[primaryball,{2 Pi j/n[[i]] }];

  ]

  ];

(*Ball initial speed*)

primaryball=Join[primaryball, Table[0,{i, Sum[n[[j]],{j, n2}]}]];

primary=Transpose[{Join[primaryball, primarydisk]}];

coord=primary;

Solution Constructor

Time=0.5*Щ;

coord=primary;

time={0};

step=0.1;

null=10^-5;

tk={};

tfunc={};

tdist={};

tspeed={};

tall={};

thit2={};

k5[time_,step_,y0_]:=Module[{k1,k2,k3,k4,t1},

  t1=SessionTime[];

  k1=step function[y0];

  k2=step function[y0+k1/2];

  k3=step function[y0+k2/2];

  k4=step function[y0+k3];

  tk=Join[tk,{SessionTime[]-t1}];

  (k1+2 k2+2k3+k4)/6

  ];

distance[y_]:=Module[{dist, asin, i,j, p,nu, t1},

  t1=SessionTime[];

  dist={};

  For[i=1,i<=n2,i++,

  If[n[[i]]==1,

  dist=Join[dist,{{,100 null}}];

  ,

  asin=2ArcSin[сball/с[[i]]];

  nu=;

  dist=Join[dist, Table[{j+nu, y[[j+1+nu]]-y[[j+nu]]-asin},{j,1,n[[i]]-1}]];

  dist=Join[dist,{{nu+n[[i]],y[[1+nu]]+2 Pi-y[[nu+n[[i]]]]-asin}}];

  ];

  ];

  tdist=Join[tdist,{SessionTime[]-t1}];

  dist

  ];

spdist[y_]:=Module[{sp, i,j, p,nu, t1},

  t1=SessionTime[];

  sp={};

  For[i=1,i<=n2,i++,

  If[n[[i]]==1,

  sp=Join[sp,{{,0}}];

  ,

  nu=;

  sp=Join[sp, Table[{j+nu, Abs[y[[j+1+nu]]-y[[j+nu]]]},{j,1,n[[i]]-1}]];

  sp=Join[sp,{{nu+n[[i]],Abs[y[[1+nu]]-y[[nu+n[[i]]]]]}}];

  ];

  ];

  tspeed=Join[tspeed,{SessionTime[]-t1}];

  sp

  ];

hit1[y_]:=Module[{dist=Abs[y],i},

  For[i=1,i<=Length[dist],i++,

  If[dist[[i,2]]>null, dist[[i,2]]=0,dist[[i,2]]=1]

  ];

  dist

  ];

hit2[y_]:=Module[{hit, u=y, i,j, p,T, b,flag, k,t1},

  t1=SessionTime[];

  hit={{}};

  For[i=1,i<=n2,i++,

  If[n[[i]]!=1,

  p=y[[1+;;,All]];

  T=Table[If[ Mod[k, n[[i]]]==j-1,1,0],{k, n[[i]]},{j, n[[i]]}];

  While[p[[1,2]]!=0,

  p=T. p;

  ];

  b={};

  flag=False;

  For[j=1,j<=n[[i]],j++,

  If[p[[j,2]]==1,

  b=Join[b,{p[[j,1]]}];

  flag=True;

  ,

  If[flag,

  b=Join[b,{Mod[b[[-1]]+1-,n[[i]],1]+}];

  hit=Join[hit,{b}];

  b={};

  flag=False;

  ];

  ];

  ];

  If[flag,

  b=Join[b,{Mod[b[[-1]]+1-,n[[i]],1]+}];

  hit=Join[hit,{b}];

  b={};

  flag=False;

  ];

  ];

  ];

  hit=Delete[hit,1];

  thit2=Join[thit2,{SessionTime[]-t1}];

  hit

  ];

tstart=SessionTime[];

Print["Session Time = ",Dynamic[SessionTime[]-tstart, UpdateInterval->2]];

Print["Iteration time = ", Dynamic[time[[-1]],UpdateInterval->5]];

Print["pop = ", Dynamic[pop, UpdateInterval->5]];

Print["step = ", Dynamic[step, UpdateInterval->5]];

Print["posmin = ", Dynamic[posmin, UpdateInterval->5]];

Print["spmin = ", Dynamic[spmin, UpdateInterval->5]];

Print["h1 = ", Dynamic[h1,UpdateInterval->5]];

Print["h2 = ", Dynamic[h2,UpdateInterval->5]];

Print["he = ", Dynamic[he, UpdateInterval->5]];

Print["e = ", Dynamic[e, UpdateInterval->5]];

(*Print["min1 = ", Dynamic[Min[distance[yhit1[[1;;N0]]]],UpdateInterval5]];

Print["min3 = ", Dynamic[Min[distance[yhit3[[1;;N0]]]],UpdateInterval5]];*)

pop=0;

pop2=0;

nstep=1;

root={};

While[time[[-1]]<=Time,

  t2=SessionTime[];

  ynext=coord[[All,-1]]+k5[time[[-1]],step, coord[[All,-1]]];

  (*Print["ynext=  ",ynext];

  Print["step=  ",step];*)

  flag1=True;

  While[flag1,

  y1=coord[[All,-1]]+k5[time[[-1]],step/2,coord[[All,-1]]];

  y2 =y1+k5[time[[-1]],step/2,y1];

  If[Norm[ynext-y2](*+Norm[distance[ynext[[1;;N0]]][[All,2]]-distance[y2[[1;;N0]]][[All,2]]]+Norm[spdist[ynext[[1+N0;;2N0]]][[All,2]]-spdist[y2[[1+N0;;2N0]]][[All,2]]]*)<null,

  flag1=False;

  ,

  ynext=y1;

  step=step/2;

  nstep=0;

  ];

  ];

  (*Print["disnatce1  ",distance[coord[[All,-1]]+k5[time[[-1]],step, coord[[All,-1]]]]];*)

  dball=distance[ynext[[1;;N0]]];

  e=Min[dball[[All,2]]];

  If[e<0,

  h1=0;

  h2=step;

  yhit2=ynext;

  While[(Abs[e]>=null),

  yhit1=coord[[All,-1]]+k5[time[[-1]],h1,coord[[All,-1]]];

  yhit3=coord[[All,-1]]+k5[time[[-1]],h2,coord[[All,-1]]];

  pop++;

  hx=h1-(h2-h1)Min[distance[yhit1[[1;;N0]]][[All,2]]]/(Min[distance[yhit3[[1;;N0]]][[All,2]]]-Min[distance[yhit1[[1;;N0]]][[All,2]]]);

  he=hx-h1;

  yhit2=coord[[All,-1]]+k5[time[[-1]],h1+he, coord[[All,-1]]];

  dball=distance[yhit2[[1;;N0]]];

  e=Min[dball[[All,2]]];

  posmin=Position[dball[[All,2]],e][[1,1]];

  spmin=spdist[yhit2[[1+N0;;2N0]]][[posmin,2]];

  If[e>0,

  h1=h1+he;

  ,

  h2=h1+he;

  ];

  ];

  If[e>0,

  h1=h1-he;

  ,

  h2=h1-he;

  ];

  h1ball=hit1[dball];

  h2ball=hit2[h1ball];

  (*yhit1 - точка перед соударением, но очень близкая к нему. Скорости как до соударения*)

  yhit1=coord[[All,-1]]+k5[time[[-1]],h1,coord[[All,-1]]];

  If[Min[distance[yhit2[[1;;]]][[All,2]]]<=0,

  For[i=1,i<=Length[h2ball],i++,

  For[j=1,j<=Length[h2ball[[i]]],j++,

  yhit2[[h2ball[[i, j]]]]=yhit1[[h2ball[[i, j]]]];

  ];

  ];

  ];

  If[Min[distance[yhit2[[1;;N0]]][[All,2]]]<=0,

  For[i=1,i<=N0,i++,

  yhit2[[i]]=yhit1[[i]];

  ];

  ];

  newy=yhit2;

  For[i=1,i<=Length[h2ball],i++,

  For[j=1,j<=Length[h2ball[[i]]],j++,

  newy[[h2ball[[i, j]]+N0]]=yhit2[[h2ball[[i,-j]]+N0]];

  ];

  ];

  (*yhit3 - точка после соударения, но очень близкая к нему. В ней уже измененные скорости.*)

  yhit3=newy+k5[time[[-1]],he, newy];

  If[Min[distance[yhit3[[1;;]]][[All,2]]]<=0,

  For[i=1,i<=Length[h2ball],i++,

  For[j=1,j<=Length[h2ball[[i]]],j++,

  yhit3[[h2ball[[i, j]]]]=newy[[h2ball[[i, j]]]];

  ];

  ];

  ];

  If[Min[distance[yhit3[[1;;N0]]][[All,2]]]<=0,

  For[i=1,i<=N0,i++,

  yhit3[[i]]=newy[[i]];

  ];

  ];

  (* yhit2 - точка в момент соударения. Скорости соударяющихся шаров равны нулю.*)

  For[i=1,i<=Length[h2ball],i++,

  For[j=1,j<=Length[h2ball[[i]]],j++,

  yhit2[[h2ball[[i, j]]+Sum[n[[k]],{k,1,n2}]]]=0;

  (*Print[h2ball[[i, j]]+Sum[n[[k]],{k,1,n2}]];*)

  ];

  ];

  If[h1!=0,

  coord=Join[coord, Transpose[{yhit1}],2];

  time=Join[time,{time[[-1]]+h1}];

  ];

  coord=Join[coord, Transpose[{yhit2}],2];

  time=Join[time,{time[[-1]]+he}];

  root=Join[root,{time[[-1]]}];

  coord=Join[coord, Transpose[{yhit3}],2];

  time=Join[time,{time[[-1]]+he}];

  step=h1+he;

  ,

  coord=Join[coord, Transpose[{ynext}],2];

  time=Join[time,{time[[-1]]+step}];

  If[(nstep==10),step=2 step;pop2++; nstep=0; , nstep++];

  (*Print[time];*)

  ];

  tall=Join[tall,{SessionTime[]-t2}];

  ];

Print["Work Time1 = ", SessionTime[]-tstart]

f1={};

f2={};

f3={};

For[i=1,i<=Length[time],i++,

  f1=Join[f1,{{{time[[i]]},coord[[-6,i]],coord[[-3,i]]}}];

  f2=Join[f2,{{{time[[i]]},coord[[-5,i]],coord[[-2,i]]}}];

  f3=Join[f3,{{{time[[i]]},coord[[-4,i]],coord[[-1,i]]}}];

  ];

Solution={{x->Interpolation[f1],y->Interpolation[f2],и->Interpolation[f3]}};

f4={};

N0=;

For[j=1,j<=n2,j++,

  N1=;

  For[j2=1,j2<=n[[j]],j2++,

  f4={};

  For[i=1,i<=Length[time],i++,

  f4=Join[f4,{{{time[[i]]},coord[[j2+N1,i]],coord[[j2+N1+N0,i]]}}];

  ];

  Solution={Join[Solution[[1]],{шj, j2->Interpolation[f4]}]};

  ];

  ];

Print["Hit Number = ",Length[root]];

Solution=Solution[[1]]

Print["Work Time2 = ", SessionTime[]-tstart]