Реферат: Расчет затвердевания плоской отливки
Министерство образования РоссийскойФедерации
Сибирский государственный индустриальный университет
Кафедра литейного производства
Расчет затвердевания плоскойотливки
в массивной форме
Выполнили: ст. гр. МЛА-97
ЗлобинаС. А.
КарпинскийА. В.
КиринаЛ. В.
ТимаревскийА. В.
ТокарА. Н.
Проверил: доцент, к.т.н.
Передернин Л.В.
Новокузнецк2001
СодержаниеСодержание. 2
Задание. 3
Постановказадачи. 4
1. Графическое представление. 4
2. Математическая формулировка задачи. 5
Методрасчета. 7
Схемаапроксимации. 8
Алгоритмрасчета. 11
Идентификаторы… 13
Блок-схема. 14
Программа. 17
Сравнениес инженерными методами расчета. 20
Результатырасчета. 21
/>/>ЗаданиеОтливка в виде бесконечной плиты толщиной2Lo=30 мм
Сплав: Латунь (10% Zn).
Форма: Песчано-глинистая объемная сырая(ПГФ).
Индексы: 1-Метв, 2- Меж,4-форма.
а1=3,6×10-5 м2/с
а2=2,1×10-5 м2/с
l1=195Вт/м×К
l2=101Вт/м×К
r1=8600кг/м3
r2=8000кг/м3
L=221000 Дж/кг
b4=1300Вт×с1/2/(м2×К)
Tф=293К
Ts=1312,5К
Tн=1345К
N=100
et=0,01 c
eТ=0,01 oC
/>/>Постановка задачи1. /> <td/> />Графическое представление
Принимаем следующие условия:
Отливка в виде бесконечной плиты толщиной2Lo затвердевает в объемной массивной песчано-глинистойформе. Принимаем, что теплофизические характеристики формы и металла постоянныи одинаковы по всему объему, системы сосредоточенные, геометрическая ось совпадаеттепловой и поэтому можно рассматривать только половину отливки. Lo<<Lф — форма массивная, т.е.форма за все время охлаждения не прогревается до конца, Тпов=Тнач;такая форма называется бесконечной
Вектор плотности теплового потока(удельный тепловой поток) имеет направление перпендикулярное к поверхностираздела отливка-форма в любой момент времени tk;
Нестационарное температурное поле –одномерное, Т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с);
Отливка и форма имеют плотный контакт впериод всего процесса затвердевания (что реально для ПГФ);
теплоотдача на границе отливка – форма подчиняется закону Ньютона(-Рихтмона): q1(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).
/>/>Метод расчетаБудем использовать МКР – метод конечныхразностей.
Сформулированную краевую задачудискретизируем на сетке.
/> />
/>=/> - шаг по пространствупостоянный;/> - шаг по времени переменный
Для аппроксимации задачи на выбраннойсетке можно использовать разные методы – шаблоны. Наиболее известные из них дляданного типа задач четырех точечный конечно разностный шаблон явный и неявный.
Явный четырех точечныйшаблон Неявный четырех точечный шаблон
/> <td/> />Использование явного шаблона для каждого временного шага получаем n+1 уравнение с n неизвестными исистема решается методом Гауса, но сходимость решения только при очень малыхшагах.
Использование неявного шаблонаобеспечивает абсолютную сходимость, но каждое из уравнений имеет 3 неизвестных,обычным методом их решить невозможно.
По явному:
/> (10)
По неявному:
/> (11)
Сходимость обеспечивается при:
/>при явном шаблоне (12)
/>-точностьаппроксимации
/> (13)
/>/>Схема апроксимацииАппроксимируем задачу 1-9 на четырехточечном неявном шаблоне
Начальные условия:
/> (14)
/> (15)
/> /> (16)
/> /> (17)
/> (18)
Граничные условия:
/> (19)
/> (20)
/> (21a)
/>=> /> (21)
Условие идеального контакта на границеотливка форма
/> (22)
Расчет временного шага />:
Величина />-var рассчитывается из условия, что за промежуток времени /> фронт перейдет из точки nf в точку nf+1
Расчет ведут итерационными (пошаговыми)методами
Строим процедуру расчета следующимобразом:
Вычисляем нулевое приближенное />для каждого шага,
За шаг итерации примем S,
Нулевое приближение S=0.
/> (23)
Уточняем шаг: S+1
/> (24)
d – параметритерации от 0 до 1
для расчета возьмем d=0.
Число S итерацийопределяется заданной точностью:
Временного шага/> (25)
И по температуре/> (26)
et<sub/> и eT – заданныеточности по времени и температуре
et=0,01c, eT=0,1°C
DtI=0,01c – время за которое образовалась корочка.
Описанный итерационный процесс называют''Ловлей фазового фронта в узел''.
Можно задать Dх, DtK=const, тогда неизвестно будет положение фронта, при помощилинейной интерполяции.
Расчет температурных полей:
Метод «прогонки»:
Считается наиболее эффективным для неявнозаданных конечно-разностных задач.
Суть метода:
Запишем в общем виде неявно заданноеконечноразностное уравнение второго порядка (14) в общем виде:
AiTi-1– BiTi + CiTi+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<sub/>всоответствии с (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 GOTO34
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 + 1THEN GOTO 33
KV1 = KV1 + KV
GOSUB 777
33 GOTO 105
103 PRINT «РЕЗУЛЬТАТЫ РАСЧЕТА»: CLS:GOSUB 777
105 IF K < N THEN GOTO101
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Е=(ТЕ-Тф)
/>