Równania różniczkowe w modelowaniu procesów technicznych cz. 1 Wyróżnione

Nie ma wątpliwości, że analiza matematyczna jest podstawowym narzędziem matematycznym do opisywania zjawisk fizycznych zachodzących w technice. Spośród narzędzi analizy matematycznej prym wiodą równania różniczkowe, ponieważ pozwalają one na stworzenie dynamicznego modelu zjawiska, a tym samym na przeprowadzenie prognozy działania danego systemu.

Jednym z elementarnych, wręcz kanonicznych modeli jest oscylator harmoniczny. Stare powiedzenie głosi, że „fizyka dotyczy tej części przyrody, którą można opisać modelem oscylatora harmonicznego”. Przypomnijmy zatem jego równanie w kanonicznej formie:

\frac{d^{2}x(t)}{dt^{2}}+\omega^{2}x(t)=0

Wszędzie, gdzie w technice spotykamy drgania, w pierwszym przybliżeniu możemy spróbować modelować je za pomocą oscylatora harmonicznego. Przykładem takiego układu jest chociażby tradycyjny zegarek ręczny.

Nie zawsze jednak taki model wystarcza. W równaniu oscylatora harmonicznego \omega oznacza stałą częstość kołową. Dla ruchu postępowego \omega^2=k/m, gdzie k jest sztywnością sprężyny, a m masą. W modelu skrętnym, takim jak układ balansu i włosa zegarka, analogiczną rolę odgrywają sztywność skrętna i moment bezwładności. Samo powołanie się na drugą zasadę dynamiki nie uzasadnia pominięcia oporów ruchu, wymuszenia ani zmian parametrów.

Dla ustalenia uwagi rozważmy uproszczony model, w którym bezwładność jest stała, opory i wymuszenie pomijamy, natomiast sprężystość układu stopniowo maleje. Jedną z możliwych hipotez jest wykładniczy spadek kwadratu częstości kołowej:

\omega^2(t)=\omega_0^2 e^{-\alpha t},\qquad \omega_0>0,\quad\alpha>0,\quad t\geq0.

Parametr \omega_0 jest stały i równy początkowej częstości kołowej, a nie funkcją czasu. Parametr \alpha wyznacza tempo zmian; obydwa mają wymiar odwrotności czasu. Wykładnik \alpha t jest więc bezwymiarowy. Wybór funkcji wykładniczej jest założeniem modelu, które trzeba sprawdzić doświadczalnie — nie wynika z samego faktu zużywania się elementu. Ponieważ maleje tu sprężystość, a nie występuje siła proporcjonalna do prędkości, nie jest to model tłumienia lepkiego.

Otrzymujemy równanie:

\frac{d^2x}{dt^2}+\omega_0^2 e^{-\alpha t}x=0.

Poszukamy takiej zmiennej bezwymiarowej, która pozwoli rozpoznać znane równanie. Naturalną pierwszą próbą jest \tau=e^{-\alpha t}. Jeżeli x(t)=f(\tau(t)), reguła łańcuchowa daje:

\frac{d^2x}{dt^2}=\left(\frac{d\tau}{dt}\right)^2\frac{d^2f}{d\tau^2}+\frac{d^2\tau}{dt^2}\frac{df}{d\tau},

\alpha^2\tau^2\frac{d^2f}{d\tau^2}+\alpha^2\tau\frac{df}{d\tau}+\omega_0^2\tau f=0.

To poprawne przekształcenie, lecz ostatni człon zawiera \tau zamiast kwadratu nowej zmiennej. Aby otrzymać standardową postać równania Bessela, wprowadzamy pierwiastek z \tau oraz odpowiedni stały czynnik:

z=\frac{2\omega_0}{\alpha}\sqrt{\tau}=\frac{2\omega_0}{\alpha}e^{-\alpha t/2},\qquad x(t)=X(z(t)).

Teraz:

\frac{dz}{dt}=-\frac{\alpha}{2}z,\qquad\frac{d^2z}{dt^2}=\frac{\alpha^2}{4}z,

\frac{dx}{dt}=-\frac{\alpha}{2}z\frac{dX}{dz},\qquad\frac{d^2x}{dt^2}=\frac{\alpha^2}{4}\left(z^2\frac{d^2X}{dz^2}+z\frac{dX}{dz}\right).

Ponadto \omega_0^2e^{-\alpha t}=\alpha^2z^2/4. Po podstawieniu i podzieleniu przez \alpha^2/4 otrzymujemy równanie Bessela zerowego rzędu:

z^2\frac{d^2X}{dz^2}+z\frac{dX}{dz}+z^2X=0.

Dla dodatniego argumentu jego ogólne rozwiązanie jest kombinacją dwóch funkcji Bessela:

X(z)=A J_0(z)+B Y_0(z),\qquad x(t)=A J_0(z(t))+B Y_0(z(t)).

Funkcja Y_0(z) jest osobliwa w zerze, ale z(t)>0 dla każdego skończonego czasu. Punkt z=0 odpowiada dopiero granicy t\to\infty. Nie ma więc podstaw, aby automatycznie odrzucać składnik z Y_0. Obydwie stałe wyznaczamy z warunków początkowych.

Przyjmijmy, że oscylator został wychylony i puszczony bez prędkości początkowej:

x(0)=x_0,\qquad\dot x(0)=0,\qquad z_0=z(0)=\frac{2\omega_0}{\alpha}.

Korzystając z tożsamości J_0'(z)=-J_1(z) oraz Y_0'(z)=-Y_1(z), dostajemy:

\dot x(t)=\frac{\alpha z(t)}{2}\left[A J_1(z(t))+B Y_1(z(t))\right].

Warunki początkowe prowadzą zatem do układu:

\begin{cases}A J_0(z_0)+B Y_0(z_0)=x_0,\\ A J_1(z_0)+B Y_1(z_0)=0.\end{cases}

Tożsamość Wrońskiego J_1(z_0)Y_0(z_0)-J_0(z_0)Y_1(z_0)=2/(\pi z_0) daje jednoznacznie:

A=-\frac{\pi z_0x_0}{2}Y_1(z_0),\qquad B=\frac{\pi z_0x_0}{2}J_1(z_0).

Rozwiązanie spełniające oba warunki początkowe ma więc postać:

x(t)=\frac{\pi z_0x_0}{2}\left[J_1(z_0)Y_0(z(t))-Y_1(z_0)J_0(z(t))\right].

W chwili zerowej tożsamość Wrońskiego daje dokładnie x(0)=x_0, natomiast w pochodnej dwa iloczyny wzajemnie się znoszą, więc \dot x(0)=0.

Dla konkretnego przykładu wybierzmy \alpha=2\omega_0, czyli z_0=1, i mierzmy wychylenie w jednostkach niezerowego x_0. Wówczas:

\frac{x(t)}{x_0}\approx 1{,}22712623\,J_0\!\left(e^{-\alpha t/2}\right)+0{,}69122984\,Y_0\!\left(e^{-\alpha t/2}\right).

Współczynniki w tym zapisie są zaokrąglone; dokładny wzór powyżej spełnia warunki x(0)/x_0=1 oraz \dot x(0)/x_0=0. W ogólnym przypadku początkowy argument wynosi 2\omega_0/\alpha, a nie 1. Wartości funkcji Bessela i przebieg otrzymanej funkcji można badać numerycznie, na przykład za pomocą Computator.NET.

Malejąca sprężystość nie gwarantuje, że wychylenie pozostanie ograniczone. W tym modelu, gdy t\to\infty, siła przywracająca zanika, a prędkość dąży do -\alpha B/\pi. Gdy B\neq0, położenie ma więc asymptotycznie składnik liniowy w czasie. Nie jest to osobliwość w skończonej chwili, lecz konsekwencja przyjętego równania; jego przydatność dla rzeczywistego zegarka trzeba ocenić w zakresie, w którym założenia modelu są spełnione.

Równanie Bessela, jego rozwiązania i użyte tożsamości podaje NIST DLMF, §10.2, §10.5, wzór 10.5.2 oraz §10.6, wzór 10.6.3. Zachowanie przy małym argumencie opisuje §10.7.

Na tym prostym przykładzie prześledziliśmy proces poszukiwania rozwiązania równania różniczkowego poprzez sprowadzenie go do znanej postaci. Ze względu na prostotę wywodu rachunki zostały przeprowadzone w sposób szkicowy, jednak każdy czytelnik powinien bez problemu móc przeprowadzić je samodzielnie. Kolejne ciekawe analizy pojawią się w kolejnych artykułach.

Dziękuję za lekturę.

Kordian Czyżewski