Адаптивная математическая модель тепловых процессов трехжильного силового кабеля

УДК 621.315

АДАПТИВНАЯ МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ТЕПЛОВЫХ ПРОЦЕССОВ ТРЕХЖИЛЬНОГО СИЛОВОГО КАБЕЛЯ

Д. И. ЗАЛИЗНЫЙ, О. Г. ШИРОКОВ

Учреждение образования «Гомельский государственный технический университет имени П. О. Сухого»,

Республика Беларусь

Введение

В большинстве научных публикаций расчет тепловых процессов в силовых кабелях выполняют сеточными методами с использованием программных пакетов ANSYS, MATLAB и др. Как правило, рассматривают стационарные тепловые процессы, например, в [1]—[3], но имеются публикации и по нестационарным процессам, такие, как [4], [5]. Основной недостаток подобных методов - необходимость в подробных сведениях обо всех геометрических параметрах кабельной линии, которые могут варьироваться в процессе ее прокладки и, соответственно, влиять на результаты расчета, а также зависимость от постоянно изменяющихся условий окружающей среды. Это не критично на этапе проектирования, однако при контроле температур в работающей кабельной линии возможны высокие погрешности. Кроме того, сеточные методы в совокупности с расчетом нестационарных тепловых процессов требуют значительных вычислительных возможностей, что ведет к усложнению и повышению стоимости аппаратной и программной составляющих систем защиты и диагностирования силовых кабелей.

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

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

Математическая модель четвертого порядка для тепловых процессов

трехжильного силового кабеля

На рис. 1 показан чертеж трехжильного кабеля с секторными жилами и обозначениями геометрических параметров в соответствии с ГОСТ Р МЭК 60287-2-1 [7]. Далее все расчеты будем вести на примере этого типа кабелей.

На рис. 1 обозначено: r1 - радиус окружности, описанной вокруг жил; t - толщина изоляции между жилами; da - наружный диаметр поясной изоляции; Da - наружный диаметр брони; ts - толщина защитного покрытия.

Рассмотрим трехжильный кабель как систему из следующих однородных тепловых тел: каждой из жил вместе с их изоляцией; защитной оболочки; внутренней охлаждающей среды (температура которой зависит от температуры поверхности кабеля [6]); внешней охлаждающей среды.

Рис. 1. Поперечный разрез трехжильного кабеля с секторными жилами

В ГОСТ Р МЭК 60287-2-1 тепловыми сопротивлениями металлических элементов пренебрегают, так как значения этих сопротивлений намного ниже, чем значения тепловых сопротивлений изоляции. Например, для алюминия удельная теплопро-

Вт Вт

водность составляет 209 --------, а для кабельной бумаги - 0,167 -------. Поэтому и

м • °С

м • °С

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

Теплоемкостями металлических элементов пренебрегать нельзя, поскольку

Дж

удельная теплоемкость для алюминия составляет 922

Дж

кг • °С

-, а для кабельной

бумаги - 2300

кг • "С

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

Для кабеля, проложенного в воздухе, теплоемкостью внутренней охлаждающей

среды будем пренебрегать. Хотя удельная теплоемкость воздуха равна 1005 —Дж ,

кг • °С

его плотность на несколько порядков ниже, чем плотность элементов кабеля.

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

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

Рис. 2. Структурная схема тепловых процессов для трехжильного кабеля

Систему дифференциальных уравнений, описывающую тепловые процессы в соответствии со структурной схемой на рис. 2, можно представить в виде

С ^1,1 + ^1,1 _ ©1,2 + ^1,1 _ ©1,3 + ^1,1 “02 _ р.

у-!- "Т" ""Т -*1 ч

л к

к

к

ЛВ, 2 0,, — 0,2 0,2 — 013 0,2 — 02

С1—12 + 1,1 „ 1,2 + 1,2 ^ 1,3 + 2 _р2.

л к, к,

к2

С Л01,3 +01,1 — 01,3 , 01,2 “ 01,3 , 01,3 “ 02

1 л к

+

к1

+

С ^02 , 01,1 — 02 , 01,2 “ 02 , 01,3 “ 02 , 02 —03 _ 0-

к2

22 й©3 02 — 03 03 — 04

к2

к3

С ^ +-^1 + ~^4 _ 0,

л к3 к4

где 011, 012, 013, 02, 03, 04 - температуры соответствующих однородных тел, °С;

° С

R1 - тепловое сопротивление между жилами, --------; R2 - тепловое сопротивление

Вт

° С

между любой из жил и защитной оболочкой, ; R3 - тепловое сопротивление ме-

° С

жду защитной оболочкой и внутренней окружающей средой, -------------------; R4 - тепловое

Вт

° С

сопротивление между внутренней и внешней окружающими средами, -----------; С, С2,

Вт

С3 - соответственно, теплоемкости жил и их изоляции, защитной оболочки и внут-

Вт • с п п п

ренней окружающей среды, --------; р, Р2, Р3 - потери активной мощности в соот-

С

ветствующих жилах, Вт.

В системе (1) параметры однородных тел «жила и ее изоляция» приняты идентичными.

Рассмотрим расчет тепловых параметров.

Тепловое сопротивление Я1 характеризует тепловой поток между центрами соседних жил, фактически идущий через прямоугольный параллелепипед, длина которого равна длине кабеля L, ширина равна г1, а толщина равна * (см. рис. 1). Тогда Я1 рассчитывается по аналогии с электрическим сопротивлением по формуле

« = % •

°С • м

где рг. - удельное тепловое сопротивление изоляции, -------.

Вт

Тепловое сопротивление R2 характеризует тепловой поток от центра жилы к внешней поверхности ее изоляции. Пренебрегая тепловым сопротивлением жил, для трехжильного кабеля с поясной изоляцией значение R2 можно рассчитать в соответствии с [7], [8]:

1 Р‘3

Г,. 3* V. da

где * - толщина изоляции между жилами (м); dx - диаметр круглой жилы с такой же площадью поперечного сечения и степенью уплотнения, что и фасонная жила (м); da - наружный диаметр поясной изоляции (м); г1 - радиус окружности, описанной вокруг жил (м).

Тепловое сопротивление R3 характеризует тепловой поток от поверхности изоляции жил через броню и наружное покрытие, т. е. через защитную оболочку. Пренебрегая тепловым сопротивлением брони, значение R3 также определим по [7], [8]:

рт 1п

R3 =---------------------------------------------------------------------4 г ау ,

3 2лЬ

°С • м

где рт - удельное тепловое сопротивление защитной оболочки, ----; * -толщина

Вт

защитной оболочки, м; Да - наружный диаметр брони, м.

Тепловое сопротивление R4• характеризующее тепловой поток между поверхностью кабеля и внешней окружающей средой, для кабеля, проложенного в воздухе, будем рассчитывать в соответствии с [9, с. 248]:

^=лДт •

где Д0 - наружный диаметр кабеля.

Теплоемкости всех однородных тел будем рассчитывать как произведения их удельных теплоемкостей на их массы, учитывая и диэлектрики, и металлы по аналогии с зависимостями, приведенными в [6].

Понижение порядка математической модели

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

Рассмотрим жилы кабеля и их изоляцию как единое однородное тело с потерями активной мощности, рассчитываемыми по формуле

Р =

где 11, 12, 13 - соответственно, токи жил кабеля; Кэ - электрическое сопротивление одной жилы кабеля.

Структурная схема тепловых процессов трехжильного кабеля, описываемых полученной упрощенной математической моделью, приведена на рис. 3.

Рис. 3. Структурная схема тепловых процессов трехжильного кабеля в соответствии с упрощенной математической моделью

Систему уравнений, соответствующих структурной схеме на рис. 3, для кабеля, проложенного в воздухе, можно записать в виде:

^01 | 0! -02 _ р.

1 Л 1 р

г

03 _ ^ 0*.

г3

В системе (7) теплоемкость С1 увеличена в 3 раза из-за суммирования теплоемкостей трех жил, а тепловое сопротивление ^ уменьшено в 3 раза из-за увеличения площади теплоотдачи полученного однородного тела «жилы и их изоляция» приблизительно во столько же раз (см. рис. 1). Как отмечено выше, теплоемкостью внутренней окружающей среды (С3) для кабеля, проложенного в воздухе, будем пренебрегать.

Тепловая схема замещения, синтезированная на основе системы (7), приведена на рис. 4.

0! 0,33^2 02 Rз 03 R4

Рис. 4. Тепловая схема замещения трехжильного кабеля по упрощенной математической модели

Таким образом, порядок математической модели понижается с четвертого до второго.

Расчетные выражения для температуры жил 01 и температуры поверхности кабеля 03 в операторной форме можно представить в виде

0 ( р ) ((Д2 (Дз + Д4 У )Р + Д + Дз + Д  )Р (^) + 04 (Р) .

'ч> (д;(Дз + д.)сс)р- + ((Дз + д.)(с;+с2)+Д2с;)р+1’

03 (Р ) =

д. р

__________з 1 2_______________________________

(Д2 (Дз + д. )С1С2 )р 2 +((Дз + д. )(с;+С2)+Д2с; )р+1

где

Д2=Д-; с;=з• с,.

Используя правила обратного преобразования Лапласа, алгоритм расчета температуры поверхности кабеля на основе соотношения (9) запишем следующим образом:

wi =

2,- =

— Д.

Щії’

- ((Д2 ДзСс )р,- + (Дз (с;+с,)+д;с; )р, +1).

8,,, = ’»'Р, +(

Нр)

-Дґ Р,) ^т;

Л,,, = 2І04,, +(Л,.; -1 - 2І04,, ) е ч;

-Д?

т,

0з, =Дз04.,

"з,,

+

¿к, +Л,., )>

Д + Д ^ V ^, ^ ^

“з ^“4 ,=1

где і - номер экспоненциальной составляющей теплового процесса; , - номер отсчета; Н(р) - характеристическое уравнение; рі - корни характеристического уравнения; т, - постоянные времени экспоненциальных составляющих; wi и г, -коэффициенты, рассчитываемые в соответствии с правилами обратного преобразования Лапласа; , и л, , - соответственно, составляющие температуры от активных

потерь в жиле и от температуры внешней охлаждающей среды по каждой экспоненциальной составляющей; Дt - интервал времени между j-м и j +1 отсчетами (период дискретизации).

Аналогично, на основе соотношения (8) можно записать алгоритм расчета температуры жил кабеля.

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

Рассчитаем температуры элементов трехжильного кабеля по исходной математической модели, решая систему (1) методом Рунге-Кутта четвертого порядка [10, с. 707], а также по модели второго порядка на основе (8)-(11).

В качестве кабеля выберем трехжильный кабель с секторными жилами и бумажной изоляцией ААШв 3х240-10. Геометрические параметры этого кабеля приведены в табл. 1. Результаты расчетов тепловых параметров даны в табл. 2.

Таблица 1

Геометрические параметры кабеля ААШв 3х240-10

Параметр Обозначение Значение

Наружный диаметр поясной изоляции ^а 43 мм

Радиус окружности, описанной вокруг жил Г1 18 мм

Толщина изоляции между жилами t 6 мм

Диаметр эквивалентной круглой жилы (Іх 17,5 мм

Наружный диаметр брони Da 48 мм

Толщина защитного покрытия 3 мм

Таблица 2

Результаты расчетов тепловых параметров кабеля ААШв 3х240-10

Параметр Rl R2 Rз R4 Сх С2

°С °С °С °С Вт • с Вт • с

Размерность Вт Вт Вт Вт °С °С

Значение 1,33 0,38 0,074 0,39 1629 2888

Рассмотрим наихудший случай, когда 11 = 270 А, 12 = 13 = 0 А и 04 = 40 °С. На рис. 5 показаны зависимости для процесса нагрева рассматриваемого кабеля в наихудшем случае.

Отклонение температур жил кабеля, рассчитанных по модели второго порядка, от этих же температур, рассчитанных по модели четвертого порядка, не превышает 2,1 °С. Для температуры поверхности кабеля это отклонение составляет не более 0,2 °С.

Таким образом, модель второго порядка является достаточно адекватной по отношению к модели четвертого порядка.

0, °С

ґ, ч

Рис. 5. Расчетные зависимости для процесса нагрева кабеля ААШв 3х240-10 при 11 =270 А, 12 = 13 = 0 А и 04 = 40 °С:

1, 2, 3 - температуры соответствующих жил кабеля по модели четвертого порядка; 4 - температура жилы кабеля по модели второго порядка;

5 - температура поверхности кабеля по модели четвертого порядка;

6 - температура поверхности кабеля по модели второго порядка

Методика косвенного измерения потерь активной мощности в кабеле

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

Для этого рассмотрим три последних уравнения алгоритма (11). Сгруппировав коэффициенты при параметре Р в третьем уравнении, получим:

( _Д Л 1 - е Ті

)

А?

Просуммировав соотношение (12) по і, можно записать:

2 2 ( >\,= Р, 2>,

+ >^

i, ,-1

Подставив (13) в последнее уравнение системы (11) и выразив Р, получим:

_Д/

2 R 0

0 ->, .е -і _>л - ^3°4,,

Х]' > г,j-1 >Лг,j R3 + R4

,=1

2 ( -Д Л

' 1 - е Хі

>

і=1

Ж-

V

)

Соотношение (14) позволяет рассчитывать неизвестные значения потерь активной мощности в трехжильном кабеле при известных значениях температуры его поверхности и температуры внешней окружающей среды. При этом полученные значения параметра Р^ будут соответствовать условным, эквивалентным потерям в

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

г=1

г=1

г=1

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

Значения эквивалентных потерь Р,, полученные на основе (14), можно использовать для расчета температуры жил кабеля.

Постановка эксперимента

На кафедре «Электроснабжение» ГГТУ им. П. О. Сухого на базе мощного понижающего трансформатора была собрана лабораторная установка для нагрева током одножильных и трехжильных силовых кабелей. Схема для проведения экспериментальных исследований тепловых процессов трехжильного кабеля и проверки адекватности разработанной адаптивной математической модели на основе этой установки показана на рис. 6.

ТА 1 - ТА 3

ип 1

и,

13 3

1Ы к Да

ІІ2 Д2

1Ь2 Д3

ш

Рис. 7. Схема для проведения эксперимента

Схема работает следующим образом. Фазные напряжения ип- иьъ подаются на входы трехфазного ЛАТРа Тг. Его выходы подключены к первичным обмоткам трехфазного трансформатора Т2 типа ТСЗВ мощностью 74 кВА, имеющего свободные окна между каждым ярмом сердечника и обмоткой. В эти окна проведены провода сечением около 100 мм2, соединенные с каждой из жил силового кабеля таким образом, что получились короткозамкнутые витки, рассчитанные на длительное пропускание тока значением до 300 А и на кратковременный ток - до 1500 А.

В качестве исследуемого образца выбран отрез кабеля ААШв 3х240-10 длиной 1 м. В просверленном до жилы отверстии размещен и плотно прижат к жиле терморезистор Кп, а само отверстие заделано высверленным материалом. Терморезистор 2 установлен на поверхность кабеля и закрыт сверху пробкой, заполненной ватой. На расстоянии 10 см от поверхности кабеля зафиксирован терморезистор 3.

С помощью экранированных сигнальных кабелей все терморезисторы подключены к соответствующим входам цифрового регистратора D2, разработанного и собранного на указанной кафедре. Вторичные обмотки измерительных трансформаторов тока ТАХ - ТА3, включенных в цепи короткозамкнутых витков, также соединены с

входами цифрового регистратора D2. Интерфейсные выводы этого регистратора подключены к СОМ - порту компьютера для обмена данными по интерфейсу RS232.

Измеренные величины, а также значение периода дискретизации на каждом шаге измерения записывались в компьютер.

Учет теплообмена через торцы кабеля

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

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

а аналогичное тепловое сопротивление брони - как R2 2. Считая, что тепловой поток

через эти тепловые сопротивления идет от центра кабеля к его торцам в обе стороны (а это эквивалентно снижению теплового сопротивления в 2 раза), можно записать:

0,5Рс — 0,5рс — р г

_ _ , Ус 2 _ » Ус 2 _Рс—

К2,1 Т7 72 _1 ,

я- с

где рс - удельное тепловое сопротивление жилы; L - длина кабеля; Fc - площадь сечения жилы; dc - диаметр жилы.

Аналогично для теплового сопротивления R2 2:

0,5Р б і 0-5Рб і Рб і

= ~^ = , „2 2.2Ч = ^ .2'

б

2 2 Fб ж{D; d; ї ЦрГ^Ї)

71 —----—

V4 4 у

где рб - удельное тепловое сопротивление брони; Da - наружный диаметр брони; da - наружный диаметр поясной изоляции.

Исходя из введенных выше определений для математической модели второго порядка, тепловые сопротивления R21 и R22 должны войти в состав однородного тела «жилы и их оболочка». Обозначим общее тепловое сопротивление этого однородного тела как К2. Исходя из принципа сложения тепловых проводимостей, для К2 можно записать:

К _п—1—-,

К + &1 + К,2

где К2 - исходное тепловое сопротивление однородного тела «жилы и их оболочка», рассчитываемое по формуле (10).

Результаты экспериментальных исследований

Было проведено 4 опыта. В табл. 3 показаны условия этих опытов.

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

На рис. 8, 9 приведены экспериментальные зависимости для опыта № 2, полученные по непосредственным измерениям, и расчетные зависимости на основе выражений (8)-(10) и (14) с учетом соотношений (15)-(17).

Таблица 3

Условия проведения опытов

Опыт Датчик температуры воздуха Обстановка в лаборатории Погодные условия

№ 1 Не закрыт фольгой Окна зашторены, дверь закрыта Переменная облачность

№2 Закрыт фольгой Окна зашторены, дверь закрыта Солнечно

№3 Закрыт фольгой Окна открыты, дверь закрыта Солнечно

№4 Закрыт фольгой Окна не зашторены, дверь закрыта Пасмурно

Рис. 8. Экспериментальная зависимость для действующих значений тока фазы «А» в опыте № 2

Рис. 9. Экспериментальные и расчетные зависимости для температур кабеля ААШв 3х240-10 в опыте № 2:

1 - экспериментальная диаграмма температуры воздуха;

2 - экспериментальная диаграмма температуры поверхности кабеля; 3 - экспериментальная диаграмма температуры жилы кабеля;

4 - расчетная диаграмма температуры жилы кабеля

Показатели качества математической модели, полученные путем сравнения измеренной и расчетной температуры жилы кабеля, даны в табл. 4.

Таблица 4

Показатели качества математической модели

Опыт Максимальная абсолютная погрешность, °С Коэффициент корреляции

№ 1 -2,9 0,998

№2 -2,4 0,993

№3 -3,2 0,998

№4 -2,2 0,996

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

Адаптивные свойства предложенной математической модели распространяются на все тепловые параметры кроме теплового сопротивления внутренней окружающей среды R4, непосредственно влияющего на тепловой поток с поверхности кабеля. При расчетах это сопротивление определялось по формуле (5), т. е. не учитывалась его зависимость от температуры. Очевидно, что в реальных условиях такая зависимость присутствует, так как на значение R4 влияют процессы конвекции и излучения [1].

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

Заключение

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

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

Литература

1. Schmidt, H-P. Efficient Simulation of Thermal and Electrical Behaviour of Industrial Cables / H-P Schmidt // Modelling and simulation. - June. - 2008. - P. 523-532.

2. Failure Modes of Underground MV Cables: Electrical and Thermal Modelling / P. Wallace ^t al.] // Proceedings of the COMSOL Conference. - 2009. - Milan ; ebookily.net.

3.Mamdooh, S. Adaptive Neuro-Fuzzy Inference System for Thermal Field Evaluation of Underground Cable System / S. Mamdooh // Journal of Energy and Power Engineering. - 2012. - № 6. - P. 1642-1650.

4. Millar, R. J. Real-Time Transient Temperature Computation of Power Cables Including Moisture Migration Modelling / R. J. Millar, М. Lehtonen // PSCC, Session 39. -Liege. - 2005. - P. 4.

5. Расчет переходных тепловых режимов одножильного силового кабеля с пластмассовой изоляцией, проложенного в воздухе / Э. Т. Ларина [и др.] // Электротехника. -1991.-№ 10.-С. 39-42.

6. Зализный, Д. И. Математическая модель тепловых процессов одножильного силового кабеля / Д. И. Зализный, С. Н. Прохоренко // Изв. высш. учеб. заведений и энергетич. объединений СНГ. - Энергетика. - 2012. - № 5. - С. 25-34.

7. ГОСТ Р МЭК 60287-2-1. Кабели электрические. Расчет номинальной токовой нагрузки. Ч. 2-1. Тепловое сопротивление. Расчет теплового сопротивления. - М. : Стандартинформ, 2009. - 19 с.

8. Капанский, А. А. Автоматизация расчета теплового сопротивления кабеля / А. А. Ка-панский, Д. И. Зализный // Вестн. Гомел. гос. техн. ун-та им. П. О. Сухого. - 2013. -№ 4. - С. 87-94.

9. Ларина, Э. Т. Силовые кабели и высоковольтные кабельные линии : учеб. для вузов / Э. Т. Ларина. - М. : Энергоатомиздат, 1996. - 464 с.

10. Корн Г. Справочник по математике (для научных работников и инженеров) / Г. Корн. - М. : Наука, 1974. - 832 с.

Получено 17.02.2014 г.