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

Соловей Георгий Олегович

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

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

Общественные науки в целом

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d364c5f1be77c40d58bef

UUID: 9828d027-70db-4b3b-a6df-5ac77e04a322

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 20

Соловей Георгий Олегович

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 2316424 bytes

Поделиться работой

САНКТ-ПЕТЕРБУРГСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ КАФЕДРА КОСМИЧЕСКИХ ТЕХНОЛОГИЙ И ПРИКЛАДНОЙ АСТРОДИНАМИКИ Соловей Георгий Олегович Дипломная работа Математическое моделирование конвективных течений в водоемах Научный руководитель, кандидат физ.-мат. наук, доцент Старков В.Н. Санкт-Петербург 2016

Содержание Введение. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 04 1. О тепловой конвекции. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 04 2. Уравнения Навье-Стокса в приближении Буссинеска. . . . . . . . . . . 07 3. Системы уравнений для функции тока. . . . . . . . . . . . . . . . . . . . . . . 08 4. Приближения для неглубоких водоемов. . . . . . . . . . . . . . . . . . . . . . 09 Глава1. Конвективные течения в атмосфере. . . . . . . . . . . . . . . . . . . . . . . . . 11 1.1 Ветры Земли. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.2 Описание морских бризов. . . . . . . . . . . . . . . . . . . . . .. . . . . . . . . . . 13 1.3 О природе термиков. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.4 Интегрирование уравнения общего вида. . . . . . . . . . . . . . . . . . . . 19 1.5 Симметричный случай. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.6 Несимметричный случай потока тепла сверху. . . . . . . . . . . . . . . . 22 1.7 Несимметричный случай температуры снизу. . . . . . . . . . . . . . . . 24 1.8 Выводы по первой главе. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Глава 2. Слой воды с источником тепла. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.1 Тепловое загрязнение. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2 Интегрирование уравнения общего вида. . . . . . . . . . . . . . . . . . . . 28 2.3 Трехслойная схема распределения источником тепла в объеме воды. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.4 Случай полной симметрии. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.5 Случай сильной асимметрии. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.6 Модель, учитывающая наличие источника тепла и стока тепла. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2

2.7 Выводы по второй главе. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Глава 3. Исследование конвективных течений в озере с учетом таяния льдины. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.1 Расчет таяния льдины в озере. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2 Однофазная задача. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.3 Двухфазная задача. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.4 Льдина на поверхности. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.5 Симметричный случай. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.6 Случай с двумя льдинами. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.7 Выводы по третьей главе. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4. Выводы. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 5. Литература. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 6. Приложение. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3

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

Рэлей [3] теоретически исследовал устойчивость равновесия в горизонтальном слое и определил крайние значения параметров конвекции для модельного случая слоя с обеими свободными границами. Дальнейшее развитие теории продвигалось весьма медленно из-за значительных вычислительных трудностей. В ряде работ рассматривались лишь некоторые усложнения задачи о горизонтальном слое, связанные с различными условиями на ограничивающих плоскостях. Г. А. Остроумов теоретически и экспериментально исследовал условия возникновения конвекции в вертикальном круговом канале [4]. По этим проблемам в [5] представлен обширный обзор. Конвекция (от лат. сonvectio – принесение, доставка) — движение жидкости или газа в поле тяжести под влиянием потока теплоты, идущего снизу (иногда сверху). Движущей (подъёмной) силой является сила Архимеда F𝐴 = 𝑔 ∙ ∆𝜌 ∙ 𝑉. Разность плотностей ∆𝜌 поднимающегося объёма V и окружающей среды (рис.2) зависит от различия их температур; вещество в объёме V должно быть горячее окружающей среды. Рис. 2. Схема, поясняющая возникновение конвекции Условия образования конвекции: температура T1 в глубине конвективного слоя выше, чем на его поверхности T2 , температура поднимающегося элемента объёма выше, а плотность ниже, чем у окружающей среды. Давление 𝑝 внутри и снаружи одинаково. Подъёмная сила на 1 см3 равнаF𝐴 = 𝑔 ∙ ∆𝜌, 𝑔 - ускорение свободного падения. 5

Конвекция широко распространена в природе: она происходит в нижнем слое атмосферы Земли (тропосфере) и в атмосферах некоторых других планет, во внешних слоях Солнца толщиной 20-30% его радиуса, в центральных частях массивных звёзд, в жидком ядре планет (рис. 3). Рис. 3. Конвективные течения внутри Солнца и в жидком ядре Земли. Большое значение конвекция имеет в технологических процессах [6]. Неравномерный нагрев жидкой среды обеспечивает ее перемешивание. Экспериментально показано наличие конвективных течений в неглубоких водоемах при неравномерном нагреве сверху (рис.4). Рис. 4. Возникновение конвективных течений при нагреве сверху. В атмосфере возникновение конвективных течений связано как с неоднородностью подстилающей поверхности, так и с неравномерностью солнечного обогрева. Солнце — основной источник тепловой энергии, поступающей на Землю. Солнечные лучи достигают земной атмосферы и на верхней ее границе отдают каждому квадратному сантиметру горизонтальной поверхности в минуту около 2 калорий тепла. Эту величину 6

тепла называют солнечной постоянной. Но сквозь атмосферу до поверхности Земли доходит далеко не все тепло, мешает облачность. 2.Уравнения Навье-Стокса в приближении Бусcинеска Рассмотрим некоторые процессы, связанные с конвективным течением среды, которую будем характеризовать изменяющимися во времени и пространстве скоростью 𝑣⃗(𝑥, 𝑦, 𝑧, 𝑡), плотностью𝜌⃗(𝑥, 𝑦, 𝑧, 𝑡) и температурой 𝑇(𝑥, 𝑦, 𝑧, 𝑡). Изменение этих параметров во времени и в пространстве описывается системой уравнений тепловой конвекции в приближении ОбербекаБуссинеска [7,8]: 𝜕𝑣⃗ 𝜌 ( + (𝑣⃗∇)𝑣⃗) = −∇𝑝 + 𝜇∇2 𝑣⃗ + 𝑝𝑔⃗ 𝜕𝑡 ∂𝜌 + ∇(𝜌𝑣⃗) = 0 ∂𝑡 𝜕𝑇 + 𝛻(𝑇𝑣⃗) = 𝑘𝛻 2 𝑇 𝜕𝑡 где 𝛻 = 𝑖⃗ 𝜕 𝜕𝑥 + 𝑗⃗ 𝜕 𝜕𝑦 ⃗⃗ +𝑘 𝜕 𝜕𝑧 ⃗⃗ ) – орты – векторный оператор Гамильтона, ( 𝑖⃗, 𝑗⃗, 𝑘 прямоугольной системы координат. Здесь вектор 𝑣⃗(𝑥, 𝑦, 𝑧, 𝑡) - скорость жидкости, 𝑝(𝑥, 𝑦, 𝑧) — поле давлений, 𝜇 — коэффициент вязкости, 𝑘 — коэффициент температуропроводности, 𝜌 — плотность жидкости; 𝑔 — ускорение силы тяжести. Для замыкания системы уравнений надо написать уравнение состояния вещества 𝜌 = 𝑓(𝑝, 𝑇). Предполагаем линейную зависимость плотности от температуры. Она получается при разложении уравнения состояния среды ρ = ρ0 (1 − β(T − 𝑇̂)). Здесь β = − 1 𝜌0 𝜕𝜌 (𝜕𝑇 ) — коэффициент теплового расширения. 𝑝 7

В жидкости, которую принято называть несжимаемой, при изменении температуры плотность, тем не менее, меняется. Это обстоятельство является причиной конвективных движений. Однако при исследовании конвекции в несжимаемой жидкости переменность плотности учитывается только при вычислении силы плавучести, а в уравнениях движения и неразрывности плотность принимают постоянной. Несмотря на то, что неоднородность плотности учитывается только в уравнении движения, приближение Буссинеска достаточно хорошо отражает важнейшие особенности тепловой конвекции. Будем рассматриваться плоское течение несжимаемой жидкости. Уравнение Навье–Стокса для несжимаемой жидкости в проекциях на оси декартовых координат (𝑥, 𝑦) имеют вид [9] 𝜌( 𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑝 +𝑢 +𝑣 )= + 𝜇∆𝑢 𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜌( 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑝 +𝑢 + 𝑣 ) = −𝜌𝑔 − + 𝜇∆𝑢 𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑢 𝜕𝑣 + =0 𝜕𝑥 𝜕𝑦 где ∆= 𝜕2 𝜕2 𝜕𝑥 𝜕𝑦 2 + 2 – оператор Лапласа, (𝑢, 𝑣)– компоненты вектора скорости, 𝜌 – плотность жидкости, 𝜇 – ее вязкость. Уравнение энергии для жидкости с постоянными параметрами без источников и без диссипации имеет вид 𝜕𝑇 𝜕𝑇 𝜕𝑇 𝜕2𝑇 𝜕2𝑇 𝜌𝑐𝑝 ( + 𝑢 + 𝑣 ) = 𝜆 ( 2 + 2 ) + 𝑞(𝑥, 𝑦) 𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 где 𝑇 – температура жидкости, 𝑐𝑝 – удельная теплоемкость, λ –коэффициент теплопроводности, 𝑞(𝑥, 𝑦) – мощность источников тепла в среде. 3. Системы уравнений для функции тока При решении в двумерной постановке удобно из уравнений движения исключить давление, введя функцию тока 𝜓 по формулам 8

𝑢= 𝜕𝜓 𝜕𝜓 ,𝑣 = − . 𝜕𝑦 𝜕𝑥 Тогда уравнения запишутся в виде 𝜕∆𝜓 1 𝜕𝜓 𝜕∆𝜓 𝜕𝜓 𝜕∆𝜓 𝜕𝑇 + − ( ) = ∆∆𝜓 − 𝑅𝑎 𝜕𝑡 𝑃𝑟 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝑃𝑟 𝜕𝑇 𝜕𝜓 𝜕𝑇 𝜕𝜓 𝜕𝑇 +( − ) = ∆𝑇 𝜕𝑡 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦 Здесь – число Прандтля 𝑃𝑟 = 𝑅𝑎 = 𝑔𝛽𝐻 3 ∆𝑇 𝜇 𝑐𝑝 𝜇 𝜆 , температурное число Релея . За единицу длины принята высота слоя, для температуры — разность между максимальным и минимальным значением температуры на границе. Для стационарного случая в линейном приближении (отбросив произведения функций и производных) получим 𝜕𝑇 , ∆𝑇 = 0 𝜕𝑥 с соответствующими краевыми условиями для температуры (в каждой ∆∆𝜓 = 𝑅𝑎 модели свои) и для функции тока. Верхнюю границу слоя предположим свободной, причем деформацией границы, вызванной конвекцией, пренебрежем, а нижнюю границу будем считать твердой. Обозначив объем жидкости, протекающей в направлении оси 𝑥 за единицу времени, через 𝜓0 (𝜓0 > 0), имеем следующие граничные условия для функции тока: 𝜕2𝜓 𝜕𝜓 (𝑥, 0) = 0. 𝜓(𝑥, 𝐻) = 𝜓0 , 2 (𝑥, 𝐻) = 0, 𝜓(𝑥, 0) = 𝜕𝑦 𝜕𝑦 4. Приближения для неглубоких водоемов Эту систему уравнений можно еще упростить, если рассматривать области, поперечный размер которых значительно меньше продольного размера. 9

Положим 𝑞(𝑥, 𝑦) = 0 , т.е. источник тепла в слое отсутствует. Уравнение для температуры Будем считать, что 𝜕2 𝑇 𝜕𝑥 2 ≪ 𝜕2 𝑇 𝜕2 𝑇 𝜕𝑥 𝜕𝑦 2 + 2 𝜕2 𝑇 𝜕𝑦 2 = 0. . Это можно показать, переходя в уравнении к безразмерным переменным. Пусть 𝑇 = 𝑇0 𝑇̂, где 𝑇0 – характерное значение температуры, а 𝑇̂ - безразмерная температура. Пусть характерные размеры по оси 𝑥 − 𝐿, а по оси 𝑦 − 𝐻. Сделаем замену: 𝑥 = 𝐿𝑥̂, 𝑦 = 𝐻𝑦̂, где 𝑥̂, 𝑦̂ – безразмерные величины. В правой части получим 𝜕 2 𝑇 𝜕 2 𝑇 𝑇0 𝜕 2 𝑇̂ 𝑇0 𝜕 2 𝑇̂ 𝑇0 𝜕 2 𝑇̂ 𝐻2 𝜕 2 𝑇̂ + = + = + ( ) 𝜕𝑦 2 𝜕𝑥 2 𝐻2 𝜕𝑦̂ 2 𝐿2 𝜕𝑥̂ 2 𝐻2 𝜕𝑦̂ 2 𝐿2 𝜕𝑥̂ 2 Для узких каналов считаем, что 𝜕2 𝑇 𝜕𝑥 2 ≪ 𝜕2 𝑇 𝜕𝑦 2 𝐻 𝐿 ≪ 1, тогда 𝐻2 𝐿2 ≪ 1. Следовательно, . Окончательно уравнения в приближении Буссинеска, описывающие конвекцию в плоском стационарном случае в безразмерных переменных значительно упростятся: уравнения для температуры 𝜕2 𝑇 𝜕𝑦 2 =0 и уравнения для функции тока 𝜓 с учетом (1) 𝜕4 𝜓 𝜕𝑥 4 𝜕4 𝜓 𝜕𝑦 4 ≪ 𝜕4 𝜓 𝜕𝑦 4 𝜕𝑇 = 𝑅𝑎 𝜕𝑥 , (2) где Ra — температурное число Рэлея. Уравнения (1) и (2) с различными граничными условиями будут основой всевозможных моделей большого количества процессов в жидких средах. 10

Глава 1. Конвективные течения в атмосфере Являются актуальными теоретические подходы и экспериментальные исследования конвекции в различных геофизических процессах [10, 11]. 1.1.Ветры Земли Рис. 1.1. Экваториальная конвекция — причина ветров [12]. Причина переноса воздушных масс — конвекция, подъем теплого легкого воздуха, замещение его снизу воздухом холодным. Сильнее всего за день прогреваются тропические области, где солнечные лучи падают на Землю почти отвесно. Воздух подымается вблизи экватора, приподнимая верхнюю границу тропосферы. Высота тропосферы в тропиках около 17 км, вдвое выше ее высоты у полюсов. Но куда же деваться этому воздуху? Легко понять — на больших высотах он растекается от экватора: северный воздух — на север, а южный — на юг (рис. 1.1). Вертикальные конвекционные потоки переходят в горизонтальные, теплый воздух в верхней тропосфере отчасти охлаждается, отдавая тепловое излучение космосу, затем в средних широтах опускается вниз и устремляется обратно к экватору. Мощность солнечных лучей, попадающих на Землю, за вычетом мощности отраженного планетой света. 1017 Вт. 11

Ветры, рожденные у экватора, расходятся по всей планете, хотя и довольно сложным образом. Надо учитывать вращение Земли и отклонение ветров силами Кориолиса. Характерное время переноса кинетической энергии по порядку величины есть время, через которое воздушная масса перемещается на расстояние земного радиуса 6•105 с ~ 1 неделя. Кстати, неделя есть характерное время изменения погоды. Объяснение возникновения атмосферной циркуляции от экватора до полюсов было дано еще в 1735 году английским ученым Хэдли. В его честь тропический круговорот воздуха называется ячейкой Хэдли. Рис. 1.2. Меридиональная циркуляция в тропосфере. Широкими стрелками отмечен вихревой перенос углового момента. Направление зонального ветра у поверхности земли указано внизу [13] Имеется и еще одна, обратная по направлению, ячейка атмосферной циркуляции, ячейка Феррела. Английский школьный учитель, в 1856 году модифицировавший схему Хэдли. В ячейках Хэдли воздух поднимается вверх вблизи экватора, опускается в обе стороны у 30° северной и южной широт (рис. 1.2). Между 30 и 60° широты в обоих полушариях также имеются меридиональные атмосферные циркуляции, называемые ячейками Феррела. В направление движения обратное по сравнению с ячейками Хэдли [14]. 12 них

Рис. 1.3. Система конвективных ячеек циркуляции воздуха в атмосфере Земли, соответствующие ячейкам приповерхностные ветры (черные стрелки) и направление силы трения атмосферы о поверхность планеты (красные стрелки) [15]. 1.2. Описание морских бризов Как только начинается дневное повышение температуры, воздух над сушей прогревается быстрее, чем над морем. Восходящие над сушей столбы воздуха создают градиент давления, направленный в верхних слоях с суши на море и вызывающий возникновение воздушных потоков в сторону моря. Тогда над сушей приземное давление начинает падать и внизу устанавливается градиент давления, направленный с моря на сушу. Так возникает морской бриз. Морской бриз продолжает усиливаться до тех пор, пока градиент температуры между сушей и морем не уменьшается вследствие перемещения более холодного воздуха вглубь суши и оседания воздуха над поверхностью моря. С ослаблением инсоляции в послеполуденное время, когда солнце начинает садиться или когда облака, образующиеся в результате морского бриза, закрывают землю от солнечных лучей, морской бриз затихает, а затем и совсем исчезает. Ночью все происходит в обратном порядке и приводит к образованию берегового бриза. Береговые бризы редко достигают такой же силы, как и морские. Скорость береговых бризов обычно составляет 3-8 км/час, а морских бризов 15-30 км/час. 13

Смена берегового бриза на морской происходит незадолго до полудня, а морского на береговой — вечером. Слой воздуха, охваченный бризом, может распространяться на высоту до нескольких сот метров, а выше, как правило, отмечается движение воздуха в обратном направлении — антибриз. Антибризы вместе с бризами образуют замкнутую циркуляцию (рис. 1.4). Рис. 1.4. Схема бризов Пусть изменение температуры вдоль поверхности земли (по оси x ) происходит по законам: днем, когда T0 суши < температуры воды 𝑇1 , T(x) = T0 + T1 − T0 1 + e−kx а ночью, когда температура воды > температуры суши T0 > 𝑇1 , T(x) = T1 − T1 − T0 1 + ekx Рис.1.5. Распределение температуры днем (слева) и ночью (справа). 𝑥 < 0 — море, 𝑥 < 0, — берег. Бризы могут образовываться по берегам не только морей, но и больших озер, крупных водохранилищ, а также на некоторых больших реках, на опушке леса, на окраине города и могут проникать на сушу от береговой черты на десятки километров. Бризы особенно часты летом при ясной и тихой (антициклональной) погоде. Они наблюдаются также на Балтийском море, когда долго стоит ясная и жаркая погода. 14

Соответствующие производные 𝑑𝑇(𝑥) 𝑑𝑥 = 𝑘(𝑇1 −𝑇0 )ekx (1+ekx )2 dT(x) dx = k(T1 −T0 )e−kx (1+e−kx )2 (днем) и (ночью). Согласно приближению Буссинеска, 𝜓 𝐼𝑉 = 𝐴 = 𝑑𝑇(𝑥, 𝑦) 𝑅𝑎 𝑑𝑥 𝜓 ′′′ = 𝐴𝑦 + 𝐵 𝜓 ′′ = 𝐴 𝑦2 + 𝐵𝑦 + 𝐶 2 𝑦3 𝑦2 𝜓 = 𝐴 + 𝐵 + 𝐶𝑦 + 𝐷 6 2 ′ 𝑦4 𝑦3 𝑦2 𝜓 = 𝐴 + 𝐵 + 𝐶 + 𝐷𝑦 + 𝐸 24 6 2 с граничными условиями 𝜓|𝑦=0 = 0, 𝜓′|𝑦=0 = 0, 𝜓|𝑦=ℎ = 𝜓0 , 𝜓′′|𝑦=ℎ = 0. Тогда из первой пары условий 𝐷 = 0, 𝐸 = 0. Из второй пары: ℎ4 ℎ3 ℎ2 𝐴 +𝐵 +𝐶 = 𝜓0 24 6 2 𝑦2 { 𝐴 2 + 𝐵𝑦 + 𝐶 = 0 Решением этой системы является 3 5𝐴ℎ4 𝐵 = − 3 [𝜓0 + ] ℎ 24 3𝜓0 𝐴ℎ2 𝐶= 2 + ℎ 8 { Следовательно 𝐴 𝜓0 𝑦 3 𝑦2 4 3 2 2 (2𝑦 − 5ℎ𝑦 + 3ℎ 𝑦 ) − 𝜓= ( − 3 2) 48 2 ℎ3 ℎ Построим график функции тока в пакете Maple [16]. 15

Рис. 1.6 Распределение функции тока ψ(x, y) днем (слева) и вечером (справа). Рис. 1.7. График функции тока. Слева – днем. Справа – ночью. Для некоторых случаев постараемся найти положения центров ячеек, где жидкость покоится, т. е. 𝑢 = 𝑣 = 0, пользуясь формулами 𝑢= 𝜕𝜓 𝜕𝜓 ,𝑣 = − 𝜕𝑦 𝜕𝑥 𝜕𝜓 𝐴′ 2 =0= 𝑦 (2𝑦 2 − 5ℎ𝑦 + 3ℎ2 ) = 0 𝜕𝑥 48 2𝑦 2 − 5ℎ𝑦 + 3ℎ2 = 2𝑦(𝑦 − ℎ) − 3ℎ(𝑦 − ℎ) = 0 𝑦=0 𝑦=ℎ [ 3 𝑦= ℎ 2 𝜕𝜓 𝑦3 𝑦2 = 0 = 𝐴 + 𝐵 + 𝐶𝑦 = 0 𝜕𝑦 6 2 𝑦2 𝑦 𝑦 (𝐴 + 𝐵 + 𝐶) = 0 6 2 𝐴𝑦 2 + 3𝐵𝑦 + 6𝐶 = 0 −3𝐵 ± √9𝐵2 − 24𝐴𝐶 𝑦= 2𝐴 16

Предположив, что 𝜓0 = 0, что означает отсутствие суммарного потока среды по всему слою толщиной h и выражается в замкнутости потока. Тогда решение упростится 5 𝐵 = − 𝐴ℎ 8 { 𝐴ℎ2 𝐶= 8 15 9 ∙ 25 2 2 𝐴ℎ ± √ 𝐴 ℎ − 3𝐴2 ℎ2 15𝐴ℎ ± 𝐴ℎ√9 ∙ 25 − 3 ∙ 64 8 64 𝑦= = 2𝐴 16𝐴 𝑦=ℎ 15 ± √33 16 Из всех пяти значений y лишь одно находится внутри ячеек: 𝑦 = ℎ 15−√33 . 16 𝑑2𝑇 𝑘𝑒 𝑘𝑥 (1 + 𝑒 𝑘𝑥 )2 − 2𝑒 𝑘𝑥 (1 + 𝑒 𝑘𝑥 )𝑘𝑒 𝑘𝑥 𝐴 = 2 = 𝑘(𝑇1 − 𝑇0 ) =0 (1 + 𝑒 𝑘𝑥 )2 𝑑𝑥 ′ 1 + 𝑒 𝑘𝑥 − 2𝑒 𝑘𝑥 = 0 𝑒 𝑘𝑥 = 1 𝑥=0 Следовательно, максимальное значение функции тока находится в 𝑥0 = 0 точке { 15−√33 , (смотри рис. 1.7) 𝑦0 = ℎ 16 1.3. О природе термиков У верхней неустойчивым, и границы возникает ламинарного совокупность подслоя поток случайных становится беспорядочных завихрений. Вихревой приземный слой (рис. 1.8), может распространяться до высот около 50 м. В этом слое переносы вещества осуществляются вынужденной и свободной конвекцией. 17

Рис. 1.8 . Ламинарное и вихревое течения у земной поверхности. Нагретый у земли воздух становится легким и устремляется вверх, неся с собой тепло: происходит конвекция. Рис. 1.9: Возникновение восходящих потоков воздуха Большая группа восходящих движений воздуха, вызываемых тепловыми (термическими) причинами (рис. 1.9), носит название термики [17]. Возникновение так называемых равнинных термиков, (рис. 1.9, c) объясняется неравномерным нагреванием почвы лучами солнца. Воздух нагревается главным образом от почвы, поглощающей солнечное тепло. Различные участки земной поверхности: озеро, лес, каменистая или песчаная поверхность, луг, нагреваются и охлаждаются поразному [18]. Так, над сухой дорогой температура может быть на 5 и более градусов выше, чем над близлежащим лугом. Картина конвекции над травяным покровом, пустыней или морем утром, днем и вечером будут сильно различаться. Над сухим участком образуется слой все более и более нагревающегося воздуха. Этот слой, вследствие своей большей легкости, отрывается и поднимается вверх. Над менее нагретыми участками воздух, наоборот, опускается. 18

Рис. 1.10: Конвективные потоки воздуха. Более нагретые от земной поверхности и (менее плотные) массы воздуха поднимаются, тогда как холодные (и более плотные) — опускаются. Вечером, после прекращения нагревания почвы, направление термического потока меняется (рис. 1.9), так как сырые участки почвы сохраняют тепло дольше, чем сухие, и поэтому над сухим участком воздух теперь опускается [19]. 1.4 Интегрирование уравнения общего вида Рассмотрим возможные варианты конвективных потоков. Далее во всем параграфе будем считать, что в слое нет источников тепла. Для уравнения 𝑑2𝑇 𝑑𝑦 2 = 0 примем следующие граничные условия. Пусть заданы распределения температур или теплового потока на границах слоя 𝑇|𝑦=0 = 𝑇0 (𝑥), 𝑇|𝑦=ℎ = 𝑇ℎ (𝑥) либо 𝑇|𝑦=0 = 𝑇0 (𝑥), Проинтегрируем 𝑑𝑇 𝑑𝑦 𝑑𝑇 | 𝑑𝑦 𝑦=ℎ = 𝑞(𝑥). = 𝐶1 и 𝑇 = 𝐶1 𝑦 + 𝐶2 . Тогда 𝑇(𝑥, 𝑦) = 𝑇0 (𝑥) + 𝑞(𝑥)𝑦. Согласно приближению Буссинеска: 𝛹′′′′ = 𝐴 = 𝜕𝑇 = 𝑇′0 (𝑥) + 𝑞′(𝑥)𝑦, 𝜕𝑥 𝛹′′′ = 𝑇′0 𝑦 + 𝛹′′ = 𝑇′0 𝛹′ = 𝑇′0 𝛹 = 𝑻′𝟎 𝒚𝟒 𝟐𝟒 𝑦3 + 6 𝑦2 2 + 𝒚𝟓 𝟏𝟐𝟎 + 𝑦4 24 𝑦3 6 𝑦2 2 𝑞 ′ + 𝐵, 𝑞 ′ + 𝐵𝑦 + 𝐶, 𝑞′ + 𝐵 𝒒′ + 𝑩 19 𝒚𝟑 𝟔 𝑦2 2 + 𝐶𝑦 + 𝐷, +𝑪 𝒚𝟐 𝟐 + 𝑫𝒚 + 𝑬.

c граничными условиями 𝑦 = 0: 𝛹 = 0, 𝛹 ′ = 0 𝑦 = ℎ: 𝛹 = 𝛹0 , 𝛹 ′′ = 0 Из первой пары граничных условий следует 𝐷 = 0, 𝐸 = 0. Тогда ℎ4 ℎ5 ′ ℎ3 ℎ2 𝑇′0 + 𝑞 +𝐵 +𝐶 = 𝛹0 24 120 6 2 ℎ2 ℎ3 ′ 𝑇′0 + 𝑞 + 𝐵ℎ + 𝐶 = 0 { 2 6 Решая эту систему, получим: 3 5 3 ′ 5 4 + 𝑇′ ℎ + 𝑞ℎ ] [𝛹 0 0 ℎ3 24 40 3 𝑇′0 ℎ2 11 ′ 3 𝐶 = 2 𝛹0 + − 𝑞ℎ ℎ 8 40 𝐵=− { Следовательно 𝑦4 𝑦5 ′ 𝑦3 5 3 𝛹 = 𝑇′0 + 𝑞 − 𝐵 3 (𝛹0 + 𝑇 ′ 0 ℎ4 + 𝑞 ′ ℎ5 ) 24 120 2ℎ 24 40 𝑦2 3 𝑇 ′ 0 ℎ2 11 ′ 3 + ( 2 𝛹0 + − 𝑞ℎ )= 2 ℎ 8 40 𝑦4 5 3 1 2 2 𝑦5 3 11 ′ = 𝑇 0( − 𝑦 ℎ+ 𝑦 ℎ )+𝑞 ( − 𝑦 3 ℎ2 − 𝑦 2 ℎ3 ) 24 48 16 120 80 80 ′ 3 𝑦2 1 𝑦3 + 𝛹0 ( 2 − ) 2ℎ 2 ℎ3 1.5 Симметричный случай В данном пункте предположим симметричное распределение как нагрева так и на верхней границе, так и снизу. Предположим, что 𝑞(𝑥) = 𝐸 − 𝑑 ∙ 𝑐𝑜𝑠(𝑥) 2 𝑇0 (𝑥) = 𝑇0 + 𝑎 ∙ 𝑒 −𝛼(𝑥−𝑥1) + 𝑏 ∙ 𝑒 −𝛽(𝑥−𝑥2) 20 2

Рис. 1.11: Распределение температур в приземном слое атмосферы. На верхней границе слоя – слева поток тепла. Справа – распределение температуры на нижней границе слоя. Тогда 2 𝑇(𝑥, 𝑦) = 𝑇0 + 𝑎 ∙ 𝑒 −𝛼(𝑥−𝑥1) + 𝑏 ∙ 𝑒 −𝛽(𝑥−𝑥2) 2 2 2 (𝐸 − 𝑑 ∙ 𝑐𝑜 𝑠(𝑥) − 𝑇0 − 𝑎 ∙ 𝑒 −𝛼(𝑥−𝑥1 ) − 𝑏 ∙ 𝑒 −𝛽(𝑥−𝑥2 ) )𝑦 + ℎ Рис. 1.12: Объемное распределение температуры 𝑇(𝑥,𝑦) 𝜓(𝑥, 𝑦) = 1 2 2 (−2𝑎𝛼(𝑥 − 𝑥1 )𝑒 −𝛼(𝑥−𝑥1 ) − 2𝑏𝛽(𝑥 − 𝑥2 )𝑒 −𝛽(𝑥−𝑥2) )𝑦 4 24 1 1 2 + ((𝑑 ∙ 𝑠𝑖𝑛(𝑥) + 2𝑎𝛼(𝑥 − 𝑥1 )𝑒 −𝛼(𝑥−𝑥1) 120 ℎ 1 1 2 + 2𝑏𝛽(𝑥 − 𝑥2 )𝑒 −𝛽(𝑥−𝑥2 ) )𝑦 5 ) + 𝐵𝑦 3 + 𝐶𝑦 2 6 2 21

Рис. 1.13: Объемный график для функции тока 𝜓(𝑥,𝑦). Рис. 1.14: Симметричные конвективные ячейки для различных значений функции тока. 1.6 Несимметричный случай потока тепла сверху Изменим картину следующим образом: пусть нагрев от земли остается симметричным, а поток на верхней границе - неравномерый. 1 𝑞(𝑥) = 𝐸 − 𝑑 ∙ 𝑐𝑜 𝑠(𝑥) − 𝑒 ∙ 𝑠𝑖𝑛 ( 𝑥) 2 2 𝑇0 (𝑥) = 𝑇0 + 𝑎 ∙ 𝑒 −𝛼(𝑥−𝑥1) + 𝑏 ∙ 𝑒 −𝛽(𝑥−𝑥2) 22 2

Рис. 1.15: Распределение температур в приземном слое атмосферы. На верхней границе слоя – слева поток тепла. Справа – распределение температуры на нижней границе слоя. Тогда 2 𝑇(𝑥, 𝑦) = 𝑇0 + 𝑎 ∙ 𝑒 −𝛼(𝑥−𝑥1) + 𝑏 ∙ 𝑒 −𝛽(𝑥−𝑥2) + 2 2 2 1 (𝐸 − 𝑑 ∙ 𝑐𝑜 𝑠(𝑥) − 𝑒 ∙ 𝑠𝑖 𝑛 (2 𝑥) − 𝑇0 − 𝑎 ∙ 𝑒 −𝛼(𝑥−𝑥1) − 𝑏 ∙ 𝑒 −𝛽(𝑥−𝑥2) ) 𝑦 ℎ Рис. 1.16: Объемное распределение температуры 𝑇(𝑥,𝑦). ψ(x, y) = + 1 2 2 (−2𝑎𝛼(𝑥 − 𝑥1 )𝑒 −𝛼(𝑥−𝑥1) − 2𝑏𝛽(𝑥 − 𝑥2 )𝑒 −𝛽(𝑥−𝑥2) )𝑦 4 24 1 1 1 1 ((2𝑑 ∙ sin(𝑥) cos(𝑥) − 𝑒 ∙ cos ( 𝑥) 120 ℎ 2 2 1 2 2 + 2𝑎𝛼(𝑥 − 𝑥1 )𝑒 −𝛼(𝑥−𝑥1 ) + 2𝑏𝛽(𝑥 − 𝑥2 )𝑒 −𝛽(𝑥−𝑥2) ) 𝑦 5 ) + 𝐵𝑦 3 6 1 + 𝐶𝑦 2 2 23

Рис. 1.17: Объемный график для функции тока 𝜓(𝑥,𝑦). Рис. 1.18: Конвективные ячейки для различных значений функции тока. 1.7 Несимметричный случай температуры снизу В данном пункте предположим равномерный поток на верхней границе и несимметричный поток тепла от поверхности. 𝑞(𝑥) = 𝐸 − 𝑑 ∙ 𝑐𝑜 𝑠(𝑥) 2 𝑇0 (𝑥) = 𝑇0 + 𝑎 ∙ 𝑒 −𝛼(𝑥−𝑥1) + 𝑏 ∙ 𝑒 −𝛽(𝑥−𝑥2) 𝑎≪𝑏 24 2

Рис. 1.19: Распределение температур в приземном слое атмосферы. На верхней границе слоя – слева поток тепла. Справа – распределение температуры на нижней границе слоя. 2 2 𝑇(𝑥, 𝑦) = 𝑇0 + 𝑎 ∙ 𝑒 −𝛼(𝑥−𝑥1) + 𝑏 ∙ 𝑒 −𝛽(𝑥−𝑥2 ) + 2 2 1 (𝐸 − 𝑑 ∙ 𝑐𝑜 𝑠(𝑥) − 𝑒 ∙ 𝑠𝑖 𝑛 (2 𝑥) − 𝑇0 − 𝑎 ∙ 𝑒 −𝛼(𝑥−𝑥1) − 𝑏 ∙ 𝑒 −𝛽(𝑥−𝑥2) ) 𝑦 ℎ Рис. 1.20: Объемное распределение температуры 𝑇(𝑥,𝑦). 𝜓(𝑥, 𝑦) = 1 2 2 (−2𝑎𝛼(𝑥 − 𝑥1 )𝑒 −𝛼(𝑥−𝑥1 ) − 2𝑏𝛽(𝑥 − 𝑥2 )𝑒 −𝛽(𝑥−𝑥2) )𝑦 4 24 1 1 2 + ((𝑑 ∙ 𝑠𝑖𝑛(𝑥) + 2𝑎𝛼(𝑥 − 𝑥1 )𝑒 −𝛼(𝑥−𝑥1) 120 ℎ 1 1 2 + 2𝑏𝛽(𝑥 − 𝑥2 )𝑒 −𝛽(𝑥−𝑥2 ) )𝑦 5 ) + 𝐵𝑦 3 + 𝐶𝑦 2 6 2 25

Рис. 1.21: Объемный график для функции тока 𝜓(𝑥,𝑦). Рис. 1.22: Конвективные ячейки для различных значений функции тока. 1.8 Выводы по первой главе. На основе упрощенных уравнений Буссинеска были построены модели конвективных течений в атмосфере. Рассмотрены различные варианты краевых условий. Видим, что все модели достаточно адекватно описывают такие процессы в приземном слое атмосферы как термики и бризы. 26

Глава 2 Слой воды с источниками тепла 2.1 Тепловое загрязнение Тепловое загрязнение водоемов и прибрежных морских акваторий возникает в результате сброса нагретых сточных вод электростанциями и некоторыми промышленными производствами. Сброс нагретых вод во многих случаях обуславливает повышение температуры воды в водоемах на 6-8 градусов Цельсия. Повышение температуры воды усиливает скорость обмена веществ у рыб и водных беспозвоночных. В свою очередь это повышает их потребность в кислороде. В то же самое в результате повышения температуры воды содержание в ней кислорода падает, тогда как потребность в нем живых организмов возрастает. Возросшая потребность в кислороде, его нехватка вызывают жестокий физиологический стресс и даже смерть. В летнее время повышение температуры воды всего на несколько градусов может вызвать 100%-ную гибель рыб и беспозвоночных, особенно тех, которые обитают у южных границ температурного интервала. Повышение температуры воды способно нарушить структуру растительного мира водоемов. Характерные для холодной воды водоросли заменяются более теплолюбивыми и, наконец, при высоких температурах полностью ими вытесняются. Все ᴨперечисленные выше последствия теплового загрязнения водоемов наносят огромный вред природным экосистемам и приводят к пагубному изменению среды обитания человека [20-22]. Рассмотрим плоское течение в бесконечном горизонтальном слое, имеющем разные температуры на границах: 𝑇|𝑦=0 = 𝑇0 и 𝑇|𝑦=ℎ = 𝑇ℎ . Если рассматривать только процесс теплопроводности и только в поперечном направлении (вдоль оси 𝑦), и учесть наличие локальных источников тепла в 27

слое, то уравнение энергии имеет очень простой вид 𝑑2𝑇 𝑑𝑦 2 = −𝑞(𝑥, 𝑦) . Его решение будет зависеть от вида функции 𝑞(𝑥, 𝑦). Пусть в слое жидкости имеются области с различными 𝑞(𝑥, 𝑦): в одних областях не будет источников тепла, в других – надо их учитывать. Уравнения в приближении Буссинеска‚ описывающие конвекцию, в плоском стационарном случае в безразмерных переменных имеют вид 𝜕4𝜓 𝜕𝜓 = 𝑅𝑎 𝜕𝑦 4 𝜕𝑥 𝑑2𝑇 = −𝑞(𝑥, 𝑦) 𝑑𝑦 2 где 𝜓 — функция тока; 𝑇 — температура; 𝑅𝑎 — тепловое число Рэлея. За единицу длины принята высота слоя ℎ , температуры — разность между максимальным и минимальным значением температуры на границе, функции тока — ν (значение кинематического коэффициента вязкости). 2.2. Интегрирование уравнения общего вида Сначала рассмотрим уравнение в неглубоком водоеме вида 𝑑2𝑇 = −𝑞(𝑥, 𝑦) 𝑑𝑦 2 с граничными условиями T|𝑦=0 = 𝑇0 (𝑥), 𝑇|𝑦=ℎ = 𝑇ℎ (𝑥). (2.1) Тогда T ′ 𝑦 = −q(x)y + 𝐶1 T = −q(x) 𝑦2 2 + 𝐶1 𝑦 + 𝐶2 , и, подставляя (2.1) в (2.2) получим 𝐶2 = 𝑇0 (𝑥) 28 (2.2)

ℎ2 𝑇ℎ (𝑥) = −q(x) + 𝐶1 ℎ + 𝑇0 (𝑥). 2 𝐶1 = ℎ2 − 𝑇0 (𝑥) 𝑇ℎ (𝑥) − 𝑇0 (𝑥) q(x) ℎ2 2 = + , ℎ ℎ ℎ 2 𝑇ℎ (𝑥) + q(x) и, следовательно, 𝑦 𝑦 ℎ2 𝑦2 T = 𝑇0 (𝑥) + (𝑇ℎ (𝑥) − 𝑇0 (𝑥)) + [ q(x) − q(x) ] ℎ ℎ 2 2 𝑦 ℎ2 𝑦 = 𝑇0 (𝑥) + [𝑇ℎ (𝑥) − 𝑇0 (𝑥) + q(x)(1 − )] ℎ 2 ℎ 𝑦 ℎ2 𝑦 𝑦2 = 𝑇0 (𝑥) + (𝑇ℎ (𝑥) − 𝑇0 (𝑥)) + q(x)( − 2 ) ℎ 2 ℎ ℎ Решим уравнение для функции тока ψ𝑦 IV = 𝜕𝑇(𝑥,𝑦) 𝜕𝑥 , т.е. IV 𝜕𝑇(𝑥, 𝑦) 𝑦 ℎ2 𝑦 𝑦2 = = 𝑇′0 (𝑥) + (𝑇′ℎ (𝑥) − 𝑇′0 (𝑥)) + q′(x)( − 2 ) 𝜕𝑥 ℎ 2 ℎ ℎ ′′′ 𝑦 2 ℎ2 ′ 𝑦2 𝑦3 = 𝑇 0 (𝑥)𝑦 + (𝑇 ℎ (𝑥) − 𝑇 0 (𝑥)) + q (x) ( − 2 ) + 𝐵 2ℎ 2 2ℎ 3ℎ ψ ψ ′ ′ ′ 𝑦2 𝑦 3 ℎ2 ′ 𝑦3 𝑦4 ′ ′ ψ = 𝑇 0 (𝑥) + (𝑇 ℎ (𝑥) − 𝑇 0 (𝑥)) + q (x) ( − ) + 𝐵𝑦 + 𝐶 2 6ℎ 2 6ℎ 12ℎ2 ′′ ′ ψ′ = 𝑇 ′ 0 (𝑥) 𝑦3 𝑦4 ℎ2 𝑦4 𝑦5 𝑦2 + (𝑇 ′ ℎ (𝑥) − 𝑇 ′ 0 (𝑥)) + q′ (x) ( − + 𝐵 ) 6 24ℎ 2 24ℎ 60ℎ2 2 + 𝐶𝑦 + 𝐷 𝑦4 𝑦5 ℎ2 ′ 𝑦5 𝑦6 𝑦3 ′ ′ ψ = 𝑇 0 (𝑥) + (𝑇 ℎ (𝑥) − 𝑇 0 (𝑥)) + q (x) ( − )+𝐵 24 120ℎ 2 120ℎ 360ℎ2 6 ′ 𝑦2 + 𝐶 + 𝐷𝑦 + 𝐸. 2 Краевые условия выглядят следующим образом: ψ|𝑦=0 = 0, ψ|𝑦=ℎ = ψ0 ψ′|𝑦=0 = 0, ψ′′|𝑦=ℎ = 0 Очевидным образом D = E = 0. Из условий на границе y = h следует 29

ℎ4 ℎ5 ℎ2 ℎ5 ℎ6 ℎ3 ℎ2 + (𝑇 ′ ℎ (𝑥) − 𝑇 ′ 0 (𝑥)) + q′ (x) ( − ) + 𝐵 + 𝐶 24 120ℎ 2 120ℎ 360ℎ2 6 2 2 3 2 3 4 ℎ ℎ ℎ ℎ ℎ 0 = 𝑇 ′ 0 (𝑥) + (𝑇 ′ ℎ (𝑥) − 𝑇 ′ 0 (𝑥)) + q′ (x) ( − ) + 𝐵ℎ + 𝐶 2 6ℎ 2 6ℎ 12ℎ2 ψ0 = 𝑇 ′ 0 (𝑥) { Решая эту систему, найдем C и B в следующем виде: 3 3ℎ4 ′ 2ℎ4 ′ 13ℎ6 ′ (𝑥) (𝑥) 𝐵 = − 3 [ψ0 + 𝑇 + 𝑇 + q (x)] ℎ 40 ℎ 15 0 720 ℎ2 ′ 1 2 ′ ℎ4 ′ (𝑥) (𝑥) 𝐶 = −ℎ𝐵 − 𝑇 ℎ − ℎ 𝑇0 − q (x) { 6 3 24 И, в конечном итоге 𝑦4 𝑦5 ℎ2 ′ 𝑦5 𝑦6 𝑦3 ′ ′ ψ = 𝑇 0 (𝑥) + (𝑇 ℎ (𝑥) − 𝑇 0 (𝑥)) + q (x) ( − )+𝐵 24 120ℎ 2 120ℎ 360ℎ2 6 ′ 𝑦2 +𝐶 . 2 2.3. Трехслойная схема распределения источников тепла в объеме воды. Рис. 2.1: Распределение температур 𝑇(𝑥,𝑦) по слоям. Пусть в среднем слое существует поток тепла 2 2 𝑞01 (𝑥) = 𝑎𝑒 −𝛼(𝑥−𝑥1) + 𝑏𝑒 −𝛽(𝑥−𝑥2) с двумя источниками с центрами в (x1 , 𝑦0 +𝑦1 2 ) и (x2 , 𝑦0 +𝑦1 2 ). В верхнем и нижнем слое положим отсутствие источников, поток нулевой: 30

d2 𝑇 dy2 = 0 для [0, 𝑦0 ] и [𝑦1 , h] Тогда интегрирование для [0, y0 ] дает: 𝑇00 = 𝐶10 𝑦 + 𝐶20 при y = 0 𝑇00 = 𝑇0 → 𝑇0 = 𝐶10 0 + 𝐶20 → 𝐶20 = 𝑇0 → 𝑇00 = 𝐶10 𝑦 + 𝑇0 Интегрирование для [y0 , ℎ]: 𝑇1ℎ = 𝐶1ℎ 𝑦 + 𝐶2ℎ 𝑦 = ℎ → 𝐶1ℎ ℎ + 𝐶2ℎ = 𝑇ℎ → 𝐶2ℎ = 𝑇ℎ − 𝐶1ℎ ℎ 𝑇1ℎ = 𝐶1ℎ (𝑦 − ℎ) + 𝑇ℎ . Интегрирование для среднего слоя [y0 , y1 ] : d𝑇01 𝑦2 = −𝑞01 (𝑥)𝑦 + 𝐶3 , 𝑇01 = −𝑞01 (𝑥) + 𝐶3 𝑦 + 𝐶4 . dy 2 Приравнивая на границах слоев функции и их производные, получим систему для определения 𝐶10, 𝐶1ℎ , 𝐶3 , 𝐶4 : 𝑦02 𝐶10 𝑦0 + 𝑇0 = −𝑞01 (𝑥) + 𝐶3 𝑦0 + 𝐶4 𝑦 = 𝑦0 2 𝐶10 = −𝑞01 (𝑥)𝑦0 + 𝐶3 𝑦 = 𝑦1 𝑦12 𝐶1ℎ (𝑦1 − ℎ) + 𝑇ℎ = −𝑞01 (𝑥) + 𝐶3 𝑦1 + 𝐶4 2 𝐶1ℎ = −𝑞01 (𝑥)𝑦1 + 𝐶3 { Решая эту систему, получаем искомые константы: 1 𝑦12 − 𝑦02 𝐶3 = [𝑇ℎ − 𝑞01 (𝑥) + 𝑞01 (𝑥)𝑦1 ℎ − 𝑇0 ] ℎ 2 𝑦02 𝐶4 = 𝑇0 − 𝑞01 (𝑥) 2 𝐶10 = 𝐶3 − 𝑞01 (𝑥)𝑦0 𝐶1h = 𝐶3 − 𝑞01 (𝑥)𝑦1 { Имея все четыре константы, можно записать функцию 𝑇(𝑥,𝑦) в едином виде с помощью функции Хевисайда следующим образом: 31

T(x, y) = (𝑇0 + 𝐶10 𝑦)(η(y) − η(y − y0 )) + (𝐶4 + 𝐶3 y − 𝑞01 (𝑥) 𝑦2 ) (η(y − y0 ) − η(y − y1 )) + (𝑇ℎ + 𝐶1h (y 2 − h))(η(y − y1) − η(y − h)) Найдем далее функцию тока из условия ∂4 Ψ ∂y4 = 𝜕𝑇(𝑥,𝑦) 𝜕𝑥 на каждом участке [0, y0 ], [y0 , y1 ]и [y0 , h]: 𝑇0 + 𝐶10 𝑦 𝑦2 2 { 𝑇ℎ + 𝐶1h (y − h). T(x, y) = 𝐶4 + 𝐶3 y − 𝑞01 (𝑥) Так 𝐶10 ≠ 𝐶10 (𝑦), 𝐶1h ≠ 𝐶1h (𝑦), 𝐶3 ≠ 𝐶3 (𝑦), 𝐶4 ≠ 𝐶4 (𝑦), 𝑞01 ≠ как 𝑞01 (𝑦), то 𝐶′10 𝑦 ∂T(x, y) 𝑦2 = 𝐶′4 + 𝐶′3 y − 𝑞′01 (𝑥) ∂x 2 𝐶 ′1h (y − h). { Тогда: Ψ1′′′′ = 𝐶′10 𝑦, Ψ1′′′ = 𝐶′10 Ψ1′′ = 𝐶′10 Ψ1′ = 𝐶′10 𝚿𝟏 = 𝑪′𝟏𝟎 Ψ2′′′′ Ψ2′′′ Ψ2′′ Ψ2′ 𝒚𝟓 𝟏𝟐𝟎 𝑦4 24 𝑦3 6 𝒚𝟑 𝟔 2 + 𝐵1 , + 𝐵1 y + 𝐶1 , + 𝐵1 + 𝑩𝟏 𝑦2 𝑦2 2 + 𝐶1 𝑦 + 𝐷1 , + 𝑪𝟏 𝒚𝟐 𝟐 + 𝑫𝟏 𝒚 + 𝑬𝟏 . 𝑦2 = 𝐶′4 + 𝐶′3 y − 𝑞 01 , 2 ′ 𝑦2 𝑦3 ′ = 𝐶′4 y + 𝐶′3 − 𝑞 01 + 𝐵2 , 2 6 𝑦2 𝑦3 𝑦4 ′ = 𝐶′4 + 𝐶′3 − 𝑞 01 + 𝐵2 𝑦 + 𝐶2 , 2 6 24 𝑦3 𝑦4 𝑦5 𝑦2 ′ = 𝐶′4 + 𝐶′3 − 𝑞 01 + 𝐵2 + 𝐶2 𝑦 + 𝐷2 , 6 24 120 2 32

𝒚𝟔 𝒚𝟓 𝒚𝟒 𝒚𝟑 𝒚𝟐 ′ ′ 𝚿𝟐 = −𝒒 𝟎𝟏 +𝑪 𝟑 +𝑪 𝟒 + 𝑩𝟐 + 𝑪𝟐 + 𝑫𝟐 𝒚 + 𝑬𝟐 . 𝟕𝟐𝟎 𝟏𝟐𝟎 𝟐𝟒 𝟔 𝟐 ′ Ψ3′′′′ = 𝐶 ′1h y − 𝐶 ′1h h, Ψ3′′′ = 𝐶 ′1h Ψ3′′ 𝑦2 − 𝐶 ′1h hy + 𝐵3 , 2 𝑦3 𝑦2 ′ = 𝐶 1h − 𝐶 1h h + 𝐵3 𝑦 + 𝐺, 6 2 ′ 𝑦4 𝑦3 𝑦2 ′ = 𝐶 1h − 𝐶 1h h + 𝐵3 + 𝐺𝑦 + 𝐷3 , 24 6 2 Ψ3′ ′ 𝒚𝟓 𝒚𝟒 𝒚𝟑 𝒚𝟐 ′ 𝚿𝟑 = 𝑪 𝟏𝐡 − 𝑪 𝟏𝐡 𝐡 + 𝑩 𝟑 + 𝑮 + 𝑫𝟑 𝒚 + 𝑬𝟑 . 𝟏𝟐𝟎 𝟐𝟒 𝟔 𝟐 ′ Константы B1 , C1 , D1 , E1 , B2 , C2 , D2 , E2 , B3 , E3 , G, D3 в уравнениях выше отыщем из граничных условий Ψ1 |𝑦=0 = 0, Ψ1 ′|𝑦=0 = 0, Ψ1 |𝑦=𝑦0 = Ψ2 |𝑦=𝑦0 , Ψ2 |𝑦=𝑦1 = Ψ3 |𝑦=𝑦1 , Ψ1 ′|𝑦=𝑦0 = Ψ2 ′|𝑦=𝑦0 , Ψ2 ′|𝑦=𝑦1 = Ψ3 ′|𝑦=𝑦1 ′′ Ψ1 |𝑦=𝑦0 = Ψ2 ′′|𝑦=𝑦0 , Ψ2 ′′|𝑦=𝑦1 = Ψ3 ′′|𝑦=𝑦1 𝑦0 𝑦1 − 𝑦0 Ψ1 |𝑦=𝑦0 = Ψ0 , Ψ | = Ψ ℎ 2 𝑦=𝑦1 0 ℎ ℎ − 𝑦1 ′′ Ψ3 |𝑦=ℎ = 0. Ψ3 |𝑦=𝑦1 = Ψ0 ℎ Из первых двух условий очевидно, что D1 = 0 и E1 = 0. Составим систему для поиска недостающих переменных: 𝑦 5 𝑦0 3 0 𝐶′10 120 + 𝐵1 𝑦 6 5 𝑦 1 1 −𝑞′ 01 720 + 𝐶 ′ 3 120 + 𝑦 5 + 𝐶1 24 1 −𝑞′ 01 120 + 𝐶 ′3 −𝑞′ 01 + 𝐵1 𝑦1 4 24 𝑦1 4 24 𝑦0 𝑦0 2 2 2 𝑦 + 𝐶 ′4 + 𝐶 ′3 𝑦 + 𝐶1 𝑦0 = 𝐶′4 2 𝑦1 6 𝑦1 6 6 3 3 𝑦 5 + 𝑦 𝐶 ′ ℎ5 1h 120 6 + 𝐶′3 𝑦0 4 24 + 𝐵2 5 𝑦0 3 6 + 𝐶2 2 24 − 𝑞′ 01 4 3 𝑦1 𝑦 + 𝐶2 𝑦1 + 𝐷2 = 𝐶 ′1h 24 − 𝐶 ′1h h 61 2 𝑦 2 𝑦 3 𝑦 2 𝐶 ′ 4 21 + 𝐵2 𝑦1 + 𝐶2 = 𝐶 ′1h 61 − 𝐶 ′1h h 21 𝑦 5 𝑦 3 𝑦 2 𝑦 𝐶′10 0 + 𝐵1 0 + 𝐶1 0 = Ψ0 0 120 6 2 ℎ + 𝐵2 𝑦1 𝑦0 3 𝑦0 4 2 5 𝑦1 4 1 1 −𝑞 ′ 01 720 + 𝐶 ′ 3 120 + 𝐶 ′4 { 6 0 0 = −𝑞′ 01 720 + 𝐶 ′ 3 120 + 𝐶 ′4 𝑦0 2 2 + 𝐷2 𝑦0 + 𝐸2 𝑦0 𝑦 + 𝐵2 0 + 𝐶2 𝑦0 + 𝐷2 120 2 𝑦0 3 𝑦0 2 𝑦0 3 𝑦0 4 ′ 𝐶′10 6 + 𝐵1 𝑦0 + 𝐶1 = 𝐶′4 2 + 𝐶′3 6 − 𝑞 01 24 + 𝐵2 𝑦0 + 𝐶2 𝑦1 4 𝑦 3 𝑦 2 𝑦1 5 𝑦1 4 𝑦 3 𝑦 2 𝐶 ′ 4 24 + 𝐵2 61 + 𝐶2 21 + 𝐷2 𝑦1 + 𝐸2 = 𝐶 ′1h 120 − 𝐶 ′1h h 24 + 𝐵3 61 + 𝐺 21 𝐶′10 𝑦0 4 6 −𝐶 ′ 𝐶 ℎ4 1h h 24 ′ ℎ3 1h 6 24 + 𝐵2 + 𝐵3 −𝐶 ′ ℎ3 6 𝑦1 3 6 + 𝐶2 +𝐺 ℎ2 1h h 2 ℎ2 2 𝑦1 2 2 + 𝐷3 ℎ + 𝐸3 = Ψ0 𝑦1 2 ℎ−𝑦1 ℎ 𝑦1 −𝑦0 ℎ + 𝐷3 𝑦1 + 𝐸3 + 𝐺𝑦1 + 𝐷3 + 𝐵3 𝑦1 + 𝐺 + 𝐷2 𝑦1 + 𝐸2 = Ψ0 + 𝐵3 ℎ + 𝐺 = 0. 33 + 𝐵3 2

Решение ее выполнено в среде Maple. Можно получить общий вид функции тока при помощи функции Хевисайда: 𝛹 = 𝛹1 (𝜂(𝑦) − 𝜂(𝑦 − 𝑦0 )) + 𝛹2 (𝜂(𝑦 − 𝑦0 ) − 𝜂(𝑦 − 𝑦1 )) + +𝛹3 (𝜂(𝑦 − 𝑦1 ) − 𝜂(𝑦 − ℎ)) Подставляя же конкретные значения, можно получить различные виды конвективных ячеек, которые мы рассмотрим далее. 2.4. Случай полной симметрии Пусть нагрев на дне и на поверхности постоянен, два источника в середине слоя располагаются симметрично друг друга и одинаковы. Пусть 2 2 𝑞01 (𝑥) = 𝑎𝑒 −𝛼(𝑥−𝑥1) + 𝑏𝑒 −𝛽(𝑥−𝑥2) Рис. 2.2: Вид функция 𝑞01 (𝑥) в среднем слое. 𝑇(𝑥, 𝑦) = (𝑇0 + 𝐶10 𝑦)(𝜂(𝑦) − 𝜂(𝑦 − 𝑦0 )) + (𝐶4 + 𝐶3 𝑦 − (𝑎𝑒 −𝛼(𝑥−𝑥1 )2 + 𝑏𝑒 −𝛽(𝑥−𝑥2 )2 𝑦2 ) ) (𝜂(𝑦 − 𝑦0 ) 2 − 𝜂(𝑦 − 𝑦1 )) + (𝑇ℎ + 𝐶1ℎ (𝑦 − ℎ))(𝜂(𝑦 − 𝑦1) − 𝜂(𝑦 − ℎ)) 34

Рис. 2.3: Объемное распределение температуры 𝑇(𝑥,𝑦). Рис. 2.4: Объемный график для функции тока 𝜓(𝑥,𝑦). 𝛹 = 𝛹1 (𝜂(𝑦) − 𝜂(𝑦 − 𝑦0 )) + 𝛹2 (𝜂(𝑦 − 𝑦0 ) − 𝜂(𝑦 − 𝑦1 )) + +𝛹3 (𝜂(𝑦 − 𝑦1 ) − 𝜂(𝑦 − ℎ)) Рис. 2.5: Конвективные ячейки для различных значений функции тока. 35

Рис. 2.6 Вид конвективных ячеек в случае ψ0 = 0 (отсутствует течение в водоеме). 2.5 Случай сильной асимметрии Пусть нагрев на дне и на поверхности постоянен, два источника в середине слоя располагаются симметрично, и правый значительно теплее левого. 2 q01 (𝑥) = 𝑎𝑒 −𝛼(𝑥−𝑥1 ) + 𝑏𝑒 −𝛽(𝑥−𝑥2) 2 𝑎≪𝑏 Тогда Рис. 2.7: Вид функция 𝑞01 (𝑥) в среднем слое. 𝑇(𝑥, 𝑦) = (𝑇0 + 𝐶10 𝑦)(𝜂(𝑦) − 𝜂(𝑦 − 𝑦0 )) + (𝐶4 + 𝐶3 𝑦 − (𝑎𝑒 −𝛼(𝑥−𝑥1 )2 + 𝑏𝑒 −𝛽(𝑥−𝑥2 )2 𝑦2 ) ) (𝜂(𝑦 − 𝑦0 ) 2 − 𝜂(𝑦 − 𝑦1 )) + (𝑇ℎ + 𝐶1ℎ (𝑦 − ℎ))(𝜂(𝑦 − 𝑦1) − 𝜂(𝑦 − ℎ)) 36

Рис. 2.8: Объемное распределение температуры 𝑇(𝑥,𝑦). Рис. 2.9: Объемный график для функции тока 𝜓(𝑥,𝑦). 𝛹 = 𝛹1 (𝜂(𝑦) − 𝜂(𝑦 − 𝑦0 )) + 𝛹2 (𝜂(𝑦 − 𝑦0 ) − 𝜂(𝑦 − 𝑦1 )) + +𝛹3 (𝜂(𝑦 − 𝑦1 ) − 𝜂(𝑦 − ℎ)) Рис. 2.10: Конвективные ячейки для различных значений функции тока. 37

Рис. 2.11 Вид конвективных ячеек в случае ψ0 = 0 (отсутствует течение в водоеме). 2.6. Модель, учитывающая наличие источника тепла и стока тепла Пусть нагрев на дне и на поверхности постоянен, два источника в середине слоя располагаются симметрично друг другу, и один из источников имеет отрицательное значение(отток тепла). 2 q01 (𝑥) = 𝑎𝑒 −𝛼(𝑥−𝑥1 ) + 𝑏𝑒 −𝛽(𝑥−𝑥2) 2 𝑎 < 0, 𝑏 > 0 Тогда Рис. 2.12: Вид функция 𝑞01 (𝑥) в среднем слое. 𝑇(𝑥, 𝑦) = (𝑇0 + 𝐶10 𝑦)(𝜂(𝑦) − 𝜂(𝑦 − 𝑦0 )) + (𝐶4 + 𝐶3 𝑦 − (𝑎𝑒 −𝛼(𝑥−𝑥1 )2 + 𝑏𝑒 −𝛽(𝑥−𝑥2 )2 𝑦2 ) ) (𝜂(𝑦 − 𝑦0 ) 2 − 𝜂(𝑦 − 𝑦1 )) + (𝑇ℎ + 𝐶1ℎ (𝑦 − ℎ))(𝜂(𝑦 − 𝑦1) − 𝜂(𝑦 − ℎ)) 38

Рис. 2.13: Объемное распределение температуры 𝑇(𝑥,𝑦). Рис. 2.14: Объемный график для функции тока 𝜓(𝑥,𝑦). 𝛹 = 𝛹1 (𝜂(𝑦) − 𝜂(𝑦 − 𝑦0 )) + 𝛹2 (𝜂(𝑦 − 𝑦0 ) − 𝜂(𝑦 − 𝑦1 )) + +𝛹3 (𝜂(𝑦 − 𝑦1 ) − 𝜂(𝑦 − ℎ)) Рис. 2.15: Конвективные ячейки для различных значений функции тока. 39

Рис. 2.16: Вид конвективных ячеек в случае ψ0 = 0 (отсутствует течение в водоеме). 2.7 Выводы по второй главе Во второй главе были рассмотрены тепловые и динамические процессы, возникающие в жидкой среде в условиях т.н. теплового загрязнения. На основе упрощения уравнений Буссинеска исследованы температурные поля и функции тока, построены соответствующие графики. Видим, что все модели достаточно адекватно описывают упомянутые процессы. 40

Глава 3. Исследование конвективных течений в озере с учетом таяния льдины В этой главе рассматриваются также конвективные течения в водоемах. Неравномерность температурных полей обусловлена наличием холодных областей (льдина) и более теплых водных пространств. Прежде чем использовать уравнения для функции тока и температуры, была построена модель нарастания и таяния льдины. Этот процесс описывается так называемой задачей Стефана. Для решения уравнений теплопроводности в областях, форма которых изменяется во времени, применяются как численные методы, так и приближенные методы математической физики [23]. Рассмотрим однофазный и двухфазный подходы для определения изменяющейся границы льдины. Показано, что для изучения конвективных течений достаточно рассматривать однофазную постановку. 3.1. Расчет таяния льдины в озере Рассматривается модельная задача о таянии в круглом озере плоской льдины в форме однородного круга. Предполагается, что температуры обеих фаз (жидкой и твердой) зависят только от радиальной координаты и времени. Для упрощения задачи будем также полагать, что соответствующие свойства обеих фаз постоянны. Наличие фазовых переходов при таянии (промерзании) льда (воды) усложняет задачу. Рассчитывается изменение фронта фазового перехода в зависимости от времени, а также время, за которое лед полностью растает или полностью закроет озеро. 41

Рис. 3.1. Изменение радиуса 𝑦(𝑡) круглой льдины Будем считать, что на границе льдины 𝑟 = 𝑦(𝑡) происходит «мгновенный» переход из твердого состояния в жидкое (учет наличия промежуточной промерзшей фазы составляет особую задачу). Возмущениями жидкой фазы передвигающимся фронтом фазового перехода будем пренебрегать и считать жидкую фазу неподвижной. 3.2. Однофазная задача Стефана Пусть по центру круглого озера располагается круглая льдина. Это дает возможность рассматривать центрально-симметричную теплопроводности для льда. Рис. 3.2. Геометрия круглой льдины Изменение температуры в твердой фазе описывается уравнением 𝜌𝑠 𝑐𝑠 𝜕𝑇𝑠 𝜆𝑠 𝜕 𝜕𝑇𝑠 = (𝑟 ) 𝜕𝑡 𝑟 𝜕𝑟 𝜕𝑟 42 задачу

и условиями 𝜕𝑇𝑠 | 𝜕𝑟 𝑟=0 =0 (3.1) (условие симметрии в центре льдины), 𝑇𝑠 |𝑟=𝑦(𝑡) = 𝑇𝑚 (3.2) температура фазового перехода (таяния или затвердевания), λ𝑠 𝜕T𝑠 𝜕𝑟 − 𝑞𝑤 (𝑡) = ρ𝑠 𝐿 𝑑𝑦 (3.3) 𝑑𝑡 условие Стефана для однофазной задачи. Здесь 𝑞𝑤 (𝑡) выражает приток или отток тепла к льдине (или от нее) от окружающей среды и берега. Числовые значения параметров [24]: Плотность воды изменяется с температурой сравнительно незначительно. Поэтому в большинстве случаев в практических расчетах ее значение может быть принято постоянным: ρ = 1000 кг/м3. Температура кристаллизации (замерзания) воды при нормальном атмосферном давлении принимается равной T𝑚 = 0℃ = 273𝐾. Коэффициент теплопроводности льда принимают в среднем равным λ𝑠 = 2,09 Вт м∙К = 2,09 м∙кг сек3 ∙К . Удельная теплоемкость воды слабо зависит от температуры, поэтому в практических расчетах ее значение может быть принято постоянным, равным c𝑙 = 4200 Дж кг∙К = 4200 м2 сек2 ∙К . Удельная теплота кристаллизации воды или энергия фазового перехода 𝐿 — это количество теплоты, которое выделяется при кристаллизации 1кг воды при постоянной температуре. Для воды она равна L = 334 ∙ 103 334 ∙ 103 Дж кг = м2 сек2 Учитывая, что при t=0°С плотность льда ρ = 917 кг/м3, а удельная теплоемкость его 𝑐𝑠 = 2100 Дж кг∙К = 2100 43 м2 сек2 ∙К получаем коэффициент

температуропроводности льда при нормальных условиях a = λ 𝑐ρ = 2.24 2.12∙971 = 4.1 ∙ 10−3 м2 /ч. Важной задачей является нахождение границы таяния (замерзания) и время этих процессов. Будем использовать интегральный подход [25]. Выразим приближенно температуру фазы так, чтобы она удовлетворяла условиям (3.2), (3.3). T𝑠 (𝑟, 𝑡) = 𝑇0 + (𝑇𝑚 − 𝑇0 ) 𝑟2 𝑦2 для 0 ≤ r ≤ y(t) (3.4) Проверим (3.1): 𝜕𝑇𝑠 2𝑟(𝑇𝑚 − 𝑇0 ) = =0 | | 𝜕𝑟 𝑟=0 𝑦(𝑡)2 𝑟=0 да, условию (3.1) удовлетворяет Проверим (3.2): 𝑇𝑠 |𝑟=𝑦(𝑡) = 𝑇𝑚 да, условию (3.2) удовлетворяет. Интегрируем уравнение теплопроводности 𝜌𝑠 𝑐𝑠 𝜕𝑇𝑠 𝜕𝑡 = 𝜆𝑠 𝜕 𝑟 𝜕𝑟 (𝑟 𝜕𝑇𝑠 𝜕𝑟 ) , умножив обе части на r: 𝑦(𝑡) 𝜕T𝑠 𝜕T𝑠 𝑦(𝑡) ρ𝑠 𝑐𝑠 ∫ 𝑟𝑑𝑟 = λ𝑠 𝑟 | 𝜕𝑡 𝜕𝑟 0 0 𝜕T𝑠 Из (3.4) имеем −2 ( 𝜕𝑡 𝑑𝑦 𝑑𝑡 𝑦(𝑡)3 −2𝑟 2 = (T𝑚 − 𝑇0 ) ( ) — подставляем, получаем: 𝑑𝑦 𝑦(𝑡) 𝑑𝑡 ) ∫ 𝑟 3 𝑑𝑟 = 𝜆𝑠 𝑦(𝑡) 𝜕𝑇𝑠 | 𝑦(𝑡)3 𝜌𝑠 𝑐𝑠 𝜕𝑟 𝑦(𝑡) (𝑇𝑚 − 𝑇0 ) 0 (𝑇𝑚 −𝑇0 ) −2 ( 4 𝜆𝑠 𝑑𝑦 𝑑𝑡 𝜆𝑠 𝜕𝑇𝑠 ) 𝑦(𝑡) = 𝜌 𝑐 𝑦(𝑡) 𝜕𝑟 | 𝑠 𝑠 𝜕𝑇𝑠 | 𝜕𝑟 𝑦(𝑡) =− 𝜌𝑠 𝑐𝑠 (𝑇𝑚 −𝑇0 ) 𝑦(𝑡) 𝑑𝑦 𝑑𝑡 𝜕T𝑠 2 = 𝜆𝑠 𝜕𝑇𝑠 | 𝜕𝑟 𝑦(𝑡) 44 − q𝑤 (𝑡) = ρ𝑠 𝐿 𝑑𝑦 𝑑𝑡 | 𝜌𝑠 𝑐𝑠 𝜕𝑟 𝑦(𝑡) — поток в твердой фазе. 2 Подставляем условие Стефана λ𝑠 , следует 𝑑𝑦 𝑑𝑡 (𝑇𝑚 −𝑇0 )

𝑑𝑦 ρ𝑠 𝑐𝑠 (T𝑚 − 𝑇0 ) 𝑑𝑦 𝑑𝑡 ρ𝑠 𝐿 + q𝑤 (𝑡) = − 𝑑𝑡 2 𝑑𝑦 (ρ𝑠 𝐿 + 𝑑𝑡 𝑑𝑦 (−ρ𝑠 𝐿 + 𝑑𝑡 ρ𝑠 𝑐𝑠 (T𝑚 −𝑇0 ) 2 ρ𝑠 𝑐𝑠 (T𝑚 −𝑇0 ) 2 ) = −q𝑤 (𝑡) для таяния ) = −q𝑤 (𝑡) для замерзания Получаем дифференциальное уравнение изменения размера льдины: 𝑑𝑦 −𝑞𝑤 (𝑡) = 𝑑𝑡 𝜌𝑠 𝑐𝑠 (𝑇𝑚 − 𝑇0 ) ± 𝜌𝑠 𝐿 2 При 𝑦(𝑡) = 0 можно найти время полного исчезновения льдины. Зная функцию 𝑦(𝑡), можно рассчитать распределение температуры в разных фазах, а также время таяния льдин различных размеров. 3.3. Двухфазная задача. Двухфазная задача таяния льдины описывается системой уравнений ρ𝑠 𝑐𝑠 𝜕T𝑠 λ𝑠 𝜕 𝜕T𝑠 = (𝑟 ) , 0 < 𝑟 < 𝑦(𝑡) 𝜕𝑡 𝑟 𝜕𝑟 𝜕𝑟 ρ𝑙 𝑐𝑙 𝜕T𝑙 λ𝑙 𝜕 𝜕T𝑙 = (𝑟 ) , 𝑦(𝑡) < 𝑟 < 𝑅 𝜕𝑡 𝑟 𝜕𝑟 𝜕𝑟 (s – твердая фаза), (l – жидкая фаза) , 𝜕T𝑠 | 𝜕𝑟 𝑟=0 =0 (3.5) (условие симметрии в центре льдины), T𝑠 |𝑟=𝑦(𝑡) = T𝑙 |𝑟=𝑦(𝑡) 𝑇𝑚 (3.6) (температура фазового перехода (таяния или замерзания)). T𝑚 = 0℃ = 273𝐾 На внешней границе жидкой фазы, соприкасающейся с берегом 𝑟 = 𝑅, либо температура является заданной функцией времени 𝑇𝑙 (𝑅, 𝑡) = 𝑇𝑤 (𝑡), либо задается внешний поток тепла от берега 45

λ𝑙 λ𝑠 𝜕T𝑠 𝜕T𝑙 | 𝜕𝑟 𝑟=𝑅 | 𝜕𝑟 𝑦(𝑡) − λ𝑙 = 𝑞𝑤 (𝑡). 𝜕T𝑙 | 𝜕𝑟 𝑦(𝑡) = ρ𝑠 𝐿 (3.7) 𝑑𝑦 𝑑𝑡 (3.8) условие Стефана. Здесь 𝑦(𝑡) – граница фазового перехода, R – радиус круглого озера, s – твердая фаза (лед). l – жидкая фаза (вода). λ𝑙 = 2,21 Вт м∙К = 2,21 м∙кг сек3 ∙К – коэффициент теплопроводности жидкой фазы (вода). Выразим приближенно температуры обеих фаз, так чтобы они удовлетворяли условиям (3.6), (3.7), (3.7). Предположим, что в твердой фазе установился параболический профиль температуры, удовлетворяющий граничным условиям (3.5) (3.6): 1 𝑟2 2 𝑦2 T𝑠 (𝑟, 𝑡) = 𝑇𝑚 (1 + ) для 0 ≤ r ≤ y(t) (3.9) Проверим (3.5): ∂T𝑠 | ∂r 𝑟=0 = 𝑇𝑚 𝑟 | 𝑦(𝑡)2 𝑟=0 = 0 – да, удовлетворяет. Проверим (3.6): T𝑠 |𝑟=𝑦(𝑡) = 𝑇𝑚 – да, удовлетворяет. Для области, занятой жидкой фазой, примем линейный профиль температуры 𝑇𝑙 (𝑟, 𝑡) = 𝐴𝑟 + 𝐵 Найдем А и В. 𝑇𝑙 |𝑟=𝑦(𝑡) = 𝐴𝑦 + 𝐵 = 𝑇𝑚 , 𝐵 = 𝑇𝑚 − 𝐴𝑦 𝜆𝑙 𝜕𝑇𝑙 𝑞𝑤 (𝑡) = 𝜆𝑙 𝐴 = 𝑞𝑤 (𝑡), 𝐴 = | 𝜕𝑟 𝑟=𝑅 𝜆𝑙 Окончательно: 𝑇𝑙 (𝑟, 𝑡) = 𝑞𝑤 (𝑡) 𝜆𝑙 (𝑟 − 𝑦) + 𝑇𝑚 (3.10) Теперь профили температуры выражены через одну неизвестную функцию 𝑦(𝑡). 46

Для широкого диапазона свойств материалов характерна так называемая тепловая квазистационарность, когда скорость продвижения фронта фазового перехода мала по сравнения со скоростью протекания тепловых процессов. В этих условиях можно считать, что в обеих фазах устанавливаются некоторые стационарные профили температуры. Это приближение достаточно часть применяется при решении двухфазных задач Стефана [26, 27, 32]. Выведем уравнение для ее нахождения, следуя интегральному методу теплового баланса, который считается в настоящее время наиболее эффективным в задачах плавления и затвердевания [25, 28]. Потребуем, чтобы профили температуры в каждой фазе удовлетворяли также и осредненным по слоям уравнениям теплопроводности. Интегрируем уравнение теплопроводности 𝜌𝑠 𝑐𝑠 𝜕𝑇𝑠 𝜕𝑡 = 𝜆𝑠 𝜕 𝑟 𝜕𝑟 (𝑟 𝜕𝑇𝑠 𝜕𝑟 ) , умножив обе части на r: 𝑦(𝑡) 𝜕𝑇𝑠 𝜕𝑇𝑠 𝑦(𝑡) 𝜌𝑠 𝑐𝑠 ∫ 𝑟𝑑𝑟 = 𝜆𝑠 𝑟 | 𝜕𝑡 𝜕𝑟 0 0 Проинтегрируем уравнения теплопроводности для твердой фазы: Из (3.9) имеем 𝜕𝑇𝑠 𝜕𝑡 𝑑𝑦 𝑑𝑡 𝑦(𝑡)3 −𝑟 2 = 𝑇𝑚 ( ). Подставляем, получаем: 𝑑𝑦 𝑦(𝑡) 𝑑𝑡 ∫ 𝑟 3 𝑑𝑟 = λ𝑠 𝑦(𝑡) 𝜕𝑇𝑠 | − 𝑦(𝑡)3 ρ𝑠 𝑐𝑠 𝜕𝑟 𝑦(𝑡) T𝑚 0 − 𝑑𝑦 𝑑𝑦 T𝑚 λ 𝜕𝑇 𝑠 𝑠 𝑑𝑡 𝑦(𝑡) = 𝑑𝑡 = λ𝑠 𝜕𝑇𝑠 | 𝑦(𝑡) , следует − | 4 ρ𝑠 𝑐𝑠 𝜕𝑟 𝑦(𝑡) 4 ρ𝑠 𝑐𝑠 𝜕𝑟 𝑦(𝑡) T𝑚 Т.е. 𝜆𝑠 𝜕𝑇𝑠 | 𝜕𝑟 𝑦(𝑡) =− 47 𝜌𝑠 𝑐𝑠 𝑇𝑚 4 𝑑𝑦 𝑑𝑡 (3.11)

поток тепла для твердой фазы на границе затвердевания. Для жидкой фазы: ρl 𝑐𝑙 𝜕T𝑙 𝜕𝑡 = λ𝑙 𝜕 𝑟 𝜕𝑟 (𝑟 𝜕T𝑙 𝜕𝑟 ), умножим обе части на r и интегрируем от y(t) до R: R 𝜕T𝑙 𝜕T𝑙 𝑅 ρl 𝑐𝑙 ∫ 𝑟𝑑𝑟 = λ𝑙 𝑟 | 𝜕𝑡 𝜕𝑟 𝑦(𝑡) 𝑦(𝑡) Из (3.10) имеем 𝜕T𝑙 1 𝑑𝑞𝑤 𝑑𝑦 = ( (𝑟 − 𝑦(𝑡)) − 𝑞𝑤 ) 𝜕𝑡 λ𝑙 𝑑𝑡 𝑑𝑡 Подставляем, получаем: R ∫ 𝑦(𝑡) 1 𝑑𝑞𝑤 2 𝑑𝑦 λ𝑙 𝜕T𝑙 𝜕T𝑙 (𝑟 − 𝑦(𝑡)𝑟) − 𝑞𝑤 𝑟) 𝑑𝑟 = − 𝑦(𝑡) ( (R | | ) λ𝑙 𝑑𝑡 𝑑𝑡 ρl 𝑐𝑙 𝜕𝑟 𝑟=𝑅 𝜕𝑡 𝑟=𝑦(𝑡) Подставляем λ𝑙 𝜕T𝑙 | 𝜕𝑟 𝑟=𝑅 = 𝑞𝑤 (𝑡), получаем: 𝑅 1 𝑑𝑞𝑤 𝑟 3 𝑟2 𝑑𝑦 𝑟 2 λ𝑙 𝑅 λ𝑙 𝜕T𝑙 = 𝑞𝑤 − 𝑦(𝑡) ( ( − 𝑦(𝑡) ) − 𝑞𝑤 )| | λ𝑙 𝑑𝑡 3 2 𝑑𝑡 2 𝑦(𝑡) ρl 𝑐𝑙 λ𝑙 ρl 𝑐𝑙 𝜕𝑡 𝑟=𝑦(𝑡) 1 𝑑𝑞𝑤 𝑅 3 − 𝑦(𝑡)3 𝑅2 − 𝑦(𝑡)2 𝑑𝑦 𝑅2 − 𝑦(𝑡)2 − 𝑦(𝑡) ( ( ) − 𝑞𝑤 ) λ𝑙 𝑑𝑡 3 2 𝑑𝑡 2 = λ𝑙 𝑅 λ𝑙 𝜕T𝑙 𝑞𝑤 − 𝑦(𝑡) | ρl 𝑐𝑙 λ𝑙 ρl 𝑐𝑙 𝜕𝑡 𝑟=𝑦(𝑡) Получаем поток тепла для жидкой фазы: 𝑦(𝑡)λ𝑙 − 𝜕T𝑙 | 𝜕𝑡 𝑟=𝑦(𝑡) ρl 𝑐𝑙 𝑑𝑞𝑤 𝑅 3 −𝑦(𝑡)3 λ𝑙 ( 𝑑𝑡 ( 3 = − 𝑦(𝑡) 𝑅 2−𝑦(𝑡)2 2 𝑑𝑦 𝑅 2 −𝑦(𝑡)2 ) − 𝑞𝑤 𝑑𝑡 2 ) + 𝑅𝑞𝑤 (3.12) Подставляем потоки (3.11), (3.12) в уравнение Стефана (3.8) − 𝑑𝑦 3 3 2 2 𝑑𝑡 + ρl 𝑐𝑙 (𝑑𝑞𝑤 𝑅 − 𝑦(𝑡) − 𝑅 − 𝑦(𝑡) (𝑑𝑞𝑤 𝑦(𝑡) + 𝑞 𝑑𝑦)) − 𝑅𝑞 𝑤 𝑤 4 λ𝑙 𝑑𝑡 3 2 𝑑𝑡 𝑑𝑡 ρ𝑠 𝑐𝑠 T𝑚 = ρ𝑠 𝐿 𝑑𝑦 𝑑𝑡 48

Получаем дифференциальное уравнение изменения радиуса льдины: 𝑑𝑦 ρ𝑠 𝑐𝑠 T𝑚 ρl 𝑐𝑙 𝑅2 − 𝑦(𝑡)2 + 𝑞𝑤 ) (ρ 𝐿 + 𝑑𝑡 𝑠 4 λ𝑙 2 ρl 𝑐𝑙 𝑑𝑞𝑤 𝑅3 − 𝑦(𝑡)3 ρl 𝑐𝑙 𝑅2 − 𝑦(𝑡)2 𝑑𝑞𝑤 = − 𝑦(𝑡) − 𝑅𝑞𝑤 λ𝑙 𝑑𝑡 3 λ𝑙 2 𝑑𝑡 Начальное условие: y(0) = R для таяния или y(0) = 0 для замерзания. Затем вычисляем температуру. Видим, что для однофазной, что для двухфазной постановок уравнения изменения размеров льдины имеют одну структуру. Можно рассчитать все параметры льдины. Однако для расчета конвективных течений при наличии льдины мы опять приближаемся к схеме, описанной в первых главах. Дело в том, что время изменения размеров льдины очень велико(длятся днями), а характерное время для конвективных движений созразмерно часам. В дальнейшем будем учитывать только влияние размера льдин на конвективные течения. 3.4. Льдина на поверхности [30, 31] Рассмотрим водоем с берегами в L1 и L2 соответственно. Пусть на поверхности водоема находится льдина с границами в x1 и x2 . Температура льдины T2 = 0℃ = 273𝐾 . Температуры берегов T0𝐿1 и T0𝐿2 соответственно, температура на дне постоянна и равна T0 . Глубина озера h. Рис. 3.3 Геометрия симметричной льдины 49

Пусть температура на поверхности открытых участков воды между берегами и льдиной выражается уравнением четвертой степени: T1 (𝑥) = 𝑎(𝑥 − 𝑥1 )4 + 𝑇2 T3 (𝑥) = 𝑏(𝑥 − 𝑥2 )4 + 𝑇2 так, чтобы, T1 (L1 ) = T0𝐿1 и T2 (L2 ) = T0𝐿2 . Очевидно, что T1 (𝑥1 ) = T2 = T2 (𝑥2 ). Тогда температуру на поверхности легко выразить в едином виде при помощи функции Хэвисайда: Tℎ = T1 (η(x − L1 ) − η(x − 𝑥1 )) + T2 (η(x − 𝑥1 ) − η(x − 𝑥2 )) + T3 (η(x − 𝑥2 ) − η(x − L2 )) Температуру по всей глубине водяного слоя при отсутствии источников тепла выразим в виде T(x, y) = 𝑇0 (𝑥) + (𝑇ℎ (𝑥) − 𝑇0 (𝑥)) 𝑦 ℎ и построить ее график. Воспользуемся выкладками первой и второй главы и найдем функции тока ψ1 , ψ2 , ψ3 на каждом из участков T1 , T2 , T3 . Для этого отыщем коэффициенты B1 , 𝐶1 , B1 , 𝐶2 , B2 , 𝐶2 , B3 , 𝐶3 для каждого из них: 3 3ℎ4 ′ 2ℎ4 ′ 𝐵1 = − 3 [ψ0 + 𝑇 (𝑥) + 𝑇 (𝑥)] ℎ 40 1 15 0 ℎ2 ′ 1 𝐶1 = −ℎ𝐵1 − 𝑇 1 (𝑥) − ℎ2 𝑇 ′ 0 (𝑥) 6 3 3 3ℎ4 ′ 2ℎ4 ′ 𝐵2 = − 3 [ψ0 + 𝑇 (𝑥) + 𝑇 (𝑥)] ℎ 40 2 15 0 ℎ2 ′ 1 𝐶2 = −ℎ𝐵2 − 𝑇 2 (𝑥) − ℎ2 𝑇 ′ 0 (𝑥) 6 3 3 3ℎ4 ′ 2ℎ4 ′ 𝐵3 = − 3 [ψ0 + 𝑇 (𝑥) + 𝑇 (𝑥)] ℎ 40 3 15 0 ℎ2 ′ 1 𝐶3 = −ℎ𝐵3 − 𝑇 3 (𝑥) − ℎ2 𝑇 ′ 0 (𝑥) 6 3 Тогда 50

𝑦4 𝑦5 𝑦3 𝑦2 ′ ′ ψ1 = 𝑇 0 (𝑥) + (𝑇 1 (𝑥) − 𝑇 0 (𝑥)) + 𝐵1 + 𝐶1 24 120ℎ 6 2 ′ 𝑦4 𝑦5 𝑦3 𝑦2 ′ ′ ψ2 = 𝑇 0 (𝑥) + (𝑇 2 (𝑥) − 𝑇 0 (𝑥)) + 𝐵2 + 𝐶2 24 120ℎ 6 2 ′ 𝑦4 𝑦5 𝑦3 𝑦2 ′ ′ ψ3 = 𝑇 0 (𝑥) + (𝑇 3 (𝑥) − 𝑇 0 (𝑥)) + 𝐵3 + 𝐶3 24 120ℎ 6 2 После чего с помощью функции Хэвисайда построим поток в едином ′ виде: ψ = ψ1 (η(x − L1 ) − η(x − 𝑥1 )) + ψ2 (η(x − 𝑥1 ) − η(x − 𝑥2 )) + ψ3 (η(x − 𝑥2 ) − η(x − L2 )) Можно построить ее график, и, задавая различные его значения, получать виды конвективных ячеек. Задавая различные значения температур на обоих берегах и границы льдины, можно получать различные виды ячеек. Задав размеры льдины в зависимости от времени (см. предыдущие пункты) можно получить различные виды этих самых ячеек в зависимости от размеров льдины, которые, в свою очередь, будут зависеть от времени. 3.5. Симметричный случай Предположим симметричный нагрев от обоих берегов и льдину с центром в центре водоема. Тогда Tℎ = (𝑎(𝑥 − 𝑥1 )4 + 𝑇2 )(η(x − L1 ) − η(x − 𝑥1 )) + T2 (η(x − 𝑥1 ) − η(x − 𝑥2 )) + (𝑎(𝑥 − 𝑥1 )4 + 𝑇2 )(η(x − 𝑥2 ) − η(x − L2 )) Рис. 3.4 Температура на поверхности водоема 51

𝑇(𝑥, 𝑦) = 𝑇0 (𝑥) 1 + ((𝑎(𝑥 − 𝑥1 )4 + 𝑇2 )(𝜂(𝑥 − 𝐿1 ) − 𝜂(𝑥 − 𝑥1 )) ℎ + 𝑇2 (𝜂(𝑥 − 𝑥1 ) − 𝜂(𝑥 − 𝑥2 )) + (𝑏(𝑥 − 𝑥2 )4 + 𝑇2 )(𝜂(𝑥 − 𝑥2 ) − 𝜂(𝑥 − 𝐿2 ))𝑦) Рис. 3.5 Распределение температуры по всей глубине. 1 (𝑎(𝑥 − 𝑥1 )3 𝑦 5 ) 1 1 𝜓=( + 𝐵1 𝑦 3 + 𝐶1 𝑦 2 ) (𝜂(𝑥 − 𝐿1 ) − 𝜂(𝑥 − 𝑥1 )) 30 ℎ 6 2 1 1 + ( 𝐵2 𝑦 3 + 𝐶2 𝑦 2 ) (𝜂(𝑥 − 𝑥1 ) − 𝜂(𝑥 − 𝑥2 )) 6 2 1 (𝑏(𝑥 − 𝑥2 )3 𝑦 5 ) 1 1 +( + 𝐵3 𝑦 3 + 𝐶3 𝑦 2 ) (𝜂(𝑥 − 𝑥2 ) 30 ℎ 6 2 − 𝜂(𝑥 − 𝐿2 )) Рис. 3.6 Пространственный график функции тока. 52

Рис. 3.7 Конвективные ячейки для двух различных размеров льдины. 3.6. Случай с двумя льдинами Рассмотрим случай, когда на поверхности находится не одна льдина, но две с границами x1 , x2 и x3 , x2 соответственно. Температура на открытом участке между ними равна T2𝑚 . Нагрев от берегов и температуру льдин возьмем такими же, как и в предыдущем пункте. Пусть 𝑇2𝑚 = 𝑠𝑖𝑛(𝑐(𝑥 − 𝑥2 ))2 + 𝑇2 Выразим температуру в едином виде: 𝑇ℎ = 𝑇1 (𝜂(𝑥 − 𝐿1 ) − 𝜂(𝑥 − 𝑥1 )) + 𝑇2 (𝜂(𝑥 − 𝑥1 ) − 𝜂(𝑥 − 𝑥2 )) + 𝑇2𝑚 (𝜂(𝑥 − 𝑥2 ) − 𝜂(𝑥 − 𝑥3 )) + 𝑇2 (𝜂(𝑥 − 𝑥3 ) − 𝜂(𝑥 − 𝑥4 )) + 𝑇3 (𝜂(𝑥 − 𝑥4 ) − 𝜂(𝑥 − 𝐿2 )) График ее будет иметь тогда вид Рис. 3.8 Распределение температуры на поверхности водоема при, наличии двух льдин 53

𝑇(𝑥, 𝑦) = 𝑇0 (𝑥) 1 + ((𝑎(𝑥 − 𝑥1 )4 + 𝑇2 )(η(x − L1 ) − η(x − 𝑥1 )) ℎ + T2 (η(x − 𝑥1 ) − η(x − 𝑥2 )) + (sin(c(x − x2 ))2 + 𝑇2 )(η(x − 𝑥2 ) − η(x − 𝑥3 )) + T2 (η(x − 𝑥3 ) − η(x − 𝑥4 )) + (𝑏(𝑥 − 𝑥2 )4 + 𝑇2 )(η(x − 𝑥4 ) − η(x − L2 ))𝑦) Рис. 3.9 Распределение температуры по всей глубине. Функцию тока найдем так же, как и в предыдущем пункте: ψ = ψ1 (η(x − L1 ) − η(x − 𝑥1 )) + ψ2 (η(x − 𝑥1 ) − η(x − 𝑥2 )) + ψ2𝑚 (η(x − 𝑥2 ) − η(x − 𝑥3 )) + ψ2 (η(x − 𝑥3 ) − η(x − 𝑥4 )) + ψ3 (η(x − 𝑥4 ) − η(x − L2 )) График ее будет иметь вид: Рис. 3.10 Пространственный график функции тока. 54

Фиксируя значения функции тока, получаем сечения этого графика в виде конвективных ячеек вида Рис. 3.11 Конвективные ячейки 3.7 Выводы по третьей главе В третьей главе исследованы конвективные течения, возникающие в водоемах под льдиной. Рассмотрен также случай для двух льдин на поверхности озера. 55

4.Выводы В дипломной работе были представлены математические модели конвективных движений в различных средах на основе упрощенных уравнениях динамики сплошной среды. Они, тем не менее, достаточно адекватно описывают некоторые конвективные процессы в объеме. Очевидна актуальность теоретических подходов в исследования конвекции в различных геофизических процессах. На основе упрощенных уравнений Буссинеска построены модели конвективных течений в атмосфере, конвективных течений в водоемах при наличии теплового загрязнения и конвективных течений в водоемах при наличии льдин на поверхности. Рассмотрены различные варианты краевых условий и описаний источников тепла. Модели подобного рода могут найти применение при исследовании технологических и природных процессов. 56

Литература 1. Benard H. Les tourbillons cellulaires dans une nappe liquid// Revue generale des Sciences, pures et appliquees, 1900, .v.12, 1261; 1309. 2. Benard H. Les tourbillons cellulaires dans une nappe liquide transportant de la chaleur par convection en regime permanent// Ann. Chim. Phys., 1901, v.7, 23, 62 3. Rayleigh, On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side// Phil. Mag., 1916, v.6, 32, 529 4. Остроумов Г.А. Свободная конвекция в условиях внутренней задачи.-М.-Л.: Гостехиздат, 1952. 5. Гершуни Г.З., Жуховицкий Е.М. Конвективная неустойчивость несжимаемой жидкости.-М.: Наука, 1972.- 302 с. 6. Лыкосов В.М., Осипков Л.П., Попов В.В., Старков В.Н. Моделирование процессов производства стекла в ванной печи// Сб. Математическая теория управления техническими объектами.- Л., изд. ЛГУ, 1982, с. 154-163. 7. Boussinesq J. Theorie analitique de la chaleur. T.2.-Paris: GauthierVillars, 1903.-625 p. 8. Oberbeck A. Uber die Warmeleitung der Flussigkeiten bei der Berucksichtigung der Stromungen infolge von Temperaturdifferenzen// Annal. Phys. Chem. - 1879. - Bd. 7, №6, S. 271-292 9. Лойцянский Л.Г. Механика жидкости и газа, М.: Наука, 1973, 848 с. 10. Голицын Г.С. Теоретические и экспериментальное исследование конвекции с геофизическими приложениями.-Л.: Гидрометеоиздат, 1980.- 56 с. 11. Голицын Г.С. Введение в динамику планетных атмосфер: Гидрометеоиздат, 1973.- 104с. 12. Бялко А.В. Наша планета- Земля.-М.: Наука, с.(библиотечка «Квант», вып. 29). 13. Планета Земля. Под ред. Д.Р.Бейтса.-М.: ИЛ, 1961.-339 с. 57 1983.-208

14. Витлицкий Г.Н. Циркуляция атмосферы в тропиках.-Л.: Гидрометеоиздат, 1971. 15. Швед Г.М. Циркуляция атмосферы//Соросовский образовательный журнал, №3, 1997, с.75-81 16. Дьяконов В. Maple 10, 11, 12, 13, 14 в математических расчетах. – СПб : Питер, 2011, 800 с. 17. Шютт К. Введение в физику полета.-М.-Л.: ОНТИ, 1938 18. Скорер Р. Аэрогидродинамика окржающей среды.-М.: Мир, 1980.550 с. 19. Качинский Н.А. Почва, ее свойства и жизнь.-М.: Наука, 1975. 20. Болбас, М. М. Основы промышленной экологии / Под ред. М. М. Болбас. - М., 1993. 21. Владимиров, А. М. Охрана окружающей среды / А. М. Владимиров и др. - СПб., 2001. 22. Протасов, В. Ф. и др. Экология, здоровье и природопользование в России / Под ред. В. Ф. Протасова. - М., 1995. 23. Гринберг Г.А. Избранные вопросы математической теории электрических и магнитных явлений.-М.-Л.: изд. АН СССР, 1948.-728 с. 24. Бонд Дж., Уотсон К., Уэлч Дж. Физическая теория газовой динамики- М.:Мир, 1968.-556с. 25. Беляев Н.М., Рядно А. А. Методы нестационарной теплопроводности .-М.: Высшая школа, 1978. 26. Саломатов В.В., Горбунов А. Д. О решении задач теплопереноса с подвижными границами методом теплового квазистационарного приближения// Физика и химия обработки материалов, 1973, №5, с. 37-44. 27. Самойлович Ю.А. Решение двумерной задачи Стефана в квазистационарном приближении/ Теплофизика высоких температур, 1975, т.13,№1, с.137-145 58

28. Гудмен Т. Применение интегральных методов в нелинейных задачах нестационарного теплообмена.- В сб. Проблемы теплообмена.-М.: Атомиздат, 1967. 29. Рубинштейн Л.И. Проблема Стефана .-Рига: Звайгзне, 1967. 30.1. Котляков В.М. Снег и лед в природе Земли. М.: Наука, 1986.-155 с. 31. Пивоваров А.А. Термика Замерзающих Водоемов. М.: изд. МГУ, 1972.-140с. 32. Колоднер И. Краевая задача с искомой границей для уравнения теплопроводности приложениями к задачам фазовых превращений//сб. переводов «Механика»,М.: ИЛ, 1957, №1, с.37-66]. 59

Приложение А К рисункам: 1.5,1.6,1.7 > with(plots); > T0d := 280; T1d := 285; k := 1/5; Psi0 := 0; h := 10; T0n := 285; T1n := 280; > T1x := T0d+(T1d-T0d)/(1+exp(-k*x)); > plot(T1x, x = -20 .. 20); 1.5 > T2x := T1n-(T1n-T0n)/(1+exp(k*x)); > plot(T2x,x=-20..20) ; 1.5 > A1:=diff(T1x,x); ; > A2 := diff(T2x, x); > Psi1 := (1/48)*A1*(2*y^4-5*h*y^3+3*h^2*y^2)-(1/2)*Psi0*(y^3/h^3-3*y^2/h^2); > Psi2 := (1/48)*A2*(2*y^4-5*h*y^3+3*h^2*y^2)-(1/2)*Psi0*(y^3/h^3-3*y^2/h^2); > plot3d(Psi1, x = -20 .. 20, y = 0 .. h, numpoints = 10000); 1.6 > plot3d(Psi2, x = -20 .. 20, y = 0 .. h, numpoints = 10000); 1.6 > y0 := (1/16)*h*(15-sqrt(33)); x0 := 0; > implicitplot([Psi1 = 5, Psi1 = 8, Psi1 = 11, (x-x0)^2+(y-y0)^2 = 0.1e-1], x = -11 .. 11, y = 0 .. h, color = [green, green, green, blue], thickness = [0, 0, 0, 10], numpoints = 100000, scaling = constrained); 1.7 > implicitplot([Psi2 = -5, Psi2 = -8, Psi2 = -11, (x-x0)^2+(y-y0)^2 = 0.1e-1], x = -11 .. 11, y = 0 .. h, color = [red, red, red, blue], thickness = [0, 0, 0, 10], numpoints = 100000, scaling = constrained); 60

1.7 Приложение Б К рисункам 1.11,1.12,1.13,1.14,1.15,1.16,1.17,1.18,1.19,1.20,1.21,1.22 > with(plots); > d := 2; e := 0; a := 2; b := 2; alpha := 5; beta := 5; T0 := 280; x1 := -3; x2 := 3; h := 10; Psi0 := 10; E := 290; q := 0; > C1 := E-d*cos(x)^2*0-d*cos(x); > C2 := T0+a*exp(-alpha*(x-x1)^2)+b*exp(-beta*(x-x2)^2); > plot(C1, x = -8 .. 8); 1.11 > plot(C2, x = -8 .. 8); 1.11 > Txy := C2+(C1-C2)*y/h+(1/2)*h^2*q*(y/h-y^3/h^3)*0; > plot3d(Txy, x = -8 .. 8, y = 0 .. h, numpoints = 10000); 1.12 > B := -3*(Psi0+(3*h^4*(1/40))*(diff(C1, x))+(2/15)*h^4*(diff(C2, x))+(13*h^6*(1/720))*(diff(q, x))*0)/h^3; > C := -h*B-(1/6)*h^2*(diff(C1, x))-(1/3)*h^2*(diff(C2, x))-(1/24)*h^4*(diff(q, x))*0; > `Ψxy` := (1/24)*(diff(C2, x))*y^4+(diff(C1, x)-(diff(C2, x)))*y^5/(120*h)+(1/2)*h^2*(diff(q, x))*(y^5/(120*h)-y^6/(360*h^2))*0+(1/6)*B*y^3+(1/2)*C*y^2; > plot3d(Psixy,x=-10..10,y=0..h,numpoints=10000) ; 1.13 61

> implicitplot([`Ψxy` = 38, `Ψxy` = 50, `Ψxy` = -30, `Ψxy` = -42], x = -10 .. 10, y = 0 .. h, color = [green, green, red, red], numpoints = 100000, scaling = constrained); 1.14 > d := 2; e := 5; a := 3; b := 3; alpha := 5; beta := 5; T0 := 280; x1 := -2; x2 := 2; h := 10; Psi0 := 10; E := 290; > C1:=E-d*cos^(2)(x)-e*sin^()(1/(2)* x); > C2 := T0+a*exp(-alpha*(x-x1)^2)+b*exp(-beta*(x-x2)^2); > plot(C1, x = -8 .. 8); 1.15 > plot(C2, x = -8 .. 8); 1.15 > Txy := C2+(C1-C2)*y/h+(1/2)*h^2*q*(y/h-y^3/h^3); > plot3d(Txy, x = -8 .. 8, y = 0 .. h, numpoints = 10000); 1.16 > B := -3*(Psi0+(3*h^4*(1/40))*(diff(C1, x))+(2/15)*h^4*(diff(C2, x))+(13*h^6*(1/720))*(diff(q, x)))/h^3; > C := -h*B-(1/6)*h^2*(diff(C1, x))-(1/3)*h^2*(diff(C2, x))-(1/24)*h^4*(diff(q, x)); > `Ψxy` := (1/24)*(diff(C2, x))*y^4+(diff(C1, x)-(diff(C2, x)))*y^5/(120*h)+(1/2)*h^2*(diff(q, x))*(y^5/(120*h)-y^6/(360*h^2))+(1/6)*B*y^3+(1/2)*C*y^2; > plot3d(Psixy,x=-10..10,y=0..h,numpoints=10000) ; 1.17 62

> implicitplot([`Ψxy` = 30, `Ψxy` = 70, `Ψxy` = -38, `Ψxy` = -90], x = -10 .. 10, y = 0 .. h, color = [green, green, red, red], numpoints = 100000, scaling = constrained); 1.18 > d := 2; e := 0; a := 3; b := 15; alpha := 5; beta := 5; T0 := 280; x1 := -2; x2 := 2; h := 10; Psi0 := 10; E := 290; > C1:=E-d*cos^(2)(x)-e*sin^()(1/(2)* x); > C2 := T0+a*exp(-alpha*(x-x1)^2)+b*exp(-beta*(x-x2)^2); > plot(C1, x = -8 .. 8); 1.19 > plot(C2, x = -8 .. 8); 1.19 > Txy := C2+(C1-C2)*y/h+(1/2)*h^2*q*(y/h-y^3/h^3); > plot3d(Txy, x = -8 .. 8, y = 0 .. h, numpoints = 10000); 1.20 > B := -3*(Psi0+(3*h^4*(1/40))*(diff(C1, x))+(2/15)*h^4*(diff(C2, x))+(13*h^6*(1/720))*(diff(q, x)))/h^3; > C := -h*B-(1/6)*h^2*(diff(C1, x))-(1/3)*h^2*(diff(C2, x))-(1/24)*h^4*(diff(q, x)); > `Ψxy` := (1/24)*(diff(C2, x))*y^4+(diff(C1, x)-(diff(C2, x)))*y^5/(120*h)+(1/2)*h^2*(diff(q, x))*(y^5/(120*h)-y^6/(360*h^2))+(1/6)*B*y^3+(1/2)*C*y^2; > plot3d(Psixy,x=-10..10,y=0..h,numpoints=10000) ; 1.21 63

> implicitplot([`Ψxy` = 50, `Ψxy` = 63, `Ψxy` = 300, `Ψxy` = -50, `Ψxy` = -55, `Ψxy` = -400], x = -10 .. 10, y = 0 .. h, color = [green, green, green, red, red, red], numpoints = 100000, scaling = constrained); 1.22 Приложение В К рисункам 2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,2.10,2.11,2.12,2.13,2.14,2.15,2.16 > with(plots); > a := 3; b := 3; alpha := 1/2; beta := 1/2; x1 := -3; x2 := 3; T0 := 281; Th := 281; h := 10; y0 := (1/3)*h; y1 := (2/3)*h; Psi0 := 10; > q01 := a*exp(-alpha*(x-x1)^2)+b*exp(-beta*(x-x2)^2); > C4 := T0-(1/2)*q01*y0^2; > C3 := (Th-(1/2)*q01*y1^2+q01*y1*h-T0+(1/2)*q01*y0^2)/h; > C10 := -q01*y0+C3; > C1h := C3-q01*y1; > Txy := (C10*y+T0)*(Heaviside(y)-Heaviside(y-y0))+(C4+C3*y-(1/2)*q01*y^2)*(Heaviside(y-y0)Heaviside(y-y1))+(C1h*(y-h)+Th)*(Heaviside(y-y1)-Heaviside(y-h)); > plot3d(Txy, x = -10 .. 10, y = 0 .. 10, numpoints = 10000); 2.3 > plot(q01, x = -10 .. 10, y = 0 .. 10); 2.2 > Txyd := (diff(C10*y+T0, x))*(Heaviside(y)-Heaviside(y-y0))+(diff(C4+C3*y-(1/2)*q01*y^2, x))*(Heaviside(y-y0)-Heaviside(y-y1))+(diff(C1h*(y-h)+Th, x))*(Heaviside(y-y1)-Heaviside(y-h)); > equations:={(diff(C10,x)*y0^(5))/(120)+(B1*y0^(3))/(6)+(C1*y0^(2))/(2)=- (diff(q01,x)*y0^(6))/(720)+(diff(C3,x)*y0^(5))/(120)+(diff(C4,x)*y0^(4))/(24)+(B2*y0^(3))/(6)+(C2*y0^(2) )/(2)+D2*y0+E2,(diff(C10,x)*y0^(4))/(24)+(B1*y0^(2))/(2)+C1*y0=(diff(q01,x)*y0^(5))/(120)+(diff(C3,x)*y0^(4))/(24)+(diff(C4,x)*y0^(3))/(6)+(B2*y0^(2))/(2)+C2*y0+D2,(di ff(C10,x)*y0^(3))/(6)+B1*y0+C1=(diff(q01,x)*y0^(4))/(24)+(diff(C3,x)*y0^(3))/(6)+(diff(C4,x)*y0^(2))/(2)+B2*y0+C2,64

(diff(q01,x)*y1^(6))/(720)+(diff(C3,x)*y1^(5))/(120)+(diff(C4,x)*y1^(4))/(24)+(B2*y1^(3))/(6)+(C2*y1^(2) )/(2)+D2*y1+E2=(diff(C1h,x)*y1^(5))/(120)(diff(C1h,x)*h*y1^(4))/(24)+(B3*y1^(3))/(6)+(G*y1^(2))/(2)+D3*y1+E3,(diff(q01,x)*y1^(5))/(120)+(diff(C3,x)*y1^(4))/(24)+(diff(C4,x)*y1^(3))/(6)+(B2*y1^(2))/(2)+C2*y1+D2=(di ff(C1h,x)*y1^(4))/(24)-(diff(C1h,x)*h*y1^(3))/(6)+(B3*y1^(2))/(2)+G*y1+D3,(diff(q01,x)*y1^(4))/(24)+(diff(C3,x)*y1^(3))/(6)+(diff(C4,x)*y1^(2))/(2)+B2*y1+C2=(diff(C1h,x)*y1^(3))/( 6)(diff(C1h,x)*h*y1^(2))/(2)+B3*y1^()+G,(diff(C10,x)*y0^(5))/(120)+(B1*y0^(3))/(6)+(C1*y0^(2))/(2)=(Psi0 *y0)/(h),(diff(q01,x)*y1^(6))/(720)+(diff(C3,x)*y1^(5))/(120)+(diff(C4,x)*y1^(4))/(24)+(B2*y1^(3))/(6)+(C2*y1^(2) )/(2)+D2*y1+E2=(Psi0*(y1-y0))/(h),(diff(C1h,x)*h^(5))/(120)(diff(C1h,x)*h*h^(4))/(24)+(B3*h^(3))/(6)+(G*h^(2))/(2)+D3*h+E3=(Psi0*(hy1))/(h),(diff(C1h,x)*h^(3))/(6)-(diff(C1h,x)*h*h^(2))/(2)+B3*h+G=0}; > solutions := solve(equations, {B1, B2, B3, C1, C2, D2, D3, E2, E3, G}); > Psi1:=(diff(C10,x)*y^(5))/(120)+(eval(B1,solutions)*y^(3))/(6)+(eval(C1,solutions)*y^(2))/(2); > Psi2:=- (diff(q01,x)*y^(6))/(720)+(diff(C3,x)*y^(5))/(120)+(diff(C4,x)*y^(4))/(24)+(eval(B2,solutions)*y^(3))/(6)+( eval(C2,solutions)*y^(2))/(2)+eval(D2,solutions)*y+eval(E2,solutions); > Psi3:=(diff(C1h,x)*y^(5))/(120)- (diff(C1h,x)*h*y^(4))/(24)+(eval(B3,solutions)*y^(3))/(6)+(eval(G,solutions)*y^(2))/(2)+eval(D3,solutions )*y+eval(E3,solutions); > Psi4 := Psi1*(Heaviside(y)-Heaviside(y-y0))+Psi2*(Heaviside(y-y0)-Heaviside(y- y1))+Psi3*(Heaviside(y-y1)-Heaviside(y-h)); > plot3d(Psi4, x = -10 .. 10, y = 0 .. h, numpoints = 10000); 2.4 > implicitplot([Psi4 = 4, Psi4 = 4.5, Psi4 = 2.5, Psi4 = 2], x = -10 .. 10, y = 0 .. h, color = [green, green, red, red], numpoints = 100000, scaling = constrained); 2.5 > Psi0 = 0; 65

>implicitplot([Psi4 = .1, Psi4 = .5, Psi4 = 1, Psi4 = 4, Psi4 = -.1, Psi4 = -.5, Psi4 = -1, Psi4 = -4], x = 10 .. 10, y = 0 .. h, color = [green, green, green, green, red, red, red, red], numpoints = 100000, scaling = constrained); 2.6 > a := 1; b := 5; alpha := 1/2; beta := 1/2; x1 := -3; x2 := 3; T0 := 281; Th := 281; h := 10; y0 := (1/3)*h; y1 := (2/3)*h; Psi0 := 10; > q01 := a*exp(-alpha*(x-x1)^2)+b*exp(-beta*(x-x2)^2); > C4 := T0-(1/2)*q01*y0^2; > C3 := (Th-(1/2)*q01*y1^2+q01*y1*h-T0+(1/2)*q01*y0^2)/h; > C10 := -q01*y0+C3; > C1h := C3-q01*y1; > Txy := (C10*y+T0)*(Heaviside(y)-Heaviside(y-y0))+(C4+C3*y-(1/2)*q01*y^2)*(Heaviside(y-y0)Heaviside(y-y1))+(C1h*(y-h)+Th)*(Heaviside(y-y1)-Heaviside(y-h)); > plot3d(Txy, x = -10 .. 10, y = 0 .. 10, numpoints = 10000); 2.8 > plot(q01, x = -10 .. 10, y = 0 .. 10); 2.7 > Txyd := (diff(C10*y+T0, x))*(Heaviside(y)-Heaviside(y-y0))+(diff(C4+C3*y-(1/2)*q01*y^2, x))*(Heaviside(y-y0)-Heaviside(y-y1))+(diff(C1h*(y-h)+Th, x))*(Heaviside(y-y1)-Heaviside(y-h)); > equations:={(diff(C10,x)*y0^(5))/(120)+(B1*y0^(3))/(6)+(C1*y0^(2))/(2)=- (diff(q01,x)*y0^(6))/(720)+(diff(C3,x)*y0^(5))/(120)+(diff(C4,x)*y0^(4))/(24)+(B2*y0^(3))/(6)+(C2*y0^(2) )/(2)+D2*y0+E2,(diff(C10,x)*y0^(4))/(24)+(B1*y0^(2))/(2)+C1*y0=(diff(q01,x)*y0^(5))/(120)+(diff(C3,x)*y0^(4))/(24)+(diff(C4,x)*y0^(3))/(6)+(B2*y0^(2))/(2)+C2*y0+D2,(di ff(C10,x)*y0^(3))/(6)+B1*y0+C1=(diff(q01,x)*y0^(4))/(24)+(diff(C3,x)*y0^(3))/(6)+(diff(C4,x)*y0^(2))/(2)+B2*y0+C2,(diff(q01,x)*y1^(6))/(720)+(diff(C3,x)*y1^(5))/(120)+(diff(C4,x)*y1^(4))/(24)+(B2*y1^(3))/(6)+(C2*y1^(2) )/(2)+D2*y1+E2=(diff(C1h,x)*y1^(5))/(120)(diff(C1h,x)*h*y1^(4))/(24)+(B3*y1^(3))/(6)+(G*y1^(2))/(2)+D3*y1+E3,66

(diff(q01,x)*y1^(5))/(120)+(diff(C3,x)*y1^(4))/(24)+(diff(C4,x)*y1^(3))/(6)+(B2*y1^(2))/(2)+C2*y1+D2=(di ff(C1h,x)*y1^(4))/(24)-(diff(C1h,x)*h*y1^(3))/(6)+(B3*y1^(2))/(2)+G*y1+D3,(diff(q01,x)*y1^(4))/(24)+(diff(C3,x)*y1^(3))/(6)+(diff(C4,x)*y1^(2))/(2)+B2*y1+C2=(diff(C1h,x)*y1^(3))/( 6)(diff(C1h,x)*h*y1^(2))/(2)+B3*y1^()+G,(diff(C10,x)*y0^(5))/(120)+(B1*y0^(3))/(6)+(C1*y0^(2))/(2)=(Psi0 *y0)/(h),(diff(q01,x)*y1^(6))/(720)+(diff(C3,x)*y1^(5))/(120)+(diff(C4,x)*y1^(4))/(24)+(B2*y1^(3))/(6)+(C2*y1^(2) )/(2)+D2*y1+E2=(Psi0*(y1-y0))/(h),(diff(C1h,x)*h^(5))/(120)(diff(C1h,x)*h*h^(4))/(24)+(B3*h^(3))/(6)+(G*h^(2))/(2)+D3*h+E3=(Psi0*(hy1))/(h),(diff(C1h,x)*h^(3))/(6)-(diff(C1h,x)*h*h^(2))/(2)+B3*h+G=0}; > solutions := solve(equations, {B1, B2, B3, C1, C2, D2, D3, E2, E3, G}); > Psi1:=(diff(C10,x)*y^(5))/(120)+(eval(B1,solutions)*y^(3))/(6)+(eval(C1,solutions)*y^(2))/(2); > Psi2:=- (diff(q01,x)*y^(6))/(720)+(diff(C3,x)*y^(5))/(120)+(diff(C4,x)*y^(4))/(24)+(eval(B2,solutions)*y^(3))/(6)+( eval(C2,solutions)*y^(2))/(2)+eval(D2,solutions)*y+eval(E2,solutions); > Psi3:=(diff(C1h,x)*y^(5))/(120)- (diff(C1h,x)*h*y^(4))/(24)+(eval(B3,solutions)*y^(3))/(6)+(eval(G,solutions)*y^(2))/(2)+eval(D3,solutions )*y+eval(E3,solutions); > > Psi4 := Psi1*(Heaviside(y)-Heaviside(y-y0))+Psi2*(Heaviside(y-y0)-Heaviside(y- y1))+Psi3*(Heaviside(y-y1)-Heaviside(y-h)); > plot3d(Psi4, x = -10 .. 10, y = 0 .. h, numpoints = 10000); 2.9 > implicitplot([Psi4 = 4, Psi4 = 5, Psi4 = 10, Psi4 = 1.8, Psi4 = 2.5, Psi4 = -5], x = -10 .. 10, y = 0 .. h, color = [green, green, green, red, red, red], numpoints = 100000, scaling = constrained); 2.10 > Psi0 = 0; >implicitplot([Psi4 = .1, Psi4 = 1, Psi4 = 5, Psi4 = -.1, Psi4 = -1, Psi4 = -5], x = -10 .. 10, y = 0 .. h, color = [green, green, green, red, red, red], numpoints = 100000, scaling = constrained) 67

2.11 > a := -1; b := 1; alpha := 1/2; beta := 1/2; x1 := -3; x2 := 3; T0 := 281; Th := 281; h := 10; y0 := (1/3)*h; y1 := (2/3)*h; Psi0 := 10; > q01 := a*exp(-alpha*(x-x1)^2)+b*exp(-beta*(x-x2)^2); > C4 := T0-(1/2)*q01*y0^2; > C3 := (Th-(1/2)*q01*y1^2+q01*y1*h-T0+(1/2)*q01*y0^2)/h; > C10 := -q01*y0+C3; > C1h := C3-q01*y1; > Txy := (C10*y+T0)*(Heaviside(y)-Heaviside(y-y0))+(C4+C3*y-(1/2)*q01*y^2)*(Heaviside(y-y0)Heaviside(y-y1))+(C1h*(y-h)+Th)*(Heaviside(y-y1)-Heaviside(y-h)); > plot3d(Txy, x = -10 .. 10, y = 0 .. 10, numpoints = 10000); 2.13 > plot(q01, x = -10 .. 10, y = -10 .. 10); 2.12 >Txyd := (diff(C10*y+T0, x))*(Heaviside(y)-Heaviside(y-y0))+(diff(C4+C3*y-(1/2)*q01*y^2, x))*(Heaviside(y-y0)-Heaviside(y-y1))+(diff(C1h*(y-h)+Th, x))*(Heaviside(y-y1)-Heaviside(y-h)); >equations:={(diff(C10,x)*y0^(5))/(120)+(B1*y0^(3))/(6)+(C1*y0^(2))/(2)=(diff(q01,x)*y0^(6))/(720)+(diff(C3,x)*y0^(5))/(120)+(diff(C4,x)*y0^(4))/(24)+(B2*y0^(3))/(6)+(C2*y0^(2) )/(2)+D2*y0+E2,(diff(C10,x)*y0^(4))/(24)+(B1*y0^(2))/(2)+C1*y0=(diff(q01,x)*y0^(5))/(120)+(diff(C3,x)*y0^(4))/(24)+(diff(C4,x)*y0^(3))/(6)+(B2*y0^(2))/(2)+C2*y0+D2,(di ff(C10,x)*y0^(3))/(6)+B1*y0+C1=(diff(q01,x)*y0^(4))/(24)+(diff(C3,x)*y0^(3))/(6)+(diff(C4,x)*y0^(2))/(2)+B2*y0+C2,(diff(q01,x)*y1^(6))/(720)+(diff(C3,x)*y1^(5))/(120)+(diff(C4,x)*y1^(4))/(24)+(B2*y1^(3))/(6)+(C2*y1^(2) )/(2)+D2*y1+E2=(diff(C1h,x)*y1^(5))/(120)(diff(C1h,x)*h*y1^(4))/(24)+(B3*y1^(3))/(6)+(G*y1^(2))/(2)+D3*y1+E3,(diff(q01,x)*y1^(5))/(120)+(diff(C3,x)*y1^(4))/(24)+(diff(C4,x)*y1^(3))/(6)+(B2*y1^(2))/(2)+C2*y1+D2=(di 68

ff(C1h,x)*y1^(4))/(24)-(diff(C1h,x)*h*y1^(3))/(6)+(B3*y1^(2))/(2)+G*y1+D3,(diff(q01,x)*y1^(4))/(24)+(diff(C3,x)*y1^(3))/(6)+(diff(C4,x)*y1^(2))/(2)+B2*y1+C2=(diff(C1h,x)*y1^(3))/( 6)(diff(C1h,x)*h*y1^(2))/(2)+B3*y1^()+G,(diff(C10,x)*y0^(5))/(120)+(B1*y0^(3))/(6)+(C1*y0^(2))/(2)=(Psi0 *y0)/(h),(diff(q01,x)*y1^(6))/(720)+(diff(C3,x)*y1^(5))/(120)+(diff(C4,x)*y1^(4))/(24)+(B2*y1^(3))/(6)+(C2*y1^(2) )/(2)+D2*y1+E2=(Psi0*(y1-y0))/(h),(diff(C1h,x)*h^(5))/(120)(diff(C1h,x)*h*h^(4))/(24)+(B3*h^(3))/(6)+(G*h^(2))/(2)+D3*h+E3=(Psi0*(hy1))/(h),(diff(C1h,x)*h^(3))/(6)-(diff(C1h,x)*h*h^(2))/(2)+B3*h+G=0}; > solutions := solve(equations, {B1, B2, B3, C1, C2, D2, D3, E2, E3, G}); > Psi1:=(diff(C10,x)*y^(5))/(120)+(eval(B1,solutions)*y^(3))/(6)+(eval(C1,solutions)*y^(2))/(2); > Psi2:=- (diff(q01,x)*y^(6))/(720)+(diff(C3,x)*y^(5))/(120)+(diff(C4,x)*y^(4))/(24)+(eval(B2,solutions)*y^(3))/(6)+( eval(C2,solutions)*y^(2))/(2)+eval(D2,solutions)*y+eval(E2,solutions); > Psi3:=(diff(C1h,x)*y^(5))/(120)- (diff(C1h,x)*h*y^(4))/(24)+(eval(B3,solutions)*y^(3))/(6)+(eval(G,solutions)*y^(2))/(2)+eval(D3,solutions )*y+eval(E3,solutions); > Psi4 := Psi1*(Heaviside(y)-Heaviside(y-y0))+Psi2*(Heaviside(y-y0)-Heaviside(y- y1))+Psi3*(Heaviside(y-y1)-Heaviside(y-h)); > plot3d(Psi4, x = -10 .. 10, y = 0 .. h, numpoints = 10000); 2.14 > implicitplot([Psi4 = 3.8, Psi4 = 4, Psi4 = 4.5, Psi4 = 3.1, Psi4 = 2.8, Psi4 = 2], x = -10 .. 10, y = 0 .. h, color = [green, green, green, red, red, red], numpoints = 100000, scaling = constrained); 2.15 > Psi0 = 0; >implicitplot([Psi4 = .1, Psi4 = 1, Psi4 = 1.5, Psi4 = -.1, Psi4 = -1, Psi4 = -1.5], x = -10 .. 10, y = 0 .. h, color = [green, green, green, red, red, red], numpoints = 100000, scaling = constrained) 2.16 69

Приложение Г К рисункам 3.6,3.10,3.11,3.12,3.13 > with(plots); > L1 := -10; L2 := 10; x1 := -3; x2 := 3; T0l1 := 280; T0l2 := 280; T2 := 273; h := 10; T0 := 277; a := (T0l1-T2)/(L1-x1)^4; b := (T0l2-T2)/(L2-x2)^4; Psi0 := 10; q := 0; > T1 := a*(x-x1)^4+T2; > T3 := b*(x-x2)^4+T2; > T := T1*(Heaviside(x-L1)-Heaviside(x-x1))+T2*(Heaviside(x-x1)-Heaviside(x- x2))+T3*(Heaviside(x-x2)-Heaviside(x-L2)); > plot(T); 3.6,3.10 > dT := (diff(T1, x))*(Heaviside(x-L1)-Heaviside(x-x1))+(diff(T2, x))*(Heaviside(x-x1)-Heaviside(xx2))+(diff(T3, x))*(Heaviside(x-x2)-Heaviside(x-L2)); > Txy := T0+(T-T0)*y/h; > plot3d(Txy, x = -10 .. 10, y = 0 .. h, numpoints = 10000); 3.11 > B1 := -3*(Psi0+(3*h^4*(1/40))*(diff(T1, x))+(2/15)*h^4*(diff(T0, x)))/h^3; > C1 := -h*B1-(1/6)*h^2*(diff(T1, x))-(1/3)*h^2*(diff(T0, x)); > B2 := -3*(Psi0+(3*h^4*(1/40))*(diff(T2, x))+(2/15)*h^4*(diff(T0, x)))/h^3; > C2 := -h*B2-(1/6)*h^2*(diff(T2, x))-(1/3)*h^2*(diff(T0, x))-(1/24)*h^4*(diff(q, x)); > B3 := -3*(Psi0+(3*h^4*(1/40))*(diff(T3, x))+(2/15)*h^4*(diff(T0, x)))/h^3; > C3 := -h*B3-(1/6)*h^2*(diff(T3, x))-(1/3)*h^2*(diff(T0, x)); > `Ψxy` := ((1/24)*(diff(T0, x))*y^4+(diff(T1, x)-(diff(T0, x)))*y^5/(120*h)+(1/6)*B1*y^3+(1/2)*C1*y^2)*(Heaviside(x-L1)-Heaviside(x-x1))+((1/24)*(diff(T0, x))*y^4+(diff(T2, x)-(diff(T0, Heaviside(x-x2))+((1/24)*(diff(T0, x)))*y^5/(120*h)+(1/6)*B2*y^3+(1/2)*C2*y^2)*(Heaviside(x-x1)x))*y^4+(diff(T3, x)))*y^5/(120*h)+(1/6)*B3*y^3+(1/2)*C3*y^2)*(Heaviside(x-x2)-Heaviside(x-L2)); > plot3d(`Ψxy`, x = -10 .. 10, y = 0 .. h, numpoints = 10000); 70 x)-(diff(T0,

3.12 > implicitplot([`Ψxy` = 10.2, `Ψxy` = 20, `Ψxy` = 50, `Ψxy` = -.4, `Ψxy` = -20, `Ψxy` = -50], x = -11 .. 11, y = 0 .. h, color = [green, green, green, red, red, red], numpoints = 100000, scaling = constrained); 3.13 >x1 := -6; x2 := 6; > implicitplot([`Ψxy` = 10.2, `Ψxy` = 20, `Ψxy` = 50, `Ψxy` = -.4, `Ψxy` = -20, `Ψxy` = -50], x = -11 .. 11, y = 0 .. h, color = [green, green, green, red, red, red], numpoints = 100000, scaling = constrained); 3.13 Приложение Д К рисунками 3.14,3.15,3.16 > with(plots); > L1 := -10; L2 := 10; x1 := -3; x2 := -2; x3 := 2; x4 := 3; T0l1 := 280; T0l2 := 280; T2 := 273; h := 10; T0 := 277; a := (T0l1-T2)/(L1-x1)^4; b := (T0l2-T2)/(L2-x4)^4; Psi0 := 10; q := 0; > c := Pi/(x3-x2); > T1 := a*(x-x1)^4+T2; > T3 := b*(x-x4)^4+T2; > T2m := sin(c*(x-x2))^2+T2; > T := T1*(Heaviside(x-L1)-Heaviside(x-x1))+T2*(Heaviside(x-x1)-Heaviside(x- x2))+T2m*(Heaviside(x-x2)-Heaviside(x-x3))+T2*(Heaviside(x-x3)-Heaviside(x-x4))+T3*(Heaviside(x-x4)Heaviside(x-L2)); > plot(T, x = L1 .. L2); 3.14 > Txy := T0+(T-T0)*y/h; 71

> plot3d(Txy, x = -10 .. 10, y = 0 .. h, numpoints = 10000); 3.15 > B1 := -3*(Psi0+(3*h^4*(1/40))*(diff(T1, x))+(2/15)*h^4*(diff(T0, x)))/h^3; > C1 := -h*B1-(1/6)*h^2*(diff(T1, x))-(1/3)*h^2*(diff(T0, x)); > B2 := -3*(Psi0+(3*h^4*(1/40))*(diff(T2, x))+(2/15)*h^4*(diff(T0, x)))/h^3; > C2 := -h*B2-(1/6)*h^2*(diff(T2, x))-(1/3)*h^2*(diff(T0, x))-(1/24)*h^4*(diff(q, x)); > B3 := -3*(Psi0+(3*h^4*(1/40))*(diff(T3, x))+(2/15)*h^4*(diff(T0, x)))/h^3; > C3 := -h*B3-(1/6)*h^2*(diff(T3, x))-(1/3)*h^2*(diff(T0, x)); > B2m := -3*(Psi0+(3*h^4*(1/40))*(diff(T2m, x))+(2/15)*h^4*(diff(T0, x)))/h^3; > C2m := -h*B2m-(1/6)*h^2*(diff(T2m, x))-(1/3)*h^2*(diff(T0, x)); > `Ψxy` := ((1/24)*(diff(T0, x))*y^4+(diff(T1, x)-(diff(T0, x)))*y^5/(120*h)+(1/6)*B1*y^3+(1/2)*C1*y^2)*(Heaviside(x-L1)-Heaviside(x-x1))+((1/24)*(diff(T0, x))*y^4+(diff(T2, x)-(diff(T0, Heaviside(x-x2))+((1/24)*(diff(T0, x)))*y^5/(120*h)+(1/6)*B2*y^3+(1/2)*C2*y^2)*(Heaviside(x-x1)x))*y^4+(diff(T2m, x)-(diff(T0, x)))*y^5/(120*h)+(1/6)*B2m*y^3+(1/2)*C2m*y^2)*(Heaviside(x-x2)-Heaviside(x-x3))+((1/24)*(diff(T0, x))*y^4+(diff(T2, x)-(diff(T0, Heaviside(x-x4))+((1/24)*(diff(T0, x)))*y^5/(120*h)+(1/6)*B2*y^3+(1/2)*C2*y^2)*(Heaviside(x-x3)x))*y^4+(diff(T3, x)-(diff(T0, x)))*y^5/(120*h)+(1/6)*B3*y^3+(1/2)*C3*y^2)*(Heaviside(x-x4)-Heaviside(x-L2)); > plot3d(`Ψxy`, x = -10 .. 10, y = 0 .. h, numpoints = 10000); > implicitplot([`Ψxy` = 10.2, `Ψxy` = 15, `Ψxy` = 50, `Ψxy` = -.4, `Ψxy` = -15, `Ψxy` = -50], x = -11 .. 11, y = 0 .. h, color = [green, green, green, red, red, red], numpoints = 100000, scaling = constrained); 3.16 72

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв