Моделирование процессов переработки пластмасс (стр. 5 из 5)
В общем случае, когда температура зависит от всех трех координат (х, у, г), дифференциальное уравнение теплопроводности конечного цилиндра имеет вид
; (3.13) 
4 СОСТАВЛЕНИЕ АЛГОРИТМА
Для решения дифференциального уравнения теплопроводности бесконечного цилиндра воспользуемся методом сеток, суть которого заключается в разбиении координатной плоскости на равные части и вычислении значения искомой функции в узлах образуемой сетки. Используя значения функции в крайних точках можно последовательно вычислить её значение в любой части координатной плоскости.
; (4.1)Заменим частный дифференциал разностным отношением:
; (4.2)Осуществим следующее преобразование функции:
; (4.3) 
; (4.4) 
; (4.5) 
(4.6) 
; (4.7) 
; (4.8)Подготовим уравнение (4.8) для рекуррентного вычисления в MatLabV6.0
Произведём переобозначения:
; (4.9) 
; (4.10) 
; (4.11) ; (4.12) 
; (4.13)Имеем формулу:
T(i+1,j+1)=T(i,j+1)+(a*dt/dr)*(((T(i,j+2)-2*T(i,j+1)+T(i,j))/dr)+((1/r)*(T(i,j+2)-T(i,j+1)))); (4.14)
В результате последовательных вычислений можно получить массив Tхарактеризующий температурное поле неограниченного цилиндра в любой момент времени.
1.Программа начинается cзадание переменных: начального и конечного момента времени, числа дискретных отсчётов по времени, радиус цилиндра и число его разбиений, констант характеризующих тепло-физические свойства полимера.
2.Следующим этапом является вычисление шага аргументов, по которым будет вычисляться исходная функция.
3.Краевые условия: значения искомой функции в начальный момент времени t0 = 0 в зависимости от радиуса, и температуры стенки литникового канала в любой момент времени задаются циклом For.
4.Каждому элементу вектора характеризующего температурное поле в начальный момент времени присваивается значение температуры, вычисленное как значение функции распределения вложенной в цикл. Число циклов присвоения значений вектору увеличивают на два так-так его элементов на один должно быть больше чем число интервалов разбиений и на одно значение больше, чтобы было возможным вычисление значения массива в центре цилиндра после перехода от внутреннего цикла к внешнему.
5.Каждому элементу вектора характеризующего температуру стенки канала в любой момент времени присваивается постоянное значение температуры Число циклов присвоения значений вектору увеличивают на один, так-так его элементов на один должно быть больше чем число интервалов разбиений.
6.Для вычисления матрицы определяющей температуру цилиндра по радиусу в любой момент времени используем два вложенных цикла For. Во внутреннем цикле предусмотрено изменение радиуса цилиндра, и вычисление температурного поля в заданный момент времени.
7.При переходе к внешнему циклу отсчёт по времени увеличивается на единицу. Значение производной температуры по радиусу в любой момент времени равно нулю и поэтому, чтобы учесть ещё одно краевое условие при переходе от внешнего цикла к внутреннему значение последней температуры копируется два раза.
8.После получения матрицы температур надо построить график. Чтобы координатные оси были проградуированные удобно для использования в матрице температур переставляют столбцы. Осуществляется это с использованием двух вложенных циклов.
9.Далее следует вывод графика и градуировка его осей.
5 СОСТАВЛЕНИЕ ПРОГРАММЫ
Программа для MatLabv6.0 R12 начинается очищения переменных графических окон функций и окна вывода результата. Осуществляют это с помощью: clear, clc, clf, clg
Чтобы программа была легка в использовании и проста в конфигурировании под любые задачи разработаем её используя понятные обозначения:
Задаём переменные:
начальный момент времени выбираем как t0=0;
конечный момент времени tk=120;
число дискретных отсчётов времени nt=120;
температура стенки Tc=30;
максимальная температура материала в середине цилиндра Tpol=170;
число дискретных отсчетов длинны цилиндра nR=10;
радиус цилиндра R=0.01 м;
температуропроводность полистирола a = 0.00000056 град/м с
Рассчитаем интервалы изменения температуры и радиуса
dr=R/(nR-1);
dt=(tk-t0)/(nt-1);
Присвоим начальные значения температуры стенки в цикле For:
for i=1:nt+1
T(i,1)=Tc;
end
Присвоим начальные значения температурного поля полимера в цикле:
for j=1:nR+2
T(1,j)=Tpol*exp(-2000*(R-dr*(j-1))^2);
end
Рассчитаем матрицу температурного поля Tво вложенном цикле For:
for i=1:nt
for j=1:nR
r=R-dr*(j-1)+0.0001*dr;
T(i+1,j+1)=T(i,j+1)+(a*dt/dr)*(((T(i,j+2)-2*T(i,j+1)+T(i,j))/dr)+((1/r)*(T(i,j+2)-T(i,j+1))));
end
T(i+1,nR+1)=T(i+1,nR);
T(i+1,nR+2)=T(i+1,nR);
end
Изменим порядок расположения столбцов обработав массив в двойном цикле For :
for i=1:nt
for j=1:nR
TT(i,j)=T(i,nR-j+1);
end
end
Построим поверхность описывающую полученную функциональную зависимость T(t,r):
figure(1)
mesh(TT)
Подпишем координатные оси
xlabel('R, MM')
ylabel('t, cek')
zlabel('TC')
6 АНАЛИЗ МОДЕЛИРОВАНИЯ И РАСЧЁТОВ
В результате численного решения дифференциального уравнения с помощью составленной программы получены данные, хорошо согласующиеся с аналитическим решением дифференциального уравнения приведенным во второй главе данной пояснительной записки.
Результаты получаемые с помощью данной программы можно использовать для моделирований реальных технологических процессов связанных с охлаждением и нагреванием цилиндрических каналов.
СПИСОК ИСПОЛЬЗОВАНОЙ ЛИТЕРАТУРЫ
1. Лыков А. В. Теория теплопроводности. М., ГИТТЛ, 1952. 391 с.
2. Карслоу Г., Егер Д. Теплопроводность твердых тел. М., «Наука», 1964. 487 с.
3. Кирпичев М. В., Михеев М. А. Моделирование тепловых устройств. М.,изд-во АН СССР, 1936. 255 с.
4. Тябин Н. В. и др. В кн.: Теплообмен. 1974. Советские исследования. М., «Наука», 1975, с. 195—198.
5. Торнер «Технология переработки пластмасс», Москва, Московский политехи, ин-т, 1965, № 1, с. 138—143.
ПРИЛОЖЕНИЕ1
clear, clc, clf, clgt0=0;
tk=120;
nt=120;
Tc=30;
Tpol=170;
nR=10;
R=0.01;
dr=R/(nR-1);
dt=(tk-t0)/(nt-1);
a=0.00000056;
for i=1:nt+1
T(i,1)=Tc;
end
for j=1:nR+2
T(1,j)=Tpol*exp(-2000*(R-dr*(j-1))^2);
end
for i=1:nt
for j=1:nR
r=R-dr*(j-1)+0.0001*dr;
T(i+1,j+1)=T(i,j+1)+(a*dt/dr)*(((T(i,j+2)-2*T(i,j+1)+T(i,j))/dr)+((1/r)*(T(i,j+2)-T(i,j+1))));
end
T(i+1,nR+1)=T(i+1,nR);
T(i+1,nR+2)=T(i+1,nR);
end
for i=1:nt
for j=1:nR
TT(i,j)=T(i,nR-j+1);
end
end
figure(1)
mesh(TT)
xlabel('R, MM')
ylabel('t, cek')
zlabel('T C')
ПРИЛОЖЕНИЕ2
