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

[Kristina_Kretsu]

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

Рассмотрим колебательную систему с одной степенью свободы.

Тело массы m соединяется с неподвижным основанием посредством упругой пружины. Его движение происходит в среде с вязким сопротивлением под действием внешней силы.

Рис.1.1. Расчетная схема системы:
масса, пружина, демпфер, внешняя сила, направление перемещения x.

Для начального момента времени заданы:

• положение тела в пространстве;
• начальная скорость.

Требуется для каждого момента времени в интервале от t₀ до tконеч определить:

• перемещение тела;
• скорость тела;
• ускорение тела.

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

В примере будем считать, что:

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

Рис.1.2. Зависимости сил:
а) сила упругости пружины;
б) сила вязкого сопротивления;
в) синусоидальное внешнее воздействие.

---

Уравнение движения

Для получения дифференциального уравнения движения используем второй закон Ньютона.

Сумма сил, действующих на тело, определяет его движение:

(1.1)

или в форме равновесия:
(1.2)

Рис.1.3. Выбранные направления усилий

где:

(1.3)

Здесь:

x — перемещение тела;
v — скорость тела;
a — ускорение тела;
k — коэффициент жесткости пружины;
μ — коэффициент вязкого сопротивления;
Q — амплитуда внешней силы;
T — период внешнего воздействия.

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

(1.4)
Также:
(1.5)

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

(1.6)

Это уравнение описывает движение тела во времени.

---

Переход к численному решению

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

Временная ось представляется как последовательность точек:

t0, t1, t2, ..., ti, ti+1, ..., tn

Рис.1.4. Временная ось и шаги интегрирования

На каждом шаге времени нужно определить:

xi, vi, ai

где:

xi — перемещение в момент времени ti;
vi — скорость в момент времени ti;
ai — ускорение в момент времени ti.

Шаг интегрирования обозначается:

Δti = ti - ti-1

Интегрирование выполняется последовательно. Выбор очередного шага зависит от:

• требуемой точности;
• результатов расчёта на предыдущих шагах;
• локальной погрешности;
• устойчивости решения.

---

Формулы метода интегрирования

В рассматриваемом примере используются формулы неявного одношагового метода.

Связь между значениями на соседних шагах задаётся так:

(1.7)
где

t0 - начальное время,
x0 ,v0 - начальные значения перемещения и скорости,
tn - конечное время.

Для начального момента времени должны быть известны:

x0, v0.

Первый шаг интегрирования

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

x1, v1, a1

для момента времени:

t1 = t0 + Δt1

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

(1.8)

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

(1.9)

Чтобы решить её, выразим неизвестные x1 и a1 через v1 :

(1.10)

Получаем:
(1.11)

После подстановки получаем одно нелинейное уравнение относительно скорости v1:
(1.12)
где

(1.12)

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

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

Именно такое уравнение необходимо решать на каждом шаге времени.

---

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

Полученное уравнение решается методом Ньютона.

Общий вид нелинейного уравнения:

(1.13)
где:

z — неизвестная величина.

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

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

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

Рис.1.5.

(1.14)

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

(1.15)

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

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

(1.16)

Геометрически решение уравнения (1.13) означает поиск точки, где кривая f(z) пересекает ось z.
В методе Ньютона на каждой j-й итерации вместо самой кривой рассматривают касательную к ней в точке z(j−1) и находят, где эта касательная пересекает ось z.

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

(1.17)

или

где
(1.18)

Для метода Ньютона нужна производная:

(1.19)

---

Исходные данные примера

Зададим исходные данные:

k = 20000 Н/м
μ = 1000 Нс²/м²
m = 0.1 кг
Q = 1000
T = 0.2π
F = 1000 sin(10t)

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

x0 = 0
v0 = 0
Δt1 = 0.001 c

Тогда коэффициенты уравнения:

α = 1000
β = (20000 · 0.001) / 2 + 0.1 / 0.001 = 110
γ = -10

Таким образом,
(1.20)
(1.21)

Зададим допустимые погрешности:

(1.22)
Максимальное число итераций:

jmax = 5

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

z0 = 0

---

Итерации метода Ньютона

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

Проверка:

Условия завершения не выполнены, переходим к следующей итерации.

---

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

Проверка:

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

---

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

Проверка:

Одно из условий уже выполнено, но изменение решения ещё больше допустимого.

---

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

Проверка:

Оба условия (1.15) выполнены, ограничение (1.16) по числу итераций не превышено.

Вычислив значение скорости для момента времени t1, получается

v1 = 0.05913 м/с

---

Определение ускорения и перемещения

После нахождения скорости вычисляем ускорение и перемещение по формулам (1.10):

Таким образом, для момента времени t1 получено:

t1 = 0.001 c
x1 = 2.96e-5 м
v1 = 0.05913 м/с
a1 = 59.13 м/с²

---

Второй шаг интегрирования

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

Для второго шага:

Формулы связи:

(1.23)

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

где

(1.24)

Предварительно принимаем:

Δt2 = Δt1 = 0.001 c

Тогда:

Получаем нелинейное уравнение:

(1.25)
которое решается методом Ньютона.

---

Выбор начального приближения

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

За начальное приближение принимаем такое значение скорости, которое тело имело бы в момент времени t2, если бы ускорение с момента времени t1 не изменилось:

(1.26)

Это называется явным шагом или прогнозом, когда в формулу для скорости входит уже известное ускорение.

Это значение используется:

• как начальное приближение в методе Ньютона;
• для оценки локальной погрешности шага.

После итераций метода Ньютона получаем:

Решение достигнуто:

v2 = 0.11159 м/с

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

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

Для этого сравниваются:

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

Используется формула:

(1.27)

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

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

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

Рис.1.6. Явный прогноз и скорректированное решение
В момент времени мы находимся в точке . Если для вычисления значения мы воспользуемся явной формулой (1.28), то точка будет лежать на касательной, проведенной к кривой в точке , поскольку есть тангенс угла наклона этой касательной к оси абсцисс.

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

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

Расчёт локальной погрешности нужен не только для оценки точности решения. Основная его задача — определить:

• является ли выполненный шаг допустимым;
• каким должен быть следующий шаг интегрирования.

Для этого задаётся допустимое значение локальной погрешности:

После выполнения очередного шага сравниваются:

• допустимая локальная погрешность;
• фактически полученная локальная погрешность.

Если:

то шаг считается успешным.

В этом случае выполняется переход к следующему шагу времени, а рекомендуемое значение нового шага вычисляется по формуле:

(1.30)

где:

— выполненный шаг интегрирования;
— рекомендуемый следующий шаг;
c — поправочный коэффициент, c < 1.

Если же:

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

В этом случае:

• расчёт текущего шага повторяется;
• используется уменьшенное значение шага интегрирования.

Для выбора нового значения шага также используется формула:

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

---

Для второго шага:

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

Допустимая локальная погрешность:

Так как:

то выполненный шаг:

не обеспечивает требуемой точности. Расчёт второго шага необходимо повторить с уменьшенным шагом.

---

Расчёт нового шага

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

где:

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

В примере:

c = 0.8

Δt2 = 0.438e-3 c

Повторный расчёт выполняется уже с шагом:

---

Повторный расчёт второго шага

После повторного расчёта получаем:

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

Так как:

второй шаг считается успешным.

Теперь можно вычислить ускорение и перемещение по формуле (1.7):

---

Полученное решение для двух точек

К этому моменту получены значения:

t0 = 0
x0 = 0
v0 = 0

t1 = 0.001 c
x1 = 2.96e-5 м
v1 = 0.05913 м/с
a1 = 59.13 м/с²

t2 = 0.001438 c
x2 = 6.12e-5 м
v2 = 0.08509 м/с
a2 = 59.21 м/с²

Рис.1.7. Численное решение для двух точек временной оси

Смысл локальной погрешности

Явный прогноз и скорректированное решение показывают, насколько выбранный шаг соответствует реальному изменению решения.

Если разница между ними велика, значит:

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

Если разница мала, шаг можно считать приемлемым.

Рис.1.6. Сравнение явного прогноза и скорректированного решения

Уменьшение шага времени уменьшает расхождение между прогнозом и скорректированным решением, а значит уменьшает локальную погрешность.

---

Интегральные кривые

При численном решении дифференциального уравнения фактически на каждом шаге происходит переход с одной интегральной кривой на другую.

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

Рис.1.8. Численное решение как переход между интегральными кривыми

Основные итоги примера

Последовательность действий в рассмотренном примере следующая.

1. Сформировали дифференциальное уравнение

На основе второго закона Ньютона получили уравнение движения:

---

2. Представили уравнение в форме без явных производных

Записали отдельно:

---

3. Заменили дифференциальные связи алгебраическими

Для выбранного метода интегрирования использовали:

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

- значения x, v и a при .

Это позволило перейти от непрерывной функции к вычислению значений в отдельных точках времени.

---

4. Сводили систему к одному нелинейному уравнению

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

или:

где:

---

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

На каждой итерации вычислялись:

После этого уточнялось значение решения.

---

6. Оценивали локальную погрешность

Точность интегрирования оценивалась сравнением:

явного прогноза и скорректированного решения

Если локальная погрешность была больше допустимой, шаг повторялся с меньшим значением Δt.

---

7. При допустимой погрешности шаг считался успешным

После успешного шага вычислялись:

• скорость;
• ускорение;
• перемещение.

---

8. Следующий шаг выбирался автоматически

Величина очередного шага определялась на основе соотношения допустимой и фактической локальной погрешности.

---

Вывод

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

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

Такой подход позволяет PRADIS выполнять расчёт переходных процессов и автоматически управлять шагом интегрирования в зависимости от поведения модели.

---