Дискретно-аналитический метод решения задачи теплопроводности

Математическая постановка задачи имеет вид:

(4.37)

где

– координата по толщине стены, ; – координата по времени, ; – значение температуры в точке во время ; – коэффициент температуропроводности материала; – функция, характеризующая мощность возможного источника тепла.

Задача (4.37) определена в пространственно-временной области W:

. (4.38)

Отметим, что

; . (4.39)

Заметим, что поскольку задача (4.37) содержит начальные условия по времени, то она является задачей Коши.

Ниже рассмотрим дискретно-аналитический метод решения задачи, который состоит в следующем: по оси x осуществляется конечно-разностная аппроксимация, а по времени t рассматривается непрерывная (континуальная) задача.

Пусть , – координаты точек разбиения, причем и – граничные точки (в которых заданы краевые условия). Таким образом, искомыми будут являться функции , во внутренних узлах сетки. Схема аппроксимации пространственно-временной области в данном случае условно показана на рис. 5.5.1.

Рис. 4.1. Пространственно-временная область:

Во всех внутренних точках узлах уравнение теплопроводности в (4.37) примет вид:

(4.40)

при этом пусть

, . (4.41)

В соответствии с краевыми условиями из (4.37) для граничных точек, в свою очередь, можем записать:

, , . (4.42)

Следовательно, уравнения теплопроводности для узлов с номерами и имеют соответственно вид:

; (4.43)

. (4.44)

Введем обозначения:

; , (4.45)

где

; . (4.46)

Получаем матричную формулировку разрешающей системы уравнений:

(4.47)

где

; ; . (4.48)

Согласно (4.31)-(4.32) общее решение задачи (4.47) имеет вид:

. (4.49)

Если не зависит от , переходим к формуле

.

Выполняем интегрирование:

,

откуда

. (4.50)

Реализация формулы (4.50) предполагает вычисление экспоненты от матрицы , для выполнения которого следует воспользоваться результатами предыдущего параграфа (см. формулу (4.17)). Имеем:

, (4.51)

где

– матрица собственных векторов матрицы ;

– обратная матрица к матрице ;

; (4.52)

– собственные числа матрицы , .

Аналогично можем вычислить

, где . (4.53)

Варианты задания.

– функция, характеризующая мощность источника тепла;

– коэффициент температуропроводности материала стены;

– краевые условия;

– начальные условия;

– толщина стены; – номер группы, – номер студента по журналу.

Текст М-файла

function teplo_1_expm

s=input('s=');

g=input('g=');

n=input('n=');

l=input('l=');

alpha=input('alpha=');

h=l/n;

c=alpha/h^2;

m=n-1; a1=ones(m-1,1);

A=diag(a1,-1)-2*eye(m)+diag(a1,1),A=c*A;

u0=g; ul=s;

F=c*[u0;zeros(m-2,1);ul];

xi=(0:h:l)';x=xi(2:n);

U0=g+(g+3*s)*x-2*(g+s)*x.^2;

t=[0 0.15 1.5];

nt=length(t);res=zeros(nt,n+1);

fprintf('\n значения функции температуры U(x,t)\n')

for i=1:nt

res(i,:)=[u0 ut(t(i),F,U0,A)' ul];

fprintf('U(%4.2f):',t(i)),fprintf('%6.2f',res(i,:)),fprintf('\n')

end

hold on

plot(xi,res(1,:),'.-')

plot(xi,res(2,:),'o-.r')

plot(xi,res(nt,:),'*:g')

grid on

s1=sprintf('t=%2.0f',t(1));

s2=sprintf('t=%4.2f',t(2));

s3=sprintf('t=%4.2f',t(nt));

legend(s1,s2,s3,0)

title(sprintf('U(x,t)=exp(At)*U0-inv(A)*(E-exp(At))*F\n%s %s %s', s1,s2,s3))

function U=ut(t,F,U0,A)

m=size(A);E=eye(m);

eAt=expm(t*A);

U=eAt*U0-A\(E-eAt)*F;

Замечание. Здесь вычисление функции от матрицы реализуется с использованием стандартной функции expm(At), в которой используется алгоритм

Результаты счета

s=12

g=3

n=8

l=1

alpha=1

A =

-2 1 0 0 0 0 0

1 -2 1 0 0 0 0

0 1 -2 1 0 0 0

0 0 1 -2 1 0 0

0 0 0 1 -2 1 0

0 0 0 0 1 -2 1

0 0 0 0 0 1 -2

значения функции температуры U(x,t)

U(0.00): 3.00 7.41 10.88 13.41 15.00 15.66 15.38 14.16 12.00

U(0.15): 3.00 4.81 6.52 8.03 9.29 10.28 11.02 11.56 12.00

U(1.50): 3.00 4.13 5.25 6.38 7.50 8.63 9.75 10.88 12.00

>>