Математическое моделирование свободной конвекции несжимаемой жидкости в двумерных областях с фиксированными и подвижными границами тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат физико-математических наук Чеблакова, Елена Анатольевна

• Чеблакова, Елена Анатольевна
• кандидат физико-математических наук
• 2000

Оглавление диссертации кандидат физико-математических наук Чеблакова, Елена Анатольевна

Обозначения

Введение.

Глава 1. Разностные схемы для решения уравнений тепловой конвекции.

§ 1. Уравнения тепловой конвекции в области с неподвижными границами.

1.1. Уравнения в приближении Обербека-Буссинеска.

1.2. Уравнения в переменных "функция тока - вихрь скорости".

1.3. Безразмерный вид уравнений тепловой конвекции в переменных у/, со.

1.4. Формулы Тома и Вудса приближенного граничного условия для вихря на стенке.

§2. Уравнения тепловой конвекции в области со свободной границей.

2.1. Кинематическое и динамические граничные условия на свободной поверхности.

2.2. Условия для вихря и функции тока на свободной поверхности.

2.3. Уравнение для касательной составляющей скорости на свободной поверхности.

§3. Метод экспоненциальной подгонки для уравнений тепловой конвекции.

3.1. Аппроксимация дифференциальных операторов.

3.2. Модифицированный метод неполной факторизации Булеева и схемы расщепления.

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

Выводы по главе 1.

Глава 2. Исследование конвекции и переноса тепла в жидкости со свободной покоящейся границей.„.

§ 1. Физико-математическая модель.

1.1. Уравнения движения.

1.2. Постановка краевых условий.

§2. Метод решения.

§3. Анализ результатов.

Выводы по главе 2.

Глава 3. Исследование конвекции и переноса тепла в жидкости с криволинейными подвижными границами.

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

§2. Исследование конвекции и переноса тепла в жидкости с подвижной свободной границей жидкость - газ.

2.1. Физико-математическая модель.

2.1.1. Уравнения движения.

2.1.2. Постановка начальных и краевых условий задачи.

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

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

§3. Решение двумерной задачи Стефана о фазовом переходе.

3.1. Физико-математическая модель.

3.1.1. Уравнения движения.

3.1.2. Постановка начальных и краевых условий задачи.

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

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

Выводы по главе 3.

Введение диссертации (часть автореферата) на тему «Математическое моделирование свободной конвекции несжимаемой жидкости в двумерных областях с фиксированными и подвижными границами»

Понятие конвекции произошло от латинского слова сотесйо — принесение, доставка; и обозначает перемещение макроскопических частей среды (газа, жидкости), приводящее к переносу массы, теплоты и других физических величин [1]. Различают естественную (свободную) конвекцию, вызванную неоднородностью среды (градиентами температуры и плотности), и вынужденную конвекцию, вызванную внешними механическими воздействиями на среду. Данная работа посвящена исследованию свободной конвекции, вызванной неоднородностями плотности за счет перепадов температуры, с помощью методов математического моделирования. Предполагается, что жидкость вязкая и несжимаемая, что ее плотность р и коэффициент поверхностного натяжения <т зависят только от температуры и не зависят от давления; что другие физические свойства вещества (теплопроводность, теплоемкость, вязкость) постоянны для рассматриваемого объема жидкости. Эти предположения наряду с предположением о малости отклонений системы от состояния гидростатического равновесия составляют основу приближения Обербека-Буссинеска [2], [3]. Несмотря на эти предположения, уравнения в приближении Обербека-Буссинеска хорошо описывают широкий спектр явлений свободной конвекции [3], [4].

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

Задачи расчета конвективных течений возникают при рассмотрении многих прикладных задач. Это, например, задачи:

• динамики ледяного покрова водоемов и водотоков [5] - [7];

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

Чохральского, Бриджмена, бестигельной зонной плавки (БЗП) [8] -[32];

• исследования течений в расплавах при пониженной гравитации [33]

42];

• расчета поля температур в многослойной среде с разной теплопроводностью слоев (например, нефть, грунт, вода) [10], [15], [19] -[26], [32];

• моделирования процессов глубинной геодинамики, таких как образование рельефа Земли, движение материков, формирование структуры недр [43] - [45]. Конвективные течения в этих задачах вызываются эффектами плавучести (конвекция Грасгофа) и термокапиллярности (конвекция Марангони). В большинстве ситуаций термокапиллярная конвекция достаточно мала по сравнению с естественной и часто опускается, но ее вклад резко возрастает при слабой гравитации (например, на борту космических станций) или при расчетах, проводимых для моделирования роста кристаллов. Поэтому в таких задачах необходимо исследовать влияние обоих типов конвекции на распределение температур и скоростей течения.

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

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

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

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

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

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

При моделировании роста кристаллов необходимо численно решать уравнения тепловой конвекции. Основной трудностью при решении этих уравнений для процессов кристаллизации является наличие большого параметра в правой части уравнения для вихря. Таким параметром является число Релея или число Грасгофа, которое из физической постановки задачи может достигать значения 106 и более. Если поделить все части уравнения на слагаемое, содержащее число Релея или Грасгофа, то мы получим уравнение с малым параметром при старшей производной. Такие уравнения имеют как быстро, так и медленно меняющиеся в пространстве компоненты [46]. При их интегрировании быстро меняющиеся компоненты решений могут привести к возникновению трудных вычислительных задач, т.к. в этом случае возникают пограничные слои или внутренние переходные слои с большими градиентами искомых физических величин. Это приводит к снижению порядка сходимости и неравномерной сходимости на равномерных сетках. Для преодоления описанных сложностей используют различные подходы: повышают порядок точности разностной схемы; в области быстрого изменения решения используют сетку с мелкими шагами. Однако эффективность этих способов уменьшается с увеличением жесткости уравнения. Другой подход состоит в применении метода экспоненциальной подгонки (или метода интегральных тождеств) [46] - [50]. Его суть состоит в том, чтобы в аппроксимирующую разностную схему ввести подгоночные коэффициенты, которые обеспечивали бы равномерную сходимость на равномерной сетке.

В настоящей работе рассматриваются сеточные операторы, аппроксимирующие исходные дифференциальные операторы первой и второй производных на трехточечных шаблонах узлов сетки. В практике расчетов к настоящему времени нашли применение три таких оператора [51]. Два из них получены при раздельном приближении первой и второй производных. При этом первая производная аппроксимируется направленными против потока разностями и центральными разностями; а вторая производная приближается традиционным способом на трехточечном шаблоне [52], [53], [100]. Третий способ состоит в совместной аппроксимации первой и второй производных (подробнее алгоритм можно посмотреть в §3 гл.1). Прежде всего уравнение приводится к дивергентному виду и вводится понятие суммарного потока, который складывается из конвективного и диффузионного потоков. Расчетная область разбивается на некоторое число непересекающихся контрольных объемов так, что каждая узловая точка содержится в одном из них. Дифференциальное уравнение интегрируется по каждому контрольному объему. Для вычисления интегралов используют кусочные профили, которые описывают изменение соответствующих функций между узловыми точками. В результате получается дискретный аналог дифференциального уравнения, в который входят значения искомой функции в нескольких узловых точках. Полученный подобным образом дискретный аналог выражает закон сохранения для конечного контрольного объема точно так же, как дифференциальное уравнение выражает закон сохранения для бесконечно малого контрольного объема. Одним из важных свойств метода интегральных тождеств является то, что в нем заложено точное интегральное сохранение таких величин, как масса, количество движения и энергия на любой группе контрольных объемов и, следовательно, на всей расчетной области. Это свойство проявляется при любом числе узловых точек, а не только в предельном случае очень большого их числа. Таким образом, даже решение на грубой сетке удовлетворяет точным интегральным балансам [50].

В результате применения метода интегральных тождеств получаются конечно-разностные уравнения с экспоненциальными выражениями для коэффициентов. Из анализа коэффициентов следует, что полученная схема монотонна [54]. В работах [47], [49] такая аппроксимация была получена для уравнения диффузии. Впервые в литературе подобная аппроксимация встречается в работе Аллена Д.Н. и Саусвелла Р.В. [55], где она была получена иным путем. В работе [55] с помощью конформного отображения осуществляется преобразование координат. При этом конечно-разностное уравнение для вихря в новых переменных содержит экспоненциальные коэффициенты. С помощью этой техники определяется установившееся ламинарное течение вязкой несжимаемой жидкости за неподвижным цилиндром.

Если дифференциальные уравнения аппроксимировать с помощью центральных и направленных разностей, то решения соответствующих разностных уравнений сходятся к решению дифференциального уравнения при условии, что шаг сетки стремится к нулю. В часто встречающихся задачах, где параметр при второй производной (в общем случае, при старшей производной) е достаточно мал, для достижения необходимой точности требуется чрезвычайно малая (практически недостижимая) величина шага сетки [46]. Известно [56], что замена первых производных центральными разностями, которая дает достаточно хорошую аппроксимацию второго порядка, непригодна для уравнений с малым параметром при старшей производной, если оператор Лапласа приближать обычным способом (на трехточечном шаблоне равномерной сетки). В [56] показано, что в этом случае при е «к (к - шаг сетки) решение разностного уравнения может не иметь ничего общего с решением дифференциального уравнения. Известный способ избавления от подобных трудностей состоит в использовании направленных против потока разностей для аппроксимации первой производной. Такой способ дает возможность численно решать краевую задачу для любых малых 8. Однако, как показано в [56], он также обладает существенными недостатками. Односторонние разности обеспечивают аппроксимацию лишь первого порядка на равномерной сетке. При этом для значений е, соизмеримых с к, появляется так называемая схемная вязкость, которая может превысить вязкость, описываемую дифференциальным уравнением. И, наконец, разностное уравнение с направленными против потока разностями плохо описывает пограничный слой - он оказывается размытым.

Поэтому целесообразно строить такие разностные схемы, решение которых было бы близко к решению исходной задачи для всех достаточно малых с. Пограничный слой описывается решениями экспоненциального типа, поэтому естественно предположить, что схема с экспоненциальными коэффициентами разностного уравнения будет хорошо описывать течение в пограничном слое. На примере одномерного уравнения для вихря в [51] показано, что экспоненциальная схема (при различных числах Рейнольдса Ке=1; 10; 100) привела к более точному решению, чем схемы с центральными и направленными разностями. Использование схемы с центральными разностями при Ке=100 привело к неустойчивости счета, что связано с невыполнением в этом случае условия диагонального преобладания [52]. Для уравнений, полученных методом экспоненциальной подгонки, условия диагонального преобладания выполняются всегда [51].

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

Вопросам построения разностных схем для решения уравнений с малым параметром при старшей производной посвящены также работы [57]-[67].

В данной работе метод экспоненциальной подгонки применяется для численного решения уравнений тепловой конвекции в задачах о течении жидкости в замкнутой квадратной полости при боковом подогреве, нижнем подогреве, а также в задаче о течении жидкости в квадратной полости с верхней свободной неподвижной границей при подогреве сбоку. Были проведены расчеты как со стационарными уравнениями тепловой конвекции, так и с нестационарными, которые решались методом установления [52]. Результаты расчетов согласуются с результатами [51], [54] в том смысле, что метод экспоненциальной подгонки позволяет получить более высокую точность по сравнению с традиционной аппроксимацией центральными и направленными против потока разностями. Также было обнаружено, что при решении нестационарных уравнений тепловой конвекции методом экспоненциальной подгонки можно брать достаточно большие шаги по времени, чего не позволяют обычные способы раздельной аппроксимации конвективных и диссипативных слагаемых (см. например, работы [4], [10], [11], [13], [35], [68] - [78]). Это объясняется большей устойчивостью экспоненциальной схемы. Автором данной работы также была проведена оценка порядка аппроксимации экспоненциальной схемы, которая показала, что на равномерной сетке схема экспоненциальной подгонки имеет 2 порядок аппроксимации.

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

1) модифицированный метод неполной факторизации Булеева [79];

2) схема продольно-поперечной прогонки [80];

3) схема стабилизирующей поправки [80].

Расчеты показали, что скорость сходимости модифицированного метода Булеева в 2-3 раза выше скорости сходимости методов расщепления, в зависимости от значений безразмерных параметров. Таким образом, был построен вычислительный алгоритм, состоящий из экспоненциальной схемы для аппроксимации дифференциальных уравнений и модифицированного метода факторизации для решения системы конечно-разностных уравнений, который обладает большей скоростью сходимости и большей устойчивостью по сравнению со схемами расщепления и применением раздельной аппроксимации конвективных и диссипативных слагаемых. Разработанный алгоритм реализован в пакет программ. Проведена серия тестовых расчетов на основе международного теста о течении вязкой несжимаемой жидкости в подогреваемой сбоку замкнутой квадратной полости ([81], [82]). На основании проведенных тестов был сделан вывод, что предложенный алгоритм обеспечивает получение достоверных результатов.

С помощью этого алгоритма была решена задача о конвективном течении жидкости в нагреваемой сбоку квадратной полости с верхней неподвижной свободной границей. Исследовано влияние гравитационной и термокапиллярной конвекции на поле температур и скоростей течения, начиная с условий отсутствия гравитации и затем при её увеличении. Актуальность решения задач о конвективных течениях в условиях пониженной гравитации возросла в последнее время. Развитие ракетнокосмической техники, сделавшей доступной длительную невесомость, дало возможность воздействовать на вещество в условиях малой гравитации, когда ускорение силы тяжести в тысячу и более раз меньше, чем на Земле. Особенно эффективным использование невесомости может быть в процессах, сопровождающихся фазовыми переходами. К этим процессам относится, например, получение монокристаллов и различных композиционных материалов[17].

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

В данной работе при моделировании стационарного течения жидкости в подогреваемой сбоку квадратной полости с верхней свободной границей предполагается, что относительное изменение поверхностного натяжения жидкости очень мало (т.е. капиллярное число стремится к нулю), так что свободная поверхность остается плоской. Подобную задачу для условий невесомости также решали Бердников B.C. и Гапонов В.А. [33]; Зебиб А. и др. [37]; Карпентер Б.М. и Хомси Г.М. [38]. В этих работах показано, что вдоль свободной поверхности образуется пограничный слой, а основные неустойчивости течения связаны с сильными потоками в углах полости. Результаты настоящей работы находятся в хорошем согласии с результатами [33], [37], [38]. В работах [33], [37], [38] представлена достаточно полная теория пограничного слоя Марангони для рассматриваемого течения. В данной работе приведено дальнейшее развитие этой теории. Показано, что в промежутке значений параметра Марангони 0<Ма<102 гравитационная конвекция становится преобладающей, начиная с числа Релея Ка^Ю4. Также найдены диапазоны изменения безразмерных параметров Марангони и Прандтля, в которых течение не меняется в условиях отсутствия силы тяжести. Были рассмотрены два предельных случая: Ма—>0 и /V—>оо в условиях невесомости и получены численные решения для этих случаев. Используя метод наименьших квадратов, были получены аппроксимационные выражения для чисел Нуссельта вдоль горячей стенки как функции от числа Релея при 10^<Ка<106 и от числа Марангони при 10:,<Ма<104. Построены зависимости интегральных чисел Нуссельта вдоль стенок от чисел Прандтля, Марангони и Релея, и тем самым изучена зависимость теплоотдачи от всех безразмерных параметров.

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

Zl. Расчет течения жидкости со свободной границей жидкость/газ при боковом подогреве в условиях невесомости;

Z2. Задача Стефана о плавлении твердого вещества при подогреве снизу [12], [83].

Эти задачи нелинейные; в них неизвестная граница определяется по неизвестному решению, и, наоборот, решение ищется в неизвестной заранее области. Перечисленные задачи являются составными частями задачи моделирования роста кристаллов методом бестигельной зонной плавки (БЗП). 4

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

1 с в т

81 (х)

А о х

Рис. 1. Эскиз к вертикальному варианту бестигелъной зонной плавки (слева) и преобразование координат. удерживается между твердыми частями образца. Градиент температуры на' свободной поверхности жидкой фазы меняет поверхностное натяжение, что в свою очередь приводит к термокапиллярной конвекции расплава. На рис.1 введены следующие обозначения: / - свободная поверхность жидкость/газ, gl к g2 - границы фазового перехода жидкость/твердое тело, А и С - области твердой фазы вещества, В - область жидкой фазы. Область А вращается с угловой скоростью со/ в одну сторону, область С вращается с угловой скоростью со2 в другую сторону, и вся система медленно опускается вниз. Кроме того, на расплав (область В) действует магнитная индукция, способствующая лучшему прогреву расплава и его поддержке.

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

1. Методы, в которых граница фазового перехода явно не выделяется. Это методы фиктивных областей (МФО) и основанные на них алгоритмы сквозного счета соответствующих уравнений с разрывными коэффициентами. Алгоритм был предложен Самарским A.A., Моисеенко Б.Д. [5] и развит позднее многими авторами: [4], [15], [43] и т.д. К этому типу решения задачи Стефана также можно отнести так называемый энтальпийный подход, в котором нестационарное уравнение теплопроводности записывается и решается в переменных энтальпии. Все теплофизические свойства считаются известными функциями энтальпии. Если они оказываются недостаточно гладкими, то они искусственно сглаживаются. Этот подход используется, например, в работах [20] - [22].

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

2.1. К первому относятся методы выпрямления фронтов ([10], [23], [25], [30]), использующие преобразования неортогональных областей в прямоугольные, при которых границы фазовых переходов становятся координатными линиями (т.е. выпрямляются). При этом уравнения движения принимают в новых переменных более громоздкий вид, появляются смешанные производные, однако все вычисления проводятся в прямоугольных областях.

2.2. Ко второму относятся методы, использующие разностные аппроксимации на неортогональных подвижных сетках ([26]), и методы, основанные на применении конечных элементов ([36]). В этом случае не надо выполнять преобразование уравнений движения к новым переменным, однако все вычисления проводятся в областях с криволинейными границами.

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

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

Метод БЗП численно моделировался в следующих работах. Анисютин Б.М. [27] и Кобаяши Н. [29] исследовали вопрос о зависимости размеров расплавленной зоны от мощности источника тепла для осесимметричной задачи. В [29] существенным ограничением является предположение о цилиндричности боковой поверхности. В [27] численно исследуется осесимметричная стационарная задача Стефана со свободной боковой поверхностью. В обеих работах [27], [29] предполагается отсутствие конвекции в расплаве, т.е. строится модель без гидродинамических эффектов. Дональд Д.К. [28] предложил одномерную модель зонной плавки, обладающую аналитическим решением. Лан Ц.В. и Коу С. [30] исследовали эффекты сильных вращений кристалла и затравки на перенос тепла, течение жидкости и форму поверхности раздела в процессе БЗП.

Большой трудностью при моделировании БЗП, а также течений жидкости со свободной границей жидкость/газ является то, что формулировка рассматриваемых при этом начально-краевых задач не содержит в явном виде условий на свободной границе для искомых функций (функции тока и вихря скорости; а также скоростей и давления). Граничные условия заданы в виде некоторых соотношений, выражающих непрерывность нормальной и касательной компонент вектора напряжений (так называемые динамические условия на свободной поверхности). В связи с этим возникают трудности с определением формы свободной поверхности. Некоторые авторы (например, Кобаяши Н. [29]) в своих расчетах делают предположение о цилиндричности боковой поверхности. В других работах (например, [27]) вычисления проводятся без гидродинамических эффектов, а свободная поверхность вычисляется как для неподвижной жидкости. В [7] изучается тепло- и массоперенос за счет поверхностного натяжения в расплаве в предположении, что поток, вызванный поверхностным натяжением, слаб. Авторами [31] проведено более полное моделирование процесса БЗП. Однако там также предполагается, что естественная конвекция и вращение преобладают над термокапиллярной конвекцией. Для функции тока на свободной поверхности ставится условие ¡¡/=0. Таким образом, выходит из рассмотрения кинематическое условие на свободной поверхности. В этом отношении наиболее естественным можно считать подход Овчаровой A.C. ([11], [35]), в котором на основе приемов, описанных в [84], разработан алгоритм, позволяющий получить в явном виде граничные условия для функции тока и вихря на свободной поверхности, и предложен метЬд расчета конвективных течений жидкости со свободной границей жидкость/газ в переменных у/ - со. В работе [11] численно решена плоская задача Стефана в области со свободной границей жидкость/газ. При этом осуществляется преобразование одной пространственной координаты для выпрямления свободной границы жидкость/газ. Граница фазового перехода явно не выделяется, а находится методом сквозного счета.

Представляет интерес построить модель, в которой осуществляется преобразование обеих пространственных координат для выпрямления как границы фазового перехода, так и свободной границы жидкость/газ. В данной работе сделано такое преобразование координат и получены уравнения тепловой конвекции в новых переменных. На основе полученных уравнений построена математическая модель для решения задач тепловой конвекции в областях со сложной геометрией, и численно решены перечисленные выше задачи Zl и Z2, являющиеся составными частями задачи моделирования роста кристаллов методом БЗП. Результаты расчетов сравнивались с результатами [11], [36]; а также [4], [10], [15], где подобные задачи решены другими методами, либо с использованием более простого преобразования координат. Это сравнение показало достоверность построенной модели. Модель может быть использована в дальнейшем для моделирования роста кристаллов методом бестигельной зонной плавки.

Цель диссертационной работы заключалась в:

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

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

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

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

• результаты численного моделирования течения в подогреваемой сбоку полости с неподвижной свободной границей жидкость/газ в условиях невесомости и при увеличении силы гравитации с помощью построенного вычислительного алгоритма;

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

Диссертация состоит из введения, трех глав, заключения, серии численных результатов в виде 4 таблиц и 29 рисунков, и списка цитированной литературы из 115 наименований. В первой главе приводятся основы математического аппарата, необходимого для моделирования процессов конвекции. Выписаны уравнения движения вязкой несжимаемой жидкости в приближении Обербека-Буссинеска; рассмотрен переход от естественных переменных и, v, р (скорость, давление) к переменным функция тока, вихрь скорости; приведены уравнения в безразмерном виде. В связи с отсутствием граничного условия для вихря на границе, являющейся неподвижной твердой стенкой, рассмотрены различные методы построения приближенного граничного

Заключение диссертации по теме «Математическое моделирование, численные методы и комплексы программ», Чеблакова, Елена Анатольевна

Выводы по главе 3.

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

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

2.1. Задача Стефана в области, состоящей из двух фаз: жидкой и твердой.

2.2. Задача о конвективном течении жидкости в полости со свободной подвижной границей жидкость/газ.

ЗАКЛЮЧЕНИЕ

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

2. С применением построенного алгоритма проведено численное моделирование течения вязкой несжимаемой жидкости в нагреваемой сбоку квадратной полости с верхней свободной неподвижной границей жидкость/газ в условиях невесомости и при увеличении влияния гравитации. Исследовано влияние гравитационной и термокапиллярной конвекции на поле температур и скоростей в широком диапазоне изменения безразмерных параметров - чисел Марангони: 0<Ма<104, чисел Релея: 0<Ра<106 и чисел Прандтля: 1<Рг<со. Тем самым осуществлено дальнейшее развитие теории пограничного слоя Марангони. ' о

• Показано, что в промежутке значений параметра Ма 0<Ма<10* гравитационная конвекция становится преобладающей, начиная с Яа=104.

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

• С использованием метода наименьших квадратов получены аппроксимационные выражения для чисел Нуссельта вдоль горячей стенки, как функции от числа Релея 102<Ка<106 и от числа Марангони

103 <Ма<104.

• Построены зависимости интегральных чисел Нуссельта вдоль стенок от чисел Прандля, Марангони и Релея, и тем самым изучена зависимость теплоотдачи от всех безразмерных параметров задачи.

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

3.1. Задача Стефана в области, состоящей из двух фаз: жидкой и твердой.

3.2. Задача о конвективном течении жидкости в полости со свободной подвижной границей жидкость/газ.

Список литературы диссертационного исследования кандидат физико-математических наук Чеблакова, Елена Анатольевна, 2000 год

1. Советский энциклопедический словарь. // М. "Советская энциклопедия", 1988.

2. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика. Том 6. Гидродинамика. //М. Наука. 1986.

3. Гершуни Г.З., Жуховицкий Е.М., Непомнящий A.A. Устойчивость конвективных течений. //М. Наука. 1989.

4. Тарунин Е.Л. Вычислительный эксперимент в задачах свободной конвекции. //Иркутск: Изд-во Иркут. Ун-та, 1990.

5. Самарский A.A., Моисеенко Б.Д. Экономичная схема сквозного счета для многомерной задачи Стефана. // ЖВМ и МФ, 1965, т.5, N5, с.816-827.

6. Соловьева E.H., Успенский А.Б. Схемы сквозного счета численного решения краевых задач с неизвестными границами для одномерных уравнений параболического типа. В ich.: Методы решения краевых и обратных задач теплопроводности. МГУ, 1975, с. 3-23.

7. Пригожин Л.Б., Булгач A.A. Численное решение одномерных задач Стефана в теплопроводности и диффузии. // Численные методы механики сплошных сред, 1981, т.12, N2, с. 71-83.

8. Полежаев В.И. Гидродинамика, тепло- и массообмен при. росте кристаллов. // В кн. Авдуевский B.C., Полежаев В.И. Гидромеханика и тепло-массообмен при получении материалов. М. Наука. 1990.

9. Полежаев В.И, Буне A.B., Верезуб H.A. и др. Математическое моделирование конвективного тепломассообмена на основе уравнений Навье-Стокса. //М. Наука. 1987.

10. Ю.Овчарова A.C. Метод решения двумерной многофронтовой задачи Стефана. // ПМТФ, 1995, т.36, N4, с. 110-119.

11. Н.Овчарова А.С. Численное решение стационарной задачи Стефана в области со свободной границей. // Вычислительные Технологии, 1999, т.4, N1, с.88-99.

12. Мейрманов A.M. Задача Стефана. // Наука, Новосибирск, 1986.

13. Пасконов В.М., Полежаев В.И., Чудов JI.A. Численное моделирование процессов тепло- и массообмена. // М. Наука, 1984.

14. Scharfetter D., Gummel Н. Large signal analysis of a silicon Read diode oscillator. // IEEE. Trans., 1969, ED-16, pp.64-77.

15. И.Вабищевич П.Н., Илиев О.П. Численное решение сопряженных задач тепло- и массопереноса с учетом фазового перехода. // Диф. ур., 1987, т.23, N7, с. 1127-1132.

16. Лодиз Р., Паркер Р. Рост монокристаллов. // М. Мир, 1974.

17. Плавление, кристаллизация и фазообразование в невесомости. // М. Наука, 1979.

18. Пак Н.И., Шикунов С.А. Численное решение двумерной задачи Стефана в энтальпийной формулировке в изменяющейся области на подвижных сетках. // Сб. Сильновозбужденные состояния в кристаллах, Томск: издание Томского научного центра СО АН СССР, 1991.

19. Bonnerot R., Jamet P. Numerical computation of the free boundary for the 2-dimensional Stefan problem by space-time finite elements. // J. of Computational Physics, 1977, v.25, p. 163-181.

20. Crowley A.B. Numerical solution of Stefan problems. // Int. J. Heat Mass Transfer, 1978, v.21, p.215-219.

21. Bonacina C., Comini G., Fasano A., Primicerio M. Numerical solution of phase-change problems. // Int. J. Heat Mass Transfer, 1973, v. 16, p. 18251832.

22. Solomon A. Some remarks on the Stefan problem. // Math. Сотр., 1966, v.20, p.347-360.

23. Успенский А.Б. О методе выпрямления фронтов для многофронтовых одномерных задач типа Стефана. // Доклады Академии Наук СССР, 1967, т.172, N1, с.61-64.

24. Журавлева Е.Н. Численное решение задач плавления и кристаллизации бинарного сплава. // Динамика сплошной среды. Сб. Научных трудов, вып. 113. Новосибирск, 1998, с.70-72.

25. Будак Б.М., Успенский А.Б. Разностный метод с выпрямлением фронтов для решения задач типа Стефана. // ЖВМ и МФ, 1969, т.9, N6, с.1299-1315.

26. Бреславский П.В., Мажукин В.И. Алгоритм численного решения гидродинамического варианта задачи Стефана при помощи динамически адаптирующихся сеток. // Мат. моделирование, 1991, т.З, N10, с.104-115.

27. Анисютин Б.М. Численное исследование тепловой задачи для процесса бестигельной зонной плавки. // Задачи гидромеханики и тепломассобмена со свободными границами. Новосибирск, 1987.

28. Дональд Д.К. Тепловой режим в условиях вакуумной плавки. // Приборы для научных исследований. 1961, N7, С.42-44.

29. Kobayashi N. Power required to form a floating zone and the zone shape. // J. Crystal Growth, 1978, v.43, p.417-424.

30. Lan C.W., Kou S. Effects of rotation on heat transfer, fluid flow and interfaces in normal gravity floating-zone crystal growth. // J. Crystal Growth, 1991, v.l 14, p.517-535.

31. Muhlbauer A., Muiznieks A., Virbulis J., Ludge A., Riemann H. Interface shape, heat transfer and fluid flow in the floating zone growth of large silicon crystals with the needle-eye technique. // J. Crystal Growth, 1995, v. 151, p.66-79.

32. Чеблакова Е.А. Численное решение плоской задачи Стефана в области со свободной границей. // Четвертый Сибирский конгресс по прикладной и индустриальной математике (INPRIM-2000). Тезисы докладов. Новосибирск, 2000.

33. Бердников B.C., Талонов В.А. Тепловая гравитационно-капиллярная конвекция в прямоугольных полостях. // Процессы тепломассопереноса и рост монокристаллов и тонкопленочных структур. Труды II Российского симпозиума (HT&CG'97). Обнинск, 1998.

34. Chippada S., Jue Т.С., Ramaswamy D. Finite element simulation of combined buoyancy and thermocapillary driven convection in open cavities. // Int. J. Num. Meth. Eng., 1995, vol.38, pp.335-351.

35. Zebib A., Homsy G.M., Meiburg E. High Marangoni number convection in a square cavity. // Phys. Fluids, 1985, vol.28, pp.3467-3476.

36. Carpenter B.M., Homsy G.M. High Marangoni number convection in a square cavity: Part II. // Phys. Fluids A, 1990, vol.2, pp.137-149.

37. Гадияк Г.В., Чеблакова Е.А. Моделирование конвекции и переноса тепла в жидкости со свободной границей в условиях пониженной гравитации. // Математические модели и методы их исследования. Международная конференция. Тезисы докладов. Красноярск, 1999.

38. Бабский В.Г., Копачевский Н.Д., Мышкис А.Д. и др. Гидромеханика невесомости. // М. Наука, 1976.

39. Волков П.К., Захаров Б.Г. Гидродинамика и конвективный тепломассоперенос в условиях слабой гравитации. // Доклады Академии Наук, 1998, т.361, N5, с. 616-619.

40. Коробицына Ж.Л., Тычков С.А. Численное моделирование процессов тепло- и массопереноса с учетом фазового перехода в геодинамике. // ЖВМ и МФ, 1997, т.37, N6, с. 733-741.

41. Spohn Т., Schubert G. Convective thinning og the litosphere: a mechanism for rifting and mid-plate volcanism on Earth, Venus and Mars. // Tectonophysics, 1983, v.94, p.67-90.

42. Трубицын В.П., Белавина Ю.Ф., Рыков B.B. Тепловое и механическое взаимодействие мантии с континентальной литосферой. // Физ. Земли, 1993, N11, с.3-15.

43. Дулан Э., Миллер Дж., Шилдерс У. Равномерные численные методы решения задач с пограничным слоем. // М. Мир, 1983.

44. Марчук Г.И. Методы расчета ядерных реакторов. // М. Атомиздат, 1961.

45. Марчук Г.И. Методы вычислительной математики. // Новосибирск. Наука. 1973.

46. Тихонов А.Н., Самарский A.A. Об одной наилучшей однородной разностной схеме. // Доклады Академии Наук СССР, 1959, т. 124, N4, с.779-782.

47. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. // М. Энергоатомиздат, 1984.

48. Булеев Н.И., Тимухин Г.И. О численном решении уравнений гидродинамики для плоского потока вязкой несжимаемой жидкости. // Известия СО АН СССР, 1969, N3, вып.1, с. 14-24.

49. Самарский A.A., Николаев Е.С. Методы решения сеточных уравнений. //М. Наука. 1978.

50. Пейре Р., Тейлор Т.Д. Вычислительные методы в задачах механики жидкости. //Л. Гидрометеоиздат. 1986.

51. Булеев Н.И., Тимухин Г.И. О составлении разностных уравнений гидродинамики вязкой несжимаемой жидкости. // Численные методы механики сплошных сред, 1972, т.З, N4, с. 19-26.

52. Allen D.N., Southwell R.V. Relaxation methods applied to determine the motion, in two dimensions, of a viscous fluid past a fixed cylinder. // Mech. and Appl. Math., 1955, vol. VIII, part 2, p. 129-145.

53. Ильин A.M. Разностная схема для дифференциального уравнения с малым параметром при старшей производной. // Мат. заметки, 1969, т.6, N2, с.237-248.

54. Данаев Н.Т., Смагулов Ш.С. Об одной методике численного решения уравнений Навье-Стокса в переменных (\|/,ю). // Моделирование в механике. Сб. научных трудов. 1991, т.5(22), N4, с.38-47.

55. Кузнецов Ю.А., Местиашвили К.В. Об оптимизации вариационно-разностных методов для задач с пограничным слоем. // Сб. Разностные и вариационно-разностные методы. Новосибирск, 1977, N2, с.101-110.

56. Игнатьев В.Н., Задорин А.И. Разностная схема для дифференциального уравнения с малым параметром на неравномерной сетке. // Препринт 229, ВЦ СО АН СССР, Новосибирск, 1980.

57. Игнатьев В.Н., Задорин А.И. О плохой обусловленности при численном решении уравнений с малым параметром. // Препринт 84, ВЦ СО АН СССР, Новосибирск, 1981.

58. Кочергин В.П., Щербаков А.В. О разностных схемах второго порядка аппроксимации для эллиптического уравнения с малым параметром при старших производных. // Численные методы механики сплошных сред, 1974, т.5, N1, с.88-97.

59. Марчук Г.И., Шайдуров В.В. Повышение точности решений разностных схем.//М. Наука, 1979.

60. Barrett К.Е. The numerical solution of singular-perturbation boundary-value problems. // Quart. J. Mech. and Appl. Math., 1974, v.27, N1, p. 57-68.

61. Pearson C.E. On a differential equation of boundary layer type. // J. Math. Phys., 1968, N47, p. 134-154.

62. Miller I.I. Sufficient conditions for the convergence, uniformly in г, of a three point difference scheme for singular perturbation problem. // Lect. Notes Math., 1978, N679, p.85-91.

63. Алексеевский M.B. О разностной схеме для дифференциального уравнения с малым параметром при старшей производной. // В кн.: Разностные методы мат. физики. М., 1979, с. 36-60.

64. Бахвалов Н.С. К оптимизации методов решения краевых задач при наличии пограничного слоя. // ЖВМ и МФ, 1969, т.9, N4, с.841-859.

65. Дубовик К.Г., Никитин С.А., Полежаев В.И. Гидродинамические эффекты температурного и концентрационного расслоения. // В сб. Проблемы вязких течений. ИТПМ СО РАН, Новосибирск, 1981, с. 5564.

66. Полежаев В.П., Грязнов B.JI. Метод расчета граничных условий для уравнений Навье-Стокса в переменных 'вихрь, функция тока'. // Доклады Академии Наук СССР, 1974, т. 219, N2, стр. 301-304.

67. Кускова Т.В. Численное исследование двумерных течений вязкой несжимаемой жидкости. // В сб. Некоторые применения метода сеток в газовой динамике. Вып.З. 1971.

68. Полежаев В.П., Валъциферов Ю.В. Численное исследование нестационарной тепловой конвекции в цилиндрическом сосуде при боковом подводе тепла. //В сб. Некоторые применения метода сеток в газовой динамике. Вып.З. 1971.

69. Тарунин E.JI. Анализ аппроксимационных формул для вихря скорости на твердой границе. // Уч. зап. Пермск. пед. ин-та. Гидродинамика. 1976, N 152, вып. 9, с. 167-179.

70. Мызникова Б.И., Тарунин E.JI. О граничных условиях для вихря скорости в задачах динамики вязкой жидкости. //В кн.: Конвективные течения и гидродинамическая устойчивость. Свердловск: Изд. УНЦ АН СССР, 1979, с. 90-101.

71. Бердников B.C., Кирдяшкин А.Г. Структура свободно-конвективных течений в горизонтальном слое жидкости при различных граничных условиях. // В сб. Структура пристенного пограничного слоя. Новосибирск, 1978, с. 5-45.

72. Thom A. An investigation of fluid flow in two dimensions. // Aer. Res. C.R.aM. 1928, N1194.

73. Woods L. Note on the numerical solution of fourth order differential equations. // Aero Quart. 1954, N5.

74. Грязнов B.JI., Полежаев В.И. Исследование некоторых разностных схем и аппроксимаций граничных условий для численного решения уравнений тепловой конвекции. // Препринт Ин-та проблем механ. АН СССР, 1974, N40.

75. Герасимов Б.П. Один метод расчета задачи конвекции несжимаемой жидкости. // Препринт Ин-та прикл. матем. АН СССР, 1975, N13.

76. Schneider G.E., Zedan М. A modified strongly implicit procedure for the numerical solution of field problems. // Numerical Heat Transfer, 1981, vol. 4, pp. 1-19.

77. Яненко Н.Н. Метод дробных шагов решения многомерных задач математической физики. // Новосибирск. Наука. 1967.

78. Талонов В.А. Численное решение задачи о конвективном течении в замкнутой полости методом компактных разностей повышенного порядка точности. // Препринт N 273-94.

79. Vahl Davis G. De, Jones I.P. Natural convection in a square cavity: a comparison exercise. // Int. J. for Num. Meth. in Fluids, 1983, vol.3, pp.227248.

80. Вабищевич П.Н. Численные методы решения задач со свободной границей. // Изд-во М. Ун-та, 1987.

81. Пухначев В.В. Движение вязкой жидкости со свободными границами: Учебное пособие. // НГУ, Новосибирск, 1989.

82. Cheblakova Е.А., Gadiyak G.V. Convective flow in an upright enclosed square cavity: a comparison exercise. // International Journal of Computational Fluid Dynamics, 1998, vol. 10, N2, pp. 139-149.

83. Cheblakova E.A., Gadiyak G.V. Natural convection in an enclosed square cavity: numerical solution using method of integral identities. // Proceedings of the V Int. Conference on Simulation of Devices and Technologies, Obninsk, 1996, pp.135-141.

84. Гадияк Г.В., Чеблакова E.A. Моделирование роста кристаллов методом Чохральского с учетом образования дислокаций. // Математические модели и численные методы механики сплошных сред. Международная конференция. Тезисы докладов. Новосибирск, 1996.

85. Гадияк Г.В., Чеблакова Е.А. Моделирование роста кристаллов методом Чохральского. // II Сибирский конгресс по прикладной и индустриальной математике. Тезисы докладов. Новосибирск, 1996.

86. Cheblakova E.A., Gadiyak G.V. Convective flow in an upright enclosed square cavity: a comparison exercise. // (электронная публикация) http://rii.arxiv.Org/find/phvs3cs/l/Cheblakova/0/l/0/1999/3/0

87. Том А., Эйплт К. Числовые расчеты полей в технике и физике. // M.-JI. Энергия, 1964.

88. Роуч П. Вычислительная гидродинамика. // М. Мир. 1980.

89. Калис Х.Э. О постановке граничных условий для решения системы уравнений Навье-Стокса в переменных функции тока и вихря скорости. // Сб. Проблемы вязких течений. СО АН СССР, 1971.

90. Кускова Т.В., Чудов JI.A. О приближенных граничных условиях для вихря при расчете течений вязкой несжимаемой жидкости. // Вычислительные методы и программирование, изд-во МГУ, 1968, вып. XI.

91. Воеводин А.Ф. Об устойчивости разностных граничных условий для функции вихря на твердой стенке. // ЖВМ и МФ, 1998, т.38, N5, с. 855859.

92. Pearson С.Е. A computational method for viscous flow problems. // J. Fluid Mech., 1965, v.21, part 4, p.611-622.

93. Булеев Н.И., Тимухин Г.И. Течение вязкой жидкости на входном участке плоского канала. // ПМТФ, 1967, N3, с.126-130.

94. Jensen V.G. Viscous flow round a sphere at low Reynolds numbers (<40). // Proc. Roy. Soc. London, 1959, ser. A, v.249, p.346-366.

95. Годунов С.К., Рябенький B.C. Разностные схемы. // М. Наука. 1973.

96. Булеев H.H. Численный метод решения двумерных уравнений диффузии. //Мат. сб., 1960, т.51(93), N2, с.227-238.

97. Смирнов В.И. Курс высшей математики. Т.1,2. // М. Гос. Изд-во Физико-мат. л-ры, 1962.

98. Остапенко В.В. Метод теоретической оценки дисбалансов неконсервативных разностных схем на ударной волне. // Доклады Академии Наук СССР, 1987, т.295, N2, с. 292-297.

99. Андерсон Д., Таннехилл Д., Плетчер Р. Вычислительная гидромеханика и теплообмен. // М. Мир, т. 1,2, 1990.

100. Флетчер К. Вычислительные методы в динамике жидкостей. // М. Мир, т. 1,2, 1991.

101. Шлихтинг Г. Теория пограничного слоя. // М. Наука, 1974.

102. Седов Л.И. Механика сплошной среды. // М. Наука, 1973.

103. Тарунин Е.Л. Нестационарная конвекция жидкости в замкнутой полости. // Известия АН СССР, МЖГ, 1968, N6, с.83-88.

104. Гершуни Г.З., Жуховицкий Е.М., Тарунин Е.Л. Численное исследование конвективного движения в замкнутой полости. // Известия АН СССР, МЖГ, 1966, N5, с.56-62.

105. Гершуни Г.З., Жуховицкий Е.М., Тарунин Е.Л. Численное исследование конвекции жидкости, подогреваемой снизу. // Известия АН СССР, МЖГ, 1966, N6, с.93-99.

106. Рихтмайер Р., Мортон К. Разностные методы решения краевых задач. //М. Мир, 1972.

107. Самарский А.А, Гулин A.B. Устойчивость разностных схем. // М. Наука, 1973.

108. Chorin A.J. A numerical method for solving incompressible viscous flow problems. // J. Comput. Phys., 1967, vol.2, pp. 12-26.

109. Антонцев С.Н., Кажихов А.В., Монахов А.В. Краевые задачи механики неоднородных жидкостей. // Н. Наука, 1983.

110. Gadiyak G.V., Korobitsina J.L., Kramarenko V.I. Numerical simulation of the thermal oxidation of silicon. // COMPEL, 1992, vol.11, N4, pp.419-431.