NUMERICAL AND STATISTICAL METHOD FOR SOLVING THE HEAT EXCHANGE PROBLEMS IN THE HEAT-PROTECTIVE CONSTRUCTIONS OF THE HONEYCOMB TYPE

全文:

Введение. При конструировании авиационной и космической техники важнейшее значение приобре- тает устройство тепловой защиты летательных аппа- ратов. Одним из возможных вариантов решения этой проблемы является применение теплозащитных пане- лей сотового типа [1-5]. Сотовая теплозащитная панель представляет собой конструкцию, состоящую из параллельно расположенных пластин, между кото- рыми заключен тонкий каркас в виде пчелиных сот, заполненных некоторым веществом с низкой тепло- проводностью. Основными преимуществами таких средств тепловой защиты являются хорошие тепло- физические свойства и небольшой вес при относи- тельно невысокой стоимости и достаточной прочности. Данная статья посвящена математическому модели- рованию процессов теплообмена в таких теплозащит- ных конструкциях и разработке метода расчёта для определения их теплового состояния в реальных условиях. В работе [6] для расчёта теплопереноса в сотовой панели предлагался численно-статистический метод, основанный на вероятностном представлении реше- ния краевой задачи для уравнения теплопроводности и численном решении стохастических дифференци- альных уравнений (СДУ). При этом делалось сглажи- вание разрывных коэффициентов уравнения тепло- проводности на основе интегрального усреднения, а численное решение СДУ осуществлялось методом Эйлера. С помощью метода Эйлера возможно получение хороших статистических оценок решения уравнения теплопроводности, но большая трудоёмкость такого подхода даже при хорошей возможности распаралле- ливания и применении суперкомпьютерной техники не позволяет проводить реальные расчёты для доста- точно больших интервалов полётного времени. С другой стороны, если теплофизические свойства материалов, из которых состоит сотовая панель, незначительно изменяются в рассматриваемом диапа- зоне температур и внутри панели отсутствуют источ- ники тепловой энергии, то возможна модификация алгоритма расчёта, позволяющая значительно умень- шить вычислительные затраты. Действительно, если коэффициенты температуропроводности материалов сотовой панели постоянны, а передача тепловой энер- гии внутри панели осуществляется исключительно теплопроводностью, то случайный процесс в одно- родной среде, который присутствует в вероятностном Следовательно, при моделировании траекторий винеровского процесса в подобластях с однородной средой можно двигаться не по траекториям, которые получаются методом Эйлера, а по границам шаров, находящихся в этих подобластях. В этом состоит основ- ная идея метода случайного блуждания по сферам [7]. Но при решении нестационарных задач требуется помимо точки первого выхода винеровского процесса на границу шара ещё моделировать время первого достижения границы шара. Как известно [8], распре- деление времени первого выхода винеровского про- цесса на границу шара совпадает с распределением времени первого достижения уровня, равного радиусу шара, одномерным бесселевским процессом. Но плотность времени первого достижения бесселевским процессом заданного уровня выражается сложной формулой в виде ряда, содержащего функции Бесселя первого рода и их положительные нули. По этой фор- муле очень затруднительно построить эффективный алгоритм для моделирования времени первого выхо- да. С другой стороны, вместо заданного уровня, кото- рый требуется достичь, можно подобрать некоторые кривые [9], зависящие от t, такие, что получаются достаточно простые и удобные для моделирования формулы для плотности времени первого достижения бесселевским процессом этих кривых. Этот метод используется в данной работе, и он называется блуж- данием по движущимся сферам. В качестве области изменения пространственных переменных x1 , x2 , x3 в работе принимается прямо- угольный параллелепипед G = (-l1, l1) ´ (-l2 , l2 ) ´ (0, l3 ) . При этом G есть объединение двух непересекающих- ся подмножеств: G = G1 È G2 , где G1 - подобласть, соответствующая каркасу и пластинам, ограничи- вающим панель, а G2 есть объединение подобластей, соответствующих ячейкам, содержащим теплоизоля- ционный материал. Каждая ячейка представляет собой внутренность правильной шестиугольной призмы, боковая поверхность которой параллельна оси x3 . Предполагается, что рассматриваемый процесс тепло- передачи происходит на отрезке времени [0,T] и опи- сывается следующей краевой задачей для уравнения теплопроводности: = å ¶u 3 ¶ æ ¶u ö ¶t ¶x ç b(x) ¶x ÷ , представлении, есть ни что иное, как винеровский процесс. Точка первого выхода винеровского процес- са на границу шара имеет на ней равномерное распре- деление. где i=1 i è i ø í b(x) = ì b1, x Î G1 , (1) îb2 , x Î G2 u(x, 0) = φ(x) , (2) ны можно, например, рассматривать его сглаживание в окрестности поверхностей разрыва с помощью ин- = 0, x1 =-l1 x1 =l1 = 0, (3) тегрального усреднения с бесконечно дифференци- руемым финитным ядром [11]: = 0, = 0, (4) b(m) (x) = ò ω m ρ x- y <ρm (| x - y |) b( y)dy = x2 =-l2 x2 =l2 1 æ | x - y | ö ¶u = ρ3 ò ω1 ç ρ ÷ b( y)dy. (8) ¶x λ = α1(t) (u - ψ1(t)) , (5) 3 x =0 m x- y <ρm è m ø 1 ò 1 3 -λ = α (t) (u - ψ (t)) . (6) При этом ω (| ξ |) = 0 при | ξ |> 1 ; ω (| ξ |) = 1 . |ξ|£1 2 2 x3 =l Предполагается, что ρm ® 0 при m ®¥ . Функции b(m) В (1)-(6) использованы следующие обозначения: φ - начальное распределение температуры в панели; λ - коэффициент теплопроводности углепластика; имеют все производные любого порядка [6]. Посколь- ку b Î Lq (G) ( q > 0 ), то для любой подобласти G¢ Ì G , отстоящей от границы G на расстоянии, не меньb1 , b2 - коэффициенты температуропроводности шем ρm , усреднение b(m) сходится к b в Lq (G¢) , углепластика и воздуха соответственно; α1 , α2 - 1 коэффициенты теплообмена между панелью и внешней средой; ψ1 , ψ2 - температура внешней среды т. е. b(m) - b æ ç q,G¢ £ sup ç ò | b(x - v) q öq b(x) | dx ÷÷ ® 0 у нижней и верхней сторон панели соответственно. 2 T T Сглаживание коэффициентов. Поскольку тепло- защитная панель состоит из двух материалов с раз- личными теплофизическими свойствами, то в задаче (1)-(6) мы имеем уравнение теплопроводности с разрывным коэффициентом температуропроводно- сти. Достаточно полное математическое исследование краевых задач для параболических уравнений с раз- рывными коэффициентами дано в книге [10]. Там же приведены условия существования и единственности обобщенного решения (теорема 5.1, гл. III) в пространстве |v|£ρm è G¢ ø при m ®¥ . Таким образом, аппроксимации сходятся к b в норме пространства Lq (G) . Известно также, что сходимость в Lq (G) влечет сходимость по мере, а из сходящейся последователь- ности функций по мере можно выбрать подпоследо- вательность, сходящуюся почти всюду [12]. Причем в [12] дан конструктивный метод построения такой подпоследовательности. В связи со сказанным далее будем считать, что такая подпоследовательность функций свойством V 0,1(Q ) (Q = G ´(0,T )) , обладающих T -h ò ò h-1(u(x, t + h) - u(x, t))2dxdt ® 0 при выбрана, и это есть последовательность b(m) . В результате замены в (1)-(6) коэффициента тем- пературопроводности функцией b(m) приходим к сле- 0 G дующей краевой задаче: h ® 0 , с нормой ¶u(m) (m) 3 ¶2u(m) 3 ¶b(m) ¶u(m) 1 ¶t = b (x)å ¶x2 + å ¶x ¶x , (9) æ ö2 i=1 i i=1 i i Q L (G) ò x (m) u = max T 0£t £T u(x, t) 2 + ç u2dxdt ÷ ç ÷ Q . (7) u (x, 0) = φ(x), (10) è T ø Как утверждается в [10] (теорема 4.5, гл. III), обобщенное решение краевой задачи с разрывными коэффициентами в параболическом уравнении в нор- ¶u(m) ¶x1 x1 =-l1 = 0, ¶u(m) ¶x1 x1=l1 = 0, (11) 2 T ме пространства V 0,1 (Q ) устойчиво относительно ¶u(m) ¶u(m) вариаций коэффициентов. Из этой теоремы примени- тельно к задаче (1)-(6) следует, что если при m ®¥ ¶x2 x2 =-l2 = 0, ¶x2 x2 =l2 = 0, (12) равномерно ограниченная последовательность b(m) ¶u(m) сходится почти всюду к b, тогда обобщенные решения задач с коэффициентами b(m) сходятся сильно в нор- λ ¶x 3 x3 =0 = α1(t) (u - ψ1(t)), (13) 2 T ме V 0,1 (Q ) к обобщенному решению задачи (1)-(6). -λ ¶u (m) = α (t) (u - ψ (t)). (14) Таким образом, при достаточно больших значениях ¶x 2 2 m приближенное решение u(m) задачи (1)-(6) с раз- 3 x3 =l рывным коэффициентом температуропроводности можно получить, если решать задачу, в которой в уравнении (1) этот коэффициент заменить на его гладкую аппроксимацию b(m) . В качестве такой заме- Далее в работе будут использоваться следующие обозначения: ¶G - граница области G; χ A - индика- торная функция множества A. Будем предполагать, что заданные в граничных условиях функции α1 , α2 , ψ1 , ψ2 такие, что для некоторого 0 < α < 1 сущест- Приближенные значения u(m) (x, t) можно полувует решение задачи (9)-(14) в пространстве Гельдера 3+α, 3+α H 2 (QT ) [10] (теорема 5.3, гл. IV). Доказано [13], что при выполнении условий существования и един- ственности решений задач (1)-(6) и (9)-(14) для лючить методом Монте-Карло путем статистического моделирования траекторий системы СДУ (15)-(18). Для этой цели в работах [6; 13] предлагалось исполь- зовать метод Эйлера. Моделирование траекторий методом Эйлера для СДУ (15)-(18) осуществляется бого ε > 0 существует натуральное число mε такое, по следующей схеме: что при m > mε выполняется неравенство D 1 vrai sup | u(m) - u |< ε . Следовательно, выбирая радиус Xi +1 = xi + h 2 σiξi + hai , (20) QT усреднения достаточно малым, можно с любой заданxi +1 = X D 1 + (Di +1K ) ni χ R3 \G ( X D ), (21) i + i +1 ной точностью приблизить решение u(m) к u. yi+1 = yi (1+ (α1)i χ¶G (xi )Δi+1K + Вероятностное представление решения. Для решений краевых задач для параболических уравне- + (a 2 )i χ¶G (xi )Di +1 K ), (22) ний известны вероятностные представления [14; 15]. Чтобы получить вероятностное представление реше- ния задачи (9)-(14), введём систему СДУ: r (3) zi+1 = zi + yi (α1)i (ψ1)i χ¶G (xi )Di+1K + + yi (α2 )i (ψ2 )i χ¶G (xi )Di+1K , ki+1 = ki + Di+1K , (23) (24) Xr = x + ò σvχG ( Xv )dWv + T -t где i, i + 1 - номера узлов сетки; ξi - вектор трех незавиr r симых N(0,1) случайных величин; Di+1K = d ( X D1 ) - + a χ ( X )dv + n( X )χ ( X )dK , i+ (15) ò v G v ò v ¶G v v D D T -t T -t r расстояние от из области. Xi +1 до ¶G в случае выхода Xi +1 Yr = 1+ r ò α1(v)Yvχ X3 (v)=0dKv + T -t Применение методов Эйлера и случайного блу- ждания по движущимся сферам. Поскольку в работе предполагается, что теплофизические свойства мате- + ò α2 (v)Yvχ X3 (v)=l3 dKv , T -t (16) риалов, из которых состоит сотовая панель, являются постоянными величинами, то случайный процесс Xv Zr = r r ò α1(v)ψ1(v)Yvχ X3 (v)=0dKv + T -t в однородной среде, соответствующей каркасу или ячейке, с точностью до постоянного множителя сов- падает с винеровским процессом. Известно, что точка первого выхода из шара вине- + ò α2 (v)ψ2 (v)Yvχ X3 (v)=l3 dKv , T -t r Kr = ò χ¶G ( Xv )dKv , T -t (17) (18) ровского процесса, начинающего движение из центра этого шара, имеет на границе этого шара равномерное распределение. При этом эта точка не зависит от вре- мени первого выхода на границу. В связи со сказан- ным предлагается для ускорения моделирования трагде x Î G - начальная точка для X ; W (3) - трехмеректорий СДУ (15)-(18) при их прохождении внутри v v ячейки использовать метод случайного блуждания 1 по сферам, а движение по каркасу и некоторой его ный винеровский процесс; sv = (2b(m) ( Xv ))2 ; av - окрестности моделировать с помощью метода Эйлера æ ¶b(m) ö с малым шагом. вектор с координатами ç è ¶xi ÷ , i = 1, 2, 3; n(P) - ø При решении нестационарных задач методом блуждания по сферам требуется каждый раз моделивнутренняя нормаль в точке P Î ¶G ; Kv - локальное время пребывания процесса Xv на границе, т. е. ска- лярный возрастающий процесс, который растет, толь- ко когда Kv Î ¶G . Вероятностным представлением решения задачи ровать не только точку первого выхода из шара, но и время первого выхода. Для моделирования времени первого выхода из шара радиуса R в работе исполь- зуется определение времени первого достижения уровня R бесселевским процессом, который является решением СДУ [8]: (9)-(14) является математическое ожидание t -1 u(x, t) = Ex,T -t (φ( XT )YT + ZT ) , (19) St = S0 + ò Sr 0 dr + Wt . (25) где символом Ex,T -t обозначено математическое Одномерный случайный процесс (19) получается в ожидание относительно вероятностной меры Px,T -t , результате применения формулы Ито к евклидовой норме трехмерного винеровского процесса с начальсоответствующей случайному процессу, стартующеной точкой x Î G . При этом выполняется равенство му в момент времени T-t из точки x, т. е. XT -t = x . S0 = x . Время первого выхода трехмерного винеровского процесса из шара радиуса R и время первого дости- жения случайным процессом уровня R имеют одина- ковое распределение. Известны явные формулы для плотности времени первого достижения бесселевским процессом заданного уровня для начальных точек 2 æ ö3 ç (( )÷ ) τ = ç a ÷ G 2 g ç 3 ÷ è 2 ø где H - гамма-распределенная случайная величина S0 = 0 и S0 > 0 [9]. Но эти формулы представлены 5 2 в виде рядов, содержащих функции Бесселя первого с параметрами , 2 . В ходе вычислений, после того, 3 рода и их положительные нули, и поэтому численно- статистическое моделирование по ним крайне затруд- нительно. С другой стороны, вместо заданного уровня, который требуется достичь, можно подобрать некото- рые кривые, зависящие от t, такие, что получаются достаточно простые и удобные для моделирования формулы для плотности времени первого достижения бесселевским процессом этих кривых. Каждая такая кривая определяется задаваемой положительной сигма- конечной мерой μ . Определим функцию 1 ¥ как получено значение τg , радиус сферы для определения очередной точки положения случайного про- цесса определяется по первой формуле в (29). Выбор параметра a дает возможность определить радиусы моделируемых шаров так, чтобы они не выходили за границы текущей ячейки и чтобы они не были слишком малы. В случае, если расстояние до границы ячейки менее некоторой заданной величины e , вы- числения производятся согласно методу Эйлера. Результаты расчётов. Вычисления проводились для сотовой панели, каркас которой и ограничиваюa r(t, x) º q(t, 0, x) - ò q(t, y, x)μ(dy), 0 (26) щие пластины изготовлены из углепластика, а напол- нителем служит воздух. Основные характеристики где q(t, y, x) - плотность вероятности перехода бессе- левского процесса; a - числовой параметр. Из свойства плотности q следует, что функции q и r удовлетворяют параболическому уравнению панели следующие: общая толщина панели - 0,035 м; толщина каждой из ограничивающих пластин панели - 0,001 м; толщина стенок сотового каркаса - 6 ×10-5 м; сотовыми ячейками являются правильные шести- = - ¶v 1 ¶2v ¶ æ 1 ö è ø ¶t 2 ¶x2 ¶x ç x v ÷. (27) угольники одинакового размера с длиной стороны, равной 0,0042 м; коэффициенты температуропровод- ности углепластика и воздуха равны соответственно Пусть функция g(t) такая, что на кривой x = g(t) выполняется равенство r(t,x) = 0. Тогда эту кривую можно рассматривать в качестве граничного условия к уравнению (27). Причем с учетом того, что бессе- левский процесс стартует из заданной точки, решение полученной краевой задачи определяется единствен- ным образом. Обозначим τg время достижения бесселевским процессом этой подвижной границы, т. е. τg = = inf {t > 0, St = g (t)} . Можно показать [9], что плот- 8, 0 ×10-4 и 2, 36 ×10-5 м2/с. Граничные условия для расчётов были сформиро- ваны на основании физических характеристик, полу- ченных по данным лётного эксперимента в холодном климате и соответствуют первым 150 с взлёта самоле- та. На рис. 1 приведён график коэффициента тепло- обмена на внешней стороне обшивки. На рис. 2 приведены полученные в результате рас- чётов при объёме выборки 2000 траекторий случайно- го процесса графики изменения значений температу- ры в середине панели и у края панели, обращённого ность времени достижения границы летворяет уравнению d æ g (t ) x - g(t) = 0 ö удоввнутрь салона. В табл. 1, 2 приведены для сравнения затраты вре- мени счёта предложенным комбинированным мето- дом и только методом Эйлера с различными шагами pg (t) = - çç ò r(t, x)dx ÷÷ . (28) по каркасу hf и по объёму внутри ячеек ha . dt è 0 ø Результаты расчётов позволяют сделать следую- Если положить μ(dy) = y2dy , то g(t) и деляются по формулам [9] 1 pg (t) опрещие выводы: а) предложенный метод значительно эффективнее использовавшегося ранее метода Эйлера; б) так как при увеличении шага в методе Эйлера æ æ öö2 эффективность комбинированного метода растет ç ç ç 2 ) ÷ g(t) = ç 2t ln ç a ÷÷ , по сравнению с методом Эйлера, то основная часть G ) (( ÷÷ ç 3 2 ÷ è è øø g3(t) вычислительных затрат приходится на движение по тонкому каркасу; в) затраты в средней точке панели меньше, чем pg (t) = . 2at (29) у края, поскольку в середине панели большая часть траектории приходится на движение по внутреннему В этом случае получается простая формула для моделирования τg пространству ячеек, где вычисления проводились с меньшим шагом. Рис. 1. Коэффициент теплообмена на внешней стороне обшивки Fig. 1. Heat transfer coefficient on the outside of the skin Рис. 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 Заключение. В работе предложен численно-стати- стический метод для оценки теплового состояния теплозащитных конструкций сотового типа. Данный метод основан на вероятностном представлении крае- вой задачи для параболического уравнения и комби- нированном применении методов Эйлера и случайно- го блуждания по движущимся сферам. Получено зна- чительное ускорение счёта по сравнению с моделиро- ванием траекторий только одним методом Эйлера.

作者简介

S. Gusev

Institute of Computational Mathematics and Mathematical Geophysics SB RAS; Novosibirsk State Technical University

6, Academika Lavrenjeva Av., Novosibirsk, 630090, Russian Federation; 20, K. Marksa Av., Novosibirsk, 630073, Russian Federation

V. Nikolaev

SibNIA named after S. A. Chaplygin

Email: sag@osmf.sscc.ru
21, Polzunovа Str., Novosibirsk, 630051, Russian Federation

参考

补充文件