авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 ||

«РОССИЙСКАЯ АКАДЕМИЯ НАУК СИБИРСКОЕ ОТДЕЛЕНИЕ ИНСТИТУТ ГИДРОДИНАМИКИ им. М.А. ЛАВРЕНТЬЕВА На правах рукописи ...»

-- [ Страница 2 ] --

Попытки получить численное решение по такой-же схеме для 4 варианта оказались неудачными: на сетках 30 30 и 30 60 процесс расчета закан чивался аварийным остановом. Поэтому для 4–6 вариантов задача решалась следующим образом: производилось 6 · 104 итераций по времени на сетке 60 60, затем 104 итераций на сетке 60 120 и окончательно 5 · итераций на сетке 120 120. На последнем шаге по времени значение начального приближения при решении задачи для для вариантов 4, 5, составило на сетке 120120 соответственно 2.5·104, 1.9·109 и 2.4·107.

Число итераций L при решении задачи для составило для 4, 5, 6 вари антов 20, 2 и 10 соответственно, значение по окончании решения задачи 1.1 · 107, 3.7 · 1010 и 9.1 · 109 соответственно, а величина для 5.9, 2.3 · 104 и 2.0 · 101 соответственно.

по окончании расчетов На рис. 6, а–е показаны изолинии функции тока для 1–6 вариантов соот ветственно, а на рис. 7, а–е изотермы для 1–6 вариантов соответственно.

На рис. 8, а, б, в показаны изолинии модифицированной азимутальной ком поненты скорости W соответственно для 4, 5, 6 вариантов.

В табл. 6 приведены характеристики течений для 1–6 вариантов расче та (здесь vrz модуль проекции вектора скорости расплава на плоскость а б в г д е Рис. а б в г д е Рис. а б в Рис. rOz;

rmin, zmin и rmax, zmax координаты r, z точек минимума и макси мума соответственно;

единицей измерения скорости служит v1 = 10 см/с, а температуры градус K;

началом отсчета температуры является точка плавления).

Из табл. 6 видно, что в 1, 4 и 5 вариантах расчета минимальная темпера тура расплава получилась меньше нуля. Но сама зона отрицательных темпе ратур очень мала. В 1 и 5 вариантах она примыкает к границам m, f, а в 4 варианте к границам 0, c. Анализ выходного массива T показывает, что в плоскости rOz эта зона помещается в квадрат со стороной 3·102 см, 5.4 · 103 см и 4.5 · 103 см соответственно для 1, 4 и 5 вариантов расчета.

Существование такой зоны может быть объяснено тем, что форма границы плавающей зоны взята близкой к приведенной в [56], а условия задачи в настоящей работе существенно отличны от принятых в [56] (более концен трированное тепловыделение и как следствие бльшая пондеромоторная о сила, которая, кроме того, не сносится из уравнения импульса в граничное условие для вихря, что приводит к значительно бльшим скоростям движе о ния расплава;

другая форма функции пондеромоторной силы на m ;

неучет теплового излучения индуктора и др.).

На рис. 9, а показаны распределение касательной скорости на свободной границе (s безразмерная длина дуги, отсчитываемая от нижней точки свободной границы, положительным для скорости является направление в Таблица Функция Вариант rmin zmin min rmax zmax max 1 2.94 0.043 0.233 1.82 0.175 1. 2 2.77 0.030 0.177 1.60 0.255 0. 3 2.65 0.034 0.075 0.770 0.076 0. 4 2.90 0.033 0.321 1.82 0.175 1. 5 2.95 0.015 0.193 1.74 0.189 1. 6 2.70 0.018 0.070 1.13 0.076 0. vrz 1 0 1.56 0.396 3. 2 0 3.40 0.002 3. 3 0 0.199 0.263 0. 4 0 0 0.058 9. 5 0 1.56 0.396 3. 6 0 0.343 0.134 0. T 1 0.885 0.891 0.239 1.03 0.558 39. 2 0 1.29 0.446 85. 3 0 2.57 0.395 122. 4 0.001 1.138 9.38 0.880 0.764 40. 5 0.887 0.898 0.072 1.20 0.476 42. 6 0 2.02 0.384 126. W 4, 5, 6 0.888 0.902 0.256 3.40 0 0. сторону возрастания s) и на различных расстояниях от нее для 5 варианта расчетов. Кривая 1 распределение скорости на самой свободной границе, на расстоянии Em, 3 на расстоянии 2Em и 4 на расстоянии 3Em от нее. На рис. 9, б показано распределение скорости вдоль внутренней нормали n, за единицу длины которой принята величина Em, при s = 1.5.

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

Кривая 1 распределение температуры на самой свободной границе, на расстоянии 3Em, 3 на расстоянии 6Em и 4 на расстоянии 9Em от нее. Рис. 10, б аналог рис. 9, б с заменой функции касательной скорости а б Рис. а б Рис. Рис. на функцию температуры.

Наконец, на рис. 11 приведены зависимости минимума кривые 1, 3, 5 и максимума кривые 2, 4, 6 функции тока от размерного времени соответ ственно для 4, 5, 6 вариантов расчета.

Проанализируем полученные результаты.

Движение, осуществляющееся в 1 варианте расчетов, можно назвать электроконвекцией, так как оно вызвано только силой электрической приро ды (влияние протягивания образца ввиду малой его скорости несуществен но). Максимум скорости в этом варианте достигается вблизи свободной гра ницы и составляет около 38 см/с. Характерная скорость внутри области течения (то есть за границами пограничных слоев для скорости) составляет 10–20 см/с. Этого оказывается достаточно для формирования температур ных пограничных слоев (сгущения изотерм вблизи границ расплава см.

рис. 7, а). Скорость на свободной границе (аналог кривой 1 на рис. 9, а) не имеет отрицательного всплеска вблизи точки s = 0 и положительного вблизи точки s = 3 (эти всплески у кривой 1 на рис. 9, а обусловлены термокапиллярным эффектом).

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

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

Движение, осуществляющееся в 4 и 5 вариантах расчета, можно назвать комбинированной конвекцией, так как оно обусловлено действием всех фак торов в совокупности. Но между ними есть качественное различие, вызван ное разным видом функции пондеромоторной силы. В 4 варианте на оси симметрии возникает струя с максимумом скорости 95 см/с. Радиус участ ка, где скорость превышает 30 см/с, составляет 1 мм, а высота этого участка около 1 см (как показали дополнительные расчеты, при отсутствии вра щения и прочих равных условиях этого участка в окрестности оси симмет рии нет). Течение в 4 варианте расчетов носит колебательный характер, о чем говорит высокое значение величины = 5.9 и вид кривых 1, 2 на рис. 11, изображающих зависимость минимума и максимума функции тока от времени для этого варианта. Период колебаний обеих величин составляет 0.2025 с, временной интервал вывода их на график 0.0075 с. (Под периодом здесь понимается расстояние между соседними максимумами, движение при этом может и не быть периодическим.) Размах колебаний составляет 0. для минимума и 0.0244 для максимума функции тока, а среднее по периоду колебаний значение равно -0.3092 для минимума и 1.1315 для максимума функции тока. В 5 варианте расчета максимум скорости на оси симметрии составляет 2.6 см/с, то есть движение вблизи оси симметрии очень медлен ное. Глобальный максимум скорости достигается вблизи свободной границы и составляет 36 см/с. Невысокое значение величины = 2.3 · 104 по окон чании расчетов в 5 варианте и вид кривых 3, 4 на рис. 11, изображающих зависимость минимума и максимума функции тока от времени для этого ва рианта, говорят о том, что решение в этом случае вышло на режим, близкий к стационарному. Начиная с момента t = 6.03 с до момента t = 75 с значе ние минимума функции тока монотонно возрастает от -0.20275 до -0.19273, а значение максимума функции тока при 4.92 c t 75 c монотонно убывает от 1.03515 до 1.00265.

Движение, осуществляющееся в 6 варианте расчетов, назовем термогра витационной конвекцией с вращением. Максимум скорости в меридиональ ной плоскости достигается для этого варианта вблизи оси симметрии и со ставляет 2.2 см/с. Ввиду медленности движения тепловых пограничных сло ев вблизи границ области, занятой расплавом, не образуется. Течение явля ется нестационарным и носит колебательный характер, о чем говорит до вольно высокое значение величины = 0.20 и вид кривых 5, 6 на рис. 11, изображающих зависимость минимума и максимума функции тока от вре мени для этого варианта. Минимум функции тока испытывает колебания, не видные на рисунке, ввиду малости их размаха, равного 0.002. Размах ко лебаний максимума функции тока составляет 0.018. Период в обоих случаях равен 6 с. Средние по периоду колебаний значения составляют -0.070695 для минимума и 0.040965 для максимума функции тока.

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

В зонах больших скоростей она близка к константе (для каждой зоны своя константа). Это происходит потому, что функция W = const с точностью до малых членов 2W r r + rH 2 Re x x y y (ввиду гармоничности функции r(x, y) и формул (2.4.3)) удовлетворяет уравнению (2.4.14). Это явление аналогично явлению постоянства вихря в плоской задаче в зоне больших скоростей, детально разобранному в работе [45], и явлению постоянства температуры внутри области течения в 1, и 5 вариантах, менее выраженному ввиду относительно небольшого числа Пекле. Такое-же явление должно иметь место и для модифицированного вихря в 1 и 2 вариантах расчета, в которых отсутствуют сила плавучести и вращение.

Дадим теоретическую оценку максимума скорости на свободной границе для 1 варианта, когда движение стационарно и вызвано только пондеромо торной силой. Уравнение (2.4.4) при x = 1 можно приближенно записать в виде:

r 2 (Hfnxy ) V = Eum.

Re x y y Так как при отсутствии термокапиллярной силы /y = 0 на m, первый член в левой части этого равенства исчезает. Кроме того, согласно (2.4.1), (2.4.2) 1 v = rH x (H, r в направлении x практически не меняются), а 1. Поэтому 1 3v 1 (Hfnxy ) = Eum 2.

ReH 3 x3 H y Это эквивалентно записи 1 3v fn = Eum, Re n3 s где s, n, и fn (s, n) определены в разделе 2.3. Характерное значение fn /s составляет 2/3 (значение для треугольного вида функции fn (s, 0)). А так как модуль скорости v убывает в пограничном слое примерно в 2 раза, а эффективное расстояние, на котором действует пондеромоторная сила, равно Em /2, то 3 v/n3 0.5v0 (2/Em )3, где v0 максимум модуля v.

Отсюда получаем искомую оценку:

1 v0 = ReEum Em = 4.4, то есть 44 см/с. Это всего в 1.16 раза больше полученного в расчетах для 1-го варианта максимума скорости на свободной границе. Так как при комбини рованной конвекции электроконвекция доминирует, то эта оценка годится и для 4-го и 5-го вариантов.

Рассчитаем характерное отношение вязкого V c и конвективного V k членов в уравнении для вихря в скин-слое:

y110 y Vc= V c(y)dy, Kc = Kc(y)dy, y10 y где 1 r 2 (U ) (V ) V c(y) = dx, Kc(y) = + dx, Re x2 x y 1(y) 1(y) (y) = 2Em /H(1, y), y10, y110 значения сеточной функции ym при m = 10 и m = 110 (они введены для того, чтобы отсечь большие значения функций V c(y), Kc(y) при y 0 и y yM, где M = 120). Для 5-го варианта расчетов найдем:

V c/V k = 1.264, то есть вязкие и конвективные члены имеют в скин-слое одинаковый порядок.

Получим аналогичный результат аналитически. Домножим уравнение (2.4.4) на r и перейдем к переменным касательная s нормаль n к свободной границе. Характерное значение конвективного члена равно 1 (rv) 1 v v Kc (V ) = rv = v.

H 2 y s s r n s n Так как характерное значение v равно 2 (т. е. 20 см/с в размерных перемен ных) и такое-же значение имеет характерный перепад скорости в скин-слое, эффективная толщина которого равна Em /2, то 2 1 Kc 2· · =.

0.5Em 0.5sm Em sm Характерное значение вязкого члена равно r 2 r 1 3v 1 3v 1 2 Vc = =.

Re n2 Re r n3 Re n3 Re (Em /2)3 ReEm Тогда их отношение будет Vc sm = = 1.772, Kc ReEm то есть они имеют одинаковый порядок.

Толщину возникающих динамического v и теплового T пограничных слоев в 1, 4 и 5 вариантах можно теоретически оценить по формулам:

m Rc m Rc v = = 1.74Em, T = = 16.1Em.

v1 m cpm v Сравнение теоретических v, T и фактических v, T (см. рис. 9, б, 10, б) толщин пограничных слоев показывает, что v v, T 2T, то есть согласие между ними вполне приемлемое.

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

0 rc = 2 = 0.00077, v1 l где 0 = 3.4 · 103 см2 /с коэффициент кинематической вязкости, исполь зуемый в [56]. 1 и 2 варианты в [56] соответствуют комбинированной кон векции с разбиением области течения на 7760 и 3006 треугольных элементов соответственно. 3 вариант соответствует термогравитационной конвекции с вращением. Приводимые ниже данные являются универсальными для всех трех вариантов, рассчитанных в [56]. Средние значения минимума и макси мума функции min = 0.058, max = 0.037, максимальная скорость (ее можно оценить из рис. 7 [56]) vmax = 0.1... 0.2, то есть 1... 2 см/с. Перегрев расплава 36.8...48.2 K. Период колебаний мини мума и максимума функции примерно 4 с и 7 с соответственно. Размах колебаний 0.0024 для минимума и 0.012... 0.016 для максимума функции.

Таким образом, результаты расчета, сделанного в настоящей работе и в [56] для комбинированной конвекции отличаются принципиально (за исклю чением максимума температуры, для которого имеет место случайное сов падение), а для случая термогравитационной конвекции с вращением до вольно близки (за исключением максимума температуры). Последний факт можно объяснить следующим образом. Хотя перегрев расплава в настоящей работе больше полученного в [56] в 3 раза, но используемый коэффициент T в 2.3 раза меньше, так что эффективное число Грасгофа больше толь ко на 30%. Эффективное число Пекле в обеих работах примерно одинако во, коэффициент теплопроводности (а следовательно и эффективное число Био) в настоящей работе меньше на 22%, вязкость на 5–30%, используется другая функция источников джоулева тепла, не учитывается тепловое излу чение индуктора и используются несколько другие граничные условия для вихря на твердых стенках и свободной границе. Таким образом, различия в постановке задачи хотя и имеют место, но не очень велики. Это и обес печивает близость результатов расчета термогравитационной конвекции с вращением в обеих работах (в том числе и близость изолиний функции тока см. рис. 6 е настоящей работы и рис. 7 d работы [56]). Разницу в макси мальной температуре (перегреве) расплава для рассматриваемого варианта можно объяснить только разной долей выделяющегося в плавающей зоне джоулева тепла (60% в настоящей работе и в несколько раз меньшей долей в работе [56]).

Были также произведены расчеты комбинированной конвекции (7 ва риант) и термогравитационной конвекции с вращением (8 вариант) при H0 = 203.7 Гс для второго варианта задания функции f (x), что соот ветствует мощности выделяющегося в плавающей зоне джоулева тепла, со ставляющей 25 % от мощности индуктора, при характеристиках материала, принятых в [56]:

cpm = 107 см2 /(К · с2 ), m = 6.7 · 106 г · см/(K · с3 ), T = 3.54 · 104 г/(см3 · K), (T ) m = 3.4 · 103 см2 /с.

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

7 вариант:

vrzmax = 2.164 (т. е. 21.64см/с), Tmax = 13.07K, min = 0.098, max = 0.611, V c/Kc = 2.166, максимальная скорость вне скин-слоя равна 1.3 (т. е. 13 см/с).

8 вариант:

vrzmax = 0.118 (т. е. 1.18см/с), Tmax = 34.86K, min = 0.067, max = 0.027.

Движение расплава в обоих вариантах получилось близким к стационарно му.

Так как максимум температуры в 8 варианте близок к полученному в [56] для аналогичного варианта, мощность выделяющегося в плавающей зоне джоулева тепла также должна быть близкой к полученной в [56] (она одна для всех рассматриваемых там вариантов). А поэтому можно сравнивать 7 вариант расчета с 1 вариантом [56] на предмет того, к чему приводит от каз от игнорирования конвективных членов в скин-слое при прочих близких условиях. Он приводит к возрастанию скорости вне скин-слоя на порядок, к уменьшению перегрева расплава в 3 раза, к увеличению модуля макси мума и минимума функции тока в 1.7 раза и к качественно другой картине течения (разному виду изолиний функции тока и температуры, в частности, к появлению температурных пограничных слоев).

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

ГЛАВА УСЛОВИЯ МОНОТОННОСТИ ФАКТОРИЗОВАННОЙ РАЗНОСТНОЙ СХЕМЫ ДЛЯ ПАРАБОЛИЧЕСКОГО УРАВНЕНИЯ В этой главе рассматривается проблема выбора временного шага и весового коэффициента схемы при численном интегрировании нестационарного урав нения с двумя пространственными переменными, описывающего диффузи онный и конвективный перенос субстанции в сплошной среде. Разностная схема записывается в факторизованном виде. Ее коэффициенты удовлетво ряют некоторым условиям, обеспечивающим монотонность при достаточ но малом шаге по времени. Формулируются ограничения на этот шаг, до статочные для монотонности схемы. (Эти результаты в некоторой степени аналогичны результатам Хартена [48, p. 6–7], которые выражают условия монотонности абстрактной трехточечной разностной схемы в случае одной пространственной переменной.) Рассмотрена также задача для уравнения теплопроводности с малым числом разбиений, для которой получены необ ходимые и достаточные условия монотонности схемы. Оказалось, что в этом случае достаточные условия монотонности весьма близки к необходимым и достаточным при значениях весового коэффициента, меньших 1/2, и далеки от них при 1. Предложено скорректировать (без доказатель ства) достаточные условия монотонности, чтобы снять это несоответствие.

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

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

3.1. Описание разностной схемы Рассмотрим уравнение относительно F :

(/t + L1 + L2 )F = G, (3.1.1) где 1 F L1 F = (p + U )F µ ;

rH 2 x x 1 F L2 F = (q + V )F µ ;

(3.1.2) rH 2 y y t время;

x, y независимые пространственные переменные, такие, что функция r(x, y)+iz(x, y) конформно отображает некоторый прямоугольник на заданную ограниченную односвязную область;

r, z, цилиндриче ские координаты, в которых r есть полярный радиус, полярный угол;

(r/x)2 + (r/y) H= коэффициент Ламэ;

U, V модифициро ванные компоненты скорости, удовлетворяющие уравнению неразрывности U/x + V /y = 0, и связанные с физическими компонентами скорости u, v в направлениях x, y соотношениями U = rHu, V = rHv;

G, p, q, µ заданные функции. Как было показано в главе 2, к (3.1.1) сводятся урав нения для функций, W, T, где = /r, W = wr, вихрь, w азимутальная компонента скорости, T температура в плавающей зоне.

Пусть в прямоугольнике задана неравномерная сетка 0 = x0 x1... xN 1 xN = 1, 0 = y0 y1... yM 1 yM = Y.

Обозначим xn + xn+ xn+1/2 =, hxn+1/2 = xn+1 xn, n = 0, N 1;

ym + ym+ ym+1/2 =, hym+1/2 = ym+1 ym, m = 0, M 1;

xn+1 xn1 ym+1 ym hxn =, n = 1, N 1;

hym =, m = 1, M 1.

2 Будем предполагать, что имеют место соотношения hxn+1/2 hxn1/2 = O(h xn+1/2 ), n = 1, N 1, hym+1/2 hym1/2 = O(h ym+1/2 ), m = 1, M 1.

Введем разностные аналоги 1, 2 операторов L1, L2 следующим образом (см. раздел 2.5):

1 Fnm = ax Fn+1m + cx Fnm bx Fn1m, nm nm nm 2 Fnm = ay Fnm+1 + cy Fnm by Fnm1, (3.1.2) nm nm nm где Un+1/2m + pn+1/2m µn+1/2m n+1/2m coth n+1/2m ax = +, nm rnm Hnm 2hxn hxn+1/2 hxn Un+1/2m + pn+1/2m Un1/2m + pn1/2m cx = + nm rnm Hnm 2hxn 2hxn µn+1/2m n+1/2m coth n+1/2m µn1/2m n1/2m coth n1/2m + +, hxn+1/2 hxn hxn1/2 hxn Un1/2m + pn1/2m µn1/2m n1/2m coth n1/2m bx = +, nm rnm Hnm 2hxn hxn1/2 hxn Vnm+1/2 + qnm+1/2 µnm+1/2 nm+1/2 coth nm+1/ ay = +, nm rnm Hnm 2hym hym+1/2 hym (3.1.3) Vnm+1/2 + qnm+1/2 Vnm1/2 + qnm1/ cy = + nm rnm Hnm 2hym 2hym µnm+1/2 nm+1/2 coth nm+1/2 µnm1/2 nm1/2 coth nm1/ + +, hym+1/2 hym hym1/2 hym Vnm1/2 + qnm1/2 µnm1/2 nm1/2 coth nm1/ by = +, nm rnm Hnm 2hym hym1/2 hym hxn+1/2 (Un+1/2m + pn+1/2m ) n+1/2m =, 2µn+1/2m hym+1/2 (Vnm+1/2 + qnm+1/2 ) nm+1/2 =.

2µnm+1/ k Здесь для произвольной функции f (x, y, t) принято обозначение fnm = f (xn, ym, k), где шаг по оси t, индекс k пока для простоты опущен.

Для численного решения уравнения (3.1.1) используем разностную схему k+1 k Fnm Fnm + (1 + 2 )(Fnm + (1 )Fnm ) = Gk, k+1 k (E + 1 2 ) nm n = 1, N 1;

m = 1, M 1, (3.1.4) где E тождественный оператор, [0, 1] весовой коэффициент.

На каждой из 4-х сторон прямоугольника могут быть поставлены условия 1-го, 2-го или 3-го рода, но так как, например, при x = 0, последние два могут быть сведены к 1-му добавлением фиктивного узла x1 = x1 (см.

[34, с.419]), то будем для простоты рассматривать только постановку условий 1-го рода:

F0m = g1m, FN m = g2m, Fn0 = g3n, FnM = g4n. (3.1.5) Зададим также начальное условие Fnm = Fnm. (3.1.6) Описанная схема является консервативной. Если функции U, V и µ од ного порядка, то она превращается в схему с центральными разностями для конвективных членов, которая имеет 2-й порядок аппроксимации по про странству. Если U, V µ, то данная схема близка к схеме для уравнений Эйлера с разностями против потока, которая имеет 1-й порядок аппрок симации по пространству. По времени описанная схема имеет 1-й порядок аппроксимации.

3.2. Достаточные условия монотонности в общем случае Согласно работе [17, с. 348] в случае одной пространственной переменной схе ма называется монотонной, если она сохраняет монотонность решения при переходе с k-го на k+1-й временной слой. Там же сформулирован признак монотонности: схема N k+1 k Fn = nm Fm m= монотонна тогда и только тогда, когда все коэффициенты nm 0. Так как этот признак необходимый и достаточный, его можно принять за опре деление монотонности. А последнее легко распрастраняется на случай 2-х пространственных переменных: схему N M k+ [nmij Fij + nmij Gk ] k Fnm = ij i=0 j= назовем монотонной, если все nmij 0.

Предложение 3.2.1. Пусть операторы 1, 2 определены равенства ми (3.1.2), причем коэффициенты ax, ay, bx, by nm nm nm nm неотрицательны, а коэффициенты cx, cy nm nm положительны (например, эти коэффициенты вычисляются по формулам (3.1.3)). Тогда разностная схема (3.1.4) с гра ничными условиями (3.1.5) и начальным условием (3.1.6) монотонна при max, где max = min{ 1, 2 }, 11, при = 0, 1 = 12, при 0 0.5,, при 0.5 1, 11 =, X +Y 2 (X Y )2 + (1 2)(X + Y ) (1 )(X + Y ) 12 =, 2 2 XY (1 ) (1 ) 13 = min,2, (3.2.1) 2X Y X = max{cx : 1 n N 1;

1 m M 1}, nm Y = max{cy : 0 n N ;

1 m M 1}, nm 2 = min{ 21, 22 },, если A 0, 21 = 1/A, если A 0,, если B 0, 22 = 1/B, если B 0, A = max{(ax + bx cx ) : 1 n N 1, 1 m M 1}, nm nm nm B = max{(ay + by cy ) : 1 n N 1, 1 m M 1}.

nm nm nm Доказательство. Запишем разностную схему (3.1.4) в факторизованном виде:

(E + 1 )(E + 2 )Fnm = (E (1 ) (1 + 2 ) + 2 2 1 2 )Fnm + Gk, k+1 k nm n = 1, N 1, m = 1, M 1. (3.2.2) Вычисления по схеме (3.2.2) можно разбить на 3 этапа:

1) вычисление правой части (3.2.2), 2) обращение оператора E + 1, 3) обращение оператора E + 2.

Для монотонности схемы в целом достаточно, чтобы она была монотонной на каждом из 3-х этапов. Рассмотрим этап 1.

nm = (E (1 ) (1 + 2 ) + 2 2 1 2 )Fnm + Gk = k k nm = (1 (1 ) (cx + cy ) + 2 2 cx cy )Fnm + k nm nm nm nm +(((1 ) 2 2 cy x k n+1m )anm )Fn+1m + (3.2.3) +(((1 ) 2 2 cy x k n1m )bnm )Fn1m + +(((1 ) 2 2 cx )ay )Fnm+1 + k nm nm +(((1 ) 2 2 cx )by )Fnm1 +... + Gk.

k nm nm nm Здесь не выписаны члены с заведомо неотрицательными коэффициента ми при Fn±1m±1. Пусть 0. Тогда нужно решить систему неравенств относительно :

1 (1 ) (x + y) + 2 2 xy 0, x (0, X], y (0, Y ], (3.2.4) min (1 )/( 2 X), (1 )/( 2 Y ), где X, Y определены в (3.2.1). Обозначим f (x, y) = 1 (1 ) (x + y) + 2 2 xy. Имеем: f /x = (1 ) + 2 2 y 0 при y y0 = (1 )/( 2 ) и f /y = (1 ) + 2 2 x 0 при x y0. Из 2 го неравенства (3.2.4) имеем: X y0 ;

Y y0. Следовательно, минимум функции f (x, y) достигается при x = X;

y = Y. Теперь первое неравенство (3.2.4) можно заменить на 1 (1 ) (X + Y ) + 2 2 XY 0. (3.2.5) Его решение есть 12 (см. формулы (3.2.1)). Пусть для определен ( 12 13 )2 2 XY = (1 )(X ности X Y. Составим разность 2 (X Y )2 + (1 2)(X + Y )2. Предположим, что таково, что Y) подкоренное выражение неотрицательно. Тогда мы имеем разность двух неотрицательных чисел, а последняя имеет тот – же знак, что и разность их квадратов, которая равна (1 2)(2X 2 + 2Y 2 ). Таким образом, при 1/2 12 13, при = 1/2 12 = 13, а при 1/2 12 на участке существования вещественного 12. Если же таково, что комплексное, то неравенство (3.2.5) выполняется при любом и, следова тельно, ограничение на дает только функция 13. При этом 1/2.

Определим 1 при = 0 как предельное значение 12 при 0. В результате получим формулу для 11 в (3.2.1) (это несколько более грубая оценка чем та, которая непосредственно следует из (3.2.3)).

Рассмотрим этап 2. Мы имеем серию одномерных задач:

k+1/2 k+1/ ax Fn+1m + (1 + cx )Fnm bx Fn1m = nm, n = 1, N 1, k+1/2 k nm nm nm k+1/2 k+1/ F0m = (E + 2 )g1m, FN m = (E + 2 )g2m, m = 1, M 1.

Известно (см. [11, с. 270]), что для монотонности схемы на этом этапе до статочно неотрицательности коэффициентов ax, bx, которая имеет ме nm nm сто в силу условий предложения, и свойства диагонального преобладания:

1 + cx (ax + bx ), которое можно записать в виде 21, где nm nm nm 21 вычисляется по формулам (3.2.1).

На этапе 3 опять имеем серию одномерных задач ay Fnm+1 + (1 + cy )Fnm by Fnm1 = Fnm, m = 1, M 1, k+1 k+1 k+1 k+1/ nm nm nm k+1 k+ Fn0 = g3n ;

FnM = g4n, n = 1, N 1, для которой аналогично этапу 2 получаем 22. Предложение доказано.

3.3. Необходимые и достаточные условия монотонности для задачи с малым числом разбиений Доказанное в предыдущем разделе предложение дает ограничения на шаг, выполнения которых достаточно для монотонности схемы. Возникает вопрос: насколько эти ограничения близки к необходимым и, в частности, является ли схема безусловно немонотонной при = 1 (max = 0 при = 1).

Для ответа на эти вопросы рассмотрим разностную схему для уравнения теплопроводности F/t ( 2 F/x2 + 2 F/y 2 ) = с числами разбиений N = M = 4 и равномерными шагами по осям x, y (здесь = const коэффициент температуропроводности). Положим на границе области F = 0. Исключив граничные точки, как показано в [23, с.

363], получим разностную задачу с 9-ю узлами на каждом временном слое:

k+1 k (E + 1 )(E + 2 )Fnm = Fnm, n, m = 1, 2, 3. (3.3.1) Введем обозначения:

= /h2, l = h2 /h2.

x x y Оператор E + 1 действует по столбцам и задается матрицей 1 + 2.

1 + 0 1 + Оператор E + 2 действует по строкам и задается матрицей 1 + 2l l l l.

1 + 2l 0 l 1 + 2l Оператор действует следующим образом:

F11 = AF11 + BF21 + CF12 + DF22, F21 = AF21 + B(F31 + F11 ) + CF22 + D(F12 + F32 ), F31 = AF31 + BF21 + CF32 + DF22, F12 = AF12 + BF22 + C(F11 + F13 ) + D(F21 + F23 ), F22 = AF22 + B(F32 + F12 ) + C(F21 + F23 ) + D(F11 + F33 + F13 + F31 ), F32 = AF32 + BF22 + C(F31 + F33 ) + D(F21 + F23 ), F13 = AF13 + BF23 + CF12 + DF22, F23 = AF23 + B(F33 + F13 ) + CF22 + D(F12 + F32 ), F33 = AF33 + BF23 + CF32 + DF22, где A = 1 2(1 )(1 + l) + 4 2 2 l, B = (1 ) 2 2 2 l, C = (1 )l 2 2 2 l, D = 2 2 l.

Оператор (E + 1 )1 действует по столбцам и задается матрицей Ex Fx Gx 1 Fx Hx Fx, x Gx Fx Ex где Ex = (1 + 2)2 2 2, Hx = (1 + 2)2, Fx = (1 + 2), Gx = 2 2, x = (1 + 4 + 2 2 2 )(1 + 2).

Оператор (E + 2 )1 действует по строкам и задается матрицей Ey Fy Gy 1 Fy Hy Fy, y Gy Fy Ey где Ey = (1 + 2l)2 2 2 l2, Hy = (1 + 2l)2, Fy = (1 + 2l)l, Gy = 2 2 l2 ;

y = (1 + 4l + 2 2 2 l2 )(1 + 2l).

k+ Разрешая (3.3.1) относительно Fnm, получим 3 k+1 1 1 k k Fnm = (E + 2 ) (E + 1 ) Fnm = nmij Fij.

x y i=1 j= Среди 81 коэффициентов nmij в силу симметрии задачи только 25 различ ных. Это следующие коэффициенты:

1111 = Ey (Ex A + Fx B) + Fy (Ex C + Fx D), 2111 = Ey ((Ex + Gx )B + Fx A) + Fy ((Ex + Gx )D + Fx C), 1211 = (Ey + Gy )(Ex C + Fx D) + Fy (Ex A + Fx B), 2211 = (Ey + Gy )((Ex + Gx )D + Fx C) + Fy ((Ex + Gx )B + Fx A), 3111 = Ey (Fx B + Gx A) + Fy (Fx D + Gx C), 3211 = (Ey + Gy )(Fx D + Gx C) + Fy (Fx B + Gx A), 1311 = Fy (Ex C + Fx D) + Gy (Ex A + Fx B), 2311 = Fy ((Ex + Gx )D + Fx C) + Gy ((Ex + Gx )B + Fx A), 3311 = Fy (Fx D + Gx C) + Gy (Fx B + Gx A), 1121 = Ey (Fx A + Hx B) + Fy (Fx C + Hx D), 2121 = Ey (2Fx B + Hx A) + Fy (2Fx D + Hx C), 1221 = (Ey + Gy )(Fx C + Hx D) + Fy (Fx A + Hx B), 2221 = (Ey + Gy )(2Fx D + Hx C) + Fy (2Fx B + Hx A), 1321 = Fy (Fx C + Hx D) + Gy (Fx A + Hx B), 2321 = Fy (2Fx D + Hx C) + Gy (2Fx B + Hx A), 1112 = Fy (Ex A + Fx B) + Hy (Ex C + Fx D), 2112 = Fy ((Ex + Gx )B + Fx A) + Hy ((Ex + Gx )D + Fx C), 1212 = 2Fy (Ex C + Fx D) + Hy (Ex A + Fx B), 2212 = 2Fy ((Ex + Gx )D + +Fx C) + Hy ((Ex + Gx )B + Fx A), 3112 = Fy (Fx B + Gx A) + Hy (Fx D + Gx C), 3212 = 2Fy (Fx D + Gx C) + Hy (Fx B + Gx A), 1122 = Fy (Fx A + Hx B) + Hy (Fx C + Hx D), 2122 = Fy (2Fx B + Hx A) + Hy (2Fx D + Hx C), 1222 = 2Fy (Fx C + Hx D) + Hy (Fx A + Hx B), 2222 = 2Fy (2Fx D + Hx C) + Hy (2Fx B + Hx A).

Рис. Анализ показывает, что при 0 и любых (в том числе = 1) все коэффициенты nmij 0 в некоторой окрестности точки = 0. Следова тельно, можно определить такой максимально допустимый шаг max, что все коэффициенты nmij будут неотрицательны, если max. Умножив в формулах (3.2.1) левые и правые части на /h2, получим достаточные x условия монотонности относительно : max (, l), где 2+2l, при = 0, (1)(1+l) 2 (1l)2 +(12)(1+l) max (, l) = (3.3.2), при 0 1/2, 4 2 l min{ 1, 1l }, при 1/2 1.

2 2 2 На рис. 12 представлены зависимости max и max от при различных значениях l. Так как max (, 1/l) = lmax (, l) и аналогичным свойством обладает величина max, то рассматривался только случай l 1. А так как при 1/2 величины max, max уже не зависят от l при l 1, то принята различная нумерация кривых при 0.5 и 0.5. Кривые 1, 3, 5 показывают зависимость max, кривые 2, 4, 6 зависимость max от при l = 1, 0.25 и 0.01 соответственно, 0 1/2. Кривые 7, показывают зависимость от величин max и max соответственно при 1/2 1. Из рис. 12 видно, что при 1/2 max весьма близко к max. При 1/2, 1 кривые max и max далеки (max не стремится к нулю, в отличие от max ). Цифрой 9 на рис. 12 обозначена кривая 11 max = min,, 1.

2 2l Видно, что эта кривая уже близка к кривой 7. Поэтому можно выдвинуть гипотезу, что и в общем случае в формулах (3.2.1) можно заменить функцию 13 на 1 13 = min,, X Y при этом полученный критерий останется достаточным.

3.4 Пример расчета Рассмотрим задачу 2F 2F F = 0, 0 x 1, 0 y 1, t 0, x2 y t F = 0.5, x = 0;

F = 0.5, x = 1;

F = f (x, t), y = 0;

F = f (x, t), y = 1;

(3.4.1) F = (x 0.5) 0.5, t = 0, где 2 cos(k/2) k2 2 t f (x, t) = x 0.5 + e sin(kx), k k= 0, при x 0;

(x) = 1, при x 0.

Она имеет аналитическое решение F (x, y, t) f (x, t).

Рис. Введем в квадрате 0 x 1, 0 y 1 равномерную сетку: xn = n/N, n = 0, N ;

ym = m/M, m = 0, M. Имеем: 1 Fnm = (Fn+1m + 2Fnm Fn1m )/h2, 2 Fnm = (Fnm+1 +2Fnm Fnm1 )/h2, где hx = 1/N, hy = 1/M.

x y Используем для численного решения задачи (3.4.1) разностную схему (3.1.4).

Положим N = M = 25 (тогда l = 1). Сначала произведем расчет при = 0.1, = 3max (0.1, 1)h2 (см. формулу (3.3.2)). На рис. 13 приведена x зависимость от x полученной в расчете функции F при y = 0.48, t = (непомеченная кривая) и t = 2 (кривая, помеченная кружками). При увеличении t осцилляции в решении распространяются на всю область значений переменной x.

Произведем аналогичный расчет при = 0.1, = max (0.1, 1)h2. На x рис. 14 приведена зависимость от x полученной в расчете функции F при y = 0.48, t = (непомеченная кривая) и t = 2 (кривая, помеченная Рис. кружками). Видно, что в этом случае осцилляций нет.

Определим теперь величину max как наибольшее значение, при котором решение при t =, y = 0.48 остается монотонным, и рассчитаем зависимость величины max = max /h2 от. Полученная зависимость x показана на рис. 15. (При = 0 max = 0.5, при = 0.5 max = 1.5.) Она обрывается на значении = 0.9, так как при = 0.9 max скачком увеличивается от значения 39 до 550 и затем плавно уменьшается до значения 500 при = 1. Это говорит о том, что с точки зрения уменьшения осцилляций при заданном лучше использовать схему с, близким к 1.

Рис. ГЛАВА ЗАДАЧА О ВРАЩАЮЩЕМСЯ КОНТЕЙНЕРЕ В этой главе рассматривается стационарное плоскопараллельное течение вязкой несжимаемой жидкости, частично заполняющей горизонтальный бес конечно длинный вращающийся цилиндр. Область течения предполагается односвязной, т. е. имеются две линии скользящего трехфазного контакта.

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

Корректность условий подобного рода с математической и физической то чек зрения в случае наличия подвижных линий контакта исследовалась в работах [4, 6] и литературе, цитированной в [4]. Результаты главы изложены в 5 главе монографии [61].

4.1. Постановка задачи Картина течения в плоскости, ортогональной оси цилиндра, показана на рис. 16. Здесь D область, занятая жидкостью, 1 граница раздела жидкость – газ, 2 граница раздела жидкость – твердая стенка, A и B точки трехфазного контакта. Центр декартовой системы координат xOy лежит на оси цилиндра, ось Oy направлена вертикально вверх, y = Rh соответствует среднему уровню жидкости в цилиндре.

4.1.1. Для определения течения, вообще говоря, необходимо найти форму свободной границы 1 и поле скоростей u = (u, v) в D. Задача ставится следующим образом. Внутри D считаются выполненными уравнения Навье – Стокса 2u 2u u u 1 p u +v = + +, x2 y x y x Рис. 2v 2v v v 1 p u +v = + + g0, (4.1.1) x2 y x y y u v + = 0. (4.1.2) x y На границе ставятся условия u · n = 0, (x, y) 1 2, (4.1.3) s · P · n = 0, (x, y) 1, (4.1.4) s · P · n = k(u · s R), (x, y) 2, (4.1.5) n · P · n = K, (x, y) 1, (4.1.6) s1 · s2 = cos(0 ), (x, y) 1 2. (4.1.7) Здесь u, v компоненты вектора скорости по осям Ox, Oy, p функция давления, плотность, g0 модуль ускорения свободного падения, коэффициент кинематической вязкости, P тензор напряжений, P = pI + 2S, (4.1.8) где I единичный тензор, S тензор скоростей деформаций, u 1 v + u x 2 x y S=, (4.1.9) 1 v u v + 2 x y y n единичный вектор внешней нормали к границе D, s вектор каса тельной, получающийся поворотом n по часовой стрелке на 90o, s1 и s векторы касательных, определенные по указанному выше правилу, к участ кам 1 и 2 границы D в точках контакта, 0 угол контакта, K кривизна свободной границы (K 0 соответствует выпуклости жидкости наружу), k коэффициент проскальзывания, угловая скорость вра щения цилиндра, R радиус цилиндра, коэффициент поверхностного натяжения.

В постановке задачи (4.1.1) система уравнений импульса, (4.1.2) уравнение неразрывности, (4.1.3) условие непротекания, (4.1.4) усло вие отсутствия касательных напряжений, (4.1.5) условие проскальзыва ния, (4.1.6) условие капиллярного равновесия, (4.1.7) условие, задающее угол 0 подхода жидкости к твердой стенке.

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

4.1.2. Перейдем к переменным "функция тока вихрь " по фор мулам v u u=, v=, =. (4.1.10) y x x y Тогда уравнение неразрывности выполняется автоматически, а из уравнений импульса следует =, 1 =, (x, y) D, (4.1.11) x y y x 2 = + y2.

где x В силу (4.1.3) функция тока имеет на границе D постоянное значение, которое можно положить равным нулю = 0, (x, y) 1 2. (4.1.12) В [30] показано, что условие (4.1.4) в терминах, можно записать в виде + 2K = 0, (x, y) 1. (4.1.13) n Для пояснения физического смысла этого условия заметим, что частица на 1 имела бы значение вихря, определяемое уравнением (4.1.13), если бы была частью твердого тела, вращающегося вокруг центра кривизны соот ветствующей точки 1. При наличии дополнительной составляющей вихря возникает взаимное скольжение слоев жидкости, которое и приводит, за счет действия сил вязкости, к возникновению касательных напряжений.

Вычисления, аналогичные проведенным в [30] при выводе (4.1.13), пока зывают, что на всей границе D имеет место равенство s · P · n = + 2K, n где K кривизна границы. Подставляя его в (4.1.5), получим 2 + = k R, (x, y) 2. (4.1.14) R n n В [30] доказано, что n·S·n=. (4.1.15) s n С учетом (4.1.8) и (4.1.15) условие (4.1.6) примет вид p 2 = K, (x, y) 1. (4.1.16) s n Пусть y = Rh0 + f (x) уравнение 1. Введем новую функцию давления p, исключив гидростатическую составляющую p = p g0 (y Rh0 ).

(4.1.17) При замене (4.1.17) в уравнениях (4.1.1) исчезнет член, отвечающий массо вым силам. Проецируя систему (1.1) на направление касательной s к 1, найдем, что p =. (4.1.18) s 2 s n n Проинтегрировав (4.1.18) по s, подставив полученное выражение для p в (4.1.17), а затем выражение для p в (4.1.16), выразив кривизну K через производные функции f, получим окончательную форму условия (4.1.6) s + g0 f + d s +C = 2 n n s n s f =, (x, y) 1, (4.1.19) (1 + f 2 )3/ где C константа интегрирования, подлежащая определению из закона сохранения массы, штрих означает дифференцирование по x.

Наконец, условие (4.1.7) можно представить в виде f f = sign(x) tg 0 + arccos h0 +, (x, y) 1 2. (4.1.20) R Таким образом, задача (4.1.1)–(4.1.7) эквивалентна задаче (4.1.11)– (4.1.14), (4.1.19), (4.1.20), сформулированной в терминах функций,, f.

4.1.3. Произведем анализ значений параметров. Перейдем к безразмер ным переменным u = U u, v = U v, x = R, y = R, f = Rf, p = U 2 p, x y (4.1.21) где U единица масштаба скорости. В результате, кроме 0, h0, в задаче появятся еще 4 безразмерных параметра U2 g0 R UR kR Re =, Fr =, Bo =, Al =, (4.1.22) g0 R где Re число Рейнольдса, F r число Фруда, Bo число Бонда, Al параметр проскальзывания (динамическое число Нуссельта [6]). Первое уравнение системы (4.1.11) и условия (4.1.12), (4.1.13) сохранят свою форму, во втором уравнении системы (4.1.11) множитель 1/ заменится множи телем Re, в условии (4.1.20) исчезнет множитель 1/R, условие (4.1.14) преобразуется к виду R +2 = Al, (4.1.23) n n U а условие (4.1.19) к виду 1 f (x) f (x) 3 = (x) + f0, (4.1.24) Bo [1 + (f (x))2 ] где s F r Fr (x) = d s, (4.1.25) 2 n Re n s n s f0 константа, подлежащая выбору из закона сохранения массы. (Здесь и ниже волны над безразмерными переменными опускаются.) В дальнейшем будем полагать U = R (это определяет значение неод нородного члена в (4.1.23)), но в этом пункте при оценке параметров примем U = Um, где Um = max(|u| : (x, y) D). (4.1.26) Задача была поставлена для моделирования течения расплава полупро водникового материала. Примем характеристики расплава германия: = 5.57 г/см3, = 1.35·103 см2 /c. Типичные значения других параметров сле дующие: g0 = 10 см/с2, R = 1 см, R = 30 см/с, Um = 3 см/с. Из (4.1.22) получим Re 2 · 103, F r 102, Bo 9. Но мы ограничимся скоростями = = = движения, при которых Re 100, т. е. предположим, что Um 0.135 см/с.

Тогда получим F r 1.8 · 105 и, в силу (4.1.25), такой же порядок имеет функция (x). Ввиду малости алгоритм решения задачи может быть следующим: 1) решим задачу (4.1.24), (4.1.20) при (x) 0, т. е. опреде лим форму свободной границы в условиях гидростатики;

2) используя все уравнения и граничные условия, кроме (4.1.24), (4.1.20), найдем функции, в области с фиксированной границей;

3) уточним форму свободной границы, решив задачу (4.1.24), (4.1.20) относительно f с функцией (x), рассчитанной по выражению (4.1.25).

Ниже поставленная задача решается для частного случая, когда парамет ры h0 и 0 связаны соотношением h0 = cos( 0 ), (4.1.27) обеспечивающим горизонтальный подход жидкости к стенкам. В этом слу чае решение задачи гидростатики тривиально f (x) 0, 1 h2 x 1 h2. (4.1.28) 0 В условии (4.1.13) можно тогда положить K = 0, а условие (4.1.24) линеа ризовать, отбросив малую второго порядка (f (x))2.

4.1.4. Перейдем к новой системе координат. При допущении (4.1.27) в силу (4.1.28) D является сегментом круга. Следующее преобразование кон формно отображает ее на бесконечную полосу:

+ i x + iy = ih0 + x0 coth, (4.1.29) где 1 h2.

x0 = (4.1.30) Разделяя действительную и мнимую части, получим sinh sin x = x0, y = h 0 x0. (4.1.31) cosh cos cosh cos Система координат (, ) называется биполярной. Для решения данной задачи она была предложена О. М. Лаврентьевой.

При замене (4.1.29) граница 1 перейдет на плоскости (, ) в прямую =, граница 2 в прямую = 0, где 0 = arccos h0, (4.1.32) точки контакта (x, y) = (±x0, h0 ) отобразятся, соответственно, на точки = ±.

С учетом замечаний, сделанных в предыдущем пункте, задача для опре деления, в области с фиксированной границей примет в биполярных координатах вид = H 2, (4.1.33) = Re, (, ) D, (4.1.34) (, ) = (, 0 ) = 0, (4.1.35) (, ) = 0, (4.1.36) (, 0 ) = Al + (2 Al), (4.1.37) H0 () = где d(x + iy) x H(, ) = =, (4.1.38) d( + i) cosh cos H0 () = H(, 0 ), (4.1.39) 2 = 2 + 2.

После линеаризации уравнения (4.1.24) и граничных условий в точках контакта получим следующую задачу для определения функции f и кон станты f0 :

f (x) f (x) = (x) + f0, x (x0, x0 ), (4.1.40) Bo f (x0 ) f (x0 ) f (x)|x=x0 =, f (x)|x=x0 =, (4.1.41) x0 x x f (x)dx = 0, (4.1.42) x где F r 1 Fr 2 (x()) = d.

2 H Re H H 0 = (4.1.43) (Здесь (4.1.42) закон сохранения массы жидкости.) Первая задача должна решаться численно. Для этого бесконечную поло су следует заменить прямоугольником и поставить на его торцах некоторые условия. Вторая задача решается аналитически, а входящие в решение ин тегралы от функции рассчитываются численно.

4.2. Асимптотика вихря и функции тока в окрестности точек кон такта 4.2.1. Для постановки задачи в ограниченной области построим асимпто тические разложения функций, при ||. При этом воспользуемся интегральной формой записи задачи (4.1.33)–(4.1.37).

Функция Грина для краевой задачи в области D с оператором Лапласа в биполярных координатах имеет вид 1 cosh (1 ) cos ( + 1 20 ) G(,, 1, 1 ) = ln, (4.2.1) 4 cosh (1 ) cos ( 1 ) где =. (4.2.2) Производная функции Грина по внешней нормали к 1 G g(1, 1, ) = = h(1, ), (4.2.3) H0 () H0 () = где sin ( 0 ) h(, ) =. (4.2.4) 2 (cosh cos ( 0 )) При вычислениях будут использоваться представления G и h в виде рядов ek|1 | G(,, 1, 1 ) = sin k( 0 ) sin k(1 0 ), (4.2.5) k k= ek|| sin k( 0 ).

h(, ) = (4.2.6) k= При известной функции решение задачи (4.1.33), (4.1.35) выписывается в виде G(,, 1, 1 )(1, 1 )H 2 (1, 1 )d1 d1.

(, ) = (4.2.7) D Скорость жидкости на границе 1 (, ) v() =. (4.2.8) H0 () = Считая правую часть (4.1.34) и v() известными функциями, выпишем решение задачи (4.1.34), (4.1.36), (4.1.37) относительно Al (, ) = 0 () + (2 Al)Ihv (, ) + ReIG (, ), (4.2.9) где 0 () = 2 h(1, )d1 = 2, (4.2.10) Ihv (, ) = h(1, )v(1 )d1, (4.2.11) IG (, ) = G(,, 1, 1 ) d1 d1. (4.2.12) 1 1 1 D (4.2.7)–(4.2.12) представляет собой систему интегральных уравнений отно сительно,.

4.2.2. Для вывода асимптотики докажем три вспомогательных утвержде ния.

Условимся высказывание "существует константа c 0, обеспечивающая выполнение неравенства |f1 ()| c|f2 ()| при 1 2 записывать в виде "f1 () = O(f2 ()), 1 2 а если 1 =, 2 =, то просто "f1 () = O(f2 ())".

Предложение 4.2.1. Пусть вещественные функции f () и g() та ковы, что |f ()|d, (4.2.13) c ean||, ||, f () = (4.2.14) n n= M d ebm || + gbM (), ||, g() = (4.2.15) m m= где gbM () = O(ebM || ). (4.2.16) = sign();

, a, c1, c+1, (n = Здесь M натуральное число, n n 1, ), d1, d+1, bm (m = 1, M ) вещественные числа, причем 0, a m m 0, 0 b1... bM, bm = an (m = 1, M, n = 1, );

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

Тогда функция If g () = f (1 )g(1 )d1 (4.2.17) представима на всей числовой прямой в виде N 1 M c A ean|| d Bm ebm || + IbM (), If g () = + (4.2.18) n n m n=1 m= где bM + 1, IbM () = O(ebM || ), N= (4.2.19) a (квадратные скобки означают взятие целой части), причем коэффициен ты A зависят от поведения g, коэффициенты Bm от поведения f n на всей числовой прямой.

Доказательство. Согласно (4.2.13), (4.2.14) имеем N c ean|| + faN () + f (), ||, f () = (4.2.20) n n= где faN () = O(eaN || ), (4.2.21) f () = 0, ||, |f ()|d. (4.2.22) Подставим (4.2.15), (4.2.20) в (4.2.17) и произведем интегрирование, предпо лагая, что ||. В результате получим (4.2.18), если обозначить M d d m m A = ean gbM ()d, + (4.2.23) n an + bm an bm m=1 N c c n n ebm faN ()d + ebm f ()d, Bm = + bm + an bm an n=1 (4.2.24) + IbM = faN (1 )gbM (1 )d1 + f (1 )gbM (1 )d1 +... (4.2.25) Интегралы в правых частях (4.2.23), (4.2.24) сходятся в силу выбора N (из a(N 1) bM aN следует, что an a(N 1) bM, bm bM 1 bM aN ) и интегрируемости f. В (4.2.25) выписаны старшие члены функции IbM (). Учитывая (4.2.16), (4.2.21), нетрудно получить оценку |IbM ()| cebM ||. Таким образом, справедливость (4.2.18) установлена для ||.

Доопределим функцию IbM () при || формулой (4.2.18). В силу условий предложения If g ограничена. Следовательно, IbM () = O(ebM || ) на всей числовой прямой. Предложение доказано.

Замечание. Если в условиях предложения 4.2.1 при некоторых n = n0 и m = m0 M имеет место равенство an0 = bm0, то в правую часть (4.2.18) добавится член d 0 c ||ean0 ||, а выражения для A 0, Bm изменят m n0 n свой вид. Если aN = bM, то остаточный член в (4.2.18) запишется в виде IbM = O(||eaN || ).

Предложение 4.2.2. Пусть 0 a = k (k = 1, ). Тогда sin a( ) a|| a|1 | bk ek|| sin k( 0 ), h(1, )e d1 = e + sin a( 0 ) k= (4.2.26) где bk вещественные коэффициенты.

Доказательство. Подставим (4.2.6) в (4.2.26) и вычислим интеграл по внешности отрезка [, + ]. Получим, что равенство (4.2.26) эквива лентно следующему:

sin a( ) a|| a|1 | k ()ea|| sin k( 0 ) = h(1, )e d1 + e, sin a( 0 ) k= (4.2.27) где 2ek (k cosh(a) + sinh(a)) k () =.

(k 2 2 a2 ) Так как начиная с некоторого k коэффициенты k () монотонно и равно мерно относительно стремятся к нулю при k, а последовательность частичных сумм от sin k( 0 ) ограничена, то ряд в левой части (4.2.27) сходится равномерно относительно (по признаку Дирихле), что делает законным предельный переход под знаком суммы при 0. В результа те получаем верное равенство, так как первый член в левой части (4.2.27) стремится к нулю при 0 ввиду ограниченности h, а второй член яв ляется рядом Фурье для функции, стоящей в правой части, продолженной нечетным образом на отрезок (20, ). Предложение доказано.

Предложение 4.2.3. Имеет место равенство Al G(,, 1, 1 )4x2 e2|1 | 0 (1 )d1 d1 = D sin Alx2 2|| bk ek|| sin k( 0 ), = e (4.2.28) 0 sin k= где bk вещественные коэффициенты.

Доказательство этого утверждения после подстановки (4.2.5), (4.2.10) в левую часть (4.2.28) проводится аналогично доказательству предложения 4.2.2.

4.2.3. Перейдем непосредственно к вычислению асимптотики функций,. Будем последовательно применять предложение 4.2.1 для оценки правых частей в выражениях (4.2.7), (4.2.11), (4.2.12). Старшинство получа емых членов зависит от значения параметра, который в общем случае лежит в пределах 1. Ниже предполагается, что 2. (4.2.29) Этот случай отвечает заполнению цилиндра более чем наполовину (см.

(4.1.32), (4.2.2)).

Предположим, что |(, )| const, (, ) D. Учитывая (4.1.38), найдем H 2 = O(e2|| ). Применив к (4.2.7) предложение 4.2.1 (с учетом (4.2.5)), получим = c sin ( 0 ) · e|| + O(e2|| ), = c sin ( 0 ) · e|| + O(e2|| ), = c cos ( 0 ) · e|| + O(e2|| ).

(при оценке производных равенство (4.2.7) дифференцируется, в пра вой части берется интеграл по, а затем применяется предложение 4.2.1), где = sign(), c+1, c1 константы, зависящие от поведения функции во всей области D. Из (4.2.8), с учетом (4.1.39), следует v = (2x0 )1 c e(1)|| + O(e|| ). Применяя к (4.2.11) предложения 4.2.1, 4.2.2 (с учетом (4.2.6)), будем иметь Ihv = c sin(( 1)( ))e(1)|| + O(e|| ), где c c =. (4.2.30) 2x0 sin(( 1)( 0 )) Оценим член IG. Чтобы не использовать дополнительных предположе ний о производных функции, проинтегрируем правую часть (4.2.12) по частям. Получим G G IG = d1 d1. (4.2.31) 1 1 1 D Так как / = O(e|| ), / = O(e|| ), то на основании заме чания к предложению 4.2.1 заключаем IG = O(||e|| ). (4.2.32) Подставляя оценки для Ihv и IG, получим Al 0 () + (2 Al)[c sin(( 1)( ))e(1)|| + O(e|| )]+ = +ReO(||e|| ). (4.2.33) Теперь, используя (4.2.33), проделаем описанную цепочку оценок еще раз.

В результате получим новые члены в разложениях,. (При оценке используется предложение 4.2.3.) Выпишем старшие из них = c sin ( 0 ) · e|| + 2 + (2 Al)(O(e(+1)|| ) + O(e3|| )) + O(e2|| ), (4.2.34) где sin 2 = Alx2 e2||.

(4.2.35) 0 sin Al 0 () + (2 Al)[c sin(( 1)( ))e(1)|| + 1 ]+ = +(|2 Al) + Re)O(||e|| ), (4.2.36) где Alx0 1 sin( ) || 1 = + 2 cot 20 e. (4.2.37) 2 0 sin( 0 ) При численном решении задачи для, в прямоугольнике (, ) [0, 0 ] [0, ] на границах = ±0 можно ставить "мягкие" граничные условия, получающиеся исключением из выражений (4.2.34), (4.2.36) неиз вестных констант c и c 1 посредством дифференцирования Al Al 0 (2 Al)1 + ( 1) 0 (2 Al)1 = 0, || 2 (4.2.38) ( 2 ) + ( 2 ) = 0, || = 0. (4.2.39) || В специальном случае Al = 2 условие (4.2.38) можно заменить на Al = 0, || = 0. (4.2.40) Если удовлетворяет (4.2.29), то точность условий (4.2.38)–(4.2.40) опре деляется порядком остаточных членов в (4.2.34), (4.2.36). Если 2, то выписанные явно в (4.2.34), (4.2.36) члены будут по-прежнему старшими, но порядок остаточных членов будет другим. Его можно повысить, если вы числить некоторые дополнительные слагаемые типа 2 и 1.


4.3. Аналитические решения При нулевом числе Рейнольдса задача (4.1.33)–(4.1.37) соответствует на плоскости переменных (x, y) краевой задаче в сегменте круга для бигармо нического уравнения. Получено два класса ее аналитических решений. Пер вый отвечает специальному значению Al, при котором упрощается условие (4.1.37). Второй специальному значению параметра заполнения, при кото ром область D становится полукругом. В последнем случае найден также предел решения при Al, соответствующий постановке на 2 усло вий прилипания. В пункте 4.4.2 будет показано, что для этого решения не существует функции f (x), удовлетворяющей задаче (4.1.40)–(4.1.43).

4.3.1. Пусть Al = 2. Тогда функция (4.2.10) в силу (4.2.9) является решением задачи (4.1.34), (4.1.36), (4.1.37) для. (Изолинии функции представляют собой в переменных (x, y) дуги окружностей, проходящие через точки контакта.) Для нахождения функции тока сделаем замену = 1, (4.3.1) где x2 ( ) cos + ln(cosh cos ) sin 1 =. (4.3.2) ( 0 ) cosh cos Функция 1 удовлетворяет условиям 1 = H 2 0, (, ) D, 1 (, ) = 0.

Поэтому в силу (4.1.33), (4.1.35) функция является решением задачи = 0, (, ) D, (, ) = 0, (, 0 ) = 1 (, 0 ).

Оно выписывается в виде (, ) = h(1, )1 (1, 0 )d1, где h определена в (4.2.4), или в эквивалентном виде 1 (1, 0 ) (, ) = h1 (1, ) d1, (4.3.3) где 1 ( 0 ) h1 (, ) = h(, )d = arctan tanh cot. (4.3.4) 2 Таким образом, решение задачи (4.1.33), (4.1.35) для в рассматриваемом случае дается формулами (4.3.1)–(4.3.4).

4.3.2. Рассмотрим теперь общий случай (параметры Al и произ вольны). Сведем систему (4.2.7), (4.2.9) к одному интегральному уравнению.

Подставляя (4.2.7) в (4.2.8) и учитывая (4.2.3), найдем, что g(1, 1, )(1, 1 )H 2 (1, 1 )d1 d1.

v() = (4.3.5) D Подставляя (4.3.5) в (4.2.9) (при Re = 0) и меняя порядок интегрирования, получим следующее интегральное уравнение для вихря в D :

Al F (,, 1, 1 )(1, 1 )H 2 (1, 1 )d1 d1 + (, ) = (2 Al) 0 (), D (4.3.6) где F (,, 1, 1 ) = g(,, 2 )g(1, 1, 2 )H0 (2 )d2. (4.3.7) (Если подставить (4.2.9) в (4.3.5), то получится интегральное уравнение для скорости v на 2.) Установим свойства собственных функций и собственных значений ядра F.

Пусть fk (, ) (k = 1, ) полная ортонормированная в D система гармонических в D функций, обращающихся в ноль на 1 :

fk fn H 2 dd = kn.

fk = 0, (, ) D, fk |1 = 0, (4.3.8) D Тогда имеет место разложение g(1, 1, ) = vk ()fk (1, 1 ), (4.3.9) k= где g(1, 1, )fk (1, 1 )H 2 (1, 1 )d1 d1.

vk () = (4.3.10) D Потребуем ортогональности vk на vk ()vn ()H0 ()d = k kn. (4.3.11) Подставляя (4.3.10) в (4.3.11), получим F (1, 1, 2, 2 )fn (2, 2 )H 2 (2, 2 )d2 d D D fk (1, 1 )H 2 (1, 1 )d1 d1 = k kn. (4.3.12) Сравнивая (4.3.12) с последним из равенств (4.3.8), приходим к выводу, что F (1, 1, 2, 2 )fn (2, 2 )H 2 (2, 2 )d2 d2, n fn (1, 1 ) = (4.3.13) D т. е. k собственные значения, а fk собственные функции ядра F.

Установим связь между функциями fk и vk. Имеем F = g(1, 1, )g(2, 2, )H0 ()d = k fk (1, 1 )fk (2, 2 ). (4.3.14) k= Воспользуемся свойством функции g, следующим из ее определения (4.2.3):

lim f ()g(1, 1, )H0 ()d = f (1 ) (4.3.15) 1 Переходя в (4.3.14) к пределу при 1 0, с учетом (4.3.15) получим g(2, 2, 1 ) = k fk (1, 0 )fk (2, 2 ). (4.3.16) k= Сравнивая (4.3.16) с (4.3.9), находим, что vk () = k fk (, 0 ). (4.3.17) Из (4.3.17), (4.3.11) получаем kn fk (, 0 )fn (, 0 )H0 ()d =. (4.3.18) k Таким образом, собственные функции fk ядра F представляют собой полную систему гармонических и ортонормированных в D, ортогональных на 2 обращающихся в ноль на 1 функций. Если они известны, то собственные значения k ядра F могут быть найдены из соотношений (4.3.18), и решение уравнения (4.3.6) выписывается в виде Ak k hk (, ) = fk (, ), (4.3.19) 2 1 (2 Al)k k= где hk = 2 fk (, 0 )H0 ()d. (4.3.20) Функция тока определяется по формуле Al k hk (, ) = gk (, ), (4.3.21) 2 1 (2 Al)k k= где G(,, 1, 1 )fk (1, 1 )H 2 (1, 1 )d1 d1.

gk (, ) = (4.3.22) D 4.3.3. Рассмотрим случай = 2. Тогда область D является полукругом.

Введем полярные координаты x = r cos, y = r sin. (4.3.23) Область D отобразится на прямоугольник (r, ) (0, 1) (, 0). Кривая 1 опишется соотношениями =, r [0, 1];

= 0, r [0, 1], а кривая 2 соотношениями r = 1, [, 0]. Функция Грина в полярных координатах примет вид G(r,, r1, 1 ) = (r2 r1 2rr1 cos( 1 ) + 1)(r2 2rr1 cos( + 1 ) + r1 ) 2 = ln 2 2 2= 4 (r r1 2rr1 cos( + 1 ) + 1)(r2 2rr1 cos( 1 ) + r1 ) r1 k (rr1)k sin k sin k, r r1 ;

при r k k= = (4.3.24) k k sin k sin k r (rr1), r r1.

при r1 k k= Производная функции Грина по внешней нормали к G g(r,, 1 ) = r1 r1 = 1 r2 1 r = = r2 2r cos( 1 ) + 1 r2 2r cos( + 1 ) + sin k sin k1 · rk.

= (4.3.25) k= Нетрудно видеть, что (4.3.25) является конкретным видом формулы (4.3.16).

Действительно, функции 4(k + 1) k fk (r, ) = r sin k (4.3.26) гармонические и ортонормированы в D, ортогональны на 2 и обраща ются в ноль на 1, т. е. являются собственными для ядра F. Построим решение задачи в соответствии с формулами (4.3.18)–(4.3.22). Имеем 0 k = fk (1, )d =, (4.3.27) 2(k + 1) 4(k + 1) (1 (1)k ) hk = 2 fk (1, )d = 2, (4.3.28) k gk (r, ) = G(r,, r1, 1 )fk (r1, 1 )r1 dr1 d1 = D 4(k + 1) 1 r2 rk sin k.

= (4.3.29) 4(k + 1) Для вихря и функции тока получим следующие выражения:

16Alm r2m1 sin(2m 1), (r, ) = (4.3.30) (2m 1)(4m + Al 2) m= 2Al 1 r2 r2m1 sin(2m 1).

(r, ) = (2m 1)(4m + Al 2) m= (4.3.31) При Al = 2 ряд (4.3.30) суммируется 4 2r sin (r, ) = arctan, (Al = 2). (4.3.32) 1 r Это решение соответствует функции 0 () в биполярных координатах. При Al = после суммирования рядов (4.3.30), (4.3.31) будем иметь 1 r2 2r sin (r, ) = arctan, (4.3.33) 1 r (r, ) = 1, 0, (4.3.34) r r= 8 r(1 + r ) sin 4 2r sin (r, ) = + arctan, (4.3.35) 1 r (1 r2 )2 + 4r2 sin (1, ) = 2, 0, (Al = ). (4.3.36) sin Из (4.3.34) видно, что при Al = на 2 будет выполняться условие прилипания. Из (4.3.36) ясен характер особенности у функции на 2 в окрестности точек контакта.

4.4. Метод и результаты расчета течения Численное решение задачи в общем случае (при Re 0) производится в два этапа. Сначала находятся функции, в фиксированной области, затем рассчитывается функция f (x), описывающая форму свободной границы, как об этом говорилось в п. 4.1.3, 4.1.4.

4.4.1. Расчет функций, производится методом установления по вре мени: уравнения (4.1.33), (4.1.34) заменяются на = + H 2, (4.4.1) t = Re, (, ) D. (4.4.2) t2 Причем для и поставлены некоторые начальные условия и гра ничные условия (4.1.35)–(4.1.37), (4.2.38), (4.2.39). Сначала решается первое уравнение, пока итерации по времени не сойдутся, затем второе и т. д. Для решения каждого из уравнений используется метод переменных направле ний [33]. Рассмотрим общее уравнение U = U + L3 U + L4 U + f. (4.4.3) t Первое уравнение получается из него при U =, L3 = L4 = 0, f = H 2.

Второе уравнение получается, если подставить U =, L3 = Re, L4 = Re, f = 0.

Запишем уравнение (4.4.3) в разностном виде, расщепив его по направле ниям, следующим образом:

k+1/2 k Unm Unm k k+1/ = (1 + 3 + (1 )4 )Unm + (2 + 4 )Unm + fnm, k k+1/ k+ Unm Unm k+1/2 k+ = (2 + 4 + (1 )3 )Unm + (1 + 3 )Unm + fnm, (4.4.4) k где Unm1 2Unm + Unm+1 Un1m 2Unm + Un+1m 2 Unm =, 1 Unm =, (4.4.5) 2, шаги сеток по осям, ;

весовой коэффициент, k, k, (k = 1, K) итерационные параметры, K число итераций. Число итераций при расчетах определяется эмпирически, при решении задачи для в каче стве k, k выбираются оптимальные итерационные параметры по Жор дану [35]. При решении задачи для в качестве итерационных выбираются те же параметры, умноженные на (Re + 1)1 или на другой эмпирический коэффициент.

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

n+1m n1m nm+1 nm = unm + vnm 2 n+1m 2nm + n1m nm+1 2nm + nm |unm | |vnm |, (4.4.6) 2 где nm+1 nm1 n+1m n1m unm =, vnm =.

2 Система решается методом прогонки.

Расчет ведется методом продолжения по числу Рейнольдса.

4.4.2. Найдем решение краевой задачи для определения формы свобод ной границы. Введем обозначение = Bo. (4.4.7) Представим искомые функцию f (x) и константу f0 в виде f (x) = f1 (x) + f2 (x), f0 = f10 + f20, (4.4.8) где функции fj (x) и константы fj0 (j = 1, 2) определяются из решения следующих краевых задач:

f1 (x) f (x) = (x) + f10, x0 x x0, (4.4.9) 2 f1 (x)|x=x0 = 0, f1 (x)|x=x0 = 0, (4.4.10) x f1 (x)dx = 0, (4.4.11) x f2 (x) f (x) = f20, x0 x x0, (4.4.12) 2 f2 (x)|x=x0 = a, f2 (x)|x=x0 = b, (4.4.13) x f2 (x)dx = 0, (4.4.14) x где a и b некоторые константы.

Решение задачи (4.4.9), (4.4.10) выписывается в виде x f1 (x) = G(x, s)(s)ds + f10, (4.4.15) x где cosh (s+x0 ) cosh (xx0 ), x0 s x x0, sinh 2x G(x, s) = (4.4.16) cosh (sx0 ) cosh (x+x0 ), x0 x s x sinh 2x функция Грина для этой задачи. Константа f10 определяется из условия (4.4.11) x0 x f10 = G(x, s)(s)dxds. (4.4.17) 2x x0 x Решение задачи (4.4.12), (4.4.13) дается формулой (b a) cosh x (b + a) sinh x f2 (x) = + + f20. (4.4.18) 2 sinh x0 2 cosh x Константа f20 определяется из условия (4.4.14) (b a) f20 =. (4.4.19) 2 2 x Подставляя первое из равенств (4.4.8) в (4.1.41) и учитывая (4.4.10), (4.4.13), (4.4.18), (4.4.19), получим алгебраическую систему относительно a и b, решая которую, находим x0 2 (f1 (x0 ) + f1 (x0 )) ba=, (4.4.20) 1 x0 coth x0 + 2 x2 (f1 (x0 ) f1 (x0 )) b+a=. (4.4.21) x0 tanh x Заметим, что при 0 знаменатели в правых частях (4.4.20), (4.4.21) сирого положительны.


Таким образом, решение задачи (4.1.40)–(4.1.42) дается формулами (4.4.15)–(4.4.21). При Al функция (x) имеет в точках контакта интегрируемую особенность и, вследствие этого, f тоже. Вместе с тем функции f (x) и f (x) в этих точках непрерывны. Для решения (4.3.33), (4.3.35), построенного при Al =, особенность функции (x) в точках контакта не интегрируема (это видно из (4.3.35)). Поэтому интеграл в правой части (4.4.15) расходится и решения задачи (4.1.40)–(4.1.42) не существует.

4.4.3. Приведем пример численного эксперимента при значениях пара метров 0 = 0.5 рад (соответственно h = 0.88, т. е. цилиндр заполнен почти полностью, = 1.19 (см. (4.1.32), (4.2.2)), Al = 1.5. Параметры Re и Bo будут варьироваться, а параметр F r не рассматривается, так как функция f (x) зависит от него линейно. Обратим внимание на следующие особенности расчета.

1) Среди асимптотических условий наибольшую ошибку дает условие для 103 в его выполнении следует вихря (4.2.38). Для достижения точности брать 0 10 12.

2) После введения сетки в области (, ) [0, 0 ] [0, ] соответству ющая часть области D в плоскости переменных (x, y) будет разбита на 2 + 2 H(, ), где ячейки, диаметры которых приближенно равны, шаги сеток по осям,. Функция H имеет максимум при = 0, = 0, равный 3.9. Чтобы обеспечить малость максимума диаметров ячеек, нужно либо делать очень большое число разбиений по оси (поряд ка 400, чтобы значение этого максимума было меньше 0.2 при 0 = 10), либо использовать неравномерную сетку.

3) При вычислении интегралов (4.4.15), (4.4.17) функция, определен ная формулой (4.1.43), должна быть доопределена вне отрезка || 0 с помощью асимптотических формул (4.2.34), (4.2.36). При этом неизвестные константы c+1, c1 приближенно определяются по решению внутри это го отрезка. Первый член в выражении (4.1.43) имеет порядок O(e(1)|| ), O(e(2)|| ). Интегралы от первых двух членов второй O(), третий по внешности отрезка || 0 пренебрежимо малы, а интеграл от третьего члена порядка единицы. При расчете функции f (x) он учитывался в виде аналитически вычисляемого слагаемого.

На рис. 17, а показаны линии тока, на рис. 17, б изолинии вихря, по строенные по решению задачи (4.1.33)–(4.1.37), рассчитанному при Re = 10.

На рис. 18, а, б приведены аналогичные графики для Re = 150. Функции и неположительны. Максимумы их модулей, соответственно, состав ляют 0.333 и 1.673 при Re = 10 и 0.335 и 1.671 при Re = 150. В точках контакта скорость жидкости равна нулю. Максимум скорости в D составляет 0.693 при Re = 10 и 0.764 при Re = 150. В первом случае он достигается на твердой стенке, во втором на свободной границе. Ре зультаты получены при 0 = 10 и числах разбиений по осям,, равных N = 64, N = 32. Сетка при расчете была неравномерной. Так, шаг по оси увеличивался от значения 0.03 при = 0 до значения 1. при || = 0. Для уменьшения влияния схемной вязкости разности против потока первого порядка при аппроксимации первых производных были заменены разностями против потока второго порядка точности. Все осталь ные члены в уравнениях и граничных условиях также аппроксимировались разностными формулами второго порядка. Система решалась методом пяти диагональной прогонки [34]. Тестовые расчеты на последовательности сеток показали, что численное решение имеет первый порядок сходимости.

Приведем результаты расчета свободной границы. Они были получены а б Рис. а б Рис. позже, и расчет функций,, необходимый для определения функции, производился другим методом по сравнению с описанным выше. А именно, он был сделан на равномерной сетке методом, описанным в п. 4.4.1. Числа разбиений были равны N = 192, N = 64. Расчет свободной границы про изводился по формулам п. 4.4.2 с учетом замечаний, сделанных в настоящем пункте. Тесты, проведенные на последовательности сеток, показали, что по лученные численно функции,, f имеют первый порядок сходимости.

Расчеты, проведенные при разных значениях 0, подтвердили, что функция правильно доопределена вне отрезка || 0.

Произведем замену f (x) = F (x) · F r/Re, f0 = F0 · F r/Re. (4.4.22) Из формул п. 4.4.2 и формулы (4.1.43) следует, что функция F (x) и кон станта F0 не зависят от параметра F r. На рис. 19 показаны графики функции F (x) при различных значениях параметров Re и Bo. Рис.

19, а, б соответствуют Re = 10, рис. 19, в, г Re = 150. Рис. 19, а, в построены при Bo = 9 (это соответствует реальным характеристикам рас плава германия см. п. 4.1.3), рис. 19, б, г Bo = 100. Расчеты функции F (x) при Bo = 0.5 показали, что в указанном диапазоне чисел Рейнольдса она визуально не отличается от наклоненной прямой линии. Максимум ее модуля составляет 8.47 при Re = 10 и 9.62 при Re = 150.

а б в г Рис. ЗАКЛЮЧЕНИЕ В диссертации получены следующие результаты.

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

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

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

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

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

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

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

СПИСОК ЛИТЕРАТУРЫ [1] Андреев В. К., Капцов О. В., Пухначев В. В., Радионов А. А. При менение теоретико-групповых методов в гидродинамике. Новосибирск:

Наука. 1994. 318 с.

[2] Анисютин Б. М. Численное исследование тепловой задачи для процесса бестигельной зонной плавки // Задачи гидромеханики со свободными границами: Межвуз. сб. науч. тр. / Новосиб. ун-т, Новосибирск, 1987, 128 с.

[3] Бадратинова Л. Г. О движении жидкого слоя по внутренней поверхно сти горизонтального вращающегося цилиндра // Динамика сплошной среды: Сб. науч. тр. / АН СССР. Сиб. отд-ние. Ин-т гидродинамики.

1993. Вып. 106. С. 179–184.

[4] Байокки К., Пухначев В. В. Задачи с односторонними ограничениями для уравнений Навье–Стокса и проблема динамического краевого угла // Прикладная механика и техническая физика. 1990. № 2. С. 27–40.

[5] Белова И. В., Волкова Г. Б. Термокапиллярная конвекция при бести гельной зонной плавке // Прикладная механика и техническая физика.

1993. № 6. С. 72–75.

[6] Богоряд И. Б. Динамика вязкой жидкости со свободной поверхностью.

Томск: Изд-во Том. ун-та. 1980. 102 с.

[7] Булеев Н. И. Пространственная модель турбулентного обмена. М.: На ука. 1989. 344 с.

[8] Веретенцев В. А. Построение разностной сетки в области с криволиней ными границами с помощью конформного отображения // Актуальные вопросы прикл. математики. М.: Изд-во МГУ. 1989. С. 88–93.

[9] Вигдорович В. Н. Совершенствование зонной перекристаллизации. М.:

Металлургия. 1974. 200 с.

[10] Воеводин А. Ф., Юшкова Т. В. Численный метод решения начально краевых задач для уравнений Навье–Стокса в замкнутых областях на основе метода расщепления // Сиб. журн. вычисл. математики. 1999. Т.

2. № 4. С. 321–332.

[11] Воеводин В. В., Кузнецов Ю. А. Матрицы и вычисления. М.: Наука.

1984. 320 с.

[12] Гончаров В. Л. Теория интерполирования и приближения функций. М.:

Гос. изд-во технико-теоретической лит-ры. 1954. 328 с.

[13] Дональд Д. К. Тепловой режим в условиях вакуумной плавки // При боры для научных исследований. 1961. № 7. С. 42–44.

[14] Ждан Л. А. Задача о движении вязкой жидкости во вращающемся круге в поле тяжести // Вестн. Моск. ун-та. Сер. 1. Мат., мех. 1987. № 1. С.

86–89.

[15] Зверев В. Г. Об одном итерационном алгоритме решения разностных эллиптических уравнений // Вычисл. технологии. 1999. Т. 4. № 1. С.

55–65.

[16] Зоммерфельд А. Электродинамика. М.: Изд-во иностранной лит-ры.

1958. 502 с.

[17] Калиткин Н. Н. Численные методы. М.: Наука. 1978. 512 с.

[18] Кочин Н. Е., Кибель И. А., Розе Н. В. Теоретическая гидромеханика.

Часть II. М.: Гос. изд-во физ.-мат. лит-ры. 1963. 727 с.

[19] Куликовский А. Г., Любимов Г. А. Магнитная гидродинамика. М.: Гос.

изд-во физ.-мат. лит-ры. 1962. 248 с.

[20] Лисейкин В. Д. Алгебраический метод построения разностных сеток.

Уч. пособие. Новосиб. гос. ун-т. Новосибирск. 2002. 46 с.

[21] Лисейкин В. Д. Метод алгебраической адаптации // Журн. вычисл.

математики и мат. физики. 1998. Т. 38. № 10. С. 1692–1709.

[22] Лойцянский Л. Г. Механика жидкости и газа. М.: Наука. 1973. 736 с.

[23] Марчук Г. И. Методы вычислительной математики. М.: Наука. 1977. с.

[24] Новиков И. И. Прикладная магнитная гидродинамика. М.: Атомиздат.

1969. 360 с.

[25] Овчарова А. С. Влияние эффектов плавучести и термокапиллярности на форму свободной поверхности расплава // Вычисл. технологии. 2002.

Т. 7. С. 323–328.

[26] Овчарова А. С. Метод расчета стационарных течений вязкой жидкости со свободной границей в переменных вихрь–функция тока // Приклад ная механика и техническая физика. 1998. Т. 39. № 2. С. 59–68.

[27] Овчарова А. С. Расчет жидкого моста при производстве монокристал лов методом бестигельной зонной плавки // Журн. вычисл. математики и математической физики. 2003. Т. 43. № 4. С. 1062–1071.

[28] Остапенко В. В. Метод теоретической оценки дисбаллансов неконсерва тивных разностных схем на ударной волне // Докл. АН СССР. 1987. Т.

295. № 2. С. 292–297.

[29] Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат. 1984. 150 с.

[30] Пухначев В. В. Задачи со свободной границей для уравнений Навье– Стокса. Дис.... докт. физ.-мат. наук. Новосибирск. 1974. 294 с.

[31] Пухначев В. В., Солонников В. А. К вопросу о динамическом краевом угле // Прикл. мат. и мех. 1982. Т. 46. Вып. 6. С. 961–971.

[32] Регель А. Р., Глазов В. М. Физические свойства электронных расплавов.

М.: Наука. 1980. 295 с.

[33] Роуч П. Вычислительная гидродинамика. М.: Мир. 1980. 616 с.

[34] Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений.

М.: Наука. 1978. 592 с.

[35] Самарский А. А. Теория разностных схем. М.: Наука. 1983. 616 c.

[36] Стечкин С. Б., Субботин Ю. И. Сплайны в вычислительной математике.

М.: Наука. 1976. 248 с.

[37] Тамм И. Е. Основы теории электричества. М.: Наука. 1976. 616 с.

[38] Хакимзянов Г. С., Шокин Ю. И., Барахнин В. Б., Шокина Н. Ю. Чис ленное моделирование течений жидкости с поверхностными волнами.

Новосибирск. Изд-во СО РАН. Ин-т вычисл. технологий. 2001. 394 с.

[39] Шрагер Г. Р., Козлобродов А. Н., Якутенок В. А. Моделирование гид родинамических процессов в технологии переработки полимерных ма териалов. Томск: Изд-во Томского ун-та. 1999. 230 с.

[40] Яковлев В. И. О влиянии термокапиллярных сил на начальный участок пленки расплава // Прикладная механика и техническая физика. 2001.

№ 2. С. 118–121.

[41] Яковлев В. И. О границах начального участка пленки расплава, фор мируемой при бестигельном зонном переплаве полупроводниковых ма териалов // Доклады Академии наук. 1999. Т. 369. № 6. С. 766–769.

[42] Яковлев В. И. О границах начального участка пленки расплава, фор мируемой при бестигельном зонном переплаве полупроводниковых ма териалов // Прикладная механика и техническая физика. 2000. № 3. С.

139–148.

[43] Яковлев В. И. О характере особенности переменных электромагнитных полей вблизи вершины проводящего клина // Прикладная механика и техническая физика. 1996. Т. 37. № 4. С. 3–8.

[44] Balmer R. T., Wang T. G. An experimental study of internal hidrocyts // Trans. ASME. Ser 1. J. of Fluids Eng. 1976. V. 98. № 4. P. 688–694.

[45] Betchelor G. K. On steady laminar ow with closed streamlines at large Reynolds number // J. Fluid Mech. 1956. V. 1. № 2. P. 177–190.

[46] Gans R. F., Yalisove S. M. Observations and Measurments of Flow in Partially-Filled Horizontally Rotating Cylinder // Trans. ASME. Ser.1. J.

of Fluids Eng. 1982. V. 104. № 3. P. 363–366.

[47] Greenspan H. P. On a rotational ow disturbed by gravity // J. Fluid. Mech.

1976. V. 74. Pt. 2. P. 335–352.

[48] Harten A. On a Class of High Resolution Total – Variation – Stable Finite – Dierence Schemes // SIAM / J. Numer. Analys. 1984. V. 21. № 1. P. 1–23.

[49] Kobayashi N. Power required to form a oating zone and the zone shapes // J. Crystal Growth. 1978. V. 43. P. 417–424.

[50] Lacis K., Muiznieks A., Ratnieks G. 3D mathematical model system for melt hydrodynamics in the silicon single crystal FZ-growth process with rotating magnetic eld // Magnetohydrodynamics. 2005. V. 41. № 2. P. 147–158.

[51] Lan C. W., Kou S. Eects of rotation on heat trasfer, uid ow and interfaces in normal gravity oating-zone crystal growth // J. Crystal Growth. 1991.

V. 114. P. 517–535.

[52] Lan C. W., Kou S. Heat transfer, uid ow and interface shapes in oating zone crystal growth // J. Crystal Growth. 1991. V. 108. P. 351–366.

[53] Lie K. H., Walker J. S., Riahi D. N. Free surface shape and AC electric current distribution for oat zone silicon growth with a radio frequency induction coil // J. Crystal Growth. 1990. V. 100. P. 450–458.

[54] Lie K. H., Walker J. S., Riahi D. N. Melt motion in the oat zone process with an axial magnetic eld // J. Crystal Growth. 1991. V. 109. P. 167–173.

[55] Muhlbauer A., Muiznieks A., Virbulis J. Analysis of the dopant segregation eects at the oating zone growth of large silicon crystals // J. Crystal Growth. 1997. V. 180. P. 372–380.

[56] Muhlbauer A., Muiznieks A., Virbulis J., Ludge A., Riemann H. Interface shape, heat transfer and uid ow in the oating zone growth of large silicon crystals with the needle-eye technique // J. Crystal Growth. 1995. V. 151.

P. 66–79.

[57] Raming G., Muiznieks A., Muhlbauer A. Numerical investigation of the inuence of EM-elds on uid motion and resistivity distribution during oating-zone growth of large silicon single crystals // J. Crystal Growth.

2001. V. 230. P. 108–117.

[58] Ratnieks G., Muiznieks A., Buliguns L., Raming G., Muhlbauer A. A. 3D numerical analysis of the inuence of EM and Marangoni forces on melt hydrodynamics and mass transport during FZ silicon crystal growth // Magnetohydrodynamics. 1999. V. 35. № 3. P. 223–236.

[59] Riahi D. N., Walker J. S. Float zone shape and stability with the electromagnetic body force due to a radio-frequency induction coil // J.

Crystal Growth. 1989. V. 94. P. 635–642.

[60] Surek T., Chalmers B. The direction of growth of the surface of a crystal in contact with its melt // J. Crystal Growth. 1975. V. 29. P. 1–11.

[61] Воеводин А. Ф., Остапенко В. В., Пивоваров Ю. В., Шугрин С. М. Про блемы вычислительной математики. Новосибирск: Изд-во Сиб. отд-ния РАН. 1995. 154 с.

[62] Пивоваров Ю. В. Одномерная тепловая задача о бестигельной зонной плавке в быстропеременном магнитном поле // Динамика сплошной сре ды: Сб. науч. тр. / РАН. Сиб. отд-ние. Ин-т гидродинамики. 1996. Вып.

111. С. 100–108.

[63] Пивоваров Ю. В. О построении ортогональной разностной сетки в кри волинейном четырехугольнике // Вычисл. технологии. 2003. Т. 8. № 5.

С. 94–101.

[64] Пивоваров Ю. В. Параметрический анализ задачи о бестигельной зон ной плавке в магнитном поле // Динамика сплошной среды: Сб. науч.

тр. / РАН. Сиб. отд-ние. Ин-т гидродинамики. 2000. Вып. 116. С. 142– 147.

[65] Пивоваров Ю. В. Расчет движения жидкости с переменной вязкостью в области с криволинейной границей // Вычисл. технологии. 2005. Т. 10.

№ 3. С. 87–107.

[66] Пивоваров Ю. В. Условия монотонности факторизованной разностной схемы для эволюционного уравнения с двумя пространственными пере менными // Вычисл. технологии. 2001. Т. 6. № 4. С. 81–91.

[67] Пивоваров Ю. В. Численное моделирование конвекции в плавающей зоне // Вычисл. технологии. 2006. Т. 11. № 1. С. 81–94.



Pages:     | 1 ||
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.