Здравствуйте!

В связи с вопросами читателей моей публикации [1] касательно условий возбуждения автоколебаний в механической системе, я решил описать явление возникновения и поддержания автоколебаний подробно, выделив основные области возникновения и применения автоколебаний.

В википедии автоколебания объясняют так [2]:

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

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

Автоколебания в технике

Автоколебательная система с запаздыванием (на примере электромеханического звонка)
Приведём пример электромеханического звонка:



При замыкании цепи кнопкой (К) электромагнит (Е) притягивает ударник, ударник бьёт по звонку и размыкает цепь питания электромагнита, механически связанным с ним контактом (Т) ударник (А) возвращается назад и процесс повторяется.

При рассмотрении процесса возникновения автоколебаний будем считать, что сила, действующая на боёк (А) звонка, изменяется пропорционально изменению тока в RL цепи.

Такое допущение сделано для упрощения рассмотрения, поскольку зависимость силы от тока в обмотке и зазора между бойком и полюсами значительно сложнее [3].

Ниже приведены конструкции электромеханических звонков и их упрощённая электрическая схема:



Боёк колеблется относительно установленного зазора согласно соотношению A*sin (w*t).
Решив численным методом дифференциальное уравнение RL цепи с начальными условиями
для замыкания и размыкания контакта, наложив на эти решения колебания бойка, получим:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt 
R=100;L=0.07;E=100;tm=L/R;T=0.0280;w=2*np.pi/T # параметры RL цепи и частота бойка w    
def dydt(y, t):# функция диф. уравнений для  численного решения
         return -y/tm
def dydt1(y1, t):
         return -y1/tm+E/L
plt.title('Колебательные процессы в звонке', size=12)
y = odeint(dydt, E/R, np.linspace(0,T/4,300))# ток (сила) при размыкании RL, y(0)=E/R
plt.plot(np.linspace(0,T/4,300), y,'b',linewidth=2,label='Сила Fp=Km*Ip')# график 
y1 = odeint(dydt1, 0,  np.linspace(T/2,3*T/4,300))# ток (сила) при замыкании RL, y(0)=0 
plt.plot( np.linspace(T/2,3*T/4,300), y1,'g',linewidth=2,label='Сила Fz=Km*Iz')
t2=np.linspace(0,T,300)
y2=(E/R)*np.sin(w*t2)# уравнение колебаний бойка
plt.plot(t2, y2,"--r",linewidth=2,label='Колебания бойка ')
t3=np.linspace(0,T/4,300)
t4=np.linspace(T/2,3*T/4,300)
y3=[0 for i in t3]
y4=[0 for i in t4]
plt.plot(t3, y3,"--k",linewidth=3,label='Запаздывание Fp от бойка ')
plt.plot(t4, y4,"--k",linewidth=3,label='Запаздывание Fz  от бойка')
plt.legend(loc='best')
plt.grid(True)
plt.show()

Для приближенной теории будем считать, что сила F?, выраженная последовательностью прямоугольных импульсов, которые возникает и исчезает мгновенно, но не в момент срабатывания контакта, а с запаздыванием ?=L/R. Добавим F? на график, получим:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt 
R=100;L=0.07;E=100;tm=L/R;T=0.0280;w=2*np.pi/T # параметры RL цепи и частота бойка w    
def dydt(y, t):# функция диф. уравнений для  численного решения
         return -y/tm
def dydt1(y1, t):
         return -y1/tm+E/L
plt.title('Колебательные процессы в звонке', size=12)
y = odeint(dydt, E/R, np.linspace(0,T/4,300))# ток (сила) при размыкании RL, y(0)=E/R
plt.plot(np.linspace(0,T/4,300), y,'b',linewidth=2,label='Сила Fp=Km*Ip')# график 
y1 = odeint(dydt1, 0,  np.linspace(T/2,3*T/4,300))# ток (сила) при замыкании RL, y(0)=0 
plt.plot( np.linspace(T/2,3*T/4,300), y1,'g',linewidth=2,label='Сила Fz=Km*Iz')
t2=np.linspace(0,T,300)
y2=(E/R)*np.sin(w*t2)# уравнение колебаний бойка
plt.plot(t2, y2,"--r",linewidth=2,label='Колебания бойка ')
def con_n(f):
         z=0
         if  np.sin(0)<=np.sin(f)<np.sin(T/4):                 
                  z=E/R
         elif np.sin(T/4)<= np.sin(f)<np.sin(3*T/4):                  
                  z=0
         elif np.sin(3*T/4)<=np.sin(f)<np.sin(T+tm):
                  z=E/R
         return z         
y3=[con_n(q) for q in t2]
plt.plot(t2, y3,"k",linewidth=3,label='Импульсы силы ') 
plt.legend(loc='best')
plt.grid(True)
plt.show()

Обозначим амплитуду силы F? через A?, получим разложения этой силы в ряд Фурье [4] (учитывая что x=a•sin(?•t), для первых двух членов ряда:



Будем считать, что постоянная составляющая силы A?/2 компенсируется регулировкой.
Тогда уравнение для колебаний бойка с учётом его приведенной массы m, трения r и изгибной жёсткости k примет вид:

(1)

Разделим обе части на массу бойка, введем обозначения, получим:

(2)

Для того, чтобы получить аналитические соотношения для частоты и амплитуды колебаний бойка, решим (2) приближённым методом [5]. Преобразуем (2) к виду:

(3)

Подставив в (3) при условии:
пропуская промежуточные выкладки получим соотношения для частоты и амплитуды автоколебаний:



На основании приведенных соотношений можно сделать вывод, что, при отсутствии самоиндукции, звонок работать не может, поскольку при L=0 нет запаздывания ?=0. Таким образом, при нулевом запаздывании автоколебания не возможны.

Автоколебания в измерительной технике (на примере механического резонатора вибрационных плотномеров)

Механические резонаторы в виде трубок пластин или цилиндров широко используются в вибрационных плотномерах, внешний вид которых приведен на рисунках:





Будем рассматривать резонатор c сосредоточенными эквивалентными параметрами: массой жесткостью и трением, характеризуемым коэффициентом

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

Запишем систему уравнений, описывающих движение резонатора в замкнутой системе возбуждения:



где: F- сила воздействия системы возбуждения на резонатор;
D(x)- неизвестный оператор обратной связи, подлежащий определению; Fупр — упругая восстанавливающая сила резонатора, которая в общем случае может описываться нелинейной функцией; х — поперечное смещение эквивалентной массы.

Воспользуемся выражением кубической упругой характеристики резонатора:



где ? — коэффициент, характеризующий отклонение реальной упругой характеристики от линейной.

Преобразуем записанную систему равенства к виду:



где — нелинейная составляющая упругой силы.

Структурная схема автоколебательной системы, работа которой характеризуется уравнениями, (1) приведена на рисунке:



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



Для решения задачи синтеза оптимальной системы возбуждения, воспользуемся методом гармонической линеаризации [6].

Механические резонаторы являются высокодобротными колебательными системами, которые можно рассматривать как узкополосные фильтры с выходным сигналом вида: x~A•cos?(?•?), где A— амплитуда колебаний резонатора; ? — частота колебаний, близкая к резонансной [7].

Поэтому для нелинейного элемента справедливо соотношение:



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



Рассмотрим уравнение для первой гармоники колебаний линеаризованной системы:



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



Из соотношения (2) определим смещение фазы, осуществляемое системой возбуждения:



Нетрудно установить, что частота автоколебаний не будет зависеть от трения при сдвиге фазы ?=?/2, тогда:



При этом условии из (2) следует, что система возбуждения должна быть дифференцирующим звеном D(i?)=(i*rэ* ?) т.е.



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

Система возбуждения состоит из трех элементов, D(i?)=Dп* Dу* D(в ), характеризующих частотные характеристики: приемника Dп, усилителя Dу и возбудителя D(в ) колебаний. Приемник является дифференцирующим – Dп=Kп* i*?, а возбудитель усилительным
звеном – Dв=Kв.

Для выполнения условия (5) усилитель должен иметь частотную характеристику:



Коэффициент усиления должен меняться вместе с изменением трения
Звено с переменным коэффициентом усиления можно реализовать простейшей нелинейностью типа двухпозиционного реле, имеющей частотную характеристику по первой гармонике [6]:



где — амплитуда первой гармоники на входе усилителя; — выходное напряжение усилителя, подаваемое на возбудитель колебаний.

Из (6) и (7) можно получить выражение для амплитуды установившихся автоколебаний резонатора:



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

Из изложенного можно сделать вывод, что частота автоколебаний резонатора вибрационного измерительного преобразователя не будет зависеть от трения при сдвиге фазы ?=?/2, когда система возбуждения является дифференцирующим звеном, и не будет зависит от амплитуды автоколебаний при стабилизации входного сигнала этого звена.

Автоколебания в радиотехнических генераторах (на примере решения уравнения
Ван-дер-Поля)

Обобщённая схема радиотехнического генератора автоколебаний приведена на рисунке:



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

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

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



Получим фазовые портреты (слева) и временные реализации колебаний (справа) осциллятора Ван-дер-Поля: ? =0.1, ? =1.1

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def f(y, x):
                y1, y2 = y              
                return [y2,(0.1-y2**2)*y2-y1]
x= np.linspace(0,200,601)
y0 = [0.0001,0.0001]
[y1,y2]=odeint(f, y0, x).T
plt.subplot(221)
plt .plot(y1,y2,linewidth=2)
plt .grid(True)
plt.subplot(222)
plt .plot(x,y1,linewidth=2, label='Автоколебания  \n квазигармонические')
plt.legend(loc='best')
plt .grid(True)
def f_1(y, x):
                y1, y2 = y              
                return [y2,(1.1-y2**2)*y2-y1]
x= np.linspace(0,200,601)
y0 = [0.0001,0.0001]
[y1,y2]=odeint(f_1, y0, x).T
plt.subplot(223)
plt .plot(y1,y2,linewidth=2)
plt .grid(True)
plt.subplot(224)
plt .plot(x,y1,linewidth=2, label='Автоколебания \n негармонические')
plt.legend(loc='best')
plt .grid(True)
plt.show()

Для ? =10.0

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
def f(y, x):
                y1, y2 = y              
                return [y2,(10.0-y2**2)*y2-y1]
x= np.linspace(0,200,601)
y0 = [0.0001,0.0001]
[y1,y2]=odeint(f, y0, x).T
plt.subplot(221)
plt .plot(y1,y2,linewidth=2)
plt .grid(True)
plt.subplot(222)
plt .plot(x,y1,linewidth=2, label='Автоколебания \n релаксационные')
plt.legend(loc='best')
plt .grid(True)
plt.show()

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

Химические колебания. Брюсселятор

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

К настоящему времени известно достаточно много колебательных реакций. Наиболее знаменитая из них была открыта Б.П. Белоусовым в 1950 г. и позднее детально изучена А.М. Жаботинским. Реакция Белоусова — Жаботинского (БЖ) представляет собой процесс окисления малоновой кислоты при взаимодействии в присутствии ионов в качестве катализатора.

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

Предложены различные математические модели реакции БЖ (например, модель Филда, Кереса и Нойеса — «орегонатор»), однако ни одна из них не описывает полностью все детали, наблюдаемые в эксперименте.

Мы рассмотрим более простой модельный пример: гипотетическую химическую реакцию, которая получила название Брюсселятор [8]. Уравнения этой реакции имеют вид:



Предполагается, что реагенты A и B имеются в избытке, так что их концентрации можно считать постоянными, а D и E ни в какие реакции не вступают. Составим кинетические уравнения, соответствующие реакции, которые описывают динамику концентраций реагирующих веществ.

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



Символами Y,X будем теперь обозначать соответствующие концентрации. Отметим, что из третьего уравнения системы следует, что скорость образования вещества X зависит от его концентрации, т.е. эта стадия реакции носит автокаталитический характер. Приведем уравнения (1) к безразмерному виду, содержащему минимальное число управляющих параметров. Для этого перейдём к новым переменным, Тогда уравнения (1) примут вид:



Построим фазовые портреты для: a=1.0; b=2.1; b=3.0;b=5.0

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
b=2.1
def f(y, x):
                y1, y2 = y              
                return [1-(b+1)*y1+(y1**2)*y2, b*y1-(y1**2)*y2]
x= np.linspace(0,100,1001)
y0 = [0.5,1.5]
[y1,y2]=odeint(f, y0, x).T
plt .plot(y1,y2,linewidth=2,label='a=1.0, b=%s'%b)
plt.legend(loc='best')
plt .grid(True)
plt.show()

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

Автоколебания в биосистемах (на примере модели Лотки Вольтерра –“Хищник -жертва”)

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

Рассмотрим модель взаимодействия хищников и их добычи, когда между особями одного вида нет соперничества. Пусть x и y— число жертв и хищников соответственно. Предположим, что относительный прирост жертв y'/x равен a-by, a>0, b>0, где a — скорость размножения жертв в отсутствие хищников, -by— потери от хищников.

Развитие популяции хищников зависит от количества пищи (жертв), при отсутствии пищи ( x=0 ) относительная скорость изменения популяции хищников равна y'/y =-c, c>0, наличие пищи компенсирует убывание, и при x>0 имеем y'/y =(-c +d*x), d>0.

Таким образом, система Вольтерра—Лотка имеет вид:



где a, b, c, d >0.

Рассмотрим фазовый портрет системы Вольтерра Лотка, для a=4 b=2.5, c=2, d=1 и графики ее решения с начальным условием x(0)=3, y(0)=1, построенные программой Python для численного решения системы обыкновенных дифференциальных уравнений:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
a=4;b=2.5; c=2; d=1 
def f(y, t):
                y1, y2 = y              
                return [y1*(a-b*y2),y2*( -c+d*y1)]
t = np.linspace(0,10,201)
y0 = [3, 1]
[y1,y2]=odeint(f, y0, t).T
plt.figure()
plt .plot(y1,y2,linewidth=2)
plt .grid(True)
plt.figure()
plt .plot(t,y1,linewidth=2)
plt .plot(t,y2,linewidth=2)
plt .grid(True)
plt.show()

Видно, что процесс имеет колебательный характер. При заданном начальном соотношении числа особей обоих видов 3: 1, обе популяции сначала растут. Когда число хищников достигает величины b=2.5, популяция жертв не успевает восстанавливаться и число жертв начинает убывать.

Уменьшение количества пищи через некоторое время начинает сказываться на популяции хищников и когда число жертв достигает величины x=c/d =2 (в этой точке y'=0), число хищников тоже начинает сокращаться вместе с сокращением числа жертв. Сокращение популяций происходит до тех пор, пока число хищников не достигнет величины y=a/b =1.6 (в этой точке x'=0).

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

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

Спасибо за внимание!!!

Ссылки:

1. Математическая модель вибрационного уровнемера с резонатором в виде консольной эллиптической трубки.
2. Автоколебания
3. Базовые уравнения задачи синтеза ш-образного электромагнита.
4. О классификации методов преобразования Фурье на примерах их программной реализации средствами Python.
5. Теодорчик К.Ф. Автоколебательные системы. ГИТЛ.,1952 г., 272 с.
6. Метод гармонической линеаризации средствами Python.
7. Жуков Ю.П. Вибрационные плотномеры. — М. Энергоатомиздат, 1991. —
144с: ил. — (Б-ка по автоматике; Вып 678)
8. И. Пригожин, Р. Лефевр Брюсселятор М. Наука,1968