Определение угловых коэффициентов излучения в программном комплексе MSC Marc Текст научной статьи по специальности «Математика»

Информатика, вычислительная техника и управление. Моделирование. Приборостроение. Метрология. Информационно-измерительные приборы и системы

УДК 621.33 Дульский Евгений Юрьевич,

аспирант кафедры «Электроподвижной состав», Иркутский государственный университет путей сообщения,

тел. 8-983-403-46-43 Гарев Николай Николаевич, аспирант кафедры «Электроподвижной состав», ИрГУПС, тел. 8-914-902-25-36

Доценко Наталия Сергеевна, аспирант кафедры «Электроподвижной состав», ИрГУПС, тел. 8-902-578-77-16

ОПРЕДЕЛЕНИЕ УГЛОВЫХ КОЭФФИЦИЕНТОВ ИЗЛУЧЕНИЯ В ПРОГРАММНОМ КОМПЛЕКСЕ MSC MARC

E.YU. Dulskiy, N.N. Garev, N.S. Docenko

DEFINING RADIATION VIEW FACTORS IN THE MSC MARC

SOFTWARE COMPLEX

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

Ключевые слова: изоляция, капсулирова-ние, теплообмен излучением, угловые коэффициенты.

Abstract. This article describes the main methods of determining the view ^ factors used in solving the problems of radiative heat transfer in the MSC Marc software package.

Keywords: isolation, encapsulation, heat radiation, view ^ factors.

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

Применение данного комплекса позволило смоделировать режимы ИК-энергоподвода в технологии капсулирования изоляции лобовых частей обмоток якоря ТЭД НБ-514Б электровозов серии «ЕРМАК», о которой упоминалось нами ранее в других работах [1, 2], с дальнейшей целью определения наиболее эффективных режимов.

Процесс моделирования режимов ИК-энергоподвода в технологии капсулирования изоляции лобовых частей обмоток якоря ТЭД НБ-514Б электровозов серии «ЕРМАК» в программном комплексе MSC Marc можно разделить на следующие этапы:

1) создание 3Б-модели;

2) разбиение модели на конечноэлементную

сетку;

3) предание конечным элементам свойств реальных материалов;

4) задание граничных условий и нагрузок;

5) анализ (решение);

6) графическая визуализация решения.

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

Программный комплекс MSC Marc при расчете УКИ использует следующие основные методы по их нахождению [3]:

1) метод прямой адаптивной интеграции;

2) метод Монте-Карло;

3) точечный метод полукуба.

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

1 rrcOS фcOS ф ,

dF(?j= A Jr

'-dAidAi

(1)

A,A,

nr

где А, А}- - площади поверхностей (рис. 1):

фг, ф;- - углы между излучающими поверхностями;

г - расстояние между излучающими поверхностями.

Рис. 1. Определение УКИ между двумя поверхностями

ИРКУТСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ПУТЕЙ СООБЩЕНИЯ

A • F = A- • Fi •

ном случае имеет вид

q, = A, • (Ыг — G г)

О, =

Мг — е,- Em M, — е, • En

Рг

1 — е,

где Ет■ - мощность излучения;

8, -коэффициент излучения; р - коэффициент отражения. Подставляя уравнение (7) в уравнение (4),

после перестановки слагаемых получаем:

М = E„. — ■

1 — е,

A, -е,

• qг •

A .8 , — (1 — е,) • A, • F

A- -е-

•{q- }=[ ^8 , — A, • F, ] }(9)

Количество излучаемого теплового потока, передаваемого между двумя поверхностями, будет зависеть от того, какая часть излучения от каждой поверхности падает на другую поверхность. Как показано на рис. 1, излучение, распространяющееся от поверхности А, на поверхность A-, будет равно:

cos®, • cos®, • Ai q, = M , • A, • (-^ ^, ^ ,) = M , • A, • F,, (2)

где q,j - тепловой поток между поверхностью At

и A

М- энергетическая светимость; А,, A, - площади поверхностей; г - расстояние между поверхностями; F, - УКИ между поверхностями А, и Ay'. Так как УКИ имеет исключительно геометрический характер, справедливо равенство:

(3)

Для задач теплового излучения черного тела, то есть при е, = 1, данное уравнение примет вид

{q, }= [A , 8, — A, • F, ].{E, }• ('°>

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

Тепловой поток, излучаемый поверхностью A1 на A2, вычисляется как

A , cos ® ,

q, = M, • (A, • cos) • (^-a-), (11)

Уравнение теплообмена излучением в дан-

(4)

п • г2

где фу - углы между нормалями т, , тj и осью г (рис. 1).

В уравнении (11) первым слагаемым в скобках является проекция А, нормали к соединительной линии, а второй член - телесный угол, под которым видна поверхность Ау от центра поверхности А .

УКИ в данном случае определяется как

где О - тепловой поток в обратном направлении.

Для О может быть сформировано два независимых друг от друга выражения:

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

N N

А■ • О = 1 Му • Ау • Р, = 1 Му • А • Ру , (5)

у=1 у=1

или

N

О = 1 Му • Ру , (6)

у=1

где N - число поверхностей, участвующих в вычислениях;

б) другое выражение для О:

F, = (■

cos ®, • cos ® ,

п г

) • Af

(1,)

Подставляем уравнение (12) в (11):

Чу = М • А, • Ру, (13)

При расчете данным методом происходит разбиение площадей излучаемых поверхностей на фасетные поверхности. Чем сложнее геометрия (как в нашем случае с ИК-излучателем и сегментом лобовой части), тем больше фасетных поверхностей и, соответственно, сложнее расчеты (рис. 2).

(7)

(8)

Подставляя данное уравнение в уравнение (5), получаем основное уравнение для задач теплового излучения диффузно-серого тела:

Рис. 2. Построение фасетных поверхностей

Формула по нахождению УКИ с учетом рис. 2 имеет вид

яА1 • йРф12 = + dA1б) • с?Рф12 = dA1a • (рф1а2а + dРф1a2б) +

+ ¿А1б • (<ЛРф1б2а + ^Рф1б2б ) +....., (14)

где dA1a и йА1б - соответственно площади элементарных площадок фасетной поверхностности;

Информатика, вычислительная техника и управление. Моделирование. Приборостроение. Метрология. Информационно-измерительные приборы и системы

dFtyu - средний угловой коэффициент излучения поверхности 1 на поверхность 2 (рис. 2);

<^ф1а2а, <^ф1а2а, ^ф1а2а, и ^ф1а2а - соответственно угловые коэффициенты излучения элементарных площадок фасетной поверхностной сетки 1 и 2.

Идея метода Монте-Карло заключается в задаче случайного испускания ИК-лучей от одной поверхности А\ к другой поверхности А2 [3]. Доля лучей, которые попали на поверхность А2, складывается в форм-факторе между ними. Метод Монте-Карло вычисляет N форм-факторов одновременно, обеспечивая линейное масштабирование. На самом деле чем больше число поверхностей, тем быстрее форм-факторы вычисляются по сравнению с методом прямой адаптивной интеграции. Следовательно, в MSC Marc этот метод принят для расчетов УКИ по умолчанию.

Некоторые из особенностей расчета УКИ в MSC Marc с использованием метода Монте-Карло:

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

- Трудоёмкость расчета почти линейно пропорционально количеству элементов. Это означает, что при обширных задачах трудоёмкость существенно не возрастёт.

- Сумма УКИ всегда равна 1.

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

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

1) Поверхность должна иметь возможность испускать ИК-лучи случайным образом. Направление луча должно быть задано в соответствии с косинусом угла между лучом и нормалью к поверхности, проведенной из начала луча, что обусловит рост выбросов ИК-лучей по нормали к поверхности. Источник ИК-луча должен быть случайным образом распределен по площади поверхности (рис. 3).

2) Учитывая происхождение, направление и длину луча, поверхность должна автоматически определять есть ли пересечение с данным лучом, и если есть, то рассчитывать его длину в точке пересечения.

Рис. 3. Распределение ИК-лучей в методе Монте-Карло 1 - поверхность; 2 - ИК-лучи

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

На рис. 4 изображено образование тени 4, что по существу представляет собой вычисление падающего излучения путем трассировки ИК-луча 3 от источника излучения 1 до точки падения ИК-луча на облучаемую поверхность 5. Эти ИК-лучи (теневые) не отражаются и не преломляются. Тени образуются, когда путь ИК-лучей от источника излучения до облучаемой поверхности частично или полностью преграждается другим объектом. Преломление и отражение излучения от самих источников не вычисляется.

1

Рис. 4. Образование тени в процессе трассировки ИК-лучей: 1 - источник излучения, 2 - геометрическое препятствие; 3 - ИК-лучи; 4 - тень;

5 - облучаемая поверхность

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

ИРКУТСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ПУТЕЙ СООБЩЕНИЯ

Рис. 5. Исходная модель

таг К?

ч

Рис. 6. Процесс рекурсивного разбиения модели

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

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

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

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

Наиболее быстрым и точным в расчетах УКИ является точечный метод полукуба. Данный метод является модификацией метода сферы единичного радиуса, предложенного В. Нуссельтом

[4]. Рассмотрим сначала его.

Если над элементарной площадкой dA1 (рис. 7) построить полусферу единичного радиуса

[5], то УКИ между dAi и некоторой поверхностью А2, согласно уравнению (1), будет равен:

dF12 = - ■ f cos р1 C0S в2 dA2 = - • f cos Pdv (15)

ТГ J е ТГ »

Рис. 7. Геометрическая схема для определения УКИ методом сферы единичного радиуса

Заметим, что d®1 — проекция dA2 на поверхность полусферы, так как

dAS cos В2dA2

d(o1 =-iL = dAS =-Ч^, (16)

r s

где r - радиус единичной полусферы. С учетом

этого УКИ рассчитывается как

1 Г

dF1_2 = — ■ I cosP1dAS.

(17)

Однако dAS cosp1 - проекция dAs на основание полусферы. Следовательно, интегрирование

S

Информатика, вычислительная техника и управление. Моделирование. Приборостроение. Метрология. Информационно-измерительные приборы и системы

оо8р1 дает проекцию Аъ поверхности Л$ на основание полусферы, или

1 г Л

й¥1-2 = —■ I еоБ= —. п • п

(18)

Главный недостаток данного метода -большие затруднения при проецировании поверхности А2 на полусферу. С этой целью, а также для увеличения скорости расчета УКИ в методе полукуба сфера заменяется на куб (рис. 8), для участков граней 1 которого заранее определены УКИ. При проецировании поверхности А2 на поверхность куба, находится сумма всех УКИ участков граней полукуба 2, на которые попадает проекция поверхности А2, хотя при этом для сложной геометрии теряется качество расчета. Сумма этих УКИ дает угловой коэффициент всей поверхности.

Рис. 8. Геометрическая схема для определения УКИ методом полукуба: 1 - участок грани полукуба; 2 - участок грани полукуба, участвующий в расчете УКИ

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

Основной алгоритм программы расчета УКИ по данному методу выглядит следующим образом:

а) считывание файла, созданного программой при анализе, содержащего геометрию модели;

б) излучение поверхности охватывает контролируемый пользователем угол (для осесиммет-ричной геометрии);

Рис. 9. Метод полуплоскости

в) в 3Б происходит разбиение четырехугольных поверхностей на треугольные;

г) создание карты неоднородных точек и соотношение их с определенным весом;

д) преобразование координат в локальной системе;

е) цикл по приемной поверхности, линии в 2Б, треугольники в 3Б;

ж) в 2Б происходит проецирование линий на линии, а в 3Б - треугольников на плоскости;

з) определение точек (количество и расположение);

и) суммирование веса точек;

к) расчет УКИ этой поверхности.

Когда присутствуют плоскости симметрии, принимающие поверхности удваиваются для каждой из них. Для циклической симметрии принимающие поверхности удваиваются в п раз.

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

Метод полукуба используется и при расчете УКИ в специализированном решателе М8С 8т-daRad. Этот решатель позволяет не только рассчитывать УКИ, но и визуализировать результаты расчетов (рис. 10).

а)

б)

Рис. 10. Визуализация расчета УКИ в системе «ИК-излучатель - сегмент лобовой части обмотки в решатели МЗС 8пн1аН;к1: а - крупная конечно-элементная сетка, б - мелкая конечноэлементная сетка

Результаты моделирования режимов ИК-энергоподвода технологии капсулирования изоляции лобовых частей обмоток якоря ТЭД НБ-514Б электровозов серии «ЕРМАК» в программном комплексе MSC Marc будут представлены в ближайшее время в дальнейших работах.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Худоногов А.М. Инновационная технология повышения надёжности и продления ресурса электрических машин тягового подвижного состава // Современные технологии. Системный анализ. Моделирование. 2012. № 4 (36). С. 261.

2. Дульский Е.Ю. Анализ пространственного распределение инфракрасного излучения в процессе капсулирования изоляции электрических машин тягового подвижного состава // Вестник ИрГТУ. 2013. № 7 (78). С. 314.

3. Marc®2012. Volume A : theory and user information // MSC Software, 2012. 813 р.

4. Nusselt W., Graphicshe Bestmmung des Winkel-verhaltnisses bei der Warmestrahlung, VDIZ. 72. 673 c.

5. Зигель Р., Хауэль Дж. Теплообмен излучением : пер. с англ. Москва, 1975. 936 с.