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

Материал из Department of Theoretical and Applied Mechanics
Перейти к: навигация, поиск
(→‎Упругая сила и угловая пружина)
(→‎Программа)
 
(не показано 11 промежуточных версий этого же участника)
Строка 11: Строка 11:
  
 
== Формулировка задачи ==
 
== Формулировка задачи ==
 +
[[Файл:foto11.jpg|thumb|Рис.1 Молекула углекислого газа.|400px]]
 
Смоделировать молекулы углекислого газа методом динамики частиц, проверить выполнение закона сохранения энергии и рассмотреть распределение энергии по степеням свободы.
 
Смоделировать молекулы углекислого газа методом динамики частиц, проверить выполнение закона сохранения энергии и рассмотреть распределение энергии по степеням свободы.
−
 
  
 
== Общие сведения ==
 
== Общие сведения ==
 
Для моделирования используем частицы, которые представляют собой абсолютно упругий шар. Масса углерода равна 12, а кислорода - 16 условных единиц. Для отталкивания молекул используется потенциал Морзе, внутри молекулы: Упругая сила и угловая пружина.
 
Для моделирования используем частицы, которые представляют собой абсолютно упругий шар. Масса углерода равна 12, а кислорода - 16 условных единиц. Для отталкивания молекул используется потенциал Морзе, внутри молекулы: Упругая сила и угловая пружина.
−
== Упругая сила и угловая пружина ==
+
===== Упругая сила и угловая пружина =====
 
Упругая сила определяется формулой:
 
Упругая сила определяется формулой:
 
::<math>
 
::<math>
Строка 61: Строка 61:
 
     {\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}
 
     {\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>
 
</math>
 +
==Программа==
 +
В данной программе в начальный момент времени системе задаются случайные скорости(начальная энергия,они достаточно велики, чтобы можно было пренебречь потенциальной энергией взаимодействия. Можно менять количество молекул углекислого газа и сбрасывать таймер расчета средних значений. Так же выводятся:кинетическая энергия системы в данный момент времени, средняя кинетическая энергия системы в данный момент времени, средние энергии, приходящиеся на атом углерода, первый и второй атом кислорода в молекуле.
 +
*Красные шары - атомы углерода
 +
*Синие шары - атомы кислорода
 +
{{#widget:Iframe |url=https://ailurus.ru/protos/foten/lab-k/ |width=680 |height=680 |border=0 }}
 +
 +
<div class="mw-collapsible mw-collapsed" style="width:100%" >
 +
'''Текст программы на языке JavaScript:'''
 +
<div class="mw-collapsible-content">
 +
Файл '''"Lab.js"'''
 +
<syntaxhighlight lang="javascript" line start="1" enclose="div">
 +
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);
 +
}
 +
</syntaxhighlight>
 +
Файл '''"Lab.html"'''
 +
<syntaxhighlight lang="html5" line start="1" enclose="div">
 +
<!DOCTYPE html>
 +
<html>
 +
<head>
 +
    <title>Balls</title>
 +
    <script src="Balls_v4_release.js"></script>
 +
</head>
 +
<body>
 +
    <canvas id="canvasBalls" width="640" height="480" style="border:1px solid #000000;"></canvas>
 +
<br>
 +
<div>Количество частиц:
 +
        <input type="button" style="width: 30px" name="" onclick="setNy(8); newSystem();return false;" value="8"/>
 +
        <input type="button" style="width: 30px" name="" onclick="setNy(20); newSystem();return false;" value="20"/>
 +
        <input type="button" style="width: 30px" name="" onclick="setNy(40); newSystem();return false;" value="40"/>
 +
        <input type="button" style="width: 30px" name="" onclick="setNy(70); newSystem();return false;" value="70"/>
 +
    </div>
 +
<div>Обновить таймер:
 +
        <input type="button" style="width: 60px" name="" onclick="ResetTime();return false;" value="Reset"/>
 +
    </div>
 +
<div>Кин энергия в данный момент времени: <span id="Ek"></span></div>
 +
<div>Средняя кинетическая энергия: <span id="Ekck"></span></div>
 +
<div>Средняя кинетическая энергия на один атом углерода: <span id="Ekc"></span></div>
 +
<div>Средняя кинетическая энергия на один атом кислорода(1): <span id="Eko1"></span></div>
 +
<div>Средняя кинетическая энергия на один атом кислорода(2): <span id="Eko2"></span></div>
 +
</body>
 +
</html>
 +
</syntaxhighlight>
 +
</div>
 +
</div>
 +
 +
==Результаты==
 +
Получена программа, которая моделирует молекулы углекислого газа в объеме и рассчитывает среднюю кинетическую энергию отдельных частей системы и всей системы в целом.
 +
Значение кинетической энергии системы совершает колебания вокруг среднего значения энергии, которое оказалось постоянным, что говорит о выполнении закона сохранения энергии. Средние энергии, приходящиеся на каждый отдельный атом, спустя какое-то время оказываются равными. Это говорит о том, что энергии, проходящиеся на каждую степень свободы, равны. Для идеального газа эта энергия равна:
 +
::<math>U= \frac{kT}{2} </math>
 +
где
 +
* <math>k</math> — постоянная Больцмана,
 +
* <math>T</math> — температура среды, 
 +
А нашу систему можно приблизительно считать идеальным газом в силу малости значения потенциальной энергии и абсолютно упругих соударений. Значит наши результаты соответствуют, теории.

Текущая версия на 03:07, 11 января 2017

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

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

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

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

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


Формулировка задачи[править]

Рис.1 Молекула углекислого газа.

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

Общие сведения[править]

Для моделирования используем частицы, которые представляют собой абсолютно упругий шар. Масса углерода равна 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:

Файл "Lab.js" 

  1 window.addEventListener("load", MainBalls, true);
  2 function MainBalls() {
  3 
  4     // Предварительные установки
  5 
  6     var canvas = canvasBalls;
  7     var context = canvas.getContext("2d");                  // на context происходит рисование
  8     canvas.oncontextmenu = function (e) {return false;};    // блокировка контекстного меню
  9 
 10     var Pi = 3.1415926;                   // число "пи"
 11 
 12     var m0 = 1;                           // масштаб массы
 13     var T0 = 1;                           // масштаб времени (период колебаний исходной системы)
 14     var a0 = 1;                           // масштаб расстояния (диаметр шара)
 15 
 16     var g0 = a0 / T0 / T0;                // масштаб ускорения (ускорение, при котором за T0 будет пройдено расстояние a0)
 17     var k0 = 2 * Pi / T0;                 // масштаб частоты
 18     var C0 = m0 * k0 * k0;                // масштаб жесткости
 19     var B0 = 2 * m0 * k0;                 // масштаб вязкости
 20 
 21     // *** Задание физических параметров ***
 22 
 23     var Ny = 5;                           // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)
 24     var m = 12 * m0;                       // масса C
 25 	var m1 = 16 * m0 ;                   //масса О 
 26     var Cwall = 8 * C0;                  // жесткость стен
 27     var B = 0 * B0;                   // вязкость среды
 28     var Bwall = 0 * B0;                // вязкость на стенках
 29     var Cball =  Cwall;              // жесткость между частицами            
 30     var r = 0.1 * a0;                     // радиус частицы в расчетных координатах
 31 	var K = 0.85;                         // сила взаимодействия ограничивается значением, реализующимся при r/a = K
 32     var a = 2.2 * r;                        // равновесное расстояние между частицами
 33 	var a_1 = 2.05 * r;                       
 34     var a_2 = 2.1 * a ;
 35 	var aCut =  2.05 * r ;                     // радиус обрезания
 36 	var alfa = 4 ;
 37 	var j = 20 ;
 38 	var Ek = 0 ;
 39 	var Ek1 = 0 ;
 40 	var k = 1 ;
 41 	var Ekc = 0 ;
 42 	var Ekck = 0;
 43 	var EC = 0 ;
 44 	var EO1 = 0;
 45 	var EO2 = 0;
 46 	var ECk = 0 ;
 47 	var EO1k = 0;
 48 	var EO2k = 0;
 49 
 50     // *** Задание вычислительных параметров ***
 51 
 52     var fps = 60;                         // frames per second - число кадров в секунду (качечтво отображения)
 53     var spf = 400;                        // steps per frame   - число шагов интегрирования между кадрами (скорость расчета)
 54     var dt  = 0.03 * T0 / fps;           // шаг интегрирования (качество расчета)
 55 
 56     // Выполнение программы
 57 
 58     var scale = canvas.height / Ny / a0 ;  // масштабный коэффициент для перехода от расчетных к экранным координатам
 59     var r2 = r * r ;                       // ___в целях оптимизации___
 60     var aCut2 = aCut * aCut ;              // ___в целях оптимизации___
 61     var a2 = a * a ;                       // ___в целях оптимизации___
 62     var a22 = a_2 * a_2 ;
 63     var a11 = a_1 * a_1 ;  	
 64 	var D = a2 * Cball / 72 ;              // энергия связи между частицами
 65     var LJCoeff = 12 * D / a2 ;            // коэффициент для расчета потенциала Л-Дж
 66 	var a1 = alfa / a ;
 67 	var MorzCoeff = 2 * a1 * D  ; 
 68  
 69     var Ka = K * r ;                       // ___в целях оптимизации___
 70     var K2a2 = Ka*Ka ;                // ___в целях оптимизации___
 71 
 72     var w = canvas.width / scale ;           // ширина окна в расчетных координатах
 73     var h = canvas.height / scale ;          // высота окна в расчетных координатах
 74 	
 75 	
 76 	 var dNd = null ;                         // ссылка на захваченный курсором шар (drag & drop)
 77 
 78     // Работа с мышью
 79 
 80     var mx_, my_;                           // буфер позиции мыши (для расчета скорости при отпускании шара)
 81 
 82     canvas.onmousedown = function(e) {      // функция при нажатии клавиши мыши
 83         var m = mouseCoords(e);             // получаем расчетные координаты курсора мыши
 84         // цикл в обратную сторону, чтобы захватывать шар, нарисованный "сверху"
 85         // (т.к. цикл рисования идет в обычном порядке)
 86         for (var i = balls.length - 1; i >= 0; i--) {
 87             var b = balls[i];
 88             var rx = b.x - m.x;
 89             var ry = b.y - m.y;
 90             var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
 91             if (rLen2 <= r2) {                  // курсор нажал на шар
 92                 if (e.which == 1) {             // нажата левая клавиша мыши
 93                     dNd = b;
 94                     dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
 95                     dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
 96                     mx_ = m.x;    my_ = m.y;
 97                     canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
 98                 } else if (e.which == 3) {              // нажата правая клавиша мыши
 99                     balls.splice(i, 1);                 // удалить шар
100                 }
101                 return;
102             }
103         }
104 		for (var i = C.length - 1; i >= 0; i--) {
105             var b = C[i];
106             var rx = b.x - m.x;
107             var ry = b.y - m.y;
108             var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
109             if (rLen2 <= r2) {                  // курсор нажал на шар
110                 if (e.which == 1) {             // нажата левая клавиша мыши
111                     dNd = b;
112                     dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
113                     dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
114                     mx_ = m.x;    my_ = m.y;
115                     canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
116                 } else if (e.which == 3) {              // нажата правая клавиша мыши
117                     C.splice(i, 1);                 // удалить шар
118                 }
119                 return;
120             }
121         }
122 		for (var i = O1.length - 1; i >= 0; i--) {
123             var b = O1[i];
124             var rx = b.x - m.x;
125             var ry = b.y - m.y;
126             var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
127             if (rLen2 <= r2) {                  // курсор нажал на шар
128                 if (e.which == 1) {             // нажата левая клавиша мыши
129                     dNd = b;
130                     dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
131                     dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
132                     mx_ = m.x;    my_ = m.y;
133                     canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
134                 } else if (e.which == 3) {              // нажата правая клавиша мыши
135                     O1.splice(i, 1);                 // удалить шар
136                 }
137                 return;
138             }
139         }
140 		for (var i = O2.length - 1; i >= 0; i--) {
141             var b = O2[i];
142             var rx = b.x - m.x;
143             var ry = b.y - m.y;
144             var rLen2 = rx * rx + ry * ry;      // квадрат расстояния между курсором и центром шара
145             if (rLen2 <= r2) {                  // курсор нажал на шар
146                 if (e.which == 1) {             // нажата левая клавиша мыши
147                     dNd = b;
148                     dNd.xPlus = dNd.x - m.x;    // сдвиг курсора относительно центра шара по x
149                     dNd.yPlus = dNd.y - m.y;    // сдвиг курсора относительно центра шара по y
150                     mx_ = m.x;    my_ = m.y;
151                     canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
152                 } else if (e.which == 3) {              // нажата правая клавиша мыши
153                     O2.splice(i, 1);                 // удалить шар
154                 }
155                 return;
156             }
157         }
158 
159         // если не вышли по return из цикла - нажатие было вне шара, добавляем новый
160        if (e.which == 1) {
161             dNd = addNewBall(m.x, m.y);         // добавляем шар и сразу захватываем его курсором
162             if (dNd == null) return;            // если шар не добавился (из за стен или других шаров) - возвращаемся
163             dNd.xPlus = 0;  dNd.yPlus = 0;      // держим шар по центру
164             mx_ = m.x;    my_ = m.y;
165             canvas.onmousemove = mouseMove;     // пока клавиша нажата - работает функция перемещения
166         }
167 
168 	
169 	};
170 	this.setNy = function(ny) {
171         j = ny;   
172     };
173 	this.ResetTime = function() {
174         timers();   
175     };
176 	
177 	this.newSystem = function() {
178 		 balls = [];  
179 		C =[] ; 
180 		O1 = [] ;
181 		O2 = [] ;
182 		if(j == 2){
183 		r = 0.2 * a0 ;
184 		 a = 3.1 * r;                        // равновесное расстояние между частицами
185 	 a_1 = 2.1 * r;                       
186      a_2 = 2.1 * a ;
187 	 aCut =  2 * r ; 
188 	 
189 	
190 		
191 		r2 = r * r ;                      
192 		aCut2 = aCut * aCut ;              
193 		a2 = a * a ;                       
194 		a22 = a_2 * a_2 ;
195 		a11 = a_1 * a_1 ;  	
196 		D = a2 * Cball / 72 ;              
197 		LJCoeff = 12 * D / a2 ;            
198 		a1 = alfa / a ;
199 		MorzCoeff = 2 * a1 * D  ; 
200 		Ka = K * r ;                       
201 		K2a2 = Ka*Ka ;        
202 		
203 		          
204        
205 	   
206 				addNewC( w* 0.25 , h*0.2, 10 );
207 				addNewO1(w* 0.25 + a , h*0.2 , 10 );
208 				addNewO2(w* 0.25 - a, h*0.2 , 10 );  
209 		  
210 				addNewC( w* 0.75 , h*0.2, 10 );
211 				addNewO1(w* 0.75 + a ,h*0.2, 10 );
212 				addNewO2(w* 0.75 - a, h*0.2, 10 );  
213 		   
214 		
215 	}
216 		
217 		if(j == 8){
218 		r = 0.2 * a0 ;
219 		 a = 2.5 * r;                        // равновесное расстояние между частицами
220 	 a_1 = 2.1 * r;                       
221      a_2 = 2.1 * a ;
222 	 aCut =  2 * r ; 
223 	 
224 		j = j*0.5 ;
225 		
226 		r2 = r * r ;                      
227 		aCut2 = aCut * aCut ;              
228 		a2 = a * a ;                       
229 		a22 = a_2 * a_2 ;
230 		a11 = a_1 * a_1 ;  	
231 		D = a2 * Cball / 72 ;              
232 		LJCoeff = 12 * D / a2 ;            
233 		a1 = alfa / a ;
234 		MorzCoeff = 2 * a1 * D  ; 
235 		Ka = K * r ;                       
236 		K2a2 = Ka*Ka ;        
237 		
238 		for (var i = 0; i < j; i++){               
239 	   
240 				addNewC( w* 0.25 , h*0.2*(1+i), 10 );
241 				addNewO1(w* 0.25 + a , h*0.2*(1+i) , 10 );
242 				addNewO2(w* 0.25 - a, h*0.2*(1+i) , 10 );  
243 		  
244 				addNewC( w* 0.75 , h*0.2*(1+i), 10 );
245 				addNewO1(w* 0.75 + a , h*0.2*(1+i) , 10 );
246 				addNewO2(w* 0.75 - a, h*0.2*(1+i) , 10 );  
247 		   
248 		}
249 	}
250 		if(j == 20){
251 		j = j*0.25 ;
252 		r = 0.15 * a0 ;
253 		a = 2.5 * r;                        
254 	 a_1 = 2.1 * r;                       
255      a_2 = 2.1 * a ;
256 	 aCut =  2 * r ; 
257 		r2 = r * r ;                      
258 		aCut2 = aCut * aCut ;              
259 		a2 = a * a ;                       
260 		a22 = a_2 * a_2 ;
261 		a11 = a_1 * a_1 ;  	
262 		D = a2 * Cball / 72 ;              
263 		LJCoeff = 12 * D / a2 ;            
264 		a1 = alfa / a ;
265 		MorzCoeff = 2 * a1 * D  ; 
266 		Ka = K * r ;                       
267 		K2a2 = Ka*Ka ;        	 
268 		
269 		for (var i = 0; i < j; i++){               
270        
271 	   
272 				addNewC( w* 0.2 , h*0.16*(1+i), 10 );
273 				addNewO1(w* 0.2 + a , h*0.16*(1+i) , 10 );
274 				addNewO2(w* 0.2 - a, h*0.16*(1+i) , 10 );  
275 		  
276 				addNewC( w* 0.4 , h*0.16*(1+i), 10 );
277 				addNewO1(w* 0.4 + a , h*0.16*(1+i) , 10 );
278 				addNewO2(w* 0.4 - a, h*0.16*(1+i) , 10 );  
279 				
280 				addNewC( w* 0.6 , h*0.16*(1+i), 10 );
281 				addNewO1(w* 0.6 + a , h*0.16*(1+i) , 10 );
282 				addNewO2(w* 0.6 - a, h*0.16*(1+i) , 10 ); 
283 				
284 				addNewC( w* 0.8 , h*0.16*(1+i), 10 );
285 				addNewO1(w* 0.8 + a , h*0.16*(1+i) , 10 );
286 				addNewO2(w* 0.8 - a, h*0.16*(1+i) , 10 ); 
287 		   
288 		}
289 	}
290 		if(j == 40){
291 		r = 0.1 * a0 ;
292 		 a = 3.1 * r;                       
293 		a_1 = 2.1 * r;                       
294 		a_2 = 2.1 * a ;
295 		aCut =  2 * r ;  
296 		
297 		j = j*0.2 ;
298 		
299 		r2 = r * r ;                      
300 		aCut2 = aCut * aCut ;              
301 		a2 = a * a ;                       
302 		a22 = a_2 * a_2 ;
303 		a11 = a_1 * a_1 ;  	
304 		D = a2 * Cball / 72 ;              
305 		LJCoeff = 12 * D / a2 ;            
306 		a1 = alfa / a ;
307 		MorzCoeff = 2 * a1 * D  ; 
308 		Ka = K * r ;                       
309 		K2a2 = Ka*Ka ;        
310 		
311 		for (var i = 0; i < j; i++){               
312 				addNewC( w* 0.16 , h*0.11*(1+i), 10 );
313 				addNewO1(w* 0.16 + a , h*0.11*(1+i) , 10 );
314 				addNewO2(w* 0.16 - a, h*0.11*(1+i) , 10 );  
315 		  
316 				addNewC( w* 0.33 , h*0.11*(1+i), 10 );
317 				addNewO1(w* 0.33 + a , h*0.11*(1+i) , 10 );
318 				addNewO2(w* 0.33 - a, h*0.11*(1+i) , 10 );  
319 				
320 				addNewC( w* 0.50 , h*0.11*(1+i), 10 );
321 				addNewO1(w* 0.50 + a , h*0.11*(1+i) , 10 );
322 				addNewO2(w* 0.50 - a, h*0.11*(1+i) , 10 ); 
323 				
324 				addNewC( w* 0.67 , h*0.11*(1+i), 10 );
325 				addNewO1(w* 0.67 + a , h*0.11*(1+i) , 10 );
326 				addNewO2(w* 0.67 - a, h*0.11*(1+i) , 10 ); 
327 				
328 				addNewC( w* 0.84 , h*0.11*(1+i), 10 );
329 				addNewO1(w* 0.84 + a , h*0.11*(1+i) , 10 );
330 				addNewO2(w* 0.84 - a, h*0.11*(1+i) , 10 ); 
331 		   
332 		}
333 	}
334 		if(j == 70){
335 		r = 0.08 * a0 ;
336 		 a = 2.3 * r;                       
337 		a_1 = 2.1 * r;                       
338 		a_2 = 2.1 * a ;
339 		aCut =  2 * r ;  
340 		
341 		j = j/7 ;
342 		
343 		r2 = r * r ;                      
344 		aCut2 = aCut * aCut ;              
345 		a2 = a * a ;                       
346 		a22 = a_2 * a_2 ;
347 		a11 = a_1 * a_1 ;  	
348 		D = a2 * Cball / 72 ;              
349 		LJCoeff = 12 * D / a2 ;            
350 		a1 = alfa / a ;
351 		MorzCoeff = 2 * a1 * D  ; 
352 		Ka = K * r ;                       
353 		K2a2 = Ka*Ka ;        
354 		
355 		for (var i = 0; i < j; i++){               
356        
357 				addNewC( w* 0.125 , h*0.09*(1+i), 10 );
358 				addNewO1(w* 0.125 + a , h*0.09*(1+i) , 10 );
359 				addNewO2(w* 0.125 - a, h*0.09*(1+i) , 10 );  
360 		  
361 				addNewC( w* 0.25 , h*0.09*(1+i), 10 );
362 				addNewO1(w* 0.25 + a , h*0.09*(1+i) , 10 );
363 				addNewO2(w* 0.25 - a, h*0.09*(1+i) , 10 );  
364 				
365 				addNewC( w* 0.375 , h*0.09*(1+i), 10 );
366 				addNewO1(w* 0.375 + a , h*0.09*(1+i) , 10 );
367 				addNewO2(w* 0.375 - a, h*0.09*(1+i) , 10 ); 
368 				
369 				addNewC( w* 0.5 , h*0.09*(1+i), 10 );
370 				addNewO1(w* 0.5 + a , h*0.09*(1+i) , 10 );
371 				addNewO2(w* 0.5 - a, h*0.09*(1+i) , 10 ); 
372 				
373 				addNewC( w* 0.625 , h*0.09*(1+i), 10 );
374 				addNewO1(w* 0.625 + a , h*0.09*(1+i) , 10 );
375 				addNewO2(w* 0.625 - a, h*0.09*(1+i) , 10 ); 
376 				
377 				addNewC( w* 0.75 , h*0.09*(1+i), 10 );
378 				addNewO1(w* 0.75 + a , h*0.09*(1+i) , 10 );
379 				addNewO2(w* 0.75 - a, h*0.09*(1+i) , 10 ); 
380 				
381 				addNewC( w* 0.875 , h*0.09*(1+i), 10 );
382 				addNewO1(w* 0.875 + a , h*0.09*(1+i) , 10 );
383 				addNewO2(w* 0.875 - a, h*0.09*(1+i) , 10 ); 
384 				}
385 		}
386 		for (var i = 0; i < C.length; i++) 
387 				{
388 					C[i].vx = 0.4 * (1 - 2 * Math.random());
389                     C[i].vy = 0.4 * (1 - 2 * Math.random());
390 					O1[i].vx = 0.3 * (1 - 2 * Math.random());
391                     O1[i].vy = 0.3 * (1 - 2 * Math.random());
392 					O2[i].vx = 0.3 * (1 - 2 * Math.random());
393                     O2[i].vy = 0.3 * (1 - 2 * Math.random());
394 				}
395 		timers();
396     }
397 
398     document.onmouseup = function(e) {          // функция при отпускании клавиши мыши
399         canvas.onmousemove = null;              // когда клавиша отпущена - функции перемещения нету
400         dNd = null;                             // когда клавиша отпущена - захваченного курсором шара нету
401     };
402 
403     function mouseMove(e) {                     // функция при перемещении мыши, работает только с зажатой ЛКМ
404         var m = mouseCoords(e);                 // получаем расчетные координаты курсора мыши
405         dNd.x = m.x + dNd.xPlus;
406         dNd.y = m.y + dNd.yPlus;
407         dNd.vx = 0.6 * (m.x - mx_) / dt / fps;   dNd.vy = 0.6 * (m.y - my_) / dt / fps;
408         mx_ = m.x;    my_ = m.y;
409     }
410 
411     function mouseCoords(e) {                   // функция возвращает расчетные координаты курсора мыши
412         var m = [];
413         var rect = canvas.getBoundingClientRect();
414         m.x = (e.clientX - rect.left) / scale;
415         m.y = (e.clientY - rect.top) / scale;
416         return m;
417     }
418 
419     
420 
421 	var addNewBall =  function(x, y) {
422         // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
423         if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
424         for (var i = 0; i < balls.length; i++) {
425             var rx = balls[i].x - x;
426             var ry = balls[i].y - y;
427             var rLen2 = rx * rx + ry * ry;
428             if (rLen2 < 4 * r2) return null;
429         }
430         var b = [];
431 
432         b.x = x;                b.y = y;        // расчетные координаты шара
433         b.fx = 0;               b.fy = 0;      // сила, действующая на шар
434         b.vx = 0;               b.vy = 0;   		// скорость 
435 
436         balls[balls.length] = b;                // добавить элемент в конец массива
437         
438 		return b; 
439     };
440 	
441     var addNewC =  function(x, y) {
442         // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
443         if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
444         for (var i = 0; i < balls.length; i++) {
445             var rx = balls[i].x - x;
446             var ry = balls[i].y - y;
447             var rLen2 = rx * rx + ry * ry;
448             if (rLen2 < 4 * r2) return null;
449         }
450         var b = [];
451 
452         b.x = x;                b.y = y;        // расчетные координаты шара
453         b.fx = 0;               b.fy = 0;      // сила, действующая на шар
454         b.vx = 0;               b.vy = 0;   		// скорость
455 
456         C[C.length] = b;                // добавить элемент в конец массива
457         
458 		return b; 
459     };
460 
461     var addNewO1 =  function(x, y) {
462         // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
463         if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
464         for (var i = 0; i < balls.length; i++) {
465             var rx = balls[i].x - x;
466             var ry = balls[i].y - y;
467             var rLen2 = rx * rx + ry * ry;
468             if (rLen2 < 4 * r2) return null;
469         }
470         var b = [];
471 
472         b.x = x;                b.y = y;        // расчетные координаты шара
473         b.fx = 0;               b.fy = 0;      // сила, действующая на шар
474         b.vx = 0;               b.vy = 0;   		// скорость
475 		
476         O1[O1.length] = b;                // добавить элемент в конец массива
477         
478 		return b; 
479     };
480 	
481 	var addNewO2 =  function(x, y) {
482         // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
483         if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
484         for (var i = 0; i < balls.length; i++) {
485             var rx = balls[i].x - x;
486             var ry = balls[i].y - y;
487             var rLen2 = rx * rx + ry * ry;
488             if (rLen2 < 4 * r2) return null;
489         }
490         var b = [];
491 
492         b.x = x;                b.y = y;        // расчетные координаты шара
493         b.fx = 0;               b.fy = 0;      // сила, действующая на шар
494         b.vx = 0;               b.vy = 0;   		// скорость
495 
496         O2[O2.length] = b;                // добавить элемент в конец массива
497         
498 		return b; 
499     };
500 	
501 	// Основной цикл программы
502 
503     function control() {
504         physics();
505         draw();
506     }
507 	function timers(){
508 		clearInterval(id1);
509 		clearInterval(id2);
510 		clearInterval(id3);
511 		clearInterval(id4);
512 		clearInterval(id5);
513 		k = 1;
514 		Ekc = 0 ;
515 		Ekck = 0;
516 		EC= 0;
517 		EO1 = 0 ;
518 		EO2 = 0;
519 		ECk = 0;
520 		EO1k =0 ; 
521 		EO2k = 0 ;
522 		Ek = 0;
523 		Ek1 = 0 ;
524 		
525 		var id1 = setInterval(function(){document.getElementById('Ekck').innerHTML = Ekck;}, 1000 / 2);
526 		var id2 =setInterval(function(){document.getElementById('Ek').innerHTML = Ek;}, 1000 / 2);
527 		var id3 =setInterval(function(){document.getElementById('Ekc').innerHTML = ECk;}, 1000 / 2);
528 		var id4 =setInterval(function(){document.getElementById('Eko1').innerHTML = EO1k;}, 1000 / 2);
529 		var id5 =setInterval(function(){document.getElementById('Eko2').innerHTML = EO2k;}, 1000 / 2);
530 		
531 	}
532     // Расчетная часть программы
533 
534 	function powers(b ,b2 , hkk , jk) {
535 		var rx = b.x - b2.x;   var ry = b.y - b2.y;         // вектор смотрит на первый шар (b)
536         var r2 = rx * rx + ry * ry;                         // квадрат расстояния между шарами
537                   	
538 					if ( r2 > aCut2 ) return null;	
539 					
540 					var rLen = (Math.sqrt(r2));
541        
542 					if (r2 < K2a2) {
543                         if (rLen > 0.00001) {                           // проверка, чтобы избежать деления на 0
544                             rx = rx / rLen * Ka;
545                             ry = ry / rLen * Ka;
546                         }
547                         r2 = K2a2;
548                         rLen = Ka;                                      // корень K2a2
549                     }
550 
551                     // сила взаимодействия
552                    
553 					var u = Math.exp( -a1 * ( rLen - hkk  )) ;
554 					var F = jk  * MorzCoeff * u* ( u - 1 ) /rLen ;
555 					
556                     var Fx = F * rx;        var Fy = F * ry;
557                     b.fx += Fx;             b.fy += Fy;
558                     b2.fx -= Fx;            b2.fy -= Fy;	
559 	}
560 	
561 	
562     function physics() {                        // то, что происходит каждый шаг времени
563         for (var s = 1; s <= spf; s++) {
564 
565             // пересчет сил идет отдельным массивом, т.к. далее будут добавляться силы взаимодействия между шарами
566            
567 		   for (var i0 = 0; i0 < balls.length; i0++) {
568                 balls[i0].fx = - B * balls[i0].vx;
569                 balls[i0].fy =  - B * balls[i0].vy;
570             }
571 			for (var i0 = 0; i0 < C.length; i0++) {
572                 C[i0].fx = - B * C[i0].vx;
573                 C[i0].fy =  - B * C[i0].vy;
574             }
575 			for (var i0 = 0; i0 < O1.length; i0++) {
576                 O1[i0].fx = - B * O1[i0].vx;
577                 O1[i0].fy =  - B * O1[i0].vy;
578             }
579 			for (var i0 = 0; i0 < O2.length; i0++) {
580                 O2[i0].fx = - B * O2[i0].vx;
581                 O2[i0].fy =  - B * O2[i0].vy;
582             }
583 
584 			
585             for (var i = 0; i < balls.length; i++) { // dlya Balls
586                 
587 				var b = balls[i];
588                 
589 				for (var j = i + 1; j < balls.length; j++) {
590                     var b2 = balls[j];
591 					powers(b,b2,a_1,0.3) ;
592 				}
593 				 
594 				 for (var j = 0; j < C.length; j++) {
595                     var b2 = C[j];
596 					powers(b,b2,a_1,0.3) ;
597 				
598 				}
599 				 for (var j = 0; j < O1.length; j++) {
600                     var b2 = O1[j];
601 					powers(b,b2,a_1,0.3) ; 
602 				}
603 				
604 				 for (var j = 0; j < O2.length; j++) {
605                     var b2 = O2[j];
606 					powers( b, b2, a_1 ,0.3) ;
607 				}
608 				
609 				if (b == dNd) continue;  // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
610 
611                 if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
612                 if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
613                 if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
614                 if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }
615 
616                 b.vx += b.fx / m * dt;        b.vy += b.fy / m * dt;
617                 b.x += b.vx * dt;             b.y += b.vy * dt;
618 				Ek1 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy);
619             }
620 			
621 			for (var i = 0; i < C.length; i++) {
622 				for (var j = i ; j < C.length; j++) {
623 					
624 					if ( i == j){
625 						var r1x = O1[i].x - C[i].x ; var r1y = O1[i].y - C[i].y ;
626 						var r2x = O2[i].x - C[i].x ; var r2y = O2[i].y - C[i].y ; 
627 					
628 						var r1 = Math.sqrt(r1x*r1x +r1y* r1y);
629 						var r2 = Math.sqrt(r2x*r2x +r2y* r2y);
630 						var xx =(r1x*r2x + r1y*r2y)/(r1*r2) ;
631 					
632 						var F1 = -0.2*Cball*(r1-a)/r1;
633 						var F2 = -0.2*Cball*(r2-a)/r2;
634 						var f = 0 ;
635 						if(xx*xx > 0.99999999999){
636 							f = -0.3*Cball;	
637 							}
638 						
639 						else{
640 						f = -0.3*Cball*(Pi - Math.acos(xx))/(Math.sqrt(1-xx*xx));
641 						}
642 						ff = f ;
643 						
644 						var f1x = f*(r2x/(r1*r2) - xx*r1x/(r1*r1)); var f1y = f*(r2y/(r1*r2) -xx*r1y/(r1*r1));
645 						var f2x = f*(r1x/(r1*r2) - xx*r2x/(r2*r2)); var f2y = f*(r1y/(r1*r2) -xx*r2y/(r2*r2));
646 					
647 						var F1x = F1 * r1x   ;        var F1y = F1 * r1y ;
648 						var F2x = F2 * r2x   ;        var F2y = F2 * r2y ;
649 						
650 						O1[i].fx += (F1x + f1x) ;  O1[i].fy += (F1y + f1y) ;
651 						O2[i].fx += (F2x + f2x) ;  O2[i].fy += (F2y + f2y) ; 
652 						C[i].fx -= (F1x + F2x + f1x + f2x) ;  C[i].fy -= (F1y + F2y + f1y + f2y) ; 
653 						
654 						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) ;
655 					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) ;
656 					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) ;
657 						
658 						}
659 					
660 					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) ;
661 					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) ;
662 					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) ;
663 					
664 				}
665 			}
666 			
667 			
668 			
669 			for (var i = 0; i < C.length; i++) { // dlya C
670                  
671 				var b = C[i];
672                 
673 				if (b == dNd) continue;  // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
674 
675                 if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
676                 if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
677                 if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
678                 if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }
679 
680              
681                 b.vx += b.fx / m * dt;        b.vy += b.fy / m * dt;
682                 b.x += b.vx * dt;             b.y += b.vy * dt;
683 				Ek1 += m* 0.5 *(b.vx*b.vx + b.vy*b.vy) ;
684 				EC += m* 0.5 *(b.vx*b.vx + b.vy*b.vy)/j ;
685             }
686 			
687 			for (var i = 0; i < C.length; i++) { // dlya O1
688               
689 				var b = O1[i];
690      
691 				if (b == dNd) continue;  // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
692 
693                 if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
694                 if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
695                 if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
696                 if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }
697 
698                 b.vx += b.fx / m1 * dt;        b.vy += b.fy / m1 * dt;
699                 b.x += b.vx * dt;             b.y += b.vy * dt;
700 				Ek1 += m1 * 0.5 *(b.vx*b.vx + b.vy * b.vy) ;
701 				EO1 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy)/j ;
702 			
703                 
704             }
705 			
706 			for (var i = 0; i < C.length; i++) { // dlya O2
707                 
708 				var b = O2[i];
709                 
710 				 
711 				if (b == dNd) continue;  // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
712 
713                 if (b.y + r > h) { b.fy += -Cwall * (b.y + r - h) - Bwall * b.vy; }
714                 if (b.y - r < 0) { b.fy += -Cwall * (b.y - r) - Bwall * b.vy;}
715                 if (b.x + r > w) { b.fx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
716                 if (b.x - r < 0) { b.fx += -Cwall * (b.x - r) - Bwall * b.vx; }
717 
718                 b.vx += b.fx / m1 * dt;        b.vy += b.fy / m1 * dt;
719                 b.x += b.vx * dt;             b.y += b.vy * dt;
720 				Ek1 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy) ;
721 				EO2 += m1* 0.5 *(b.vx*b.vx + b.vy*b.vy)/j ;
722 				
723 				
724             
725 			}
726 			 
727         }
728 		Ekc +=Ek1 ;
729 		
730 		if(Ek1 > 1){
731 			k++;
732 			}
733 		Ekck = Ekc/k;
734 		ECk = EC/k ;
735 		EO1k = EO1/k ; 
736 		EO2k = EO2/k ;
737 		Ek = Ek1;
738 		Ek1 = 0 ;
739     }
740 	
741    
742    
743 
744    // Рисование
745 
746     var rScale13 = r * scale * 1.3;         // ___в целях оптимизации___
747     var rScaleShift = r * scale / 5;        // ___в целях оптимизации___
748     function draw() {
749         context.clearRect(0, 0, w * scale, h * scale);      // очистить экран
750         for (var i = 0; i < balls.length; i++){
751             var xS = balls[i].x * scale;           var yS = balls[i].y * scale;
752             
753             context.fillStyle = "#FFD700";
754 
755             context.beginPath();
756             context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
757             context.closePath();
758             context.fill();
759         }
760     
761 		for (var i = 0; i < C.length; i++){
762             var xS = C[i].x * scale;           var yS = C[i].y * scale;
763             
764             context.fillStyle = "#f08080";
765 
766             context.beginPath();
767             context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
768             context.closePath();
769             context.fill();
770         }
771 		for (var i = 0; i < O1.length; i++){
772             var xS = O1[i].x * scale;           var yS = O1[i].y * scale;
773             
774             context.fillStyle = "#3070d0";
775 
776             context.beginPath();
777             context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
778             context.closePath();
779             context.fill();
780         }
781 		for (var i = 0; i < O2.length; i++){
782             var xS = O2[i].x * scale;           var yS = O2[i].y * scale;
783             
784             context.fillStyle = "#3070d0";
785 
786             context.beginPath();
787             context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
788             context.closePath();
789             context.fill();
790         }
791 	}
792 
793     // Запуск системы
794 	
795 	this.newSystem();
796 	var id1 = setInterval(function(){document.getElementById('Ekck').innerHTML = Ekck;}, 1000 / 2);
797 	var id2 =setInterval(function(){document.getElementById('Ek').innerHTML = Ek;}, 1000 / 2);
798 	var id3 =setInterval(function(){document.getElementById('Ekc').innerHTML = ECk;}, 1000 / 2);
799 	var id4 =setInterval(function(){document.getElementById('Eko1').innerHTML = EO1k;}, 1000 / 2);
800 	var id5 =setInterval(function(){document.getElementById('Eko2').innerHTML = EO2k;}, 1000 / 2);
801 	
802     setInterval(control, 1000 / fps);
803 }

Файл "Lab.html" 

 1 <!DOCTYPE html>
 2 <html>
 3 <head>
 4     <title>Balls</title>
 5     <script src="Balls_v4_release.js"></script>
 6 </head>
 7 <body>
 8     <canvas id="canvasBalls" width="640" height="480" style="border:1px solid #000000;"></canvas>
 9 	<br>
10 	 <div>Количество частиц:
11         <input type="button" style="width: 30px" name="" onclick="setNy(8); newSystem();return false;" value="8"/>
12         <input type="button" style="width: 30px" name="" onclick="setNy(20); newSystem();return false;" value="20"/>
13         <input type="button" style="width: 30px" name="" onclick="setNy(40); newSystem();return false;" value="40"/>
14         <input type="button" style="width: 30px" name="" onclick="setNy(70); newSystem();return false;" value="70"/>
15     </div>
16 	 <div>Обновить таймер:
17         <input type="button" style="width: 60px" name="" onclick="ResetTime();return false;" value="Reset"/>
18     </div>
19 	<div>Кин энергия в данный момент времени: <span id="Ek"></span></div>
20 	<div>Средняя кинетическая энергия: <span id="Ekck"></span></div>
21 	<div>Средняя кинетическая энергия на один атом углерода: <span id="Ekc"></span></div>
22 	<div>Средняя кинетическая энергия на один атом кислорода(1): <span id="Eko1"></span></div>
23 	<div>Средняя кинетическая энергия на один атом кислорода(2): <span id="Eko2"></span></div>
24 </body>
25 </html>

Результаты[править]

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

[math]U= \frac{kT}{2} [/math]

где

  • [math]k[/math] — постоянная Больцмана,
  • [math]T[/math] — температура среды,

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

Источник — «http://tm.spbstu.ru/?title=Моделирование_углекислого_газа_методом_динамики_частиц&oldid=45780»