Обратный каскад энергии (двумерная турбулентность) — различия между версиями
(→Ссылки) |
(→Анализ) |
||
(не показано 5 промежуточных версий этого же участника) | |||
Строка 1: | Строка 1: | ||
+ | |||
+ | |||
+ | |||
'''''Курсовой проект по [[Механика дискретных сред|Механике дискретных сред]]''''' | '''''Курсовой проект по [[Механика дискретных сред|Механике дискретных сред]]''''' | ||
Строка 49: | Строка 52: | ||
Частички расположены впритык друг к другу и образуют прямоугольную решетку. С помощью написанной программы можно изменять количество шаров, коэффициент вязкости и температуру в среде(термостат). | Частички расположены впритык друг к другу и образуют прямоугольную решетку. С помощью написанной программы можно изменять количество шаров, коэффициент вязкости и температуру в среде(термостат). | ||
− | {{#widget:Iframe |url=http://tm.spbstu.ru/htmlets/Stepanov/inverse% | + | {{#widget:Iframe |url=http://tm.spbstu.ru/htmlets/Stepanov/inverse%20cascade4_19.12.16/newVar.html |width=950 |height=1900 |border=0 }} |
<div class="mw-collapsible mw-collapsed"> | <div class="mw-collapsible mw-collapsed"> | ||
'''Текст программы на языке JavaScript:''' <div class="mw-collapsible-content"> | '''Текст программы на языке JavaScript:''' <div class="mw-collapsible-content"> | ||
Строка 56: | Строка 59: | ||
// Версия 16.12.16 | // Версия 16.12.16 | ||
+ | |||
+ | |||
+ | // Версия Total 19.12.16 | ||
+ | //Stepanov MD | ||
function MainBalls(canvas, slider_01, text_01, slider_02, text_02) { | function MainBalls(canvas, slider_01, text_01, slider_02, text_02) { | ||
Строка 72: | Строка 79: | ||
var a0 = 1; // масштаб расстояния (диаметр шара) | var a0 = 1; // масштаб расстояния (диаметр шара) | ||
+ | var g0 = a0 / T0 / T0; // масштаб ускорения (ускорение, при котором за T0 будет пройдено расстояние a0) | ||
+ | var k0 = 2 * Pi / T0; // масштаб частоты | ||
var C0 = m0 * k0 * k0; // масштаб жесткости | var C0 = m0 * k0 * k0; // масштаб жесткости | ||
var B0 = 2 * m0 * k0; // масштаб вязкости | var B0 = 2 * m0 * k0; // масштаб вязкости | ||
− | + | ||
− | + | var L = 0; | |
− | + | var G = 0; | |
− | + | var G1 = 0; | |
+ | var G2 = 0; | ||
// *** Задание физических параметров *** | // *** Задание физических параметров *** | ||
− | var Ny = | + | var Ny = 20; // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна) |
var m = 1 * m0; // масса | var m = 1 * m0; // масса | ||
var Cwall = 10 * C0; // жесткость стен | var Cwall = 10 * C0; // жесткость стен | ||
var Cball = 0.1 * Cwall; // жесткость между частицами | var Cball = 0.1 * Cwall; // жесткость между частицами | ||
+ | |||
+ | /*var B = 0.000 * B0; // вязкость среды*/ | ||
var Bball = 0.01 * B0; // вязкость между частицами | var Bball = 0.01 * B0; // вязкость между частицами | ||
var Bwall = 0.03 * B0; // вязкость на стенках | var Bwall = 0.03 * B0; // вязкость на стенках | ||
+ | var mg = 0.25 * m * g0; // сила тяжести | ||
var r = 0.5 * a0; // радиус частицы в расчетных координатах | var r = 0.5 * a0; // радиус частицы в расчетных координатах | ||
var K = 0.7; // все силы, зависящие от радиуса, ограничиваются значением, реализующимся при r/a = K | var K = 0.7; // все силы, зависящие от радиуса, ограничиваются значением, реализующимся при r/a = K | ||
Строка 98: | Строка 111: | ||
var spf = 100; // steps per frame - число шагов интегрирования между кадрами (скорость расчета) | var spf = 100; // steps per frame - число шагов интегрирования между кадрами (скорость расчета) | ||
var dt = 0.045 * T0 / fps; // шаг интегрирования (качество расчета) | var dt = 0.045 * T0 / fps; // шаг интегрирования (качество расчета) | ||
− | + | ||
+ | |||
+ | var steps = 0; | ||
+ | |||
+ | |||
// Выполнение программы | // Выполнение программы | ||
Строка 105: | Строка 122: | ||
var a2 = a * a; // ___в целях оптимизации___ | var a2 = a * a; // ___в целях оптимизации___ | ||
var D = a2 * Cball / 72; // энергия связи между частицами | var D = a2 * Cball / 72; // энергия связи между частицами | ||
− | + | var Temp = 0; | |
− | + | Temp0 = 0.12 * D; | |
− | + | var B1 = 12*Math.sqrt((2*D)/a2); | |
− | + | var B = 0.026 * B1; | |
− | + | var betta = 3.5*B; | |
− | + | ||
− | + | var LJCoeff = 6 * D / a2; // коэффициент для расчета потенциала Л-Дж | |
− | |||
− | |||
var Ka = K * a; // ___в целях оптимизации___ | var Ka = K * a; // ___в целях оптимизации___ | ||
var K2a2 = K * K * a2; // ___в целях оптимизации___ | var K2a2 = K * K * a2; // ___в целях оптимизации___ | ||
Строка 120: | Строка 135: | ||
var grad; // должен ли работать градиент (регулируется в функции setNy()) | var grad; // должен ли работать градиент (регулируется в функции setNy()) | ||
− | + | this.set_01 = function(p) {betta = p * B;}; | |
− | + | this.set_02 = function(pT) {Temp0 = pT *D;}; | |
− | + | ||
this.setNy = function(ny) { | this.setNy = function(ny) { | ||
Ny = ny; | Ny = ny; | ||
Строка 145: | Строка 160: | ||
}; | }; | ||
− | + | // настройка слайдеров и текстовых полей | |
slider_01.min = 0.1; slider_01.max = 20; | slider_01.min = 0.1; slider_01.max = 20; | ||
slider_01.step = 0.1; | slider_01.step = 0.1; | ||
Строка 151: | Строка 166: | ||
text_01.value = betta / B; | text_01.value = betta / B; | ||
− | + | // настройка слайдеров и текстовых полей | |
slider_02.min = 0.1; slider_02.max = 1; | slider_02.min = 0.1; slider_02.max = 1; | ||
slider_02.step = 0.1; | slider_02.step = 0.1; | ||
slider_02.value = Temp0 / D; // начальное значение ползунка должно задаваться после min и max | slider_02.value = Temp0 / D; // начальное значение ползунка должно задаваться после min и max | ||
text_02.value = Temp0 / D; | text_02.value = Temp0 / D; | ||
− | + | ||
− | + | // график | |
− | + | var vGraph1 = new TM_graph( // определить график | |
− | + | "#vGraph1", // на html-элементе #vGraph | |
− | + | 20000, // сколько шагов по оси "x" отображается | |
− | + | 0, 1,0.1); // мин. значение оси Y, макс. значение оси Y, шаг по оси Y | |
− | + | // график | |
+ | var vGraph2 = new TM_graph( // определить график | ||
+ | "#vGraph2", // на html-элементе #vGraph | ||
+ | 20000, // сколько шагов по оси "x" отображается | ||
+ | 0, 19000,1000); // мин. значение оси Y, макс. значение оси Y, шаг по оси Y | ||
+ | |||
// Работа с массивом | // Работа с массивом | ||
Строка 177: | Строка 197: | ||
} | } | ||
} | } | ||
− | + | ||
var b = []; | var b = []; | ||
Строка 185: | Строка 205: | ||
// balls[balls.length] = b; // добавить элемент в конец массива | // balls[balls.length] = b; // добавить элемент в конец массива | ||
− | + | balls.push(b); | |
return b; | return b; | ||
}; | }; | ||
− | + | ||
− | + | this.setMy = function() { | |
− | + | balls = []; | |
for (var j = 0; j < Ny; j++) | for (var j = 0; j < Ny; j++) | ||
for (var i = 0; i < Ny; i++) | for (var i = 0; i < Ny; i++) | ||
− | + | addNewBall(w*j/Ny + r, h*i/Ny + r , true); // задаем 50x50 | |
− | + | document.getElementById('ballsNum').innerHTML = balls.length; | |
− | } | + | vGraph1.vArray = new Array(); |
− | + | vGraph2.vArray = new Array(); | |
+ | steps = 0; | ||
+ | G1 = 0; | ||
+ | G2 = 0; | ||
+ | } | ||
+ | |||
// Основной цикл программы | // Основной цикл программы | ||
Строка 209: | Строка 234: | ||
for (var s = 1; s <= spf; s++) { | for (var s = 1; s <= spf; s++) { | ||
− | + | L = 0; | |
+ | |||
+ | // пересчет сил идет отдельным массивом, т.к. далее будут добавляться силы взаимодействия между шарами | ||
+ | // for (var i0 = 0; i0 < balls.length; i0++) { | ||
+ | // balls[i0].vx += - B * balls[i0].vx*dt; | ||
+ | // balls[i0].vy += -B * balls[i0].vy*dt; | ||
+ | // } | ||
+ | for (var i0 = 0; i0 < balls.length; i0++) { | ||
Temp = 1/2 * (balls[i0].vx * balls[i0].vx + balls[i0].vy * balls[i0].vy) / Ny; | Temp = 1/2 * (balls[i0].vx * balls[i0].vx + balls[i0].vy * balls[i0].vy) / Ny; | ||
− | + | var stat = Math.sqrt(Temp0/Temp); | |
− | + | balls[i0].vx*=stat; | |
− | + | balls[i0].vy*=stat; | |
− | } | + | } |
+ | |||
for (var i = 0; i < balls.length; i++) | for (var i = 0; i < balls.length; i++) | ||
− | + | { | |
// расчет взаимодействия производится со всеми следующими шарами в массиве, | // расчет взаимодействия производится со всеми следующими шарами в массиве, | ||
// чтобы не считать каждое взаимодействие дважды | // чтобы не считать каждое взаимодействие дважды | ||
− | var b = balls[i]; | + | var b = balls[i]; |
+ | |||
for (var j = i + 1; j < balls.length; j++) { | for (var j = i + 1; j < balls.length; j++) { | ||
var b2 = balls[j]; | var b2 = balls[j]; | ||
var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) | var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b) | ||
− | var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами | + | var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами |
if (r2 > aCut2) continue; // проверка на радиус обрезания | if (r2 > aCut2) continue; // проверка на радиус обрезания | ||
var rLen = (Math.sqrt(r2)); | var rLen = (Math.sqrt(r2)); | ||
Строка 238: | Строка 272: | ||
// сила взаимодействия | // сила взаимодействия | ||
var s2 = a2 / r2; var s4 = s2 * s2; // ___в целях оптимизации___ | var s2 = a2 / r2; var s4 = s2 * s2; // ___в целях оптимизации___ | ||
− | + | var s8 = s4 * s4; | |
− | + | var F = LJCoeff * s8 * s2 * ((aCut - rLen)/(aCut - a)); // сила взаимодействия Леннарда-Джонса | |
− | + | ||
− | + | var vx21 = b.vx - b2.vx; var vy21 = b.vy - b2.vy; // вектор смотрит на первый шар (b) | |
− | + | var ex = rx / rLen; var ey = ry / rLen; | |
− | + | var v = vx21 * ex + vy21 * ey; | |
− | + | F -= betta * ((aCut - rLen)/(aCut - a)) * v * (1/rLen); | |
− | + | ||
// суммируем силы | // суммируем силы | ||
var Fx = F * rx*dt; var Fy = F * ry*dt; | var Fx = F * rx*dt; var Fy = F * ry*dt; | ||
− | + | ||
− | + | b.vx += Fx; b.vy += Fy; | |
− | + | b2.vx -= Fx; b2.vy -= Fy; | |
− | + | G2 += b.vx * b2.vx + b.vy * b2.vy; | |
− | + | G1 += 1; | |
− | + | ||
} | } | ||
Строка 263: | Строка 297: | ||
} | } | ||
− | + | for (var i0 = 0; i0 < balls.length; i0++) { | |
balls[i0].x += dt * balls[i0].vx; | balls[i0].x += dt * balls[i0].vx; | ||
balls[i0].y += dt * balls[i0].vy; | balls[i0].y += dt * balls[i0].vy; | ||
+ | L += Math.abs((balls[i0].x - w/2) * balls[i0].vy) + Math.abs((balls[i0].y - h/2) * balls[i0].vx); | ||
} | } | ||
− | + | ||
− | + | G = G2 / G1 / 2/ Temp0/Ny; | |
− | + | ||
− | + | steps++; | |
if (steps % 200 == 0) { | if (steps % 200 == 0) { | ||
− | + | vGraph1.graphIter(steps, G); | |
+ | vGraph2.graphIter(steps, L); | ||
+ | } | ||
} | } | ||
} | } | ||
Строка 298: | Строка 335: | ||
// Запуск системы | // Запуск системы | ||
this.newSystem(); | this.newSystem(); | ||
− | + | /*for (var i0 = 0; i0 < balls.length; i0++) { // задаем частицам случайные скорости | |
balls[i0].vx = (3 * (1 - 2 * Math.random())); | balls[i0].vx = (3 * (1 - 2 * Math.random())); | ||
balls[i0].vy = (3 * (1 - 2 * Math.random())); | balls[i0].vy = (3 * (1 - 2 * Math.random())); | ||
}*/ | }*/ | ||
− | + | if(!window.requestAnimationFrame){ | |
− | + | window.requestAnimationFrame = (function(){ | |
− | + | return window.webkitRequestAnimationFrame | |
− | + | || window.mozRequestAnimationFrame | |
− | + | || window.oRequestAnimationFrame | |
− | + | || window.msRequestAnimationFrame | |
− | + | || function(callback, element){window.setTimeout(callback, 1000 / fps); | |
− | + | }; | |
+ | })(); | ||
} | } | ||
//setInterval(control, 1000 / fps); | //setInterval(control, 1000 / fps); | ||
− | + | function animate(){ | |
− | + | requestAnimationFrame(animate); | |
− | + | control(); | |
− | + | } | |
− | + | animate(); | |
+ | // след. функция обновляет информацию о количестве частиц на поле | ||
+ | //setInterval(function(){document.getElementById('ballsNum').innerHTML = balls.length;}, 1000 / 20); | ||
+ | |||
} | } | ||
</syntaxhighlight> | </syntaxhighlight> | ||
Строка 323: | Строка 364: | ||
==Анализ== | ==Анализ== | ||
− | С помощью графика можно наблюдать за зависимостью корреляции скоростей <big><math> Г </math></big> от входных параметров (вязкость, температура, кол-во частиц). | + | С помощью графика можно наблюдать за зависимостью корреляции скоростей <big><math> Г </math></big> поведением кинетического моменты системы <big><math> L </math></big> от входных параметров (вязкость, температура, кол-во частиц). |
Уравнение для корреляции скоростей : | Уравнение для корреляции скоростей : | ||
Строка 334: | Строка 375: | ||
\end{cases} | \end{cases} | ||
</math></big> | </math></big> | ||
+ | |||
Например, для начальной конфигурации на 60 000 шаге по времени будет наблюдаться график 1. | Например, для начальной конфигурации на 60 000 шаге по времени будет наблюдаться график 1. | ||
Строка 366: | Строка 408: | ||
+ | |||
+ | |||
+ | |||
+ | Уравнение для кинетического момента системы : | ||
+ | |||
+ | <big><math>L = \sum_{i} \pmb r_i\times m_i \pmb V_i </math></big> | ||
==Ссылки== | ==Ссылки== | ||
*Vitaly A. Kuzkin, Anton M. Krivtsov. Interscale energy transport and velocity correlations in thermostated dissipative soft disc system. | *Vitaly A. Kuzkin, Anton M. Krivtsov. Interscale energy transport and velocity correlations in thermostated dissipative soft disc system. |
Текущая версия на 21:25, 19 декабря 2016
Курсовой проект по Механике дискретных сред
Исполнитель: Степанов Матвей
Группа: 09 (43604)
Семестр: осень 2016
Переход энергии с микро на макро уровень (и обратно) - одно из фундаментальных физических явлений. В данной работе,на примере двумерных турбулентных вихрей, рассматривается случай перехода с мелкомасштабного механического движения на крупномасштабное. В литературе это явление обычно упоминается как “обратный каскад” энергии. Построенная модель позволяет исследовать переход энергии с микро на макро уровень и корреляции скоростей в дискретной среде.
Модель представлена набором взаимодействующих частиц с случайными начальными скоростями. Частицы взаимодействуют за счет отталкивающих потенциальных и диссипативных сил. Динамика взаимодействия описана набором уравнений движения Ньютона:
где m, r - масса и радиус вектор i-ой частицы,
- отталкивающая потенциальная и диссипативная силы соответственно.Выражения для потенциальной и диссипативной сил:
где
- радиус обрезания, - энергетические параметры системы, - коэффициент вязкости.Диссипативные силы уменьшают полную энергию системы. Таким образом, без внешнего добавления энергии система перейдет в состояние равновесия. В данной модели энергия добавляется при помощи термостата Берендсена. Скорости частиц на каждом шаге по времени перемножаются на параметр
:
где
- температура средыПрограмма[править]
Частички расположены впритык друг к другу и образуют прямоугольную решетку. С помощью написанной программы можно изменять количество шаров, коэффициент вязкости и температуру в среде(термостат).
Файл "Inverse cascade.js"
1 // Версия 16.12.16
2
3
4 // Версия Total 19.12.16
5 //Stepanov MD
6
7 function MainBalls(canvas, slider_01, text_01, slider_02, text_02) {
8
9 canvas.onselectstart = function () {return false;}; // запрет выделения canvas
10
11 // Предварительные установки
12
13 var context = canvas.getContext("2d"); // на context происходит рисование
14 canvas.oncontextmenu = function (e) {return false;}; // блокировка контекстного меню
15
16 var Pi = 3.1415926; // число "пи"
17
18 var m0 = 1; // масштаб массы
19 var T0 = 1; // масштаб времени (период колебаний исходной системы)
20 var a0 = 1; // масштаб расстояния (диаметр шара)
21
22 var g0 = a0 / T0 / T0; // масштаб ускорения (ускорение, при котором за T0 будет пройдено расстояние a0)
23 var k0 = 2 * Pi / T0; // масштаб частоты
24 var C0 = m0 * k0 * k0; // масштаб жесткости
25 var B0 = 2 * m0 * k0; // масштаб вязкости
26
27 var L = 0;
28 var G = 0;
29 var G1 = 0;
30 var G2 = 0;
31
32 // *** Задание физических параметров ***
33
34 var Ny = 20; // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)
35 var m = 1 * m0; // масса
36 var Cwall = 10 * C0; // жесткость стен
37 var Cball = 0.1 * Cwall; // жесткость между частицами
38
39 /*var B = 0.000 * B0; // вязкость среды*/
40
41 var Bball = 0.01 * B0; // вязкость между частицами
42 var Bwall = 0.03 * B0; // вязкость на стенках
43 var mg = 0.25 * m * g0; // сила тяжести
44 var r = 0.5 * a0; // радиус частицы в расчетных координатах
45 var K = 0.7; // все силы, зависящие от радиуса, ограничиваются значением, реализующимся при r/a = K
46 var a = 2 * r; // равновесное расстояние между частицами
47 var aCut = 2 * a; // радиус обрезания
48
49 // *** Задание вычислительных параметров ***
50
51 var fps = 50; // frames per second - число кадров в секунду (качечтво отображения)
52 var spf = 100; // steps per frame - число шагов интегрирования между кадрами (скорость расчета)
53 var dt = 0.045 * T0 / fps; // шаг интегрирования (качество расчета)
54
55
56 var steps = 0;
57
58
59 // Выполнение программы
60
61 var r2 = r * r; // ___в целях оптимизации___
62 var aCut2 = aCut * aCut; // ___в целях оптимизации___
63 var a2 = a * a; // ___в целях оптимизации___
64 var D = a2 * Cball / 72; // энергия связи между частицами
65 var Temp = 0;
66 Temp0 = 0.12 * D;
67 var B1 = 12*Math.sqrt((2*D)/a2);
68 var B = 0.026 * B1;
69 var betta = 3.5*B;
70
71 var LJCoeff = 6 * D / a2; // коэффициент для расчета потенциала Л-Дж
72 var Ka = K * a; // ___в целях оптимизации___
73 var K2a2 = K * K * a2; // ___в целях оптимизации___
74
75 var dNd = null; // ссылка на захваченный курсором шар (drag & drop)
76 var grad; // должен ли работать градиент (регулируется в функции setNy())
77
78 this.set_01 = function(p) {betta = p * B;};
79 this.set_02 = function(pT) {Temp0 = pT *D;};
80
81 this.setNy = function(ny) {
82 Ny = ny;
83 context.fillStyle = "#3070d0"; // цвет, шара
84 };
85 this.setNy(Ny); // запускаем с уже присвоенным значением, чтобы обновились настройки градиента
86
87 // Запуск новой системы
88
89 // следующие переменные должны пересчитываться каждый раз, когда мы изменяем значение Ny
90 var scale, w, h;
91 var rScale13, rScaleShift;
92 this.newSystem = function() {
93 scale = canvas.height / Ny / a0; // масштабный коэффициент для перехода от расчетных к экранным координатам
94 w = canvas.width / scale; // ширина окна в расчетных координатах
95 h = canvas.height / scale; // высота окна в расчетных координатах
96
97 rScale13 = r * scale * 1.3; // ___в целях оптимизации___
98 rScaleShift = r * scale / 5; // ___в целях оптимизации___
99
100 this.setMy();
101 };
102
103 // настройка слайдеров и текстовых полей
104 slider_01.min = 0.1; slider_01.max = 20;
105 slider_01.step = 0.1;
106 slider_01.value = betta / B; // начальное значение ползунка должно задаваться после min и max
107 text_01.value = betta / B;
108
109 // настройка слайдеров и текстовых полей
110 slider_02.min = 0.1; slider_02.max = 1;
111 slider_02.step = 0.1;
112 slider_02.value = Temp0 / D; // начальное значение ползунка должно задаваться после min и max
113 text_02.value = Temp0 / D;
114
115 // график
116 var vGraph1 = new TM_graph( // определить график
117 "#vGraph1", // на html-элементе #vGraph
118 20000, // сколько шагов по оси "x" отображается
119 0, 1,0.1); // мин. значение оси Y, макс. значение оси Y, шаг по оси Y
120 // график
121 var vGraph2 = new TM_graph( // определить график
122 "#vGraph2", // на html-элементе #vGraph
123 20000, // сколько шагов по оси "x" отображается
124 0, 19000,1000); // мин. значение оси Y, макс. значение оси Y, шаг по оси Y
125
126 // Работа с массивом
127
128 var balls = []; // массив шаров
129 var addNewBall = function(x, y, check) {
130 // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
131 if (check) {
132 if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
133 for (var i = 0; i < balls.length; i++) {
134 var rx = balls[i].x - x;
135 var ry = balls[i].y - y;
136 var rLen2 = rx * rx + ry * ry;
137 if (rLen2 < 4 * r2) return null;
138 }
139 }
140
141 var b = [];
142
143 b.x = x; b.y = y; // расчетные координаты шара
144 b.fx = 0; b.fy = 0; // сила, действующая на шар
145 b.vx = (3 * (1 - 2 * Math.random())); b.vy = (3 * (1 - 2 * Math.random())); // скорость
146
147 // balls[balls.length] = b; // добавить элемент в конец массива
148 balls.push(b);
149 return b;
150 };
151
152 this.setMy = function() {
153 balls = [];
154 for (var j = 0; j < Ny; j++)
155 for (var i = 0; i < Ny; i++)
156 addNewBall(w*j/Ny + r, h*i/Ny + r , true); // задаем 50x50
157 document.getElementById('ballsNum').innerHTML = balls.length;
158 vGraph1.vArray = new Array();
159 vGraph2.vArray = new Array();
160 steps = 0;
161 G1 = 0;
162 G2 = 0;
163 }
164
165 // Основной цикл программы
166
167 function control() {
168 physics();
169 draw();
170 }
171
172 // Расчетная часть программы
173
174 function physics() { // то, что происходит каждый шаг времени
175 for (var s = 1; s <= spf; s++) {
176
177 L = 0;
178
179 // пересчет сил идет отдельным массивом, т.к. далее будут добавляться силы взаимодействия между шарами
180 // for (var i0 = 0; i0 < balls.length; i0++) {
181 // balls[i0].vx += - B * balls[i0].vx*dt;
182 // balls[i0].vy += -B * balls[i0].vy*dt;
183 // }
184 for (var i0 = 0; i0 < balls.length; i0++) {
185 Temp = 1/2 * (balls[i0].vx * balls[i0].vx + balls[i0].vy * balls[i0].vy) / Ny;
186 var stat = Math.sqrt(Temp0/Temp);
187 balls[i0].vx*=stat;
188 balls[i0].vy*=stat;
189 }
190
191 for (var i = 0; i < balls.length; i++)
192 {
193 // расчет взаимодействия производится со всеми следующими шарами в массиве,
194 // чтобы не считать каждое взаимодействие дважды
195 var b = balls[i];
196
197 for (var j = i + 1; j < balls.length; j++) {
198 var b2 = balls[j];
199 var rx = b.x - b2.x; var ry = b.y - b2.y; // вектор смотрит на первый шар (b)
200 var r2 = rx * rx + ry * ry; // квадрат расстояния между шарами
201 if (r2 > aCut2) continue; // проверка на радиус обрезания
202 var rLen = (Math.sqrt(r2));
203
204 // если расстояние между частицами мало, силы будут посчитаны для K * a
205 if (r2 < K2a2) {
206 if (rLen > 0.00001) { // проверка, чтобы избежать деления на 0
207 rx = rx / rLen * Ka;
208 ry = ry / rLen * Ka;
209 }
210 r2 = K2a2;
211 rLen = Ka; // корень K2a2
212 }
213 // сила взаимодействия
214 var s2 = a2 / r2; var s4 = s2 * s2; // ___в целях оптимизации___
215 var s8 = s4 * s4;
216 var F = LJCoeff * s8 * s2 * ((aCut - rLen)/(aCut - a)); // сила взаимодействия Леннарда-Джонса
217
218 var vx21 = b.vx - b2.vx; var vy21 = b.vy - b2.vy; // вектор смотрит на первый шар (b)
219 var ex = rx / rLen; var ey = ry / rLen;
220 var v = vx21 * ex + vy21 * ey;
221 F -= betta * ((aCut - rLen)/(aCut - a)) * v * (1/rLen);
222
223 // суммируем силы
224 var Fx = F * rx*dt; var Fy = F * ry*dt;
225
226 b.vx += Fx; b.vy += Fy;
227 b2.vx -= Fx; b2.vy -= Fy;
228
229 G2 += b.vx * b2.vx + b.vy * b2.vy;
230 G1 += 1;
231
232 }
233
234 if (b.y + r > h) { b.vy += (-Cwall * (b.y + r - h) - Bwall * b.vy)*dt; }
235 if (b.y - r < 0) { b.vy += -Cwall * (b.y - r) - Bwall * b.vy;}
236 if (b.x + r > w) { b.vx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
237 if (b.x - r < 0) { b.vx += -Cwall * (b.x - r) - Bwall * b.vx; }
238
239 }
240 for (var i0 = 0; i0 < balls.length; i0++) {
241 balls[i0].x += dt * balls[i0].vx;
242 balls[i0].y += dt * balls[i0].vy;
243 L += Math.abs((balls[i0].x - w/2) * balls[i0].vy) + Math.abs((balls[i0].y - h/2) * balls[i0].vx);
244 }
245
246 G = G2 / G1 / 2/ Temp0/Ny;
247
248 steps++;
249 if (steps % 200 == 0) {
250 vGraph1.graphIter(steps, G);
251 vGraph2.graphIter(steps, L);
252 }
253 }
254 }
255
256 // Рисование
257 function draw() {
258 context.clearRect(0, 0, w * scale, h * scale); // очистить экран
259 for (var i = 0; i < balls.length; i++){
260 var xS = balls[i].x * scale; var yS = balls[i].y * scale;
261 if (grad) {
262 // расчет градиента нужно проводить для каждого шара
263 var gradient = context.createRadialGradient(xS, yS, rScale13, xS - rScaleShift, yS + rScaleShift, 0);
264 gradient.addColorStop(0, "#0000bb");
265 gradient.addColorStop(1, "#44ddff");
266 context.fillStyle = gradient;
267 }
268
269 context.beginPath();
270 context.arc(xS, yS, r*scale, 0, 2 * Math.PI, false);
271 context.closePath();
272 context.fill();
273 }
274 }
275
276 // Запуск системы
277 this.newSystem();
278 /*for (var i0 = 0; i0 < balls.length; i0++) { // задаем частицам случайные скорости
279 balls[i0].vx = (3 * (1 - 2 * Math.random()));
280 balls[i0].vy = (3 * (1 - 2 * Math.random()));
281 }*/
282 if(!window.requestAnimationFrame){
283 window.requestAnimationFrame = (function(){
284 return window.webkitRequestAnimationFrame
285 || window.mozRequestAnimationFrame
286 || window.oRequestAnimationFrame
287 || window.msRequestAnimationFrame
288 || function(callback, element){window.setTimeout(callback, 1000 / fps);
289 };
290 })();
291 }
292 //setInterval(control, 1000 / fps);
293
294 function animate(){
295 requestAnimationFrame(animate);
296 control();
297 }
298 animate();
299 // след. функция обновляет информацию о количестве частиц на поле
300 //setInterval(function(){document.getElementById('ballsNum').innerHTML = balls.length;}, 1000 / 20);
301
302 }
Анализ[править]
С помощью графика можно наблюдать за зависимостью корреляции скоростей
поведением кинетического моменты системы от входных параметров (вязкость, температура, кол-во частиц).Уравнение для корреляции скоростей :
Например, для начальной конфигурации на 60 000 шаге по времени будет наблюдаться график 1.
А для измененных входных параметров график 2 и 3.
Уравнение для кинетического момента системы :
Ссылки[править]
- Vitaly A. Kuzkin, Anton M. Krivtsov. Interscale energy transport and velocity correlations in thermostated dissipative soft disc system.