Министерство образования Российской Федерации

Сибирский государственный индустриальный университет

Кафедра литейного производства

Расчет затвердевания плоской отливки

в массивной форме

 

Выполнили: ст. гр. МЛА-97

Злобина С. А.

Карпинский А. В.

Кирина Л. В.

Тимаревский А. В.

Токар А. Н.

Проверил: доцент, к.т.н.

Передернин Л.В.

Новокузнецк 2001

Содержание

Содержание. 2

Задание. 3

Постановка задачи. 4

1.   Графическое представление. 4

2.   Математическая формулировка задачи. 5

Метод расчета. 7

Схема апроксимации. 8

Алгоритм расчета. 11

Идентификаторы.. 13

Блок-схема. 14

Программа. 17

Сравнение с инженерными методами расчета. 20

Результаты расчета. 21

Задание

Отливка в виде бесконечной плиты толщиной 2L_o=30 мм

Сплав: Латунь (10% Zn).

Форма: Песчано-глинистая объемная сырая (ПГФ).

Индексы: 1-Ме_тв, 2- Ме_ж, 4-форма.

а₁=3,6×10^-5 м²/с

а₂=2,1×10^-5 м²/с

l₁=195 Вт/м×К

l₂=101 Вт/м×К

r₁=8600 кг/м³

r₂=8000 кг/м³

L=221000 Дж/кг

b₄=1300 Вт×с^1/2/(м²×К)

T_ф=293 К

T_s=1312,5 К

T_н=1345 К

N=100

e_t=0,01 c

e_Т=0,01 ^oC

Постановка задачи

1. 
Графическое представление

Принимаем следующие условия:

Отливка в виде бесконечной плиты толщиной 2Lo затвердевает в объемной массивной песчано-глинистой форме. Принимаем, что теплофизические характеристики формы и металла постоянны и одинаковы по всему объему, системы сосредоточенные, геометрическая ось совпадает тепловой и поэтому можно рассматривать только половину отливки. Lo<<Lф - форма массивная, т.е. форма за все время охлаждения не прогревается до конца, Т_пов=Т_нач; такая форма называется бесконечной

Вектор плотности теплового потока (удельный тепловой поток) имеет направление перпендикулярное к поверхности раздела отливка-форма в любой момент времени t_k;

Нестационарное температурное поле – одномерное, Тj(х, t_k), j=1,2,4;

Температура затвердевания принимается постоянной, равной Ts;

Теплофизические характеристики сред, a_j=l_j/c_jr_j, j=1,2,4;

Теплоаккумулирующую способность формы примем постоянной, b_ф==const;

C,l,r - теплофизические характеристики формы;

Переохлаждение не учитываем;

Удельная теплота кристаллизации L(Дж/кг) выделяется только на фронте затвердевания (nf) - условие Стефана;

Не учитывается диффузия химических элементов – квазиравновесное условие;

Перенос тепла за счет теплопроводности и конвекции учитывается введением коэффициента эффективной электропроводности:

для жидкой среды l₂=n*l₀, где l₀ – теплопроводность неподвижного жидкого металла; n=10;

Не учитывается усадка металла при переходе из жидкого состояния в твердое;

Передача тепла в жидком и твердом металле происходит за счет теплопроводности и описывается законом Фурье:

q = - l_jgradT, плотность теплового потока,Дж/(м²с);

Отливка и форма имеют плотный контакт в период всего процесса затвердевания (что реально для ПГФ);

теплоотдача на границе отливка – форма  подчиняется закону Ньютона(-Рихтмона): q₁(t^k)=a(T₁^к - T^ф) – для каждого момента времени t_к, где a - коэффициент теплоотдачи, для установившегося режима (автомодельного) a=;

Полученная таким образом содержательная модель и ее графическая интерпретация затвердевания плоской отливки в объемной массивной форме, упрощает формулировку математической модели и достаточно хорошо отражает затвердевание на тепловом уровне, т.е. позволяет получить закон T=f(x;t).

2.  Математическая формулировка задачи

Математическая модель формулируется в виде краевой задачи, которая включает следующие положения:

а) Математическое выражение уравнения распределения теплоты в изучаемых средах.

Дифференциальное уравнение теплопроводности Фурье, которое имеет смысл связи, между временным изменением температуры и ее пространственным распределением:

Или в соответствии с условием 5 запишем:

  ; xÌ[0,l_o], j=                                                      (1)

б) Условия однозначности:

1. Теплофизические характеристики сред

r_j, l_j, c_j, b_j, a_j, T_L, T_S

2. Начальные условия

2.1 Считаем, что заливка происходит мгновенно и мгновенно же образуется тончайшая корка твердого металла.

T₁^н(x, t^н)= T_S(E)                                                                                             (2)

2.2 Положение фронта затвердевания

t=t^н_задан. ,x=0, y(t^н)=0                                                                                (3)

2.3 Температура металла в отливке

T_j,i^н=T_н ; j=2, iÌ(2,n)                                                                                (4)

2.4 Температура на внешней поверхности формы (контакт форма - атмосфера) и температура формы.

T₄^н=T_ф                                                                                                                                                                          (5)

3. Граничные условия

3.1 Условия сопряжения на фронте затвердевания (условия Стефана) i=nf

                                                                (6)

3.2 Температура на фронте затвердевания

                                                                      (7)

3.3 Условие теплоотдачи на границе отливка-форма

                                             (8)

 - граничное условие третьего рода

3.4 Условие на оси симметрии

                                                                                   (9)

Задача, сформулированная в выражениях (1-9) есть краевая задача, которая решается численным методом.

Аппроксимировав на сетке методом конечных разностей (МКР), получим дискретное сеточное решение.

T_i=f(x_i;t^k).

Метод расчета

Будем использовать МКР – метод конечных разностей.

Сформулированную краевую задачу дискретизируем на сетке.

        

= - шаг по пространству постоянный; - шаг по времени переменный

Для аппроксимации задачи на выбранной сетке можно использовать разные методы – шаблоны. Наиболее известные из них для данного типа задач четырех точечный конечно разностный шаблон явный и неявный.

Явный четырех точечный шаблон            Неявный четырех точечный шаблон

Использование явного шаблона для каждого временного шага получаем 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)

e_t  и e_T – заданные точности по времени и температуре

e_t=0,01c,   e_T=0,1°C

Dt^I=0,01c – время за которое образовалась корочка.

Описанный итерационный процесс называют ''Ловлей фазового фронта в узел''.

Можно задать Dх, Dt^K=const, тогда неизвестно будет положение фронта, при помощи линейной интерполяции.

Расчет температурных полей:

Метод «прогонки»:

Считается наиболее эффективным для неявно заданных конечно-разностных задач.

Суть метода:

Запишем в общем виде неявно заданное конечноразностное уравнение второго порядка (14) в общем виде:

A_iT_i-1 – B_iT_i + C_iT_i+1 + D_i = 0 ; i = 2, 3, 4, …n-1                                       (27)

действительно для всех j и k.

и краевые условия для него:         

T₁ = p₂T₂ + q₂                                                                                                                                                       (28 а)

T_n = p_nT_m-1 + q_n                                                                                                                                                   (28 б)

T_i = f(A_i; X_i; t_k)   - сеточное решение.

A_i, B_i, C_i, D_i – известные коэффициенты, определенные их условий однозначности и дискретизации задачи.

Решение уравнения (27) – ищем в том же виде, в котором задано краевое условие (28 а)

T_i = а_i+1T_i+1 + b_i+1 ; i = 2, 3, 4, …n-1                                                                                                              (29)

A_i+1, b_i+1 – пока не определенные «прогоночные» коэффициенты (или коэффициенты разностной факторизации)

Запишем уравнение (29) с шагом назад:

T_i-1 = а_iT_i + b_i                                                                                                                                                        (30)

Подставим уравнение (30) в уравнение (27):

A_i(a_iT_i + b_i) – B_iT_i + C_iT_i+1 + D_i = 0

Решение нужно получить в виде (29):

                                                                           (31)

Найдем метод расчета прогоночных коэффициентов.

Сравним уравнение (29) и (31):

                                                                                           (32)

                                                                                             (33)

(32),(33)– рекуррентные прогоночные отношения позволяющие вычислить прогоночные коэффициенты точке  (i+1) если известны их значения в точке i.

Процедура определения коэффициентов а_i+1 и b_i+1 называется прямой прогонкой или прогонкой вперед.

Зная коэффициенты конечных точек и температуру в конечной точке Т_i+1 можно вычислить все Т_i.

Процедура расчета температур называется обратной прогонкой. То есть, чтобы вычислить все Т поля для любого t_k нужно вычислить процедуры прямой и обратной прогонки.

Чтобы определить начальные а₂и b₂, сравним уравнение (29) и уравнение (28 а):

a₂ = p₂;  b₂ = q₂

Запишем уравнение 29 с шагом назад:

T_n = p_nT_n-1 + q_n

T_n-1 = q_nT_n + b_n      

                                                                                          (34)

Новая задача определить p_n , q_n

Вывод расчетных формул:

Преобразуем конечноразностное уравнение (14) в виде (27)

, j=1,2                                       (35)

относиться к моменту времени k

Из (35) => A_i=C_i= B_i=2A_i+ D_i=                                          (36)

Определим значения коэффициентов для граничных условий:

на границе раздела отливка-форма

                                                                                         (37)

приведем это выражение к виду (28 а)

 отсюда                                                                              (38)

b₂=q₂= a₂=p₂=1                                                                                    (39)

на границе раздела Me_тв - Ме_ж

из (29), T_nf=T_n=> a_nf+1=0, b_nf+1=T_s                                                           (40)

условие на оси симметрии

T_n-1=T_n в соответствии с (21)

p_n=1, q_n=0                                                                                                (41)

подставив (41) в (34) получим

                                                                                                  (42)

Алгоритм расчета

1)   Определить теплофизические характеристики сред, участвующих в тепловом взаимодействии λ₁, λ₂, ρ₁, ρ₂, L, а₁, а₂, Т_s, Т_н, Т_ф.

2)   Определить размеры отливки, параметры дискретизации и точность расчета

2l₀=30 мм, l₀=R=15 мм=0,015 м

n=100,

первый шаг по времени: Δt¹=0,01 с, t=t+Δt

е_t=0,01 с, e_t=0,1 ^оC

3)   Принять, что на первом временном шаге к=1, t₁=Δt¹, nf=1, Т₁=Т₃, Т_i=Т_н, , i=2,…,n, Т₄=Т_ф

4)   Величина плотности теплового потока на границе раздела отливка – форма

                                                                       (43)

, s=0, (нулевое приближение)

к=2,                                                                     (44)

5)   Найти нулевое приближение Δt^{к, 0} на к-том шаге

переход  nf → i → i+1 по формуле (23)

6)   Найти коэффициенты A_i, С_i, В_i, D_i по соответствующим формулам для сред Ме_тв. и Ме_ж. В нулевом приближении при s=0

7)   Рассчитать прогоночные коэффициенты a_i+1, b_i+1 для Ме_тв. и Ме_ж., s=0 с учетом что Т_nf=Т_з.

Т₁=р₂Т₂+g₂

Т_i=а₂Т₂+в₂

Найти а₂ и в₂:

а₂=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) и задаем температуру в первой точке

- [Температурное поле для первого шага по времени

- [Делаем шаг по времени

- [Вычисляем плотность теплового потока

- [Шаг по времени в нулевом приближении

- [Начальные прогоночные коэффициенты

- [Шаг по итерации

- [Вычисляем коэффициенты B_j для подстановки в (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_зал=T_L+70, Т_н=1/2(T_зал+Т_L)

Расчет:

с

с

c

Скорость затвердевания во времени характеризуется следующим выражением:

, где u_Е=(Т_Е-Т_ф)