Описание методов PRADIS на примере
Рассмотрим колебательную систему с одной степенью свободы. Тело массы \(m\) соединяется с неподвижным основанием посредством упругой пружины. Его движение происходит в среде с вязким сопротивлением под действием внешней силы (рис. 1.1).
Рис. 1.1
Для начального момента времени заданы положение тела в пространстве и его начальная скорость. Требуется для каждого момента времени в интервале от \(t_0\) до \(t_{\text{конеч}}\) определить перемещение, скорость и ускорение тела.
Поскольку необходимо проиллюстрировать решение нелинейной задачи, будем считать, что сила вязкого сопротивления пропорциональна квадрату относительной скорости концов демпфера (рис. 1.2б). Сила упругости пружины линейно зависит от смещения (рис. 1.2а), воздействующая сила имеет синусоидальный характер (рис. 1.2в).
Рис. 1.2
Зависимости, позволяющие получить дифференциальное уравнение движения, в соответствии со вторым законом Ньютона имеют вид:
или
где с учетом выбранных направлений (рис. 1.3):
Рис. 1.3
Здесь \(x\), \(\nu\) и \(a\) — соответственно перемещение, скорость и ускорение тела.
Подстановка соотношений (1.3) в уравнение (1.2) дает:
Учитывая, что
получаем дифференциальное уравнение движения тела:
Численное интегрирование
Использование численного подхода к интегрированию уравнения (1.6) предполагает нахождение приближенного решения для определенных моментов времени. Временная ось представляется совокупностью точек \(t_0, t_1, t_2, \dots, t_i, t_{i+1}, \dots, t_n\) (рис. 1.4), в каждой из которых ищется приближенное решение.
Рис. 1.4
Интегрирование осуществляется последовательно; выбор величины очередного шага \(\Delta t_i\) зависит как от требуемых показателей точности, так и от результатов интегрирования на предыдущих шагах.
Таким образом, использование численного подхода позволяет перейти от непрерывных значений \(x\), \(\nu\), \(a\) на всем промежутке времени от \(t_0\) до \(t_{\text{конеч}}\) к совокупности дискретных значений \(x_i\), \(\nu_i\), \(a_i\) для определенных моментов времени \(t_i\). При этом алгебраические формулы выбранного метода интегрирования заменяют дифференциальные соотношения (1.5).
Формулы неявного одношагового метода Штермера устанавливают следующую зависимость для переменных \(x_i\), \(\nu_i\) по известным с предыдущего шага значениям \(x_{i-1}\), \(\nu_{i-1}\) [1]:
где
\(t_0\) — начальное время,
\(x_0\), \(\nu_0\) — начальные значения перемещения и скорости,
\(t_n\) — конечное время.
Для начального момента времени \(t_0\) значения \(x_0\) и \(\nu_0\) должны быть известны. Определим значения \(x_1\), \(\nu_1\), \(a_1\) для момента времени \(t_1 = t_0 + \Delta t_1\).
Уравнение (1.4) для момента времени \(t_i\) имеет вид:
Дополняя это соотношение формулами метода интегрирования (1.7), получаем для момента времени \(t_1\) замкнутую систему уравнений:
Приведем полученную систему к одному уравнению, выразив неизвестные \(x_1\) и \(a_1\) через \(\nu_1\):
Получаем:
Группируя сомножители при одинаковых степенях неизвестного \(\nu_1\), соотношение можно записать в виде:
где
Соотношение (1.12) сохраняет свой вид для любого момента времени \(t_i\) при соответствующей замене подстрочных индексов (\(1 \to i\), \(0 \to i-1\)).
Итак, использование формул метода интегрирования позволяет преобразовать исходное дифференциальное уравнение (1.6) к нелинейному уравнению вида (1.12), которое решается на каждом шаге по времени.
Решение нелинейного уравнения методом Ньютона
Уравнение (1.12) решается методом Ньютона. Рассмотрим последовательность действий при решении нелинейного уравнения.
Рассматривается уравнение вида:
где \(f(z)\) — нелинейная функция относительно неизвестного \(z\).
Рис. 1.5
Алгоритм численного решения включает следующие шаги:
Выбор начального приближения к решению — величины \(z^0\).
Организация последовательности итераций, на каждой из которых уточняется полученное на предыдущей итерации значение \(z\) по схеме (рис. 1.5а):
\[\begin{split}\begin{cases} z^{j} = z^{j-1} + \Delta z^{j} \\[6pt] \Delta z^{j} = -\dfrac{f(z^{j-1})}{f'(z^{j-1})} \end{cases} \qquad (1.14)\end{split}\]где \(f(z^{j-1})\) — значение функции \(f(z)\) при \(z = z^{j-1}\), \(f'(z^{j-1})\) — значение производной \(df/dz\) при \(z = z^{j-1}\).
Проверка на каждой итерации условия прекращения итераций (рис. 1.5б):
\[|z^{j} - z^{j-1}| = |\Delta z^{j}| \leq \delta_z \qquad (1.15)\]\[|f(z^{j})| \leq \delta_f\]где \(\delta_f\) — допустимая невязка (отклонение от нуля) правой части уравнения (1.13); \(\delta_z\) — допустимая величина отличия решения на двух соседних итерациях.
Проверка ограничения на максимально допустимое количество итераций:
\[j \leq j_{\text{max}} \qquad (1.16)\]
Геометрически решение уравнения (1.13) сводится к отысканию абсциссы точки пересечения с осью \(z\) кривой \(f(z)\). На каждой \(j\)-й итерации метода Ньютона решение заменяется отысканием точки пересечения касательной к кривой \(f(z)\) с осью \(z\); касательная строится для \(z = z^{j-1}\).
Первый шаг: численный расчет
Возвращаемся к численному решению уравнения (1.12). Обозначив \(z = \nu_1\), имеем:
или
где
Для решения уравнения (1.17) методом Ньютона потребуется выражение для производной функции \(f(z)\):
Зададимся исходными данными:
\(k = 20000\) Н/м,
\(\mu = 1000\) Н·с²/м²,
\(m = 0.1\) кг,
\(Q = 1000\),
\(T = 0.2\pi\),
\(F = 1000 \sin(10t)\).
Начальные условия и шаг интегрирования:
Тогда, согласно (1.12):
Таким образом,
Зададим значения допустимых погрешностей:
Максимально допустимое количество итераций примем равным \(5\).
Начальное приближение к решению выберем \(z^0 = 0\).
Первая итерация.
Проверка завершения итераций:
Переход на следующую итерацию.
Вторая итерация.
Проверка:
Третья итерация.
Четвертая итерация.
Оба условия (1.15) удовлетворены, ограничение (1.16) не превышено. Решение достигнуто. Мы вычислили значение скорости для момента времени \(t_1\):
Воспользовавшись формулами (1.10), определим значения ускорения и перемещения:
Второй шаг: выбор величины шага
Решение для момента времени \(t_1\) получено. Сделаем еще один шаг по времени, чтобы проиллюстрировать выбор величины шага. Уравнения (1.9)–(1.12) справедливы для любого момента времени с учетом соответствующей замены подстрочных индексов. Для второго шага имеем:
Приводим эту систему к одному уравнению относительно скорости:
где
Предварительно величину шага \(\Delta t_2\) принимаем равной \(\Delta t_1\), т.е. \(0.001\) с. Тогда, с учетом исходных данных и полученного на первом шаге решения, подсчитаем коэффициенты:
Опять имеем нелинейное уравнение:
которое решаем методом Ньютона.
В этом месте необходимо обратить особое внимание на выбор начального приближения к решению в алгоритме метода Ньютона.
За начальное приближение примем такое значение скорости, которое тело имело бы в момент времени \(t_2\), если бы ускорение с момента времени \(t_1\) не изменилось:
Это так называемый явный шаг (прогноз). Величина скорости, полученная явным шагом, будет использована не только как начальное приближение в методе Ньютона, но и при оценке величины выбранного шага по времени.
Итак, начальное приближение (прогноз):
Опуская подробные выкладки (они аналогичны приведенным для первого шага), итерации метода Ньютона приводят к следующей последовательности:
начальное приближение: \(\nu_2^0 = 0.11826\)
первая итерация: \(\nu_2^1 = 0.11172\)
вторая итерация: \(\nu_2^2 = 0.11159\) — решение достигнуто.
Таким образом, при величине шага \(\Delta t_2 = 0.001\) с скорость для момента времени \(t_2\) составляет:
Оценка локальной погрешности
Погрешность метода интегрирования на \(i\)-м шаге, называемую локальной погрешностью, будем оценивать по формуле:
где \(\nu_i^p\) — явный прогноз величины скорости на \(i\)-м шаге, определяемый формулой
\(\nu_i^c\) — значение скорости, полученное в результате итерационного решения с использованием неявной формулы
Рис. 1.6
Соотношение (1.28) уже применялось при выборе начального приближения (1.26). Смысл оценки локальной погрешности поясняет рис. 1.6.
В момент времени \(t_{i-1}\) мы находимся в точке \(\nu_{i-1}\). Если для вычисления \(\nu_i\) воспользоваться явной формулой (1.28), то точка \(\nu_i = \nu_i^p\) будет лежать на касательной к кривой \(\nu(t)\) в точке \(t_{i-1}\), поскольку \(a_{i-1}\) есть тангенс угла наклона этой касательной.
При вычислении \(\nu_i\) с использованием формулы (1.29) требуется значение \(a_i\), т.е. тангенса угла наклона касательной к кривой \(\nu(t)\) в точке \(t_i\). Так как в момент времени \(t_{i-1}\) мы ничего не знаем о поведении функции \(\nu(t)\) при \(t = t_i\), то \(\nu_i = \nu_i^c\) вычисляется путем совместного решения системы уравнений (1.7), куда входит соотношение (1.29). При этом требуется несколько ньютоновских итераций.
Рис. 1.6 показывает, что явный прогноз \(\nu_i^p\) и скорректированное решение \(\nu_i^c\) лежат по разные стороны от кривой \(\nu(t)\), проходящей через точку \(t_{i-1}\). Чем больше разница между \(\nu_i^c\) и \(\nu_i^p\), тем сильнее на текущем шаге график скорости отличается от прямой и тем выше погрешность интегрирования. Уменьшение величины шага \(\Delta t_i\) приводит к уменьшению локальной погрешности, оцениваемой по формуле (1.27).
Управление величиной шага
Расчет локальной погрешности позволяет оценить приемлемость сделанного шага по времени и рекомендовать величину следующего шага.
Задается значение предельно допустимой локальной погрешности на шаге интегрирования \(\delta_l\). По результатам очередного \(i\)-го шага сравниваются значения \(\delta_l\) и фактически полученной локальной погрешности \(lp_i\).
Если \(lp_i \leq \delta_l\), то сделанный шаг признается успешным. Осуществляется переход на следующий шаг; величина его для одношаговых методов интегрирования первого порядка точности (к которым относятся формулы (1.7)) выбирается по зависимости:
\[\Delta t_{i+1} = c \Delta t_i \sqrt{\frac{\delta_l}{lp_i}} \qquad (1.30)\]где \(\Delta t_i\) — величина совершенного шага, \(\Delta t_{i+1}\) — рекомендуемое значение следующего шага, \(c\) — поправочный коэффициент (\(c < 1\)).
Если \(lp_i > \delta_l\), то значение шага \(\Delta t_i\) слишком велико и не обеспечивает требуемой точности. Необходимо выполнить расчет на \(i\)-м шаге повторно с уменьшенным значением \(\Delta t_i\). Для выбора величины шага также используется формула (1.30), но полученное значение применяется для повторного расчета на текущем шаге.
Вернемся к примеру. Для второго шага интегрирования имеем:
Локальная погрешность:
Приняв допустимую погрешность \(\delta_l = 0.001\), констатируем, что \(lp_2 > \delta_l\). Выполненный шаг с \(\Delta t_2 = 0.001\) с не обеспечивает требуемой точности; необходимо повторить расчет на втором шаге с уменьшенным значением \(\Delta t_2\).
Рекомендуемое значение для повторного расчета определим по формуле (1.30) с коэффициентом \(c = 0.8\):
Результаты повторного расчета с шагом \(\Delta t_2 = 0.438 \times 10^{-3}\) с:
Поскольку \(lp_2 < \delta_l\), второй шаг считается успешным. Дополним значение скорости \(\nu_2 = 0.08509\) м/с значениями ускорения и перемещения, используя формулы связи (1.7):
Геометрическая интерпретация
В момент времени \(t_{i-1}\) мы находимся в точке \(\nu_{i-1}\). Через нее проходит кривая \(\nu(t)\) — так называемая интегральная кривая для момента времени \(t_{i-1}\), т.е. график скорости, соответствующий точному решению уравнения (1.6) при начальном условии
Поскольку уравнение (1.6) решается приближенно, на каждом \(i\)-м шаге численного интегрирования происходит переход с одной интегральной кривой, удовлетворяющей начальному условию (1.31), на другую интегральную кривую, представляющую собой точное решение уравнения (1.6) при начальном условии
На рис. 1.6 интегральная кривая для \(t = t_i\) обозначена пунктирной линией.
Результатом численного решения служит ломаная, проходящая через совокупность интегральных кривых, каждая из которых является точным решением уравнения (1.6) при начальных условиях, определяемых численным решением на текущем шаге (рис. 1.8).
Рис. 1.8
Основные этапы численного решения
Подытожим основные моменты, существенные с точки зрения численного анализа рассмотренного примера:
Формирование дифференциального уравнения, описывающего поведение системы:
\[kx + \mu \frac{dx}{dt} \left| \frac{dx}{dt} \right| + m \frac{d^2 x}{dt^2} - Q \sin \frac{2\pi}{T} t = 0\]При формировании уравнения использован второй закон Ньютона.
Представление уравнения в форме, не содержащей явно дифференциальных соотношений:
\[\begin{split}\begin{cases} kx + \mu \nu |\nu| + ma - Q \sin \dfrac{2\pi}{T} t = 0 \\[8pt] \nu = \dfrac{dx}{dt} \\[8pt] a = \dfrac{d\nu}{dt} = \dfrac{d^2 x}{dt^2} \end{cases}\end{split}\]Замена дифференциальной связи алгебраическими уравнениями выбранного метода интегрирования:
\[\begin{split}\begin{cases} k x_i + \mu \nu_i |\nu_i| + m a_i - Q \sin \dfrac{2\pi}{T} t_i = 0 \\[8pt] x_i = x_{i-1} + \nu_{i-1} \Delta t_i + a_i \dfrac{\Delta t_i^2}{2} \\[8pt] \nu_i = \nu_{i-1} + a_i \Delta t_i \end{cases}\end{split}\]Приведение системы к одному уравнению относительно \(\nu_i\):
\[\alpha \nu_i |\nu_i| + \beta \nu_i + \gamma = 0\]Решение нелинейного уравнения методом Ньютона. Начальное приближение определяется формулой явного прогноза.
Оценка точности интегрирования путем контроля локальной погрешности. При неудовлетворительной величине локальной погрешности расчет на текущем шаге повторяется с уменьшенным шагом.
Вычисление ускорения и перемещения по уравнениям связи после успешного завершения шага.
Выбор величины следующего шага на основе соотношения допустимой и фактической локальной погрешности.
Автоматическое формирование математической модели
На рассмотренном примере мы обозначили канву численного решения, которой придерживается алгоритм вычислительного ядра PRADIS. Однако формирование дифференциального уравнения проводилось «вручную», а некоторые вопросы решались неформально (например, аналитическое определение производной \(df(z)/dz\)). Поэтому продолжим рассмотрение методов PRADIS с выяснения принципов автоматического формирования системы дифференциальных уравнений.
Вернемся к системе на рис. 1.1. Уравнение равновесия:
Для каждого \(i\)-го момента времени уравнение может быть записано в виде:
где
(Обратите внимание: знаки отличаются от (1.3), поскольку в PRADIS при рассмотрении условий равновесия суммируются усилия, действующие со стороны системы на элементы.)
Связь между \(x_i\), \(\nu_i\), \(a_i\) задается формулами метода интегрирования (1.7), которые воспроизведем еще раз:
Соотношение (2.2) представляет собой алгебраическое нелинейное уравнение, которое можно решать методом Ньютона. Имеем:
где
В качестве переменной \(z\) может быть выбрана любая из компонент \(x_i\), \(\nu_i\) или \(a_i\), поскольку они связаны соотношениями (2.4). Примем, как и ранее, \(z = \nu_i\).
На каждой итерации необходимо вычислить \(f(z)\) и \(df(z)/dz\).
Вычисление \(f(z)\) сводится к суммированию текущих значений сил при текущих значениях \(x_i\), \(\nu_i\), \(a_i\) (\(i\) — номер шага по времени, \(j\) — номер итерации по Ньютону). Вычисляемое значение \(f(z)\) представляет собой погрешность выполнения условия равновесия, которую необходимо уменьшить до допустимых пределов.
Производная \(df(z)/dz\):
Каждую из производных следует представить как производную сложной функции. Для силы \(F_i^c\):
Поскольку \(z = \nu_i\):
Аналогично для остальных сил (2.7a).
Используем уравнения связи (2.4) для получения зависимостей \(a_i\) и \(x_i\) от \(\nu_i\):
Дифференцируем (2.8) по \(\nu_i\):
Вычисление частных производных
С использованием зависимостей (2.3) вычислим частные производные, входящие в выражения (2.7) и (2.7а).
Для силы \(F_i^c\):
поскольку \(F_i^c\) не зависит от перемещения, скорости и ускорения тела (2.10).
Для силы \(F_i^y\):
Для силы \(F_i^в\):
Для силы \(F_i^u\):
Подставляя найденные значения частных производных и коэффициенты из (2.9) в формулы (2.7) и (2.7а), получаем:
Суммируя слагаемые в формуле (2.6), приходим к выражению:
Полученный результат структурно аналогичен ранее выведенной формуле (1.19), коэффициенты которой определяются соотношениями (1.12).
Преимущества автоматического подхода
Теперь, имея возможность вычислять \(f(z)\) и \(df(z)/dz\), продолжить расчет по описанному алгоритму не представляет труда. Какими же преимуществами окупаются дополнительные выкладки?
Дифференциальное уравнение движения выписывать не пришлось.
Развернутая форма нелинейного уравнения вида (1.12) оказалась невостребованной.
Функциональная зависимость для производной \(df(z)/dz\) не потребовалась.
Рис. 2.1
Для формирования математической модели была использована следующая информация (далее — «перечень»):
Сведения о стыковке элементов схемы (рис. 2.1).
Условие равновесия сил для \(i\)-го момента времени — нелинейное алгебраическое уравнение относительно \(x_i\), \(\nu_i\), \(a_i\):
\[\sum_{k=1}^{4} F_i^{(k)} = 0 \qquad (2.16)\](топологическое уравнение, определяемое структурой связей).
Компонентные уравнения — выражения, позволяющие определить усилия в каждом элементе как функции перемещения, скорости, ускорения и времени (2.3).
Частные производные усилий по перемещению, скорости, ускорению (2.10)–(2.13).
Алгебраические уравнения связи \(x\), \(\nu\), \(a\) — формулы метода интегрирования (2.4).
Выбор базисной переменной \(z\) (из \(x\), \(\nu\) или \(a\)) для решения нелинейного уравнения.
Этой информации достаточно для реализации машинных алгоритмов формирования математической модели объекта.
Модели элементов
Пусть имеется техническая система, процессы в которой требуют анализа. Расчетная схема соответствует рис. 1.1.
Систему необходимо представить в виде совокупности элементов, стыкующихся по общим степеням свободы (узлам). Количество степеней свободы у каждого элемента определяется его типом (рис. 2.2).
Рис. 2.2
Существует понятие модели элемента. Пользователь собирает модель системы из моделей отдельных элементов, подобно детскому конструктору. Задача пользователя — правильность сборки; остальные вопросы формирования математической модели решаются программным обеспечением.
Для расчетной схемы (рис. 1.1) пользователь находит модели элементов в библиотеке и описывает структуру:
$ FRAGMENT: Пример
# BASE: 1
# STRUCTURE:
Пружина' K (1 2; Коэффициент жесткости)
Нелинейный демпфер ' MUNL (1 2; Коэффициент вязкости)
Масса ' M (2; Масса тела)
Воздействие ' FSIN (2 1; Q, T, начальная фаза)
Этим описанием пользователь сообщает программному комплексу PRADIS всю необходимую информацию по стыковке элементов схемы (пункт 1 перечня). Выбраны модели K, MUNL, M, FSIN; они соединены надлежащим образом; узел 1 описан как неподвижный.
В процессе обработки описания структуры определяется размерность системы уравнений — количество узлов, в которых должны выполняться условия равновесия. В рассматриваемом примере два узла, один из которых закреплен. На этапе формирования математической модели структура данных готовится по обоим узлам, но на этапе расчета уравнение для закрепленного узла исключается, а его кинематические характеристики обнуляются.
Этап численного интегрирования представляет собой последовательность шагов по времени, каждый из которых сводится к решению нелинейного уравнения равновесия вида (2.16). Для этого необходима информация, приведенная в пунктах 3–6 перечня.
Входными данными для любой модели элемента являются:
неизменный список параметров модели;
текущие значения перемещений, скоростей и ускорений узлов, с которыми элемент соединен.
По этим данным модель элемента обязана для текущего момента времени вычислить:
усилия, действующие со стороны системы на элементы (вектор усилий элемента) — пункт 3 перечня;
частные производные усилий по перемещениям, скоростям и ускорениям узлов элемента (матрицу Якоби элемента) — пункт 4 перечня.
Если элемент имеет \(N\) степеней свободы, длина вектора усилий элемента равна \(N\), а якобиан элемента имеет размер \(N \times N \times 3\).
Рис. 2.3
Например, разработчик двухузловой модели одномерной упругой пружины реализовал следующие зависимости (рис. 2.3):
Для любого момента времени модель элемента по переданным значениям перемещений, скоростей и ускорений вычисляет усилия на концах пружины и частные производные. Вектор усилий состоит из двух элементов, якобиан — из 12.
Поскольку узел 1 пружины закреплен, востребуется только информация, связанная с незакрепленным вторым узлом:
Эта информация позволяет учесть вклад пружины при решении нелинейного уравнения (2.5). Вклад остальных элементов учитывается аналогично.
Алгоритм выполнения шага интегрирования
Продолжим ранее начатое интегрирование, сделав третий шаг по времени, используя формальный алгоритм, базирующийся на формулах (2.5)–(2.15).
По результатам второго шага (\(t_2 = 1.438 \times 10^{-3}\) с):
Величина шага \(\Delta t_2 = 0.438 \times 10^{-3}\) с, локальная погрешность \(lp_2 = 0.000018\).
Рекомендуемое значение \(\Delta t_3\) определяем по формуле (1.30) с \(c = 0.8\) и \(\delta_l = 0.001\):
Действуем по схеме, представленной на рис. 2.4а–2.4в.
Рис. 2.4а
Рис. 2.4б
Рис. 2.4в
Перед третьим шагом известны значения \(t_2\), \(x_2\), \(\nu_2\), \(a_2\), \(\Delta t_3\).
Определение коэффициентов приведения якобиана (2.9):
\[\frac{dx_3}{d\nu_3} = \frac{\Delta t_3}{2} = \frac{2.63 \times 10^{-3}}{2} = 1.31 \times 10^{-3}\]\[\frac{d\nu_3}{d\nu_3} = 1,\quad \frac{da_3}{d\nu_3} = \frac{1}{\Delta t_3} = \frac{1}{2.63 \times 10^{-3}} = 380.2\]Начальное приближение (явный прогноз по формуле (1.28)):
\[\nu_3^0 = \nu_2 + a_2 \Delta t_3 = 0.08509 + 59.21 \cdot 2.63 \times 10^{-3} = 0.24081\]Значения \(x_3^0\) и \(a_3^0\), необходимые для расчета в моделях элементов:
\[a_3^0 = a_2 = 59.21\]\[x_3^0 = x_2 + \nu_2 \Delta t_3 + a_3^0 \frac{\Delta t_3^2}{2} = 6.12 \times 10^{-5} + 0.08509 \cdot 2.63 \times 10^{-3} + 59.21 \cdot \frac{(2.63 \times 10^{-3})^2}{2} = 46.5 \times 10^{-5}\]Первая итерация Ньютона (см. рис. 2.4в).
Обращение к моделям элементов. Вычисление вектора сил и якобиана каждого элемента по текущим значениям \(x_3^0\), \(\nu_3^0\), \(a_3^0\) для незакрепленного узла:
Пружина:
\[F^y = k x = 20000 \cdot 46.5 \times 10^{-5} = 9.3\]\[\frac{\partial F^y}{\partial x} = k = 20000,\quad \frac{\partial F^y}{\partial \nu} = \frac{\partial F^y}{\partial a} = 0\]Демпфер:
\[F^в = \mu \nu |\nu| = 1000 \cdot 0.24081 \cdot |0.24081| = 58.0\]\[\frac{\partial F^в}{\partial x} = 0,\quad \frac{\partial F^в}{\partial \nu} = 2\mu |\nu| = 2 \cdot 1000 \cdot 0.24081 = 481.6,\quad \frac{\partial F^в}{\partial a} = 0\]Точечная масса:
\[F^u = m a = 0.1 \cdot 59.21 = 5.9\]\[\frac{\partial F^u}{\partial x} = \frac{\partial F^u}{\partial \nu} = 0,\quad \frac{\partial F^u}{\partial a} = m = 0.1\]Прикладываемая сила:
\[F^c = Q \sin \frac{2\pi}{T} t = -1000 \sin \frac{2\pi}{0.2\pi} (1.438 \times 10^{-3} + 2.63 \times 10^{-3}) = -40.6\]\[\frac{\partial F^c}{\partial x} = \frac{\partial F^c}{\partial \nu} = \frac{\partial F^c}{\partial a} = 0\]Суммирование сил (значение \(f(z)\) на текущей итерации):
\[\sum F = F^y + F^в + F^u + F^c = -40.6 + 9.3 + 58.0 + 5.9 = 32.6\]Вычисление \(df(z)/dz\) по формулам (2.6), (2.7), (2.7а):
\[\frac{dF^c}{dz} = 0 \cdot 1.31 \times 10^{-3} + 0 \cdot 1 + 0 \cdot 380.2 = 0\]\[\frac{dF^y}{dz} = 20000 \cdot 1.31 \times 10^{-3} + 0 \cdot 1 + 0 \cdot 380.2 = 26.2\]\[\frac{dF^в}{dz} = 0 \cdot 1.31 \times 10^{-3} + 481.6 \cdot 1 + 0 \cdot 380.2 = 481.6\]\[\frac{dF^u}{dz} = 0 \cdot 1.31 \times 10^{-3} + 0 \cdot 1 + 0.1 \cdot 380.2 = 38.0\]\[\frac{df(z)}{dz} = 0 + 26.2 + 481.6 + 38.0 = 545.8\]Приращение \(\Delta z^1\):
\[\Delta z^1 = -\frac{f(z^0)}{f'(z^0)} = -\frac{32.6}{545.8} = -0.05973\]Очередное приближение:
\[z^1 = z^0 + \Delta z^1 = 0.24081 - 0.05973 = 0.18108\]Уточнение \(x_3^1\), \(\nu_3^1\), \(a_3^1\) по формулам (2.8):
\[\nu_3^1 = z^1 = 0.18108\]\[a_3^1 = \frac{\nu_3^1 - \nu_2}{\Delta t_3} = \frac{0.18108 - 0.08509}{2.63 \times 10^{-3}} = 36.5\]\[x_3^1 = x_2 + \frac{\nu_2 + \nu_3^1}{2} \Delta t_3 = 6.12 \times 10^{-5} + \frac{0.08509 + 0.18108}{2} \cdot 2.63 \times 10^{-3} = 41.1 \times 10^{-5}\]Проверка условий завершения итераций:
\[|f(z^0)| = 32.6 > \delta_f = 0.1,\quad |\Delta z^1| = 0.05973 > \delta_z = 0.001\]Итерации необходимо продолжить.
Проверка лимита итераций:
\[j = 1 < j_{\text{max}} = 5\]Увеличение счетчика итераций: \(j = 2\).
Вторая итерация дает:
\[f(z^1) = -3.6,\quad \Delta z^2 = -0.00858\]\[x_3^2 = 39.8 \times 10^{-5},\quad \nu_3^2 = 0.17250,\quad a_3^2 = 33.2\]Условия завершения не выполнены.
Третья итерация (последняя):
\[|f(z^2)| = 0.07 < \delta_f,\quad |\Delta z^3| = 0.00018 < \delta_z\]\[x_3^3 = 39.8 \times 10^{-5},\quad \nu_3^3 = 0.17232,\quad a_3^3 = 33.2\]Оценка локальной погрешности:
\[lp_3 = \left| \frac{\nu_3^p - \nu_3^c}{2} \right| = \left| \frac{0.24081 - 0.17232}{2} \right| = 0.034\]Корректировка правила выбора шага. Практика расчетов показала, что формула (1.30) приемлема только при соотношении \(\delta_l/lp_i\) вблизи единицы. При значительных отклонениях рекомендуемое значение шага часто завышено. Используем уточненное правило:
\[\begin{split}\Delta t_{\text{рек}} = \begin{cases} c \cdot \Delta t_i \cdot \dfrac{\delta_l}{lp_i}, & \text{если } \dfrac{\delta_l}{lp_i} < 0.25 \\[12pt] c \cdot \Delta t_i \cdot \sqrt[4]{\dfrac{\delta_l}{lp_i}}, & \text{если } \dfrac{\delta_l}{lp_i} > 7 \\[12pt] c \cdot \Delta t_i \cdot \sqrt{\dfrac{\delta_l}{lp_i}}, & \text{если } 0.25 \leq \dfrac{\delta_l}{lp_i} \leq 7 \end{cases} \qquad (2.23)\end{split}\]Для нашего случая:
\[\frac{\delta_l}{lp_3} = \frac{0.001}{0.034} = 0.03 < 0.25\]\[\Delta t_{\text{рек}} = 0.8 \cdot 2.63 \times 10^{-3} \cdot \frac{0.001}{0.034} = 0.061 \times 10^{-3}\ \text{с}\]Повторный расчет с \(\Delta t_3 = 0.061 \times 10^{-3}\) с дает для момента времени \(t_3 = t_2 + \Delta t_3 = 1.499 \times 10^{-3}\) с:
\[x_3 = 6.64 \times 10^{-5}\ \text{м},\quad \nu_3 = 0.08862\ \text{м/с},\quad a_3 = 58.1\ \text{м/с}^2\]Локальная погрешность в пределах нормы. Рекомендуемое значение для следующего шага: \(\Delta t_{\text{рек}} = 0.264 \times 10^{-3}\) с.
Расчет на третьем шаге завершен.
Разделение функций в вычислительном ядре
Основное, на что следует обратить внимание, — это разделение функций между программой интегрирования и программами реализации моделей элементов.
Программе интегрирования, работающей по алгоритму рис. 2.4а–2.4в, не важно, какие процессы интегрируются. Ее зависимость от моделей элементов сводится только к своевременному получению векторов сил и матриц якобианов. Какие свойства отдельных элементов отражает эта информация, программу интегрирования не касается.
Модели элементов, в свою очередь, имеют свой уровень независимости при четко очерченных обязанностях перед «верхами». Физические свойства отдельного элемента отражаются в компонентных уравнениях на уровне модели элемента, а программа интегрирования работает на уровне уравнений равновесия потоков.
Такое разграничение функций определяет универсальность вычислительного ядра PRADIS — принципиальную возможность расчета любых объектов, процессы в которых подчиняются законам равновесия потоковых переменных (равновесие сил, электрических и тепловых потоков, расходов жидкости и газа).
Угловые степени свободы в пространственных элементах PRADIS
Согласно теореме Эйлера, твердое тело из одного углового положения в другое можно перевести одним поворотом вокруг некоторой оси, называемой осью конечного вращения. Обозначим \(e_1, e_2, e_3\) — направляющие косинусы оси конечного вращения, \(\Phi\) — угол конечного вращения. Введем четыре кинематических параметра, описывающих угловое движение твердого тела [1,2]:
Эти параметры связаны уравнением:
В отличие от любой совокупности трех кинематических параметров (например, углов Эйлера), указанные четыре параметра не вырождаются при любом положении твердого тела.
Угловые степени свободы, принятые в пространственных элементах PRADIS, выражаются через кинематические параметры (1) следующим образом:
где
Первые три степени свободы являются внешними для моделей элементов, четвертая — внутренняя, скрытая от пользователя. Начальное значение потенциальной переменной, соответствующей внутренней степени свободы, устанавливается в моделях элементов равным 1.
Потоковыми переменными для первых трех степеней свободы являются моменты по глобальным осям X, Y, Z. Четвертая (внутренняя) потоковая переменная сдерживает изменение во времени величины \(L_q\):
где \(M_u\) — коэффициент пропорциональности, одинаковый для всех степеней свободы такого рода и принимаемый в моделях элементов равным
Для пользователя корректны практически все приемы, характерные для поступательного движения:
можно запрещать движение по выбранным угловым степеням свободы (базируя соответствующие узлы), что равносильно сокращению размерности вектора, направленного по оси конечного вращения;
связь по вращению между точками стыкуемых элементов можно осуществлять не по всем трем степеням свободы, а только по выбранным.
Следует проявлять осторожность: в отличие от плоского вращения, первая и вторая производные от потенциальных переменных (3) не являются угловой скоростью и угловым ускорением соответственно. Начальные условия, задаваемые моделью VN, не определяют, в общем случае, начальную угловую скорость. ПРВП типов V и A выводят текущие значения первой и второй производной от потенциальной переменной. Значения угловых скоростей и ускорений доступны из рабочего вектора некоторых моделей элементов, в частности J3O.
Литература
Виттенбург Й. Динамика систем твердых тел. — М.: Мир, 1980.