Моделирование углекислого газа методом динамики частиц
Материал из Department of Theoretical and Applied Mechanics
Курсовые работы 2016-2017 учебного года > Моделирование углекислого газа методом динамики частиц
Курсовой проект по Механике дискретных сред
Исполнитель: Смирнов Александр
Группа: 09 (43604/1)
Семестр: осень 2016
Содержание
Формулировка задачи
Смоделировать молекулы углекислого газа методом динамики частиц, проверить выполнение закона сохранения энергии и рассмотреть распределение энергии по степеням свободы.
Общие сведения
Для моделирования используем частицы, которые представляют собой абсолютно упругий шар. Масса углерода равна 12, а кислорода - 16 условных единиц. Для отталкивания молекул используется потенциал Морзе, внутри молекулы: Упругая сила и угловая пружина.
Упругая сила и угловая пружина
Упругая сила определяется формулой:
где
- — жесткость,
- — отклонение от положения равновесия.
Потенциал пружины равен:
При этом
где
- —угловая жесткость,
- — радиус вектор от атома углерода к первому или второму кислороду соответственно.
А сила вычисляется как:
Потенциал Морзе
Парный силовой потенциал взаимодействия. Определяется формулой:
где
- — энергия связи,
- — длина связи,
- — параметр, характеризующий ширину потенциальной ямы.
Сила, соответствующая потенциалу Морзе, вычисляется по формуле:
Или в векторной форме:
Программа
В данной программе в начальный момент времени системе задаются случайные скорости(начальная энергия,они достаточно велики, чтобы можно было пренебречь потенциальной энергией взаимодействия. Можно менять количество молекул углекислого газа и сбрасывать таймер расчета средних значений. Так же выводятся:кинетическая энергия системы в данный момент времени, средняя кинетическая энергия системы в данный момент времени, средние энергии, приходящиеся на атом углерода, первый и второй атом кислорода в молекуле.
- Красные шары - атомы углерода
- Синие шары - атомы кислорода
Текст программы на языке 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>
Результаты
Получена программа, которая моделирует молекулы углекислого газа в объеме и рассчитывает среднюю кинетическую энергию отдельных частей системы и всей системы в целом. Значение кинетической энергии системы совершает колебания вокруг среднего значения энергии, которое оказалось постоянным, что говорит о выполнении закона сохранения энергии. Средние энергии, приходящиеся на каждый отдельный атом, спустя какое-то время оказываются равными. Это говорит о том, что энергии, проходящиеся на каждую степень свободы, равны. Для идеального газа эта энергия равна:
где
- — постоянная Больцмана,
- — температура среды,
А нашу систему можно приблизительно считать идеальным газом в силу малости значения потенциальной энергии и абсолютно упругих соударений. Значит наши результаты соответствуют, теории.
