|
Схема нарисована без учёта сил трения и моментов.
Функция виброхода.
{
//---------------
// Счёт времени.
//---------------
t2 = t2 + dt; //Время процесса через внутренний счётчик функции.
t = j*dt; //Время процесса через внешний счётчик.
/*Счёт времени в модели выполняется двумя счётчиками. Счётчик t2 – считает абсолютное время независимое от работы процедур вывода результата. Счётчик j – обнуляется каждый цикл вывода результата. Он обеспечивает синхронизацию модели и процедур вывода графики. По умолчанию dt = 10 мкс, т. е. 100 кГц, что вполне достаточно как по времени счёта, так и по качеству.*/
//----------------------------------------------
// Управление трением колёс виброхода об опору.
//----------------------------------------------
/*Направление силы трения - Ftr всегда противоположно направлению движения корпуса - v1, относительно опоры - vo. Используем упрощенный линейный закон изменения трения от скорости.*/
vo = op*vopx; //Скорость движения опоры подвеса.
if(v1-vo<0) {Ftr = Ftrm - ktr*(v1-vo);};
if(v1-vo>0) {Ftr = -Ftrm - ktr*(v1-vo);};
if(v1-vo==0) {Ftr = 0;};
/*Сила трения покоя не может ускорять в реальном движении. Как только относительная скорость
скользящих тел стала равной нулю и тела оказались в относительном покое, происходит их
объединение, через сцепление трением. Эту ситуацию мы будем обрабатывать, считая для
простоты, что энергия сцепления равна нулю (не требуется разрыв механической связи
между телами),а «квант скорости» относительного движения автоматически включает
нужное направление силы трения. Статическое сцепление заменяется в этом случае
на динамическое. Т.е, относительное удержание тел происходит за счёт переключения направления
силы трения – «дребезга» (если нет попадания в ноль). Энергия «дребезга» мала и не нарушает баланса энергии. Фактически это следящий процесс. На графике силы трения это отображается в виде сплошных вертикальных полос. Это режимы работы, когда в системе явно не хватает кинетической энергии для непрерывного поступательного движения и происходит полное торможение. Не смотря на то, что статическое и динамическое сцепление принципиально отличаются,
это не приводит к недопустимому нарушению правильного взаимодействия тел (приятный презент от творца нашего мира). Для тех, кто хочет увидеть, что трение обрабатывается правильно, предусмотрен режим «тест трения», где (корпус + грузы без вращения) взаимодействуют с индикаторной тележкой через трение,
демонстрируя корректное поведение. Проблема в том, что от величины трения зависит
математическая формула закона взаимодействия, что и усложняет понимание этого процесса
на уровне простого логического анализа, а вот логика программы без проблем обрабатывает
этот кошмар, выдавая на-гора результат, который заставляет чесать «репу».*/
//Сила трения не может совпадать с направлением скорости.
if((v1-vo)*Ftr>0) {Ftr = 0;}; //Ограничивает спад трения из-за
//коэффициента Ktr нулём.
/*Используя коэффициент ktr Ньютон/(м/сек) можно задавать спад или нарастание трения относительно начального уровня Ftrm в зависимости от скорости, моделируя приблизительно сухое или гидродинамическое трение.*/
//--------------------
// Движение качелей.
//--------------------
sinop = sopx/hop; //Синус угла отклонения от вертикали.
cosop = sqrt(1 - sinop*sinop); //Косинус угла отклонения от вертикали.
tgop = sinop/cosop; //Тангенс угла отклонения от вертикали.
aopx = -((m + mop)*(g0 + aopy)*tgop + Ftr)/mop; //Ускорение качелей по
//горизонтали.
vopx = vopx + aopx*dt; //Скорость качелей по горизонтали.
sopx = sopx + vopx*dt; //Перемещение качелей по горизонтали.
vopy = vopx*tgop; //Скорость вертикального движения.
dhop = hop*(1.0 - cosop); //Высота подъёма.
aopy = (vopy - vopy1)/dt; //Вертикальное ускорение подвеса
vopy1 = vopy; //изменяет вес инерцоида.
//Энергия подвеса для контроля.
Aop = (mop + m)*g0*dhop + (mop + m)*vopy*vopy/2 + mop*vopx*vopx/2; //=0
//Сумма потенциальной энергии; вертикальной и горизонтальной энергии
//движения.
//-----------------------------
//-----------------------------
// Начало медленного полутакта.
//-----------------------------
if(sy > 0 && n==2);
{n = 1;
t15 = t2 - t16; //Интервал времени быстрого полутакта.
t14 = t2; //Момент времени начала медленного полутакта.
Ms = Mm; //Силовой момент в период медленного полутакта.
E01 = m2*vy*vy/2; //Конечная энергия быстрого полутакта.
if(kd == 1.0){vy = vy0min;}; //Включение свободного хода
//(нет тактирования).
E0 = m2*vy*vy/2; //Начальная энергия в медленном полутакте.
Ei = Ei + E01 - E0; //Учёт энергии тактирования грузов.
us0 = vy0min/R; //Начальная угловая скорость.
};
//Синхронизация работы медленного полутакта. Обеспечивает учёт энергии
//между полутактами.
//----------------------------
// Начало быстрого полутакта.
//----------------------------
if(sy < 0 && n==1)
{n = 2;
t13 = t2 - t14; //Время медленного полутакта.
t16 = t2; //Момент начала быстрого полутакта.
t12 = t2 - t1; t1 = t2; //Время такта.
ds = s1 - s11; s11 = s1; //Перемещение за такт через движение корпуса.
v12 = ds/t12; //Шаговая скорость по быстрому полутакту.
sc = sv - sv1; sv1 = sv; //Перемещение центра масс за такт.
Ms = Mb; //Силовой момент в период быстрого полутакта.
E01 = m2*vy*vy/2; //Конечная энергия медленного полутакта.
if(kd == 1.0){vy = vy0max;}; //Задание скорости вращения в быстром
//полутакте.
E0 = m2*vy*vy/2; //Начальная энергия в быстром полутакте.
Ei = Ei + E01 - E0; //Учёт энергии тактирования грузов.
us0 = -vy0max/R; //Начальная угловая скорость.
nc = nc++; //Счётчик тактов по быстрому полутакту.
};
//Синхронизация работы быстрого полутакта. Обеспечивает учёт энергии между //полутактами.
//---------------------------------------
// Кинематика через численные интегралы.
//---------------------------------------
v1 = v1 + a1*dt; //Скорость корпуса по оси х.
v2 = v2 + a2*dt; //Скорость грузов в направлении движения.
s1 = s1 + v1*dt; //Перемещение корпуса.
s2 = s2 + v2*dt; //Перемещение грузов в направлении движения.
vy = vy + ay*dt; //Поперечная скорость грузов.
sy = sy + vy*dt; //Поперечное перемещение грузов.
//---------------------------------------
//-------------------------------------------
// Деформация рычага и сила упругости в нём.
//-------------------------------------------
xd = s2 - s1; //Расстояние между корпусом и грузами в направлении
//движения.
L = sqrt(xd*xd + sy*sy); //Длинна деформированного стержня.
dr = L - R; //Величина деформация стержня.
dr2 = dr - dr1; //Скорость деформации стержня умноженная на dt.
dr4 = dr3 - dr2; //Ускорение деформации стержня на два диф. времени.
dr_ = dr + dr2 + dr4/2; //Скорректированная деформация стержня.
Fr = -kr*dr_; //Упругая сила при деформации стержня.
dr1 = dr;
dr3 = dr2;
/*Выражение для деформации (dr + dr2 + dr4/2) позволяет более гладко замкнуть вычислительный цикл с учётом предсказания процесса и уменьшить переходные процессы и обеспечить устойчивость работы модели.*/
//-------------------------------------------
//----------------
// функции угла.
//----------------
cosr = xd/L; //Косинус угла поворота грузов.
sinr = sy/L; //Соответственно синус.
//----------------
//----------------
// Динамика.
//----------------
Fgs = m*au; //Сила скатывания прибора на уклоне.
Fx = Fr*cosr; //Сила растяжения стержня вдоль направления движения.
Fy = Fr*sinr; //Сила растяжения стержня поперёк направления движения.
F1 = -Fx + Ftr + Ms*sinr/L + m1/m*Fgs; //Сила на корпус по оси х.
F2 = Fx - Ms*sinr/L + m2/m*Fgs; //Сила на грузы по оси х.
Fry = Fy + Ms*cosr/L; //Сила на грузы поперёк движения.
a1 = F1/m1; //Ускорение корпуса по оси х.
a2 = F2/m2; //Ускорение грузов по направлению движения.
ay = Fry/m2; //Ускорение грузов поперек движения, т.е. по оси y.
/*--------------------------------------------------
Конец алгоритма обработки процесса движения,
далее
только вспомогательные и проверочные функции.
//------------------------------------------------*/
//----------------------------
// Вычисление отката корпуса.
//----------------------------
if(v1>0 && nt == 0) {sot1 = s1; nt = 1;};
if(v1<0 && nt == 1) {sot = s1 - sot1; nt = 0;};
//----------------------------------------------
//---------------------------------------------------
// Вычисление кинематических параметров центра масс.
//---------------------------------------------------
a = (m1*a1 + m2*a2)/m; //Ускорение центра масс.
v = (m1*v1 + m2*v2)/m; //Скорость движения центра масс.
s = (m1*s1 + m2*s2)/m; //Перемещение центра масс.
//---------------------------
// Работа против сил трения.
//---------------------------
A = A - Ftr*v1*dt; //Потери энергии корпуса на торможение силой трения.
//-------------------------------------
//---------------------------------------------------
// Работа против сил тяжести (для опытов на горке).
//---------------------------------------------------
A1 = A1 - Fgs*v*dt; //Работа силы тяжести на уклоне.
//----------------------------------------------------
//-----------------------------------------
// Упругая энергия связи корпуса и грузов.
//-----------------------------------------
Ep = kr*dr_*dr_/2; //Можно смело пренебречь при выбранных параметрах
//жесткости стержня 5 тон/мм.
//-----------------------------------------
//-----------------------------------------------------------------------
// Работа момента сил (для опытов по выявлению влияния момента на работу
// виброхода.
//-----------------------------------------------------------------------
vvr = sqrt((v1 - v2)*(v1 - v2) + vy*vy); //Модуль скорости вращения.
us = vvr/L; //Угловая скорость вращения грузов.
Am = Am - Ms*us*dt; //Работа момента сил.
ug = (a1 - a2)*sinr/R + ay*cosr/R; //Угловое ускорение грузов.
//----------------------------------
//--------------------------------------------------------------
// Суммарная энергия системы. Проверка закона сохр. энергии.
//--------------------------------------------------------------
E = m1*v1*v1/2 + m2*v2*v2/2 + m2*vy*vy/2 + A + Am + Ep + A1 + Ei;
//=const;
//--------------------------------------------------------------
/*Всегда константа, как и положено, при учёте всех кинетических энергий и работ и конечно начальных и конечных условий в каждом полутакте (Ei).*/
//----------------------------
// Проверка импульса системы.
//----------------------------
Ft = Ft + Ftr*dt; //Импульс от силы трения.
p1 = m1*v1 + m2*v2 - Ft; //=0; Всегда ноль; m*v==Ft;
//---------------------------
//Только трение даёт внешний импульс, как и положено по классической теории виброхода.
//--------------------------------------------------------
// Проверка аналитического выражения для угловой скорости
// при движении вибратора без трения и силового момента.
//--------------------------------------------------------
k = m2/m;
f1 = us - us0/sqrt(1.0 - k*sinr*sinr));
//--------------------------------------------------------
//----------------------------------------------------------
// Проверка аналитического выражения для углового ускорения
// при движении вибратора без трения и силового момента.
//----------------------------------------------------------
f2 = ug - k*us0*us0*sinr*cosr/
(1.0 - k*sinr*sinr))/(1.0 - k*sinr*sinr));
//----------------------------------------------------------
};
Назад
|
