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

кандидата технических наук
Леонов, Виктор Витальевич
город
Москва
год
2011
специальность ВАК РФ
05.13.18
Диссертация по информатике, вычислительной технике и управлению на тему «Разработка математической модели концентратора высокотемпературной солнечной энергоустановки космического назначения»

Автореферат диссертации по теме "Разработка математической модели концентратора высокотемпературной солнечной энергоустановки космического назначения"

ЛЕОНОВ ВИКТОР ВИТАЛЬЕВИЧ

РАЗРАБОТКА МАТЕМАТИЧЕСКОЙ МОДЕЛИ КОНЦЕНТРАТОРА ВЫСОКОТЕМПЕРАТУРНОЙ СОЛНЕЧНОЙ ЭНЕРГОУСТАНОВКИ КОСМИЧЕСКОГО НАЗНАЧЕНИЯ

Специальность 05.13.18 - Математическое моделирование, численные методы и комплексы программ

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата технических наук

-Б ОКТ 2011

Москва-2011

4855264

Работа выполнена на кафедре прикладной математики Московского государственного технического университета имени Н.Э. Баумана

Научный руководитель:

доктор технических наук, профессор Кувыркин Георгий Николаевич

Официальные оппоненты:

доктор технических наук Градов Владимир Михайлович

кандидат технических наук Гавриш Сергей Викторович

Ведущая организация:

ФГУП «НПО им. С .А. Лавочкина»

Защита состоится « 18 » октября 2011 года в 13 часов 00 минут на заседании диссертационного совета Д 212.141.15 при Московском государственном техническом университете имени Н.Э. Баумана по адресу: 105082, г. Москва, Рубцовская наб., д. 2/18, ауд. 1006.

Отзывы на автореферат в двух экземплярах, заверенные печатью, просим отправлять по адресу: 105005, г. Москва, 2-я Бауманская ул., д. 5, МГТУ им. Н.Э. Баумана, учёному секретарю совета Д 212.141.15.

С диссертацией можно ознакомиться в библиотеке Московского государственного технического университета имени Н.Э. Баумана.

Автореферат диссертации разослан « /Ф » Cxtru<_»eS^iSi 2011г.

Ученый секретарь диссертационного совета, кандидат технических наук, старший научный сотрудник, доцент ^ Атгетков A.B.

Общая характеристика работы

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

В настоящее время солнечная энергетика продолжает активно развиваться. Если на конец 2009 года она занимала в структуре мирового производства электроэнергии около 1%, то к середине XXI века, по прогнозам экспертов Международного энергетического агентства (1ЕА), при сохранении современной динамики развития, может достигнуть 25%. В космосе более 98% космических аппаратов используют низкотемпературные СЭУ, работающие с применением фотоэлектрических преобразователей. Разработка новых проектов по созданию космических солнечных электростанций, космических аппаратов с применением солнечных и электрореактивных ракетных двигателей, систем для освещения приполярных областей, энергосистем для космических станций и баз как на орбите, так и поверхности других планет, требует разработки и создания СЭУ большей мощности по сравнению с существующими. Как показывает опыт эксплуатации СЭУ, в основном солнечных электростанций, для создания систем большой мощности наиболее эффективно использование высокотемпературных СЭУ (ВТСЭУ) с применением высокопотенциальных концентраторов солнечной энергии (параболических, сферических, параболоцилиндрических) или зеркальных концентрирующих систем (ЗКС), позволяющих значительно повысить плотность солнечной энергии в рабочей зоне.

Проектирование ВТСЭУ космического назначения ведётся с конца 70-х годов прошлого века. За пройденный период был разработан ряд проектов как по созданию подобных систем в целом, так и её элементов; основное внимание уделялось проектированию и созданию концентраторов солнечной энергии. Был разработан ряд методов по определению характеристик подобных систем. При этом, в отличие от задач и методов классической оптики, направленной на получение качественного изображения объекта, в задачах, связанных с проектированием ВТСЭУ, в первую очередь интересуют энергетические характеристики сконцентрированного излучения (плотность потока, распределение температуры по поверхности концентратора и приёмника) и связанная с ними эффективность системы концентратор-приёмник (СКП). Из-за невозможности в то время создать крупногабаритные концентраторы с приемлемым соотношением точностных и массовых характеристик подобные системы в космосе не были реализованы.

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

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

Экспериментальное определение радиационных характеристик или исследуемых параметров требует проведения сложных экспериментов, имеющих значительные ограничения на интерпретацию результатов (зависимость от расстояния до Солнца, спектра излучения, точности изготовления поверхности зеркал). Поэтому особый интерес представляет разработка математической модели, позволяющей рассчитывать характеристики сложных зеркальных систем с учётом влияния условий эксплуатации, конструктивных особенностей, шероховатости поверхности. Математическое моделирование позволит обеспечить значительное уменьшение материальных затрат, связанных с проектированием и отработкой подобных систем, по сравнению с традиционными экспериментальными методами. А возможность рассмотреть большое число вариантов в результате математического моделирования приведёт к разработке рациональной конструкции при обеспечении максимальной энергетической и массовой эффективности ВТСЭУ, что очень актуально для космической техники.

Целью работы является разработка математической модели, позволяющей рассчитывать радиационные характеристики ЗКС, работающей в составе ВТСЭУ космического назначения.

Для достижения этой цели были поставлены и решены следующие задачи.

1. Разработка математической модели, позволяющую исследовать распределение радиационных потоков в ЗКС различной конфигурации с учётом влияния возможных погрешностей, характерных для рабочей поверхности.

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

3. Проведение верификации разрабатываемой математической модели, алгоритмов и программного комплекса на тестовых задачах.

4. Исследование влияния погрешностей рабочей поверхности зеркала на радиационные и энергетические характеристики СКП, работающей в составе ВТСЭУ.

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

Научная новизна. В диссертации получены следующие основные результаты, которые выносятся на защиту.

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

2. Комплекс прикладных программ «Tracer», реализующий разработанную математическую модель.

3. Результаты анализа влияния погрешностей рабочих поверхностей на энергетические характеристики СКП, работающей в составе ВТСЭУ.

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

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

Апробация результатов работы. Основные положения и результаты диссертационной работы докладывались и обсуждались на Общероссийском форуме «Использование космоса в мирных целях» (Москва, 2007), VI Международной научной конференции «Ракетно-космическая техника» (Москва, 2007), XXXII -XXXV Академических чтениях по космонавтике (Москва, 2008-2011), Всероссийской конференции молодых учёных и специалистов «Будущее машиностроения России» (Москва, 2008), II Международной научно-технической конференции «Аэрокосмические технологии», посвященной 95-летию со дня рождения академика В.Н. Челомея, (Реутов-Москва, 2009), Международной научной конференции UKSEDS (Гилфорд, 2009; Манчестер, 2011), Международном научном семинаре «Разработка космических систем: взгляд из EPFL и МГТУ», (Лозанна, 2009-2010),

XVII Школе-семинаре молодых учёных и специалистов под руководством академика РАН А.И. Леонтьева «Проблемы газодинамики и тепломассообмена в аэрокосмических технологиях» (Жуковский, 2009), X-XI Всероссийской научно-технической конференции и школе молодых учёных, аспирантов и студентов «Научные исследования в области транспортных, авиационных и космических систем» (Воронеж, 2009; Таруса, 2010), Международном научном семинаре Space Station Design Workshop, (Штутгарт, 2009), IV Всероссийской молодежной научно-инновационной школе «Математика и математическое моделирование», (Саров, 2010), V Российской национальной конференции по теплообмену (Москва, 2010),

XVIII Школе-семинаре молодых учёных и специалистов под руководством академика РАН А.И. Леонтьева «Проблемы газодинамики и тепломассообмена в новых энергетических технологиях» (Звенигород, 2011).

Публикации. Основные положения и результаты диссертационной работы опубликованы в 17 печатных и электронных изданиях, из них в трёх статьях из Перечня рецензируемых ведущих научных журналов и изданий. Получено 1 свидетельство о государственной регистрации программы на ЭВМ.

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

Структура и объём работы. Диссертация состоит из введения, 4 глав, выводов и списка литературы. Работа изложена на 150 страницах, содержит 96 иллюстраций. Библиография включает 106 наименований.

Содержание работы

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

В первой главе проведён анализ конструкций ВТСЭУ и её элементов, их основных характеристик, методов их определения, применявшихся ранее и используемых в настоящее время.

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

т]кл = Л/т/Л/с.

Величина зависит от геометрических и оптических параметров всех элементов СКП, а также учитывает потери тепловой мощности приёмником из-за излучения в результате нагрева.

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

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

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

2

падающее излучение (О , Д, X, 3)

отраженное излучение (0,,!!г вгГ{^К Т. Ьш)

У

,- отражающая поверхность {Т, Иш, тш)

Рис. 1. Моделирование отражения

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

отраженного излучения к энергии а!£( = /; • со^^, п) йш падающего излучения в пределах телесного угла с1ш

сИ^.ы^) РУОч.^Г) = 77-^-}-—,

где 5 = (в,Р) - направление излучения (траектория пучка) (рис. 1), Я - длина волны, Т - температура поверхности, <5 - угловой размер источника излучения, I -интенсивность излучения, <о - телесный угол, - среднеквадратичная высота шероховатости, шш - шаг шероховатости. Индексами /' и г обозначены величины, относящиеся соответственно к падающему и отражённому излучениям.

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

Анализ методов расчёта распределения радиационных тепловых потоков в ЗКС показал, что необходим метод, позволяющий учитывать влияние локальных дефектов и температурных полей, особенностей распространения излучения (зависимость от длины волны, направления падения и отражения излучения на каждом участке рабочей поверхности). Этим требованием будет полностью удовлетворять «пучковая» модель радиационного теплообмена, построенная на принципах статистического моделирования с использованием методов Монте-Карло.

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

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

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

В соответствии с основной идеей метода Монте-Карло для приближенного вычисления величины а необходимо «придумать» такую случайную величину X, чтобы математическое ожидание МХ = а. При этом сама величина X может быть функцией скалярной или векторной случайной величины, или даже функционалом от случайного процесса.

Случайные величины нужного распределения обычно моделируются с помощью преобразований одного или нескольких независимых значений случайного числа К, равномерно распределённого в интервале от 0 до 1. Поэтому первоочередной задачей при использовании метода Монте-Карло, от решения которой в значительной степени зависят результаты расчета и их достоверность, является задача моделирования (генерирования) таких случайных величин.

В результате анализа используемых генераторов случайных чисел для решения поставленных задач был выбран генератор псевдослучайных чисел, разработанный Макото Мацумото и Такудзи Нисимура и получивший название «Вихрь Мерсенна». Вихрь Мерсенна является витковым регистром сдвига с обобщённой отдачей, обеспечивает быструю генерацию высококачественных псевдослучайных чисел (в частности, в 2-3 раза быстрее линейных конгруэнтных генераторов) с равномерным распределение генерируемых чисел в интервале от 0 до 1 в 623 измерениях (для линейных конгруэнтных генераторов оно ограничено 5-ю измерениями) и периодом

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

где Я - длина волны; Г - температура поверхности; й - постоянная Планка, с0 -скорость света в вакууме; к - постоянная Больцмана; <т - постоянная Стефана-Больцмана.

Кумулятивная функция К (Л, Г) может принимать любое значение на отрезке [0;1] и определяет вероятность события, происходящего в пределах от 0 до Я. При этом Я задаётся как число, генерируемое случайным образом, а соответствующее значение Я определяется с использованием функционального соотношения К (Я, Г) при заданной температуре.

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

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

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

равным, 219937-1.

л

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

Рис. 2. Отклонение нормали (п -> п') к поверхности треугольного элемента в точке М, моделирующее грань £¿5 шероховатой поверхности

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

где 1Г и - направляющие векторы отражённого и падающего луча соответственно, п - вектор нормали. А изменение направления отражения луча из-за шероховатости можно учитывать добавлением значения среднеквадратичного угла наклона шероховатости к нормали поверхности в точке падения пучка излучения в соответствии с параметрами шероховатости (рис. 2). Данный подход аналогичен методу гауссовой нормали, предложенному В.А. Грилихесом. Изменение отражательной способности поверхности в зависимости от углов падения и отражения излучения учитывается с помощью соотношения Г. Дэвиса

где рЧб;) - зеркальная отражательная способность шероховатой поверхности, Ро(Ю - отражательная способность идеальной поверхности того же материала, кш/Х -отношение среднеквадратичной высоты шероховатости к длине волны излучения, в, -угол падения. Это соотношение было основано на статистическом подходе к рассмотрению отражения электромагнитных волн от шероховатой проводящей поверхности и предполагает, что вещество является идеальным проводником, распределение высот неровностей поверхности является гауссовым около среднего значения, а отражения являются однократными. Поэтому оно пригодно для материалов, у которых отражение происходит очень близко от поверхности, например для металлов, и не применимо для диэлектриков, у которых отражение происходит не только на поверхности, но также и внутри материала.

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

3

2

1Г = 1£ — 2п(п ■

КХХ — 1Х + Куу — 1у + К„ — 1Ж + А(Г - Г») + <7 = О, (2)

Уравнение стационарной теплопроводности в твёрдом теле имеет вид д2Т д2Т д2Т

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

дТ дТ дТ

где /г - коэффициент теплообмена; Т - температура на граничных поверхностях тела;

- температура окружающей среды; 1Х, 1у и \г - направляющие косинусы вектора нормали к поверхности; д - плотность теплового потока.

Решение уравнения (1) с граничными условиями (2) для СКП произвольной формы удобно искать, переходя к МКЭ. Возможность использования МКЭ закладывалась при создании геометрической модели СКП и поэтому данный переход не требует дополнительных вычислительных затрат.

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

В данном случае для решения квазидвумерной задачи удобно использовать треугольные трёхузловые конечные элементы (КЭ). При этом граничные условия (2) задаются на передней и задней гранях КЭ. Совместное решение получаемых в МКЭ матричных уравнений, записанных для каждого КЭ, позволяет определить распределение температур по рабочим поверхностям СКП.

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

Таким образом, алгоритм моделирования радиационного теплообмена состоит из следующих шагов.

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

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

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

3. Задание начальных и граничных условий в зависимости от целей исследования.

4. Генерация пучков излучения. Для этого в соответствии с подходом, присущим методам Монте-Карло, производится независимая выборка из равномерно распределенного ряда чисел от нуля до единицы случайных чисел Л, обеспечивающих заданный закон распределения спектральных и энергетических характеристик излучателя.

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

5.1. Определение точки пересечения траектории пучка с геометрической моделью.

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

- Если это число удовлетворяет условию

О < Яа < а, то падающая порция энергии поглощается.

- Если случайное число Ка удовлетворяет условию

а < Яа < 1,

то рассматриваемая порция энергии будет отражена поверхностью в точке падения.

5.3. Определение направления отражения в зависимости от оптических свойств поверхности в точке падения:

- если поверхность зеркальная, то отражение происходит по закону Декарта-Снеллиуса;

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

- для случая диффузного или направленного отражения выбирается соответствующая распределению кумулятивная функция. Например, для случая диффузного отражения

бш2 в = Дв и /? = 2пЯр.

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

7. Решение тепловой задачи с целью определение полей температур производится МКЭ в квазидвумерной постановке по данным, полученным на предыдущем шаге.

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

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

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

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

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

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

Для проверки правильности работы блока, реализующего математическую модель радиационного теплообмена, были рассмотрены диффузные угловые коэффициенты (перпендикулярные пластины, концентрические сферы, параболоид), расчёт результирующего теплового потока при заданных температурах (концентрические сферы), расчёт температуры селективной поверхности при заданном тепловом потоке (параболоид). Было проведено сравнение полученных в КПП значений индикатрис отражения с экспериментальными данными и результатами расчётов, приведённых в работах Дж. Хауэлла.

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

Результаты тестирования показывают, что разработанная математическая модель и реализующий её КПП «Tracer» являются работоспособными, в результаты их применения - достоверны.

В четвёртой главе проведено исследование влияния возможных дефектов отражающей поверхности на её радиационные характеристики. Исследованы радиационные и энергетические характеристики СКП ВТСЭУ в зависимости от геометрических и оптических характеристик отражающей поверхности с учётом влияния дефектов.

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

Рис. 3. Индикатрисы отражения для плоской пластинки: а) при различных углах падения излучения (Лш = 3.2 мкм); б) при различных значениях коэффициента отражения р (кш = 2.0 мкм)

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

Рис. 4. Распространение потока падающего излучения (0; = 0) в полости микрометеоритного кратера (а) и индикатриса отражения от её поверхности (б)

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

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

1,0 0,8

о-

0,4 0,2 0,0

-0,06 -0,04 -0,02 0 0,02 0,04 0,06

Г = Г /гк

Рис. 6. Распределение относительной плотности потока излучение по плоскости входного отверстия приёмника: а) 2-х мерное (метод Апариси); б) 3-х мерное (моделирование)

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

Рис. 5. Геометрическая модель (а) и схема (б) СКП, состоящей из параболического концентратора и сферического полостного приёмника: /- фокусное расстояние, «м - угол полураскрытия концентратора, гк, г0 и г*. -соответственно радиусы концентратора, входного отверстия приёмника и фокального пятна, 8 - угловой размер источника излучения (для Земли 5 = 32'), 8' - угловой размер источника излучения с учётом погрешности зеркала

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

а) 12 .. . б) ц/Ч"

н = 60°

Е0р БШ2 цм'

ЛА = 3 • N

можно определить ЛА. Здесь р - коэффициент отражения рабочей поверхности (зеркала) концентратора, Е0 - солнечная постоянная, ик ~ угол полураскрытия концентратора (рис. 5,6).

Такой подход позволит с малыми вычислительными затратами приближённо определять энергетические характеристики СКП различных габаритов при постоянной интегральной точности (мере точности) зеркала концентратора.

По полученному распределению тепловых потоков в рассматриваемой геометрической модели было вычислено поле температур на концентраторе (рис. 7).

-Э457

= 3469

= 1225 223

Рис. 7. Поле температур на концентраторе: а) при изолированной задней стороне без учёта теплопроводности; б) без изоляции с учётом теплопроводности

Разработка ВТСЭУ высокой мощности (МВт и более) потребует создания крупногабаритных концентраторов. Создание подобных конструкций в космосе возможно принципиально двумя способами:

1) сборка (производство) из привезённых элементов или материалов;

2) развёртывание (раскрытие) конструкции из транспортного положения в рабочее.

Реализацию второго способа можно обеспечить применением многозвенных

(стержневых, лепестковых, зонтичных), бескаркасных (роторных) и надувных конструкций.

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

Был проведён анализ влияния изменения геометрии элементов рассматриваемой СКП на её геометрический КПД т)г. Анализ показал, что рационально использовать концентраторы с углом полураскрытия из диапазона от 60° до 70° (рис. 8,6). Это позволит обеспечить высокий геометрический КПД СКП при отсутствии переотражения излучения концентратором на себя. Это утверждение согласуется с результатами экспериментов с различными СКП и опыта эксплуатации ВТСЭУ, в которых было показано, что наилучшими характеристиками обладают концентраторы с углом полураскрытия 60-65°.

0,01 0,02 0,03 Л) = го/''к

0,04 0,05 10 20 30 40 50

60 70 80 90

Рис. 8. Зависимость геометрического КПД: а) от относительного радиуса входного отверстия приёмника; б) от угла полураскрытия концентратора

Кроме потерь, вызванных несовершенством геометрической формы отражателя и дефектов отражающей поверхности, уменьшение эффективности СКП происходит за счёт энергетических потерь, вызванных нагревом и переизлучением (вторичным излучением) её элементов, в первую очередь приёмником (рис. 9). Эти потери позволяет учесть энергетический КПД (рис. 10).

: ------------------------ - . -

0.9 ....... ........................„..... •

0,02

0,05

0,06

Рис. 9. Собственное излучение приёмника

0,03 0,01

'Ъ = 'о/г*

Рис. 10. Зависимость энергетического КПД от относительного радиуса входного отверстия приёмника

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

Таким образом, можно утверждать, что для исследуемой СКП существует оптимальное соотношение габаритов концентратора (его радиуса) и приёмника (радиуса его входного отверстия). Причём если при рабочей температуре Т = 1300 К это соотношение можно выбирать из диапазона г0 = 0.02... 0.03, ввода дополнительные инженерные и технологические ограничения без ощутимых потерь в КПД, то при Т = 1900 К этот диапазон сужается уже до г0 = 0.02 ... 0.022.

Основные результаты работы

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

2. Разработан комплекс прикладных программ «Tracer», реализующий разработанные математическую модель и алгоритмы.

3. Проведена верификация математической модели и комплекса прикладных программ, подтвердившая их работоспособность.

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

5. Проведён анализ влияния геометрических параметров и различных дефектов поверхности на радиационные характеристики зеркальных концентрирующих систем на примере параболического концентратора.

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

Основные результаты диссертации опубликованы в работах

1. Свидетельство о государственной регистрации программы для ЭВМ №2011611669. Tracer - определение радиационных и энергетических характеристик зеркальных концентрирующих систем / В.В. Леонов. Зарегистрировано в Реестре программ для ЭВМ 18.02.2011.

2. Леонов В.В. Перспективы использования солнечных электростанций для освоения Луны // Актуальные проблемы Российской космонавтики: Труды XXXII Академических чтений по космонавтике. М., 2008. С. 385-387.

3. Леонов В.В. Проектирование и расчёт характеристик крупногабаритных концентраторов солнечной энергии высокотемпературных солнечных энергоустановок космического назначения // Будущее машиностроения России: Сборник трудов Всероссийской конференции молодых ученных и специалистов. М„ 2008. С. 230-231.

4. Леонов В.В. Проектирование и расчёт характеристик концентраторов высокотемпературных солнечных энергоустановок космического назначения //Актуальные проблемы Российской космонавтики: Труды XXXIII Академических чтений по космонавтике. М., 2009. С. 420-422.

5. Леонов В.В. Статистическое моделирование процесса теплообмена излучением в системе концентратор-приёмник солнечной энергии // Актуальные проблемы фундаментальных наук: Сборник трудов. М., 2009. С. 43-45.

6. Леонов В.В. Моделирование процесса теплообмена излучением в системе концентратор-приёмник солнечной энергии // Аэрокосмические технологии: Научные материалы II Международной научно-технической конференции, посвящённой 95-летию со дня рождения академика В.Н. Челомея. М., 2009. С. 138-139.

7. Леонов В.В. Статистическое моделирование процесса теплообмена излучением в системе концентратор-приёмник солнечной энергии // Проблемы газодинамики и тепломассообмена в аэрокосмических технологиях: Труды XVII Школы-семинара молодых учёных и специалистов под руководством академика РАН А.И. Леонтьева. М., 2009. Т. 2. С. 221-224.

8. Леонов В.В. Моделирование процесса теплообмена излучением в системе концентратор-приёмник солнечной энергии //Тепловые процессы в технике. 2009. Т. 1, № 8. С. 340-342.

9. Кувыркин Г.Н., Леонов В.В. Моделирование процесса теплообмена излучением в зеркальных концентрирующих системах // Научные исследования в области транспортных, авиационных и космических систем "АКТ-2009": Труды X Всероссийской научно-технической конференции и школы молодых учёных, аспирантов и студентов. Воронеж, 2009. С. 86-91.

Ю.Леонов В.В. Моделирование процесса теплообмена излучением в системе концентратор-приёмник солнечной энергии высокотемпературной солнечной энергоустановки космического назначения // Актуальные проблемы Российской космонавтики: Труды XXXIV Академических чтений по космонавтике. М., 2010. С. 415-417.

П.Леонов В.В. Расчёт характеристик радиационного теплообмена в системе концентратор-приёмник солнечной энергии // Актуальные проблемы фундаментальных наук: Сборник трудов. М., 2010. С. 24-26.

12. Леонов В.В. Математическое моделирование радиационного теплообмена в системе концентратор-приёмник солнечной энергоустановки космического назначения // Математика и математическое моделирование: Сборник трудов IV Всероссийской молодежной научно-инновационной школы. Саров, 2010. С. 81-83.

13.Leonov V.V. Mathematical Modeling of a Radiative Heat Transfer Process in the System of the Solar Energy Concentrator-receiver for Space High-temperature Solar Power Plant // Proc. 14th Int. Heat Transfer Conf. Washington DC, 2010 (Электронное издание).

14. Леонов В.В. Проектирование зеркальных концентрирующих систем для солнечных энергоустановок космического назначения // Научные исследования в области транспортных, авиационных и космических систем: Труды XI Всероссийской научно-технической конференции и школы молодых учёных, аспирантов и студентов. Воронеж, 2010. С. 175-180.

15. Леонов В.В. Моделирование радиационного теплообмена в системе концентратор-приёмник солнечной энергии // Труды V Российской национальной конференции по теплообмену. М., 2010. Т. 6. С. 231-234.

16. Леонов В.В. Определение радиационных характеристик системы концентратор-приёмник солнечной энергии высокотемпературной солнечной энергоустановки // Актуальные проблемы Российской космонавтики: Труды XXXV Академических чтений по космонавтике. М., 2011. С. 98-99.

17. Кувыркин Г.Н., Зарубин B.C., Леонов В.В. Математическая модель системы концентратор-приёмник высокотемпературной солнечной энергоустановки космического назначения //Вестник МГТУ. Серия «Машиностроение». 2011. № 1. С. 82-91.

18. Леонов В.В. Моделирование радиационного теплообмена в системе концентратор-приёмник солнечной энергии //Тепловые процессы в технике. 2011. Т. 1,№ 5. С. 228-233.

Подписано в печать 13.09.2011г.

Усл.п.л. - 1.0 Заказ №05956 Тираж: ЮОэкз.

Копицентр «ЧЕРТЕЖ.ру» ИНН 7701723201 107023, Москва, ул.Б.Семеновская 11, стр.12 (495) 542-7389 www.chertez.ru

Оглавление автор диссертации — кандидата технических наук Леонов, Виктор Витальевич

Список сокращений.

Список обозначений.

Введение.

Глава 1. Обзор литературы. Анализ существующих методик.

1.1. Тенденции в проектировании и создании солнечных энергоустановок.

1.2. Основные принципы построения высокотемпературных солнечных энергоустановок космического назначения.

1.3. Зеркальные концентрирующие системы.

1.4. Основные характеристики зеркальных концентрирующих систем.

1.5. Характеристики отражённого излучения. Выбор основных расчётных величин.

1.5.1. Рассеяние излучения в среде.

1.5.2. Рассеяние излучения на границе двух сред.

1.6. Основные методики определения распределения облучённости в системе концентратор-приёмник.

1.7. Выводы по 1-й главе.

Глава 2. Разработка математических моделей.

2.1. Статистическое моделирование и методы Монте-Карло.

2.1.1. Генерирование случайных чисел.

2.2. Геометрическая модель.

2.2.1. Параметры шероховатости.

2.2.2. Учёт влияния шероховатости поверхности в задачах теплообмена.

2.3. Модель излучения.

2.3.1. Генерация пучков излучения.

2.3.2. Распространение пучков излучения.

2.3.3. Взаимодействие пучков с поверхностью.

2.4. Моделирование полей температуры.

2.5. Общий алгоритм моделирования радиационного теплообмена.

2.6. Выводы по 2-й главе.

Глава 3. Комплекс прикладных программ Tracer.

3.1. Общая функциональная структура КПП Tracer.

3.1.1. Программы подготовки данных.

3.1.2. Программы основного вычислительного блока.

3.1.3. Программы отображения результатов расчета.

3.2. Верификация КПП.

3.2.1. Угловые коэффициенты.

3.2.2. Индикатрисы отражения.

3.2.3. Радиационный теплообмен.

3.2.4. Кондуктивный теплообмен.

3.3. Выводы по 3-й главе.

Глава 4. Исследование зеркальных концентрирующих систем.

4.1. Индикатриса отражения.

4.1.1. Малые неровности (шероховатость).

4.1.2. Глубокие полости.

4.2. Исследование радиационных характеристик СКП.

4.2.1. Плотность потока сконцентрированного излучения.

4.2.2. Поле температур.

4.2.3. Геометрический КПД СКП.

4.2.4. Энергетический КПД СКП.

4.2.5. Доконцентраторы.

Заключение диссертация на тему "Разработка математической модели концентратора высокотемпературной солнечной энергоустановки космического назначения"

2.6. Выводы по 2-й главе

Результаты» расчета по разработанным и использованным известным моделям» и их достоверность непосредственно зависят от характеристик используемого ГСЧ. Поэтому основным моментом в реализации? математических моделей, построенных на принципах статистического моделирования с использованием методов Монте-Карло, является генерация последовательностей случайных чисел. В результате анализа различных ГСЧ для генерации последовательности случайных чисел, равномерно распределенных в интервале от 0 до 1, решено использовать вихрь Меерсона в виде алгоритма МТ19937.

На основе полученных случайных чисел или их. модификации» (для получения распределения^ отличного от равномерного) и построена модель-распространения радиационных тепловых потоков. При этом процесс распространения каждого конкретного пучка излучения был рассмотрен как случайная марковская цепочка. Каждая цепочка описывает историю движения: конкретного пучка излучения от момента его зарождения (генерации) до тех пор, пока пучок не будет поглощён или не покинет исследуемую систему.

Развитие метода гауссовой нормали применительно к пучковой модели излучения позволило отказаться от диффузной модели отражения, излучения от поверхности, рассмотренной- применительно к методу Монте-Карло в работах [96, 102]. Это дало возможность реалистичней описать процесс отражения« и привело к большей гибкости в описании свойств отражающих (рабочих) поверхностей ЗКС.

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

Глава 3. Комплекс прикладных программ Tracer

Для реализации описанных в предыдущей главе математических моделей и алгоритмов был разработан комплекс прикладных программ (КПП) Tracer. Приведём общее описание этого КПП и верификацию его работы на тестовых задачах.

3.1. Общая функциональная структура КПП Tracer

Препроцессор

Комплекс программ подготовки исходных данных

Комплекс программ подготовки графической модели

Процессор

Комплекс программ по определению распределения тепловых потоков

Комплекс программ по определению полей температур а е и и О

SJ о S. в н и о в

Рис. 3.1.

Общая структурная схема КПП Tracer

На рис. 3.1 представлена общая структурная схема КПП Tracer. Комплекс состоит из трех больших групп программ: препроцессора, процессора и постпроцессора. Препроцессор подготавливает графическую модель и всю совокупность данных для решения задачи радиационного теплообмена в заданной системе. Блок процессора состоит из двух частей. Первую часть составляют программы для решения задач распространения радиационных тепловых потоков в исследуемой системе, а вторую - программы, предназначенные для решения задач определения полей температур. Использование программ постпроцессора происходит каждый раз, когда возникает необходимость представления результатов обработки данных как на стадии их подготовки, так и на этапах численного решения задач. Все программы* комплекса написаны на языке программирования Delphi (Изначально Бе1рЫ являлся средой объектного программирования, основанной на языке Object Pascal, который в свою очередь основывается на процедурном языке программирования Pascal фирмы Borland, но после перехода прав на среду Delphi от фирмы Borland к CodeGear в 2007 году, за языком» окончательно закрепилось название Delphi.). Общий объем программного кода составляет около 18000 строк.

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

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

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

3.1.1. Программы подготовки данных

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

Построение геометрической модели

Основным геометрическим примитивом, используемым в данном комплексе, является треугольный элемент (см. п. 2.2). Группируя связкютаких элементов; получаем поверхности произвольной формы (плоскости; цилиндры-, сферы, параболоиды и др.). При этом имеется возможность как поэлементного построения; задавая координаты вершин каждого треугольного элемента, так и. воспользоваться библиотекой геометрических объектов dll graf.

Задание материалов и радиационных свойств поверхностей

Задание материалов может производиться как непосредственно при создании геометрической модели; так и впоследствии при её редактировании. Термомеханические свойства задаваемых материалов хранятся в подключаемой базе данных. Так же есть возможность задавать/редактировать материалы на стадии подготовки КЭ модели.

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

Построение конечноэлементной модели и задание граничных условий

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

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

В качестве примера на рис. 3.2 приведена геометрическая модель, покрытая сеткой КЭ, для системы Солнце-концентратор-приёмник, построенная при помощи препроцессора КПП Tracer.

Рис. 3.2.

Геометрическая модель системы Солнце-концентратор-приёмник, покрытая сеткой треугольных КЭ

3.1:2. Программы основного вычислительного блока

Основным в структуре КПП Tracer является программный блок, входящий» в состав процессора. Программы этого блока реализуют все математические модели, изложенные в предшествующих главах. Ядро блока-построено по модульному принципу И" имеет две главные составляющие -модули Trace и Therm. Первый из них предназначен для моделирования распространения радиационных тепловых, потоков, а второй для решения задач теплопроводности и вычисления полей температур.

Работу ядра процессора и, его модулей обеспечивают две динамические библиотеки: dlltrace и dllfem. Первая из них содержит все процедуры и функции необходимые для реализации математических моделей радиационного теплообмена: генерации излучения, его распространения, взаимодействия с геометрической^ моделью: Вторая обеспечивает реализацию МКЭ: формирование рабочих матриц МКЭ, решение систем линейных алгебраических уравнений (СЛАУ), коррекцию параметров в соответствии с найденными узловыми значениями искомых функций, а также выполнение всех необходимых промежуточных процедур.

Основные процедуры при моделировании радиационного теплообмена

Последовательность, работы модуля. Trace, реализующего процесс моделирования радиационного теплообмена, представлена на укрупнённой структурной схеме рис: 3.3, и состоит из нескольких этапов.

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

Второй представляет собой макроцикл, реализующий непосредственно алгоритм моделирования радиационного теплообмена, приведённый в п. 2.5. Выполнение условия Пу > N (см. рис. 3.3) означает, что текущий номер nj пучка излучения превысил общее количество N пучков и завершение макроцикла. Макроцикл состоит из трех шагов: первый-это генерация пучка излучения и прослеживание его траектории, второй — моделирование свойств поверхности в точке падения пучка и третий-моделирование взаимодействия пучка и поверхности. Второй и третий этапы повторяются до тех пор, пока пучок либо не будет поглощен (Ra < от), либо не покинет исследуемую систему (flagi = false).

Рис. 3.3.

Общая структурная схема модуля Trace

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

Основные процедуры при решении задач теплопроводности

На рис. 3.4 представлена укрупнённая структурная схема модуля Therm, реализующего алгоритмы численного решения задач теплопроводности. Процесс численного решения этих задач состоит из нескольких этапов.

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

Рис. 3.4.

Общая структурная схема КПП Therm

Вторым этапом является формирование глобальной матрицы СЛАУ и матрицы-столбца тепловых нагрузок. Их формирование осуществляется в процессе макроцикла интегрирования по КЭ: Выполнение условия nj > N (см. рис. 3.4) означает, что текущий номер п7- КЭ превысил общее количество N элементов в сетке КЭ и завершение макроцикла. Макроцикл состоит из трех шагов: первый - это подготовка данных для фиксированного КЭ, второй -собственно численное интегрирование, позволяющее сформировать локальные матрицы. КЭ и его граней, если эти грани соответствуют границам и, наконец, третий - объединение локальных матриц в глобальные. Особенностью использования квазидвумерных конечных элементов является то, что все элементы, как правило, являются граничными:

Решение нелинейных задач теплопроводности связано с процессом последовательных приближений [23], что составляет существо третьего этапа. Основным методом решения СЛАУ, сформированной на втором этапе, является метод Зейделя [7, 96]. При этом подготовку данных и численное интегрирование приходится проводить на каждой итерации с учетом результатов, полученных на предыдущей итерации. Ясно, что в случае линейной» стационарной задачи сформированную СЛАУ достаточно решить один раз, а в случае линейной нестационарной задачи-один раз на каждом шаге по времени.

На «четвертом* этапе происходит анализ сходимости, численного решения, нелинейной задачи и/или проверка условия достижения заданного предельного времени (в случае нестационарной задачи). Если сходимость не достигнута и/или рассматриваемый момент времени не является предельным, то следует переход к описанным выше процедурам формирования, глобальных матриц СЛАУ.

3.1.3. Программы отображения результатов расчета

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

Работу постпроцессора обеспечивают две библиотеки: и

111роз1ргос. Первая отвечает за графические построения, а вторая за обработку результатов моделирования.

3.2. Верификация КПП

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

3.2.1. Угловые коэффициенты

Одной из основных математических трудностей в исследовании, теплообмена излучением между поверхностями является учёт геометрических соотношений, определяющих под каким углом поверхности видны друг другу; Для учёта этих величин вводятся угловые коэффициенты [9].

Рассмотрим угловой коэффициент для случая, когда излучение, испускаемое изотермической поверхностью попадает на абсолютно чёрную поверхность 52 (рис. 3.5). По определению [9, 102], угловой коэффициент есть доля энергии испускаемого излучения, падающая на Б2 стТ^Б1 где аТ*^-поток излучения, испускаемый поверхностью

Рис. 3.5.

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

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

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

Сфера 1

Рис. 3.6. Расчётная схема

Всё излучение, испускаемое с поверхности внутренней сферы попадает на поверхность тогда

1-2 = 1 •

Используя соотношение взаимности для угловых коэффициентов [102], получаем р ^1^1-2 £1 1 $2 $2

Часть энергии, испускаемой поверхностью 52, попадает на эту же поверхность, тогда, учитывая, что для замкнутой системы, состоящей из N поверхностей N %

IV 1

7=1 получаем

Тогда

Р2-г + Р2-2 ~ 1 ■

2 —

Р2-2 =-с-■ 2

Проведём исследование подобной системы с помощью разрабатываемого программного комплекса. В качестве примера возьмём сферы с радиусами = 1 м, Я2 = 3 м. Тогда расчётное значение геометрических коэффициентов

1-2 = 1, = \ = 0.1111, 8

Р2-2 = д = 0.8889.

Значения, получаемые в результате математического эксперимента, показаны для Р2г на рис. 3.7 и для Р2-г и Р2-2 совместно на рис. 3.8 в зависимо от числа рассмотренных историй.

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

Рис. 3.9. Расчётная схема

Для этой системы угловой коэффициент можно определить из соотношения [102] Лт [п■ ага§ (¿) + Я- аг*§ - V"2 + ^ ал* (-===) + 1 -1п 4 где

7Г ■ IV

1 + Ж2)(1 + Я2)] [ \У2(1 + IV2 + Я2) я

IV2

1 + Ж2 + Я2 11 (1 + И^2)(Ж2 + Я2)

Я2(1 + Ж2+Я2) н2м

1 + Я2)(Ж2 + Я2) I я=т,

IV — Г

На рис. 3.10 показана зависимость углового коэффициента от геометрических размеров элементов системы.

Диффузный угловой коэффициент Р1„2 можно определить из соотношения [22]

Р±2 |

К^Ч'-1 где г и Н - соответственно радиус и высота параболической оболочки (см. рис. 3.11).

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

Рг-2 = 1 — Р\-г

В таблице 3 приведено сравнение результатов моделирования при различном числе историй с результатом аналитического решения (для Н/г = 0.5, = 0.82).

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

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

2. Разработан комплекс прикладных программ «Tracer», реализующий разработанные математическую модель и алгоритмы.

3. Проведена верификация математической модели и комплекса прикладных программ, подтвердившая их работоспособность.

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

5. Проведён анализ влияние геометрических параметров и различных дефектов поверхности на радиационные характеристики зеркальных концентрирующих систем на примере параболического концентратора.

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

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

Библиография Леонов, Виктор Витальевич, диссертация по теме Математическое моделирование, численные методы и комплексы программ

1. Агабабов С.Г. Влияние шероховатости поверхности твёрдого тела на его радиационные свойства и методы их экспериментального определения // Теплофизика экспериментальных температур. 1968. Т.8, № 11. С. 78-88.

2. Адрианов В.Н. Основы радиационного и сложного теплообмена. М: Энергия, 1972. 464 с.

3. Алфёров Ж.И., Андреев В.М., Румянцев В.Д. Тенденции и перспективы развития солнечной фотоэнергетики // Физика и техника полупроводников. 2004. Т. 38, № 8. С. 937-948.

4. Апариси Р.Р. Концентрация солнечной энергии в гелиотехнических сооружениях. Автореф. канд. дисс. М., 1955. 16 с.

5. Апариси Р.Р. Экспериментальная установка для получения высоких температур // Использование солнечной энергии. М.: Изд-во АН СССР, 1961. № 1.С. 151-162.

6. БаумИ.В: Формирование поля облучённости приёмника в "точных" и "неточных" гелиоконцентраторах // Солнечные энергетические установки. М: ЭНИН, 1974. С. 213-234.

7. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: БИНОМ, 2006. 636 с.

8. БирбекР.С. Влияние шероховатости поверхности на общую полусферическую и зеркальную отражательную способность металлических поверхностей // Труды амер. о-ва инж.-мех. сер. С. 1964. Т. 2, № 74. С. 24-36.

9. Блох А.Г. Теплообмен излучением: Справочник. М: Энергоатомиздат, 1991.432 с.

10. Ю.Бугер П. Оптический трактат о градации света. М.: Изд-во АН СССР, 1950. 484 с.

11. П.Гершун A.A. Световое поле. Л.: ОНТИ, 1936. 178 с.

12. Грилихес В.А. Солнечные высокотемпературные источники тепла для космических аппаратов. М: Машиностроение, 1975. 248 с.

13. Грилихес В.А. Солнечные космические электростанции. Л.: Наука, 1986. 182 с.

14. Грилихес В.А., Захидов P.A. К выводу формулы распределения облучённости в фокальной плоскости параболических концентраторов //Гелиотехника. 1971. Т. 4. С. 9-13.

15. ГрилихесВ.А. Алгоритм статистического решения задачи о распределении лучистых потоков в приёмниках солнечных установок с параболоидными концентраторами // Гелиотехника. 1966. Т. 4. С. 25- 34.

16. Грилихес В.А. Концентратор солнечной энергии: Сборник статей. JI: Энергия, 1972.248 с.

17. Грилихес В.А. Метод расчёта распределения лучистых потоков в приёмниках солнечных установок//Гелиотехника. 1966. Т. 1. С. 3-11.

18. ГэтландК. Космическая техника. М.: Мир, 1986. 295 с.

19. Дрессер Д. Основы проектирования солнечных энергетических устройств // Использование солнечной энергии при космических исследованиях. М.: Мир, 1964. С. 137-204.

20. Жозе П. Распределение плотности потока энергии в фокальном изображении солнечной печи // Солнечные высокотемпературные печи. М.: ИЛ, 1960. С. 229-238.

21. Зарубин B.C. Кувыркин Г.Н. Математические модели механики и электродинамики сплошной среды. М.: Изд-во МГТУ им. Н.Э. Баумана, 2008. 512 с.22.3арубин B.C. Температурные поля в конструкции летательных аппаратов. М.: Машиностроение, 1978.184 с.

22. Изменение морфологии поверхности металлов^ при сверхзвуковом соударении» О.Н. Никитушкина, Л.И. Иванов, С.А. Бедняков, Л.С. Новиков // ФХОМ, 2001. № 1. С. 48-51.

23. Использование солнечной энергии при космических исследованиях: сб. статей. М.: Мир, 1964. 415 с.

24. Кабан Ф., Вен Л.Ф: Распределение плотности, потока энергии в, фокальном изображении солнечной печи // Солнечные высокотемпературные печи. М.: ИЛ, 1960. С. 239-245.

25. Канатников А.Н., Кршценко А.П, Четвериков В.Н. Дифференциальное исчисление функций многих переменных. / Сер. Математика в техническом университете; Вып. V. М.: Изд-во МГТУ им. Н.Э. Баумана, 2000. Т. 5.456 с.

26. КарякинН.А. Световые приборы. М.: Высшая школа, 1975. 335 с.

27. Красина Е.А., Невежин O.A., Рубанович И.М. Оценка точности, отражающей поверхности параболического концентратора при изменении углового размера источника излучения // Гелиотехника. 1974. Т. 3. С. 31-36.

28. Кудрин И.О. Солнечные высокотемпературные космические энергодвигательные установки. М: Машиностроение, 1987. 248 с.

29. Лыков A.B. Теория теплопроводности. М.: Высшая школа, 1967. 600 с.

30. Ляшков В.И., Кузьмин С.Н. Нетрадициошше и возобновляемые источники энергии / Учебное пособие. Тамбов: Изд-во ТГТУ, 2003. 96 с:38;Майзель< C.Oi Сборник; трудов« по* светотехнике. Mi:. Изд-во АЕй СССР* 19381240 с.

31. Мак-Келланд Д. Концентраторы- солнечного излучения! для высокотемпературных энергетических установок космических летательных аппаратов^ // Энергетические установки; для* космических аппаратов; М : Мир, 1964; С. 95-110:

32. Международная космическая станция // Центр» подготовки космонавтов им:Гагарина: сайт.- URb: http://www.gctc.ru/iss/index.litml1 (дата обращения 12.04:2011);

33. Мешков В.В. Основы светотехники. М.: Госэнергоиздат, 1957. 352 с.42 .Михайлов Г. А., Войтишек А.В. Численное статистическое; моделирование: Методы Монте-Карло. М.: Академия, 2006; 368 с.

34. Модель космоса: Воздействие; космической^ среды; на? материалы; и оборудование космических аппаратов! М.: КДУ, 2007. Т. 2.1144 с.

35. Модель космоса: Физические условия в космическом пространстве: М.: КДУ, 2007. Т. 1. 1000 с.

36. Основы идентификации и проектирования?тепловых процессов и систем ОМ1 Алифанов- Е.Н: Вабишевич, В:В: Михайлов» Hi др:.| / Учебное пособие: Mi: Логос, 20011 400 с.46.0цисикМ.Н. Сложный теплообмен / перевод с англ. М: Мир, 1976. 616 с.

37. ПланкМ: Теория теплового излучения: Mi: Гостехиздат, 1938:204'с:.48;Поляк Г.Л. Исследование теплообмена1 излучением между диффузными поверхностями! //Журнал Технической« Физики. 1935. Т. 1, № 5. С. 555-590:

38. Рубанович И.М. Теоретическое исследование влияния точности? наведения- параболического концентратора на Солнце на работу гелиоустановки //Гелиотехника: 1966; Т. 7. С. 18-30:

39. Семенов Н.Н. Об энергетике будущего // Наука и общество: 1973: С. 109-130^

40. Симон А. Расчёт концентрации солнечной энергии в фокальном изображении? параболоидного отражателя; // Солнечные высокотемпературные:печи. М:: ИЛ, 1960. С. 252-263.

41. Скребушевский Б.С. Космические энергетические установки с преобразованием солнечной энергии. М: Машиностроение, 1992. 224 с.

42. Слюсарев Г.Г. Расчёт оптических систем. М.: Машиностроение, 1975. 639 с.

43. Солнечная? энергетика: В.И. Виссарионов; Г.В: Дерюгина, В.А. Кузнецова, Н.К. Малинин М.: Издательский дом МЭИ, 2011. 276 с.

44. Справочник технолога-машиностроителя: М.: Машиностроение, 1986. Т. 1. 656 с.

45. Спэрроу Э.М., Сэсс Р.Д. Теплообмен излучением. Л.: Энергия, 1971. 294 с.

46. Структура микрократеров на поверхности металлических образцов, экспонировавшихся в открытом космосе О.Н: Никитушкина, Л.И. Иванов, С.А. Бедняков и др. // ФХОМ; 2002. № 2. С. 21-25.

47. Твердович Э.В., Мадаев В.В. Испытание параболических концентраторов наоберографе системы Леонова//Гелиотехника. 1974. Т. 3. С. 28-33.

48. Теория теплообмена. Сборник рекомендуемых терминов. М: Наука, 197Г. 80 с.

49. Тепляков Д.И. Влияние центрального затенения на энергетические характеристики параболоидных зеркал // Гелиотехника; 1967. Т. 6. С. 18-24.

50. Тепляков Д.И. Исследование структуры поля излучения в гелиоустановках с отражающими концентраторами // ИФЖ. 1958. Т. 4. С. 31-39.

51. Тепляков Д.И. Особенности переноса и распределения энергии в высокотемпературных солнечных,установках // Гелиотехника. 1967. Т. 3. С. 13-19.

52. Тепляков Д.И., Аннаев А. Исследование оптимальной геометрии приёмников для точных зеркальных гелиоустановок. М.: Изд-во АН СССР, 1961. С. 21-30.

53. Титов В.М., Фалеенко Ю.И. Сквозное пробивание при метеоритном ударе //Космические исследования^ 1972. Т. 10. №4. С. 589-595.

54. Торренс К.С. Двухпараметрическая отражательная способность непроводника электричества как функция длины волны и шероховатости поверхности // Труды амер. о-ва инж.-мех. сер. С. 1965. Т. 2. № 145. С. 56-62.

55. Торренс К.С. Незеркальные пики в пространственном распределении отражённого теплового излучения // Труды амер. о-ва инж.-мех. сер. С. 1966. Т. 2. №81. С. 12-28.

56. Фаворский О.Н., Фишгойт В.В., Литовский Е.И. Основы теории космических электрореактивных двигательных установок. М.: Высшая школа, 1970. 488 с.

57. Физические величины: Справочник. М.: Энергоатомиздат, 1991. 1232 с.

58. Хююоо Мин Проблемы проектирования солнечных печей // Солнечные высокотемпературные печи. М.: ИЛ, 1960. С. 326-339.

59. Центробежные бескаркасные крупногабаритные космические конструкции. Г.Г. Райкунов, В.А. Комков, В.М. Мельников, Б.Н. Харлов М.: ФИЗМАТЛИТ, 2009. 448 с.

60. Чиколев В:Н. Осветительная способность прожекторов электрического света. СПБ, 1892. 224 с.

61. Шмальц Г. Качество поверхности. М.: Машгиз, 1941. 333 с.

62. Abbety G. The Sun. New York: MacMillan, 1957. 336 p.

63. Bliss R.W. Holloman Air Development Center Technical Memorandum. Alamogordo (NM): HDW-TW, 1958. V. 1. 24 p.

64. Corlett R.C. Direct numerical simulation of thermal radiation in vacuum // J. Heat Transfer. 1966. V. 88. P. 376-382.

65. DaviesH. The reflection of electromagnetic waves from a rough surface // In: Proc. of the Institution of Electrical Engineers. 1954. P. 209-214.

66. Drolshagen G. Meteoroid. Debris impact analysis application to LDEF, EURECA and COLUMBUS // In: Proc. of the 1st Europ. Conf. on Spacc Debris, Darmstadt, Germany. 1993. P. 515-522.

67. Eco-towns. Living a greener future. London: Department for Communities and Local Government, 2008. 57 p.

68. Energy technology perspectives 2010. Paris: International Energy Agency, 2010. 20 p.

69. Glaser P.E. Power from the Sun: ifs Future // Science. 1968. V. 168. P. 857-861.

70. Hopf E. Mathematical Problems of Radiative Equilibrium. London: Cambridge University Press, 1934.110 p.

71. Howell J.R., Perlmutter M. Monte Carlo solution of radiant heat transfer through radiant media between gray walls // Heat Transfer. 1964. P. 116-122.

72. Howell J.R., Strite M.K., Renkell H. Heat transfer analasis of rocket nozzles using veiy hight temperature propellants // AALA J. 1965. V. 3. P. 669-673.95 .Lide D.R. CRC Handbook of Chemistry and Physics. Boca Raton, FL: CRC Press, 2005.2660 p.

73. Mahan J.R. Radiation heat transfer: a statistical approach. New York: Wiley & Sons, 2002. 500 p.

74. Matsumoto M., NishimuraT. Mersenne twister: A 623-dimensionally equidistributed uniform pseudorandom number generator //ACM Trans, on Modeling and Computer Simulations. 1998. V. 8, № 1. P. 3-30.

75. Palik E.D. Handbook of Optical Constants of Solids. Academic Press, 1998. V. 3. 999 p.

76. Perlmutter M., Siegel R. Effect of specularly reflecting gray surface on thermal radiation through a tube and from its heated wall // J. Heat Transfer. 1963. P. 55-62.

77. Polgar L.G., Howell J.R. Directional thermal-radiative properties of conical cavities. NASA, 1965.158 p.

78. SafwatH.H. Effect of Surface Roughness on Specular and Diffuse Reflectance Components of Carbon Steel in visible Wevelength Range // ASME. 1969. Paper 60-WA/HT-42. P. 45-58.

79. Siegel R., Howell J.R. Thermal radiation heat transfer. New York: Taylor & Francis, 2002. 857 p.

80. Soules J. Design and Fabrication of a Dielectric Total Internal Reflecting Solar Concentrator and Associated Flux Extractor for Extreme High Temperature (2500K) Applications. San Diego: NASA Contractor Report, 1997. 17 p.

81. Space-Based Solar Power as an Opportunity for Strategic Security // Report to the Director, National Security Space Office. Arlington, VA, 2007. 75 p.

82. W.E.F. Smithsonian Physical Tables. Washington: Smithsonian Institution, 1954. V. 9. 722 p.

83. Weiner M.M., Tindall J.W., Candell M. Radiative interchange factors by Monte Carlo // ASME Paper 65-WA/HT-51,1965. P. 34-49.