Продолжаем серию статей по автоматизации выполнения фигур пилотажа на малом ДПЛА. Настоящая статья имеет, прежде всего, образовательную цель: здесь мы покажем, как можно создать простейшую систему автоматического управления (САУ) на примере задачи выполнения фигуры пилотажа «бочка» при управлении летательным аппаратом только элеронами. Статья является второй в цикле публикаций «Пилотажный ДПЛА», рассказывающем о процессе постройки аппаратной и программной частей САУ в обучающей форме.
Введение
Итак, мы решили реализовать «бочку» в автоматическом режиме. Очевидно, что для автоматического выполнения фигуры необходимо сформулировать соответствующий закон управления. Процесс придумывания будет гораздо безболезненнее и быстрее, если использовать математическую модель движения ЛА. Проверка закона управления в лётном эксперименте хоть и возможна, но требует гораздо большего времени, к тому же может оказаться гораздо дороже в случае потери или повреждения аппарата.Так как при небольших углах атаки и скольжения самолёта его движение по крену практически не связано с движением в двух других каналах: путевом и продольном - для выполнения простой «бочки» достаточно будет построить модель движения только вокруг одной оси - оси ОХ связанной СК. По этой же причине, закон управления элеронами не будет существенно изменяться, когда дело дойдет до создания полной системы управления.
Модель движения
Уравнение движения ЛА вокруг продольной оси ОХ связанной СК крайне простое:Где - момент инерции относительно оси ОХ , а момент состоит из нескольких составляющих, из которых для реалистичного описания движения нашего самолёта достаточно рассмотреть только две:
Где - момент, обусловленный вращением ЛА вокруг оси ОХ (демпфирующий момент), - момент, обусловленный отклонением элеронов (управляющий момент). Последнее выражение записано в линеаризованной форме: момент крена линейно зависит от угловой скорости и угла отклонения элеронов с постоянными коэффициентами пропорциональности и соответственно.
Как известно (например, из Вики), линейному дифференциальному уравнению
Соответствует апериодическое звено первого порядка
Где - передаточная функция, - оператор дифференцирования, - постоянная времени, а - коэффициент усиления.
Как перейти от дифференциального уравнения к передаточной функции?
В нашем случае от параметров уравнения к параметрам передаточной функции можно перейти следующим образом (зная, что производная отрицательная):
Для апериодического звена постоянная времени равна времени, за которое выходная величина при единичном ступенчатом воздействии входной величины принимает значение, отличающееся от установившегося на величину ~5%, а коэффициент усиления численно равен установившемуся значению выходной величины при единичном ступенчатом воздействии:
В построенной модели движения есть два неизвестных параметра: коэффициент усиления и постоянная времени . Эти параметры выражаются через характеристики физической системы: момент инерции , а также производные момента крена и :
Таким образом, если известен момент инерции , то, определив параметры модели, по ним можно восстановить параметры системы.
Параметры модели. Момент инерции
Наш летательный аппарат состоит из следующих частей: крыло, фюзеляж с оперением, двигатель, аккумулятор (АКБ) и БРЭО:
К БРЭО относятся: плата автопилота, плата приёмника СНС, плата радиомодема, плата приёмника сигнала от управляющей аппаратуры, два регулятора напряжения, регулятор оборотов двигателя, а также соединительные провода.
В силу малого веса БРЭО, его вкладом в общий момент инерции можно пренебречь.
Как проводилась оценка величины момента инерции?
Оценку момента инерции можно провести следующим образом. Посмотрим на самолёт вдоль оси ОХ :
А затем представим его в виде следующей упрощённой модели:
Схема для расчёта момента инерции . Слева вверху - аккумулятор, справа внизу - двигатель. Двигатель и аккумулятор располагаются на оси фюзеляжа
Видно, что для создания модели были отброшены: киль, горизонтальное оперение, винт, а также БРЭО. При этом остались: фюзеляж, крыло, аккумулятор, двигатель. Измерив массы и характерные размеры каждой части, можно вычислить моменты инерции каждой части относительно продольной оси фюзеляжа:
Общее значение момента инерции ЛА относительно оси ОХ получим сложением моментов инерции частей:
Оценив вклад каждой из частей ЛА в общий момент инерции , получилось следующее:
- крыло - 96.3%,
- фюзеляж - 1.6 %,
- двигатель и аккумулятор - 2%,
Поэтому, несмотря на скромный вес (около 20% от общей взлётной массы ЛА), крыло имеет значительный момент инерции.
Параметры модели. Производные момента крена и
Вычисление производных момента крена – довольно трудная задача, связанная с расчётом аэродинамических характеристик ЛА численными методами или с помощью инженерных методик. Применение первых и вторых требует значительных временных, интеллектуальных и вычислительных затрат, которые оправданы при разработке систем управления большими самолётами, где стоимость ошибки всё же превышает затраты на построение хорошей модели. Для задачи управления БПЛА, масса которого не превосходит 2 кг, такой подход вряд ли оправдан. Другой способ вычисления этих производных - лётный эксперимент. Учитывая дешевизну нашего самолёта, а также близость подходящего поля для проведения такого эксперимента, для нас выбор был очевиден.Записав в автопилот прошивку для ручного управления и регистрации параметров, мы собрали самолёт и подготовили его к испытаниям:
В лётном эксперименте удалось получить данные по углу отклонения элеронов и угловой скорости вращения ЛА. Пилот управлял самолётом в ручном режиме, выполняя полёт по кругу, развороты и «бочки», а бортовое оборудование регистрировало и отправляло необходимую информацию на наземную станцию. В результате были получены необходимые зависимости: (град/с) и (б/р). Величина представляет собой нормированный угол отклонения элеронов: значение 1 соответствует максимальному отклонению, а значение −1 - минимальному:
Как теперь определить и из полученных данных? Ответ - измерить параметры переходного процесса по графикам и .
Как определялись коэффициенты k и T?
Коэффициент усиления определялся путём отнесения величины установившегося значения угловой скорости к величине отклонения элерона:
Зависимости угла отклонения элеронов и угловой скорости крена от времени, полученные в лётном эксперименте
На предыдущем рисунке участкам установившегося значения угловой скорости приблизительно соответствуют, например, отрезки вблизи моментов времени 422, 425 и 438 с (отмечены тёмно-красным на рисунке).
Постоянная времени определялась из тех же графиков. Для этого найдены участки резкого изменения угла отклонения элерона, а затем измерено время, за которое угловая скорость принимает значение, отличающееся от установившегося значения на 5%.
Результат определения значений постоянной времени и коэффициента усиления следующий: , . Этим значениям коэффициентов при известном значении момента инерции соответствуют следующие значения производных момента крена:
Верификация модели
Итак, построив модель, основой которой является апериодическое звеноМожно провести её верификацию, подав на вход сигнал , полученный из лётного эксперимента, и сравнить выходной сигнал модели с величиной , также полученной в эксперименте.
Как проводилось моделирование?
Инструмент для проведения моделирования мы выбирали, прежде всего, основываясь на возможности повторения результатов широким кругом читателей: это прежде всего означает, что программа должна быть в общественном доступе. В принципе задачу моделирования поведения апериодического звена первого порядка можно решить и создав свой собственный инструмент с нуля. Но так как в дальнейшем модель будет усложняться, то создание своего инструмента может отвлечь от основной задачи - создания САУ пилотажного ДПЛА. С учётом принципа открытости инструмента мы выбрали JSBsim .
В предыдущем разделе мы получили значения коэффициентов и . Используем их для моделирования движения самолета. Из прошлой статьи мы помним, что конфигурация модели самолета в JSBsim
задается при помощи XML
файла. Создадим собственную модель:
Так как мы строим модель движения аппарата только по крену, оставим многие из секций файла пустыми. В файле модели последовательно задаются следующие характеристики.
Геометрические размеры самолёта задаются в секции metrics : площадь крыла, размах, длина средней аэродинамической хорды, площадь горизонтального оперения, плечо горизонтального оперения, площадь вертикального оперения, плечо вертикального оперения, положение аэродинамического фокуса.
Массовые характеристики самолёта задаются в секции mass_balance : тензор инерции самолета, вес пустого самолета, положение центра масс.
Стоит отметить, что абсолютные положения аэродинамического фокуса и центра масс самолета не участвуют в расчете динамики аппарата, важно их относительное расположение.
Далее следуют секции, описывающие характеристики шасси самолета и его силовой установки.
В следующей секции, ответственной за систему управления , заполним канал, ответственный за управление по крену: укажем единственный вход fcs/aileron-cmd-norm , величина которого будет нормирована от -1 до 1.
Аэродинамические характеристики задаются в секции aerodynamics : силы задаются в скоростной системе координат, моменты - в связанной. Нас интересует момент по крену. В секции axis name=«ROLL» задаются функции, которые определяют момент сил от различных составляющих проекции момента аэродинамических сил на ось OX связанной системы координат. В нашей модели таких составляющих две. Первая составляющая - демпфирующий момент, который равен произведению угловой скорости на определенный ранее коэффициент . Вторая составляющая - это момент от элеронов при фиксированной скорости полёта: он равен произведению ранее определенного коэффициента на величину отклонения элеронов.
Стоит отметить, что при определении коэффициента была использована размерная величина . В наших полетных данные величина угловой скорости измерялась в градусах в секунду, тогда как в JSBSim используются радианы в секунду, поэтому коэффициент должен быть приведен к нужной нам размерности, т. е. разделен на 180 градусов и умножен на радиан. Записываем эти составляющие момента аэродинамических сил внутри функций произведения product . При моделировании результат выполнения всех функций суммируется и получается значение проекции аэродинамического момента на соответствующую ось.
Проверить созданную модель можно на экспериментальных данных, полученных при лётных испытаниях. Для этого создаем скрипт следующего содержания:
где точками обозначены пропущенные данные. В знакомом нам по предыдущей статье файле скрипта появился новый тип события («Time Notif»
), который позволяет задавать непрерывное изменение параметра по времени. Зависимость параметра от времени задана табличной функцией. JSBSim
линейно интерполирует значение функции между табличными данными. Процедура верификации модели движения по крену заключается в исполнении данного скрипта на созданной модели и сравнения результатов с экспериментальными.
Результат верификации показан на рисунке:
Как видно из рисунка, совпадение модели с реальностью чуть менее, чем полное.
Синтез управления для выполнения «бочки»
Получив модель, несложно определить, на какую величину и как долго необходимо отклонять элероны, чтобы выполнить «бочку». Одним из вариантов является следующий алгоритм отклонения:
Наличие отрезков длительностью по 0.1 с в начале и в конце алгоритма отклонения элеронов моделирует инертность сервопривода, который не может отклонить поверхности мгновенно. Модель показывает, что при таком законе отклонения элеронов ЛА должен выполнить один полный оборот вокруг оси ОХ , проверим?
Лётный эксперимент
Полученный закон управления элеронами был запрограммирован в автопилот, установленный на самолёт. Идея эксперимента проста: вывести самолёт в горизонтальный полёт, после чего задействовать полученный закон управления. В случае если реальное движение ЛА по крену соответствует созданной модели, то самолёт должен выполнить «бочку» - один полный оборот на 360 градусов.Отдельно выражаем благодарность нашему верному пилоту за труд, профессионализм и удобный багажник на приоре-универсал!
В процессе проведения эксперимента стало ясно, что модель движения по крену была построена удачно - самолёт выполнял одну «бочку» за другой, как только пилот задействовал запрограммированный закон управления. На следующем рисунке показана угловая скорость , записанная в процессе эксперимента, и полученная по результатам моделирования, а также угол крена и тангажа из лётного эксперимента:
А на следующем рисунке показаны зарегистрированные в лётном эксперименте сигналы на элероны, руль высоты (РВ) и руль направления (РН):
Вертикальными линиями обозначены моменты начала и окончания выполнения «бочки». Из рисунков видно, что в процессе выполнения «бочки» пилот не вмешивается в управление рулём высоты и рулём направления, также видно, что угол тангажа неизменно стремится уменьшиться при выполнении «бочки» - самолёт затягивает в пикирование, как это и было предсказано по результатам моделирования в авиасимуляторе (см. статью «Пилотажный ДПЛА. Как правильно сделать бочку»). Если внимательно рассмотреть предыдущие графики, то станет видно, что третья «бочка» даже не была закончена, потому что пилот вмешался в управление, чтобы вывести самолёт из пикирования: настолько сильно изменяется угол тангажа при выполнении «бочки» одними элеронами.
Замечания
Выводы
В результате проведённой работы мы показали один из способов создания модели движения ДПЛА по угловой скорости . В лётном эксперименте было доказано, что созданная модель движения вполне соответствует моделируемому объекту. На основании разработанной модели получен закон программного управления, позволяющий выполнять «бочку» в автоматическом режиме. Также мы убедились, что выполнить правильную «бочку» одними элеронами не получится, а также наглядно это продемонстрировали.Следующим этапом будет доработка закона управления путём добавления обратной связи, а также включение в управление руля высоты. Последнее потребует создание модели продольного движения нашего самолёта. По результатам работы выйдет следующая публикация.
Раздел программы: “Формализация и моделирование”.
Тема урока: “Моделирование движения”.
Тип урока: урок изучения нового материала.
Вид урока: комбинированный.
Технология: личностно-ориентированная.
Время проведения: второй урок по теме “Моделирование графических объектов”.
Цели урока:
- развитие представлений о моделировании как методе познания;
- формирование системно-информационного подхода к анализу окружающего мира;
- формирование общеучебных и общенаучных навыков работы с информацией.
Задачи урока:
- Воспитательная – развитие познавательного интереса, воспитание информационной культуры, воспитание умения четко организовать самостоятельную работу.
- Учебная – изучить и закрепить прием моделирования динамических объектов.
- Развивающая – развитие системно-конструктивного мышления, расширение кругозора.
Методы: словесные, наглядные, практические.
Организационные формы работы: фронтальные, индивидуальные.
Материально-техническая база:
- презентация “Моделирование движение”;
- комплекс: демонстрационный экран и компьютер с ОС Windows-9x с установленным MS Office 2000;
- компьютеры с программной средой Turbo Pascal 7.0.
Межпредметная связь : математика.
1. Подготовка к уроку
Для урока подготовлена презентация с помощью Power Point с целью визуализации информации по ходу объяснения нового материала. (Приложение1.ppt)
План урока:
| Содержание этапа урока | Вид и формы работы |
| 1. Организационный момент | Приветствие |
| 2. Мотивационное начало урока | Постановка цели урока. Фронтальный опрос |
| 3. Изучение нового материала | Использование слайдов, работа в тетради |
| 4. Этап закрепления, проверки полученных знаний | Практическая работа: компьютерный эксперимент по проверке программы |
| 5. Этап систематизации, обобщения изученного | Самостоятельная работа за компьютером:
компьютерный эксперимент по исследованию
модели. Работа в тетради |
| 6. Подведение итогов, домашнее задание | Работа в тетради |
Ход урока
2. Организационный момент
3. Мотивационное начало урока. Постановка цели урока
Учитель: На прошлом занятии мы строили статичное изображение.
Вопрос: Какая модель называется статической? Какая модель называется динамической?
Ответ: Модель, описывающая состояние объекта, называется статической. Модель, описывающая поведение объекта, называется динамической.
Учитель: Сегодня продолжим тему построение изображений, но уже в динамике, т.е. объект будет изменять свое положение на плоскости во времени. Начну с демонстрации имеющейся у меня копилки программ, которые хорошо иллюстрируют тему сегодняшнего урока. (Начинается показ через запуск программ на ЯП Паскаль “Хаотичное движение” , “Полет в космосе” , “Движение колеса” (Приложение2.pas, Приложение3.pas, Приложение4.pas). Изучение модели движения мы и посвятим сегодняшний урок.
В классе на экране тема урока “Моделирование движения”.
Запишите тему сегодняшнего урока.
Учитель: Условие задачи зафиксируйте в тетради.
Для решения задачи смоделируем процесс движения сначала через описательную модель, затем формализованную и, наконец, компьютерную, чтоб можно было реализовать модель на компьютере.
Для начала давайте обсудим вопрос, что значит создать анимацию (иллюзию движения какого-либо объекта)?
Обсуждение. Заслушивание всех вариантов ответов, вплоть до невозможных.
Предполагаемый ответ: Если это как в мультипликации, то, наверное, это должно быть в виде набора статичных изображений сменяющих друг друга через какое-то время.
Учитель: Хорошо.
4. Изучение нового материала
Словесную описательную модель нашей задачи можно сформулировать так:
Учитель вслух комментирует описательную модель, просит учащихся зафиксировать ее в тетради.
Учитель: Перейдем к формализованной модели, и раз это изображение, то воспользуемся системой координат компьютера и схематично изобразим, как это должно выглядеть.
Учащиеся фиксируют эту модель в тетрадь.
Учитель: А вот как это будет выглядеть на экране (слайд выполнен с анимацией, круг проделывает движение слева на право).
Учащиеся наблюдают.
Учитель: Запишем словесный алгоритм реализации нашей модели. Ясно, что для повторения многократного изображений круга каждый раз в новой точке экрана понадобится цикл.
Вопрос: Какой цикл лучше использовать?
Ответ: For-To-Do.
Вопрос: Какая процедура поможет нам нарисовать круг белого цвета? Черного цвета?
Ответ: SetColor(15) и Circle(X,Y,R), затем SetColor(0) и Circle(X, Y, R).
Вопрос: Как осуществить задержку времени на пример на 100 м/сек?
Ответ: Delay(100).
Учитель: Правильно.
Демонстрируем слайды с 8 по 10. Учащиеся сверяют свои ответы с правильными.
Учитель: А теперь запишите всю программу целиком у себя в тетради.
Выдерживаем паузу 5–7 минут. Затем даем возможность свериться с образцом.
Движение автомобиля рассматривается как плоскопараллельное движение твердого тела по горизонтальной поверхности (рис. 1). В общем случае движение автомобиля описывается следующей системой дифференциальных уравнений:
где - вектор ускорения центра масс автомобиля; m - масса автомобиля; fi - вектор силы сопротивления прямолинейному движению i-го колеса; i - вектор взаимодействия с грунтом i-го колеса; w - вектор силы сопротивления воздуха; J z - момент инерции автомобиля относительно оси z; M nki - момент сопротивления повороту i-го колеса.
Ускорение определяется как
где dV/dt - относительная производная скорости центра масс автомобиля. Проекции скоростей в координатах x`, y`, z`:
Учитывая, что:
можно записать следующую систему уравнений:
Данную систему уравнений решим с помощью пакета DEE (Differential Equation Editor) входящего в состав Simulink. Для этого записываем уравнения в нормальной форме Коши и настраиваем входные данные:
Рисунок 6. Решатель систем дифференциальных уравнений
Входными данными будут являться выходы с предыдущих блоков. Общий вид модели представлен на следующем рисунке:
Рисунок 7. Модель транспортного средства с колесной формулой 4х4
Результаты моделирования представим графически:
Рисунок 8. Траектория движения автомобиля
Результаты моделирования представляют собой траекторию движения автомобиля в форме окружности, что говорит об адекватности данной модели. Данная работа может послужить фундаментом для дальнейших перспективных исследований в области разработки систем автоматического управления движением автомобиля, в том числе систем активной безопасности.
Сокращенный перевод-пересказ статьи: .
Введение
В статье описаны основы подхода к моделированию движения объектов, который удобно применять в компьютерных играх. Этот подход прост, реализующие его программы работают быстро и достаточно устойчиво. Кроме того, для понимания его основ не требуется особых математических знаний (хотя сам подход имеет под собой твердые математические основания). С его помощью можно моделировать движение тканей, мягких и твердых тел, а также тел с учетом связей.
Физическое моделирование, т. е. моделирование движения персонажей, основанное на законах физики (а точнее - механики), изучается достаточно давно. В литературе (см. и др.) предлагаются различные подходы, много усилий вложено в создание точных и надежных алгоритмов. Точные методы моделирования движения известны в физике уже давно. Однако для игр и систем виртуальной реальности, точность - не самое главное достоинство (хотя хорошо, когда она есть). Гораздо важнее правдоподобие (программист может искажать модель реальности сколько угодно, лишь бы при этом ему удалось увлечь игрока) и скорость выполнения (на выполнение расчетов по моделированию движения отводится лишь часть времени, которое длится кадр анимации). В случае физического моделирования, термин "правдоподобие" подразумевает также устойчивость: нельзя признать удачным метод моделирования, при котором тела проникают сквозь препятствия или подпрыгивают, когда должны лежать неподвижно. Методы, описанные в данной работе, разрабатывались в первую очередь для достижения правдоподобия и скорости расчета. Они имеют высокую производительность и достаточно просты в реализации (по крайней мере, в сравнении с другими методами, решающими те же задачи).
Рассматриваемый метод является итеративным, так что, начиная с определенного шага, он может быть остановлен в любой момент. Это дает возможность выбирать между точностью расчетов и затраченным временем: если некоторая величина погрешности считается приемлемой, коду можно разрешить работать быстрее; причем величина погрешности может подбираться адаптивно во время выполнения. Метод также обрабатывает столкновения и контакты покоя (resting contact) и справляется с моделированием стопки сложенных друг на друга тел, что для многих физических движков является проблемой.
Успех применения метода складывается из правильного комбинирования и использования преимуществ нескольких техник:
- метода численного интегрирования Верле;
- обработки столкновений и проникновения тел при помощи проецирования (by projection);
- простого решателя связей (constraint solver), использующего релаксацию;
- аппроксимации квадратного корня, повышающей скорость вычислений;
- моделирования твердых тел, как частиц, соединенных связями.
К каждой из указанных техник мы дадим короткое пояснение. При написании этого документа, автор старался сделать его доступным для максимально широкой аудитории, не потеряв при этом существенной информации, необходимой для реализации. Это означает, что математические объяснения и обоснования сведены к минимуму, если только они не имеют решающего значения для понимания предмета. Цель работы состоит в том, чтобы продемонстрировать возможности реализации довольно продвинутого и устойчивого метода физического моделирования, не утонув в математических тонкостях.
Содержание организовано следующим образом. В разделе 2, описано представление системы частиц без использования скорости. Такое представление имеет ряд преимуществ, из которых наиболее существенными являются устойчивость и простота реализации связей и других ограничений (constraints). В раздел 3 описано, как происходит обработка столкновений. Затем, в разделе 4, система частиц дополняется связями, позволяющими моделировать движение ткани. Раздел 5 объясняет, как настроить систему связанных между собой частиц для моделирования твердого тела. Далее, в разделе 6, показано как реализовать соединения между телами (в частности, шарниры). Раздел 7 содержит различные короткие заметки и некоторый опыт по реализации трения.
В дальнейшем, векторы обозначаются буквами со стрелочками, а их компоненты - нижними индексами: \(\vec{x}=(x_1,x_2,x_3)\) .
Метод интегрирования Верле
Сердцем симуляции (т. е. имитации физического процесса при помощи компьютерной системы) является система частиц. Обычно, при реализации такой системы предполагается, что каждая частица имеет две основные характеристики: координату (положение, position) \(\vec{x}\) и скорость \(\vec{v}\) . Тогда, новые значения координат \(\vec{x}^\prime\) и скорости \(\vec{v}^\prime\) вычисляются по формулам
$$ \begin{aligned} \vec{x}^\prime &= \vec{x} + \vec{v} \Delta t, \\ \vec{v}^\prime &= \vec{v} + \vec{a} \Delta t, \end{aligned} $$
где \(\Delta t\) - шаг по времени, \(\vec{a}\) - ускорение, вычисленное в соответствии со 2-м законом Ньютона \(\vec{f}=m \vec{a}\) (где \(\vec{f}\) - суммарная сила, действующая на частицу). Приведенные формулы реализуют простейший метод численного интегрирования - метод Эйлера.
Мы рассмотрим другое описание частицы, в котором скорость не используется: вместо хранения положения и скорости каждой частицы, мы будем сохранять текущее положение частицы \(\vec{x}\) и ее положение на предыдущем шаге интегрирования \(\vec{x}^{*}\) . Предполагая шаг интегрирования постоянным, получим следующие формулы для вычисления новых значений:
$$ \begin{aligned} \vec{x}^\prime &= 2\vec{x} - \vec{x}^{*} + \vec{a} \Delta t^2, \\ \vec{x}^{*} &= \vec{x}. \end{aligned} $$
Этот способ численного интегрирования называется методом Верле (см. ) и активно используется в молекулярной динамике.
Метод Верле опирается на приближенную формулу вычисления второй производной
$$ \frac{\Delta^2 \vec{x}}{\Delta t^2} = \frac{ \frac{\vec{x}^\prime - \vec{x}}{\Delta t} - \frac{\vec{x} - \vec{x}^{*}}{\Delta t} }{\Delta t} = \frac{\vec{x}^\prime - 2\vec{x} + \vec{x}^{*}}{\Delta t^2} = \vec{a} $$
Такое приближение не является самым точным (есть и более совершенные методы численного интегрирования), зато оно устойчиво и работает быстро. Уменьшая коэффициент 2 до, скажем, 1.99, мы тем самым вводим силу сопротивления, рассеивающую энергию системы. Отметим также, что \(\vec{x}-\vec{x}^{*}\) - это расстояние, пройденное за последний шаг интегрирования (\(\vec{v}\Delta t\) ).
В конце шага интегрирования текущее положение каждой частицы \(\vec{x}\) сохраняется в соответствующей переменной \(\vec{x}^{*}\) для использования на следующем шаге. Если частиц в системе много, то вместо копирования их координат удобно использовать перенаправление указателей.
Код, реализующий описанные выше идеи, может выглядеть так (класс Vector3 содержит все необходимые операции над векторами)
class ParticleSystem { Vector3 m_x [ NUM_PARTICLES ]; // Текущее положение Vector3 m_oldx [ NUM_PARTICLES ]; // Предыдущее положение Vector3 m_a [ NUM_PARTICLES ]; // Суммарная сила (ускорение) Vector3 m_vGravity ; // Гравитация float m_fTimeStep ; public : void TimeStep (); private : void Verlet (); void SatisfyConstraints (); void AccumulateForces (); // (конструкторы, инициализацию полей и т.п. опустим) }; // шаг интегрирования методом Верле void ParticleSystem :: Verlet () { for (int i = 0 ; i < NUM_PARTICLES ; i ++ ) { Vector3 & x = m_x [ i ]; Vector3 temp = x ; Vector3 & oldx = m_oldx [ i ]; Vector3 & a = m_a [ i ]; x += x - oldx + a * fTimeStep * fTimeStep ; oldx = temp ; } } // суммирование сил, действующих на каждую частицу void ParticleSystem :: AccumulateForces () { // Все частицы находятся под действием гравитации for (int i = 0 ; i < NUM_PARTICLES ; i ++ ) m_a [ i ] = m_vGravity ; } // проверка соблюдения наложенных связей void ParticleSystem :: SatisfyConstraints () { // Сейчас нам не важно, как это реализовано. } // шаг расчета void ParticleSystem :: TimeStep () { AccumulateForces (); Verlet (); SatisfyConstraints (); }
Пока все описанное выше выглядит не очень-то впечатляюще. Преимущества этого подхода станут ясны, как только мы перейдем к использованию связей и к описанию твердых тел.
Попробуйте задать \(\vec{a}=(0,0,1)\) и начальные условия \(\vec{x}=(1,0,0)\) , \(\vec{x}^*=(0,0,0)\) . Вычислите вручную несколько шагов и посмотрите, что получится.
Столкновения и обработка контактов при помощи проецирования
Способы обработки контактов между телами, основанные на использовании штрафных функций (penalty-based schemes), предполагают, что в месте контакта, где возможно проникновение тел друг в друга, нужно вставить пружину для моделирования этого контакта. Такой подход прост в реализации, но приводит к ряду серьезных проблем. В частности, очень трудно подобрать жесткость пружины так, чтобы, с одной стороны, объекты не проникали друг в друга слишком глубоко, а с другой - чтобы система не потеряла устойчивость из-за слишком большой жесткости пружин. Еще один подход к обработке столкновений заключается в том, что при обнаружении столкновения время "отматывается" назад, вплоть до точного момента контакта тел (например, с помощью бинарного поиска), затем корректируются положения и скорости тел (по известным из курса физики формулам для столкновений), после чего расчет начинается заново с этого момента времени. И так - для каждого столкновения. Не слишком экономный подход, если предполагается моделировать в реальном времени движение множества тел.
Здесь мы рассмотрим другой подход. Проникшие в препятствие частицы мы будем проецировать за пределы препятствия. Под проецированием мы понимаем перемещение частицы, настолько малое, чтобы только освободить ее от препятствия. Как правило, это предполагает перемещение частицы по направлению нормали к поверхности контакта (препятствия) - отсюда и происхождение термина "проецирование".
Рассмотрим следующий пример. Пусть наш "мир" представляет собой внутренность куба размером (0,0,0)--(1000,1000,1000) и, кроме того, коэффициент восстановления (restitution coefficient) частиц равен нулю (т. е. столкнувшиеся с поверхностью куба частицы не отражаются от нее). Чтобы координаты частиц оставались внутри куба, запишем следующий код, реализующий проецирование:
// Заставляет частицы оставаться внутри куба void ParticleSystem :: SatisfyConstraints () { for (int i = 0 ; i < NUM_PARTICLES ; i ++ ) { // Для всех частиц Vector3 & x = m_x [ i ]; x = vmin (vmax (x , Vector3 (0 , 0 , 0 )), Vector3 (1000 , 1000 , 1000 )); } }
(vmax представляет собой покомпонентную операцию вычисления максимума, а vmin - аналогичное вычисление минимума). Это код позволяет обработать как столкновения, так и неподвижные контакты (resting contact) (т. е. случаи, когда точка покоится на поверхности куба), и сохраняет положение всех частиц внутри куба. Прелесть метода Верле в том, что соответствующие изменения в значения скоростей вносятся автоматически. В последующих вызовах TimeStep() скорость уже будет скорректирована так, чтобы не содержать составляющей, перпендикулярной поверхности куба (что соответствует нулевому значению коэффициента восстановления). См. рис.1.
Попробуйте сами проделать эти вычисления - и вы уведите, что не нужно обнулять скорость в направлении, перпендикулярном стенке куба - это происходит "само собой". Описанное может показаться тривиальным, если ограничится моделированием частиц, но основные преимущества метода Верле проявятся, как только мы перейдем к рассмотрению связей и связанных твердых тел. Т. е., прямо сейчас.
Обработка нескольких одновременно наложенных связей методом релаксации
Модель ткани обычно представляет собой систему частиц, соединенных пружинами. Дифференциальные уравнения такой системы построить несложно. Но одно дело построить, и совсем другое - решить. При этом всплывают все те проблемы, что мы имели при использовании штрафных функций: слишком жесткие пружины приводят к тому, что система уравнений сама становиться "жесткой" (stiff system), а это приводит к неустойчивости, если используются простейшие методы численного интегрирования или к медленной работе, если используются методы более совершенные - в обоих случаях головная боль обеспечена. И наоборот, слишком мягкие пружины приводят к тому, что ткань будет выглядеть нереалистично упругой.
Однако самое интересное произойдет, если устремить жесткость пружин к бесконечности: система вдруг становится разрешима даже для весьма простого (и быстрого) метода интегрирования, оставаясь при этом устойчивой. Но прежде чем мы продолжим разговор о ткани, давайте вернемся к предыдущему примеру. Куб, с которым мы имели дело, можно рассматривать как совокупность односторонних (unilateral) связей (т. е. связей, записываемых в форме неравенств) - по одной для каждой стороны куба - которые должны выполняться все время моделирования.
\begin{equation} x_i \geq 0 \ \text{and}\ x_i \leq 1000 \quad (i=1,2,3). \label{eq:C1} \end{equation}
В рассмотренном примере, для того чтобы соблюсти ограничения, накладываемые связями (чтобы частицы оставались внутри куба), достаточно просто спроектировать координаты "вылезших" частиц на поверхность куба. Эта идея описывается следующим псевдокодом:
// Псевдокод, позволяющий выполнить ограничения (1) for i = 1 , 2 , 3 set xi = min { max { xi , 0 }, 1000 }
Это можно представить себе так, будто частица и поверхность препятствия соединены бесконечно жесткой пружиной, которая, в случае удлинения, мгновенно возвращается к своей нормальной длине, равной нулю.
Расширим нашу модель, добавив к ней стержень длиной 100. Для этого нам понадобится задать две частицы (\(\vec{x}_1\) и \(\vec{x}_2\) ) и потребовать, чтобы расстояние между ними всегда было равно 100. Математическая запись этой двухсторонней (bilateral) связи имеет вид:
\begin{equation} |\vec{x}_2-\vec{x}_1| = 100. \label{eq:C2} \end{equation}
Даже если в начальный момент времени положения частиц удовлетворяют условиям \eqref{eq:C2}, то уже в следующий момент эти условия, скорее всего, выполняться не будут. Для того чтобы получить корректное значение расстояния, переместим частицы, проецируя их на множество решений, описанных \eqref{eq:C2}. Для этого частицы либо отодвигаются друг от друга, либо подтягиваются ближе, в зависимости от того, мало или велико расстояние, полученное численным интегрированием. См. рис.2.
Перемещение частиц для исправления расстояния, не удовлетворяющего ограничению \eqref{eq:C2}
Псевдокод, реализующий выполнение условий \eqref{eq:C2}:
Delta = x2 - x1 ; deltalength = sqrt (delta * delta ); diff = (deltalength - restlength ) / deltalength ; x1 -= delta * 0.5 * diff ; x2 += delta * 0.5 * diff ;
Заметим, что delta - вектор, а delta*delta - скалярное произведение. Этот псевдокод будет раздвигать или сдвигать частицы так, чтобы добиться требуемого расстояния между ними. И вновь мы можем рассматривать это как бесконечно жесткую пружину, мгновенно возвращающую себе нормальную длину, равную 100.
Теперь предположим, что, помимо условия \eqref{eq:C2}, должно выполняться и условие \eqref{eq:C1} (частицы обязаны находиться внутри куба). Может оказаться, что при попытке соблюсти условие \eqref{eq:C2}, какая-то из частиц стержня нарушит требования \eqref{eq:C1} (стержень будет торчать из куба). Можно, конечно, снова спроектировать частицу-нарушителя на поверхность куба, выполняя \eqref{eq:C1}, но тогда будет нарушено уже \eqref{eq:C2}.
Чтобы удовлетворить одновременно требованиям \eqref{eq:C1} и \eqref{eq:C2} нам нужно решить систему уравнений. Мы это и сделаем, но не напрямую: просто будем повторять два фрагмента псевдокода друг за другом какое-то количество раз, в надежде, что результат окажется полезным. Такой подход реализован в следующем коде:
// реализация моделирования стержня внутри куба void ParticleSystem :: SatisfyConstraints () { for (int j = 0 ; j < NUM_ITERATIONS ; j ++ ) { // Сначала выполним условия (1) for (int i = 0 ; i < NUM_PARTICLES ; i ++ ) { // Для всех частиц Vector3 & x = m_x [ i ]; x = vmin (vmax (x , Vector3 (0 , 0 , 0 )), Vector3 (1000 , 1000 , 1000 )); } // Теперь удовлетворим (2) Vector3 & x1 = m_x [ 0 ]; Vector3 & x2 = m_x [ 1 ]; Vector3 delta = x2 - x1 ; float deltalength = sqrt (delta * delta ); float diff = (deltalength - restlength ) / deltalength ; x1 -= delta * 0.5 * diff ; x2 += delta * 0.5 * diff ; } }
(здесь опущена инициализация частиц). Хотя такой способ "тупого" повторения и может показаться несколько наивным, тем не менее, он сходится к решению, которое мы ищем! В математике он называется методом релаксации (или Якоби, или Гаусса-Зейделя - в зависимости от того, как именно вы это делаете, см. ). Он работает, последовательно удовлетворяя отдельные ограничения, и сходится к глобальной конфигурации, которая удовлетворяет всем ограничениям одновременно. Этот метод очень полезен в ситуациях, когда должны одновременно выполняться несколько независимых ограничений.
Число необходимых итераций зависит от моделируемой системы и характера движения. Можно сделать выбор этого числа адаптируемым, измеряя изменение, произошедшее относительно предыдущей итерации. Если мы остановим итерации слишком рано, результат будет недостаточно точным, но, благодаря методу Верле, в следующем кадре он, вероятно, будет чуть лучше, а в следующем кадре - еще лучше, и т. д. Это означает, что преждевременная остановка релаксации не уничтожит анимацию полностью, но сделает картинку более дерганой.
Моделирование ткани
Тот факт, что связь типа стержня можно рассматривать как очень жесткую пружину, позволяет использовать этот вид связей для моделирования тканей. Предположим, например, что ткань представляется шестиугольной сеткой, состоящей из треугольников. Каждый узел сетки представляет собой частицу, а каждая грань - связь типа стержня, соединяющую частицы (нормальная длина стержня равна расстоянию между соединяемыми им узлами).
Функция HandleConstraints() , отвечающая за обработку связей, использует релаксацию по всем ограничениям. Цикл релаксации может повторятся несколько раз. Однако, чтобы получить хорошо выглядящую анимацию, в большинстве случаев достаточно всего одной итерации. Это означает, что расход времени в симуляции ткани зависит в основном от того, как долго выполняются \(N\) операций вычисления квадратного корня и \(N\) делений (где \(N\) - число ребер в сетке, моделирующей ткань). Ниже мы покажем один трюк, позволяющий избавиться от вычисления квадратного корня. Но сначала рассмотрим как выглядит обработка ограничений.
// Реализация моделирования тканей struct Constraint { int particleA , particleB ; float restlength ; }; // Предположим, что массив ограничений m_constraints уже существует void ParticleSystem :: SatisfyConstraints () { for (int j = 0 ; j < NUM_ITERATIONS ; j ++ ) { for (int i = 0 ; i < NUM_CONSTRAINTS ; i ++ ) { Constraint & c = m_constraints [ i ]; Vector3 & x1 = m_x [ c . particleA ]; Vector3 & x2 = m_x [ c . particleB ]; Vector3 delta = x2 - x1 ; float deltalength = sqrt (delta * delta ); float diff = (deltalength - c . restlength ) / deltalength ; x1 -= delta * 0.5 * diff ; x2 += delta * 0.5 * diff ; } // Прикрепим одну из частиц, составляющих ткань, к началу координат m_x [ 0 ] = Vector3 (0 , 0 , 0 ); } }
Теперь обсудим, как избавиться от вычисления квадратного корня. Если все ограничения соблюдены (ну, или почти соблюдены), то, как мы уже знаем, результат вычисления квадратного корня стремится к \(r\) - нормальной длине связи (стержня). Мы используем этот факт чтобы получить приближенное выражение для функции квадратного корня. Заменим функцию \(\sqrt{x}\) членом 1-го порядка из ее разложения в ряд Тейлора в окрестности длины \(r\) (это эквивалентно одной итерации метода Ньютона-Рафсона с начальным приближением \(r\) ). После некоторых преобразований, получим следующий псевдокод:
// Псевдокод для соблюдения ограничений (2), использующий приближение sqrt delta = x2 - x1 ; delta *= restlength * restlength / (delta * delta + restlength * restlength ) - 0.5 ; x1 -= delta ; x2 += delta ;
Обратите внимание, что если расстояние уже удовлетворяет ограничениям (т. е., если |delta|=restlength), то мы получим delta равное (0,0,0) и никаких изменений не произойдет.
Теперь при обработке каждой связи мы больше не используем квадратные корни. Кроме того, квадрат значения restlength * restlength можно вычислить заранее. Трудоемкие операции сокращены до выполнения \(N\) делений за кадр (и доступа к соответствующей памяти) - трудно придумать что-то, работающее существенно быстрее.
Ограничения не обязательно будут удовлетворены за одну итерацию, но, благодаря методу Верле, система быстро сходится к правильному состоянию (когда все ограничения соблюдены) - буквально за несколько кадров. На самом деле, использование только одной итерации и аппроксимации квадратного корня снимает проблему жесткости системы уравнений, которая обязательно проявилась бы у системы с абсолютно жесткими стержнями.
Размещая связи-стержни между парами соседних вершин, можно распространить алгоритм моделирования ткани на моделирование растений.
Код и уравнения, рассмотренные в этом разделе предполагают, что все частицы имеют одинаковую массу. Тем же способом можно моделировать и частицы с разными массами, но полученные уравнения будут немного сложнее.
Так, соблюдение ограничения \eqref{eq:C2} для частиц с разными массами реализует следующий псевдокод:
// Псевдокод для соблюдения ограничений (2) delta = x2 - x1 ; deltalength = sqrt (delta * delta ); diff = (deltalength - restlength ) / (deltalength * (invmass1 + invmass2 )); x1 -= invmass1 * delta * diff ; x2 += invmass2 * delta * diff ;
Здесь invmass1 и invmass2 хранят обратные массы частиц \(\vec{x}_1\) и \(\vec{x}_2\) . Если мы хотим, чтобы частица оставаться неподвижной, нужно установить для нее invmass = 0 , что соответствует бесконечной массе. Как и выше, для ускорения расчетов можно использовать приближенное вычисление квадратного корня.
Твердые тела
Уравнения движения твердых тел были предложены задолго до изобретения современных компьютеров. Для того, чтобы в те времена получить какие-то полезные результаты, математики должны были выполнять преобразования формул. Это привело к появлению таких полезных понятий и инструментов, как тензор инерции, момент импульса, момент сил, кватернионы для представления ориентации и т. п. Между тем, имеющиеся сейчас возможности обрабатывать огромные объемы данных в цифровой форме позволяют проводить расчеты для более простых элементов, а в некоторых случаях даже делают такие расчеты более выгодными. В случае трехмерных твердых тел, это означает, что может оказаться удобным моделировать твердое тело с помощью четырех частиц и шести связей (дающих правильное количество степеней свободы: \(4 \cdot 3 - 6 = 6\) ). Это упрощает множество вещей, и именно этим мы займемся далее.
Рассмотрим тетраэдр, в каждую из четырех вершин которого помещена частица. Кроме того, каждое из шести ребер тетраэдра представляет собой ограничение типа стержня, рассмотренное в предыдущем разделе. Этого вполне достаточно, чтобы имитировать твердое тело. Тетраэдр можно поместить внутрь куба, рассмотренного выше, и интегратор Верле обеспечит его правильное движение. Функция SatisfyConstraints() должна позаботиться о двух вещах: 1) чтобы частицы оставались внутри куба, и 2) чтобы были соблюдены шесть ограничений-стержней. Сделать это, как и раньше, можно с использованием релаксации: обычно достаточно 3-х--4-х итераций. Не забывайте также об эффективном вычислении квадратного корня.
Однако ясно, что при столкновениях твердые тела будут вести себя не так как "скелетные" тетраэдры. Существует и другая проблема: до сих пор мы обнаруживали факт столкновения между твердым телом и окружающим "миром" только на основе информации о вершинах: если вершина оказывалась вне куба, она снова проектировалась вовнутрь. Это прекрасно работает, пока внутренняя часть "мира" выпукла. Если же это не так, то тетраэдр сможет проникнуть сквозь границу "мира" даже когда ни одна из его вершин эту границу не пересекала (см. рис.3, где треугольник представляет собой плоский аналог тетраэдра). Рассмотрим, как решается эта проблема.
Сначала разберем более простой вариант задачи. Возьмем стержень, поместим его в кубический "мир" и предположим, что у куба есть небольшой выступ, направленный внутрь. Теперь стержень может пересечь границы "мира", хотя обе частицы на его концах остаются внутри куба (рис.4). Мы не будем вдаваться в тонкости разработки механизма обнаружения столкновений (collision detection), так как это целая отдельная наука. Вместо этого предположим, что подсистема обнаружения столкновений уже существует и делает свое дело: позволяет определить глубину проникновения и координаты точек проникновения для каждого из двух сталкивающихся объектов. Одно из определений точек проникновения и глубины проникновения звучит так: глубина проникновения \(d_p\) - это кратчайшее расстояние, на которое нужно развести два объекта в подходящем направлении, чтобы избежать их столкновения. Точки проникновения - это точки на каждом из объектов, которыми объекты касаются друг друга после того, как упомянутый выше перенос состоялся.
Взгляните еще раз на рис.4. Здесь, после этапа численного интегрирования, стержень проник через границу. Детектор столкновений определил две точки проникновения: \(\vec{p}\) и \(\vec{q}\) . На рис.4а, точка \(\vec{p}\) фактически совпадает с одной из концевых частиц: \(\vec{p}=\vec{x}_1\) . На рис.4б, \(\vec{p}\) лежит между \(\vec{x}_1\) и \(\vec{x}_2\) на расстоянии 1/4 длины стержня от \(\vec{x}_1\) . В обоих случаях, точка \(\vec{p}\) лежит на стержне и, следовательно, ее координаты могут быть выражены в виде линейной комбинации координат точек \(\vec{x}_1\) и \(\vec{x}_2\) : \(\vec{p} = c_1\vec{x}_1 + c_2\vec{x}_2\) такой, что \(c_1 + c_2 = 1\) . В первом случае \(c_1 = 1\) , \(c_2 = 0\) , а во втором - \(c_1 = 0.75\) и \(c_2 = 0.25\) . Эти значения говорят нам, на какое расстояние нужно передвинуть соответствующие частицы.
Чтобы скорректировать положение стержня, переместим его так, чтобы точка \(\vec{p}\) совпала с \(\vec{q}\) . Для этого передвинем частицы \(\vec{x}_1\) и \(\vec{x}_2\) в направлении, заданном вектором, соединяющим \(\vec{p}\) и \(\vec{q}\) : .
В первом случае (рис.4а), мы просто cпроектируем \(\vec{x}_1\) за пределы области, где ей "запрещено" находиться, также как это делали раньше (в направлении \(\vec{q}\) ). Этого будет достаточно, а координаты \(\vec{x}_2\) вообще не нужно изменять. Во втором случае (рис.4б) точку \(\vec{x}_1\) также нужно перенести на большее расстояние, чем \(\vec{x}_2\) , поскольку точка \(\vec{p}\) расположена ближе к \(\vec{x}_1\) (действительно, так как \(\vec{p} = 0.75\vec{x}_1 + 0.25\vec{x}_2\) , то всякий раз перемещая \(\vec{x}_1\) на величину 0.75, мы перемещаем \(\vec{x}_2\) только на 0.25). Другими словами, новые положения частиц \(\vec{x}_1^\prime\) и \(\vec{x}_2^\prime\) задаются соотношениями
\begin{equation} \begin{aligned} \vec{x}_1^\prime &= \vec{x}_1 + 0.75\lambda\cdot\vec{\Delta}, \\ \vec{x}_2^\prime &= \vec{x}_2 + 0.25\lambda\cdot\vec{\Delta}, \end{aligned} \label{eq:x_new} \end{equation}
где \(\lambda\) - неизвестная величина. Новое положение частицы \(\vec{p}\) - \(\vec{p}^\prime\) - вычисляется по формуле
$$ \vec{p}^\prime = c_1\vec{x}_1^\prime + c_2\vec{x}_2^\prime . $$
Вспомним, что мы хотим добиться, чтобы \(\vec{p}^\prime = \vec{q}\) , т. е. должны выбрать \(\lambda\) в точности таким, чтобы \(\vec{p}^\prime\) в результате совпало с \(\vec{q}\) . Так как мы перемещаем частицы только в направлении \(\vec{\Delta}\) , то \(\vec{p}\) также перемещается в направлении \(\vec{\Delta}\) и, следовательно, решение уравнения \(\vec{p}^\prime = \vec{q}\) можно найти, выразив \(\lambda\) из
\begin{equation} \vec{p}^\prime\cdot\vec{\Delta} = \vec{q}\cdot\vec{\Delta} . \label{eq:pq} \end{equation}
Расписывая выражение, стоящее в левой части равенства, получим
$$ \begin{aligned} \vec{p}^\prime\cdot\vec{\Delta} &= (0.75\vec{x}_1^\prime + 0.25\vec{x}_2^\prime) \cdot\vec{\Delta} \\ &= (0.75 (\vec{x}_1 + 0.75\lambda\cdot\vec{\Delta}) + 0.25 (\vec{x}_2 + 0.25\lambda\cdot\vec{\Delta})) \cdot\vec{\Delta} \\ &= ((0.75\vec{x}_1 + 0.25\vec{x}_2)\cdot\vec{\Delta} + \lambda(0.75^2 + 0.25^2)\cdot\Delta^2 \\ &= \vec{p} \cdot\vec{\Delta} + \lambda(0.75^2 + 0.25^2)\cdot\Delta^2 , \end{aligned} $$
что, с учетом правой части \eqref{eq:pq}, дает
$$ \lambda = \frac{(\vec{p}-\vec{q}) \cdot\vec{\Delta}}{(0.75^2 + 0.25^2)\cdot\Delta^2} . $$
Подставляя найденное \(\lambda\) в \eqref{eq:x_new}, получим скорректированные положения частиц \(\vec{x}_1\) и \(\vec{x}_2\) , при которых \(\vec{p}^\prime\) совпадет с \(\vec{q}\) .
На рис.5 показано положение, возникшее после перемещения частиц. Взаимного проникновения объектов теперь нет, но зато нарушено требование неизменности длины стержня. Чтобы это исправить, сделаем еще одну итерацию цикла релаксации (или даже несколько), после чего завершаем цикл исправлений положений частиц.
В случае тетраэдра описанная выше стратегия будет работать аналогично. Сначала находятся точки взаимопроникновения \(\vec{p}\) и \(\vec{q}\) (они также могут находится внутри треугольника) и \(\vec{p}\) представляется линейной комбинацией четырех частиц \(\vec{p}=c_1\vec{x}_1+c_2\vec{x}_2+c_3\vec{x}_3+c_4\vec{x}_4\) таких, что \(c_1+c_2+c_3+c_4=1\) (это потребует решения небольшой системы линейных уравнений). После того, как будет найден \(\vec{\Delta} = \vec{q}-\vec{p}\) , можно будет найти значение \(\lambda\) по формуле
$$ \lambda = \frac{(\vec{p}-\vec{q}) \cdot\vec{\Delta}}{(c_1^2 + c_2^2 + c_3^2 + c_4^2)\cdot\Delta^2} , $$
а исправленные положения частиц определяются как
$$ \begin{align*} \vec{x}_1^\prime &= \vec{x}_1 + c_1\lambda\cdot\vec{\Delta}, \\ \vec{x}_2^\prime &= \vec{x}_2 + c_2\lambda\cdot\vec{\Delta}, \\ \vec{x}_3^\prime &= \vec{x}_3 + c_3\lambda\cdot\vec{\Delta}, \\ \vec{x}_4^\prime &= \vec{x}_4 + c_4\lambda\cdot\vec{\Delta}. \end{align*} $$
Итак, мы рассмотрели столкновение одного твердого тела с неподвижным "миром". Описанный выше метод можно легко обобщить для обработки столкновений нескольких твердых тел. При этом столкновения обрабатываются для одной пары тел в один момент времени и, вместо того, чтобы перемещать только \(\vec{p}\) , понадобится перемещать \(\vec{p}\) и \(\vec{q}\) по направлению друг к другу.
И вновь, после корректировки положений частиц во избежание взаимного проникновения тел, необходимо позаботиться о выполнении еще шести ограничений - неизменности расстояний между частицами, составляющими твердое тело. С помощью этого метода, тетраэдр можно даже вложить внутрь другого объекта, который удобнее использовать вместо самого тетраэдра при обработке столкновений. На рис.6 показан тетраэдр, помещенный в куб.
Во-первых, куб должен быть каким-то образом прикреплен к тетраэдру. Один из подходов состоит в том, что в качестве центра куба выбирается центр масс тетраэдра \(0.25\cdot (\vec{x}_1 + \vec{x}_2 + \vec{x}_3 + \vec{x}_4)\) , а затем по текущим координатам тетраэдра вычислить координаты вершин куба. При обнаружении столкновения, точка контакта \(\vec{p}\) (которая теперь расположена на кубе) обрабатывается также, как и выше. Аналогично вычисляются и обновленные значения координат частиц. Для ускорения расчетов можно заранее вычислить коэффициенты \(c_1\) --\(c_4\) для всех вершин куба. Если \(\vec{p}\) окажется вершиной, то значения \(c_1\) --\(c_4\) могут быть найдены и использованы непосредственно. В противном случае, \(\vec{p}\) лежит внутри треугольника или на одной из его сторон, и значения \(c_1\) --\(c_4\) можно получить из предварительно вычисленных значений вершин треугольника при помощи интерполяции.
Как правило, для обработки столкновений достаточно 3--4 итераций. Если релаксацию остановить слишком рано, то тела не будут вести себя как абсолютно твердые. Но это даже хорошо, ведь абсолютно твердых тел в природе не существует. Кроме того, это делает систему более устойчивой.
При перестановке положений частиц, составляющих тетраэдр, физические свойства тела должны быть изменены соответствующим образом (математически это означает, что изменяется тензор инерции тела, зависящий от положений и масс частиц).
По тому же принципу, что и тетраэдр, можно задать другую подобную конфигурацию частиц и связей, расположив частицы в точках с координатами \((0,0,0)\) , \((1,0,0)\) , \((0,1,0)\) и \((0,0,1)\) . Пусть \(\vec{a}\) , \(\vec{b}\) , и \(\vec{c}\) - векторы, направленные из частицы 1 к частицами 2, 3 и 4, соответственно. Ограничим положения частиц, требованием, чтобы векторы \(\vec{a}\) , \(\vec{b}\) , и \(\vec{c}\) имели единичную длину, и угол между каждой из трех пар векторов был равен \(90^\circ\) (соответствующие скалярные произведения должны быть равны нулю). Заметим, что это снова, как и в случае с тетраэдром, даст 4 частицы и 6 связей.
Сочлененные тела
Теперь мы можем соединять несколько твердых тел при помощи шарниров (цилиндрических, сферических и т. п.). Предположив, что два тела имеют одну общую частицу, мы получим сферический шарнир (pin joint), а если общими являются две частицы, то получим цилиндрический шарнир (hinge) (рис.7). Таким же образом можно связать два тела при помощью стержня или любого другого вида связи - надо только не забыть добавить код для обработки нового вида связи в цикл релаксации.
Такой подход позволяет построить полную модель сочлененного человеческого тела. Реалистичность увеличится, если дополнительно реализовать ограничения на угловые перемещения в шарнирах. Существует несколько способов реализации таких ограничений. Простейший способ предполагает использование ограничения типа стержня, которые срабатывает только тогда, когда расстояние между двумя частицами станет ниже некоторого порогового значения (в данном случае мы имеем дело с односторонней связь вида \(|\vec{x}_2 - \vec{x}_1| > 100\) ). Как следствие этого, обе частицы никогда не смогут слишком приблизиться друг к другу (рис.8).
Другой метод создания ограничений на угловые перемещения требует соблюдения следующего условия для скалярного произведения
$$ (\vec{x}_2 - \vec{x}_0)\cdot (\vec{x}_1 - \vec{x}_0) < \alpha . $$
Можно также ограничить перемещение частицы определенной плоскостью. И вновь, положения частиц, не удовлетворяющие заданным ограничениям, должны быть скорректированы. Делается это аналогично случаю стержня, хотя соответствующие формулы будут немного сложнее.
В транспортном потоке каждый автомобиль движется либо под влиянием со стороны других участников движения, либо в отсутствии их влияния. Движение автомобиля будем называть свободным, если ни один из участников дорожного движения не оказывает влияния на движение этого автомобиля, а также на мнение водителя о дорожно-транспортной обстановке, в результате изменения которого он мог бы изменить режим движения своего автомобиля. Скорость движения такого автомобиля будем называть скоростью свободного движения на данном участке дороги.
Основанием для моделирования свободного движения автомобиля явились:
Уравнения теории эксплуатационных свойств автомобиля:
а) тягово-скоростные и тормозные свойства автомобиля;
б) уравнения криволинейного движения и устойчивости автомобиля.
Натурные наблюдения за параметрами движения автомобилей на двухполосных дорогах.
Математическая модель свободного движения автомобиля имеет следующую концептуальную основу:
Мера воздействия на органы управления автомобиля, а также мнение водителя о ДТС могут измениться только при наступлении одного из ситуаций, перечисленных ниже.
На каждом участке дороги водитель стремится поддерживать оптимальную с его точки зрения (базовую) скорость движения, которая зависит от цели и дальности поездки, вида перевозимого груза (количества пассажиров), состояния здоровья и степени утомления водителя и других факторов. Базовая скорость в модели задается случайным законом распределения, полученным в результате натурных наблюдений.
Если базовая скорость движения автомобиля на следующем участке дороги отличается от базовой скорости на текущем участке, то водитель заранее изменяет скорость движения автомобиля таким образом, чтобы к моменту въезда на новый участок скорость движения достигла величины базовой скорости на новом участке.
Водитель может через органы управления автомобиля воздействовать на параметры движения следующими способами:
а) изменить скорость движения и ускорение нажатием на педаль тормоза или акселератора (выдвижением рейки);
б) изменить передаточное число КПП, что позволяет изменить диапазон значений скорости движения автомобиля;
в) изменить направление движения автомобиля, вращением рулевого колеса.
Кроме перечисленных действий водитель может включить стоп-сигналы (нажатием педали тормоза) или сигналы поворота, что может служить причиной изменения режима движения других автомобилей.
С точки зрения обеспечения базовой скорости движения автомобиля в конкретных дорожных условиях могут возникнуть следующие характерные обстоятельства:
возможность водителя увеличить скорость движения автомобиля до базовой скорости ограничена тягово-динамическими характеристиками автомобиля;
возможность водителя уменьшить скорость движения автомобиля в режиме торможения (экстренное торможение) ограничена коэффициентом сцепления шины с дорогой и/или тормозными характеристиками автомобиля;
возможность водителя изменить скорость движения автомобиля до базовой скорости не ограничено ни тягово-динамическими или тормозными характеристиками автомобиля, ни сцепными качествами поверхности дороги.
Рассмотрим подробнее, каким образом моделируется движение автомобиля в перечисленных выше случаях.
В первом случае движение автомобиля моделируется на основе известных в теории автомобиля дифференциальных уравнений, полученными на основе уравнения силового баланса автомобиля:
P т = P п + P к + P в + P и, (2.5)
где P т - тяговая сила при установившейся скорости автомобиля;
P п - сила сопротивления подъему;
P к - сила сопротивления качению;
P в - сила сопротивления воздуху;
P и - сила сопротивления разгону (приведенная сила инерции).
Существуют различные зависимости, аппроксимирующие внешние характеристики двигателя. В рассматриваемой модели дифференциальные уравнения движения автомобиля получены на основе аппроксимации внешней характеристики двигателя, приведенной в работе :
где N e , N max - соответственно, мощность и максимальная мощность двигателя, квт;
M k - крутящий момент двигателя, Нм;
M kN - крутящий момент двигателя, при максимальной мощности, Нм;
a, b, c - постоянные коэффициенты для данного двигателя;
n - угловая скорость коленчатого вала двигателя, об/мин;
n N - угловая скорость коленчатого вала при максимальной мощности двигателя, об/мин.
После замены всех членов уравнения (2.5) на соответствующие значения и некоторых преобразований получается:
Где m a - масса автомобиля, кг;
m 0 - масса автомобиля, с номинальной нагрузкой, кг;
u k i - передаточное число коробки передач;
v - скорость движения автомобиля, м/с;
тр - коэффициент полезного действия трансмиссии;
k p - коэффициент коррекции двигателя;
Номинальная скорость движения автомобиля при i-й передаче, м/c;
G a - сила тяжести, действующий на автомобиль, Н;
k f - параметр учета влияния скорости движения на коэффициент сопротивления качению колеса;
W - фактор обтекаемости автомобиля, кг/м;
f 0 - коэффициент сопротивления качению при малой скорости движения;
б - продольный уклон дороги.
Уравнение (2.8) определяет ускорение автомобиля в зависимости от скорости движения. Для рассматриваемой имитационной модели зависимости вида “ускорение - скорость”, “путь разгона - скорость” и т. д. непригодны, поскольку при пересчете векторов-координат автомобилей через интервал времени t min (см. блок 12 на рис 2.16) возникает необходимость определения этих параметров в зависимости от времени.
С целью определения зависимости скорости движения от времени при полной подаче топлива можно проинтегрировать выражение (2.8). Пусть начальным условием будет v = v 0 при t=0. Тогда после интегрирования получим:
Проинтегрируем (2.13) еще раз, при начальных условиях t=0 и s=s 0 . Получим:
где v 0 - начальная скорость движения автомобиля;
s 0 - начальное положение автомобиля;
v 1 и v 2 - корни уравнения.
Для того чтобы получить зависимость a = a(t) нужно найти производную выражения (2. 13) по времени. Получим:
Выражения (2.13) - (2.15) позволяют пересчитать параметры движения автомобиля через произвольный промежуток времени t min , в условиях ограничения параметров движения тягово-динамическими характеристиками автомобиля.
Во втором случае моделирование движения автомобиля осуществляется при следующих допущениях:
силы реакции R x достигают максимального значения одновременно для всех колес;
коэффициенты сцепления x всех колес с дорогой, а следовательно, и ускорение автомобиля j з остаются неизменными за весь период установившегося замедления.
При таких допущениях процесс торможения может быть описан тормозной диаграммой j з = j(t) (рис. 2.3) . Весь процесс торможения с момента обнаружения опасности до полной остановки автомобиля состоит из следующих этапов:
время реакции водителя t рв;
время запаздывания t з;
время нарастания тормозного усилия t н;
время установившегося замедления t уст;
время растормаживания t р.
В случае торможения при полном использовании сил сцепления (экстренное торможение) j уст зависит только от коэффициента сцепления шин с дорогой и продольного уклона дороги, а значение ускорения можно считать постоянным:
Скорость движения и пройденный путь автомобиля в произвольный момент времени t легко определяются интегрированием выражений (2.16) и (2.17):
Рис. 2.3 Тормозная диаграмма автомобиля
И, наконец, в третьем случае водитель имеет возможность придавать автомобилю тот режим движения, который по его мнению является наиболее безопасным и целесообразным в сложившейся дорожно-транспортной ситуации. В этом случае значение ускорения находится в интервале
j уст < a < a max , (2. 20)
и определяется водителем.
С целью упрощения расчетов предполагается, что в таком режиме автомобиль движется равноускоренно до следующего особого состояния. Скорость движения и пройденный путь автомобиля через промежуток времени t в таком случае определяются следующим образом:
Теперь рассмотрим движение автомобиля в условиях изменения направления движения.
В модели предусмотрены следующие виды траекторий движения автомобилей:
прямолинейное движение;
круговое движение (угол поворота управляемых колес не изменяется);
криволинейное движение при постоянной угловой скорости поворота управляемых колес;
Маневрирование автомобилей (обгоны, перестраивания по полосам движения, повороты и т.д.) являются опасными, но в то же время неотъемлемыми элементами движения автомобилей. С изменением направления движения автомобилей связаны многочисленные конфликтные ситуации и ДТП, поскольку такие маневры, часто совершаемые внезапно для других участников дорожного движения, создают возмущения в транспортном потоке. Несмотря на это, в существующих имитационных моделях дорожного движения различных отечественных и зарубежных исследователей упрощенно описываются движение автомобиля и траектории отдельных его точек при изменении направления движения. Часто положение автомобиля на дороге задается только продольной координатой; при этом подразумевается, что в поперечном направлении автомобиль не движется, а изменение полосы движения при обгонах или маневрах перестроения происходит скачкообразно. В тех моделях, где учитывается поперечное движение, в лучшем случае движение автомобиля рассматривается как плоскопараллельное движение, при котором продольная ось автомобиля остается параллельной относительно продольной оси дороги. Такое упрощение оправдано при решении некоторых задач, таких, как определение скорости сообщения, определение пропускной способности участка дороги, при решении проблем экологии и т. д., потому что значительно упрощает модель и уменьшает объем расчетов. Однако при решении задач оценки уровня безопасности движения такое упрощение не оправдано.
Криволинейное движение автомобиля в рассматриваемой имитационной модели определяется:
координатами расчетной точки автомобиля относительно неподвижной системе координат;
курсовым углом движения автомобиля;
углом поворота управляемых колес автомобиля;
угловой скоростью поворота управляемых колес.
Расчетной точкой, относительно которого производятся все расчеты по определению координат автомобиля, в модели принята середина задней оси автомобиля. В таком случае уравнения, определяющие координаты имеют наименее громоздкий вид.
В модели движение автомобиля рассматривается как чередование прямолинейного движения, кругового движения и криволинейного движения переменного радиуса. Первые два достаточно хорошо описываются сравнительно простыми аналитическими выражениями. Рассмотрим подробней криволинейное движение (рис. 2.4). В дальнейшем будем принимать следующие допущения:
углы поворота обоих управляемых колес автомобиля равны между собой;
у колес автомобиля отсутствует боковой увод;
угловая скорость поворота управляемых колес постоянна;
расчетная точка автомобиля движется с постоянным ускорением;
автомобиль движется на плоскости (плоское движение);
отсутствует буксование колес;
крен автомобиля не влияет на траекторию.
Для оценки уровня безопасности движения первые два допущения мало значимы. Если в процессе движения значение угловой скорости поворота управляемых колес либо ускорение расчетной точки автомобиля существенно изменяется, то траектория автомобиля разбивается
на несколько участков, на каждом из которых значения указанных величин принимается постоянными.
Рис. 2.4
Таким образом, в пределах приемлемой точности вышеуказанные допущения значительно упрощают расчеты.
В работах приводятся уравнения, которые позволяют определить траекторию движения автомобиля. Так, в работе приведены следующие формулы, определяющие курсовой угол и координаты движения середины задней оси автомобиля:
где x в, y в - координаты середины задней оси автомобиля;
C 1 и C 2 - постоянные, определяемые начальными условиями;
Угол поворота управляемых колес автомобиля;
v o - скорость движения автомобиля;
k - угловая скорость поворота управляемых колес автомобиля;
L 0 - база автомобиля;
k p - режимный параметр, характеризующий режим криволинейного движения:
Аналогичные формулы, но относительно движения центра масс автомобиля приведены в работе:
где x ц.м. , y ц.м. - координаты центра масс автомобиля;
C 1 , C 2 , C 3 - постоянные, определяемые начальными условиями;
v a - скорость движения автомобиля;
v y - скорость бокового смещения центра масс автомобиля;
a - угловая скорость продольной оси автомобиля в горизонтальной плоскости.
Траектории точек автомобиля, находящихся на продольной оси, могут быть заданы зависимостью ее кривизны от времени. Однако, моделирование криволинейного движения автомобиля в неподвижной системе координат через x в, и y в, а также угла гораздо удобнее.
В уравнениях (2.23) - (2.26) скорость движения автомобиля предполагается постоянной величиной. В реальном процессе движения автомобиля в транспортном потоке часто сочетаются криволинейное движение с изменением скорости движения. Это происходит, в частности, в следующих ситуациях:
перед остановкой автомобиль одновременно уменьшает скорость движения и криволинейным движением приближается к правой кромке проезжей части;
после трогания с места автомобиль увеличивает скорость движения и по криволинейной траектории приближается к середине полосы движения;
во время обгона, особенно в первой фазе обгона с ожиданием, автомобиль меняет полосу движения и одновременно увеличивает скорость движения;
на перекрестках, после запрещающего сигнала светофора, или после пропуска помехи справа, поворачивающие автомобили двигаются криволинейно и с ускорением и т.д.
Сочетание криволинейного движения с изменением скорости движения является потенциальным источником конфликтных ситуаций, т.к. оно часто заставляет других участников движения изменить режим движения. Криволинейное движение автомобиля обязательно необходимо моделировать с учетом ускорения.
Уравнения, приведенные в разных исследованиях для определения траектории криволинейного движения, чаще всего выражают параметры движения автомобиля в виде дифференциальных уравнений или интегралов. Как правило, подынтегральные выражения имеют сложный аналитический вид и в явном виде не интегрируются. Поэтому, либо такие уравнения решаются численными, приближенными методами с применением компьютеров, либо траектория движения строится графическим или графоаналитическим методом. При решении многих конкретных задач эти методы оказываются целесообразными, особенно первый из них, при помощи которого можно рассчитать траекторию движения практически с любой точностью. Однако это достигается за счет больших затрат машинного времени счета. В течение имитационных экспериментов на задаваемом участке дороги одновременно могут находиться сотни автомобилей, каждый из которых совершает различные маневры, в том числе и криволинейное движение с ускорением. При использовании итерационных методов с целью расчета траектории движения большие затраты машинного времени неизбежны. Это резко замедляет работу имитационной модели в реальном масштабе времени. Поэтому здесь более целесообразным представляется разложение подынтегральных выражений в степенной ряд с требуемой точностью. Далее, после интегрирования полученных многочленов можно получить искомые параметры криволинейного движения в аналитическом виде.
Пусть автомобиль движется с постоянным ускорением a. Тогда выражение (2.23) принимает следующий вид:
Разложим в степенной ряд подынтегральное выражение, содержащееся в правой части (2.30) и возьмем столько членов, сколько требуется для обеспечения необходимой точности. Обозначим
Тогда (2.30) можно переписать в следующем виде:
Приближенные аналитические решения уравнений (2.24) и (2.25) также будем искать, разлагая в степенной ряд функции f()=cos(()) и g()=sin(()), где курсовой угол определяется из выражения (2.31). Производные от упомянутых функций первых нескольких порядков приведены в табл. 2.1.
Пусть в начале криволинейного движения курсовой угол автомобиля 0 =0. Тогда:
Теперь можно определить траекторию движения точки B автомобиля (рис. 2.4). Выражения (2.24) и (2.32) позволяют после интегрирования (2.24) получить
Аналогично, из выражений (2.25) и (2.33) после интегрирования (2.25) можно получить
Выражения (2.34) и (2.35) вполне приемлемы для расчета траектории движения автомобиля, если за время криволинейного движения курсовой угол автомобиля изменяется не более чем на 90°. В реальных дорожных условиях такое изменение практически никогда не превышает
Таблица 2.1
Производные от функций f(и) и g(и) первых нескольких порядков
|
Порядок производной n |
f(n)(и), при и=0 |
g(n)(и), при и=0 |
|
105 B4-588 B2-896 D2 |
420 B3-272 B+1120 B D2 |
|
|
2520 B3 D+8160 B D |
8064 B2 D-2240 D3+2176 D |
|
|
6300 B4-18960 B2+25200 B2 D2-31104 D2 |
945 B5+16380 B3-7936 B+57120 B D2 |
Таблица 2.2
Значения функций f(и) и g(и), описывающие процесс поворота автомобиля ГАЗ-24 при v 0 =10 м/с, k =0,05 рад/с, L 0 =2,8 м, 45°.
|
г(и), по формуле (2.23), град |
f(и), по формуле (2.23) |
f(и), по формуле (2.32) |
g(и), по формуле (2.23) |
g(и), по формуле (2.33) |
|
В табл. 2.2 приведены значения функций f() и g(), которые описывают процесс поворота автомобиля ГАЗ-24 при v 0 =10 м/с, k =0,05 рад/с, L 0 =2,8 м. Вычисления проводились с одной стороны с применением выражения (2.23), а с другой стороны - выражений (2.32) и (2.33). При <90° расхождения в значениях функций не превышает 0,1%, что вполне обеспечивает требуемую точность вычислений.
Координаты любой другой точки E автомобиля (рис. 2.4), находящейся от точки B на расстояниях a и b, соответственно по продольной и поперечной оси автомобиля, можно определить следующим образом:
Рассмотрим, каким образом происходит процесс поворота в имитационной модели в режиме свободного движения. До начала криволинейного движения водитель уменьшает скорость движения до значения базовой скорости свободного движения v св на участке поворота. С такой скоростью он завершает поворот, а после поворота при необходимости опять увеличивает скорость движения. Значение скорости движения v св определяется на основе натурных наблюдений.
Модель процесса поворота представляет собой последовательность трех этапов (рис. 2.5):
водитель с постоянной угловой скоростью р1 вращает рулевое колесо по/против направления часовой стрелки, при правом/левом повороте. Управляемые колеса совершают поворот с угловой скоростью k1 . Автомобиль совершает криволинейное движение. Радиус кривизны траектории движения уменьшается от + до R п (участок E 1);
водитель держит рулевое колесо в неподвижном состоянии. Автомобиль двигается по круговой траектории с постоянным радиусом поворота R п относительно центра поворота С. Этот этап в процессе поворота может отсутствовать (участок E 2);
водитель вращает рулевое колесо в обратном направлении с постоянной угловой скоростью р2 . Управляемые колеса поворачиваются с угловой скоростью k2 . Автомобиль снова движется по криволинейной траектории. Радиус кривизны траектории движения увеличивается от R п до + (участок E 3).
Первая часть неравенства (2.38) объясняется стремлением водителя удержать автомобиль на своей полосе движения. При < S 0min водитель не в состоянии предотвратить выезд автомобиля за пределы своей полосы движения, если ускорение не изменять (если ускорение изменится, то S 0min тоже изменится).
Рис. 2.5 Схема к описанию модели процесса поворота автомобиля
Вторая часть неравенства объясняется тем, что, если в начале поворота > S 0max , и в дальнейшем водитель не изменит параметры движения автомобиля (k и a), то автомобиль либо выедет за пределы проезжей части дороги, либо пересечет биссектрису OC угла ц под острым углом, т.е. в дальнейшем опять-таки выедет за пределы полосы движения. Поэтому, для предотвращения этого, водитель вынужден приближаться к повороту до тех пор, пока не выполнится второе условие неравенства (2.38) и только после этого начинать поворот.
Водитель начинает поворот на некотором расстоянии до центра O перекрестка. является случайной величиной. При принятых выше ограничениях значение величины находится в пределах:
S 0min < < S 0max (2.38)
Угловая скорость поворота управляемых колес k , которая определяется как отношение угловой скорости поворота рулевого колеса р на передаточное число рулевого механизма, рассматривается как случайная величина и определяется следующим образом:
определяется минимальное значение kmin угловой скорости поворота управляемых колес, при котором автомобиль в процессе поворота не выедет за пределы полосы движения с противоположной стороны от центра поворота С (рис. 2.6а);
определяется максимальное значение kmax угловой скорости поворота управляемых колес, при котором автомобиль в процессе поворота не выедет за пределы полосы движения со стороны центра поворота (рис. 2.6б);
определяются значения р min и р max угловой скорости поворота рулевого колеса, соответствующие значениям kmin и kmax ;
между значениями р min и р max генератором случайных чисел по заранее заданному закону распределения разыгрывается случайное число k ", которое и является значением угловой скорости поворота управляемых колес на первом этапе поворота.
Рис. 2.6 Схема к вычислению угловой скорости поворота управляемых колес автомобиля в процессе поворота
Закон распределения случайной величины k ", определяется в результате натурных наблюдений (см. п. 3.2, рис. 3.12). На третьем этапе значение угловой скорости поворота управляемых колес k "" определяется аналогично.