Численное решение уравнения теплопроводности и волнового уравнения
Материал из Department of Theoretical and Applied Mechanics
Виртуальная лаборатория > Численное решение уравнения теплопроводности и волнового уравнения
Скачать программу: Equation v6-9-mini release.zip
Текст программы на языке JavaScript (разработчики Цветков Денис, Кривцов Антон):
Файл "Equations.js"
1 window.addEventListener("load", main_equations, false);
2 function main_equations() {
3
4 // Предварительные установки
5
6 var a0 = 1; // масштаб расстояния
7 var t0 = 1; // масштаб времени
8 var m0 = 1; // масштаб массы
9 var T0 = 1; // масштаб температуры
10
11 var k0 = 2 * Math.PI / t0; // масштаб частоты
12 var C0 = m0 * k0 * k0; // масштаб жесткости
13
14 var m = 1 * m0; // масса
15 var C = 1 * C0; // жесткость
16
17 var omega = Math.sqrt(C / m);
18 var c = a0 * omega;
19
20 var tw = CONFIG.N_Wave * a0 / c / 4;
21 var beta = Math.PI * c * c * Math.sqrt(tw);
22
23 var dx = 1 * a0; // шаг сетки по оси x
24
25 var spf = CONFIG.spf; // steps per frame - сколько расчетов проходит за каждый кадр отображения
26 var dt = CONFIG.dt * t0; // шаг интегрирования по времени
27
28 var c3 = init_canvas(Equations_canvas); // инициализация canvas
29
30 // слайдер
31 slider_dt.max = dt * 2;
32 slider_dt.min = dt / 50;
33 slider_dt.step = dt / 100;
34 slider_dt.value = dt;
35 slider_dt.oninput = function() {dt = parseFloat(slider_dt.value);};
36
37 var heat, wave;
38 function start_new_system() {
39 heat = new Heat_model(T0, dx, dt, beta);
40 wave = new Wave_model(a0, k0, m, c, dx, dt);
41 }
42 function start_new_system_stair() {
43 var func = function(M, M_max, N, val) {return IC_stair(M, M_max, N, val,
44 parseFloat(stair_x1_number.value), parseFloat(stair_x2_number.value))};
45 CONFIG.IC_Heat = func;
46 CONFIG.IC_Wave = func;
47 start_new_system()
48 }
49 function start_new_system_Gauss() {
50 CONFIG.IC_Heat = IC_Gauss;
51 CONFIG.IC_Wave = IC_Gauss;
52 start_new_system()
53 }
54
55
56
57 button_restart_stair.onclick = start_new_system_stair;
58 button_restart_Gauss.onclick = start_new_system_Gauss;
59
60 start_new_system();
61 control();
62
63 function control() {
64 heat.calculate_steps(spf, dt);
65 wave.calculate_steps(spf, dt);
66 draw_interface(c3);
67 if (checkbox_graph_01.checked) draw_graph(c3, heat.get_D(), heat.get_div(), 1, "#ff0000", "u");
68 if (checkbox_graph_02.checked) draw_graph(c3, wave.get_D(), wave.get_div(), 1, "#00ff00", "u");
69 requestAnimationFrame(control);
70 }
71
72 function draw_graph(c, D, D_div, mult_y, color, val) {
73 // график
74 c.ctx.strokeStyle = color;
75 c.ctx.beginPath();
76 c.ctx.moveTo(0, c.h - D[0].u * (c.h - 20));
77 for (var i = 1; i < D.length - 1; i++) {
78 c.ctx.lineTo(i / (D.length - 1) * c.w, c.h - D[i][val] / D_div * mult_y * (c.h - 20));
79 }
80 c.ctx.stroke();
81 }
82
83 function draw_interface(c) {
84 c.ctx.clearRect(0, 0, c.w, c.h);
85
86 // вертикальная линия
87 c.ctx.strokeStyle = "#bbbbbb";
88 c.ctx.beginPath();
89 c.ctx.moveTo(c.w / 2, 0);
90 c.ctx.lineTo(c.w / 2, c.h);
91 c.ctx.stroke();
92 }
93
94 function init_canvas(canvas) {
95 canvas.onselectstart = function () {return false;}; // запрет выделения canvas
96
97 var canv_obj = {};
98 canv_obj.ctx = canvas.getContext("2d"); // на context происходит рисование
99 canv_obj.w = canvas.width; // ширина окна в расчетных координатах
100 canv_obj.h = canvas.height; // высота окна в расчетных координатах
101
102 return canv_obj;
103 }
104 }
105
106
107 function Heat_model(T0, dx, dt, X) {
108 var N = CONFIG.N_Heat + 2; // количество узлов по оси x + 2 для ГУ
109 var T_start = 1 * T0; // начальная температура
110
111 var T = make_mass(N);
112 var IC = CONFIG.IC_Heat;
113 IC(T, T_start, N, "u"); // начальные условия
114 BC_mirrow(T, N); // граничные условия
115
116 var div = 1 * T0; // делитель для функции рисования
117 this.get_div = function() {return div};
118 var t = dt / 2; // суммарное время с начала расчета
119
120 this.get_D = function() {return T};
121 this.calculate_steps = calc_implicit_Euler_method;
122
123
124 var alpha = new Array(N + 2), beta = new Array(N + 2);
125 function calc_implicit_Euler_method(spf, dt) {
126 for (var s = 0; s < spf; s++) {
127 var A = X * dt / (dx * dx + 2 * X * dt);
128 var B = dx * dx / (dx * dx + 2 * X * dt);
129
130 alpha[N - 2] = 0;
131 beta[N - 2] = T[N - 2].u;
132 for (var i = N - 3; i >= 1; i--) {
133 alpha[i] = A / (1 - alpha[i + 1] * A);
134 beta[i] = (B * T[i + 1].u + beta[i + 1] * A) / (1 - alpha[i + 1] * A);
135 }
136
137 T[1].u = (B * T[1].u + beta[1] * A) / (1 - alpha[1] * A - A);
138 for (var i = 2; i < N - 1; i++) {
139 T[i].u = (A * T[i - 1].u + B * T[i].u + beta[i] * A) / (1 - alpha[i] * A);
140 }
141
142 t += dt;
143 }
144 }
145
146 return this;
147 }
148 function Wave_model(a0, k0, m, c, dx, dt) {
149 var N = CONFIG.N_Wave + 2; // количество узлов по оси x + 2 для ГУ
150 var U_start = a0 * k0 / (2 * Math.PI);
151
152 var U = make_mass(N);
153 var IC = CONFIG.IC_Wave;
154 IC(U, U_start, N, "u"); // начальные условия
155 IC_zero(U, U_start, N, "v"); // начальные условия для скоростей
156 BC_mirrow(U, N); // граничные условия
157
158 var div = 1 * a0 * k0 / (2 * Math.PI); // делитель для функции рисования
159 this.get_div = function() {return div};
160 var t = dt / 2; // суммарное время с начала расчета
161
162 this.get_D = function() {
163 if (CONFIG.smooth_wave) return get_smoothed(U, N, "u", CONFIG.averaging_quality_Wave);
164 else return U;
165 };
166 this.calculate_steps = function(spf, dt) {
167 var koeff1 = c * c / (dx * dx) / m;
168 for (var s = 0; s < spf; s++) {
169 for (var i = 1; i < N - 1; i++) {
170 U[i].v += (U[i + 1].u - 2 * U[i].u + U[i - 1].u) * koeff1 * dt;
171 }
172 for (var i = 1; i < N - 1; i++) {
173 U[i].u += U[i].v * dt;
174 }
175 t += dt;
176 }
177 };
178 return this;
179 }
180
181 function get_smoothed(D0, N, val, avl_quality) {
182 var D = [];
183 var avl = (N - 2) / avl_quality;
184 var s = 1;
185 var avl_half = avl / 2;
186 for (var i = avl_half + 1; i < N - 1 ; i += avl) {
187 var av = 0;
188 for (var j = i - avl_half; j < i + avl_half; j++) {
189 av += D0[j][val];
190 }
191 D[s] = {};
192 D[s][val] = av / avl;
193 s++;
194 }
195 // для первой точки добавляется объект, чтобы график располагался правильно по центру
196 // но сама первая точка не добавляется - она не усреднена
197 D[0] = {};
198 return D;
199 }
Файл "Config.js"
1 var CONFIG = {
2 fps: 60 // количество кадров в секунду
3 ,spf: 200 // сколько раз происходит вычисление между кадрами
4
5 ,dt: 0.005
6
7 // начальные условия
8 // IC_zero: все значения - нули
9 // IC_stair: ступенька посередине экрана
10 // IC_half_stair: половина ступеньки
11 // IC_stair_random: ступенька, составленная из случайных значений
12 // IC_half_stair_random: половина ступеньки из случайных значений
13 // IC_Gauss: распределение Гаусса
14 // ,IC_Heat: IC_Gauss
15 // ,IC_Wave: IC_Gauss
16 ,IC_Heat: function(M, M_max, N, val) {return IC_stair(M, M_max, N, val, 0.33, 0.66)}
17 ,IC_Wave: function(M, M_max, N, val) {return IC_stair(M, M_max, N, val, 0.33, 0.66)}
18
19 ,N_Heat: 1000 // количество частиц по оси x
20 ,N_Wave: 1000 // количество частиц по оси x
21
22 ,averaging_quality_Wave: 200 // на сколько отрезков усреднения поделить N (число должно быть делителем N)
23 ,smooth_wave: true // сглаживать волновое уравнение
24 };
Файл "IC_BC_LIB_v1.js"
1 // Аргументы в функции IC и BC:
2 // M - массив (сформированный с помощью функции make_mass данной библиотеки)
3 // M_max - максимальное значение, которое должна задать функция (иногда данное значение передается,
4 // но не испозьзуется, сделано для того, чтобы у всех фунций был одинаковый набор аргументов,
5 // и можно было бы быстро менять используемую функцию (например, из файла настроек))
6 // N - Размер массива
7 // val - ключ, по которому будет организован доступ к данным. В программе можно будет получить данные через него,
8 // например, так: M[0].val
9 // Чтобы лучше разобраться - просмотрите код какой либо программы, использующей эту библиотеку.
10
11 /**
12 * Создает специальный массив для последующего задания в нем начальных и граничных условий
13 * N - Длина массива (нужно учитывать, что первый и последний элементы будут использованы для граничных условий))
14 * */
15 function make_mass(N) {
16 var M = [];
17 for (var i = 1; i < N - 1; i++)
18 M[i] = {};
19 return M;
20 }
21
22 function IC_zero(M, M_max, N, val) { // Заполняет массив нулями
23 for (var i = 1; i < N - 1; i++) {
24 M[i][val] = 0;
25 }
26 }
27
28 // 0 <= start, end <= 1
29 // start <= end
30 // чтобы получить функцию вида IC_stair(M, M_max, N, val), используйте замыкание
31 // например: var func = function(M, M_max, N, val) {return IC_stair(M, M_max, N, val, 0.33, 0.66)}
32 function IC_stair(M, M_max, N, val, start, end) { // Задает "ступеньку" высоты M_max посередине экрана
33 for (var i = 1; i < N - 1; i++) {
34 if (i >= N * start && i < N * end) M[i][val] = M_max;
35 else M[i][val] = 0;
36 }
37 }
38 function IC_half_stair(M, M_max, N, val) { // Задает мгновенный подъем от нуля до максимального значения
39 for (var i = 1; i < N - 1; i++) {
40 if (i < N / 2) M[i][val] = 0;
41 else M[i][val] = M_max;
42 }
43 }
44 function IC_stair_random(M, M_max, N, val) { // Задает "ступеньку" из случайных значений от -M_max до M_max
45 for (var i = 1; i < N - 1; i++) {
46 if (i >= N / 3 && i < 2 * N / 3) M[i][val] = (2 * Math.random() - 1) * M_max;
47 else M[i][val] = 0;
48 }
49 }
50 function IC_half_stair_random(M, M_max, N, val) {// Правая половина экрана будет заполнена случайными значениями от -M_max до M_max
51 for (var i = 1; i < N - 1; i++) {
52 if (i < N / 2) M[i][val] = 0;
53 else M[i][val] = (2 * Math.random() - 1) * M_max;
54 }
55 }
56 function IC_Gauss(M, M_max, N, val) { // Задает распределение Гаусса.
57 for (var i = 1; i < N - 1; i++) {
58 var x = (i - 1) / (N - 2);
59 M[i][val] = Math.exp(-Math.pow(x - 0.5, 2) * 50);
60 }
61 }
62
63 function BC_mirrow(M, N) { // Зеркальные граничные условия
64 M[0] = M[1];
65 M[N - 1] = M[N - 2];
66 }
67 function BC_periodic(M, N) { // Периодические граничные условия.
68 M[0] = M[N - 2];
69 M[N - 1] = M[1];
70 }
Файл "Equations.html"
1 <!DOCTYPE html>
2 <html>
3 <head>
4 <meta charset="UTF-8" />
5 <title>Equations</title>
6 <script src="IC_BC_LIB_v1.js"></script>
7 <script src="Config.js"></script>
8 <script src="Equations.js"></script>
9 </head>
10 <body>
11 <canvas id="Equations_canvas" width="1000" height="500" style="border:1px solid #000000;"></canvas><br>
12 Скорость расчета: <input type="range" id="slider_dt" style="width: 150px;"><br><br>
13 <input type="checkbox" id="checkbox_graph_01" checked/><font color="#ff0000" size="5"><B>—</B></font> Уравнение теплопроводности<br>
14 <input type="checkbox" id="checkbox_graph_02" checked/><font color="#00ff00" size="5"><B>—</B></font> Волновое уравнение<br><br>
15 <input type="button" id="button_restart_Gauss" value="Рестарт Гаусс"/><br><br>
16 <input type="button" id="button_restart_stair" value="Рестарт ступенька"/> (0 <= <b>x1</b>, <b>x2</b> <= 1; <b>x1</b> < <b>x2</b>)<br>
17 x1 = <input type="number" id="stair_x1_number" min=0 max=1 step=0.001 value="0.33"/><br>
18 x2 = <input type="number" id="stair_x2_number" min=0 max=1 step=0.001 value="0.66"/>
19 </body>
20 </html>
