Интернет-версия справочника "Теплоэнергетика и теплотехника"

Web version of the handbook "Power Engineering and Thermal Engineering" (in construction)

 

 

Электронный курс Тепломассообмен в энергетических установках

 

А.П.Солодов

(English Version)

 Cover    Cover

А.П.Солодов, В.Ф.Очков «Mathcad / Дифференциальные модели»

Москва, Издательство МЭИ (publish@admin.mpei.ac.ru тел. (095) 361-16-81), 2002 г., – 239 с.: илл.

ISBN 5-7046-0823-x

A. Solodov, V. Ochkov «Differential Models»

In English, Springer Publishing House (Springer-Verlag, spring, 2004)

Заказ книги в издательстве - www.springer.de

Mathcad-документы книги

WebSheets (Mathcad Application Server)

Предисловие

Глава 1.      Дифференциальные математические модели

1.1.      Введение

1.2.      Законы в дифференциальной форме

1.3.      Модели роста

1.4.      Законы сохранения

1.5.      Одномерные стационарные модели. ТВЭЛ ядерного реактора.

1.6.      Заключение

Глава 2.      Интегрируемые в квадратурах дифференциальные уравнения

2.1.      Перечень интегрируемых уравнений

2.2.      Линейное уравнение первого порядка

2.3.      Линейные однородные уравнения с постоянными коэффициентами

2.4.      Линейные неоднородные уравнения с постоянными коэффициентами

2.5.      Уравнения с разделяющимися переменными

2.6.      Однородные дифференциальные уравнения

2.7.      Уравнения, допускающие понижение порядка

Глава 3.      Динамическая модель системы с тепловыделением

3.1.      Введение

3.2.      Математическая модель

3.3.      Фазовый портрет системы. Устойчивые и неустойчивые состояния равновесия

3.4.      Представление множества равновесных состояний

3.5.      Построение бифуркационной диаграммы

3.6.      Трехмерное представление равновесных состояний в форме «катастрофы сборки»

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

3.8.      Временная эволюция системы с тепловыделением

3.9.      Заключение

Глава 4.      Жесткие дифференциальные уравнения

4.1.      Введение

4.2.      Метод rkfixed. Неустойчивость решения

4.3.      Метод rkadapt

4.4.      Mетод stiffr. Решение модельной жесткой задачи

4.5.      Метод stiffr. Решение уравнений химической кинетики.

4.6.      Явные и неявные методы

4.7.      Заключение

Глава 5.      Теплопередача в окрестности критической точки поперечно обтекаемой трубы

5.1.      Введение

5.2.      Интегральное уравнение теплового пограничного слоя

5.3.      Математическое описание задачи

5.4.      Распределение скорости внешнего потока по окружности трубы

5.5.      Вывод соотношений для критической точки

5.6.      Приведение математического описания к безразмерному виду

5.7.      Представление правой части дифференциального уравнения в форме алгоритма оптимизации

5.8.      Численное интегрирование с помощью встроенной функции Odesolve

5.9.      Заключение

Глава 6.      Уравнение Фолкнера-Скэн. Трение и теплообмен в пограничном слое

6.1.      Введение

6.2.      Математическая формулировка

6.3.      Сведение краевой задачи к начальной задаче методом  sbval

6.4.      Решение начальной задачи методом rkfixed

6.5.      Построение поля течения

6.6.      Пограничный слой на проницаемой поверхности.

6.7.      Уравнение теплового пограничного слоя

6.8.      Закон теплообмена

6.9.      Неприятности при решении краевой задачи посредством функции Odesolve

6.10.        Заключение

Глава 7.      Уравнение Рэлея. Гидродинамическая неустойчивость

7.1.      Введение

7.2.      Уравнения гидродинамики для свободного сдвигового потока

7.3.      Метод малых возмущений. Линеаризация.

7.4.      Переход в комплексную область

7.5.      Численное интегрирование в комплексной области. Программа Euler

7.6.      Интегрирование и поиск собственных значений

7.7.      Возвращение в действительную область

7.8.      Заключение

Глава 8.      Кинематические волны концентрации в ионитовом фильтре

8.1.      Введение

8.2.      Уравнение сохранения концентрации

8.3.      Волновое уравнение для концентрации

8.4.      Уравнение сохранения в безразмерной форме

8.5.      Уравнение изотермы адсорбции

8.6.      Решение волнового уравнения методом характеристик (Пример 1 анимации, Пример 2 анимации)

8.7.      Заключение

Глава 9.      Кинематические ударные волны

9.1.      Введение

9.2.      Уравнение сохранения в конечно-разностной форме.

9.3.      Разрывные решения. Ударные волны

9.4.      Метод Мак-Кормака. Вычислительная программа McCrm

9.5.      Ударные волны концентрации примесей в фильтре

9.6.      Ударные волны на автостраде

9.7.      Гравитационные пузырьковые течения. Ударные волны паросодержания

9.8.      Заключение

Глава 10.    Численное моделирование температурного поля платы компьютера

10.1.        Введение

10.2.        Встроенные функции для решения уравнений в частных производных

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

10.4.        Метод Гаусса-Зайделя. Программа Plate

10.5.        Тепловая модель платы компьютера

10.6.        Задача об орбитальной платформе. Функция bvalfit

10.7.        Заключение

Глава 11.    Температурные волны

11.1.        Введение

11.2.        Математическая формулировка задачи

11.3.        Дискретное представление

11.4.        Метод прогонки. Вычислительные программы Coef и SysTRD

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

11.6.        Заключение

Литература

Приложение. Встроенные функции для решения обыкновенных дифференциальных уравнений

Предметный указатель

Предисловие

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

q              Постановка задачи и формулировка физической модели

q              Формулировка дифференциальной математической модели, т.е. модели в форме дифференциальных уравнений

q              Интегрирование дифференциальных уравнений

q              Визуализация результатов.

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

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

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

Книга состоит из одиннадцати глав и приложения.

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

 

Контрольный объем для двухмерной задачи

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

Логистическая модель с учетом случайных колебаний фактора роста

 

Глава вторая «Интегрируемые в квадратурах дифференциальные уравнения» содержит перечень и краткий рецептурный справочник для уравнений указанного типа. Показано, как можно применить символьный процессор Mathcad’а при математических выкладках в этом случае.

 

Глава третья «Динамическая модель системы с тепловыделением» посвящена системам, которые могут развиваться по катастрофическим сценариям – с воспламенением и взрывом.

 

Поверхность равновесных состояний (вид сбоку и вид сверху)

 

Эволюция системы в случае трех равновесных состояний

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

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

В главе четвертой «Жесткие дифференциальные уравнения» подробно рассмотрена проблема численного интегрирования этого особого типа уравнений. Несмотря на практическую важность вопроса, студенту или инженеру трудно найти доступное описание сущности феномена. Mathcad имеет эффективную встроенную функцию для интегрирования жестких уравнений, но его справочная система также почти ничего не содержит по обсуждаемой проблеме. Материалы данной главы восполняют в некоторой степени этот пробел.

 

Интегрирование жестких уравнений химической кинетики.

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

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

           

Идеальное обтекание цилиндра                                                                                            

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

 

     

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

 

Глава шестая «Уравнение Фолкнера-Скэн. Трение и теплообмен в пограничном слое» посвящена численному анализу фундаментальной задачи гидродинамики и тепломассообмена. Совместно с разделом 1.4 «Законы сохранения» шестая глава образует краткий теоретический курс конвективного теплообмена. С точки зрения овладения техникой работы в Mathcad’е, главной темой является численное решение краевой двухточечной задачи для системы обыкновенных дифференциальных уравнений с применением встроенной функции, позволяющей свести краевую проблему к задаче с начальными условиями.

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

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

Глава восьмая «Кинематические волны концентрации в ионитовом фильтре» развивает тему интегрирования уравнений в частных производных, открытую в первой главе задачей об ударных волнах на автостраде. Сформулирована физическая и математическая модель ионитовых фильтров – устройств, обеспечивающих получение очень чистой воды для паросиловых установок тепловых и атомных электростанций. Методом характеристик построено решение уравнения в частных производных, описывающее пространственно-временную эволюцию поля концентраций. Показано, как нелинейность дифференциальной модели приводит к появлению разрывных неоднозначных решений.

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

Глава девятая «Кинематические ударные волны» продолжает тему волновых уравнений. Разработана Mathcad-программа для численного интегрирования дифференциальных уравнений в частных производных, обеспечивающая эффективное воспроизведение ударных волн. Благодаря этому доводится до конца исследование задач об ударных волнах на автострадах и ударных волнах в фильтрах.

 

 

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

 

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

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

 

Температурное поле компьютерной платы

Температурное поле в грунте вокруг трубопровода

 

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

Температурные волны

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

Авторы стремились показать в этой книге на реальных примерах, как можно эффективно использовать Mathcad на всех этапах разработки инженерного проекта:

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

·        при получении аналитических решений на этапах, где это возможно (или где это позволяют возможности символьного процессора Mathcad)

·        при численном решении, когда аналитические решения невозможны или неэффективны

·        при презентации результатов с помощью многообразных и простых в использовании средств визуализации.

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

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

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

В конце книги приведен список литературы. Кроме прямых ссылок на источники, он включает перечень общей литературы по численному анализу и компьютерному моделированию, по-видимому, далеко не полный. Он содержит те относящиеся к нашей теме книги, которые были прочитаны авторами, показались им интересными и, следовательно, в явной или неявной форме были использованы при подготовке пособия. «Цель расчетов – не числа, а понимание» - этот девиз  книги Р.Хемминга [Ошибка! Источник ссылки не найден.] мы хотели бы сделать главной идеей и нашей книги.

Глава 11.    Температурные волны

11.1.   Введение

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

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

Например, в проблеме Управляемого Термоядерного Синтеза плотность теплового потока на тепловоспринимающих твердых поверхностях может достигать 108 Вт/м2.

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

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

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

Мы рассмотрим далее в качестве модели описанных выше процессов одномерную нестационарную задачу теплопроводности с внутренними источниками теплоты (Рис. 11.1).

 

Рис. 11.1. Одномерная нестационарная задача теплопроводности

 

Предполагается, что пространственные изменения происходят только вдоль оси координат x. Боковая поверхность считается адиабатически изолированной. Если же требуется учесть теплоообмен на боковой поверхности, то теплосъем следует имитировать внутренним стоком, точно так же как в предыдущей главе.

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

11.2.   Математическая формулировка задачи

В качестве исходной формулировки можно принять уравнение энергии, в котором следует исключить оператор конвективного переноса и учесть одномерность задачи, t = t(x, τ). В результате получится:

( 11.1)

На левом и правом торце объекта (см. Рис. 11.1) происходит тепловое взаимодействие со средой, и здесь необходимо задать граничные условия. Универсальным способом описать самые различные ситуации будет применение граничных условий третьего рода на левом (x = 0) и правом (x = L)  торцах объекта:

( 11.2)

В этих соотношениях приравнены значения плотности теплового потока,

q       поступающего из окружающей среды и вычисленного по уравнению Ньютона-Рихмана (правые части)

q       и отводимого внутрь тела посредством теплопроводности и вычисленного по закону Фурье (левые части).

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

11.3.   Дискретное представление

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

Рис. 11.2. Контрольный объем и потоки энергии

Итак, мы собираемся записать баланс энергии для выделенного контрольного объема (Рис. 11.2), сформированного в окрестности узла P. Последний находится сейчас в фокусе нашего внимания, но выведенные далее соотношения будут годиться для любого внутреннего узла.

От соседних, западного W и восточного E узлов поступают потоки теплоты за счет теплопроводности. Внутри контрольного объема выделяется теплота при действии внутреннего источника. В результате произойдет увеличение тепловой энергии в рассматриваемом объеме, что будет обнаружено по повышению температуры – от T0P до TP за время Δτ :

( 11.3)

Величина баланса bal должна быть нулевой. Большими греческими «лямбда» Λ обозначены относительные значения коэффициента теплопроводности, вычисленные для контрольных поверхностей, проходящих между узлами. Например, для восточной поверхности:

,

(11.4)

где величина λ без индекса означает некоторое характерное значение теплопроводности.

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

Уравнение ( 11.3) записано для внутренних узлов. Ниже будет получено аналогичное уравнение  для граничных узлов ( 11.8) .

Неявная схема. Подчеркнем, что температуры в других узлах, т.е. TW и TE – неизвестные величины из «будущего», такие же как TP. Поэтому ( 11.3) представляет собой уравнение с тремя неизвестными. Всего таких уравнений можно записать столько же, сколько имеется неизвестных температур в узлах сетки. Для нахождения температур необходимо решить систему уравнений.

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

Явная схема. Если в уравнении ( 11.3) разности температур во втором и третьем операторах (теплопроводности) записать для «прошлого» и снабдить обозначением «0», то получится явная схема: в каждом уравнении содержится одно единственное неизвестное значение TP (в первом операторе, описывающем изменение энергии во времени). Программа и вычисления будут очень простыми, но если временной шаг превысит некоторое значение, при счете возникнут прогрессирующие паразитные осцилляции. Ограничения на шаг довольно обременительны, поэтому в настоящее время предпочитают неявные схемы.

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

Вернемся к анализу уравнения баланса ( 11.3). Далее приводятся два варианта записи этого уравнения. Первое – в мнемонической форме, второе – в индексной, необходимой для составления программы. Неизвестные величины в этих уравнениях – температуры в узлах для следующего момента времени. Остальные, известные величины входят в коэффициенты уравнения. Среди этих известных величин – значение температуры в рассматриваемом узле T0P, взятое с предыдущего шага по времени, либо из начальных условий.

( 11.5)

Сравнивая выражения ( 11.3) и ( 11.5), получают значения коэффициентов A, B, C, D. Эти операции можно выполнить в Mathcad’е с помощью оператора collect, как показано ниже. Предварительно вводится обозначение для безразмерного шага по времени, т.н. сеточного числа Фурье:

 

 

( 11.6)

Структура получившейся линейной системы уравнений с коэффициентами A, B, C, D оказывается трехдиагональной, как в демонстрационном примере для сетки из пяти узлов:

 

( 11.7)

 

При большом числе узлов матрица коэффициентов окажется почти пустой. Например, для ста улов всего три процента ячеек будет занято числами, остальные останутся нулевыми. Ясно, что при использовании классического Гауссова исключения в основном будут перемножаться и складываться нули. Однако существует специальный вариант метода исключения, учитывающий трехдиагональную структуру матрицы и называемый методом прогонки [Ошибка! Источник ссылки не найден., Ошибка! Источник ссылки не найден., Ошибка! Источник ссылки не найден., Ошибка! Источник ссылки не найден., Ошибка! Источник ссылки не найден.]. Mathcad-программа этого метода приведена на Рис. 11.5 .

 

 

Рис. 11.3. Баланс энергии для поверхностного узла

 

Для граничных узлов контрольные объемы оказываются половинного размера, как показано на Рис. 11.3. Кроме того, для контрольной грани, попадающей на поверхность объекта, необходимо специальным образом записать тепловой поток:

( 11.8)

Это уравнение является дискретным аналогом граничных условий, заданных выше соотношениями ( 11.2). Индекс «inner» означает ближайший внутренний узел: (n-1) правого торца, 2 – для левого торца. Мы не будем повторять вычисления для коэффициентов в уравнениях для поверхностных узлов. Они полностью аналогичны тем, что выполнены для внутренних узлов. Окончательные выражения можно прочитать в тексте программы Coef (Рис. 11.4).

11.4.   Метод прогонки. Вычислительные программы Coef и SysTRD

Mathcad-программа Coef

Эта подпрограмма (Рис. 11.4) вычисляет коэффициенты A, B, C, D для системы уравнений, структура которой иллюстрируется соотношением ( 11.7).

Величина T0 в списке формальных параметров означает вектор начальных значений температуры. Индексация начинается с единицы, т.е. значение переменной Origin среды Mathcad следует установить в единицу. Размерность вектора T0 равна числу узлов сетки.

Величины Tf и  Bi – векторы из двух элементов, задающие соответственно температуры среды и числа Био на левом и правом торце объекта (Рис. 11.1). Число Био – это безразмерный коэффициент теплоотдачи: Bi = α·Δx/λ.

Цикл for обеспечивает вычисления во внутренних узлах по формулам ( 11.6). Отдельно рассчитываются коэффициенты для граничных узлов, исходя из соотношения ( 11.7).

Параметр iTime сообщает подпрограмме, на каком шаге по времени находится процесс вычислений. Это важно для вычисления температуры среды, контактирующей с правым торцом. Задавая в векторе Pulse ненулевую амплитуду Ampl и некоторое значение частоты ν, мы сможем имитировать периодические тепловые воздействия на объект.

Программа Coef возвращает массив, собранный посредством функции augment из векторов A, B, C, D.

 

 

Рис. 11.4. Подпрограмма для вычисления матрицы коэффициентов

 

Mathcad-программа SysTRD

 

 

Рис. 11.5. Метод прогонки

 

Коэффициенты A, B, C, D для трехдиагональной системы уравнений подготовлены подпрограммой Coef. В теле процедуры значения коэффициента D замещаются вычисленными значениями температуры в узлах сетки. Итак, процедура SysTRD возвращает вектор решений системы.

Mathcad-программа TimeHistory

Программа TimeHistory организует вычисления по временным шагам. Параметр nTime задает число шагов по времени. Функция TimeHistory возвращает распределения температуры по координате x для последовательных моментов времени iTime как столбцы матрицы F.

 

 

Рис. 11.6. Главная программа

 

На каждом шаге по времени вызывается процедура Coef для вычисления коэффициентов A, B, C, D. Затем происходит обращение к решателю системы уравнений методом прогонки SysTRD. Обновляется вектор значений температуры T0 на предыдущем временном шаге и добавляется столбец к матрице F, хранящей распределения температуры для последовательных моментов времени.

 

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

Этап подготовки данных для расчета показан на Рис. 11.7. Предполагается исследовать температурные колебания в латунном стержне длиной 39 мм, если на правом торце задана температура, пульсирующая около среднего значения 800˚C  c амплитудой 320˚C , а на левом торце – постоянная нулевая температура. Эти граничные условия имитируются заданием соответствующих температур жидкости и таких значений сеточного числа Био, которые соответствуют очень большим величинам коэффициентов теплоотдачи на правом и левом торцах объекта.

Подготовительные вычисления на Рис. 11.7 понятны без комментариев. В конце этого фрагмента c с помощью функции matrix формируется вектор начальных температур T0 , задающий равномерное начальное распределение 100˚C, а также вектор значений относительного коэффициента теплопроводности Λ.

В данном примере величина Λ принимается постоянной. Если теплопроводность зависит от координат и/или температуры, то понадобится дополнительная процедура или вставка в программу Coef. Лучшим способом осреднения коэффициента теплопроводности будет вычисление среднегармонических значений по формуле (11.4). Мы оставляем эти усовершенствования программы как упражнение для читателей.

 

 

Рис. 11.7. Пример расчета: исходные данные

 

 

 

Рис. 11.8. Результаты расчета: температурные волны

 

Результаты вычислений показаны на Рис. 11.8. Правый график представляет серию распределений температуры вдоль оси  x для различных моментов времени. Более наглядным является трехмерное представление на левой диаграмме. Плоскость в основании графика построена на продольной координате x (отметки 1..40) и временной координате (отметки номера временного шага в пределах 1..200). По вертикальной оси отложена температура. Отчетливо видны вынужденные колебания температуры на одном из торцов ( отметка 40 по координатной оси) и уменьшение амплитуды колебаний по мере приближения к противоположному торцу.

Эффектное представление нестационарного температурного поля можно получить с помощью анимации, как это описано в главе 9 .

11.6.   Заключение

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

  1. Исследовать влияние коэффициента теплоотдачи на температуру горячего торца ( на ее максимальное и среднее значение)
  2. Исследовать влияние свойств материала на максимальную температуру поверхности и теплоотвод через стержень
  3. Исследовать глубину проникания температурных волн, взяв длинный стержень и задав адиабатические условия на левом торце
  4. Исследовать теплопередачу через стенку, если коэффициент теплоотдачи пульсирует во времени (для этого понадобится несколько модифицировать программу Coef, обеспечив вариации числа Био так же, как сейчас это сделано с температурой жидкости)
  5. Исследовать температурные режимы стенок цилиндра двигателя внутреннего сгорания с воздушным охлаждением
  6. Исследовать температурные режимы стен здания при погодных и сезонных изменениях температуры
  7. И т.д.

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

Последние три главы книги посвящены реализации в Mathcad’е численных методов для дифференциальных уравнений в частных производных. Достоинство разработанных на Mathcad’е простых программ в том, что они демонстрируют основные алгоритмы численного анализа, будучи доступными и открытыми для модификации и экспериментов. Более широкое представление о численном анализе как фундаменте современных вычислительных программ дают включенные в список литературы книги: [Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.,Ошибка! Источник ссылки не найден.]

Приложение. Встроенные функции для решения обыкновенных дифференциальных уравнений

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

odesolve (x,b,step)

Интегрирование одного дифференциального уравнения (первого, второго и т.д. порядка) при решении начальной или краевой задачи.

f = odesolve (x,b,step)

Решение возвращается как функция. Можно строить график функции f(x) или вычислять ее значения, например, f(6) (см. рисунок).

Аргументы:

Параметры x, b, step должны быть действительными (не комплексными). Рекомендуется поэскпериментировать со значениями параметра  step , чтобы оценить точность решения.  Можно выбрать применяемый при работе odesolve метод интегрирования: с фиксированным шагом (fixed step) или с адаптирующимся шагом (аdaptive step size) , дважды щелкнув мышкой на слове odesolve и отмечая Fixed или Adaptive в контекстном меню.

 

 

Применение odesolve для решения дифференциального уравнения третьего порядка

 

rkfixed (IC, x1, x2, npoints, D)

Решение с фиксированным шагом интегрирования.

 

S = rkfixed (IC, x1, x2, npoints, D)

S -  матрица, содержащая независимую переменную  x  в первом столбце, а в остальных столбцах - искомую функцию и ее первые n-1 производные.

Аргументы: 

 

 

Применение rkfixed для решения системы дифференциальных уравнений, эквивалентной уравнению Фолкнера-Скэн

 

rkadapt (IC, x1, x2, acc, D, kmax, s)

Решение с адаптирующимся шагом интегрирования.

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

S = rkadapt (IC, x1, x2, acc, D, kmax, s)

Аргументы:

 

 

Применение rkadapt для решения системы уравнений, эквивалентной уравнению Фолкнер-Скан

 

Rkadapt (y, x1, x2, npoints, D)

Аргументы Rkadapt такие же, как у функции rkfixed

Rkadapt  возвращает решение в равноотстоящих точках (всего в  1+ npoints  точках). На каждом таком интервале Rkadapt вызывает функцию rkadapt, для которой параметр точности acc  равен значению переменной среды Mathcad TOL.

 

 

Применение Rkadapt для решения системы  дифференциальных уравнений, эквивалентной уравнению Фолкнер-Скан

 

bulstoer (y,x1,x2,acc,D,kmax,s)

Bulstoer(y, x1, x2, npoints, D)

Решение гладких систем.

Аргументы и формат возвращаемого значения для этих функций такие же, как для rkadapt  и Rkadapt соответственно.

В справочной системе Mathcad функция  bulstoer  рекомендована для гладких (smooth) систем,  а  rkadapt - для медленно меняющихся систем (slowly varying systems).  Обычный пользователь вряд ли заметит разницу в применении этих функций. Сопоставление можно найти в книге [Ошибка! Источник ссылки не найден.].  Принципиальные различия методов интегрирования rkfixed, rkadapt, stiffr обсуждаются на вычислительных примерах в главе 4.

 

Stiffr (Y, x1, x2, npoints, D, J)

Stiffb(y, x1, x2, npoints, D, J)

stiffr (y,x1,x2,acc,D,J,kmax,s)

stiffb(y,x1,x2,acc,D,J,kmax,s)

Эти интеграторы специально разработаны для жестких систем.

Список аргументов аналогичен тому, какой имеют функции rkadapt (y,x1,x2,acc,D,kmax,s) или  Rkadapt (y, x1, x2, npoints, D).  Однако появляется дополнительный параметр J.

J есть функция, возвращающая  n·(n+1)-матрицу,  первый столбец которой содержит вектор производных от правой части системы  по независимой переменной, а остальные столбцы образуют матрицу Якоби для системы дифференциальных уравнений.

Stiffb использует метод Булирша - Штоера  для жестких систем, Stiffr - метод Розенброка [Ошибка! Источник ссылки не найден.].

Решение жестких систем подробно рассматривается в главе 4. Там же подробно обсуждаются принципиальные различия методов интегрирования rkfixed, rkadapt, stiffr.

 

sbval (v, x1, x2, D, load, score)

Решение двухточечной краевой задачи.

sbval  сводит двухточечную краевую задачу к задаче с начальными данными. Возвращает недостающие начальные условия  на левой границе x1.

Аргументы:

 

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

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

 

 

 

Решение двухточечной краевой задачи для уравнения Фолкнера-Скэн методом sbval

 

bvalfit (v1, v2, x1, x2, xin, D, load1, load2, score)

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

Аргументы:

 

 

Решение двухточечной краевой задачи для уравнения Фолкнера-Скан методом bvalfit

Применение функции bvalfit подробно рассмотрено в главе 10.

 

            Предметный указатель

A

augment, 83, 221

B

break, 201

Break, 201

bulstoer, 234, 235

bvalfit, 129, 206, 207, 209, 211, 237, 238

C

collect, 133, 199, 219

complex, 11, 137

Contour Plot, 112, 172

cspline, 105, 122

E

eigenvals, 78

F

factor, 92

Find, 141

for, 167, 201, 221

G

Given, 39, 45, 141

I

if, 40, 90, 91

interp, 105, 122, 169

L

last, 78

linterp, 112

M

matrix, 105, 111, 223

Minimize, 96, 100

multigrid, 28, 196, 197

O

Odesolve, 39, 44, 45, 96, 97, 105, 121, 122, 127, 128, 230, 231

Origin, 220

P

Parametric surface plots, 34, 158, 172

pspline, 169

Pspline, 169

R

relax, 28, 196, 197, 201, 202, 203, 204, 205, 212

rkadapt, 17, 60, 64, 71, 72, 73, 74, 75, 81, 96, 232, 233, 234, 235

Rkadapt, 234, 235

rkfixed, 64, 65, 66, 68, 69, 70, 71, 72, 81, 96, 105, 109, 110, 127, 137, 138, 231, 232, 233, 234, 235

rnorm, 18

root, 105, 122

S

sbval, 105, 108, 109, 110, 127, 129, 207, 211, 235, 236, 237

simplify, 136

solve, 11, 15, 80

stack, 105, 111, 168

stiffr, 64, 65, 73, 74, 75, 77, 82, 83, 235

Stiffr, 235

submatrix, 168

substitute, 92

V

Vector Field Plot, 90, 111, 112

Vectorize, 178

А

Автомодельные переменные, 106

Адиабатическое приближение, 77

Анимация, 173, 174, 224

Аррениуса закон, 51

Б

Баланса уравнение, 20

Био число, 200, 201, 221, 223, 225

Бифуркация, 56

В

Вдув, 86, 114, 115, 116, 120, 132

Вечная мерзлота, 211

Волновая скорость, 32, 33, 35, 36, 131, 149, 156, 158, 175, 178, 179, 187, 190, 191

Волновое уравнение, 31, 36, 151, 156, 161, 162, 187

Волновое число, 131, 135

Вывод дифференциального уравнения, 27

Вязкого трения закон, 24

Г

Газопровод, 212

Гармонического осциллятора уравнение, 10, 11, 47

Гаусса-Зайделя метод итераций, 200

Гистерезис, 60

Граничные условия, 20, 28, 84, 89, 92, 101, 104, 108, 120, 197, 199, 200, 203, 205, 206, 209, 212, 213, 216, 220, 223, 236

Д

Двухточечная краевая задача, 6, 105, 108, 127, 129, 207, 208, 211, 235, 236, 237, 238

Дивергенция, 21, 37

Ж

Жесткие уравнения, 6, 64, 65, 67, 68, 69, 71, 72, 73, 74, 75, 77, 78, 79, 81, 82, 148, 235

И

Интегральный метод, 88, 228

Источник, 7, 8, 20, 21, 22, 23, 29, 30, 37, 43, 44, 51, 87, 161, 185, 196, 197, 201, 207, 215, 217, 225

К

Катастрофа сборки, 5, 51, 56, 57, 58

Комплексная амплитуда, 131

Комплексная фазовая скорость, 135

Комплексная экспонента, 46, 134

Комплексное решение, 47, 131, 143

Контрольный объем, 20, 21, 22, 24, 26, 27, 28, 37, 44, 84, 147, 161, 162, 164, 166, 185, 187, 198, 217

Концентрация

автомобилей, 30

примеси, 146

пузырьков пара, 183

Коэффициент теплоотдачи, 51, 52, 55, 57, 59, 88, 89, 92, 93, 98, 99, 120, 121, 124, 126, 200, 201, 202, 205, 221, 224

Куранта число, 164, 166, 168, 175

Куэтта поток, 12, 13

Л

Ленгмюра изотерма, 150, 151, 152, 154, 156, 157, 169, 170, 173

Линеаризация, 82, 130

Логистическая модель, 15, 19, 153, 154

М

Мак-Кормака метод, 163, 164, 166, 168, 169, 171, 179

Метод малых возмущений, 133

Н

Неустойчивая стационарная точка, 53

Неустойчивость численная, 72, 81, 211

Неявная схема, 79, 81, 218

Нуссельта число, 121

Ньютона-Рихмана уравнение, 51, 92, 120, 216

О

Опрокидывание волны, 35, 158

Отсос, 84, 102, 117, 119, 129

П

Параметр в дифференциальном уравнении, 63

Паросодержание, 166, 184

Пограничный слой, 6, 26, 41, 42, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 97, 98, 101, 102, 103, 104, 105, 107, 109, 111, 113, 114, 115, 117, 120, 122, 126, 128, 129, 227, 228

Прандтля число, 93, 121, 122, 123, 124

Преобразование подобия, 104

Пристрелки (стрельбы) метод, 108, 109, 129, 211, 236

Прогонки метод, 7, 129, 216, 220, 222

Программа

Coef, 220, 221, 222, 223, 225

Euler, 138

Plate, 200, 201, 202, 203, 204, 205, 206, 212

SysTRD, 220, 221, 222

TimeHistory, 172, 179, 222

Пуассона уравнение, 28, 29, 196, 197, 202

Р

Радиационный теплоотвод, 207

Разрывная правая часть, 207

Рейнольдса число, 91, 103, 108, 126, 133, 136

Релаксации параметр, 69, 148

Релаксационное уравнение, 148

Рэлея уравнение, 6, 26, 41, 130, 131, 136

С

Семенова критерий, 57

Среднеквадратичные пульсации, 143, 144

Стандартный закон теплообмена, 89

Стефана-Больцмана постоянная, 207

Стохастические модели, 16

Т

Теплового взрыва теория, 51

Тепловой пограничный слой, 84, 85, 86, 87, 88, 89, 90, 93, 95, 102, 103, 120, 122, 126, 129

Тепловыделяющий элемент ядерного реактора, 38, 212, 213

Теплопроводности закон, 23

Толщина потери энтальпии, 85, 86, 87, 89, 94, 95, 96, 98, 126

Трехдиагональная система, 219, 222

У

Ударные волны, 6, 7, 28, 33, 36, 40, 159, 161, 163, 165, 169, 173, 175, 179, 181, 182, 192, 193, 224

Уоллиса модель, 187, 229

Устойчивые стационарные точки, 53

Ф

Фазовый портрет, 53, 56

Фолкнера-Скэн уравнение, 6, 26, 41, 42, 84, 101, 102, 104, 105, 106, 108, 109, 127, 129, 211, 230, 232, 237, 238

Функция тока, 101, 104, 105, 106, 107, 108, 111, 112, 114, 116, 118, 119, 121, 122

Фурье число, 23, 217, 219

Х

Характеристики, 6, 32, 34, 35, 36, 156, 157, 158, 159, 162, 169, 173, 178, 179, 187, 192, 193

Характеристическое уравнение, 46, 47, 48

Химической кинетики уравнения, 75, 76, 77, 78

Э

Экспоненциальный рост, 13, 14, 15, 16, 43, 135, 137, 141, 152

Я

Явная схема, 79, 218

Якоби матрица, 74, 76, 78, 82, 83, 235