Моделирование углекислого газа методом динамики частиц

Материал из Department of Theoretical and Applied Mechanics
Версия от 01:57, 11 января 2017; Foten (обсуждение | вклад) (→‎Программа)

Перейти к: навигация, поиск
Курсовые работы 2016-2017 учебного года > Моделирование углекислого газа методом динамики частиц

Курсовой проект по Механике дискретных сред

Исполнитель: Смирнов Александр

Группа: 09 (43604/1)

Семестр: осень 2016


Формулировка задачи

Смоделировать молекулы углекислого газа методом динамики частиц, проверить выполнение закона сохранения энергии и рассмотреть распределение энергии по степеням свободы.


Общие сведения

Для моделирования используем частицы, которые представляют собой абсолютно упругий шар. Масса углерода равна 12, а кислорода - 16 условных единиц. Для отталкивания молекул используется потенциал Морзе, внутри молекулы: Упругая сила и угловая пружина.

Упругая сила и угловая пружина

Упругая сила определяется формулой:

[math] F( r) = -kr, [/math]

где

  • [math]k[/math] — жесткость,
  • [math] r[/math] — отклонение от положения равновесия.

Потенциал пружины равен:

[math] \varPi(r) = С\frac{(\pi-φ)^2}{2} , [/math]

При этом

[math]φ =\frac{ {\bf r_1}·{\bf r_2}}{r_1 r_2} [/math]

где

  • [math]С[/math] —угловая жесткость,
  • [math] {\bf r_i} [/math] — радиус вектор от атома углерода к первому или второму кислороду соответственно.

А сила вычисляется как:

[math]F(r) = -\varPi'(r) , [/math]


Потенциал Морзе

Парный силовой потенциал взаимодействия. Определяется формулой:

[math] \varPi(r) = D\left[e^{-2\alpha(r-a)}-2e^{-\alpha(r-a)}\right], [/math]

где

  • [math]D[/math] — энергия связи,
  • [math]a[/math] — длина связи,
  • [math]\alpha[/math] — параметр, характеризующий ширину потенциальной ямы.

Сила, соответствующая потенциалу Морзе, вычисляется по формуле:

[math] F(r) = -\varPi'(r) = 2\alpha D\left[e^{-2\alpha(r-a)}-e^{-\alpha(r-a)}\right]. [/math]

Или в векторной форме:

[math] {\bf F}({\bf r})= -\nabla\varPi(r) = 2\alpha D\left[e^{-2\alpha(r-a)}-e^{-\alpha(r-a)}\right]\frac{{\bf r}}{r} [/math]

Программа

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

Текст программы на языке JavaScript:

window.addEventListener("load", MainBalls, true); function MainBalls() { 

   // Предварительные установки

   var canvas = canvasBalls;
   var context = canvas.getContext("2d");                  // на context происходит рисование
   canvas.oncontextmenu = function (e) {return false;};    // блокировка контекстного меню

   var Pi = 3.1415926;                   // число "пи"

   var m0 = 1;                           // масштаб массы
   var T0 = 1;                           // масштаб времени (период колебаний исходной системы)
   var a0 = 1;                           // масштаб расстояния (диаметр шара)

   var g0 = a0 / T0 / T0;                // масштаб ускорения (ускорение, при котором за T0 будет пройдено расстояние a0)
   var k0 = 2 * Pi / T0;                 // масштаб частоты
   var C0 = m0 * k0 * k0;                // масштаб жесткости
   var B0 = 2 * m0 * k0;                 // масштаб вязкости

   // *** Задание физических параметров ***

   var Ny = 5;                           // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)
   var m = 12 * m0;                       // масса C

var m1 = 16 * m0 ; //масса О 

   var Cwall = 8 * C0;                  // жесткость стен
   var B = 0 * B0;                   // вязкость среды
   var Bwall = 0 * B0;                // вязкость на стенках
   var Cball =  Cwall;              // жесткость между частицами            
   var r = 0.1 * a0;                     // радиус частицы в расчетных координатах

var K = 0.85; // сила взаимодействия ограничивается значением, реализующимся при r/a = K 

   var a = 2.2 * r;                        // равновесное расстояние между частицами

var a_1 = 2.05 * r; 

   var a_2 = 2.1 * a ;

var aCut = 2.05 * r ; // радиус обрезания var alfa = 4 ; var j = 20 ; var Ek = 0 ; var Ek1 = 0 ; var k = 1 ; var Ekc = 0 ; var Ekck = 0; var EC = 0 ; var EO1 = 0; var EO2 = 0; var ECk = 0 ; var EO1k = 0; var EO2k = 0; 

   // *** Задание вычислительных параметров ***

   var fps = 60;                         // frames per second - число кадров в секунду (качечтво отображения)
   var spf = 400;                        // steps per frame   - число шагов интегрирования между кадрами (скорость расчета)
   var dt  = 0.03 * T0 / fps;           // шаг интегрирования (качество расчета)

   // Выполнение программы

   var scale = canvas.height / Ny / a0 ;  // масштабный коэффициент для перехода от расчетных к экранным координатам
   var r2 = r * r ;                       // ___в целях оптимизации___
   var aCut2 = aCut * aCut ;              // ___в целях оптимизации___
   var a2 = a * a ;                       // ___в целях оптимизации___
   var a22 = a_2 * a_2 ;
   var a11 = a_1 * a_1 ;

var D = a2 * Cball / 72 ; // энергия связи между частицами 

   var LJCoeff = 12 * D / a2 ;            // коэффициент для расчета потенциала Л-Дж

var a1 = alfa / a ; var MorzCoeff = 2 * a1 * D  ; 

   var Ka = K * r ;                       // ___в целях оптимизации___
   var K2a2 = Ka*Ka ;                // ___в целях оптимизации___

   var w = canvas.width / scale ;           // ширина окна в расчетных координатах
   var h = canvas.height / scale ;          // высота окна в расчетных координатах


var dNd = null ; // ссылка на захваченный курсором шар (drag & drop) 

   // Работа с мышью

   var mx_, my_;                           // буфер позиции мыши (для расчета скорости при отпускании шара)

   canvas.onmousedown = function(e) {      // функция при нажатии клавиши мыши
       var m = mouseCoords(e);             // получаем расчетные координаты курсора мыши
       // цикл в обратную сторону, чтобы захватывать шар, нарисованный "сверху"
       // (т.к. цикл рисования идет в обычном порядке)
       for (var i = balls.length - 1; i >= 0; i--) {
           var b = balls[i];
           var rx = b.x - m.x;
           var ry = b.y - m.y;
           var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
           if (rLen2 <= r2) {                  // курсор нажал на шар
               if (e.which == 1) {             // нажата левая клавиша мыши
                   dNd = b;
                   dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
                   dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
                   mx_ = m.x;    my_ = m.y;
                   canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
               } else if (e.which == 3) {              // нажата правая клавиша мыши
                   balls.splice(i, 1);                 // удалить шар
               }
               return;
           }
       }

for (var i = C.length - 1; i >= 0; i--) { 

           var b = C[i];
           var rx = b.x - m.x;
           var ry = b.y - m.y;
           var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
           if (rLen2 <= r2) {                  // курсор нажал на шар
               if (e.which == 1) {             // нажата левая клавиша мыши
                   dNd = b;
                   dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
                   dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
                   mx_ = m.x;    my_ = m.y;
                   canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
               } else if (e.which == 3) {              // нажата правая клавиша мыши
                   C.splice(i, 1);                 // удалить шар
               }
               return;
           }
       }

for (var i = O1.length - 1; i >= 0; i--) { 

           var b = O1[i];
           var rx = b.x - m.x;
           var ry = b.y - m.y;
           var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
           if (rLen2 <= r2) {                  // курсор нажал на шар
               if (e.which == 1) {             // нажата левая клавиша мыши
                   dNd = b;
                   dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
                   dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
                   mx_ = m.x;    my_ = m.y;
                   canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
               } else if (e.which == 3) {              // нажата правая клавиша мыши
                   O1.splice(i, 1);                 // удалить шар
               }
               return;
           }
       }

for (var i = O2.length - 1; i >= 0; i--) { 

           var b = O2[i];
           var rx = b.x - m.x;
           var ry = b.y - m.y;
           var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
           if (rLen2 <= r2) {                  // курсор нажал на шар
               if (e.which == 1) {             // нажата левая клавиша мыши
                   dNd = b;
                   dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
                   dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
                   mx_ = m.x;    my_ = m.y;
                   canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
               } else if (e.which == 3) {              // нажата правая клавиша мыши
                   O2.splice(i, 1);                 // удалить шар
               }
               return;
           }
       }

       // если не вышли по return из цикла - нажатие было вне шара, добавляем новый
      if (e.which == 1) {
           dNd = addNewBall(m.x, m.y);         // добавляем шар и сразу захватываем его курсором
           if (dNd == null) return;            // если шар не добавился (из за стен или других шаров) - возвращаемся
           dNd.xPlus = 0;  dNd.yPlus = 0;      // держим шар по центру
           mx_ = m.x;    my_ = m.y;
           canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
       }


}; this.setNy = function(ny) { 

       j = ny;   
   };

this.ResetTime = function() { 

       timers();   
   };

this.newSystem = function() { balls = []; C =[] ; O1 = [] ; O2 = [] ; if(j == 2){ r = 0.2 * a0 ; a = 3.1 * r; // равновесное расстояние между частицами a_1 = 2.1 * r; 

    a_2 = 2.1 * a ;

aCut = 2 * r ;


r2 = r * r ; aCut2 = aCut * aCut ; a2 = a * a ; a22 = a_2 * a_2 ; a11 = a_1 * a_1 ; D = a2 * Cball / 72 ; LJCoeff = 12 * D / a2 ; a1 = alfa / a ; MorzCoeff = 2 * a1 * D  ; Ka = K * r ; K2a2 = Ka*Ka ;



addNewC( w* 0.25 , h*0.2, 10 ); addNewO1(w* 0.25 + a , h*0.2 , 10 ); addNewO2(w* 0.25 - a, h*0.2 , 10 );

addNewC( w* 0.75 , h*0.2, 10 ); addNewO1(w* 0.75 + a ,h*0.2, 10 ); addNewO2(w* 0.75 - a, h*0.2, 10 );


}

if(j == 8){ r = 0.2 * a0 ; a = 2.5 * r; // равновесное расстояние между частицами a_1 = 2.1 * r; 

    a_2 = 2.1 * a ;

aCut = 2 * r ;

j = j*0.5 ;

r2 = r * r ; aCut2 = aCut * aCut ; a2 = a * a ; a22 = a_2 * a_2 ; a11 = a_1 * a_1 ; D = a2 * Cball / 72 ; LJCoeff = 12 * D / a2 ; a1 = alfa / a ; MorzCoeff = 2 * a1 * D  ; Ka = K * r ; K2a2 = Ka*Ka ;

for (var i = 0; i < j; i++){

addNewC( w* 0.25 , h*0.2*(1+i), 10 ); addNewO1(w* 0.25 + a , h*0.2*(1+i) , 10 ); addNewO2(w* 0.25 - a, h*0.2*(1+i) , 10 );

addNewC( w* 0.75 , h*0.2*(1+i), 10 ); addNewO1(w* 0.75 + a , h*0.2*(1+i) , 10 ); addNewO2(w* 0.75 - a, h*0.2*(1+i) , 10 );

} } if(j == 20){ j = j*0.25 ; r = 0.15 * a0 ; a = 2.5 * r; a_1 = 2.1 * r; 

    a_2 = 2.1 * a ;

aCut = 2 * r ; r2 = r * r ; aCut2 = aCut * aCut ; a2 = a * a ; a22 = a_2 * a_2 ; a11 = a_1 * a_1 ; D = a2 * Cball / 72 ; LJCoeff = 12 * D / a2 ; a1 = alfa / a ; MorzCoeff = 2 * a1 * D  ; Ka = K * r ; K2a2 = Ka*Ka ;

for (var i = 0; i < j; i++){


addNewC( w* 0.2 , h*0.16*(1+i), 10 ); addNewO1(w* 0.2 + a , h*0.16*(1+i) , 10 ); addNewO2(w* 0.2 - a, h*0.16*(1+i) , 10 );

addNewC( w* 0.4 , h*0.16*(1+i), 10 ); addNewO1(w* 0.4 + a , h*0.16*(1+i) , 10 ); addNewO2(w* 0.4 - a, h*0.16*(1+i) , 10 );

addNewC( w* 0.6 , h*0.16*(1+i), 10 ); addNewO1(w* 0.6 + a , h*0.16*(1+i) , 10 ); addNewO2(w* 0.6 - a, h*0.16*(1+i) , 10 );

addNewC( w* 0.8 , h*0.16*(1+i), 10 ); addNewO1(w* 0.8 + a , h*0.16*(1+i) , 10 ); addNewO2(w* 0.8 - a, h*0.16*(1+i) , 10 );

} } if(j == 40){ r = 0.1 * a0 ; a = 3.1 * r; a_1 = 2.1 * r; a_2 = 2.1 * a ; aCut = 2 * r ;

j = j*0.2 ;

r2 = r * r ; aCut2 = aCut * aCut ; a2 = a * a ; a22 = a_2 * a_2 ; a11 = a_1 * a_1 ; D = a2 * Cball / 72 ; LJCoeff = 12 * D / a2 ; a1 = alfa / a ; MorzCoeff = 2 * a1 * D  ; Ka = K * r ; K2a2 = Ka*Ka ;

for (var i = 0; i < j; i++){ addNewC( w* 0.16 , h*0.11*(1+i), 10 ); addNewO1(w* 0.16 + a , h*0.11*(1+i) , 10 ); addNewO2(w* 0.16 - a, h*0.11*(1+i) , 10 );

addNewC( w* 0.33 , h*0.11*(1+i), 10 ); addNewO1(w* 0.33 + a , h*0.11*(1+i) , 10 ); addNewO2(w* 0.33 - a, h*0.11*(1+i) , 10 );

addNewC( w* 0.50 , h*0.11*(1+i), 10 ); addNewO1(w* 0.50 + a , h*0.11*(1+i) , 10 ); addNewO2(w* 0.50 - a, h*0.11*(1+i) , 10 );

addNewC( w* 0.67 , h*0.11*(1+i), 10 ); addNewO1(w* 0.67 + a , h*0.11*(1+i) , 10 ); addNewO2(w* 0.67 - a, h*0.11*(1+i) , 10 );

addNewC( w* 0.84 , h*0.11*(1+i), 10 ); addNewO1(w* 0.84 + a , h*0.11*(1+i) , 10 ); addNewO2(w* 0.84 - a, h*0.11*(1+i) , 10 );

} } if(j == 70){ r = 0.08 * a0 ; a = 2.3 * r; a_1 = 2.1 * r; a_2 = 2.1 * a ; aCut = 2 * r ;

j = j/7 ;

r2 = r * r ; aCut2 = aCut * aCut ; a2 = a * a ; a22 = a_2 * a_2 ; a11 = a_1 * a_1 ; D = a2 * Cball / 72 ; LJCoeff = 12 * D / a2 ; a1 = alfa / a ; MorzCoeff = 2 * a1 * D  ; Ka = K * r ; K2a2 = Ka*Ka ;

for (var i = 0; i < j; i++){

addNewC( w* 0.125 , h*0.09*(1+i), 10 ); addNewO1(w* 0.125 + a , h*0.09*(1+i) , 10 ); addNewO2(w* 0.125 - a, h*0.09*(1+i) , 10 );

addNewC( w* 0.25 , h*0.09*(1+i), 10 ); addNewO1(w* 0.25 + a , h*0.09*(1+i) , 10 ); addNewO2(w* 0.25 - a, h*0.09*(1+i) , 10 );

addNewC( w* 0.375 , h*0.09*(1+i), 10 ); addNewO1(w* 0.375 + a , h*0.09*(1+i) , 10 ); addNewO2(w* 0.375 - a, h*0.09*(1+i) , 10 );

addNewC( w* 0.5 , h*0.09*(1+i), 10 ); addNewO1(w* 0.5 + a , h*0.09*(1+i) , 10 ); addNewO2(w* 0.5 - a, h*0.09*(1+i) , 10 );

addNewC( w* 0.625 , h*0.09*(1+i), 10 ); addNewO1(w* 0.625 + a , h*0.09*(1+i) , 10 ); addNewO2(w* 0.625 - a, h*0.09*(1+i) , 10 );

addNewC( w* 0.75 , h*0.09*(1+i), 10 ); addNewO1(w* 0.75 + a , h*0.09*(1+i) , 10 ); addNewO2(w* 0.75 - a, h*0.09*(1+i) , 10 );

addNewC( w* 0.875 , h*0.09*(1+i), 10 ); addNewO1(w* 0.875 + a , h*0.09*(1+i) , 10 ); addNewO2(w* 0.875 - a, h*0.09*(1+i) , 10 ); } } for (var i = 0; i < C.length; i++) { C[i].vx = 0.4 * (1 - 2 * Math.random()); 

                   C[i].vy = 0.4 * (1 - 2 * Math.random());

O1[i].vx = 0.3 * (1 - 2 * Math.random()); 

                   O1[i].vy = 0.3 * (1 - 2 * Math.random());

O2[i].vx = 0.3 * (1 - 2 * Math.random()); 

                   O2[i].vy = 0.3 * (1 - 2 * Math.random());

} timers(); 

   }

   document.onmouseup = function(e) {          // функция при отпускании клавиши мыши
       canvas.onmousemove = null;              // когда клавиша отпущена - функции перемещения нету
       dNd = null;                             // когда клавиша отпущена - захваченного курсором шара нету
   };

   function mouseMove(e) {                     // функция при перемещении мыши, работает только с зажатой ЛКМ
       var m = mouseCoords(e);                 // получаем расчетные координаты курсора мыши
       dNd.x = m.x + dNd.xPlus;
       dNd.y = m.y + dNd.yPlus;
       dNd.vx = 0.6 * (m.x - mx_) / dt / fps;   dNd.vy = 0.6 * (m.y - my_) / dt / fps;
       mx_ = m.x;    my_ = m.y;
   }

   function mouseCoords(e) {                   // функция возвращает расчетные координаты курсора мыши
       var m = [];
       var rect = canvas.getBoundingClientRect();
       m.x = (e.clientX - rect.left) / scale;
       m.y = (e.clientY - rect.top) / scale;
       return m;
   }


var addNewBall = function(x, y) { 

       // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
       if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
       for (var i = 0; i < balls.length; i++) {
           var rx = balls[i].x - x;
           var ry = balls[i].y - y;
           var rLen2 = rx * rx + ry * ry;
           if (rLen2 < 4 * r2) return null;
       }
       var b = [];

       b.x = x;                b.y = y;        // расчетные координаты шара
       b.fx = 0;               b.fy = 0;      // сила, действующая на шар
       b.vx = 0;               b.vy = 0;   		// скорость 

       balls[balls.length] = b;                // добавить элемент в конец массива

return b; 

   };

   var addNewC =  function(x, y) {
       // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
       if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
       for (var i = 0; i < balls.length; i++) {
           var rx = balls[i].x - x;
           var ry = balls[i].y - y;
           var rLen2 = rx * rx + ry * ry;
           if (rLen2 < 4 * r2) return null;
       }
       var b = [];

       b.x = x;                b.y = y;        // расчетные координаты шара
       b.fx = 0;               b.fy = 0;      // сила, действующая на шар
       b.vx = 0;               b.vy = 0;   		// скорость

       C[C.length] = b;                // добавить элемент в конец массива

return b; 

   };

   var addNewO1 =  function(x, y) {
       // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
       if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
       for (var i = 0; i < balls.length; i++) {
           var rx = balls[i].x - x;
           var ry = balls[i].y - y;
           var rLen2 = rx * rx + ry * ry;
           if (rLen2 < 4 * r2) return null;
       }
       var b = [];

       b.x = x;                b.y = y;        // расчетные координаты шара
       b.fx = 0;               b.fy = 0;      // сила, действующая на шар
       b.vx = 0;               b.vy = 0;   		// скорость

       O1[O1.length] = b;                // добавить элемент в конец массива

return b; 

   };

var addNewO2 = function(x, y) { 

       // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
       if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
       for (var i = 0; i < balls.length; i++) {
           var rx = balls[i].x - x;
           var ry = balls[i].y - y;
           var rLen2 = rx * rx + ry * ry;
           if (rLen2 < 4 * r2) return null;
       }
       var b = [];

       b.x = x;                b.y = y;        // расчетные координаты шара
       b.fx = 0;               b.fy = 0;      // сила, действующая на шар
       b.vx = 0;               b.vy = 0;   		// скорость

       O2[O2.length] = b;                // добавить элемент в конец массива

return b; 

   };

// Основной цикл программы 

   function control() {
       physics();
       draw();
   }

function timers(){ clearInterval(id1); clearInterval(id2); clearInterval(id3); clearInterval(id4); clearInterval(id5); k = 1; Ekc = 0 ; Ekck = 0; EC= 0; EO1 = 0 ; EO2 = 0; ECk = 0; EO1k =0 ; EO2k = 0 ; Ek = 0; Ek1 = 0 ;

var id1 = setInterval(function(){document.getElementById('Ekck').innerHTML = Ekck;}, 1000 / 2); var id2 =setInterval(function(){document.getElementById('Ek').innerHTML = Ek;}, 1000 / 2); var id3 =setInterval(function(){document.getElementById('Ekc').innerHTML = ECk;}, 1000 / 2); var id4 =setInterval(function(){document.getElementById('Eko1').innerHTML = EO1k;}, 1000 / 2); var id5 =setInterval(function(){document.getElementById('Eko2').innerHTML = EO2k;}, 1000 / 2);

} 

   // Расчетная часть программы

function powers(b ,b2 , hkk , jk) { var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) 

       var r2 = rx * rx + ry * ry;                         // квадрат расстояния между шарами

if ( r2 > aCut2 ) return null;

var rLen = (Math.sqrt(r2));

if (r2 < K2a2) { 

                       if (rLen > 0.00001) {                           // проверка, чтобы избежать деления на 0
                           rx = rx / rLen * Ka;
                           ry = ry / rLen * Ka;
                       }
                       r2 = K2a2;
                       rLen = Ka;                                      // корень K2a2
                   }

                   // сила взаимодействия

var u = Math.exp( -a1 * ( rLen - hkk )) ; var F = jk * MorzCoeff * u* ( u - 1 ) /rLen ; 

                   var Fx = F * rx;        var Fy = F * ry;
                   b.fx += Fx;             b.fy += Fy;
                   b2.fx -= Fx;            b2.fy -= Fy;

} 


   function physics() {                        // то, что происходит каждый шаг времени
       for (var s = 1; s <= spf; s++) {

           // пересчет сил идет отдельным массивом, т.к. далее будут добавляться силы взаимодействия между шарами

for (var i0 = 0; i0 < balls.length; i0++) { 

               balls[i0].fx = - B * balls[i0].vx;
               balls[i0].fy =  - B * balls[i0].vy;
           }

for (var i0 = 0; i0 < C.length; i0++) { 

               C[i0].fx = - B * C[i0].vx;
               C[i0].fy =  - B * C[i0].vy;
           }

for (var i0 = 0; i0 < O1.length; i0++) { 

               O1[i0].fx = - B * O1[i0].vx;
               O1[i0].fy =  - B * O1[i0].vy;
           }

for (var i0 = 0; i0 < O2.length; i0++) { 

               O2[i0].fx = - B * O2[i0].vx;
               O2[i0].fy =  - B * O2[i0].vy;
           }



           for (var i = 0; i < balls.length; i++) { // dlya Balls

var b = balls[i];

for (var j = i + 1; j < balls.length; j++) { 

                   var b2 = balls[j];

powers(b,b2,a_1,0.3) ; }

for (var j = 0; j < C.length; j++) { 

                   var b2 = C[j];

powers(b,b2,a_1,0.3) ;

} for (var j = 0; j < O1.length; j++) { 

                   var b2 = O1[j];

powers(b,b2,a_1,0.3) ; }

for (var j = 0; j < O2.length; j++) { 

                   var b2 = O2[j];

powers( b, b2, a_1 ,0.3) ; }

if (b == dNd) continue; // если шар схвачен курсором - его вз. со стенами и перемещение не считаем 

               if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
               if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
               if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
               if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }

               b.vx += b.fx / m * dt;        b.vy += b.fy / m * dt;
               b.x += b.vx * dt;             b.y += b.vy * dt;

Ek1 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy); 

           }

for (var i = 0; i < C.length; i++) { for (var j = i ; j < C.length; j++) {

if ( i == j){ var r1x = O1[i].x - C[i].x ; var r1y = O1[i].y - C[i].y ; var r2x = O2[i].x - C[i].x ; var r2y = O2[i].y - C[i].y ;

var r1 = Math.sqrt(r1x*r1x +r1y* r1y); var r2 = Math.sqrt(r2x*r2x +r2y* r2y); var xx =(r1x*r2x + r1y*r2y)/(r1*r2) ;

var F1 = -0.2*Cball*(r1-a)/r1; var F2 = -0.2*Cball*(r2-a)/r2; var f = 0 ; if(xx*xx > 0.99999999999){ f = -0.3*Cball; }

else{ f = -0.3*Cball*(Pi - Math.acos(xx))/(Math.sqrt(1-xx*xx)); } ff = f ;

var f1x = f*(r2x/(r1*r2) - xx*r1x/(r1*r1)); var f1y = f*(r2y/(r1*r2) -xx*r1y/(r1*r1)); var f2x = f*(r1x/(r1*r2) - xx*r2x/(r2*r2)); var f2y = f*(r1y/(r1*r2) -xx*r2y/(r2*r2));

var F1x = F1 * r1x  ; var F1y = F1 * r1y ; var F2x = F2 * r2x  ; var F2y = F2 * r2y ;

O1[i].fx += (F1x + f1x) ; O1[i].fy += (F1y + f1y) ; O2[i].fx += (F2x + f2x) ; O2[i].fy += (F2y + f2y) ; C[i].fx -= (F1x + F2x + f1x + f2x) ; C[i].fy -= (F1y + F2y + f1y + f2y) ;

powers(C[i],C[j],a_1, 0.3) ; powers(C[i],O1[j],a_1, 0.3) ; powers(C[i],O2[j],a_1, 0.3) ; powers(O1[i],C[j],a_1, 0.3) ; powers(O1[i],O1[j],a_1, 0.3) ; powers(O1[i],O2[j],a_1, 0.3) ; powers(O2[i],C[j],a_1, 0.3) ; powers(O2[i],O1[j],a_1, 0.3) ; powers(O2[i],O2[j],a_1, 0.3) ;

}

powers(C[i],C[j],a_1, 0.3) ; powers(C[i],O1[j],a_1, 0.3) ; powers(C[i],O2[j],a_1, 0.3) ; powers(O1[i],C[j],a_1, 0.3) ; powers(O1[i],O1[j],a_1, 0.3) ; powers(O1[i],O2[j],a_1, 0.3) ; powers(O2[i],C[j],a_1, 0.3) ; powers(O2[i],O1[j],a_1, 0.3) ; powers(O2[i],O2[j],a_1, 0.3) ;

} }


for (var i = 0; i < C.length; i++) { // dlya C

var b = C[i];

if (b == dNd) continue; // если шар схвачен курсором - его вз. со стенами и перемещение не считаем 

               if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
               if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
               if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
               if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }



               b.vx += b.fx / m * dt;        b.vy += b.fy / m * dt;
               b.x += b.vx * dt;             b.y += b.vy * dt;

Ek1 += m* 0.5 *(b.vx*b.vx + b.vy*b.vy) ; EC += m* 0.5 *(b.vx*b.vx + b.vy*b.vy)/j ; 

           }

for (var i = 0; i < C.length; i++) { // dlya O1

var b = O1[i];

if (b == dNd) continue; // если шар схвачен курсором - его вз. со стенами и перемещение не считаем 

               if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
               if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
               if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
               if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }

               b.vx += b.fx / m1 * dt;        b.vy += b.fy / m1 * dt;
               b.x += b.vx * dt;             b.y += b.vy * dt;

Ek1 += m1 * 0.5 *(b.vx*b.vx + b.vy * b.vy) ; EO1 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy)/j ; 


           }

for (var i = 0; i < C.length; i++) { // dlya O2

var b = O2[i];


if (b == dNd) continue; // если шар схвачен курсором - его вз. со стенами и перемещение не считаем 

               if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
               if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
               if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
               if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }

               b.vx += b.fx / m1 * dt;        b.vy += b.fy / m1 * dt;
               b.x += b.vx * dt;             b.y += b.vy * dt;

Ek1 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy) ; EO2 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy)/j ;


} 

       }

Ekc +=Ek1 ;

if(Ek1 > 1){ k++; } Ekck = Ekc/k; ECk = EC/k ; EO1k = EO1/k ; EO2k = EO2/k ; Ek = Ek1; Ek1 = 0 ; 

   }




  // Рисование

   var rScale13 = r * scale * 1.3;         // ___в целях оптимизации___
   var rScaleShift = r * scale / 5;        // ___в целях оптимизации___
   function draw() {
       context.clearRect(0, 0, w * scale, h * scale);      // очистить экран
       for (var i = 0; i < balls.length; i++){
           var xS = balls[i].x * scale;           var yS = balls[i].y * scale;
           
           context.fillStyle = "#FFD700";

           context.beginPath();
           context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
           context.closePath();
           context.fill();
       }

for (var i = 0; i < C.length; i++){ 

           var xS = C[i].x * scale;           var yS = C[i].y * scale;
           
           context.fillStyle = "#f08080";

           context.beginPath();
           context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
           context.closePath();
           context.fill();
       }

for (var i = 0; i < O1.length; i++){ 

           var xS = O1[i].x * scale;           var yS = O1[i].y * scale;
           
           context.fillStyle = "#3070d0";

           context.beginPath();
           context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
           context.closePath();
           context.fill();
       }

for (var i = 0; i < O2.length; i++){ 

           var xS = O2[i].x * scale;           var yS = O2[i].y * scale;
           
           context.fillStyle = "#3070d0";

           context.beginPath();
           context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
           context.closePath();
           context.fill();
       }

} 

   // Запуск системы

this.newSystem(); var id1 = setInterval(function(){document.getElementById('Ekck').innerHTML = Ekck;}, 1000 / 2); var id2 =setInterval(function(){document.getElementById('Ek').innerHTML = Ek;}, 1000 / 2); var id3 =setInterval(function(){document.getElementById('Ekc').innerHTML = ECk;}, 1000 / 2); var id4 =setInterval(function(){document.getElementById('Eko1').innerHTML = EO1k;}, 1000 / 2); var id5 =setInterval(function(){document.getElementById('Eko2').innerHTML = EO2k;}, 1000 / 2); 

   setInterval(control, 1000 / fps);

}

Результаты

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

[math]U= \frac{kT}{2} [/math]
а нашу систему можно приблизительно считать идеальным газом в силу малости значения потенциальной энергии и абсолютно упругих соударений.
Источник — «http://tm.spbstu.ru/?title=Моделирование_углекислого_газа_методом_динамики_частиц&oldid=45771»