Расчет одноконтурной автоматической системы регулирования температуры печи котельного агрегата
Федеральное агентство по рыболовству
Федеральное государственное образовательное учреждение
высшего профессионального образования
«Астраханский государственный технический университет»
«УТВЕРЖДАЮ»
заведующий кафедрой АТП
Есауленко В.Н.
«_____» _________ 2011 г.
Курсовой проект
«Расчет одноконтурной автоматической системы регулирования температуры печи котельного агрегата»
Пояснительная записка КП 220301.072676.2011
проект выполнил
ст.гр. ДИА-41 Югов С.Г.
руководитель к.т.н., доц. Кокуев А. Г.
Астрахань 2011 г.
Астраханский Государственный Технический Университет
Кафедра «Автоматизация технологических процессов»
Дисциплина «Теория автоматического управления»
Специальность «Автоматизация технологических процессов и производств»
Курс 4 Группа ИА-41 Семестр 7
ЗАДАНИЕ
на курсовой проект (работу) студента Югова С.Г.
Тема проекта (работы): расчет одноконтурной автоматической системы регулирования.
Срок сдачи студентом законченного проекта:
Исходные данные к проекту: кривая разгона объекта регулирования по каналу регулирующего воздействия.
Содержание расчетно-пояснительной записки (перечень подлежащих разработке вопросов): 1. Передаточная функция объекта управления; 2. Расчет и построение переходной функции; 3. Получение аппроксимирующих передаточных функций для объекта управления; 4. Расчет и построение частотных характеристик объекта управления; 5. Расчет оптимальных настроек регулятора методом расширенных частотных характеристик; 6. Расчет и построение КЧХ замкнутой системы; 7. Построение переходного процесса по возмущению в системе регулирования приближенными методами; 8. Оценка качества переходного процесса по возмущению.
Перечень графического материала: частотные характеристики объекта управления; переходные процессы по каналам управления и возмущения; схема моделирования; кривая настроек регулятора; КЧХ замкнутой системы регулирования.
Дата выдачи задания:
Руководитель к.т.н., доц. Кокуев А.Г.
Выполнил Югов С.Г.
Календарный план
Наименование этапов курсового проекта (работы) |
Срок выполнения этапов проекта (работы) |
Примечание |
|
1 |
Получение задания |
6.09.2010 |
|
2 |
Введение |
12.09.2010 |
|
3 |
Краткие теоретические сведения |
16.09.2010 |
|
4 |
Расчет одноконтурной АСР |
8.10.2010 |
|
5 |
Заключение |
22.10.2010 |
|
6 |
Список литературы |
25.11.2010 |
|
7 |
Графическая часть |
26.01.2011 |
|
Студент
Руководитель
Введение
Современные автоматизированные системы управления техническими процессами требуют значительного количества и разнообразия средств измерений, обеспечивающих выработку сигналов измерительной информации в форме, удобной для дистанционной передачи, сбора, дальнейшего преобразования, обработки и передачи.
В настоящее время существует большое число различных по своему назначению систем автоматического регулирования. Одни из них поддерживают заданную температуру, давление, расход жидкости или газов в объектах регулирования, другие изменяют эти параметры по различным законам.
Автоматические системы применяют и для управления скоростью вращения гидравлических и паровых турбин, дизелей, регулирования напряжения на электростанциях. Их используют также для регулирования мощности в ядерных энергетических реакторах, удержания электронного пучка в линейных ускорителях, регулирование тока в физических установках.
Автоматизация - это применение комплекса средств, позволяющих осуществлять производственные процессы без непосредственного участия человека, но под его контролем. Автоматизация производственных процессов приводит к увеличению выпуска, снижению себестоимости и улучшению качества продукции, уменьшает численность обслуживающего персонала, повышает надежность и долговечность машин, дает экономию материалов, улучшает условия труда и техники безопасности.
Автоматизация агрегатов включает в себя автоматическое регулирование, дистанционное управление, технологическую защиту, теплотехнический контроль, технологические блокировки и сигнализацию.
Автоматическое регулирование обеспечивает ход непрерывно протекающих процессов в парогенераторе (питание водой, горение, перегрев пара и др.)
Системы управления современными химико-технологическими процессами характеризуется большим числом регулируемых параметров. Так, что контуров регулирования сложных химико-технологических комплексов может достигать нескольких сотен.
В качестве регуляторов в подавляющем большинстве систем управления в нефтехимии, нефтепереработке, энергетике, металлургии и др. отраслях промышленности России и зарубежных стран в основном используются так называемые типовые промышленные регуляторы П-, ПИ- и ПИД-законы регулирования. Широкий диапазон изменения настроечных параметров типовых регуляторов позволяет использовать их для управления процессами с различной инерционностью, обеспечивает их взаимозаменяемость, удобство в эксплуатации и, в конечном счете, надежность систем управления.
Несмотря на развитие теории оптимального управления, разработку серийных регуляторов с переменной структурой, типовые законы регулирования по-прежнему составляют значительное большинство в системах управления промышленными процессами. В автоматизированных системах управления технологическими процессами (АСУТП), реализованных на основе мини- и микро-ЭВМ, доля типовых законов (алгоритмов) все еще велика, в особенности на нижних уровнях управления. В числе алгоритмов управления, реализуемых микропроцессорными контроллерами типа «Ремиконт», ПИ- и ПИД-законы являются одними из основных. В современных микропроцессорных системах Микро-Z и МОД-300 управление на нижнем уровне в значительной степени также осуществляется по типовым законам, реализованным цифровым способом.
Указанными обстоятельствами объясняется внимание, уделяемое проблеме расчета настроечных параметров типовых регуляторов учебными программами вузов соответствующих специальностей. Значительная часть существующих в настоящее время промышленных АСР являются одноконтурными. Это объясняется целым рядом причин. Но последнюю роль играют в этом отсутствие надежных технических средств и сложность алгоритмов расчета, требующих большого объема информации. В то же время широко применяемые на практике каскадные АСР, системы с дифференцированием промежуточной переменной и другие в силу специфики их динамических свойств приводятся к одноконтурным и рассчитываются в 2 этапа как одноконтурные.
Это же относится и к многомерным системам, расчет настроек которых во многих случаях сводится к многократному расчету приведенных одноконтурных. Таким образом, методы и алгоритмы расчета настроек одноконтурных АСР, рассматриваемые в пособии, являются основой расчета систем более сложной структуры. К настоящему времени разработано достаточно большое число приближенных и точных методов и методик расчета настроек.
Отдельную группу составляют приближенные методы (экспресс методы), позволяющие по минимуму информации о динамике процесса определить параметры настройки регулятора.
Точные методы, использующие полную информацию о динамике процесса, требуют определенных вычислительных затрат и, как правило, обеспечивают минимум некоторых критериев оптимальности.
Понятие о котельной установке
Паровым котлом называется комплекс агрегатов, предназначенных для получения водяного пара. Этот комплекс состоит из ряда теплообменных устройств, связанных между собой и служащих для передачи тепла от продуктов сгорания топлива к воде и пару. Исходным носителем энергии, наличие которого необходимо для образования пар из воды, служит топливо.
Основными элементами рабочего процесса, осуществляемого в котельной установке, являются:
1) процесс горения топлива,
2) процесс теплообмена между продуктами сгорания или самим горящим топливом с водой,
3) процесс парообразования, состоящий из нагрева воды, ее испарения и нагрева полученного пара.
Во время работы в котлоагрегатах образуются два взаимодействующих друг с другом потока: поток рабочего тела и поток образующегося в топке теплоносителя. В результате этого взаимодействия на выходе объекта получается пар заданного давления и температуры.
Одной из основных задач, возникающей при эксплуатации котельного агрегата, является обеспечение равенства между производимой и потребляемой энергией. В свою очередь процессы парообразования и передачи энергии в котлоагрегате однозначно связаны с количеством вещества в потоках рабочего тела и теплоносителя.
Горение топлива является сплошным физико-химическим процессом. Химическая сторона горения представляет собой процесс окисления его горючих элементов кислородом, проходящий при определенной температуре и сопровождающийся выделением тепла. Интенсивность горения, а так же экономичность и устойчивость процесса горения топлива зависят от способа подвода и распределения воздуха между частицами топлива. Условно принято процесс сжигания топлива делить на три стадии: зажигание, горение и дожигание. Эти стадии в основном протекают последовательно во времени, частично накладываются одна на другую.
Расчет процесса горения обычно сводится к определению количества воздуха в м3, необходимого для сгорания единицы массы или объема топлива количества и состава теплового баланса и определению температуры горения.
Значение теплоотдачи заключается в теплопередаче тепловой энергии, выделяющейся при сжигании топлива, воде, из которой необходимо получить пар, или пару, если необходимо повысить его температуру выше температуры насыщения. Процесс теплообмена в котле идет через водогазонепроницаемые теплопроводные стенки, называющиеся поверхностью нагрева. Поверхности нагрева выполняются в виде труб. Внутри труб происходит непрерывная циркуляция воды, а снаружи они омываются горячими топочными газами или воспринимают тепловую энергию лучеиспусканием. Таким образом, в котлоагрегате имеют место все виды теплопередачи: теплопроводность, конвекция и лучеиспускание. Соответственно поверхность нагрева подразделяется на конвективные и радиационные. Количество тепла, передаваемое через единицу площади нагрева в единицу времени носит название теплового напряжения поверхности нагрева. Величина напряжения ограничена, во-первых, свойствами материала поверхности нагрева, во-вторых, максимально возможной интенсивностью теплопередачи от горячего теплоносителя к поверхности, от поверхности нагрева к холодному теплоносителю.
Интенсивность коэффициента теплопередачи тем выше, чем выше разности температур теплоносителей, скорость их перемещения относительно поверхности нагрева и чем выше чистота поверхности.
Образование пара в котлоагрегатах протекает с определенной последовательностью. Уже в экранных трубах начинается образование пара. Этот процесс протекает при больших температуре и давлении. Явление испарения заключается в том, что отдельные молекулы жидкости, находящиеся у ее поверхности и обладающие высокими скоростями, а следовательно, и большей по сравнению с другими молекулами кинетической энергией, преодолевая силовые воздействия соседних молекул, создающее поверхностное натяжение, вылетают в окружающее пространство. С увеличением температуры интенсивность испарения возрастает. Процесс обратный парообразованию называют конденсацией. Жидкость, образующуюся при конденсации, называют конденсатом. Она используется для охлаждения поверхностей металла в пароперегревателях.
Пар, образуемый в котлоагрегате, подразделяется на насыщенный и перегретый. Насыщенный пар в свою очередь делится на сухой и влажный. Так как на теплоэлектростанциях требуется перегретый пар, то для его перегрева устанавливается пароперегреватель, в которых для перегрева пара используется тепло, полученное в результате сгорания топлива и отходящих газов. Полученный перегретый пар идет на технологические нужды.
Автоматическое регулирование котельных установок
Система автоматического регулирования котельных установок обеспечивает изменение производительности установки при сохранении заданных параметров (давления и температуры пара) и максимального КПД установки. Кроме того, повышает безопасность, надежность и экономичность работы котла, сокращает количество обслуживающего персонала и облегчает условия его труда. Автоматическое регулирование котла включает регулирование подачи воды, температуры перегретого пара и процесса горения. При регулировании питания котла обеспечивается соответствие между расходами воды, подаваемой в котел, и вырабатываемого пара, что характеризуется постоянством уровня воды в барабане.
Регулирование питания котлов малой производительности обычно осуществляется одноимпульсными регуляторами, управляемыми датчиками изменения уровня воды в барабане. В котлах средней и большой паропроизводительности с малым водяным объемом применяются двухимпульсные регуляторы питания котла по уровню воды и расходу пара, а также трехимпульсные. Управляющие питанием котла по уровню воды, расходу пара и перепаду давлений на регулирующем клапане.
Регулирование температуры пара осуществляется регулятором, управляемым датчиками изменения температуры перегретого пара на выходе из пароперегревателя, изменения температуры пара в промежуточном коллекторе пароперегревателя и изменения температуры газов в газоходе пароперегревателя, а иногда еще датчиком изменения давления пара.
Регулирование процесса горения в топке котла (в соответствии с расходом пара) осуществляется регуляторами подачи топлива II, воздуха III и регулятором тяги IV (см. рис 3.22). Регуляторы подачи топлива II и воздуха III управляются датчиком изменения давления перегретого пара I, а регулятор тяги IV – датчиком изменения разрежения в топке 7 котла.
Рисунок 2. Схема автоматического регулирования котельной установки
1 — бункер угля;
2—шаровая мельница;
3 — сепаратор;
4 — циклон;
5 — бункер пыли:
6 -мельничный вентилятор;
7 — топка котла;
8 — барабан котла;
9—пароперегреватель;
10 — пароохладитель;
11 — экономайзер;
12-воздухоподогреватель;
13 — вентилятор;
14 — дымосос; I — датчик измерения давления перегретого пара:
II — регулятор топлива;
III — регулятор воздуха;
IV — регулятор тяги;
V — регулятор загрузки мельницы;
VI — регулятор температуры мельницы.
В котельных установках, работающих на пылевидном топливе, осуществляется также регулирование работы пылеприготовительной системы регулятором V загрузки мельниц, обеспечивающим постоянство загрузки шаровых барабанных мельниц и регулятором VI температуры пылевоздушной смеси за мельницей. Для предупреждения персонала о недопустимости отклонения важнейших параметров котельной установки от заданных служат звуковые и световые сигнализаторы предельных уровней воды в барабане, предельных температур перегретого пара и низшего давления питательной воды.
Для обеспечения правильной последовательности операций при пуске и остановке механизмов применяется блокировка. Так, при аварийном отключении дымососов отключаются дутьевые вентиляторы и прекращается подача топлива в топку.
Работа котельных установок должна быть надежной, экономичной и безопасной для обслуживающего персонала. Для выполнения этих требований котельные установки эксплуатируются в соответствии с правилами устройства и безопасной эксплуатации паровых котлов и рабочими инструкциями, составленными на основе правил Госгортехнадзора с учетом местных условий и особенностей оборудования.
Котел должен быть оборудован необходимым количеством контрольно-измерительных приборов, автоматической системой регулирования важнейших параметров котла, защитными устройствами, блокировкой и сигнализацией. Режимы работы котла должны соответствовать режимной карте, в которой указываются рекомендуемые технологические и экономические показатели его работы: параметры пара и питательной воды, содержание СO>2> в газах, температура и разрежение по газовому тракту, коэффициент избытка воздуха и т.п. Большинство современных котельных установок полностью автоматизированы.
При нарушении нормальной работы котла в следствие неисправностей, которые могут привести к аварии, он должен быть немедленно остановлен. Капитальный ремонт котлов производится через каждые два-три года. Котел периодически подвергается техническому освидетельствованию по трем видам:
наружный осмотр (не реже одного раза в год);
внутренний осмотр (не реже одного раза в четыре года);
гидравлическое испытание (не реже одного раза в восемь лет).
Регулирование температуры перегретого пара
Регулирование t>ne> осуществляется пароохладителями .Обычно используются пароохладители впрыскивающего типа (рис. 1.2). Они представляют собой участок паропровода, внутрь которого при помощи форсунок впрыскивается охлаждающая среда. Снижение температуры происходит за счет испарения капель среды. На участке, где происходит испарение, внутри трубопровода есть защитная металлическая рубашка – нужна, чтобы капли не попадали на
стенки труб, т.к. они работают под давлением и на них будут образовываться трещины. Проще всего было бы впрыскивать питательную воду в необходимом количестве (так делают на западе).
В России экономят на системах очистки воды, поэтому питательная вода грязная и соль будет оседать на стенках. В России впрыскивают конденсат, который образуется из-за охлаждения части насыщенного пара в водяном экономайзере (конденсационная установка, рис. 1.3). На котлах низкого давления используются поверхностные пароохладители.
Таким образом, пароохладитель позволяет изменить в сторону уменьшения.
Рис.1.2. Пароохладитель впрыскивающего типа.
Рис. 1.3. Конденсационная установка.
Обычно, на котлах имеются системы автоматического регулирования, в том числе и t>ne>. В данном случае регулирование осуществляется изменением расхода конденсата на впрыск.
Регулирование температуры перегрева нужно:
1. Современные котлы работают в широком диапазоне нагрузок. Производительность по пару меняется от 100% до 50% (и даже до 30%). При снижении производительности в два раза, котел сжигает примерно в два раза меньше топлива, расход продуктов сгорания также снижается в двое. Т.к. конструкция неизменна, то скорость газов снизится в два раза, а коэффициент теплопередачи k убывает. Тепловосприятие ПП ( H – площадь поверхности нагрева), т.к. конструкция не меняется? а температурный напор уменьшается при снижении нагрузки. В связи с этим и будет снижаться. При проектировании поверхность нагрева ПП выбирают такой, чтобы поддерживалась на минимальной нагрузке, но тогда на номинальной нагрузке 100% эта поверхность будет избыточной и ее лишнее тепловосприятие компенсируется при помощи пароохладителей.
2. Котел – сложная динамическая система, которая непрерывно находится под влиянием внешних воздействий:
а)
Плановые и аварийные остановки мельниц
b)
Изменение режимов работы тягодутьевых
машин
c) Переход котла в процессе работы
с одного топлива на другое (даже качество
одного топлива непрерывно меняется).
3.
Знания и расчетные методики несовершенны
(возможно даже есть ошибки). Пароохладитель
позволяет компенсировать.
4. Пароохладитель устанавливается перед теми ступенями ПП, которые работают в наиболее тяжелых температурных условиях. Снижение температуры перед ними повышает их надежность.
Динамические свойства объекта регулирования по каналу регулирующего воздействия определены кривой разгона.
Регулируемая величина — температура (С) нагревательной печи. Кривая разгона получена при скачкообразном изменении подачи топлива (газа) - x>вх> (% х.р.о.). Максимальный расход газа 40000 нм3/час.
Безразмерная кривая разгона и вариант задания представлен в таблицах 1 и 2.
Таблица 1
0 |
0 |
1 |
3 |
7 |
11 |
15 |
18 |
21 |
23 |
26 |
27 |
28 |
29 |
29,5 |
30 |
30 |
|
0 |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
8 |
9 |
10 |
11 |
12 |
13 |
14 |
15 |
16 |
Таблица 2
Вариант |
x>вх> имп, % х.р.о. |
a>t>, сек |
а>>, С |
Регулятор |
Метод расчета |
2 |
25 |
25 |
1,5 |
ПИ |
М=1,62 |
1. Построить кривую разгона в размерном виде и определить по ней динамические параметры объекта – k, T, τ, ε, ρ. Записать передаточную функцию, аппроксимирующие объект сочетанием инерционного звена первого порядка с запаздыванием — .
Для того чтобы получить кривую разгона в размерном виде умножим ее безразмерные значения на соответствующие масштабные коэффициенты, результат запишем в таблице 3.
Таблица 3
0 |
0 |
1,5 |
4,5 |
10,5 |
16,5 |
22,5 |
27 |
31,5 |
34,5 |
39 |
40,5 |
42 |
43,5 |
44,25 |
45 |
45 |
|
0 |
25 |
50 |
75 |
100 |
125 |
150 |
175 |
200 |
225 |
250 |
275 |
300 |
325 |
350 |
375 |
400 |
Определение динамических параметров объекта:
Т=170 с — постоянная времени объекта
— коэффициент усиления объекта;
— степень самовыравнивания;
— скорость разгона;
=70 c — полное запаздывание объекта;
>т>=35 c — транспортное запаздывание объекта.
2. Определить по кривой разгона методом интегральных площадей (Симою) передаточную функцию регулируемого объекта -
Передаточную функцию объекта представим в виде:
где τ – транспортное запаздывание; k – коэффициент усиления; T>3>, T>2>, T>1> – постоянные времени.
Постоянные времени T>3>, T>2>, T>1> определим по следующим формулам:
, где dt =30 – шаг дискретизации;
, где α = dt/T>1>;
.
– значение переходной характеристики в i-й момент времени
Перепишем передаточную функцию объекта, подставляя найденные коэффициенты:
3. Решить полученное из передаточной функции, дифференциальное уравнение при заданном значении входа x>вх>, построить расчетную кривую разгона и сопоставить ее с заданной. Представить динамическую модель объекта соединением типовых динамических звеньев, смоделировать объект на базе имитационного моделирования (Simulink), получить на модели кривую разгона и сравнить ее с рассчитанной
Из передаточной функции объекта получаем дифференциальное уравнение:
Решим это дифференциальное уравнение без учета запаздывания.
Подставим в это уравнение известные значения:
Начальные условия:
Характеристическое уравнение имеет вид:
Корни характеристического уравнения:
Общее решение имеет вид:
Частное решение:
Для нахождения коэффициентов С>1>, С>2>, С>3>> >воспользуемся начальными условиями:
Решим эту систему уравнений матричным методом
Таким образом, уравнение кривой разгона, учитывая запаздывание, имеет вид:
, С
Объект управления можно представить в следующем виде:
– усилительное звено;
– звено чистого запаздывания;
– апериодическое звено первого порядка;
– колебательное звено;
где ;
;
, С
4. Рассчитать и построить частотные характеристики объекта регулирования (КЧХ, АЧХ и ФЧХ).
Имеем объект управления:
Для того, чтобы построить частотные характеристики объекта, перейдем от преобразований Лапласа к преобразованиям Фурье. Для этого сделаем замену
Разложим знаменатель на множители:
Амплитудно-частотная характеристика (АЧХ) имеет вид:
A(), C
, c-1
Фазо-частотная характеристика (ФЧХ):
(), рад
, c-1
Комплексная частотная характеристика (КЧХ):
5. Для заданного регулятора заданным методом рассчитать оптимальные значения параметров настройки, обеспечивающих заданный запас устойчивости системы регулирования и качество переходных процессов. Расчет параметров ведется по расширенным частотным характеристикам.
Передаточная функция ПИ-регулятора имеет вид
Рассчитаем оптимальные настроечные параметры методом расширенных частотных характеристик для М=1,62
М=(1+)/2м
м=0,35
Для этого в выражение передаточной функции подставим :
Разложим знаменатель на множители:
Расширенная амплитудно-частотная характеристика (АЧХ):
Расширенная фазо-частотная характеристика (ФЧХ):
Воспользуемся критерием Найквиста
Отсюда:
Выбор параметров настройки: точек С>1> и С>0>
Выберем точку правее максимума:
C>1>=0.2;
C>0>=0.0012;
6. Составить структурную схему системы регулирования (при найденных оптимальных настройках регулятора). Получить передаточную функцию замкнутой системы относительно внешнего возмущающего воздействия. В качестве передаточной функции объекта относительно возмущающего воздействия взять передаточную функцию , определенную в п.1.
— передаточная функция замкнутой системы по возмущающему воздействию.
7. Рассчитать и построить КЧХ замкнутой системы относительно возмущающего воздействия.
Сделаем замену
Рассмотрим 2-ю скобку знаменателя:
Ее реальная часть:
Ее мнимая часть:
Запишем это комплексное число в показательной форме в виде:
,
где и .
Аналогично запишем первую скобку знаменателя:
, где и .
Числитель можно записать в виде:
По полученной аналитически комплексной частотной характеристике строим годограф.
8. Методом Акульшина построить переходный процесс в системе регулирования при единичном скачкообразном возмущающем воздействии. Смоделировать систему в программе Matlab, получить переходной процесс при заданном возмущающем воздействии и сравнить его с расчетным.
Простейшим методом приближенного построения переходных процессов по частотным характеристикам является метод Акульшина. Исследуемый переходный процесс представляется в виде элемента некоторой периодической функции, период которой равен удвоенному значению длительности процесса. Тогда входная величина в случае скачкообразного возмущающего воздействия может быть представлена следующим рядом Фурье:
Так как рассматриваемая система линейна, то колебания выходной (регулируемой) величины могут быть представлены только единственным образом в виде ряда Фурье:
;
;
Отсюда:
9. Оценить качество полученного переходного процесса
Время регулирования t>р>=2000 с;
Перерегулирование ;
Время достижения первого максимума: t>max>=500 с;
Частота переходного процесса с-1
5. Коэфф. затухания за 1 период:
10. Вывод
В ходе выполнения курсовой работы исследовали заданную систему:
- построили кривую разгона, в результате получили передаточную функцию методом Симою;
- построили частотные характеристики объекта регулирования;
- рассчитали оптимальные настройки ПИ-регулятора, обеспечивающие заданный запас устойчивости и качество переходных процессов;
- построили переходной процесс в системе регулирования при скачкообразном возмущающем воздействии;
- определили качество полученного процесса, из которого:
, по методу Акульшина
Так как , значит метод Акульшина дает достаточно точный результат
11. Список литературы:
1. «Методы Классической и Современной Теории Автоматического Управления», т.1,2 под.ред. Н.Д.Егупова, М., изд. МГТУ им.Баумана, 2004.
2. «Теория Автоматического Управления», т.1 под. ред. А.А. Воронова, М., «Высшая Школа», 1986
3. Б.Р.Андриевский, А.Л.Фрадков «Элементы Математического Моделирования», С.-П. «Наука», 2001