Некоторые вопросы численного моделирования нестационарных газодинамических течений тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

РГб од

РОССИЙСКАЯ АКАДЕМИЯ НАУК ВЫЧИСЛИТЕЛЬНЫЙ ЦЕНТР

на правах рукописи УДК 519.6:537.84

ЧАРАХЧЬЯН АЛЕКСАНДР АГАСИЕВИЧ

НЕКОТОРЫЕ ВОПРОСЫ ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ НЕСТАЦИОНАРНЫХ ГАЗОДИНАМИЧЕСКИХ ТЕЧЕНИЙ

01.01.07 - вычислительная математика 01.02.05 - механика жидкости, газа и плазмы

Автореферат диссертации на соискание ученой степени доктора физико-математических наук

Москва 1993

Работа выполнена в Вычислительной центре РАН

Официальные оппоненты:

доктор физико-математических паук, профессор ТИ1КИН В.Ф.,

доктор физико-математических наук, профессор ФАВОРСКИЙ А.П.,

доктор физико-математических наук, профессор ХОЛОДОВ A.C.

Ведущая организация: Институт прикладной математики РАН им. М.В.Келдыша.

Защита состоится << 6 >> ^ часов на заседании специализированного совета Д 002.32.01 при Вычислительном центре РАН по адресу: 117967, г. Москва, ул. Вавилова, 40..

С диссертацией можно ознакомиться в библиотеке МИ РАН и ВЦ РАБ по адресу: ул. Вавилова, 42.

Автореферат разослан ■. 2 >> _ 199 года.

Ученый секретарь специализированного совета доктор физ.-ыат. наук

Е.Д.Терентьев

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

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

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

В вычислительной гидродинамике имеется три направления, в рамках которых развивались метода численного решения задач с областью течения из нескольких сред. К первому направлению относятся методы сквозного счета через границу раздела сред, начало которым положил известный метод частиц в ячейке Харлоу. Здесь используется либо техника смешанных ячеек, состоящих из нескольких сред (метод частиц в ячейке; метод концентраций,

- ? -

Бахрах С.М. и др., 1981), либо специальная процедура переинтерполяции данных с лагранжевых частиц на эйлерову сетку и обратно (метод индивидуальных частиц, Агурейкин В.А., Криков Б.П., 1986). Другое направление образуют активно развивающиеся в последнее время свободно-лагрвнжевы методы, когда шаблон разностной схемы определяется на кавдом шаге по времени в зависимости от расположения узлов лагранжевой сетки. Несомненным достоинством методов из этих двух направлений является их универсальность, т.е. возможность сравнительно быстрого переключения соответствующих пакетов программ с одной задачи на другую. Вместе с тем естественным продолжением этого достоинства является следующий недостаток. В методах первого направления форма границы раздела сред вообще не определяется, а в методах второго направления граница раздела сред как правило представляется ломаной линией. В то же время границы раздела сред в решении уравнений Эйлера являются гладкими кривыми. В результате использование методов из первых двух направлений для расчета течений с интенсивным перемещением одной из сред вдоль границы раздела приводит к "схемному1? трению, т.е. к торможению одной среда и возмущению в другой, которых нет в решении уравнений Эйлера. Поэтому для получения надежных результатов эти методы требуют очень подробных сеток в окрестности границ раздела.

Настоящая диссертация посвящена развитию третьего из указанных выше направлений, основы которого заложены известной монографией О.К.Годунова и др. (1976). Граница раздела сред в каждый момент времени представляется гладкой кривой, вдоль которой расставляются узлы некоторой линии сетки. "Схемного" тре-

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

В настоящее время для расчета эволюционных течений используются как консервативные, так и неконсервативные схемы. За исключением полностью консервативных схем, о которых особо будет сказано ниже, расчет внутренней энергии е в консервативных схемах происходит следующим образом. Из разностного уравнения энергии определяется рвзностная функция полной энергии е=е+|и|2/2, где вектор скорости и определяется из разностного

уравнения движения, а значение е определяется соответствующим вычитанием. Неконсервативные схемы строятся на базе недивергентного уравнения энергии, из которого исключена кинетическая энергия. В этом случае значение е определяется непосредственно из разностного уравнения энергии. Поскольку разрывные решения уравнений газовой динамики определяются законами сохранения, естественно предположить, что разностное решение, полученное по неконсервативной схеме с невыполнящимися точно законами сохранения, при наличии разрывов будет сходиться к функциям, отличным от точного решения. В диссертации приведен такой пример, связанный с неконсервативным вариантом схемы Годунова, построенном на недивергентной форме уравнения энергии. Аналогичные примеры есть и в работах других авторов. Тем не менее анализ литературы, как отечественной, так и зарубекной, показывает, что, несмотря на наличие указанной выше методической ошибки, неконсервативные схемы отнюдь не являются менее распространенными по сравнению с консервативными, по крайней мере при расчете эволюционных течений. Причина заключается в потенциальной возможности потери точности расчета внутренней энергии в консервативных схемах. Такой недостаток отсутствует в полностью консервативных схемах (Самарский A.A., Попов Ю.П.; Гольдин В.Я., Ионкин Н.И., Калиткин H.H.), полученных путем выделения консервативных схем из множества схем, аппроксимирующих некоторым естественным образом недивергентное уравнение энергии. Эти схемы построены на сетках с разнесенными узлами, когда термодинамические переменные и вектор скорости определены в разных узлах сетки. В то же время имеетоя много схем, в частности схе-

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

По сравнении с невязкими течениями, расчет вязких течений представляет из себя значительно более сложную проблему вычислительной гидродинамики. Эта сложность связана не столько с появлением дополнительных слагаемых в уравнениях, сколько с появлением нового масштаба - толщины погранслоя, где вязкое течение существенно отличается от Иевязкого. В упоминавшейся выше задаче о сжатии дейтерия в конических мишенях схлбпываю-щийся на оси симметрии алюминий гененировал высокоскоростные струи дейтерия. Необходимо было проверить гипотезу о значительном увеличении температуры плазмы вдоль границы струи за счет вязкого нагрэва. Линейныэ размеры области дейтерия при сжатии уменьшались в 5+10 раз, а толщина погрвнслоя в дейтерии была еще в 50+100 раз меньше. Наряду с дозвуковым течением в погрзн-слое, в основной части течения возникали сильные ударные волны, поровденные высокоскоростными струями. Такая задача, даже оставаясь двумерной, по своей сложности значительно превосходит задачи, когда-либо рассмотренные как в отечественной так и в зарубежной литературе по вычислительной гидродинамике.

Значительный прогресс в области численного моделирования вязких течений был достигнут в связи с необходимостью расчета вязкого нагрева при входе тел в атмосферу планет. Из отечественных работ отметим работы О.М.Белоцерковского, Ю.П.Головачева, В.М.Ковени, ¿.И.Толстых, а также работы М.Я.Иванова, В.Г.Нрупы и Р.3.Нигматуллина по расчету внутренних вязких течений. Здесь, по-видимому, впервые в реальных многомерных задачах стали использоватся сетки с сильным сгущением, необходимым для расчета течения в погранслое. Тем не менее применение развитых здесь методов для расчета упомянутой выше задачи о сжатии дейтерия в конической мишени связано с серьезными трудностями. И дело не столько в том, что в этой задаче граница области вязкого течения подвижна, а используемые сетки неортогональны. Основная трудность заключается в следующем.

Обозначим через шаг по времени, который следует из условия Куранта для ячейки сетки (I,./) и определяется размером ячейки и местной скоростью звука. При наличии погранслоев шаг по времени т°=т1п(т;°^) оказывается очень маленьким из-за небольших размеров ячеек в погранслое. К - счастью течение в погранслое существенно дозвуковое, и поэтому скорость изменения искомых функций там много меньше скорости звука, что делает возможным расчет эволюционного течения с шагом по времени т»т°. В то же время задача о входе тел в атмосферу планет, по крайней мере в связи с расчетом вязкого нагрева, является стационарной задачей. Поэтому развитые здесь методы являются неявными схемами переменных направлений,-которые очень эффективны для расчета стационарных течений методом установления. Однако использование

таких схем для расчета эволюционных течений с шагом по времени т>т° может приводить к существенной потере точности. Это было наглядно продемонстрировано в работе Г.М.Махвиладзе и С.В.Щербака (препринт *113, ИПМехан АН СССР, 1978) на примере задачи конвекции. В этой задаче течение везде дозвуковое, что позволило авторам получить надежные результаты на равномерных сетках. Было показано, что точность расчета эволюционной задачи по некоторой двухслойной схеме переменных направлений заметно падает уже при тм4т° при отстутствии каких-либо признаков нарушения устойчивости расчета. Авторами была Предложена трехслойная схема переменных направлений, для которой эффект понижения точности отсутствовал вплоть до т«8т°. Тем не менее при т«1бг° точность расчета опять заметно падала. Заметим, что для прикладных задач, решаемых на пределе возможностей ЭВМ, диагностика такой потери точности является серьезной проблемой.

Причина потери точности заключается в следующем. Дозвуковое течение является одновременно и слабосжймаемым, т.е. в уравнении неразрывности р-1 (а/5£+1Г7)р=-<Цлт левая часть много меньше чем" каждая из пространственных производных, входящих в сНти. В результате при т»т° погрешность аппроксимации в основном определяется той своей частью, которая вносится в схему процедурой расщепления по направлениям. Поэтому схемы для дозвуковых течений следует строить аналогично схемам для несжимаемой жидкости: оба слагаемых, входящих в сИти следует аппроксимировать на одном и том же временном слое. Хотя на каждом временном слое приходится решать двумерную систему алгебраических уравнений, допустимый шаг по времени практически перестает за-

- 8 - .

висеть от условия Куранта. Это было продемонстрировано в работе В.А.Гончарова, В.М.Кривцова и автора (1988), где построенная таким образом схема позволила проводить расчет упомянутой выше задачи конвекции с шагом по времени т«64т° без заметного падения точности.

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

Весьма актуальной является рассмотренная в диссертации задача численного исследования сжатия дейтерия ударниками в конической мишени при наличии сильного кумулятивного аффекта. Результат соответствущего натурного эксперимента до сих пор не получил удовлетворительного объяснения. Неожиданность этого результата в сравнении с другими аналогичными экспериментами заключается в заметном (<*10б) нейтронном выходе при сравнительно небольшой (<«5.4 км/с) скорости ударника. Хотя этот эксперимент не был продолжен после 1980 г. и, по-видимому, не был никем повторен, авторы эксперимента по-прежнему уверены в надежности своего результата, неизменно включая его в свои обзорные статьи, последняя из которых опубликована в 1992 г.

Целью работы является: - Значительно поднять уровень сложности задач с несколькими

средами, доступных для численного моделирования в рамках методики с представлением границ раздела сред гладкими кривыми..

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

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

Научная новизна. Остановимся вначале на результатах, полученных в связи с разработкой новой методики расчета задач с несколькими средами. Методика развивает известный подход к таким задачам, когда граница раздела сред представляется гладкой кривой. Выше уже отмечались трудности, связанные с применением такого подхода к задачам со сложной геометрией. Первая из них заключается в том, что в случае нескольких пересекающихся границ раздела их явное выделение, по крайней мере в случае сеток фиксированной структуры, может оказаться невозможным. Эта трудность была преодолена путем компиляции в единую методику известных алгоритмов. Использование расщепления на лагранжев этап и этап пересчета поля течения с лагранжевой сетки на подвижную эйлерову позволило в рамках одной разностной схемы проводить расчет как с явном выделением границы раздела, так и "сквозным" методом концентраций, в котором граница раздела сред моделируется полосой смешанных ячеек. Это существенно увеличило гибкость методики. Для задач с несколькими границами раздела сред

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

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

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

ности разностной схемы (использовалась схема Годунова). Цри сохранении полной энергии изменение внутренней энергии на лаг-ранжевом этапе в некоторых ячейках сетки оказывается противоположным по знаку тому изменении, которое следует из недивергентной формы уравнения энергии. Отталкиваясь от этого наблюдения, было найдено такое преобразование консервативной схемы в переменных Лагранжа, которое, с одной стороны, давало схему, близкую к консервативной, а, с другой стороны, гарантировало отсутствие потери точности расчета внутренней энергии. Это было под-твервдено тестовыми расчетами. В отличие от неконсервативной схемы, новая схема, так же как и консервативная, сходилась к некоторому разрывному решению уравнений газовой динамики. В то же время эффект ложной кумуляции, в отличие от консервативной схемы, для новой схемы отсутствовал. Найденное преобразование консервативной схемы применимо к широкому классу схем, как двумерных, так и трехмерных, построенных на сетках произвольной структуры, в частности к схеме Годунова и ее модификациям высокого порядка точности - схемам В.П.Колгана, Ван Леера, А.В.Родионова, схеме с параболической интерполяцией (Colella and Woodward, 1984), по крайней мере к двуслойным МО (essentially non-oacillatory) схемам, и др. Не видно никаких причин, по которым это преобразование может ухудшить такие свойства исходной консервативной схемы, как порядок аппроксимации, устойчивость и монотонность. Рассмотрены также схемы без расщепления на лаг-ранжев этап и этап пересчета. Найдэно аналогичное преобразование широкого класса таких консервативных схем на подвижных сетках произвольной структуры.

- 12 -

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

Помимо использованного в методике расщепления на лагранжев этап и этап пересчета, достаточно распространенным является другое близкое расщепление, используемое в методах "крупных частиц" и "жидкости в ячейке". В диссертации выполнен сравнительный анализ обеих расщбплений. На примере схем Харлоу и Годунова показано, что разностные схемы с использованием первого расщепления более устойчивы и более монотонны по сравнению с ' аналогичными схемами, построенными на втором расщеплении.

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

переменных направлений для расчета эволюционных течений не годятся, а формальный отказ от расщепления по направлениям приводит к очень неэкономичным схемам. Развитая в диссертации схема является схемой расщепления по физичесикм процессам: на первом этапе аппроксимируются уравнения Эйлера в переменных Лагранжа, на втором этапе учитываются слагаемые, описывающие вязкость и теплопроводность, а третий этап является алгоритмом пересчета шля течения с лагранжевой сетки на заданную, отвечающую очередному моменту времени. На каждом этапе строится монотонная схема не ниже второго порядка точности по пространственным переменным.

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

Следующий результат является центральным с точки зрения создания схемы, позволяющей вести расчет одновременно и дозвукового течения в погранслое, и сильных ударных волн в основной части течения. Каждый из этих типов течения предъявляет свои, в значительной степени противоположные требования к разностным схемам. Для расчета дозвуковых течений нужна неявная' схема без расщепления по направлениям, позоляющая вести счет с шагом х»10,. При расчете ударных волн одно из основных требований к

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

На защиту выносится также новый монотонный алгоритм пересчета высокого* порядка точности между неблизкими сетками, который используется на последнем этапе расщепления. Алгоритм состоит из экономичной процедуры поиска узлов старой сетки для каждого узла новой и двумерной интерполяционной процедуры на неортогональной сетке. Используются одномерные монотонные интерполяции полиномами Эрмита (De Boor and Swartz, 1977; Fritch. and Carlson, 1980), модифицированные таким образом, чтобы не допускать падения точности в окрестности локальных экстремумов. Алгоритм можно рассматривать в качестве монотонной явной разностной схемы третьего порядка точности для некоторого уравнения переноса, которая имеет повышенную устойчивость (если скорость переноса постоянна, то схема абсолютно устойчива). Выпол-

нэн ¿.я . тестовых расчетов, как одномерных, так и двумерных, в частности в сравнении с известной неявной TVD схемой второго порядка точности. Во всех тестах построенный алгоритм давал более точные результаты.

Остановимся теперь на задаче о сжатии дейтерия в конических мишенях. Подобные задачи в различных постановках рассматривались в ряде работ, однако в большинстве из них кумулятивный эффект не оказывал заметного влияния на течение. Начальная стадия образования кумулятивной струи численно исследовалась в работе В.В.Рассказовой, В.Г.Рогачева, Н.Ф.Свидинской (1985). Задача, рассмотренная в настоящей диссертации, поставлена в связи с экспериментом, опубликованном в работе С.И.Анисимова и др. (1980). На возможность появления в этом эксперименте кольцевой кумулятивной струи было впервые указано, по-видимому, В.Я.Терновым (1984). А.В.Бушманом и др. (198Т) впервые была численно решена двумерная задача, поставленная в связи с этим экспериментом, в которой дейтерий был заменен вакуумом. Автору принадлежит приоритет в решении полной задачи, когда мишень заполнена дейтерием. Численно решена задача в рамках уравнений Эйлера с реальными уравнениями состояния металлов. Численно иследован вызванный высокоскоростной струей вязкий нагрев дей-териевой плазмы в погранслое у границы с алшинием. В рамках уравнений двухтемпературной магнитной гидродинамики численно исследовано влияние на течение спонтанного электромагнитного поля.

Практическая значимость. Создан пакет программ, позволяющий численно моделировать широкий класс невязких и вязких тече-

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

Апробация. Результаты диссертации докладывались на VII Всесоюзном семинаре "Теоретические основы и конструирование численных алгоритмов решения задач математической физики" (Кемерово, 1988), на II Всесоюзном совещании по методам построения сеток для задач математической физики (Кемерово, 1988), на III Российско-японском симпозиуме по вычислительной гидродинамике (Владивосток, 1992), на семинаре под руководством В.Е.Фортова в ИВТАНе и на семинаре по методам решения задач математической физики в ВЦ РАН.

Структура и объем диссертации. Диссертация состоит из введения, трех глав и заключения. Объем диссертации 211 стр., включая 62 рис., 1 табл. и библиографию из 170 наименований.

СОДЕРЖАНИЕ РАБОТЫ В кратком введении обосновывается актуальность темы и да-

- 17 -

ется информация о структуре диссертации.

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

В §1 представлен новый метод построения регулярных сеток в областях с криволинейной границей. Рассматривается следующая задача. На плоскости (г,у) нужно построить сетку С=Цх,у) , (=1,2...гс, /=1,2,...т) при заданных координатах граничных узлов (х,у)(т, (х,у)1>?, .(х,у) . Приводится пример, для которого некоторые известные алгоритмы генерируют сетки с самопэ-ресекающимися ячейками. Вводится вспомогательная плоскость (Е.т]) и на ней - квадратная сетка (£,т]).., (=1,2,...,гг, /=1,2,...,т, такая, что т|,, =3. Для нахождения отобраше-

ния х=х(£,т]), у=у(£,,т)), определящего искомую сетку, ранее использовались нелинейные эллиптические уравнения, которые получаются после' обращения переменных в уравнениях Д£=0, Ат)=0. Предлагается аппроксимировать не эти уравнения, а функционал

(1) I = = / е 71^ 6 ад, ^х^ -х^.

Рассмотрим ячейку сетки ((+1/2^+1/2) с вершинами ((,/), ((>/+1). (1+1 ,/+1), (1+1 „Я, которые занумеруем внутри ячейки от 1 до 4. Пусть ^ - удвоенная площадь треугольника с вершиной к, полученного при разбиении ячейки диагональю, не проходящей через эту вершину, функционал (1 ) заменяется функцией

(2) Г 1

{=1 ¿ = 1 I- 4 +1/2.¿+1/2

где каждое из четырех слагаемых аппроксимирует подинтегральную

функцию на одном из треугольников, например, <р123=[ (г,-г2)г+ + (У, -Уг )2+ )2+ (2/3"Уе)2]Л7н ■ М30®90™ сеток Б, сос;--ящих из выпуклых четырехугольников, определяется набором неравенств

1/г.^+1/2>а' и1>2.....л-1; /=1,2,...т-1:

Через дБ обозначим границу множества В, состоящее из сеток, для которых хотя бы одно из этих неравенств обращается в равенство. Функция I71 обладает следующим свойством. Если Се!>, т.е.

если хотя бы одна из величин Jb для какой-то ячейки стремится к нулю, оставаясь при этом положительной, то Это свойс-

тво позволяет реализовать следующий способ построения сетки. Вначале с помощью отдельного алгоритма находится некоторая сетка а затем используется какой-либо метод безусловной минимизации функции (2), каждый шаг которого выбирается так, что СеО. Далее в §1 описывается конкретный алгоритм и приводя;гя примеры расчетов, демонстрирующие надежность этого алгоритм?, по сравнению с известными. В последнем пункте приводится обобщенно функции (2), позволяющее эффективно управлять расположением узлов.

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

изводных не используется. На первом этапе применяется какая-либо схема в переменных Лагранка, а на втором этапе поле течения с использованием законов сохранения пересчитывается со сдвинутой лагранжевой сетки обратно на заданную. Оба подхода приводят к внешне похожим разностным схемам. В качестве первого шага второй подход формулируется в виде схемы расщепления для интегродафференциальных уравнений, полученных интегрированием исходных уравнений по произвольному объему. Затем делается переход к схеме расщепления для уравнений в частных производных. После линеаризации схема расщепления первого подхода оказывается несимметричной, а второго подхода - симметричной. Исследованы свойства двух конкретных разностных схем - Годунова и Хвр-лоу, как с симметричным, так и с несимметричным расщеплением. Показано, что схемы с несимметричным расщеплением менее устойчивы и менее монотонны.

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

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

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

(3) р<11^^+я?р=0> р(1е/(11+с!1ури=0, й/<Н=дШ+(ия), е=е+|и|г/2. Здесь приведены уравнения движения и уравнение энергии, р -плотность, и- вектор скорости, 1 - время, р - давление, е -полная энергия, е - внутренняя энергия. Умножая уравнение движения скалярно на и и вычитая из уравнения энергии, получим уравнение энергии в недивергентной форме

(4) рйЕ/<^+рсИта=0.

Пусть имеется произвольная, возможно, неструктурированная сетка в трехмерном пространстве, разбивающая область течения на ячейки, некоторым образом пронумерованные. Рассмотрим множество консервативных схем, аппроксимирующих (3) на лагранжевой сетке,

*4(и;-и{) + т §Ма = 0, + 1§Р{11с1в) = О,

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

уравнение скалярно на вектор й4=(и^+и{)/2 и вычитая из второго, получим разностное уравнение, которое приводится к виду

(5) ¥{(е|-е{) + хр^Щйа) = 0, р = #Р(срйа), )/§Щйа). Нетрудно убедиться, что <р является осреднящей вектор-функцией, т. е. £(фйа)=1. Разностное уравнение (5) аппроксимирует (4) для любой осредняющей функции <р. Новые схемы, которые были названы почти консервативными, вводятся дополнительным условием

(6) (фп)^О,

где п - вектор внешней нормали к поверхности ячейки. Из (б) очевидным образом следует, что р^О при РЮ. Если функция <р, вычисленная из консервативной схемы, удовлетворяет (6) на всей поверхности, то она и используется в (5) для определения внутренней энергии, т.е. для такой ячейки счет ведется по консервативной схеме. Если условие (6) где-либо нарушено, то функция <р модифицируется следупцим образом Г 0, (Срт)^],

(7) Ф = ]

I ф/11— ; (фйа)], (<ргс)>0 .

СсрпХО

Нетрудно убедиться, что <р удовлетворяет неравенству (6); если ф удовлетворяет (б), то ф=ф; #(фйа)=1. Не видно■ никаких причин для того, чтобы преобразование (7) могло ухудшить порядок точности схемы, ее монотонность или устойчивость. Для приведенного во введении к §4 примера такая схема, так же, как и консервативная, обеспечивала сходимость к точному решению. В случае постоянства функций Р и II на каждой грани ячейки, р{ в (5) вычисляется по интерполяционной формуле р{=ЗзьРь, Х)зй=1, где суммирование ведется по всем граням ячейки. Для консервативной

схемы X (срс2а), где ф определено в (5), а интегрирование выполняется по грани. Условие (6) переходит в условие о^^О, а (7) переходит в -соответствущее преобразование коэффициентов ак.

Далее в п. 4.1 рассмотрены схемы для этапа пересчета, где предлагается ограничить консервативные формулы пересчета принципом максимума для внутренней энергии, который имеет место в случае неконсервативных формул. В п. 4.2 приведен пример расчета, демонстрирующий появление чисто счетного эффекта ложной кумуляции при использовании консервативной схемы Годунова и отсутствие этого эффекта для почти консервативной схемы. В п. 4.3 рассмотрено широкое множество консервативных схем без расщепления на произвольной подвижной сетке, аппроксимирующих дивергентную форму уравнений в переменных Эйлера, для которых аналогичным образом удается построить почти консервативные схемы.

В §5 представлена комбинированная методика расчета течений с областью из нескольких сред, допускающая' использование как явного выделения границ раздела сред в виде линий сетки, так и сквозного расчета методом концентраций. Описан вариант этого метода и приведен пример его использования, демонстрирующий остутствие грубых ошибок в определении положения границы раздела. Гибкость методики в целом демонстрируется на примере начальной стадии расчета сжатия дейтерия в конической свинцовой мишени, закрытой сверху алюминиевой крышкой. Использование метода концентраций для верхней границы между алюминием и свинцом, не оказывающей существенного влияния на динамику сжатия дейтерия, делает возможным явное выделение всех остальных гра-

ниц раздела.

Глава 2 посвящена проблеме построения монотонной схемы второго порядка точности для уравнений Навье-Стокса на сильно неравномерных сетках, позолящей одновременно вести расчет дозвукового течения в узком погранслое и интенсивных ударных волн в основной части течения. Введением к этой главе является §6, в котором, в частности, приведен пример расчета эволюционной задачи конвекции из работы Г.М.Махвиладзе и С.Б.Щербака, демонстрирующий значительную потерю точности при использовании неявной двухслойной схемы переменных направлений с шагом по времени т>т0=т1л(т° ), где 1°., определяется условием Куранта для ячейки ((,./). Для создания- экономичной схемы с нужными свойствами предлагается использовать расщепеление по "физическим процессам", строя на каждом этапе монотонную схему не ниже второго порядка точности.

В §7 рассмотрена новая неявная схема для уравнений невязкого газа в переменных Лагранжа на неортогональной сетке. Схема состоит из двух шагов. На первом шаге строится неявная схема на сетке с разнесенными узлами: термодинамические переменные определены в серединах ячеек, а в серединах граней ячеек определены нормальные компоненты скорости, которые на нижнем временном слое интерполируются по значениям вектора скорости в серединах ячеек. Для такой схемы алгебраические преобразования позоляют исключить скорость и получить уравнения относительно только давления на 9-ти точечном шаблоне. На втором шаге найденное поле давления используется для аппроксимации уравнений движения в серединах ячеек. Приводятся результаты тестовых расчетов за-

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

В §8 рассмотрен новый принцип объединения двух схем в одну составную на неравномерных сетках. Вначале рассматриваются уравнения невязкого газа в переменных Лагранжа. Пусть задан некоторый шаг по времени т. Переход на верхний временной слой состоит из двух этапов. На первом используется какая-либо явная монотонная схаМа второго порядка точности (был выбран вариант известной схемы А.В.Родионова) с шагом по времени т '

V

=min(%,%°j), своим для каждой ячейки. На втором этапе используется неявная схема-из §7 с шагом по времени z1.2.5 = т-г'1.до-

t j tj

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

Далее в §8 рассматриваются параболические уравнения, возникающие на этапе учета вязкости в уравнениях движения. После аппроксимации пространственных производных получается система обыкновенных дифференциальных уравнений на 9-ти точечном шаблоне

г 1 1

«Ъ -2 I 1

t+s.J+ff

ß=1 а=-1 s=-1

а=1,2; i=1,2J=1,2,.

f • • • 9

9 ■ • ■ 9

м

- 25 - •

где u^j - г-компонента вектора скорости в ячейке IJ, и^ -z-компонента, - коэффициенты схемы, зависящие от коорди-

нат узлов сетки. Значения и™., иа . ., и™ , uf .,

О J n+1.J tO t,m+1

i=0,1____,п+1, J=0,1,...,m+1, a=1,2 предполагаются заданными.

Рассмотрим две схемы для (8): какой-либо аналог схемы переменных направлений и полностью неявную схему, когда все слагаемые в правой части (8) взяты с верхнего временного слоя. Первая схема, реализуемая с помощью прогонок для двухкомпонентных векторов, экономична, но ответить на вопрос о ее монотонности и даже устойчивости представляется весьма затруднительным. Вторая схема абсолютно устойчива и, по-видимому, монотонна, но не экономична. Промежуточными свойствами обладает следующая составная схема. На первом этапе используется схема переменных направлений с шагом го времени т'1.где

т° = е/ max s>0.

J a.p.e.s

На втором этапе используется неявная схема с шагом i\z,). При е=>0 схема стремится к неявной, при е=*» - к схеме переменных направлений. Если положить £=0.1+0.5, то оператор перехода первого этапа схемы будет не сильно отличаться от единичного, что дает обоснованную надежду на устойчивость этого этапа. Так как правая часть (8) аппроксимирует некоторый эллиптический дифференциальный оператор второго порядка, коэффициенты bfP3S будут малы для крупных ячеек и велики для мелких. Поэтому значение е можно выбрать так, чтобы для многих крупных ячеек т^}=0 и их можно было исключить из расчета по неявной схеме. В то же время для мелких ячеек, где составная схема будет близка к

неявной. Далее в §8 рассмотрен этап учета двухтемпературной теплопроводности, где также строится экономичная составная схема.

Параграф 9 посвящен новому алгоритму пересчета высокого порядка точности между неблизкими сетками. Одномерный аналог алгоритма как разностная схема для уравнения переноса рассматривается в п. 9.1. Используется известная монотонная аппроксимация кубическими полиномами Эрмита, модифицированная таким образом, чтобы избежать потери точности в окрестности локальных экстремумов. Если скорость в уравнении переноса постоянна, то алгоритм является явной безусловно устойчивой монотонной схемой третьего порядка точности. Двумерный алгоритм рассмотрен в п. 9.2. В обоих пунктах представлены результаты расчета тестовых задач для уравнения переноса в сравнении с известной неявной TVD схемой, описанной в приложении А. Поскольку алгоритм неконсервативен, полученное с его использованием численное решение задач с ударными волнами сходится к функциям, отличным от точного решения. В конце параграфа приводится пример, дапций представление о величине такой методической ошибки. Для ее устранения предлагается строить составные- алгоритмы, которые, аналогично составным схемам из §8, объединяют одношаговые консервативные алгоритмы для близких сеток и алгоритм из настоящего параграфа. Применение предлагаемого алгоритма к сеткам с сильным сгущением вблизи подвижной криволинейной границы требует специального корректирующего алгоритма для граничных узлов. Этот алгоритм рассмотрен в Приложении Б.

Глава 3 посвящена численному исследованию сжатия дейтерия

в конических мишенях при наличии сильного кумулятивного эффекта. Глава начинается с введения, где дан краткий обзор экспериментальных и теоретических работ. . Описаны результаты работы А.В.Бушмана и др., которые рассматривали близкую задачу с заменой дейтерия на вакуум.

В §10 представлены результаты численного исследования задачи в рамках уравнений Эйлера. Используются широкодиапазонные уравнения состояния металлов (А.В.Бушман, В.Е.Фортов). Первые результаты относятся к задаче без дейтерия. Рассматриваются две свинцовые мишени с углами раствора конуса 60° и 30°. Радиус отверстия на поверхности. мишени 1 мм. Толщина ударника 2 мм, начальная скорость 5.6 км/с. Показано, что для мишени с углом раствора конуса 30° кольцевая кумулятивная струя не возникает, а с углом 60° - возникает и схлопывается на оси симметрии. Далее для мишени с углом раствора конуса 60° приводятся результаты расчета полной задачи, когда мишень заполнена дейтерием, который полагался совершенным полностью ионизованным газом. Схлопыванщийся алюминий делит дейтериевый объем на две части, генерируя в каждой из них высокоскоростную струю («180+200 км/с). Далее обе части дейтерия продолжают сжиматься вторичными кумулятивными струями алюминия. Верхняя часть при.этом вдавливается в окружающий алюминий, а нижняя часть вдавливается в свинец. С ростом давления в дейтерии его дальнейшее сжатие прекращается, а струи дейтерия и вызванные ими ударные волны затухают. Температура дейтерия оказывается весьма умеренной: 30+40 эВ. Максимальные значения температуры, вызванные появлением струй, также невелики: 60+70 эВ. В результате использование

известной формулы для расчета числа образующихся нейтронов D-D реакции дает почти нулевой результат, в то время как в эксперименте стабильно регистрировалось около 10б нейтронов. Приводится грубая нижняя оценка температуры плазмы в'эксперименте, которая дает «200 эВ для равномерно нагретой плазмы. Имеются два возможных объяснения такого несоответствия между экспериментом и расчетом. Либо в эксперименте интенсивная D-D реакция идет в отсутствие высокотемпературной плазмы, либо имеется эффективный механизм нагрева плазмы, которого нет в уравнениях Эйлера и который, следовательно, связан не только со сжатием плазмы.

В §11 представлены результаты двух попыток автора путем численного моделирования найти тот механизм нагрева плазмы, который мог бы объяснить результаты эксперимента. Обе эти попытки связаны с высокоскоростными струями дейтерия. Цель первой задачи - проверить роль вязкого нагрева, возникающего на границе такой струи., Течение в дейтерии описывается уравнениями Навье-Стокса для двухтемпературной полностью ионизованной неза-магниченной плазмы с дополнительными условиями на границах с металлами (непрерывность касательной составляющей скорости и отсутствие потока тепла). Вязкий нагрев в погранслое приводит к заметному увеличению ионной температуры (*100 эВ). Тем не менее такого нагрева явно недостаточно, чтобы объяснить результаты эксперимента.

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

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

В заключении сформулированы основные результаты диссертации. Приведем эти формулировки в сокращенном виде.

В области вычислительной математики.

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

2. Выполнен сравнительный анализ двух широко применяющихся типов разностных схем расщепления по физическим процессам.

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

4. Построены многошаговые алгоритмы пересчета с одной сетки на другую, не ограниченные условием близости сеток.

5. Построена новая методика расчета задач с областью течения из нескольких сред с различными уравнениями состояния.

6. Построена экономичная неявная схема без расщепления по направлениям, для уравнений газовой динамики в переменных Лаг-ранжа на черырехугольных неортогональных сетках.

7. Предложен новый прицип объединения двух схем в одну составную для эволюционных уравнений на неравномерных сетках.

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

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

В области механики жидкости, газа и плазмы.

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

2. С использованием известной формулы для скорости D-D реакции показано, что решение уравнений Эйлера не может объяснить появление значительного числа нейтронов D-D реакции в эксперименте .

3. В рамках уравнений Навье-Стокса двухтемпературной гидродинамики численно исследован вызванный высокоскоростной струей вязкий нагрев двйтериевой плазмы в погранслое у границы с алюминием.

4. В рамках уравнений двухтемпературной магнитной гидродинамики численно исследовано влияние на течение спонтанного элект-

ромагнитного поля.

Эти результаты опубликованы в статьях:

1. Модификация схемы С.К.Годунова в переменных Эйлера // ЖВМ и МФ, 23, JS5 (1983).

2: Об одной схеме расщепления для уравнений газовой динамики // ЖВМ И МФ, 25; Д12 (1985).

3. Об одном подходе к расчету газодинамических течений с поверхностями раздела сред // Проблемы прикладной математики и информатики. М.: Наука, 1987.

4-. Алгоритм построения криволинейных сеток из выпуклых' четырехугольников// ДАН, 295. Ж? (1987) (совместно с С.А.Иваненко ).

5. Эффект ложной кумуляции для схемы O.K. Годунова на подвиж-

(

ных сетках// ЖВМ и МФ, 28, JS1 (1988).

6. Неконсервативная разностная схема для уравнений газовой динамики на базе схемы С.К.Годунова// ВАНТ. Сер. Методики и программы числ. решения задач матем. физ. 1988. Вып. 2.

7. Криволинейные сетки из выпуклых четырехугольников// ЖВМ и МФ, 28, JS4 (1988) (совместно с О.А.Иваненко).

8. Численное моделирование сжатия газа в конических твердотельных мишенях. М.: ВЦ АН ССОР, 1988.

9. Application of moving regular grids to computation of gasdynamic flows with interfaces// Modern Problems in Comput. Aerogasdynamics, CRC Press, Boca Ration, 1992.

10. Расчет вязкого нагрева плазмы, сжимаемой ударником в конической мишени// ЖВМ и МФ, 30, J65 (1990).

11. О симметричной и .несимметричной схемах расщепления для

- 32 -

уравнений газовой динамики// ЖВМ и МФ, 31, (1991). 12. Расчет сжатия дейтерия в конической мишени в рамках уравнений Навье-Стокса для двухтемпературной магнитной гидродинамики// ЖВМ и МФ,- 33, *5 (1993).

A.A.Чарахчьян Некоторые вопросы численного моделирования нестационарных газодинамических течений

Подписано в печать 29.10.93 Формат бумаги 60x84 1/16 Тираж 100 вкз. Заказ 66. Бесплатно

Отпечатано на ротапринтах в ВЦ РАН 117967 Москва, ул. Вавилова, 40