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




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

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

'''Группа:''' [[Группа 09|09]] (43604)

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



[[Файл:ObratnyKasakad2vyhrya.jpg|thumb|Рис.1 Взаимодействие двух вихрей.|450px]]

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

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

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

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

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

<big><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} > a_{c}.
\end{cases}
</math></big>

<big><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} > a_{c}.
\end{cases}
</math></big>

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

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

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

где <big><math> T </math></big> - температура среды

==Программа==

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

{{#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">
'''Текст программы на языке JavaScript:''' <div class="mw-collapsible-content">
Файл '''"Inverse cascade.js"'''
<syntaxhighlight lang="javascript" line start="1" enclose="div">

// Версия 16.12.16


// Версия Total 19.12.16
//Stepanov MD

function MainBalls(canvas, slider_01, text_01, slider_02, text_02) {

    canvas.onselectstart = function () {return false;};     // запрет выделения canvas

    // Предварительные установки

    var context = canvas.getContext("2d");                  // на context происходит рисование
    canvas.oncontextmenu = function (e) {return false;};    // блокировка контекстного меню

    var Pi = 3.1415926;                 // число "пи"

    var m0 = 1;                         // масштаб массы
    var T0 = 1;                         // масштаб времени (период колебаний исходной системы)
    var a0 = 1;                         // масштаб расстояния (диаметр шара)

    var g0 = a0 / T0 / T0;              // масштаб ускорения (ускорение, при котором за T0 будет пройдено расстояние a0)
    var k0 = 2 * Pi / T0;               // масштаб частоты
    var C0 = m0 * k0 * k0;              // масштаб жесткости
    var B0 = 2 * m0 * k0;               // масштаб вязкости
  
  var L = 0;
  var G = 0;
  var G1 = 0;
  var G2 = 0;

    // *** Задание физических параметров ***

    var Ny = 20;                         // число шаров, помещающихся по вертикали в окно (задает размер шара относительно размера окна)
    var m = 1 * m0;                     // масса
    var Cwall = 10 * C0;                // жесткость стен
    var Cball = 0.1 * Cwall;            // жесткость между частицами

    /*var B = 0.000 * B0;                 // вязкость среды*/
 
    var Bball = 0.01 * B0;              // вязкость между частицами
    var Bwall = 0.03 * B0;              // вязкость на стенках
    var mg = 0.25 * m * g0;             // сила тяжести
    var r = 0.5 * a0;                   // радиус частицы в расчетных координатах
    var K = 0.7;                        // все силы, зависящие от радиуса, ограничиваются значением, реализующимся при r/a = K
    var a = 2 * r;                      // равновесное расстояние между частицами
    var aCut = 2 * a;                   // радиус обрезания

    // *** Задание вычислительных параметров ***

    var fps = 50;                       // frames per second - число кадров в секунду (качечтво отображения)
    var spf = 100;                      // steps per frame   - число шагов интегрирования между кадрами (скорость расчета)
    var dt  = 0.045 * T0 / fps;         // шаг интегрирования (качество расчета)
                        
                        
                        var steps = 0;
                        
                        
    // Выполнение программы

    var r2 = r * r;                     // ___в целях оптимизации___
    var aCut2 = aCut * aCut;            // ___в целях оптимизации___
    var a2 = a * a;                     // ___в целях оптимизации___
    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 K2a2 = K * K * a2;              // ___в целях оптимизации___

    var dNd = null;                     // ссылка на захваченный курсором шар (drag & drop)
    var grad;                           // должен ли работать градиент (регулируется в функции setNy())

  this.set_01 = function(p) {betta = p * B;};
  this.set_02 = function(pT) {Temp0 = pT *D;};
  
    this.setNy = function(ny) {
        Ny = ny;
        context.fillStyle = "#3070d0";  // цвет, шара
    };
    this.setNy(Ny);                         // запускаем с уже присвоенным значением, чтобы обновились настройки градиента

    // Запуск новой системы

    // следующие переменные должны пересчитываться каждый раз, когда мы изменяем значение Ny
    var scale, w, h;
    var rScale13, rScaleShift;
    this.newSystem = function() {
        scale = canvas.height / Ny / a0;    // масштабный коэффициент для перехода от расчетных к экранным координатам
        w = canvas.width / scale;           // ширина окна в расчетных координатах
        h = canvas.height / scale;          // высота окна в расчетных координатах

        rScale13 = r * scale * 1.3;         // ___в целях оптимизации___
        rScaleShift = r * scale / 5;        // ___в целях оптимизации___

      this.setMy();
    };

   // настройка слайдеров и текстовых полей
    slider_01.min = 0.1;               slider_01.max = 20;
    slider_01.step = 0.1;
    slider_01.value = betta / B;          // начальное значение ползунка должно задаваться после min и max
    text_01.value = betta / B;

// настройка слайдеров и текстовых полей
    slider_02.min = 0.1;               slider_02.max = 1;
    slider_02.step = 0.1;
    slider_02.value = Temp0 / D;          // начальное значение ползунка должно задаваться после min и max
    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
  
    // Работа с массивом

    var balls = [];                             // массив шаров
    var addNewBall =  function(x, y, check) {
        // проверка - не пересекается ли новый шар со стенами или уже существующими шарами
        if (check) {
            if (x - r < 0 || x + r > w || y - r < 0 || y + r > h) return null;
            for (var i = 0; i < balls.length; i++) {
                var rx = balls[i].x - x;
                var ry = balls[i].y - y;
                var rLen2 = rx * rx + ry * ry;
                if (rLen2 < 4 * r2) return null;
            }
        }
    
        var b = [];

        b.x = x;                b.y = y;        // расчетные координаты шара
        b.fx = 0;               b.fy = 0;      // сила, действующая на шар
        b.vx = (3 * (1 - 2 * Math.random()));               b.vy = (3 * (1 - 2 * Math.random()));       // скорость

       // balls[balls.length] = b;                // добавить элемент в конец массива
    balls.push(b);
        return b;
    };
  
  this.setMy = function() {
    balls = [];
        for (var j = 0; j < Ny; j++)
            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;
  }
  
    // Основной цикл программы

    function control() {
        physics();
        draw();
    }

    // Расчетная часть программы

    function physics() {                        // то, что происходит каждый шаг времени
        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;
        var stat = Math.sqrt(Temp0/Temp);
        balls[i0].vx*=stat;
        balls[i0].vy*=stat;
            }
      
            for (var i = 0; i < balls.length; i++) 
      {
                // расчет взаимодействия производится со всеми следующими шарами в массиве,
                // чтобы не считать каждое взаимодействие дважды
                var b = balls[i];          
        
                for (var j = i + 1; j < balls.length; j++) {
                    var b2 = balls[j];
                    var rx = b.x - b2.x;   var ry = b.y - b2.y;         // вектор смотрит на первый шар (b)
                    var r2 = rx * rx + ry * ry;                         // квадрат расстояния между шарами          
                    if (r2 > aCut2) continue;                           // проверка на радиус обрезания
                    var rLen = (Math.sqrt(r2));

                    // если расстояние между частицами мало, силы будут посчитаны для K * a
                    if (r2 < K2a2) {
                        if (rLen > 0.00001) {                           // проверка, чтобы избежать деления на 0
                            rx = rx / rLen * Ka;
                            ry = ry / rLen * Ka;
                        }
                        r2 = K2a2;
                        rLen = Ka;                                      // корень K2a2
                    }
                    // сила взаимодействия
                    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;
          
          b.vx += Fx;        b.vy += Fy;
          b2.vx -= Fx;        b2.vy -= Fy;

          G2 += b.vx * b2.vx + b.vy * b2.vy;
          G1 += 1;
          
                }

                if (b.y + r > h) { b.vy += (-Cwall * (b.y + r - h) - Bwall * b.vy)*dt; }
                if (b.y - r < 0) { b.vy += -Cwall * (b.y - r) - Bwall * b.vy;}
                if (b.x + r > w) { b.vx += -Cwall * (b.x + r - w) - Bwall * b.vx; }
                if (b.x - r < 0) { b.vx += -Cwall * (b.x - r) - Bwall * b.vx; }

            }  
      for (var i0 = 0; i0 < balls.length; i0++) {
                balls[i0].x += dt * balls[i0].vx;
                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) {
        vGraph1.graphIter(steps, G);
		vGraph2.graphIter(steps, L);
      }
        }
    }

    // Рисование
    function draw() {
        context.clearRect(0, 0, w * scale, h * scale);      // очистить экран
        for (var i = 0; i < balls.length; i++){
            var xS = balls[i].x * scale;           var yS = balls[i].y * scale;
            if (grad) {
                // расчет градиента нужно проводить для каждого шара
                var gradient = context.createRadialGradient(xS, yS, rScale13, xS - rScaleShift, yS + rScaleShift, 0);
                gradient.addColorStop(0, "#0000bb");
                gradient.addColorStop(1, "#44ddff");
                context.fillStyle = gradient;
            }

            context.beginPath();
            context.arc(xS, yS, r*scale, 0, 2 * Math.PI, false);
            context.closePath();
            context.fill();
        }
    }

    // Запуск системы
    this.newSystem();
  /*for (var i0 = 0; i0 < balls.length; i0++) {     // задаем частицам случайные скорости
        balls[i0].vx = (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);

  function animate(){
    requestAnimationFrame(animate);
    control();
  }
  animate();
    // след. функция обновляет информацию о количестве частиц на поле
   //setInterval(function(){document.getElementById('ballsNum').innerHTML = balls.length;}, 1000 / 20);

}
</syntaxhighlight>
</div>

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

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

<big><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 > a_{c}.
\end{cases}
</math></big>

Например, для начальной конфигурации на 60 000 шаге по времени будет наблюдаться график 1.
[[Файл:1cascade.png|thumb|График 1. Начальная конфигурация:'''betta''' = 3.5 B; '''Temp0''' = 0.12 D; '''Количество частиц''':900 .|центр|450px]]
А для измененных входных параметров график 2 и 3.
[[Файл:3cascade.png|thumb|График 2. Конфигурация:'''betta''' = 5 B; '''Temp0''' = 0.12 D; '''Количество частиц''':900 .|слева|450px]]
[[Файл:2cascade.png|thumb|График 3. Конфигурация:'''betta''' = 1.5 B; '''Temp0''' = 0.12 D; '''Количество частиц''':900 .|справа|450px]]



























==Ссылки==
*Vitaly A. Kuzkin, Anton M. Krivtsov. Interscale energy transport and velocity correlations in thermostated dissipative soft disc system.