PRADIS

PRADIS

Методы формирования и численного расчета математических моделей переходных процессов

ПРОГРАММНЫЙ КОМПЛЕКС ДЛЯ АВТОМАТИЗАЦИИ МОДЕЛИРОВАНИЯ НЕСТАЦИОНАРНЫХ ПРОЦЕССОВ В МЕХАНИЧЕСКИХ СИСТЕМАХ И СИСТЕМАХ ИНОЙ ФИЗИЧЕСКОЙ ПРИРОДЫ

ВЕРСИЯ 4.2

Содержание

1. Простой пример. 1

2. Механизм формирования математической модели. 15

3. Кратко об угловых степенях свободы, используемых в пространственных элементах PRADIS 32

1. Простой пример

Рис. 1.1

Рассмотрим колебательную систему с одной степенью свободы. Тело массы m соединяется с неподвижным основанием посредством упругой пружины. Его движение происходит в среде с вязким сопротивлением под действием внешней силы (рис. 1.1.). Для начального момента времени заданы положение тела в пространстве и его начальная скорость. Требуется для каждого момента времени в интервале от t₀ до t_конеч определить перемещение, скорость и ускорение тела.

Рис. 1.2

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

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

(1.1)

или

, (1.2)

Рис. 1.3.

где с учетом выбранных направлений (рис.1.3):

(1.3)

Здесь x, v и a - соответственно, перемещение, скорость и ускорение тела.

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

(1.4)

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

и , (1.5)

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

(1.6)

Использование численного подхода к интегрированию уравнения (1.6) предполагает нахождение приближенного решения для определенных моментов времени, т.е., временная ось представляется совокупностью точек t₀, t₁, t₂, ... t_i, t_i+₁, ... t_n (рис. 1.4), в каждой из которых ищется приближенное решение уравнения(1.6).

Рис. 1.4.

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

Таким образом, использование численного подхода к решению уравнения (1.6) позволяет перейти от непрерывных значений x, v, a на всем промежутке времени от t=0 до t=t_конеч. к совокупности дискретных значений x_i, v_i, a_i для определенных моментов времени t_i. При этом алгебраические формулы выбранного метода интегрирования заменяют дифференциальные соотношения (1.5). Так, формулы неявного одношагового метода Штермера устанавливают следующую зависимость для переменных x_i, v_i по известным с предыдущего шага значениям x_i-1, v_i-₁ [1]:

(1.7)

где ,

t₀ - начальное время,

x₀,v₀ - начальные значения перемещения и скорости,

t_n - конечное время.

Для начального момента времени t₀ значения x₀ и v₀ должны быть известны. Оставляя пока в стороне вопрос выбора величины шага интегрирования t_i, определим значения x₁, v₁, a₁ для момента времени t₁= t₀ +D t₁.

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

(1.8)

Дополняем это соотношение формулами выбранного метода интегрирования (1.7) и получаем для момента времени t₁ замкнутую систему уравнений:

(1.9)

Приведем полученную систему к одному уравнению, выразив неизвестные x₁ и a₁ через v₁:

(1.10)

Получаем:

(1.11)

Группируя сомножители при одинаковых степенях неизвестного v₁, соотношение (1.11) можно записать в виде:

, (1.12)

где

(1.12а)

Заметим, что соотношение (1.12) сохраняет свой вид для любого момента времени ti при соответствующей замене подстрочных индексов (1 на i, 0 на i-1).

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

Уравнение (1.12) решается методом Ньютона. Позволим себе напомнить последовательность действий при решении нелинейного уравнения методом Ньютона.

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

, (1.13)

где f(z) - нелинейная функция относительно неизвестного z.

Рис. 1.5

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

1) выбор начального приближения к решению - величины z⁰;

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

(1.14)

где f(z^j-¹) - значение функции f(z) при z=z^j-¹, f '(z^j-1) - значение производной f(z)/dz при при z=z^j-1;

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

(1.15)

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

- допустимая величина отличия решения на двух соседних итерациях;

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

(1.16)

Геометрически решение уравнения (1.13) сводится к отысканию абсциссы точки пересечения с осью z кривой f(z). На каждой j-й итерации метода Ньютона решение этой задачи заменяется отысканием точки пересечения касательной к кривой f(z) с осью z, при этом касательная строится для z=z^j-¹.

Возвращаемся к численному решению уравнения (1.12). Обозначив z=v₁, имеем:

(1.17)

или

где

(1.18)

Для решения уравнения (1.17) методом Ньютона нам потребуется выражение для производной функции f(z):

(1.19)

Зададимся исходными данными, чтобы подсчитать значения коэффициентов a , b, g уравнения (1.17):

k = 20000 Н/м,

m = 0.1 кг,

Q = 1000, T = 0.2p, F = 1000sin 10t,

Начальные условия и шаг интегрирования:

Тогда, согласно (1.12)

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

(1.20)

(1.21)

Зададим значения допустимых погрешностей для проверки условий (1.15):

(1.22)

Максимально допустимое количество итераций примем равным 5.

Начальное приближение к решению выберем

z⁰ = 0

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

Проверка завершения итераций:

Переход на следующую итерацию.

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

Проверка завершения итераций:

Переход на следующую итерацию.

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

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

Оба условия (1.15) удовлетворены, ограничение (1.16) не превышено. Решение достигнуто. Мы вычислили значение скорости для момента времени t₁, получив

Воспользовавшись формулами (1.10), определим значения ускорения и перемещения для этого же момента времени.

Решение для момента времени t₁ получено. Сделаем еще один шаг по времени для того, чтобы проиллюстрировать теперь выбор величины шага. Уравнения (1.9) - (1.12) справедливы для любого момента времени с учетом соответствующей замены подстрочных индексов. Для 2-го шага по времени имеем:

(1.23)

Так же, как и на первом шаге, приводим эту систему к одному уравнению относительно скорости:

где

(1.24)

Предварительно величину шага Dt₂ принимаем равной Dt₁, т.е. 0.001 с. Тогда, с учетом исходных данных и полученного на первом шаге решения, можно подсчитать коэффициенты a , b , g:

Опять имеем нелинейное уравнение:

, (1.25)

которое решаем методом Ньютона.

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

За начальное приближение к решению примем такое значение скорости, которое тело имело бы в момент времени t₂, если бы ускорение тела с момента времени t₂ не изменилось, т.е., считаем, что

(1.26)

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

Итак, начальное приближение (прогноз):

Опуская подробные выкладки (они аналогичны приведенным для первого шага по времени), итерации метода Ньютона приводят к следующей последовательности значений:

начальное приближение:

первая итерация:

вторая итерация: , решение достигнуто.

Мы получили, что при величине шага Dt₂=0.001, скорость для момента времени t₂

Пришло время оценить погрешность сделанного шага по времени.

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

, (1.27)

где v_i^p - явный прогноз величины скорости на i-м шаге, определяемый формулой

, (1.28)

v_i^c - значение скорости, которое мы получили в результате итерационного решения, используя неявную формулу

(1.29)

Рис. 1.6

Заметим, что соотношение (1.28) уже применялось нами при выборе начального приближения к решению в алгоритме метода Ньютона (смотри зависимость (1.26)).

Вычисление скоростей по формулам (1.28) и (1.29) и суть оценки локальной погрешности по формуле (1.27) поясняет рис. 1.6.

В момент времени t_i-₁ мы находимся в точке v_i-₁. Если для вычисления значения v_i мы воспользуемся явной формулой (1.28), то точка v_i=v_i^p будет лежать на касательной, проведенной к кривой v(t) в точке t_i-₁, поскольку a_i-₁ есть тангенс угла наклона этой касательной к оси абсцисс.

При вычислении v_i с использованием формулы (1.29) нам требуется значение a_i, т.е. тангенса угла наклона касательной, проведенной к кривой v(t) уже в точке t_i. Так как в момент времени t_i-₁ мы ничего не знаем о поведении функции v(t) при t=t_i и касательную к кривой v(t) в точке t_i тоже провести не можем, то мы вычисляем v_i =v_i^c не непосредственно по формуле (1.29), а путем совместного решения системы уравнений (1.7), куда входит и соотношение (1.29). При этом нам приходится последовательно приближаться к решению (т.е., к v_i^c) за несколько Ньютоновских итераций.

Рисунок 1.6. показывает, что явный прогноз v_i^p и скорректированное решение v_i^c лежат по разные стороны от кривой v(t), проходящей через точку t_i-₁. Чем больше разница между v_i^c и v_i^p, тем сильнее на текущем шаге отличается график скорости от прямой и тем выше погрешность интегрирования на шаге. Рисунок также позволяет понять, что уменьшение величины шага t_i приводит к уменьшению локальной погрешности, оцениваемой по формуле (1.27), поскольку уменьшается расхождение значений v_i^p и v_i^c.

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

Механизм определения величины шага, исходя из критерия локальной погрешности, достаточно прост. Задается значение предельно допустимой локальной погрешности на шаге интегрирования . По результатам очередного i-го шага сравниваются значения допустимой ()и фактически полученной локальной погрешности(). Если , то сделанный шаг признается успешным. Осуществляется переход на следующий шаг по времени; величина его для одношаговых методов интегрирования первого порядка точности, которым соответствуют используемые нами формулы (1.7), выбирается по зависимости:

, (1.30)

где Dt_i - величина совершенного шага по времени,

Dt_i+₁ - рекомендуемое значение следующего шага,

c - поправочный коэффициент, c < 1.

Если же , то значение сделанного шага Dt_i слишком велико и не обеспечивает требуемой точности. Поэтому необходимо провести расчет на i-м шаге еще раз, используя с уменьшенным значением Dt_i. В этом случае для выбора величины Dt_i также используется формула (1.30), только получаемое по ней значение шага используется не для следующего (i+1) -го шага, а для повторного расчета на текущем i-м шаге.

Возвращаясь к разбираемому примеру численного решения уравнения (1.6), проведем для 2-го шага интегрирования оценку локальной погрешности и величины шага.

В ходе расчета мы получили значения:

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

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

мы вынуждены констатировать, что , т.е. выполненный шаг со значением Dt₂=0.001 с не обеспечивает требуемой точности результатов и необходимо повторить расчет на 2-м шаге с уменьшенным значением Dt₂. Рекомендуемое значение Dt₂ для повторного расчета определим с помощью формулы (1.30), используя коэффициент c=0.8:

Результаты повторного расчета с шагом Dt₂ = 0.438e-3 с дают следующие значения прогноза скорости, скорректированного решения и локальной погрешности:

Поскольку полученное теперь значение , то второй шаг можно считать успешным с точки зрения заданной точности решения. Теперь дополним вычисленное значение скорости v₂=0.08509 м/с значениями ускорения a₂ и перемещения x₂, воспользовавшись формулами связи (1.7):

К настоящему моменту мы получили численное решение для двух точек временной оси:

Рис. 1.7

Следуя приведенному алгоритму, можно продолжить расчет и получить решение для всего интервала времени, который интересует исследователя.

Прежде чем подводить первые итоги, хотелось бы вернуться к рис. 1.6. для некоторых пояснений. В момент времени t_i-₁ мы находимся в точке v_i-₁. Через нее проходит кривая v(t). Она представляет собой так называемую интегральную кривую для момента времени t_i-₁, т.е., график скорости, соответствующий точному решению уравнения (1.6) при начальном условии

(1.31)

Поскольку уравнение (1.6) мы решаем приближенно, то фактически на каждом i-м шаге численного интегрирования по времени мы переходим с одной интегральной кривой, удовлетворяющей начальному условию (1.31), на другую интегральную кривую, которая уже представляет собой точное решение уравнения (1.6) при начальном условии

(1.32)

(на рис. 1.6. интегральная кривая для t=t_i обозначена пунктирной линией).

Поэтому результатом численного решения служит ломанная, проходящая через совокупность интегральных кривых, каждая из которых является точным решением уравнения (1.6) при начальных условиях, определяемых численным решением на текущем шаге (рис. 1.8).

Рис. 1.8

Подытожим основные моменты, существенные с точки зрения численного анализа рассмотренного примера.

Последовательность наших действий сводилась к следующему:

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

При формировании уравнения использовали 2-й закон Ньютона, являющийся одним из способов записи условия динамического равновесия.

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

3. Заменили дифференциальную связь между x, v и a алгебраическими уравнениями связи, справедливыми для выбранного метода интегрирования, сведя тем самым задачу получения решения в виде непрерывных функций к задаче отыскания совокупности значений неизвестной функции в отдельных точках временной оси:

где Dt_i - величина i-го шага по времени

(Dt_i = t_i - t_i-₁_);

x_i, v_i, a_i - значения x, v и a при t=t_i.

4. Полученную систему привели к одному уравнению, выразив x_i, a_i через v_i:

Таким образом, на каждом шаге по времени расчет сводился к решению нелинейного алгебраического уравнения вида f(z)=0, где z=v_i.

5. Решение нелинейного уравнения проводили методом Ньютона. Это итерационный численный метод (решение получаем приближенное, с заранее заданной точностью, за несколько проходов). Для получения решения на каждом проходе необходимо вычислять значения функции f(z) и ее производной df(z)/dz. Начальное приближение к решению определяем, используя формулу явного прогноза.

6. Точность численного интегрирования по времени оценивали путем контроля локальной погрешности на шаге интегрирования, зависящей от разности явного и неявного решения. При неудовлетворительной величине локальной погрешности повторяли расчет на текущем шаге с уменьшенным значением шага Dt_i.

7. Если локальная погрешность на шаге лежала в пределах допустимой, то считали шаг успешным и, используя вычисленное значение скорости v_i, подсчитывали ускорение a_i и перемещение x_i по уравнениям связи, справедливым для выбранного метода интегрирования.

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

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

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

2. Механизм формирования математической модели

Вернемся к рассмотрению системы на рис.1.1.

Перепишем еще раз уравнение равновесия:

(2.1)

Поскольку при численном интегрировании мы получаем решение в отдельных точках временной оси, то для каждого i-го момента времени уравнение (2.1) может быть записано в виде:

, (2.2)

где

, (2.3)

Считаем также, что для i-го момента времени значения x_i, v_i, a_i связаны уравнениями (1.7), которые ввиду их использования в наших дальнейших выкладках воспроизведем еще раз:

(2.4)

Обратите внимание, что уравнения (2.3) отличаются от аналогичных выражений (1.3) знаком. Это связано с тем, что в PRADIS при рассмотрении условий равновесия суммируются усилия, действующие со стороны системы на элементы, а не усилия со стороны элементов, как это было принято нами при выборе положительного направления для сил в соответствии с рис. 1.3.

Подстановкой (2.3) в (2.2) можно было бы получить уравнение вида (1.8), но мы этого делать не будем, поскольку нам нужно формировать и анализировать математическую модель по универсальному алгоритму. Итак, имеем достаточно универсальное уравнение равновесия вида (2.2), справедливое для каждого i-го момента времени. Заметим, что переход от записи уравнения равновесия в форме (2.1) к записи в форме (2.2) ознаменовал собой качественное изменение типа уравнения. Если соотношение (2.1) представляет собой дифференциальное уравнение (поскольку входящие в него зависимости для сил используют производные перемещения по времени), то соотношение (2.2) есть уже просто алгебраическое нелинейное уравнение (поскольку связь между x_i, v_i, a_iопределяется алгебраическими уравнениями (2.4)). А с нелинейным уравнением в форме (2.2) можно поступать также, как мы поступаем с уравнением (1.12), а именно: решать его методом Ньютона.

Имеем

, (2.5)

где

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

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

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

Теперь распишем производную df(z)/dz.

(2.6)

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

Так как z=v_i, то

(2.7)

Аналогично:

(2.7a)

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

, (2.8)

Продифференцируем выражения (2.8) по v_i

(2.9)

Теперь вычислим частные производные в выражениях (2.7), (2.7а), пользуясь зависимостями (2.3):

Для силы F_i^с:

(F_i^с не зависит от перемещения тела)

(F_i^с не зависит от скорости тела) (2.10)

(F_i^с не зависит от ускорения тела)

Для силы F_i^у:

(F_i^у не зависит от скорости тела) (2.11)

(F_i^у не зависит от ускорения тела)

Для силы F_i^в:

(F_i^в не зависит от перемещения тела)

(2.12)

(F_i^в не зависит от ускорения тела)

Для силы F_i^и:

(F_i^и не зависит от перемещения тела)

(F_i^и не зависит от скорости тела) (2.13)

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

, (2.14)

Суммируем слагаемые формулы (2.6):

(2.15)

Сравнивая результат с ранее полученной формулой (1.19), для которой коэффициенты берутся из соотношений (1.12), можно констатировать, что они похожи.

Теперь, имея возможность вычислять f(z) и df(z)/dz, продолжить расчет по описанному ранее алгоритму не представляет труда. Однако, как говорил один их героев Толкиена, "ситуация в настоящий момент может представляться требующей некоторых разъяснений ". Продравшись через частокол полных и частных производных, мы получили тот же результат, что и ранее, но терпения у читателя эти выкладки наверняка изрядно поубавили. Какими же приобретениями окупаются эти трудозатраты?

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

2. Развернутая форма нелинейного уравнения вида (1.12), к решению которого сводится расчет на каждом шаге по времени, тоже оказалась невостребованной.

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

Рис. 2.1

Если внимательно просмотреть наши рассуждения, то мы использовали следующую информацию (приводимый ниже список необходимой для формирования математической модели информации будем далее называть "перечнем"):

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

2. Условие равновесия сил, записанное для i-го момента времени, и представляющее собой нелинейное алгебраическое уравнение относительно x_i, v_i, a_i вида:

, (2.16)

где 4 - количество сил, сходящихся в узле стыковки (количество стыкующихся ветвей элементов).

Уравнение (2.16) еще называют топологическим, поскольку оно определяется топологией, т.е., структурой связей в схеме.

3. Выражения, позволяющие определить усилия в каждом элементе как функцию перемещения, скорости, ускорения и времени - (2.3). Это так называемые компонентные уравнения, описывающие поведение отдельной компоненты (элемента) схемы.

4. Выражения для частных производных от усилий, действующих в элементе, по перемещению, скорости, ускорению - см. зависимости (2.10)-(2.14).

5. Алгебраические уравнения связи x , v , a для текущего момента времени - формулы метода интегрирования (2.4).

6. Указание, относительно какой из переменных - x, v или a - вести решение нелинейного уравнения, т.е., что выбирается в качестве переменной z функции f(z) в формуле (2.5).

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

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

Итак, пусть имеется техническая система, процессы в которой требуют анализа. Расчетная схема соответствует рис.1.1.

Систему необходимо представить в виде совокупности элементов,стыкующихся по общим степеням свободы (узлам). Количество степеней свободы (узлов) у каждого элемента определяется разновидностью элемента (рис. 2.2).

Рис. 2.2

Существует понятие модели элемента. Пользователь собирает модель системы из моделей отдельных элементов, как игрушку в детском конструкторе. Его (пользователя) должна заботить только правильность сборки, остальные вопросы формирования математической модели - головная боль разработчиков программного обеспечения. Имея перед собой расчетную схему (рис. 1.1) и вычленив из нее элементы (рис. 2.2), пользователь находит модели соответствующих элементов в библиотеке моделей элементов программного комплекса и описывает структуру исследуемой схемы.

Описание структуры заключается в соединении элементов по общим степеням свободы и указании закрепленных узлов. Для полноты картины приведем кусок текста описания структуры на входном языке PRADIS:

$ FRAGMENT: Пример

# BASE: 1

# STRUCTURE:

Пружина' K (1 2; Коэффициент жесткости)

Нелинейный демпфер ' MUNL (1 2; Коэффициент вязкости)

Масса ' M (2; Масса тела)

Воздействие ' FSIN (2 1; Q, T, начальная фаза)

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

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

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

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

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

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

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

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

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

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

Если элемент имеет N степеней свободы, то длина вектора усилий элемента также N, а якобиан элемента имеет длину N*N*3.

Рис. 2.3

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

реализовал следующие зависимости для усилий и якобиана элемента (рис. 2.3):

(2.17)

, (2.18)

(2.19)

(2.20)

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

Поскольку в рассматриваемой расчетной схеме объекта узел 1 пружины закреплен, то в данном конкретном случае из всей информации, вычисляемой моделью и передаваемой "наверх", востребуется только та часть, которая связана с незакрепленным вторым узлом:

(2.21)

(2.22)

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

Чтобы конкретизировать сказанное, продолжим ранее начатое интегрирование примера, сделав очередной 3-й шаг по времени. При этом будем использовать формальный алгоритм, базирующийся на последовательности вычислений по формулам (2.5)-(2.15).

Напомним, что по результатам 2-го шага по времени (при t₂=1.438e-3 с) мы получили следующие значения неизвестных:

Величина шага равнялась ∆t₂=0.438 e-3 с, полученное на шаге значение локальной погрешности составило lp₂=0.000018.

Рекомендуемое значение ∆t₃ для следующего шага определяем по формуле (1.30) с учетом c=0.8 и δ_l = 0.001:

Дальше действуем по схеме, представленной на рисунках 2.4а-2.4в. Из рис. 2.4.а следует, что перед i-м шагом по времени должны быть известны значения t_i-₁, x_i-₁, v_i-₁, a_i-₁, ∆t_i. Можно убедиться, что перед началом 3-го шага мы действительно располагаем информацией о значениях t₂, x₂, v₂, a₂, ∆t₃.

Подробности алгоритма выполнения отдельного шага почерпнем из рис. 2.4.б.

1. Определяем значения коэффициентов приведения якобиана - , , (см. формулы (2.7), (2.7а)), зависящие от величины шага. Поскольку при расчетах на первых двух шагах за базисную переменную мы уже приняли скорость (т.е., z=v_i), то суммирование якобиана проводится по формулам (2.7), (2.7а), для которых коэффициенты приведения вычисляются по зависимостям (2.9):

Напомним, что соотношения (2.9) определяются формулами метода интегрирования (2.4).

2.Вычисляем начальное приближение к решению по формуле явного прогноза (1.28):

Обратите внимание, что в качестве начального приближения должно быть вычислено не только значение v_i⁰, но и значения x_i⁰, a_i⁰, необходимые для расчета в моделях элементов.

Поэтому:

3. Установив счетчик итераций равным 1, реализуем последовательность действий на 1-й итерации Ньютона (см. рис. 2.4в).

4. Обращение к моделям элементов. Вычисление вектора сил и якобиана каждого элемента по текущим значениям x₃⁰, v₃⁰, a₃⁰.

В настоящий момент ограничимся анализом информации только по незакрепленному узлу:

Рис. 2.4а

Рис. 2.4б

Рис. 2.4в

Пружина:

Демпфер:

Точечная масса:

Прикладываемая сила:

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

Полученная сумма является значением функции f(z) на текущей итерации (см. выражение (2.5)).

Вычислим теперь df(z)/dz, используя соотношения (2.6), (2.7), (2.7а):

6.Определяем приращение :

7. Вычисляем очередное приближение к решению

8. По полученному значению уточняем текущие значения x₃¹, v₃¹, a₃¹, используя формулы (2.81):

Вычисления на первой итерации Ньютона закончены.

9. Проверяем условия завершения ньютоновских итераций. Напомним, что ранее мы приняли следующие значения допустимых погрешностей для проверки условий (1.15):

Исходя из этих значений, заключаем, что

Таким образом, ньютоновские итерации на текущем шаге по времени необходимо продолжить.

10. Прежде, чем перейти на следующую итерацию, проверяем, не исчерпано ли предельное количество итераций:

j = 1 < j_max = 5

11. Увеличиваем счетчик номера итераций:

j = j + 1 = 1 + 1 = 2

12. Проверяем последовательность действий на рис. 2.4в для второй итерации метода Ньютона. Эти действия приведут нас к следующему решению:

Проверка условий окончания итераций покажет, что итерации еще не закончены:

Проверка:

j = 2 < j_max = 5

устраняет последние преграды с пути выполнения очередной, третьей итерации.

13. Третья итерация окажется последней. Будет получен следующий результат:

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

15. Вычисляем рекомендуемое по критерию локальной погрешности значение следующего шага.

Здесь следует вставить одну ремарку. Практика расчетов показала, что формула (1.30) приемлема только в определенном диапазоне соотношений δl / lp_i, а именно - вблизи единицы. При значительных отличиях δl / lp_i от единицы рекомендуемое формулой (1.30) значение шага, как правило, завышено и приводит к неоправданной потере шагов из-за неудовлетворения требованиям точности интегрирования. В этом мы убедимся уже сейчас, поскольку для выбора величины текущего шага использовали формулу (1.30) при соотношении δl / lpi = 0.001/0.000018 = 55.5. Как следствие этого, сделанный шаг с рекомендуемым формулой (1.30) значением шага ∆t₃=2.63e-3 привел нас к результату, когда сравнение фактически полученной (lp₃ = 0.034) и предельно допустимой (δl=0.001) локальных погрешностей определяет необходимость повторения расчетов на 3-м шаге с уменьшенным значением шага.

Скорректируем правило выбора шага по критерию локальной погрешности. Оно выглядит следующим образом:

(2.23)

Тогда, продолжая обсуждение алгоритма с прерванного места, рекомендуемое значение шага для повторного расчета на 3-м шаге определим с учетом (2.23):

Поскольку δl / lp₃ < 0.25, то

16. Устанавливаем ∆t₃= 0.061*e-3 с и повторяем вычисления на 3-м шаге, начиная с пункта 1.

Повторный расчет с этим значением шага приводит к следующим результатам для момента времени t₃=t₂+ t₃ = 1.499*e-3 с:

Локальная погрешность в пределах нормы. Рекомендуемое значение для следующего шага ∆t_рек=0.264*e-3с.

Расчет на третьем шаге по времени закончен.

Основное, на что хотелось бы обратить внимание по завершении разбора примера, - это разделение функций между собственно программой интегрирования и программами реализации моделей элементов. Программе интегрирования, работающей по алгоритму рис. 2.4а-2.4в, вообще говоря, все равно, какие процессы интегрировать. Ее зависимость от моделей элементов сводится только к своевременному получению от них векторов сил и матриц якобианов. А какие свойства отдельных элементов эта информация отражает, программу интегрирования это касаться не должно. Модели элементов, в свою очередь, имеют свой уровень независимости информации при вполне очерченных обязанностях перед "верхами". То есть, физические свойства отдельного элемента объекта отражаются в компонентных уравнениях на уровне модели элемента, а программа интегрирования работает на уровне уравнений равновесия потоков, не касаясь, по каким соотношениям составляющие этих потоков вычисляются. Такое разграничение функций определяет универсальность вычислительного ядра PRADIS, т.е. принципиальную возможность расчета любых объектов, процессы в которых подчиняются законам равновесия потоковых переменных (равновесие сил, электрических и тепловых потоков, расходов жидкости и газа).

3. Кратко об угловых степенях свободы, используемых в пространственных элементах PRADIS

Известно, что твердое тело из одного углового положения в другое можно перевести одним поворотом вокруг некоторой оси, называемой осью конечного вращения (теорема Эйлера). Обозначим e1, e2, e3 - направляющие косинусы оси конечного вращения, Fi - угол конечного вращения. Тогда можно ввести четыре кинематических параметра, описывающих угловое движение твердого тела [1,2]:

x1 = e1 * sin(Fi/2),

x2 = e2 * sin(Fi/2), (1)

x3 = e3 * sin(Fi/2),

x4 = cos(Fi/2),

и одно уравнение связи для этих параметров:

x1*x1 + x2*x2 + x3*x3 + x4*x4 = 1 (2)

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

Угловые степени свободы, принятые в пространственных элементах PRADIS, выражаются через кинематические параметры (1) следующим образом:

q1 = x1 * Lq,

q2 = x2 * Lq, (3)

q3 = x3 * Lq,

q4 = x4 * Lq,

где

Lq = sqrt(q1*q1+q2*q2+q3*q3+q4*q4). (4)

Первые три степени свободы являются внешними для моделей элементов, четвертая - внутренняя, скрытая от пользователя. Начальное значение потенциальной переменной, соответствующей внутренней степени свободы, устанавливается в моделях элементов равным 1.

Потоковыми переменными для первых трех степеней свободы являются моменты по глобальным осям X, Y, Z. Четвертая (внутренняя) потоковая переменная сдерживает изменение во времени величины Lq (4):

i4 = Mu * d(Lq)/dT, (5)

где Mu - коэффициент пропорциональности, одинаковый для всех степеней свободы такого рода и принимаемый в моделях элементов равным

Mu = DABSI/sqrt(MSHEPS). (6)

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

В частности:

· можно запрещать (базируя соответствующие узлы) движение по выбранным угловым степеням свободы, что равносильно сокращению размерности вектора, направленного по оси конечного вращения (например, при двух закрепленных угловых степенях свободы, точка может вращаться только вокруг оси, соответствующей незакрепленному узлу);

· связь по вращению между точками стыкуемых элементов тоже можно осуществлять (если это необходимо) не по всем трем степеням свободы, а только по выбранным.

С чем нужно быть осторожнее? В отличие от плоского вращения, первая и вторая производные от потенциальных переменных (3) не будут являться угловой скоростью и угловым ускорением соответственно. Поэтому, например, начальные условия, задаваемые моделью VN, не будут, в общем случае, определять начальную угловую скорость. Естественно, что и ПРВП типа V и A будут выводить не угловую скорость, а текущие значения первой и второй производной от потенциальной переменной. Значения же угловых скоростей и ускорений доступны из рабочего вектора некоторых моделей элементов, в частности J3O.

4. Литература:

1. Виттенбург Й. Динамика систем твердых тел. - М.: Мир, 1980.