CC BY
13ISSN 1814-1196 Научный вестник НГТУ том 71, № 2, 2018, с. 23-38
http://journals.nstu.ru/vestnik Science Bulletin of the NSTU Vol. 71, No. 2, 2018, pp. 23-38
ИНФОРМАТИКА, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И УПРАВЛЕНИЕ
INFORMATICS, COMPPUTER ENGINEERIN AND CONTROL
УДК 629.7.042.2.001.24:622.998 DOI: 10.17212/1814-1196-2018-2-23-38
Математическое моделирование теплового состояния негерметизированного отсека беспилотного летательного аппарата
С.А. ГУСЕВ12, ВН. НИКОЛАЕВ3
1 630090, пр. академика Лаврентьева, 6, Институт вычислительной математики и математической геофизики Сибирского отделения Российской академии наук, доктор физико-математических наук, старший научный сотрудник. E-mail: sag@ osmf.sscc.ru
630073, РФ, г. Новосибирск, пр. Карла Маркса, 20, Новосибирский государственный технический университет, доктор физико-математических наук, профессор. E-mail: sag@osmf.sscc.ru
3 630051, РФ, г. Новосибирск ул. Ползунова, 21, ФГУП «Сибирский научно-исследовательский институт авиации имени С.А. Чаплыгина», доктор технических наук, начальник отдела. E-mail: nikvla50@mail.ru
Разработана математическая модель теплового состояния отсека летательного аппарата, описывающая теплообмен обшивки с сотовыми конструкциями из углепластика, процесса теплопередачи бортового оборудования и воздушной среды. Рассматриваемый процесс переноса тепла в неоднородной среде описывается краевой задачей для уравнения теплопроводности с граничными условиями третьего рода.
Для решения прямой задачи теплового состояния сотовой конструкции фюзеляжа использован метод Монте-Карло на основе вероятностного представления решения в виде математического ожидания функционала от диффузионного процесса.
При этом коэффициенты температуропроводности материалов, из которых состоит панель, являются постоянными величинами, а процесс передачи тепла осуществляется только посредством теплопроводности. Поэтому в вероятностном представлении решения краевой задачи соответствующий случайный процесс с точностью до постоянного множителя совпадает с винеровским процессом без сноса в подобластях, в которых среда однородна. В работе предлагается для моделирования винеровского процесса внутри ячеек сотовой панели на некотором заданном удалении от каркаса использовать метод случайного блуждания по движущимся сферам, а для расчетов по каркасу панели и в его окрестности - метод Эйлера. Такой подход дает значительное ускорение счета по сравнению с использованием для моделирования траекторий только одного метода Эйлера.
Обратная задачи теплообмена сотовой конструкции решена путем минимизации функции взвешенной суммы квадратов невязок с помощью итерационного стохастического квазиградиентного алгоритма.
Статья получена 30 марта 2018 г.
Разработанная математическая модель теплового состояния отсека была применена для оптимизации температуры и расхода воздуха системы обеспечения теплового режима приборного продуваемого теплоизолированного отсека летательного аппарата.
Ключевые слова: математическая модель, тепловое состояние, сотовая конструкция, гетерогенные структуры, численное решение, параболическая краевая задача, метод Монте-Карло, разрывные коэффициенты
ВВЕДЕНИЕ
При проектировании отсеков летательного аппарата с заданной надежностью его работы необходимо исследование процессов теплообмена, происходящих во время полета, между внешней средой и находящимся внутри оборудованием. Для этой цели наиболее удобным и признанным средством является метод математического моделирования.
Тепловое состояние приборных отсеков летательного аппарата формируется внешними и внутренними факторами. К внешним факторам относятся теплообмен между наружной поверхностью обшивки и тепловым и аэродинамическим пограничными слоями, радиационный теплообмен наружной поверхности обшивки между Солнцем, Землей и другими поверхностями. Внутри отсека тепловая энергия выделяется бортовым оборудованием и отводится или подводится системой обеспечения теплового режима [1].
Модель должна учитывать нестационарный теплообмен теплоизолированной обшивки, перенос воздушным потоком тепловой энергии из одной части отсека в другую, конвективный теплообмен элементов отсеков, и в данной работе она представляет собой систему дифференциальных уравнений в частных производных и обыкновенных дифференциальных уравнений. Для реального отсека количество уравнений модели может быть очень большим и достигать нескольких сотен единиц. Поэтому возникает необходимость в эффективных численных методах решения прямой и обратной задач по исследованию теплообмена элементов отсеков, оценивания погрешности параметрической идентификации.
1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ТЕПЛОВОГО СОСТОЯНИЯ ОТСЕКОВ ЛЕТАТЕЛЬНОГО АППАРАТА
Обшивка рассматриваемого типа летательных аппаратов представляет собой сложную конструкцию, которая включает в качестве основного элемента теплозащиты сотовые панели с каркасом из углепластика, заполненным воздухом в качестве основного наполнителя. Процесс теплообмена в сотовых панелях описывается краевой задачей для уравнения теплопроводности с разрывными коэффициентами. Для решения этой краевой задачи используется метод Монте-Карло на основе стохастических дифференциальных уравнений в сочетании с методом блуждания по движущимся сферам [2]. В общем виде процесс теплопередачи в многослойной обшивке описывается следующими уравнениями [1, 3]:
^ (х) Т^ = (х) х)х, 0 < х < I, 0 < X < Хк; (1)
Хcv (x) Fcv Tcv,x acv,out
(t) Fcv (Tcv (t, x) - Te (t)) + Qcv out, x = 0; (2) Хcv (x) Fcv Tcv,x = acv,in (t) Fcv (Tair (t) — Tcv (t,x)) + Qcv,in, x = l; (3)
Tcv (0, x) = T0(x), 0 < x < l, (4)
„ , ч \Ccompo, x e comPo , гДе Ccv (x) = 1
I Cair, x e air ,
л . . \hompo , x e comPo ,
Kv (x) = 1
I Kir, x e air,
т. е. коэффициенты Ccv, Xcv зависят от того, в каком слое рассматривается перенос тепла.
В уравнениях (1) - (4) использованы следующие обозначения: T^ (t, x) - температура сотовой конструкции; Tcvt - первая производная T^
по времени t; Tcvx - первая производная T^ по x; Tcv xx - вторая производная Tcv по x; Ccv (x) - объемная теплоемкость сотовой конструкции обшивки, определяемая теплоемкостью композита Се0тр0 и теплоемкостью воздуха Сс^г ; (l) - теплопроводность сотовой конструкции, определяемая теплопроводностью композита Xc0mp0 и теплопроводностью воздуха *kair; acvout - коэффициент теплоотдачи наружной поверхности конструкции; acvin - коэффициент теплоотдачи внутренней поверхности конструкции; Fcv - площадь конструкции при наружном и внутреннем теплообмене; Qcvcmt - тепловая энергия внешних источников; Qcvin - тепловая энергия внутренних источников; Te - температура восстановления; t - время; Tair -температура воздушной среды в отсеке или в части отсека; Tj - температура
j-го элемента отсека; l - толщина сотовой конструкции.
Процесс теплопередачи бортового оборудования представим в виде обыкновенного дифференциального уравнения, описывающего конвективно-радиационный теплообмен с окружающими конструкциями
Tm,t = aair ,m (t) Fair ,m ^ Cm (Tair (t) — Tm ) +
+ S S j,m ICm Tj (t) Tms — c0^m Fm ICmTm + Qm ^ Cm, (5)
ml /
где Tm - температура m-го бортового оборудования; Tmt - первая производная Tm по t; aair m - коэффициент теплоотдачи m-го бортового оборудования; Fair m - площадь m-го бортового оборудования при конвективном теплообмене; Cm - теплоемкость i-го бортового оборудования; gj m - коэффи-
циент радиационного теплообмена системы «у-й элемент отсека - т-й блок бортового оборудования»; 8т - коэффициент черноты излучения т-го блока; Qm - энергия тепловыделения или теплопоглощения т-м бортовым оборудованием от системы кондиционирования и преобразованная из электрической энергии.
Процесс теплообмена воздушной среды в нетеплоизолированном негерметичном непродуваемом отсеке представим в виде обыкновенного дифференциального уравнения, описывающего конвективный теплообмен между внутренней поверхностью обшивки и бортовым оборудованием т:
Ттг,ырг/ _ атг,т () ^су,т,ирг ^ Сагг (Тсу,1п,ырг (1су,ырг, ^) _ Тагг,ырг ) ^
Уравнение теплообмена воздушной среды в негерметичном продуваемом теплоизолированном или нетеплоизолированном приборном отсеке представим в виде обыкновенного дифференциального уравнения, описывающего конвективный теплообмен внутренней поверхности обшивки, элементов отсека и перенос энтальпии из одной части отсека в другую:
где Таг к_, Таг к - температуры воздушного потока соответственно в (к-1)-й и к-й частях отсека; Jajr к - массовая скорость воздушного потока в к-й части отсека; Рк - суммарная площадь воздушных каналов в к-й части отсека; Ср - удельная теплоемкость воздуха; Са^ к - теплоемкость воздуха в к-й части отсека.
Суммирование в уравнении (7) ведется по у-му элементу, входящему в к-ю часть отсека.
Теплоемкость воздуха Са¿г к определяется по выражению
где рагг к - плотность воздуха в к-й части отсека; п - скорость воздуха на входе в отсек; еп{ - площадь воздушных каналов на входе в первую часть отсека; & - интервал дискретизации времени при решении системы дифференциальных уравнений; ¥а^г к - объем воздуха в к-й части отсека.
(6)
(7)
(8)
Коэффициенты теплоотдачи поверхностей асг,,омг, асг;п, аа;г,т в уравнениях (2)-(7) будем вычислять с помощью методик, описанных в работах [3, 4].
Коэффициенты радиационного теплообмена в уравнении (5) определяются методом Монте-Карло [5].
2. ПРИМЕНЕНИЕ МЕТОДА МОНТЕ-КАРЛО ДЛЯ РЕШЕНИЯ ПРЯМОЙ ЗАДАЧИ ТЕПЛОВОГО СОСТОЯНИЯ СОТОВОЙ КОНСТРУКЦИИ ФЮЗЕЛЯЖА
В том случае, когда сотовая панель (рис. 1) в составе многослойной обшивки рассматривается как однородная среда с усредненными коэффициентами объемной теплоемкости и теплопроводности, теплоперенос через обшивку описывается уравнениями (1)-(4). Однако известно, что усредненные теплофизические свойства неоднородных сред могут меняться при изменении направления теплового потока [6]. По этой причине мы также рассматриваем определение теплового состояния сотовой панели как решение трехмерной краевой задачи для уравнения теплопроводности с разрывным коэффициентом температуропроводности. Здесь в силу особенностей применяемого метода делается допущение, что коэффициенты температуропроводности углепластика и воздуха являются величинами постоянными.
Рис. 1. Сотовая конструкция фюзеляжа
Ниже дано описание рассматриваемой краевой задачи. Областью, в которой определяется рассматриваемая краевая задача, является прямоугольный параллелепипед О = (-/1, /1) х (-/2, /2) х (0, /3). При этом О есть объединение двух непересекающихся подмножеств: О = и О2, где О - подобласть, соответствующая каркасу и пластинам, ограничивающим панель, а О2 есть объединение подобластей, соответствующих ячейкам, содержащим воздух. Рассматриваемый процесс теплопередачи происходит на отрезке времени [0, Т] и описывается следующей краевой задачей для уравнения теплопроводности:
¡Тс* д г
з д ■=£ дд-
1=1д X ч
( дТ ^ а (х)
д X;
(9)
где
а( х) =
I асотро, х е О1,
х е Ол х е О2;
(10)
д Тс
д XI
= 0,
С1=-/1
д Тс
д XI
= 0;
С1 =1
(11)
д Тс
д х2
= 0,
С2 =-/2
д Тс
д х2
= 0;
С2 =/2
(12)
х дТсу "А,
д х3
и )(Т^ - " Тагг ,ои и));
х3 =1
(13)
х, дТ~
д х3
асу,1п () (Тсу Та1глп ())-
х3 =1
(14)
В уравнениях (9)-(14) использованы следующие обозначения: асотро, аа1г - коэффициенты температуропроводности композита и воздуха соответственно; асу оШ, асу п - коэффициенты теплообмена поверхности панели и воздушной средой снаружи и внутри бортового оборудования соответственно; Та{г оШ, Та{г п - температура воздушной среды у наружной стороны панели и внутренней соответственно.
В работе [7] доказано существование обобщенных решений краевых задач с разрывными коэффициентами. Причем эти решения могут быть аппроксимированы решениями краевых задач, у которых коэффициенты являются приближениями исходных разрывных коэффициентов. Например, можно получить приближенное решение исходной задачи, решая задачу со сглаженными коэффициентами на основе интегрального усреднения [8]. В данной работе предлагается определять приближенное решение задачи как задачи со сглаженными коэффициентами методом Монте-Карло на основе вероятностного представления решения в виде математического ожидания функционала от диффузионного процесса. Первоначально в работе [9] оценки математического ожидания этого функционала определялись на основе численного решения стохастических дифференциальных уравнений методом Эйлера. Недостатком такого метода является его большая трудоемкость. Значительное ускорение счета удалось получить с помощью предложенного в работе [2] комбинированного метода, в котором расчет траекторий диффузионного процесса в ячейках, заполненных воздухом (область G2), осуществлялся методом блуждания по движущимся сферам, а по каркасу, ограничивающим пластинам (область G|) и в близкой их окрестности - методом Эйлера. Заметим, что применение комбинированного метода возможно только в случае постоянных теплофизических свойств веществ, из которых состоит сотовая панель. Подробное описание комбинированного метода дано в работе [2].
3. АЛГОРИТМ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ТЕПЛОВОГО СОСТОЯНИЯ СОТОВОЙ КОНСТРУКЦИИ ФЮЗЕЛЯЖА
Для определения вектора коэффициентов 0 модели теплового состояния сотовой конструкции будем определять минимум функции Ф (0) взвешенной суммы квадратов невязок [10] с помощью итерационного алгоритма минимизации, использующего производные функции Ф (0) . Для этой цели предлагается использовать вариант стохастического квазиградиентного алгоритма с переменной метрикой [8], в котором приближения к точке минимума строятся по правилу
0к+1 = 0к - рк НкVк Ф, к = 0,1,..., (15)
где Нк - случайная квадратная матрица размером I х I; Vк Ф - градиент целевой функции в точке 0к; Рк - параметр шага.
Последовательность матриц Нк вычисляется по схеме
Н0 = I, Нк +1 +(I - цк Vк+1Ф(Дк+10)Г )Нк, Дк +10 = 0к+1 - 0к. (16)
Параметр цк выбирается из равенства цк = Ц /(^к+1Ф | -| Дк+101), где ц - константа, такая что 0 < ц <1.
На каждом шаге алгоритма происходит автоматическая подстройка параметра шага Рк. Если Ф(0к+1) > Ф(0к), то Рк+1 = Рк / У, где у > 1 - фиксированный параметр. Если Ф(0к+1) <Ф(0к), то выполняется следующая последовательность действий: Рк / = Рк'У, 0к+1,г = 0к - Ркг НкVk Ф, и вычисление Ф(0к+1,г'), г = 0, 1,. Эти действия производятся до тех пор, пока убывает значение функции Ф и при этом выполняются условия: Ртт - Ркг - Ртах (Ртт, Ртах - минимальная и максимальная длина шага соответственно) и г < гтах (гтах - заданное максимальное число итераций по
увеличению шага). Значения 0к+1, Рк+1 полагаются равными соответствен-
пк+1, г
но значениям 0 ' , Рк г , которые равны минимальному из полученных таким образом значений Ф (0).
Параметрическую идентификацию математической модели теплового состояния блоков оборудования в отсеках предлагается проводить по композиции метода наискорейшего спуска, метода Ньютона и квазиньютоновского метода Бройдена-Флетчера-Гольдфарба-Шэнно [11]. При решении жесткой системы обыкновенных дифференциальных уравнений предлагается использовать неявный метод Розенброка второго порядка [12].
4. ДОВЕРИТЕЛЬНЫЕ ОБЛАСТИ РЕЗУЛЬТАТОВ ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ
При определении погрешностей оценок вектора искомых коэффициентов модели © = [©д необходимо использовать совместную доверительную область оценок вектора искомых коэффициентов модели. Форма совместной доверительной области, как и вид ковариационной матрицы Р(0)
погрешностей оценок коэффициентов 0 , зависит от особенностей модели, вида и количества параметров теплообмена математической модели и других значимых факторов.
Для рассматриваемой нами нелинейной модели процессов теплообмена
У = Р(0,и,X); X е (0,хепе1); У(0,и,0) = У0; Р,У е Ят; 0е Яг
(17)
эта задача в настоящее время не имеет достаточно конструктивного решения.
В уравнении (17) использованы следующие обозначения: У = [У1 ]т=1 -вектор параметров теплообмена; и; - вектор параметров режима полета.
Приближенное уравнение совместной доверительной области будет определяться для следующего функционала, соответствующего рассматриваемой математической модели:
N
Ф(0) = - Р (0,и, ) Я-1 [Г; - Р (0,и] )] , j=1
(18)
где N - число дискретных значений параметров теплообмена; Я; - ковариационная матрица погрешностей параметров теплообмена.
Будем предполагать, что путем минимизации Ф(0) по 0 с использованием одного из методов [10 - 12] найдены оптимальные оценки 0 , достаточно близкие истинным значениям 00 вектора искомых параметров 0 . Для этого функционал (18) С.В. Епифановым и С.А. Каплуном [13] был преобразован к следующему виду:
Ф(0, © )=£
! ;=1
У; Р (0, и; ) + Р (0, и; ]-
;
;
)-РI0, и;)_
я; 1 *
У;-Р (^ и; )+ Р ^ и; )- Р (©, и;)
= Ф1 (00 ) + Ф2 (0, 0 ) + Фз (0, (I),
(19)
где
Ф
N
(0 )=?,
; =1
-г
У;-Р (©, и; ) + Р ^ и; ^ Д,-1 - Р (©, ^ )] ;
Ф2 (©, ©) = I [Р(©, и, )-Р(©, и,) ]=1
*71
р (4 и, )-р ( и})
Ф
N
.(©, © )=!
у ' }=1
7,- р (©,
)]Ч-1 [Р (© и,)-Р ( ид
ха-
Если в уравнении (3) положить © = ©0, то функционал Ф2 (©0, ©) , рактеризующий отклонения получаемых оценок © от ©0 , будет равен
Ф2 (©0, ©) = Ф(©0, ©)-Ф (©)-Фз (©0, ©). (20)
Будем также полагать, что погрешности в 7 являются центрированными случайными величинами, распределенными по нормальному закону. Тогда функционалы ф(©0, ©) и Ф1 (©) являются независимыми случайными
величинами, распределенными по закону х2 со степенями свободы т и (т - г) соответственно [10]. Следовательно, их разность также распределена по закону х2 с г степенями свободы. Кроме того, используя приведенную в [10] линеаризацию модели (17) в виде
7, = Н]
(©-© 0),
(21)
можно показать, что Ф31 ©0
Ф3 (©0, ©) = 0, где Н, - (т х г) - матрица коэффициентов влияния. Тогда, исходя из (19), с учетом (21) и приведенных выше выводов по функционалам ф(©0, ©), Ф (©) и Ф3 (©0, ©) методами дисперсионного анализа [10] можно получить следующее приближенное выражение для границы совместной доверительной области в пространстве исходных
коэффициентов © :
N
1[Р (©, и, )-Р (©0, и, -1 [Р (©, и, )-Р (©0, и, )] = х2-а (г), (22)
3=1
22 где Х1-а (г) - квантиль х -распределения, соответствующий вероятности 1 - а.
Предполагая достаточную близость © и ©0 , подставим в выражение (22) линеаризованную модель (21) и получим следующее более простое выражение для совместной доверительной области, описывающее в пространстве © г-мерный эллипсоид с центром в ©0 :
(©-©0 )т 4©-©0 )=х2-а (г);
(23)
А = [а; ] =
2
N т к«; N т к«; к,-,г ^ ^ «^ ^ Ч1 Чг
;=1 ««=1
;=1 «=1
N т к;г к;-1 N т Л« ^ ^ «;1 ^ ^ «;г
;=1 «=1
;=1«=1 о«-
(24)
где к«;д - коэффициент влияния д-го искомого коэффициента на «-й параметр теплообмена в ;-й момент времени; о2 - дисперсия «-го параметра теплообмена.
Очевидно, что А, являясь матрицей Грама векторов к; =
«;д о«
, адек-
ватна обратной ковариационной матрице Р (0) погрешностей их оценок [10].
Для случая большой размерности вектора 0 (г >2) использование совместной доверительной области в выражении (23) сопряжено с вычислительными трудностями. Поэтому можно ввести некоторые условные совместные
доверительные интервалы А0^ каждого из искомых параметров в виде проекций совместной доверительной области на соответствующие координатные
оси пространства 0 , что эквивалентно замене эллиптической области на описанный вокруг нее параллелепипед. В работе [10] для них получено следующее выражение:
«=1
где
А0? =±,
XI2- а (г )
Р1 АРд
'д-1 1 ед +1
• ег ]
= 1, — , г,
Т Е =-с-1 о
Т
0д = [а1д ^а(д- 1)д а(д +1)д ••• агд ] ;
(25)
сд =
а11 ^ 1( д -1) а1( д +1) ^ а1г
а(д-1)1 ^а(д-1)(д-1) а(д-1)(д+1) • • • а(д- 1)г а(д +1)1 • а(д +1)(д-1) а(д +1)(д +1) • • • а(д +1)г
аг1 ••• аг(д-1) аг(д +1) ••• агг
Здесь ^ - вектор размером г х 1; Еу, - векторы размером (г -1) х 1; Су — (г — 1) х (г — 1) - матрица, которая получается из матрицы А вычеркиванием у-го столбца у-й строки.
Элементы матрицы А определяются как особенностями параметров
теплообмена и математической моделью, так и значимыми факторами вектора параметров режима полета и .
5. ОЦЕНИВАНИЕ ПАРАМЕТРОВ МАТЕМАТИЧЕСКОЙ МОДЕЛИ БЕСПИЛОТНОГО ЛЕТАТЕЛЬНОГО АППАРАТА
Проверка предложенного теоретического метода была выполнена для приборного отсека беспилотного летательного аппарата, который представляет собой негерметизированный отсек с сотовой конструкцией, продуваемый воздухом из системы обеспечения теплового режима. Воздух охлаждает расположенное в отсеке оборудование. Блоки бортового оборудования разделены воздушными прослойками. При этом была проведена оптимизация температуры и расхода воздуха системы обеспечения теплового режима приборного отсека.
Компоновка приборного отсека негерметизированного отсека летательного аппарата с сотовыми конструкциями представлена на рис. 2.
Рис. 2. Компоновка негерметизированного продуваемого теплоизолированного отсека с сотовыми конструкциями:
I—VIII - части отсека; 1-22 - блоки бортового оборудования; 23-31 - части обшивки; x, y, z - координаты; pV - плотность воздушной среды за бортом;
Vair out - воздушная скорость полета; Tair out - температура воздушной среды
за бортом; Tstm - температура воздуха на выходе системы обеспечения теплового
режима; Gstm - расход воздуха на выходе системы обеспечения теплового режима
_О
Толщина обшивки I-V частей отсека равна /¡_у = 2 • 10 м, VI-VIII ча-
_3
стей - /vi_ viii = 4 •lO м. Толщина теплоизоляции обшивки равна
_2
lins = 2 • 10 м.
Коэффициент теплопроводности обшивки I-V частей отсека равен VV = 2,11 •lO2 Вт/(м ■ К), а VI-VIII частей - A,VI_VIII = 1,63 • 10_1 Вт/(м ■ К). Коэффициент теплопроводности сотовой конструкции обшивки равен Xins = 8 -10_2 Вт/(м ■ К).
Адекватность математической модели теплового состояния негерметичного продуваемого теплоизолированного отсека была подтверждена для холодного и теплового типов климата [14].
Основным критерием при оптимизации температуры и расхода воздуха системы обеспечения теплового режима является температура воздуха в отсеке с регулируемой температурой не ниже 253,15 К и не выше 328,15 К [15]. При этом значения температуры и расхода воздуха могут находиться в пределах соответственно 283,15...293,15 К и 0,5...1,0 кг/c. Вектор коэффициентов модели
© = [ , Gstm ]T (26)
включает в себя необходимые характеристики (значения температуры воздуха Tstm в кельвинах (К) и расход воздуха Gstm в килограммах в секунду (кг/с)).
Оценки коэффициентов модели для холодного ® c и теплового © v типов климата для температуры и расхода воздуха соответственно равны
© c = [290,4 0,89 ]T,
© v = [281,6 0,96]T.
*
Совместные доверительные интервалы А© неопределенностей оценок коэффициентов (25) для холодного © c и теплового типов климата © v при доверительной вероятности ß = 0,99 соответственно равны
* от
А© c = [7,0 6,4-10_2]T ,
* о т
А© v = [7,8 5,1 • 10 ]T.
*
Совместные доверительные интервалы А© каждого из искомых коэффициентов получены по вышеизложенному методу.
ВЫВОДЫ
Предложен метод определения теплового состояния отсеков и систем беспилотного летательного аппарата на основе метода математического моделирования. Характерной особенностью рассматриваемого типа летательных аппаратов является наличие в обшивке фюзеляжа сотовых теплозащитных панелей. Расчет теплопереноса в такой панели проводится методом Монте-Карло на основе вероятностного представления решения параболического уравнения. При этом используется предложенный авторами комбинированный метод, в котором расчет траекторий диффузионного процесса в ячейках, заполненных воздухом, осуществляется методом блуждания по движущимся сферам, а по каркасу, ограничивающим пластинам и в их окрестности - методом Эйлера.
При параметрической идентификации математической модели теплового состояния сотовой конструкции фюзеляжа используется квазиградиентный алгоритм с переменной метрикой.
Предложенный метод позволяет оптимизировать теплофизические параметры отсеков и систем при проектировании летательных аппаратов.
СПИСОК ЛИТЕРАТУРЫ
1. Воронин Г.И. Системы кондиционирования на летательных аппаратах. - М.: Машиностроение, 1973. - 443 с.
2. Гусев С.А., Николаев В.Н. Численно-статистический метод для решения задач теплообмена в теплозащитных конструкциях сотового типа // Сибирский журнал науки и технологий. - 2017. - Т. 18, № 4. - С. 719-726.
3. Дульнев Г.Н., Тарновский Н.Н. Тепловые режимы электронной аппаратуры. - Л.: Энергия, 1971. - 248 с.
4. Дульнев Г.Н., Польщиков Б.В., Потягайло А.Ю. Алгоритмы иерархического моделирования процессов теплообмена в сложных радиоэлектронных комплексах // Радиоэлектроника. - 1979. - № 11. - С. 49-54.
5. НиколаевВ.Н., Гусев С.А., Махоткин О.А. Математическая модель конвективно-лучистого теплообмена продуваемого теплоизолированного негерметичного отсека летательного аппарата // Прочность летательных аппаратов. Расчет на прочность элементов авиационных конструкций. - Новосибирск: СибНИА, 1996. - Вып. 1. - C. 98-108.
6. Миснар А. Теплопроводность твердых тел, жидкостей, газов и их композиций. - М.: Мир, 1968. - 460 с.
7. Ладыженская О.А., Солонников В.А., Уральцева Н.Н. Линейные и квазилинейные уравнения параболического типа. - М.: Наука, 1967. - 736 с.
8. Соболев С.Л. Некоторые применения функционального анализа в математической физике. - Изд. 3-е. - М.: Наука, 1988. - 336 с.
9. Gusev S.A. Application of SDEs to estimating solutions to heat conduction équations with discontinuous coefficients // Numerical Analysis and Applications. - 2015. - Vol. 8, N 2. -P. 122-134.
10. ХиммельблауД. Анализ процессов статистическими методами: пер. с англ. - М.: Мир, 1973. - 957 с.
11. Gill P., Murray E. Quasi-Newton methods for unconstrained optimization // IMA Journal of Applied Mathematics. - 1971. - Vol. 9, N 1. - P. 91-108.
12. Артемьев С.С., Демидов Г.В., НовиковЕ.А. Минимизация овражных функций численным методом для решения жестких систем уравнений. - Новосибирск, 1980. - 13 с. - (Препринт / ВЦ СО АН СССР; № 74).
13. НиколаевВ.Н., СимбирскийД.Ф. Доверительные области результатов параметрической идентификации процессов теплообмена бортового оборудования самолета // Методы и средства исследования внешних воздействующих факторов на бортовое оборудование летательных аппаратов. - Новосибирск: СибНИА, 1991. - Вып. 2. - С. 11-15.
14. Николаев В.Н. Экспериментально-теоретические исследования теплового состояния приборного отсека фоторазведчика // Журнал Сибирского федерального университета. Техника и технологии. - 2011. - Т. 4, № 3. - Р. 337-347.
15. ГОСТ РВ 20.39.304-98. Комплексная система общих технических требований. Аппаратура, приборы, устройства и оборудование военного назначения. Требования стойкости к внешним воздействующим факторам. - М.: Госстандарт России, 1998. - 55 с.
Гусев Сергей Анатольевич, доктор физико-математических наук, старший научный сотрудник Института вычислительной математики и математической геофизики СО РАН, профессор Новосибирского государственного технического университета. Основные научные интересы - разработка методов решения прикладных задач с использованием численного решения стохастических дифференциальных уравнений и оптимизация параметров динамических систем. Автор свыше 80 научных трудов. E-mail: sag@ osmf.sscc.ru.
Николаев Владимир Николаевич, доктор технических наук, начальник отдела ФГУП «Сибирский научно-исследовательский институт авиации им. С.А. Чаплыгина». Основные научные интересы - математическое моделирование теплового состояния отсеков и систем самолета при проектировании и летных испытаниях, электромагнитная совместимость радиоэлектронного оборудования самолета. Автор свыше 100 научных трудов. E-mail: nikvla50@mail.ru.
DOI: 10.17212/1814-1196-2018-2-23-38
Mathematical simulation of the thermal state of a pilotless vehicle unpressurized compartment
S.A. GUSEV1'2, V.N. NIKOLAEV3
1 ICMMG SB RAS, 6, Lavrent'eva Prospekt, Novosibirsk, 630090, Russian Federation, D. Sc. (Phys. & Math.), senior researcher in ICMMG SB RAS. E-mail: sag@osmf.sscc.ru Novosibirsk State Technical University, 20, K. Marx Prospekt, Novosibirsk, 630073, Russian Federation, D. Sc. (Phys. & Math.), professor. E-mail: sag@osmf.sscc.ru 3 FSUE S.A. Chaplygin Siberian Aeronautical Research Institute, 21, Polzunov St., Novosibirsk, 630051, Russian Federation, D. Sc. (Eng.), head of the FSUE department. E-mail: nikvla50@ mail.ru
A mathematical model of the thermal state of the aircraft compartment is developed which describes the heat exchange of the honeycomb structure skin made of carbon fiber and the process of heat transfer of the on-board equipment and air environment. Heat transfer in the honeycomb panel is described by the boundary value problem for a parabolic equation with a discontinuous diffusion coefficient and the boundary conditions of the third kind.
To solve the direct problem of the fuselage honeycomb structure thermal state, the Monte Carlo method is used on the basis of the probabilistic representation of the solution in the form of a mathematical expectation of a diffusion process functional.
Smoothing this coefficient by integral averaging is used in the proposed method. An approximate solution of the problem with a smoothed coefficient is obtained using the probabilistic representation of its solution. This representation is an expectation of a diffusion process corresponding to the boundary value problem. To obtain an approximate solution of the problem it is necessary to simulate a huge number of paths of the diffusion process in the area under consideration. The Euler method was earlier used for this purpose. But this approach requires significant computational costs. In this paper the method is modified by using the method of
*
Received 30 March 2018.
random walk on moving spheres in subareas corresponding to cells of the panel. The new approach allows us to significantly reduce computation costs.
The inverse problem of the honeycomb structure heat exchange is solved by minimizing the weighted sum function of the squared residuals using an iterative stochastic quasigradient algorithm.
The developed mathematical model of the compartment thermal state was used for optimizing the temperature and airflow of the thermal control system of the instrumental ventilated thermal-insulated compartment of the aircraft.
Keywords: mathematical model, thermal state, honeycomb structure, heterogeneous structures, numerical solution, parabolic boundary value problem, Monte Carlo method, discontinuous coefficients
REFERENCES
1. Voronin G.I. Sistemy konditsionirovaniya na letatel'nykh apparatakh [Air-conditioning systems onboard the aircrafts]. Moscow, Mashinostroenie Publ., 1973. 443 p.
2. Gusev S.A., Nikolaev V.N. Chislenno-statisticheskii metod dlya resheniya zadach teploo-bmena v tep-lozashchitnykh konstruktsiyakh sotovogo tipa [Numerical and statistical method for solving the heat exchange problems in the heat-protective constructions of the honeycomb type]. Sibirskii zhurnal nauki i tekhnologii - Scientific Journal of Science and Technology, 2017, vol. 18, no. 4, pp. 719-726.
3. Dul'nev G.N., Tarnovskii N.N. Teplovye rezhimy elektronnoi apparatury [Thermal conditions of electronics]. Leningrad, Energiya Publ., 1971. 248 p.
4. Dul'nev G.N., Pol'shchikov B.V., Potyagailo A.Yu. Algoritmy ierarkhicheskogo modeliro-vaniya protsessov teploobmena v slozhnykh radioelektronnykh kompleksakh [Algorithms for hierarchical modeling of heat transfer processes in complex electronic systems]. Radioelektronika - Electronics, 1979, no. 11, pp. 49-54.
5. Nikolaev V.N., Gusev S.A., Makhotkin O.A. [Mathematical model of the convective radiant heat exchange of the venting heat-insulated unpressurized bay of the aircraft]. Prochnost' letatel'nykh apparatov. Raschet na prochnost' elementov aviatsionnykh konstruktsii [Strength calculation of the airframe elements. A series of the aircraft strength]. Novosibirsk, SibNIA Publ., 1996, iss. 1, pp. 98108. (In Russuan).
6. Missenard A. Conductivite' thermique des solides, liquides, gaz et de leurs melanges [The thermal conductivity of solids, liquids, gases and their compositions]. Paris, Eyrolles, 1965 (Russ. ed.: Misnar A. Teploprovodnost' tverdykh tel, zhidkostei, gazov i ikh kompozitsii. Moscow, Mir Publ., 1968. 460 p.).
7. Ladyzhenskaya O.A., Solonnikov V.A., Ural'tseva N.N. Lineinye i kvazilineinye uravneniya pa-rabolicheskogo tipa [Linear and quasilinear equations of parabolic type]. Moscow, Nauka Publ., 1967. 736 p.
8. Sobolev S.L. Nekotorye primeneniya funktsional'nogo analiza v matematicheskoi fizike [Applications of functional analysis in mathematical physics]. 3rd ed. Moscow, Nauka Publ., 1988. 336 p.
9. Gusev S.A. Application of SDEs to estimating solutions to heat conduction equations with discontinuous coefficients. Numerical Analysis and Applications, 2015, vol. 8, no. 2, pp. 122-134.
10. Himmelblau D. Process analysis by statistical methods. New York, Wiley, 1970. 463 p. (Russ. ed.: Khimmel'blau D. Analiz protsessov statisticheskimi metodami. Translated from English. Moscow, Mir Publ., 1973. 957 p.).
11. Gill P., Murray E. Quasi-Newton methods for unconstrained optimization. IMA Journal of Applied Mathematics, 1971, vol. 9, no. 1, pp. 91-108.
12. Artem'ev S.S., Demidov G.V., Novikov E.A. Minimizatsiya ovrazhnykh funktsii chislennym metodom dlya resheniya zhestkikh sistem uravnenii [Minimization of ravine functions by numerical method for the stiff sets of equations solving]. Preprint no. 74. Data center of Siberian Department of the Academy of Sciences of the USSR. Novosibirsk, 1980. 13 p.
13. Nikolaev V.N., Simbirskii D.F. [The confidence regions of the parametric identification results of the heat exchange processes of the aircraft airborne equipment]. Metody i sredstva issledo-vaniya vneshnikh vozdeistvuyushchikh faktorov na bortovoe oborudovanie letatel'nykh apparatov [Methods and means of study of the external factors effecting on the aircraft airborne equipment]. Novosibirsk, SibNIA Publ., 1991, iss. 2, pp. 11-15. (In Russuan).
14. Nikolaev V.N. Eksperimental'no-teoreticheskie issledovaniya teplovogo sostoyaniya pribor-nogo otseka fotorazvedchika [Theoretical and experimental investigation of the photographic reconnaissance plane instrument bay thermal state]. ZhurnalSibirskogofederal'nogo universiteta. Tekhnika i tekhnologii - Journal of Siberian Federal University. Engineering & Technologies, 2011, vol. 4, no. 3, pp. 337-347.
15. GOST RV 20.39.304-98. Kompleksnaya sistema obshchikh tekhnicheskikh trebovanii. Ap-paratura, pribory, ustroistva i oborudovanie voennogo naznacheniya. Trebovaniya stoikosti k vnesh-nim vozdeistvuyushchim faktoram [State Military Standard 20.39.304-98. The integrated system of general technical requirements. Military equipment, instruments, devices and apparatus. Requirements of resistance to the effecting external factors]. Moscow, Gosstandart Rossii Publ., 1998. 55 p.).
Для цитирования:
Гусев С.А., Николаев В.Н. Математическое моделирование теплового состояния негерме-тизированного отсека беспилотного летательного аппарата // Научный вестник НГТУ. - 2018. -№ 2 (71). - С. 23-38. - doi: 10.17212/1814-1196-2018-2-23-38.
For citation:
Gusev S.A., Nikolaev V.N. Matematicheskoe modelirovanie teplovogo sostoyaniya negermet-izirovannogo otseka bespilotnogo letatel'nogo apparata [Mathematical simulation of the thermal state of a pilotless vehicle unpressurized compartment]. Nauchnyi vestnik Novosibirskogo gosudarstven-nogo tekhnicheskogo universiteta - Science bulletin of the Novosibirsk state technical university, 2018, no. 2 (71), pp. 23-38. doi: 10.17212/1814-1196-2018-2-23-38.
ISSN 1814-1196, http://journals.nstu.ru/vestnik Science Bulletin of the NSTU Vol. 71, No 2, 2018, pp. 23-38
