3.2. Расчет нестационарного распределения температуры

Задание: рассчитать изменение распределения температуры в области ABCD, приведенной на рисунке 1, в течение 10 временных шагов. Отрезок AЕ – контакт K1, CF – контакт K2. Начальное распределение температуры нулевое. В начальный момент времени температура контакта К1 мгновенно повышается до значения Т1, контакта К2 – до Т2. Температуропроводность материала – а, шаг пространственной дискретизации — h, временной дискретизации – t, число временных шагов – k_T. Для обеспечения высокой точности решения задачи (чтобы получить хорошую консервативность) СЛАУ, возникающие на каждом временном шаге, решаются с помощью функции LinearSolve. Программа, решающая эту задачу:

(1) (2) (3) (4) (5)

ktmod=10; temp1=40.; temp2=10.; a=1.; h=1.; t=1.; xb=6.; ib=IntegerPart[xb/h]+1; ya=5.; ja=IntegerPart[ya/h]+1; yc=2.; jc=IntegerPart[yc/h]+1; yd=3., jd=IntegerPart[yd/h]+1; neizv=ib*ja; temp=Table[0,{i, neizv}]; temp=Partition[temp, ib];

30

(6) (7) (8) (9) (10) (11) (12) (13) (14) (15) (16) (17) (18) (19) (20) (21) (22) (23) (24) (25) (26) (27) (28) (29) (30) (31) (32) (33) (34) (35) (36) (37) (38) (39) (40) (41) (42) (43) (44) (45) (46)

For[j=1, jjc, j++, temp[[j,1]]=temp1]; For[j=jd, jja, j++, temp[[j,ib]]=temp2]; nurvn=neizv-jc+jd-ja-1; m=Table[0,{i, nurvn*nurvn}]; m=Partition[m,nurvn]; b=Table[0,{i, nurvn}]; x=Table[0,{i, nurvn}]; For[k=1, kktmod, k++, nur=0; For[i=1, iib, i++, For[j=1, jja, j++, If[(j==1)&&(i>1)&&(i<ib), nur++; m[[nur,nur]]=0.5*h*h+2*a*t; m[[nur,nur+ja]]=-0.5*a*t; m[[nur,nur+1]]=-a*t; b[[nur]]=0.5*h*h*temp[[j,i]; If[i>2, m[[nur,nur-ja]]=-0.5*a*t]; If[i==2, b[[nur]]= b[[nur]]+0.5*a*t*temp1]; Continue[]]; If[(i==ib)&&(j==1), nur++; m[[nur,nur]]=0.25*h*h+a*t; m[[nur,nur-ja]]=-0.5*a*t; m[[nur,nur+1]]=-0.5*a*t; b[[nur]]=0.25*h*h*temp[[j,i]]; Continue[]]; If[(i==1)&&(j>jc)&&(j<ja), nur++; m[[nur,nur]]=0.5*h*h+2*a*t; m[[nur,nur+1]]=-0.5*a*t; m[[nur,nur+ja]]=-a*t; b[[nur]]=0.5*h*h*temp[[j,i]]; If[j>jc+1, m[[nur,nur-1]]=-0.5*a*t; If[j==jc+1, b[[nur]]=b[[nur]]+0.5*a*t*temp1]; Continue[]]; If[(i==1)&&(j==ja), nur++; m[[nur,nur]]=0.25*h*h+a*t; m[[nur,nur-1]]=-0.5*a*t; m[[nur,nur+ja]]=-0.5*a*t; b[[nur]]=0.25*h*h*temp[[j,i]]; Continue[]]; If[(i==ib)&&(j>1)&&(j<jd), nur++; m[[nur,nur]]=0.5*h*h+2*a*t; m[[nur,nur-ja]]=-a*t; m[[nur,nur-1]]=-0.5*a*t; b[[nur]]=0.5*h*h*temp[[j,i]; If[j<jd-1, m[[nur,nur+1]]=-0.5*a*t]; If[i==jd-1, b[[nur]]= b[[nur]]+0.5*a*t*temp2]; Continue[]]; If[(j==ja)&&(i>1)&&(i<ib), nur++; m[[nur,nur]]=0.5*h*h+2*a*t; m[[nur,nur-ja]]=-0.5*a*t; m[[nur,nur-1]]=-a*t; b[[nur]]=0.5*h*h*temp[[j,i]; If[i<ib-1, m[[nur,nur+ja]]=-0.5*a*t]; If[i==ib-1, b[[nur]]= b[[nur]]+0.5*a*t*temp2]; Continue[]]; If[(i>1)&&(i<ib)&&(j>1)&&(j<ja), nur++;

31

(47) (48) (49) (50) (51) (52) (53) (54) (55) (56) (57) (58) (59)

m[[nur,nur]]=h*h+4*a*t; m[[nur,nur-1]]=-a*t; m[[nur,nur+1]]=-a*t; b[[nur]]=h*h*temp[[j,i]; If[(i>2)||(j>jc), m[[nur,nur-ja]]=-a*t]; If[(i==2)&&(jjc), b[[nur]]= b[[nur]]+a*t*temp1]; If[(i<ib-1)||(j<jd), m[[nur,nur+ja]]=-a*t]; If[(i==ib-1)&&(jjd), b[[nur]]= b[[nur]]+a*t*temp2];] ]]; x=LinearSolve[m,b]; nur=0; For[i=1, iib, i++, For[j=1, jja, j++, If[!(((j==1)&&(j1)&&(jjc))||((i==ib)&&(jjd)&&(jja))), nur++; temp[[j,i]]=x[[nur]]]; ]]; ListPlot3D[temp,ViewPoint->{1.2,1.2,1.2}]];

Подготовка к решению задачи аналогична той, что была в предыдущем пункте. Поэтому поясним только смысл переменных: ktmod – число временных шагов; temp1, temp2 – температура 1-го и 2-го контактов; xb – длина отрезка AB; ya, yc, yd – длина отрезков AD, AC и AD; ib, ja, jc,jd – соответствующие счетчики; neizv – число узлов разностной сетки; temp – двумерный массив, содержащий рассчитываемое распределение; nurvn – размерность решаемой СЛАУ; m – матрица СЛАУ; b – вектор-столбец правых частей; x – вектор-столбец решений.

Искомое распределение будет размещаться в двумерном массиве temp. Для обращения к элементу этого массива будем использовать два индекса i и j, которые являются и номерами узлов сетки по осям (i – по оси x, j – по оси y). Однако при обращении к функции LinearSolve искомое решение (значения температуры в узлах сетки) будет размещено в одномерном массиве x. Поэтому, в дополнение к упомянутым индексам, введем еще один, nur, который будем использовать при сквозной нумерации узлов по мере их обхода. Порядок обхода следующий. Начинаем с левой границы области, исключая контакт (температура в узлах сетки, попадающих на этот контакт, известна). Это отрезок EB. Его проходим снизу вверх. Далее переходим к ближайшей вертикальной линии сетки, ее проходим тоже снизу вверх (это линия x₂ на рис. 2). Аналогичным образом обрабатываются остальные узлы сетки. Завершаем эту процедуру, пройдя по узлам, попадающим на отрезок DE.

В строке (12) – начало цикла выполнения всех шагов по времени, в строке (13) – начало двойного цикла прохода по всем узлам разностной сетки. Для каждого узла рассчитываются и заносятся в массивы m и b соответствующие коэффициенты. Здесь переменная nur – текущий номер уравнения.

32

Строкой (53) заканчивается проход по всем узлам сетки, СЛАУ сформирована и решается (строка (54)). Полученное решение переносится в массив temp (строки (55)(58)) и отображается на экране (строка (59)).

Эта реализация достаточно компактная по объему текста, но не экономичная по затратам процессорного времени. Дело в том, что временные изменения отражаются только в свободных членах СЛАУ, а матрица СЛАУ содержит постоянные коэффициенты. Поэтому здесь следовало бы один раз рассчитать коэффициенты матрицы (массив m), а в цикле по времени пересчитывать только свободные члены (массив b).