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

Рис. 5.4. Электрическая цепь постоянного тока.

Используя законы Кирхгофа составим систему уравнений (рис.01). Подставим вместо сопротивлений и ЭДС конкретные значения, перенесем все члены в левую часть: Для решения системы уравнений с n неизвестными используется программа ПР - 5.3. На экран выводятся значения переменных и осществляется проверка: в исходную систему подставляются найденные значения xi, i=1, 2, ... ,n. При этом левые части уравнений обращаются в ноль, что подтверждает правильность решения. Итак, полученные значения токов: I1=0,481 А, I2=0,812 А, I3=-0,331 А, I4=-0,086 А, I5=-0,395 А, I6=0,726 А. Знак минус означает, что реальное направление тока в соответствующей ветви противоположно предполагаемому.

Задача 5.5. Цепь состоит из источника переменной ЭДС и трех ветвей, в каждой из которых резистор, конденсатор и катушка индуктивности. Вторая и третья ветви соединены параллельно между собой, последовательно с ними включена первая ветвь. Рассчитайте все токи и напряжения, полную, активную и реактивную мощности. Постройте векторную диаграмму. Определите действующие значения всех токов и напряжений.

Рис. 5.5. Однофазная цепь переменного тока.

Допустим, цепь имеет параметры: R1=10 Ом, L1=0.08 Гн, C1=1200 мкФ, R2=8 Ом, L2=0.42 Гн, C2=1500 мкФ, R3=14 Ом, L3=0.14 Гн, C3=1270 мкФ. Для вычисления импедансов, токов и напряжений для каждой ветви используются формулы: Представленная ниже программа ПР - 5.5 содержит процедуры Sum, Razn, Proizv, Delen, Modul, позволяющие определять сумму, разность, произведение, отношение двух комплексных чисел и находить модуль комплексного числа. С помощью этих процедур осуществляются расчеты комплексов тока и напряжения, комплекса мощности.

uses crt, graph; { ПР - 5.5 } const n=4; type komplex = record Re : real; Im : real; end; var I1,Iba, Icb, Iac, E,EB, EC, z1,z2,z3,S, R,P, D,U, Z,C : komplex; var i, ii, j,jj, k : integer; w, m:real; aa, a: array[1..N, 1..N+1] of komplex; x, bb, b: array[1..N] of komplex; Procedure Sum(Z1, Z2: komplex); begin S. Re:=Z1.Re+Z2.Re; S. Im:=Z1.Im+Z2.Im; end; Procedure Razn(Z1, Z2: komplex); begin R. Re:=Z1.Re-Z2.Re; R. Im:=Z1.Im-Z2.Im; end; Procedure Proizv(Z1, Z2: komplex); begin P. Re:=Z1.Re*Z2.Re-Z1.Im*Z2.Im; P. Im:=Z1.Re*Z2.Im+Z2.Re*Z1.Im; end; Procedure Delen(Z1, Z2: komplex); begin D. Re:=(Z1.Re*Z2.Re+Z1.Im*Z2.Im)/(Z2.Re*Z2.Re+Z2.Im*Z2.Im); D. Im:=(Z1.Im*Z2.Re-Z1.Re*Z2.Im)/(Z2.Re*Z2.Re+Z2.Im*Z2.Im); end; Procedure Modul(Z1: komplex); begin M:=sqrt(Z1.Re*Z1.Re+Z1.Im*Z1.Im); end; BEGIN clrscr; w:=2*3.1415926*50; Z1.Re:=10; Z1.Im:=w*0./(w*1200); Z2.Re:=8; Z2.Im:=w*0./(w*1500); Z3.Re:=14; Z3.Im:=w*0./(w*1270); writeln('Z1= ', Z1.Re:2:3,' ',Z1.Im:2:3,'j'); writeln('Z1= ', Z2.Re:2:3,' ',Z2.Im:2:3,'j'); writeln('Z1= ', Z3.Re:2:3,' ',Z3.Im:2:3,'j'); Proizv(Z2,Z3); Sum(Z2,Z3); Delen(P, S); Sum(Z1,D); Modul(S); writeln('Z = ',S. Re:2:3,' ',S. Im:2:3,'j ', M:2:3); E. Re:=127; E. Im:=0; Delen(E, S); I1:=D; Modul(D); writeln('I1= ', D. Re:2:3,' ',D. Im:2:3,'j ', M:2:3); Proizv(D, Z1); Modul(P); writeln('U1= ', P. Re:2:3,' ',P. Im:2:3,'j ', M:2:3); Razn(E, P); Modul(R); writeln('U2= ', R. Re:2:3,' ',R. Im:2:3,'j ', M:2:3); Delen(R, Z2); Modul(D); writeln('I2= ', D. Re:2:3,' ',D. Im:2:3,'j ', M:2:3); Delen(R, Z3); Modul(D); writeln('I3= ', D. Re:2:3,' ',D. Im:2:3,'j ', M:2:3); I1.Im:=-I1.Im; Proizv(E, I1); writeln('S = ', P. Re:2:3,' ', P. Im:2:3,'j'); ReadKey; END.

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

Рис. 5.6. Трехфазная цепь.

Изучение трехфазных цепей предполагает нахождение комплексных значений токов, что требует решения системы уравнений в комплексных числах. По законам Кирхгофа составим систему уравнений в комплексных числах. Цепь имеет три независимых контура и два узла O и O', всего 4 ветви; получаем систему из 4 независимых уравнений: Эта системы уравнений может быть решена путем приведения матрицы коэффициентов к треугольному виду. В используемой программе ПР - 5.6 операции сложения, вычитания, умножения и деления комплексных чисел осуществляются в процедурах Sum, Razn, Proizv, Delen. Сначала, чтобы исключить деление на 0, складываются второе, третье и последующие уравнения с первым. Делят элементы каждой строки матрицы на первый элемент, после чего вычитают из второй, третьей и т. д. строк первую. В результате первые элементы второй, третьей и последующих строк оказываются равными 0. Аналогичным образом приводят к 0 вторые элементы третьей, четвертой и последующих строк и т. д. Эти операции повторяют до тех пор, пока матрица не будет приведена к треугольному виду. Это позволяет вычислить искомые значения комплексов токов и осуществить проверку решения. Эта программа может быть использована для решения системы уравнений в действительных числах.

	uses crt, graph; { ПР - 5.6 } const n=4; type komplex = record Re : real; Im : real; end; var z1,z2,z3,S, R,P, D,U, Z,C :komplex; var i, ii, j,jj, k : integer; m:real; aa, a: array[1..N, 1..N+1] of komplex; x, bb, b: array[1..N] of komplex; Procedure Sum(Z1, Z2: komplex); begin S. Re:=Z1.Re+Z2.Re; S. Im:=Z1.Im+Z2.Im; end; Procedure Razn(Z1, Z2: komplex); begin R. Re:=Z1.Re-Z2.Re; R. Im:=Z1.Im-Z2.Im; end; Procedure Proizv(Z1, Z2: komplex); begin P. Re:=Z1.Re*Z2.Re-Z1.Im*Z2.Im; P. Im:=Z1.Re*Z2.Im+Z2.Re*Z1.Im; end; Procedure Delen(Z1, Z2: komplex); begin D. Re:=(Z1.Re*Z2.Re+Z1.Im*Z2.Im)/(Z2.Re*Z2.Re+Z2.Im*Z2.Im); D. Im:=(Z1.Im*Z2.Re-Z1.Re*Z2.Im)/(Z2.Re*Z2.Re+Z2.Im*Z2.Im); end; Procedure Modul(Z1: komplex); begin M:=sqrt(Z1.Re*Z1.Re+Z1.Im*Z1.Im); end; BEGIN clrscr; a[1,1].Re:=-1;a[1,2].Re:=1; a[1,3].Re:=1; a[1,4].Re:=1; a[1,5].Re:=0; a[2,1].Re:=10;a[2,2].Re:=100; a[2,3].Re:=0; a[2,4].Re:=0; a[2,5].Re:=-220; a[3,1].Re:=10;a[3,2].Re:=0; a[3,3].Re:=0; a[3,4].Re:=0; a[3,5].Re:=110; a[4,1].Re:=10;a[4,2].Re:=0; a[4,3].Re:=0; a[4,4].Re:=0; a[4,5].Re:=110; a[1,1].Im:=0; a[1,2].Im:=0; a[1,3].Im:=0; a[1,4].Im:=0; a[1,5].Im:=0; a[2,1].Im:=0; a[2,2].Im:=0; a[2,3].Im:=0; a[2,4].Im:=0; a[2,5].Im:=0; a[3,1].Im:=0; a[3,2].Im:=0; a[3,3].Im:=-120; a[3,4].Im:=0; a[3,5].Im:=-110*1.73; a[4,1].Im:=0; a[4,2].Im:=0; a[4,3].Im:=0; a[4,4].Im:=80; a[4,5].Im:=+110*1.73; For i:=1 to N do For j:=1 to N+1 do aa[i, j]:=a[i, j]; For i:=2 to N do For j:=1 to N+1 do begin Sum(a[i, j],a[1,j]); a[i, j]:=S; end; For i:=3 to N do For j:=1 to N+1 do begin Sum(a[i, j],a[2,j]); a[i, j]:=S; end; For i:=4 to N do For j:=1 to N+1 do begin Sum(a[i, j],a[3,j]); a[i, j]:=S; end; For k:=1 to N do begin For i:=k to N do For j:=1 to N+1 do begin Delen(a[i, N+2-j],a[i, k]); a[i, N+2-j]:=D; end; For i:=k+1 to N do For j:=1 to N+1 do begin Razn(a[i, j],a[k, j]); a[i, j]:=R; end; end; For i:=1 to N do begin For j:=1 to N+1 do write(a[i, j].Re:2:3, ' '); writeln; end; writeln; For i:=1 to N do begin For j:=1 to N+1 do write(a[i, j].Im:2:3, ' '); writeln; end; writeln; For ii:=1 to N do begin i:=N-ii+1; x[i].Re:=-a[i, N+1].Re; x[i].Im:=-a[i, N+1].Im; For j:=i+1 to N do begin Proizv(a[i, j],x[j]); Razn(x[i],P); x[i]:=R; end; end; writeln(' ОТВЕТ '); For i:=1 to N do begin Modul(x[i]); writeln('x[',i,']= ',x[i].Re,' ',x[i].Im, M); end; writeln(' ПРОВЕРКА '); For i:=1 to N do begin Proizv(aa[i,1],x[1]); C:=P; Proizv(aa[i,2],x[2]); Sum(P, C); C:=S; Proizv(aa[i,3],x[3]); Sum(P, C); C:=S; Proizv(aa[i,4],x[4]); Sum(P, C); C:=S; Sum(aa[i,5],C); writeln(S. Re,' ',S. Im); end; ReadKey; END.
6. РАСЧЕТ ЭЛЕКТРИЧЕСКОГО И МАГНИТНОГО ПОЛЕЙ

Задача 6.1. Три точечных заряда q1, q2, q3 имеют координаты (x1, y1), (x2, y2), (x3, y3). Рассчитайте потенциал электростатического поля во всех точках плоскости, содержащей эти заряды.

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

uses crt, graph; { ПР - 6.1.1 } var r1,r2,r3,x1,x2,x3,y1,y2,y3,fi, q1,q2,q3: real; Gd, Gm, i,j :integer; BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi\'); q1:=15; q2:=-20; q3:=10; q1:=8; q2:=46; x1:=180.5; y1:=150.5; x2:=400.5; y2:=180.5; x3:=240.5; y3:=340.5; For i:=0 to 640 do For j:=0 to 480 do begin r1:=sqrt(sqr(x1-i)+sqr(y1-j)); r2:=sqrt(sqr(x2-i)+sqr(y2-j)); r3:=sqrt(sqr(x3-i)+sqr(y3-j)); fi:=q1/r1+q2/r2+q3/r3; putpixel(i, j,round((fi+1)*13)); {If round(fi/0.05)-fi/0.05<0.002 then putpixel(i, j,15);} end; Repeat until Keypressed; CloseGraph; END.

Рис. 6.1.1. Расчет потенциала поля точечных зарядов.

Рассмотренный метод может быть использован для расчета полей, созданных протяженными заряженными телами. Например, чтобы промоделировать электростатическое поле плоской пластины можно вместо нее расположить 10-20 точечных зарядов, лежащих на прямой. Программа ПР - 6.1.2 рассчитывает поле пластины и расположенного рядом с ней точечного заряда. Результат вычислений приведен на рис. 6.1.2.

uses crt, graph; const N=40; var fi : real; Gd, Gm, i,j, k :integer; q, x, y, r : array[1..N]of real; BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi\'); For k:=1 to 20 do begin q[k]:=10; y[k]:=100; x[k]:=200+20*k; end; {For k:=21 to 40 do begin q[k]:=-10; y[k]:=200; x[k]:=100+20*(k-20); end;} q[k]:=-40; y[k]:=200; x[k]:=400; For i:=0 to 640 do For j:=0 to 480 do begin fi:=0; For k:=1 to N do begin r[k]:=sqrt(sqr(x[k]-i)+sqr(y[k]-j)); If r[k]<>0 then fi:=fi+q[k]/r[k]; end; putpixel(i, j,round((fi+1)*13)); end; Repeat until Keypressed; CloseGraph; END.

Рис. 6.1.2. Расчет потенциала поля точечных зарядов.

Задача 6.2. Две пластины имеют потенциалы 100 В и -100 В. Между ними вблизи краев находится металлический шарик с потенциалом 20 В. Рассчитайте распределение потенциала.

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

{$N+} uses crt, graph; { ПР - 6.2 } const N=100; M=100; NN=200; h=0.2; dt=0.1; SizeOfReal=6; type RealPoint= ^real; var t: single; i, j,ii, jj, DV, MV, EC: integer; a: array[1..NN] of pointer; vv, rho, kk, k : real; Function fi(i, j: word) : RealPoint; begin fi:=Ptr(Seg(a[i]^),ofs(a[i]^)+(j-i)*SizeOfReal); end; Procedure Raschet; begin vv:=1; k:=0.1; For i:=2 to N-1 do for j:=2 to M-1 do begin If j<M then begin If (i=80)and(j>30) then fi(i, j)^:=100; If (i=20)and(j>30) then fi(i, j)^:=-100; end; If sqr(i-60)+sqr(j-70)<100 then fi(i, j)^:=20; fi(i+n, j)^:=(fi(i-1,j)^+fi(i+1,j)^+ fi(i, j-1)^+fi(i, j+1)^+rho*h*h)/4; end; For j:=2 to M-1 do for i:=2 to N-1 do fi(i, j)^:=fi(i+n, j)^; For j:=1 to M do fi(1,j)^:=0; For j:=1 to M do fi(N, j)^:=0; For i:=1 to N do fi(i,2)^:=0; For i:=1 to N do fi(i, M)^:=fi(i, M-1)^; end; Procedure Draw; begin for i:=2 to N-1 do for j:=2 to M-1 do begin setcolor(round(fi(i, j)^/8)+1); If sqr(i-60)+sqr(j-70)<100 then setcolor(15); rectangle(11+i*4,450-j*4,13+i*4,453-j*4); rectangle(10+i*4,450-j*4,12+i*4,452-j*4); end; end; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); For i:=1 to NN do begin GetMem(a[i], M*SizeOfReal); For j:=1 to M do fi(i, j)^:=0; end; Repeat t:=t+dt; Raschet; kk:=kk+1; If kk=round(kk/25)*25 then Draw; until KeyPressed; Release(HeapOrg); CloseGraph; END.

Рис. 6.2. Моделирование электростатического поля.

Задача 6.3. Постройте силовые линии магнитного поля, созданного: 1) одной обмоткой с током; 2) двумя обмотками, по которым текут токи в одном направлении или в противоположных направлениях.

Рассмотрим виток с током радиусом r, лежащий в плоскости XOY. Разобъем его на элементы и введем вектора сонаправленные с током. Расположение этого элемента витка и точки A, в которой определяется индукция магнитного поля, задаются векторами: По закону Био-Савара-Лапласа для вектора индукции и его проекций можно записать: При этом используются формулы: Используется программа ПР - 6.3. Для построения силовой линии задают начальную точку A0, определяют в ней направление индукции магнитного поля, строят небольшой отрезок в направлении вектора индукции, переходя в новую точку A1. Затем повторяют все снова. Результаты вычислений -- на рис. 6.3.1, 6.3.2, 6.3.3.

Uses crt, graph; { ПР - 6.3 } Const N=200; pi=3.1415; Var S, k,kk, nn, ii, jj, i,j, DV, MV, r,r1,TOK: integer; Rx, Ry, Rz, BBB: real; Bx, By, Bz, dfi, deltax, deltay, deltaz, deltaR, dlx, dly, dlz, alpha: real; Procedure Raschet; begin For i:=1 to N do begin dfi:=2*pi/N; dlx:=TOK*r*sin(dfi*i)*dfi; dly:=-TOK*r*cos(dfi*i)*dfi; deltax:=r*cos(dfi*i)-Rx; deltay:=r*sin(dfi*i)-Ry; deltaR:=sqrt(deltax*deltax+deltay*deltay+deltaz*deltaz); Bz:=(dlx*deltay-dly*deltax)/(deltaR*deltaR*deltaR)+Bz; By:=(dlz*deltax-dlx*deltaz)/(deltaR*deltaR*deltaR)+By; end; end; Procedure Postroenie; begin Repeat Bz:=0; By:=0; r:=80; For k:=0 to 3 do begin deltaz:=50*k-Rz; TOK:=1; Raschet; end; r:=50; For kk:=0 to 1 do begin deltaz:=30*kk-Rz+300; TOK:=-1; Raschet; end; BBB:=sqrt(By*By+Bz*Bz); circle(320+round(Ry),450-round(Rz),1); Ry:=Ry+By/BBB; Rz:=Rz+Bz/BBB; nn:=nn+1; For k:=0 to 3 do begin circle(320+80,450-50*k,3); circle(320-80,450-50*k,3); end; For kk:=0 to 1 do begin circle(320+50,450-30*kk-300,3); circle(320-50,450-30*kk-300,3); end; until (KeyPressed)or(nn>300); end; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); Rx:=0; For ii:=-4 to 5 do begin Ry:=-5+10*ii; Rz:=20; nn:=0; Postroenie; end; For ii:=-3 to 4 do begin Ry:=-5+10*ii; Rz:=350; nn:=0; Postroenie; end; Repeat until KeyPressed; CloseGraph; END.

Рис. 6.3.1. Магнитное поле обмотки.

Рис. 6.3.2. Магнитное поле двух обмоток. Токи сонаправлены.

Рис. 6.3.3. Магнитное поле двух обмоток. Токи противоположнонаправлены.

Задача 6.4. Постройте силовые линии магнитного поля, созданного тремя параллельными проводниками с током.

Рис. 6.4.1. Расчет магнитного поля двух токов.

Допустим, два проводника с токами I1 и I2 пересекают плоскость XOY в точках с координатами (x1,y1) и (x2,y2). Индукция магнитного поля в точке A(x, y) может быть найдена по формулам: Используемая программа ПР-1 случайным образом выбирает точку A0, определяет проекции вектора индукции и в его направлении строит отрезок единичной длины, находя точку A1. Затем все повторяется снова. Выполнив 10000 шагов, программа снова случайно выбирает точку A0 и начинает строить вторую силовую линию и т. д.

Uses crt, graph; Var DrV, MV, i,x1,y1,x2,y2,x3,y3,k : integer; cosa1,sina1,cosa2,sina2,cosa3,sina3,I1,I2,I3, B, B1,B2,B3,Bx, By, r1,r2,r3,x, y,pi: real; BEGIN DrV:=Detect; InitGraph(DrV, MV,'c:\bp\bgi'); pi:=arctan(1)*4; x1:=150; y1:=100; x2:=400; y2:=100; x3:=300; y3:=300; circle(x1,y1,2); circle(x2,y2,2); circle(x3,y3,2); x:=300; I1:=1; I2:=-2; I3:=2; Randomize; For i:=1 to 80 do begin k:=0; x:=round(random(64))*10; y:=round(random(48))*10; Repeat k:=k+1; r1:=sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1)); r2:=sqrt((x-x2)*(x-x2)+(y-y2)*(y-y2)); r3:=sqrt((x-x3)*(x-x3)+(y-y3)*(y-y3)); B1:=I1/r1; cosa1:=(x-x1)/r1; sina1:=(y-y1)/r1; B2:=I2/r2; cosa2:=(x-x2)/r2; sina2:=(y-y2)/r2; B3:=I3/r3; cosa3:=(x-x3)/r3; sina3:=(y-y3)/r3; Bx:=B1*sina1+B2*sina2+B3*sina3; By:=B1*cosa1+B2*cosa2+B3*cosa3; B:=sqrt(Bx*Bx+By*By); x:=x+0.1*Bx/B; y:=y-0.1*By/B; putpixel(round(x),round(y),15); until (k>10000)or KeyPressed; end; Repeat until KeyPressed; CloseGraph; END.

Рис. 6.4.2. Магнитное поле трех проводников с током.

Задача 6.5. Создайте компьютерную модель ферромагнетика Изинга. Промоделируйте процесс образования областей спонтанной намагниченности (доменов).

Рассмотрим двумерную модель магнетика. Представим себе сетку, в каждом из узлов которой находится атом, спин (собственный магнитный момент) которого может быть направлен от нас, к нам или лежать в плоскости рисунка. Энергия взаимодействия соседних атомов друг с другом - Jsi, jsi+1,j, атомы, расположенные далеко друг от друга, не взаимодействуют. Внешнее магнитное поле способствует повороту спина атома и изменяет энергию системы на H или - H. Энергия всей системы равна: Процесс моделирования состоит в следующем: 1) задается начальная ориентация всех спинов случайным образом; 2) определяется энергия системы; 3) случайно выбирается один атом и его спин изменяется на противоположный, определяется энергия системы в новой конфигурации; 4) если в результате переориентации спина энергия уменьшилась, то такое изменение конфигурации принимается, в противном случае -- отклоняется; 5) результат выводится на экран; 6) возврат к операции 2. Для решения этой задачи используется программа ПР - 6.5, при этом моделируется ситуация, в которой в левом верхнем углу индукция магнитного поля направлена от нас, а в правом нижнем углу -- к нам. Используется не квадратная, а шестиугольная сетка, результат вычислений -- на рис. 6.5.

{$N+} Uses crt, graph; { ПР - 6.5 } Const N=50; M=50; J_vzaim=1; Var i, j,k, l,u, d,Gd, Gm : integer; p, H,E, EE, dE : real; s : array [0..N+1,0..M+1] of single; Energiya : string; Label metka; Procedure Energy; Var i, j : integer; begin E:=0; dE:=0; For i:=1 to N do For j:=1 to M do begin If j mod 2=0 then E:=E-J_vzaim*(s[i, j]*s[i+1,j]+s[i, j]*s[i-1,j]+s[i, j]*s[i, j+1]+ s[i, j]*s[i, j-1]+s[i, j]*s[i-1,j+1]+s[i, j]*s[i-1,j-1]); If j mod 2<>0 then E:=E-J_vzaim*(s[i, j]*s[i+1,j]+s[i, j]*s[i-1,j]+s[i, j]*s[i, j+1]+ s[i, j]*s[i, j-1]+s[i, j]*s[i+1,j+1]+s[i, j]*s[i+1,j-1]); end; For i:=1 to N do For j:=1 to M do begin h:=0; If (i>0.5*N)and(j>0.5*M) then H:=-2.5; If (i<0.4*N)and(j<0.4*M) then H:=1.5; dE:=dE-H*s[i, j]; end; E:=E+dE; end; Procedure Draw; Var i, j : integer; begin cleardevice; setbkcolor(white); For i:=1 to N do For j:=1 to M do begin If j mod 2=0 then d:=0 else d:=3; If s[i, j]=1 then setcolor(12); If s[i, j]=-1 then setcolor(9); If s[i, j]=0 then setcolor(15); circle(20+6*i+d,20+6*j,1); circle(20+6*i+d,20+6*j,2); circle(20+6*i+d,20+6*j,3); end; Str(E, Energiya); OuttextXY(10,400,Energiya); end; BEGIN Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi'); If GraphResult <> grOk then Halt(1); Randomize; For i:=1 to N do For j:=1 to M do begin p:=Random(100)/100; s[i, j]:=0; If p<0.33 then s[i, j]:=-1; If p>0.66 then s[i, j]:=1; end; Draw; Repeat For k:=1 to round(N*M/5) do begin Energy; EE:=E; i:=round(random(N)); j:=round(random(M)); If s[i, j]=1 then begin u:=1; s[i, j]:=0; goto metka; end; If s[i, j]=-1 then begin u:=-1; s[i, j]:=0; goto metka; end; If (s[i, j]=0)and(random(100)>50) then s[i, j]:=-1 else s[i, j]:=1; u:=0; metka: Energy; If (E>EE) then s[i, j]:=u; end; Energy; Draw; until KeyPressed; Repeat until KeyPressed; CloseGraph; END.

Рис. 6.5. Модель ферромагнетика Изинга.

7. МОДЕЛИРОВАНИЕ ТЕПЛОПРОВОДНОСТИ

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

Уравнение теплопроводности для двумерной среды в конечных разностях имеет вид: Для расчета распределения температуры по поверхности пластины используется программа ПР - 7.1. Она содержит цикл по времени с вложенными в него двумя циклами по i и по j, в которых перебираются все элементы пластины и вычисляются их температуры в последующие моменты времени. Результат решения задачи представлен на рис. 7.1, -- на нем разными цветами изображены области с различными температурами. Изотермы (границы разноцветных областей) перепендикулярны теплоизолированным краям пластины и параллельны краям, температура которых поддерживается постоянной.

Uses crt, graph; { ПР - 7.1. } Const n=75; m=70; h=1; dt=0.1; Var kk, i, j, DV, MV : integer; tt, t : array[1..N, 1..M] of real; q, a : real; uslovie : boolean; Procedure Istoch; begin If ((i>20)and(i<35))and((j>20)and(j<30)) then q:=20 else q:=0; If ((i>45)and(i<65))and((j>50)and(j<60)) then q:=-10; end; Procedure Nach_uslov; begin For i:=1 to N do For j:=1 to M do t[i, j]:=1; end; Procedure Raschet; begin Istoch; If (abs(i-40)<10)and(abs(j-40)<10) then a:=0 else a:=1; tt[i, j]:=t[i, j]+a*(t[i, j+1]-2*t[i, j]+t[i, j-1])*dt/(h*h)+ a*(t[i+1,j]-2*t[i, j]+t[i-1,j])*dt/(h*h)+q; end; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); Nach_uslov; Repeat kk:=kk+1; For i:=2 to N-1 do For j:=2 to M-1 do begin t[N, j]:=t[N-1,j]; Raschet; end; For i:=2 to N-1 do For j:=2 to M-1 do t[i, j]:=tt[i, j]; If kk/25=round(kk/25) then For i:=2 to N-1 do For j:=2 to M-1 do begin setcolor(round(t[i, j]/300+1)); rectangle(i*5+50,j*5,i*5+54,j*5+4); rectangle(i*5+51,j*5+1,i*5+53,j*5+3); end; until KeyPressed; CloseGraph; END.

Рис. 7.1. Распределение температуры: двумерная среда.

Задача 7.2. Пластина состоит из трех полосок с различными коэффициентами теплопроводности. Нижняя сторона пластины теплоизолирована, остальные поддерживаются при постоянной температуре. В разных местах пластины расположены протяженный источник тепла и источник холода. Необходимо рассчитать распределение температуры в последовательные моменты времени.

Заменяя частные производные их конечно-разностными аппроксимациями, запишем уравнение теплопроводности для неоднородной двумерной среды в конечных разностях: Для решения задачи используется программа ПР - 7.2. В ней последовательно перебираются узлы двумерной сетки по строкам и по столбцам. Когда переменная naprav равна 0, то элементы перебираются по столбцам сверху вниз и снизу вверх. Когда переменная naprav равна 1, то элементы перебираются по строкам слева направо и справа налево. Нижникй край пластины теплоизолирован, это задается циклом: For i:=1 to N do t[i, M]:=t[i, M-1]; Все остальные края пластины поддерживаются при постоянной температуре. Используется программа ПР - 7.2, результат ее работы представлен на рис. 7.2.

Uses crt, graph; { ПР - 7.2. } Const n=72; m=72; h=1; dt=0.05; Var naprav, ii, jj, kk, i,j, DV, MV, EC : integer; k, t : array[1..N, 1..M] of real; q, a,b, bb : real; uslovie : boolean; Procedure Raschet; begin q:=0; If (i>20)and(i<50)and(j>50)and(j<60) then q:=20; If (i>50)and(i<60)and(j>20)and(j<30) then q:=-30; If naprav=1 then t[i, j]:=t[i, j]+(k[i, j+1]-k[i, j-1])*(t[i, j+1]-t[i, j-1])*dt/(4*h*h)+ k[i, j]*(t[i, j-1]-2*t[i, j]+t[i, j+1])*dt/(h*h)+q*dt else t[i, j]:=t[i, j]+(k[i+1,j]-k[i-1,j])*(t[i+1,j]-t[i-1,j])*dt/(4*h*h)+ k[i, j]*(t[i-1,j]-2*t[i, j]+t[i+1,j])*dt/(h*h)+q*dt; end; BEGIN DV:=Detect; InitGraph(DV, MV,'c:\bp\bgi'); For i:=1 to N do For j:=1 to M do begin k[i, j]:=1; If (i>55) then k[i, j]:=2; If (i<30) then k[i, j]:=0.3; end; Repeat kk:=kk+1; naprav:=0; For i:=2 to N-1 do For j:=2 to M-1 do Raschet; For i:=2 to N-1 do For j:=M-1 downto 2 do Raschet; naprav:=1; For j:=2 to M-1 do For i:=2 to N-1 do Raschet; For j:=2 to M-1 do For i:=N-1 downto 2 do Raschet; For i:=1 to N do t[i, M]:=t[i, M-1]; If kk/10=round(kk/10) then For i:=2 to N-1 do For j:=2 to M-1 do begin setcolor(round(t[i, j]/100)+1); rectangle(i*5+50,j*5,i*5+54,j*5+4); rectangle(i*5+51,j*5+1,i*5+53,j*5+3); end; Line(55*5+52,0,55*5+52,480); Line(30*5+52,0,30*5+52,480); until KeyPressed; CloseGraph; END.