Задача Стефана о фазовом переходе

Размещено на http://www.allbest.ru/

Курсовая работа

по дисциплине "Уравнения в частных производных"

на тему: "Задача Стефана о фазовом переходе"

Содержание

Введение

1. Автомодельное решение классической задачи Стефана

2. Численные методы, применяемые для решения задачи Стефана

2.1 Метод ловли фазового фронта в узел сетки

2.2 Метод выпрямления фронтов

2.3 Метод сглаживания коэффициентов

2.4 О выборе параметра сглаживания

2.5 Разностные схемы сквозного счета

Заключение

Список используемой литературы

Введение

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

математический модель задача стефан

1. Автомодельное решение классической задачи Стефана

Классической задачей Стефана называют простейшую одномерную задачу промерзания (оттаивания), кристаллизации (плавления), когда теплофизические характеристики, начальные и граничные условия принимаются постоянными. Рассмотрим процесс промерзания грунта. Координатную ось 0х направим вглубь грунта. Пусть начальное распределение температуры постоянно и равно С>0. Если на поверхности х=0 температура мгновенно изменяется и все время поддерживается постоянной, отличной по знаку от начальной температуры, и равной <0, то граница промерзания х=?(t) будет со временем проникать вглубь грунта. Задача о распределении температуры при наличии фазового перехода и о скорости движения границы раздела фаз сводится к решению уравнений:

(1.1)

(1.2)

С дополнительными условиями:

(1.3)

(1.4)

а на границе фазового перехода заданы условия:

(1.5)

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

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

Отсюда, в частности, следует, что движение нулевой изотермы Т=0 будет описываться уравнением:

где ? - значение аргумента, при котором f(?)=0.

Введя новую переменную

вместо задачи (3.10- (3.6) получим новую задачу

(1.7)

(1.8)

(1.9)

Решения обыкновенных дифференциальных уравнений (1.7), (1.8) имеют следующий вид:

Постоянные определяются с помощью условий (1.9):

где - интеграл ошибок.

Условие (1.10) преобразуется в следующее трансцендентное уравнение относительно ?:

(1.11)

Так как при x>0 erf(x)> 0 и при x>? erf(x)> 1, то левая часть (1.11) для при изменении ? от 0 до ? изменяется от -? до +?, а правая часть от 0 до -?. Отсюда следует существование хотя бы одного корня уравнения (1.11). Тогда глубина промерзания определится по формуле:

В связи с тем, что решение трансцендентного уравнения (1.11) представляет некоторые трудности, для ориентировочных расчетов часто применяется формула, известная в литературе как формула Стефана.

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

Интегрируя по t от 0 до некоторого , получим формулу Стефана для определения глубины промерзания

Обобщением этой формулы является формула Лейбензона, которая получается, если распределения температуры в талой и мерзлой зонах задаются в виде:

Очевидно, что выбранные функции удовлетворяют уравнениям (1.1), (1.2) и условиям (1.3), (1.4), а условие Стефана (1.6) преобразуется к виду

Полагая здесь приходим к квадратному уравнению относительно ?:

определив его корень и подставив значение ?, получаем значение глубины промерзания

В частном случае при с=0 получается формула Стефана.

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

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

2. Численные методы, применяемые для решения задачи Стефана

Наиболее универсальным методом решения задач типа Стефана являются численные методы, которые начали разрабатываться с 50-х годов 20 века.

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

2.1 Метод ловли фазового фронта в узел сетки

Метод пригоден для решения одномерных задач типа Стефана. Рассмотрим вначале однофазную задачу типа Стефана. Для простоты примем коэффициенты уравнения теплопроводности равными 1. Пусть требуется определить функцию Т(x,t) и ?(t) из условий

(1.12)

(1.13)

(1.14)

(1.15)

Предположим, что функции µ(t)<0, ?(x)?0, а функция ?(t) монотонно возрастающая. Для численного решения данной задачи применим разностный метод ловли фронта в узел сетки. Характерная особенности этого метода заключается в специальном способе построения разностной сетки. На отрезке [0,l], где l=?(t*), t* - конечный момент времени, строим неравномерную пространственную сетку, состоящую из узлов, так чтобы точка совпадала с одним из узлов

На отрезке [0,t*] также строится неравномерная сетка

Шаг сетки по времени выбирается таким образом, чтобы за каждый шаг по времени фронт фазового перехода перемещался по координате х ровно на один шаг, т.е. ?()-?()=. Задачу (1.12)-(1.15) заменим следующей разностной задачей для определения и

(1.16)

(1.17)

(1.18)

(1.19)

(1.20)

Для простоты мы рассматриваем здесь чисто неявную разностную схему, которая имеет первый порядок аппроксимации по ?, а условие Стефана аппроксимировано разностным уравнением (1.19) с первым порядком по h и ?. Аналогично могут быть рассмотрены разностные схемы более высокого порядка аппроксимации. Индекс «v» над функцией обозначает ее значение на предыдущем моменте времени. Система (1.16)-(1.20) относительно неизвестных , на j-м временном слое представляет собой нелинейную систему алгебраических уравнений и ее решение может быть найдено итерационным методом. Предполагая, что искомые величины для всех k=1,2,...,j-1, i=0,1,...,m+k найдены, будем определять их значение на j-м временном слое по следующей итерационной схеме.

(1.21) (1.22)

(1.23)

Если известно , то из системы (1.21), (1.22) методом прогонки определяются , а из (1.23) - следующая итерация и т.д. Начальную итерацию можно принять за

Конец итерационного процесса проверяется условиями:



где - достаточно малые числа. Теперь рассмотрим применение данного метода к решению двухфазной задачи Стефана, которую мы сформулировали в главе 1:

(1.24)

(1.25)

(1.26)

(1.27)

Предположим, что решения уравнений удовлетворяют следующим начальному и граничным условиям:

Левая часть (1.26) обозначает разность значений функции слева и справа от фронта фазового перехода. Аналогично предыдущему строится разностная сетка с шагами , . Заметим, что узел, лежащий на фронте фазового перехода имеет индексы (j,j). Это будем иметь в виду в дальнейшем при записи разностных уравнений в талой и мерзлой частях области. Для численного решения задачи (1.24) - (1.30) возьмем также чисто неявную разностную схему, которую представим в виде:

(1.31)

(1.32)

(1.33)

(1.34)

(1.35)

Граничному условию (1.30) поставим в соответствие его дискретный аналог, имеющий второй порядок аппроксимации по h:

(1.36)

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

j-ом слое, следующую s+1 итерацию определим из решения системы уравнений

(1.37)

 (1.38)

(1.39)

(1.40)

(1.41)

(1.42)

Если известно , то решая систему (1.37), (1.38), найдем s+1 итерацию в узлах i=1,2,...j-1 и из системы (1.39) - (1.41) определим значения для i=j+1,j+2,...,n. На последнем этапе из уравнения (1.42) найдем следующую s+1 итерацию. Процесс итерации продолжаем до тех пор, пока не выполнится условие

где - заданные достаточно малые числа. Решения систем (1.37), (1.38) и (1.39) - (1.41) легко находятся с помощью метода прогонки. Для однофазной задачи Стефана рассмотрим модификацию метода ловли фронта.

Если учесть монотонное возрастание функции ?(t) в задаче (1.12) - (1.15), то можно построить разностную сетку иначе, чем мы делали выше. Вначале построим фиксированные сетки на отрезках [0, ] и [0,t*] с переменными шагами , :



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

Разностная задача на построенной пространственно-временной сетке будет иметь такой же вид, как (1.16) - (1.20). Но в отличие от предыдущего, в силу выбора шага , численная реализация ее проводится по другому алгоритму. Предположим, что все искомые величины ,, до j-1 временного слоя найдены. Рассмотрим определение их на j-м слое. Система уравнений (1.16) на j-м слое состоит из m+j-1 уравнений и последнее из них ( при i=m+j-1) содержит неизвестный шаг . Поэтому возьмем систему без последнего уравнения и решение ее ищем по методу прогонки

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



Определение по формуле прогонки начинаем с i=m+j-2 до i=0. Для того чтобы начинать вычисления необходимо иметь значение . Это неизвестное значение определяется из совместного решения двух уравнений, первое из которых получается из системы (1.16) при i=m+j-1, а второе - из формулы прогонки при i=m+j-2. Таким путем получим следующее равенство:

(1.43)

При выводе этого равенства учтено равенство нулю значений в узлах, лежащих на фронте фазового перехода. Для определения воспользуемся первым условием из (1.19). Подставив значение из (1.43), получаем относительно искомого шага z= кубическое уравнение

Таким образом, алгоритм модифицированного метода ловли фронта в узел сетки состоит из следующих этапов:

- вычисляются прогоночные коэффициенты ,;

- по формуле вычисляется неизвестный шаг пространственной сетки ;

-на последнем этапе по формуле прогонки определяются значения ,

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

2.2 Метод выпрямления фронтов

Идея метода состоит в преобразовании задачи типа Стефана в области с криволинейной границей, в другую задачу, определенной в прямоугольной области. Рассмотрим применение его к однофазной задаче типа Стефана.

Пусть требуется найти функции T(x,t), ?(t), удовлетворяющие условиям (1.12)-(1.15). Данную задачу преобразуем, введя новую пространственную переменную

Тогда отрезок [0, ?(t)] преобразуется в единичный отрезок [0,1] и область переходит в прямоугольную область . Функция T(x,t)=T(?,t) и с учетом зависимости ? от x и t для производных функции T получаем следующие выражения:

Уравнения (1.12)-(1.15) преобразуются в следующие:

, (1.45)

(1.46)

, t>0 (1.47)

, (1.48)

Заметим, что относительно неизвестных функций T(x,t), ?(t) получили нелинейную систему, хотя первоначальное уравнение (1.12) было линейным. Преимущество такого преобразования заключается в том, что задача теперь решается в прямоугольной области с известными границами. Это дает возможность построить разностную задачу на фиксированной пространственно-временной сетке.

Введем в области (,) сетку, для простоты равномерную:

Задачу (1.45)-(1.48) заменим следующей разностной задачей, записанной в итерационной форме.

(1.50)

(1.51)

Для определения s+1 итерации неизвестных , при известных их значениях на j-1 слое и s-й итерации сначала находим s+1 итерацию из уравнения (1.52), затем из системы (1.49)-(1.51), применяя метод прогонки, определяем s+1 итерацию . Процесс итераций продолжается до тех пор,пока не достигнута заданная точность

где - заданные достаточно малые числа.

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

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

2.3 Метод сглаживания коэффициентов

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

Пусть G есть p-мерная область пространства с границей Г, обозначим через цилиндр с основанием G: ={G [0,t*]}. Боковую поверхность этого цилиндра обозначим через S= {Г [0,t*]}.

Требуется найти функции Т, Ф из следующих условий:

(1.53)

(1.54)

(1.55)

(1.56)

(1.57)

(1.58)

Введем новую функцию - удельное теплосодержание

где ?(х) - дельта-функция Дирака.

Подставив ее в уравнение энергии

получим

 (1.59)

Покажем, что уравнение (1.59) включает уравнения (1.53), (1.54) и условие (1.56) на фазовой поверхности. Первая часть утверждения следует немедленно, если учесть свойство ? - функции ?(x)=0 при х?0. Для доказательства второй части рассмотрим точку Р на поверхности фазового перехода Ф=0 (см.рис.1.1)

Рис

Проведем нормаль к поверхности и касательную в этой точке и построим цилиндр достаточно малого объема с осью совпадающей с нормалью, симметричной относительно касательной и образующей, параллельной к нормали. Пусть - боковая поверхность, - нижнее, - верхнее основание цилиндра; их площадь обозначим соответственно ||, ||=||=|?2|. Проинтегрируем уравнение (1.59) по объему цилиндра и устремим ||>0, а затем и ||>0. Тогда интеграл от С??Т??t в пределе обратится в нуль. Объемные интегралы

где

преобразуем в поверхностные по поперечному сечению цилиндра. Предположим, что для функции Ф(Т)=0 выполнено условие существования обратной функции dФ/dТ>0. Тогда ?(T-)=?(Ф) и ??(Ф)?t=?(Ф)(?Ф??t). Элемент объема dV=d?dn, где d? - элемент площади плоских участков, параллельных. Так как gradФ направлен по нормали к поверхности Ф(Т)=0, то интегрирование по нормали можно заменить интегрированием по Ф; при этом, принимая во внимание следующие соотношения

при ||>0, получим

где - среднее поперечное сечение цилиндра.

Для преобразования воспользуемся формулой Остроградского-Гаусса

при ||>0, интеграл по обратится в нуль, а интегралы по ||,|| перейдут в интеграл по



Подставим сюда выражение

и заменим единичный вектор нормали из равенства

Тогда для получим следующее выражение:

Вычислив значение подынтегрального выражения в точке Р, окончательно имеем

Из уравнения (1.59) следует=, т.е. получаем условие (1.56).

Таким образом, решение задачи типа Стефана (1.53)-(1.58) сводится к решению уравнения (1.59) с дополнительными условиями (1.55), (1.57), (1.58). Левая часть уравнения (1.59) содержит сосредоточенную теплоемкость L??(Т-) на поверхности фазового перехода Т=, т.е. она обращается в нуль при Т?. Теплофизический смысл этого члена заключается в том, что теплота фазового перехода L? выделяется на фазовом фронте. Если построить для этого уравнения разностное уравнение, то коэффициент левой части уравнения (1.59) будет вычисляться в узлах сетки, а меняющийся фронт фазового перехода не всегда совпадает узлом сетки. Отсюда следует, что разностная схема не всегда будет учитывать теплоту фазового перехода, т.е. она не будет обладать свойством консервативности, Возникновение такой ситуации приводит к необходимости сглаживания коэффициентов уравнения (1.59). Для этого дельта-функция приближенно заменяется дельтообразной, или размазанной, дельта-функцией ?(Т-,?)?0, где ? - величина полуинтервала, на котором отлична от нуля ?(Т-,?).

Таким образом, вводится сглаженная или эффективная теплоемкость

?)

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

Т<-?,

Т>+?

2. Изменение энтальпии на интервале (-?,+?) сохраняется

(+?)- (-?)=H(+?)-H(-?)

На этом же интервале (-?,+?) проводится сглаживание коэффициента теплопроводности ?.

В результате вместо задачи (1.59), (1.55), (1.57), (1.58) получается задача для уравнения теплопроводности со сглаженными коэффициентами

(1.60)

(1.61)

(1.62)

На практике сглаженную теплоемкость выбирают по разному.

Рассмотрим некоторые часто применяемые варианты сглаживания.

Пусть коэффициенты не зависят от Т.

1. Сглаженная теплоемкость на интервале сглаживания постоянна. Тогда из второго условия следует

2.

2?+/2

и уравнение (1.60) решается с коэффициентом

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

Подставив во второе условие, имеем

2b?=+(+)?

а равенство скачков дает

(-?)+= (+?)+

Таким образом, для этого варианта сглаживания уравнение (1.60) имеет следующий коэффициент:

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

2.4 О выборе параметра сглаживания

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

Коэффициенты разностного уравнения вычисляются в зависимости от интервала температуры по разным выражениям. Теплота фазового перехода, выделяющаяся на фронте, учитывается в выражении , определенном в интервале сглаживания. Если в указанном интервале не попадает ни одно узловое значение температуры, то температурное поле будет определено без учета выделения теплоты на фронте. Из сказанного следует условие выбора длины интервала сглаживания: параметр сглаживания ? должен быть выбран так, чтобы интервал сглаживания содержал в себе на каждом шаге по времени интервал, определяемый значениями температуры хотя бы в двух соседних узлах, между которыми находится фронт фазового перехода. Если фронт находится между узлами , то параметр ? должен удовлетворять условию

2? ? ||

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

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

С?=0,84*, L?=0,334*

На первой строке расположены значения температуры в узлах сетки при t=16,5 ч, на второй - при t=21,5 ч после начала замораживания

Рис

В первых двух вариантах ?=1. В первом варианте при наличии свободной воды решение имеет колебательный характер во времени. Это показывает, что тепловыделение на фронте при выбранной величине ? учитывается в расчетах не во всех моментах времени: разность значений температуры в двух соседних узлах намного превосходит 2?, т.е. указанное условие выбора параметра ? не выполнено. Когда отсутствует свободная вода (вариант 2), в выражении коэффициента слагаемое, содержащее теплоту фазового перехода обращается в нуль и решение задачи при любом ? имеет монотонный характер во времени. В вариантах 3 и 4 величина тепловыделения за счет фазового перехода воды превосходит на порядок коэффициент удельной теплоемкости материала. Когда параметр сглаживания ?=3, происходит незначительное колебание значений температуры во времени. При большой величине ? (?=10) эти колебания исчезают и идет монотонное промерзание. Таким образом, правильный учет тепловыделения на фронте фазового перехода зависит от выбора длины интервала сглаживания.

2.5 Разностные схемы сквозного счета

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

(1.63)

(1.64)

 (1.66)

(1.67)

(1.68)

(1.69)

На отрезке [0,l] введем квазиравномерную сетку основных и потоковых узлов:

область [0,l] потоковыми узлами разбивается на ячейки

Введем сетку по времени, также неравномерную

Для построения разностной задачи уравнения (1.63), (1.64) при t= проинтегрируем на каждой ячейке i с учетом условий (1.66) - (1.69):



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

(1.70)

Система (1.70) нелинейная и ее решение проводится методом итераций.

Итерационный процесс решения системы (1.70), состоит из следующих этапов (для сокращения записи итерационную схему не приводим):

1. Задаем начальное приближение , например, по формуле

Здесь при ?=0 фронт неподвижен, а при ?=1 фронт продвигается по линейному закону.

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

3. По найденной итерации вычисляется следующее приближение Процесс итерации повторяется до тех пор, пока не выполнятся условия

где заданные малые положительные числа.

Усвоив идею метода из изложенного материала, распространение данного метода на другие постановки задач типа Стефана читатель сможет выполнить самостоятельно.

Заключение

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

Список используемой литературы

1.А. Д. Полянин, В. Ф. Зайцев, Справочник по нелинейным уравнениям математической физики, М.: ФИЗМАТЛИТ, 2002.

2. Б. М. Будак, А. А. Самарский А. Н. Тихонов, Сборник задач по математической физике-2001г.

3. Пикулин В.В., Похожаев С.И Практический курс по уравнениям математической физики-2002г.

4. Гладков С.О, Сборник задач по теоретической и математической физике-2001г.

5. Несис Е.И.Методы математической физики-2003г.

6. Тихонов А.Н., Самарский А.А., Уравнения математической физики-2002г.

Размещено на Allbest.ru