Министерство образования Российской Федерации
Сибирский государственный индустриальный университет
Кафедра литейного производства
Расчет затвердевания плоской отливки
в массивной форме
Выполнили: ст. гр. МЛА-97
Злобина С. А.
Карпинский А. В.
Кирина Л. В.
Тимаревский А. В.
Токар А. Н.
Проверил: доцент, к.т.н.
Передернин Л.В.
Новокузнецк 2001
Графическое представление
Принимаем следующие условия:
Отливка в виде бесконечной плиты толщиной 2Lo затвердевает в объемной массивной
песчано-глинистой форме. Принимаем, что теплофизические характеристики формы и
металла постоянны и одинаковы по всему объему, системы сосредоточенные,
геометрическая ось совпадает тепловой и поэтому можно рассматривать только
половину отливки. Lo<<Lф - форма массивная, т.е. форма за все время
охлаждения не прогревается до конца, Тпов=Тнач; такая
форма называется бесконечной
Вектор плотности теплового потока (удельный тепловой поток) имеет направление
перпендикулярное к поверхности раздела отливка-форма в любой момент времени t
k;
Нестационарное температурное поле – одномерное, Тj(х, tk), j=1,2,4;
Температура затвердевания принимается постоянной, равной Ts;
Теплофизические характеристики сред, aj=lj/cjrj, j=1,2,4;
Теплоаккумулирующую способность формы примем постоянной, bф==const;
C,l,r - теплофизические характеристики формы;
Переохлаждение не учитываем;
Удельная теплота кристаллизации L(Дж/кг) выделяется только на фронте
затвердевания (nf) - условие Стефана;
Не учитывается диффузия химических элементов – квазиравновесное условие;
Перенос тепла за счет теплопроводности и конвекции учитывается введением
коэффициента эффективной электропроводности:
для жидкой среды l2=n*l0, где l0 –
теплопроводность неподвижного жидкого металла; n=10;
Не учитывается усадка металла при переходе из жидкого состояния в твердое;
Передача тепла в жидком и твердом металле происходит за счет теплопроводности
и описывается законом Фурье:
q = - ljgradT, плотность теплового потока,Дж/(м2с);
Отливка и форма имеют плотный контакт в период всего процесса затвердевания
(что реально для ПГФ);
теплоотдача на границе отливка – форма подчиняется закону Ньютона(-Рихтмона): q
1(tk)=a(T1к - Tф) – для
каждого момента времени tк, где a - коэффициент теплоотдачи, для
установившегося режима (автомодельного) a=
;
Полученная таким образом содержательная модель и ее графическая интерпретация
затвердевания плоской отливки в объемной массивной форме, упрощает
формулировку математической модели и достаточно хорошо отражает затвердевание
на тепловом уровне, т.е. позволяет получить закон T=f(x;t).
2. Математическая формулировка задачи
Математическая модель формулируется в виде краевой задачи, которая включает
следующие положения:
а) Математическое выражение уравнения распределения теплоты в изучаемых средах.
Дифференциальное уравнение теплопроводности Фурье, которое имеет смысл связи,
между временным изменением температуры и ее пространственным распределением:
Или в соответствии с условием 5 запишем:
; xÌ[0,lo], j= (1)
б) Условия однозначности:
1. Теплофизические характеристики сред
rj, lj, cj, bj, aj, TL, TS
2. Начальные условия
2.1 Считаем, что заливка происходит мгновенно и мгновенно же образуется
тончайшая корка твердого металла.
T1н(x, tн)= TS(E)
(2)
2.2 Положение фронта затвердевания
t=tнзадан. ,x=0, y(tн)=0
(3)
2.3 Температура металла в отливке
Tj,iн=Tн ; j=2, iÌ(2,n)
(4)
2.4 Температура на внешней поверхности формы (контакт форма - атмосфера) и
температура формы.
T4н=Tф
(5)
3. Граничные условия
3.1 Условия сопряжения на фронте затвердевания (условия Стефана) i=nf
(6)
3.2 Температура на фронте затвердевания
(7)
3.3 Условие теплоотдачи на границе отливка-форма
(8)
- граничное условие третьего рода
3.4 Условие на оси симметрии
(9)
Задача, сформулированная в выражениях (1-9) есть краевая задача, которая
решается численным методом.
Аппроксимировав на сетке методом конечных разностей (МКР), получим дискретное
сеточное решение.
Ti=f(xi;tk).
Метод расчета
Будем использовать МКР – метод конечных разностей.
Сформулированную краевую задачу дискретизируем на сетке.
= - шаг по пространству постоянный; - шаг по времени переменный
Для аппроксимации задачи на выбранной сетке можно использовать разные методы
– шаблоны. Наиболее известные из них для данного типа задач четырех точечный
конечно разностный шаблон явный и неявный.
Явный четырех точечный шаблон Неявный четырех точечный шаблон
Использование явного шаблона для каждого временного шага получаем n+1
уравнение с n неизвестными и система решается методом Гауса, но сходимость
решения только при очень малых шагах.
Использование неявного шаблона обеспечивает абсолютную сходимость, но каждое
из уравнений имеет 3 неизвестных, обычным методом их решить невозможно.
По явному:
(10)
По неявному:
(11)
Сходимость обеспечивается при:
при явном шаблоне
(12)
-точность аппроксимации
(13)
Схема апроксимации
Аппроксимируем задачу 1-9 на четырех точечном неявном шаблоне
Начальные условия: (14)
(15)
(16)
(17)
(18)
Граничные условия: (19)
(20)
(21 a)
=>
(21)
Условие идеального контакта на границе отливка форма
(22)
Расчет временного шага :
Величина -var
рассчитывается из условия, что за промежуток времени
фронт перейдет из точки nf в точку nf+1
Расчет ведут итерационными (пошаговыми) методами
Строим процедуру расчета следующим образом:
Вычисляем нулевое приближенное для каждого шага,
За шаг итерации примем S,
Нулевое приближение S=0.
(23)
Уточняем шаг: S+1
(24)
d – параметр итерации от 0 до 1
для расчета возьмем d=0.
Число S итераций определяется заданной точностью:
Временного шага
(25)
И по температуре (26)
et и eT – заданные точности по времени и температуре
et=0,01c, eT=0,1°C
DtI=0,01c – время за которое образовалась корочка.
Описанный итерационный процесс называют ''Ловлей фазового фронта в узел''.
Можно задать Dх, DtK=const, тогда неизвестно будет положение фронта,
при помощи линейной интерполяции.
Расчет температурных полей:
Метод «прогонки»:
Считается наиболее эффективным для неявно заданных конечно-разностных задач.
Суть метода:
Запишем в общем виде неявно заданное конечноразностное уравнение второго
порядка (14) в общем виде:
AiTi-1 – BiTi + CiT
i+1 + Di = 0 ; i = 2, 3, 4, .n-1
(27)
действительно для всех j и k.
и краевые условия для него:
T1 = p2T2 + q2
(28 а)
Tn = pnTm-1 + qn
(28 б)
Ti = f(Ai; Xi; tk) - сеточное решение.
Ai, Bi, Ci, Di – известные
коэффициенты, определенные их условий однозначности и дискретизации задачи.
Решение уравнения (27) – ищем в том же виде, в котором задано краевое условие
(28 а)
Ti = аi+1Ti+1 + bi+1 ; i = 2, 3, 4, .n-1 (29)
Ai+1, bi+1 – пока не определенные
«прогоночные» коэффициенты (или коэффициенты разностной факторизации)
Запишем уравнение (29) с шагом назад:
Ti-1 = аiTi + bi
(30)
Подставим уравнение (30) в уравнение (27):
Ai(aiTi + bi) – BiTi + CiTi+1 + Di = 0
Решение нужно получить в виде (29):
(31)
Найдем метод расчета прогоночных коэффициентов.
Сравним уравнение (29) и (31):
(32)
(33)
(32),(33)– рекуррентные прогоночные отношения позволяющие вычислить
прогоночные коэффициенты точке (i+1) если известны их значения в точке i.
Процедура определения коэффициентов аi+1 и bi+1
называется прямой прогонкой или прогонкой вперед.
Зная коэффициенты конечных точек и температуру в конечной точке Тi+1
можно вычислить все Тi.
Процедура расчета температур называется обратной прогонкой. То есть, чтобы
вычислить все Т поля для любого tk нужно вычислить процедуры прямой
и обратной прогонки.
Чтобы определить начальные а2и b2, сравним уравнение (29) и уравнение (28 а):
a2 = p2; b2 = q2
Запишем уравнение 29 с шагом назад:
Tn = pnTn-1 + qn
Tn-1 = qnTn + bn
(34)
Новая задача определить pn , qnВывод расчетных формул:
Преобразуем конечноразностное уравнение (14) в виде (27)
, j=1,2 (35)
относиться к моменту времени k
Из (35) => Ai=Ci= Bi=2Ai+ Di= (36)
Определим значения коэффициентов для граничных условий:
на границе раздела отливка-форма
(37)
приведем это выражение к виду (28 а)
отсюда
(38)
b2=q2=
a2=p2=1
(39)
на границе раздела Meтв - Меж
из (29), Tnf=Tn=> anf+1=0, bnf+1
=Ts (40)
условие на оси симметрии
Tn-1=Tn в соответствии с (21)
pn=1, qn=0
(41)
подставив (41) в (34) получим
(42)
Алгоритм расчета
1) Определить теплофизические характеристики сред, участвующих в тепловом
взаимодействии λ1, λ2, ρ1,
ρ2, L, а1, а2, Тs, Тн
, Тф.
2) Определить размеры отливки, параметры дискретизации и точность расчета
2l0=30 мм, l0=R=15 мм=0,015 м
n=100,
первый шаг по времени: Δt1=0,01 с, t=t+Δt
еt=0,01 с, et=0,1 оC
3) Принять, что на первом временном шаге к=1, t1=Δt1
, nf=1, Т1=Т3, Тi=Тн, , i=2,.,n, Т
4=Тф
4) Величина плотности теплового потока на границе раздела отливка – форма
(43)
, s=0, (нулевое приближение)
к=2, (44)
5) Найти нулевое приближение Δtк, 0 на к-том шаге
переход nf → i → i+1 по формуле (23)
6) Найти коэффициенты Ai, Сi, Вi, Di
по соответствующим формулам для сред Метв. и Меж. В
нулевом приближении при s=0
7) Рассчитать прогоночные коэффициенты ai+1, bi+1
для Метв. и Меж., s=0 с учетом что Тnf=Тз.
Т1=р2Т2+g2
Тi=а2Т2+в2
Найти а2 и в2:
а2=1, (45)
(46)
8) Рассчитать температуру на оси симметрии
(47)
9) Рассчитать температурное поле жидкого и твердого металла
(48)
10) Пересчитать значения ∆tк по итерационному процессу (24)
d – параметр итерации (d=0.1)
проверяем точность;
11) Скорость охлаждения в каждом узле i рассчитать по формуле:
, оС/с
(50)
12) Скорость затвердевания на каждом временном шаге:
, м/с
(51)
13) Средняя скорость охлаждения на оси отливки:
14) Положение фронта затвердевания по отношению к поверхности отливки
, к – шаг по времени (52)
15) Полное время затвердевания
, к′ -
последний шаг (53)
16) Средняя скорость затвердевания отливки
(54)
Идентификаторы
Блок-схема
- [Вводим исходные данные
- [Вычисляем шаг по пространству
- [Вычисляем коэффициенты Аj, Сj для подстановки в (32),
(33) и задаем температуру в первой точке
- [Температурное поле для первого шага по времени
- [Делаем шаг по времени
- [Вычисляем плотность теплового потока
- [Шаг по времени в нулевом приближении
- [Начальные прогоночные коэффициенты
- [Шаг по итерации
- [Вычисляем коэффициенты Bj для подстановки в (32), (33)
- [Вычисляем прогоночные коэффициенты по твердому металлу
- [Прогоночные коэффициенты для фронта
- [Вычисляем прогоночные коэффициенты по жидкому металлу
- [Температура на оси симметрии
- [Расчет температурного поля
- [Ищем максимальный температурный шаг
- [Уточняем Dt
- [Точность временного шага
- [Проверка точности
- [Расчет времени
- [Скорость охлаждения в каждом узле
- [Скорость затвердевания и положение фронта
- [Вывод результатов
- [Проверка достижения фронтом центра отливки
- [Расчет полного времени, ср. скорости затвердевания ср. скорости охлаждения
на оси отливки
Вывод результатов
- [Конец.
Программа
CLEAR , , 2000
DIM T(1000), T1(1000), AP(1000), BP(1000), Vox(1000), N$(50)
2 CLS
N = 100: KV = 50: N9 = 5: L = .015
TM = 293: TI = 1345: TS = 1312.5
BM = 1300: a1 = .000036: a2 = .000021
TA0 = .01: ETA = .01: E = .01
l1 = 195: l2 = 101
R0 = 8600: LS = 221000
AF = 0: Pi = 3.14159265359#
3 PRINT "Число шагов N, штук"; N
PRINT "Длина отливки L, м"; L
PRINT "Температура формы Tf, К"; TM
PRINT "Начальная температура сплава Tн, К"; TI
PRINT "Температура затвердевания Tz, К"; TS
PRINT "Bф "; BM
PRINT "Первый шаг по времени, Tk0 "; TA0
PRINT "Точность по времени, Еt "; ETA
PRINT "Точность по температуре, ЕТ "; E
PRINT "Температуропроводность Ме твердого, а1 "; a1
PRINT "Температуропроводность Ме жидкого, а2 "; a2
PRINT "LS= "; LS
PRINT "Коэф. теплопроводности, l1 "; l1
PRINT "Коэф. теплопроводности, l2"; l2
PRINT "Плотность Ме твердого, р1 "; R0
INPUT "Изменить данные <y/n>"; QV$
IF QV$ = "Y" THEN GOSUB 222
48 N1 = N - 1
DX = L / (N - 1)
A = a1 / DX ^ 2
B1 = 2 * A
RL = R0 * LS * DX
NF = 1
B2 = l1 / DX
KV1 = 1
AL = a2 / DX ^ 2
BL1 = 2 * AL
BL2 = l2 / DX
T(1) = TS
T1(1) = TS
FOR i = 2 TO N
T(i) = TI
T1(i) = TI
NEXT i
TA = TA0
K = 1
dta = .01
GOTO 103
101 K = K + 1
NF = NF + 1
B3 = SQR(Pi * TA)
q = BM * (T(1) - TM) / B3
dta = RL / (AF + q)
B5 = BM * TM / B3
B3 = BM / B3
B4 = B2 + B3
AP(1) = B2 / B4
BP(1) = B5 / B4
T(NF) = TS
NF1 = NF - 1
NF2 = NF + 1
K1 = 0
102 K1 = K1 + 1
Et = 0
B3 = SQR(Pi * (TA + dta))
q = BM * (T(1) - TM) / B3
B5 = BM * TM / B3
B3 = BM / B3
B4 = B2 + B3
AP(1) = B2 / B4
BP(1) = B5 / B4
DTA1 = 1 / dta
IF NF1 = 1 THEN GOTO 23
FOR i = 2 TO NF1
B = B1 + DTA1
f = DTA1 * T1(i)
B4 = B - A * AP(i - 1)
AP(i) = A / B4
BP(i) = (A * BP(i - 1) + f) / B4
NEXT i
23 FOR i = NF1 TO 1 STEP -1
TC = AP(i) * T(i + 1) + BP(i)
B = ABS(TC - T(i)) / TC
IF B > Et THEN Et = B
T(i) = TC
NEXT i
AP(NF) = 0
BP(NF) = TS
B = BL1 + DTA1
FOR i = NF2 TO N
f = DTA1 * T1(i)
B4 = B - AL * AP(i - 1)
AP(i) = AL / B4
BP(i) = (AL * BP(i - 1) + f) / B4
NEXT i
IF NF = N THEN GOTO 34
TC = BP(N) / (1 - AP(N))
B = ABS(TC - T(N)) / TC
T(N) = TC
IF B > Et THEN Et = B
IF NF >= N1 THEN GOTO 34
FOR i = N1 TO NF2 STEP -1
TC = AP(i) * T(i + 1) + BP(i)
B = ABS(TC - T(i)) / TC
IF B > Et THEN Et = B
T(i) = TC
NEXT i
34 P = AF + q
P1 = 1 / P
TM2 = BL2 * (T(NF2) - TS)
IF NF = N THEN GOTO 80
TM1 = B2 * (TS - T(NF1))
DTF = P1 * (RL + dta * (TM2 - TM1 + P))
P3 = ABS(DTF - dta) / DTF
dta = DTF
IF (P3 > ETA) OR (Et > E) THEN GOTO 102
80 TA = TA + dta
IF NF = 1 THEN dta = TA0
Vox = (T1(NF) - TS) / dta
FOR i = 1 TO N
Vox(i) = (T1(i) - T(i)) / dta
T1(i) = T(i)
NEXT i
VS = DX / dta
Xf = (K - 1) * DX
IF K <> KV1 + 1 THEN GOTO 33
KV1 = KV1 + KV
GOSUB 777
33 GOTO 105
103 PRINT "РЕЗУЛЬТАТЫ РАСЧЕТА": CLS : GOSUB 777
105 IF K < N THEN GOTO 101
GOSUB 777
Vz = 1000 * L / TA
Voxl = (TI - TS) / TA
PRINT "Полное время затв. отл. TA="; TA; "с."
PRINT "Ср. скорость охл. на оси отл. Voxl="; Voxl; " K/с"
PRINT "Ср. скорость затв. отл. Vz="; Vz; " мм/с"
END
777 PRINT "К="; K; " DTA="; dta; "VS="; VS * 1000; " мм/с XF="; Xf; " мм"
PRINT "T="; T(1); : FOR i = 1 TO 10: PRINT T(i * 10); : NEXT i: PRINT "K"
PRINT "Vox="; Vox(1); : FOR i = 1 TO 10: PRINT Vox(i * 10); : NEXT i: PRINT
"K/c"
RETURN
222 CLS
INPUT "Число шагов N, штук"; N
INPUT "Длина отливки L, м"; L
INPUT "Температура формы Tf, К"; TM
INPUT "Начальная температура сплава Tн, К"; TI
INPUT "Температура затвердевания Tz, К"; TS
INPUT "Bф "; BM
INPUT "Первый шаг по времени, Tk0 "; TA0
INPUT "Точность по времени, Еt "; ETA
INPUT "Точность по температуре, ЕТ "; E
INPUT "Температуропроводность Ме твердого, а1 "; a1
INPUT "Температуропроводность Ме жидкого, а2 "; a2
INPUT "LS= "; LS
INPUT "Коэф. теплопроводности, l1 "; l1
INPUT "Коэф. теплопроводности, l2"; l2
INPUT "Плотность Ме твердого, р1 "; R0
CLS
GOTO 3
RETURN
Сравнение с инженерными методами расчета
Г. Ф. Баландин для расчета продолжительности затвердевания отливки
эвтектического сплава предложил следующие выражения:
-время заливки
-время снятия перегрева
-время затвердевания
Принимаем Tзал=TL+70, Тн=1/2(Tзал+ТL)
Расчет:
с
с
c
Скорость затвердевания во времени характеризуется следующим выражением:
, где uЕ=(ТЕ-Тф)