Обратный каскад энергии (двумерная турбулентность)

Материал из Department of Theoretical and Applied Mechanics
Перейти к: навигация, поиск


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

Исполнитель: Степанов Матвей

Группа: 09 (43604)

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


Рис.1 Взаимодействие двух вихрей.

Переход энергии с микро на макро уровень (и обратно) - одно из фундаментальных физических явлений. В данной работе,на примере двумерных турбулентных вихрей, рассматривается случай перехода с мелкомасштабного механического движения на крупномасштабное. В литературе это явление обычно упоминается как “обратный каскад” энергии. Построенная модель позволяет исследовать переход энергии с микро на макро уровень и корреляции скоростей в дискретной среде.

Модель представлена набором взаимодействующих частиц с случайными начальными скоростями. Частицы взаимодействуют за счет отталкивающих потенциальных и диссипативных сил. Динамика взаимодействия описана набором уравнений движения Ньютона:

[math]m\pmb{\ddot r}= \sum_{j}\left({{}^p}{\pmb F_{ij}} + {{}^d}{\pmb F_{ij}}\right),\qquad i,j = 1...N[/math]

где m, r - масса и радиус вектор i-ой частицы, [math]{{}^p}{\pmb F_{ij}}, {{}^d}{\pmb F_{ij}}[/math] - отталкивающая потенциальная и диссипативная силы соответственно.

Выражения для потенциальной и диссипативной сил:

[math] {{}^p}{\pmb F_{ij}} = \begin{cases} \frac{6D}{a^{2}}{\left(\frac{a_{c} - \pmb r_{ij}}{a_{c} - a}\right)}^{2}{\left(\frac{a}{\pmb r_{ij}}\right)}^{10}, & \text{если} \; \pmb r_{ij} \leqslant a_{c}; \\ 0, & \text{если} \; \pmb r_{ij} \gt a_{c}. \end{cases} [/math]

[math] {{}^d}{\pmb F_{ij}} = \begin{cases} \beta {\left(\frac{a_{c} - \pmb r_{ij}}{a_{c} - a}\right)}^{2} \left( \pmb v_{j} - \pmb v_{i}\right)\cdot \pmb e_{ij} \pmb e_{ij} , & \text{если} \; \pmb r_{ij} \leqslant a_{c}; \\ 0, & \text{если} \; \pmb r_{ij} \gt a_{c}. \end{cases} [/math]

где [math]\pmb r_{ij} =\pmb r_{i} - \pmb r_{j},\; \pmb e_{ij} = \frac{\pmb r_{ij}}{r_{ij}},\; a_{c}[/math] - радиус обрезания, [math]a, D[/math] - энергетические параметры системы, [math] \beta [/math] - коэффициент вязкости.

Диссипативные силы уменьшают полную энергию системы. Таким образом, без внешнего добавления энергии система перейдет в состояние равновесия. В данной модели энергия добавляется при помощи термостата Берендсена. Скорости частиц на каждом шаге по времени перемножаются на параметр [math] s [/math] :

[math] s = \sqrt{\frac{T_{0}}{T}},\qquad T =\frac{1}{2}\sum_{i} m {\pmb v_{i}}^{2} [/math]

где [math] T [/math] - температура среды

Программа

Частички расположены впритык друг к другу и образуют прямоугольную решетку. С помощью написанной программы можно изменять количество шаров, коэффициент вязкости и температуру в среде(термостат).

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

Файл "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 }

Анализ

С помощью графика можно наблюдать за зависимостью корреляции скоростей [math] Г [/math] поведением кинетического моменты системы [math] L [/math] от входных параметров (вязкость, температура, кол-во частиц).

Уравнение для корреляции скоростей :

[math] Г= \frac{\sum_{i,j\neq i} m \pmb v_{i} \cdot \pmb v_{j} w\left(\pmb r_{ij}\right)}{2\sum_{i,j\neq i} w\left(\pmb r_{ij}\right) },\qquad w\left(r\right)= \begin{cases} 1, & \text{если} \; r \leqslant a_{c}; \\ 0, & \text{если} \; r \gt a_{c}. \end{cases} [/math]


Например, для начальной конфигурации на 60 000 шаге по времени будет наблюдаться график 1.

График 1. Начальная конфигурация:betta = 3.5 B; Temp0 = 0.12 D; Количество частиц:900 .

А для измененных входных параметров график 2 и 3.

График 2. Конфигурация:betta = 5 B; Temp0 = 0.12 D; Количество частиц:900 .
График 3. Конфигурация:betta = 1.5 B; Temp0 = 0.12 D; Количество частиц:900 .















Уравнение для кинетического момента системы :

[math]L = \sum_{i} \pmb r_i\times m_i \pmb V_i [/math]

Ссылки

  • Vitaly A. Kuzkin, Anton M. Krivtsov. Interscale energy transport and velocity correlations in thermostated dissipative soft disc system.