автореферат диссертации по информатике, вычислительной технике и управлению, 05.13.16, диссертация на тему:Численное решение нестационарных теплофизических задач с фазовым переходом в интервале температур

Гб ОА

1998

российская академия паук сибирское отделе!iii!: институт вычислительной математики и математической геофизики

На правах рукописи

ПОПОВ Владимир Николаевич

ЧИСЛЕННОЕ РЕШЕНИЕ НЕСТАЦИОНАРНЫХ ТЕПЛОФИЗИЧЕСКИХ ЗАДАЧ С ФАЗОВЫМ ПЕРЕХОДОМ В ИНТЕРВАЛЕ ТЕМПЕРАТУР

Специальность 05.13.16 - применение вычислительной техники, математического моделирования и математических методов в научных исследованиях

АВТОРЕФЕРАТ диссертации на соискание ученой степени доктора физико-математических наук

ц ДЕК

Новосибирск - 1998

Работа выполнена в Институте теоретической и прикладной механики СО РАН

Научные консультанты:

доктор физико-математических наук, пофессор В.П. Ильин, доктор физико-математических наук, с.н.с. А.Н. Черепанов

Официальные оппоненты:

доктор физико-математических наук В.И. Дробышевич доктор физико-математических наук Ю.Е. Воскобойников доктор технических наук В.В. Саломатов

Ведущая организация:

Институт физики прочности и материаловедения СО РАН (г.Томск)

Защита состоится января 1998 г. в 15 ч. 00 мин. на заседанк специализированного совета Д 002.10.02 при Институт вычислительной математики и математической геофизики С РАН по адресу: 630090, Новосибирск-90, пр. ак. Лаврентьева, 6

С диссертацией можно ознакомиться в читальном за)

ИВМиМГ СО РАН.

Автореферат разослан 1998

г.

Ученый секретарь специализированного совета кандидат технических наук

Г.И.Забиняко

I. ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность темы. Возрастающее требования к современной тех->логии получения металла выдвигают задачи исследования процессов [твердевания. Метода моделирования с использованием быстродейст--вдих ЭВМ часто оказывается единственным способом определения :тимальных режимов, при которых происходит тот или иной процесс» к как возможности экспериментального изучения этих явлений огра-чены методическими к техническими трудностями.

Теория магматогенных гидротермальных рудных месторождений язана с исследованием динамики дифференциации (перераспределе-я) магмы в камерах земной коры. С момента возникновения пэтроло-и развивается и совершенствуется модельное прогнозирование при иске руд, пространственно сопряженных с проявлением маги в форме трузивных массивов и вулканов. Решение проблемы количественного исания процессов перераспределения пока продвинулось мало, поэму представляет теоретический и прикладной интерес разработка едельных количественных моделей дифференциации магм и получение ответствупцих оценок.

Большое количество задач в металлургии и петрологии, при мольном исследовании которых используются уравнения теплоиереноса, лают актуальным поиск эффективных численных алгоритмов решения.

Целью диссертационной работы является: разработка и развитие математических моделей нестационарных тешюфизических процессов, сопровождавшихся затвердеванием расплавов в интервале температур;

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

изучение итерационных алгоритмов решения линейных и нелинейных алгебраических систем, получаемых при неявных аппроксимациях уравнений тенлопереноса;

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

моделирование петрологических процессов и проведение вычислительных экспериментов по исследованию формирования магматических

интрузивов на основе решения нестационарной модели.

Методика исследования. Рассмотрение разностных алгоритмов, получение условий сходимости и оценок погрешности решения в равномерной норме проводилось с использованием свойств монотонных матриц и при помощи методических расчетов модельных задач. Изучение процессов теплопереноса в практических и исследовательских задачах основывалось на использовании результатов численного моделирования

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

- исследованы неявные разностные алгоритмы численного решения одномерных и двухмерных задач затвердевания с фазовым переходом в интервале температур, описываемых квазилинейными параболическими уравнениями с краевыми условиями 1-Ш родов, и на основе свойств Ы-матриц получена оценка погрешности решения в равномерной норме;

- построены разностные схемы повышенного порядка точности для многомерных уравнений теплопроводности и определены соотношения между пространственными и временными шагами, при которых они являются монотонными;

- рассмотрена эффективность применения итерационного метода неполной факторизации с ускорением по методу сопряженных градиентов для решения линейных и нелинейных систем, получаемых при неявной аппроксимации с различным порядком точности многомерных задач теплопереноса;

- проведено -численное моделирование процесса получения круглого стального слитка в машине непрерывного литья и определена зависимость, связывающая значения коэффициентов теплоотдачи в секциях зоны вторичного охлаждения с удельным расходом вода установки

- решена практическая задача исследования процессов теплопереноса в изложнице сложной геометрической формы и получена. оценка возникающих в конструкции термонапряжений;

- проведено численное моделирование процессов интрузии и определена параметра, характеризующие режима формирования магматических систем.

Научная и практическая ценность. Диссертация посвящена разработке и применен«» математических моделей тешюфизичесних процессов в металлургии и петрологии, связанных с затвердеванием раснла-

юв, а также развитию математического аппарата решения соответст-гувщих задач. Выполненные численные эксперименты при моделировании [роцессов теплопереноса и структурообразовапия служат в качестве >боснования достоверности теоретических исследований. Проведено апробирование метода неволной факторизации с ускорением по методу »пряженных градиентов для решения систем уравнений, получаемых [ри неявной аппроксимации модельных задач. Полученные в диссерта-[ии результаты нашли практическое применение в усовершенствовании ¡уществувдих и разработке новых технологических процессов на йпадно-Сибирском металлургическом комбинате (г. Новокузнецк), [овосибирском заводе хммконцентратов (г. Новосибирск). Математи-:еская модель процесса формирования интрузивов используется в 1б1,единенном институте геологии, геофизики и минералогии СО РАН г. Новосибирск) в виде пакета программ для ЭВМ.

Апробация работы. Основные результаты работы докладывались ■ на Международных конференциях "Численные метода и приложения" София, 1989), "Высокоазотистые стали - 89" (Варна, 1989), "При-шение компьютеров в литейном производстве" (София, 1990), "Крис-аллизация и компьютерные модели" (Ижевск, 1991, 1992, 1994); на 111-ом Российском совещании по экспериментальной минералогии Черноголовка, 1995); обсуждались на семинарах Института вычисли-елъной математики и математической геофизики и Института теорети-еской и прикладной механики СО РАН.

Публикации и личный вклад автора. Основные результаты, издоенные в диссертации, опубликованы в работах 11-201. Все теорети-еские, экспериментальные и расчетные исследования, вошедшие в аботу, выполнены лично автором.

Структура и объем работы. Диссертация состоит из введения, рех глав, заключения и списка цитируемой литературы. Работа изло-ена на 232 страницах и включает 39 рисунков, 26 таблиц, список итераторы из 192 наименований.

II. СОДЕРЖАНКЕ РАБОТЫ

Во введении дается общая характеристика работы, определящая е актуальность, формулируется задача исследования - математичее-ое моделирование нестационарных теплофизических процессов» сопро-ождаицихся затвердеванием расплавов в интервале температур, крат-о излагается структура диссертации и ее содержание по главам.

з

Процесс затвердевания, связанный с тешюотводом, определяется теплофизическими свойствами расплава: теплопроводиостью. теплоемкость!) и скрытой теплотой кристаллизации. Скрытая теплота кристаллизации физически соответствует энергии, которую необходимо отвести от металла при переходе от неупорядоченной структуры жидкости к упорядоченной структуре твердой фазы.

Традиционная постановка задачи о затвердевании жидкости, сформулированная австрийским математиком Стефаном более ста лет назад, основана на представлении о том, что: а) кристаллизация происходит на достаточно гладкой границе раздела фаз; б) фазовый переход происходит при постоянной температуре; в) скорость движения границы раздела фаз определяется условием отвода скрытой теплоты фазового перехода. Однако в настоящее время приближение Стефана используется при решении узкого круга задач: затвердевание чистых металлов, рост монокристаллов, кристаллизация с гладкой границей фазового перехода в контролируемых условиях.

Изучение динамики затвердевания сплавов при использовании различных экспериментальных методов позволило сделать вывод о том, что фронт кристаллизации существует в виде пространственно-распределенной двухфазной зоны и лишь в очень редких случаях как совершенно гладкая поверхность раздела фаз. Затвердевание сплавов происходит в интервале температур ликвидуса и солидуса ?5. Скрытая теплота кристаллизации выделяется в пределах двухфазной зоны.

В последнее время для решения задач кристаллизации металлических сплавов широко и наиболее успешно используется теория "квазиравновесной" двухфазной зоны» предложенная и развитая В.Т. Борисовым, согласно которой затвердевание полностью завершается при температуре, соответствующей особой точке диаграммы состояния. Чаще всего ?5=Гц, ТЕ - где температура эвтектики. При этом сечение жидкой фазы меняется скачком от значения до нуля. При численном моделировании окончание затвердевания определяется заданием величины /гв=о.025-0.05 в конце двухфазной зоны.

В предположении того, что является функцией температуры и вполне определяется видом диаграммы состояния, в пределах двухфазной зоны уравнение теплопереноса имеет следующий вид

дТ э/.

-л7 ^^двгаЩ*)!-«^ —1 . и.у.а,)^.

где зед - удельная скрытая теплота кристаллизации, р£ - плотность твердой фазы материала. Последнее слагаемое отлично от нуля только в интервале температур кристаллизации. Индекс 25 характеризует теплофизические параметры в области После использования подстановки

°Л - 8Л

34 дТ 31

уравнение теплопроводности в двухфазной зоне переписывается в виде

81

где сЭф называется эффективной удельной теплоемкостью сплава в двухфазном его состоянии и определяется выражением

«ЧГ - -

а для (рс)г

(рс)15=(рс)г/1+(рс)3(1~/г)

сэф(7)=(рс)гя+3еор5 — ,

обычно используются выражения вида

где индексы I, б характеризует теплофизические параметры в областях жидкого и твердого состояний вещества соответственно.

Таким образом, учет выделения тепла кристаллизации в двухфазной зоне сводится к соответствующему заданию зависимости эффективной удельной теплоемкости сплава от температуры, а уравнение теп-лопереноса приобретает вид, характерный для квазилинейных параболических уравнений.

На границе жидкой и двухфазных зон т)^ выполняется условия

-О'

Т^+0'

д?

Х,(?) — 1 дп

дТ да

а на границе двухфазной зоны и твердой фазы г/,

д1

Т

ть-<Г

справедливо дТ

т, -о~МП — "я и 5 зп

V0'

Далее описав процесс теплопереноса в области перегретого расплава Пг и ого затвердевшей части 05, задавая начальные и граничные условия получаем замкнутую систему для определения закона

затвердевания и изменения температурных полей.

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

В настоящее время при решении задач теплопереноса широко используется метод конечных элементов, вместе с тем конечно-разностные алгоритмы сохраняют свою привлекательность и находят широкое применение. При практических расчетах используются, как правило, наиболее известные разностные алгоритмы численного решения линейных или нелинейных задач теплопереноса, существенный вклад в развитие которых внесли. A.A. Самарский, Г.И. Марчук, H.H. Яненко, J. Douglas. Основными показателями качества разностного уравнения являются устойчивость и сходимость» порядок погрешности, экономичность. Однако теоретическое исследование разностных алгоритмов для решения нелинейных задач приводит к громоздким вычислениям, а получаемые оценки весьма грубы и. дают не вполне правильное представление об условиях применимости рассматриваемых схем.

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

Параграф I.I посвящен исследованию алгоритма решения одномерной задачи затвердевания металла в форме, математическая модель которой сводится к нелинейному параболическому уравнению с разрывными коэффициентами и включает основные сильноменявдиеся тешюфи-зические параметры - эффективную теплоемкость и коэффициент теплоотдачи. Процесс кристаллизации в двухфазной зоне описывается с использованием равновесного приближения. Задача формулируется в следухщем виде:

ду д2у

<р(у> — = й. —0$r^rn-Q, dt 1 аг2 0

ау

= а„

<Г У

01 с «Эх"

"1

дх Щ (СМ)

дх <9у(1,г) дх

= К

Зу(10+-0,1) дх

(^[у (Х^О,* >-у<Х0-0,г ) ] ,

о.

о<г <т.

у(х,0) = У10(х), 0$хао-0, у(х,0) = У20(х), х0+СКт$1.

Индексы 1 и 2 относятся к разным средам - отливке и изложнице соответственно; теплофизические параметры а^, , а1, а2 -положительные константы, а ус - температура внешней среда. Функция моделирует эффективную теплоемкость с учетом выделения тепла при кристаллизации, а функция а(у}>0 определяется в соответствии с законом Стефана-Волъимана- Предполагается, что функции У» ф(у), а(у) - непрерывны, ограничены, обладат достаточной гладкостью, и для них выполняется условия

ОУ ОФ1V)

* й2*

дг

а и

дц

д<х{*с\)

дц

-«><7)<со.

где М1 > 0, 1=1,2,3.

Разностная аппроксимация задачи проводится интегро-интерполя-ционным методом на равномерной сетке в области изменения х с шагом Л таким, чтобы точка хо попала в узел. В результате получается нелинейная система разностных уравнений, аппроксимирующая задачу с порядком 0{1+Ъ?), которая записывается в виде векторного уравнения

В

ип-ип~1

■о

+ - ВЭД1 =

1

ТА0П+

1

У31.

Здесь 13 - вектор с компонентами, соответствующими значениям решения системы в узлах сетки в момент времени гп; А - трехдиатональ-ная, а и Т)7^ - диагональные матрицы; вектор.

Для решения системы нелинейных уравнений применяется метод простой итерации с использованием в качестве начального приблюке-

л

%

1Ш1 значения с предыдущего временного шага:

цП.З-ИпП-Ч 1 « 1

В^.з - + _ дп.зуп.а+т = _ Аип,з+1+ _-уп,а

и % к 1 Нг К

„П.О = уП-1^

который сходится при выполнении условия

1 < 1тахШ1М2,К5)Г1,

Отсвда, если при достаточно больших Бп

рп.зп^п.вг.-!^ ^ е

и ап<в(И1М2)-1, 0<9<1» то для вектора погрешности разностного решения гп в каждый момент времени (Х^? справедливо

где ЬгИ^а-бГ1, 1®|с=1®1с+мчй2е» погрешность аппроксимации

и |в|с=0(1.Л2).

Принимая во внимание равенство нулю начальной ошибки и за счет выбора е (например, при условии £=Иг), можно получить следующую оценку

¡2% « с0|®§с,

где «§вс=0{1+?I2), С0=ТетрСТЬ).

В § 1.2 рассматривается неявные алгоритмы для решения задач описываемых многомерными уравнениями теплопроводности. В первом разделе параграфа, на основе свойств монотонных матриц, проводится исследование неявной разностной схемы, используемой для решения задачи в следующей постановке

ди <3, ди . д г ди . с {и) — = — Ми) — + — Х(и) — |-нр(и), О^хО, 0<у<1, 51 дх дх ) ду 1 ду }

ди

— = о, х=0,

дх

ди

— = а^иИи -и), х^, 0<у<.1, дх с

и(г,у,0)=иои,у>, 0<у<1, и(.1,0,0)=^, и(х,1,0)=и2(О), СХх$1, 1=0.

Функции и, с(и), \(и), ф(и>, а(и) - непрерывны, обладает лопаточной гладкостью, и для них выполняется условия .

0<с1<с(и)<а>, Ос^сАДи )<оо, |ф(и)Кф0<™, а(и)> 0,

¡си|. 1^1,^Мсу. ¡О • ¡VI

'де с,, Я.,, ф0> q, И - константы.

Аппроксимация балансных соотношений, получаемых после интег-(ирования исходных уравнений, проводится на пятиточечном шаблоне : шагами кх, по пространственным переменным. Полученная нели-юйяая система неявных разностных уравнений записывается в век-тарном виде

Вц —- + + = Р11.

5десь Vй - вектор с компонентами, соответствующими значениям решетя системы в узлах сетки в момент времени - диагональная атрица; А1 - блочно-даагональная матрица, каждый блок которой )сть квадратная трехдиагональная матрица; Л^ - блочно-трехдиаго-шльная матрица, где блоки - диагональные матрицы; Рп - вектор. Полученное уравнение приводится к виду

"де матрица [ВиШ^1 (А"+А^)] имеет строгое диагональное преобладаем и является монотонной при любых соотношениях 1, Лх, йу и Из этого следует, что при неотрицательных [ачальных данных и правой части, значения решений на всех времен-шх шагах неотрицательны. При выполнении условий

шг^.уос, /М2, Шшах)На^) (ишах-ис)^0.

а<0с1/Ш(1+5Ы)], Ь=М(1+5М)/1с1(1-6)],

для вектора погрешности разностного решения времени справедлива оценка

где 1®|с = 0<и-й|+?ф.

Дальнейшие части § 1.2 посвящены рассмотрению "компактных." разностных схем повышенного порядка точности (с использованием узлов ближайших к рассматриваемому узлу сетки) для решения много-гомерных уравнений параболического типа. В разделах параграфа исследуются способы аппроксимаций повышенной точности для уравнений теплопереноса, а также анализируются условия на соотношения между пространственными и временными шагами, обеспечивающие монотонность разностных схем. Наряду с этим рассматриваются аппроксимации порядка для уравнения теплопроводности в цилиндрической системе координат, которое имеет свои особенности на оси симметрии.

В случае использования неявных разностных схем при решении нестацонарных многомерных задач тешхопереноса возникает проблема реализации расчетов на каждом временном шаге. К наиболее эффективным методам решения алгебраических систем высокого порядка относятся итерационные алгоритмы неполной факторизации с использованием ускорения при помощи сопряженных градиентов, успешно применяемые для эллиптических уравнений. Существенный вклад в развитие этих методов, которые, в силу своих достоинств, открывают новые возможности решения пяти- и девятиточечных сеточных уравнений, внесли Н.И. Булеев, о. Ахе1ззоп, В.П. Ильин.

Параграф 1.3 посвящзн анализу применения метода неполной факторизации с ускорением по методу сопряженных градиентов (МНФСГ) при итерационном решении систем алгебраических уравнений, возникающих в случае неявной аппроксимации нестационарной двухмерной задачи теплопроводности. Рассматриваются схемы - чисто неявная, Кранка-Николсова и повышенной точности, аппроксимирующая исходную задачу с порядком . Проводится сравнение МНФСГ с други-

ми методами решения линейных и нелинейных систем.

Задача рассматривается в следующей формулировке:

и =ш{и(1 ,у,о:>,

3 у,1

о<е<1

2П в каждый момент

ю

ди

— = v(A.vu)-нp(^,J/). (х,у)<5№,

ог

Ы = В(Х,у), (Х,у>еЗМ, 0<1$Г«в,

и(Х,у,0)=и0(Х,у), (Х,у)еИ, 1=0,

?дв ЧИШ - двухмерная область определения функции и, I - линейный шератор, дЧ1 - область определения граничного условия, Л>0 - заданная функция координат.

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

дп_дп-1= _ "У (Аип-?п), Я0={и0(х±,у5)}.

'дй и11 - искомый вектор с компонентами, соответствующими значе-1иям решения системы в узлах сетки в момент времени и11-1, Рп -гевестные векторы, а и Ь - величины временного и некоторого характерного шага по пространственной переменной. Блочно-трехдкаго-

[альная матрица А=Б-Ь-и ОНГ , МВ^, и^Ш^, 1=1,----1)

гредао латается стильтьесовской, то есть симметричной М-матрицей.

Для решения на текущем временном шаге системы линейных •равнений

кг _ ?12 _ , •де В= —Е+А, 7-и , + —и , рассматривается метод сопряжен-1 а

ых градиентов:

г° = е-Б7°. г1^1 = г5с-а5{Вр5с,

р° =

0,1 ~ {Врк,рк)> ук+1 = 75с+акрк,

П

(гк+1 ик+1}

Рк+1 = Рк= - (ГК>? .

с предобуславливазощей матрицей вида

К =

Здесь 0={С1) - блочно-диагоналъная матрица, определяемая ре-курентными соотношеними

С

ь2

г

Л2

где е - вектор с единичными элементами, 5 - диагональная матрица, а (С1) ) _ »трехдиагональная часть" матрицы С-1, 6 - числовой итерационный параметр, .

Исходя из естественного предположения, что решение исходной дифференциальной задачи - достаточно гладко, так что ип=ип~1+0(а), в качестве начального приближения для системы на п-ом временном шаге используется У°=ип~1, либо результаты интерполяции: V0-

=гип~1-пп_2, т&г- ?°=зи1Х~1-зип~г+ип~3, п^З.

Число итераций, необходимое для уменьшения нормы вектора ошибки в б-1 раз, оценивается неравенством

11гг(5)! Г ?г , К2 .0.5 Г? Н(б) « -, В £ 2) - |1+ - } + -

Если итерации заканчиваются при величине невязки »г5ев= то для вектора 7,п={и(1.,у^Л„справедливо

1 1 (Е+—о = +—5- Лг|", Кг Ь2

где Ш - вектор погрешность аппроксимации. Откуда, в силу свойств матрицы А, следует

п Ь2

т.е. величина погрешности итераций накапливается линейно.

За счет выбора е. которые могут быть сделаны сколь угодно (алыми, при условии е=71г«фи можно подучить выражение вида

иг"» ^ «ФИ.

При начальная ошибка »У0и по крайней мере 0(1). Если

ФИ~0(т+712), то в этом случае вук«/йу°в~Лг.

Для задачи

ди. О2 и д2и

— = Л —д- +А. —5 , 0^1, Огда, 0<*г?1,

дг ахг ду

■ точным решением вида

)=(ау+Ь)-ехрШ+сх}.

ри \=с~г и полученными из него граничными условиями II рода были роведены расчеты с использованием метода блочной верхней релакса-ии и мнфсг.

На этой же задаче, но с граничными условиями I рода, рассмат-ивалось применение МНФСГ для схем, которые имеет погрешности ап-роксимации ), В последнем слу-

а у х у х у

ае получается девятиточечное сеточные выражение для которого тре-уется выполнение условия по соотношении шагов сетки

ЗХк-А^. у х

На основе проведенных экспериментов можно сделать следующее включение:

1. В случае использования метода неполной факторизации с старением по методу сопряженных градиентов итерации сходятся при вачениях итерационного параметра ве[0,1] с наилучшим при измене-ки количества узлов 6=1. Количество итераций при фиксированном т ронорционально Л-0-5, что соответствует теоретической оценке, риводимой в работе.

2. Метод неполной факторизации с сопряженными градиентами аачительно превосходит по экономичности метод блочной верхней злаксации и лишен недостатка последнего - подбора оптимального герационного параметра, связанного с й и т. Выигрыш увеличива-гся с уменьшением Л. При этом точность расчетов при использовании

обоих методов фактически одинаковая.

3. Связь между ростом временного шага и увеличением количества итераций в случае фиксированного к, определяется приблизительно как ч0"^, что близко теоретической оценке.

4. Использование в качестве начального приближения значений полученных по результатам интерполяции решений с предшествующих временных шагов, позволяет сократить общий объем вычислений.

5. Применение метода неполной факторизации с ускорением по методу сопряженных градиентов для схем повышенной точности приводит к значительному уменьшению величины погрешности вычислений, что дает возможность сокращения количества выполняемых операций за счет увеличения г и к.

6. При решении нелинейных задач тенлонереноса использование МНФСГ в сочетании с методом простой итерации ведет к снижению объема вычислений, особенно с уменьшением Л. Точность, приемлимая в решении практических задач, и сокращение времени счета позволяют считать такой алгоритмический подход очень эффективным. Ограничением является сходимость метода простой итерации при больших т.

7. Результаты численных расчетов подтверждают полученную теоретическую оценку устойчивости решения исходной задачи к погрешностям итерационного процесса на временных шагах.

8. Высокая эффективность метода неполной факторизации с ускорением по методу сопряженных градиентов позволяет использовать неявные разностные схемы как наиболее устойчивые при численном решений многомерных уравнений параболического типа.

В последующих параграфах рассматриваются результаты решения задач, подученные с использованием изложенного в первой главе алгоритмического подхода.

Во второй главе рассматривается две задачи связанные с металлургическим производством.

Использование непрерывного литья стали для получения заготовок дает возможность высокоэффективного производства при экономии энергии и увеличении выхода годного металла. В связи с повышением требований к скорости разливки, качеству и расширению марочного состава сталей, исследование тепловых условий формирования непрерывного слитка является актуальной практической задачей. Рассмотрение вопросов теплофизики, тепло- и массопереноса в тесной связи с технологией определяет все дальнейшие стороны интенсификации и

оптимизации технологических процессов. Наиболее существенный вклад в модельное исследование процессов получения непрерывного слитка внесли В.А. Емельянов, Е.М. Китаев, Ю.А. Самойлович, В.А.Куравлев.

В & 2.1 предлагается математическая модель, а также алгоритм решения использовавшийся при проведении вычислительных экспериментов по отработке рациональных режимов затвердевания круглого стального слитка в машине непрерывного литья заготовок и определении зависимости, связывающей величины коэффициентов теплоотдачи с удельным расходом воды в зоне вторичного охлаждения установки.

В большинстве известных работ рассматриваются установившиеся процессы затвердевания при заданных значениях коэффициентов теплоотдачи в секциях зоны вторичного охлаждения, без определения связи их величины с расходом воды. В параграфе рассматривается нестационарная математическая модель, описываемая квазилинейным параболическим уравнением, которая позволяет рассматривать все стадии формирования круглой стальной заготовки в машине непрерывного литья: от момента, когда затравка находится в кристаллизаторе и происходит "намерзание" металла, до выхода процесса кристаллизации на установившийся режим. Это дает возможность исследовать переходные процессы, когда по тем или иным технологическим причинам при литье меняются скорость вытягивания слитка или режим вторичного охлаждения.

При моделировании процесса кристаллизации стального слитка используется приближенно квазиравновесной двухфазной зоны для многокомпонентного сплава. Влияние конвекции расплава в кристаллизаторе предусматривается введением эффективного коэффициента теплопроводности. Учитывается зависимость теплсфизических параметров стали от температуры. Предполагается аксиальная симметрия процессов охлаждения и затвердевания.

Соотношение, связывание температуру в двухфазной зоне с объемной долей жидкой фазы, имеет следукщий вид:

Н Н

?=?А- I РгФГ^1'^- где I

1=1

- температура ликвидуса стали при исходном химическом составе Сг=С°. С учетом принятых допущений уравнение теплопереноса представляется в виде

'Эф|

ат ат

— —

01

дг

ч 1 д г дТ ч д г дТ .

= _ — [гЛ. — | + — [Л — |, > г дг дг > дг дг >

где сдф - эффективная теплоемкость стали, определяемая соотношениями

сэф~

^ п -(г+8-в к.) п

45 ¿=1

(ср),. /.,-0,

(ср)гз=(ср)г/г+(ср)5/5, 4=1-/г

Здесь индексами 1,5 обозначены параметра для жидкой и твердой фаз: ге^ - теплота кристаллизации при фазовом переходе; - объемная доля жидкой фазы; С^, - концентрация, коэффициент распреде-

ления и модуль углового коэффициента наклона линии ликвидуса для 1-того растворенното примесного компонента; е1=р3/рг, £=^-1.

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

На рис. I представлены некоторые результаты расчетов температуры' поверхности круглого стального слитка (Ст. 20) в зоне вторичного водовоздушного охлаждения (при работе трех секций), которые получили подтверждение опытными измерениями на промышленной установке Западно-Сибирского металлургического комбината.

Для зоны вторичного охлаждения, используя опытные данные в качестве реперных точек, из результатов расчетов были определены средние по длине секций значения коэффициентов теплоотдачи а, и устанавлена их зависимость от удельного расхода воды # в виде

Х=1,2,3,4, ^=<42=312, А3=А4= 211.

Графики на рис.2 иллюстрирую установленные экспериментально-расчетным путем зависимости коэффициентов теплоотдачи ак от удельного расхода воды по секциям, которые располагается в области данных полученных другими авторами.

т, °с

1000-

и 11IнI|IIм|11|||1111| 1005 10 15 20 25 0.

а, Вт/ы2К

д, ы /м ч

1—ГТТТПТ]-1 III 11П| п

1 1 10

о

Рис.1. Изменение температуры на Рис.2. Экспериментальные и рас-повергности (сплошная линия) и четные (сплошные линии) завкси-оси (штриховая линия) стального мости коэффициентов теплоотдачи слитка (Ст.20) в установившемся от удельного расхода воды, реааше, при скорости вытягива- -ния 0.57 м/мия (1 ), 0.67 и/иин (2) (Л, о - опытные данные).

Результаты численного моделирования позволили' определить, оптимальные условия формирования отливок, температурные поля в затвердевшей корочке металла, глубину жидкой души, локальные значения скоростей охлаждения и кристаллизации, а.также оценить характерные размеры дендритной структуры получаемых слитнов.

В металлургии йолътоо внимание уделяется проблемам затвердевания расплавов в формах, и основная масса работ посвящена прогнозированию свойств затвердеващэго металла, а проблемы изложниц рассматривается только с точки зрения оценки в них теплового потока и поглощения тепла. Главные требования к изложницам состоят в обеспечении заданного ресурса работы и требуемой скорости оборачиваемости. Процессы, отвечахщие за их эксплуатационные свойства, основываются на законах теплопроводности и тармомеханики поверхностных явлений. Уровень термоупругих и термопластических напряжений изложницы определяются температурными условиями ее

эксплуатации.

Целью рассматриваемой в § 2.2 задачи было исследование изменения температурного поля изложницы сложной геометрической формы с затвердевающими и остыващими в ней пятью слитками для последующей оценки величин возникающих напряжений- При решении были приняты следующие упрощающие"допущения:

- в виду того, что время заливки много меньше как периода затвердевания, так и последующего охлаждения, заполнение считается мгновенным;

- на поверхности сквозного отверствия в центральной части изложницы задаются адиабатические условия теплообмена, на верхней торцевой плоскости и боковой поверхности - условия теплообмена излучением и свободной конвекцией, а на нижнем торце - кондуктив-ного теплообмена с поддоном:

- в соответствии с технологическими требованиями, охлаждение изложницы со слитками осуществляется вначале излучением в вакууми-роваяном объеме печи, а после извлечения из нее - на воздухе.

В основу математической модели положено нестационарное уравнение теплолереноса в трехмерном пространстве с учетом выделения теплоты кристаллизации па границе фазового перехода при затвердевании слитков в приближении Стефана. Теплофизические параметры для затвердевавдего металла (урана) считаются зависящими от температуры и задается согласно справочным данным, а для материала изложницы - постоянными, равными их средним значениям в рассматриваемом интервале температур. В момент заполнения изложница находится в печи и начальное распределение температурного поля задается в ней согласно данным измерений.

На поверхность каждого из пяти рабочих каналов изложницы нанесено теплозащитное покрытие заданной толщины. Между поверхностями рабочего канала изложницы, покрытого слоем краски, и слитка выполняется условие непрерывности потоков тепла.

Ввиду равномерного в азимутальном направлении расположения заливочных каналов в изложнице, решение задачи осуществляется для сектора (см. рис. 3). На границах смежных с соседними секторами выполняется условия симметрии.

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

три системы цилиндрических координат: одна связана с центром слитка (0 ). другая - с центром изложницы (0), третья расположена в плоскости проходящей через линии симметрии между двумя секциями (0 ). В соответствии с этим рассматривается решение задачи в четырех сопряженных областях Л.В,С,0, для каждой из которых строилась своя пространственная сетка.

Рис. 3. Схема узлов сетки в поперечном сечении фланца (а) и изложницы (б) { - - лингси сетки;---- - рабочий канал).

После заливки металла в рабочие каналы происходит разогрев изложницы за счет перегрева и теплоты кристаллизации. Слиток затвердевает в течение 4-5 минут. Наибольшая температура системы наблюдается в нижней плоскости фланца, совпадающей с поверхностью зеркала металла. Разогрев в первоначальный период имеет неравномерный характер из-за сложной геометрической формы моделируемого объекта. Последующее выравнивание по сечению наступает спустя 20-30 минут, а окончательное перераспрэленда заввршается к 60-ой минуте рассматриваемого цикла, после чего охлаждение система происходит вследствии теплоотдачи от боковой поверхности и дна изложницы.

В течении всего технологического цикла наиболее разогретыми являются части изложницы в области А, прилегащие к рабочим каналам, на внешней границе которой замкнутая изотерма имеет форму, близкую к аксиальной симметрии. В этой цилиндрической области возникает плоская осесимметричная деформация. Для приближенной количественной оценки соответствующих термических упругих напряжений используется полученное решение уравнения теплопроводности. Вычисленные значения радиальной (ог) и азимутальной (с^) составляющих напряжений, радиального перемещения (ц^) приводятся в таблице.

Таблица.

И, м аг, Ша а , МПа мкм 7, °С

0.103 0.00 -694.92 -326.36 798.4

0.108 -31.38 -661.45 -275.46 797.0

0.113 -58.61 -632.61 -227.63 795.9

0.120 -91.15 -598.29 -164.87 794.6

0.127 -118.32 -569.70 -106.52 793.8

0.134 -141 .31 -545.55 -51.72 792.8

0.141 -160.90 -525.02 0.00 792.2

йз результатов следует, что хотя рассчитанные значения су и о не превосходят предела прочности на сжатие серых чугунов (800 Ша), возникающие в рассматриваемой области значительные азимутальные сжимащие напряжения способствуют появлению поверхностных трещин, которые в процессе эксплуатации (вследствии термоцик-лирования) сливаясь, могут привести к разрушению рабочих каналов.

Третья глава посвящена исследованию переходных процессов теп-лоперецоса при формировании интрузивных магматических массивов.

В настоящее время при моделировании этой задачи в основном ограничивайся заданием простейших условий - "мгновенное" заполнение камеры, однородное распределение в ней температур и фаз, рассмотрение процесса затвердевания в приближении Стефана, стационарный расход жидкости. Однако состояния в реальных магматических системах во время и после прекращения напорного течения расплава далеки от такой идеализированной картины.

Количественное описание заполнения интрузивной камеры может быть самосогласованно выполнено только в рамках общей модели дина-

мики развития магматической системы с момента генерации расплава в очаге и заканчиваться его затвердеванием. Последовательно замкнутое решение этой задачи излагается п работах Ю.П. Желтова. М.Р. Ryan, D.D. Pollard, А.Н. Черепанова и В.Н. Шарапова. При этом проблема распадается па две взаимосвязанные, но представляющие самостоятельный интерес, задачи: иптруднрование магмы в трещинные проводники и заполнение магматических камер различной формы.

Для анализа динамики охлаждения расплава при интрузии используется простейшая структура магматической системы и, с учетом упрощений в постановке задачи, рассматривается следующая физическая модель.

1) Давление в магматическом очаге, представляющим собой некий резервуар с базитовым расплавом, избыточное против литостатичес-кого из-за наличия стресса и разности плотностей жидкости и вышележащих ультрабазитовых пород.

2) Под действием этого давления над очагом зарождается однородный плоский канал, в который начинает поступать магматическая жидкость. Головная часть потока представляет собой горизонтальную плоскость контакта между расплавом и вмещающими породами, на котором жидкость обменивается теплом с вмещающей средой по сложному закону, а интенсивность теплообмена определяется коэффициентом теплоотдачи.

3) Избыточное давление расходится на преодоление гидравлического сопротивления, которое зависит от ширины канала и может меняться со временем из-за намораживания изверженных пород, развития гетерофазной зоны, кипения магмы, а также на поддержание некоторой ширины трещины, не связанной с процессом фазовых переходов в потоке расплава. Гидравлические потери в головной части потока определяет скорость движения магмы по каналу.

4) Интрузивная камера возникает над проводником в земной коре, подвергнутой рифтогенезу. и процесс образования резервуара продолжается до тех пор, пока давление в магме не сравняется с литоста-тическим или у ее верхней границы не наморозится достаточно прочная корка. При формировании камеры предполагается только "всплыва-ние" ее кровли и не предусматривается разрушение толщи пород.

Нужно отметить, что такая схема является модификацией модели Райяна, в которой "плавучесть" расплава является элементом не только процесса интрудирования магмы по проводнику, но и заполне-

ния камеры.

Обще потери давления в потоке магмы по проводнику учитываются введением некоторого эффективного коэффициента трения кото-

рый считается постоянным, поскольку не известен закон его изменения. Это упрощающее предположение позволяет учесть при решении тепловой задачи изменение средней скорости течения магма в проводнике за счет "срабатывания" давления и уменьшения вплоть до нуля "плавучести" расплава.

В потоке, двигающемуся по проводнику, из-за охлаждения развивается гетерофазная зона и намораживается корка изверженных пород. При заполнении камеры происходит сходный процесс охлаждения расплава, особенности которого зависят от характера его растекания в возникающей полоста.

В рамках сформулированной во введении третьей главы модели, в § 3.1 осуществляется математическая постановка задачи. Процесс кристаллизации магмы в двухфазной зоне описывается в рамках квазиравновесного приближения для многокомпонентного сплава.

Для определения изменения температурного поля во вмещашей породе и магме применялся алгоритм основанный на балансной разностной схеме. Узлы сетки вдоль длины проводника располагаются таким образом, чтобы все пространство заполняемой полости делилось на равные объемы. Первоначально просчитывается гидродинамическая задача и находится время, за которое верхняя граница потока магмы достигает уровня следующего узла сетки. При этом в канал поступает фиксированный объем расплава, а его незатвердевшая часть смещается.

В рамках сформулированной модели оказалось возможным получить оценки высоты подъема магматической колонны по проводнику к определить изменение средней скорости движения потока по каналу в зависимости от значений эффективных коэффициентов трения при заданной величине избыточного давления в очаге над гидростатическим. Правдоподобность выбора параметров была оценена сравниванием расчетных величин скоростей течения магмы на поверхность земли с наблюдаемыми значениями для магматической системы вулкана Килауэа.

Сопоставление рассчитанных и наблюдаемых величин как скоростей движения магматических "фронтов", так и содержания объемной доли твердой фазы в базитовых расплавах при извержениях позволяет предположить, что в реальных магматических системах диапазон варь-

рования Л.Эф лежит где-то в интервале от нескольких десятков до ¡00-500.

Неизотермичность процесса ивтрудирования базитового расплава :з очага генерации накладывает ряд существенных ограничений на со-1Тношения в качественной модели возникновения и развития нагмати-:еских систем. Одно из них состоит в том. что скорости движения [агмы не могут быть ниже некоторого критического уровня, поскольку : противном случае происходит "перемораживание" проводника.

Впервые с использованием квазиравновесного приближения для »писания развития гетерофазной зоны в § 3.2 получены параметры, ^растеризующие нестационарный процесс охлаждения лавы в магмати-:еской камере в процессе развития системы. Ввиду того, что реаль-не магматические тела имеют довольно разнообразные геометрические срмы, для оценки возможных режимов затвердевания рассматривалось ;аполнение камер трех морфологических типов, предельных в отноше-ми особеностей охлаждения расплава: воронковидных (V), лякколито-'ИДНЫХ (А), силлообразных (а).

По результатам расчетов оказалось, что можно ожидать минимум за динамических режима охлаждения и солидификации расплава как 1ри заполнении камеры, так и после прекращения напорного движения ¡агмы: 1} при скоростях течения в проводнике выше критических, :) при докритических скоростях течения в канале. Объемы тел, кото-ие могут сформироваться при докритических скоростях течения 1асплава в питаквдх проводниках шириной не более 5 м, составляют 1.03-0.2 км3. Это весьма существенное петрогенетическое следствие, ©тения нестационарной задачи интрудирования базитовых магм.

Получены также особенности охлаждения и солидификации распла-:а для закритических режимов заполнения камер. Характер солидафи-ации в воронковидной камере связан в первую очередь с тем, что [ри заполнении полости приток магмы постоянно падает, а площадь юверхности теплосъема увеличивается. Отсюда следует существенный х5ъбм затвердевшего расплава и большое разнообразие в режимах »хлаадения. В случае лакколитовидного тела ситуация прямо противо-Еоложная, так как при заполнении такой камеры отношение товерхнос-и теплосъема к объему расплава уменьшается. В силлоподобных или ¡токовидных интрузивных телах особенности теплосъема мало меняются :о временем. Поэтому для этого случая получается качественно сходке картины характера солидификации расплава при различных скорое-

тях заполнения.

Можно предположить, что для реальных магматических систем подученные оценки скорее всего следует относить к категории "верхних", поскольку в расчетах использована система, для которой принята линейная аппроксимация поверхности ликвидуса, а потому с температурами фазовых переходов на 100-150 градусов выше тех, что характерны для базальтовых расплавов.

Сравнение результов вычислений и данных полученных из наблюдений по величине содержания твердой фазы, динамике изменения температур в лавах при квазистационарных извержениях свидетельствуй об их совпадении. Таким образом определенные из численных экспериментов параметры интрудирования и солидафикации расплавов правдоподобны не только качественно. Мх значения близки к реальным. Это позволяет сформулировать некоторые петрогенетические следствия, определяющиеся динамикой интрудирования базитовых мата.

В настоящее время петрологи предполагает возможность проявления двух основных случаев солидафикации и дифференциации базитовых. расплавов в магматических телах после заполнения интрузивной камеры: 1) гравитационная, с сопутствующим комплексом механизмов фазовых переходов и тепло-массопереноса при содержании твердой фракции в объеме менее 30%; 2) затвердевание магмы в жидкотекучем ядре массива, содержащем 30-60 % твердой фазы, когда развиваются токи усадочной конвекции, фильтр-прессинг и ретроградное кипение в остаточной жидкости. Очевидно, что первый из вариантов возможен только при скоростях движения магмы по проводнику выше критических. Но и тогда он зависит от характера охлаждения магмы при заполнении, которая определяется формой камеры и ее размерами.

Для случая воронковидной геометрии можно предполагать достаточно жесткие границы проявления гравитационной солидафикации -приток расплава должен быть значительным и составлять не менее нескольких десятков кубических метров в секунду. То есть заполнение должно быть весьма быстрым. По оценкам это десятки или сотни часов. При этом даже если расплав был несколько перегрет в магматическом теле на стенках намораживается довольно значительный слой изверженных пород порядка нескольких десятков метров, около которого развивается гетерсфазная область примерно такой же мощности.

При заполнении камер другой формы такого рода осложнений, то есть значительных масштабов затвердевания расплава при скоростях

течения близких к критическим, видимо, не возникает. Расчеты показали, что в лакколитовидных или силложздобньтх магматических телах скорее реализуется первый из отмоченных выше режимов солидификации. Заложение таких камер при перегреве расплава характеризуется относительно небольшими размерами затвердевшей и гетерофазной зон.

Из результатов реиения нестационарной задачи динамики интру-дирования следует, что охлаждение базитового расплава в процессе заполнения камер существенно различается в зависимости от их формы. Следовательно и дальнейшие явления солидификации интрузивных тел, а также разделение компонентов в процессе фазовых переходов могут протекать различным образом как в локальных частях магматических тел, так и во всем объеме массива.

III. ОСНОВНЫЕ РЕЗУЛЬТАТЫ РАБОТЫ

1. На основе свойств монотонных матриц исследованы неявные разностные алгоритмы численного решения одномерных и двухмерных задач затвердевания сплавов с фазовым переходом в интервале температур, описываемых квазилинейными параболическими уравнениями с краевыми условиями 1-Ш родов, и получены оценки погрешности в равномерной норме.

2. Построены "'компактные" разностные схемы повышенного порядка точности для многомерных уравнений теплопроводности и определены соотношения между пространственными и временными шагами, при которых они является монотонными. Показано, что схема Кранка-Николсона с погрешностью порядка 0(а2+И} для решения задачи, описываемой трехмерным уравнением, абсолзотно немонотонна.

3. Исследована возможность применения итерационного метода неполной факторизации с ускорением по методу сопряженных градиентов для решения линейных и нелинейных систем, получаемых при неявной аппроксимации многомерных задач теплопереноса, и показана высокая эффективность такого алгоритмического подхода.

4. Предложена нестационарная модель, описывающая получение круглого стального слитка в машине непрерывного литья, и численно исследованы температурные режимы переходных процессов формирования заготовки: от момента запуска установки до выхода процесса затвердевания металла на установившийся режим. По результатам расчетов и опытных данных определены зависимости, связывание значения

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

5. 11а основе численного моделирования проведено исследование нестационарных теплофизических процессов в многоканальной чугунной изложнице при затвердевании в ней пяти цилиндрических слитков. Согласно полученному решению задачи рассмотрен характер изменения температурных полей системы и проведена оценка возникающих термонапряжений в изложнице.

6. С помощью численного решения нестационарной задачи динамики интрудирования базитовых расплавов по щелевидаому проводнику, получены оценки средних скоростей движения магмы и высоты ее подъема, определены зависимости критических скоростей от размеров проводника. Эти результаты, с учетом размеров дайковых проводников и базитовых плутонов» дают возможность очертить диапазон расходов расплава при формировании магматических камер и вероятные времена заполнения.

7. Решение нестационарной модели, описывающей динамику формирование интрузивов, показывает существование различных режимов затвердевания магма, определяемых, ее движением и геометрической формой заполняемых камер. Из этого можно сделать вывод, что и последующая солидафшсация» а также разделение компонентов в процессе фазовых переходов могут протекать различным образом как в локальных частях, так и во всем объеме интрузивного массива.

IV. ПУБЛИКАЦИЙ ПО ТЕМЕ ДИССЕРТАЦИИ

1. Попов В.Н., Черепанов А.Н., Манолов В.Н. Сравнение результатов численного исследования процесса затвердевания стального слитка с помощью равновесной и квазиравновесной математических моделей // Численные методы и приложения: Сб. научн. тр. - София, 1339. - с. 385-389.

2. Ильин В.П.. Попов В.Н. Сходимость разностного решения для одномерной задачи кристаллизации металла в форме // Численные метода механики сплошной среды: Сб. научн. тр. - Новосибирск, 1989. - Т. 3(20), № 4. - С. 73-82.

3. Попов В.Н. Сходимость разностной схемы для двухмерного квазилинейного уравнения параболического типа // Вариационно-разностные методы в задачах численного анализа: Сб. научн. тр. -

Новосибирск. 1991. - с. 124-133.

4. Ильин В.П., Попов В.Н. О применении метода неполной факторизации для решения двухмерного параболического уравнения // Вариационно-разностные методы в задачах численного анализа: Сб. научн. тр. - Новосибирск, 1991. - С. 61-75.

5. Попов В.Н. Применение метода неполной факторизации для решения задач теплопереноса с использованием неявных разностных схем.

- Новосибирск. 1997. - 15 с. - (Препр./ РАН. Сиб. отд-ние. йн-т теорет. и прикл. механ.; > 4-97).

6. Черепанов А.Н., Попов В.Н., Айзатулов P.C. и др. Расчетно-экспериментальное исследование тепловых режимов формирования непрерывного слитка стали // Изв. ВУЗов. Черная металлургия.

- 1997. -» 8. - С. 43-47.

7. Черепанов А.Н.. Манолов В.К., Попов В.Н., Рашев П.В., Млиев О.П. Математическое моделирование затвердевающей стали с учетом образования двухфазной зоны // Высокоазотистые стаж-SS: Труды Междунар- научно-техн. конф. - НРБ, Варна, 1989.

8. Черепанов А.Н.. Попов В.Н.. Караник Ю.А., Григорьева Г.М. Оптимизация процесса формирования отливки с помощью вычислительного эксперимента // Кристаллизация и компьютерные модели: Труды IV Всесошной конф.- 1991 г.- Ижевск, 1991.- С. I07-II5.

9. Черепанов А.Н., Попов В.Н., Георгиев Г.. Караник Ю.А. Математическое моделирование теплофизических и термоупругих процессов при 'затвердевании стальной отливки // Литейное производство. -1994. -ЯЗ.

10. Черепанов А.Н.» Попов В.Н. Компьютерная модель формирования профильного изделия, получаемого вытягиванием вниз из расплава // Кристаллизация и компьютерные модели: Труда V Междунар. конф., охт. 1992 г. - Ижевск. 1994. - С. 69-74.

11. Черепанов А.Н., Казаков A.A., Плаксия С.М., Попов В.Н. Компьютерное исследование формирования химических соединений переменного состава при кристаллизации /У Кристаллизация и компьютерные модели: Труда V Междунар. ковф., окт. 1992 г. - Ижевск, 1994. - С. 75-84.

12. Черепанов А.Н., Попов В.Н., Караник Ю.А., Георгиев Г. Численные исследования тепловых термоупругих процессов формирования цилиндрической отливки // Применение компьютеров в литейном производстве: Материалы III Междунар. конф., июль 1990 г. -

НРБ. София, 1990.

13. Черепанов А.Н.» Попов В.Н., Черепанова В.К. Моделирование динамики тепловых процессов в изложнице и затвердевающих слитках металла // Металлы. - 1998. - Jí 2.

14. Попов В.Н., Черепанов А.Н., Аброськин А.Е. и др. Компьютерное моделирование процессов затвердевания металла и тепловых полей в многоочзсовой чугунной изложнице // Кристаллизация: компьютерные модели, эксперимент, технология : Тезисы VI Междунар. конф.. 1ЭЭ4 г. - Ижевск, 1994. - С. 50-51.

15. Sharapov V.U., Cherepanov A.M., Popov V.N., Xobov A.G. Dynamic solidification oí basic melt at enplacement oí magmatic chambers of different shapes // Experiment in geosciencea. 1996. - Л 2.

16. Шарапов B.H., Черепанов А.Н., Попов В.Н.. Лобов А.Г. Динамика охлаждения базитового расплава при заполнении воронковидаой интрузивной камеры // Петрология. - 1997. - Т. 5, » I.

17. Шарапов В.Н., Черепанов А.Н., Попов В.Н. Термодинамические условия развития и вырождения структурных фаций магматических тел // Геология и геофизика.- 1995.- Т. 36, » 12. - С. 80-98.

18. Шарапов В.Н., Черепанов А.Н., Попов В.Н. Динамика развития и формирования лополитов // Докл. РАН. - 1997. - Т. 352. & 5.

19. Черепанов А.Н., Шарапов В.Н., Попов В.Н. Динамика охлаждения базитового расплава при заполнении магматических камер // Геология и геофизика. - 1996. - Т. 37, Л 7. - С. 47-59.

20. Шарапов В.Н., Черепанов А.Н., Попов В.Н. Динамика охлаждения базитового расплава при заполнении воронковидаой интрузивной камеры // XIII-е Российское совещание по экспериментальной минералогии: Тез. докл., сент. 1995 г. - Черноголовка, 1995.

7 л ■

О) а&ез-омл/ог

СИБИРСКОЕ ОТДЕЛЕНИЕ РОССИЙСКОЙ АКАДЕМИЙ НАУК ИНСТИТУТ ТЕОРЕТИЧЕСКОЙ И ПРИКЛАДНОЙ МЕХАНИКИ

На правах рукописи

ЙФ

аре

мир

[лаевич

•"5 'А

'НЫХ^ШШОФЙБЙЧЕСКЙХ ЗАДАЧ ТЕМПЕРАТУР

ОБ^З^ёг^ применение вычислительной техники, математического моделирования и математических методов в научных исследованиях

Диссертация на соискание ученой степени доктора физико-математических наук

Научные консультанты: доктор физико-математических наук, профессор Ильин В.П. доктор физико-математических наук, с.н.с. Черепанов А.Н.

Новосибирск, 1998

СОДЕРЖАНИЕ

ВВЕДЕНИЕ 5

ГЛАВА I. НЕЯВНЫЕ РАЗНОСТНЫЕ СХЕШ ДЛЯ РЕШЕНИЯ ЗАДАЧ

ТЕПЛОПРОВОДНОСТИ 23

§ 1.1. Алгоритм численного решения одномерной

задачи затвердевания металла в изложнице 30

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

1.1.2. Описание алгоритма решения 31

1.1.3. Исследование сходимости итерационного процесса 34

1.1.4. Оценка погрешности разностного решения 34

1.1.5. Выводы 39

§ 1.2. Неявные разностные схемы для многомерных

уравнений теплопроводности 40

1.2.1. Сходимость разностной схемы для двухмерного квазилинейного уравнения параболического типа

с краевыми условиями 1-Ш родов 40

1.2.2. Исследование условий монотонности разностных схем повышенной точности для решения нестационарных двухмерных и трехмерных

уравнений теплопроводности 49

1.2.3. Схема повышенного порядка точности для уравнения тешюпереноса в цилиндрической

системе координат 58

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

1.2.5. Выводы 68

§ 1.3. Использование метода неполной факторизации с ускорением по методу сопряженных градиентов для решения нестационарных двухмерных задач теплопереноса 74

1.3.1. Описание и свойства итерационного алгоритма 74

1.3.2. Сравнивание МНФСГ и метода блочной верхней релаксации по результатам численных расчетов

модельной задачи 78

1.3.3. Результаты -численных расчетов при использовании МНФСГ для схем повышенной точности 85

1.3.4. Результаты расчетов нелинейной модельной задачи 96

1.3.5. Выводы 106

ГЛАВА 2. ЧИСЛЕННЫЕ РЕШЕНИЯ ЗАДАЧ ТЕПЯОПЕРЕНОСА В МЕТАЛЛУРГИЙ 114

§ 2.1. Исследование тепловых режимов формирования

непрерывного слитка стаж 119

2.1.1. Постановка задачи и основные уравнения 119

2.1.2. Алгоритм решения задачи 123

2.1.3. Определение зависимости коэффициентов теплоотдачи от расходом воды при теплообмене в зоне вторичного охлаждения 124

2.1.4. Оценка структуры получаемого слитка 128

2.1.5. Выводы 130

§ 2.2. Исследование тепловых процессов в пятиканальной

изложнице при затвердевании металла 138

2.2.1. Математическая постановка задачи 138

2.2.2. Алгоритм численного решения 141

2.2.3. Результаты расчетов 145

2.2.4. Анализ термонапряженного состояния изложницы 148

2.2.5. Выводы 151

ГЛАВА 3. ИССЛЕДОВАНИЕ ДИНАМИКИ ОХЛАЖДЕНИЯ ВАЗИТ0В0Г0 РАСПЛАВА.

ПРИ ЗАПОЛНЕНИИ ПЛОСКИХ МАГМАТИЧЕСКИХ КАНАЛОВ И КАМЕР 162

§ 3.1. Численные исследования процесса интрудирования

базитовой магмы по щелевому каналу 170

3.1.1. Математическое описание модели 170

3.1.2. Алгоритм решения и значения принятых

в расчетах физических параметров 176

3.1.3. Динамика интрудирования базитовой магмы

по проводнику 179

3.1.4. Выводы 181

§ 3.2. Динамика охлаждения базитового расплава при

заполнении магматических камер 190

3.2.1. Модель заложения магматической камеры 190

3.2.2. Алгоритм решения задачи 191

3.2.3. Условия формирования магматических тел 1%

3.2.4. Результаты численных исследований динамики охлаждения расплава в магматических камерах 195

3.2.5. Выводы 199

ЗАКЛЮЧЕНИЕ 210

ЛИТЕРАТУРА 212

СПРАВКА 0 ВНЕДРЕНИИ РЕЗУЛЬТАТОВ 232

ВВЕДЕНИЕ

Настоящая работа посвящена разработке и применению математических моделей описывающих теплофизические процессы в металлургии и петрологии, связанные с затвердеванием расплавов в интервале температур, а также развитию математического аппарата решения возникающих задач.

Возрастающие требования к современной технологии получения металла выдвигают задачи исследования процессов затвердевания сплавов с помощью построения адекватных математических моделей [1-7,46,47,501. Следует отметить, что методы моделирования с использованием быстродействующих ЭВМ часто оказываются единственным способом определения оптимальных режимов, при которых происходит тот или иной процесс, так как возможности экспериментального изучения этих явлений ограничены методическими и техническими трудностями.

При фундаментальном теоретическом подходе к исследованию затвердевания металла требуется анализ всего комплекса происходящих процессов, т.е. следует учитывать тепловые, диффузионные и гидродинамические явления, протекающие одновременно. Однако при этом математическая модель получается сложной и трудно применимой в численных исследованиях [83. Теплофизические процессы в связанной динамической системе развиваются с неодинаковыми скоростями, т.е. характерные времена их протекания неодинаковы. Это дает возможность оценивать основные особенности изучаемых явлений и тем самым при разработке математической модели и проведении последующих расчетов учитывать только главные "движущие силы". В случае необходимости этот подход может реализовыватся не только для всего рассматриваемого временного интервала, но и для отдельных его этапов. Поэтому на практике нередко используются упрощенные вари-

анты моделей для анализа и исследования тепловых режимов, что вполне оправдано.

При затвердевании расплавов сложного состава, т.е. переходе из жидкого состояния в твердое, ведущая роль принадлежит процессам тепломассообмена. На основе теории теплопроводности решается множество задач, основные из которых - определение динамики нарастания объема твердой фазы, продвижения фазового фронта и продолжительности процессов затвердевания.

Традиционная постановка задачи о затвердевании жидкости, сформулированная австрийским математиком Стефаном более ста лет назад 191, основана на представлении о существовании поверхности фронта кристаллизации, где выполняются следующие условия:

а) кристаллизация (плавление) происходит на достаточно гладкой границе раздела фаз, определяемой уравнением )=0, которое можно переписать в виде

г=Фи,у>г), <р(*»у.о)=Ф0и,уь <х»у,Оеа=аг-на5>

где аз - области жидкого и твердого состояний вещества;

б) фазовый переход происходит при постоянной температуре кристаллизации

т = т

1 z кр*

в) скорость движения границы раздела фаз определяется условием отвода скрытой теплоты фазового перехода

дТ дТ

э^р V = к — - кг —, -У 5 П 5 дп I. дп

где А. - коэффициенты теплопроводности вещества в жидком и твердом состояниях, р3 - плотность твердой фазы, э^ - удельная

скрытая теплота фазового перехода.

После дополнения математической формулировки задачи уравнениями нестационарной теплопроводности для жидкой и твердой фаз, соответствующими начальными и граничными условиями, можно определить закон затвердевания расплава (движение границы раздела фаз), а также распределение температурных полей в системе с течением времени.

Однако в настоящее время приближение Стефана используется при решении узкого круга задач: затвердевание чистых металлов, рост монокристаллов, кристаллизация с гладкой границей фазового перехода в контролируемых условиях.

Изучение динамики затвердевания сплавов при использовании различных экспериментальных методов позволило сделать вывод о том, что фронт кристаллизации существует в виде пространственно-распределенной двухфазной зоны и лишь в очень редких случаях как совершенно гладкая поверхность раздела фаз. Именно в переходной зоне закладываются многие структурные, химические и механические свойства металла [1,10-173. Влияние каким либо доступным способом на процессы происходящие в переходной зоне, дает возможность получения слитков или готовых изделий с наперед заданными характеристиками [48,491. Поэтому особый интерес представляет определение зависимости изменения температурных полей от типа сплава, начальных условий заливки и внешних условий охлаждения во время процесса затвердевания.

Формирование слитка обусловливается совместным протеканием процессов затвердевания и кристаллизации, которые являются двумя неразрывными сторонами сложного процесса перехода расплава из жидкого состяния в твердое. Кристаллизацией называется процесс формирования структуры металла с учетом всех сопровождающих его явле-

ний. Под затвердеванием понимают процессы теплопереноса в условиях фазового превращения. В любой момент затвердевание и кристаллиза-зация характеризуются одщой и той же скорость протекания, но роль их в процессе перехода из жидкого состояния в твердое неодинакова. Процесс кристаллизации связан условиями возникновения и роста кристаллических зародышей, а процесс затвердевания, связанный прежде всего с теплоотводом, определяет окончательное формирование. твердой фазы слитка £181.

Таким образом, теплоперенос в условиях фазового превращения определяет обе стороны процесса перехода металла из жидкого состояния в твердое. Тепловые условия в слитке (в частности температурное поле) являются главными факторами.

Процесс затвердевания, связанный с теплоотводом, определяется теплофизическими свойствами расплава: теплопроводностью, теплоемкостью и скрытой теплотой кристаллизации. Скрытая теплота кристаллизации физически соответствует энергии, которую необходимо отвести от металла при переходе от неупорядоченной структуры жидкости к упорядоченной структуре твердой фазы. При этом энергия тепловых колебаний атомов не изменяется, что и соответствует постоянству температуры при этом переходе. Затвердевание сплавов происходит в интервале температур ликвидуса и солидуса Эти температуры определяются диаграммой состояния сплавов. Скрытая теплота кристаллизации выделяется в пределах двухфазной зоны в интервале

В последнее время для решения задач затвердевания металлических сплавов широко и наиболее успешно используется теория "квазиравновесной" двухфазной зоны, предложенная и развитая В.Т. Борисовым £19-223. Эта теория не учитывает концентрационное (или диффузионное) переохлаждение расплава. Так как его величины в реальных

условиях малы (по результатам экспериментов не более 1°), при математическом моделировании этим параметром можно пренебречь. В указанной теории выявлена высокая скорость роста твердой фазы при малых отклонениях от равновесных условий. Это позволяет предполагать, что в каждом элементарном объеме двухфазной области жидкость и твердая фаза находятся в равновесии, а концентрация жидкости связана с температурой условием равновесия, т. е. уравнением линии ликвидуса расплава: где С - концентация растворенного

элемента.

В систему уравнений теории квазиравновесной двухфазной зоны входят уравнения теплопроводности и переноса массы в жидкой части двухфазной области. Решение этой системы в соответствии с началь-начальными и граничными условиями позволяет с достаточно высокой точностью определить температурное поле в слитке, поле концентраций и объемной доли жидкой фазы.

Ввиду малой диффузионной подвижности компонентов в твердой фазе, распределение вещества в растущих кристаллах неоднородно и полное равновесие в системе не достигается. Считается, что оно имеет место только на фазовой поверхности. Поэтому кристаллизация завершается не при температуре равновесного солидуса, соответствующей исходной концентрации компонентов в расплаве С)Э> а при более низкой - в связи с обогащением остаточной жидкости примесными компонентами. Считается, что затвердевание полностью завершается при температуре, соответствующей особой точке диаграммы состояния: эвтектике ?Е> перитектике ? и т. д. Чаще всего используется При этом сечение жидкой фазы меняется скачком от значения до нуля. В практических расчетах при численном моделировании окончание затвердевания определяется заданием величины 0.025-0.05 в конце двухфазной зоны. При исследовании чисто тепло-

вого режима кристаллизации диффузией компонентов в жидкости пренебрегают, то есть не учитывают макроскопическую дифференциацию. Из предположения того, что Си? связаны, следует: сечение жидкой и твердой / фаз, а также спектральная теплота кристаллизации эе^р З/^/дУ являются функциями температуры и вполне определяются видом диаграммы состояния.

В пределах двухфазной зоны уравнение теплопереноса имеет следующий вид [223

дТ а/,

{Ре)18 ^ =^СЛ.258гай(?)3-ае)Эр5 , (х.у.й,)^.

Последнее слагаемое отлично от нуля только в интервале температур кристаллизации. Индекс 1з характеризует соответствующие теплофизи-ческие параметры.

Функция относительного количества жидкой фазы, кристаллизующейся в двухфазной зоне отливки, определяется по диаграмме состояния, согласно различным методикам. Например, для линейной диаграммы состояния по правилу равновесного рычага [53:

Г тго~т 1п

где показатель степени п варьируется для разных сплавов и зависит также от концентрации примесей в материалах; либо по правилу неравновесного рычага для бинарных сплавов [21,233:

/г

1/(1-6,)

либо, для многокомпонентны! сплавов, выражением £24}

N

связывающим температуру в двухфазной зоне с сечением жидкой фазы Здесь Тк - температура плавления основного компонента сплава, ?20 - равновесная температура ликвидуса, - исходная

концентрация, коэффициент распределения и модуль углового коэффициента наклона линии ликвидуса для I-того растворенного примесного компонента; 81=р5/рг. Существуют и другие способы определения , см. £4,6,18,22,25] и цитируемые в них работы. Выбор того или иного приближения сказывается на точности определения ширины двухфазной зоны, времени затвердевания, оценке размеров структурных параметров сплавов £26,511.

В связи с тем, что определяется как функция температуры, то после использования подстановки

щ д^ дт_ ог дт м *

уравнение теплопроводности в двухфазной зоне переписывается в виде

дТ

где едф называется эффективной удельной теплоемкостью сплава в двухфазном его состоянии и определяется выражением

а для (рс) обычно используются выражения вида

(рс)г5=(рс)г/г+(рс)5/5.

Таким образом, учет выделения тепла кристаллизации в двухфазной зоне сводится к соответствующему заданию зависимости эффективной удельной теплоемкости сплава от температуры, а уравнение теплопереноса приобретает вид, характерный для квазилинейных параболических уравнений.

На границе жидкой и двухфазных зон т^ выполняется условия

-о

дТ

Т)г+0>

а для фронта конца кристаллизации в двухфазной зоне на границе т)3 справедливо

=Т

дТ

М?)

дТ

5 дп

%+0>

Используя уравнения теплопроводности для области перегретого расплава 0-,

Ч

дТ

(ре), — =<Илг[А.7(Г)згас1(Г)3, (г,у, 3,)^!,,

ь и с

И

(рс),

дТ

дг

для затвердевшей части 05, и задав начальные и граничные условия Ти,у,2,0)=Т0(х,у,г), (х,у,г)е0=аг+0г5+а5.

(Х,у,2)еШ,

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

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

0-

дТ

дТ

=М*>

V0" 5

V3,

иж алгоритмы сквозного счета для решения задачи Стефана [27-393. й хотя существуют подходы позволяющие при этом точно фиксировать фазовую границу 1341, наиболее распространенным является жбо сглаживание разрывных коэффициентов задачи £29,351, жбо введение функции теплосодержания. Последнее позволяет свести искомую задачу к нежнейному уравнению для разрывной функции, описывающему тепло-перенос во всей области включая фазовую границу. Так как функция имеет разрыв первого рода, то при использовании аппроксимации, она вместе с теплофизическими константами подвергается сглаживанию на некотором конечном интервале температур [27,29,323.

Границы фазового превращения в этом случае явно не выделяются. йх положение, когда это необходимо, устанавливается после расчета температурного поля сжтка £1,333.

Метода математического моделирования и алгоритмы численных исследований в металлургии нашж свое применение и в других областях науки. Теория магматических и магматогенных гидротермальных рудных месторождений связана с разработкой динамики дифференциации (перераспределения) магмы в камерах земной коры, которая интенсивно развивается три последних десятилетия [40-453. Магматогенные

процессы являются главными генераторами большинства рудных месторождений. Поэтому с момента возникновения петрологии развивается и совершенствуется основанное на модельных исследованиях прогнозирование при поиске руд, пространственно сопряженных с проявлением магм в форме интрузивных массивов и вужанов. Центральное место здесь занимает проблема описания процессов формирования рудонос