2. Механизм формирования математической модели
В этом разделе разберём, как PRADIS переходит от расчетной схемы к математической модели, которую затем решает solver.
Продолжим рассматривать систему из простого примера: масса соединена с основанием через пружину и демпфер, а к телу приложена внешняя сила.
2.1. Уравнение равновесия в дискретный момент времени
Для системы можно записать уравнение равновесия сил:
(2.1)
При численном интегрировании решение ищется не непрерывно, а в отдельных точках времени. Поэтому для каждого i-го момента времени уравнение записывается так:
(2.2)
где:
Здесь xi, vi, ai — перемещение, скорость и ускорение тела в момент времени ti.
2.2. Связь перемещения, скорости и ускорения
Для i-го момента времени значения xi, vi, ai связаны формулами метода интегрирования:
Важно: в PRADIS при записи равновесия суммируются усилия, действующие со стороны системы на элементы. Поэтому знаки в выражениях для сил могут отличаться от тех, которые используются при ручном выводе через выбранное положительное направление.
После подстановки выражений для сил и формул интегрирования дифференциальная задача превращается в нелинейное алгебраическое уравнение для текущего шага времени.
2.3. Функция невязки
Для решения вводится функция:
(2.5)
где:
Переменной z можно выбрать одну из связанных величин:
В данном примере выбираем:
Тогда значение f(z) — это невязка равновесия. Иными словами, это сумма сил в узле при текущем приближении решения. Метод Ньютона должен уменьшить эту невязку до допустимого значения.
2.4. Производная функции невязки
Для метода Ньютона требуется не только значение f(z), но и производная:
Так как функция складывается из сил, производная записывается как сумма производных:
Каждая сила может зависеть от перемещения, скорости и ускорения. Поэтому производная каждой силы записывается по правилу сложной функции:
Так как в примере выбрано:
то далее используются производные по 
Аналогично:
2.5. Зависимости xi и ai от vi
Из формул метода интегрирования получаем:
(2.8)
Продифференцируем выражения по 
:
Эти коэффициенты используются при приведении якобианов элементов к выбранной переменной решения.
2.6. Частные производные сил
Теперь вычислим частные производные для каждой силы.
Внешняя сила
Внешняя сила 
зависит только от времени, поэтому:
Сила упругости
Сила упругости 
зависит от перемещения:
Сила вязкого сопротивления
Сила вязкого сопротивления 
зависит от скорости:
Сила инерции
Сила инерции 
зависит от ускорения:
2.7. Производная для метода Ньютона
Подставим частные производные и коэффициенты связи:
Получаем:
Суммарно:
(2.15)
Эта формула совпадает по смыслу с результатом, полученным при ручном выводе. Но теперь видно главное: PRADIS не обязан вручную составлять полное дифференциальное уравнение движения. Достаточно знать, как элементы соединены, какие силы они создают и какие производные имеют.
2.8. Какая информация нужна PRADIS для формирования модели
Для автоматического формирования математической модели PRADIS использует следующий набор данных.
1. Сведения о стыковке элементов
Нужно знать, какие элементы соединены между собой и через какие узлы.
Рис.2.1. Схема стыковки элементов: пружина, демпфер, масса, воздействие, узел стыковки
2. Условие равновесия
Для каждого узла записывается условие равновесия потоковых переменных. В механической системе это равновесие сил:
(2.16)
где 4 - количество сил, сходящихся в узле стыковки (количество стыкующихся ветвей элементов).
Такое уравнение называют топологическим, потому что оно определяется структурой связей в схеме.
3. Компонентные уравнения элементов
Для каждого элемента должны быть известны зависимости, по которым вычисляются усилия:
Эти уравнения описывают физику конкретного элемента.
4. Частные производные усилий
Для метода Ньютона нужны производные усилий по:
перемещению,скорости,ускорению.
Эти производные формируют якобиан элемента.
5. Формулы метода интегрирования
Используются алгебраические связи между x, v, a на текущем шаге:
6. Выбор переменной решения
Solver должен знать, относительно какой величины решается нелинейное уравнение:
x, v или a
В примере используется:
2.9. Представление системы через элементы PRADIS
Система в PRADIS собирается как совокупность элементов, соединённых по общим степеням свободы.
Каждый элемент имеет определённое число узлов и степеней свободы.
Рис.2.2. Модели элементов: пружина, демпфер, точечная масса, сосредоточенная сила
Пользователь собирает модель системы из готовых моделей элементов. Его задача — корректно выбрать элементы, задать параметры и соединить их между собой. Формирование математической модели выполняет PRADIS.
Пример описания структуры:
$ FRAGMENT: Пример
# BASE: 1
# STRUCTURE:
Пружина' K (1 2; Коэффициент жесткости)
Нелинейный демпфер ' MUNL (1 2; Коэффициент вязкости)
Масса ' M (2; Масса тела)
Воздействие ' FSIN (2 1; Q, T, начальная фаза)
В этом описании:
K — модель пружины;
MUNL — модель нелинейного демпфера;
M — модель массы;
FSIN — модель синусоидального воздействия.
Узел 2 связан с массой, пружиной, демпфером и внешним воздействием. Узел 1 закреплён, то есть свободные концы пружины и демпфера неподвижны.
2.10. Что делает программа при обработке структуры
При обработке структуры PRADIS определяет размерность системы уравнений.
В данном примере есть два узла, но один из них закреплён. Поэтому на этапе расчёта уравнение для закреплённого узла исключается, а его кинематические характеристики принимаются равными нулю:
x = 0
v = 0
a = 0
Далее расчёт представляет собой последовательность шагов по времени. На каждом шаге решается нелинейное уравнение равновесия.
2.11. Роль моделей элементов
Для любой модели элемента входными данными являются:
параметры элемента;
текущие перемещения узлов;
текущие скорости узлов;
текущие ускорения узлов.
По этим данным модель элемента должна вычислить:
1. Вектор усилий
Это усилия, действующие со стороны системы на элемент.
2. Якобиан элемента
Это частные производные усилий по перемещениям, скоростям и ускорениям узлов.
Если элемент имеет N степеней свободы, то:
длина вектора усилий = N
а размер якобиана:
N × N × 3
2.12. Пример элемента: пружина
Для двухузловой одномерной идеально упругой пружины:
Рис.2.3. Двухузловая модель пружины
Усилия:
Производные по перемещениям:
Производные по скоростям:
(2.19)
Производные по ускорениям:
(2.20)
Так как узел 1 закреплён, в расчёте используется часть, относящаяся ко второму узлу:
(2.21)
Аналогично учитывается вклад остальных элементов: демпфера, массы и внешнего воздействия.
2.13. Третий шаг интегрирования
Продолжим расчёт и выполним третий шаг по времени.
После второго шага были получены:
Рекомендуемый третий шаг с учетом параметров:
Рис.2.4. Общий алгоритм перехода между временными шагами
Рис.2.4б. Алгоритм выполнения одного временного шага
Рис.2.4в. Алгоритм одной итерации метода Ньютона
2.13.1. Коэффициенты приведения якобиана
2.13.2. Начальное приближение
По явному прогнозу:
Также:
2.14. Первая итерация Ньютона
На первой итерации модели элементов вычисляют силы и производные по текущим значениям:
x_{3}^{0}, \quad v_{3}^{0}, \quad a_{3}^{0}
Пружина
F_y = kx = 20000 \cdot 46.5 \times 10^{-5} = 9.3
\frac{\partial F_y}{\partial x} = 20000
\frac{\partial F_y}{\partial v} = 0
\frac{\partial F_y}{\partial a} = 0
Демпфер
F_{\text{б}} = \mu v |v| = 1000 \cdot 0.24081 \cdot |0.24081| = 58.0
\frac{\partial F_{\text{б}}}{\partial x} = 0
\frac{\partial F_{\text{б}}}{\partial v} = 2 \cdot 1000 \cdot |0.24081| = 481.6
\frac{\partial F_{\text{б}}}{\partial a} = 0
Точечная масса
F_{\text{и}} = ma = 0.1 \cdot 59.21 = 5.9
\frac{\partial F_{\text{и}}}{\partial x} = 0
\frac{\partial F_{\text{и}}}{\partial v} = 0
\frac{\partial F_{\text{и}}}{\partial a} = 0.1
Внешняя сила
F_c = -1000 \sin\!\left(\frac{2\pi}{0.2\pi} \cdot (1.438 \times 10^{-3} + 2.63 \times 10^{-3})\right)
F_c = -40.6
\frac{\partial F_c}{\partial x} = \frac{\partial F_c}{\partial v} = \frac{\partial F_c}{\partial a} = 0
Сумма сил
\Sigma F = F_y + F_{\text{б}} + F_{\text{и}} + F_c
\Sigma F = -40.6 + 9.3 + 58.0 + 5.9 = 32.6
Следовательно:
f(z^0) = 32.6
Производная
\frac{df(z)}{dz} = \frac{dF_c}{dz} + \frac{dF_y}{dz} + \frac{dF_{\text{б}}}{dz} + \frac{dF_{\text{и}}}{dz}
\frac{dF_c}{dz} = 0
\frac{dF_y}{dz} = 20000 \cdot 1.31 \times 10^{-3} = 26.2
\frac{dF_{\text{б}}}{dz} = 481.6
\frac{dF_{\text{и}}}{dz} = 0.1 \cdot 380.2 = 38.0
\frac{df(z)}{dz} = 0 + 26.2 + 481.6 + 38.0 = 545.8
Поправка и новое приближение
\Delta z^1 = -\frac{f(z^0)}{f'(z^0)}
\Delta z^1 = -\frac{32.6}{545.8} = -0.05973
z^1 = z^0 + \Delta z^1 = 0.24081 - 0.05973 = 0.18108
Уточняем значения:
v_3^1 = z^1 = 0.18108
a_3^1 = \frac{v_3^1 - v_2}{\Delta t_3} = \frac{0.18108 - 0.08509}{2.63 \times 10^{-3}} = 36.5
x_3^1 = x_2 + \left(\frac{v_2 + v_3^1}{2}\right)\Delta t_3
x_3^1 = 41.1 \times 10^{-5}
Проверяем условия завершения:
|f(z^0)| > \delta_f
|\Delta z^1| > \delta_z
Итерации продолжаются.
2.15. Вторая и третья итерации Ньютона
Для второй итерации получаем:
f(z^1) = -3.6
\Delta z^2 = -0.00858
x_3^2 = 39.8 \times 10^{-5}
v_3^2 = 0.17250
a_3^2 = 33.2
Проверка показывает, что условия завершения ещё не выполнены:
|f(z^1)| > \delta_f
|\Delta z^2| > \delta_z
Переходим к третьей итерации.
На третьей итерации:
f(z^2) = -0.07
|f(z^2)| < \delta_f
\Delta z^3 = -0.00018
|\Delta z^3| < \delta_z
x_3^3 = 39.8 \times 10^{-5}
v_3^3 = 0.17232
a_3^3 = 33.2
Оба критерия завершения выполнены.
2.16. Локальная погрешность третьего шага
После завершения итераций оцениваем локальную погрешность:
l_{p3} = \left|\frac{v_{3p} - v_{3c}}{2}\right|
l_{p3} = \left|\frac{0.24081 - 0.17232}{2}\right| = 0.034
Допустимая локальная погрешность:
\delta_l = 0.001
\text{Так как:}
l_{p3} > \delta_l
шаг слишком большой и расчёт необходимо повторить с меньшим значением \Delta t_3.
2.17. Коррекция шага интегрирования
Обычная формула выбора шага хорошо работает, когда отношение:
\dfrac{\delta_l}{l_{pi}}
близко к единице.
Если это отношение сильно отличается от единицы, шаг может получиться завышенным. Поэтому используется скорректированное правило:
\text{если } \dfrac{\delta_l}{l_{pi}} < 0.25:
\Delta t_{\text{рек}} = c \Delta t_i \left(\dfrac{\delta_l}{l_{pi}}\right)
\text{если } \dfrac{\delta_l}{l_{pi}} > 7:
\Delta t_{\text{рек}} = 4c \Delta t_i
\text{если } 0.25 \le \dfrac{\delta_l}{l_{pi}} \le 7:
\Delta t_{\text{рек}} = c \Delta t_i \sqrt{\dfrac{\delta_l}{l_{pi}}}
В данном случае:
\dfrac{\delta_l}{l_{p3}} = \dfrac{0.001}{0.034} = 0.03
Так как:
0.03 < 0.25
используем:
\Delta t_{\text{рек}} = c \Delta t_3 \left(\dfrac{\delta_l}{l_{p3}}\right)
\Delta t_{\text{рек}} = 0.8 \cdot 2.63 \times 10^{-3} \cdot \left(\dfrac{0.001}{0.034}\right)
\Delta t_{\text{рек}} = 0.061 \times 10^{-3}\ \text{с}
2.18. Повторный расчёт третьего шага
Устанавливаем:
\Delta t_3 = 0.061 \times 10^{-3}\ \text{с}
и повторяем расчёт третьего шага.
После повторного расчёта получаем:
t_3 = 1.499 \times 10^{-3}\ \text{с}
x_3 = 6.64 \times 10^{-5}\ \text{м}
v_3 = 0.08862\ \text{м/с}
a_3 = 58.1\ \text{м/с}^2
Локальная погрешность находится в пределах нормы.
Рекомендуемое значение следующего шага:
\Delta t_{\text{рек}} = 0.264 \times 10^{-3}\ \text{с}
Расчёт третьего шага завершён.
2.19. Основной вывод
Главная идея этого раздела — разделение функций между программой интегрирования и моделями элементов.
Программа интегрирования:
- выполняет шаги по времени;
- организует итерации метода Ньютона;
- собирает уравнения равновесия;
- использует векторы сил и якобианы элементов;
- контролирует локальную погрешность;
- выбирает следующий шаг.
Модели элементов:
- описывают физические свойства конкретных элементов;
- вычисляют усилия;
- вычисляют частные производные;
- передают эти данные вычислительному ядру.
Программа интегрирования не должна знать, по каким физическим законам работает конкретный элемент. Она работает только с равновесием потоковых переменных.
Именно такое разделение обеспечивает универсальность вычислительного ядра PRADIS.
Один и тот же алгоритм может использоваться для объектов разной физической природы:
механических;
гидравлических;
газодинамических;
тепловых;
электрических.
Главное условие — процессы должны подчиняться законам равновесия потоковых переменных:
равновесие сил;
равновесие расходов жидкости и газа;
равновесие тепловых потоков;
равновесие электрических потоков.
3. Кратко об угловых степенях свободы, используемых в пространственных элементах PRADIS
При моделировании пространственных механизмов и конструкций необходимо описывать не только поступательные, но и вращательные движения тел. Для этого в пространственных элементах PRADIS используются специальные угловые степени свободы, основанные на параметрах конечного вращения.
Такой подход позволяет избежать проблем, возникающих при использовании традиционных углов Эйлера, и обеспечивает устойчивое описание пространственного движения при любых положениях тела.
Описание конечного вращения
Согласно теореме Эйлера, любое изменение ориентации твердого тела может быть представлено одним поворотом вокруг некоторой оси, называемой осью конечного вращения.
Обозначим:
- e_1,\ e_2,\ e_3 — направляющие косинусы оси конечного вращения;
- F_i — угол конечного вращения.
Тогда вводятся четыре кинематических параметра:
x_1 = e_1 \cdot \sin\left(\frac{F_i}{2}\right)
x_2 = e_2 \cdot \sin\left(\frac{F_i}{2}\right)
x_3 = e_3 \cdot \sin\left(\frac{F_i}{2}\right)
x_4 = \cos\left(\frac{F_i}{2}\right)
Для этих параметров выполняется уравнение связи:
x_1^2 + x_2^2 + x_3^2 + x_4^2 = 1
Данная форма представления широко применяется в задачах пространственной механики, поскольку не содержит особых положений и не приводит к вырождению параметров при вращении тела.
Почему используются четыре параметра
При использовании углов Эйлера существуют положения, в которых происходит потеря одной степени свободы и возникают проблемы при численном расчете.
Параметры конечного вращения лишены этого недостатка:
- не обращаются в бесконечность;
- их производные остаются конечными;
- одинаково хорошо работают при любых ориентациях тела.
Именно поэтому данный подход принят в пространственных моделях PRADIS.
Угловые степени свободы в PRADIS
Угловые степени свободы пространственных элементов определяются через параметры конечного вращения следующим образом:
q_1 = x_1 \cdot L_q
q_2 = x_2 \cdot L_q
q_3 = x_3 \cdot L_q
q_4 = x_4 \cdot L_q
где
L_q = \sqrt{q_1^2 + q_2^2 + q_3^2 + q_4^2}
Первые три степени свободы являются внешними и доступны пользователю при построении модели.
Четвертая степень свободы является внутренней и используется только внутри математических моделей элементов.
Начальное значение потенциальной переменной, соответствующей внутренней степени свободы, принимается равным:
q_4 = 1
Потоковые переменные
Для первых трех угловых степеней свободы потоковыми переменными являются моменты относительно глобальных осей координат:
M_x,\ M_y,\ M_z
Четвертая потоковая переменная используется для поддержания нормировки параметров вращения.
Она определяется выражением:
i_4 = M_u \cdot \frac{dL_q}{dT}
где коэффициент
M_u = \frac{DABSI}{\sqrt{MSHEPS}}
одинаков для всех элементов такого типа.
Физически эта переменная препятствует накоплению численных ошибок и сохраняет условие нормировки параметров вращения.
Какие операции допустимы
Для пользователя работа с внешними угловыми степенями свободы практически не отличается от работы с обычными поступательными координатами.
В частности, допускается:
Базирование вращения
Можно запрещать вращение по выбранным угловым степеням свободы.
Например:
- запретить вращение вокруг оси X;
- разрешить вращение только вокруг оси Z;
- полностью зафиксировать ориентацию тела.
Такое ограничение эквивалентно уменьшению размерности вектора конечного вращения.
Создание вращательных связей
Связь между элементами можно задавать не по всем угловым степеням свободы, а только по необходимым.
Например:
- передавать вращение только вокруг одной оси;
- разрешить свободное вращение вокруг двух осей;
- реализовать различные варианты шарниров и кинематических соединений.
На что необходимо обратить внимание
Несмотря на внешнее сходство с поступательными степенями свободы, существует важное отличие.
Для пространственных элементов первая и вторая производные потенциальных переменных не являются непосредственно угловой скоростью и угловым ускорением.
То есть в общем случае неверны соотношения:
\frac{dq}{dt} = \omega
\frac{d^2 q}{dt^2} = \varepsilon
где:
- ω — угловая скорость;
- ε — угловое ускорение.
Это принципиальное отличие пространственной постановки от плоского вращения.
Особенности задания начальных условий
По этой причине при использовании моделей начальных условий, например модели VN, задаваемые значения не всегда соответствуют реальной начальной угловой скорости тела.
Аналогично при выводе результатов через ПРВП типа:
V или A
пользователь получает:
- первую производную потенциальной переменной;
- вторую производную потенциальной переменной.
Но не непосредственно угловую скорость и угловое ускорение.
Где получить реальные угловые скорости
Для анализа реальных угловых скоростей и ускорений следует использовать рабочие переменные соответствующих пространственных элементов.
Например, необходимые данные доступны в рабочем векторе модели:
J3O
и ряда других пространственных элементов.
Именно эти значения следует использовать при анализе динамики вращения, построении графиков и сравнении результатов расчета с экспериментальными данными.
Выводы
Угловые степени свободы пространственных элементов PRADIS основаны на параметрах конечного вращения и обеспечивают устойчивое описание вращательного движения без вырождений, характерных для углов Эйлера.
При работе с ними важно помнить:
- используются четыре параметра вращения;
- первые три степени свободы являются внешними;
- четвертая степень свободы является внутренней;
- потоковыми переменными являются моменты по глобальным осям;
- производные угловых координат не являются напрямую угловой скоростью и угловым ускорением;
- реальные значения угловых скоростей следует получать из рабочих переменных пространственных элементов.
Такой подход обеспечивает надежное и универсальное моделирование пространственных механических систем в PRADIS.