Двумерное уравнение теплопроводности + волновое уравнение
Материал из Department of Theoretical and Applied Mechanics
Виртуальная лаборатория > Двумерное уравнение теплопроводности + волновое уравнение
Численное решение двумерного уравнения теплопроводности и двумерного волнового уравнения.
Управление:
- Левая кнопка мыши - нагреть область
- Средняя кнопка мыши - придать области нормальную температуру
- Правая кнопка мыши - охладить область
Скачать Wave_Heat_2D_v2-4_release.zip.
Текст программы на языке JavaScript (разработчики Цветков Денис, Кривцов Антон):
Файл "Wave_Heat_2D.js"
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 }
Файл "Wave_Heat_2D.html"
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, который будет задавать двумерные ГУ системы. Задавать можно следующие значения: M[i][j].u - температура/перемещение; M[i][j].v - первая производная (имеет смысл задавать только для волнового уравнения). Также можно использовать следующие переменные: nx, ny - количество частиц по горизонтали и по вертикали; M_max, M_min - максимальное и минимальное значение, возможное в массиве (на эти значения настроен цветовой диапазон программы, значения, выходящие за пределы, могут отображаться некорректно). Превратить участок в ГУ можно с помощью M[i][j].BC = true;">Начальные условия:</abbr>
31 </div>
32 for (var i = 0; i < nx; i++) {<br>
33 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 }<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 © 2014
48 </td></tr>
49
50 </table>
51 </body>
52 </html>
