Методы PRADIS: руководство по численному моделированию

Введение

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

Основная задача документа - показать:

• как математическая модель преобразуется в расчетную схему;

• как выполняется шаг интегрирования;

• как формируется система нелинейных алгебраических уравнений;

• как применяется метод Ньютона;

• какую роль играют модели элементов;

• как обеспечивается универсальность вычислительного ядра PRADIS.

Весь материал изложен на сквозном числовом примере с конкретными итерациями, коэффициентами и критериями останова.

Обозначения

Основные обозначения, используемые в документе

Символ	Единица измерения (в примере)	Описание
\(x\)	м	Перемещение (координата) тела
\(\nu\)	м/с	Скорость тела
\(a\)	м/с²	Ускорение тела
\(z\)	м/с	Базисная переменная (в методе Ньютона — скорость \(\nu_i\))
\(F^{\text{с}}\)	Н	Внешняя синусоидальная сила
\(F^{\text{у}}\)	Н	Сила упругости пружины
\(F^{\text{в}}\)	Н	Сила вязкого сопротивления
\(F^{\text{и}}\)	Н	Сила инерции (\(-ma\))
\(\delta_f\)	Н	Критерий останова по невязке уравнения
\(\delta_z\)	м/с	Критерий останова по приращению переменной
\(\delta_l\)	м/с	Допустимая локальная погрешность
\(lp_i\)	м/с	Оценка локальной погрешности на шаге \(i\)

1. Постановка задачи

На примере колебательной системы с одной степенью свободы показана ключевая идея алгоритма вычисленного ядра PRADIS: как исходная элементная схема (пружина, демпфер, масса, внешняя сила) автоматически преобразуется на каждом шаге по времени в систему нелинейных алгебраических уравнений относительно скорости, и как эта система решается методом Ньютона с автоматическим управлением шагом.

Весь материал изложен на сквозном числовом примере, который последовательно проводится через три шага интегрирования:

• Шаг 1 (раздел 3.3) — первый шаг интегрирования от \(t_0\) до \(t_1\). Демонстрируется решение нелинейного уравнения скорости методом Ньютона без управления шагом.

• Шаг 2 (раздел 4.3) — второй шаг интегрирования от \(t_1\) до \(t_2\). Показывается оценка локальной погрешности и автоматическое уменьшение шага при превышении допустимой погрешности.

• Шаг 3 (раздел 7) — третий шаг интегрирования от \(t_2\) до \(t_3\). Иллюстрируется повторное применение управления шагом в более сложной ситуации.

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

1.1 Описание физической системы

Рассматривается колебательная система с одной степенью свободы: тело массы \(m\), соединённое с неподвижным основанием упругой пружиной. Движение происходит в среде с вязким сопротивлением под действием внешней силы (рис. 1.1). Начальные условия: заданы положение тела и его начальная скорость.

Рис. 1.1 — Колебательная система

Требуется определить перемещение \(x(t)\), скорость \(\nu(t)\) и ускорение \(a(t)\) на интервале \([t_0, t_{\text{конеч}}]\).

1.2 Характеристики модели

В качестве иллюстрации решения нелинейной задачи примем следующие модели элементов: демпфер — с квадратичным трением (рис. 1.2.б), пружина — линейная (рис. 1.2.а), внешняя нагрузка — гармоническая (рис. 1.2.в). Их характеристики приведены на рис.1.2

Рис. 1.2 — Характеристики элементов: а) упругая сила, б) сила вязкого сопротивления, в) внешняя сила

Каждый из этих элементов создаёт свою силу:

• Сила упругости пружины (рис. 1.2, а) — линейна: \(F^{\text{у}} = -kx\)

• Сила вязкого сопротивления (рис. 1.2, б) пропорциональна квадрату относительной скорости: \(F^{\text{в}} = -\mu \nu |\nu|\)

• Внешняя сила (рис. 1.2, в) — синусоидальная: \(F^{\text{с}} = Q \sin(2\pi t / T)\)

В соответствии с принципом Даламбера, добавим силу инерции:

\[F^{\text{и}} = -ma\]

Таким образом, все силы, действующие на систему, можно записать в виде системы уравнений:

\[\begin{split}\begin{cases} F^{\text{с}} = Q \sin \dfrac{2\pi}{T} t \\[6pt] F^{\text{у}} = -kx \\[6pt] F^{\text{в}} = -\mu \nu \left| \nu \right| \\[6pt] F^{\text{и}} = -ma \end{cases} \qquad (1.1)\end{split}\]

1.3 Дифференциальное уравнение движения

Зависимости, позволяющие получить дифференциальное уравнение движения, в соответствии со вторым законом Ньютона, имеют вид:

\[ma = F^{\text{с}} + F^{\text{у}} + F^{\text{в}} \qquad (1.2)\]

Удобной для численного решения является форма равновесия, где сумма всех сил (включая силу инерции) равна нулю:

\[F^{\text{с}} + F^{\text{у}} + F^{\text{в}} + F^{\text{и}} = 0,\quad F^{\text{и}} = -ma \qquad (1.3)\]

На рис. 1.3 показаны положительные направления действующих сил.

Рис. 1.3 — Направление сил

Подстановка соотношений для сил (раздел 1.2) в формулу равновесия (1.3):

\[kx + \mu \nu |\nu| + ma - Q \sin \frac{2\pi}{T} t = 0 \qquad (1.4)\]

Учитывая, что

\[\nu = \frac{dx}{dt}, \qquad a = \frac{d^{2}x}{dt^{2}}, \qquad (1.5)\]

получаем дифференциальное уравнение движения тела:

\[kx + \mu \frac{dx}{dt}\left|\frac{dx}{dt}\right| + m\frac{d^2x}{dt^2} - Q \sin \frac{2\pi}{T} t = 0 \qquad (1.6)\]

2. Численное интегрирование

2.1 Дискретизация времени

Использование численного подхода к интегрированию уравнения (1.5) предполагает нахождение приближенного решения для определенных моментов времени. Временная ось представляется совокупностью точек \(t_0, t_1, t_2, \dots, t_i, t_{i+1}, \dots, t_n\) (рис. 2.1), в каждой из которых ищется приближенное решение.

Рис. 2.1 — Дискретизация временной оси

Интегрирование осуществляется последовательно; выбор величины очередного шага \(\Delta t_i\) зависит как от требуемых показателей точности, так и от результатов интегрирования на предыдущих шагах.

Таким образом, использование численного подхода позволяет перейти от непрерывных значений \(x\), \(\nu\), \(a\) на всем промежутке времени от \(t_0\) до \(t_{\text{конеч}}\) к совокупности дискретных значений \(x_i\), \(\nu_i\), \(a_i\) для определенных моментов времени \(t_i\). При этом алгебраические формулы выбранного метода интегрирования заменяют дифференциальные соотношения (1.4).

2.2 Неявная одношаговая схема интегрирования уравнений движения

В вычислительном ядре PRADIS для интегрирования уравнений движения используется семейство неявных одношаговых методов. Общая идея таких методов заключается в том, что значения переменных в новый момент времени \(t_i\) вычисляются с использованием информации о состоянии системы в тот же момент \(t_i\), что и делает схему неявной.

В отличие от явных методов, где новое состояние вычисляется непосредственно по известным значениям с предыдущего шага, в неявных методах интересующие нас величины \(x_i\), \(\nu_i\), \(a_i\) оказываются связанными между собой и с ранее известными значениями \(x_{i-1}\), \(\nu_{i-1}\), \(a_{i-1}\) через систему алгебраических соотношений. Эти соотношения зависят от выбранного метода интегрирования и величины шага \(\Delta t_i\).

В качестве основной схемы в настоящем руководстве рассматривается метод Штермера. Этот метод является одним из наиболее распространённых в задачах структурной динамики благодаря своей устойчивости и хорошим энергетическим свойствам. Для системы с одной степенью свободы метод Штермера задаёт связь между переменными на новом и предыдущем шагах следующим образом:

\[\begin{split}\begin{cases} x_i = x_{i-1} + \nu_{i-1}\Delta t_i + a_i \dfrac{\Delta t_i^2}{2} \\[10pt] \nu_i = \nu_{i-1} + a_i \Delta t_i \end{cases} \qquad \text{(2.1)}\end{split}\]

Метод является неявным, так как ускорение \(a_i\) в момент \(t_i\) неизвестно заранее и входит в правую часть обоих уравнений. Подстановка (2.1) в уравнение движения приводит к нелинейному алгебраическому уравнению относительно \(a_i\) (или, что эквивалентно, относительно \(\nu_i\)), которое решается итерационно.

В состав вычислительного ядра PRADIS также включены и другие неявные разностные схемы, такие как неявный метод Эйлера, метод трапеций и метод Ньюмарка. Их применение может быть предпочтительным в зависимости от специфики решаемой задачи (например, для систем с сильным демпфированием или для задач, где важна сохранность энергии). Справочные формулы для этих методов, а также соответствующие коэффициенты приведения производных, приведены в Приложении А.

2.3 Реализация нелинейных систем

Для любого из перечисленных методов интегрирования последовательность действий при численном решении остаётся единой:

1. Уравнение равновесия в момент \(t_i\) дополняется выбранными формулами интегрирования, связывающими \(x_i\), \(\nu_i\) и \(a_i\) с известными значениями с предыдущего шага.

2. С использованием этих формул все неизвестные переменные выражаются через одну базисную переменную (как правило, это скорость \(\nu_i\) или ускорение \(a_i\)). В результате формируется одно нелинейное алгебраическое уравнение относительно выбранной базисной переменной.

3. Полученное уравнение решается итерационно (например, методом Ньютона–Рафсона) до достижения сходимости с заданной точностью. На каждой итерации вычисляются значение функции (невязка уравнения) и её производная по базисной переменной.

4. После нахождения базисной переменной по формулам интегрирования восстанавливаются остальные кинематические характеристики (\(x_i\) и \(a_i\)), после чего выполняется оценка локальной погрешности и, при необходимости, корректировка шага интегрирования.

В дальнейшем при подробном численном разборе в качестве примера выбран метод Штермера. Для остальных методов, перечисленных в Приложении А, порядок действий и логика работы вычислительного ядра полностью аналогичны; отличие заключается лишь в конкретном виде формул связи (2.1) и, как следствие, в коэффициентах приведения производных.

2.4 Система уравнений для одного шага

Определим значения \(x_1\), \(\nu_1\), \(a_1\) для момента времени \(t_1 = t_0 + \Delta t_1\).

Уравнение (1.4) для момента времени \(t_i\) имеет вид:

\[k x_i + \mu \nu_i \left| \nu_i \right| + m a_i - Q \sin \frac{2\pi}{T} t_i = 0 \qquad (2.2)\]

Дополняя это соотношение формулами метода интегрирования (2.1), получаем для момента времени \(t_1\) замкнутую систему уравнений:

\[\begin{split}\begin{cases} k x_1 + \mu \nu_1 \left| \nu_1 \right| + m a_1 - Q \sin \dfrac{2\pi}{T} t_1 = 0 \\[10pt] x_1 = x_0 + \nu_0 \Delta t_1 + a_1 \dfrac{\Delta t_1^{2}}{2} \\[10pt] \nu_1 = \nu_0 + a_1 \Delta t_1 \end{cases} \qquad (2.3)\end{split}\]

2.5 Сведение к одному уравнению относительно скорости

Приведем полученную систему (2.3) к одному уравнению, выразив неизвестные \(x_1\) и \(a_1\) через \(\nu_1\) с помощью формул метода Штермера:

\[\begin{split}\begin{cases} a_1 = \dfrac{\nu_1 - \nu_0}{\Delta t_1} \\[10pt] x_1 = x_0 + \dfrac{\nu_0 + \nu_1}{2} \Delta t_1 \end{cases} \qquad (2.4)\end{split}\]

Подставляя выражения для \(a_1\) и \(x_1\) из (2.1) в исходное уравнение движения (1.4), получаем:

\[k \left( x_0 + \frac{\nu_0 + \nu_1}{2} \Delta t_1 \right) + \mu \nu_1 \left| \nu_1 \right| + m \frac{\nu_1 - \nu_0}{\Delta t_1} - Q \sin \frac{2\pi}{T} t_1 = 0\]

Сгруппируем слагаемые относительно \(\nu_1\) и \(|\nu_1|\):

\[\mu \nu_1 |\nu_1| + \left( \frac{k \Delta t_1}{2} + \frac{m}{\Delta t_1} \right) \nu_1 + \left( k x_0 + \frac{k \nu_0 \Delta t_1}{2} - \frac{m \nu_0}{\Delta t_1} - Q \sin \frac{2\pi}{T} t_1 \right) = 0\]

Введем обозначения:

\[\alpha = \mu, \qquad \beta = \frac{k \Delta t_1}{2} + \frac{m}{\Delta t_1}, \qquad \gamma = k x_0 + \frac{k \nu_0 \Delta t_1}{2} - \frac{m \nu_0}{\Delta t_1} - Q \sin \frac{2\pi}{T} t_1 \qquad (2.5)\]

В результате дискретизации получено нелинейное алгебраическое уравнение относительно неизвестной скорости \(\nu_1\) на первом шаге интегрирования, которое назовём уравнением скорости:

\[\alpha \nu_1 |\nu_1| + \beta \nu_1 + \gamma = 0 \qquad (2.6)\]

Уравнение скорости (2.6) сохраняет свою форму для любого момента времени \(t_i\) при замене индексов (\(1 \to i\), \(0 \to i-1\)), при этом коэффициенты alpha, beta, gamma зависят от шага и известных с предыдущего шага величин.

Таким образом, применение описанного метода интегрирования сводит исходное дифференциальное уравнение (2.2) к последовательности уравнений скорости вида (2.6), решаемых на каждом шаге по времени.

3. Решение нелинейного уравнения методом Ньютона

Перейдем к решению уравнения (2.6). Для этого используем метод Ньютона.

3.1 Постановка задачи и геометрическая интерпретация

Рассматривается уравнение вида:

\[f(z) = 0 \qquad (3.1)\]

где \(f(z)\) — нелинейная функция относительно неизвестного \(z\). Корнем уравнения (3.1) называется такое значение \(z = z^*\), при котором функция \(f(z)\) обращается в ноль.

Геометрически корень \(z^*\) соответствует точке пересечения графика функции \(f(z)\) с горизонтальной осью \(z\). Таким образом, решение уравнения (3.1) сводится к отысканию абсциссы точки пересечения кривой \(f(z)\) с осью \(z\).

Для численного нахождения этого корня применяется метод Ньютона. На каждой \(j\)-й итерации метода решение заменяется отысканием точки пересечения касательной к кривой \(f(z)\) с осью \(z\); касательная строится для \(z = z^{j-1}\).

На рисунке 3.1 приведена геометрическая интерпретация описанного выше итерационного процесса.

Рис. 3.1 — Геометрическая интерпретация метода Ньютона

3.2 Алгоритм метода Ньютона

Алгоритм численного решения включает следующие шаги:

1. Выбор начального приближения к решению — величины \(z^0\).

2. Организация последовательности итераций, на каждой из которых уточняется полученное на предыдущей итерации значение \(z\) по схеме (рис. 3.1, а):

   \[\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 (3.2)\end{split}\]

   где \(f(z^{j-1})\) — значение функции \(f(z)\) при \(z = z^{j-1}\), \(f'(z^{j-1})\) — значение производной \(df/dz\) при \(z = z^{j-1}\).

3. Проверка на каждой итерации условия прекращения итераций (рис. 3.1, б):

   \[|\Delta z^{j}| \leq \delta_z,\qquad |f(z^{j})| \leq \delta_f \qquad (3.3)\]

   где \(\delta_f\) — допустимая невязка (отклонение от нуля) правой части уравнения (3.1); \(\delta_z\) — допустимая величина отличия решения на двух соседних итерациях.

4. Проверка ограничения на максимально допустимое количество итераций:

   \[j \leq j_{\text{max}} \qquad (3.4)\]

3.3 Пример численного расчёта (Шаг 1)

Выполним первый шаг интегрирования — от \(t_0\) до \(t_1\). Все обозначения соответствуют \(i = 1\).

Применим описанный выше метод Ньютона для решения уравнения (2.6). Для этого перепишем его в виде (3.1), приняв \(z = \nu_1\):

\[f(z) = \alpha z|z| + \beta z + \gamma \qquad (3.5)\]

Производная функции (3.5):

\[f'(z) = 2\alpha|z| + \beta \qquad (3.6)\]

Коэффициенты \(\alpha, \beta, \gamma\) определяются через параметры механической системы согласно (2.1) и (2.6).

Исходные данные для шага 1:

\[k = 20000\ \text{Н/м},\quad \mu = 1000\ \text{Н·с}^2/\text{м}^2,\quad m = 0.1\ \text{кг},\quad Q = 1000,\quad T = 0.2\pi\]

\[x_0 = 0,\quad \nu_0 = 0,\quad \Delta t_1 = 0.001\ \text{с}\]

Тогда, согласно (2.1):

\[\alpha = \mu = 1000\]

\[\beta = \frac{k \Delta t_1}{2} + \frac{m}{\Delta t_1} = \frac{20000 \cdot 0.001}{2} + \frac{0.1}{0.001} = 10 + 100 = 110\]

\[\gamma = k x_0 + \frac{k \nu_0 \Delta t_1}{2} - \frac{m \nu_0}{\Delta t_1} - Q \sin \frac{2\pi}{T} t_1\]

Подставляя значения:

\[t_1 = \Delta t_1 = 0.001\ \text{с},\qquad \frac{2\pi}{T} = \frac{2\pi}{0.2\pi} = 10,\qquad \frac{2\pi}{T} t_1 = 10 \cdot 0.001 = 0.01\ \text{рад}\]

Таким образом:

\[\gamma = 20000 \cdot 0 + \frac{20000 \cdot 0 \cdot 0.001}{2} - \frac{0.1 \cdot 0}{0.001} - 1000 \sin(0.01) \approx -10\]

Итак, коэффициенты уравнения скорости:

\[\alpha = 1000,\quad \beta = 110,\quad \gamma = -10\]

Следовательно:

\[f(z) = 1000\, z\, |z| + 110\, z - 10,\qquad f'(z) = 2000\, |z| + 110\]

Зададим параметры метода Ньютона:

\[\delta_z = 0.001,\quad \delta_f = 0.1,\quad j_{\text{max}} = 5,\quad z^0 = 0\]

Первая итерация.

\[\Delta z^1 = -\frac{f(z^0)}{f'(z^0)} = -\frac{1000 \cdot 0 \cdot 0 + 110 \cdot 0 - 10}{2000 \cdot |0| + 110} = \frac{10}{110} = 0.09091\]

\[z^1 = 0 + 0.09091 = 0.09091\]

Проверка критериев остановки:

\[|\Delta z^1| = 0.09091 > \delta_z,\qquad f(z^1) = 1000 \cdot 0.09091 \cdot 0.09091 + 110 \cdot 0.09091 - 10 = 8.26 > \delta_f\]

Продолжаем итерации.

Вторая итерация.

\[\Delta z^2 = -\frac{f(z^1)}{f'(z^1)} = -\frac{8.26}{2000 \cdot 0.09091 + 110} = -0.02831\]

\[z^2 = 0.09091 - 0.02831 = 0.06260\]

Проверка: \(|\Delta z^2| = 0.02831 > \delta_z\), \(|f(z^2)| = 0.80 > \delta_f\).

Третья итерация.

\[\Delta z^3 = -0.00341,\quad z^3 = 0.05918\]

\[|\Delta z^3| = 0.00341 > \delta_z,\qquad |f(z^3)| = 0.012 < \delta_f\]

Четвёртая итерация.

\[\Delta z^4 = -0.00005,\quad z^4 = 0.05913\]

\[|\Delta z^4| = 0.00005 < \delta_z,\qquad |f(z^4)| = 0.0006 < \delta_f\]

Оба условия (3.3) удовлетворены, ограничение (3.4) не превышено. Решение достигнуто. Скорость для момента времени \(t_1\):

\[\nu_1 = 0.05913\ \text{м/с}\]

Используя формулы (2.1), определим ускорение и перемещение:

\[a_1 = \frac{\nu_1 - \nu_0}{\Delta t_1} = \frac{0.05913 - 0}{0.001} = 59.13\ \text{м/с}^2\]

\[x_1 = x_0 + \frac{\nu_0 + \nu_1}{2} \Delta t_1 = 0 + \frac{0 + 0.05913}{2} \cdot 0.001 = 2.96 \times 10^{-5}\ \text{м}\]

Результат шага 1:

\[x_1 = 2.96 \times 10^{-5}\ \text{м},\quad \nu_1 = 0.05913\ \text{м/с},\quad a_1 = 59.13\ \text{м/с}^2\]

Примечание

Полученные значения \(x_1\), \(\nu_1\), \(a_1\) будут использованы как начальные условия для шага 2 (раздел 4.3).

4. Управление шагом интегрирования

При использовании описанного выше метода интегрирования (раздел 2) возникает необходимость автоматического выбора шага по времени для обеспечения заданной точности. Для этого используется оценка локальной погрешности.

4.1 Локальная погрешность

Погрешность, возникающую при интегрировании на текущем \(i\)-м шаге (её называют локальной), мы будем оценивать по формуле:

\[lp_i = \left| \frac{\nu_i^p - \nu_i^c}{2} \right| \qquad (4.1)\]

Расшифруем новые обозначения:

1. \(\nu_i^p\) — это прогноз (предсказание) скорости, полученный по явной формуле, знакомой нам по выбору начального приближения (4.1):

   \[\nu_i^p = \nu_{i-1} + a_{i-1} \Delta t_i \qquad (4.2)\]

2. \(\nu_i^c\) — это скорректированное значение скорости, которое находится из неявной формулы с помощью итераций (например, методом Ньютона):

   \[\nu_i^c = \nu_{i-1} + a_i \Delta t_i \qquad (4.3)\]

Обратимся к рис. 3.2, который поясняет эту оценку.

Рис. 4.4

Допустим, в момент \(t_{i-1}\) мы знаем точное значение \(\nu_{i-1}\). Если мы воспользуемся явной формулой (4.2), мы просто продолжим движение по касательной, проведенной к графику \(\nu(t)\) в точке \(t_{i-1}\). Угол наклона этой касательной задаётся величиной \(a_{i-1}\).

Для расчёта по неявной формуле (4.3) нам понадобится значение \(a_i\) — то есть наклон касательной уже в конечной точке шага \(t_i\), который нам пока неизвестен. Поэтому мы решаем систему уравнений (2.1) и (4.3) итеративно.

Как видно из рисунка, истинная кривая \(\nu(t)\) проходит между двумя полученными точками: явным прогнозом \(\nu_i^p\) и уточнённым решением \(\nu_i^c\). Разность между ними показывает, насколько сильно график скорости изогнут на этом шаге. Чем больше эта разность, тем выше ошибка интегрирования. Уменьшая шаг \(\Delta t_i\), мы уменьшаем и локальную погрешность, рассчитанную по формуле (4.1).

4.2 Правило выбора шага

Если \(lp_i \leq \delta_l\) — шаг успешен. Если \(lp_i > \delta_l\) — шаг повторяется с уменьшенным \(\Delta t_i\).

Рекомендуемое значение шага:

\[\Delta t_{i+1} = c\,\Delta t_i \sqrt{\frac{\delta_l}{lp_i}},\quad c < 1 \qquad (4.4)\]

4.3 Пример расчёта с управлением шага (Шаг 2)

Выполним второй шаг интегрирования — от \(t_1\) до \(t_2\). Начальные условия — результат шага 1 (раздел 3.3).

Исходные данные (после шага 1):

\[x_1 = 2.96 \times 10^{-5}\ \text{м},\quad \nu_1 = 0.05913\ \text{м/с},\quad a_1 = 59.13\ \text{м/с}^2\]

\[\Delta t_2 = 0.001\ \text{с},\quad \delta_l = 0.001\]

1. Прогноз скорости (явная формула):

\[\nu_2^p = \nu_1 + a_1 \Delta t_2 = 0.05913 + 59.13 \cdot 0.001 = 0.11826\ \text{м/с}\]

2. Итерационное уточнение (неявная схема):

После решения методом Ньютона (аналогично шагу 1) получаем:

\[\nu_2^c = 0.11159\ \text{м/с}\]

3. Оценка локальной погрешности (см. раздел 4.1):

\[lp_2 = \left| \dfrac{\nu_2^p - \nu_2^c}{2} \right| = \left| \dfrac{0.11826 - 0.11159}{2} \right| = 0.00333\]

4. Сравнение с допустимой погрешностью:

\[lp_2 = 0.00333 > \delta_l = 0.001 \quad \Rightarrow\quad \text{шаг необходимо уменьшить}\]

5. Расчёт нового шага (см. раздел 4.2, \(c = 0.8\)):

\[\Delta t_2^{\text{(новый)}} = 0.8 \cdot 0.001 \cdot \sqrt{\frac{0.001}{0.00333}} = 0.438 \times 10^{-3}\ \text{с}\]

6. Повторный расчёт шага 2 с уменьшенным шагом:

После повторного решения методом Ньютона:

\[\nu_2 = 0.08509\ \text{м/с},\qquad a_2 = \frac{\nu_2 - \nu_1}{\Delta t_2} = 59.21\ \text{м/с}^2,\qquad x_2 = x_1 + \frac{\nu_1 + \nu_2}{2} \Delta t_2 = 6.12 \times 10^{-5}\ \text{м}\]

\[lp_2 = 0.00064 \leq \delta_l \quad \Rightarrow\quad \text{шаг принят}\]

Результат шага 2:

\[x_2 = 6.12 \times 10^{-5}\ \text{м},\quad \nu_2 = 0.08509\ \text{м/с},\quad a_2 = 59.21\ \text{м/с}^2,\quad lp_2 = 0.00064\]

Примечание

Эти данные будут использованы как начальные условия для шага 3 (раздел 7).

4.4 Основные этапы численного решения

Алгоритм численного интегрирования может быть представлен в виде следующей последовательности действий:

1. Формирование дифференциального уравнения, описывающего поведение системы:

   \[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\]

   При формировании уравнения использован второй закон Ньютона.

2. Дискретизация — замена производных конечно-разностными соотношениями на сетке \(t_i\). При использовании неявной схемы второго порядка уравнения связи записываются непосредственно для текущего шага:

   \[\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}\]

   Тем самым исходная задача сводится к последовательности алгебраических систем.

3. Исключение переменных — приведение системы к одному нелинейному уравнению относительно \(\nu_i\):

   \[\alpha \nu_i |\nu_i| + \beta \nu_i + \gamma = 0\]

4. Решение нелинейного уравнения методом Ньютона с явным прогнозом в качестве начального приближения.

5. Адаптивный контроль шага:

   • оценка локальной погрешности интегрирования;

   • при превышении допустимой погрешности — повтор шага с уменьшенным \(\Delta t\).

6. Обновление переменных — вычисление \(a_i\) и \(x_i\) по найденному \(\nu_i\) с использованием уравнений связи.

7. Коррекция шага для следующей итерации на основе отношения допустимой и фактической погрешности.

5. Автоматическое формирование математической модели

5.1 Автоматическое построение нелинейного уравнения и его линеаризация

На рассмотренном примере мы обозначили канву численного решения, которой придерживается алгоритм вычислительного ядра PRADIS. Однако формирование дифференциального уравнения проводилось «вручную», а некоторые вопросы решались неформально (например, аналитическое определение производной \(df(z)/dz\)). Поэтому продолжим рассмотрение методов PRADIS с выяснения принципов автоматического формирования системы дифференциальных уравнений.

Вернемся к системе на рис. 1.1. Уравнение равновесия:

\[F^c + F^y + F^в + F^u = 0 \qquad (5.1)\]

Для каждого \(i\)-го момента времени уравнение может быть записано в виде:

\[F_i^c + F_i^y + F_i^в + F_i^u = 0 \qquad (5.2)\]

где

\[\begin{split}\begin{cases} F_i^c = -Q \sin \dfrac{2\pi}{T} t \\[8pt] F_i^y = k x_i \\[8pt] F_i^в = \mu \nu_i |\nu_i| \\[8pt] F_i^u = m a_i \end{cases} \qquad (5.3)\end{split}\]

(Обратите внимание: знаки отличаются от (1.1), поскольку в PRADIS при рассмотрении условий равновесия суммируются усилия, действующие со стороны системы на элементы.)

Связь между \(x_i\), \(\nu_i\), \(a_i\) задается формулами метода интегрирования (1.7), которые воспроизведем еще раз:

\[x_i = x_{i-1} + \nu_{i-1} \Delta t_i + a_i \frac{\Delta t_i^2}{2} \qquad (5.4)\]

\[\nu_i = \nu_{i-1} + a_i \Delta t_i \qquad (5.5)\]

Соотношение (5.2) представляет собой алгебраическое нелинейное уравнение, которое можно решать методом Ньютона. Имеем:

\[f(z) = 0, \qquad (5.6)\]

где

\[f(z) = F_i^c + F_i^y + F_i^в + F_i^u\]

В качестве переменной \(z\) может быть выбрана любая из компонент \(x_i\), \(\nu_i\) или \(a_i\), поскольку они связаны соотношениями (5.4). Примем, как и ранее, \(z = \nu_i\).

На каждой итерации необходимо вычислить \(f(z)\) и \(df(z)/dz\).

Вычисление \(f(z)\) сводится к суммированию текущих значений сил при текущих значениях \(x_i\), \(\nu_i\), \(a_i\) (\(i\) — номер шага по времени, \(j\) — номер итерации по Ньютону). Вычисляемое значение \(f(z)\) представляет собой погрешность выполнения условия равновесия, которую необходимо уменьшить до допустимых пределов.

Производная \(df(z)/dz\):

\[\frac{df(z)}{dz} = \frac{dF_i^c}{dz} + \frac{dF_i^y}{dz} + \frac{dF_i^в}{dz} + \frac{dF_i^u}{dz} \qquad (5.7)\]

Каждую из производных следует представить как производную сложной функции. Для силы \(F_i^c\):

\[\frac{dF_i^c}{dz} = \frac{\partial F_i^c}{\partial x_i} \frac{dx_i}{dz} + \frac{\partial F_i^c}{\partial \nu_i} \frac{d\nu_i}{dz} + \frac{\partial F_i^c}{\partial a_i} \frac{da_i}{dz}\]

Поскольку \(z = \nu_i\):

\[\frac{dF_i^c}{dz} = \frac{\partial F_i^c}{\partial x_i} \frac{dx_i}{d\nu_i} + \frac{\partial F_i^c}{\partial \nu_i} \frac{d\nu_i}{d\nu_i} + \frac{\partial F_i^c}{\partial a_i} \frac{da_i}{d\nu_i} \qquad (5.8)\]

Аналогично для остальных сил (5.9):

\[\frac{dF_i^y}{dz} = \frac{\partial F_i^y}{\partial x_i} \frac{dx_i}{d\nu_i} + \frac{\partial F_i^y}{\partial \nu_i} \frac{d\nu_i}{d\nu_i} + \frac{\partial F_i^y}{\partial a_i} \frac{da_i}{d\nu_i}\]

\[\frac{dF_i^в}{dz} = \frac{\partial F_i^в}{\partial x_i} \frac{dx_i}{d\nu_i} + \frac{\partial F_i^в}{\partial \nu_i} \frac{d\nu_i}{d\nu_i} + \frac{\partial F_i^в}{\partial a_i} \frac{da_i}{d\nu_i}\]

\[\frac{dF_i^u}{dz} = \frac{\partial F_i^u}{\partial x_i} \frac{dx_i}{d\nu_i} + \frac{\partial F_i^u}{\partial \nu_i} \frac{d\nu_i}{d\nu_i} + \frac{\partial F_i^u}{\partial a_i} \frac{da_i}{d\nu_i}\]

Используем уравнения связи (5.4) для получения зависимостей \(a_i\) и \(x_i\) от \(\nu_i\):

\[a_i = \frac{\nu_i - \nu_{i-1}}{\Delta t_i}\]

\[x_i = x_{i-1} + \frac{\nu_i + \nu_{i-1}}{2} \Delta t_i \qquad (5.10)\]

Дифференцируем (5.10) по \(\nu_i\):

\[\frac{dx_i}{d\nu_i} = \frac{\Delta t_i}{2},\qquad \frac{d\nu_i}{d\nu_i} = 1,\qquad \frac{da_i}{d\nu_i} = \frac{1}{\Delta t_i} \qquad (5.11)\]

5.2 Вычисление частных производных

С использованием зависимостей (5.3) вычислим частные производные, входящие в выражения (5.8) и (5.9).

Для силы \(F_i^c\):

\[\frac{\partial F_i^c}{\partial x_i} = 0,\quad \frac{\partial F_i^c}{\partial \nu_i} = 0,\quad \frac{\partial F_i^c}{\partial a_i} = 0, \qquad (5.12)\]

поскольку \(F_i^c\) не зависит от перемещения, скорости и ускорения тела (5.12).

Для силы \(F_i^y\):

\[\frac{\partial F_i^y}{\partial x_i} = k,\quad \frac{\partial F_i^y}{\partial \nu_i} = 0,\quad \frac{\partial F_i^y}{\partial a_i} = 0 \qquad (5.13)\]

Для силы \(F_i^в\):

\[\frac{\partial F_i^в}{\partial x_i} = 0,\quad \frac{\partial F_i^в}{\partial \nu_i} = 2\mu |\nu_i|,\quad \frac{\partial F_i^в}{\partial a_i} = 0 \qquad (5.14)\]

Для силы \(F_i^u\):

\[\frac{\partial F_i^u}{\partial x_i} = 0,\quad \frac{\partial F_i^u}{\partial \nu_i} = 0,\quad \frac{\partial F_i^u}{\partial a_i} = m \qquad (5.15)\]

Подставляя найденные значения частных производных и коэффициенты из (5.11) в формулы (5.8) и (5.9), получаем:

\[\frac{\partial F_i^c}{\partial z} = 0,\quad \frac{\partial F_i^y}{\partial z} = \frac{k \Delta t}{2},\quad \frac{\partial F_i^в}{\partial z} = 2\mu |\nu_i|,\quad \frac{\partial F_i^u}{\partial z} = \frac{m}{\Delta t_i} \qquad (5.16)\]

Суммируя слагаемые в формуле (5.7), приходим к выражению:

\[\frac{df(z)}{dz} = \frac{k \Delta t_i}{2} + 2\mu |\nu_i| + \frac{m}{\Delta t_i} \qquad (5.17)\]

Полученный результат структурно аналогичен ранее выведенной формуле (3.6), коэффициенты которой определяются соотношениями (2.6).

5.3 Преимущества автоматического подхода

Теперь, имея возможность вычислять \(f(z)\) и \(df(z)/dz\), продолжить расчет по описанному алгоритму не представляет труда. Какими же преимуществами окупаются дополнительные выкладки?

1. Дифференциальное уравнение движения выписывать не пришлось.

2. Развернутая форма нелинейного уравнения вида (2.6) оказалась невостребованной.

3. Функциональная зависимость для производной \(df(z)/dz\) не потребовалась.

Рис. 5.1

Для формирования математической модели была использована следующая информация (далее — «перечень»):

1. Сведения о стыковке элементов схемы (рис. 5.1).

2. Условие равновесия сил для \(i\)-го момента времени — нелинейное алгебраическое уравнение относительно \(x_i\), \(\nu_i\), \(a_i\):

   \[\sum_{k=1}^{4} F_i^{(k)} = 0 \qquad (5.18)\]

   (топологическое уравнение, определяемое структурой связей).

3. Компонентные уравнения — выражения, позволяющие определить усилия в каждом элементе как функции перемещения, скорости, ускорения и времени (5.3).

4. Частные производные усилий по перемещению, скорости, ускорению (5.12)–(5.15).

5. Алгебраические уравнения связи \(x\), \(\nu\), \(a\) — формулы метода интегрирования (5.4).

6. Выбор базисной переменной \(z\) (из \(x\), \(\nu\) или \(a\)) для решения нелинейного уравнения.

Этой информации достаточно для реализации машинных алгоритмов формирования математической модели объекта.

6. Модели элементов

В качестве объекта моделирования рассматривается техническая система, расчетная схема которой изображена на рис. 1.1.

В программном комплексе Pradis модель системы строится в виде набора элементов, объединяемых в узлах. Набор степеней свободы, приписываемых элементу, зависит от его типа (рис. 6.1).

Рис. 6.1

Работа строится на использовании моделей элементов. Из них, как из деталей конструктора, пользователь формирует общую модель системы. Его задача — обеспечить корректную стыковку компонентов; формирование математического описания полностью берет на себя ПО.

Для построения расчетной схемы, приведенной на рис. 1.1, используются готовые модели элементов из библиотеки base («Базовые компоненты»).

Рис. 6.2

Этим описанием пользователь сообщает программному комплексу PRADIS всю необходимую информацию по стыковке элементов схемы (пункт 1 перечня п.5.3). Выбраны модели K, MUNL, M, FSIN; они соединены надлежащим образом; узел 1 описан как неподвижный.

В процессе обработки описания структуры определяется размерность системы уравнений — количество узлов, в которых должны выполняться условия равновесия. В рассматриваемом примере два узла, один из которых закреплен. На этапе формирования математической модели структура данных готовится по обоим узлам, но на этапе расчета уравнение для закрепленного узла исключается, а его кинематические характеристики обнуляются.

Этап численного интегрирования представляет собой последовательность шагов по времени, каждый из которых сводится к решению нелинейного уравнения равновесия вида (5.18). Для этого необходима информация, приведенная в пунктах 3–6 перечня (раздел 5.3).

Входными данными для любой модели элемента являются:

• неизменный список параметров модели;

• текущие значения перемещений, скоростей и ускорений узлов, с которыми элемент соединен.

По этим данным модель элемента обязана для текущего момента времени вычислить:

1. усилия, действующие со стороны системы на элементы (вектор усилий элемента) — пункт 3 перечня (раздел 5.3);

2. частные производные усилий по перемещениям, скоростям и ускорениям узлов элемента (матрицу Якоби элемента) — пункт 4 перечня (раздел 5.3).

Если элемент имеет \(N\) степеней свободы, длина вектора усилий элемента равна \(N\), а якобиан элемента имеет размер \(N \times N \times 3\).

Рис. 6.3

Например, разработчик двухузловой модели одномерной упругой пружины реализовал следующие зависимости (рис. 6.3):

\[F_1 = k(x_1 - x_2),\quad F_2 = k(x_2 - x_1) \qquad (6.1)\]

\[\frac{\partial F_1}{\partial x_1} = k,\; \frac{\partial F_1}{\partial x_2} = -k,\; \frac{\partial F_2}{\partial x_1} = -k,\; \frac{\partial F_2}{\partial x_2} = k \qquad (6.2)\]

\[\frac{\partial F_1}{\partial \nu_1} = \frac{\partial F_1}{\partial \nu_2} = \frac{\partial F_2}{\partial \nu_1} = \frac{\partial F_2}{\partial \nu_2} = 0 \qquad (6.3)\]

\[\frac{\partial F_1}{\partial a_1} = \frac{\partial F_1}{\partial a_2} = \frac{\partial F_2}{\partial a_1} = \frac{\partial F_2}{\partial a_2} = 0 \qquad (6.4)\]

Для любого момента времени модель элемента по переданным значениям перемещений, скоростей и ускорений вычисляет усилия на концах пружины и частные производные. Вектор усилий состоит из двух элементов, якобиан — из 12.

Поскольку узел 1 пружины закреплен, востребуется только информация, связанная с незакрепленным вторым узлом:

\[F_2 = k(x_2 - x_1) \qquad (6.5)\]

\[\frac{\partial F_2}{\partial x_2} = k,\quad \frac{\partial F_2}{\partial \nu_2} = \frac{\partial F_2}{\partial a_2} = 0 \qquad (6.6)\]

Эта информация позволяет учесть вклад пружины при решении нелинейного уравнения (5.6). Вклад остальных элементов учитывается аналогично.

7. Алгоритм выполнения шага интегрирования

Продолжим ранее начатое интегрирование, сделав третий шаг по времени, используя формальный алгоритм, базирующийся на формулах (5.6)–(5.17).

По результатам второго шага (\(t_2 = 1.438 \times 10^{-3}\) с):

\[x_2 = 6.12 \times 10^{-5}\ \text{м},\quad \nu_2 = 0.08509\ \text{м/с},\quad a_2 = 59.21\ \text{м/с}^2\]

Величина шага \(\Delta t_2 = 0.438 \times 10^{-3}\) с, локальная погрешность \(lp_2 = 0.000018\).

Рекомендуемое значение \(\Delta t_3\) определяем по формуле (4.4) с \(c = 0.8\) и \(\delta_l = 0.001\):

\[\Delta t_3 = 0.8 \cdot 0.438 \times 10^{-3} \cdot \sqrt{\frac{0.001}{0.000018}} = 2.63 \times 10^{-3}\ \text{с}\]

Действуем по схеме, представленной на рис. 7.1–7.2.

Рис. 7.1

Рис. 7.2

Перед третьим шагом известны значения \(t_2\), \(x_2\), \(\nu_2\), \(a_2\), \(\Delta t_3\).

1. Определение коэффициентов приведения якобиана (5.11):

   \[\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\]

2. Начальное приближение (явный прогноз по формуле (4.2)):

   \[\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} = 48.97 \times 10^{-5}\]

3. Первая итерация Ньютона (см. рис. 7.1в).

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\]

5. Суммирование сил (значение \(f(z)\) на текущей итерации):

   \[\sum F = F^y + F^в + F^u + F^c = -40.6 + 9.3 + 58.0 + 5.9 = 32.6\]

6. Вычисление \(df(z)/dz\) по формулам (5.7), (5.8), (5.9):

   \[\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\]

7. Приращение \(\Delta z^1\):

   \[\Delta z^1 = -\frac{f(z^0)}{f'(z^0)} = -\frac{32.6}{545.8} = -0.05973\]

8. Очередное приближение:

   \[z^1 = z^0 + \Delta z^1 = 0.24081 - 0.05973 = 0.18108\]

9. Уточнение \(x_3^1\), \(\nu_3^1\), \(a_3^1\) по формулам (5.10):

   \[\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}\]

10. Проверка условий завершения итераций:

    \[|f(z^0)| = 32.6 > \delta_f = 0.1,\quad |\Delta z^1| = 0.05973 > \delta_z = 0.001\]

    Итерации необходимо продолжить.

11. Проверка лимита итераций:

    \[j = 1 < j_{\text{max}} = 5\]

12. Увеличение счетчика итераций: \(j = 2\).

13. Вторая итерация дает:

    \[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\]

    Условия завершения не выполнены.

14. Третья итерация (последняя):

    \[|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\]

15. Оценка локальной погрешности:

    \[lp_3 = \left| \frac{\nu_3^p - \nu_3^c}{2} \right| = \left| \frac{0.24081 - 0.17232}{2} \right| = 0.034\]

16. Корректировка правила выбора шага. Практика расчетов показала, что формула (4.4) приемлема только при соотношении \(\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 (7.1)\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{с}\]

17. Повторный расчет с \(\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}\) с.

Расчет на третьем шаге завершен.

8. Разделение функций в вычислительном ядре

Основное, на что следует обратить внимание, — это разделение функций между программой интегрирования и программами реализации моделей элементов.

Программе интегрирования, работающей по алгоритму рис. 2.4а–2.4в, не важно, какие процессы интегрируются. Ее зависимость от моделей элементов сводится только к своевременному получению векторов сил и матриц якобианов. Какие свойства отдельных элементов отражает эта информация, программу интегрирования не касается.

Модели элементов, в свою очередь, имеют свой уровень независимости при четко очерченных обязанностях перед «верхами». Физические свойства отдельного элемента отражаются в компонентных уравнениях на уровне модели элемента, а программа интегрирования работает на уровне уравнений равновесия потоков.

Такое разграничение функций определяет универсальность вычислительного ядра PRADIS — принципиальную возможность расчета любых объектов, процессы в которых подчиняются законам равновесия потоковых переменных (равновесие сил, электрических и тепловых потоков, расходов жидкости и газа).

ПРИЛОЖЕНИЕ А — формулы методов интегрирования и коэффициенты приведения якобиана

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

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

На шаге интегрирования известны значения перемещения, скорости и ускорения на предыдущем шаге:

\[x_{i-1},\ \nu_{i-1},\ a_{i-1}\]

и требуется определить значения на текущем шаге:

\[x_i,\ \nu_i,\ a_i\]

где:

• \(x\) — перемещение или обобщенная потенциальная переменная;

• \(\nu\) — первая производная потенциальной переменной по времени;

• \(a\) — вторая производная потенциальной переменной по времени;

• \(\Delta t_i = t_i - t_{i-1}\) — величина текущего шага интегрирования.

В PRADIS формулы интегрирования задают алгебраические связи между \(x_i\), \(\nu_i\), \(a_i\) и значениями, известными с предыдущего шага. Эти связи позволяют на каждом шаге интегрирования перейти от системы нелинейных ОДУ к системе нелинейных алгебраических уравнений.

А.1. Метод Штермера

Для метода Штермера используются соотношения:

\[x_i = x_{i-1} + \nu_{i-1} \Delta t_i + a_i \frac{\Delta t_i^2}{2}\]

\[\nu_i = \nu_{i-1} + a_i \Delta t_i\]

Если в качестве базисной неизвестной на текущем шаге выбрана скорость \(z = \nu_i\), то:

\[a_i = \frac{\nu_i - \nu_{i-1}}{\Delta t_i}\]

\[x_i = x_{i-1} + \frac{\nu_{i-1} + \nu_i}{2} \Delta t_i\]

Коэффициенты приведения производных:

\[C_x = \frac{dx_i}{d\nu_i} = \frac{\Delta t_i}{2}\]

\[C_\nu = \frac{d\nu_i}{d\nu_i} = 1\]

\[C_a = \frac{da_i}{d\nu_i} = \frac{1}{\Delta t_i}\]

А.2. Метод Ньюмарка

В обобщенной форме метод Ньюмарка записывается через параметры \(\beta_N\) и \(\gamma_N\):

\[x_i = x_{i-1} + \nu_{i-1} \Delta t_i + \left( \frac{1}{2} - \beta_N \right) a_{i-1} \Delta t_i^2 + \beta_N a_i \Delta t_i^2\]

\[\nu_i = \nu_{i-1} + (1 - \gamma_N) a_{i-1} \Delta t_i + \gamma_N a_i \Delta t_i\]

При выборе базисной неизвестной \(z = \nu_i\):

\[a_i = \frac{\nu_i - \nu_{i-1} - (1 - \gamma_N) a_{i-1} \Delta t_i}{\gamma_N \Delta t_i}\]

\[x_i = x_{i-1} + \nu_{i-1} \Delta t_i + \left( \frac{1}{2} - \beta_N \right) a_{i-1} \Delta t_i^2 + \beta_N a_i \Delta t_i^2\]

Коэффициенты приведения производных:

\[C_x = \frac{dx_i}{d\nu_i} = \frac{\beta_N \Delta t_i}{\gamma_N}\]

\[C_\nu = 1\]

\[C_a = \frac{da_i}{d\nu_i} = \frac{1}{\gamma_N \Delta t_i}\]

А.3. Неявный метод Эйлера

Для неявного метода Эйлера используются соотношения:

\[x_i = x_{i-1} + \nu_i \Delta t_i\]

\[\nu_i = \nu_{i-1} + a_i \Delta t_i\]

При выборе базисной неизвестной \(z = \nu_i\):

\[a_i = \frac{\nu_i - \nu_{i-1}}{\Delta t_i}\]

\[x_i = x_{i-1} + \nu_i \Delta t_i\]

Коэффициенты приведения производных:

\[C_x = \Delta t_i\]

\[C_\nu = 1\]

\[C_a = \frac{1}{\Delta t_i}\]

А.4. Метод трапеций

Для метода трапеций используются соотношения:

\[x_i = x_{i-1} + \frac{\nu_{i-1} + \nu_i}{2} \Delta t_i\]

\[\nu_i = \nu_{i-1} + \frac{a_{i-1} + a_i}{2} \Delta t_i\]

При выборе базисной неизвестной \(z = \nu_i\):

\[a_i = 2 \frac{\nu_i - \nu_{i-1}}{\Delta t_i} - a_{i-1}\]

\[x_i = x_{i-1} + \frac{\nu_{i-1} + \nu_i}{2} \Delta t_i\]

Коэффициенты приведения производных:

\[C_x = \frac{\Delta t_i}{2}\]

\[C_\nu = 1\]

\[C_a = \frac{2}{\Delta t_i}\]

А.5. Использование коэффициентов в якобиане

Модель элемента возвращает частные производные потоковой переменной \(F\) по перемещению, скорости и ускорению:

\[\frac{\partial F}{\partial x},\ \frac{\partial F}{\partial \nu},\ \frac{\partial F}{\partial a}\]

Программа интегрирования преобразует эти производные к производной по выбранной базисной неизвестной \(z\). Если \(z = \nu_i\), то:

\[\frac{dF}{dz} = \frac{\partial F}{\partial x} C_x + \frac{\partial F}{\partial \nu} C_\nu + \frac{\partial F}{\partial a} C_a\]

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

В многомерной задаче такая операция выполняется для всех активных степеней свободы и всех моделей элементов, после чего локальные вклады собираются в глобальную матрицу якобиана расчетной системы.