The aircraft hydraulic system units and pipelines heat exchange parameters study

Full Text

Введение

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

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

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

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

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

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

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

$T_{cov,t}={\vartheta }_1{\text{\hspace{0.17em}ρ}}_V\left(t\right)V_{air,out}\left(t\right)\left[T_e\right(t)-T_{cov}(t\left)\right]-{\vartheta }_2{\text{\hspace{0.17em}β}}^{{\vartheta }_{17}}\left(t\right)\left[T_{cov}\right(t)-T_{air}(t\left)\right]-$

$-{\vartheta }_3\left\{{\left[\frac{T_{cov}\left(t\right)}{100}\right]}^4-{\left[\frac{T_{eq}\left(t\right)}{100}\right]}^4\right\}-{\vartheta }_4\left\{{\left[\frac{T_{cov}\left(t\right)}{100}\right]}^4-{\left[\frac{T_{eq,ef}\left(t\right)}{100}\right]}^4\right\}+$                               (1)

$+{\vartheta }_{18}{\text{\hspace{0.17em}β}}^{{\vartheta }_{17}}\left(t\right)\left[T_e\right(t)-T_{cov}(t\left)\right]+{\vartheta }_{19}{\left[\frac{p_V\left(t\right)}{T_{cv}\left(t\right)}\frac{d{p}_V\left(t\right)}{dt}\right]}^{{\vartheta }_{17}}\times \left[T_e\right(t)-T_{cov}(t\left)\right]+{\vartheta }_{21}{Q}_{cov,out}.$

$T_{air,t}={\vartheta }_5{\text{\hspace{0.17em}β}}^{{\vartheta }_{17}}\left(t\right)\left[T_{cov}\right(t)-T_{air}(t\left)\right]+{\vartheta }_6\text{\hspace{0.17em}β}\left(t\right)\left[T_{cov}\right(t)-T_{air}(t\left)\right]+$

$+{\vartheta }_{20}\frac{p_V\left(t\right)}{T_{cov}\left(t\right)}\frac{d{p}_V\left(t\right)}{dt}\left[T_{cov}\right(t)-T_{air}(t\left)\right]-{\vartheta }_7{\text{\hspace{0.17em}β}}^{{\vartheta }_{17}}\left(t\right)\left[T_{air}\right(t)-T_{eq}(t\left)\right]-$

$-{\vartheta }_8{\text{\hspace{0.17em}β}}^{{\vartheta }_{17}}\left(t\right)\left[T_{air}\right(t)-T_{eq,ef}(t\left)\right]+{\vartheta }_{22}.$                                                                                                       (2)

$T_{eq,t}={\vartheta }_9{\text{\hspace{0.17em}β}}^{{\vartheta }_{17}}\left(t\right)\left[T_{air}\right(t)-T_{eq}(t\left)\right]+{\vartheta }_{10}\left\{{\left[\frac{T_{cov}\left(t\right)}{100}\right]}^4-{\left[\frac{T_{eq}\left(t\right)}{100}\right]}^4\right\}+$

$+{\vartheta }_{11}\left\{{\left[\frac{T_{eq,ef}\left(t\right)}{100}\right]}^4-{\left[\frac{T_{eq}\left(t\right)}{100}\right]}^4\right\}+{\vartheta }_{12}(T_{lq}-T_{eq}).$                                                                            (3)

$T_{eq,ef,t}={\vartheta }_{13}{\text{\hspace{0.17em}β}}^{{\vartheta }_{17}}\left(t\right)\left[T_{air}\right(t)-T_{eq,ef}(t\left)\right]+{\vartheta }_{14}\left\{{\left[\frac{T_{cov}\left(t\right)}{100}\right]}^4-{\left[\frac{T_{eq,ef}\left(t\right)}{100}\right]}^4\right\}-$

$-{\vartheta }_{15}\left\{{\left[\frac{T_{eq}{T}_{eq,ef}\left(t\right)}{100}\right]}^4-{\left[\frac{T_{eq,ef}\left(t\right)}{100}\right]}^4\right\}+{\vartheta }_{16}.$                                                                                      (4)

где $T_{lq}$ – температура гидрожидкости, протекающей в агрегате или трубопроводе; $T_{cov},T_{air},T_{eq},T_{eq,ef}$ в выражениях (1)–(4) с индексом t означает дифференцирование по времени t температуры обшивки, воздушной среды в отсеке, поверхности теплоизоляции агрегата или трубопровода, эффективного оборудования, окружающего исследуемый агрегат или трубопровод соответственно; $p_V$ – плотность воздуха за бортом; $\Theta ={\left[{\text{υ}}_{\text{1}}{\text{,υ}}_{\text{2}}{\text{,\hspace{0.17em}\hspace{0.33em}…\hspace{0.17em},υ}}_{\text{17}}\right]}^T$ – вектор ко- эффициентов модели; T – верхний индекс, обозначающий операцию транспонирования.

Вместо расчётных значений ${\text{α}}_{eq}$ коэффициента теплоотдачи конвекцией исследуемого агрегата или трубопровода введено произведение измеряемых плотности воздуха за бортом ${\text{ρ}}_V$ и числа Маха M полёта летательного аппарата при $\text{M}<\text{1}$ и произведение ${\text{ρ}}_V{\text{\hspace{0.05em}M}}^{\text{2}}$ при $\text{M}\ge 1$.
В общем виде процесс теплопередачи в углепластиковой теплоизоляции описывается уравнениями [5; 6]:

$C_{cv}\left(x\right)T_{cv,t}={\left({\lambda }_{cv}\right(x\left)T_{cv,x}\right)}_x,0<x<l,0<t\le t_k;$                                                                        (5)

${\lambda }_{cv}\left(x\right)F_{cv}{T}_{cv,x}={\text{α}}_{cv,out}\left(t\right)F_{cv}\left(T_{cv}\right(t,x)-T_{air}(t\left)\right)+Q_{cv,out}^{},x=0;$                                             (6)

${\lambda }_{cv}\left(x\right)F_{cv}{T}_{cv,x}={\alpha }_{cv,in}\left(t\right)F_{cv}\left(T_{cv,in}\right(t)-T_{cv}(t,x\left)\right)+Q_{cv,in}^{},x=l;$                                                          (7)

$T_{cv}\left(\mathrm{0,}x\right)=T_0\left(x\right), 0<x<l,$                                                                                                                (8)

$\mathrm{где} C_{cv}\left(x\right)=\left\{\begin{array}{c}{C}_{compo},x\in compo,\\ C_{air},x\in air,\end{array}\right.$

${\lambda }_{cv}\left(x\right)=\left\{\begin{array}{c}{\lambda }_{compo},x\in compo,\\ {\lambda }_{air},x\in air,\end{array}\right.$

где коэффициенты $C_{cv }, {\lambda }_{cv}$  определяются рассматриваемым слоем переноса тепла.

В уравнениях (5)–(8) использованы следующие обозначения:

$C_{cv}\left(x\right)$ – объемная теплоемкость теплоизоляции, определяемая теплоемкостью композита ${С}_{compo}$ и теплоемкостью воздуха ${С}_{air}$ ; ${\lambda }_{cv}\left(l\right)$ – теплопроводность теплоизоляции, определяемая теплопроводностью композита ${\lambda }_{compo}$ и теплопроводностью воздуха ${\lambda }_{air}$ ; ${\alpha }_{cv,out}$ – коэффициент теплоотдачи наружной поверхности агрегата или трубопровода; ${\alpha }_{cv,in}$– коэффициент теплоотдачи внутренней поверхности агрегата или трубопровода; $F_{cv}$ – площадь агрегата или трубопровода при наружном и внутреннем теплообмене; $Q_{cv,out}$ – тепловая энергия внешних источников;$Q_{cv,in}$– тепловая энергия внутренних источников; $l$ – толщина теплоизоляции.

Решение прямой задачи теплового состояния элементов гидросистемы

Обыкновенные дифференциальные уравнения (1)–(4) составляют следующую систему уравнений:

$Y_t=F\left(Y\right(t,\Theta \left)\right),t\in \left(\mathrm{0,}{t}_t\right);Y_t=Y_{\Theta },Y\in R^s;\Theta \in R^r,$                                                            (9)

$\mathrm{где} \mathrm{Y}={[{\mathrm{T}}_{\mathrm{cov}},{\mathrm{T}}_{\mathrm{air}},{\mathrm{T}}_{\mathrm{eq}},{\mathrm{T}}_{\mathrm{eq},\mathrm{ef}}]}^{\mathrm{T}}$– вектор параметров теплового состояния гидросистемы;  – вектор первых производных Y по t.

Решение жёсткой системы (9) обыкновенных дифференциальных уравнений предлагается проводить методом Розенброка по следующей схеме [7]:

${\overrightarrow{Y}}_{n+1}={\overrightarrow{Y}}_n+\alpha K_1+(1-\alpha )K_2;$                                                                                                        (10)

$K_1=h{(I-\alpha h{F}_Y({\overrightarrow{Y}}_n,t_n,\Theta \left)\right)}^{-1}F({\overrightarrow{Y}}_n,t_n+\alpha h,\Theta );$                                                                               (11)

$K_2=h{(I-\alpha h{F}_Y({\overrightarrow{Y}}_n,t_n,\Theta \left)\right)}^{-1}F({\overrightarrow{Y}}_n,t_n+\alpha K_1,t_n+2\alpha h,\Theta );\alpha =1-1/\sqrt{2}$                                     (12)

где ${\overrightarrow{Y}}_n$, ${\overrightarrow{Y}}_{n+1}$ – решение системы, полученной на n-й и (n+1)-й итерациях, соответственно; F – правая часть системы; F_y – матрица Якоби; I – единичная матрица; h – шаг интегрирования.

Решение системы дифференциальных уравнений с частными производными (5)–(8) предлагается определять методом Монте-Карло со сглаженными коэффициентами в виде математического ожидания функционала от диффузионного процесса [8–10].

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

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

Определение вектора коэффициентов Θ модели теплового состояния элементов гидросистемы будем определять минимумом функции Φ(Θ) взвешенной суммы квадратов невязок [11] с помощью итерационного алгоритма минимизации:

$\Phi \left(\Theta \right)=\sum _{j=1}^N{[{\overrightarrow{Y}}_j-\overrightarrow{F}(\Theta ,{\overrightarrow{U}}_j\left)\right]}^T{R}_j^{-1}[{\overrightarrow{Y}}_j-\overrightarrow{F}(\Theta ,{\overrightarrow{U}}_j)],$                                                                          (13)

где N – число точек по времени; ${\overrightarrow{U}}_j$ – вектор управления; $R_j$ – ковариационная матрица неопределённостей параметров теплообмена.

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

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

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

Критерием при оптимизации коэффициентов модели гидросистемы является температура гидрожидкости до 333 K. Температура воздуха в отсеке до 400 К.

 

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

p_air,out – давление воздуха за бортом за пределами теплового пограничного слоя; M – число Маха на внешней границе пограничного слоя; T_air,out – температура воздуха за бортом за пределами теплового пограничного слоя

Fig. 1. Parameters of the flight mode and the air environment overboard of the aircraft application model: p_air,out – air pressure outside the thermal boundary layer; – Mach number at the outer boundary of the boundary layer; T_air,out – air temperature outside the thermal boundary layer

 

Коэффициент теплопроводности углепластиковой теплоизоляции равен  λ_cv = 8∙10^–2  Вт/(м∙К). Толщина углепластиковой теплоизоляции  агрегата приняли за 2∙10^–2 м.

Полученные оценки коэффициентов модели Θ агрегата имеют следующие значения:

$\mathrm{\Theta }$ = [1,2584∙10–4	4,1542∙10–1	5,3417∙10–2	1,2215∙10–2	5,3456∙10–3
2,0357∙10–3	9,2045∙10–3	1,1904∙10–1	3,9123∙10–2	3,5162∙10–1
2,6877∙10–3	2,0979∙10–2	3,2077∙10–4	2,0343∙10–4	1,6850∙10–1
1,3344∙10–4	5,1202∙10–1]T.

Cовместные доверительные интервалы искомых коэффициентов модели

При большой размерности вектора коэффициентов модели  классическое использование совместной доверительной области [13] (рис. 2) связано со значительными проблемами.

 

Рис. 2. Линии равного уровня функции невязки Φ(Θ) = 2734 вблизи действительных значенийкоэффициентов модели ${\mathbf{\text{υ}}}_{\mathbf{\text{1}}}{\mathbf{\text{,υ}}}_{\mathbf{\text{2}}}{\mathbf{\text{,\hspace{0.17em}\hspace{0.33em}…\hspace{0.17em},υ}}}_{\mathbf{\text{5}}}$ 

Fig. 2. Equal level lines of the residual function Φ(Θ) = 2734 near the actual values of the coefficients of the model ${\mathbf{\text{υ}}}_{\mathbf{\text{1}}}{\mathbf{\text{,υ}}}_{\mathbf{\text{2}}}{\mathbf{\text{,\hspace{0.17em}\hspace{0.33em}…\hspace{0.17em},υ}}}_{\mathbf{\text{5}}}$

 

В связи с этим было предложено использовать совместные доверительные интервалы  коэффициентов в виде проекций совместной доверительной области на координатные оси пространства Θ. Это соответствует замене гиперэллиптической области на описанный вокруг неё гиперпараллелепипед. Было получено для $\Delta {\Theta }_q^{\ast }$ следующее выражение [14; 15]:

$\Delta {\Theta }_q^{\ast }=\pm \sqrt{\frac{{\text{χ}}_{1-\alpha }^2\left(r\right)}{{\overrightarrow{F}}_q^{T}A{\overrightarrow{F}}_q}},q=\mathrm{1,}\dots ,r$                                                                                                       (14)

${\overrightarrow{F}}_q={\left[e_1\dots e_{q-1}1e_{q+1}\dots e_r\right]}^T;{\overrightarrow{E}}_q=-C_q^{-1}{\overrightarrow{D}}_q,$

${\overrightarrow{D}}_q={\left[a_{1q}\dots a_{(q-1)q}{a}_{(q+1)q}\dots a_{rq}\right]}^T.$                                                                                                 (15)

${\text{C}}_q=\left[\begin{array}{c}{a}_{11}\dots {}_{1(q-1)}{a}_{1(q+1)}\dots a_{1r}\\ \dots \dots \dots \\ a_{(q-1)1}\dots a_{(q-1)(q-1)}{a}_{(q-1)(q+1)}\dots a_{(q-1)r}\\ a_{(q+1)1}\dots a_{(q+1)(q-1)}{a}_{(q+1)(q+1)}\dots a_{(q+1)r}\\ \dots \dots \dots \\ a_{r1}\dots a_{r(q-1)}{a}_{r(q+1)}\dots a_{rr}\end{array}\right],$                                                                         (16)

где ${\text{χ}}_{1-\text{α}}^2$ – распределение; ${\overrightarrow{F}}_q-(r\times 1)$ – вектор;  ${\overrightarrow{E}}_q,{\overrightarrow{D}}_q-(r-1)\times 1$  – векторы; $C_q-(r-1)\times (r-1)$ – матрица, которая получается из матрицы A вычеркиванием q -го столбца q -й строки.

При этом оценки доверительных интервалов $\Delta {\Theta }^{*}$ полученных коэффициентов агрегата при доверительной вероятности β = 0,95 имеют следующие величины:

${\mathrm{\Delta \Theta }}^{*}$ = [2,4339∙10–6	4,5678∙10–3	2,7423∙10–5	6,5632∙10–3	3,2545∙10–5
4,5579∙10–5	2,2781∙10–6	1,9907∙10–3	8,5377∙10–4	3,4863∙10–3
6,3324∙10–6	4,6360∙10–4	4,6578∙10–7	3,2878∙10–6	5,5735∙10–4
3,7469∙10–6	7,4590∙10–5]T

Заключение

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

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

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

Построена математическая модель теплового состояния агрегата гидросистемы в негерметическом отсеке летательного аппарата и оценены доверительные интервалы каждого из искомых коэффициентов модели с использованием χ²_1-α  распределения при доверительной вероятности β = 0,95.

About the authors

Vladimir N. Nikolaev

S. A. Chaplygin Siberian Aeronautical Research Institute

Author for correspondence.
Email: nikvla50@mail.ru

Dr. Sc., Chief Researcher

Russian Federation, 21, Polzunov St., Novosibirsk, 630051

References

Supplementary files