Введение к работе
Актуальность проблемы. Методам расчета эволюционных газодинамических течений посвящено очень большое число работ. Имеется много монографий, как отечественных, так и зарубежных. Тем не менее поток публикаций по этой теме не уменьшается, поскольку как сложность задач, так и требования к точности расчета непрерывно возрастают.
Диссертация посвящена двумерным эволюционным течениям. Хотя некоторые отечественные и многие зарубежные публикации по вычислительной гидродинамике посвящены трехмерным течениям, разработка новых методов расчета двумерных течений с областью из нескольких сред с различными уравнениями состояния остается весьма актуальной. Сложность таких задач состоит в том, что с течением времени сильно меняются такие геометрические характеристики течения, как форма границ раздела или размеры области, занимаемой той или иной средой. Могут появиться новые границы раздела. Односвязная область, занимаемая какой-либо средой, может превратиться в двусвязную. Такой задачей является рассмотренное в диссертации сжатие дейтерия в конической мишени при наличии сильного кумулятивного эффекта.
В вычислительной гидродинамике имеется три направления, в рамках которых развивались методы численного решения задач с областью течения из нескольких сред. К первому направлению относятся методы сквозного счета через границу раздела сред, начало которым положил известный метод частиц в ячейке Харлоу. Здесь используется либо техника смешанных ячеек, состоящих из нескольких сред (метод частиц в ячейке; метод концентраций,
- г -
Баграх СМ. и др., 1981), либо специальная процедура переинтер-поляции данных с лагранжевых частиц на эйлерову сетку и обратно (метод индивидуальных частиц, Агурейкин В.А., Криков Б.П., 1986). Другое направление образуют активно развивающиеся в последнее время свободно-лагранжевы методы, когда шаблон разностной схемы определяется на каждом шаге по времени в зависимости от расположения узлов лагранжевой сетки. Несомненным достоинством методов из этих двух направлений является их универсальность, т.е. возможность сравнительно быстрого переключения соответствующих пакетов программ с одной задачи на другую. Вместе с тем естественным продолжением этого достоинства является следующий недостаток. В методах первого направления форма границы раздела сред вообще не определяется, а в методах второго направления граница раздела сред как правило представляется ломаной линией. В то же время границы раздела сред в решении уравнений Эйлера являются гладкими кривыми. В результате использование методов из первых двух направлений для расчета течений с интенсивным перемещением одной из сред вдоль границы раздела приводит к "схемному1? трению, т.е. к торможению одной среды и возмущению в другой, которых нет в решении уравнений Эйлера. Поэтому для получения надежных результатов эти методы требуют очень подробных сеток в окрестности границ раздела.
Настоящая диссертация посвящена развитию третьего из указанных выше направлений, основы которого заложены известной монографией О.К.Годунова и др. (1976). Граница раздела сред в каждый момент времени представляется гладкой кривой, вдоль которой расставляются узлы некоторой линии сетки. "Схемного" тре-
- З -ния вдоль границы раздела здесь нет, однако до работ автора такой подход удавалось использовать .только для относительно узкого круга задач. Сюда не попадали задачи с несколькими пересекающимися границами раздела, когда выделить явно все эти границы не представляется возможным. Кроме того границы раздела могут с течением времени сильно искривляться, а удовлетворительного алгоритма расчета сеток с заданной структурой в областях с сильно искривленными границами фактически не было. При решении задачи о сжатии дейтерия автором было опробовано несколько известных алгоритмов генерации сеток. В некоторый момент времени все эти алгоритмы начинали генерировать непригодные для проведения расчетов сетки, часть ячеек которых были самопересекающимися. Поэтому приходилось периодически изменять структуру сетки и проводить переинтерполяцию решения на новую сетку, что значительно снижало точность расчета, не говоря уже о заметном увеличении сложности всего алгоритма. Все эти трудности привели к тому, что до работ автора задачи описанного выше класса решались в основном по методикам, построенным либо на методах сквозного счета через границу раздела сред, либо на свободно-лагранжевых методах.
В настоящее время для расчета эволюционных течений используются как консервативные, так и неконсервативные схемы. За исключением полностью консервативных схем, о которых особо будет сказано ниже, расчет внутренней гнергии є в консервативных схемах происходит следующим образом. Из разностного уравнения энергии определяется разностная функция полной энергии е=Е+|и|2/2, где вектор скорости и определяется из разностного
- 4 -уравнения движения, а значение є определяется соответствующим вычитанием. Неконсервативные схемы строятся на базе недивергентного уравнения энергии, из которого исключена кинетическая энергия. В этом случае значение є определяется непосредственно из разностного уравнения энергии. Поскольку разрывные решения уравнений газовой динамики определяются законами сохранения, естественно предположить, что разностное решение, полученное по неконсервативной схеме с невыполнящимися точно законами сохранения, при наличии разрывов будет сходиться к функциям, отличным от точного решения. В диссертации приведен такой пример, связанный с неконсервативным вариантом схемы Годунова, построенном на недивергентной форме уравнения энергии. Аналогичные примеры есть и в работах других авторов. Тем не менее анализ литературы, как отечественной, так и зарубекной, показывает, что, несмотря на наличие указанной выше методической ошибки, неконсервативные схемы отнюдь не являются менее распространенными по сравнению с консервативными, по крайней мере при расчете эволюционных течений. Причина заключается в потенциальной возможности потери точности расчета внутренней энергии в консервативных схемах. Такой недостаток отсутствует в полностью консервативных схемах (Самарский А.А., Попов Ю.П.; Гольдин В.Я., Ионкин Н.И., Калиткин Н.Н.), полученных путем выделения консервативных схем из множества схем, вшпхжсимирущих некоторым естественным образом недивергентное уравнение энергии. Эти схемы построены на сетках с разнесенными узлами, когда термодинамические переменные и вектор скорости определены в разных узлах сетки. В то же время имеетоя много схем, в частности схе-
- 5 -ма Годунова и все ее модификации высокого порядка точности, построенных на сетках с совмещенными узлами, когда все переменные определены в одних и тех же узлах. Поэтому актуальной проблемой вычислительной-гидродинамики является построение схем на сетках с совмещенными узлами, которые, с одной-стороны, были бы близки к консервативным схемам настолько, чтобы избежать методической ошибки при расчете разрывных решений, и, с другой стороны, гарантировали бы отсутствие потери точности расчета внутренней анергии.
По сравнении с невязкими течениями, расчет вязких течений представляет из себя значительно более сложнуи проблему вычислительной гидродинамики. Эта сложность связана не столько с появлением дополнительных слагаемых в уравнениях, сколько с появлением нового масштаба - толщины погранслоя, где вязкое течение существенно отличается от невязкого. В упоминавшейся выше задаче о сжатии дейтерия в конических мишенях схлбпываю-щийся на оси симметрии алюминий гененировал высокоскоростные струи дейтерия. Необходимо было проверить гипотезу о значительном увеличении температуры плазмы вдоль границы струи за счет вязкого нагрева. Линейные размеры области дейтерия при сжатии уменьшались в 5+10 раз, а толщина погранслоя в дейтерии была еще в 50+100 раз меньше. Наряду с дозвуковым течением в погран-слое, в основной части течения возникали сильные ударные волны, порожденные высокоскоростными струями. Такая задача, даже оставаясь двумерной, по своей сложности значительно превосходит задачи, когда-либо рассмотренные как в отечественной так и в зарубежной литературе по вычислительной гидродинамике.
- б -
Значительный прогресе в области численного моделирования вязких течений был достигнут в связи с необходимостью расчета вязкого нагрева при входе тел в атмосферу планет. Из отечественных работ отметим работы О.М.Белоцерковского, Ю.П.Головачева, В.М.Ковени, А.И.Толстых, а также работы М.Я.Иванова, В.Г.Крупы и Р.З.Нигматуллина по расчету внутренних вязких течений. Здесь, по-видимому, впервые в реальных многомерных задачах стали использоватся сетки с сильным сгущением, необходимым для расчета течения в погранслое. Тем не менее применение развитых здесь методов для расчета упомянутой выше задачи о сжатии дейтерия в конической мишени связано с серьезными трудностями. И дело не столько в том, что в етой задаче граница области вязкого течения подвижна, а используемые сетки неортогональны. Основная трудность заключается в следующем.
Обозначим через а шаг по времени, который следует из условия Куранта для ячейки сетки (IJ) и определяется размером ячейки и местной скоростью звука. При наличии погранслоев шаг по времени T0=mln(T ) оказывается очень маленьким из-за небольших размеров ячеек в погранслое. К - счастью течение в погранслое существенно дозвуковое, и поэтому скорость изменения искомых функций там много меньше скорости звука, что делает возможным расчет эволюционного течения с шагом по времени т»т. В то же время задача о входе тел в атмосферу планет, по крайней мере в связи с расчетом вязкого нагрева, является стационарной задачей. Поэтому развитые здесь методы являются неявными схемами переменных направлений,-которые очень эффективны для расчета стационарных течений методом установления. Однако использование
- 7 -таких схем для расчета эволюционных течений с шагом по времени т»т может приводить к существенной потере точности. Это было наглядно продемонстрировано в работе Г.М.Махвиладзе и С.В.Щербака (препринт *113, ИПМехан АН СССР, 1978) на примере задачи конвекции. В этой задаче течение везде дозвуковое, что позволило авторам получить надежные результаты на равномерных сетках. Было показано, что точность расчета эволюционной задачи по некоторой двухслойной схеме переменных направлений заметно падает уже при тм4т при отстутствии каких-либо признаков нарушения устойчивости расчета. Авторами была предложена трехслойная схема переменных направлений, для которой эффект понижения точности отсутствовал вплоть до т«8т. Тем не менее при т«1бт точность расчета опять заметно падала. Заметим, что для прикладных задач, решаемых на пределе возможностей ЭВМ, диагностика такой потери точности является серьезной проблемой.
Причина потери точности заключается в следующем. Дозвуковое течение является одновременно и слабосжймаемым, т.е. в уравнении неразрывности p_1 (3/flt+tr7)p=-dlvu левая часть много меньше чем" каждая из пространственных производных, входящих в divu. В результате при т»т погрешность аппроксимации в основном определяется той своей частью, которая вносится в схему процедурой расщепления по направлениям. Поэтому схемы для дозвуковых течений следует строить аналогично схемам для несжимаемой жидкости: оба слагаемых, входящих в divu следует аппроксимировать на одном и том же временном слое. Хотя на каждом временном слое приходится решать двумерную систему алгебраических уравнений, допустимый шаг по времени практически перестает за-
-8-. висеть от условия Куранта. Это было продемонстрировано в работе В.А.Гончарова, В.М.Кривцова и автора (1988), где построенная таким образом схема позволила проводить расчет упомянутой выше задачи конвекции с шагом по времени x«64t без заметного падения точности.
Таким образом, методы, развитые в связи с расчетом вязкого нагрева при входе тел в атмосферу планет, мокно использовать для расчета эволюционных течений с шагом по времени т»т;0 только отказавшись формальна от расщепления по направлениям. Однако такая схема оказывается очень неэкономичной из-за большой размерности системы алгебраических уравнений на каждом временном слое, в которую входят значения всех неизвестных функций - давления, плотности и компонент вектора скорости.
Весьма актуальной является рассмотренная в диссертации задача численного исследования сжатия дейтерия ударниками в конической мишени при наличии сильного кумулятивного эффекта. Результат соответствующего натурного эксперимента до сих пор не получил удовлетворительного объяснения. Неокиданность этого результата в сравнении с другими аналогичными экспериментами заключается в заметном («10б) нейтронном выходе при сравнительно небольшой (*5.4 км/с) скорости ударника. Хотя этот эксперимент не был продолжен после 1980 г. и, по-видимому, не был никем повторен, авторы эксперимента по-преннему уверены в надежности своего результата, неизменно включая его в свои обзорные статьи, последняя из которых опубликована в 1992 г.
Целью работы является: - Значительно поднять уровень сложности задач с несколькими
- 9 -средами, доступных для численного моделирования в рамках методики с представлением границ раздела сред гладкими кривыми..
Построить методику расчета эволюционных вязких снимаемых течений в областях с подвижной границей, имея в виду возможность одновременного расчета дозвукового течения в узком по-гранслоэ и сильных ударных волн в основной части течения.
Численно исследовать сжатие дейтерия в конических мишенях при наличии сильного кумулятивного эффекта с учетом возможного влияния на течение вязкого нагрева плазмы и спонтанного электромагнитного поля.
Научная новизна. Остановимся вначале на результатах, полученных в связи с разработкой новой методики расчета задач с несколькими средами. Методика развивает известный подход к таким задачам, когда граница раздела сред представляется гладкой кривой. Выше уже отмечались трудности, связанные с применением такого подхода к задачам со сложной геометрией. Первая из них заключается в том, что в случае нескольких пересекающихся границ раздела их явное выделение, по крайней мере в случае сеток фиксированной структуры, может оказаться невозможным. Эта трудность была преодолена путем компиляции в единую методику известных алгоритмов. Использование расщепления на лагранжев этап и этап пересчета поля течения с лагранжевой сетки на подвикную эйлерову позволило в рамках одной разностной схемы проводить расчет как с явном выделением границы раздела, так и "сквозным" методом концентраций, в котором граница раздела сред моделируется полосой смешанных ячеек. Это существенно увеличило гибкость методики. Для задач с несколькими границами раздела сред
- 10 -появилась возможность оставаться в рамках сеток с фиксированной структурой, "заплатив" зв это сквозным расчетом некоторых границ раздела, не оказывающих значительного влияния на ту область течения, которая представляет основной интерес.
Другая трудность, о которой уже упоминалось выше, - возможное сильное искривление границ раздела. Для построения сеток в областях с криволинейными границами широко используются нелинейные эллиптические уравнения, которые определяют отображение области на прямоугольник во вспомогательном пространстве. Автор совместно с О.А.Иваненко предложил новый метод построения сеток, в котором аппроксимируются не сами уравнения, а тот функционал, минимум которого достигается на решении этих уравнений. Удалось найти такую аппроксимацию функционала, минимизация которой гарантирует выпуклость всех четурехугольных ячеек сетки практически при любом искривлении границ области, что было подтверждено многочисленными расчетами.
Следующий результат заключается в построении нового класса схем для уравнений газовой динамики. Стимулом для этого исследования послужили численные эксперименты на различных сетках, которые привели к обнаружению так называемого эффекта ложной кумуляции - возникающего при расчете на подвижных сетках чисто вычислительного эффекта существенного завышения интенсивности кумулятивной струи, качественным образом менявшего характер течения. Первоначально эффект был обнаружен при сравнении различных вариантов расчета задачи о сжатии в конических мишенях, а затем наглядно воспроизведен на достаточно простой задаче. Как оказвлось, причина этого эффекта заключается в консерватив-
ности разностной схемы (использовалась схема Годунова). При сохранении полной энергии изменение внутренней энергии на лаг-ранжевом этапе в некоторых ячейках сетки оказывается противоположным по знаку тому изменению, которое следует из недивергентной формы уравнения энергии. Отталкиваясь от этого наблюдения, было найдено такое преобразование консервативной схемы в переменных Лагранкв, которое, с одной стороны, давало схему, близкую к консервативной, а, с другой стороны, гарантировало отсутствие потери точности расчета внутренней энергии. Это было подтверждено тестовыми расчетами. В отличие от неконсервативной схемы, новая схема, так же как и консервативная, сходилась к некоторому разрывному решению уравнений газовой динамики. В то же время эффект ложной кумуляции, в отличие от консервативной схемы, для новой схемы отсутствовал. Найденное преобразование консервативной схемы применимо к широкому классу схем, как двумерных, так и трехмерных, построенных на сетках произвольной структуры, в частности к схеме Годунова и ее модификациям высокого порядка точности - схемам В.П.Колгана, Ван Леера, А.В.Родионова, схеме с параболической интерполяцией (Colella and Woodward, 1984), по крайней мере к двуслойным МО (essentially non-oaclllatory) схемам, и др. Не видно никаких причин, по которым это преобразование может ухудшить такие свойства исходной консервативной схемы, как порядок аппроксимации, устойчивость и монотонность. Рассмотрены также схемы без расщепления на лаг-ранжев этап и этап пересчета. Найдено аналогичное преобразование широкого класса таких консервативных схем на подвижных сетках произвольной структуры.
- 12 -В состав методики входит алгоритм пересчета поля течения с лагранжевой сетки на заданную подвижную. Как правило эти сетки
- близки, и поэтому используются быстрые алгоритмы, в которых вещество из ячейки старой сетки передается только в ближайшие по значениям индексов ячейки новой сетки. Применимость таких алгоритмов ограничена условием близости сеток. Вместе с тем для сложных течений часто возникает необходимость в расчетах с различными законами расстановки узлов выделенных линий сетки или с различными алгоритмами расчета внутренних узлов. Смена алгоритма происходит в некоторый момент времени. В этом случае условие близости сеток оказывается нарушенным. Поэтому были разработаны многошаговые алгоритмы, которые анализируют условие близости и, в случав необходимости, вводят промежуточные сетки.
Помимо использованного в методике расщепления на лагранжев этап и этап пересчета, достаточно распространенным является другое близкое расщепление, используемое в методах "крупных частиц" и "жидкости в ячейке". В диссертации выполнен сравнительный анализ обеих расщеплений. На примере схем Харлоу и Годунова показано, что разностные схемы с использованием первого расщепления более устойчивы и более монотонны по сравнению с
' аналогичными схемами, построенными на втором расщеплении.
Перейдем теперь к результатам, полученным в связи с расчетом эволюционых вязких сжимаемых течений. Проблема заключается в построении монотонной (т.е. без заметных численных осцилляции) схемы второго порядка точности для уравнений Навье-Стокса на подвижных неортогональных сильно неравномерных сетках. Как уже указывалось выше, развитые до работ автора неявные схемы
- 13 -переменных направлений для расчета эволюционных течений не годятся, а формальный отказ от расщепления по направлениям приводит к очень неэкономичным схемам. Развитая в диссертации схема является схемой расщепления по физичесикм процессам: на первом этапе аппроксимируются уравнения Эйлера в переменных Лагранжа, на втором этапе учитываются слагаемые, описывающие вязкость и теплопроводность, а третий этап является алгоритмом пересчета поля течения с лагранжевой сетки на заданную, йтвечающую очередному моменту времени. На каждом этапе строится монотонная схема не ниже второго порядка точности по пространственным переменным.
Первый результат заключается в построении новой неявной схемы второго порядка точности на неортогональной сетке для уравнений Эйлера в переменных Лагранжа. Отсутствие расщепления по направлениям делает схему пригодной для расчета эволюционных течений с шагом по времени т>т. В то же время схема достаточно экономична, поскольку получающаяся система алгебраических уравнений содержит в качестве неизвестных значения только одной скалярной функции.
Следующий результат является центральным с точки зрения создания схемы, позволяющей вести расчет одновременно и дозвукового течения в погранслое, и сильных ударных волн в основной части течения. Каждый из этих типов течения предъявляет свои, в значительной степени противоположные требования к разностным схемам. Для расчета дозвуковых течений нужна неявная' схема без расщепления по направлениям, позоляющая вести счет с шагом х»10,. При расчете ударных волн одно из основных требований к
- 14 -разностной схеме это отсутствие заметных численных колебаний. Решение.этой проблемы было найдено в виде новогс принципа объединения двух разных схем на лагранхевом етапе расщепления в одну составную схему. Такая схема в погранслое переходит в неявную схему, в основной части течения переходит в явную монотонную схему второго порядка точности, а плавность перехода от одной схемы к другой определяется плавностью перехода от мелкой сетки к крупной. Такой ке принцип построения схемы удается эффективно использовать и для параболических уравнений, возникающих при учете вязкости и теплопроводности. Здесь на неортогональных сетках получена схема, которая по своей экономичности мало отличается от схемы переменных направлений, и одновременно оказывается близкой к чисто неявной схеме, абсолютно устойчивой и монотонной.
На защиту выносится также новый монотонный алгоритм пересчета высокого* порядка точности между неблизкими сетками, который используется на последнем этапе расщепления. Алгоритм состоит из экономичной процедуры поиска узлов старой сетки для каждого узла новой и двумерной интерполяционной процедуры на неортогональной сетке. Используются одномерные монотонные интерполяции полиномами Эрмита (De Боог and Swartz, 1977; Pritch. and Carlson, 1980), модифицированные таким образом, чтобы не допускать падения точности в окрестности локальных экстремумов. Алгоритм можно рассматривать в качестве монотонной явной разностной схемы третьего порядка точности для некоторого уравнения переноса, которая имеет повышенную устойчивость (если скорость переноса постоянна, то схема абсолютно устойчива). Выпол-
- 15 -нен [..я . тестовых расчетов, как одномерных, так и двумерных, в частности в сравнении с известной неявной TVD схемой второго порядка точности. Во всех тестах построенный алгоритм давал более точные результаты.
Остановимся теперь на задаче о сжатии дейтерия в конических мишенях. Подобные задачи в различных постановках рассматривались в ряде работ, однако в большинстве из них кумулятивный эффект не оказывал заметного влияния на течение. Начальная стадия образования кумулятивной струи численно исследовалась в работе В.В.Рассказовой, В.Г.Рогачева, Н.Ф.Свидинской (1985). Задача, рассмотренная в настоящей диссертации, поставлена в связи с экспериментом, опубликованном в работе С.И.Анисимова и др. (1980). На возможность появления в этом эксперименте кольцевой кумулятивной струи было впервые указано, по-видимому, В.Я.Терновым (1984). А.В.Бушманом и др. (198Т) впервые была численно решена двумерная задача, поставленная в связи с этим экспериментом, в которой дейтерий был заменен вакуумом. Автору принадлежит приоритет в решении полной задачи, когда мишень заполнена дейтерием. Численно решена задача в рамках уравнений Эйлера с реальными уравнениями состояния металлов. Численно ислэдован вызванный высокоскоростной струей вязкий нагрев дей-териевой плазмы в погранслое у границы с алтинием. В рамках уравнений двухтемпературной магнитной гидродинамики численно исследовано влияние на течение спонтанного электромагнитного поля.
Практическая значимость. Создан пакет программ, позволяющий численно моделировать широкий класс невязких и вязких тече-
ний. Развитые новые алгоритмы могут использоваться в других пакетах программ. Реализация алгоритма построения.сеток в криволинейных областях сдана в Государственный фонд алгоритмов и программ. Найденное преобразование консервативной схемы в почти консервативную позволяет сравнительно легко модифицировать имеющиеся пакеты программ, базирующиеся на схеме Годунова, и ее модификациях высокого порядка точности, что. сделает невозможным потерю точности расчета внутренней анергии. Дана легко выполнимая рекомендация по модификации имеющихся пакетов программ, основанных на методах "крупных частиц" и "кидкости в ячейке", которая сделает получаемые численные решения более устойчивыми и более монотонными.
Апробация. Результаты диссертации докладывались на VII Всесоюзном семинаре "Теоретические основы и конструирование численных алгоритмов решения задач математической физики" (Кемерово, 1988), на II Всесоюзном совещании по методам построения сеток для задач математической физики (Кемерово, 1988), на III Российско-японском симпозиуме по вычислительной гидродинамике (Владивосток, 1992), на семинаре под руководством В.Е.Фортова в ИВТАНе и на семинаре по методам решения задач математической физики в ВЦ РАН.
Структура и объем диссертации. Диссертация состоит из введения, трех глав и заключения. Объем диссертации 211 стр., включая 62 рис., 1 табл. и библиографию из 170 наименований.
Похожие диссертации на Некоторые вопросы численного моделирования нестационарных газодинамических течений
