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

Рис. 2.2. К доказательству брахистохронных свойств циклоиды.

Задача 2.3. Промоделируйте скольжение частицы по сферической поверхности и ее движение после отрыва. Найдите точку отрыва.

Допустим, шайба скользит по шероховатой сферической поверхности и отрывается от нее на некоторой высоте h. Необходимо рассчитать траекторию шайбы и время движения. Из второго закона Ньютона можно вычислить ускорение, скорость и координаты точки в последовательные моменты времени: В момент отрыва шайба движется с ускорением g, при этом ее нормальное ускорение равно: После отрыва от сферической поверхности шайба движется под действием силы тяжести: Используется программа ПР - 2.3, в ней определяется высота h отрыва шайбы и общее время движения, которое выводится на экран. Результат моделирования приведен на рис. 2.3.

uses dos, crt, graph; {ПР - 2.3} const r=200; m=1; mu=0.1; g=10; dt=0.001; pi=3.1415926; var a, at, v,vx, vy, x,y, t : real; Gd, Gm, n : integer; tt : string; BEGIN Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); Line(20,400,600,400); Circle(120,400,199); a:=0.3; x:=r*sin(a); y:=r*cos(a); Repeat inc(n); t:=t+dt; If v<sqrt(g*r*cos(a)) then begin a:=arctan(x/y); at:=g*(sin(a)-mu*cos(a)); v:=v+at*dt; x:=x+v*cos(a)*dt; y:=y-v*sin(a)*dt; vx:=v*cos(a); vy:=-v*sin(a); end else begin vx:=v*cos(a); vy:=vy-g*dt; x:=x+vx*dt; y:=y+vy*dt; Circle(250,400-round(y),1); end; Circle(120,400,2); If n mod 300=0 then Circle(120+round(x),400-round(y),3); until (KeyPressed)or(y<-1); Str(round(t*1000),tt); OutTextXY(10,10,tt); Readkey; END.

Рис. 2.3. Скольжение частице по сферической поверхности.

Задача 2.4. Изучите движение лыжника, скатывающегося с горы и прыгающего с трамплина на склон. Профили горки и склона заданы.

Пусть профиль горки и склона, на который прыгает лыжник, задаются уравнениями: Для расчета движения лыжника используются уравнения: Программа ПР - 2.4 содержит цикл по времени, в котором пересчитываются координаты и скорость лыжника в последовательные моменты времени, а результаты выводятся на экран (рис. 2).

uses dos, crt, graph; { ПР - 2.4 } const r=200; m=1; mu=0.2; g=10; dt=0.001; dx=0.1; pi=3.1415926; var a, b,x1,x2,y1,y2,at, v,vx, vy, x,y, t : real; Gd, Gm, n : integer; tt : string; BEGIN Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); b:=0.015; x:=-160; y:=b*x*x; Repeat inc(n); t:=t+dt; delay(1); If x<30 then begin x1:=x; y1:=b*x1*x1; x2:=x+dx; y2:=b*x2*x2; a:=arctan((y1-y2)/(x2-x1)); at:=g*(sin(a)-mu*cos(a)); v:=v+at*dt; x:=x+v*cos(a)*dt; y:=y-v*sin(a)*dt; vx:=v*cos(a); vy:=-v*sin(a); end else begin vy:=vy-g*dt; x:=x+vx*dt; y:=y+vy*dt; end; If x<30 then Circle(120+round(0.5*x),280-round(0.5*y),1); If n mod 300=0 then Circle(120+round(0.5*x),280-round(0.5*y),3); until (KeyPressed)or((y<0)and(y<-sqrt(50*x-3000))); x:=60; Repeat x:=x+1; y:=-sqrt(50*x-3000); Circle(120+round(0.5*x),280-round(0.5*y),1); until x>1000; Str(round(t*1000),tt); OutTextXY(40,330,tt); Readkey; END.

Рис. 2.4. Движение лыжника.

Задача 2.5. Из пушки вылетает снаряд с известной скоростью. Определить угол под которым необходимо произвести выстрел для того, чтобы снаряд попал в цель. Сила сопротивления, действующая на снаряд, пропорциональна его скорости.

При движении снаряда на него со стороны воздуха действует сила вязкого трения. Запишем второй закон Ньютона в конечных разностях: Для решения задачи используется метод стрельбы: в программе ПР - 2.5 определяется траектория движения снаряда при заданной начальной скорости и произвольном угле α между стволом орудия и горизонталью. Если снаряд не поразил цель, то угол α увеличивают на некоторый шаг Δα и повторяют расчет, затем еще раз и еще раз до тех пор, пока снаряд не попадет в цель. Искомое значение угла в радианах выводится на экран компьютера. Задача имеет два решения. На рис. 2.5 показан метод нахождения одного из них.

uses crt, graph; {Метод стрельбы} { ПР - 2.5 } var x, y,vx, vy, ax, ay, v,Fx, Fy, alfa : real; Gd, Gm, i: integer; ugol: string; Label metka; const m=100; dt=0.005; r=1.5; xc=490; yc=30; BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); Line(10,0,10,450); Line(0,450,640,450); Circle(xc+10,450-yc,2); v:=10; For i:=1 to 40 do begin x:=0; y:=0; alfa:=0.6+0.02*i; str(alfa, ugol); vx:=v*cos(alfa); vy:=v*sin(alfa); Repeat Fy:=-3; Fx:=0; ax:=(Fx-r*vx)/m; ay:=(Fy-r*vy)/m; vx:=vx+ax*dt; vy:=vy+ay*dt; x:=x+vx*dt; y:=y+vy*dt; Circle(round(x)+10,450-round(y),1); If sqr(x-xc)+sqr(y-yc)<16 then goto metka; until (y<0)or(KeyPressed); end; metka: OutTextXY(10,10,ugol); Repeat until KeyPressed; CloseGraph; END.

Рис. 2.5. Определение направления ствола методом стрельбы.

Задача 2.6. Ракета движется вокруг Земли по круговой орбите радиуса r. При однократном включении двигателя в течение заданного интервала времени скорость ракеты увеличивается в k раз. Промоделируйте переход спутника на круговую орбиту с большим радиусом R.

Свяжем с центром Земли инерциальную систему отсчета. На спутник действует сила гравитационного притяжения, проекции его ускорения, скорости, а также его координаты в последовательные моменты времени вычисляются по формулам: В программе ПР - 2.6 моделируется переход спутника с одной круговой орбиты на другую более высокую орбиту. Для этого в моменты времени t=20000 и t=40000 следует на небольшое время включить двигатели так, чтобы скорость спутника возрасла в k раз. Результат расчета траектории движения спутника представлен на рис. 2.6.

uses crt, graph; { ПР - 2.6 } var v, B, q, F, Fx, Fy : real; r, x, y, vx, vy, ax, ay, t : real; Gd, Gm, i: integer; const M=1500; mm=200; dt=0.005; k1=1.2; k2=1.22; BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); If GraphResult <> grOk then Halt(1); Line(320,240,640,240); Line(320,240,320,0); Circle(320,240,5); x:=0; y:=80; vx:=4.1; vy:=0; Repeat begin t:=t+1; r:=sqrt(x*x+y*y); F:=M*mm/(r*r); Fx:=-F*x/r; Fy:=-F*y/r; ax:=Fx/mm; ay:=Fy/mm; If t=20000 then begin vx:=k1*vx; vy:=k1*vy; Circle(round(x)+320,240-round(y),3); end; If t=40000 then begin vx:=k2*vx; vy:=k2*vy; Circle(round(x)+320,240-round(y),3); end; vx:=vx+ax*dt; vy:=vy+ay*dt; x:=x+vx*dt; y:=y+vy*dt; Circle(round(x)+320,240-round(y),1); end; until KeyPressed; CloseGraph; END.

Рис. 2.6. Переход спутника на более высокую орбиту.

3. МЕХАНИКА СИСТЕМЫ МАТЕРИАЛЬНЫХ ТОЧЕК

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

Рис. 3.1. Движение тележки с грузом при ударе (вязкое трение).

При взаимодействии тележки с пружиной на тележку действует сила упругости пружины, груз начинает скользить и со стороны него на тележку действует сила вязкого трения. Из второго закона Ньютона можно определить ускорения тележки и груза, а затем вычислить их скорости и координаты: Для моделирования этого явления используется программа ПР - 3.1, содержащая цикл по времени, в котором вычисляются скорости и координаты тележки и груза. На экране появляется изображение движущихся тележки и груза до и после взаимодействия с пружиной.

uses crt, graph; { ПР - 3.1 } var F, Ftr, x1,x2,v1,v2,a1,a2,xx1,xx2,t: real; Gd, Gm, i: integer; const m1=1; m2=1.2; dt=0.005; k=0.1; r=0.1; BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); line(220,0,220,480); x1:=250; v1:=-5; x2:=250; v2:=-5; Repeat If x1<0 then F:=-k*x1 else F:=0; Ftr:=r*abs(v2-v1); a1:=(F-Ftr)/m1; v1:=v1+a1*dt; x1:=x1+v1*dt; a2:=Ftr/m2; v2:=v2+a2*dt; x2:=x2+v2*dt; setcolor(8); Rectangle(round(xx1)+220,240,round(xx1)+300,250); Circle(round(xx2)+290,230,3); setcolor(15); Rectangle(round(x1)+220,240,round(x1)+300,250); Circle(round(x2)+290,230,3); delay(5); xx1:=x1; xx2:=x2; until KeyPressed; CloseGraph; END.

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

Рис. 3.2. Движение тележки с грузом при ударе (сухое трение).

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

uses crt, graph; { ПР - 3.2 } var F, Ftr, x1,x2,xx1,xx2,v1,v2,a1,a2,t: real; Gd, Gm, i: integer; const m1=1; m2=1; dt=0.005; k=0.1; g=10; mu=0.03; BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); line(220,0,220,480); x1:=250; v1:=-7; x2:=250; v2:=-7; Repeat If x1<0 then F:=-k*x1 else F:=0; If (v1=v2) then If (abs(a2)+abs(a1)<=mu*g) then Ftr:=m2*(abs(a2)+abs(a1)) else Ftr:=-m2*(abs(a2)+abs(a1)); If v1<>v2 then begin If (v1>0)and(v2<0) then Ftr:=mu*m2*g; If (v1<0)and(v2>0) then Ftr:=-mu*m2*g; If (v1>0)and(v2>0)and(abs(v2)>abs(v1)) then Ftr:=-mu*m2*g; If (v1>0)and(v2>0)and(abs(v2)<abs(v1)) then Ftr:=mu*m2*g; end; a1:=(F-Ftr)/m1; v1:=v1+a1*dt; x1:=x1+v1*dt; a2:=Ftr/m2; v2:=v2+a2*dt; x2:=x2+v2*dt; setcolor(8); Rectangle(round(xx1)+220,240,round(xx1)+400,248); Circle(round(xx2)+400,230,3); setcolor(15); Rectangle(round(x1)+220,240,round(x1)+400,248); Circle(round(x2)+400,230,3); delay(5); xx1:=x1; xx2:=x2; until KeyPressed; CloseGraph; END.

Задача 3.3. Тележка с телом движется с некоторой скоростью. В тело попадает пуля и застревает в нем (или проходит сквозь него). Как изменяются скорости тележки и пули в процессе взаимодействия?

Рис. 3.3.1. Движение тела после попадания пули.

Запишем формулы, позволяющие рассчитать ускорение пули и тела, их скорости и координаты в следующий момент времени t+1: Программа ПР - 3.3 содержит цикл по времени, в котором вычисляются сила взаимодействия пули и тела, их ускорения, скорости и координаты. Получающиеся графики зависимости скоростей пули и тела от времени представлены на рис. 3.3.2: слева -- пуля и тело двигались навстречу друг другу, пуля прошла на вылет; справа -- пуля догнала тело и застряла в нем.

uses crt, graph; { ПР - 3.3 } var Ftr, x1,x2,v1,v2,a1,a2,xx1,xx2,t: real; Gd, Gm, i: integer; const m1=4; m2=1; dt=0.005; r=0.1; BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); line(220,0,220,40); line(0,350,640,350); x1:=250; v1:=-4; x2:=0; v2:=10; Repeat If abs(x2-x1)<40 then Ftr:=r*abs(v2-v1) else Ftr:=0; a1:=Ftr/m1; v1:=v1+a1*dt; x1:=x1+v1*dt; t:=t+dt; a2:=-Ftr/m2; v2:=v2+a2*dt; x2:=x2+v2*dt; setcolor(8); Rectangle(round(xx1)+220-40,10,round(xx1)+260,30); Circle(round(xx2)+220,20,3); setcolor(15); Rectangle(round(x1)+220-40,10,round(x1)+260,30); Circle(round(x2)+220,20,3); delay(10); xx1:=x1; xx2:=x2; Circle(round(t*10),350-round(15*v1),1); Circle(round(t*10),350-round(15*v2),1); until KeyPressed; CloseGraph; END.

Рис. 3.3.2. Графики зависимостей скрости пули и тележки от времени.

Задача 3.4. Создайте двумерную модель газа, рассматривая молекулы как частицы, между которыми действуют силы отталкивания.

Рассмотрим систему материальных точек, движущихся по плоскости, на которые действуют внешние силы Fi (например, сила тяжести) и внутренние силы F'i. Из второго закона Ньютона следует: Программа, моделирующее движение молекул, должна содержать цикл по времени t, а также цикл, в котором перебираются все молекулы, вычисляются проекции действующих на сил, определяются ускорение, скорость и координата в следующий момент времени: Используется программа ПР - 3.4.

{$N+} { ПР - 3.4 } uses dos, crt, graph; const N=200; dt=0.0002; m=0.1; var Fx, Fy, x,y, vx, vy: array[1..N] of single; t, Gd, Gm, i,j: integer; k, ax, ay, F,l : single; hh: string; Procedure Sila; label Metka; begin For i:=1 to N do begin Fx[i]:=0; Fy[i]:=0; end; For i:=1 to N do For j:=1 to N do begin If j=i then goto Metka; l:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); If (l>0)and(l<10) then F:=40000/l else F:=0; Fx[i]:=Fx[i]+F*(x[i]-x[j])/(l+0.0001); Fy[i]:=Fy[i]+F*(y[i]-y[j])/(l+0.0001); Metka: end; end; Procedure Nach_uslov; begin Randomize; For i:=1 to N do begin x[i]:=10+random(390); y[i]:=10+random(390); vy[i]:=random(400)/10-20; vx[i]:=random(400)/10-20; end; end; BEGIN Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); Nach_uslov; Repeat Sila; For i:=1 to N do begin vx[i]:=vx[i]+Fx[i]/m*dt; vy[i]:=vy[i]+Fy[i]/m*dt; x[i]:=x[i]+vx[i]*dt; y[i]:=y[i]+vy[i]*dt; If y[i]>400 then y[i]:=10; If y[i]<10 then y[i]:=400; If x[i]<10 then x[i]:=400; If x[i]>400 then x[i]:=10; end; If t mod 20=0 then begin cleardevice; For i:=1 to N do circle(round(x[i]),round(y[i]),2); end; t:=t+1; str(t, hh); OuttextXY(550,150,hh); until KeyPressed; Repeat until keypressed; CloseGraph; END.

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

Будем моделировать жидкость как совокупность частиц (шариков), между которыми действуют силы притяжения (при удалении) и отталкивания (при сближении). Движение каждого шарика подчиняется второму закону Ньютона и рассчитывается как в задаче 3.4. На систему частиц действует сила тяжести, она как единое целое падает на выступ или горизонтальную поверхность. Используется программа ПР - 3.5, результат моделирования представлен на рис. 3.5.

{$N+} { ПР - 3.5 } uses dos, crt, graph; const N=250; dt=0.003; m=1; var Fx, Fy, x,y, vx, vy, xx, yy: array[1..N] of single; t, Gd, Gm, i,j: integer; k, ax, ay, F,l : single; hh: string; Procedure Sila; label Metka; begin For i:=1 to N do begin Fx[i]:=0; Fy[i]:=0; end; For i:=1 to N do For j:=1 to N do begin If j=i then goto Metka; l:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); If (l>9)and(l<30) then F:=-40000/l/l else F:=0; If l<8 then F:=4000000/l/l/l; Fx[i]:=Fx[i]+F*(x[i]-x[j])/(l+0.001); Fy[i]:=Fy[i]+F*(y[i]-y[j])/(l+0.001)+m*10; Metka: end; end; Procedure Nach_uslov; begin Randomize; For j:=0 to 24 do For i:=1 to 10 do begin x[i+10*j]:=8*i+170; y[i+10*j]:=8*j+150; vy[i]:=0; vx[i]:=0; end; end; BEGIN Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); Nach_uslov; Repeat Sila; For i:=1 to N do begin k:=0.9; xx[i]:=x[i]; yy[i]:=y[i]; ax:=Fx[i]/m; ay:=Fy[i]/m; vx[i]:=k*vx[i]+ax*dt; vy[i]:=k*vy[i]+ay*dt; x[i]:=x[i]+vx[i]*dt; y[i]:=y[i]+vy[i]*dt; If (x[i]<220)and(y[i]>352) then begin vx[i]:=-k*vx[i]; x[i]:=221; end; If (x[i]<220)and(y[i]>350) then begin vy[i]:=-k*vy[i]; y[i]:=349; end; If (y[i]>400) then begin vy[i]:=-k*vy[i]; y[i]:=399; end; end; cleardevice; setcolor(15); line(0,402,640,402); rectangle(0,402,318,353); For i:=1 to N do begin Circle(100+round(x[i]),round(y[i]),3); Circle(100+round(x[i]),round(y[i]),1); end; t:=t+1; str(t, hh); OuttextXY(150,150,hh); until KeyPressed; Repeat until keypressed; CloseGraph; END.

Рис. 3.5. Падение столба жидкости на выступ.

Задача 3.6. Используя решение предыдущей задачи создайте компьютерную модель диффузии двух газов. Считайте молекулы газов твердыми шарами, между которыми действуют силы отталкивания.

Для решения задачи следует промоделировать движение нескольких десятков молекул двух сортов внутри прямоугольного сосуда. Чтобы учесть отражение молекул от стенок в результате абсолютно упругого удара будем использовать операторы: If (x[i]>500)or(x[i]<100) then vx[i]:=-vx[i]; If (y[i]>200)or(y[i]<100) then vy[i]:=-vy[i]; Результат работы программы ПР - 3.6, моделирующей диффузию, представлен на рис. 3.6.

{$N+} { ПР - 3.6 } uses dos, crt, graph; const N=100; dt=0.00005; m=0.01; var Fx, Fy, x,y, vx, vy: array[1..N] of single; t, Gd, Gm, i,j: integer; k, ax, ay, F,l : single; hh: string; Procedure Sila; label Metka; begin For i:=1 to N do begin Fx[i]:=0; Fy[i]:=0; end; For i:=1 to N do For j:=1 to N do begin If j=i then goto Metka; l:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); If (l>0)and(l<10) then F:=40000/l else F:=0; Fx[i]:=Fx[i]+F*(x[i]-x[j])/(l+0.0001); Fy[i]:=Fy[i]+F*(y[i]-y[j])/(l+0.0001); Metka: end; end; Procedure Nach_uslov; begin Randomize; For i:=1 to N do begin y[i]:=102+random(98); If i<N/2 then x[i]:=102+random(200) else x[i]:=302+random(198); vx[i]:=random(800)/10-40; vy[i]:=random(800)/10-40; end; end; BEGIN Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); Nach_uslov; Repeat Sila; inc(t); For i:=1 to N do begin vx[i]:=vx[i]+Fx[i]/m*dt; vy[i]:=vy[i]+Fy[i]/m*dt; x[i]:=x[i]+vx[i]*dt; y[i]:=y[i]+vy[i]*dt; If (x[i]>500)or(x[i]<100) then vx[i]:=-vx[i]; If (y[i]>200)or(y[i]<100) then vy[i]:=-vy[i]; end; If t mod 10=0 then begin cleardevice; setbkcolor(white); setcolor(8); rectangle(98,98,502,202); For i:=1 to N do If i<N/2 then begin setcolor(9); circle(round(x[i]),round(y[i]),2); end else begin setcolor(2); circle(round(x[i]),round(y[i]),4); end; end; until KeyPressed; Repeat until keypressed; CloseGraph; END.

Рис. 3.6. Моделирование диффузии двух газов.

Задача 3.7. Промоделируйте обтекание пластины газом (сжимаемой жидкостью). Газ следует представлять в виде совокупности частиц-маркеров, между которыми действуют силы притяжения и отталкивания.

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

{$N+} { ПР - 3.7 } uses dos, crt, graph; const N=250; dt=0.01; var Fx, Fy, x,y, vx, vy : array[1..N] of real; m, t,Gd, Gm, i,j : integer; k, F,l : real; Procedure Sila; begin For i:=1 to N do begin Fx[i]:=0; Fy[i]:=0; end; For i:=1 to N do For j:=1 to N do begin l:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); If (l>0)and(l<100) then F:=-700/(l) else F:=0; If (l>0)and(l<10) then F:=25000/l; Fx[i]:=Fx[i]+F*(x[i]-x[j])/(l+0.001); Fy[i]:=Fy[i]+F*(y[i]-y[j])/(l+0.001)+m*20; end; end; Procedure Nach_uslov; begin Randomize; m:=1; For j:=0 to 24 do For i:=1 to 10 do begin x[i+10*j]:=10*i+155; y[i+10*j]:=10*j+120; end; end; BEGIN Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); setbkcolor(15); Nach_uslov; Repeat Sila; For i:=1 to N do begin k:=0.2; vx[i]:=k*vx[i]+Fx[i]/m*dt; vy[i]:=k*vy[i]+Fy[i]/m*dt; x[i]:=x[i]+vx[i]*dt; y[i]:=y[i]+vy[i]*dt; If (x[i]<160) then begin vx[i]:=-vx[i]; x[i]:=161; end; If (x[i]>260) then begin vx[i]:=-vx[i]; x[i]:=259; end; If (x[i]>180)and(x[i]<240)and(y[i]>320)and(y[i]<324) then begin vy[i]:=-vy[i]; y[i]:=319; end; If (y[i]>410) then begin vy[i]:=1; y[i]:=240; x[i]:=random(90)+165; end; end; delay(50); If t mod 2=0 then begin cleardevice; setcolor(12); For i:=1 to N do begin circle(round(2*x[i]),round(2*y[i])-400,3); end; setcolor(8); rectangle(2*180,2*320-400,2*240,2*324-400); line(2*160,0,2*160,480); line(2*260,0,2*260,480); end; inc(t); until KeyPressed; Repeat until keypressed; CloseGraph; END.

Рис. 3.7. Обтекание пластины газом.

Задача 3.8. Промоделируйте конвекцию газа в прямоугольном сосуде с помощью частиц-маркеров, между которыми действуют силы отталкивания.

Допустим, квадратный сосуд заполнен двумерным газом, который моделируется совокупностью частиц-маркеров. Эти частицы соответствуют элементарным объемам газа, они имеют равные массы, координаты x[i] и y[i], скорости vx[i] и vy[i], и температуры t[i]. Между частицами-маркерами действуют силы отталкивания, вычисляемые в процедуре Sila. Левая часть сосуда нагревается, температура частиц с x[i]<50 равна t[i]:=1, на них действует сила, направленная вверх. Правая часть сосуда охлаждается, температура частиц с x[i]>150 составляет t[i]:=-1, на них действует сила, направленная вниз. Для того, чтобы промоделировать конвекцию газа, необходимо учесть, что давление в точках, лежащих на одной горизонтали, примерно одинаково. То есть молекулы (а значит и частицы-маркеры), поднимающиеся вблизи левой стенки вверх, не должны скапливаться в левом верхнем углу сосуда. Для учета этого фактора сосуд разбивается на 16 квадратов. В процедуре Davlenie рассчитывается число частиц-маркеров в каждом квадрате, считается, что оно пропорционально давлению газа. Если в левом квадрате давление меньше чем в данном, то на частицы данного квадрата действует сила, направленная влево. Если в правом квадрате давление меньше, то на частицы действует сила, направленная вправо. Это учитывается в процедуре Sila. Используется программа PR-1. На рис. 3.8 показаны траектории частиц-маркеров за небольшой промежуток времени. Левый рисунок соответсвует ситуации, когда частицы движутся хаотично, конвекция отсутствует.

{$N+} { ПР - 3.8 } uses dos, crt, graph; const N=80; dt=0.0005; k=0.95; m=0.1; var Fx, Fy, x,y, vx, vy, xx, yy : array[1..N] of single; t: array[1..N]of single; p: array[1..10,1..10] of integer; FFy, FFx, Gd, Gm, i,j, kk, v,w, v1,w1,time: integer; r, F,l: single; Procedure Sila; begin For i:=1 to N do begin Fx[i]:=0; Fy[i]:=0; end; For i:=1 to N do For j:=1 to N do begin l:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); If (l>0)and(l<30) then F:=300000/l else F:=0; v1:=round(x[i]/50); w1:=round(y[i]/50); FFx:=0; If ((p[v1,w1]>p[v1-1,w1])and(v1>2)) then FFx:=70; If ((p[v1,w1]>p[v1+1,w1])and(v1<3)) then FFx:=-70; Fx[i]:=Fx[i]+F*(x[i]-x[j])/(l+0.001)+FFx; Fy[i]:=Fy[i]+F*(y[i]-y[j])/(l+0.001)-100*t[i]; end; end; Procedure Davlenie; begin For v:=1 to 4 do For w:=1 to 4 do p[v, w]:=0; For i:=1 to N do begin For v:=1 to 4 do For w:=1 to 4 do If (50*(v-1)<x[i])and(x[i]<50*v)and(50*(w-1)<y[i])and(y[i]<50*w) then inc(p[v, w]); end; end; BEGIN Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); For i:=1 to N do begin x[i]:=random(198); y[i]:=random(198); end; Repeat Davlenie; Sila; For i:=1 to N do begin vx[i]:=k*vx[i]+Fx[i]/m*dt; vy[i]:=k*vy[i]+Fy[i]/m*dt; x[i]:=x[i]+vx[i]*dt; y[i]:=y[i]+vy[i]*dt; t[i]:=0.8*t[i]; If x[i]>150 then t[i]:=-1; If x[i]<50 then t[i]:=1; If x[i]<0 then begin x[i]:=1; vx[i]:=-vx[i]; end; If x[i]>200 then begin x[i]:=199; vx[i]:=-vx[i]; end; If y[i]<0 then begin y[i]:=1; vy[i]:=-vy[i]; end; If y[i]>200 then begin y[i]:=199; vy[i]:=-vy[i]; end; end; delay(10); inc(time); rectangle(198,198,402,402); If time mod 500=0 then cleardevice; If time mod 20<>0 then For i:=1 to N do circle(200+round(x[i]),200+round(y[i]),1); until KeyPressed; Repeat until keypressed; CloseGraph; END.