Создание модели возникновения Солнечной системы из межзвездного газа на базе численного моделирования с учетом гравитационного взаимодействия частиц
Курсовой проект - Авиация, Астрономия, Космонавтика
Другие курсовые по предмету Авиация, Астрономия, Космонавтика
vxcum,vycum:real;:string;:text;:char;
begin(число частиц=);(N);(размеры ящика Lx и Ly=);(Lx,Ly);(шаг по времени=);
readln(dt);:=dt*dt;
write(число шагов по времени между усреднениями=);(nave);(количество наборов усреднения=);(nset);(новая конфигурация?(y/n));
readln(newconf);(newconf=y) then
write(число частиц в ряду=);(nrow);(максимальная скорость=);
readln(vmax);:=N div nrow;:=Ly/nrow;:=Lx/ncol;:=0;icol:=1 to ncol doirow:=1 to nrow do
i:=i+1;[i]:=ay*(irow-0.5);
if((irow mod 2)=0) then[i]:=ax*(icol-0.25)[i]:=ax*(icol-0.75);[i]:=random*vmax;[i]:=random*vmax;:=vxcum+vx[i];:=vycum+vy[i];;:=vxcum/N;:=vycum/N;i:=1 to N do
vx[i]:=vx[i]-vxcum;[i]:=vy[i]-vycum;
end;
{write(имя файла с конфигурацией);
readln(fname);}(относительное изменение скорости=);
readln(vscale);(datain,novij.txt);(datain);(datain,N,Mx,My);:=Lx/Mx;:=Ly/My;i:=1 to N do(datain,x[i],y[i],vx[i],vy[i]);[i]:=x[i]*xscale;[i]:=y[i]*yscale;
vx[i]:=vx[i]*vscale;[i]:=vy[i]*vscale;
end;(datain);;;
Рис. 5.4 Начальные условия, использованные в задаче 5.1 с параметрами, приведенными в подпрограмме initiаl
Затем мы реализуем алгоритм Верле для численного решения уравнений движения Ньютона. Обратите внимание на то, что скорость частично обновляется с использованием старого ускорения. Далее вызывается подпрограмма ассеl, которая вычисляет ускорение, используя новую координату, и скорость еще раз изменяется. Подпрограмма Vеrlеt обращается к подпрограмме реriоdiс, которая в свою очередь обеспечивает рассмотрение частиц только в центральной ячейке.
procedure Verlef(var x,y,vx,vy,ax,ay:component;
N:integer;,Ly,dt,dt2:real;
var virial,xflux,yflux,pe,ke:real);
i:integer;,ynew:real;i:=1 to N do:=x[i]+vx[i]*dt+0.5*ax[i]*dt2;
ynew:=y[i]+vy[i]*dt+0.5*ay[i]*dt2;[i]:=vx[i]+0.5*ax[i]*dt;[i]:=vy[i]+0.5*ay[i]*dt;
periodic(xnew,ynew,xflux,yflux,vx[i],vy[i],Lx,Ly);[i]:=xnew;[i]:=ynew;;
accel(x,y,ax,ay,N,Lx,Ly,pe);
for i:=1 to N do
vx[i]:=vx[i]+0.5*ax[i]*dt;[i]:=vy[i]+0.5*ay[i]*dt;:=ke+0.5*(vx[i]*vx[i]+vy[i]*vy[i]);:=virial+x[i]*ax[i]+y[i]*ay[i];
end;;
В подпрограмме ассеl для нахождения полной силы, действующей и каждую частицу, используется третий закон Ньютона. (Напомним, что нашей системе единиц масса частицы равна единице.) Обращение к подпрограмме sераrаtiоn обеспечивает, что расстояние между частицами никогда не будет больше Lх/2 в х-направлении и Lу/2 в у-направлении.
procedure accel(var x,y,ax,ay:component;:integer;,Ly:real;pe:real);
var,j:integer;,dy,r,force,potential:real;i:=1 to N do[i]:=0.0;[i]:=0.0;;i:=1 to (N-1) doj:=(i+1) to N do
dx:=x[i]-x[j];:=y[i]-y[j];
separation(dx,dy,Lx,Ly);
r:=sqrt(dx*dx+dy*dy);
f(r,force,potential);[i]:=ax[i]+force*dx;[i]:=ay[i]+force*dy;[j]:=ax[j]-force*dx;[j]:=ay[j]-force*dy;:=pe+potential;;;
Потенциал и сила вычисляются в подпрограмме f.
procedure f (var force,potential:real;:real);,ri3,ri6,g:real;
r:=1;:=1.0/r;:=ri*ri*ri;:=ri3*ri3;:=24.0*ri*ri6*(2.0*ri6-1.0);
end;
Координаты частиц выдаются на экран в подпрограмме snapshot. Частота, с которой координаты частиц должны обновляться на экране, зависит от ?t.
procedure results (N,nave:integer;,Ly,dt:real;virial,ke,pe,xflux,yflux,time:real);
pflux,pvirial,E,T:real;
begintime=0.0 then
writeln(время,T,E,pflux,pviral);
time:=time+dt*nave;:=ke/nave;
pe:=pe/nave;:=(pe+ke)/N;:=ke/N;:=0.0;:=0.0;:=((xflux/(2.0*lx))+(yflux/(2.0*Ly)))/(dt*nave);:=0.0;:=0.0;:=(N*T)/(Lx*Ly)+0.5*virial/(nave*Lx*Ly);:=0.0;(time:9:3,T:9:4,E:9:4,pflux:9:4,pvirial:9:4);
end;
5.4 Измерение макроскопических величин
Макросостояние газа мы характеризовали числом частиц в левой половине ящика. Известно, что равновесное макросостояние можно характеризовать и другими параметрами, такими как абсолютная температура T, среднее давление Р, а также объемом м полной энергией.
Кинетическое определение температуры вытекает из теоремы о равнораспределении: каждый квадратичный член, входящий в выражение для энергии классической системы, находящейся в равновесии при температуре Т, имеет среднее значение , где kB -постоянная Больцмана. Отсюда мы можем определить температуру Т системы в d-мерном пространстве соотношением
(5.4)
где сумма берется по всем N частицам системы и d компонентам скорости. Скобки обозначают усреднение по времени. Выражение (5.4) представляет собой пример связи макроскопической величины, в данном случае температуры, с временным средним по траекториям частиц. Заметим, что соотношение (6.4) справедливо только в том случае, если движение центра масс системы равно нулю-не хватает еще, чтобы движение центра масс резервуара изменяло температуру! В системе СИ температура Т измеряется в градусах Кельвина (К) и постоянная kB= 1.38*10-23 Дж/К. В дальнейшем мы будем измерять температуру в единицах ?/kB.
Еще одной тепловой характеристикой является теплоемкость при постоянном объеме
СV = (?Е/?Т)V,
где Е - полная энергия. СV есть мера количества тепла, необходимого, чтобы произвести изменение температуры. Поскольку теплоемкость зависит от размеров системы, удобно определить удельную теплоемкость на частицу, а именно сV = СV/N. Легче всего получить сV путем нахождения средней потенциальной энергии при соседних температурах Т и Т + ?Т; сV складывается из температурной зависимости потенциальной энергии и удельной теплоемкости, связанной с кинетической энергией, т.е. (d/2)kB.
Чтобы определить среднее давление, предположим что частицы наxoдятся в резервуаре с жесткими стенками. Как известно, столкновения частиц со стенками резервуара приводят к тому, что на каждый элемент стенки действует средняя результирующая сила. Средняя сила, действующая на единичную площадку, и есть давление P газа. Давление можно найти, связывая силу в горизонтальном направлении со скоростью изменения соответствующей компоненты импульса частиц, ударяющихся о стенку.
Из тех же соображений можно получить среднее давление в отсутствие стенок. Поскольку давление в равновесном состоянии во всех направлениях одинаково, мы можем связать давление с общим переносом количества движения через элемент поверхности в любом месте системы. Рассмотри?/p>