Динамика взаимодействующих частиц
Материал из Department of Theoretical and Applied Mechanics
Виртуальная лаборатория > Динамика взаимодействующих частиц
Здесь расположена последняя версия программы, моделирующей динамику взаимодуйствующих частиц. Каждая частица представляет из себя вязкоупругий шар. Взаимодействие между шарами описывается потенциалом Леннарда-Джонса
Гравитация:
mg =
⋅ m ⋅ g0
Сколько шаров помещается по вертикали:
Конфигурация:
Short Lennard-Jones potential
Термостат:
T ⋅ T0 =
Разгон случайными скоростями
Термостат действует на: Внешнее трение Внутреннее трение
Разгон случайными скоростями
Термостат действует на: Внешнее трение Внутреннее трение
T ≈
Количество частиц:
Скачать программу: Balls_v6_release.zip
Текст программы на языке JavaScript (разработчики Кривцов Антон, Цветков Денис):
Файл "Balls_v6_release.js"
1 function MainBalls(canvas, slider_01, text_01, slider_02, text_02) {
2
3 canvas.onselectstart = function () {return false;}; // запрет выделения canvas
4
5 // Предварительные установки
6
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 = 1 * m0; // масса
25 var CWall = 10 * C0; // жесткость стен
26 var CBall = 0.1 * CWall; // жесткость между частицами
27 var BVisc = 0.008 * B0; // вязкость среды
28 var BInternal = 0.01 * B0; // внутреннее трение
29 var BWall = 0.03 * B0; // вязкость на стенках
30 var mg = 0.25 * m * g0; // сила тяжести
31 var r = 0.5 * a0; // радиус частицы в расчетных координатах
32 var K = 0.7; // все силы, зависящие от радиуса, ограничиваются значением, реализующимся при r/a = K
33 var a = 2 * r; // равновесное расстояние между частицами
34 var aCut = 2 * a; // радиус обрезания
35 var TGoalK = 2; // целевая температура системы равна TGoalK * D
36 var TActualMaxK = 200; // макимальная температура, при которой работает термостат равна TActualMaxK * D
37
38 // *** Задание вычислительных параметров ***
39
40 var fps = 60; // frames per second - число кадров в секунду (качечтво отображения)
41 var spf = 100; // steps per frame - число шагов интегрирования между кадрами (скорость расчета)
42 var dt = 0.04 * t0 / fps; // шаг интегрирования (качество расчета)
43
44 // Выполнение программы
45
46 var r2 = r * r; // ___в целях оптимизации___
47 var a2 = a * a; // ___в целях оптимизации___
48 var D = a2 * CBall / 72; // энергия связи между частицами
49 var LJCoeff = 12 * D / a2; // коэффициент для расчета потенциала Л-Дж
50 var b = Math.pow(13 / 7, 6) * a; // коэффициент для SLJ потенциала
51 var b2 = b * b; // ___в целях оптимизации___
52 var SLJDenominator = 1 / (aCut * aCut - b2); // знаменатель для расчета SLJ потенциала
53
54 var thermostatEnabled = document.getElementById('checkbox_02').checked; // термостат применяется к вязкости среды
55 var addRandomV = document.getElementById('checkbox_03').checked; // случайные скорости для разгона
56 var T0 = 1 * D; // масштаб температуры
57 var TGoal = TGoalK * T0; // целевая температура системы
58 var TActualMax = TActualMaxK * T0; // макимальная температура, при которой работает термостатс (для избежания беск. скоростей)
59 var TActual = 0; // актуальная температура
60 var k = 1; // постоянную Больцмана примем за единицу
61 var Tk = m / k; // ___в целях оптимизации___
62 var viscFrictionTh = document.getElementById('checkbox_04').checked; // термостат применяется к вязкости среды
63 var internalFrictionTh = document.getElementById('checkbox_05').checked; // термостат применяется к внутреннему трению
64 var TempIntervalID;
65
66 var Ka = K * a; // ___в целях оптимизации___
67 var K2a2 = K * K * a2; // ___в целях оптимизации___
68
69 var dNd = null; // ссылка на захваченный курсором шар (drag & drop)
70 var grad; // должен ли работать градиент (регулируется в функции setNy())
71 var SLJEnabled = document.getElementById('checkbox_01').checked;
72
73 this.setSlider_01 = function(c) {mg = c * m * g0;}; // функция для слайдера гравитации
74 this.setSlider_02 = function(c) {TGoal = c;}; // функция для слайдера термостата
75 this.setNy = function(ny) {
76 Ny = ny;
77 if (Ny > 8) {
78 grad = false; // градиент не работает, если Ny > 8
79 context.fillStyle = "#3070d0"; // цвет, шара
80 } else
81 grad = true;
82 };
83 this.setNy(Ny); // запускаем с уже присвоенным значением, чтобы обновились настройки градиента
84 this.setCheckbox_01 = function(bool) {SLJEnabled = bool;};
85 this.setCheckbox_02 = function(bool) {
86 thermostatEnabled = bool;
87 document.getElementById('checkbox_03').disabled = !bool;
88 document.getElementById('checkbox_04').disabled = !bool;
89 document.getElementById('checkbox_05').disabled = !bool;
90 document.getElementById('slider_02').disabled = !bool;
91 document.getElementById('text_02').disabled = !bool;
92 if (bool) {
93 TempIntervalID = setInterval( // обновление информации о температуре
94 function(){document.getElementById('Temperature').innerHTML = TActual.toFixed(3);}, 1000 / 3);
95 }
96 else {
97 clearInterval(TempIntervalID); // температура больше не подсчитывается - удаляем обновление информации о ней
98 document.getElementById('Temperature').innerHTML = "???"
99 }
100 };
101 this.setCheckbox_02(thermostatEnabled); // запускаем сразу, чтобы обновить состояния элементов интерфейса
102 this.setCheckbox_03 = function(bool) {addRandomV = bool;};
103 this.setCheckbox_04 = function(bool) {viscFrictionTh = bool;};
104 this.setCheckbox_05 = function(bool) {internalFrictionTh = bool;};
105
106 // Настройка интерфейса
107
108 slider_01.min = 0; slider_01.max = 5;
109 slider_01.step = 0.05;
110 slider_01.value = mg / m / g0; // начальное значение ползунка должно задаваться после min и max
111 text_01.value = mg / m / g0;
112 slider_02.min = 0; slider_02.max = 5;
113 slider_02.step = 0.05;
114 slider_02.value = TGoal; // начальное значение ползунка должно задаваться после min и max
115 text_02.value = TGoal.toFixed(1);
116
117 // Запуск новой системы
118
119 // следующие переменные должны пересчитываться каждый раз, когда мы изменяем значение Ny
120 var scale, w, h;
121 var rScale13, rScaleShift;
122 this.newSystem = function() {
123 scale = canvas.height / Ny / a0; // масштабный коэффициент для перехода от расчетных к экранным координатам
124 w = canvas.width / scale; // ширина окна в расчетных координатах
125 h = canvas.height / scale; // высота окна в расчетных координатах
126
127 rScale13 = r * scale * 1.3; // ___в целях оптимизации___
128 rScaleShift = r * scale / 5; // ___в целях оптимизации___
129
130 this.setRandom(); // задаем случайную конфигурацию
131 };
132
133 // Работа с мышью
134
135 var mx_, my_; // буфер позиции мыши (для расчета скорости при отпускании шара)
136
137 canvas.onmousedown = function(e) { // функция при нажатии клавиши мыши
138 var m = mouseCoords(e); // получаем расчетные координаты курсора мыши
139 // цикл в обратную сторону, чтобы захватывать шар, нарисованный "сверху"
140 // (т.к. цикл рисования идет в обычном порядке)
141 for (var i = balls.length - 1; i >= 0; i--) {
142 var b = balls[i];
143 var rx = b.x - m.x;
144 var ry = b.y - m.y;
145 var rLen2 = rx * rx + ry * ry; // квадрат расстояния между курсором и центром шара
146 if (rLen2 <= r2) { // курсор нажал на шар
147 if (e.which == 1) { // нажата левая клавиша мыши
148 dNd = b;
149 dNd.xPlus = dNd.x - m.x; // сдвиг курсора относительно центра шара по x
150 dNd.yPlus = dNd.y - m.y; // сдвиг курсора относительно центра шара по y
151 mx_ = m.x; my_ = m.y;
152 canvas.onmousemove = mouseMove; // пока клавиша нажата - работает функция перемещения
153 } else if (e.which == 3) // нажата правая клавиша мыши
154 balls.splice(i, 1); // удалить шар
155 return;
156 }
157 }
158
159 // если не вышли по return из цикла - нажатие было вне шара, добавляем новый
160 if (e.which == 1) {
161 dNd = addNewBall(m.x, m.y, true); // добавляем шар и сразу захватываем его курсором
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 document.onmouseup = function(e) { // функция при отпускании клавиши мыши
170 canvas.onmousemove = null; // когда клавиша отпущена - функции перемещения нету
171 dNd = null; // когда клавиша отпущена - захваченного курсором шара нету
172 };
173
174 function mouseMove(e) { // функция при перемещении мыши, работает только с зажатой ЛКМ
175 var m = mouseCoords(e); // получаем расчетные координаты курсора мыши
176 dNd.x = m.x + dNd.xPlus;
177 dNd.y = m.y + dNd.yPlus;
178 dNd.vx = 0.6 * (m.x - mx_) / dt / fps; dNd.vy = 0.6 * (m.y - my_) / dt / fps;
179 mx_ = m.x; my_ = m.y;
180 }
181
182 function mouseCoords(e) { // функция возвращает расчетные координаты курсора мыши
183 var m = [];
184 var rect = canvas.getBoundingClientRect();
185 m.x = (e.clientX - rect.left) / scale;
186 m.y = (e.clientY - rect.top) / scale;
187 return m;
188 }
189
190 // Работа с массивом
191
192 var balls = []; // массив шаров
193 var addNewBall = function(x, y, check) {
194 // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
195 if (check) {
196 if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
197 for (var i = 0; i < balls.length; i++) {
198 var rx = balls[i].x - x;
199 var ry = balls[i].y - y;
200 var rLen2 = rx * rx + ry * ry;
201 if (rLen2 < 4 * r2) return null;
202 }
203 }
204
205 var b = [];
206
207 b.x = x; b.y = y; // расчетные координаты шара
208 b.fx = 0; b.fy = mg; // сила, действующая на шар
209 b.vx = 0; b.vy = 0; // скорость
210
211 balls[balls.length] = b; // добавить элемент в конец массива
212 return b;
213 };
214
215 this.setEmpty = function() {balls = [];}; // пустое поле
216
217 this.setRandom = function() { // случайная конфигурация
218 balls = [];
219 for (var i = 0; i < 1000; i++)
220 addNewBall(Math.random() * w, Math.random() * h, true);
221 };
222
223 var sqrt3 = Math.sqrt(3);
224 this.setTriangularLattice = function() { // задать на поле треугольную решетку
225 balls = [];
226 var center = (w - Math.floor(w / r) * r) / 2; // сдвиг, решетка будет появляться по середине по горизонтали
227 for (var j = 0; j < Math.floor(h / (sqrt3 * r)); j++)
228 for (var i = 0; i < Math.floor(w / r) - 1; i++)
229 if ((i + j) % 2 == 0) addNewBall(r * (i + 1) + center, h - r * (1 + sqrt3 * j), false);
230 };
231
232 // Основной цикл программы
233
234 function control() {
235 physics();
236 draw();
237 }
238
239 // Расчетная часть программы
240
241 function physics() { // то, что происходит каждый шаг времени
242 for (var s = 1; s <= spf; s++) {
243
244 var BViscTh = BVisc;
245 var BInternalTh = BInternal;
246 // работа термостата
247 if (thermostatEnabled) {
248 if (balls.length > 0) {
249 var v2Sum = 0;
250 for (var i1 = 0; i1 < balls.length; i1++)
251 v2Sum += balls[i1].vx * balls[i1].vx + balls[i1].vy * balls[i1].vy;
252 var v2Average = v2Sum / balls.length;
253 TActual = Tk * v2Average;
254
255 if (addRandomV) { // случайные скорости, если температура слишком мала
256 if (TGoal > 0.15 && TActual < 0.1) {
257 for (var i2 = 0; i2 < balls.length; i2++) {
258 balls[i2].vx += 0.3 * (1 - 2 * Math.random());
259 balls[i2].vy += 0.3 * (1 - 2 * Math.random());
260 }
261 }
262 }
263
264 if (TActual < TActualMax) { // из за того, что мышкой можно задать шарам запредельную скорость
265 if (viscFrictionTh) BViscTh = BVisc * (TActual - TGoal); // действие термостата
266 if (internalFrictionTh) BInternalTh = BInternal * (TActual - TGoal); // действие термостата
267 }
268 } else
269 TActual = 0; // для датчика температуры на странице
270 }
271
272 // пересчет сил идет отдельным массивом, т.к. далее будут добавляться силы взаимодействия между шарами
273 for (var i0 = 0; i0 < balls.length; i0++) {
274 balls[i0].fx = - BViscTh * balls[i0].vx;
275 balls[i0].fy = mg - BViscTh * balls[i0].vy;
276 }
277
278 for (var i = 0; i < balls.length; i++) {
279 // расчет взаимодействия производится со всеми следующими шарами в массиве,
280 // чтобы не считать каждое взаимодействие дважды
281 var b = balls[i];
282 for (var j = i + 1; j < balls.length; j++) {
283 var b2 = balls[j];
284 var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b)
285 var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
286 var rLen = (Math.sqrt(r2));
287 if (rLen > aCut) continue; // проверка на радиус обрезания
288
289 // если расстояние между частицами мало, силы будут посчитаны для K * a
290 if (rLen < Ka) {
291 if (rLen > 0.00001) { // проверка, чтобы избежать деления на 0
292 rx = rx / rLen * Ka;
293 ry = ry / rLen * Ka;
294 }
295 r2 = K2a2;
296 rLen = Ka; // корень K2a2
297 }
298
299 // сила взаимодействия
300 var s2 = a2 / r2; var s4 = s2 * s2; // ___в целях оптимизации___
301 var F = LJCoeff * s4 * s4 * (s4 * s2 - 1); // сила взаимодействия Леннарда-Джонса
302 if (SLJEnabled) {
303 var kSLJ; // k(r) - сглаживающий коэффициент SLJ потенциала
304 if (r <= b) kSLJ = 1;
305 else {
306 var brackets = (r2 - b2) * SLJDenominator;
307 kSLJ = 1 - brackets * brackets;
308 } // случай rLen > aCut обработан выше
309 F *= kSLJ;
310 }
311
312 // сила внутреннего трения между частицами
313 if (r2 < a2) {
314 var vx21 = b.vx - b2.vx; var vy21 = b.vy - b2.vy; // вектор смотрит на первый шар (b)
315 var ex = rx / rLen; var ey = ry / rLen;
316 var v = vx21 * ex + vy21 * ey;
317 F -= F * BInternalTh / rLen * v;
318 }
319
320 // суммируем силы
321 var Fx = F * rx; var Fy = F * ry;
322
323 b.fx += Fx; b.fy += Fy;
324 b2.fx -= Fx; b2.fy -= Fy;
325 }
326
327 if (b == dNd) continue; // если шар схвачен курсором - его вз. со стенами и перемещение не считаем
328
329 if (b.y + r > h) { b.fy += -CWall * (b.y + r - h) - BWall * b.vy; }
330 if (b.y - r < 0) { b.fy += -CWall * (b.y - r) - BWall * b.vy;}
331 if (b.x + r > w) { b.fx += -CWall * (b.x + r - w) - BWall * b.vx; }
332 if (b.x - r < 0) { b.fx += -CWall * (b.x - r) - BWall * b.vx; }
333
334 b.vx += b.fx / m * dt; b.vy += b.fy / m * dt;
335 b.x += b.vx * dt; b.y += b.vy * dt;
336 }
337 }
338 }
339
340 // Рисование
341 function draw() {
342 context.clearRect(0, 0, w * scale, h * scale); // очистить экран
343 for (var i = 0; i < balls.length; i++){
344 var xS = balls[i].x * scale; var yS = balls[i].y * scale;
345 if (grad) {
346 // расчет градиента нужно проводить для каждого шара
347 var gradient = context.createRadialGradient(xS, yS, rScale13, xS - rScaleShift, yS + rScaleShift, 0);
348 gradient.addColorStop(0, "#0000bb");
349 gradient.addColorStop(1, "#44ddff");
350 context.fillStyle = gradient;
351 }
352
353 context.beginPath();
354 context.arc(xS, yS, r * scale, 0, 2 * Math.PI, false);
355 context.closePath();
356 context.fill();
357 }
358 }
359
360 // Запуск системы
361 this.newSystem();
362 setInterval(control, 1000 / fps);
363 // след. функция обновляет информацию о количестве частиц на поле
364 setInterval(function(){document.getElementById('ballsNum').innerHTML = balls.length;}, 1000 / 20);
365 }
Файл "Balls_v6_release.html"
1 <!DOCTYPE html>
2 <html>
3 <head>
4 <meta charset="UTF-8" />
5 <title>Balls</title>
6 <script src="Balls_v6_release.js"></script>
7 </head>
8 <body>
9 <canvas id="canvasBalls" width="800" height="600" style="border:1px solid #000000;"></canvas>
10 <br>
11 <div>Гравитация:
12 <input type="range" id="slider_01" style="width: 150px;" oninput="app.setSlider_01(this.value); document.getElementById('text_01').value = this.value;">
13 mg =
14 <input id="text_01" style="width: 5ex;" required pattern="[-+]?([0-9]*\.[0-9]+|[0-9]+)" oninput="
15 // если введено не число - строка не пройдет валидацию по паттерну выше, и checkValidity() вернет false
16 if (!this.checkValidity()) return;
17 app.setSlider_01(this.value);
18 document.getElementById('slider_01').value = this.value;
19 ">
20 ⋅ m ⋅ g0</div><br>
21
22 <div>Сколько шаров помещается по вертикали:
23 <input type="button" style="width: 30px" name="" onclick="app.setNy(3); app.newSystem();return false;" value="3"/>
24 <input type="button" style="width: 30px" name="" onclick="app.setNy(4); app.newSystem();return false;" value="4"/>
25 <input type="button" style="width: 30px" name="" onclick="app.setNy(5); app.newSystem();return false;" value="5"/>
26 <input type="button" style="width: 30px" name="" onclick="app.setNy(7); app.newSystem();return false;" value="7"/>
27 <!--знак пробела-->
28 <input type="button" style="width: 30px" name="" onclick="app.setNy(9); app.newSystem();return false;" value="9"/>
29 <input type="button" style="width: 30px" name="" onclick="app.setNy(12); app.newSystem();return false;" value="12"/>
30 <input type="button" style="width: 30px" name="" onclick="app.setNy(15); app.newSystem();return false;" value="15"/>
31 </div><br>
32
33 <div>Конфигурация:
34 <input type="button" name="" onclick="app.setTriangularLattice(); return false;" value="Треугольная решетка"/>
35 <input type="button" name="" onclick="app.setRandom(); return false;" value="Как попало"/>
36 <input type="button" name="" onclick="app.setEmpty(); return false;" value="Пустое поле"/>
37 </div><br>
38
39 <div>
40 <input type="checkbox" id="checkbox_01" name="" onchange="app.setCheckbox_01(this.checked);"/>
41 <a href="/SLJ" title="SLJ" class="mw-redirect">Short Lennard-Jones</a> potential
42 </div><br>
43
44 <div>
45 <input type="checkbox" id="checkbox_02" name="" onchange="app.setCheckbox_02(this.checked);"/>
46 Термостат:
47 <input type="range" id="slider_02" style="width: 150px;" oninput="app.setSlider_02(this.value); document.getElementById('text_02').value = this.value;">
48 T ⋅ T0 = <input id="text_02" style="width: 5ex;" required pattern="[-+]?([0-9]*\.[0-9]+|[0-9]+)" oninput="
49 // если введено не число - строка не пройдет валидацию по паттерну выше, и checkValidity() вернет false
50 if (!this.checkValidity()) return;
51 app.setSlider_02(this.value);
52 document.getElementById('slider_02').value = this.value;
53 ">
54 <br>
55 <input type="checkbox" checked id="checkbox_03" name="" onchange="app.setCheckbox_03(this.checked);"/>Разгон случайными скоростями
56 <br>
57 Термостат действует на:
58 <input type="checkbox" checked id="checkbox_04" name="" onchange="app.setCheckbox_04(this.checked);"/>Внешнее трение
59 <input type="checkbox" checked id="checkbox_05" name="" onchange="app.setCheckbox_05(this.checked);"/>Внутреннее трение
60 <div>T ≈ <span id="Temperature"></span></div>
61 </div><br>
62
63 <div>Количество частиц: <span id="ballsNum"></span></div>
64
65 <script type="text/javascript">var app = new MainBalls(
66 document.getElementById('canvasBalls'),
67 document.getElementById('slider_01'),
68 document.getElementById('text_01'),
69 document.getElementById('slider_02'),
70 document.getElementById('text_02')
71 );</script>
72 </body>
73 </html>
Здесь вы можете найти предыдущие версии программы.
Предлагаемые направления развития стенда
- Добавить возможность задавать температуру при помощи термостата Хувера.
- Исследовать скорость работы программы, оптимизировать.
- Применить принцип копирования.
- Добавить различные потенциалы.
- Добавить возможность задавать упруговязкие стенки в любом месте.
- Добавить возможность включения периодической системы.
