Двумерное уравнение теплопроводности + волновое уравнение — различия между версиями

  1 window.addEventListener("load", main_equations, false);
  2 function main_equations() {
  3 
  4     // Предварительные установки
  5     var canv = init_canvas(Wave_Heat_2D_canvas);
  6 
  7     var a0 = 1;                             // масштаб расстояния
  8     var t0 = 1;                             // масштаб времени
  9     var m0 = 1;                             // масштаб массы
 10     var T0 = 1;                             // масштаб температуры
 11 
 12     var dx = 1 * a0;                        // шаг сетки по оси x
 13     var dy = 1 * a0;                        // шаг сетки по оси y
 14 
 15     var spf = 2;                            // steps per frame - сколько расчетов проходит за каждый кадр отображения
 16     var fps = 15;                           // frames per second - количество кадров в секунду
 17     var dt = 0.005 * t0;                    // шаг интегрирования по времени
 18     var t = 0;
 19 
 20     var p0 = m0 / (a0 * a0 * a0);           // единица плотности, кг/м3
 21     var c0 = a0 * a0 / (t0 * t0 * T0);      // единица удельной теплоемкости, Дж/(кг * К)    (Дж = кг * м2 / с2)
 22     var kap0 = m0 * a0 / (t0 * t0 * t0 * T0);// единица теплопроводности (каппа)
 23 
 24     var p = 0.15 * p0;
 25     var c = 0.2 * c0;
 26     var k = 1 * kap0;
 27     var X = k / (c * p);                    // температуропроводность (хи)
 28 
 29     var Nx = 50 + 2;                        // количество узлов по оси x + 2 для ГУ
 30     var Ny = 50 + 2;                        // количество узлов по оси y + 2 для ГУ
 31 
 32     var cell_w = canv.w / (Nx - 2);         // ширина клетки
 33     var cell_h = canv.h / (Ny - 2);         // высота клетки
 34 
 35     var T_max = 10 * T0;
 36     var T_min = - T_max;
 37 
 38     var k0 = 2 * Math.PI / t0;              // масштаб частоты
 39     var C0 = m0 * k0 * k0;                  // масштаб жесткости
 40 
 41     var m = 1 * m0;                 	    // масса
 42     var C = 1 * C0;                 	    // жесткость
 43 
 44     var omega = Math.sqrt(C / m);
 45     var cc = 2.3 * a0 * omega;
 46 
 47     var pause = false;
 48     button_pause.onclick = function () {pause = !pause;};
 49     button_clear.onclick = clear_button_func;
 50     function clear_button_func() {
 51         clear_all();
 52         draw();
 53     }
 54     function clear_all() {
 55         for (var i = 1; i < Nx - 1; i++)
 56         for (var j = 1; j < Ny - 1; j++) {
 57                 T[i][j].u = 0;
 58                 T[i][j].v = 0;
 59                 T[i][j].BC = false;
 60         }
 61         BC2 = [];
 62     }
 63 
 64     var koeff1;
 65 
 66     // начальные условия - распределение Гаусса
 67     var T = [];
 68     for (var i = 0; i < Nx; i++) {
 69         T[i] = [];
 70         for (var j = 0; j < Ny; j++) {
 71             T[i][j] = {};
 72             var x = i / (Nx - 1);
 73             var y = j / (Ny - 1);
 74 
 75             var width = 0.05;
 76             var power = 15;
 77             var minus = 0;
 78 
 79             T[i][j].u = Math.exp(-Math.pow(x - 0.5, 2) / width) * Math.exp(-Math.pow(y - 0.5, 2) / width) * power - minus;
 80             T[i][j].v = 0;
 81             T[i][j].BC = false;
 82         }
 83     }
 84 
 85     function set_calc_func() {
 86         if (radio_heat.checked) {
 87             heat_equation = true;
 88             slider_damping_power.disabled = true;
 89             damping_span.style.color="#888888";
 90             koeff1 = X / (dx * dy);
 91         } else {
 92             heat_equation = false;
 93             slider_damping_power.disabled = false;
 94             damping_span.style.color="#000000";
 95             koeff1 = cc * cc / (dx * dy);
 96         }
 97     }
 98     var heat_equation;
 99     set_calc_func();
100     radio_heat.onchange = set_calc_func;
101     radio_wave.onchange = set_calc_func;
102 
103     function set_BC_func() {
104         if (radio_BC_mirror.checked) BC_mirror();
105         else if (radio_BC_periodic.checked) BC_periodic();
106     }
107     radio_BC_mirror.onclick = set_BC_func;
108     radio_BC_periodic.onclick = set_BC_func;
109 
110     // функции и переменные, объявленные через this, можно будет вызывать по сокращенному названию (например, sin(x), PI)
111     var math_methods = Object.getOwnPropertyNames(Math);
112     for (var i in math_methods)
113         this[math_methods[i]] = Math[math_methods[i]];
114 
115     var damping_power = 0;
116     slider_damping_power.min = 0;
117     slider_damping_power.max = 5;
118     slider_damping_power.step = 5 / 100;
119     slider_damping_power.value = damping_power;
120     slider_damping_power.oninput = function() {damping_power = slider_damping_power.value};
121 
122     apply_IC_button.onclick = apply_IC;
123     function apply_IC() {
124         var script = document.createElement('script');
125         script.innerHTML = "function IC_function(nx, ny, M, M_max, M_min, t, add_BC2) {" +
126             "for (var i = 0; i < nx; i++) {\r\n" +
127             "    for (var j = 0; j < ny; j++) {\r\n" +
128             IC_text.value + "}}};";
129         document.getElementsByTagName('head')[0].appendChild(script);
130         clear_all();
131         IC_function(nx, ny, M, T_max, T_min, t, add_BC2);
132         draw();
133     }
134     clear_IC_window_button.onclick =  clear_IC_window;
135     function clear_IC_window() {
136         IC_text.value = "";
137 //        IC_text.value =
138 //            "for (var i = 0; i < nx; i++) {\r\n" +
139 //            "    for (var j = 0; j < ny; j++) {\r\n" +
140 //            "       \r\n" +
141 //            "    }\r\n" +
142 //            "}";
143     }
144     var functions_for_buttons = [
145         "M[i][j].u = M_max * (sin(i / (nx - 1) * 2 * PI));\r\n"
146        ,"M[i][j].v = M_max * (sin(i / (nx - 1) * 2 * PI));\r\n"
147        ,"M[i][j].u = M_max * pow(sin(i / (nx - 1) * 2 * PI), 25);\r\n"
148        ,"M[i][j].u = M_max * (sin((i + j) / (2*(nx - 1)) * 2 * PI));\r\n"
149        ,"var lx = 0.5 - i / (nx - 1);\r\n" +
150         "var ly = 0.5 - j / (ny - 1);\r\n" +
151         "M[i][j].u = sin(sqrt(lx * lx + ly * ly) * 12 * PI) * 6;\r\n"
152        ,"var lx = 0.5 - i / (nx - 1);\r\n" +
153         "var ly = 0.5 - j / (ny - 1);\r\n" +
154         "M[i][j].u = 1 / max(sqrt(lx * lx + ly * ly), 0.01);    // max - чтобы избежать деления на 0\r\n"
155        ,"M[i][j].u = M_max * sin(i / (nx - 1) * 4 * PI) * sin(j / (ny - 1) * 4 * PI);\r\n"
156        ,"M[i][j].u = M_max * sin(i / (nx - 1) * nx * PI) * sin(j / (ny - 1) * ny * PI);\r\n"
157        ,"M[i][j].u = M_max * (min(sin(i / (nx - 1) * 9 * PI), sin(j / (ny - 1) * 9 * PI)));\r\n"
158        ,"M[i][j].u = M_max * (2 * random() - 1);\r\n"
159        ,"M[i][j].u = M_max * pow(2 * random() - 1, 201);\r\n"
160        ,"M[i][j].u = - M_max * (sin(i / (nx - 1) * PI) * random());\r\n"
161        ,"M[i][j].u = 2 * M_max * (i / (nx - 1)) - M_max;\r\n"
162        ,"M[i][j].u = i / (nx - 1) * j / (ny - 1) * M_max;\r\n"
163        ,"if (i == 0) {\r\n" +
164         "    M[i][j].u = M_max;\r\n" +
165         "    M[i][j].BC = true;\r\n" +
166         "} else if (i == nx - 1) {\r\n" +
167         "    M[i][j].u = M_min;\r\n" +
168         "    M[i][j].BC = true;\r\n" +
169         "}\r\n"
170        ,"if (i == 0) {\r\n" +
171         "    if ((j >= 35 && j <= ny)) M[i][j].u = M_max;\r\n" +
172         "    else M[i][j].u = M_min;\r\n" +
173         "    M[i][j].BC = true;\r\n" +
174         "} else if (i == nx - 1) {\r\n" +
175         "    if (Math.round(j - j % 10) % 8 == 0) M[i][j].u = M_min;\r\n" +
176         "    else M[i][j].u = M_max;\r\n" +
177         "M[i][j].BC = true;\r\n" +
178         "}\r\n"
179        ,"if (j >= 20 && j <= 29) {\r\n" +
180         "    if (i == nx / 2 - 1) {\r\n" +
181         "        M[i][j].u = M_max;\r\n" +
182         "        M[i][j].BC = true;\r\n" +
183         "    } else if (i == nx / 2) {\r\n" +
184         "        M[i][j].u = M_min;\r\n" +
185         "        M[i][j].BC = true;\r\n" +
186         "    }\r\n" +
187         "}\r\n"
188        ,"if (j >= 24 && j <= 25 && (i == nx / 2 - 1 || i == nx / 2))\r\n" +
189         "    add_BC2(function(t){return sin(t * 10) * 10;}, i, j);\r\n"
190 //       ,"        if (((i == 20 || i == 30) && j >= 40) || ((j == ny - 1 || j == 40) && i > 20 && i < 30)) {\r\n" +
191 //        "            M[i][j].u = 0;\r\n" +
192 //        "            M[i][j].BC = true;\r\n" +
193 //        "        }\r\n" +
194 //        "        //M[25][40].BC = false;\r\n" +
195 //        "        M[25][40].u = M_min;\r\n" +
196 //        "        if (i > 20 && i < 30 && j > 40 && j < 49) {\r\n" +
197 //        "            M[i][j].u = M_min;\r\n" +
198 //        "            M[i][j].BC = true;\r\n" +
199 //        "        }\r\n"
200     ];
201     function generate_function_buttons() {
202         for (var i = 0; i < functions_for_buttons.length; i++) {
203             var btn = document.createElement('input');
204             btn.type = "button";
205             btn.value = i + 1;
206             btn.onclick = function(i) {
207                 return function() {
208                 IC_text.value =
209 //                    "for (var i = 0; i < nx; i++) {\r\n" +
210 //                    "    for (var j = 0; j < ny; j++) {\r\n" +
211                         functions_for_buttons[i];// +
212 //                    "    }\r\n" +
213 //                    "}";
214                     apply_IC();
215                 }
216             }(i);
217             generated_buttons_span.appendChild(btn);
218         }
219     }
220     generate_function_buttons();
221     IC_text.value = functions_for_buttons[0];
222 
223 
224     // граничные условия
225     function BC_periodic() {
226         for (i = 1; i < Nx - 1; i++) {
227             T[i][0] = T[i][Ny - 2];
228             T[i][Ny - 1] = T[i][1];
229         }
230         for (j = 0; j < Ny; j++) {
231             T[0][j] = T[Nx - 2][j];
232             T[Nx - 1][j] = T[1][j];
233         }
234     }
235     function BC_mirror() {
236         for (i = 1; i < Nx - 1; i++) {
237             T[i][0] = T[i][1];
238             T[i][Ny - 1] = T[i][Ny - 2];
239         }
240         for (j = 0; j < Ny; j++) {
241             T[0][j] = T[1][j];
242             T[Nx - 1][j] = T[Nx - 2][j];
243         }
244     }
245     BC_mirror();
246 
247     // дублирующий массив без ГУ и переменные, для внешнего интерфейса
248     var M = [];
249     for (var i = 0; i < Nx - 2; i++) {
250         M[i] = [];
251         for (var j = 0; j < Ny - 2; j++) {
252             M[i][j] = T[i + 1][j + 1];
253         }
254     }
255     var nx = Nx - 2;
256     var ny = Ny - 2;
257 
258     // массив для граничных условий Коши 2-го рода
259     var BC2 = [];
260     function add_BC2(func, i, j) {
261         BC2.push({F:func, i:i, j:j});
262     }
263 
264     var color_N = 120;                           // цветов не больше, чем color_N, саму переменную color_N в расчетах лучше не использовать
265     var colors = prepare_colors(color_N);
266     var cell_pics = prepare_cell_pics(colors);
267 
268     function control() {
269         if (!pause) {
270             calculate_steps(spf, dt);
271             draw();
272         }
273     }
274     setInterval(control, 1000 / fps);                       // Запуск системы
275 
276     // Работа с мышью
277 
278     Wave_Heat_2D_canvas.onmousedown = function(e){           // функция при нажатии клавиши мыши
279         var T1;
280         if (e.which == 1) T1 = T_max * 0.3;        // при нажатии левой клавиши мыши клетка нагревается
281         else if (e.which == 2) T1 = 0;        // при нажатии средней клавиши мыши клетка нагревается
282         else if (e.which == 3) T1 = T_min * 0.3;              // при нажатии правой клавиши мыши клетка остывает
283         else return false;
284 
285         var area = (e.which != 2);
286         set_cell(e, T1, area);
287         Wave_Heat_2D_canvas.onmousemove = function(e) {set_cell(e, T1, area);};   // функция, выполняющаяся при перемещении курсора мыши
288         return false;
289     };
290 
291     document.onmouseup = function(e){           // функция при отпускании клавиши мыши
292         Wave_Heat_2D_canvas.onmousemove = null;              // когда клавиша отпущена - функции перемещения нету
293     };
294 
295     function set_cell(e, T1, area){       // придать клетке определенное состояние с нажатия клавиши мыши
296         var m = mouse_coords(e);                // обновить координаты в переменных mouseX, mouseY
297         if (m.x < 0 || m.x >= canv.w || m.y < 0 || m.y >= canv.h) return;         // проверка на ошибочные координаты
298         var i = Math.floor(m.x / cell_w) + 1;       // получаем ячейку по горизонтали
299         var j = Math.floor(m.y / cell_h) + 1;       // получаем ячейку по вертикали
300         // везде прибавляем 1 - из за периодических условий массив сдвинут на 1
301         if (!area) {
302             if (!T[i][j].BC) {
303                 T[i][j].u = T1;
304                 draw_cell(i, j);
305             }
306         }
307         else for (var shift_i = -3; shift_i <= 3; shift_i++) {
308             var ii;
309             if (i + shift_i < 1) {
310                 if (!radio_BC_periodic.checked) continue;
311                 ii = i + shift_i + Nx - 2;
312             }
313             else if (i + shift_i > Nx - 2) {
314                 if (!radio_BC_periodic.checked) continue;
315                 ii = i + shift_i - Nx + 2;
316             }
317             else ii = shift_i + i;
318             for (var shift_j = -3; shift_j <= 3; shift_j++) {
319                 var jj;
320                 if (j + shift_j < 1) {
321                     if (!radio_BC_periodic.checked) continue;
322                     jj = j + shift_j + Ny - 2;
323                 }
324                 else if (j + shift_j > Ny - 2) {
325                     if (!radio_BC_periodic.checked) continue;
326                     jj = j + shift_j - Ny + 2;
327                 }
328                 else jj = shift_j + j;
329 
330                 if (T[ii][jj].BC) continue;
331 
332                 var r = Math.abs(shift_i) + Math.abs(shift_j);
333 
334                 var T_new;
335                 if (r == 0) T_new = T1;
336                 else if (r == 1) T_new = T1 / 20 * 15;
337                 else if (r == 2) T_new = T1 / 20 * 6;
338                 else if (r == 3) T_new = T1 / 20;
339                 else continue;
340                 
341                 T[ii][jj].u += T_new;
342                 if (T[ii][jj].u > T_max) T[ii][jj].u = T_max;
343                 if (T[ii][jj].u < T_min) T[ii][jj].u = T_min;
344                 draw_cell(ii, jj);
345             }
346         }
347     }
348 
349     function mouse_coords(e) {                  // функция возвращает экранные координаты курсора мыши
350         var m = [];
351         var rect = Wave_Heat_2D_canvas.getBoundingClientRect();
352         m.x = e.clientX - rect.left;
353         m.y = e.clientY - rect.top;
354         return m;
355     }
356 
357     function calculate_steps(spf, dt) {
358         for (var s = 0; s < spf; s++) {
359             for (var i = 1; i < Nx - 1; i++) {
360                 for (var j = 1; j < Ny - 1; j++) {
361                     if (T[i][j].BC) continue;
362                     if (heat_equation) T[i][j].v = (T[i + 1][j].u + T[i][j + 1].u - 4 * T[i][j].u + T[i - 1][j].u + T[i][j - 1].u) * koeff1;
363                     else T[i][j].v += ((T[i + 1][j].u + T[i][j + 1].u - 4 * T[i][j].u + T[i - 1][j].u + T[i][j - 1].u) * koeff1 - T[i][j].v * damping_power) * dt;
364                 }
365             }
366 
367             for (var i = 1; i < Nx - 1; i++) {
368                 for (var j = 1; j < Ny - 1; j++) {
369                     T[i][j].u += T[i][j].v * dt;
370                 }
371             }
372             for (var i = 0; i < BC2.length; i++) {
373                 var f = BC2[i];
374                 M[f.i][f.j].u = f.F(t);
375             }
376             t += dt;
377         }
378     }
379 
380     function draw() {
381         canv.ctx.clearRect(0, 0, canv.w, canv.h);               // очистить экран
382         for (var i = 1; i <= Nx - 1; i++) {
383             for (var j = 1; j <= Ny - 1; j++) {
384                 var col = Math.round((T[i][j].u - T_min) * (colors.length - 1) / (T_max - T_min));            // 0 <= col <= color_N
385                 if (col < 0) col = 0;
386                 if (col >  (colors.length - 1)) col =  (colors.length - 1);
387                 canv.ctx.drawImage(cell_pics[col], (i - 1) * cell_w, (j - 1) * cell_h);
388             }
389         }
390     }
391     function draw_cell(i, j) {
392         var col = Math.round((T[i][j].u - T_min) * (colors.length - 1) / (T_max - T_min));            // 0 < col < color_N
393         canv.ctx.drawImage(cell_pics[col], (i - 1) * cell_w, (j - 1) * cell_h);
394     }
395 
396     function RGBtoHEX(r, g, b) {
397         return "#" + ((1 << 24) + (r << 16) + (g << 8) + b).toString(16).slice(1);
398     }
399     function prepare_colors(N) {
400         var col = [];
401         var n = N / 6 - 1;
402         for (var i = 0; i < 256; i += Math.floor(255 / n))
403             col.push(RGBtoHEX(255 - i, 255, 255));
404         for (var i = 0; i < 256; i += Math.floor(255 / n))
405             col.push(RGBtoHEX(0, 255 - i, 255));
406         for (var i = 0; i < 256; i += Math.floor(255 / n))
407             col.push(RGBtoHEX(0, 0, 255 - i));
408         for (var i = 0; i < 256; i += Math.floor(255 / n))
409             col.push(RGBtoHEX(i, 0, 0));
410         for (var i = 0; i < 256; i += Math.floor(255 / n))
411             col.push(RGBtoHEX(255, i, 0));
412         for (var i = 0; i < 256; i += Math.floor(255 / n))
413             col.push(RGBtoHEX(255, 255, i));
414         return col;
415     }
416     function prepare_cell_pics(col) {
417         var pics = [];
418         for (var i = 0; i < col.length; i++) {
419             var cell = document.createElement('canvas');
420             var cell_ctx = cell.getContext('2d');
421 
422             cell_ctx.fillStyle = col[i];
423             cell_ctx.beginPath();
424             cell_ctx.rect(0, 0, cell_w, cell_h);
425             cell_ctx.fill();
426             pics.push(cell);
427         }
428         return pics;
429     }
430 
431     function init_canvas(canvas) {
432         canvas.onselectstart = function () {return false;};     // запрет выделения canvas
433         canvas.oncontextmenu = function () {return false;};      // блокировка контекстного меню
434 
435         var canv_obj = {};
436         canv_obj.ctx = canvas.getContext("2d");                  // на context происходит рисование
437         canv_obj.w = canvas.width;          // ширина окна в расчетных координатах
438         canv_obj.h = canvas.height;         // высота окна в расчетных координатах
439 
440         return canv_obj;
441     }
442 
443 }

 1 <!DOCTYPE html>
 2 <html>
 3 <head>
 4     <meta charset="UTF-8" />
 5     <title>Wave 2D</title>
 6     <script src="Wave_Heat_2D.js"></script>
 7 </head>
 8 <body>
 9     <table align=center>
10         <tr>
11             <td><canvas id="Wave_Heat_2D_canvas" width="650" height="650" style="border:1px solid #000000; border-radius:6px; margin-bottom: -5px"></canvas><br>
12             <td style="font:bold 12px sans-serif; color:#003366; background-color: #bbbbdd; border-radius:10px; padding: 8px; padding-right: 12px; vertical-align:top">
13                 <div style="text-align: center; font: bold italic 20px Arial, 'Helvetica Neue', Helvetica, sans-serif;">Опции</div>
14                 <hr>
15                 <input type="radio" id="radio_heat" name="eq_radio" /> Уравнение теплопроводности<br>
16                 <input type="radio" id="radio_wave" name="eq_radio" checked/> Волновое уравнение<br><br>
17                 <span id="damping_span">Сила затухания:</span><br><input type="range" id="slider_damping_power" style="width: 150px;"><br><br>
18                 <input type="button" id="button_pause" value="Пауза/Старт"/><br>
19                 <input type="button" id="button_clear" value="Очистить"/><br><br>
20                 <hr>
21                 <span style="color: #000000">Граничные условия:</span><br>
22                 <input type="radio" id="radio_BC_mirror" name="BC_radio" checked/>Зеркальные<br>
23                 <input type="radio" id="radio_BC_periodic" name="BC_radio" />Периодические<br>
24             </td>
25         </tr>
26         <tr>
27             <td colspan="2" style="font:bold 12px sans-serif; color:#333333; background-color: #cccccc; border-radius:10px; padding: 8px; padding-right: 12px; vertical-align:top">
28                 <div style=" padding-bottom: 6px">
29                 <abbr style="border-bottom: 1px dotted black;"
30                       title="Сюда можно вписать код на языке JavaScript, который будет задавать двумерные ГУ системы.&#10;Задавать можно следующие значения:&#10;    M[i][j].u - температура/перемещение;&#10;    M[i][j].v - первая производная (имеет смысл задавать только для волнового уравнения).&#10;Также можно использовать следующие переменные:&#10;    nx, ny - количество частиц по горизонтали и по вертикали;&#10;    M_max, M_min - максимальное и минимальное значение, возможное в массиве (на эти значения настроен цветовой диапазон программы, значения, выходящие за пределы, могут отображаться некорректно).&#10;Превратить участок в ГУ можно с помощью M[i][j].BC = true;">Начальные условия:</abbr>
31                 </div>
32                 for (var i = 0; i < nx; i++) {<br>
33                 &nbsp;&nbsp;&nbsp;&nbsp;for (var j = 0; j < ny; j++) {
34                 <textarea id="IC_text" cols="40" rows="5" wrap=off style="width:95%; height:5.3em; resize: vertical;"></textarea><br>
35                 &nbsp;&nbsp;&nbsp;&nbsp;}<br>
36                 }<br>
37                 <input type="button" id="apply_IC_button" value="Применить" style="margin-top: 8px"/>
38                 <input type="button" id="clear_IC_window_button" value="Очистить поле для ввода"/><br><br>
39                 <span style="padding: 0 4px 0 0;">Примеры:</span><span id="generated_buttons_span"></span>
40             </td>
41             <td></td>
42         </tr>
43         <tr><td colspan="2" style="font:italic bold 12px Georgia, sans-serif; color:#ffffff; background-color: #888888; border-radius:6px; padding: 8px;">
44             <a href="http://tm.spbstu.ru/%D0%A6%D0%B2%D0%B5%D1%82%D0%BA%D0%BE%D0%B2_%D0%94%D0%B5%D0%BD%D0%B8%D1%81_%D0%92%D0%B0%D0%BB%D0%B5%D1%80%D1%8C%D0%B5%D0%B2%D0%B8%D1%87">Denis Tcvetkov</a>
45             &
46             <a href="http://tm.spbstu.ru/%D0%90.%D0%9C._%D0%9A%D1%80%D0%B8%D0%B2%D1%86%D0%BE%D0%B2">Anton Krivtsov</a>
47             &copy 2014
48         </td></tr>
49 
50     </table>
51 </body>
52 </html>