МЕТОД БОЛЬШИХ ЧАСТИЦ

(МОДЕЛИРОВАНИЕ НА ПЭВМ)

Майер Р.В.

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

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

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

Здесь и далее расстояние между частицами г и сила взаимодействия F заданы в условных единицах.
При неупругих деформациях тела, когда i-тая частица смещается относительно соседней j-той частицы, на i-тую частицу действует сила вязкого трения. Она направлена в сторону противоположную относительному движению. Чтобы рассчитать ее проекции на координатные оси, на каждом временном шаге t нужно сформировать матрицу S, элементами которой являются расстояния     между i-ой и j-ой частицами в момент t-1. Тогда сила внутреннего трения будет равна:

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

Рис. 1. Падение неупругого тела на горизонтальную поверхность.

program Dvigen_tela_metod_chastic;
{$N+} uses crt, graphs-
const N=120; m=l; r=25E+4; dt=0.001; b=10;
var Fx,Fy,x,y,vx,vy:array[1..N] of single;
s:array[0..N,0..N]of single;
Gd,Gm,i,j,d,t: integer; k,ax,ay,F,l: single;
Procedure Sila; label Metka;
{PR-1}
{ Free Pascal }
begin For i:=l 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>b then goto Metka;
If (l>0)and(K=b) then F: =1000* (8-1) else F:=0;
If 1<6 then begin F:=10000; end;
Fx[i] :=Fx[i] + (P-r*(l-s[ifj]))*(x[i]-x[j])/(l+0.001) ;
Fy[i]:=Fy[i]+(F-r*(l-s[i,j]))*(y[i]-y[j])/(l+0.001)+m*2;
Metka: end; end;
Procedure Nach_uslov;
begin Randomize; d:=6;
For j:=l to 20 do For i:=l to d do begin
If j mod 2=0 then x[i+d*(j-1)]:=8*i
else x[i+d*(j-1)]:=8*i+4;
у[i+d*(j-1)]:=7*j+140; end;
For i:=1 to N do For j:=1 to N do begin
s[i,j]:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); end;
For i:=l to N do begin vx[i]:=100; vy[i]:=60; end; end;
BEGIN Gd:=Detect;
InitGraph(Gd, Gm, f с:\bp\bgi f) ; Nach_uslov;
Repeat Sila;
For i:=1 to N do For j:=1 to N do
begin s[i,j]:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); end;
For i:=l to N do begin ax:=Fx[i]/m; ay:=Fy[i]/m;
vx[i]:=vx[i]+ax*dt; vy[i]:=vy[i]+ay*dt;
x[i]:=x[i]+vx[i]*dt; y[i]:=y[i]+vy[i]*dt;
If (y[i]>400) then begin vy[i]:=0; vx[i]:=0;
y[i]:=399; end; end; line(0,402,940,402);
circle(10+round(x[l]),round(y[1]),1);
circle(10+round(x[N]),round(y[N]), 1) ;
If t mod 1000=0 then begin t:=0;
For i:=l to N do circle(10+round(x[i]),round(y[i]),3);
end; inc(t);
until KeyPressed; CloseGraph;
END.

Задача 2
Тело кубической формы находится в невесомости и взаимодействует с частицей известной массы так:
1) частица отлетает от одной из его вершин, двигаясь вверх;
2) движущаяся частица соударяется с одной из вершин куба.
Считая, что движение кубика плоское, необходимо рассчитать состояние системы в последующие моменты времени.
Задача решается аналогично. При моделировании взаимодействия частицы с первоначально покоящимся кубом получаются результаты, представленные на рис. 2.1. Если частица отлетает от одной из его вершин, двигаясь вверх, то в соответствии с законами сохранения импульса и момента импульса, центр масс куба начнет двигаться вниз, а сам куб - вращаться. Когда движущаяся частица
соударяется с одной из вершин куба, в общем случае происходит нецентральный удар; частица и куб разлетаются (рис. 2.2). На рисунках также показаны траектории движения двух максимально удаленных точек тела. Чтобы куб сохранял свою форму, необходимо увеличить коэффициенты в выражении для силы упругости и вязкого трения. При упругом соударении частицы куба с горизонтальной поверхностью изменяется знак вертикальной проекции скорости.

Рис. 2. Частица отлетает от одной из вершин куба, Частица сталкивается с кубом; куб отскакивает от поверхности.
Задача 3
Горизонтальная балка закреплена за левый конец, а к правому концу подвешен груз известной массы. Центральная часть балки имеет дефекты, уменьшающие ее прочность. Промоделируйте деформацию балки, разрушение и движение ее частей.
Все точки балки перемещаются параллельно вертикальной плоскости (движение плоское).

Рис. 3. Моделирование разрушения балки.
Балку представим в виде совокупности частиц шарообразной формы, расположенных в плоскости XOY и взаимодействующих друг с другом с силами притяжения и отталкивания. Кроме того, при относительном движении частиц возникает сила вязкого трения. Чтобы учесть дефекты в центральной части балки, можно  If (j>40)and(j<50)then kk:=300 else kk:=1000;
Правый конец балки нагружен. Это задается так:
If i>N-2*d then FF:=5 else FF:=0;
В нашем случае d = 6. Используется программа PR-2, результаты моделирования представлены на рис. 3. Видно, что под действием силы тяжести балка деформируется, разламывается на две части, одна остается на месте, а другая падает вниз и соударяется с горизонтальной поверхностью.

program Razrushenie_balki; { PR-2 }
{$N+} uses crt, graph; { Free Pascal }
const N=120; m=l; r=5E+4; dt=0.001; b=10;
var Fx,Fy,x,y,vx,vy,xx,yy: array[l..N] of single;
s:array[0..N,0..N]of single;
Gd,Gm,i,j, d,t: integer; FF,k,kk,ax,ay,F,l: single;
Procedure Sila; label Metka;
begin For i:=l 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>b then goto Metka;
If (j>40)and(j<50)then kk:=300 else kk:=1000;
If (l>0)and(K=b) then F:=kk*(8-1) else F:=0;
If 1<6 then begin F:=10000; end;
If i>N-2*d then FF:=5 else FF:=0;
Fx[i]:=Fx[i]+(P-r*(l-s[ifj]))*(x[i]-x[j])/(l+0.001);
Fy[i]:=Fy[i]+(F-r*(l-s[i,j]))*(y[i]-y[j])/(1+0.001)+m*0.3+FF;
Metka: end; end;
Procedure Nach_uslov;
begin Randomize; d:=6;
For j:=l to 20 do For i:=l to d do begin
If j mod 2=0 then у[i+d*(j-1)]:=8*i+10
else y[i+d*(j-l)]:=8*i+4+10; x[i+d*(j-1)]:=7*j+400; end;
For i:=1 to N do For j:=1 to N do
s[i,j]:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); end;
BEGIN Gd:=Detect; InitGraph(Gd,Gm,f с:\bp\bgi f) ;
Nach_uslov; Randomize;
Repeat Sila;
For i:=1 to N do For j:=1 to N do
s[i,j]:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j]));
For i:=d+l to N do begin xx[i]:=x[i]; yy[i]:=y[i];
ax:=Fx[i]/m; ay:=Fy[i]/m; vx[i]:=vx[i]+ax*dt;
vy[i]:=vy[i]+ay*dt; x[i]:=x[i]+vx[i]*dt; у[i]:=y[i]+vy[i]*dt;
If y[i]>400 then begin
y[i]:=400; vy[i]:=0; vx[i]:=0.5*vx[i]; end; end;
circle(10+round(x[N]),round(y[N]),1);
If t mod 1000=0 then begin t:=0; cleardevice;
line(0,402,940,402); line(410,0,410,80);
For i:=l to N do circle(10+round(x[i]),round(y[i]),3); end;
inc(t);
until KeyPressed; CloseGraph;
END.