Министерство образования и науки Российской Федерации
Федеральное государственное бюджетное образовательное учреяадение
высшего профессионального образования
«Тамбовский государственный технический университет»
ВЛ. Ляшков
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
И АЛГОРИТМИЗАЦИЯ ЗАДАЧ
ТЕПЛОЭНЕРГЕТИКИ
Утверждено Ученым советом университета в качестве
учебного пособия для магистрантов первого года обучения
по направлению 140100.68 «Теплоэнергетика и теплотехника»
(магистерская программа 140100.01 «Технология производства
электрической и тепловой энергии»)
Тамбов
Издательство ФГБОУ ВПО «ТГТУ»
2012
УДК 621.1.001.575(075)
ББК ^3в631я73
Л993
Рецензенты:
Кандидат педагогических наук, доцент кафедры информатики и
информационных технологий ФГБОУ ВПО «ТГУ им. Г.Р. Державина»
Е.В. Клыгина
Кандидат технических наук, доцент кафедры «Информационные
процессы и управление» ФГБОУ ВПО «Тамбовский
государственный технический университет»
С.И. Татаренко
Ляшков, В.И.
Л993 Математическое моделирование и алгоритмизация за-
дач теплоэнергетики [Электронный ресурс] : учебное посо-
бие / В.И. Ляшков. - Тамбов : Изд-во ФГБОУ ВПО «ТГТУ»,
2012.-140 с.
Рассмотрены вопросы разработки математических моделей
различных теплоэнергетических процессов, методы реализации их
с помощью ПК, приемы решения типовых задач, возникающих при
математическом моделировании. При этом особенно выделены
объекты теплотехнического назначения.
Предназначено для магистрантов первого года обучения по
направлению 140100.68 «Теплоэнергетика и теплотехника» (маги-
стерская программа 140100.01 «Технология производства электри-
ческой и тепловой энергии»).
УДК 621.1.001.575(075)
ББК ^3в631я73
© Федеральное государственное бюджетное
образовательное учреждение высшего
профессионального образования
«Тамбовский государственный технический
университет» (ФГБОУ ВПО «ТГТУ»), 2012
ВВЕДЕНИЕ
Прошедшие два десятилетия наглядно показали, что при сохране-
нии существующих условий развития общества и экономики России в
ближайшие годы нам предстоит опуститься до уровня бывших коло-
ниальных, а ныне так называемых развивающихся стран Африки и
Латинской Америки. Это наконец-то дошло и до правящей элиты.
В последние год-два все громче и громче заговорили о необходимости
существенной модернизации в области промышленного производства,
сельского хозяйства, образования и науки, политической организации
общества, регулирования рыночных отношений и др. Восстановление
и развитие промышленности должно основываться на применении
новых методов организации производства, на внедрении новейших
технологий и оборудования. Только это позволит обеспечить выпуск
продукции, конкурентоспособной на мировом рынке.
Энергетика является одной из ключевых отраслей экономики лю-
бой страны. Именно энергетика обеспечивает энерговооруженность, а
значит, и производительность труда, создавая предпосылки для высо-
кой экономической мощи страны и высокого благосостояния ее насе-
ления. К сожалению, поспешно осуществленное реформирование
энергетики не принесло обещанных результатов (об этом неоднократ-
но указывал действующий Президент страны), и это вызывает неуве-
ренность в осуществимости планов, зафиксированных в «Энергетиче-
ской стратегии России на период до 2020 года», принятой в 2003 г.
В этом правительственном документе установлено, что для пре-
одоления проблем, тормозящих развитие производства, ежегодный
темп роста потребления (а значит, и производства) энергии должен
составлять почти 5 %. Однако за истекшие 7 лет таких темпов роста ни
разу не было достигнуто. В принятом в 2009 г. новом документе, уточ-
няющем стратегию до 2030 г., невыполнение намеченных ранее рубе-
жей объясняется мировым экономическим кризисом, а все новые гори-
зонты прописаны общими направлениями с минимальными цифровы-
ми прогнозами, дабы не краснеть потом за провалы в работе.
В настоящее время по производству энергии мы вышли только на
уровень 1990 г., но это уже не способно обеспечить все возрастающие
потребности страны. Сегодня практически в каждом четвертом субъ-
екте Российской Федерации ощущается большой дефицит энергии,
обусловленный практической выработкой ресурсов работоспособно-
сти генерирующих установок, предельным износом электрических и
теплофикационных сетей и неурегулированными экономическими от-
ношениями между субъектами вновь созданной инфраструктуры,
предназначенной для энергообеспечения потребителей энергии.
3
Будущие технические специалисты, умом и руками которых
должны быть решены задачи модернизации в области энергетики и
энергомашиностроительного производства, должны знать и уметь ис-
пользовать в своей практической деятельности современные эффек-
тивные методы научных исследований, поскольку это позволяет со-
кратить продолжительность, уменьшить стоимость выполняемых ра-
бот и обеспечить высокий качественный уровень принимаемых реше-
ний. К таким эффективным методам в первую очередь следует отнести
математизацию исследований, как теоретических, так и эксперимен-
тальных. Это предполагает обязательное использование математиче-
ских моделей, достаточно точно (адекватно) описывающих объект ис-
следования. При этом возможны 2 варианта получения таких моделей.
Теоретический подход, основанный на использовании известных зако-
нов природы и вытекающих из них взаимосвязей между отдельными
параметрами объекта в виде алгебраических или более сложных урав-
нений или формул. Другой, экспериментальный подход, основан на
предварительном планировании экспериментов, проведении их в соот-
ветствии с разработанным оптимальным планом и специальной мате-
матической обработке результатов, позволяющей получить достаточно
простую и пригодную для неширокого диапазона изменения входных
параметров объекта математическую модель.
Отдавая на суд читателя свой скромный труд, автор согревает се-
бя надеждой, что его работа, хотя бы в малой мере, поможет студенту
понять основные подходы и особенности разработки и реализации ма-
тематических моделей в области теплоэнергетики.
4
1. ОБЩИЕ ПОНЯТИЯ И ПОДХОДЫ
1.1. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
КАК СРЕДСТВО РЕШЕНИЯ ПРАКТИЧЕСКИХ ЗАДАЧ
Моделирование присуще человеческой природе и возникает во
всех видах творческой активности людей. Каждый из нас всегда стара-
ется смоделировать предстоящую ситуацию с тем, чтобы попытаться
предвидеть ее результаты. Великому Шекспиру приписывают крыла-
тую фразу: «Весь мир - театр, и все мы в нем актеры», которая как
нельзя лучше демонстрирует приверженность человека к моделирова-
нию. Если подробнее разбираться с моделированием, то можно выде-
лить ситуационное моделирование (сейчас я играю роль профессора, а
вы — роли прилежных студентов), физическое моделирование, связан-
ное с созданием физических моделей устройств и процессов (подробно
такое моделирование мы изучали при знакомстве с теорией подобия) и
математическое моделирование, в настоящее время переходящее все
больше и больше в информационное. С ним-то и предстоит нам более
подробно познакомиться.
Так почему же моделирование получило такое широкое распро-
странение? Да потому, что моделирование — это изучение свойств и по-
ведения объекта исследований без использования для этого самого объ-
екта, без внедрения в него. Работа с моделью позволяет относительно
быстро и дешево исследовать свойства объекта, его поведение в любых
мыслимых (и даже нереальных) ситуациях. С помощью математической
модели мы можем имитировать процессы и явления, не реализуемые в
реальной жизни (например, авиационный тренажер, на котором отраба-
тывается аварийная ситуация, способная при неправильных действиях
экипажа привести к катастрофе, или такой режим работы химического
реактора, который непременно приведет к взрыву аппарата и т.п.).
Математическое моделирование появилось с внедрением матема-
тики во все области человеческой деятельности. Это проникновение в
свое время позволило великому философу и математику И. Канту ска-
зать: «В каждой дисциплине науки столько, сколько в ней математи-
ки». Так что математическое моделирование появилось с возникнове-
нием точных наук. Сегодня практически все инженерные дисциплины
(для примера вспомним любимые «Теоретические основы теплотехни-
ки») в той или иной мере посвящены разработке методик инженерного
расчета отдельных процессов, устройств, конструкций. И каждая такая
методика включает некие аналитические связи, полученные на основе
известных законов природы с учетом конкретных обстоятельств (это и
есть математическая модель) и описания последовательности приме-
нения этих связей для того, чтобы получить конечный результат (а это
и есть описание алгоритма расчетов).
5
За последние пятьдесят лет роль математического моделирования
исключительно высокими темпами возрастает, поскольку человек по-
лучает для этого все более и более совершенный инструмент в виде
персонального компьютера (ПК), обеспечивающего обработку гигант-
ского количества информации с невообразимо высокой скоростью.
Это позволяет разрабатывать математические модели для все более
сложных объектов, точно и всесторонне учитывать влияние на объект
различных факторов, которые ранее для упрощения модели просто не
учитывались и это отражалось на точности и достоверности результа-
тов моделирования. Сегодня технические, экономические и другие
объекты достигли такой степени сложности, что уже не поддаются
обычным методам теоретического исследования, позволяющим полу-
чать какие-то аналитические решения в виде неких алгебраических
(или более сложных) формул. Поэтому переход к математическому
моделированию является одним из векторов, составляющих научно-
технический прогресс.
При математическом моделировании реализуются взаимодейст-
вие и взаимовлияние, отражаемые схемой на рис. 1.1. Разработка ма-
тематической модели ведется в итерационном режиме выполнением
нескольких последовательных этапов. Сначала формулируется перво-
начальный «математический эквивалент», описывающий в математи-
ческой форме важнейшие свойства объекта: геометрические, физиче-
ские, уравнения, отражающие основные законы природы (закон сохра-
нения массы, сохранения энергии, сохранения количества движения,
закон Ньютона и др.). Как правило, этот «математический прообраз»
объекта представляется системой дифференциальных уравнений (обыч-
ных или в частных производных) и набором условий однозначности.
6
После этого проводится разработка алгоритма, т.е. такой последова-
тельности действий, реализация которых приведет к решению постав-
ленной при моделировании задачи. Для реализации алгоритма на ПК в
соответствии с ним разрабатывается программа для ПК на языках вы-
сокого уровня. При этом широко используются «стандартные» проце-
дуры применения тех или иных вычислительных алгоритмов. Триада
модель - алгоритм - программа в процессе ее разработки требует не-
однократного согласования и уточнения с реальным объектом. В ре-
зультате этого в отдельную часть ее вносятся определенные корректи-
вы, которые должны быть учтены в двух других частях триады. В ре-
зультате такой работы и возникает работоспособная математическая
модель, которая подвергается окончательной доработке, и это принято
называть отладкой модели и программы.
На отлаженной программе и проводится численное исследование
свойств и поведения объекта путем расчетов интересующих исследо-
вателя параметров при различных условиях однозначности. Отдельные
результаты таких расчетов затем сравнивают с результатами экспери-
ментальных исследований при таких же условиях проведения опытов.
Если опытные и расчетные результаты удовлетворительно совпадают,
то говорят, что разработанная модель адекватна, т.е. с приемлемой
точностью описывает реальную действительность. В силу большой
сложности исследуемых сегодня объектов добиться полной адекватно-
сти часто не удается. Но даже тогда, когда разработанная модель обес-
печивает хотя бы принципиальное соответствие с экспериментальны-
ми результатами по тенденциям, по влиянию входных факторов на
результаты опыта и т.п., считается, что полученная модель вполне
приемлема для практики, поскольку введением в нее специальных
экспериментально определенных поправочных множителей или спе-
циальных функциональных зависимостей для расчета этих поправок
позволяет обеспечить адекватность модели.
Основные идеи и термины, положенные в основу математическо-
го моделирования, наглядно представляются при знакомстве со схемой
на рис. 1.2.
Рис. 1.2. Схема анализа объекта при разработке его математической модели
7
Будем называть объектом исследования некоторую часть пред-
метной области, поведение и свойства которой нас интересуют.
Все параметры, через которые происходит взаимодействие между
объектом и окружающей средой, принято делить на 4 группы. Входные
параметры X - это различные потоки (вещества, энергии, информа-
ции) из окружающей среды к объекту. Считается, что эти параметры
мы можем измерить, но изменить их состав и интенсивность при за-
данных значениях выходных параметров Y мы не можем (например,
массовый расход воды, подогреваемой в теплообменнике до заданной
температуры).
Величины выходных параметров Y определяются величинами
всех остальных параметров. Управляющие параметры U - это тоже
входные параметры, но изменением их мы можем воздействовать и
изменять величины выходных параметров (например, температура
греющего теплоносителя в теплообменнике, изменением которой
можно влиять на конечную температуру подогреваемого теплоносите-
ля). Возмущающие параметры (или шумы) Z способны влиять на вы-
ходные параметры, но это влияние происходит без нашего участия,
они возникают и изменяются по величине случайным образом. Если
влияние этих параметров на Y велико и даже определяет поведение
объекта, то объект называют стохастическим и для его описания ис-
пользуют вероятностные характеристики. Если величины Y очень мало
или совсем не зависят от возмущающих параметров, то объект и его
математическую модель называют детерминированной.
Зависимость, которая в математической форме однозначно опи-
сывает связь между всеми входными и выходными параметрами объ-
екта, называют математической моделью объекта:
Y= ф(^, U, Z).
На рисунке 1.3 приведена классификация математических моде-
лей по пространственным и временным признакам.
Рис. 1.3. Классификация математических моделей
8
Если параметры объекта изменяются только по времени, а в про-
странстве параметры одинаковы, то говорят, что объект «с сосредото-
ченными параметрами». В другом, более частом случае, когда пара-
метры меняются и по времени, и в пространстве, объект называют «с
распределенными параметрами». Конечно же, для объектов с распре-
деленными параметрами математическое описание получается гораздо
более сложным. Чтобы в какой-то мере упростить задачу, широко
применяется промежуточная (ячеистая) модель, в которой весь объем
объекта делится на несколько частей и считается, что в пределах каж-
дой части параметры меняются только по времени, а вот между частя-
ми имеет место неравномерность распределения параметров в про-
странстве. Такая же идея используется в квазистационарных моделях:
весь временной интервал разбивается условно на несколько частей и
считают, что в течение каждой такой части процессы в системе уста-
новившиеся, стационарные, а на концах этих временных интервалов
происходит скачок значений параметров объекта.
Приведенная на рис. 1.4 схема наглядно отражает состав матема-
тического описания объекта, составляющего его математическую мо-
дель. Конечно же, основу математической модели обычно составляют
материальные и энергетические балансы по соответствующим потокам
вещества и энергии. В качестве составляющих таких балансов обычно
используются отдельные теоретические и эмпирические соотношения,
описывающие взаимосвязи между отдельными внешними и внутрен-
ними параметрами объекта. В математическое описание включаются
также математические зависимости между параметрами и известные
закономерности протекания «элементарных» процессов (например,
теплопроводности, теплопередачи и т.п.), из суммы которых и склады-
вается взаимодействие между объектом и окружающим миром.
Рис. 1.4. Состав математического описания объекта
9
Эти закономерности описываются известными физическими за-
конами (закон Фурье для теплопроводности, закон Фика для массооб-
менных процессов и т.п.). Особое значение играют ограничения на
переменные параметры объекта. При выходе за такие предельные ог-
раничения может теряться физический смысл результата расчетов,
получаются такие значения выходных параметров (например, прибли-
жающуюся к нулю производительность), которые не представляют
практической ценности или требуются такие управляющие параметры,
которые не реализуемы на практике.
Одно из достоинств математического моделирования состоит в
том, что различные по характеру процессы могут иметь сходное мате-
матическое описание, что упрощает и облегчает алгоритмизацию мо-
дели (вспомните, например, электро-, тепло-, гидроаналогию, о кото-
рой говорилось в «Теоретических основах теплотехники»). Более того,
даже для очень трудно формализуемых объектов удается построить,
используя математические описания отдельных физических явлений,
математическую модель, позволяющую делать очень важные для
практики выводы. Для примера назову одну такую задачу из [1]: «Ди-
намика распределения власти в иерархии», где в качестве основного
уравнения, описывающего взаимодействие власти и общества, исполь-
зовано хорошо известное нам дифференциальное уравнение теплопро-
водности, переформулированное через специфические для исследуе-
мого объекта (властные структуры в обществе) параметры, такие как
уровень властных полномочий, оговоренных Конституцией и законами
страны (некий аналог температуры), поток власти, осуществляемый
через пресловутую вертикаль власти (аналог теплового потока), ско-
рость изменения величины власти, определяемая уровнем коррупции в
субъекте федерации (аналог коэффициента температуропроводности),
реакция общества в виде зачатков гражданского общества (аналог
внутренних источников или стоков теплоты) и др. Такая модель по-
зволила установить условия стабильного существования власти и об-
щества без тенденций к переходу в диктатуру или анархию, прогнози-
ровать кризисы власти различного характера, изучать реакцию граж-
данского общества на изменения в законодательстве или его право-
применения и др.
Общее требование к математическим моделям: они должны со-
держать разрешимые уравнения и системы уравнений. При этом полу-
чить алгебраические решения (в виде алгебраических формул) удается
только для самых простых объектов, и тогда говорят, что получено
решение «в квадратурах». Для более сложных задач, когда математи-
ческое описание включает трансцендентные уравнения или системы
дифференциальных уравнений (и особенно нелинейных или в частных
производных), решение может быть получено только численным спо-
10
собом. Разработкой методов и алгоритмов решения таких задач зани-
мается специальный раздел математики, так и называемый «Числен-
ные методы» [2]. При этом наибольшее применение находят методы
последовательных приближений, метод сеток, метод конечных эле-
ментов и вариационные методы [3, 4]. С некоторыми из них нам с вами
предстоит еще познакомиться.
Средством реализации математической модели является ПК, реа-
лизующий соответствующую программу расчетов, разработанную в
соответствии с выбранным алгоритмом решения поставленной задачи.
Напомним, что алгоритмом принято называть ту последовательность
действий, которая непременно приведет к решению поставленной за-
дачи.
Часто задачей математического моделирования является исследо-
вание особенностей поведения объекта при взаимодействии его с
внешней средой. Очень широко математические модели используются
для оптимизации в каком-либо смысле отдельных характеристик (на-
пример, экономической эффективности) объекта. Под оптимизацией
понимают определение таких режимных, конструктивных или других
управляющих факторов, при которых объект получает наиболее жела-
тельные качества. В рыночной экономике этим качеством является,
как вы понимаете, величина прибыли на единицу затрат, необходимых
для функционирования объекта.
1.2. ПРИЕМЫ РАЗРАБОТКИ МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ
Важным приемом, который широко используется при разработке
математических моделей, является декомпозиция сложных техниче-
ских объектов [5]. Чтобы исследовать сложный технический объект,
его мысленно разбивают на отдельные, в меру самостоятельные части
(элементы), определенным образом взаимодействующие между сбой.
Обычно такая декомпозиция ведется в соответствии с внутренней
структурой объекта. Такой прием позволяет построить модель объекта
как некую совокупность моделей отдельных его элементов. Исследо-
вание отдельных элементов позволяет выявить «слабое звено» среди
выделенных элементов и, в первую очередь, предложить мероприятия
для увеличения эффективности именно этого звена, а затем и объекта в
целом.
1.	В простейшем случае и применительно к энерготехнологиче-
ским объектам, при анализе сложную структуру представляют в виде
ряда простых структур с последовательным или параллельным соеди-
нением элементов. При этом, упрощая реальную действительность,
принимается допущение, что изменения КПД одного элемента не
влияют на величину КПД всех остальных элементов.
11
Для примера рассмотрим задачу разработки математической мо-
дели, описывающей рабочий процесс простейшей газотурбинной уста-
новки (ГТУ), состоящей из последовательно включенных компрессора,
камеры сгорания и непосредственно газовой турбины. При последова-
тельном соединении трех элементов (рис. 1.5), в которых происходит
преобразование различных видов энергии Е с различными КПД: rji, г|2
и т|3, энергия, передаваемая во 2-й элемент, определится по формуле:
ЕА -^вхП 1 •
Аналогично рассчитаем энергию, передаваемую в 3-й элемент:
Ев = Eav[2 = E^viz.
Энергия на выходе из группы элементов будет:
Явых = ЕвПз = Евхг[ 1П2Пз •
КПД группы последовательно включенных элементов найдется
как отношение энергии на выходе к энергии на входе для всей группы:
Ппосл ^вых^вх = П1П2'ПЗ-
Вывод: при последовательном соединении элементов КПД всей
группы определяется величиной произведения КПД всех элементов, и
увеличение одного из КПД во столько же раз увеличивает КПД всей
группы элементов.
Возьмем другой пример: на крупной ТЭЦ работает паровой котел
общей мощностью 300 МВт. Свежий пар от этого котла подается на
3 паровых турбины мощностью 150, 100 и 50 МВт соответственно. Вы-
рабатываемая электроэнергия всех турбин отдается в единую сеть. Об-
щая схема объекта соответствует параллельному включению элементов.
При параллельном включении элементов все они работают с об-
щим входом и общим выходом (см. рис. 1.6). Здесь также каждый эле-
мент может иметь свое значение КПД (т|ь т|2, т|3) и увеличение одного из
них приводит к увеличению КПД всей группы, но влияют они на общий
КПД по-другому. Связь между Евых и Евх можно получить, если ввести
весовые коэффициенты Ki9 показывающие, какая доля Евх подводится к
данному элементу. Величины К; должны быть известны по постановке
задачи или определены предварительно специальным расчетом.
Рис. 1.5. Последовательное включение элементов
12
Рис. 1.6. Параллельное включение элементов
В этом случае можем записать: Ел = ЕВХК]Г\], Ев = ЕвхК2Г\г, Ес =
= £'вх^3т]з. Тогда величина ЕВЪ1Х определится суммой Евых = ЕВХ(Е]Г|] +
+ К2Х\2 + ^зПз) и КПД параллельной схемы будет т|пар = Евых/Евх = Е]Г| i +
+ К2^\2 + Е3г|з. Из формулы видно, что наибольшее влияние на величи-
ну Л пар будет оказывать тот элемент, у которого произведение
наибольшее, а самым проблемным является элемент с наименьшим
значением т|(. Именно с него и следует начинать разработку мероприя-
тий, направленных на повышение КПД всей схемы.
2.	Кроме рассмотренных простейших структур, на практике дос-
таточно широко встречаются более сложные схемы соединения от-
дельных частей объекта, получившие название схем с обратной свя-
зью. Такие схемы включают байпасирование или рециркуляцию пото-
ков. Типичные два варианта таких схем приведены на рис. 1.7.
б)
Л(у)
Рис. 1.7. Схемы структур с обратной связью:
а — отрицательная обратная связь; б — положительная обратная связь
13
При отрицательной обратной связи какая-то часть входного пара-
метра Xi, подвергаясь преобразованиям в дополнительном устройстве,
суммируется с выходом основного блока и они формируют общий вы-
ходной сигнал у. При положительной обратной связи все происходит
наоборот: малая часть выходного параметра направляется в дополни-
тельное устройство, где подвергается некоторому преобразованию
(например, усилению), и результат этого преобразования суммируется
с величиной входного параметра хь
Рассмотренные схемы с обратной связью обычно используются
для организации автоматического регулирования и поддержания на
заданном уровне режимов работы энергетического оборудования. Ко-
нечно же, обратные связи существенно усложняют математическое
описание объекта. Однако же, если известны статические характери-
стики и основного, и дополнительного блоков (известны виды функ-
ций f_ и/2), то статическая характеристика объекта с обратной связью
может быть определена. Если в простейшем случае это линейные за-
висимости /1(х) = к\х и^(х) = к?х, то для объекта с отрицательной об-
ратной связью получаем: у =/i(xi — Х2) +/2(^2) = #i(xi — Х2) + к^х^ = к\Х\ +
+ (к2 -
3.	В восьмидесятые годы прошлого века стала широко приме-
няться теория графов, позволившая расширить возможности метода
декомпозиции. Граф - это совокупность отдельных вершин и соеди-
няющих их дуг. С помощью графа наглядно отображаются связи меж-
ду отдельными элементами, выделенными при декомпозиции объекта,
его топологические свойства, как правило, отражающие особенности
технологической схемы.
Так, систему основных элементов и связей, моделирующих теп-
лосиловую установку, показанную на рис. 1.8, можно представить в
виде графа, у которого каждому элементу соответствует вершина, а
связи между элементами отражены соответствующими дугами. Чтобы
облегчить перенос такого графа и его последующий анализ на ПК, его
описывают с помощью матрицы соединений вершин графа. Ее назы-
вают структурной матрицей.
Рис. 1.8. Упрощенная тепловая схема паросиловой установки
14
На следующем рисунке (рис. 1.9) приведен граф, отражающий
структуру и связи между элементами паросиловой установки. Узлы
графа пронумерованы римскими цифрами; I - паровой котел, II - ци-
линдр высокого давления, III — регенеративный подогреватель конден-
сата, IV - цилиндр низкого давления, V - конденсатор, VI - электроге-
нератор. Каждая связь тоже пронумерована. Расшифруем приведенные
связи: 1, 2 и 3- потоки топлива, воздуха и дымовых газов, соответст-
венно; 4 - поток перегретого пара в ЦВД; 5 - поток пара в ЦНД; 6 и
7 — потоки механической энергии, передаваемой валом турбины элек-
трогенератору; 8 - поток электроэнергии, отдаваемый потребителю;
9 - поток отработанного пара в конденсатор установки; 10 и 11 — пото-
ки подводимой и отводимой охлаждающей воды в конденсаторе; 12 -
поток конденсата в регенеративный подогреватель; 13 - поток пара
для подогрева конденсата; 14 — поток подогретого конденсата в котел.
Чтобы формализовать описание объекта при таком подходе, ис-
пользуются две матрицы. В структурной матрице (табл. 1.1) для каж-
дого элемента указывают входящие связи числом +1, а выходящие
связи числом -1 (или только знаками + и -). Как видно из таблицы,
каждая внутренняя связь является одновременно входной и выходной
для разных узлов.
Кроме структурной матрицы (или матрицы соединений) делается
еще матрица видов связей (табл. 1.2), в которой некоторым числовым
признаком указывается вид энергоносителя (например: 1 — перегретый
пар, 2, — влажный пар, 3 - конденсат, 4 — питательная вода, 5 - топли-
во, 6 - дымовые газы, 7 - механическая энергия, 8 - электрическая
энергия, 9 — вода из системы охлаждения конденсатора, 10 - воздух).
Из специальной базы данных по таким числовым кодам компьютер
может самостоятельно выбрать для каждого элемента и каждого пото-
ка всю необходимую исходную информацию для моделирования объ-
екта: геометрические размеры, расчетные формулы, теплофизические
характеристики и многие другие сведения.
Рис. 1.9. Граф предыдущей схемы
15
1.1. Структурная матрица
№ дуги	Узлы (элементы схемы)					
	I	II	III	IV		VI
1	+1					
2	+1					
3	-1					
4	-1	+1				
5		-1		+1		
6		-1				+1
7				-1		+1
8						-1
9				-1	+1	
10					+1	
11					-1	
12			+1		-1	
13		-1	+1			
14	+1		-1			
1.2. Матрица видов связей
№ дуги	Узлы (элементы схемы)					
	I	II	III	IV	V	VI
1	10					
2	5					
3	6					
4	1	1				
5		1		1		
6		7				7
7				7		7
8						8
9				2	2	
10					9	
11					9	
12			3		3	
13		1	1			
14	3		3			
16
1.3.	СОДЕРЖАНИЕ И СОСТАВ ТЕОРЕТИЧЕСКОЙ
МАТЕМАТИЧЕСКОЙ МОДЕЛИ
Теоретической (или неформальной) моделью называют матема-
тическую модель, которая получена в результате использования ос-
новных законов природы для описания свойств объекта. Антиподом
такой модели является формальная модель (обычно для очень слож-
ных объектов), которую получают на основании экспериментальных
исследований объекта. Отражая те же самые законы природы, фор-
мальная модель не содержит математического описания этих законов,
а представляется часто некоторыми эмпирическими зависимостями,
аппроксимирующими экспериментальные данные.
Основу формальной математической модели составляет система
балансовых уравнений для объекта в целом и для каждого из узлов
графа. При установившихся процессах (когда параметры объекта не
изменяются по времени) для каждого узла можно записать следующие
уравнения:
1.	Уравнение баланса энергии
#	м
£(Y1GA),.-£(Y2P)y=0,
Z=1	J=1
где у] и у2 - коэффициенты, учитывающие потери энергии в окружаю-
щую среду; G - расход энергоносителя; h - энтальпия теплоносителя;
Р - мощность механической или электрической связи; N и М — число
элементов с тепловыми связями и число элементов со связями механи-
ческими и электрическими, соответственно.
2.	Уравнение баланса расходов теплоносителей
N
YGi=Q-
И=1
3.	Уравнение гидравлического баланса для каждого узла
(Рвх — ДР — Рвых) “ 0 •
4.	Уравнение изменения энтальпии
(й^ -kh- Авых) = 0 (или + ДА - Авых) = 0 при подводе теплоты).
В общем случае изменения энтальпии Ай и давления Др зависят
от совокупности параметров связей Лпс, а также от совокупности кон-
структивных параметров Ак. Например, изменение энтальпии находит-
ся так:
ДА (Авх Авых)г|, где т| У(ЛПС, у4к).
17
Термодинамические, расходные и конструктивные параметры Лпс
и Лк при моделировании не могут приниматься произвольно, поэтому
математическая модель должна быть дополнена соответствующими
ограничениями:
^^ncmin < ^пс < ^псгпах ® 4стт < ^ктах ?
однозначно определяющими область существования решений. Область
существования решений может быть ограничена и более сложными
зависимостями типа у < (p(%i, Х2).
В зависимости от исследуемого объекта или процесса, их матема-
тическая модель может выражаться в виде некоторой линейной или
нелинейной функции, в форме замкнутой системы линейных уравне-
ний или в виде системы нелинейных уравнений.
При нестационарных процессах или для объектов с распределен-
ными параметрами обычно математическая модель представляется
системой линейных или нелинейных дифференциальных уравнений
или уравнений в частных производных, отражающих вышеназванные
законы в дифференциальной форме, совокупно с соответствующими
условиями однозначности. Решение таких систем может представлять
большие сложности, поэтому необходимо познакомиться с методами
реализации различных математических моделей на ПК при моделиро-
вании и оптимизации объектов, чему будет посвящено несколько по-
следующих глав.
1.4.	ХАРАКТЕРИСТИКИ ОБЪЕКТОВ МОДЕЛИРОВАНИЯ
Как уже говорилось, математическая модель позволяет изучать
свойства и поведение объекта при изменении его входных параметров.
Обычно зависимости выходных параметров от величин входных пара-
метров принято называть характеристиками объекта. Для большинства
объектов теплоэнергетики и число входных управляющих параметров,
и число выходных параметров всегда достаточно большое. Поэтому
характеристики рассчитывают по каждому из выходных параметров:
У; =fixi, Х2, Хз, ..., хт), При этом если число входных параметров т > 2,
говорят, что координатные оси Xi, хт образуют фазовое гиперпро-
странство, в котором и отражаются все возможные состояния парамет-
ра у(. При некоторых зафиксированных численных значениях величин
Xi, Х2, ... в соответствии с функцией f мы получаем точку в этом фазо-
вом пространстве. Все множество таких состояний образует некую
поверхность фазовых состояний. При изменении только одного из
входных параметров Xj на этой поверхности образуется некая изолиния
(например, изотерма), а при изменении двух и более входных парамет-
ров — линия, которую называют фазовой траекторией. Для того, чтобы
18
наглядно отразить зависимость величины выходного параметра от ве-
личины входного используется проекция фазовой траектории на одну
из координатных плоскостей и именно такие кривые принято называть
отдельными характеристиками объекта. В курсе «Тепловые двигатели
и нагнетатели» мы широко использовали такие характеристики
(вспомните скоростные и нагрузочные характеристики ДВС или ха-
рактеристики поршневых компрессоров).
Если исследуется объект, функционирующий при установившем-
ся режиме, то характеристики его называют статическими. Упомяну-
тые выше характеристики — статические.
Однако в отдельных случаях объектом исследования являются
устройства или процессы, работающие при неустановившемся режиме
(процесс от запуска установки до выхода на установившийся режим,
задачи регулирования режимов работы и многие другие). Такие харак-
теристики, показывающие как изменяется анализируемый выходной
фактор по времени, называют динамическими. Для практики очень
важны два вида динамических характеристик:
1. Это так называемые кривые разгона, при определении которых
изменение входного управляющего параметра осуществляется в виде
некоторого скачка (рис. 1.10). Снятые экспериментально кривые раз-
гона позволяют подобрать математическую модель объекта в виде
дифференциального уравнения первого или второго порядка по време-
ни.
2. Характеристики переходных процессов, в которых изменение
входного управляющего параметра осуществляется по некоторому
заранее заданному закону (рис. 1.11).
Математическое моделирование позволяет рассчитывать статиче-
ские и динамические характеристики исследуемых объектов (решать
прямую задачу), а также, используя специальные методы оптимиза-
ции, выбирать такие управляющие воздействия, при которых фазовая
траектория процесса будет наилучшей в некотором заранее заданном
смысле (в этом случае говорят о решении обратной задачи моделиро-
вания или об оптимизации переходного процесса).
Рис. 1.10. Кривая разгона объекта
Рис. 1.11. Кривая переходного
процесса
19
1.5.	ВОПРОСЫ И ЗАДАНИЯ ДЛЯ САМОПРОВЕРКИ
(могут использоваться в качестве вопросов при проведении зачета)
1.	Какие виды моделирования Вы знаете? Чем они схожи и в
чем их принципиальное различие?
2.	Какие преимущества обеспечивает математическое модели-
рование по сравнению с исследованием свойств объекта непосредст-
венно на самом объекте?
3.	Опишите самую общую схему создания математической мо-
дели объекта.
4.	Какие математические зависимости лежат в основе математи-
ческой модели объекта?
5.	Что называют алгоритмом реализации математической моде-
ли? В чем его отличие от программы для ПК?
6.	В чем суть метода последовательной детализации при разра-
ботке алгоритма реализации математической модели?
7.	Как вы понимаете термин «адекватность»? В частности, что
означают понятия: адекватная модель объекта, адекватное поведение
гражданина, адекватный ответ на провокацию?
8.	Как проверяется адекватность математической модели?
9.	Как принято подразделять отдельные свойства объекта при
анализе на предмет разработки математической модели этого объекта?
10.	В чем сходство и в чем различие между входными и управ-
ляющими параметрами объекта исследования?
11.	Какие параметры объекта исследования называют выходны-
ми? Зависит ли число выходных параметров объекта от постановки
задачи исследования?
12.	Что называют возмущающими параметрами? Расскажите об
особенностях этих параметров.
13.	Как учитывают влияние возмущающих параметров на резуль-
таты моделирования?
14.	Какие объекты принято называть детерминированными, а ка-
кие стохастическими?
15.	Что принято называть математической моделью объекта?
16.	Как принято классифицировать модели по пространственным
особенностям?
17.	Как принято классифицировать объекты и их модели по осо-
бенностям их поведения во времени?
18.	Расскажите об особенностях ячеистых моделей.
19.	Какие модели называют квазистационарными?
20.	Перечислите типовой состав математического описания объ-
екта исследований.
20
21.	Перечислите наиболее общие требования к математическим
моделям объектов.
22.	Для чего используются математические модели?
23.	Какие приемы используются при разработке математических
моделей сложных объектов?
24.	Как определяется КПД объекта при последовательном соеди-
нении элементов сложного объекта?
25.	Как оцениваются свойства объекта при параллельном соеди-
нении отдельных элементов?
26.	Расскажите о соединении элементов объекта в структуры с
обратной связью. Как организуются положительная и отрицательная
обратные связи? В каких случаях используются такие структуры?
27.	Какие преимущества обеспечивает метод графов при деком-
позиции сложных объектов?
28.	Как формируется структурная матрица графа декомпозиции
сложного объекта?
29.	Как заполняется матрица видов связей графа декомпозиции
сложного объекта?
30.	Какие ограничения обязательно должны содержаться в мате-
матической модели?
31.	В чем сходство и разница между математическим и имитаци-
онным моделированием?
32.	Что принято называть характеристиками объекта моделиро-
вания? Как принято отображать такие характеристики? Как Вы пони-
маете термины «гиперпространство», «фазовое состояние», «фазовая
траектория»?
33.	В чем особенность статических характеристик объекта?
34.	Какие динамические характеристики используются для опи-
сания объектов, работающих на неустановившихся режимах?
35.	В чем различие между кривыми разгона и характеристиками
переходных процессов?
21
2. МЕТОДЫ АНАЛИЗА ПРОСТЫХ
МАТЕМАТИЧЕСКИХ МОДЕЛЕЙ
2.1.	РЕШЕНИЕ НЕЛИНЕЙНЫХ УРАВНЕНИЙ ВИДА F(x) = О
Достаточно часто математическая модель объекта может быть
представлена в виде нелинейного уравнения с одним входным пара-
метром. Для примера рассмотрим тепловые волны в полуограничен-
ном теле, возникающие при действии на его поверхности гармониче-
ского температурного возмущения по закону
0пов= Oo-cos(g)t),
где 0о _ амплитуда колебаний температуры; со = 2я/- угловая частота;
f= 1/Т; Т- период колебаний. Распределение температуры внутри тела
в этом случае описывается формулой
где а - коэффициент температуропроводности материала; х - расстоя-
ние от поверхности тела; т - текущее время [6].
Если при анализе этой математической модели ставится задача
определить, на каком расстоянии х в момент времени т = тзад амплиту-
да колебания температуры 0 будет равна 0,3 0о, то мы приходим к не-
обходимости решить следующее нелинейное уравнение, записав пре-
дыдущую формулу в виде
0о[ехр(-2а)cos(an;3a7j -х^(^2а)- 0,3] = 0 .
Математики разработали множество численных методов решения
таких уравнений, различающихся по своей эффективности с точки
зрения количества операций для получения результата или объема
оперативной памяти, занимаемой при решении. Правда, для современ-
ных ПК, отличающихся очень большим быстродействием и огромной
оперативной памятью, эти различия становятся непринципиальными.
Поэтому для общего знакомства рассмотрим только 3 простых метода:
метод сканирования, метод половинного деления и метод хорд.
1.	Метод сканирования с автоматическим уменьшением шага.
Чтобы понять суть этого метода, изобразим, как меняется функция
F(x) с ростом х (см. рис. 2.1). Отметим, что корень уравнения находит-
ся в точке пересечения кривой с осью абсцисс (именно здесь F(x) = 0).
Чтобы найти корень уравнения, от некоторой начальной точки хнач с
шагом Ах вычисляют значение функции F(x + Ах) до тех пор, пока
предыдущее и последующее значения этой функции не окажутся с
разными знаками. В такой ситуации величина шага уменьшается в не-
22
сколько раз и меняется на противоположный знак этого шага. Умень-
шение и изменение знака Дх продолжается до тех пор, пока величина
Дх не станет менее некоторого наперед заданного значения е, опреде-
ляющего точность получаемого решения. Описанный в общем алго-
ритм решения поставленной задачи более подробно представлен на
рис. 2.2.
Рис. 2.2. Блок-схема алгоритма решения уравнения F(x) = О
методом сканирования
23
2.	Метод половинного деления отрезка заключается в том, что на
оси х выделяется участок ab, содержащий корень уравнения (это можно
выявить путем сканирования с большим шагом Ах). Алгоритм решения
хорошо иллюстрируется рис. 2.3. Чтобы найти корни уравнения, вычис-
лив значения функции на концах интервала F(a) и F(b), рассчитывают
абсциссу положения середины отрезка с = (а + Ь)/2 и величину F(c). Да-
лее сопоставляют знаки F(a) и F(c). Если знаки этих величин одинаковы,
то на участке ас корня нет и тогда начало участка (точку а) переносят в
точку с и повторяют деление и сопоставление функций. Если же знаки
F(a) и F(c) противоположны (как для функции, изображенной на рис. 2.3
штриховой линией), то на отрезке ас корень есть и тогда в точку с пере-
носят правую границу интервала b и повторяют деление отрезка, вычис-
ление функции F(c) и сопоставление знаков до тех пор, пока
abs(a - b) станет меньше наперед заданной величины е. Описанный ал-
горитм наглядно отображен блок-схемой, приведенной на рис. 2.4.
F(b)
Рис. 2.3. Иллюстрация метода половинного деления
Рис. 2.4. Блок-схема алгоритма метода половинного деления
24
3.	Метод хорд также используется
для нахождения корня на [а, 6]. Прибли-
женное значение корня определяется
здесь точкой пересечения хорды, соеди-
няющей величины F(a) и F(b), с осью
абсцисс (см. рис. 2.5). При этом каждое
последующее приближение корня (или
положение точки с) рассчитывается по
итерационной формуле
Рис. 2.5. Иллюстрация
метода хорд
с = а + F(a)
b-a
+ F(b)
При этом величины а и F(a) заменяются на с и F(c), если F(c) и
F(b) имели одинаковые знаки, или же b и F(b) заменяются на с и F(c),
если F(c) и F(b) имели разные знаки. После замены цикл расчетов с и
F(c) повторяется и снова проводится анализ знаков (штриховой линией
на рисунке показан результат второй итерации, показывающий на-
сколько точка пересечения новой хорды приближается к значению
корня). Итерации прекращаются, когда абсолютное значение разницы
между с и а (или с и Ь) не станет меньше s. На рисунке 2.6 приведена
блок-схема описанного алгоритма.
Рис. 2.6. Блок-схема алгоритма метода хорд
25
Приведенные алгоритмы широко используются для решения
уравнений рассмотренного типа в различных математических пакетах,
предлагаемых на рынке: MathCAD, MathLAB, Maple и др. Для реше-
ния таких уравнений в электронной таблице MS Excel имеется специ-
альная надстройка «Подбор параметра», методику работы с которой
мы освоим на практических занятиях, решая приведенную ниже зада-
чу, сводящуюся к системе из двух нелинейных уравнений. Здесь же мы
познакомимся с одним из возможных приемов, позволяющих решать
такие системы.
Продемонстрируем использование рассмотренных алгоритмов для
решения одной из непростых задач термодинамики влажного воздуха.
Итак, состояние влажного воздуха задано параметрами ср (относи-
тельная влажность, %) и Н (удельная энтальпия влажного воздуха,
кДж/кг). Требуется определить два других параметра: t (температура
влажного воздуха, °C) и d (влагосодержание влажного воздуха,
кгвл/кгсв) [7, 8].
Изучая «Теоретические основы теплотехники» в разделе «Влаж-
ный воздух и другие газовые смеси» в результате специального анали-
за мы получили 2 уравнения, определяющие связи между всеми че-
тырьмя параметрами влажного воздуха (/, J, ср и Н):
d = 0,622 ФЛ| ,
Я-ФРн
Н = 1,005? + (2500 + 1,97?)J.
(2-1)
(2.2)
Правда, в формулу (2.1) величина t входит опосредствованно, че-
рез значение давления насыщенного водяного пара рк при температуре t.
Зная /, величину /?н можно найти по таблицам насыщения или рассчи-
тать по эмпирической формуле, приведенной в учебном пособии [9]:
(t+ 273,15)'
1000
+ 62 ln(Z + 273,15)
(2.3)
В этой формуле коэффициенты bi имеют следующие значения:
Ь_! = -7,82154; bQ= 82,8657; bx= 10,28; Ь2 = -11,4878.
С учетом этого формулу (2.1) можно представить в виде функции
от параметров (р и t:
d - 0,622
z=—1
1
1000
+ £>2ln(Z+273,15)
(2.4)
z=—1
/ +273,15
1000
+ Z>2ln(^ + 273,15)
26
Точно также формулу (2.2) путем простых алгебраических преоб-
разований представим в виде зависимости d от остальных параметров
(Яи?):
Я-1,005?
2500 +1,97?
(2.5)
Теперь мы можем решить задачу, приводя ее к виду F(x) = 0. Дей-
ствительно, разница значений d, рассчитанных по формулам (2.4) и
(2.5), равна нулю (d - d = 0). В результате получаем нелинейное урав-
нение вида F(?) = 0:
0,622
ф-ехр
+ 621п(? + 273,15)
В - ф • ехр
+ 621п(? + 273,15)
Я-1,005?
2500 +1,97?
Для решения этого нелинейного уравнения можно использовать
любой из рассмотренных нами алгоритмов. Не прибегая к написанию
программы для ПК, мы получим решение с помощью надстройки «Под-
бор параметра» пакета MS Excel, поскольку с этим пакетом хорошо зна-
ком каждый студент после изучения дисциплин «Компьютерное дело-
производство» и «Основы информатики и программирование». Фраг-
менты программы с открытыми ячейками, содержащими расчетные
формулы, приведены на рис. 2.7. Фрагмент первого шага решения и его
результаты приведены на рис. 2.8 (значение ? в ячейке Е56, значение d -
ячейке G53, где приведена полусумма значений d, рассчитанных по
формулам (2.4) и (2.5) при найденной температуре ?). Для сравнения
приведем результаты определения названных параметров по известной
Н- d-диаграмме Рамзина: при Я = 80 кДж/кг и ф = 50 % (при расчетах ф
берется в долях ф = 0,5) d = 18 г/кг, ? = 35 °C. Как это можно видеть по
результатам расчетов, они прекрасно совпадают с приведенными циф-
рами, хотя результаты графического определения содержат естественно
большую погрешность, чем результаты расчетов. Заметим, что при рас-
четах величину d мы получаем в системе СИ (кг/кг).
Менее точно, но вполне сопоставимо с результатами, полученны-
ми по Я - J-диаграмме, нашу задачу можно решить и по-другому,
прибегая к графической интерпретации формул (2.4) и (2.5). Для этого,
задавая ряд значений температур, например с шагом через 5 °C, будем
рассчитывать величины d по формулам (2.4) и (2.5) и отображать их на
графике d=flt). Точка пересечения двух кривых и будет соответство-
вать решению системы уравнений (поскольку здесь d = d). Абсцисса
27
этой точки даст нам значение t, а ордината - величину d. Выполнять
описанные расчеты очень удобно с помощью MS Excel. Ниже на
рис. 2.9 и 2.10 показаны приемы создания программы и результаты
таких расчетов в виде соответствующего графика, иллюстрирующего
данные табличных расчетов.
с Microsoft Excel - Задача 1 с вл воз
®Файл Правка Вид Вставка Формат Сервис Дэн-ые Qw Огравка______________________________________________
100% у о .|| Arial Суг
F55	~_г1	-|^E<P(SF»371000/(E56^273,15)]-»-SFH^ff5‘((E5^-273,15yiOOO)-»->Ft6-Ln(E56<73,15))
В	С	D Е F 1 G | К | I ] J
162
60 Определи 0016391
05 Определи 30
153	Задаю Н>
54	Задано фи=
"351	Расчёт Pw=f(t| по формуле (2.3)
I ffi	t нач
67
53	Расчет d по формул* (2.4) 4е
I 69	Расчет d по формуле (2.5) й=
60	' F(i>
HJiXWE
30
0013302
О,СП 948
0006178
L Microsoft Excel - Задача 1 с вл воз
®фаил Оэаека Вид Встда^а Формат Сервис Дам-ью Q<no £п
□	-п с. . % £ л *И1
F58 /d »| =(0j62y$DS4*F55MIMS-ffD6rF56)) _______________
B~f С | D | Е | F | G Н
62
60
61
Задано Н=
Задаче фи=
Расчёт Рн=ЭД по фермуле (23)
t нач
Расчет d по формуле (2.4) 4=
Расчет d по формуле (2.5) d=
F(x)-
60 Определи' 0,016391
0.5 Определи1 30
0004242
30
I 00133021
001946
0006178
F60
J ! =(ФЕ653-1.ОО5-Е5В)Г(25ОО*107>Е56)	'
I С | D Т I F I G [
Задано О	60 Определи	0016391
Задано фи=	05 Определи	30
Расчёт Рй=ф) по формуле (2.3)	0ДО4242
1 нач	30
Расист dno формул* (2.4) 4е
Расчет d по формуле (2.5)
Р(я)=
0013302
I QJ01948I
ОД06178
| =F59T68
С	D
53
54
5Б
66
57
Я
S
ДР I
Задаю №
Задано фи=
Расчёт по формуле (2.3)
t нач
Расчет d по формуле (2.4) 4s
Расчет d по формуле (2.5) d=
F(x)=
Ш Определи 0Д16391
05 Определи 30
ОРО4242
30
0013302
001948
I ojocsizei
В
G
Рис. 2.7. Фрагменты таблицы Excel с открытыми окнами
расчетных формул
28
Microsoft E>cel - Задача 1 с вл воз
DfiaAn qpbu аж Вст«*а Фадазт едв< 1Х<> Охмвд
» <№ГЮ
С D
к
I
N
Задмфм*
ПОпрмаиод d * OpiWrt
ОЗОодмытч 1 *
ОДОЙ IQ
ГЬдбср парам
51
одпытв
Microccft E>cei - Задача 1 с гл во*
в®йп В^ека frvi йсгдака Форм* Сервис Д»*м Оно Qrpsf
о огя«й х *•<* • - - •’ ж л м n:tlw
го «I « -го го
"7---5”“^----Е-------------Е *	Т~5------
вЭ;
Йми №	ЮО^щиь d. * QJJ1ZO
Л^*н> ф**“	Он^йДФЛЛъ i t> □ЧДМХ
Рк** ВДОДод +w*« рл 0Ш67Г»
1одч ДОЖ
Родчп d nt	*
Pww d од *ютм« 1ДО1
FW-
ОДОИСПЕ
nntrei
Рис. 2.8. Фрагменты таблицы Excel с диалоговым окном надстройки
«Подбор параметра» и результатом расчета
- Гн r<xoft Excel Зэдачл 1 с ол гиз
Ойвйл Вши Ь® Вставь Формат Сервис Д»*«е 0*о фтхва
D	X	-	-  < Ж A il II В4 1004 » (В. А>ШО
ftl2 _d •	ЫЧГ15ТСМ2*2П1€>П1ХА»**Т^1^?Э.^)
; *	~ »- 1 С О	« . и i
1	РашиIIи	од gy щмнн| одрод«фа* МОДОД^О одвдм
2	.	КмффмфМы феодjjm (J)
?	ИД**	W	Ч’г	>АИМ
4 .	2-д«ф>*	ОД	ЬДО	ЮЖ?
О	0.1013	ЬЛ>»	wje
с	К> -нда
? .
.9 .	йдмм»	80 Снодламгъ d
>	Змм ♦*"	05 Онршмтъ I
W.
11 1	fw 14*} <нт	4^4
01	< ftjHHWJwl одйутош «доло1№ ^дазйв
U подавит 0ДМПИБ 1ДШ21«И ОДЕБ
ricrozoft Lxccl ’ Задача 1 с од зоз
Q файл Омска м встде Формат Эдис Дамм Оно СТ
о * а «й ал? х ъж* • < х * и л а 4
______________________________________________________
а в ______________с 1	в f o'
Р*Щ»М *мни ОД рушодим ОДОДВГфф* одпвп воыуод
К< IЦI fill Ilu foytfjMiP)
1ЛДАДА И-	М Ц7.O1U
Хдм 9^	05	bflr МР
Бфшгр«одстдмм«од О* О JOO Ц1>» WJB
11,46^0
Змм^
И|м>г
00 Спрвдвжти d
0> Споддшть I
АходИЯ йя^м)4од (3,11 ОДМ
и олтзэтадиВЯРД ашлимт *дзж
v opt№75 6ДОШ4Ц дедок» одаь
пи jJ ’ *до№1да;л1^ря&.1УА1а}____________________
—. *	— в 1 с	О 1 Е f о
1	ОДодв яМ** • одЦОДиа MfMMFtpui 1Ммп1 одщяв
2	1й»фф«одты ферм/лы fty
W	-7521Ы
оз га> «да?
0,013	1020
гаг и.мж
ВО ОфМ*** 4
ОД Оподдеоть L
Рн ОД (ДЛ б ОД (2Л <-ОД -£Д?1 _OKd
N ОДО1222ЛО OXBWb<dTBCTg*f| wm
П 0ДОЮТК 0ДЮ1« tJBEiljr ода»
Рис. 2.9. Фрагменты программы расчетов d в пакете Excel
с открытыми окнами формул для каждой ячейки
О
29
Рис. 2.10. График результатов расчета
2.2.	РЕШЕНИЕ СИСТЕМЫ ЛИНЕЙНЫХ
УРАВНЕНИЙ МЕТОДОМ ГАУССА
Очень часто, особенно когда математическая модель составляется
из балансовых уравнений для стационарного объекта или установив-
шегося режима, описание связей между параметрами объекта пред-
ставляется в виде замкнутой системы линейных уравнений. Как пра-
вило, эту систему записывают в следующей форме
r«n%i + ли %2 + «13 *з + ...+ ain хп = Ь\,
«21 %1 + «22 %2+ «23 Хз + ...+ «2л Хп = Ь^,
«л1 ^1 “г ««2 Х% “г «лЗ *3 ' • ’ • ' «лл
хп= Ь
П9
где «у - коэффициенты при неизвестных Xj (первый индекс — номер
строки, второй индекс - номер столбца); bi (i = 1, 2, ..., п) - свободные
члены уравнений. Приведенную систему можно представить и в со-
кращенном виде
или в матричной форме
А хХ=В,
где А - матрица коэффициентов «^; X - вектор-столбец неизвестных
В - столбец свободных членов bj (i - номер строки; j - номер столбца):
а11> а12> •••
а21> а22> •" а2п
• • •
^л1 ’ @п2> ••• ®пп
30
Система называется несовместимой, если у нее нет ни одного
решения. В противном случае система совместима. Совместимую сис-
тему называют определенной, если она имеет единственное решение и
неопределенной, если имеется несколько решений. В теории линейных
систем сформулированы условия, необходимые для совместимости
системы и для ее определенности, и этот анализ не требует решения
системы. Метод Гаусса позволяет решить задачу, одновременно он же
дает ответы и на вопрос о совместимости и определенности.
И еще одно замечание. Доказано, что решение систем не изме-
нится, если поменять местами любые две строки системы; если по-
членно умножить правую и левую часть любого уравнения на любое
число, не равное нулю, и если любое уравнение заменить суммой или
разностью его с другим уравнением системы. Именно эти 3 преобразо-
вания используются в методе Гаусса.
Основная идея этого метода заключается в последовательном ис-
ключении неизвестных, начиная с первого, до тех пор, пока не полу-
чим уравнение с одной неизвестной хпп. Чтобы лучше понять суть ме-
тода Гаусса, рассмотрим систему из трех уравнений
ЛцХ1+ «12 %2+ «13 Х3 = Ь1,
Л «21 Xi + «22 *2 + «23 ХЗ = Ь2,
«31 *1 + «32 х2 + «„3 Х3 = Ь3;
(2.6)
(2.7)
(2.8)
Чтобы исключить X] из уравнения (2.6), его умножают на дробь
«21/«ц и вычитают из уравнения (2.7). При этом, если «ц = 0, то урав-
нения меняют местами, так, чтобы первым было уравнение с «п 0.
Выполняем описанные действия:
(«21/«11) X «цХ1+ («21/«11) X «12 Х2+ («21/«11) X «1з Х3= («21/«11) X ЬЪ
(«21 Х1 + «22 Х2+ «23 Х3) - («21/«11) X «цХ1 - («21/«11) X в12 Х2~
~ («21/«11) X «13 Х3 = Z?2-(«21/«ll) X
Слагаемые, содержащие Xi, взаимно уничтожаются, а слагаемые с
х2 и х3 группируем следующим образом
[«22 - («21/«11)Х«12] Х2 + [«23 - («21/«11)Х«13] Х3 = \Ь2 - («21/«11)х&1].
Для упрощения обозначим выражения, стоящие в квадратных
скобках, через А22, А23 и р2, соответственно:
А22Х2 + А23Х3 — р2-
(2.9)
Итак, формулу (2.7) мы заменили на формулу (2.9), не содержа-
щую Х1.
Точно так же исключим Xi из уравнения (2.8):
(«31/«11)х«цХ1 + («31/«11)х«12 х2+ («31/«11)Х«13 Х3= («31/«ll)xZ>i,
31
(«31 xi+1?32 х2+ а33 х3) - (а^/а^ха^Х! - (а31/ап)ха12 х2 -
- (а31/а31)ха13 х3 = b3-(a31la^xbx.
Группируем подобные члены:
[«32- («3i/«ii)x«i2] х2+ [а33-(а31/ап)ха13] х3= [63 - (a3i/«n)xZ>i]
и переписываем полученное уравнение в сокращенной форме
^32^2 + Л33Х3 — Рз.	(2.10)
Теперь вместо уравнения (2.8) мы имеем уравнение (2.10), экви-
валентное замененному. В итоге осталось 2 уравнения с неизвестными
х2 и х3.
По тому же алгоритму исключим неизвестную х2 из уравнений
(2.9) и (2.10). Переписываем (2.9), умножив его почленно на дробь
(А 32/А22)\
(А32/А22)хА22х2 + (А32/А22)хА23х3 = (А32/А22)хр2.
Вычитаем почленно из (2.10):
А32х2 + А33х3 — (А32/А22)хА22х2 — (А32/А22)хА23х3 = р3 — (А32/А22)хр2.
После исключения х2 группируем подобные члены:
[Лзз - (А32/А22)хА23]х3 - р3 - (А32/А22)хр2,
(2.П)
откуда находим:
Формулу (2.11) для простоты записи представим так:
#ззГз = Уз
и отметим, что исходная система (2.6) - (2.8) превращена нами в тре-
угольную:
+ Я12 х2 + 6?13 Х3 = Ь\,
5 0 + А22х2 + А23х3 = р2;
ч 0	0	В33х3 = у3.
Значения неизвестных х2 и х3 теперь легко находятся, если рас-
чет вести снизу вверх, рассчитав сначала х3, затем, подставив получен-
ное значение в формулу (2.9) или (2.10), найдем х2, после чего, под-
ставляя значения х3 и х2 в формулу (2.6), находим %].
В процессе последовательных исключений неизвестных возмож-
на ситуация, когда все коэффициенты преобразованного уравнения
32
превратились в нули, а свободный член не равен нулю. Это означает,
что система не совместима, т.е. система решения не имеет. Если на
каком-то шаге исключений все коэффициенты становятся нулевыми и
свободный член тоже равен нулю, это означает, что анализируемое
уравнение является линейной комбинацией другого уравнения систе-
мы и тогда анализируемое уравнение выпадает из системы, она стано-
вится не замкнутой. В такой ситуации система является совместимой,
но неопределенной, т.е. имеет несколько решений.
Благодаря достаточно простому итеративному алгоритму, для
решения систем методом Гаусса разработаны и опубликованы стан-
дартные процедуры на языках программирования высокого уровня
(Паскаль, СИ-н- и др.). В настоящее время такие процедуры включены
в популярные интегрированные пакеты прикладных программ, пред-
лагаемые на рынке математического обеспечения ПК (MS Excel,
MathCAD, MatLAB, Maple и др.). На практических занятиях мы осво-
им методику решения таких систем.
2.3.	ИТЕРАЦИОННЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ
ЛИНЕЙНЫХ УРАВНЕНИЙ
Алгебраический метод Гаусса в отдельных случаях не может
обеспечить приемлемой точности решения задачи. Дело в том, что при
определенных значениях величин и bi задача становится «плохо
обусловленной»: даже при небольших, лежащих в пределах приемле-
мой точности, изменениях любой из названных величин результаты
решения изменяются очень существенно. Для примера возьмем систе-
му с двумя неизвестными [2]:
Г300х1 + 400 х2 = 700;
[1 OOxi + 133 х2 = 233.
Точное решение такой системы Xi = 1, х2 = 1. Представим, что
свободные члены 700 и 233 были определены нами в результате экспе-
риментальных исследований и записывались с дисплея цифрового
прибора, где третий знак всегда дает погрешность ±1. Тогда с такой же
вероятностью число 233 может быть заменено на 232. При этом реше-
ние системы будет совершенно другим: Xi = —3, х2 = 4. В итоге получа-
ем: при погрешности коэффициента b2 А62 = 1/233-100 = 0,43% по-
грешность неизвестной х2 получается Дх2= (4 - 1)/1 • 100 = 300%.
Численные итерационные методы позволяют находить прибли-
женное решение отдельных плохо обустроенных систем с заранее за-
данной точностью и именно поэтому широко используются на практи-
ке, хотя, конечно же, требуют гораздо больших вычислительных ре-
сурсов (больше времени для расчетов, иногда больше оперативной
памяти и т.п.).
33
Чтобы понять основную идею, заложенную в рассматриваемый
подход, запишем сначала простейшую систему:
«цХ1 + «12 %2 + «13 -*3 = Ь\,
< «21 *1 + «22 *2 + «23 *3 = Ьъ
«31 *1 + «32 *2 + «33 *3 = Ь3
и перепишем ее следующим образом, выразив по очереди каждую из
неизвестных:
rxi = -«12/«п Л-2 - «1з/«11 *3 + Z>i/«n;
< Х2 = ”«21/«22 *1 - «2з/«22 *3 + Ьо/а22,
Лз = -«31/«33 *1 - «2з/«зз Х2 + Ь3/а33.
(2.12)
(2.13)
(2.14)
В начальном, нулевом приближении зададим некоторые значения
неизвестных xi, х3, например так: Xi = х2 = х3 = 0 или xi = х2 = х3 = 1
или как-либо по-другому. Далее последовательно, используя формулы
(2.12), (2.13) и (2.14), рассчитывают величины xi, х2 и хз в первом при-
ближении, затем, подставляя полученные значения в правые части тех
же формул, находят Xi, х2 и хз во втором приближении. Итерации по-
вторяются до тех пор, пока каждая из величин xi, х2 и хз в последую-
щем приближении будет отличаться от значений предыдущего при-
ближения не больше, чем на величину наперед заданной точности е.
Если в общем случае исходную систему уравнений записать так
(напомним: i — номер строки, j - номер столбца):
^aijxj=bi, 1 = 1, 2,п
7=1
то переписанная система с выделенными неизвестными в общем виде
—, /= 1,2, i*j.
(2.15)
При этом реализуется следующий итерационный процесс: каждое
последующее {к + 1)-е приближение рассчитывается по формуле
(2.15), причем в правую часть формулы ставятся величины х7 преды-
дущего, £-го приближения. В этом и состоит суть метода простых ите-
раций.
Доказано, что процесс итерации сходится к точному решению,
если по абсолютному значению любой множитель перед х7 меньше или
равняется единице (K«j/««)| < 1) и хотя бы один из них удовлетворял
только неравенству.
34
Чтобы система удовлетворяла приведенному условию сходимо-
сти, необходимо, чтобы ее диагональные элементы преобладали по
величине над другими. Для этого во многих случаях приходится про-
водить известные преобразования системы, не влияющие на результат
решения, в частности это перестановка уравнений и замена одного из
уравнений суммой или разностью его с другим уравнением системы.
Разберемся с этим на примере. Имеется система линейных неод-
нородных уравнений:
-2,8xi + Х2 + 4х3 = 60;
5 IOX1-X2+ 8х3= 10;
-Xi + 2х2 - 0,6хз= 60.
Поскольку «13/^14 = 4/(-2,8) = -1,43, и это по абсолютному значе-
нию больше 1,0, то в представленном виде итеративный процесс будет
расходящимся. Обратим внимание: диагональный коэффициент в пер-
вой строке по абсолютному значению меньше, чем сумма двух других
коэффициентов (2,8 < 1 + 4).
Переставим уравнения в системе следующим образом:
IOX1-X2+ 8x3= 10;
X -Xi + 2x2 - 0,6хз= 60;
-2,8xi + х2 + 4х3 = 60.
Теперь диагональные элементы преобладают в каждой строке:
10>1 + 8;2>1+ 0,6; 4 > 2,8 + 1. Поэтому соблюдаются все условия
сходимости: |aiy а1Д| = 1/10 < 1,0; |а13/ au| = 8/10 < 1,0; |a2,i/	= 1/2 < 1,0;
|^2,У #2,2!= 0,6/2 < 1,0; |а3д/	| = 2,8/4 < 1,0; |я3>2/ <*з,з1= 2,8/4 < 1,0.
Итерационный процесс сходится быстрее, если при расчете по-
следующего приближения в правую часть формулы (2.15) при i > 1
подставлять значения xi, Х2 и т.д., рассчитанные уже в последующем
приближении в предыдущих расчетах. Например, при расчете х3 по
формуле (2.14) в правую часть ее подставлять значения xi и хг, рассчи-
танные перед этим по формулам (2.12) и (2.13). Такой метод называют
методом Гаусса-Зейделя и итерационная формула в общем виде здесь
представляется так:
Достоинством итерационных методов является простота алгорит-
ма и минимальная память для хранения коэффициентов матрицы. Осо-
бенно удобно применять этот метод для сильно разреженных матриц,
где многие коэффициенты а у равны нулю.
35
2.4.	МЕТОДЫ РЕШЕНИЯ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИИ
Рассмотрим сначала решение систем, которые принято называть
«квазилинейными». Так называют системы уравнений, у которых ве-
личины коэффициентов перед неизвестными не являются постоянны-
ми, а зависят от значений неизвестных хц Х2, хз,..., т.е.
&ij	^2, -*3, • • •)*
Такие системы очень широко встречаются при разработке мате-
матических моделей в теплоэнергетике, где многие теплофизические
характеристики зависят от величины температур, которые и выступа-
ют в качестве искомых неизвестных. Для примера запишем математи-
ческую модель для анализа процесса стационарной теплопередачи че-
рез плоскую стенку (см. рис. 2.11), не прибегая к обычным упрощаю-
щим предпосылкам, т.е. учитывая, что коэффициент теплопроводности
X и коэффициенты теплоотдачи ai и 0С2 прямо или косвенно зависят от
величины температуры. Основу модели составляют следующие тепло-
балансовые уравнения
q СС1(/Ж1 41),
<q = X(4i - 4г)/8;
q — 0,2(42 - АкД
где ai =/1(/ж1, 41), a2=/2(4<2, 4г), =/з(4ь 4з). Величины ?жЬ ^и5 счи-
таются заданными, а в качестве неизвестных выступают параметры q,
4i и 42*
Квазистатическая система формально имеет прежний вид:
п
7 = 1,2,...,«.
j=i
Рис. 2.11. Теплопередача через плоскую неограниченную стенку
36
Метод простых итераций здесь состоит в том, что перед каждым
-	£_i_1	£+1
расчетом последующего приближения неизвестных х ,	х2?
х3,... вычисляются значения коэффициентов а у по величинам
хх , х2? х3?..., полученным при расчетах предыдущего приближения.
В результате на этом этапе получается обычная замкнутая система
уравнений с постоянными коэффициентами, решая которую находят
значения неизвестных /+1 х2/+1 х3,.... Далее повторяют описан-
ныи алгоритм до тех пор, пока разница между последующими и пре-
дыдущими приближениями по каждой из неизвестных не станет
меньше заданных величин щ, 82 и т.д.
Второй подход более сложен, но обеспечивает, как правило, более
быструю сходимость итерационного процесса, а в некоторых случаях
позволяет решить задачу, неразрешимую другими способами. Для
лучшего понимания рассмотрим простейшую квазилинейную систему
из двух уравнений
(2.16)
Предположим, что найдены значения неизвестных х} и х2 пре-
дыдущего приближения (начальное приближение при к = 0 принима-
ется так же, как это описано ранее). Тогда для последующего прибли-
жения неизвестные можно представить так
2 >
(2.18)
(2.19)
где Ajq и Дх2 - уточняющие добавки для перехода от предыдущего к
последующему приближению. В рассматриваемом методе (его назы-
вают методом Ньютона) записывается система уравнений для этих
уточняющих добавок. Чтобы составить такую систему, разложим
функцию Яц(Х1, Х2) в ряд Тейлора в точке jq и ограничимся первы-
ми двумя членами ряда:
ия +
«ц -г
ЙХ1
—dx2.
х2
Переходя к конечным приращениям Axj и Ах2 (dx заменяется на
Ах), перепишем приведенную формулу:
37
Совершенно аналогично можем получить разложение коэффици-
ента <Я12<Л1, х2):
Запишем теперь первое уравнение системы (2.16) для последую-
щего приближения, учитывая формулы (2.18) и (2.19):
а11( *1 +^xi)+	<212 (
+ Лг2) = ^,
и перепишем ее, подставляя приближенные значения коэффициентов
fc+l fc+1
аи и а\2 > приведенные выше:
Раскроем скобки и сгруппируем подобные члены, отбрасывая, как
бесконечно малые второго порядка слагаемые, содержащие произве-
дения ( Дх?! х Axj ), ( Лх2 х Лх2 ) и (Лх1 х Лх2 ):
Перепишем это уравнение, вынося за скобки величины ЛХ] и
Лх2 и перенося вправо слагаемые, содержащие коэффициенты преды-
дущего приближения:
к ыи > , и W..
L41 1 I	Ли П	Ло
dx2 ,
= Ь]—кацкх1 —kdi2kx2.
к .	ц12 *
«12 -I
dx>
^а12 k
Лп
dx2
1
•2 “
38
Совершенно аналогично получим вместо (2.17) следующее урав-
нение
В итоге мы получили систему из двух уравнений, которую пред-
ставим в краткой форме, обозначая содержимое в круглых скобках
перед Axj и Лх2 через правые части уравнений через В\ и В2, а
Axj и Дх2 через и
Решив систему обычным способом, находим последующие при-
ближения неизвестных. Продолжая итерации, мы будем уменьшать
уточняющие добавки Axj и Дх2 до тех пор, пока они не превратятся в
нули (а точнее — не станут меньше наперед заданной малой величины с).
В заключение отметим, что оба рассмотренные метода в опреде-
ленных случаях могут быть использованы и для решения нелинейных
алгебраических уравнений - таких уравнений, в которых неизвестные xi,
х2, ... входят в уравнения системы в виде некоторых нелинейных алгеб-
раических зависимостей или содержат трансцендентные функции.
Для примера запишем такую замкнутую нелинейную систему:
ГлцХ12 + tfi2 x2=Z>i;
|a2i xi + я22 х22 = Z>2,	(2.20)
методам решения которой при заданных постоянных численных зна-
чениях коэффициентов мы посвятим дальнейшее теоретическое
рассмотрение и одно из ближайших практических занятий за компью-
тером.
2.5.	ОСОБЕННОСТИ РЕШЕНИЯ СИСТЕМ
НЕЛИНЕЙНЫХ УРАВНЕНИЙ
Решать такие системы часто бывает необходимым при поиске
безусловного экстремума в задачах оптимизации, при использовании
явных схем решения систем обыкновенных дифференциальных урав-
нений и во многих других случаях.
В общем поставленная задача может быть записана так:
39
/i(xbx2, ...,х„) = 0;
Л(хьх2, ...,х„) = 0;
(2.21)
Л(хьх2, ...,х„) = 0,
где /1 -fn - некоторые нелинейные функции. Решение «в квадратурах»
(т.е. с получением алгебраических формул для расчета любой неиз-
вестной X/) таких систем, как правило, невозможно. Поэтому решают
их численным способом.
Известны несколько методов решения таких нелинейных систем.
Самый простой подход — это метод простых итераций, знакомый нам
по решениям одномерных задач вида х =fix), а также рассмотренный
выше для решения системы линейных уравнений. Чтобы реализовать
этот метод, систему (2.21) путем простых алгебраических преобразо-
ваний приводим к виду:
Х1 Ф1(Х1, х2, ..., х„),
Х2 = ф2(Х1, х2, ..., х„);	(2.22)
Теперь, задавшись начальными приближенными значениями не-
известных °х1,°х2 , ...,°хи (о том, как определить эти величины, пого-
ворим ниже), можем составить следующие итерационные формулы:
Xi =Ф1(ах1,^Х2,	*хл);
^2 — Фг( *4» ^2> •••?
которые последовательно реализуются, позволяя рассчитать любые
- £+1 £+1
последующие приближения	х2 и т.д., если известны значения
кхх,кх2 и т.д. предыдущего приближения. После каждой реализации
очередной формулы принято рассчитывать абсолютное значение раз-
ности между последующим и предыдущим приближением
Лг- = abs(*+1xf-\-) и выбирать среди них максимальное значение
тахЛ;, i = 1,2, ..., п. Когда это максимальное значение среди Л/ станет
меньше, чем задаваемая точность s (например, £ = 0,001), итерации
jt+i &+1
прекращаются и в качестве решения принимают значения Xj, х2,
..., хп последней итерации. Мы описали алгоритм так называемого
параллельного итерирования.
40
Приближенные значения для начала итераций 0x1,°x2j
можно найти графическим способом. Для этого правую часть формул
(2.22) обозначим через yi, у2 и т.д., запишем таблицу с рядом положи-
тельных значений xi и такую же таблицу с значениями х2 (и так же для
всех остальных х(), а в следующие таблицы будем заносить расчетные
значения yt. Это позволяет (правда, это только при малой размерности
системы) в координатах xi - х2 (для примера) изобразить по соответст-
вующим значениям уг и у2 точки, объединяя которые найдем график
каждой из зависимостей, а приближенные решения будут лежать в
точках пересечения этих графиков.
Кроме определения начальных приближений решения, для систе-
мы (2.22) необходимо установить, в окрестностях приближенного ре-
шения эта система обладает или нет свойством сходимости? Вспом-
ним, что для уравнения вида х = ф(х) таким условием является нера-
венство
<1
dx
Без доказательства (но оно существует в учебниках по численным
методам) для системы (2.22) запишем аналогичные условия
max max
J i
d<Pt
dx;
«X
Другими словами, каждая из частных производных функций ф, не
должна превышать величины 1,0.
Ускорить процесс решения (уменьшить число итераций) можно,
применив последовательное итерирование, суть такого изменения со-
стоит в том, что при каждом итеративном расчете каждой неизвестной
с большим номером в уравнения вместо предыдущего значения в пра-
вой части формул (2.22) ставят определенные выше значения после-
дующего приближения:
&+1
*1 =Ф1( *1, х2,.
к+1
*2 =<р2(
Аг 4-1 v _	/Ан-1	&4-1 k \.
^3 Фз( ^1’	’ ’’’’
i+1x„ = <p„(*+V+1x2,.... Ч) •
Этот уже знакомый нам алгоритм называют методом Зейделя.
Чтобы освоить на практике эти методы с помощью пакета MS Ex-
cel найдем решения одной-двух систем нелинейных уравнений. Пер-
41
вым примером будет решение приведенной выше нелинейной системы
(2.20), представленное на рис. 2.12. Приведенная на нем методика ре-
шения достаточно проста и наглядна, не требующая дополнительных
комментариев.
файл драма Вид Вставка Федат Сдоис Дажые (>мо Орема
а у в 4л л ’ л z я н -
Ariel Суг	-10 - Ж Z Н
НЮ ▼ fit =ЕСЛИ(Е10>0 JD5; "Продолжай итерации’; “Конец итераций-)
Ё I F 1 G Г"-Н I 1
1
2
Э
У
6
7
8
9
10
11
12
13
14
15
1Ь
17
18
19
20
НЕЛИНЕЙНАЯ СИСТЕМА
апх
Коэффициенты при х
а11= 4	а21= -2
а 12= 3	а22= 5
Свободные члены
Ь1= 4
Ь2= 6
Преобразуем уравнения
Приближения
1 х1 пред х2 пред
1
x1=((b1-a1Zx2)/a1iyO0
к2-((Ь2 а21*х1)/о22)Ю£
х1 послед х2 послед dell
0,5 1,204911
del? Условие выхода
0,5 0,20491итерации
х1 пред х2 пред
0,3 1,264911
а
х1 послед х2 послед dell
0,226532 1,133216 0,273460
а	а
del?
О ДЗ1393 Продолжай итерации
3 х1 пред х2 пред х1 послед х2 послед dell deQ
и/лы/ i,iui2ib ид±А41 i.ijuubi u.iuwi и Д4/1 bbi 1родолжай итерации
4 х1 пред х2 пред х1 послед х2 последнем del?
0,335541 1,136051 0384658 1,155083 0049116 0019032 Конец итераций
2
Охкрьиа лпсйкл и раичсдаий фиржулий ддм л! ииыдсдридцсс
РЮ » fit =((ЩЗ-$Е$4*С10)/ВЕ*ЗуО,5
“А I В I Cl D I Ё Г
I G I H
6
7_
8
9
10
11
Преобразуем уравнения
Приближения
1 х1 пред х2 пред
1
х1=((Ь1-а12*х2уа11уО0
х2=((Ь2-а2Гх1)/а22уО0
х1 послед х2 послед dell
ll	0.б] 1264911
del? Условие выхода i
0,5 0264911 Продол
>у-
й итерации
Открыта ячейка с расчетной формулой для х2 последующее
ЕЮ
Л =((HW4Gn*B10)/$G»4y0.5
Преобразуем уравнения
D I Е |
х1=((Ь1-а12*х2уа11)Ч),5
Х2=	82 Гх 1 у aZ2fll ,Ь
G
н
11
Приближения
1 х1 пред х2 пред
0,5
х1 послед х2 послед dell del? Условие выхода I
0,5 0264911 Продолжай итераци^
I
Рис. 2.12. Решение системы нелинейных уравнений в табличном
процессоре MS Excel с фрагментами открытых ячеек для расчета
последующих значений неизвестных X] и
42
D31	▼	=((В31*(СЗИ5>1У2УОД
А В C D | Е F | G | H
23 НЕЛИНЕЙНАЯ неоднородная система
24 2Xj3 - Х1Х2 - 5Xi+1=0	Преобразуем уравнения
2Г Xi+3lgXi • х2*=О	Х|=КОРЕНЬ(ХГ(Х2+6)-1)/2)
26	X2=KOPEHb(Xi<rigXi)
I  I-! - «	II*	Т	*	«
27
26 Применяем метод Зейделя: при расчете Xj послед используем рассчитанное Х1
29 Приближения
х2 пред xi послед х2 послед dell
х1 пред
3
х1 пред
3,162278
х1 пред
3,276296
х1 пред
3,359781
х1 пред
3,410676
х1 пред
3,441431

deQ Условие
0,10508 Продол?
х2 пред х1 послед х2 послед dell deQ
2,10508 3,276296 2,196007 0,114019 О Д90926 Продол?
х2 пред xi послед х2 послед dell deQ
2,196007 3,359781 2,222322 0,083464 0Д26Э15 Продол?
х2 пред xi послед х2 послед dell deQ
2,222322 3,410676 2,238124 0,050895 ОД15802 Продол?
х2 пред х1 послед х2 послед dell deQ
2,238124 3.441431 2247587 0,030755’ О £09464 Продол?
х2 пред xi послед х2 послед dell deQ
2,247587 3,459918 2,253246 0,018488 0Д05658 Конец hi
Фрагмент программы с опфытым оптом Х2 последующее
Е31
| А | ЁГ"
29 Приближения
30 1 х1 пред
зГ]
• « *
fi. =(D31+31QG10(D31))«0.5
I С | D | E ~I
F | G | H
x2 пред х1 послед x2 послед dell deQ Условие
2 3.16227ВГТ15923П 0.162278 0,159231 Прадоли
Фрагмент программы с отхрытым окном "Условие выхода*
К31	▼ А =ЕСЛИ(Е31 >0,02; "Продолжай итерации”;"Конец итераций*
I А | В | С [	Ё | F |	6	| Н
29 Приближения
ЯП |1	х1 пред Y? пряд «1 ПЛГ Пвд Y? ЛАГПАД doll КяО Угплйив
Щ	3	2 3,162278 2,10508 0,162278 0.1050б1Продшп
Рис. 2.13. Решение системы нелинейных уравнений итерационным
методом Зейделя в MS Excel и фрагменты таблицы с открытыми
расчетными формулами
Другой пример: найти действительные и неотрицательные корни
системы уравнений
43
= 0.
В соответствии с изложенным выше, приводим уравнения систе-
мы к виду (2.22):
рч = 7[Х1(х2+5)-1]/2 ;
I х2 +3^%! .
Далее для расчетов будем использовать итерационную схему
Зейделя, в частности при вычислении величины последующего при-
ближения неизвестной х2 будем использовать значение %], уже полу-
ченное в последующем приближении. Программа расчетов легко фор-
мируется путем простого копирования расчета второго приближения в
последующие строки. Результаты расчетов приведены на рис. 2.13.
2.6.	РЕШЕНИЕ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ
МЕТОДОМ «ПРИСТРЕЛКИ ПО ЦЕЛИ»
С этим весьма специфическим методом познакомимся на примере
решения задачи по определению параметров перегретого водяного
пара, которая очень широко встречается в практической работе любого
теплоэнергетика и в докомпьютерную эру решалась обычно с помо-
щью известной h - ^-диаграммы. Сегодня же для расчета параметров
перегретого пара в технической литературе [9] рекомендуются упро-
щенные уравнения, обеспечивающие весьма высокую точность расче-
тов, превышающую определения по h - х-диаграмме.
Уравнение состояния для перегретого пара, определяющее связь
между давлением р, температурой t и удельным объемом у приводится
в виде
Ат
у = —- + В + Ср,	(2.23)
Р
где параметр т = (t + 273,15)/1 ООО, а величины В и С вычисляются как
следующие функции этого параметра т:
В =	+ 62т2 + Ьз/(х - b4)2, С = со + С1Т-8 + с2т-14.
Значения коэффициентов bi и с( рассчитаны на основании сведе-
ний о сжимаемости перегретого пара и будут приведены нами позже.
Для расчетов энтальпии h пара там же рекомендуется следующая
зависимость
h = /?о + (ВъР + С^?2/2)-1000,
(2.24)
44
в которой величина /?о и коэффициенты Bh и С* также рассчитываются
по сложным нелинейным зависимостям:
Ло = <Яо +	+ ^2^2+ «З’ЬэТ,
Bh = b0 + 3 Z>2t2 + 3 йз/(т - Z>4)2 + 2 Ь3 - brf,
C*= Cq+ 9 cit"8+ 15-сгт-14.
Для расчетов энтропии перегретого пара приводится следующая
формула
s = sq + Bs'P + Cs-p2/l - TMn(lOOOp),	(2.25)
где для расчета величин и Cs следует использовать следующие
формулы
= <яг1пт + 2л2т - аз/т + а^ЮОО,
Bs = -bi + 2Z>2t-3 + 263/(т - b4)3,
Cs = 8cit-9+ 14 c2t-13.
Заметим, что все приведенные формулы содержат те же самые
константы bi и с, и дополнительно еще константы ai9 значения которых
выписаны из [9] и приводятся в табл. 2.1.
В приведенном математическом описании значения параметров
перегретого пара задаются или получаются в следующих размерно-
стях:
-	давление перегретого пара р - в МПа;
-	удельный объем пара у - в м3/кг;
-	температура перегретого пара t - в °C;
-	удельная энтальпия пара h - в кДж/кг;
-	удельная энтропия пара - в кДж/(кг-К),
-	газовая постоянная R — в кДж/(кг-К) (для воды и водяного пара
R = 0,46151 кДж/(кг-К)).
2.1. Значения коэффициентов ah Ь; и С{
для формул (2.23), (2.24) и (2.25)
1	л.	bi	Ci
0	2125,05	3,23740^	5,608440^
1	1482,85	2,540^	-2,599340^
2	379,026	-1,1354-10’3	-1,2604-1О’8
3	46,174	-4,38140^	
4	10816,1	0,21	
45
Из формул (2.23) — (2.25) видно, что они легко реализуются толь-
ко при так называемой стандартной постановке задачи, когда задаются
параметры р и Z, а рассчитать необходимо остальные 3 параметра: v, h
и s. При всех других сочетаниях задаваемых параметров эти формулы
выступают уже в виде нелинейных уравнений, содержащих опреде-
ляемые параметры. В зависимости от поставленной задачи это приво-
дит к необходимости численным методом решить или одно из приве-
денных уравнений или даже систему из двух таких нелинейных урав-
нений. Подробный анализ этих ситуаций и предлагаемые алгоритмы
решения задач приведены в [10]. Мы же остановимся только на одном
из возможных вариантов нестандартной постановки задачи, решение
которой дает возможность познакомиться с еще одним оригинальным
методом решения систем нелинейных уравнений.
Давайте вспомним с вами изученную в курсе «Тепловые двигате-
ли и нагнетатели» методику расчета числа ступеней паровой турбины.
В качестве исходных данных здесь задаются параметры свежего пара
Pq и и давление в конденсаторе паросиловой установки рк. При рас-
чете числа ступеней задаются некоторым значением степени пониже-
ния давления в каждой из ступеней pi = р( = pz (обычно для активных
турбин р( = 0,65...0,75) и находят сначала приближенно число ступе-
ней Z (из формулы рк/pQ = pf ):
• _ 1п(а/^о)
П₽ In Pf
После этого величина Znp округляется до ближайшего целого и
уточняется принятое ранее значение р( (р, = (рк//?о)12чел)- При таком
подходе перепады давлений в каждой ступени получаются одинако-
выми, а перепады энтальпий (располагаемые теплоперепады AZzp) -
различными.
Далее для каждой из ступеней через величины коэффициентов
скорости, учитывающих потери энергии на преодоление внутреннего
трения при течении пара в соплах и в каналах рабочего колеса, находят
эти потери энергии, а также потери работоспособности в результате
частичного дросселирования потока при переходе от одной ступени к
другой. После этого строят изображение условного необратимого про-
цесса расширения, определяя при этом величины действительных теп-
лоперепадов А/?д, через которые и определяется действительная ско-
рость истечение пара из сопла и>д. Для наглядности на рис. 2.14 приве-
ден участок h - ^-диаграммы с изображением этих процессов.
46
Рис. 2.14. Н -s-диаграмма рабочего процесса многоступенчатой турбины
На лекциях или в учебной литературе по этой теме обращается
внимание на то, что такой подход не является наилучшим и использует-
ся только на первоначальных этапах проектирования. В дальнейшем
перепады давлений, а значит, и величины р( корректируются так, что для
каждой ступени степень уменьшения давления отличается от соседних.
Давайте представим себе, что имея математическую модель, по-
зволяющую определять параметры перегретого пара и рассчитывать
по ним важнейшие характеристики для любой ступени паровой турби-
ны, мы поставили себе задачу поискать более оптимальное распреде-
ление перепадов давления в каждой ступени, чем это получается при
изложенном выше подходе. Одной из рассматриваемых при этом ги-
потез может быть и такая: пусть перепады давлений в каждой из сту-
пеней выбираются из условия, что в любой ступени турбины срабаты-
вается одинаковый действительный теплоперепад при одинаковых
внутренних потерях энергии (рис. 2.15).
Рис. 2.15. Рабочий процесс в турбине
при одинаковых теплоперепадах &h в каждой ступени
47
В этом случае для любой ступени турбины мы можем найти значе-
ния параметров hi и по формулам пропорциональности для i = 1,2,Z:
Теперь возникает та самая сложная математическая задача, о ко-
торой говорилось выше: зная параметры h и s, необходимо найти соот-
ветствующие им значения t, р и у, используя приведенную выше мате-
матическую модель, описанную формулами (2.23) - (2.25). Задача бу-
дет решена, если нам удастся отыскать такие значения t и р, при кото-
рых разности между заданными значениями h и s и значениями Лр и £р,
рассчитанными по формулам (2.24) и (2.25), станут равными нулю или
эти параметры будут отличаться друг от друга на величину а.
Понятно, что попытки выразить из названных формул величины р
и т (чтобы потом, зная т, найти t) абсолютно неперспективны. Поэтому
задача может быть решена только численным методом. Для этого за-
дадимся некоторыми приближенными значениями искомых параметров
4тр и /?пр, величины которых можно рассчитать, например, приняв в пер-
вом приближении перегретый водяной пар за идеальный газ. В этом
случае h = cpm(tnp + 273), откуда tnp = h!cpm - 273 (срь перегретого водяного
пара приближенно оценивается величиной 1,54 кДж/(кг-К) [9]). Из фор-
мулы для расчета энтропии идеального газа
5 -С In—
р Т. р.
можно определить и второй приближенный параметр /?пр.
Теперь мы можем, зная значения /пр и рпр, по формулам (2.24) и
(2.25) вычислить приближенное значение величин h и s и сравнить их с
заданными значениями. Конечно же, такое сравнение покажет замет-
ную разницу, и наша задача - уменьшить эту разницу до приемлемой
величины. Чтобы определить стратегию такого уменьшения названной
разницы, зададимся вопросом: а как изменятся результаты расчета,
если величины /пр и рпр изменить на dt и dpi Поскольку и h, и s являют-
ся функциями двух аргументов (t и р), то соответствующие прираще-
ния, а это полные дифференциалы, будут определяться суммой част-
ных дифференциалов:
dh -
^OnpjPnp) ,
-------—dt
dt
dp
dp
(2.26)

48
^(^пр’/^пр) -	^(^npj/^np)
(2.27)
В результате новые значения расчетных параметров найдутся как
суммы предыдущих значений и дополнительных слагаемых dh и ds:
^Опрfi^j Т^прdp) h(tnp> Рпр)dh9
+ dt, p^ + dp) = s(fa> Pnp) + ds.
Найдем интересующие нас разности между заданными и расчет-
ными величинами h и s:
АЛ h [й(^,рпр) + dh], As s [s(^np,Т^пр) "I- ^].
Приравняв нулю эти разности, мы получаем систему уравнений
относительно dp и dt, решив которую мы найдем такие величины dp и
dt, при которых разницы Ай и А? равны нулю. Выполняем намеченное:
Л [й(4ф, /?пр) + dh] 0, s' [*?(^пр, /?пр) + fife] О*
Перепишем эти уравнения в виде системы в канонической форме,
подставляя вместо dh и ds их значения по формулам (2.26) и (2.27) и
отмечая искомые теперь величины специальным образом:
dKt Рпр)^+дК( Рпр)^ = h _ А(	(2 28)
dt	др
=5 (2 29)
Мы получили систему линейных уравнений вида
«и
dt + «i2fi^ = b},
a2\dt + «22^ - b
2 ?
решение которой легко находим обычным способом последовательно-
го исключения неизвестных (метод Гаусса):
-^dp,
«11
«11
49
Откуда после деления на находим
12“21
Если теперь подставить это значение в формулу для eft , то полу-
чим
_ д1Д~^21^2
Заметим, что коэффициенты а у - это значения частных производ-
ных, вычисленные при Znp и /?пр, а величины свободных членов bj - это
разности, стоящие в правых частях уравнений (2.28) и (2.29).
И хотя значения частных производных можно определить, диф-
ференцируя зависимости (2.24) и (2.25) с учетом всех их дополнитель-
ных усложнений, проще эти производные определить численным спо-
собом, задаваясь некоторыми малыми приращениями Л4 и A/?i (на-
пример, Afi = 5 °C и Api = 0,01 МПа) и проводя расчеты по формулам:
_ ^Опр’-Атр) _ ^(^пр + ^1’—	/?пр)
^(/пр ? /?пр )
ф
^(/пр’.Рпр+Да) ^^пр’Рпр)
Api
^Опр > Рпр У *^0lip "I"	5 Рпр ) ^(^пр J Рпр )
6Z71 —	—	j
dt
„ _ ^Опрэ-Рпр) _ ^(/пр’ Рпр +	— ^(Лтр’/^пр)
др	Ар}
Далее, когда мы рассчитали значения dt и dp, находим после-
дующее приближение искомых параметров:
^npl — ^пр +	Рпр1 — -Рпр + df) •
Теперь можно при новых приближенных значениях 4^1 и /?пр1 по
формулам (2.24) и (2.25) снова рассчитать значения h и л1. Поскольку
наш теоретический анализ велся для бесконечно малых изменений, а
численные расчеты были проведены для небольших, но конечных
50
приращении, конечно же, при новых значениях Znpi и /?пр] разницы ме-
жду заданными и расчетными значениями h и s не будут в точности
нулевыми, однако несовпадение будет намного меньше, чем было при
первом приближении.
Принимая результаты полученного приближения за исходную
точку и повторяя описанный выше алгоритм расчетов, можно в таком
итеративном цикле уточнять величины tnpj и рпр/- до тех пор, пока раз-
ницы не приблизятся к нулю (станут меньше а). Чтобы уменьшить
число итераций, а значит и суммарное время счета, обычно использу-
ется тактика «экономной стрельбы». Суть здесь в том, что после рас-
чета величин dt и dp организуется специальный дополнительный
цикл определения значений 4^1 и рпр\ с шагом, равным некоторой доли
(например, по 0,36) от найденных dt и dp .
Чтобы лучше понять ход наших рассуждений, разобьем проведен-
ный нами анализ на 2 этапа и проиллюстрируем эти этапы графически.
На рисунке 2.16 на h-s-диаграмме показан первый этап наших рассуж-
дений. Точка 1 соответствует заданным значениям величин h и s, при
которых следует найти значения параметров р и t. Точка 1пр определена
по приближенным значениям /?пр и Znp. Далее по формулам (2.26) и (2.27)
найдены значения dh и ds, что позволило определить положение точки
1™в и найти величины Лй и As, а также вектор Ар Чтобы точки 1 и 1™в
совпали, необходимо величину 7?i устремить к нулю.
Это может быть достигнуто только тогда, когда Ай = 0 и As = 0.
Эти два полученных равенства с учетом формул (2.28) и (2.29) приво-
дят нас к системе простых линейных алгебраических уравнений, ре-
шением которой находим величины dt и <$?. Это позволяет опреде-
лить такие новые значения /пр] и /?прь при которых теоретически равен-
ства Ай = 0 и As = 0 соблюдаются и вектор R равен нулю. Именно этот
этап расчетов и отражен на рис. 2.17.
Рис. 2.16. Первый этап решения задачи
51
Рис. 2.17. Второй этап решения задачи
В действительности же вектор R оказывается не нулевым, поэто-
му шагая вдоль него малыми шагами, мы можем приблизиться (или
даже немного перешагнуть) к точке 1. Как только мы перешагнули
точку 1, прежние шаги заменяются на вдвое меньшие и меняются зна-
ки на противоположные. В результате такого алгоритма можно при-
ближаться к точке 1 до тех пор, пока АТ? не станет < s.
2.7.	ВОПРОСЫ И ЗАДАНИЯ ДЛЯ САМОПРОВЕРКИ
1.	Перечислите методы решения нелинейных уравнений вида
F(x) = 0 или х = ф(х).
2.	В чем суть метода сканирования при решении уравнения вида
F(x) = О?
3.	В чем суть метода половинного деления при решении уравне-
ния вида F(x) = О?
4.	В чем суть метода хорд при решении уравнения вида F(x) = О?
5.	Расскажите о трех формах записи системы линейных алгеб-
раических уравнений.
6.	Какие системы линейных уравнений называют несовмести-
мыми?
7.	По какому признаку совместимые системы алгебраических
уравнений подразделяют на определенные и неопределенные?
8.	Перечислите три основные свойства совместимых систем ли-
нейных уравнений.
9.	В чем суть решения системы линейных уравнений методом
Г аусса?
10.	Как конкретно проводится последовательное исключение
неизвестных в методе Гаусса? Какие свойства системы при этом ис-
пользуются?
52
11.	Перечислите те «пять шагов к победе», которые следует вы-
полнить при матричном решении системы линейных уравнений.
12.	Почему возникает необходимость и как осуществляется ре-
шение системы линейных уравнений итерационным методом? Можно
ли в этом случае получить точное решение системы?
13.	Сформулируйте условие сходимости итерационного процес-
са при решении системы уравнений итерационным методом.
14.	Какой особенностью отличается метод Гаусса — Зейделя от
обычного итерационного метода решения системы уравнений?
15.	Каковы особенности применения метода простых итераций
для решения систем нелинейных уравнений при квазилинейной поста-
новке задачи?
16.	В чем суть метода Ньютона, используемого для итерацион-
ного решения системы линейных или нелинейных алгебраических
уравнений?
17.	Каковы особенности применения метода простых итераций
для решения систем нелинейных уравнений?
18.	При каком условии итерационное решение системы нели-
нейных уравнений сходится?
19.	Чем отличаются алгоритмы параллельного и последователь-
ного итерирования при решении систем нелинейных уравнений?
20.	В чем особенность метода «пристрелки» при решении сис-
темы нелинейных уравнений?
21.	Как ставится задача определения неизвестных в методе
«пристрелки»?
22.	Как организуется минимизация невязки между величинами
известными и найденными приблизительно этих параметров?
23.	В чем особенность «экономной стрельбы» при решении сис-
темы нелинейных уравнений?
53
3.	МАТЕМАТИЧЕСКИЕ МОДЕЛИ НЕУСТАНОВИВШИХСЯ
ПРОЦЕССОВ И МЕТОДЫ ИХ РЕАЛИЗАЦИИ
3.1.	ЗАДАЧА ОБ ИСТЕЧЕНИИ ГАЗА ИЗ ОГРАНИЧЕННОГО ОБЪЕМА
Напомним предварительно, что установившимся режимом приня-
то называть такое функционирование исследуемого объекта, когда его
параметры не изменяются с течением времени. Если же с течением
времени параметры объекта изменяют свое значение, то говорят, что
реализуется неустановившийся, нестационарный процесс.
При математическом моделировании инженерных задач, включая
и задачи, возникающие в теплоэнергетике, очень часто встречаются
неустановившиеся процессы. Как правило, их удается описать некото-
рыми дифференциальными уравнениями (или системой дифференци-
альных уравнений). Для объектов с сосредоточенными параметрами
обычно получаются обыкновенные дифференциальные уравнения, а
для систем с распределенными параметрами - дифференциальные
уравнения в частных производных.
Естественно, что знакомство с приемами разработки математиче-
ских моделей для неустановившихся процессов начнем с более про-
стых задач, т.е. с объектов с сосредоточенными параметрами.
В этом случае обычно задача формулируется так: следует опреде-
лить, как меняются отдельные параметры объекта в течение некоторо-
го промежутка времени, если известны состояние его в начальный мо-
мент времени и величина этого промежутка. Математики такую поста-
новку задачи называют задачей Коши. Основные подходы к разработ-
ке математической модели и последующему решению задачи Коши
продемонстрируем на ряде конкретных примеров.
Сначала рассмотрим такую задачу: рассчитать продолжитель-
ность истечения воздуха из баллона объемом V = 3 м3 и давление в
баллоне через 60 с после начала процесса, если это истечение происхо-
дит в атмосферу (В = 0,1 МПа) через отверстие диаметром d = 0,01 м.
Начальное давление в баллоне />нач = 0,18 МПа, температура ^нач =
= 18 °C. Изменением температуры воздуха в баллоне в процессе исте-
чения пренебрегать, считать также, что давление в баллоне в течение
процесса распределяется равномерно.
Найдем отношение давлений р = В/р^ = 0,1/0,18 = 0,556. По-
скольку р > ркр (для воздуха ркр = 0,528), то течение газа будет дозву-
ковым. По мере истечения воздуха давление в баллоне будет умень-
шаться до тех пор, пока не достигнет барометрического, после чего
процесс прекратится (см. рис. 3.1).
54
Рис. 3.1. Изменение давления в баллоне
В некоторый момент времени т давление в баллоне будет р, а
мгновенный массовый расход воздуха будет М(р). За элементарно ма-
лый промежуток времени dx при этом из баллона будет вытекать dm кг
воздуха:
-dm = M(p)dt.
Приведенную формулу можно рассматривать как дифференци-
альное уравнение, описывающее зависимость т =fiy). Перепишем его
пока так:
dm
dx
(3.1)

Массовый расход воздуха М(р) через отверстие определяется ве-
личиной поперечного сечения отверстия и перепадом давлений р-В.
Сжатие струи и другие препятствующие эффекты принято учитывать
величиной коэффициента истечения ц (примем ц = О, 9). Из термоди-
намики потока известно, что величину расхода легко рассчитать по
формуле [3]:
, ч то/2
(3.2)
где v = V/m = удельный объем воздуха; к — показатель адиабаты.
Величину т можно связать с давлением в баллоне и по-другому.
Для этого запишем уравнение состояния для воздуха в баллоне, спра-
ведливое для любого момента времени т
pV - mRT
и продифференцируем его по времени, учитывая, что во время истече-
ния воздух в баллоне расширяется очень постепенно и это не приводит
55
к неоднородности давлений по объему баллона, а также что темпера-
тура воздуха в баллоне практически не уменьшается в результате теп-
лопритока через стенки баллона:
= RT—.
dt dt
Выразим отсюда производную
dm _ V dp
dt RT dt
и подставим значение производной в формулу (3.1):
V dp
RT dt

Разнеся переменные, мы получили дифференциальное уравнение,
описывающее искомую зависимость т =fip\.
dt --------
RTM(p)
(3.3)
Интегрируя это уравнение, получим
^кон
О Рнач
RTM(p)
К сожалению, из-за сложности функции М(р) проинтегрировать
правую часть в квадратурах не удается. Более того, если бы смогли это
сделать, то за бортом нашего анализа осталась бы зависимость р =flt),
поскольку для нее мы имели бы только две точки: давление в начале
процесса />нач при т = 0 и давление В при т = ткон. А как это следует из
поставленной задачи, на практике важно знать зависимость р = Дт),
чтобы определять давление р в баллоне в различные моменты времени
т. Поэтому формулу (3.3) перепишем по-другому:
RTM (р)
dn —---------dt.	(3.4)
Решить такое уравнение аналитически тоже невозможно, поэтому
решать будем численным методом, заменяя бесконечно малые прира-
щения dp и dt малыми конечными приращениями:
RTM(p)
56
При такой замене мы мысленно как бы заменяем реальный не-
прерывный процесс некоторым условным скачкообразным процессом,
в котором в течение интервала Ат давление в баллоне не изменяется, а
уменьшение его происходит скачком в начале следующего интервала,
как показано это на рис. 3.1. Для такого скачкообразного процесса по-
лучается, что в течение интервала Ат по формуле (3.2) легко рассчи-
тать величину расхода М(р) и тогда для расчета искомой функции
можно использовать следующую итерационную формулу
ИТ
к+1р=кр +\кр=кр - у-М(кр)^
(3.6)
где символом кр обозначено предыдущее приближенное значение
давления /?, зависящее от величины принятого шага интегрирования
Лт, а к+1р - последующее его значение.
В результате нашего анализа мы пришли к так называемой итера-
ционной схеме Эйлера. Это простейшая явная схема, поскольку значе-
ние к+1р вычисляется по явной формуле, не содержащей в правой
части неизвестной к+1р. Можно построить и другой алгебраический
аналог дифференциального уравнения (3.4), при этом величина к+1р
будет входить и в правую часть формулы (3.6). Тогда для определе-
ния к+1р придется решить нелинейное уравнение, например, приводя
его к виду F(x) = 0. В этом случае говорят о неявной схеме расчетов.
Подробнее с этим подходом мы познакомимся позже.
Чем меньше будет величина Ат, тем ближе полученное давление
кр будет совпадать с действительным давлением />, которое получа-
лось бы при аналитическом решении уравнения (3.4). Чтобы различать
эти результаты, вводится понятие о погрешности численного решения,
которая определится разностью
&i = p- р-
Если при уменьшении величины Ат эта погрешность уменьшается
(а это во многом зависит от вида искомой функции), то говорят, что
применяемая разностная схема обладает свойством сходимости. Точ-
нее признаком сходимости является
lim^o S/
= 0.
Чтобы разностная схема обладала таким свойством необходимо и
достаточно выполнение двух других условий: аппроксимируемости и
устойчивости. Об этих понятиях поговорим в дальнейшем.
57
Когда за время dx вытекает dm кг газа, то оставшаяся часть
(т - dm) расширяется, заполняя весь объем V. При получении диффе-
ренциального уравнения (3.4), алгоритм численного интегрирования
которого приведен выше, мы рассмотрели один из предельных случа-
ев, характеризующейся высокой интенсивностью теплообмена между
газом и внешней средой. Настолько высокой, что уменьшение темпе-
ратуры в результате расширения газа полностью компенсируется ее
увеличением за счет притока теплоты извне, в результате чего процесс
идет при Т = const.
Однако можно рассмотреть и другой предельный случай, когда
теплообмен между газом и окружающей средой практически исключен
(баллон теплоизолирован или изготовлен из материала с очень малой
теплопроводностью). Тогда в процессе вытекания газа оставшаяся в
баллоне часть газа будет расширяться адиабатически и температура
его будет уменьшаться по мере уменьшения давления р. В этом случае
в итерационной формуле (3.6) величину Т тоже нужно будет опреде-
лять в зависимости от кр :
мр =кр +Ькр =кр - RT^ м^р)^.
Чтобы найти зависимость Т(кр), вспомним основы термодина-
мики и запишем известную связь между параметрами для адиабатного
процесса [3]:
х(к-1)/к
^/?нач )
где Тър- температура и давление газа в любой точке процесса; Тнач и
Лич - эти же параметры в начальной точке процесса, а через к обозна-
чен показатель адиабаты (для двухатомных газов к = 1,4). Из приве-
денной формулы находим требуемую зависимость
х(к-1)/к
которую следует включить в итерационную формулу для расчета ве-
личины T(kp).
Поскольку действительные давления р нам неизвестны, то и ве-
личину 8j на каждом шаге интегрирования мы не знаем. Поэтому при
реализации формулы (3.6) (или ее усложненного аналога) величину
шага Лт принято определять автоматически, реализуя следующий ите-
58
рационный алгоритм: сначала задают достаточно большой шаг Ат и
рассчитывают величину д, вычисляя правую часть формулы (3.6) по
Ро - Агач • Затем величина шага делится пополам (Ati = Ат/2) и снова
рассчитывают давление рх в конце шага Ат, выполняя 2 уменьшенных
шага Ati (см. рис. 3.2).
Рис. 3.2. Автоматический подбор шага Дт;
Г Вход Л
I с вводом I
I ЛТнач, pjt |
\R,T,VhmJ
____ I
Расчёт Дт:=ДтНЕН,
Аон по формуле (3.6)
J -
Расчёт Дт:=Дт/2 *
О
да
Рис. 3.3. Блок-схема процедуры автоматического выбора
шага интегрирования
59
Рис. 3.4. Укрупненная блок-схема алгоритма численного интегрирования
обыкновенного дифференциального уравнения методом Эйлера
Если полученный новый результат с заданной точностью (напри-
мер, с точностью до 0,001 МПа) совпадает с предыдущим, то с этим
новым шагом рассчитывается следующая итерация. Если же требуемо-
го совпадения результатов не было получено, то шаг Atj еще раз де-
лится пополам и опять проводится сопоставление последующего и
предыдущего результатов расчета рг для начального шага Ат. Такое
деление шага будем продолжать до тех пор, пока заданное совпадение
не произойдет, и только после этого сделаем очередной шаг по време-
ни и начнем рассчитывать р2 с шагом Атвыбр, полученным после всех
его предыдущих делений. На рисунке 3.2 приведено графическое ото-
бражение алгоритма автоматического подбора величины шага Атвыбр, а
на рис. 3.3 приведена блок-схема этого алгоритма.
Заканчивая рассмотрение задачи, отобразим изложенный алго-
ритм решения нашей задачи в виде блок-схемы на рис. 3.4, содержа-
щей укрупненный блок (процедуру) автоматического определения оп-
тимального шага Атвыб.
3.2.	ИССЛЕДОВАНИЕ РЕЖИМОВ РАБОТЫ СМЕСИТЕЛЬНОГО
ТЕПЛООБМЕННИКА
Для закрепления методики разработки математических моделей
для неустановившихся процессов рассмотрим еще одну простую зада-
чу. В смесительном пароводяном теплообменнике барботажного типа
находится 2 м3 воды при температуре /нач = 20 °C. В начале рабочего
дня для приготовления горячей воды на технологические нужды от-
60
крывают вентиль и подают в теплообменник сухой насыщенный водя-
ной пар с температурой насыщения ?нас = 120 °C. При полном открытии
вентиля расход пара составляет Мп = 0,2 кг/с. Считая, что высота стол-
ба воды в смесителе настолько велика, а отверстия в барботажных
трубках настолько малы, что поднимающиеся вверх пузырьки пара
успевают полностью сконденсироваться, определите, через сколько
времени вода в теплообменнике нагреется до температуры /кон = 85 °C.
Расчет провести для двух вариантов:
-	без учета теплопотерь от наружной поверхности теплообмен-
ника в окружающую среду;
-	с учетом теплопотерь при условии, что коэффициент теплопе-
редачи от воды к окружающему воздуху к = 15 Вт/(м2-К). При этом
температура воздуха 4ОЗ =12 °C, а общая поверхность теплообмена
5*то = 5 м .
Основу математической модели рассматриваемой задачи состав-
ляют соответствующие тепловые балансы. Разберемся с логикой и ме-
тодикой их составления.
Для неустановившихся режимов всегда рассуждают так: пусть в
некоторый момент времени т температура воды в аппарате t. Тогда за
малый промежуток времени dt источники теплоты отдадут dQ\ тепло-
ты, а приемники теплоты получат dQ2 теплоты. При этом температура
воды за время dt изменится на величину dt. Равенство dQ\ = dQ^ и вы-
ражает собой тепловой баланс за этот промежуток времени t/u.
В нашем случае источником теплоты является конденсирующий-
ся водяной пар. Мы знаем, что при конденсации каждый килограмм
пара выделяет теплоту парообразования г и после полной конденсации
образуется 1 кг воды с температурой /нас. Часть этой теплоты тоже
пойдет на нагревание воды, находящейся в аппарате. Кроме того, сле-
дует учитывать, что в течение процесса разогрева масса воды в тепло-
обменнике увеличивается за счет образующегося конденсата, причем
за время t/т увеличение массы будет dm = Mndx. Итак, количество теп-
лоты, отданное на нагрев воды за время dt
dQi M^dtr “И Л/пб/тСВод(^нас /),
где первое слагаемое описывает количество теплоты, отданное паром,
а второе - отданное конденсатом на нагрев воды при температуре ее t
в момент времени т.
Теперь рассчитываем вторую часть баланса (пока без учета теп-
лопотерь с поверхности теплообменника):
t/^2 (^вод ~+" ^/г/О^воД^э
где произведение Мпт описывает количество конденсата, образовавше-
гося из пара к моменту времени т, а свод - удельная теплоемкость воды.
61
Приравнивая правые части приведенных формул и разделив ре-
зультат на dx, получаем
водЧ*нас
- 0] = (^вод + MnT)cWJSdt/dx
или в канонической форме
dt _ ^П[г + Свод(^нас 0]
dx (^вод "I" ^пт)свод
(3.7)
Обозначим правую часть формулы (3.7) через F(t,x) и запишем
наше дифференциальное уравнение в общем виде
Начальные условия задачи нам известны: при т = 0 /нач = 20 °C,
т.е. мы имеем обычную задачу Коши, решаем которую самым простым
способом - по явной схеме метода Эйлера:
4+1 / =^ + Р(^,т)Ат.
Итерации будем повторять до тех пор, пока i+1Z > ZK0H (ZK0H =
= 85 °C), при этом на каждом шаге итераций фиксируется текущее
время *+1 т =\ +Ат.
Блок-схема программы для ПК во многом аналогична предыду-
щей (см. рис. 3.5).
Если учитывать и теплопотери с поверхности аппарата, то вели-
чину dQz нужно будет дополнить еще одним слагаемым. Это второе
слагаемое находится в соответствии с основным уравнением теплопе-
редачи и описывает величину внешних теплопотерь за интервал вре-
мени dx:
dQ^ (^вод -Л^пОСвод^ ~+" t^S^dx.'
При этом соответствующим образом изменится формула для рас-
чета F(/,t), числитель которой будет уменьшен на величину А(/ - 4оз)^то-
Алгоритм решения останется прежним, но поскольку величина F
уменьшилась, термограмма разогрева будет проходить заметно положе
и время на разогрев до /кон соответственно увеличится.
На приведенных примерах мы познакомились с простейшим ме-
тодом интегрирования обыкновенных дифференциальных уравнений.
Рассмотрим теперь этот вопрос в общем случае.
62
Рис. 3.5. Блок-схема алгоритма
Итак, необходимо решить задачу Коши, сформулированную, так
же, как в рассмотренных примерах:
— = /(w),
max
Ут = о=Уо.
Записав для (£ + 1)-го момента времени соотношение
к+ху-ку + Ду, разложим искомое решение у	в ряд Тейлора:
/(т + б/т) = /(т) +
Если отбросить слагаемые более высоких порядков малости и пе-
рейти к малым конечным приращениям Лт, то получим
Дт+Ат) = /(т) + -у- Ат
ат
или, заменяя fir) на соответствующие значения у и учитывая, что
dx dx
- F(t, у) , можем записать
63
ми=ки+/(тЛг/)Дт.
(3.8)
Здесь вместо переменной у введено ее приближенное значение U,
при этом понятно, что разница между истинным значением у и вели-
чиной U существенно зависит от величины шага Ат, о чем упомина-
лось выше.
Формула (3.8) описывает реализацию метода Эйлера по так назы-
ваемой явной схеме, когда по значению функции на предыдущем шаге
находится ее значение на последующем по времени шаге.
Более сложный, но гораздо более точный подход реализуется при
так называемой неявной схеме, когда значение функции/рассчитыва-
ется через последующее приближение k+1U величины £7, т.е. по фор-
муле
ми=ки + /(т,4+1 £/)Ат.	(3.9)
При этом на каждом шаге по т приходится решать, как правило,
нелинейное уравнение относительно k+1U. Доказано, что такая неяв-
ная схема всегда устойчива, вне зависимости от величины шага Ат, в
то время как явная схема устойчива только при очень малых Ат, удов-
летворяющих специальному условию (например, для задач нагрева -
охлаждения это темп охлаждения пт. Ат < 2/т).
Схемы Эйлера содержали только один из членов ряда Тейлора,
поэтому их называют еще схемой первого порядка аппроксимации.
Схемы более высоких порядков будут обеспечивать более высокую
точность интегрирования, поскольку будут включать большее число
слагаемых. Правда, при этом приходится вычислять вторые, третьи
и т.д. производные от функции/(т, у), что усложняет расчеты. Поэтому
возникает потребность в разработке методов, обеспечивающих высо-
кий порядок точности, но не требующих расчета производных второго,
третьего и других порядков.
3.3.	МЕТОДЫ СИМПСОНА, РУНГЕ-КУТТА, АДАМСА
Запишем выражение для приращения функции у =Дт) на некото-
ром интервале времени от до ад
Л+1 - У, =	= I.
Напомним, что величина интеграла I определяется площадью под
кривой, изображающей анализируемую функцию. Эту площадь числен-
ным способом приближенно можно определить с разной степенью точ-
ности. Простейшим (но наименее точным) из таких подходов является
64
метод прямоугольников (см. рис. 3.6, схемы а и Ь). Для его реализации
достаточно вычислить значение функции fir, у) только один раз в начале
или конце выделенного интервала (при или t,+i). При этом заштрихо-
ванные горизонтально площадки характеризуют величину погрешности
такого метода. Несколько повысить точность определения величины
интеграла I можно, если рассчитать значение функции в середине ин-
тервала и это значение принять за высоту прямоугольника (рис. 3.7).
Более точным является метод трапеций, представленный на рис. 3.8.
Здесь вычисление значения функции проводится дважды, на краях
интервала [т„ т,+1], а площадь под кривой находится как площадь за-
штрихованной трапеции. Как это видно из рисунка, площадь трапеции
гораздо более точно отражает величину площади под кривой.
Еще более точным является метод Симпсона. В простейшем случае
расчетный интервал [т(, т(+1] делится пополам и вычисляют 3 значения
функции - на краях интервала и в его середине, т.е. при тс = т, + (тГ| -
- т/)/2. Далее через полученные точки проводится квадратичная парабола
/1(т) = Лт2 + Вт + С.
Рис. 3.6. Определение величины интеграла
на отрезке г/+)] методом
прямоугольника по крайним точкам
Рис. 3.7. Определение
величины интеграла на
отрезке	т ;+1] методом
прямоугольника
по средней точке
Рис. 3.8. Графическое
представление метода трапеций
Рис. 3.9. Графическое
представление метода Симпсона
65
Записав эту формулу для каждого из трех моментов времени (т<, тс
и т,+1), получаем замкнутую систему из трех уравнений, позволяющую
рассчитать значения параметров Л, В и С на рассматриваемом интер-
вале. Если теперь проинтегрировать новую функцию fi, то мы получим
формулу для расчета приближенного значения величины I. Результаты
такого решения и соответствующая расчетная формула приведены в
учебной и справочной литературе. Графическая интерпретация этого
метода показана на рис. 3.9. Конечно же, в качестве интерполяционно-
го полинома можно использовать полином более высокой степени и
получать более точные результаты. Однако это связано с увеличением
вычислений, поэтому, как правило, используется только полином вто-
рой степени.
Понятно, что приведенные формулы не приемлемы для интегри-
рования дифференциальных уравнений при решении задачи Коши,
поскольку точное значение искомой функции мы знаем только при
т = 0. Схемы Эйлера первого порядка можно представлять получен-
ными путем замены величины интеграла I простейшими формулами
tty или 7«At-F(t/+i, Ц+i),
где через Ц обозначено приближенное значение искомой функции при
известной величине ее производной F.
В схемах Рунге - Кутта используются приемы уточнения резуль-
татов интегрирования, рассмотренные выше. Здесь тоже исследуемый
интервал делится на две части и значение производной F рассчитыва-
ется 2 раза: для тж и для Тр После этого сначала по формуле Эйлера
находят промежуточное приближенное значение функции
Гпром "I- F(y, Т/)Ат.
После этого уточняют значение интеграла I по формуле трапеций
I & [0,5F(C4 т() + ОЛ^С/пром, т/+1)]Ат
и определяют более точное значение функции
ui+1= Ui+I* Ui+ 0,5M[F(Ui, тг-) + F(CZnp0M, тж)].
Рассмотренный алгоритм называют схемой второго порядка, гра-
фическая интерпретация этого подхода приведена на рис. 3.10.
В настоящее время самое широкое применение находит метод
Рунге — Кутта четвертого порядка. В таком методе интервал [т(, тж]
сначала делится на 4 равных части. Далее рассчитывают значения про-
изводных на границах каждого участка и находят по формулам Эйлера
приближенные значения функции U в этих точках. По результатам
66
Рис. 3.10. Определение величины Ui+1 на интервале -с;+1]
методом Рунге - Кутта второго порядка
этих расчетов для каждого интервала Лт/2 находят приближенное зна-
чение интеграла по квадратурной формуле Симпсона и, суммируя ре-
зультаты, рассчитывают величину Ц+ь В итоге получена следующая
итерационная формула
Ui+1 = Ui + Ax/6-(Fi + 2F2 + 2F3 + F4),
где функции Fi, F2, F3 и F4 вычисляются по следующим формулам:
F, = F(yh Uj),
F2 = F(t,- + bx/2,Uf + At/2Fj),
F3 = F(t; + At/2, Uj + At/2F2),
F4 = F(ri+1, Ui + AxF3).
Приведенный алгоритм запрограммирован в разных названных
выше математических пакетах, опубликован в соответствующей лите-
ратуре (например, на Паскале), его можно просто позаимствовать из
Интернета.
Известны также и другие подходы, в частности методы Адамса, в
которых для расчета величины Ui+i используют результаты расчетов
на нескольких предыдущих интервалах времени, т.е. значения Ц, Ц_1,
Ui-г, Ui.3 и Uj.4 (в методах 4-го порядка).
В заключение несколько слов о решении систем обыкновенных
дифференциальных уравнений
^1	17/	А
—Г=^(У,У1,У2)>
dx
= Р2(иУ1’У2^
67
maxj JVlT-Q-JOb У2т = 0“№«
Разностные схемы для системы уравнений записываются так же,
как и для одного уравнения. Например, по явной схеме Эйлера это бу-
дет система уравнений:
k+lU. =kU, +\iFl{xikU1 ,kU2),
MU2=kU2+\TF2(Ti,kUxkU1),
где верхними индексами отмечены предыдущие и последующие при-
ближенные значения искомых функций, рассчитываемые в итерацион-
ном процессе: т/ - начало очередного интервала времени Ат. Решения
каждого уравнения обычно следуют друг за другом.
Некоторые сложности могут возникать при интегрировании так
называемых жестких систем, но об этом нужно читать специальную
литературу.
Обыкновенные дифференциальные уравнения более высоких по-
рядков предварительно введением новых переменных сводятся к сис-
теме обыкновенных дифференциальных уравнений первого порядка,
каждое из которых решают методом Эйлера (или Рунге - Кутта), на-
чиная с самого нижнего.
Пусть, например, нужно проинтегрировать дифференциальное
уравнение
В качестве новой переменной примем значение первой производ-
dy _
нои по у: ср =	. Тогда уравнение второго порядка через эту пере-
dx
менную можно записать как уравнение первого порядка:
dx
= F{y, т).
Решив его, мы найдем значения производной ф, в каждый момент
времени т(. Поле этого переходим к решению таким же способом диф-
ференциального уравнения
dx
= ф
и находим значения у на всем интервале интегрирования от 0 до ткон.
68
3.4.	ВОПРОСЫ И ЗАДАНИЯ ДЛЯ САМОПРОВЕРКИ
1.	Расскажите об особенностях формулирования и математиче-
ского описания нестационарных задач в области теплоэнергетики.
2.	Какие упрощающие предпосылки вводятся при моделирова-
нии процесса истечения воздуха из баллона через ограниченное отвер-
стие?
3.	В чем суть итерационной явной схемы Эйлера, используемой
для численного интегрирования простых дифференциальных уравне-
ний?
4.	Что называют сходимостью разностных схем численного ин-
тегрирования дифференциальных уравнений?
5.	Опишите алгоритм автоматического определения величины
шага интегрирования, обеспечивающего сходимость решения.
6.	В чем суть итерационной неявной схемы Эйлера, используе-
мой для численного интегрирования простых дифференциальных
уравнений?
7.	Как при численном интегрировании дифференциальных урав-
нений реализуется метод Адамса?
8.	Что положено в основу методов Рунге — Кутта, используемых
при численном интегрировании дифференциальных уравнений?
9.	Расскажите об особенностях численного интегрирования диф-
ференциальных уравнений второго и более высоких порядков.
69
4.	ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ ТЕПЛООБМЕНА
4.1.	ЧИСЛЕННОЕ РЕШЕНИЕ ЗАДАЧ СТАЦИОНАРНОЙ
ТЕПЛОПРОВОДНОСТИ
Вопросам математического описания процессов теплообмена по-
священа отдельная область научных исследований, получившая назва-
ние «Теория тепло- и массообмена» или просто «Теплопередача».
В этой науке подробно рассматриваются основные законы природы,
определяющие интенсивность теплообмена, разрабатываются инже-
нерные методики расчета типовых процессов, выводятся дифференци-
альные уравнения, описывающие различные классы явлений, связан-
ных с тепло- или массообменом. Поскольку углубленное изучение од-
ноименной дисциплины предусмотрено государственным стандартом
образования на предыдущем этапе вашего обучения в стенах универ-
ситета, в дальнейшем мы будем широко ссылаться на полученные ра-
нее знания.
Аналитическое решение задач теплопроводности возможно лишь
для тел простой геометрической формы и при простейших граничных
условиях. На практике же часто возникает необходимость определить
температурное поле в телах более сложной формы или при таких ус-
ловиях однозначности, когда температура или условия теплообмена на
поверхности тела изменяются в пространстве или по времени, или ко-
гда величина X существенно и нелинейно зависит от t, или когда тело
неоднородно и величина X различна в разных точках тела и по разным
направлениям. Во всех этих случаях задача может быть решена только
численным методом.
Чтобы перейти к численным методам, исследуемое тело мыслен-
но разделяют на небольшие объемы простой формы (чаще всего пря-
моугольной). При этом считают, что в пределах каждого такого объема
свойства вещества, мощность внутренних источников и температура
остаются постоянными, а изменение этих параметров происходит
скачками на границах каждого объема. Другими словами, непрерыв-
ный процесс теплопроводности заменяется некоторым дискретным в
пространстве и по времени процессом.
Центральные точки выделенных объемов (их называют узлами)
образуют внутри тела пространственную сетку. Для любого узла такой
сетки на основе теплобалансовых уравнений или путем замены диф-
ференциального уравнения теплопроводности его конечно-разностным
аналогом (от бесконечно малых приращений переходят к малым ко-
нечным приращениям) можно получить алгебраические соотношения,
в совокупности составляющие замкнутую систему уравнений, для ре-
шения которой используются стандартные или специально разрабо-
70
тайные методы. Изложенный подход называют методом конечных
разностей или методом сеток.
В качестве примера рассмотрим двумерное температурное поле,
возникающее в однородной теплопроводной пластине толщиной 6
(толщина пластины намного меньше остальных ее размеров), когда
температуры на боковых поверхностях ее различны (см. рис. 4.1, по-
верхности А, В, С и D). Разделим пластину на элементарные прямо-
угольники со сторонами Дх и Ду, выделим центры масс прямоугольни-
ков и будем считать, что в пределах выделенного элемента температу-
ра его постоянна, а изменяется она скачками на границах каждого из
элементов.
Изобразим теперь на рис. 4.2 расчетную сетку с шагом Дх по оси х
и Ду по оси у. Промаркируем узлы сетки вокруг одного из элементов,
лежащего внутри пластины, обозначая их номера вдоль оси х индексом
/, а вдоль оси у - индексом J.
Температурное поле такой пластины будет плоским, и дифференци-
альное уравнение для него при стационарных условиях принимает вид:
(4.1)
тт	д2*
Производную —у
дх
представим в виде
и заменим беско-
нечно малые приращения dtwdx малыми конечными величинами Д/ и Дх:
Рис. 4.1. Пластина, разбитая
на прямоугольные элементы
Рис. 4.2. Расчетная сетка
71
Здесь отношение Az/Ax - величина, близкая к величине проекции
температурного градиента на ось х, а выражение А( Az/Ax) представляет
собою разницу между величинами (Az/Ax) справа и слева от анализи-
руемого узла. Выпишем значения этих отношений:
Значит
d2Z ~	+
дх2	Ах2
и далее
Аналогичные рассуждения для оси у позволяют получить
d2t
В итоге дискретный аналог дифференциального уравнения (4.1)
представляется в виде
(4.2)
Это уравнение описывает связь между температурами соседних
элементов для любого из внутренних узлов сетки. Для узлов, принад-
лежащих внешним границам пластины (всего их пгр^ температура или
известна (если заданы ГУ-1), или может быть определена из уравне-
ний, описывающих граничные условия. Таким образом, для всех п(п =
= ij - и^) неизвестных температур формула (4.2) дает п линейных ал-
гебраических уравнений. Эти уравнения заметно упрощаются, если
принять одинаковыми величины шагов по х и у: Ах = Ау:
^+1, j + *i-i, j + j+i + *i, j-i "Ч, j о.
(4.3)
Прежде чем говорить о методах решения полученной системы, от-
метим три важнейших свойства разностных схем: аппроксимируемость,
устойчивость и сходимость решения. Первое означает, что при Ах О
и Ау 0, т.е. решение системы алгебраических уравнений стремится к
72
решению исходного дифференциального уравнения. Устойчивой на-
зывается схема, для которой ошибки округления при уменьшении ша-
гов Ах и Ду не приводят к большим искажениям решения. Сходимость
означает, что по мере уменьшения Дх и Ду решение системы становит-
ся все точнее, все лучше приближается к истинному решению. Сходи-
мость выступает как следствие аппроксимируемости и устойчивости.
Анализ различных конечно-разностных схем на устойчивость и схо-
димость приведен в специальной литературе.
Полученную методом сеток систему уравнений решают обычно
на ЭВМ методом последовательного исключения неизвестных (метод
Гаусса) или другими методами, о которых речь пойдет позже. Теперь
же рассмотрим другой оригинальный метод - метод релаксаций. Суть
этого метода в следующем. Сначала в узлах сетки записывают ожи-
даемые, интуитивно выбранные значения температур. Конечно же, они
не будут удовлетворять уравнению (4.3) и вместо равенства нулю в
каждом узле сетки мы будем получать некоторый остаток Rig:
RiJ ~ ^+l,j + h-l,j + h, j+l + h, j-l ~ ^i,j •
Величина говорит о том, насколько удачно и правильно были
выбраны значения температур в окрестностях каждого узла в первом
приближении.
Найдем значения остатков Rig для всех узлов. Там, где величина
Rij окажется наибольшей, температуры были выбраны наименее удач-
но, и именно для этого узла нужно их скорректировать. Для этого наи-
большее по абсолютному значению делим на четыре части и резуль-
тат добавляем к температуре этого узла, а также к остаткам четырех
соседних узлов.
Как это видно из приведенной выше формулы (4.3), в результате
этого значение R^- в выбранном узле станет равно нулю. Но из-за из-
менения температуры в узле с номером z, j настолько же изменятся
остатки в соседних узлах. Чтобы понять это, запишем остаток для узла,
лежащего слева от остатка с наибольшим значением Rig:
В приведенной формуле только величина t# стала больше на чет-
вертую часть остатка Я^/4. Поэтому и остаток в этом узле станет
больше именно на такую же величину. Теперь понятно, почему алго-
ритм сведения остатков к нулю включает и увеличение остатков в со-
седних узлах на величину RiJ^.
Изменив температуру и остатки в соседних узлах, вновь про-
сматриваем все остатки, снова выбираем узел, где остаток наибольший и
повторяем процедуру сглаживания остатков в этом узле. Повторяя такое
73
сглаживание до тех пор, пока все остатки не станут равными нулю (точ-
нее - некоторой относительно небольшой величине), приходим к реше-
нию задачи.
Модификацией этого метода, наиболее удобной для реализации
на ЭВМ, является метод Зейделя, где уменьшение остатков ведется не
для узлов с наибольшим R, а поочередно, от первого узла к последне-
му. При этом для расчета температуры в последующем приближении
используют значение температуры в том же узле, но рассчитанное в
предыдущем приближении. Более подробно описание этого метода
можно найти в соответствующей литературе.
4.2.	ЧИСЛЕННОЕ РЕШЕНИЕ НЕСТАЦИОНАРНЫХ
ЗАДАЧ ТЕПЛОПРОВОДНОСТИ
Еще раз отметим, что все численные методы основаны на допу-
щении возможности без внесения существенных погрешностей заме-
нять непрерывный в пространстве и во времени процесс некоторым
дискретным, скачкообразным процессом. При расчетах нестационар-
ных процессов, в частности, считают, что температура в любой точке
тела в течение некоторого промежутка времени остается неизменной, а
в начале каждого следующего промежутка меняется скачкообразно,
принимая новое значение.
Для решения задач обычно и здесь используется метод сеток, те-
перь уже пространственно-временных. В качестве примера рассмот-
рим процесс нестационарной теплопроводности неограниченной пла-
стины при ГУ-1 с обеих сторон. В этом случае температурное поле
является одномерным: температура меняется только по толщине пла-
стины. В общем случае принимаем, что температуры поверхностей гс1
и tC2 не постоянны, а заданы как некоторые функции времени:
Мысленно разделим пластину на несколько тонких слоев толщи-
ной Аг и будем считать, что в пределах каждого слоя температура ха-
рактеризуется величиной tifk, где i - номер слоя (/ = О, 1,2, ..., п), к - но-
мер текущего интервала времени Ат (к = 0, 1, 2, ..., т). Для некоторого
момента времени т = Ат • к действительное распределение температу-
ры в теле заменим ступенчатым распределением ее по отдельным сло-
ям (см. рис. 4.3). Тогда температурное поле, представляющее совокуп-
ность всех значений температур tifk, можно отразить при помощи про-
странственно-временной сетки с шагами Ат и Ат, представленной на
рис. 4.4.
74
Рис. 4.3. Температурное поле
при сложных ГУ в момент времени г
Рис. 4.4. Пространственно-
временная сетка
Значения температур в некоторых узлах этой сетки известны по
условиям однозначности.
ГУ слева:	= /ДЛтА:), к= 0, 1, 2, 3,..., т.
ГУ справа: tnk - ,к = 0, 1,2,3,т.
НУ:	tiQ = /3(Ах/), i = 1, 2,	— 1.
Дифференциальное уравнение теплопроводности одномерного
температурного поля при переходе от бесконечно малых к малым ко-
нечным приращениям принимает вид
Здесь а - коэффициент температуропроводности материала пла-
стины;
В результате алгебраический аналог дифференциального уравне-
ния будет
h, £+1 h, к
Ат
(4.4)
75
Приведенная формула связывает между собой четыре соседние в
пространстве температуры по схеме, приведенной на рис. 4.5, которую
называют расчетным шаблоном явной схемы. Если начинать расчет с
i = 1 и к = 0, то эта схема (или формула (4.4), соответственно) будет
содержать только одну неизвестную /1д. Определив ее, шаблон сдви-
гают вправо на шаг и рассчитывают следующую температуру 2 и т.д.
В итоге последовательно определяются все температуры сначала пер-
вого временного слоя (к = 1), затем второго (к = 2) и т.д. Недостатком
явной схемы является необходимость определенным образом соразме-
рить шаги Лт и Дх, поскольку доказано, что решение бывает устойчи-
вым только при условии
Конечно-разностный аналог можно построить и для другого рас-
четного шаблона (неявная схема), приведенного на рис. 4.6. В такой
схеме используются значения температур в соседних точках, но не для
одного, а для двух соседних интервалов времени. Алгебраический ана-
лог дифференциального уравнения при этом принимает вид
Ат
(4.5)
Отметим, что формула (4.5) по структуре идентична формуле
(4.4), но правая часть ее записана для (#+1)-го интервала времени. Для
решения задачи формулу (4.5) записывают последовательно для всех
узлов сетки, в результате чего получается замкнутая система алгебраи-
ческих уравнений, которую решают обычно методом прогонки.
Чтобы понять суть и особенности этого метода, проведем сначала
несложные преобразования, записав формулу (4.5) в виде
$i,k+l $i,k (А+1Л+1
4,k+l
4-i,a+i 4,a+i	4+i,z±i
о—
4-1Л
e-------о
4л 4+1Л
Рис. 4.5. Расчетный шаблон
явной схемы
Рис. 4.6. Шаблон неявной схемы
и сгруппировав подобные члены по возрастанию номера /:
4А-1, &4-1 +	£+1 + Q^z+1, к+1	’
(4.6)
где ради сокращения записей введены обозначения
с - п -
’ '-'i ~ ~ 2 5	~ к •
Ах
В качестве примера по соотношению (4.6) выпишем систему
уравнений для временного слоя при к = 1 и п = 5, расставляя темпера-
туры с одинаковыми индексами друг под другом и дополняя соответ-
ствующие уравнения нулями:
i - 4:	0+	0+	0+	Л4/3д + B4t4l +
= А;
= D2,
= d3;
= d4.
Отметим, что величины Zo,i и г5д известны из граничных условий,
известны и значения О1? D2, D$, D4, поскольку температуры Zo,o, 4),ь to,2?
to,39 и 4),4 известны из начальных условий. Слои с номерами i = 0 и i = п
неизвестных не содержат, поскольку эти температуры известны из
граничных условий. В результате убеждаемся, что записанные четыре
уравнения содержат четыре неизвестных: ^д, f2.i, 6,1 и ^4,ь Аналогичная
трехдиагональная система может быть получена для любого значения
п, а последовательное решение таких систем для к = 1, 2, ... позволяет
найти значения температур во всех узлах сетки. Доказано, что неявная
схема всегда устойчива и это позволяет делать расчеты с достаточно
крупными шагами Ах и Ат.
Метод прогонки основан на допущении, что между соседними уз-
лами температура меняется по линейному закону. Тогда для двух со-
седних (в пространстве) точек сетки можно записать
,к 4-1 Ц 4-1, к 4-1 ”1”	>
(4.7)
где Ej и Fj — некоторые неизвестные пока константы (их называют про-
гоночными коэффициентами), свои для каждого слоя с номером /. Для
предыдущего по порядку слоя с номером i - 1 формула (4.7) запишется
в виде
t'i—1, к4-1	к4-1 "I" 7^—1 •
Подставим теперь это значение температуры в формулу (4.6):
77
(A^f-1	+ ^ih,k+l *" ^4+1, Л+1 = &i
и сгруппируем подобные члены, переписав предыдущую формулу так:
^л+1(А+ АА-1) - “GA+u+i + А _ 4Л-1 •
Отсюда находим
Сопоставляя теперь формулы (4.7) и (4.8), видим, что прогоноч-
ные коэффициенты можно рассчитать по рекуррентным формулам:
£.=-Cz./(Bf + 4£f_1);
^.=(d/-4F/_1)/(4 + 4A-i)-
Заметим, что значения Е\ и F\ легко найти, записав формулу (4.6)
для i = 1 и к = 1:
и преобразовав ее аналогично предыдущему:
Из последней формулы понятно, что
£i = -с} /4 и д = (рг — 4/од)/4 •
Подчеркнем, что все величины, входящие в правые части этих
формул, известны и это позволяет рассчитать численные значения ко-
эффициентов £i и Fi.
Далее по рекуррентным соотношениям рассчитывают значения
прогоночных коэффициентов £( и F, для первого временного слоя
(прямая прогонка: i = 2, 3,..., п), а после этого по формуле
^i-\,k ~ ^i-l^i,k+l + ^1-1
находят температуры во всех узлах этого временного слоя, начиная от
предпоследнего (обратная прогонка: i = и, п - 1, ..., 2, 1). В начале об-
ратной прогонки приведенная формула кроме рассчитанных уже прого-
ночных коэффициентов содержит известную из граничных условий
температуру при этом на каждом последующем шаге используется
результат предыдущего расчета. Закончив расчет первого временного
78
слоя (к = 1), начинают прямую прогонку для следующего (к = 2), а рас-
считав значения прогоночных коэффициентов - обратной прогонкой
рассчитывают температуры и для этого слоя. Повторяя описанный
алгоритм, полностью решают задачу. Итеративный алгоритм метода
прогонки ориентирован на широкое применение ЭВМ для выполнения
расчетов.
Заметим, что метод прогонки - это модификация метода Гаусса
применительно к трехдиагональной системе линейных алгебраических
уравнений.
В теории разностных схем рассматривается еще и линейная ком-
бинация явной и неявной схем с некоторыми весовыми коэффициен-
тами а и (1 — а):
h,k+l ^,к	. 2	2^-£+1 + д+i) + (1 a)(/f+u	+ ^--i £)].
При а = 0,5 такую схему называют схемой Кранка - Николсона.
При получаем неявную схему, при о = 0 - явную. Эта схема об-
ладает определенными преимуществами, в частности более высоким
порядком точности по сравнению даже с неявной схемой, однако при
увеличении шага Лт могут возникать колебания разностных решений.
4.3.	РЕШЕНИЕ МНОГОМЕРНЫХ
ЗАДАЧ ТЕПЛОПРОВОДНОСТИ
В этом случае тело разбивается на элементы простой геометриче-
ской формы. В простейшем случае это прямоугольные параллелепипе-
ды или кубы для трехмерных тел или прямоугольники и квадраты -
для плоских фигур.
Среди множества более сложных подходов наибольшее распро-
странение находит так называемая локально-одномерная схема. При
этом протекание многомерного процесса распространения теплоты
представляется условно как двух- или трехстадийный процесс: на пер-
вой стадии за время Лт теплота распространяется только вдоль оси х, а
на второй стадии в течение этого же промежутка времени - вдоль оси
у (для плоской задачи). Этот прием, называемый расщеплением реше-
ния, обеспечивает достаточно точное решение, если применять метод
прогонки для решения по каждому из направлений.
На рисунке 4.7 приведена схема, наглядно иллюстрирующая суть
метода расщепления решения. Считая скорость распространения теп-
лоты очень высокой, представим, что на первой стадии между гори-
зонтальными слоями тела установлены теплоизолирующие прослойки,
в результате чего тепловой поток распространяется только по направ-
лению х до самой наружной границы.
79
Рис. 4.7. Физическая интерпретация метода расщепления решения
для двумерных задач
После этого такие же прослойки укладываются между вертикаль-
ными слоями, и поток проходит по другому направлению. Из-за боль-
шой скорости оба эти процесса успевают закончиться за время Ат, а
результирующая температура в любом узле состоит из суммы рассчи-
танных по каждому из направлений значений температур.
По сути дела при расщеплении решения дифференциальное урав-
нение теплопроводности (для плоской задачи)
заменяется эквивалентной системой из двух уравнений в частных про-
изводных
Алгебраические аналоги для направлений х и у будут иметь сле-
дующий вид:
аЛт
Д/
((/+1Л+1
^уД+1 +G-l,i+l) >
О, £+!
где z, j - индексы, определяющие номер узла сетки по координатам х и
у, соответственно. На рис. 4.8 приведен расчетный шеститочечный
шаблон, используемый для решения рассматриваемой задачи.
$0
^/+1,£+1
О
tij,k
Рис. 4.8. Расчетный шаблон неявной схемы для двумерной задачи
Решается эта система методом прогонки поочередно по каждому
из направлений (/ = О, 1, 2, ..., n‘,j = 0, 1, 2, ..., т) для каждого времен-
ного слоя (к = 1, 2, 3,	I). Определенные сложности возникают толь-
ко при формулировании и использовании начальных и граничных ус-
ловий. Обычно известны начальные температуры и это позволяет
методом прогонки рассчитать все температуры ^дд и ^д, а затем,
сдвинув шаблон, и все остальные температуры в узлах пространствен-
но-временной сетки.
4.4.	РЕШЕНИЕ НЕЛИНЕЙНЫХ ЗАДАЧ ТЕПЛОПРОВОДНОСТИ
Известно, что у многих материалов, особенно капиллярно-
пористых, величины X, р и ср заметно зависят от температуры. Обычно
эти зависимости находят опытным путем, а опытные данные аппрок-
симируют простыми расчетными формулами. Например, величину
коэффициента теплопроводности обычно описывают линейной зави-
симостью:
~I-
где Хо - теплопроводность материала при t = 0 °C, а коэффициент а
показывает, насколько изменяется величина X при изменении темпера-
туры на один градус.
Нелинейность обязательно приходится учитывать при расчетах
задач, в которых в течение процесса температура тела меняется в
очень больших пределах. Например, теплопроводность хромоникеле-
вых нержавеющих сталей, используемых в криогенной технике, при
изменении температуры от 5 К до 300 К изменяется от 1 до
15 Вт/(м-К).
Если при расчетах учитывать эти зависимости, то дифференци-
альное уравнение теплопроводности становится нелинейным. Анали-
тическое решение нелинейных дифференциальных уравнений в част-
81
ных производных даже для одномерных задач практически невозмож-
но или (при определенных допущениях) отличается большой сложно-
стью. Поэтому эти задачи обычно решают численным способом.
Система алгебраических уравнений, получаемая на основе алгеб-
раического аналога, в этом случае становится тоже нелинейной. Для
решения используется тот же подход, который разработан для обыч-
ных нелинейных уравнений. Наиболее распространенный прием со-
стоит в том, что реализуется квазилинейная модель, суть которой за-
ключается в том, что значения величин X, р и ср и других (например qv)
при расчетах каждого последующего временного слоя рассчиты-
ваются предварительно по известным уже температурам для преды-
дущего момента времени. Это дает возможность рассчитать величину
в первом приближении, а затем, используя полученную темпера-
туру, рассчитать значения X, р и другие во втором приближении и да-
лее - во втором приближении. Такой внутренний цикл итераций
повторяют до тех пор, пока предыдущие и последующие приближения
температуры практически перестанут отличаться. После этого все
повторяется на следующем пространственном шаге (при i + 1) до тех
пор, пока i < п.
Чтобы лучше понять изложенный алгоритм решения квазилиней-
ных задач, на рис. 4.9 приведена его укрупненная схема для расчета
одномерного нестационарного температурного поля (без внутреннего
цикла уточнения температуры
Ввод исходных данных (начальные и
граничные условия)
Расчёты л, р и с по температуре tirk.}
*
Расчёт прогоночных
коэффициентов Е: и
I
Обратная прогонка,
расчёт температур tn_iik
нет
Рис. 4.9. Блок-схема алгоритма решения по квазилинейной модели
82
При более сложных ситуациях на каждом временном шаге Ат ме-
тодом последовательных приближений с наперед заданной точностью
решается система нелинейных алгебраических уравнений (основы та-
кого решения рассмотрены нами ранее).
4.5.	МЕТОД КОНЕЧНЫХ ЭЛЕМЕНТОВ
Выше мы достаточно подробно познакомились с одним из подхо-
дов, широко применяемых при численном решении задач теплообме-
на. Это так называемый метод сеток, суть которого состоит в том, что
исследуемое пространство мысленно представляется как совокупность
большого числа элементарных параллелепипедов с весьма малыми
сторонами Ах, Ау и Az, а непрерывный в пространстве и по времени
процесс распространения теплоты внутри тела заменяется некоторым
скачкообразным. При этом в течение отрезка времени Ат температура
внутри элемента считается неизменной, а изменение температуры
происходит в момент перехода к следующему отрезку Ат. Изменение
температуры в пространстве также происходит скачками на гранях
между соседними элементами. Температуру каждого элемента условно
сосредотачивают в центре элементов, и соединяя эти центры прямыми,
получают прямоугольные сетки в проекции на соответствующие коор-
динатные плоскости. Для поиска решения в дифференциальных урав-
нениях, описывающих процесс, элементарно малые приращения dx, dy,
dz, dT и dx также заменяют выбранными конечными приращениями Ах,
Ау и т.д., тем самым превращая дифференциальные уравнения в алгеб-
раические, для решение которых известны регулярные общие методы.
Однако при таком подходе определенные трудности возникают
при описании граничных условий, а также с проблемой сходимости,
устойчивости и аппроксимируемости результатов численного решения.
Действительно, при неудачном выборе размеров Ах, Ay, Az и отрезка
времени Ат результаты расчетов могут совершенно не соответствовать
действительному распределению температуры в теле или при уменьше-
нии названных величин результаты расчетов начинают то приближать-
ся, то отдаляться от действительного решения. Конечно же, разработаны
отдельные рекомендации, обеспечивающие получение правильных ре-
зультатов, но и при этом остаются определенные, а иногда и непреодо-
лимые трудности при реализации названного подхода.
В настоящее время очень широкое применение находят вариаци-
онные методы решения задач, суть которых заключается в том, что
решается задача минимизации некоторого специально разработанного
функционала 7[/(х)]. В качестве такого функционала обычно выступает
величина определенного интеграла, где подынтегральной функцией
является искомая зависимость fix). В отличие от функции, дающей
непрерывное множество значений, функционал отличается тем, что
83
результат определения функционала всегда только одно число. Вели-
чина этого числа зависит от вида функции fix) и пределов интегриро-
вания. Именно это свойство функционала и позволяет использовать
его в качестве критерия оптимизации, например.
В какой-то мере подход на основе вариационного метода мы
встречали при описании метода релаксаций, используемого для реше-
ния задач стационарной теплопроводности. Там для некоторой пло-
ской фигуры записывалось известное дифференциальное уравнение
стационарной теплопроводности
82Т д2Т
дх2 dv2
после чего тело было разбито на прямоугольники со сторонами Дх и
Ду, и для некоторого центрального узла с координатами i и j был полу-
чен алгебраический аналог дифференциального уравнения в виде сле-
дующего соотношения между температурами выделенной точки и
температурами в соседних с ней узлах (при Дх = Ду):
Ti+ij+ TiAj+ TiJ+i + Tjj.i - 4Tjj- 0.
Далее для каждого узла сетки задавали приближенные значения
температур Т, при этом, конечно же, приведенное равенство не соблю-
далось, поэтому для каждого узла находились величины остатков
-4Г
а затем специальным алгоритмом минимизировались величины 7?^.
Если просуммировать все R^, то полученная сумма может быть
принята за функционал I, и тогда можно будет решать задачу опреде-
ления минимума обычным образом, т.е. записывая и решая уравнения
д!
дх
поскольку это позволяет определить экстремальное (в нашем случае
min) значение суммы остатков и тем самым найти искомые температу-
ры Tij.
Главным преимуществом метода конечных элементов является
то, что здесь сняты отдельные ограничения на принципы построения
сеток. Это позволяет более точно описывать наружную конфигурацию
исследуемых тел, в необходимых местах увеличивая или уменьшая
густоту узлов. Упрощается решение нелинейных задач, когда физиче-
ские свойства тел или граничные условия оказываются зависимыми от
координат расположения соответствующих узлов пространственной
сетки или от номера временного интервала. Проще решаются и вопро-
сы устойчивости и сходимости решений.
84
Чтобы разобраться с основами метода конечных элементов, рас-
смотрим задачу о стационарном температурном поле в тонкой пласти-
не сложной геометрической формы при наличии внутренних источни-
ков теплоты и граничных условиях третьего рода (рис. 4.10). Поверх-
ность пластины обозначена через D, а длина по контуру - L. Измене-
нием температуры по толщине пластины пренебрегаем, считая темпе-
ратурное поле двумерным.
Дифференциальное уравнение теплопроводности при наличии
постоянно действующих внутренних источников тепла производи-
тельностью qv запишется в виде трех слагаемых (через Т обозначена
избыточная температура: Т = t-
а граничные условия третьего рода описываются известным диффе-
ренциальным уравнением, отражающим закон Ньютона-Рихмана:
Величины X, а и qv могут быть функциями от координат хи у.
Задача формулируется как обычно: необходимо найти такую
функцию Т = fix, у), которая будет наилучшим образом удовлетворять
приведенным дифференциальным уравнениям в любой точке, при лю-
бых хи у, принадлежащих плоскости D.
Опуская пока технологию получения подынтегральных выраже-
ний (в принципе она аналогична описанной нами выше), запишем
функционал, минимизация которого и позволит решить поставленную
задачу:
4л*,х>]=
(4.9)
Чтобы решить задачу численным способом, пластину разбивают
на множество элементов простой геометрической формы. Чаще всего -
это треугольные симплекс-элементы, хотя в качестве таких элементов
могут выступать четырехугольники (необязательно прямоугольные) и
фигуры более сложной формы. Выработаны некоторые рекомендации
для того, чтобы на самых первых этапах способствовать упрощению и
облегчению будущего решения. В частности желательно, чтобы при
разбивке поверхности D треугольные симплек-элементы были близки
по форме к равносторонним треугольникам. Пример такого разбиения
показан на рис. 4.11. В произвольном порядке пронумеруем все эле-
менты от 1 до N. При этом мы получили Музловых точек (М= N + 2).
85
Рис. 4.10. Плоская стенка
сложной формы
X
Рис. 4.11. Пример разбивки тела
на треугольные элементы
На рисунке 4.12 показан один из типовых симплекс-элементов с
номером п и глобальными номерами узловых точек тп, т + 1 и т + 2
(на предыдущем рисунке он выделен штриховкой). Каждый такой тре-
угольник содержит 3 узла, и именно в этих узлах будем определять
температуру Т(х, у), добиваясь того, чтобы эти температуры мини-
мально отличались от истинных, полностью соответствующих приве-
денному математическому описанию температурного поля.
Порядок глобальной нумерации узлов рекомендуется проводить
так, чтобы максимальная разница R = №тах - №min между наибольшим
номером узла первого элемента и наименьшим номером того же эле-
мента была минимальна. На рисунке 4.13 приведены 2 варианта нуме-
рации узлов.
При нумерации вдоль наименьшей из сторон элемента (рис. 4.13, а)
величина R = 7, а при другой последовательности нумерации узлов
(рис. 4.13, б) 7? = 21. Первый вариант предпочтительнее, поскольку это
впоследствии дает более удобную для расчетов матрицу.
86
Рис. 4.13. Примеры разбиения фигуры на симплекс-элементы
Для каждой узловой точки искомую температуру Т(х, у) предста-
вим в виде линейного многочлена
м
T{x,y) = ^Amfm{x,y),
m=l
(4.Ю)
где Am — некоторые неизвестные пока постоянные коэффициенты, а
fm(x, у) - некоторые функции пространственных координат (функции
формы), которые ниже мы научимся определять для любого узла по-
лученной сетки.
Чтобы определить эти функции fm(x, у), рассмотрим выделенный
нами треугольник (см. рис. 4.12), и специальным образом (всегда про-
тив часовой стрелки) промаркируем принадлежащие этому треуголь-
нику узлы локальными номерами z, j и к. Узел i называют главной уз-
ловой точкой симплекс-элемента. На рисунке приведен и глобальный
номер тп выделенного узла. При этом выбор начальной точки i произ-
волен.
Для треугольного элемента функция формы любого узла с прием-
лемой точностью может быть представлена простой линейной зависи-
мостью
87
fm{x, y) = a+ bx + cy,
где постоянные коэффициенты a, b и с легко находятся через коорди-
наты узлов к с учетом того, что в узловой точке i значение коэффи-
циента формы принимается равным единице fm(x, y)j = 1, а для других
двух точек - нулю :/„(%, у), = 0,fm(x,y)k = 0. Такое распределение значе-
ний функции формы аналогично выделению в методе релаксации тем-
пературы tij, которая последовательно уточняется с учетом температур
ti+\j и др. в соседних точках. Два другие узла этого элемента вы-
ступят в качестве главных узловых точек в соседних треугольных эле-
ментах, соприкасающихся с выделенным.
Итак, запишем функции формы для главной узловой точки i вы-
деленного элемента с номером п и для двух других узлов (локальные
номера у, к)\
Clm "I" b^i "И СтУj 1,
(4.П)
"1“ ЬпрСк "И СтУк 0*
Поскольку координаты Xj, хк и yi9 yj, ук нам известны по резуль-
татам разбиения поверхности D на симплекс-элементы, мы видим, что
получили замкнутую систему линейных уравнений, содержащую 3
неизвестных: ап, Ьп и сп, которые найдем, решая систему (4.11):
ат = (Х]Ук ~ Xjyj)/2S‘
Ьщ (yj У
(4.12)
Cm (Хк Xj)/2S9
где 5 — площадь выделенного симплекс-элемента. Эту площадь легко
рассчитать через координаты вершин треугольника, вспомнив сле-
дующую формулу из аналитической геометрии:
Sm = (ху* - X&J + xyj - хук + Xjyt - Xjy^/2.
Определив описанным способом константы ат, Ьт и ст для всех М
узлов области D, по формуле (4.10) заменим функцию Т(х, у) в функ-
ционале (4.9). Тогда получим новый функционал, зависящий только от
констант Ат. При этом, учитывая, что величины fm(x, у) равны 1 для
узла с номером тп и равны 0 для всех других узлов, можно считать, что
Ат ~ Т(х, у).
Тогда задача минимизации нового функционала решается обыч-
ным способом: его производную приравниваем нулю и, решая полу-
ченное уравнение, находим экстремальное значение (в нашем случае -
минимальное) функционала. Поскольку наш функционал зависит от
нескольких величин Ат, то следует рассматривать частные производ-
ные:
88
di
дЛг
= 0,
После дифференцирования мы получим М уравнений, решив ко-
торые можно определить все коэффициенты Ат, а значит и определить
все температуры Т(х, у) в узлах сетки, т.е. решить задачу теплопровод-
ности.
Поняв суть излагаемого метода, попробуем разобраться с задачей
подробнее на простейшем примере. Возьмем фигуру, состоящую всего
из трех симплекс-элементов, приведенную на рис. 4.14. Пронумеруем
элементы так: (1), (2), (3) и узлы 1, 2, 3, 4, 5. Кроме того, нанесем ло-
кальные номера узлов в каждом элементе (/, у, к). Для перехода к гло-
бальным координатам составим следующую таблицу (в скобках при-
ведены координаты х и у узлов каждого из элементов).
0	12	3	4
Рис. 4.14. Фигура из трех симплекс-элементов
Таблица 4.1
№ элемента	Номера узлов		
	Локальные		
	/(/', к)	Лк, 0	k(i,f)
	Глобальный номер узла, координаты х; у узла		
(1)	1(1,3; 1)	2(2,9; 0,5)	3(2,4; 1,7)
(2)	2(2,9; 0,5)	4(3,6; 1,7)	3(2,4; 1,7)
(3)	4(3,6; 1,7)	5(3,2; 2,9)	3(2,4; 1,7)
Для элемента (1), записывая систему (4.12) для каждого узла, как
для узловой точки (поскольку выбор узловой точки произволен), полу-
чим три следующих системы:
^(i),z = (*<1)уУ(1)л - ^(1),^(1)Л/2^;	Л(1),/К*(1)л J'O),* - X(i)jZ- у(ц *)/2S;
= (y(i)j -у(1)Л)/25;	fyi) j =О(1)Л - у(1)9/)/25;
C(l),z = (X(l)J - Х(1)Д/25,	C(i)j =(x(l),z “ x(i),
89
°(i), t= (x(i),/ Ут j - xmjym^/2S;
b(D,k= (ym,i ~ ymj)/2S'^
C(i), *= (Xi),; - ^(i),<)/25,
Совершенно аналогично можно записать по 3 системы для эле-
мента (2) и для элемента (3), получив соответствующие формулы для
а<2),ь	b(2),h b(3),iи т.д. В результате мы получили переопределенную
систему, содержащую 3x3x3 = 27 уравнений. А для пяти точек доста-
точно 3x5 = 15 уравнений, чтобы система была замкнутой.
Обратим внимание на рис. 4.14. Узел 2 принадлежит двум разным
элементам (1) и (2). В первом он имеет локальный номер у, а во втором — /.
Чтобы при расчетах по элементу (1) и по элементу (2) получать одина-
ковую температуру, в узле 2 коэффициенты ат, Ьт, ст и а(2)Л 6(2),/,
C(2)jZ- должны быть равными друг другу. Складывая правые и левые час-
ти соответствующих формул и разделив сумму пополам, из 6-ти урав-
нений получим только 3. Совершенно аналогично уменьшим и число
уравнений для узла 4. А для узла 3 придется сложить по 3 правых и
левых части формул для	^(2), ъ с<2), к и ^(3), ь
С(з)^, и суммы разделить на 3. Тогда для каждого узла действительно
получаем по 3 уравнения для вычисления коэффициентов для функции
формы. Действительно, складывая формулы для коэффициентов а, Ъ и
с в узлах 2 и 4 мы уменьшим число уравнений на 6. Делая то же самое
для узла 3, мы еще на 6 уравнений уменьшим полученную систему. В
результате останется замкнутая система, содержащая 15 уравнений.
Если посмотреть на рис. 4.13 (схема а), то легко найти такие узлы
(один из них выделен жирно), которые принадлежат четырем, пяти,
шести или даже семи симплекс-элементам. Поэтому после суммирова-
ния соответствующих формул и перехода к канонической форме сис-
темы уравнений мы можем встретить в литературе множители, стоя-
щие перед соответствующими суммами 1/4, 1/5, 1/6 или 1/7 и др.
Теперь нам осталось только вместо локальных номеров i,j, к под-
ставить глобальные номера из приведенной выше таблицы и значения
координат соответствующих точек. После этого следует привести сис-
тему к канонической форме. Полученную замкнутую систему обычно
представляют как матричное уравнение
677= Ф,
где G - глобальная матрица теплопроводности размером М х М, со-
держащая значения частных производных функционала I по темпера-
турам в каждом узле; U - вектор-столбец искомых значений темпера-
тур в узлах 1.. .М\ Ф - глобальный вектор-столбец тепловых потоков
90
размерностью состоящий из значений коэффициентов формы
для каждого узла.
Для решения этого уравнения можно применять и стандартный
метод (мы его осваивали в пакете Excel), а также специальные методы,
ускоряющие расчеты.
В заключение отметим, что несмотря на определенную сложность
для понимания, изложенная методика хорошо формализована, и это
позволяет легко перекладывать все описанные действия на плечи ПК.
В настоящее время на рынке имеются предложения ряда хороших па-
кетов прикладных программ, реализующих автоматически и разбиение
заданной области на конечные элементы, и формирование соответст-
вующих вектора-столбца Ф и матрицы G, и решение полученного мат-
ричного уравнения как для плоских стационарных и нестационарных
задач, так и для трехмерных задач. Примером таких пакетов могут
служить пакет Elcut для плоских задач и пакет Ansys для трехмерных
объектов, а также пакет COMSOL Multiphysics.
4.6.	КОНЕЧНО-РАЗНОСТНЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ
КОНВЕКТИВНОГО ТЕПЛООБМЕНА
Чтобы описать и найти температурное поле в движущейся жидко-
сти аналогично дифференциальному уравнению теплопроводности на
основе закона сохранения энергии выводится специальное дифферен-
циальное уравнение - дифференциальное уравнение энергии [8]. Это
уравнение учитывает и перенос тепла теплопроводностью, и накопле-
ние тепла в элементарно малом объеме в результате изменения его
теплосодержания при протекании через него потока теплоносителя.
По форме оно похоже на дифференциальное уравнение Фурье
однако содержит дополнительное слагаемое в левой части - конвек-
тивную составляющую, отражающую вклад конвекции в общий тепло-
вой баланс. Из уравнения видно, что температурное поле в движущей-
ся жидкости (или газе) зависит от ее скорости, и при решении тепло-
вой задачи необходимо одновременно решить задачу гидродинамиче-
скую, т.е. найти скорость жидкости w в любой точке потока.
Если к элементарно малому объему движущейся жидкости при-
менить известный из механики принцип Даламбера (сумма всех сил,
действующих на тело, включая и силу инерции, взятую с обратным
знаком, равна нулю), то можно получить дифференциальное уравне-
ние движения, которое в символьной форме можно представить в виде
91
dw dw 1	„2
—+ w—= g— Vp + vV w,
CT on p
где Vp - оператор Гамильтона; g - ускорение силы тяжести; V2w -
оператор Лапласа. Левая часть этого уравнения отражает собою силу
инерции (сумма локальной и конвективной составляющих), в правой
части величина g - отражает силу веса, слагаемое с Vp - силу давле-
ния, а последнее слагаемое — силу трения. Как видим, уравнение дви-
жения содержит два неизвестных параметра: w и р. Поэтому уравнение
движения всегда рассматривается совместно с другим уравнением
гидродинамики - уравнением неразрывности.
Дифференциальное уравнение неразрывности отражает закон со-
хранения массы применительно к бесконечно малому объему движу-
щейся жидкости. В общем случае оно имеет вид
d(pw) = 0
дп
Для несжимаемых жидкостей (р = const) это уравнение упрощается:
dw „ dwv	dw
— = 0 или —— ч-----— ч---— = 0.
dn	dx dy dz
В дополнение к приведенным дифференциальным уравнением
необходимо добавить еще дифференциальное уравнение теплоотдачи,
описывающее особенности теплообмена на границе жидкость - стенка:
= а(Г-Гс).
В итоге для несжимаемых жидкостей приведенные дифференци-
альные уравнения составляют замкнутую систему, содержащую четы-
ре неизвестных: a, t, w и р. Для газов и сжимаемых жидкостей величи-
на р тоже войдет в список неизвестных и система уравнений должна
быть дополнена еще уравнением состояния:
pv - ZRT или р / р = ZRT,
где Z — общий коэффициент сжимаемости; R — газовая постоянная.
Приведенная система уравнений описывает весь класс явлений
конвективного теплообмена. Чтобы решить некоторую конкретную
задачу, необходимо проинтегрировать уравнения, учитывая еще и ус-
ловия однозначности этой конкретной задачи. Формулирование этих
условий гораздо сложнее, чем в задачах теплопроводности. Так, началь-
92
ные и граничные условия, например, должны быть заданы для каждого
неизвестного параметра, а не только для температуры.
Из-за сильной зависимости вязкости от температуры уравнение
движения является нелинейным. Другие дифференциальные уравнения
тоже достаточно сложные. Поэтому аналитическое решение задачи пу-
тем интегрирования системы дифференциальных уравнений в общем
случае невозможно. Даже для самых простых задач, чтобы получить
аналитическое решение, приходится вводить массу упрощающих пред-
посылок (например, что кинематическая вязкость v = const, движение
равномерное и т.д.), которые в итоге делают полученное решение при-
ближенным и малодостоверным.
Поэтому приведенные дифференциальные уравнения обычно ис-
пользуют для численного решения задач конвективного теплообмена.
Именно на их основе строятся конечно-разностные аналоги для расче-
тов методом сеток. Большинство же важнейших для практики задач
решены на основании экспериментальных исследований с привлече-
нием для организации опытов и обработки результатов этих экспери-
ментов основ теории подобия.
Чтобы понять, как осуществляется переход от дифференциальных
уравнений к их алгебраическим аналогам, используемым при числен-
ном интегрировании, рассмотрим простейшее установившееся одно-
мерное течение вдоль оси х. Будем считать, что тепло физические
свойства с, 1 и р жидкости в течение процесса остаются постоянными
в любой точке потока, скорость w любой элементарной струйки посто-
янна и известна (например, определена экспериментально) и темпера-
турное поле стационарно.
В этом случае дифференциальное уравнение энергии примет вид:
(4.13)
Для построения алгебраического аналога этого уравнения приме-
ним другой подход, называемый методом баланса. Суть метода опре-
деляется тем положением, что при стационарном температурном поле
в элементарном объеме не может происходить накопление теплоты,
иначе температура его увеличивалась бы с течением времени.
Рассмотрим элементарную струйку, текущую слева направо, вы-
делив в ней элементарный объем длиной Дх (см. рис. 4.15). Как и
прежде, будем считать, что температура любого элемента сосредото-
чена в его центре и равняется 7}. Вдоль оси х теплота передается двумя
способами: конвективным путем при течении жидкости и за счет теп-
лопроводности вдоль выделенного слоя.
93
Рис. 4.15. Расчетный элемент X;
Накопление теплоты в элементе от конвекции определится разни-
цей энтальпий на выходе и входе в элемент:
^вых (через сечение при xt+ Ах/2) - <?*х (через сечение при х~ Ах/2).
Накопление теплоты за счет теплопроводности найдется анало-
гично
^вых (через сечение при х,+ Дх/2) - ^х (через сечение при х,— Дх/2).
В методе баланса для стационарного температурного поля сумму
двух приведенных выражений следует приравнять нулю:
( ^ых - ^ ) + ( <7вых - <7вх ) = 0.	(4.14)
Аппроксимацию потоков д*ых и можно сделать разными спо-
собами (через Ui обозначено приближенное значение температуры 7}):
- со сдвигом по ходу потока
*7вых —	j *?вх —	(4.15)
- по средней температуре
^Кых=сри<Ц.+1+Ц)/2;	= cp</z. +	) / 2;	(4.16)
- со сдвигом против хода потока
^вых = cpwVi q*x = cpwUt_x.	(4.17)
На основании закона Фурье можно записать алгебраические ана-
логи для тепловых потоков, передаваемых теплопроводностью
Аппроксимация (4.15) самая неудачная, поскольку предполагает,
что температура в точке i не зависит от температуры в предыдущем
объеме, а определяется только температурой следующего по ходу по-
94
тока объема. Это противоречит физической сути явления. Ведь буду-
щее всегда туманно, а прошлое хорошо известно. Именно поэтому при
использовании аппроксимации (4.15) часто получаются абсурдные
результаты.
Второй подход не вызывает противоречий, но может плохо рабо-
тать при больших скоростях w. Доказано, что наилучшей является ап-
проксимация (4.17).
12	~
Будем считать, что величина производной — является линеи-
дх
ным многочленом приближенных температур U}:
дТ
дх
+ a2Ut + a^U м.
— щи щ
Эта формула обобщает все три изложенных выше подхода. Тогда
алгебраический аналог дифференциального уравнения (4.13) будет
cpw(a1t//+1 + a2Ut + a3t/M) = X

или (поскольку правую часть приведенной формулы мы уже рассчи-
тывали при анализе дифференциального уравнения теплопроводности)
ср^а1им +a2Ut +а3и^) = K-^(UM. (4.18)
Ах
При аппроксимации центральной разностью щ = Ах/2, а2 = 0 и
аз = -Ах/2 и полученная формула после простых преобразований при-
нимает вид
и _ им + Ц-1 cpwAx
'	2 4л
В общем случае температура Т вдоль потока может совершенно
не меняться, может возрастать, а может и уменьшаться. Это зависит от
интенсивности поперечного теплообмена, который в нашем рассмот-
рении совершенно исключен. Представим, что в результате попереч-
ного подвода теплоты Т монотонно возрастает. Тогда для трех сосед-
них температур должно соблюдаться соотношение:
и^си^и^.
Однако из формулы видно, что при большой скорости w величина
Ui может оказаться меньше что противоречит приведенному не-
равенству. Поэтому, чтобы решение получалось правильным, нужно
95
ввести такое ограничение, при котором требуемое соотношение со-
блюдается. Таким ограничением является условие:
4Л
откуда следует
cpwAx
При аппроксимации разностью назад а\ = 0, аг = Лх/2 и а3 = -Лх/2
и тогда алгебраический аналог после преобразований формулы (4.18)
получаем в виде
им + Ц 1
cpwAx
X(2 + cpwAx/X)
1 + cpwAx/X
Из формулы видно, что она не требует никаких ограничений и
для рассмотренного выше случая условие возрастания температур со-
блюдается. Именно эта формула чаще всего и используется при расче-
тах.
Значение i последовательно изменяется от 1 до /кон, что дает замк-
нутую систему алгебраических уравнений, решая которую находят все
значения Ц.
Пока мы рассмотрели только одну элементарную струйку (точнее —
некое множество таких струек) внутри ламинарного потока при упро-
щающем предположении о том, что скорость w является постоянной и
не зависящей от любых обстоятельств.
Естественно, что при решении реальных задач необходимо учи-
тывать действительные граничные и начальные условия, распределе-
ние скорости по сечению и другие факторы.
Для примера рассмотрим ламинарное установившееся течение
несжимаемой жидкости внутри трубы (см. рис. 4.16).
^ст
Рис. 4.16. Течение жидкости в трубе
96
Пусть задана постоянная температура стенки Тст и коэффициент
теплоотдачи а (см. дифференциальное уравнение теплоотдачи). В силу
радиальной симметрии картины можно записать еще одно граничное
условие
Если считать, что ТФС, включая и вязкость жидкости, не зависят
от температуры и работа сил давления пренебрежимо мала, то интег-
рируя уравнение движения, можно получить следующее распределе-
ние скорости в потоке:
(4.19)
Дифференциальное уравнение энергии, аналогичное (4.13), но за-
писанное в цилиндрической системе координат, имеет вид:
cpwr
дТ
dz
В абсолютном большинстве случаев переносом теплоты вдоль
потока за счет теплопроводности можно пренебрегать (если критерий
Пекле Ре » 1.0) и тогда приведенное уравнение упрощается, посколь-
ку производная избыточной температуры по z равна нулю

дТ
Отметим, что это уравнение аналогично уравнению нестационар-
ной теплопроводности при одномерном температурном поле, когда
Т = fix, т). Только в нашем случае Т = fiz, г). Поэтому разобьем верх-
нюю часть потока равномерной по г и по z сеткой с шагами Az, Аг
(рис. 4.17), обозначая координаты узлов сетки через i и j. Это позволя-
ет записать алгебраические аналоги уравнения энергии и уравнений
граничных условий.
Рис. 4.17. Пространственная сетка
97
Если считать, что изменение энтальпии на границе потока проис-
ходит только за счет теплообмена со стенкой, то для элемента жидко-
сти, соприкасающегося со стенкой, получаем уравнение, аналогичное
ГУ-3, отражающее тепловой баланс ^кон = ^то*
cpwr
—	- П
где F — поверхность теплообмена, F — 2 тс/? Az; 5 - сечение слоя,
5 = 2л(/?-Аг/2)Аг.
Для внутренних узлов сетки при аппроксимации против потока
методом баланса получаем алгебраический аналог уравнения энергии:
СРЧ UtJ 1 =^^[rM/2('Ui+lJ ’
где величину wr рассчитывают по формуле (4.19), а величина ттн/2 на-
ходится по формуле
Аналогично получается разностная схема для элемента, лежащего
на оси симметрии (г = 0, j = 1):
cpWr=0	(С7Л1 - ц. 0) = Хж (t/2J - UXJ )Az.
Для узлов, соприкасающихся со стенкой трубы (г = R, j = т
wr= 0), из приведенного выше уравнения получим:
Как уже отмечалось выше, приведенные уравнения решаются со-
вместно с использованием явной или неявной схем для каждого из
слоев с индексом j = При этом для каждого элемента с индексом
i приходится решать систему уравнений, содержащую одну (Ц при
явной схеме) или несколько (Ц_ь £4 Ц+i) неизвестных.
4.7. РЕШЕНИЕ СОВМЕСТНЫХ ЗАДАЧ
ПРИ КОНВЕКТИВНОМ ТЕПЛООБМЕНЕ
При расчетах систем охлаждения или обогрева различных техни-
ческих устройств часто встречаются случаи, когда приходится рас-
сматривать такие задачи, в которых требуется одновременно опреде-
лить и температурное поле внутри стенки канала, и температурное
поле внутри движущегося в этом канале потока жидкости. Подходы к
решению каждой из этих задач мы рассмотрели выше на примере за-
98
дачи теплопроводности и задачи конвективного теплообмена при за-
данных граничных условиях. Теперь же предстоит рассмотреть более
сложную задачу, в которой нужно проинтегрировать совместно урав-
нения математической модели, включающие как дифференциальное
уравнение теплопроводности для стенки, так и дифференциальные
уравнения конвективного теплообмена для потока жидкости с учетом
соответствующих условий однозначности.
Чтобы понять основы подхода к численному решению таких за-
дач, используя при этом наработанные нами выше приемы, рассмот-
рим наиболее простой вариант поставленной задачи. Пусть требуется
определить температурное поле внутри тонкостенной трубки из лату-
ни и температурное поле протекающей внутри трубки жидкости, кото-
рую обогревают, пропуская по трубке электрический ток при постоян-
ном напряжении на концах трубки (рис. 4.18).
Сформулируем задачу более точно: в канале, образованном ци-
линдрической стенкой длиной I с площадью сечения 5СТ и смоченном
периметром FCT, протекает жидкость с удельной теплоемкостью ср и
постоянным массовым расходом G. Коэффициент теплопроводности
материала стенки может зависеть от температуры Хст=/1(Тст), а в стен-
ке действует внутренний источник теплоты, для которого задана
удельная линейная мощность qi = f2(x, Тст) (функции ft и f2 известны).
Местные термические сопротивления в местах подсоединения элек-
тродов к трубке равны соответственно 7?Ti и Т?т2.
Температурные поля стенки и жидкости, как и ранее, будем счи-
тать стационарными и одномерными.
Для построения разностной схемы введем равномерную про-
странственную сетку с шагом Ах = 1/(п - 1), где п — число узлов по на-
правлению х. Аналогично предыдущей задаче следовало бы ввести
сетку и по радиусу г. Чтобы упростить понимание общего подхода, не
будем этого делать и посчитаем, что приемлемую точность решения
можно получить, если и стенку, и жидкость представить в виде одного
ряда последовательных элементов, температура которых характеризу-
ется приближенными величинами TiCT и TiyR (см. рис. 4.19).
Рис. 4.18. Течение жидкости в трубе
99
Рис. 4.19. Упрощенная пространственная сетка
Чтобы получить дифференциальное уравнение, описывающее
процесс теплопроводности в стенке, выделим внутри стенки элемен-
тарно малый объем длиной dx и запишем для него уравнение теплового
баланса: при установившемся режиме количество теплоты, входящее в
выделенный элемент, и теплоты, уходящее из элемента, одинаково:
*
*
Распишем составляющие этого баланса, учитывая, что правая
часть приведенного равенства складывается из теплоты, поступающей
за счет теплопроводности, и теплоты, выделенной внутренними источ-
никами:
*
*
тепл,х
Левая часть тоже представляется двумя слагаемыми:
*
тепл,
*
КОН •
Вспоминая законы Фурье и Ньютона — Рихмана, распишем вели-
чины соответствующих слагаемых:
- количество теплоты, вносимой теплопроводностью в выде-
ленный элемент:
*	=-Х, Sch 
тепл,х /ъст uwl ,
количество теплоты, выделившееся за время <А:
4*
бвн = qidxdx;
- количество теплоты, уходящее из элемента за счет конвектив-
ного теплообмена, полагая при этом, что он « «2 и это позволяет не
учитывать потери тепла с наружной поверхности трубки:
2кон = «2 (Тот - 4 )Pdxdz .
100
Поскольку функция 2тепл непрерывна, то воспользуемся ее раз-
ложением в ряд Тейлора для величины 2*еПл x+dx :
(?тепл, x+dx
©тепл, х
^Qtqwl, X J
------— dx,
dx
где всеми остальными слагаемыми ряда мы пренебрегаем, как величи-
нами более высоких порядков малости.
Подставляя теперь приведенные величины всех слагаемых тепло-
вого баланса в первоначальное уравнение и сокращая на произведение
dx'dx, сначала получим
*
теплых
= а2(Гст-Тж)Г,
*
а затем, подставляя величину утепл х , окончательное уравнение
A Q ^СТ
--- Лсто-----
dx dx
+ ^-а2Гст(Тст-Тж) = 0.
(4.20)
Температурное поле в жидкости в нашей задаче описывается од-
номерным стационарным дифференциальным уравнением энергии,
отражающим в нашем случае баланс между количеством теплоты, за-
траченным на нагрев жидкости в рассматриваемом элементе длиной
dx, и количеством теплоты, полученным в пределах этого элемента от
стенки за счет теплоотдачи:
cpG^ = a2TuimidX(T„-Тх).	(4.21)
dx
При формулировании граничных условий прежде всего запишем
2 уравнения, определяющих теплообмен на входе (х = 0) и выходе из
стенки, отражающих равенство тепловых потоков в зоне теплового
контакта трубки с подводящими электродами:
Температуру жидкости на входе в канал считаем заданной:
101
Перейдем теперь к формированию конечно-разностных аналогов
приведенных дифференциальных уравнений, используя разности
«против хода». Для граничной точки стенки при i = 1 из формулы
(4.20) с учетом значения производной при х = 0 получаем
-X
(h - Л)) — [<71 - “Л (71 - их)].
Дг Я1ст	2
Для внутренних точек стенки (i = 2, 3, п - 1) формула (4.20)
принимает вид:
[^ст(0 (7,+1 - tt) - ^ст(М) (7, - 7,_1)] + qt ~	- Ut) = 0 .
Для граничной точки стенки (при i = п)
[qn &nFCT(tn ^и)] •
Для внутренних узлов жидкости по уравнению (4.21) находим
В приведенных выше формулах ради сокращения записей введено
обозначение:
FrT = лх7пиЛх.
V1
Для первой точки жидкости, при i = 1
Uq— ^вх-
Записанные выше разностные уравнения образуют систему из 2п
нелинейных уравнений, решение которой позволяет определить все
неизвестные температуры и Ц. При этом итерационный процесс в
простейшем случае можно строить в квазистационарной постановке:
на каждом шаге Ах, сначала по значениям температуры _ i на преды-
дущем шаге по функциям fa и f2 рассчитываются значения Хст и qh а
также величины ср и а2, если они тоже меняются вдоль потока, завися
от температуры Ц _ ь Далее любым известным способом (например,
методом Гаусса) решаются приведенные линеаризованные уравнения
и находятся значения и Ц для рассчитываемого шага Ах,. Далее опи-
санный алгоритм повторяется до тех пор, пока i < п.
102
4.8.	ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ
ПРИ ТЕПЛООБМЕНЕ ИЗЛУЧЕНИЕМ
Во многих теплотехнических расче-
тах, особенно при анализе тепловых пото-
ков в топках и нагревательных печах, воз-
никает необходимость рассчитать лучи-
стый теплообмен между телами, разделен-
ными абсолютно прозрачной средой (это
называют диатермической оболочкой), как
показано на рис. 4.20. Введем сначала не-
сколько упрощающих предположений:
-	температурное поле каждой по-
верхности равномерное, стационарное;
Рис. 4.20. Система тел
1, 2,..., ..., п
-	поверхности тел серые, диффузионно отражающие и диффу-
зионно поглощающие;
-	все поверхности плоские, исключающие самооблучение.
Из теории теплопередачи мы знаем, что при решении такой зада-
чи требуется учесть не только собственное излучение поверхности, но
множество других тепловых потоков, действующих на поверхность и
показанных на рис. 4.21. Величина потока собственного излучения
согласно закону Стефана-Больцмана определяется так:
где - степень черноты серого тела; Со - коэффициент излучения аб-
солютно черного тела; 7) - абсолютная температура поверхности тела,
Sj - площадь поверхности.
Падающим потоком называют тепловой поток, падающий
на z-ю поверхность от всех остальных тел
с учетом многократных отражений. Этот
падающий поток частично отражается, а
частично поглощается. По закону Кирхго-
фа коэффициент поглощения равен степе-
ни черноты £(, а коэффициент отражения,
если тело не прозрачно для тепловых лу-
чей, определяется разницей (1 - £(). Тогда
величины поглощенного и отраженного
потоков можно рассчитать так:
Рис. 4.21. Тепловые потоки
при излучении стенки
егл=е;е°ад и er = а - )е°ад.
103
Эффективный поток представляет собой суммарный лучи-
стый поток, уходящий с поверхности на все другие тела. При этом учи-
тываются и собственное излучение, и суммарный отраженный поток:
Величина падающего на z-ю поверхность теплового потока опре-
деляется эффективными потоками от всех остальных поверхностей и
угловыми коэффициентами облучения фу, показывающими, какая доля
тепла, излучаемая телом у, попадает на тело /. Таким образом:
С учетом предыдущих формул можно записать следующее выра-
жение:
Записывая последовательно формулу (4.22) для всех /, мы полу-
чим замкнутую систему уравнений, содержащих п неизвестных эф-
фективных потоков , что дает возможность, применяя известные
процедуры, решить систему и определить все неизвестные.
Как правило, эту систему приводят к канонической форме мат-
ричного уравнения
АХ = В,
где матрица коэффициентов А формируется так:
#м = ~(1 - ег)фу, i,j = 1, 2,..., п, i *j,
а столбец свободных членов состоит из величин
b; Si9 i 1, 2, ..,, п.
Чтобы полнее разобраться с этим вопросом, распишем подробно
все элементы матричного уравнения для приведенной на рис. 4.20 сис-
темы стенок (при п = 5):
ап = 0
#21 = “(I “ ^2)<р21
#31= ~(1 - £з)9з1
«41 = -(1 “ ^4)<р41
#51 = “(1 - ^5)<р51
#12—(1 — £1)<Р12
#22 = 0
#32 = “(I “£з)Ф32
#42=-(1 — ^4)ф42
#52=—(1 - £5)Ф52
#13 = —(1 “£1)<Р13
«23 = “(I “ £2)<р23
«зз = 0
«43 = -(1 “ £4)943
#53 = “(1 - £5)953
104
#14= “0 — £1)914
#24 =_(1 “ £1)ф24
#	34 =-(l — ез)ф34
a44=0
#	54=—(1 - е5)ф55
#	15 = —0 “ £1)915
#	25 = -(1 “ £2)925
#	35 = -(l - £з)935
#	45 = -(1 “ £4)945
#	55 = 0
еэф
е2эф
2зф
е4эф
еэф
b5 = E5CnT52S5
Вспомним о двух свойствах, которыми обладают угловые коэф-
фициенты:
-	свойство замыкаемости состоит в том, что £12 + £21= 1;
-	свойство взаимности для любой пары поверхностей дает
£12*$1 = £21^2-
Используя эти свойства, можно заметно упростить исходную
матрицу и поиск решения.
После решения матричного уравнения и определения величин
Qf$ легко находим результирующие потоки, падающие на каждое тело:
п
вг = <2ГЛ -Qf= ~&° = е£<Р>,&ф - егОД4Хг.;
м
/ — 1,2, ..п, 9y-i— О*
Если какая-нибудь поверхность (или все они) вогнутая, то тогда
возникает самооблучение части этой поверхности (см. рис. 4.22) и для
того, чтобы учесть это, водится коэффициент самооблучения Ф 0.
Это приводит к тому, что в диагонали матрицы один или несколько
коэффициентов становятся не нулевыми.
Конечно же, нам понятно, что каждое тело в системе тел может
иметь собственную температуру 7} при установившемся режиме толь-
ко при наличии в каждом из тел внут-
ренних источников теплоты. В против-
ном случае рано или поздно в резуль-
тате теплообмена в системе наступит
тепловое равновесие при одинаковых
для всех тел температурах.
Поэтому рассмотрим задачу не-
стационарного теплообмена между
системой тел при наличии в них внут-
ренних источников теплоты удельной
производительностью qvi. Для упроще-
ния задачи будем полагать, что тела
настолько высокотеплопроводны, что
температура внутри тела практически
Рис. 4.22. Лучеиспускание
при сложной
конфигурации тел:
1, 2, 3 - лучи, падающие
на другие стенки;
4 - самооблучение
105
одинакова в каждой точке тела. Между телами движется воздух с за-
данными объемными расходами G, и заданной температурой Твоз, в
результате чего между воздухом и телами возникает еще и конвектив-
ный теплообмен, характеризуемый величинами коэффициентов тепло-
отдачи а/.
Решение задачи начинаем с записи уравнения теплового баланса
для /-го тела, приравнивая приходы и расходы теплоты за элементарно
малый промежуток времени dv.
= G™Pic»idTi+at(Ti -4,3)^K0H^ ,
J Ft I	^*4 /	I * I fJl I	I I	HvJ Z I	Z
откуда получаем систему обыкновенных дифференциальных уравне-
ний первого порядка
GiPicpi	= qviVt + Qrsr ~ а, (Г, ~ 4>з)№°Н ,i= 1,2,..., п. (4.23)
ат
В этом уравнении pi,cpi - плотность и удельная теплоемкость ма-
териала /-го тела; /\луч, 7^-кон - поверхности z-го тела, участвующие в
лучистом теплообмене, и величина аналогичной поверхности, участ-
вующей в конвективном теплообмене с воздухом (ниже мы увидим,
что величины этих поверхностей могут быть различными).
Мы уже разбирались, что такие системы решают с помощью ка-
кой-либо явной схемы, например, применяя метод Рунге-Кутта. Одна-
ко в нашем случае перед каждым очередным шагом по времени необ-
ходимо для каждого из тел рассчитать величины QfG3, для чего, как
это показано выше, сначала нужно решить систему алгебраических
уравнений и найти величины эффективных потоков Qf$. Используя
подход, аналогичный квазилинейному, это можно сделать, решив пре-
дыдущую задачу при температурах поверхностей 7}_ь найденных на
предыдущем шаге Ат по времени. При всяком другом подходе сис-
тема уравнений (4.23) практически не разрешима, поскольку каждое из
уравнений содержит по две неизвестных 7} и QfQ3.
Чтобы облегчить работу, заметим, что угловые коэффициенты фу
определяются только геометрическими характеристиками тел и не за-
висят от температуры тел. Это позволяет в самом начале расчетов
сформировать матрицу коэффициентов А, транспонировать ее и найти
обратную матрицу, которая будет использоваться потом при расчетах
теплообмена излучением при увеличении т, поскольку на каждом шаге
будут изменяться только величины свободных членов Ьр
В конечном результате мы получим зависимости 7} = /J(t) для ка-
ждого из тел, участвующих в теплообмене.
106
4.9.	РАСЧЕТ УГЛОВЫХ КОЭФФИЦИЕНТОВ
Угловой коэффициент облучения
между двумя поверхностями площа-
дью 51 и Sz рассчитывается по форму-
ле, которая получена в результате ин-
тегрирования по соответствующим
площадям формулы закона Ламберта:
1 VrcosGi cos07 Tr, Tr,
Ф12 — —	2	dS}dS2, (4.24)
J J r2
где r - расстояние между элементар-
ными площадками на поверхностях 51
и S2, 61 и 02 - углы между лучом, со-
Рис. 4.23. Взаимооблучение
элементарных площадок

единяющим элементарные площадки (рис. 4.23), и нормалями к этим
площадкам.
Таким образом, для вычисления угловых коэффициентов, как это вид-
но из приведенной формулы, нам необходимо вычислить сначала
двойной интеграл по площади Si, а затем — двойной интеграл по пло-
щади Sz, что в целом означает определение четырехкратного интегра-
ла. Такое интегрирование осуществляется достаточно просто только
для тел простой геометрической формы типа плоских прямоугольных
пластин, цилиндров или шаров. Однако в технических устройствах
очень часто встречаются более сложные поверхности, и поэтому при-
ходится прибегать к численному интегрированию.
Рассмотрим отдельные приемы вычисления двойных интегралов.
Прежде всего это метод последовательного интегрирования, суть ко-
торого в том, что двойной интеграл представляют как следующую по-
следовательность обычных интегралов (для прямоугольной области с
координатами сторон ab и cd, см. рис. 4.24):
= [F(y)cfy
(4.24)
с а
где
ь
?(У) = \f(x,y)dx •
а
При расчетах сначала по лю-
бой квадратурной формуле (на-
пример, методом трапеций) про-
Рис. 4.24. Область интегрирования
107
»*
Рис. 4.25. Разбиение
плоской поверхности
на симплекс-элементы
приближенно, заполняя
водят численное интегрирование по гори-
зонтальным прямым и определяют значе-
ния функции F(y) в горизонтальных узлах
разбиения. Затем на основе этих значений
по такой же квадратурной формуле рассчи-
тывают интегралы по вертикальным пря-
мым, а суммируя получаемые при этом
величины, находят интеграл (4.24).
Конечно же, вычисление двойного ин-
теграла усложняется, когда форма поверх-
ности более сложная. Тогда величину двой-
ного интеграла можно определить только
исследуемую площадь прямоугольниками или
другими простыми фигурами. Мы уже знаем, что более полно произволь-
ную поверхность можно учесть с помощью симплексов треугольной фор-
мы (см. рис. 4.25).
В этом случае двойной интеграл
Рис. 4.26. К расчету угловых
коэффициентов при сложных
формах поверхностей тел
I = f(x, y)dxdy
Q
рассчитывают приближенно как
сумму
п________
/=1
где Si - площадь z-го элемента
xt и yt - координаты центра тяже-
сти элемента; п - число симплекс-
элементов. Естественно, что в нача-
ле расчетов нумеруются все верши-
ны элементов, для каждой из них
записываются координаты и по
которым определяются координаты
центра тяжести каждого элемента и
только потом по заданной функции
fix, у) рассчитывают величину оче-
редного слагаемого. Этот метод
называется методом ячеек.
Трудности при расчетах угло-
вых коэффициентов заметно увели-
чиваются, когда исследуемые по-
108
верхности не плоские, а имеют выпуклости или вогнутости. Для при-
мера на рис. 4.26 приведен вид сверху двух криволинейных поверхно-
стей АВ и CD. Если рассматривать элементарную площадку при ка-
ком-либо значении х, например при х = а, то в зависимости от значе-
ния у поведение функции J[a, у) будет меняться и, как это видно из
рисунка, на интервале [с, е] значения функции будут отличаться от
нуля, а на интервале [е, d\ эта функция равна нулю, поскольку на часть
поверхности с дугой ED лучи из точки А не попадают, эта часть по-
верхности не облучается. Доля облучаемой поверхности при увеличе-
нии х будет увеличиваться, и при х = Ъ облучение будет практически
стопроцентным.
Таким образом, при вычислении интеграла вида (4.24) в процессе
численного интегрирования на каждом шаге необходимо проверять
условие облученности и учитывать только облученную поверхность.
То же самое необходимо проверять и при вычислении интеграла по
поверхности СО, при этом из точки D противоположная поверхность
практически не видна и соответствующие слагаемые интеграла равны
нулю. Это значительно усложняет алгоритм расчетов, так что разра-
ботка программы для расчета величины q>i2 — дело весьма не простое,
но общие идеи, позволяющие преодолевать все трудности, мы с вами
теперь понимаем.
Для расчета угловых коэффициентов при сложных и неплоских
поверхностях широко применяется метод Монте-Карло, использую-
щий отдельные особенности последовательности случайных величин.
Чтобы понять суть этого метода, рассмотрим его применение на при-
мере одномерного определенного интеграла
ь
I = f (x)dx.
а
Перепишем его в следующем виде, умножив и разделив на неко-
торую функцию р(х), которая обладает следующим свойством:
ь
p(x)dx -1.
а
Тогда предыдущий интеграл принимает вид
p(x)dx.
' Р(х-)
(4.25)
Познакомимся теперь с некоторыми положениями теории вероят-
ностей. Если мы имеем некоторую случайную величину X, прини-
109
мающую свои случайные значения только из промежутка \а,Ь\ на чи-
словой оси х, и имеющую функцию распределения р(х), то математи-
ческое ожидание величины X (практически это среднее арифметиче-
ское значение для исследуемой последовательности, которое будем
обозначать как Е{Х}) определяется величиной интеграла
ь
Е{Х} = xp(x)dx.
а
Пусть имеется другая случайная величина Л, связанная со слу-
чайной величиной X некоторой функциональной зависимостью
Л = y(JV).
Тогда математическое ожидание последовательности Л будет
ь
Е{Л} = Гт(Др(х)«& •	(4.26)
а
Сравнивая формулы (4.25) и (4.26), видим, что интеграл (4.25)
можно трактовать как математическое ожидание случайной величины
Л, имеющей такую же плотность распределения вероятности р(х), и
поэтому:
Л = Т(Х)
/(*)
А*)
В методе Монте-Карло путем статистического моделирования на
ПК последовательно генерируются случайные числа в интервале [а, Z>]
с заданной функцией распределения (обычно принимается закон рав-
номерного распределения случайной величины). Имея набор случай-
ных чисел xi, Х2,	х„, рассчитывают случайные величины X/ =
= J(xi)/p(xi) и далее находят математическое ожидание величины Л,
заменяя интеграл соответствующей суммой
1 п _	1 п
£{Л} = -£ Х;. = X =	) / p(Xi ) = I*I.
Случайные числа генерируются с помощью специальных стан-
дартных процедур, имеющихся во всех вычислительных пакетах. Вы-
числение кратных интегралов выполняется аналогично, только вместо
случайной величины создается случайный вектор X, два (а в нашем
случае, при расчете угловых коэффициентов, даже 4) набора случай-
ных чисел, соответственно на одной и другой поверхности S] и Sz.
ПО
4.10.	ВОПРОСЫ И ЗАДАНИЯ ДЛЯ САМОПРОВЕРКИ
1.	Как получается алгебраический аналог при численном интег-
рировании задач стационарной теплопроводности?
2.	В чем суть метода сеток, с помощью которого производится
дискретизация температурного поля?
3.	Как получают алгебраический аналог для численного решения
задач при нестационарном одномерном температурном поле?
4.	Явная схема и шаблон численного интегрирования дифферен-
циального уравнения теплопроводности для нестационарного одно-
мерного температурного поля.
5.	Неявная схема и шаблон численного интегрирования диффе-
ренциального уравнения теплопроводности для нестационарного од-
номерного температурного поля.
6.	В чем суть метода прогонки при интегрировании дифференци-
ального уравнения теплопроводности для нестационарного одномер-
ного температурного поля.
7.	Расскажите об особенностях численного интегрирования мно-
гомерных задач нестационарной теплопроводности.
8.	В чем суть метода расщепления решения, используемого для
численного интегрирования многомерных задач нестационарной теп-
лопроводности?
9.	В чем особенности численного интегрирования нелинейных
задач теплопроводности?
10.	Как интегрируют дифференциальное уравнение теплопро-
водности, если тепло физические коэффициенты исследуемого тела
заметно зависят от температуры?
11.	Как интегрируют дифференциальное уравнение теплопро-
водности, если внутренние источники (или стоки) теплоты нелинейно
зависят от температуры?
12.	В чем основная суть метода конечных элементов?
13.	Как составляется функционал /(/(%)) при решении задачи о
плоском нестационарном температурном поле?
14.	Что принимается в качестве симплекс-элемента в методе ко-
нечных элементов?
15.	Расскажите о принятом порядке глобальной и локальной ну-
мерации симплекс-элементов.
16.	Что характеризует собой функция формы симплекс-элемента»?
17.	Как принято описывать функцию формы треугольного сим-
плекс-элемента?
18.	За счет чего переопределенная система уравнений, позво-
ляющая определить коэффициенты функции формы в методе конеч-
ных элементов, упрощается до замкнутой системы уравнений?
111
19.	Какие упрощающие предпосылки вводятся при разработке
алгебраических аналогов для численного интегрирования задач кон-
вективного теплообмена?
20.	Чем метод балансов отличается от метода сеток?
21.	Расскажите о трех подходах к аппроксимации потоков теп-
лоты и массы при реализации метода балансов.
22.	Как формулируются совместные задачи при конвективном
теплообмене?
23.	Что принято называть диатермической оболочкой при расче-
тах теплообмена излучением?
24.	Как формулируется задача теплового расчета диатермиче-
ской оболочки?
25.	Какие упрощающие предпосылки вводятся при анализе теп-
лообмена излучением между замкнутой системой серых тел?
26.	Что называют эффективным тепловым потоком при лучи-
стом теплообмене?
27.	За счет каких свойств коэффициентов облученности удается
значительно упростить систему уравнений, описывающую теплообмен
излучением между телами диатермической оболочки?
28.	Как решается задача о теплообмене между системой тел при
наличии у них внутренних источников теплоты постоянной мощности
и наличии конвективных потоков между телами?
29.	Запишите уравнение теплового баланса для системы тел с
внутренними источниками теплоты, обменивающимися между собой
лучистым и конвективным способами.
30.	Какие математические методы используются для определе-
ния угловых коэффициентов облученности тел сложной формы?
112
5.	ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ
ГИДРО- И АЭРОДИНАМИКИ
5.1.	ИСПОЛЬЗОВАНИЕ МЕТОДА СЕТОК
ДЛЯ РЕШЕНИЯ ПРОСТЕЙШИХ ЗАДАЧ
До сих пор мы сосредоточивали внимание на тепловых задачах.
Однако в практической деятельности и при научных исследованиях в
области теплотехники часто возникает необходимость решать задачи
гидро- или аэромеханики. В первую очередь это задачи о волновых
явлениях при течении потоков (например, о гидравлическом ударе,
резонансных явлениях и др.). В математической физике дифференци-
альное уравнение гиперболического типа второго порядка, описываю-
щее колебательные явления в механических, электрических, гидравли-
ческих и других системах принято называть волновым уравнением,
которое в общем случае для одномерной задачи и без учета отдельных
влияющих факторов (идеализированная модель) записывается в виде:
а2ф(х,т) _ 52ф(х,т)
_	—	_	, 0<х<Д0<т<ткон,
дт2 дх2
(5.1)
где L - длина канала, струны, стержня; ткон - продолжительность про-
цесса (L и ткон - величины, заданные по условиям однозначности). Ре-
шение задачи, которое получается в результате интегрирования этого
дифференциального уравнения, должно удовлетворять как приведен-
ному уравнению, так и задаваемым граничным и начальным условиям.
В разных задачах (электрических, механических или гидравлических)
под функцией ср понимаются разные зависимости отдельных парамет-
ров, меняющихся как в пространстве, так и по времени.
Чтобы понять общий подход при численном интегрировании
уравнения (5.1), рассмотрим сначала задачу о механических колебани-
ях невесомого и абсолютно упругого стержня длиной £, концы которо-
го жестко защемлены (рис. 5.1). Начальное смещение от среднего по-
ложения Ан описывается функцией Ан = fix), а начальная скорость пе-
ремещения струны и>н - функцией wH = g(x). Тогда начально-краевая
задача (задача Коши) будет состоять из дифференциального уравнения
(5.1) и следующих заданных зависимостей:
Ф1(х, 0) =fix), ф2(х, 0) = g(x), 0 < х < L - начальные условия;
Ф1(О, т) = 0, ф2(Д т) = 0, 0 < т < ткон - граничные условия
(слева и справа).
113
Рис. 5.1. Положение стержня в момент времени г
Разобьем непрерывную функцию ф(х, т) пространственно-
временной сеткой с шагами Дх и Лт, узлы которой будем отмечать
индексами / и у, где индекс i определяет положение узла в простран-
стве, а индекс j относится к определенному временному интервалу и
отражает текущий момент времени т (рис. 5.2). Далее, заменяя вторые
частные производные центральными разностями (вывод их приводил-
ся ранее), получаем следующий алгебраический аналог формулы (5.1):
(5.2)
где через U обозначены приближенные значения функции ф в узловых
точках.
Как это видно, формула (5.2) содержит только одну неизвестную
величину Uij+i для последующего момента времени, а все остальные
слагаемые относятся к двум предыдущим моментам времени j и j — 1.
Выразим эту неизвестную:
Рис. 5.2. Пространственно-
временная сетка
В отличие от аналогичных схем для уравнения теплопроводности,
в этой формуле используется 3 временных слоя: (/ - 1)-й,у-й и рассчи-
тываемый (/ + 1)-й. Такой формуле соответствует пятиточечный шаб-
лон, получивший название «крест» (см. рис. 5.3). Это шаблон явной
схемы. Перебирая последовательно i = 1, 2, ..., п - 1, мы рассчитаем
все значения Ui и С72 для всего у + 1-го временного слоя. Эта явная
схема условно-устойчива	(при
Дт/Дх < 1), но требует использова-
ния достаточно малых шагов Лт,
чтобы получать в меру точные ре-
зультаты. О приемах, позволяю-
щих гарантировать необходимую
точность, мы говорили ранее, вы-
деляя, например, метод последова-
тельного уменьшения шага Лт.
114
Определенной особенностью
отличается начало расчетов. На-
чальные условия в разностной фор-
ме (при т = 0, у = 0 и i = 1, 2, п - 1)
записываются так:
t7i(z*Ax, 0-Ат)
0-Ат) = g(z Ах).
Тогда, применяя известный нам
трехточечный шаблон явной схемы,
мы можем рассчитать первый вре-
менной слой (см. рис. 5.4). А уже
после этого можно переходить к
расчету второго временного слоя,
используя пятиточечный шаблон.
На рисунке 5.4 наглядно демон-
стрируется изложенный выше алго-
ритм начала расчетов. Символом о
у-1 О
Рис. 5.3. Пятиточечный
шаблон явной схемы
отражены узлы, для которых значе- Рис. 5.4. Начало расчётов
ния функции ср известны из началь-	по явной схеме
ных и граничных условий, символом * - те узлы, в которых значения
функции рассчитываются по трехточечному шаблону, а символом • -
те узлы, в которых значения функции рассчитываются по пятиточеч-
ному шаблону.
Если, как это делали мы при анализе дифференциального уравне-
ния теплопроводности, правую часть уравнения (5.1) записать для
(/ + 1)-го отрезка времени Ат, то она примет вид
Ц и+1 - 2Ц.у+1 + иМг
Ах2
и тогда получим девятиточечный шаб-
лон неявной трехслойной схемы, со-
держащий 3 неизвестных:	U^+i и
Ui+ij+i. Такой шаблон приведен на
рис. 5.5. Неявная схема является устой-
чивой, но на каждом шаге по времени
приходится решать систему алгебраи-
ческих уравнений (в линейной поста-
новке — например, методом прогонки).
&------сУ-----У''
о------о------о
-м-1 у-1	Z+M-1
Рис. 5.5. Девятиточечный
шаблон неявной схемы
115
5.2.	МЕТОД ХАРАКТЕРИСТИК
Рассмотрим теперь более подробно, как решается линейная зада-
ча Коши для одномерного волнового уравнения, описывающего рас-
пространение плоских звуковых волн в неподвижной (или медленно
движущейся) среде. Математическая формулировка задачи включает
волновое уравнение (уравнение акустики) и соответствующие уравне-
ния условий однозначности, но в более сложном виде, чем это прини-
малось выше.
Заметим, что задачу всегда удобнее решать в относительной, без-
размерной форме. Поэтому перепишем волновое уравнение (5.1), вво-
дя безразмерные параметры X = x/L и безразмерное текущее время
0 = ттск/ткон (L - длина канала, ттек - текущее время, ткон - продолжи-
тельность исследуемого процесса):
^9) = а2Ф(^,е) О<Х<1;О<0<1.	(5.3)
де2 дх2
Для примера возьмем более сложные начальные условия
<рх(Х, 0) = X2 + X и <рт(Т, 0) = cos X.
Обычно задачу с такими достаточно сложными начальными усло-
виями решают методом характеристик, сводя ее сначала к системе
двух уравнений первого порядка. Для этого введем новые искомые
функции, описывающие отклонения скорости w и давления р от их
значений в невозмущенной среде, вызванные распространением зву-
ковых волн вдоль оси х:
„ P(X9)_W)
ае	дх
Тогда преобразованное к безразмерной форме уравнение (5.3)
можно представить так:
д (дд(Х,0)_ д(&p(JT,0p
или, переписывая через вновь введенные функции:
дн(Х,8) др(Х, 6)
ае ат	1'
Полагая, что функции ф, w ир непрерывны и имеют непрерывные
первую и вторую производные, проведем перекрестное дифференци-
рование формулы (5.4):
116
d2w(X,Q) _ d2p(X,Q)
8QdX ~ dXdQ
или, переписав это по-другому
д ( д(?(Х, 0)^|
50 < дх j
д Ла<р(^9е)Л
дх t 50 )
и заменяя частные производные введенными функциями, получаем:
др{X, 0) _ dw(X, 0)
50	” дХ
(5.5)
Таким образом, от исходной задачи мы действительно перешли к
системе двух дифференциальных уравнений (5.4) и (5.5) со следую-
щими условиями однозначности:
w(X, 0) = cos X, p(Xfi) =	<p(T,0) =	(X2 + X) = 2X +1.
оХ	дХ
Систему уравнений (5.4) и (5.5) можно записать в векторно-
матричной форме, вводя вектор-столбец U и матрицы коэффициентов
Ли В размерностью 2x2:
А^ + В^
50 дХ
(5.6)
с начальными условиями
U(X,0) = JC0Sj¥ |.
[2T + 1J
В уравнении (5.6) вектор U и матрицы коэффициентов следую-
щие:
Запишем теперь формулу для полного дифференциала U, как
функции двух переменных X и 0:
dU =
откуда выразим второе слагаемое
117
dU
50
<70 = dU-
dU
dX
(5.7)
Уравнение (5.6) умножим почленно на dQ:
A^-dQ + B^-de = O
50 дх
и подставим в первое слагаемое выражение по формуле (5.7):
в™
дХ
<70 = 0
или, вынося частную производную за скобки, получим другое матрич-
ное уравнение:
(в<70 - AdX)
dU
дХ
= -AdU
Будем считать, что искомые функции гладкие и имеют непрерыв-
ные первую и вторую производные, а значит, система (5.4) — (5.5) име-
ет решение U(X, 0), которое можно представить в координатах X - 0
некоторой точкой у (см. рис. 5.6). Определитель матрицы коэффициен-
тов (в<70 - AdX) может быть отличным от нуля и тогда система имеет
действительные корни. Если же определитель этой матрицы равен ну-
лю, то система имеет множество решений, удовлетворяющих ей, т.е.
функция U(X, 0) находится в этом случае неоднозначно. И тогда кри-
вую, отражающую все множество возможных решений у, называют
характеристикой матричного уравнения.
Приравнивая нулю определитель названной матрицы и, посколь-
ку <70 ф 0, разделив на <70 оба слагаемые, получаем
или В - АХ = 0
Рис. 5.6. Характеристика
матричного уравнения
где X = dX/d§. Далее распишем этот оп-
ределитель
Рассчитываем определитель (-Х)-
(-Х) - (-1 )•(-!) = X2 - 1 и решаем урав-
118
2
нение X - 1 = 0: Xi = 1, Х2 = —1. Вместо величин X запишем их значения
и
Разнесем переменные и проинтегрируем полученные уравнения:
JV— 0 + Ci и X—0 + С2
Этим решениям соответствуют
два семейства характеристик вида СЛ и
С~, проходящих под углом 45° к оси X
и параллельных друг другу (рис. 5.7).
Чтобы получить так называемые
характеристические соотношения, ум-
ножим правую и левую части формулы
(5.5) на произведение dXd&.
Рис. 5.7. Семейство характери-
стик для системы двух
линейных уравнений
др{Хdw{X
-----uXdxj —--d/Cdxj,
50	dX
откуда, выделяя полные дифференциалы, получаем
dp-dX= dwdft или
Поскольку
= X = ±1, то предыдущая формула дает
dp = dw и dp = -dw.
Интегрированием приведенных выражений получаем так назы
ваемые инварианты Римана:
p — w — const] и
которые используют при разработ-
ке численной схемы для расчетов
значений искомых функций.
Чтобы лучше понять суть под-
хода, возьмем на оси X (рис. 5.8)
отрезок [0; 1] и зададим на нем все-
го 3 точки (для получения более
подробных и точных результатов
таких точек задают гораздо боль-
ше): X = 0, X = 0,5 и
X = 1,0. Проведем через эти точки
линии характеристик и пронумеру-
р + w = const2,
Рис. 5.8. Характеристическая сетка
119
ем индексами i и j узловые точки характеристической сетки. Точки 1 о,о,
21,о и 32д соответствуют нулевому временному слою и величины ско-
рости w и давления р в них известны из начальных условий:
wi = и>(0, 0) = cosO = 1; p\=p(Q, 0) = 2Х+ 1=20+1 = 1;
w2 = w(0,5, 0) = cosO,5 = 0,8776; p2=p(0,5, 0) = 2X+ 1 = 2-0,5 + 1=2;
w3 = w(l, 0) = cos 1 = 0,5403; p2=p(l, 0) = 2X+ 1 = 2-1 + 1= 3.
Из рисунка 5.6 мы видим, что при таком подходе мы получаем
трехточечный угловой шаблон, который содержит две нижних точки с
известными к началу расчета величинами w и р, одну точку (в верши-
не), содержащую неизвестные w и р (в нашем примере это шаблоны с
узлами 1, 2 и 4, 2, 3 и 5, 4, 5 и 6.
Определим теперь неизвестные параметры для узла 4, аппрокси-
мируя значение производной вдоль характеристики из семейства С4"
(при X = 1) разностью «против хода (параметры конца минус парамет-
ры начала)»:
Аналогично вдоль другой характеристики (при X = —1):
Приравнивая полученные дроби числам 1 и -1, соответственно,
приходим к системе с двумя алгебраическими уравнениями, содержа-
щими неизвестные Х4 и 04:
А4-А] = l(04-0i);
Х2 - Х4 = -1 (02 - Од)-
Находим эти неизвестные, учитывая, что Xi = 0, Х2= 0,5, 01 = 0,02= 0.
Первое уравнение дает Х4 = 04, а второе -0,5 - Х4 = 04. С учетом
предыдущего равенства получаем 2 Х4 = 0,5, откудаХ4 = 04 = 0,25.
Полученное ранее условие Римана вдоль характеристики вида С4"
для лежащих на ней узлов 1 и 4 дает: р4- w4 = pi~ wi, а для другой
характеристики С~- р4 + w4 = р2 + w2. В итоге мы получили еще одну
систему, содержащую две неизвестных р4 и w4. Решаем ее, учитывая,
что р\ = 1, р2 = 2, wi = 1, w2 = 0,8776, и находим: из первого уравнения
р4 - w4 = 1 - 1 = 0, откуда р4 = w4. Второе уравнение дает: р4 + w4 =
= 2 + 0,8778 = 2,8778. С учетом предыдущего равенства получаем
1р4= 2,8778, откудар4 = и>4 = 1,4388.
120
Совершенно аналогично можно рассчитать теперь угловой шаб-
лон с точками 2, 3 и 5, а затем и последний шаблон с точками 4, 5 и 6,
поскольку все искомые параметры для точек 4 и 5 будут уже опреде-
лены.
Таким образом, мы подробно рассмотрели и теоретические осно-
вы, и алгоритм решения задач данного класса. Заканчивая знакомство
с решением такого типа задач, заметим лишь, что ради упрощения мы
рассматривали волновые процессы в идеальной среде, что позволило
не учитывать ее сжимаемости, плотности и других параметров. При
решении конкретных задач, например о гидравлическом ударе в тру-
бах, формулы (5.4) и (5.5) записываются более сложно (их называют
формулами Н.Е. Жуковского):
ди>(х,т) др(х,т)
др(х,т) dw(x,x)
где ро - плотность среды при атмосферном давлении; Е - приведенный
модуль упругости трубопровода:
1__I ^вн _|_]_
Е' 5 Е”
tjxq Ег - модуль упругости среды; <7ВН - внутренний диаметр трубы;
6 - толщина ее стенки; Е” - модуль упругости материала трубопрово-
да. Конечно же, введение названных констант не меняет общего под-
хода, но несколько усложняет расчеты, поскольку величины р0 и Е
будут входить соответствующими множителями в расчетные формулы
для каждого элемента треугольной сетки.
5.3.	ВОПРОСЫ И ЗАДАНИЯ ДЛЯ САМОПРОВЕРКИ
1.	Как получают алгебраический аналог для волнового уравне-
ния, описывающего явления гидро- и аэромеханики?
2.	Почему для расчетов неустановившихся процессов гидроди-
намики используется пятиточечный шаблон явной схемы?
3.	В чем суть метода характеристик, используемого для решения
задач гидродинамики?
4.	Что называют инвариантами Римана?
5.	Какие физические свойства жидкости и материала стенки тру-
бы учитываются при расчетах гидравлического удара в трубах?
121
ЗАКЛЮЧЕНИЕ
Подытоживая изложенный в настоящем учебном пособии учеб-
ный материал, можно с уверенностью сказать, что знакомство с основ-
ными подходами математического моделирования и алгоритмизации
инженерных задач в области теплоэнергетики и теплотехники воору-
жает будущего специалиста в названной предметной области новыми,
весьма важными и перспективными компетенциями, дающими воз-
можность на самом современном уровне получать требуемые практи-
кой решения.
Теплоэнергетика и теплотехника как предметная область прило-
жения полученных знаний и компетенций отличается очень большим
разнообразием используемых процессов, устройств и агрегатов. Мате-
матические модели здесь всегда разрабатываются или на основе ре-
зультатов экспериментальных исследований, или на основе известных
законов сохранения (массы, энергии, количества движения и др.). По-
этому эти модели могут представляться как в виде очень простых ли-
нейных или нелинейных уравнений, так и в виде сложных, линейных
или нелинейных дифференциальных уравнений первого или второго
порядка в частных производных и систем таких уравнений.
Умение разрабатывать математические модели объекта исследо-
вания и алгоритмы их реализации на ПК позволяет выполнять имита-
ционное моделирование, обеспечивающее выявление свойств объекта
при изменении его отдельных исходных параметров, открывая воз-
можность выбрать такие значения варьируемых величин, при которых
результат будет наилучшим (в заданном смысле) по сравнению со все-
ми другими допустимыми вариантами. Другими словами - наличие
математической модели является первым необходимым условием для
того, чтобы поставить и решить задачу оптимизации объекта. Именно
поэтому в федеральном государственном образовательном стандарте и
рабочей программе, разработанной в соответствии с этим документом,
изучению методов математического моделирования, приобретению
практических навыков реализации полученных моделей на ПК уделе-
но особое внимание.
122
СПИСОК ЛИТЕРАТУРЫ
1.	Самарский, А.А. Математическое моделирование. Идеи, мето-
ды, примеры / А.А. Самарский, А.Г. Михайлов. - М., 2001. - 320 с.
2.	Киреев, В.И. Численные методы в примерах и задачах /
В.И. Киреев, А.В. Пантелеев. - М. : Высш, шк., 2004. - 480 с.
3.	Применение ЭВМ для решения задач теплообмена /
Г.Н. Дульнев и др. - М.: Высш, шк., 1990. - 207 с.
4.	Сегерлинд, Л. Применение метода конечных элементов /
Л. Сегерлинд. - М.: Мир, 1979. - 392 с.
5.	Бояринов, А.Н. Методы оптимизации в химической техноло-
гии / А.Н. Бояринов, В.В. Кафаров. - М. : Химия, 1975.-575 с.
6.	Дульнев, Г.Н. Тепло- и массообмен в радиоэлектронной аппа-
ратуре / Г.Н. Дульнев. - М. : Высш, шк., 1984. - 247 с.
7.	Теоретические основы теплотехники. Теплотехнический экс-
перимент : справочник / под общ. ред. А.В. Клименко и В.М. Зорина. -
3-е изд - М. : Изд-во МЭИ, 2001. - 564 с.
8.	Ляшков, В.И. Теоретические основы теплотехники / В.И. Ляш-
ков. - М. : Высш, шк., 2008. - 318 с.
9.	Зубарев, В.Н. Практикум по технической термодинамике /
В.Н. Зубарев, А.А. Александров, В.С. Охотин. - М. : Высш, шк., 1986. -
314 с.
10.	Ляшков, В.И. Компьютерные расчеты в термодинамике /
В.И. Ляшков. - Тамбов : ТГТУ, 1997. - 163 с.
11.	Арутюнов, В.А. Математическое моделирование тепловой ра-
боты промышленных печей / В.А. Арутюнов, В.В. Бухмиров,
С.С. Крупенников. - М. : Металлургия, 1990. - 239 с. (Даются про-
граммы расчетов разных тепловых задач).
12.	Тарко, Л.М. Волновые процессы в трубопроводах гидромеха-
низмов / Л.М. Тарко. - М.: Машгиз, 1963. - 184 с.
12.	ELCUT. Руководство пользователя. - СПб. : Производствен-
ный кооператив ТОР, 2003. - 231 с.
123
ГЛОССАРИЙ
№
п/п
Новые понятия
Содержание
9
10
Алгоритм
решения задачи
Аналитические
методы решения
задач («решения
в квадратурах»)
Декомпозиция
объекта
исследований
Квази стацио-
нарная модель
объекта
исследований
Математическая
модель объекта
Моделирование
Неявная схема
численного
интегрирования
Объект с рас-
пределенными
параметрами
Объект с сосре-
доточенными
параметрами
Расчетный
шаблон
Последовательность действий, в результате
которых поставленная задача будет обяза-
тельно решена_________________________
Методы, позволяющие получить точное
решение в виде математических формул
Условное разбиение объекта на отдельные
части, определенным образом взаимодей-
ствующие между собой, обеспечивая нор-
мальное функционирование объекта_____
Абстрактная упрощенная модель объекта, в
которой считается, что на отдельных про-
межутках времени Ат параметры не меня-
ются по времени, а изменения их происхо-
дят в конце каждого такого промежутка
Зависимость, которая в математической
форме однозначно описывает связь между
всеми входными и выходными параметра-
ми объекта
Изучение свойств и поведения объекта ис-
следований без непосредственного исполь-
зования этого объекта, без внедрения в не-
го____________________________________
Реализация алгебраического аналога, обес-
печивающего получение для каждой точки
расчетной сетки замкнутой системы из
трех (или более) уравнений____________
Объект, физические свойства которого за-
висят от пространственных координат
Объект, физические свойства которого не
зависят от пространственных координат и
одинаковы в любой его точке___________
Графическое отображение алгебраического
аналога дифференциального уравнения,
записанного для определения величины
расчетного параметра в одной или группе
соседних точек расчетной сетки
124
Продолжение табл.
№
п/п
11
12
13
14
Новые понятия
Сходимость
численного
метода решения
У стойчивость
численного
решения
Фазовое гипер-
пространство
Хорошо
о бу с лов ленная
задача
15 Черный ящик
16
Численные
методы решения
задач
17 Явная схема
численного
интегрирования
18 Ячеистая модель
объекта
исследований
Содержание
Свойство метода, обеспечивающее при
уменьшении шага итераций стремление
результатов расчета к точному решению
Один из признаков хорошо обусловленной
задачи, выражающийся в том, что в диапа-
зоне возможного изменения входных па-
раметров численное решение не имеет раз-
рывов или не удаляется от точного реше-
ния___________________________________
Абстрактное графическое представление п-
мерного пространства (п > 3), в котором
изображается зависимость выходного па-
раметра у от ряда входных параметров хь
Х2, ..., Хп_________________________________
Задача, поставленная так, что при неболь-
ших изменениях исходных данных (напри-
мер, в пределах допустимой погрешности
измерений) результаты решения изменяют-
ся незначительно
Абстрактная модель исследуемого объекта,
характеризующаяся наличием входных и
выходных параметров, зависимости между
которыми неизвестны___________________
Методы приближенного решения задач
путем реализации алгоритмов, соответст-
вующих математическим моделям и даю-
щих только частные решения____________
Реализация алгебраического аналога, обес-
печивающего получение для каждой точки
расчетной сетки одного уравнения с одной
неизвестной___________________________
Упрощенная абстрактная модель объекта,
условно разбитого в пространстве на не-
сколько частей (ячеек), внутри которых
параметры одинаковы. Однако в разных
таких ячейках параметры различаются по
величине. Изменения параметров происхо-
дят на границах каждой ячейки
125
ПРИЛОЖЕНИЕ
ВАРИАНТЫ ЗАДАНИЙ ДЛЯ КУРСОВОЙ РАБОТЫ
Курсовая работа выполняется магистрантом самостоятельно в те-
чение всего учебного семестра и служит для практического освоения и
закрепления методик и приемов математического моделирования с
решением задач на ПК. Работа выдается на второй неделе второго се-
местра. Выполненную работу сдают и защищают в начале зачетной
недели.
При итоговом контроле в конце семестра наряду с проверкой и
устным собеседованием непосредственно по работе, итоговая оценка
определяется ещё и результатами устного собеседования по отдель-
ным теоретическим вопросам (см. вопросы для самопроверки). Работа
оценивается по общепринятой пятибалльной системе.
Курсовая работа содержит 5 типовых практических задач, пере-
чень которых формируется для каждого магистранта индивидуально
путем указания номера задачи из каждой группы. Каждая группа задач
ориентирована на отработку одного из приемов моделирования и со-
держит не менее десяти задач.
Группа 1. Решение нелинейных уравнений вида F(x) = 0 или
X =fix).
Содержание: определить действительные корни нелинейного
уравнения на заданном интервале значений аргумента х.
№ задачи	Задача	Интервал поиска решения: а <х < Ъ		Вид функции F(x) или fix)
1	Решить уравнение F(x) = 0	а = 0	b = 5	F(x) = Зх3 - 2х2 + х - 10
2	Решить уравнение F(x) = 0	а = 0	b = 10	F(x) = Зх3 - 2ехр(х2/10) + + 4х - 50
3	Решить уравнение F(x) = 0	а = 0	b =2	F(x) = Зх3 + 4ехр(-х/10) + + 4х-20
4	Решить уравнение F(x) = 0	а = -2	Ь = 2	F(x) = Зх2+ 4sin(x/10) + + 2х- 2
5	Решить уравнение X =Лх)	а = -2	Ь=2	fix) = -Зх2 - 4cos(x/5) + 6
126
Продолжение табл.
№ задачи	Задача	Интервал поиска решения: а <х < b		Вид функции F(x) или fix)
6	Решить уравнение F(x) = 0	а = 0	b =4	F(x) = sin(x)cos (х) - 0,2
7	Решить уравнение F(x) = 0	а = ОД	b = 1,1	F(x) = sinx + 2exp(3x-J°’2) - 70
8	Решить уравнение F(x) = 0	а = 0	b =2	F(x) = 0,4exp(-5x°’8) - 0,01
9	Решить уравнение F(x) = 0	а = 0	b = 2	F(x) = 3x4 - 4exp(-x/15) - - 2 Ox +10 = 0
10	Решить уравнение x=flx)	а = 0,2	Ъ = 1	fix) = 5 sinx + + 2exp(2x“°’4) - 60
Группа 2. Решить систему линейных уравнений.
№
задачи
Система уравнений
№
задачи
Система уравнений
5xi + 6x2+ 7хз = 38
2xi + 4х2+6х3 = 28
3xi - 8х2 - 4х3 = -31
5xi + 4х2 + 7х3 = 44
2xi + 6x2 + 6x3 = -26
-3xi + Зх2 - 4х3 =15
-2xi + 4х2 + 1х3 = 6
Зх1 + Зх2+2хз= 16
4xi + 2х2 + Зх3= 18
5xi + 6x2 + 7х3 = 34
2xi + 4х2 + 6х3 = 20
-3xi - 8х2 - 4х3 = -29
5xi + 6х2 + 7х3 = 54
2xi + 4х2 + 6х3 = 36
-3xi - 8х2- 4х3 = -45
2xi + 4х2 + 1х3 = 20
—3x1 Зх2 2х3 = 23
4xi + 2х2 - Зх3 = 26
4
2xi + 4х2 + 1х3= 13
Зх1 + Зх2+2хз= 15
—4х] + 2х2 + Зх3 = 21
5xi + 4х2 + Зх3= 13
3xi + 4х2 + 5х3= 15
Xi + 2х2 + Зх3= 21
Х1 + Х2 + Хз + Х4= 10
2xi + 4х2 + Х3+Х4 =17
-3xi + Зх2 + 2х3+х4 =19
4xi + 2X2 + Зхз -Х4 =21
10
2xi + 2х2 + 2х3 + 2х4 = 20
2xi + Зх2 + 4х3 + 5х4 = 40
3xi + Зх2 + Зх3 + Зх4 = 30
4xi + Зх2 + 2х3 + х4 = 20
127
Группа 3. Решить систему нелинейных уравнений.
№
задачи
Система уравнений
№
задачи
Система уравнений
2х3 - х| -1 = О
cos(xj + 0,5) - х2 = 2
sin хэ - 2х^ -1
sin(x1 +1) - х2 = 1
2хг + cos(x2) = 2
2xj + tg(x1x2) = 0
6х2 - б)2 + In х1 = О
tg(xj х2 + 0,2) = х2
+ 2х2
х1 - 1х2 +1 = 0
- х2 + 2х2 -1 = 0
10
exp(xj -1) + х2 - 2 = О
х1 + х2 - 3 = О
х2+х2-9 = 0
х^ + х2 —1 = 0
Xi3 - Хо = О
х2 + х2 - 9 = О
4
9
Группа 4. Составить математическую модель и алгоритм реше-
ния задачи, решить задачу численным способом. Записанные курсивом
значения параметров для каждого исполнителя задаются преподавате-
лем индивидуально.
Проанализировать недостатки модели, разработанной с учетом
оговоренных упрощающих предпосылок, сформулировать предложе-
ния для уточнения модели за счет более точного и полного учета со-
ставляющих теплового или материального балансов.
Задача 1
Провести численный анализ процесса разогрева электропечи для
непрерывного отжига сердечников трансформаторов. Схематически
конструкция печи представлена на рис. 1. Через окна в передней и зад-
ней стенке с постоянной скоростью на роликах перемещается цепной
конвейер, на который закрепляются обрабатываемые заготовки (кон-
вейер на схеме не показан). Разогрев печи и заготовок производится
десятью карбид-кремниевыми нагревателями, которые питаются от
сети трехфазного переменного тока напряжением 380 В. Масса нагре-
вателя /ин = 2 кг, диаметр = 0,04 м, длина /н = 1 м, удельная теплоем-
кость Сн= 0,9 кДж/(кг-К), степень черноты ен= 0,7. Стенки печи выпол-
128
Рис. 1. Устройство электропечи
йены из огнеупорного теплоизоляционного материала толщиной
5 = 0,06 м с плотностью рст= 2500 кг/м3, теплоемкостью сст= 1 кДж/(кг-К)
и коэффициентом теплопроводности Хст = 1,5 Вт/(м-К). Размеры печи
(длина, ширина и высота) указаны на рисунке, размеры входного и
выходного окон проставлены там же. В процессе разогрева конвейер с
заготовками не работает, его включают только тогда, когда температу-
ра в печи, измеряемая специальной термопарой, достигнет 810 °C, что
составляет 90% от номинальной рабочей температуры.
Рассчитать термограмму разогрева и определить время т включе-
ния конвейера, принимая средний коэффициент теплоотдачи от на-
ружных поверхностей печи в окружающий воздух а = 40 ВтДм^К),
Зависимость электрического удельного сопротивления Яуд нагревате-
лей принимать линейной: 7?уд = 2,92 + 0,001/, Ом/м. В течение разогре-
ва в печи находятся 5 заготовок массой т3 = 3 кг с удельной теплоем-
костью Сн = 0,8 кДж/(кг-К). Потерями теплоты на разогрев воздуха в
объеме печи пренебрегать. Считать, что теплообмен между электрона-
гревателями и внутренними стенками печи и потери теплоты через
окна в окружающее пространство происходят только тепловым излу-
чением. Температуру наружной поверхности печи в течение процесса
разогрева принимать средней и постоянной, равной tCT2= 35 °C. Темпе-
ратуру на внутренней поверхности стенки /CTi определять приближенно
по формуле
^Т1 = /Ст2+ 1314п(т+ 1),
где т - текущее время.
Задача 2
На Тамбовскую ТЭЦ пришла цистерна топочного мазута для пи-
кового котла. Для слива мазута к цистерне присоединен вентиль с тру-
бой, опущенной в резервную емкость (см. рис. 2). Проходное сечение
полностью открытого вентиля равно FBCHT = 0,01 м2. Расстояние от
днища цистерны до конца сливной трубы по вертикали h = 1,2 м. Цис-
терну считать цилиндром с диаметром d = 3 м и длиной I = 5 м.
129
I
Рис. 2
Определить, с какой скоростью будет опускаться уровень мазута
в цистерне, каким он будет через 170 с после начала слива и через
сколько времени весь мазут выльется из цистерны.
Упрощающие указания: возможное изменение температуры ок-
ружающей среды в течение процесса не учитывать; вязкость мазута
учитывается только коэффициентом истечения ц = 0,75', плотность
мазута р = 875 кг/м3.
Задача 3
Составить математическую модель про-
цесса сжатия дизтоплива в плунжерном насосе
высокого давления и на ее основе определить
максимальную цикловую подачу топлива при
давлении открытия выпускного клапана рот =
= 24 МПа. В рабочую полость насоса топливо
подается топливоподкачивающим насосом под
давлением /?нач =0,15 МПа. Размеры основных
деталей: диаметр втулки d = 15 мм, ход плун-
жера S = 25 мм (см. рис. 3), зазор между плун-
жером и втулкой 6 = 0,01 мм, длина зазора L =
= 30 мм. Плунжер приводится в движение с
помощью кулачка и толкателя. Кулачек спроек-
тирован так, что в пределах угла
а = 140° закон перемещения плунжера сину-
соидальный: X = Ssin(180a/140). Коэффициент
изотермической сжимаемости дизтоплива pz =
= 0,7-КГ4 см2/кг = 0,7 -1(Г3 1/МПа. Утечки через
зазор рассчитывать по формуле Бернулли:
Рис. 3
(7ут Рут^щели	Р ’
130
где Gyr - объемный расход утечек; ц - коэффициент истечения (при-
нимать ц = 0,87)', р — плотность дизтоплива (принимать р = 875 кг/м3);
Др - перепад давлений на концах щели (текущего давления в полости
над плунжером р и давления /?нач).
Построить также график зависимости изменения давления р по
углу поворота а и по времени т, учитывая, что после открытия выпу-
скного клапана давление выталкивания не меняется и остается
24 МПа. При расчетах принимать число оборотов четырехтактного
двигателя п= 750 об/мин.
Задача 4
Разработать математическую модель, описывающую динамику
процесса разогрева воды в электрокотле с пассивными электродами,
устройство которого приведено на рис. 4, а электрическая схема - на
рис. 5. С помощью этой модели определить продолжительность выхо-
да котла на установившийся режим при фазовом напряжении U= 380 В
и массовом расходе нагреваемой воды М= 20 кг/мин. Какую темпера-
туру будет иметь вода на выходе из котла через 150 с после включения
котла?
Рис. 4. Схема электрокотла и его основные конструктивные размеры:
1 - пассивный электрод; 2 - фазовый (активный) электрод; 3 — корпус котла;
4 — компенсационная емкость; 5 — разборный электроизолятор
131
Рис. 5. Распределение электрических сопротивлений:
ПЭ - пассивный электрод; АЭ — активный электрод
Указания к анализу и моделированию:
1.	Абсолютное давление воды в котле определяется суммой баро-
метрического давления В и избыточного давления, определяемого вы-
сотой столба воды Н.
2.	Расход воды через котел таков, что образование пара внутри
пассивного электрода не происходит (размер h = 0).
3.	Электрическое сопротивление проточной зоны определяется по
формуле:
D 1 . 3a2(J2-a2)3
Rn = Pt —i-ln—,
2л/п J2(JK5 6-a6)
где pz - удельное электрическое сопротивление воды, принимать
pz= 15 Омм.
4.	Для непроточной зоны (внутри пассивного электрода) электри-
ческое сопротивление определяется по-другому:
5. Рабочий фазный ток в котле
U
*n+2^
6. Мощность, выделяемая по зонам
7У=73/2ЯП и N„ = 24bi1Rw.
132
7. Коэффициент теплопередачи, определяющий величину тепло-
потерь в окружающую среду, принимать к = 28 Вт/(м2-К), коэффици-
ент теплопередачи через стенку пассивного электрода к наружной воде
кп = 185 Вт/(м2-К).
8. При расчетах принимать следующие размеры и характеристи-
ки: dK = 0,23 м, dn = 0,082 м, d3 = 0,057 м, /э = 0,6 м, Н = 3 м, hK = 0,8 м,
радиус а = 0,07 м, массы стальных деталей типэ = 0,85 кг, пц3 = 1,4 кг,
корпуса тК = 18 кг.
Задача 5
В инструкциях по эксплуатации автомобиля содержится указание
о том, что после запуска двигателя необходимо, работая на холостом
ходу, дождаться его прогрева хотя бы до температуры 60 °C и только
после этого начинать движение.
Составьте математическую модель, описывающую процесс про-
грева двигателя, полагая, что среднюю температуру его хорошо харак-
теризует температура охлаждающей воды, циркулирующей в контуре:
рубашка цилиндра - радиатор (см. рис. 6). Масса воды тв = 8 кг, масса
двигателя 7ида = 90 кг, масса радиатора тр = 10 кг. Радиатор выполнен
из латунных трубок, теплоемкость латуни ся = 0,64 кДж/(кг-К). Сред-
нюю теплоемкость деталей двигателя принимать сд = 0,47 кДж/(кг-К).
На холостом ходу удельный расход топлива составляет Вхх = 0,00022 кг/с.
Низшая теплота сгорания автомобильного бензина 0РН = 44,1 МДж/кг.
При расчетах принимать:
1.	Теплообмен между двигателем и окружающей средой осуще-
ствляется вынужденной конвекцией в окружающий воздух с коэффи-
циентом теплоотдачи анар = 250 Вт/^-К).. Поверхность теплоотдачи
принимать FHap = 1,8 м2.
2.	Коэффициент теплопередачи в радиаторе принимать завися-
щим от температуры t охлаждающей воды по формуле £рал = 35 + 0,05/.
Поверхность теплообмена Fpaa = 0,8 м2.
3.	Из теплового баланса известно, что на холостом ходу с дымо-
выми газами теряется 52% выделившейся в двигателе теплоты.
Рис. 6. Схема двигателя с радиатором
133
4.	Механический КПД двигателя т)м= 0,96.
С помощью математической модели определить продолжитель-
ность прогрева двигателя на холостом ходу от температуры гнач = 15 °C
до температуры t = 95 °C. Определить также, через сколько секунд
после запуска двигателя можно начинать движение, не нарушив заво-
дской рекомендации?
Задача 6
Составьте математическую модель, описывающую процесс вы-
равнивания давлений в двух баллонах с азотом, соединенных трубкой
диаметром d = 0,006 м и снабженной краном, практически мгновенно
позволяющим открыть трубку для прохода газа. Объемы баллонов
одинаковы (V= 250 литров). Давление в первом баллоне р\ = 15 МПа,
во втором - р2= 2 МПа. Температура газов в начале и в течение про-
цесса практически не меняется и равна t = 20 °C.
Определить продолжительность процесса до полного выравнива-
ния давлений, а также давления в баллонах через 3 с после открытия
крана на трубке. Построить график изменения давлений в баллонах от
времени.
Задача 7
Стальную цилиндрическую заготовку диаметром d= 45 мм и дли-
ной L = 60 мм необходимо прогреть для последующей закалки до
700 °C. Разогрев ведется в муфельной печи с микропроцессорным
управлением, позволяющим задавать различные режимы разогрева пу-
тем изменения напряжения, подаваемого на вход электронагревателей.
Определить продолжительность прогрева заготовки до заданной
температуры, а также температуру через 10 мин после включения
печи. При этом отрабатываемая программа приведена на рис. 7, общее
сопротивление электронагревателей печи 7?н = 10 Ом, теплопотери в
окружающую среду происходят в результате свободной конвекции, а
удельный тепловой поток теплопотерь составляет 8% от тепловыделе-
ния нагревателей в любой момент времени. Определить также, какова
будет температура на поверхности электропечи в конце разогрева, ес-
ли температура в помещении во время процесса была 23 °C, а коэффи-
циент теплоотдачи от поверхности
печи в окружающую среду анар =
= 15 Вт/(м2-К).
Задача 8
В водогрейной газовой котель-
ной для хранения резервного дизель-
ного топлива установлены два резер-
вуара диаметром D = 2 м и высотой
Рис. 7. Программа работы
электронагревателей
134
Н = 2 м с рабочим объемом каждого по V = 6 м3. Резервуары соедине-
ны между собой трубопроводом с задвижкой на dy = 30 мм. Первый
резервуар был заполнен дизтопливом плотностью р = 875 кг/м3 на весь
рабочий объем, а второй оставался пустым, когда на 5 мин открыли
разделяющую их задвижку.
Принимая коэффициент истечения при полном открытии задвиж-
ки ц = 0,88, определить уровни дизтоплива в каждом из резервуаров в
момент закрытия задвижки. Каковы были эти уровни через 2 мин по-
сле открытия задвижки? Какие объемы дизтоплива при этом находи-
лись в резервуарах? Изобразите графически изменения уровней во
времени за все 5 мин процесса перетекания дизтоплива из бака в бак.
Задача 9
На промышленном предприятии, имеющем собственную котель-
ную с паровыми котлами, в абонентском вводе одного из отдельно
стоящих производственных зданий установлен емкостный водоподог-
реватель (см. рис. 8). При запуске системы отопления здания в змеевик
этого подогревателя подают насыщенный водяной пар при давлении
рп= 0,17 МПа. Водяной объем подогревателя Кв= 4 м3, начальная тем-
пература воды в нем /нач = 15 °C. Масса стальных деталей подогревате-
ля типод = 850 кг, удельная теплоемкость стали сст = 0,46 кДж/(кг-К).
Средний коэффициент теплоотдачи от наружной поверхности подог-
ревателя окружающему воздуху анар = 8 Вт/(м2-К), поверхность тепло-
отдачи FHap = 6 м2, температуру поверхности принимать равной темпе-
ратуре воды в подогревателе, температура воздуха в помещении або-
нентского ввода ?возд= 15 °C.
На время запуска патрубок 7 для выхода нагретой воды перекры-
вается вентилем и циркуляция воды в системе отопления отсутствует.
Вентиль открывают, когда температура воды в подогревателе достиг-
нет 75 °C.
8	7	6’
Рис. 8. Емкостный водоподогреватель:
7 - корпус; 2 — вход воды из системы отопления; 3 — змеевик;
4 - коллекторы змеевика; 5 — вход пара и выход конденсата; 6 — манометр;
7 — патрубок для выхода нагретой воды в систему отопления;
8 - бобышка для предохранительного клапана
135
Разработать математическую модель и алгоритм ее реализации
для расчета кривой разгона описанного аппарата. Определить время, за
которое температура воды в подогревателе поднимется до 50 °C, если
во время разогрева расход пара М остается постоянным и равным
0,4 кг/с. Определите время, через которое следует открыть вентиль и
подать подогретую воду в систему отопления.
Задача 10
В современных котельных агрегатах, особенно в прямоточных
котлах и котлах со сверхкритическим давлением, для регулирования
температуры первичного пара широко используются впрыскивающие
пароохладители. С их помощью в поток перегретого пара впрыскива-
ется питательная вода, которая там испаряется, забирая для этого часть
теплосодержания пара. В результате температура пара несколько
уменьшается, а его расход увеличивается на величину расхода впрысну-
той воды. Одна из конструкций такого устройства приведена на рис. 9.
Рассчитать характеристики переходного процесса пароохладите-
ля, если в начальный момент времени пароохладитель работал при
следующих установившихся параметрах на входе: давление, темпера-
тура и расход пара были соответственно /?вх = 12,75 МПа (130 кг/см2),
£вх= 550 °C, Сгп.вх= 62 кг/с, расход питательной воды был GB. нач= 0,6 кг/с,
а ее температура 4 = 250 °C.
Переход на новый режим осуществлялся системой автоматиче-
ского регулирования так, что расход питательной воды увеличивался
по линейному закону в течение 5 минут от GB нач до величины
GB кон = 0,9 кг/с при постоянстве всех остальных входных параметров.
При расчетах принимать давление питательной воды одинаковым
с давлением пара, изменением теплопотерь и гидравлических сопро-
тивлений за процесс пренебрегать.
Рис. 9. Впрыскивающий пароохладитель:
1 - корпус; 2 - диффузор; 3 - бобышки для датчиков САР;
4 — отверстия для распиливания воды; 5 — сопло; 6 — компенсатор;
7 - водоподводящая трубка; 8 — защитный патрубок водоподводящей трубки
136
ОГЛАВЛЕНИЕ
ВЕДЕНИЕ................................................. 3
1.	ОБЩИЕ ПОНЯТИЯ И ПОДХОДЫ...............................5
1.1.	Математическое моделирование как средство решения
практических задач................................5
1.2.	Приемы разработки математических моделей.........11
1.3	Содержание и состав теоретической математической
модели............................................17
1.4.	Характеристики объектов моделирования........18
1.5.	Вопросы и задания для самопроверки...........20
2.	МЕТОДЫ АНАЛИЗА ПРОСТЫХ МАТЕМАТИЧЕСКИХ
МОДЕЛЕЙ.............................................22
2.1.	Решение нелинейных уравнений вида F(x) = 0...22
2.2.	Решение системы линейных уравнений методом Гаусса .. 30
2.3.	Итерационные методы решения систем линейных
уравнений.........................................33
2.4.	Методы решения систем нелинейных уравнений...36
2.5.	Особенности решения систем нелинейных уравнений..39
2.6.	Решение систем нелинейных уравнений методом
«пристрелки по цели»..............................44
2.7.	Вопросы и задания для самопроверки...........52
3.	МАТЕМАТИЧЕСКИЕ МОДЕЛИ НЕУ СТАНОВИВШИХСЯ
ПРОЦЕССОВ И МЕТОДЫ ИХ РЕАЛИЗАЦИИ....................54
3.1.	Задача об истечении газа из ограниченного объема.54
3.2.	Исследование режимов работы смесительного
теплообменника....................................60
3.3.	Методы Симпсона,	Рунге-Кутта, Адамса.........64
3.4.	Вопросы и задания	для самопроверки...........69
4.	ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ
ТЕПЛООБМЕНА.........................................70
4.1.	Численное решение задач стационарной
теплопроводности..................................70
4.2.	Численное решение нестационарных задач
теплопроводности..................................74
4.3.	Решение многомерных задач теплопроводности...79
4.4.	Решение нелинейных задач теплопроводности....81
4.5.	Метод конечных элементов.....................83
4.6.	Конечно-разностные методы решения задач
конвективного теплообмена.........................91
4.7.	Решение совместных задач при конвективном
теплообмене.......................................98
137
4.8.	Численные методы решения задач при теплообмене
излучением....................................103
4.9.	Р асчет угловых коэффициентов............107
4.10.	В опросы и задания для самопроверки.....111
5.	ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧ ГИДРО-
И АЭРО ДИНАМИКИ.................................113
5.1.	Использование метода сеток для решения простейших
задач......................................113
5.2.	Метод характеристик......................116
5.3.	Вопросы и задания для самопроверки.......121
ЗАКЛЮЧЕНИЕ.......................................122
СПИСОК ЛИТЕРАТУРЫ................................123
ГЛОССАРИЙ........................................124
ПРИЛОЖЕНИЕ...................................... 126
Учебное электронное издание
ЛЯШКОВ Василий Игнатьевич
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И
АЛГОРИТМИЗАЦИЯ ЗАДАЧ
ТЕПЛОЭНЕРГЕТИКИ
Учебное пособие
Редактор Т.М. Г линкина
Инженер по компьютерному макетированию М.С. Анурьева
Подписано в печать 28.05.2012.
Формат 60 х 84 /16. 8,14 уел. печ. л. Заказ № 251
Издательско-полиграфический центр ФГБОУ ВПО «ТГТУ»
392000, г. Тамбов, ул. Советская, д. 106, к. 14