ИДЗ_1_С прямым методом решения СЛАУ
.docxМИНОБРНАУКИ РОССИИ
Санкт-Петербургский государственный
электротехнический университет
«ЛЭТИ» им. В.И. Ульянова (Ленина)
Кафедра МНЭ
отчёт
по индивидуальному заданию №1
по дисциплине «Моделирование и проектирование микро- и наносистем»
Тема: ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ НЕСТАЦИОНАРНОГО ПРОЦЕССА ТЕПЛОПРОВОДНОСТИ В НЕОДНОРОДНОМ ТЕЛЕ С ИСПОЛЬЗОВАНИЕМ ПРЯМЫХ МЕТОДОВ РЕШЕНИЯ СЛАУ
Вариант №16
Студентка гр. 9282 |
|
Зикратова А. А. |
Преподаватель |
|
Рындин Е. А. |
Санкт-Петербург
2022
Цель работы.
Нахождение распределения температуры в многослойной структуре в каждый момент времени из заданного диапазона прямым методом решения СЛАУ.
Задание.
Рис. 1 – Общий вид неоднородной структуры
Таблица 1 – Исходные данные
№ п/п |
L, мм |
W, мм |
Коэффициент теплопроводности k, Вт/(м К) |
Плотность , кг/м3 |
Удельная теплоемкость C, Дж/(кг К) |
Начальное и граничные условия (род) gt gxmin gxmax gymin gymax, К или К/м |
Зависимость плотности мощности источников (стоков) тепла от температуры f(T), Вт/м3 |
Интервал времени t, c |
16. |
30 150 230 |
10 10 20 10 50 10 20 30 30 130 |
1000 100 200 100 300 10 10000 |
10000 2000 1000 2000 1000 2000 300
|
10000 100 1000 100 1000 100 300 |
300(1) 50(2) 50(2) 300(1) 0(2)
|
-1107 (ИДЗ №1)
-1107/lgT (ИДЗ №2,3)
|
20 |
Теоретические положения
Ⅰ. Дифференциальное уравнение теплопроводности, начальные и граничные условия в обычном виде:
- -
Начальное условие:
Г
} y = [ymin, ymax], t = (tmin, tmax]
раничные условия:
1)
2)
3
} x = (xmin, xmax), t = (tmin, tmax]
)
4)
Ⅱ. Дифференциальное уравнение теплопроводности, начальные и граничные условия в дискретном виде:
Начальное условие:
,
Граничные условия:
1
} j = 1…J, m = 2…M
} i = 2…I - 1, m = 2…M
) 3) =2) 4) =
Данное уравнение в дискретном виде решается путём представления его в матричной форме: AX = B → X = A-1B
A – матрица коэффициентов, X – вектор-столбец искомых переменных (температур), B – вектор-столбец свободных членов.
Программа в Matlab:
clear all
close all
clc
L=[30 150 230];
L=L*1e-3;
L(3)=L(3)-L(2)-L(1);
W=[10 10 20 10 50 10 20 30 30 130]*1e-3;
kL=[1000 100 200 100 300 10 10000];
rL=[10000 2000 1000 2000 1000 2000 300];
cL=[10000 100 1000 100 1000 100 300];
gt=300;
gxmin=50;
gxmax=50;
gymin=300;
gymax=0;
F=-1e7;
tmax=20;
St=5;
t=0:tmax/St:tmax;
dt=diff(t);
M=length(t);
Sx=7;
x(1)=0;
kV(1)=kL(1);
rV(1)=rL(1);
cV(1)=cL(1);
for i=1:length(kL)
x=[x max(x)+W(i)/Sx:W(i)/Sx:max(x)+W(i)];
kV=[kV ones(1,Sx).*kL(i)];
rV=[rV ones(1,Sx).*rL(i)];
cV=[cV ones(1,Sx).*cL(i)];
end
I=length(x);
dx=diff(x);
kV=kV';
rV=rV';
cV=cV';
Sy=9;
y(1)=0;
for i=1:length(L)
y=[y max(y)+L(i)/Sy:L(i)/Sy:max(y)+L(i)];
end
dy=diff(y);
J=length(y);
k=kV;
r=rV;
c=cV;
for j=2:J
k=[k kV];
r=[r rV];
c=[c cV];
end
f=zeros(I,J,M);
for m=1:M
for j=Sy+1: 2*Sy+1
for i=1:I
if x(i)>=W(9) && x(i)<=W(9)+W(8)
f(i,j,m)=F;
end
end
end
end
A=zeros(I*J*M,I*J*M);
B=zeros(I*J*M,1);
for i=1:I
for j=1:J
A(I*J*(1-1)+I*(j-1)+i,I*J*(1-1)+I*(j-1)+i)=1;
B(I*J*(1-1)+I*(j-1)+i)=gt;
end
end
for j=1:J
for m=2:M
A(I*J*(m-1)+I*(j-1)+1, I*J*(m-1)+I*(j-1)+1)=-1/dx(1);
A(I*J*(m-1)+I*(j-1)+1, I*J*(m-1)+I*(j-1)+2)=1/dx(1);
B(I*J*(m-1)+I*(j-1)+1)=gxmin;
A(I*J*(m-1)+I*(j-1)+I, I*J*(m-1)+I*(j-1)+I)=1/dx(I-1);
A(I*J*(m-1)+I*(j-1)+I, I*J*(m-1)+I*(j-1)+I-1)=-1/dx(I-1);
B(I*J*(m-1)+I*(j-1)+I)=gxmax;
end
end
for i=2:I-1
for m=2:M
A(I*J*(m-1)+I*(1-1)+i, I*J*(m-1)+I*(1-1)+i)=1;
B(I*J*(m-1)+I*(1-1)+i)=gymin;
end
end
for i=2:I-1
for m=2:M
A(I*J*(m-1)+I*(J-1)+i, I*J*(m-1)+I*(J-1-1)+i)=-1/dy(J-1);
A(I*J*(m-1)+I*(J-1)+i, I*J*(m-1)+I*(J-1)+i)=1/dy(J-1);
B(I*J*(m-1)+I*(J-1)+i)=gymax;
end
end
for i=2:I-1
for j=2:J-1
for m=2:M
A(I*J*(m-1)+I*(j-1)+i,I*J*(m-1)+I*(j-1)+i)=r(i,j)*c(i,j)/dt(m-1)-...
2/(dx(i)+dx(i-1))*(-k(i,j)/dx(i)-k(i-1,j)/dx(i-1))-...
2/(dy(j)+dy(j-1))*(-k(i,j)/dy(j)-k(i,j-1)/dy(j-1));
A(I*J*(m-1)+I*(j-1)+i,I*J*(m-1)+I*(j-1)+i+1)=-2/(dx(i)+dx(i-1))*k(i,j)/dx(i);
A(I*J*(m-1)+I*(j-1)+i,I*J*(m-1)+I*(j-1)+i-1)=-2/(dx(i)+dx(i-1))*k(i-1,j)/dx(i-1);
A(I*J*(m-1)+I*(j-1)+i,I*J*(m-1)+I*(j+1-1)+i)=-2/(dy(j)+dy(j-1))*(k(i,j)/dy(j));
A(I*J*(m-1)+I*(j-1)+i,I*J*(m-1)+I*(j-1-1)+i)=-2/(dy(j)+dy(j-1))*(k(i,j-1)/dy(j-1));
A(I*J*(m-1)+I*(j-1)+i,I*J*(m-1-1)+I*(j-1)+i)=-r(i,j)*c(i,j)/dt(m-1);
B(I*J*(m-1)+I*(j-1)+i)=f(i,j,m);
end
end
end
TV=A^(-1)*B;
for i=1:I
for j=1:J
for m=1:M
T(i,j,m)=TV(I*J*(m-1)+I*(j-1)+i);
end
end
end
NN='Graphic_';
for m=1:M
NNN=[NN num2str(m)];
figure
mesh(y.*1e3,x.*1e3,T(:,:,m)-273);
xlabel('y, mm','FontSize',19)
ylabel('x, mm','FontSize',19)
zlabel('T, ^oC','FontSize',19)
xlim([min(y.*1e3) max(y.*1e3)])
ylim([min(x.*1e3) max(x.*1e3)])
zlim([-20 100])
grid on
colormap([0 0 0])
print(gcf,'-djpeg',NNN)
end
Результаты моделирования:
Рис. 2 – Пространственное распределение температуры в неоднородном теле при t = 0 с (прямой метод)
Рис. 3 – Пространственное распределение температуры в неоднородном теле при t = 4 с (прямой метод)
Рис. 4 – Пространственное распределение температуры в неоднородном теле при t = 8 с (прямой метод)
Рис. 5 – Пространственное распределение температуры в неоднородном теле при t = 12 с (прямой метод)
Рис. 6 – Пространственное распределение температуры в неоднородном теле при t = 16 с (прямой метод)
Рис. 7 – Пространственное распределение температуры в неоднородном теле при t = 20 с (прямой метод)
Вывод: в данной работе программа выполняет численное решение уравнения теплопроводности для неоднородного тела с использованием прямого метода решения СЛАУ. С помощью этого метода можно решить задачу одной матричной операцией, но для этого затрачивается внушительный объём оперативной памяти компьютера даже при малом количестве шагов координатных и временных сеток.