Задачи, алгоритмы, программы (стр. 5 )

Рис. 7.2. Теплопроводность в неоднородной среде.

Задача 7.3. Имеется неоднородный стержень, известна начальная температура различных его точек, координаты и мощность источника тепла (холода). Один конец теплоизолирован, другой поддерживается при постоянной температуре. Необходимо рассчитать распределение температуры вдоль стержня.

Для решения задачи может быть использован алгоритм АЛ - 1. Алгоритм АЛ - 1 НАЧАЛО ПРОГРАММЫ ДЛЯ i:=1 ДО N ДЕЛАТЬ ЕСЛИ (i>18)И(i<22) ТО T[i]:=5 ИНАЧЕ T[i]:=0.1; ДЛЯ j:=1 ДО N ДЕЛАТЬ ЕСЛИ j>20 ТО k[j]:=1.8 ИНАЧЕ k[j]:=1; ПОВТОРЯТЬ ДО НАЖАТИЯ НА КЛАВИШУ {== kk:=kk+1; ДЛЯ i:=2 ДО N-1 ДЕЛАТЬ { ЕСЛИ (i>80) И (i<83) ТО q:=0.5 ИНАЧЕ q:=0; TT[i]:=T[i]+k[i]*(T[i+1]-2*T[i]+T[i-1])*dt/(h*h)+ (k[i+1]-k[i-1])*(T[i+1]-T[i-1])*dt/(4*h*h)+q*dt; } ДЛЯ i:=2 ДО N-1 ДЕЛАТЬ { T[i]:=TT[i]; } T[1]:=0.1; T[N]:=T[N-1]; ЕСЛИ kk/1000=round(kk/1000) ТО { ДЛЯ i:=2 ДО N ДЕЛАТЬ ПОСТАВИТЬ ТОЧКУ (i, T[i]); ==} КОНЕЦ ПРОГРАММЫ В нем используются два массива T[i] и TT[i], в которых сохраняются значения температуры элементов стержня в моменты t и t+1 соответственно. Расчет температуры осуществляется по выведенной выше формуле в цикле. Будем считать, что длина стержня l, его коэффициент температуропроводности k при x>0,2l равен 1,8, а при x<0,2l равен 1. Температура точек с координатами от 0,18l до 0,22l равна T1; все остальные точки имеют температуру T2<t1. Точки, координата $x$ которых лежит в интервале от 0,81l до 0,82l, нагреваются источником тепла известной мощностью q. Левый конец поддерживается при постоянной температуре, а правый --- теплоизолирован. Результат решения задачи представлен на рис. 7.3.</t

	uses crt, graph; { ПР - 7.3 } const n=100; m=1; h=1; dt=0.001; var ii, jj, kk, i,j, DV, MV, EC: integer; q, q1 :real; tt, t: array[1..N] of real; k: array[1..N] of real; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); EC:=GraphResult; If EC <> grOK then Halt(1); For i:=1 to N do If (i>18)and(i<22) then t[i]:=5 else t[i]:=0.1; For j:=1 to N do If j>20 then k[j]:=1.8 else k[j]:=1; Repeat For i:=2 to N-1 do begin If (i>80)and(i<83) then q:=0.5 else q:=0; tt[i]:=t[i]+k[i]*(t[i+1]-2*t[i]+t[i-1])*dt/(h*h)+ (k[i+1]-k[i-1])*(t[i+1]-t[i-1])*dt/(4*h*h)+q*dt; end; For i:=2 to N-1 do t[i]:=tt[i]; t[1]:=0.1; t[N]:=t[N-1]; If kk/3000=round(kk/3000) then For i:=2 to N do begin line(i*5+20,420-round(80*t[i]), (i-1)*5+20,420-round(80*t[i-1])); end; delay(1); kk:=kk+1; until KeyPressed; CloseGraph; END.
8. ВОЛНОВЫЕ И АВТОВОЛНОВЫЕ ПРОЦЕССЫ

Задача 8.1. Промоделируйте процесс распространения импульса в одномерной среде (струне). Изучите прохождение импульса через границу раздела двух сред с разными скоростями распространения возмущения.

Запишем одномерное волновое уравнение в конечных разностях: Для получения движущегося профиля волны необходимо организовать цикл по i, в котором перебираются все элементы среды, вычисляются их смещения из положения равновесия. После этого стирается предыдущая моментальная фотография волны и строится новая. Все это должно находиться внутри цикла по времени (программа ПР - 8.1). Результат моделирования прохождения волны через границу раздела двух сред приведен на рис. 8.1.

Uses crt, graph; { ПР - 8.1 } Const n=200; h=1; dt=0.02; Var i, j, DV, MV, EC : integer; eta, xi, xi1,xi2 : array[1..N] of real; t, b, vv : real; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); Repeat t:=t+dt; If t<6.28 then xi1[1]:=2*sin(t) else xi1[1]:=0; For i:=2 to N-1 do begin If i<N/2 then vv:=8 else vv:=3; xi2[i]:=2*xi1[i]-xi[i]+vv*(xi1[i-1]- 2*xi1[i]+xi1[i+1])/h/h*dt*dt; end; For i:=2 to N-1 do begin xi[i]:=xi1[i]; xi1[i]:=xi2[i]; end; For i:=1 to N do begin setcolor(black); circle(i*3-3,240-round(xi[i-1]*50),1); setcolor(white); circle(i*3-3,240-round(xi1[i-1]*50),1); end; until KeyPressed; CloseGraph; END.

Рис. 8.1. Одномерная волна: отражение и прохождение через границу раздела двух сред.

Задача 8.2. Создайте компьютерную модель распространения волны в двумерной среде (пластине). С помощью нее изучите прохождение волны от точечного источника через границу раздела двух сред.

Запишем волновое уравнение в конечных разностях: Для решения задачи используется программа ПР - 8.2.1. Для того, чтобы увеличить число узлов сетки, на которую разбита прямоугольная пластина используются указатели. В двумерном массиве xi[i, j] хранятся значения смещения различных элементов среды (xi(i, j)^) и их скорости (xi(i+N, j)^) в формате single (4 байта). Чтобы исключить влияние волн, отраженных от края пластины, искусственно создается эффект затухания. Коэффициент затухания волн вблизи краев пластины считается равным r=0,2. Директива компилятора {$N+} необходима, чтобы включить математический сопроцессор. Результат моделирования прохождения волны через границу раздела двух сред приведен на рис. 8.2.

{$N+} uses crt, graph; { ПР - 8.2.1 } const N=200; M=200; h=1; dt=0.2; Size=4; type RealPoint= ^single; var i, j,DV, MV : integer; vv, t,r : single; a: array[1..2*N] of pointer; Function xi(i, j: word): RealPoint; begin xi:=Ptr(Seg(a[i]^),ofs(a[i]^)+(j-i)*Size); end; Procedure Formula; begin vv:=4; {If i<50 then vv:=2 else vv:=4; } If (i>10)and(i<190)and(j>10)and(j<190) then r:=0 else r:=0.2; xi(i+N, j)^:=xi(i+N, j)^+vv*(xi(i, j+1)^-2*xi(i, j)^+xi(i, j-1)^+ xi(i+1,j)^-2*xi(i, j)^+xi(i-1,j)^)/h*dt-r*xi(i+N, j)^*dt; end; Procedure Raschet; begin For i:=2 to N-1 do For j:=2 to M-1 do Formula; For i:=2 to N-1 do For j:=2 to M-1 do xi(i, j)^:=xi(i, j)^+xi(i+N, j)^*dt; For i:=N-1 downto 2 do For j:=M-1 downto 2 do Formula; For j:=2 to M-1 do For i:=2 to N-1 do xi(i, j)^:=xi(i, j)^+xi(i+N, j)^*dt; end; Procedure Draw; begin For i:=2 to N-1 do For j:=2 to M-1 do begin setcolor(round(xi(i, j)^/5)); rectangle(10+i*2,450-j*2,13+i*2,453-j*2); rectangle(11+i*2,451-j*2,12+i*2,452-j*2); end; end; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); setbkcolor(15); setcolor(8); rectangle(11,451,12+200*2,450-200*2); For i:=1 to 2*N+1 do begin GetMem(a[i], M*Size); For j:=1 to M do xi(i, j)^:=0; end; Repeat t:=t+dt; xi(60,60)^:=100*sin(t*1); xi(120,120)^:=100*sin(t*1); Raschet; Draw; until KeyPressed; Release(HeapOrg); CloseGraph; END.

Рис. 8.2. Двумерная волна: прохождение через границу раздела двух сред.

В представленной ниже программе ПР - 8.2.2 использован другой прием: часть данных (скорости элементов среды eta) записывается в файл a. dat. Это позволяет увеличить число элементов массива xi[i, j].

{$N+} { ПР - 8.2.2 } uses crt, graph; const N=120; M=120; h=1; dt=0.2; var k, i, j, DV, MV : integer; vv, eta, a,t : single; xi: array[1..N,1..M] of single; var F: file of single; Procedure Formula; begin eta:=eta+vv*(xi[i, j+1]-2*xi[i, j]+xi[i, j-1]+ xi[i+1,j]-2*xi[i, j]+xi[i-1,j])/h/h*dt; end; Procedure Raschet; var aa: real; k: integer; begin Assign(F,'a. dat'); Reset(F); k:=0; For i:=2 to N-1 do For j:=2 to M-1 do begin inc(k); Seek(F, K); Read(F, eta); Formula; Seek(F, K); Write(F, eta); end; k:=0; For i:=2 to N-1 do For j:=2 to M-1 do begin inc(k); Seek(F, K); Read(F, eta); xi[i, j]:=xi[i, j]+eta*dt; end; k:=0; For i:=2 downto N-1 do For j:=2 downto M-1 do begin inc(k); Seek(F, K); Read(F, eta); Formula; Seek(F, K); Write(F, eta); end; k:=0; For i:=2 downto N-1 do For j:=2 downto M-1 do begin inc(k); Seek(F, K); Read(F, eta); xi[i, j]:=xi[i, j]+eta*dt; end; Close(F); end; Procedure Draw; begin for i:=2 to N-1 do for j:=2 to M-1 do begin setcolor(round(xi[i, j]/5)); rectangle(10+i*3,450-j*3,13+i*3,453-j*3); rectangle(11+i*3,451-j*3,12+i*3,452-j*3); end; end; Procedure FileCreate; begin Assign(F,'a. dat'); Rewrite(F); For i:=1 to N do For j:=1 to M do Write(F, a); end; BEGIN vv:=2; DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); setbkcolor(15); setcolor(8); rectangle(11,451,12+120*3,452-120*3); FileCreate; Repeat t:=t+dt; xi[20,45]:=50*sin(t*2); xi[60,65]:=50*sin(t*2); Raschet; Draw; until KeyPressed; CloseGraph; END.

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

Пусть источники имеют координаты (x1, y1), (x2, y2). Обозначим расстояния, проходимые волнами от источников до произвольной точки с координатами (x, y) через l1, l2. Результирующее смещение этой точки в момент T можно рассчитать по формулам: Задача решается с помощью программы ПР - 8.3. В ней перебираются все точки экрана и вычисляются смещения; в зависимости от их величины ставится точка соответствующего цвета.

uses crt, graph; { ПР - 8.3 } const m=500; h=1; dt=0.001; var i, j,jj, k,DV, MV, EC : integer; x1,x2,s, pi, lambda, xx, x,y : real; a1,a2,d, l1,l2,II, In1,In2,razn, fi1,fi2 : real; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); EC:=GraphResult; if EC <> grOK then Halt(1); pi:=3.1415; lambda:=6; d:=15; x1:=32-d; x2:=32+d; For i:=1 to 640 do For j:=0 to 480 do begin x:=0.1*i; y:=-10+0.1*j; l1:=sqrt((x-x1)*(x-x1)+y*y); l2:=sqrt((x-x2)*(x-x2)+y*y); If l1<>0 then A1:=1/l1; if l2<>0 then A2:=2/l2; s:=A1*sin(2*pi*l1/lambda)+A2*sin(2*pi*l2/lambda); putpixel(i, j,round(s*50)+1); end; ReadKey; CloseGraph; END.

Рис. 8.3. Интерференционная картина.

Задача 8.4. Численными методами исследуйте двухкомпанентную модель автоволнового процесса. Промоделируйте распространение автоволны в одномерной среде, аннигиляцию автоволн.

Двухкомпанентная модель исходит из того, что автоволновой процесс описывается двумя величинами: концентрацией активатора и концентрацией ингибитора. Активатор вызывает протекание химической реакции, а ингибитор препятствует этому. В случае распространения огня по полю, на котором быстро растет трава, активатором является температура T: когда она достигает критической температуры возгорания tkr трава загорается. Концентрация ингибитора u тем больше, чем больше дыма и меньше травы в данной точке плоскости. Когда концентрация ингибитора достигает порогового значения ukr (вся трава выгорает), химическая реакция прекращается. Изменение концентраций активатора и ингибитора описывается уравнением теплопроводности (диффузии). Будем использовать следующую математическую модель: Когда трава горит (t[i]>tkr и u[i]>ukr) выделяется энергия (P=2), а количество травы уменьшается (Q=-0,07). Если трава сгорела, но температура элемента активной среды еще высока (t[i]>0), температура должна быстро уменьшаться за счет охлаждения (P=-0,1). Когда элемент среды охладился (t[i]<0.2), а трава не выросла (u[i]>u0), она растет, что описывается слагаемым r*u[i]*(u0-u[i]), где r=0,003. При горении температура травы не может превышать 1,2 tkr. Для численного решения рассмотренной выше системы двух дифференциальных уравнений необходимо заменить производные их конечно-разностными аппроксимациями и выразить t[i] и u[i] в дискретный момент времени t+1. Используется программа ПР - 8.4. На рис. 8.4.1 и 2 представлены результаты моделирования распространения автоволны в одномерной среде и аннигиляции двух автоволн.

{$N+} { ПР - 8.4 } uses crt, graph; const n=320; M=40; t0=10; tkr=5; u0=3; ukr=1; h=1; dt=0.01; var time, i,j, DV, MV, EC, k : integer; a, b,r, P,Q : real; t, u: array[1..N] of single; gor: array[1..N] of integer; Procedure Raschet; --- Вычисление температуры --- begin P:=0; Q:=0; If (t[i]>tkr)and(u[i]>ukr) then gor[i]:=1 else gor[i]:=0; If gor[i]=1 then begin P:=2; Q:=-0.07; end; {трава сгорает} If (gor[i]=0)and(t[i]>0) then begin P:=-0.1; end; {охлаждается} If (gor[i]=0)and(u[i]<u0)and(t[i]1.2*tkr) then t[i]:=1.2*tkr; t[i]:=t[i]+a*(t[i+1]-2*t[i]+t[i-1])*dt/(h*h)+P*dt; {u[i]:=u[i]+Q*dt+r*u[i]*(u0-u[i])/5;} u[i]:=u[i]+b*(u[i+1]-2*u[i]+u[i-1])*dt/(h*h)+ Q*dt+r*u[i]*(u0-u[i]); end; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); a:=6; b:=3; setbkcolor(15); For i:=1 to N do begin u[i]:=u0; If (i>1)and(i<20){or(i>300)and(i100000)or(KeyPressed); Repeat until KeyPressed; CloseGraph; END. </u0)and(t[i]

Рис. 8.4.1. Распространение одномерной автоволны.

Рис. 8.4.2. Аннигиляция одномерных автоволн.

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

Рассмотрим автоволновой процесс в двумерной активной среде на примере распространения фронта огня по полю, на котором быстро растет трава. Исходя из двухкомпанентной модели запишем дифференциальные уравнения для активатора (температуры) T и количества травы u: Разобъем прямоугольную область на элементы, каждый из которых характеризуется величинами T[i, j] и u[i, j]. С целью увеличения числа элементов подключим математический сопроцессор, переменную T[i, j] объявим как single, а u[i, j] как word. Используется программа ПР - 8.5, результаты моделирования представлены на рис. 8.5.1 и 2. Для того, чтобы получить сферическую волну, элементы внутри круга небольшого радуса переводят в возбужденное состояние, то есть "поджигают", задавая их начальную темепературу достаточно высокой. Из рис. 8.5.1 видно, как образовавшаяся автоволна огибает препятствие (группу элементов черного цвета с низкой теплопроводностью). Для получения спиральной однорукавной автоволны задают фронт волны, оборванный в центре экрана (рис. 8.5.2).

{$N+} { ПР - 8.5 } uses crt, graph; const n=110; M=95; t0=10; tkr=5; u0=65000; ukr=32000; h=1; dt=0.05; var time, i,j, k,DV, MV, EC: integer; a, P,Q, uu: single; t: array[1..N,1..M] of single; u: array[1..N,1..M] of word; gor: boolean; Procedure Raschet; --- Вычисление температуры --- begin P:=0; Q:=0; If (t[i, j]>tkr)and(u[i, j]>ukr) then gor:=true else gor:=false; If gor=true then begin P:=1.9; Q:=-1600; end; {трава сгорает} If (gor=false)and(t[i, j]>0) then begin P:=-0.25; end; {быстро охлаждается} If (gor=false)and(u[i, j]<u0)and(t[i, j]2*tkr then t[i, j]:=2*tkr; If uu>=u0 then u[i, j]:=u0; If (uu>1)and(uu40)and(j<70) then a:=0.1 else a:=0.8; end; BEGIN a:=0.8; DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); setbkcolor(15); For i:=1 to N do For j:=1 to M do begin u[i, j]:=u0; t[i, j]:=0; {If sqr(i-45)+sqr(j-40)<100 then t[i, j]:=2*t0; {одиночная волна} If (i>45)and(i<55)and(j<40) then t[i, j]:=2*t0; {однорукавная волна} If (i>40)and(i<46)and(j<40) then begin u[i, j]:=round(0.3*ukr); t[i, j]:=0.9*tkr; end; end; Repeat inc(time); For i:=1 to N do For j:=1 to M do begin u[1,j]:=u[2,j]; u[N, j]:=u[N-1,j]; u[i,1]:=u[i,2]; u[i, M]:=u[i, M-1]; t[1,j]:=t[2,j]; t[N, j]:=t[N-1,j]; t[i,1]:=t[i,2]; t[i, M]:=t[i, M-1]; end; For i:=2 to N-1 do For j:=2 to M-1 do Raschet; For i:=N-1 downto 2 do For j:=M-1 downto 2 do Raschet; If time mod 100=1 then begin cleardevice; For i:=2 to N do For j:=2 to M-1 do begin setcolor(9); If u[i, j]>ukr then setcolor(2); If (t[i, j]>tkr)and(u[i, j]>ukr) then setcolor(12); If u[i, j]90)and(i<94)and(j>40)and(j<70) then setcolor(8); rectangle(i*5,j*5,i*5+4,j*5+4); rectangle(i*5+1,j*5+1,i*5+3,j*5+3); end; end; until KeyPressed; CloseGraph; END. </u0)and(t[i, j]

Рис. 8.5.1. Огибание сферической автоволнрй препятствия.

Рис. 8.5.2. Однорукавная спиральная автоволна.

Задача 8.6. Имеется модель одномерной упругой среды, состоящая из осцилляторов, связанных пружинами жесткостью q. Каждый осциллятор представляет собой материальную точку массой m, подвешенную на пружине жесткостью k. Промоделируйте распространение волны в данной среде, убедитесь в том, что фазовая скорость в общем случае не равна групповой, то есть имеет место дисперсия.

Рис. 8.6.1. Модель одномерной упругой среды.

На каждый осциллятор рассматриваемой модели одномерной упругой среды действуют: 1) сила упругости со стороны пружины осциллятора; 2) силы упругости со стороны соседних осцилляторов; 3) силы вязкого стрения. Из второго закона Ньютона следует, что ускорение θ, скорость η и смещение ξ i-ого осциллятора в момент времени t+1 могут быть найдены из формул: Ниже представлена программа ПР - 8.8.1, моделирующая распространение импульса по цепочке осцилляторов, связанных упругими связями. Она содержит цикл по времени t, и вложенный в него цикл по i, в котором последовательно перебираются все осцилляторы и вычисляются их ускорения, скорости и смещения. Тут же определяется кинетическая и потенциальная энергия колеблющихся осцилляторов. На экране строится "моментальная фотография волны" (зависимость смещения ξ от координаты x в фиксированный момент времени t) и зависимость энергии осцилляторов от координаты x. Эта компьютераная модель позволяет "пронаблюдать" распространение импульса, его отражение от края упругой среды, от "границы раздела двух сред". В случае когда k=0 импульс распространяется, не изменяя своей формы, его фазовая скорость (быстрота переноса колебаний) и групповая скорость (быстрота переноса энергии) равны (рис. 8.6.2). Если жесткость k и массу m осциллятора подобрать так, чтобы частота волны w была бы близка к его собственной частоте колебаний, то происходит дисперсия. Импульс растягивается, фазовая скорость заметно отличается от групповой (рис. 8.6.3). Ниже приведен другой вариант программы (ПР - 8.6.2) и результат ее использования (рис. 8.6.4).

{$N+} { ПР - 8.6.1 } uses dos, crt, graph; Const m=0.5; r=0.1; dt=0.0015; q=150; N=150; w=6; k=13; Var F, teta, Sum, S,t : single; i, x,Gd, Gm : integer; eta, ksi, y, a,b : array [0..N] of single; BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); setbkcolor(15); setcolor(8); line(10,250,630,250); Repeat t:=t+dt; cleardevice; For i:=1 to N-1 do y[i]:=ksi[i]; For i:=1 to N-1 do begin b[i]:=a[i]; x:=10+4*i; F:=q*(ksi[i-1]-ksi[i])+q*(ksi[i+1]-ksi[i]); teta:=(F-r*eta[i]-k*ksi[i])/m; eta[i]:=eta[i]+teta*dt; ksi[i]:=ksi[i]+eta[i]*dt; If w*t<2*3.142 then ksi[1]:=9*sin(w*t) else begin ksi[1]:=0; eta[1]:=0; end; a[i]:=(k*sqr(ksi[i])+q*sqr(ksi[i+1]-ksi[i])+m*sqr(eta[i]))/2/10; circle(x,370-round(ksi[i]*10),2); circle(x,370-round(ksi[i]*10),3); line(10+4*(i-1),369-round(ksi[i-1]*10),x,369-round(ksi[i]*10)); rectangle(x,250-round(a[i]),x+1,250); rectangle(x+2,250-round(a[i]),x+3,250); end; Sum:=0; S:=0; For i:=1 to N do begin Sum:=Sum+a[i]; S:=S+a[i]*i; end; circle(10+4*round(S/Sum),260,2); circle(10+4*round(S/Sum),260,3); delay(1000); until KeyPressed; CloseGraph; END.

Рис. 8.6.2. Дисперсия отсутствует. Групповая и фазовая скорости равны.

Рис. 8.6.3. Групповая скорость меньше фазовой. Имеет место дисперсия.

uses dos, crt, graph; { ПР - 8.6.2 } const m=0.05; r=0.015; k=12; dt=0.0005; q=100; N=120; pi=3.1415; var F, teta: Real; i, j,t, Gd, Gm: integer; eta, ksi, y,e, e0: array [1..N] of real; BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); Repeat inc(t); For i:=2 to N-1 do begin y[i]:=ksi[i]; e0[i]:=e[i]; ksi[N]:=0; If (5*pi*t*dt)<2*pi then ksi[1]:=2*sin(5*pi*t*dt) else ksi[1]:=0; F:=q*(ksi[i-1]-ksi[i])+q*(ksi[i+1]-ksi[i]); teta:=(F-r*eta[i]-k*ksi[i])/m; ksi[i]:=ksi[i]+eta[i]*dt+teta*dt*dt/2; eta[i]:=eta[i]+teta*dt; e[i]:=q*sqr(ksi[i-1]-ksi[i])/2+q*sqr(ksi[i+1]- ksi[i])+k*ksi[i]*ksi[i]/2+m*eta[i]*eta[i]/2; end; For i:=2 to N-1 do begin setcolor(8); line(5*i,400,5*i,400-round(y[i]*30)); line(5*i,250,5*i,250-round(e0[i]*2)); setcolor(15); line(5*i,400,5*i,400-round(ksi[i]*30)); line(5*i,250,5*i,250-round(e[i]*2)); end; until KeyPressed; CloseGraph; END.

Рис. 8.6.4. Перемещение фазы волны в диспергирующей среде.

9. РАСЧЕТ ТЕЧЕНИЯ ЖИДКОСТИ

Задача 9.1. Изучите установившееся потенциальное течение идеальной жидкости по трубе прямоугольного сечения. Внутри трубы имеются выступы и различные препятствия. Постройте линии тока.

Функция тока, характеризующая потенциальное течение жидкости, удовлетворяет уравнению Лапласа. Запишем его в конечных разностях: Для решения задачи используется программа ПР - 9.1. Результаты расчета функции тока представлены на рис. 9.1.

Uses crt, graph; { ПР - 9.1 } Const n=140; m=50; Var i, ii, j,jj, k,DV, MV, EC : integer; psi: array[1..N, 1..M] of real; Procedure Raschet; begin psi[i, j]:=(psi[i+1,j]+psi[i-1,j]+psi[i, j+1]+psi[i, j-1])/4; end; Procedure Gran_usl; begin For i:=1 to N do begin psi[i,2]:=-200; psi[i, M-1]:=200; psi[i,1]:=-200; psi[i, M]:=200; end; For j:=1 to M do begin psi[N-1,j]:=-204+8*j; psi[N, j]:=-204+8*j; end; For j:=1 to 25 do begin psi[2,j]:=-204+16*j; psi[N-1,j]:=-204+8*j; psi[1,j]:=-204+16*j; psi[N, j]:=-204+8*j; end; For i:=1 to N do For j:=1 to M do If (j>25)and((i-60)*(i-60)+(j-25)*(j-25)<200) then psi[i, j]:=0; For i:=1 to N do For j:=1 to M do If (abs(j-25)<15)and(abs(i-110)<5) then psi[i, j]:=0; For i:=1 to N do For j:=1 to M do If (j>25)and(i<30) then psi[i, j]:=200; end; Procedure Draw; --- Вывод на экран --- begin setcolor(round((psi[i, j]+200)/30)); If ((j>25)and((i-60)*(i-60)+(j-25)*(j-25)<200)) or ((j>25)and(i<30))or((abs(j-25)<15) and(abs(i-110)<5)) then setcolor(15); rectangle(i*4+50,j*4,i*4+53,j*4+3); rectangle(i*4+51,j*4+1,i*4+52,j*4+2); end; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); Repeat inc(k); Gran_usl; For i:=2 to N-1 do For j:=2 to M-1 do Raschet; Gran_usl; For j:=M-1 downto 2 do For i:=N-1 downto 2 do Raschet; If k/10=round(k/10) then begin cleardevice; For i:=2 to N-1 do For j:=2 to M-1 do Draw; end; until KeyPressed; CloseGraph; END.

Задача 9.2. Исследуйте установившееся безвихревое течение вязкой жидкости в длинной трубе (канале) постоянного сечения, внутри которой движется бесконечно длинное тело.

Это течение инвариантно по отношению к переносам в направлении движения. Для его расчета следует определить скорости в сечении, перпендикулярном направлению течения (оси z), то есть решить уравнение: Необходимо правильно задать граничные условия. Слои вязкой жидкости, прилегающие к поверхности твердого тела, имеют одинаковую с ним скорость. Если жидкость имеет свободную поверхность, то скорость частиц этой поверхности равна скорости частиц, расположенных слоем ниже. Результаты расчета обтекания длинного корпуса корабля и катамарана (рис. 9.2).

Uses crt, graph; { ПР - 9.2 } Const n= 90; m=58; h=1; dt=0.02; Var ii, jj, kk, i,j, DV, MV, EC : integer; vv, v: array[1..N, 1..M] of real; Procedure Gr_usl; begin For i:=1 to N do For j:=1 to M do v[i, j]:=vv[i, j]; For i:=2 to N-1 do For j:=2 to M-1 do begin If (j<5)and(abs(5-i)<5) then v[i, j]:=0; If (j<40)and(abs(i-50)<1+20*sqr(cos(j/25))) then v[i, j]:=300; end; For i:=1 to N do v[i, M]:=0; For j:=1 to M do begin v[1,j]:=0; v[N, j]:=0; end; end; Procedure Raschet; begin vv[i, j]:=v[i, j]+(v[i, j+1]-2*v[i, j]+v[i, j-1])*dt/(h*h) +(v[i+1,j]-2*v[i, j]+v[i-1,j])*dt/(h*h); end; procedure Draw; begin setcolor(round(v[i, j]/30)+1); If v[i, j]<1 then setcolor(9); rectangle(i*3+9,j*3,i*3+2,j*3+2); rectangle(i*3+10,j*3,i*3+1,j*3+1); end; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); Repeat Gr_usl; For i:=1 to N-1 do For j:=2 to M-1 do Raschet; Gr_usl; For j:=M-1 downto 2 do For i:=N-1 downto 2 do Raschet; For i:=1 to N do vv[i,1]:=vv[i,2]; kk:=kk+1; If kk/30=round(kk/30) then For i:=2 to N-1 do For j:=2 to M-1 do Draw; until KeyPressed; CloseGraph; END.