О выборе функционала для одной вариационной задачи газовой динамики

осложненных условиях путем применения магнитных устройств // Нефтяное хозяйство. - 1997. - № 7. - С. 46-47.

5. Пивоварова Н.А., Унгер Ф.Г., Туманян Б.П. Влияние обработки постоянным магнитным полем на парамагнитную активность нефтяных систем // Химия и технология топлив и масел.

- 2002. - № 1. - С. 30-32.

6. Лесин В.И., Дюнин А.Г., Хавкин А.Я. Изменение физико-химических свойств водных растворов под влиянием электромагнитного поля // Журнал физической химии. - 1993. - Т. 67. -№ 7. - С. 1561-1562.

7. Борсуцкий З.Р., Ильясов С.Е. Исследования механизма магнитной обработки нефтей на основе результатов лабораторных и промысловых испытаний // Нефтепромысловое дело. - 2002.

- № 8. - С. 28-37.

8. Пат. 2153126 РФ. МКИ5 С01Б 111/22. Устройство для защиты трубопроводов от коррозии / В.М. Кондаков, А.Н. Качуров-

ский, А.Л. Бушковский, В.А. Кольцов, Л.В. Прасс, В.Н. Лялин. Заявлено 28.09.1998; Опубл. 20.07.2000, Бюл. № 20. - 16 с.: 7 ил.

9. Спектроскопия оптического смешения и корреляция фотонов / Под ред. Г. Камминса, Э. Пайка. - М.: Мир, 1978. - 574 с.

10. Посадов И.А., Поконова Ю.В. Структура нефтяных асфальте-нов. - Л.: ЛТИ, 1977. - 76 с.

11. Лоскутова Ю.В., Юдина Н.В. Влияние постоянного магнитного поля на реологические свойства высокопарафинистых неф-тей // Коллоидный журнал. - 2003. - Т. 65. - № 4. -С. 510-515.

12. Лоскутова Ю.В., Юдина Н.В. Реологическое поведение нефтей в магнитном поле // Инженерно-физический журнал. - 2006. - Т. 79. - № 1. - С. 102-110.

13. Лоскутова Ю.В. Влияние магнитного поля на реологические свойства нефтей: Дис. ... канд. хим. наук. - Томск, 2003. -144 с.

УДК 533.6.011

О ВЫБОРЕ ФУНКЦИОНАЛА ДЛЯ ОДНОЙ ВАРИАЦИОННОЙ ЗАДАЧИ ГАЗОВОЙ ДИНАМИКИ

В.М. Галкин

Томский политехнический университет E-mail: vlg@tpu.ru

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

1. Введение

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

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

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

сверхзвуковым и возможно использование простых в реализации и быстрых в расчете маршевых схем; второй подход допускает существование дозвуковых течений и, соответственно, требуется применение гораздо более сложных и более медленных численных методов. Пусть в основу первой численной реализации положен метод характеристик [2]. Тогда очевидно, что если в некоторых контурах не будет реализовываться полностью сверхзвуковое течение, то это приведет к аварийному останову (авосту). Следствием этого будет неодносвязность области определения, поскольку функционал, использующий параметры потока на выходе из сопла, вычислить нельзя. Второй подход, учитывающий появление дозвукового течения и использующий, например, метод установления и схему Годунова [3], лишен это недостатка, однако интеллектуальные и временные затраты по сравнению с методом характеристик возрастают на порядок.

В работе [4] наряду с методом характеристик, позволяющим быстро получать решение, предложен функционал, который, в частности, использует значения параметров течения, найденные на каждой характеристике С+. Это позволяет даже при наличии авостов продолжить функционал на одно-связную область.

В данной статье, которая является развитием работы [4], предлагается вычислять функционал

более простым способом. Для оценки правильности полученных результатов использовалось решение, полученное по методике [2].

2. Постановка задачи

Дано установившееся незакрученное изоэнтро-пическое и изоэнергетическое течение идеального совершенного газа в осесимметричном сопле, рис. 1. Характеристические уравнения и условия совместности имеют вид:

^ = ^0±а),

ах

а0±-

С082

а 5"1

где у - показатель адиабаты, далее 7=1,4; х и у -продольная и поперечная координаты, отнесенные к радиусу сопла в минимальном сечении; а=агап(1/М) - угол Маха, М - число Маха; 0 - угол наклона вектора скорости (линии тока) к оси х; знак +(-) соответствует характеристике С+(С~); полагается, что во входном сечении ха=0. На стенке реализуется условие непротекания: tg(0)=/'(x), где /(х) - функция, описывающая контур сопла, а штрих обозначает производную по х. Во входном сечении ха поток однородный: 0п=0, Мп=1, где индекс «ш» соответствует входу в сопло. Далее индекс «0» соответствует выходу из сопла.

Рис. 1. Схема сопла: оа - минимальное сечение; аЬ - сопло с угловой точкой; аё - характеристика С-, принадлежащая пучку волн разрежения; ае - начальная характеристика ; к1 - характеристика С+ пт - характеристика С++; сЬ - конечная характеристика С

При вычислении параметров течения в угловой точке «а» используется формула Прандтля-Май-ера, следующая из (1):

62 +^(«2) = д1 +

у/(а) - -а + агС^

у — 1

у +1

tg

где нижние индексы 1 и 2 соответствуют параметрам до и после разворота в угловой точке; 01=0, а=ап.

Необходимо найти контур сопла (рис. 1), на выходе из которого в сечении х=хь поток должен иметь однородные поперек сопла параметры 0=0 и число Маха М0>1. Считается, что давление на выходе из сопла больше давления окружающей среды.

3. Тестовая задача

При заданном показателе адиабаты рассмотрим однопараметрическое, зависящее от М0, семейство сопел с угловой точкой и равномерной характеристикой на выходе. Известно [1], что эти сопла имеют максимальную тягу, являются кратчайшими и обеспечивают нулевые потери на рассеяние в выходном сечении. Поскольку на выходе поток будет однородный и параллельный, то указанные сопла будут удовлетворять поставленным выше условиям.

Использование методики [2] при заданных М0>1 и показателе адиабаты позволяет найти единственное сопло, принадлежащее указанному семейству и удовлетворяющее условиям задачи. Контур сопла, полученный таким образом, считается эталоном, и с ним будет проводиться сравнение прямого метода, а координаты начальной точки «а» и найденной конечной точки «Ь» будут являться исходными данными для прямого метода.

4. Прямой метод

Даны координаты точек «а», «Ь» и тангенс угла наклона сопла в точке «Ь»:

/(ха ) = Уа, /(хЬ ) = УЬ, /'(хЬ ) = У^

ха = 0, Уа = 1, У' = 0.

Последнее равенство следует из условия 0=0 на выходе из сопла.

Для аппроксимации искомого контура в качестве базисных функций используются степенные полиномы:

N+3

/(Х) = X 1 = (2Х - ХЪ - Ха )/(X - Ха ) ,

1=1

г е [-1,1], Х е[Ха, Хъ].

Так как должны выполняться условия (3), то коэффициенты ст, Сдт+2, ст и профиль сопла/(х) выражаются через N линейно независимых коэффициентов с1, ..., Сдт. Требуется найти профиль сопла, удовлетворяющий геометрическим условиям (3) и постановке задачи.

Прямой расчет. В рамках прямого метода принято называть прямым расчетом единичный расчет поля течения, который заканчивается вычислением функционала. Выбор функционала будет обсуждаться ниже. В прямом расчете система уравнений (1) решается по схеме, рис. 2, где цифрами «1» и «2» обозначены точки с известными параметрами, из которых выходят характеристики С + и С-. Цифрой «3» обозначена точка их пересечения. Записывая уравнения (1) в разностном виде, получим систему уравнений относительно неизвестных параметров а/, 0/, х/, у3 в точке «3»:

У3 - У 2

Уз - У1

= ^ё(023 -«2зХ

= 1Б(01з +«13),

соэ2а

(у +1)12 - соэ2 а

-(а]3 -а1) +

зта13 зт013

У13 СО§

соэ а

(у +1)/ 2 - соэ2 а

(^ -*2) = 0,

*) = о,

К -а2) -

э1па23 зт023

У23 СО§

здесь 7=1,2,... - номер итерации. Обозначим р=а,в,х,у, тогда Р1з=(Р1+Рзу-1)/2, Р2з=(Р2+Рз-1)/2. Полученная система решалась итерационно до выполнения условия тахр-рз^^О-8. На начальной итерации полагалось р30=(р1+р2)/2. Так как на оси 0=0, у=0, а на стенке 0=аг^(Дх)) и у=Цх), то в данной системе уравнений производились очевидные упрощения.

Рис. 2. Схема расчета: 13 - отрезок характеристики С+;

23 - отрезок характеристики С

Таким образом, от оси к стенке, по известной характеристике С+ рассчитывается следующая характеристика С ++1, где индексы I и /+1 - порядковые номера вычисляемых характеристик. Расчет производится до тех пор, пока очередная характеристика С не приходит в точку «Ь». В качестве начальных условий используются 0Ы=О и М=1,001 на начальной характеристике ае. Такие начальные условия достаточно часто применяют вместо плоской звуковой линии [5]. Отметим, что для вычисления интегралов использовалась формула трапеции.

Выбор функционала. В [4] использовался функционал, который рассчитывался вдоль характеристики (рис. 1):

I=л .

Если при вычислении характеристики С++1 происходил авост, то формула (5) заменялась выражением:

Г/ ь

I = (хь - х,) + П02й/ + |ф(х)йх,

У к а

ОI(х)|, I'(х) < о

ф( х) =

О, I'(х) > О

где первый интеграл вычислялся вдоль последней рассчитанной характеристики С+, а второй - суммировал площадь участков контура сопла с отрицательными углами наклона. Если авост произошел при расчете характеристики С+, то х=0 и 0=О вдоль начальной характеристики С2+. В результате функционал (6) приобретает вид:

ь

I = хь + |ф( х )йх.

Как видно, способ вычисления функционала (5-7), предложенный в [4], достаточно сложен. И, кроме этого, замена формулы (5) на формулу (6) приводит к тому, что функционал, сохраняя непрерывность, в общем случае может стать недиффе-ренцируемым.

В данной работе предлагается другой подход. Рассмотрим два функционала

11 =

12 = Ь1 ((а - «О )/а )2й/, Ь =| й/,

\ С+ С+

которые вычисляются вдоль произвольной характеристики С+. Выбор функционалов (8) обусловлен следующими соображениями.

Пусть сопло с известным Мо удовлетворяет условиям задачи. Тогда, как видно на рис. 3, первый функционал из (8) имеет два минимума 11=0. Первому минимуму на характеристике ае соответствует 0=О и ай=агс8т(1/Мя). Второму минимуму на характеристике сЬ соответствует 0=О и аО=аг8т(1/М0). Как следует из уравнений (1), только второй минимум является необходимым и достаточным условием решения задачи. Очевидно, что использование первого функционала не гарантирует единственности решения.

Рис. 3. Значение функционалов вдоль характеристик С+, / - номер характеристики

Второй функционал из (8), в отличие от первого, имеет единственный минимум 12=0 на характеристике сЬ (рис. 3), чему соответствует а=а0. Однако этот минимум является необходимым, но не достаточным условием решения задачи, поскольку может не выполняться условие 0=О. Комбинация I и 12 дает следующий функционал:

а

с

«о)/ а )2+02] dl.

который имеет единственный минимум /3=0 на характеристике сЬ (рис. 3). Так как этот минимум является необходимым и достаточным условием решения задачи и ему соответствует 6=0 и а=а0, то далее будет использоваться функционал (9).

При вычислении функционала (9) воспользуемся приемом, предложенным в [4], и если при вычислении характеристики С+ происходит авост, то функционал (9) вычисляется вдоль уже рассчитанной характеристики С+.

Если сопло имеет отрицательный наклон сопла Дха)<0 в точке «а», то авост происходит при вычислении характеристики С+, и функционал (9) вычисляется на начальной характеристике ае. Поскольку при вычислении параметров течения в угловой точке «а» используется формула Прандтля-Майера (2), то можно использовать эту формулу и для вычисления а2 при 62<0. В этом случае найденные значения а2 не имеют физического смысла, т. к. этому соответствует дозвуковое течение. Однако в силу своей непрерывности (рис. 4) функция а2(62) позволяет продолжить функционал (9) на односвя-зную область в случае отрицательного наклона сопла Дха)<0 и обеспечить функционалу непрерывность и дифференцируемость. Как показали расчеты, при вычислении по формуле трапеции функционала на начальной характеристике ае достаточно брать только две точки - «а» и «е».

Зависимость а2=а2

Поскольку на входе в сопло граничные условия не меняются, а параметры течения определяются только профилем сопла (4), то значение функционала (9) будет неявно зависеть от этого профиля. Таким образом, задача нахождения профиля сопла У(х), доставляющего экстремум функционалу (9), сводится к поиску точки (с1,^,сЛГ), в которой функция многих переменных /3=/3(с1,^,сж) имеет экстремум. Для нахождения минимума этой функции использовался квазиньютоновский метод Бройдена из [6].

5. Численные результаты

Число N варьируемых переменных с менялось от 1 до 10, начальное значение с=0. На характери-

стике ае задавалось 50 точек. Координаты точек «а» и «Ь»: (0; 1) и (3,576; 1,299). Этим значениям соответствует найденное из [2] эталонное сопло с угловой точкой и числом Маха на выходе М0=2.

Сравнение решения для обоих функционалов с эталонным соплом показало, что при N=10 ординаты сопел отличались в четвертом знаке после запятой, а максимальная относительная погрешность по составила 0,02 %.

На рис. 5 представлены результаты минимизации в виде окончательных значений функционалов для различного числа N. Таблица демонстрирует влияние числа коэффициентов и используемого функционала как на число прямых расчетов, так и на количество авостов. Из представленных результатов видно, что предложенный функционал (9), несмотря на простоту, по эффективности не уступает функционалу (5-7), предложенному в [4].

Рис. 5. Минимизация функционалов: 1) формулы

Таблица. Минимизация функционалов. В числителе - число прямых расчетов, в знаменателе - авостов

№ Формулы

1 13/2 17/2

2 36/3 39/4

3 63/4 69/5

4 98/4 99/5

5 135/5 115/7

6 150/7 166/7

7 199/7 203/6

8 273/6 304/7

9 385/8 320/8

10 444/7 492/9

Заключение

Проведенные численные исследования показывают, что предложенный функционал достаточно прост в реализации, при наличии авостов имеет продолжение на односвязную область и сохраняет дифференцируемость, а по эффективности не уступает ранее рассмотренному в работе [4] функционалу. Сравнение профиля найденного сопла с эталонным профилем показало, что при 10 варьируемых переменных максимальная относительная погрешность по ординате составляет 0,02 %.

СПИСОК ЛИТЕРАТУРЫ

1. Крайко А.Н. Вариационные задачи газовой динамики. - М.: Наука, 1979. - 448 с.

2. Кацкова О.Н. Расчет равновесных течений газа в сверхзвуковых соплах. - М.: ВЦ АН СССР, 1964. - 61 с.

3. Годунов С.К., Забродин А.В., Иванов М.Я. и др. Численное решение многомерных задач газовой динамики. - М.: Наука, 1976. - 400 с.

4. Галкин В.М., Волков Ю.С. Сравнение базисных функций в прямой задаче профилирования сверхзвуковой части сопла //

Сибирский журнал индустриальной математики. - 2004. - Т. 7. - № 4(20). - С. 48-58.

5. Крайко А.Н., Тилляева Н.И. Оптимальное профилирование контура сверхзвуковой части тарельчатого сопла // Механика жидкости и газа. - 2000. - № 6. - С. 172-184.

6. Дэннис Дж. мл., Шнабель Р. Численные методы безусловной оптимизации и решения нелинейных уравнений. - М.: Мир, 1988. - 440 с.

УДК 536.21+692.2:691.327:666.973.2:666.64-492.3

НЕСТАЦИОНАРНЫЙ ПРОСТРАНСТВЕННЫЙ ТЕПЛОПЕРЕНОС В НЕОДНОРОДНОЙ КЕРАМЗИТОБЕТОННОЙ СТЕНЕ

А.Н. Хуторной, А.Я. Кузин, Н.А. Цветков, Т.А. Мирошниченко, А.В. Колесникова

Томский государственный архитектурно-строительный университет E-mail: kaftgs@tsuab.ru

С помощью математического моделирования исследовано тепловое состояние неоднородного трехслойного керамзитобетон-ного фрагмента стены с фасадным утеплением на гибкой связи. Установлен характер распределения полей температуры в ке-рамзитобетонной конструкции. Выполнена оценка зоны влияния гибкой связи на температурное поле стены. Разработанная численная методика позволяет проводить тепловую экспресс-диагностику наружных неоднородных керамзитобетонных стен с различными теплофизическими и геометрическими характеристиками материалов системы в реальных условиях эксплуатации.

Вопросы энергосбережения относятся к приоритетным направлениям развития науки, технологий и техники Российской Федерации. Очевидно, что решение этих вопросов должно иметь комплексный характер. В работе [1] выполнена оценка энергоэффективности комплекса энергосберегающих мероприятий в нескольких зданиях общественного назначения г. Москвы. Суммарная экономия энергии оказалась весьма значительной (56...63 %), при этом доля утепления несветопро-зрачных ограждений составила 16.20 %. Там же [1] подчеркнута важность энергосберегающих мероприятий, связанных с утеплением несветопрозрач-ных ограждений, как обеспечивающих суммарное снижение энергопотребления здания не менее чем в 2 раза. В этой связи разработка способов повышения теплоэффективности наружных стен существующих зданий и разработка новых конструкций с повышенными теплозащитными свойствами является объективной необходимостью.

Повысить уровень теплозащиты наружных стен зданий можно за счет устройства в толще конструкции пустот, заполненных эффективным утеплителем [2]. Однако, как показано в работе [3], для ряда климатических зон такое решение может быть не достаточным. С целью обеспечения требованиям СНиП 11-3-79* «Строительная теплотехника» возникает необходимость в устройстве дополнительного утепления ограждающих конструкций, которое может быть реализовано с помощью фасадных систем утепления на гибких связях.

Физико-математическая постановка задачи. Рассмотрим нестационарный теплоперенос через трехслойную неоднородную конструкцию, состоящую из керамзитобетона - 1 с вертикальной пустотой, заполненной утеплителем - 4, утеплителя фасадной системы утепления - 2 и обшивки - 3 (рис. 1). Керамзитобетон с обшивкой соединены гибкой связью (коннектором) - 5. Форма керамзитобето-на, внутреннего утеплителя, коннектора, утеплителя фасадной системы утепления и обшивки - прямые параллелепипеды, поперечные и продольные сечения которых в общем случае представляют собой разносторонние прямоугольники. Известны те-плофизические характеристики (А;, р, с,, /=1,5) материалов системы, ее геометрические размеры, температуры наружной (^£) и внутренней (г1^) сред, коэффициенты теплоотдачи на наружной (о,) и внутренней (а0) поверхностях ограждения [4].

Необходимо рассчитать поля температур и плотностей тепловых потоков в сечениях неоднородной керамзитобетонной стены с фасадным утеплением.

Теплоперенос в неоднородном керамзитобетон-ном фрагменте стены описывается в декартовой системе координат трехмерными нелинейными нестационарными уравнениями теплопроводности

i = 1,5.