ЧИСЕЛЬНИЙ АНАЛІЗ ПЕРЕХІДНОГО ПРОЦЕСУ В СТРУНІ

Карачун В.В., Мельник В.М.

Національний технічний університет України «КПІ»

ЧИСЕЛЬНИЙ АНАЛІЗ ПЕРЕХІДНОГО ПРОЦЕСУ В СТРУНІ

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

Розглянемо підготовку задачі до розв’язання її на ЕОМ.

Перейдемо від неперервних змінних до дискретних -

(1)

Замінимо частинні похідні відповідними співвідношеннями, наприклад, таким чином:

(2)

Якщо рівняння (2) роз’язати відносно , то отримаємо розрахункове співвідношення для так званої явної схеми інтегрування

(3)

де ; ; .

Перш ніж програмувати формулу (3), визначимося з формою представлення стану струни в довільний поточний момент часу . Функцію будемо представляти у вигляді масиву Yt такої структури

		-1	0	1	2		...	m
Yt	m					...

де m - число розрахункових відрізків, на які поділяється струна, (- довжина струни, - крок вздовж координати , як вже відзначалося раніше, вираз (3)).

Початковий стан струни (за ) визначимо відповідною підпрограмою-функцією Yn. Але, крім цього ще необхідно задати швидкості точок струни (підпрограма Dyn). Скористаємося щойно згаданими функціями для обчислення значень за умови . Для цього запишемо ряд Тейлора -

(4)

обмежившись в ньому трьома першими доданками. Тут означає , а відповідає .

Значення другої похідної у формулі (4) визначаємо за умови :

. (5)

Підставляємо (5) в (4) -

. (6)

Формула (6) покладена в основу підпрограми Init, яка формує масиви (для ) і Yt (для ), а також задає .

procedure Init;

var g,s:integer;

begin

for s:=0 to m do Yp[s]:=Yn(s*Hx);

for s:=1 to m-1 do

Yt[s]:=Yp[s]+Dyn(s*Hx)*Tau+E/2*(Yp[s-1]-2*Yp[s]+Yp[s+1])+

Tau2/2*Fxt(s,0);

g:=round(Xg(0)/Hx);

if Mg>0 then Yt[g]:=Yp[g]+Dyg*Tau+

sqr(Tau)/2*(Mh*fxt(g,0)+Fg(0)+

Ts/Hx*(Yp[g-1]-2*Yp[g]+Yp[g+1]))/(Mg+Mh);

case Ngl of

1:Yt[0]:=Yl(Tau);

2:if g>0 then Yt[0]:=Yp[0]+Dym0*Tau+

sqr(Tau)/2*(Ts*(Yp[1]-Yp[0])/Hx+

fxt(0,0)*Mh/2+Fl(0))/(M0+Mh/2)

else Yt[0]:=Yp[0]+(Dyn(0)*Mh/2+

Dym0*M0+Dyg*Mg)/(M0+Mh/2+Mg)+

(Ts*(Yp[1]-Yp[0])/Hx+fxt(0,0)*Mh/2+

Fl(0)+Fg(0))*sqr(Tau)/2/(M0+Mh/2+Mg);

3: if g>0 then Yt[0]:=Yp[0]+Tau*(Dyn(0)+Mh/2+

Dym0*M0)/(M0+Mh/2)+

sqr(Tau)/2/(M0+Mh/2)*(Ts*(Yp[1]-Yp[0])/Hx+

C0*(Ylk(0)-Yp[0]))

else Yt[0]:=Yp[0]+Tau*(Dyn(0)*Mh/2+

Dym0+M0+Dyg*Mg)/(M0+Mg/2+Mg)+

sqr(Tau)/2/(M0+Mh/2+Mg)*(Ts*Yp[1]-

Yp[0]+C0*(Ylk(0)-Yp[0])+fxt(0,0)*Mh/2+Fg(0))

end;

case Ngr of

1: Yt[m]:=Yr(Tau);

2: if g<m then Yt[m]:=Yp[m]+Tau*(Mh/2*Dyn(m*Hx)+

Mm*DyMm)/(Mh/2+Mm)+

sqr(Tau)/2/(Mh/2+Mm)*(Ts*(Yp[m-1]-

Yp[m])/Hx+Fr(0)+Fxt(m*Hx,0)*Mh/2)

else Yt[m]:=Yp[m]+Tau*(Mh/2*Dyn(M*Hx)+

Mm*DyMm+Mg*Dyg)/(Mh/2+Mm+Mg)+

sqr(Tau)/2/(Mm+Mh/2+Mg)*(Ts*(Yp[m-1]-

Yp[m])/Hx+Fr(0)+Fg(0)+Fxt(m*Hx,0)*Mh/2);

3:if g<m then Yt[m]:=Yp[m]+Tau*(Mm*DYMm+

Mh/2*Dyn(M*Hx))/(Mm+Mh/2)+

sqr(Tau)/2/(Mm+Mh/2)*(Ts*(Yp[m-1]-

Yp[m])/Hx+C1*(Yrk(0)-Yp[m])+

fxt(m*Hx,0)*Mh/2+Fr(0))

else Yt[m]:=Yp[m]+Tau*(Mm*DyMm+Mh/2*Dyn(M*Hx)+

Mg*Dyg)/(Mm+Mg+Mh/2)+

sqr(Tau)/2/(Mm+Mg+Mh/2)*(Ts*(Yp[m-1]-

Yp[m])/Hx+C1*(Yrk(0)-Yp[m])+

fxt(m*Hx,0)*Mh/2+Fr(0)+Fg(0))

end;

T:=Tau

end;

Тут Type CoefR=array[-1..200]of real;

Коефіцієнти Tau, Tau2, Е, що використовуються в процедурі Init, як і ряд інших коефіцієнтів, що будуть залучені далі, обчислюються процедурою CoefB (буде наведена нижче).

Звернемося тепер до формули (3), що реалізує явну схему. Вона «обслуговує» внутрішні вузли сітки, яка покриває область , в яких відсутня маса . Якщо ж маса в момент часу, що розглядається, знаходиться у вузлі з номером тоді для цього вузла повинна виконуватися умова

(7)

де - маса відрізку струни довжини .

Перейдемо в умові (7) до дискретних змінних

звідки

або після зведення подібних -

Після введення коефіцієнтів

(8)

де , , .

Зовнішні вузли сітки з номерами і «обслуговуються» граничними умовами. Розглянемо граничні умови послідовно:

. (9)

звідки

або

Ввівши коефіцієнти , приводимо до вигляду

, (10)

де ; ; .

Співвідношення (10) слушне за умови, коли ковзаюча вздовж струни маса в даний момент часу не виявляється на лівому кінці струни. В іншому випадку ця маса приєднується до , а сила, діюча на неї, додається до і тоді: