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

Математическая модель теплового состояния герметизированных отсеков сверхзвукового пассажирского летательного аппарата описывает процесс передачи теплоты по конструкциям летательного аппарата, тепловыделение или теплопоглощение бортовым оборудованием, системами кондиционирования и вентиляции [1–4].

Параметрическая идентификация математической модели теплового состояния герметизированных отсеков относится к решению некорректных задач в механике вязкоупругих сред [5] при использовании методов численного анализа математических моделей [6].

Использование теплоизолированных сотовых конструкций фюзеляжа является перспективным направлением в области авиационной техники [1; 7; 8].

Критерии оптимизации теплофизических характеристик герметизированного отсека определены в соответствии с нормами лётной годности АП 25 [9].

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

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

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

Упомянутые уравнения имеют следующий вид [10; 11]:

C cv ( x ) T cv,t = (4 v ( x ) Tv, x ) x , 0 <  x <  l hon , 0 <  t ^ t k ;                       (1)

^λcv⁽ x ^{) F} cv,in ^T cv, x

acv,out

( t ) F cv,in ( T cv,out⁽ t , ^x ) - ^T e⁽ t )) + Q cv ,out^,

x = 0;

^λcv⁽ x ^{) F} cv,in ^T cv, x

^acv,in ⁽ t ) F cv,in ^{( T} air,pr ⁽ t )   T cv,in ⁽ t ^{, x} )) + Q cv,in^,

x = l ;

T cv (0, x ) = T _o( x ), 0 <  x <  l hon ,                                   (4)

где C cv ( x ) =

x e hon ;

x e air ,

λ cv ( x )

x e hon ;

x e air ,

т. е. коэффициенты C _cv , λ_cv зависят от того, в каком слое рассматривается перенос тепла.

В уравнениях (1)–(4) использованы следующие обозначения: T _cv( t , x ) – температура сотовой конструкции; T _cv,t – первая производная T _cv по t ; t – время; T _cv,out – первая производная T _cv по x ; T _{cv, x , x} – вторая производная T _cv по x ; T _cv,out , T _cv,in – температура наружной поверхности сотовой конструкции и внутренней; l _hon – толщина сотовой конструкции; C _cv( x ) – объёмная теплоёмкость сотовой конструкции обшивки, определяемая теплоёмкостью дюралюминия С _hon и теплоёмкостью воздуха С _air ; λ_cv( l ) – теплопроводность сотовой конструкции, определяемая теплопроводностью дюралюминия λ_hon и теплопроводностью воздуха λ_air ; α_cv,out – коэффициент теплоотдачи наружной поверхности сотовой конструкции; α_cv,in – коэффициент теплоотдачи внутренней поверхности сотовой конструкции; F _cv,in – площадь сотовой конструкции при наружном и внутреннем теплообмене; Q _cv,out – тепловая энергия внешних источников; Q _cv,in – тепловая энергия внутренних источников; T _e – температура восстановления; T _{air, pr} – температура воздушной среды в герметизированном отсеке или в части отсека; T_j – температура j -го элемента отсека; l – толщина сотовой конструкции.

Вынужденная конвекция снаружи герметизированного отсека и граничные условия третьего рода . Коэффициент теплоотдачи α_cv,out наружной поверхности обшивки будем вычислять по формуле [10]

a cv,out = 3,26(Re Out ) ^{- 0,5} (Pr *ut ) ^{- 0,666} p * c * f (5)

* 6

для ламинарного пограничного слоя при критерии Рейнольдса Reout < 1·106 и по формуле [10]

a cv,out = 0,181(lg Re Out ) — ^2,584 (Pr out ) " ² / ³ P * C * < „-                        (6)

для турбулентного пограничного слоя при критерии Рейнольдса Г10⁶ <  Re O _ut <  Г10⁹.

В уравнениях (5), (6) использованы следующие обозначения: Re^*_out , Pr_o^*_ut – критерии соответственно Рейнольдса и Прандтля при температуре воздуха T ^* ; ρ _V ^* – плотность воздуха за бортом при температуре T ^* ; c _p^* – удельная теплоёмкость воздуха при температуре T ^* ; V _air,out – воздушная скорость полёта.

Температура воздуха T ^* определяется по формуле

*

T = T air,out + 0, ^{72 (} T e ^- T air,out ) ,

где T _air,out – температура воздуха за бортом за пределами теплового пограничного слоя.

Критерий Reout вычисляется по формуле

*   **

Reout = Vair,out pV 7cv,out / ^air > где Lcv,out – характерный размер для наружной поверхности обшивки; μ*air – динамическая вязкость воздуха.

Характерный размер L _cv,out принимается равным расстоянию от носовой точки фюзеляжа до середины рассчитываемого отсека или его части.

Динамическая вязкость μ*air определяется по выражению

μ^*_air = 0,1222229∙10^–6 + 0,682674∙10^–8 T ^* – 0,313155∙10^–11 ( T ^* )².

Критерий Pro*ut вычисляется по формуле

Pr*ut = cP^ir/Cr, где λ *air – теплопроводность воздуха.

Теплопроводность λair определяется по выражению

λair = 0,141483∙10–2 + 0,896161∙10–4 T* – 0,204759∙10–7 (T*)2.(11)

Плотность воздуха ρV* за бортом вычисляется по выражению p, = 3,4826-10 3 pair,out / T ,(12)

где p _air,out – давление воздуха за бортом.

Температура восстановления T _e в уравнении (2) определяется по выражению

Te = Tair,out(1 + M2 Г (k - 1) / 2),(13)

где r – коэффициент восстановления температуры для ламинарного пограничного слоя r = V Pr_out , для турбулентного - r = 3 Pr_out ; M - число Маха на внешней границе пограничного слоя; к = c p / C v - отношение удельной теплоёмкостей газа при постоянном давлении и объёме.

Вынужденная конвекция в герметизированном отсеке и граничные условия третьего рода. Коэффициент теплоотдачи α_cv,in внутренней поверхности сотовой конструкции можно вычислить по формуле [11]

a cv,in = 0,66 / - Re in ⁵ PrF / 4v,in                             (14)

для ламинарного пограничного слоя при критерии Рейнольдса Re _in< 4∙10⁴ и по формуле [11]

a cv,in = 0,037 ^Re^Pr⁰ ’ ⁴³/ L cv,in

при критерии Рейнольдса Re in >  440 8 .

Теплопроводность λ _air и динамическая вязкость μ _in в формулах (14), (15) в критериях Re, Pr определяются соответственно по выражениям (11), (9) для температуры воздуха в отсеке T _air .

Произведение J_W плотности ρ_air и скорости W _air воздуха в отсеке или массовая скорость J_W в критерии

Rein   pair Wair Lin / Pair принимается из результатов эксперимента.

Характерный размер L _cv,in принимается по аналогии с L _cv,out равным расстоянию от начала отсека до середины рассчитываемой части отсека.

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

Сы (x) T,t = (Xbi (x) T,x) x,  0 < x < l, 0 < t < tk; (17) Xbl (x) Fbl,uprTbl, x =abl,upr (t) Fbl,upr (Tbl (t)  Tair, upr (t’ x))’  x = 0; (18) X bl (x) Fbl, prT, x = a bl, pr (t) Fbl, pr (T pr (t) — T (t, x)),  x = l; (19) Tbl (0, x) = T0(x),     0 < x < l, (20) где Cbi(x) =

^C compo ’ ^{x e com}P^o ;    , , ,

^X compo , x ^e comPo ;

Х_{Я| r}, x e air .

air

.         ^X bl ( x ) =

C air , x e air ,

В уравнениях (17)–(20) использованы следующие обозначения: Tbl(x, t) – температура сотовой перегородки; Tbl,t – первая производная Tbl по t; Tbl,x – первая производная Tbl по x ; Tbl x x — вторая производная Tbi по x ; Cbi (x) - объёмная теплоёмкость перегородки; Xbi (l) — теплопроводность перегородки; a bi, upr, a bi, pr - коэффициенты теплоотдачи поверхности пе- регородки соответственно со стороны негерметизированного и герметизированного отсеков; Fbl, upr, Fbl, pr – площадь перегородки, участвующая в конвективном теплообмене соответст- венно со стороны негерметизированного и герметизированного отсеков Tair,upr – температура воздушной среды в негерметизированном отсеке или в части отсека.

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

Tn t = a_alr m ( t ) F ir m / Cm ( T air ( t ) - Tn ) + m ,t      air, m      air, m    m air, pr        m

⁺ E g j , m / C m ^T j ⁽ t ) / ^T ms - c 0 ^e m F m / C m T m ⁺ Q m / C m ,

m где Tm – температура m-го бортового оборудования; Tm,t – первая производная Tm по t; Tms – эталонная температура; a air m - коэффициент теплоотдачи m-го бортового оборудования; Fair,m – площадь m-го бортового оборудования при конвективном теплообмене; Cm – теплоёмкость i-го бортового оборудования; c0 – постоянная Стефана – Больцмана; gj,m – коэффициент радиационного теплообмена системы «j-й элемент отсека – m-й блок бортового оборудования»; εm – коэффициент черноты излучения m-го блока; Qm – энергия тепловыделения или теплопоглощения m-м бортовым оборудованием от системы кондиционирования и преобразованная из электрической энергии.

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

Tr ,t = a air, r ( t ) F air, r / C r T air, pr ( t ) - T- ) +

+£ gj, r / Crjt) / Tms -0 e rFr / CrT + Qr / Cr, (22) r где Tr - температура r-й конструкции; Tr t - первая производная Tr по t; a air r - коэффициент теплоотдачи r-й конструкции; Fair,r – площадь r-й конструкции при конвективном теплообмене; Cr – теплоёмкость r-й конструкции; Qr – энергия тепловыделения r-й конструкции; gj,r – коэффициент радиационного теплообмена системы «j-й элемент отсека – r-я конструкция»; εr – коэффициент черноты излучения r-й конструкции.

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

T air, pr ,t    ^a cv,in ⁽ t ) F cv,in, pr / C air ^{( T} cv,in, pr ^{( l} cv, pr ’ t )    T air, pr ) ⁺

^+a bl , pr ⁽ t ) F bl , pr / C air (T bl , pr ^{( l} bl , pr , t )    T air, pr ) ⁺

⁺ E ^{( a} air, r ⁽ t ) F air, r / C air ⁽ T r ⁽ t ) - T air, pr ) ^{+ Q} r / C air ) ⁺ r

+ E ^{( a} air, p ^p

+ E ^{( a} air, m m

⁽ t) F air, p / ^C air ^{( T} p ⁽ t ) - T air, pr ) ^{+ Q} p / ^C air ) ⁺

( t ) F air, m / C air ( T n ( t ) - T air, pr ) + Q m / C air ) +

+E Kir, j (t) Fair, j / Cair( Tj (t) - Tair, pr ) + Cp Gstm / Cair(TStm — Tair, pr ), j где Tair,pr,t – первая производная Tair по t; Gstm – расход воздуха, вытекающего из системы кондиционирования и (или) из системы обеспечения теплового режима; Tstm – температура воздуха, вытекающего из системы кондиционирования и (или) из системы обеспечения теплового режима; Tcv,in, pr – температура внутренней поверхности сотовой конструкций обшивки.

Теплоёмкость воздуха C _air определяется по выражению

С air    ^c p ^{( p} air G stm ^Л t ^{+ p} air ^ air ),

где pair - плотность воздуха в отсеке; cp - удельная теплоёмкость воздуха; Лt - интервал дискретизации времени при решении системы дифференциальных уравнений; Vair – объём воздуха в отсеке.

Теплоёмкость воздуха C _air определяется по выражению (24).

Коэффициенты теплоотдачи поверхностей acv,out, acv,in, a bl, upr, а bl, pr , а air, m, aair, r в уравнениях (2), (3), (18), (19), (21) – (23) будем вычислять с помощью методов, описанных в работах [10–13].

Коэффициенты радиационного теплообмена в уравнениях (21), (22) определяется методом Монте-Карло [14].

Алгоритм решения прямой задачи теплового состояния сотовой конструкции обшивки

Для решения прямой задачи теплового состояния сотовой конструкции фюзеляжа рассмотрим математическую модель теплообмена в гетерогенной структуре (1)–(4), (17)–(20) как параболическую краевую задачу с разрывными коэффициентами. Описание и физические свойства гетерогенных структур можно найти в [15].

Теплообмен в гетерогенной структуре описывается следующим параболическим уравнением:

d T d t

d T

Ef a ij ( t , x ) 7 T = 0, ( t , x ) ^G Q t , i ^1 1 ^d x i

i

d x

' j 7

где Q_T = (0, t k ) x G ; T ( t , x )

– температура конструкции; a_ij – коэффициенты, являющиеся

Липшиц непрерывными по x функциями в подобластях G(k), к = 1,...,M . Применительно к сотовой панели коэффициенты aij представляют собой коэффициенты температуропроводности

' 4v ( x )

a ij ( t , x ) = 5 ^C cv^{( x} )’

i = J ,

0, i ^ j .

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

иЕ -2

i              i,ji

Также требуется, чтобы искомая функция удовлетворяла условию:

Tt=0 = ф( x)(28)

и граничному условию на дG:

(

Е ^aijni ч i, j

д T , _ , , — + n( t, x) T + Y( t, x) дX;

j 7

= 0.

x Ед G

В [16] доказано существование обобщенных решений краевых задач (25), (27), (25) и (29). Причём эти решения могут быть аппроксимированы решениями соответствующих краевых задач с коэффициентами, которые являются приближениями исходных разрывных коэффициентов. Например, можно получить приближённое решение исходной задачи, решая задачу со сглаженными коэффициентами. В данной работе предлагается статистически оценивать решение задачи со сглаженными коэффициентами с помощью метода, основанного на численном решении стохастических дифференциальных уравнений и, таким образом, получать оценки приближённого решения исходной краевой задачи с разрывными коэффициентами.

Для сглаживания коэффициентов предлагается использовать интегральное усреднение [16; 17]

f^(p)(x) = p ³j Mp(| ^x - y I)f⁽y)dy                            ⁽³⁰⁾

Ix - y |

p и j Ир(| ^ |)d^ = 1 . В (19) и далее р означает радиус усреднения, а символ |а| обозначает евк-

I ^р лидову норму вектора a .

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

Pешение параболического уравнения может быть представлено в виде математического ожидания функционала решения стохастических дифференциальных уравнений [17]. Этот факт может быть использован для получения статистических оценок решений параболических уравнений с помощью численного решения стохастических дифференциальных уравнений. Приближенное решение задачи таким способом может быть получено в одной или нескольких точках внутри области, что часто бывает достаточным для практического применения. При этом нет необходимости строить сетку по пространственным переменным и решать большие системы линейных алгебраических уравнений. Метод статистического моделирования легко распараллеливается с высокой эффективностью. Поэтому для решения задачи можно использовать суперкомпьютерную технику. Мы применяем этот метод [18] для оценки решения параболического уравнения со сглаженными коэффициентами, полученными на основе интегрального усреднения (30).

Приближенное решение уравнения (25) будем находить как решение следующего уравнения дT  у д (р) дT

/ aj

д t   i j i дxi Ч 3  дxj J где ai(jρ) – сглаженные коэффициенты уравнения (25) в окрестности Г. Соответственно гранич- ное условие, соответствующее сглаженным коэффициентам, имеет вид

(                                У

Е a^(p)^^ + П( t, x) T + Y( t, x)

I i, j^{j     дx}J                       7

= 0.

x ед G

Для точки (t, x) е Qt мы определим случайный процесс X•, который начинает движение из точки x в момент времени t и является решением следующего векторного стохастического дифференциального уравнения [19; 20]

Xv = x + J b (xr, r) dr + J CT (X_r, r) dWr + J nA (X_r, r) dlkrl,               (33)

T-t                    T-t                        T-t где W• – винеровский процесс [16]; nA – нормализованный вектор внутренней конормали,

v т. е. nA = An/ | An |; | kv |= J 1gg (Xr)d | kr | - неотрицательный случайный процесс, который рас- t тет только тогда, когда процесс X• достигает

2стTст = A^(p),

границы; ст - (3 x 3) - матрица такая, что,

AW = ( a    b = [£

k i=1

^dxi ’

£<, £3<'I i=1 ^dxi      i=1 ^dxi )

.

Дополнительно введем обозначения: E_t,x – математическое ожидание относительно вероятностной меры P_t,x , соответствующей случайному процессу, исходящему из точки x в момент времени t ; 1_S – функция множества S .

Таким образом, мы можем получить оценки решения задачи (31), (32), (28) путём численного моделирования траектории процесса X_•. Для этого мы используем модифицированный метод Эйлера, согласно которому приближённые траектории X_•рассчитываются по формуле [21]

xi+1 = xi+ ^hb(xi,ti) ⁺hb^ст(xi,ti^)Zi +(^Ai+1K)nA, ti+1 = ti + h, (35)

Ai+1К = [d (xi + hb⁽xi,ti) + hh^ст(xi,ti^)Zi)] , ⁽³⁶⁾ где nA - единичный вектор внутренней конормали в точке xi, если xi находится на 5G;

[ a ] = max {0, - a }; d (x) - неположительная вещественная функция, удовлетворяющая для любой точки x £ G следующему уравнению x = р( x) + d (x) nA (p( x)).                                    (3 7)

В уравнении (37) используются следующие обозначения: р(x) - проекция точки x £ G на

G^* в направлении конормали; n^A(ρ(x)) – единичный вектор внутренней конормали в ρ(x) .

Заметим, что d (x) = 0, если x e G.

Для получения приближенного решения (оценки) задачи (31), (32), (28) требуется также вычислять следующие функции, определенные на сетке по времени:

f                   1, i = 0,

yi =

zi =

( i-1

к

i exp £ П(tk, xk Яаg (xk)Ak+1К , i ^ 1,

k k=0

0, i = 0,

-

£ (Y⁽tk, ^xk )1ag^(xk^)Ak+1^K) Ук, i ^ ¹

i-1

I k=0

Оценка решения задачи (20), (21), (15) определяется по формуле T⁽t,^x) = Et_k -t,x [ф (^xN ) yN⁺zN ] ,

где N=

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

Yt = F (Y (t, 0 )), t e (0,11); Y_t= Y 0 , F, Y e RS ; 0 e Rr,            (40)

где Y = [Tcv, Tbl, Tm , Tr,... ]^T - вектор параметров теплового состояния отсека; Y _t - вектор первых производных Y по t; © = [и1,и2,...,и5]^T - вектор коэффициентов модели; T -верхний индекс, обозначающий операцию транспонирования.

Для решения системы уравнений (40) предлагается использовать следующую численную схему типа Розенброка второго порядка аппроксимации для неавтономных систем [22]

^Yn+1 = ^Yn + ^aK1 +^{(1 —}a)^K2^;

K1 = h(I - ah^_Y (Yn,t_n,©))^-1Ф(Yn,t_n + ah,©);

K2 = h(I - ahФ_у (Yn,tn,©))^-1Ф(Yn,t_n + aK1,tn + 2ah,©);

a = 1-1/V2

где, Y_n, Y_n+1 - решение системы, полученной на n-й и (n + 1)-й итерациях, соответственно; Ψ – правая часть системы; Ψ_Y – матрица Якоби; I – единичная матрица; h – шаг интегрирования.

Алгоритм параметрической идентификации математической модели теплового состояния сотовой конструкции обшивки

Решение обратной задачи, т. е. оценивание коэффициентов Θ модели (1)–(4), (17)–(20), сводится к минимизации суммы квадратов невязок между заданными по принятому критерию значениями Y и соответствующими значениями Y^ˆ(t, Θ) , полученными в ходе расчётов по уравнениям модели:

N

Ф(©) = £ (Yk - Yk (tk, ©))^T (Yk - Yk (tk, ©)),                      (44)

k=1

где t_k – моменты времени при k = 1,...,N.

Для этого применяется итерационный алгоритм минимизации, использующий производные функции Φ(Θ) . Для этой цели предлагается использовать вариант стохастического квазигра-диентного алгоритма с переменной метрикой [23], в котором приближения к точке минимума строятся по правилу

Θ^k+1=Θ^k-ρ_kΗ^k∇^kΦ, k=0,1,K,                         (45)

где H^k - случайная квадратная матрица размера l х l; V^kФ - градиент целевой функции в точке Θ^k; ρ_k – параметр шага.

Последовательность матриц Η^k вычисляется по схеме:

_Η0 _{= I, Η}k+1 _{+(I -μk∇}k+1_{Φ ⋅(Δ}k+1_Θ)T_{) Η}k_,

_Δk+1_Θ=Θk+1 _-Θk

Параметр μ_k выбирается из равенства μ_k = μ/( ∇^k+1Φ ⋅ Δ^k+1Θ , где μ – константа такая, что 0 < μ <1.

На каждом шаге алгоритма происходит автоматическая подстройка параметра шага ρ_k . Если Ф(0k+1) > Ф(0k), то рк₊1 = рк /у, где у > 1 - фиксированный параметр. Если Φ(Θ^k+1) < Φ(Θ^k), то выполняется следующая последовательность действий рki = рк Y, Θ^k+1,i = Θ^k - ρ_k,i Η^k∇^k Φ и вычисление Φ(Θ^k+1,i), i = 0,1,K.

Эти действия производятся до тех пор пока убывает значение функции Φ(Θ) и при этом выполняются условия: ρ_min≤ ρ_k,i ≤ ρ_max (ρ_min,ρ_max – минимальная, максимальная длина шага, соответственно) и i< i_max(i_max– заданное максимальное число итераций по увеличению шага). Значения Θ^k+1, ρ_k+1 полагаются равными соответственно значениям Θ^k+1,i, ρ_k,i , которые равны минимальному из полученных значений Φ(Θ) .

Алгоритм параметрической идентификации математической модели теплового состояния внутренних элементов отсеков

Тепловое состояние внутренних элементов отсеков описывают обыкновенные дифференциальные уравнения (21)–(23).

Как было отмечено в работе [23] для минимизации функции (44) целесообразно использовать композицию метода наискорейшего спуска [24], квазиньютоновского метода Бройдена – Флетчера – Гольдфарба – Шэнно [25] и метода Ньютона [26], которые реализуются в соответствии с формулой:

Θj+1 = Θj +aj ⋅S(Θj), (47) где a_j – коэффициент, характеризующий длину шага на j-й итерации; S – параметр, указывающий направление поиска вектора Θ₀действительных значений коэффициентов Θ .

Очередное направление S поиска j вектора Θ в этом алгоритме определяется из системы уравнений:

∇²Φ(Θ_j)⋅S(Θ_j)=-∇Φ(Θ_j), (48) где ∇²Φ(Θ_j) - (r × r) – матрица Гессе, представляющая собой квадратную матрицу вторых частных производных функции Φ по вектору Θ ; ∇Φ(Θ_j) – градиент функции Φ в текущей точке Θ _j .

Начальная матрица ∇ ²Φ ( Θ _k ) в уравнении (48) была принята единичной.

Для решения системы уравнений (48) матрицу ∇²Φ(Θ _j) представляют в факторизованной форме:

∇²Φ(Θ_j)=L(Θ_j)D(Θ_j)L^T(Θ_j), (49) где L(Θ _j) – нижнетреугольная матрица с единичной диагональю; D(Θ _j) – диагональная матрица.

Матрицы L(Θ _j), D(Θ _j) получают разложением Холесского матрицы ∇²Φ(Θ _j) по алгоритму, описанному в работе [25].

Для минимизации функции (44) квазиньютоновским методом Бройдена – Флетчера – Гольдфарба – Шэнно очередное направление поиска Θ _j определяется из системы уравнений (48).

Оценивание матрицы Гессе вторых частных производных в этом алгоритме проводится по формуле Бройдена – Флетчера – Гольдфарба – Шэнно [25].

В процессе минимизации с использованием квазиньютоновского алгоритма на каждой ите-2

рации требуются вычисления градиента функции ∇²Φ(Θ _j) . Компоненты градиента функции (48) вычисляются по формуле

N ∂Φ/∂ϑ1=2∑(Yk-Yk(tk,Θ))^T∂Yk(tk,Θ)/∂ϑ1,                      (50)

k=1

где H_1,k = ∂Y^ˆ_k(t_k, Θ) / ∂ϑ₁ – производные от решения уравнений (21)–(23) по ϑ₁ , которые называются функциями чувствительности [25].

Применение алгоритма (47) – (50) к системе (21) – (23) и функциям чувствительности

Y_t,Θ = f_YY_Θ + f_Θ,

YΘ(Θ)=0                                    (51)

требует вычисления на каждом шаге системы уравнений размера S × (r + 1) .

В системе (51) f_Y – матрица первых частных производных f по Y размера S × S . f_Θ– матрица первых частных производных f по Θ размера S × r . Матрица системы (51) содержит на главной диагонали блоки I - ah f_Y , а в первых S столбцах находятся элементы, включающие производные f_Y,Y и f_Y,Θ; f_Y,Y – матрица первых частных производных по Y элементов матрицы f_Y размера S × S × S ; f_Y,Θ– матрица первых частных производных по Θ элементов матрицы f_Y размеры S × S × r .

Рис. 1. Сопоставление рассчитанных и измеренных значений температуры T a r, pr в отсеке летательного аппарата:

• 1    – реализация средних арифметических рассчитанных в нескольких точках значений температуры; 2 – рассчитанные значения температуры; 3 – границы измеренных в нескольких точках значений температуры; 4 – доверительные интервалы неопределённостей измерений

Fig. 1. Comparison of calculated and measured temperature values T in the aircraft compartment: a r, pr

• 1    – implementation of arithmetic averages calculated at several points of temperature values;

• 2    – calculated temperature values; 3 – boundaries of temperature values measured at several points;

• 4 – confidence intervals of measurement uncertainties

Описанный алгоритм был исследован нами путём численного моделирования для ряда задаваемых (эталонных) значений вектора параметров режима полёта U = [ρ_V , V_air,out, T_e]^T (где ρ_V – плотность воздушной среды за бортом; V_air,out – воздушная скорость полёта; T_e – температура восстановления) и искомых коэффициентов Θ₀ , которым приводились в соответствие рассчитанные по (21)–(23) и зашумлённые величины вектора измерений Y . Изучались устойчивость, скорость и точность сходимости алгоритма в зависимости от формы и точности задания исходных условий и других факторов. При этом критериями качества являлись величины функции невязки Φ(Θ) последней итерации и величины (Θ^ˆ_j - Θ₀) / Θ₀ , определяющие неопределённости оценок коэффициентов Θ^ˆотносительно их действительных значений Θ₀ .

По результатам моделирования нами был сделан вывод об удовлетворительной сходимости предложенного алгоритма. Конечные результаты незначительно зависят от неопределённостей задания начальных значений Θ_int .

Проверка адекватности построенной модели реальным процессам была проведена по методу анализа остатков [26] для случая аддитивных случайных некоррелированных неопределённостей, распределённых по нормальному закону, так и для случая произвольных статистических свойств вектора неопределённостей.

Результаты проверки адекватности модели приведены на рис. 1.

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

Оценивание коэффициентов математической модели теплового состояния носового приборного продуваемого отсека

В качестве объекта исследования был принят сверхзвуковой летательного аппарата (рис. 2).

Параметры режима полёта и воздушной среды за бортом модели применения летательного аппарата приведены на рис. 3.

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

Характеристики размещенных в отсеке блоков бортового оборудования представлены в таблице.

При оптимизации толщины сотовой конструкции обшивки, температуры и расхода воздуха системы кондиционирования отсека вектор коэффициентов модели имеет следующий вид:

Для холодного типа климата начальная температура воздуха должна соответствовать 283,15 К, для экстремального теплого сухого – 315,15 К. Температура воздуха на выходе из системы кондиционирования воздуха в полёте не должна быть ниже 276,15 и выше 333,15 К.

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

Gstm(t) =ϑG,0+ϑG,1ρV(t) M²(t),                               (52)

где ϑ_G,0,ϑ_G,1 – коэффициенты модели, которые будут оцениваться.

Θ=[lhon, Tstm, JG,0, JG,1]^T                                      (53)

и включает в себя необходимые характеристики (значения толщины сотовой конструкции l_hon в м, температуры воздуха T_stm , T_cl в К и коэффициенты ϑ_G,0,ϑ_G,1 , описывающие расход воздуха системы кондиционирования G_stm в кг/с).

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

• Fig. 2.    Compartment of onboard equipment on a supersonic aircraft

Рис. 3. Параметры режима полёта и воздушной среды за бортом модели применения летательного аппарата: pair,out, Па;        М;          Tair,out, К; pair,out – давление воздуха за бортом за пределами теплового пограничного слоя; M – число Маха на внешней границе пограничного слоя;

T_air,out – температура воздуха за бортом за пределами теплового пограничного слоя

• Fig. 3.    Parameters of the flight mode and air environment overboard the aircraft application model: pair,out^, Pa^;         М;           T_air,out, К; p_air,out – air pressure outside the thermal boundary layer;

M – is the Mach number at the outer boundary of the boundary layer; T_air,out is the air temperature overboard outside the thermal boundary layer

Характеристики размещенных в отсеке блоков бортового оборудования

Номер части отсека	Номер блока бортового оборудования	Площадь теплоотдающей поверхности блока, Fair, i , м2	Масса блока m, кг	Энергия тепловыделения блоком Qi , Вт
I	1	0,081	1,0	0
	2	0,330	3,5	41
	3	0,230	18,0	1200
	4	0,124	2,5	180
II	5	0,488	10,0	100
	6	0,125	8,6	378
	7	0,180	8,5	27
	8	0,325	8,5	66
	9	0,450	10,0	38
	10	0,249	6,0	60
III	11	1,005	35,0	463
IV	12	0,105	3,0	50
	13	0,081	1,0	0
	14	0,054	8,0	0
	15	0,054	8,0	0
	16	0,273	4,5	146
	17	0,762	25,0	400
V	18	1,151	36,0	430
VI	19	0,984	40,0	285
VII	20	0,530	14,0	15
	21	0,390	12,0	600
VIII	22	0,118	3,3	15

Оценки коэффициентов модели для температуры и расхода хладагента соответственно равны

Θ = [0,06 275,7 0,002 0,001]^T.

Вектор P_D неопределённостей оценок коэффициентов Θ Он составляет при доверительной вероятности β = 0,99 соответственно

P_D = [0,008 3,7 0,0008 0,00007]^T .

Заключение

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

В качестве методов параметрической идентификации модели теплового состояния отсеков при проектировании предложено использовать композицию метода наискорейшего спуска, метода Ньютона и квазиньютоновского метода Бройдена – Флетчера – Гольдфарба – Шэнно. Для решения прямой задачи сотовых конструкций предложено применить метод, основанный на численном решении стохастических дифференциальных уравнений. При решении жёсткой системы обыкновенных дифференциальных уравнений, описывающих теплообмен бортового оборудования, человека, груза, кресла, конструкций, характеристик системы кондиционирования воздуха – неявный метод Розенброка второго порядка аппроксимации для неавтономных систем.

В качестве методов параметрической идентификации модели теплового состояния отсеков при проектировании предложено использовать вариант стохастического квазиградиентного алгоритма с переменной метрикой и композицию метода наискорейшего спуска, метода Ньютона и квазиньютоновского метода Бройдена – Флетчера – Гольдфарба – Шэнно.

Проведена проверки адекватности модели по анализу остатков.

Исследования проведены в соответствии с Нормами лётной годности Авиационных правил АП 25.