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

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

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

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

Следовательно, при моделировании траекторий винеровского процесса в подобластях с однородной средой можно двигаться не по траекториям, которые получаются методом Эйлера, а по границам шаров, находящихся в этих подобластях. В этом состоит основная идея метода случайного блуждания по сферам [7]. Но при решении нестационарных задач требуется помимо точки первого выхода винеровского процесса на границу шара ещё моделировать время первого достижения границы шара. Как известно [8], распределение времени первого выхода винеровского процесса на границу шара совпадает с распределением времени первого достижения уровня, равного радиусу шара, одномерным бесселевским процессом. Но плотность времени первого достижения бесселевским процессом заданного уровня выражается сложной формулой в виде ряда, содержащего функции Бесселя первого рода и их положительные нули. По этой формуле очень затруднительно построить эффективный алгоритм для моделирования времени первого выхода. С другой стороны, вместо заданного уровня, который требуется достичь, можно подобрать некоторые кривые [9], зависящие от t , такие, что получаются достаточно простые и удобные для моделирования формулы для плотности времени первого достижения бесселевским процессом этих кривых. Этот метод используется в данной работе, и он называется блужданием по движущимся сферам.

В качестве области изменения пространственных переменных x1 , x2 , x3 в работе принимается прямоугольный параллелепипед G = (-11,11) х (-12,12) х (0,13). При этом G есть объединение двух непересекающих-ся подмножеств: G = G1 и G2, где G1 - подобласть, соответствующая каркасу и пластинам, ограничивающим панель, а G2 есть объединение подобластей, соответствующих ячейкам, содержащим теплоизоляционный материал. Каждая ячейка представляет собой внутренность правильной шестиугольной призмы, боковая поверхность которой параллельна оси x3 . Предполагается, что рассматриваемый процесс теплопередачи происходит на отрезке времени [0,T] и описывается следующей краевой задачей для уравнения теплопроводности:

dи       д Г, , . ди )

= / I b ( x ) I ,

д t i=1 дx l     дx J где b (x) =

x e G 1

x e G 2’

u ( x ,0) = Ф( x ),                    (2)

д u д x₁

x i =- l 1

— 0, ^

д x ₁

— 0, ^x 1 ^— l 1

дu дx 2

— 0, u x 2 —— 1 2        ⁵ x 2

— 0, ^x 2 ^{— 1} 2

X — дx3    n

• ³ x 3 =0

. д u

—X— дx3 x.-l x

= ^ai⁽ t ) ( u - v ⁽ t ) ) ,

= ^a2 ⁽ t ) ( u ^— V 2⁽ t ) ) ^.

В (1)–(6) использованы следующие обозначения: φ – начальное распределение температуры в панели; λ – коэффициент теплопроводности углепластика; b1 , b2 – коэффициенты температуропроводности углепластика и воздуха соответственно; α1 , α2 – коэффициенты теплообмена между панелью и внешней средой; ψ1 , ψ2 – температура внешней среды у нижней и верхней сторон панели соответственно.

Сглаживание коэффициентов . Поскольку теплозащитная панель состоит из двух материалов с различными теплофизическими свойствами, то в задаче (1)–(6) мы имеем уравнение теплопроводности с разрывным коэффициентом температуропроводности. Достаточно полное математическое исследование краевых задач для параболических уравнений с разрывными коэффициентами дано в книге [10]. Там же приведены условия существования и единственности обобщенного решения (теорема 5.1, гл. III) в пространстве функций P 2 ^0,1( Q T ) ( Q T — G x (0, T ) ) , обладающих

T—h свойством J J h—1(u (x, t + h) — u (x, t))2dxdt ^ 0

0G при

h ^ 0 , с нормой

I “ k — max I| u ( x , t )IL

Q T 0< t <  T          L 2 ^{( G} )

^            ) 2

J u x dxdt .

V О т

Как утверждается в [10] (теорема 4.5, гл. III), обобщенное решение краевой задачи с разрывными коэффициентами в параболическом уравнении в норме пространства Р 2 ^0,1 ( Q_T ) устойчиво относительно вариаций коэффициентов. Из этой теоремы применительно к задаче (1)-(6) следует, что если при m ^да равномерно ограниченная последовательность b ^{( m )} сходится почти всюду к b , тогда обобщенные решения задач с коэффициентами b ^{( m )} сходятся сильно в норме k,¹ ( Q T ) к обобщенному решению задачи (1)-(6).

Таким образом, при достаточно больших значениях m приближенное решение u(m) задачи (1)–(6) с разрывным коэффициентом температуропроводности можно получить, если решать задачу, в которой в уравнении (1) этот коэффициент заменить на его гладкую аппроксимацию b(m). В качестве такой заме- ны можно, например, рассматривать его сглаживание в окрестности поверхностей разрыва с помощью интегрального усреднения с бесконечно дифференцируемым финитным ядром [11]:

^{b ( m ) (} x ) ^—    J to p m ( I x ^— y I ) b ⁽ y ) dy ^—

I ^{x —} y |

m

^—4" J to ^f1x ^{y 1 b(}y)dy.

^Pm x—y

m^{V P m}'

При этом to₁(| £ |) — 0 при | £ |>1; J to₁(| £ |) — 1.

Предполагается, что p _m ^ 0 при m ^ да . Функции b (^m) имеют все производные любого порядка [6]. Поскольку b e Lq (G) (q > 0), то для любой подобласти G 'с G , отстоящей от границы G на расстоянии, не меньшем p_m, усреднение b(^m) сходится к b в L_q (G'),

) q

■ ^ 0

т. е. |b<-) ^—bq,G'

q

< sup J | b (x — v) — b (x) | dx I^v|

m VG'

при m ^ да . Таким образом, аппроксимации b(^m) сходятся к b в норме пространства L_q (G).

Известно также, что сходимость в L_q (G) влечет сходимость по мере, а из сходящейся последовательности функций по мере можно выбрать подпоследовательность, сходящуюся почти всюду [12]. Причем в [12] дан конструктивный метод построения такой подпоследовательности. В связи со сказанным далее будем считать, что такая подпоследовательность выбрана, и это есть последовательность b^(m).

В результате замены в (1)–(6) коэффициента температуропроводности функцией b^(m) приходим к следующей краевой задаче:

д u (^m)_=b (_m) Л д²u (^m)+ Д д b (^m) д u (^m) дt " ^(x)t; дxi ⁺h дx дx_z ’

u (^m)(x ,0) — Ф( x),

д u (^m) дx₁

— 0,

^x1^——l1

д u (^m) дx₁

— 0, ^x1^—l1

д u (^m) ^дx 2

— 0,

x 2^——l2

д u (^m) ^дx 2

— 0, x 2^—l 2

, д u (m) λ дx3

^{— a}1⁽t ) ( u^—V1⁽t ) ) ,

x3 —0

, д u (^m)

—X---- дx3

^{— a}2 ⁽t)(u^—V2⁽t)).

Далее в работе будут использоваться следующие обозначения: дG - граница области G; xA — индикаторная функция множества A. Будем предполагать, что заданные в граничных условиях функции α₁ , α₂ ,

• V1, V2 такие, что для некоторого 0 < а < 1 существует решение задачи (9)–(14) в пространстве Гельдера 3+_а3+а

H , 2 (QT) [10] (теорема 5.3, гл. IV). Доказано [13], что при выполнении условий существования и единственности решений задач (1)–(6) и (9)–(14) для любого е > 0 существует натуральное число mе такое, что при m > mе выполняется неравенство vraisup | и(m) - и |< е . Следовательно, выбирая радиус QT усреднения достаточно малым, можно с любой заданной точностью приблизить решение u(m) к u.

Вероятностное представление решения. Для решений краевых задач для параболических уравнений известны вероятностные представления [14; 15]. Чтобы получить вероятностное представление решения задачи (9)–(14), введём систему СДУ:

r

Xr = x + J с_vXG (X_v)dWT +

T — t rr

+ J avXg (Xv)dv + J n(Xv )Xag (Xv)dKv,    (15)

T — t                    T—t

Y_r = 1 + J ai(v)Yvxx₃(„)=o dK_v +

T — t

r

⁺J ^a2^(v)YvXX₃(v)=1₃^dKv,              ⁽¹⁶⁾

T — t

r

Zr = J ^a1^(v)V1(^v)YvXX₃(v)=0dKv⁺

T — t

r

⁺J ^a2^(v W^v)Yv XX₃( v )=1₃^dKv,         ⁽¹⁷⁾

T — t

r

Kr = J XdG (Xv)dKv,             (18)

T — t где x e G - начальная точка для Xv; W.^ - трехмер- ный винеровский процесс; стv = (2b(m)(Xv))2; av -( 5b(m) )

вектор с координатами , i = 1, 2, 3; n(P) – ( dxi )

внутренняя нормаль в точке P e dG ; K_v - локальное время пребывания процесса X_v на границе, т. е. скалярный возрастающий процесс, который растет, только когда K_vedG .

Вероятностным представлением решения задачи (9)–(14) является математическое ожидание и (x, t) = Ex,T—t (ф(Xt )Yt + Zt ),          (19)

где символом Ex T—t обозначено математическое ожидание относительно вероятностной меры Px T—t, соответствующей случайному процессу, стартующему в момент времени T-t из точки x, т. е. XT—t = x.

Приближенные значения u^(m) (x,t) можно получить методом Монте-Карло путем статистического моделирования траекторий системы СДУ (15)–(18). Для этой цели в работах [6; 13] предлагалось использовать метод Эйлера. Моделирование траекторий методом Эйлера для СДУ (15)–(18) осуществляется по следующей схеме:

Xi+1 = xi+ h2 0,4+ ha,(20)

x+1 = X,A+1 +(A,+1 K)nxR3\G (XA1),(21)

y,+1 = y, (1^{+ (a}1)i Xag⁽xi)^Ai+1^{K +}

+ C«2)iX ag (x, )Ai+1 K),(22)

z,+1 = z,⁺y,^(a1⁾,⁽V1), X ag⁽x^)AI+1^{K +}

+У, (a2),(V2), X ag (x,)A,+1K, k+1 = k,+A,+1 K,(24)

где , + 1 – номера узлов сетки; ξ – вектор трех независимых N(0,1) случайных величин; A,₊₁K = d (X A+1) -расстояние от X^А1 до dG в случае выхода X^А1 из области.

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

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

При решении нестационарных задач методом блуждания по сферам требуется каждый раз моделировать не только точку первого выхода из шара, но и время первого выхода. Для моделирования времени первого выхода из шара радиуса R в работе используется определение времени первого достижения уровня R бесселевским процессом, который является решением СДУ [8]:

t

St = S0 +J S— dr + Wt.

Одномерный случайный процесс (19) получается в результате применения формулы Ито к евклидовой норме трехмерного винеровского процесса с начальной точкой x e G. При этом выполняется равенство S0 = Il^xll.

Время первого выхода трехмерного винеровского процесса из шара радиуса R и время первого достижения случайным процессом уровня R имеют одинаковое распределение. Известны явные формулы для плотности времени первого достижения бесселевским процессом заданного уровня для начальных точек S0 = 0 и S₀> 0 [9]. Но эти формулы представлены в виде рядов, содержащих функции Бесселя первого рода и их положительные нули, и поэтому численностатистическое моделирование по ним крайне затруднительно. С другой стороны, вместо заданного уровня, который требуется достичь, можно подобрать некоторые кривые, зависящие от t, такие, что получаются достаточно простые и удобные для моделирования формулы для плотности времени первого достижения бесселевским процессом этих кривых. Каждая такая кривая определяется задаваемой положительной сигмаконечной мерой μ . Определим функцию

τ

g

(

a

) 3

e

H

,

где H – гамма-распределенная случайная величина с параметрами ₂ , ₃ . В ходе вычислений, после того,

1 ^Л

Г (t, x) ⁼q (t, 0, x) — [ q (t, y, x )p( dy), a • a₀

где q(t, y, x) – плотность вероятности перехода бессе-левского процесса; a – числовой параметр.

Из свойства плотности q следует, что функции q и r удовлетворяют параболическому уравнению

dv 1 52v 9 f 1 ) — =--т--1 — v I.

dt 2 dx²dx V x )

Пусть функция g(t) такая, что на кривой x = g(t) выполняется равенство r(t,x) = 0. Тогда эту кривую можно рассматривать в качестве граничного условия к уравнению (27). Причем с учетом того, что бессе-левский процесс стартует из заданной точки, решение полученной краевой задачи определяется единственным образом.

Обозначим τ_g время достижения бесселевским процессом этой подвижной границы, т. е. т_g = = inf {t > 0, St = g(t)}. Можно показать [9], что плотность времени достижения границы x - g(t) = 0 удовлетворяет уравнению

d f T „ p (t) =-— I r(t,x)dx g dt

1 V 0

Если положить ц(dy) = y²dy, то g(t) и p_g (t) опре-

деляются по формулам [9]

f r

g⁽t) =

2tln

a

Pg⁽t) =

g³(t) 2at

В этом случае получается простая формула для моделирования τ_g

как получено значение τ_g , радиус сферы для определения очередной точки положения случайного процесса определяется по первой формуле в (29). Выбор параметра a дает возможность определить радиусы моделируемых шаров так, чтобы они не выходили за границы текущей ячейки и чтобы они не были слишком малы. В случае, если расстояние до границы ячейки менее некоторой заданной величины £, вычисления производятся согласно методу Эйлера.

Результаты расчётов. Вычисления проводились для сотовой панели, каркас которой и ограничивающие пластины изготовлены из углепластика, а наполнителем служит воздух. Основные характеристики панели следующие: общая толщина панели – 0,035 м; толщина каждой из ограничивающих пластин панели – 0,001 м; толщина стенок сотового каркаса - 6 -10^-5 м; сотовыми ячейками являются правильные шестиугольники одинакового размера с длиной стороны, равной 0,0042 м; коэффициенты температуропроводности углепластика и воздуха равны соответственно 8,0 -10-4 и 2,36 -10^-5 м²/с.

Граничные условия для расчётов были сформированы на основании физических характеристик, полученных по данным лётного эксперимента в холодном климате и соответствуют первым 150 с взлёта самолета. На рис. 1 приведён график коэффициента теплообмена на внешней стороне обшивки.

На рис. 2 приведены полученные в результате расчётов при объёме выборки 2000 траекторий случайного процесса графики изменения значений температуры в середине панели и у края панели, обращённого внутрь салона.

В табл. 1, 2 приведены для сравнения затраты времени счёта предложенным комбинированным методом и только методом Эйлера с различными шагами по каркасу h_f и по объёму внутри ячеек h_a .

Результаты расчётов позволяют сделать следующие выводы:

• а)    предложенный метод значительно эффективнее использовавшегося ранее метода Эйлера;

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

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

Рис. 1. Коэффициент теплообмена на внешней стороне обшивки

• Fig. 1.    Heat transfer coefficient on the outside of the skin

  
  
  
  
  
  
  
  

  М

  
  

  нижний край панели

  
  
  
  
  
  

середина панели

235 L 0

^с

Рис. 2. Значения температуры внутри панели

• Fig. 2.    Temperature values inside the panel

Таблица 1

Затраты времени счёта в средней точке панели, с

Метод	ha = 2,5-10-7, hf = 2,5-10-8	ha = 5,0-10-7, hf = 5,0-10-8	ha = 10-6, hf = 10-7
Комбинированный метод	1064	548	154
Метод Эйлера	5839	3070	1563

Таблица 2

Затраты времени счёта вблизи границы x3 = 0, c

Метод	ha = 2,5 -10-7,hf = 2,5 -10-8	ha = 5,0.10-7, hf. = 5,0.10-8	ha = 10-6, hf = 10-7
Комбинированный метод	1394	479	188
Метод Эйлера	6193	4156	1854

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

Acknowledgments. S. A. Gusev’s work was carried out with the financial support of the Russian Foundation for Basic Research, grant No. 17-01-00698 A.