Моделирование систем с одной степенью свободы

1. Задача. Имеется физическая система с одной степенью свободы, состоящая из инерционного элемента массой m, упругого элемента жесткостью k и диссипативного элемента с коэффициентом сопротивления r. Определить отклик системы x(t), а также ее первую и вторую производные v(t), a(t), на внешнее воздействие F_x = F_x (t), если известны начальные условия x(0), v(0).

2. Теория. Из второго закона Ньютона следует линейное неоднородное дифференциальное уравнение второго порядка:

Процессы, происходящие в последовательном колебательном контуре, состоящем из последовательно соединенных резистора R, конденсатора C и катушки индуктивности L, на который подано напряжение U(t), описываются уравнением:

где q -- заряд, проходящий по цепи, i = dq/dt -- сила тока. При этом реализуется электромеханическая аналогия "сила - напряжение".

Итак, неоднородное диффуравнение второго порядка описывает широкий класс задач, изучаемых в школьном и вузовском курсе физики: движение связанного с пружиной тела в вязкой среде под действием произвольной силы F_x (t), протекание тока через последовательно соединенные резистор, конденсатор и катушку индуктивности, подключенные к источнику ЭДС e(t) или источнику тока i(t).

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

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

Кроме того, физические явления, возникающие в системе, зависят от ее параметров m, r, k и начальных условий, к которым относятся координата x(0) и скорость v_x (0) в начальный момент времени.

3. Алгоритм. Дифференциальное уравнение второго порядка может быть решено методом конечных разностей Эйлера. Он состоит в том, что бесконечно малые приращения функции x(t) и ее первых двух производных v_x , a_x заменяются конечными разностями, а бесконечно малые приращения dt -- малыми, но конечными приращениями Δ t. Сначала, исходя из параметров системы m, r, k, ее координаты x и скорости v_x в момент времени t, расчитывается ее координата и скорость в следующий момент t + Δ t:

a_x (t + Δ t) = (F_x (t) - r v_x (t) - kx(t))/m,

v_x (t + Δ t) = v_x (t) + a_x (t + Δ t)Δ t.

x(t + Δ t) = x(t) + v_x (t + Δ t)Δ t.

Затем это состояние рассматривается как исходное и процедура расчета повторяется для момента времени t + 2Δ t и так далее. Одновременно с вычислениями строятся графики x = x(t), v_x = v_x (t), a_x = a_x (t).

Рассмотрим алгоритм численного решения уравнения.

1. Задают параметры физической системы m, k, r, зависимость внешнего воздействия от времени F_x (t), а также начальные условия x(0), v_x (0) и шаг по времени Δ t.

2. Начало цикла по t. Дают приращение по времени: переменной t присваивают значение t + Δ t.

3. Определяют ускорение, скорость и координату тела в момент t + Δ t:

a_x (t + Δ t) = (F_x (t) - r v_x (t) - kx(t))/m,

v_x (t + Δ t) = v_x (t) + a_x (t + Δ t)Δ t,

x(t + Δ t) = x(t) + v_x (t + Δ t)Δ t.

4. Результаты вычислений x(t + Δ t), v_x (t + Δ t), a_x (t + Δ t) выводят на экран в числовом виде либо строят соответствующие точки на координатной плоскости.

5. Возвращение к операции 2. Если цикл по t закончился, -- выход из цикла.

4. Компьютерная программа. При запуске программа рисует графики зависимостей координаты x = x(t), проекции скорости v_x = v_x (t) и проекции ускорения a_x = a_x (t).

Некоторые строчки программы заключены в скобки "(*" и "*)". Убрав скобки и активизировав соответствующие операторы, можно промоделировать различные явления.

program PROGRAMMA2;uses dos, crt, graph;Const Fm=10;w=5;m=2;r=0;k=0;Mx=20; Mv=40; Ma=8; Mf=2; Mt=100;dt=0.00006;Var x,v,a,F,t : Real;j,xx,vv,aa,FF,tt,Gd,Gm : Integer;BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); if GraphResult<>grOk then Halt(1); t:=0; v:=0; x:=-3; line(30,300,650,300); line(31,500,31,10); OutTextXY(50,20,'X, V, A');Repeat begin {Задание функции F=F(t)} t:=t+dt; (* F:=Fm*sin(w*t); *)(*If sin(w*t)<0 then F:=0; If sin(w*t)>0 then F:=Fm;*) F:=0; If t<1 then F:= Fm; If t>3 then F:=-Fm; a:=(F-r*v-k*x)/m; x:=x+v*dt; v:=v+a*dt; tt:=round(t*Mt); xx:=round(x*Mx); vv:=round(v*Mv); aa:=round(a*Ma); FF:=round(F*Mf); circle(30+tt,300-xx,1); circle(30+tt,300-vv,1); circle(30+tt,300-aa,2); end;until KeyPressed;CloseGraph;END.