ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 3, с. 481-489
УДК 519.634
ПРИМЕНЕНИЕ СОСТАВНЫХ РАЗНОСТНЫХ СХЕМ ДЛЯ РАСЧЕТА НЕСТАЦИОНАРНЫХ ТЕЧЕНИЙ С УЗКИМИ ТЕПЛОВЫМИ СЛОЯМИ1)
© 2007 г. Г. В. Беляков*, В. И. Грынь**, А. А. Чарахчьян**
(* 119334 Москва, Ленинский пр-т, 38, ИДГ РАН;
** 119991 Москва, ул. Вавилова, 40, ВЦ РАН) e-mail: chara@ccas.ru Поступила в редакцию 12.05.2006 г.
Описана методика расчета нестационарных гидродинамических течений с узкими тепловыми слоями, которая применяется к моделированию высокоскоростного нагрева трением при ударном сжатии пластины на клине. Библ. 21. Фиг. 6.
Ключевые слова: нестационарные гидродинамические течения, ударные волны, трение.
ВВЕДЕНИЕ
Для расчета гидродинамических течений с узкими погранслоями различного происхождения и ударными волнами в основной части течения необходимо использовать разностные сетки с сильным сгущением в области погранслоя. В случае стационарных течений такого типа высокую эффективность показали квазимонотонные неявные схемы переменных направлений (см., например, [1], [2]). Однако возможность эффективного применения таких схем к расчету нестационарных течений далеко не очевидна. В самом деле, в погранслое течение, как правило, дозвуковое, что в принципе дает возможность проводить расчет с шагом по времени
• K
т > тк = minii; j,,
и j
где TKj определяется условием Куранта для ячейки сетки (i, j) и зависит от размеров ячейки и местной скорости звука. В настоящей работе мы ограничиваемся случаем двумерных течений. Далее сошлемся на работу [3], где изучалось применение схем переменных направлений к расчету нестационарной задачи конвекции, поставленной для уравнений Навье-Стокса сжимаемого газа. Было показано, что нарушение условия Куранта, по крайней мере для обычных двухслойных схем расщепления, приводит к значительной потере точности. Так как устойчивость расчета при этом сохраняется, диагностировать такую потерю точности можно только сравнением расчетов с существенно разными шагами по времени. Причиной потери точности, как было показано в [4], является рост той части погрешности аппроксимации, которая вносится в схему процедурой расщепления. Переход к неявной схеме без расщепления позволил в [4] значительно увеличить шаг по времени. Хотя в [3], [4] применялись равномерные прямоугольные сетки, а течение везде было дозвуковым, нет оснований надеяться, что эффект понижения точности схем переменных направлений исчезнет в случае неравномерных сеток.
Другая трудность создания эффективных методов расчета нестационарных течений с ударными волнами и узкими пограслоями заключается в том, что отказ от расщепления по направлениям в известных квазимонотонных схемах приводит к очень неэкономичным схемам, так как на каждом временном слое приходится решать систему алгебраических уравнений относительно значений нескольких разностных функций.
В [5], [6] был предложен новый принцип объединения двух разных схем в одну составную схему, свободную от рассмотренных выше недостатков. Такие схемы в погранслое переходят в неявную схему, в основной части течения - в явную квазимонотонную схему, а плавность перехода от одной схемы к другой определяется плавностью перехода от мелкой сетки к крупной. Экономичность составной схемы определяется, во-первых, тем, что отсутствие требования квазимоно-
1) Работа выполнена при финансовой поддержке РФФИ (код проекта 04-01-00051) и Российской академии наук (программа ОМН РАН "Соврем. вычисл. и информ. технологии решения больших задач").
тонности позволяет построить неявную схему с системой алгебраических уравнений относительно значений только одной скалярной функции. Во-вторых, число уравнений примерно равно числу узлов сетки в погранслое, что значительно меньше общего числа узлов. Это позволяет построить эффективные составные схемы и для параболических уравнений, возникающих при использовании расщепления по физическим процессам на этапах учета вязкости и теплопроводности, объединяя схему переменных направлений, "работающую" в основной части течения, и чисто неявную схему, "работающую" в погранслое.
В [5], [6] составные схемы применяются для расчета вязкого нагрева дейтерия, сжимаемого внутри конической мишени. Настоящая работа посвящена применению составных схем для расчета нагрева трением, моделируемого заданием потока тепла на подвижной границе раздела сред. Так как скорость тепловой волны много меньше скорости вещества и скорости звука, тепловая волна имеет вид узкого погранслоя вблизи границы раздела.
1. ЧИСЛЕННЫЙ МЕТОД
Рассматриваются плоские течения сжимаемой теплопроводной жидкости, которые описываются уравнениями
Ш 1 1 , • ди д иу
-— = - ШУи, ШУи = -=-— + -=--- ,
р дх ду
(1)
дих = 1 др йи = 1 дР д t рдх' Ш рду'
^ = -Р&у и + &у кУТ, ш р
где t - время, х, у - пространственные координаты, и - вектор скорости с координатами их, иу, Ш/й = = д^ + ихд/дх + иуд/ду - лагранжева производная по времени, р - плотность, р - давление, £ -удельная внутренняя энергия, Г - температура, к = к(р, Г) - коэффициент теплопроводности. Уравнения (1) замыкаются уравнениями состояния р = р(р, Г), £ = £(р, Г).
В плоскости х, у имеется подвижная регулярная четырехугольная невырожденная сетка О^), некоторые линии которой являются лагранжевыми, т.е. движутся вместе с веществом, а внутренние узлы сетки, не принадлежащие выделенным линиям, вычисляются тем или иным методом построения сеток (см. [7]). Разностные функции определены в серединах ячеек сетки. Используется расщепление на лагранжев этап и этап пересчета поля течения с лагранжевой сетки О-ф на сетку О^ + т) (см. [8]).
Рассмотрим лагранжев этап. Множество всех ячеек сетки разделим на два подмножества, которые назовем множеством крупных Ц. и мелких Ц ячеек. Например, если сетка сгущается к некоторой линии, к множеству Ц можно отнести некоторое количество слоев сетки вблизи этой границы, а к множеству Ц. - остальные ячейки сетки. Шаг по времени определим формулой
т = шт тк. (2)
г, 1
Переход на верхний временной слой состоит из двух этапов. На первом этапе шаг по времени
Гт, /, 1 еЦ -,
<1 = \ к (3)
К-, 1, г, 1 е Ц,
свой для каждой ячейки сетки. На втором этапе шаг по времени
(2) Г 0, г, 1 еЦ,
Т" = \ к . . Ц (4)
[т - Я1, г, 1 е Ц,
дополняет шаг (3) до полного шага т. Шаг по времени первого этапа допускает использование явной схемы для уравнений гидродинамики. На втором этапе предполагается использование неявной схемы. Как видно из (3), (4), для крупных ячеек из множества Ц. переход на верхний вре-
к
менной слой выполняется по схеме первого этапа, для очень мелких ячеек с тг , < т основной
вклад в переход вносит неявная схема второго этапа, а плавность перехода от одной схемы к другой определяется плавностью перехода от мелкой сетки к крупной.
На первом этапе используется квазимонотонная схема типа С.К. Годунова второго порядка точности, которая применялась ранее для расчета различных течений с ударными волнами (см., например, [9]) в сочетании с неявной схемой для уравнения теплопроводности. Последняя схема, в свою очередь, строится как составная, объединяющая схему переменных направлений и схему без расщепления по направлениям (см. [5], [6]).
На втором этапе строится неявная схема без расщепления как по направлениям, так и по физическим процессам, обобщающая схему [10], [6] на уравнения (1). Применяется метод Ньютона, на каждой итерации которого возникает система линейных уравнений относительно значений двух разностных функций, которая, как и в [5], [6], решается прямым методом из [11], использующим разложение ленточной матрицы на треугольные множители. Так как множество Ц^, на котором строится неявная схема, сравнительно невелико, вычислительные затраты на ее реализацию не влияют заметно на общую эффективность метода.
Схема лагранжевого этапа тестировалась на одномерной задаче о нагреве сжимаемой жидкости постоянным потоком тепла (см. [12]). В уравнениях (1) все неизвестные зависят только от х и иу = 0. В начальный момент t = 0 область х > 0 заполнена неподвижным алюминием с давлением 60 ГПа и плотностью 3.5 г/см3, что примерно соответствует параметрам алюминия за фронтом присоединенной ударной волны при бесструйном сжатии пластины со скоростью 4 км/с (см. [13]). Соответствующая двумерная задача будет рассмотрена в следующем разделе. На границе х = 0 ставится условие симметрии их = 0 и задается постоянный поток тепла q = -кЭГ/Эх = = 30 ГВт/см2. Используется широкодиапазонное уравнение состояния алюминия (см. [14]) в табличной форме. Для коэффициента теплопроводности используется та же модель, что и в [12].
После небольшого по времени начального этапа, когда формируется ударная волна, решение задачи имеет вид двух волн: ударной волны, после которой давление повышается до некоторого значения, и тепловой волны, в которой давление не меняется и скорость распространения которой много меньше скорости ударной волны. Расчет велся в переменных Лагранжа. При t = 0 задавалась сетка по х, которая являлась комбинацией сеток с постоянным и переменным шагом, меняющимся по геометрической или арифметической прогрессии. Результаты расчета на разных сетках в момент времени, когда плотность при х = 0 уменьшилась до некоторого заданного значения, показаны на фиг. 1. Для всех трех сеток шаг вблизи границы х = 0, Н1 = 2 х 10-5 мм достаточно мал, чтобы можно было разрешить структуру тепловой волны. Расчет на самой подробной сетке велся без применения составной схемы, т.е. множество ячеек полагалось пустым. Расчет на грубых сетках велся по описанной выше составной схеме с шагом по времени т - 100тк. Видно, что расчеты на грубых сетках правильно воспроизводят как структуру тепловой волны, так и положение ударной волны.
Фиг. 1. Давление вблизи ударной волны (график (а)) и температура вблизи тепловой волны (график (б)) в одномерной задаче о нагреве постоянным потоком тепла на разных сетках: сплошная линия - шаг перед ударной волной к - 0.0005, точки - шаг к - 0.02, треугольники - шаг к - 0.05 мм.
О", 3 - 1)
и ■
з + 1,3
и + 1, з
(3 + 1,3)
1,3
■ 3 (3 - 1
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.