CC BY
25
УДК 533.6
В. М. Пасконов, М. Г. Лебедев
РАЗВИТИЕ МОДЕЛЕЙ И МЕТОДОВ ВЫЧИСЛИТЕЛЬНОЙ АЭРОГИДРОДИНАМИКИ1
(кафедра математической физики, лаборатория моделирования процессов тепломассопереноса
факультета ВМиК)
Одной из центральных проблем аэродинамики во второй половине XX в. являлось исследование сверхзвукового обтекания затупленных тел. В изучении этой проблемы решающую роль сыграли вычислительные методы, которые первоначально были разработаны на основе уравнений Эйлера, а впоследствии — на основе уравнений Навье-Стокса для вязкого сжимаемого теплопроводного газа. Работы в этом направлении в МГУ, и в частности на факультете ВМиК, инициировал академик Г. И. Петров в связи с необходимостью решения задач теплозащиты летательных аппаратов, расчета отрывных течений, аэродинамических следов и определения донного давления.
При сверхзвуковом полете летательного аппарата в атмосфере структура течения газа во всей возмущенной области достаточно сложна и в значительной степени зависит от формы тела. Построение численного решения сложной структуры для системы уравнений Навье-Стокса разностным методом в области значительной протяженности предъявляет специфические требования к методике расчета [1]. Это связано с тем, что одновременно должны рассчитываться такие разные по свойствам области резких изменений всех искомых функций, как области ударных волн различной интенсивности (головная и хвостовая ударные волны), область пограничного слоя вблизи поверхности обтекаемого тела, область возвратно-циркуляционного течения, области ближнего и дальнего следа [1-3]. Расчет таких областей предъявляет высокие требования к точности, которую могут обеспечить достаточно мелкие (подробные) разностные сетки. Структура этих сеток в какой-то мере должна согласовываться с полученным решением, что может быть достигнуто путем введения существенно неравномерных в физической плоскости сеток.
Подход к численному моделированию рассматриваемой задачи, позволяющий частично обойти указанные трудности, был основан на разделении всей области интегрирования на ряд взаимно перекрывающихся областей и последовательном расчете течения в каждой из них. Такое разделение области возможно в силу слабой передачи возмущений вверх по потоку при обтекании тел сверхзвуковым набегающим потоком вязкого газа.
Для численного исследования стационарных и нестационарных течений вязкого сжимаемого теплопроводного газа на основе уравнений Навье-Стокса были разработаны явно-неявные разностные сплайн-схемы повышенного порядка точности [4-6], итерационные процедуры, связанные с существенной нелинейностью систем дифференциальных уравнений, построение разностных сплайн-сеток для областей сложной формы.
Из работ последнего десятилетия отметим работы по численному моделированию течения в следе за пластиной конечной толщины, обтекаемой сверхзвуковым потоком, а также влияния включения поперечных сверхзвуковых вдувов с верхней и нижней боковых поверхностей вблизи кормового среза [7-10]. В этих работах показано влияние граничных условий симметрии на результаты расчета вихревых течений за кормой при умеренных числах Рейнольдса от 400 до 1200 и при числе Маха, равном 2. Приводятся распределения давления и теплового потока на корме пластины в зависимости от температуры стенки кормы.
На основании проведенных расчетов установлено, что применение условий симметрии для расчета течения только в верхней полуплоскости за пластиной конечной толщины приводит к уменьшению области возвратно-циркуляционного течения в ближнем следе и уменьшению интеграла донного давления. Интересно отметить также существенное влияние температурного фактора при изменении температуры поверхности пластины от Тш = 0,1 до Тш = 1,0 на структуру течения в ближнем следе при расчете двустороннего обтекания кормы пластины конечной толщины сверхзвуковым потоком для числа Маха М = 2 и в диапазоне изменения числа Рейнольдса от Ее = 400 до Ее = 1200.
1 Работа выполнена при финансовой поддержке РФФИ, проекты № 02-07-90407, № 04-01-00332.
Интересные, и в какой-то степени неожиданные, результаты получены в исследовании структур течений, при которых с боковых поверхностей производятся различные по модулю вдувы. В зависимости от интенсивности вдувов с боковых поверхностей и температуры пластины получены стационарные и нестационарные режимы течения. Потеря стационарности происходит в диапазоне чисел Маха вдуваемых струй 1,005 < М{ < 1,03 и при температуре поверхности пластины 0,1 < Tw < 0,4 (здесь Mi (г = 1,2) — числа Маха боковых вдувов).
Отметим также, что расчеты проводились на многопроцессорной вычислительной системе PowerXplorer фирмы Parsytec Computer Gmbh. Переход к использованию многопроцессорных вычислительных систем для моделирования течений вязких жидкостей и газов на основе модели Навье-Стокса открывает возможность расчета задач с использованием разностных методов на сетках, содержащих большое количество узлов, что необходимо для исследования течений с вихревыми структурами.
Однако проблема расчета вязких течений при больших числах Рейнольдса с помощью пусть даже существенного увеличения узлов разностной сетки решена быть не может. Это показало сравнение результатов расчетов и экспериментальных данных по сверхзвуковому обтеканию сферически затупленных конусов. Расчеты проводились в лаборатории моделирования процессов тепломассопереноса ф-та ВМиК, а эксперименты — на баллистической трассе Ленинградского физико-технического института им. А. Ф. Иоффе. Сравнение проводилось по характеру изменений донного давления с увеличением числа Рейнольдса. При ламинарном режиме течения в диапазоне изменения числа Рейнольдса от 103 до 104 результаты расчетов и экспериментальные данные хорошо соответствовали друг другу. При числах Рейнольдса больше 3 • 105 в расчетах характер изменения донного давления оставался прежним, а данные эксперимента показывали резкое изменение донного давления, что связывалось с переходом от ламинарного режима течения к турбулентному.
Если оставаться на позиции применимости уравнений Навье-Стокса для исследования течений при больших числах Рейнольдса, то причину сложности вычислительной реализации таких режимов следует отнести к потере точности в самой вычислительной модели. Построение вычислительной модели обычно связывается с дискретизацией дифференциальной задачи, в данном случае с построением разностных уравнений, аппроксимирующих соответствующую систему дифференциальных уравнений с необходимыми граничными условиями. Важным, но иногда забываемым моментом является взаимное согласование принципов и предположений, при которых получены дифференциальные и разностные уравнения. Для моделирования вязких течений этот момент может являться существенным. Заметим, что вывод основных уравнений динамики вязкой жидкости базируется фактически на дискретизации среды с последующим предельным переходом. В то же время расчеты ведутся на основе дискретной модели. Естественно, вычислительная модель должна учитывать предположения, заложенные в дифференциальной модели. Известно, что вывод уравнений Навье-Стокса основан на первой теореме Гельмгольца, описывающей скоростное поле сплошной среды в окрестности точки и дающей выражение компонент тензора скоростей через проекции вектора скорости на оси декартовой системы координат. Эти выражения получены с использованием линейного приближения представления скорости в виде ряда Тейлора (в терминологии разностных схем) с первым порядком точности относительно пространственных шагов сетки, если ячейку или их некоторую совокупность рассматривать как элементарный объем. Следовательно, аппроксимация уравнений, полученных при таком предположении, разностными схемами выше первого порядка теряет смысл, так как сам вывод уравнений предполагает линейное изменение скорости в элементарном объеме. Однако использование разностных схем первого порядка точности для решения уравнений Навье-Стокса приводит фактически к искажению физической вязкости, являющейся малым параметром при второй старшей производной, за счет искусственной "счетной" вязкости, вследствие чего нельзя определить, какому числу Рейнольдса соответствует такой расчет. В то же время соответствие результатов расчета числу Рейнольдса, при котором проводились вычисления, имеет существенное значение для анализа результатов, их интерпретации и сравнения с экспериментальными данными.
В вычислительной практике переход к разностным уравнениям, как правило, осуществляется на основе разностных схем, аппроксимирующих систему уравнений Навье-Стокса с порядком не ниже второго относительно пространственных шагов сетки. Но и в этом случае "искажение" физической вязкости присутствует при расчетах за счет линейного изменения скорости в каждой ячейке сетки в силу предположения, сделанного при выводе самих дифференциальных уравнений. В работах [11, 12] описан один подход к модификации классической модели Навье-Стокса, преследующий цель — увеличение точности самой модели. Для этого используется квадратичное приближение при разложении
скорости в ряд Тейлора в процессе выражения элементов тензора скоростей деформации через компоненты вектора скорости. При этом естественно возникает симметричный тензор третьего ранга — тензор ускорения деформации, который в результате умножения его на вектор смещений, определяемый в каждой точке пространства компонентами скорости и некоторым временем (например, временным шагом при численном решении), дает симметричный тензор второго ранга, который проводит коррекцию тензора скоростей деформации в основном реологическом соотношении модели Навье-Стокса. В результате получается модифицированная система уравнений Навье-Стокса, состоящая из трех уравнений количества движения третьего порядка, уравнения энергии второго порядка с модифицированной диссипативной функцией, уравнения неразрывности и уравнения состояния.
В качестве модельного уравнения при разработке численного метода решения краевых задач для эволюционных уравнений в частных производных третьего порядка было выбрано нелинейное уравнение типа Кортевега-де Фриса-Бюргерса с двумя малыми параметрами при старших производных. В работах [13, 14] предложен разностный метод численного решения краевой задачи на основе аппроксимации повышенного порядка точности для этого уравнения, построена итерационная разностная сплайн-схема и показана возможность получения численного решения уравнения третьего порядка при определенных соотношениях между двумя малыми параметрами и шагами по пространству и времени. Исследован механизм взаимодействия решения вырожденной задачи (при отсутствии третьей производной) с решением уравнения с двумя малыми параметрами. Разностная схема строилась на основе аппроксимаций повышенного порядка точности, полученных в виде линейных комбинаций конечно-разностных операторов от функции и ее сплайн-производных [4].
На каждом временном слое расчет уравнения производился в три этапа. На первом этапе для нахождения предварительных значений искомой функции проводилось решение уравнения второго порядка с заданными граничными условиями по явно-неявной разностной схеме с итерационным процессом типа Зейделя [2] или по неявной схеме с использованием метода прогонки.
Второй этап заключался в вычислении по полученным на первом этапе значениям функции методом прогонки значений второй производной во внутренних точках заданного отрезка. Дополнительно задавалось условие непрерывности второй производной в предграничных точках [4]. Затем решение уравнения второго порядка уточнялось путем внесения поправки. Для нахождения искомых значений функции строился итерационный процесс.
Третий этап состоял в уточнении полученных значений функции путем внесения поправки от производной третьего порядка и проведения необходимых итераций. Здесь использовалось граничное условие для второй производной от решения. Таким образом находилось решение разностного уравнения третьего порядка на каждом временном слое. Численное решение вырожденной задачи (при отсутствии третьей производной) состояло в реализации первого и второго этапов, описанных выше.
Приведем некоторые результаты расчетов, которые показывают влияние малых параметров при старших производных на характер решения.
-11=_,_I_,_I
0,4 0,6 х
Рис. 1. Численное решение уравнения третьего порядка при значениях параметров г] = 0,01, М = 0,0001, с!х = 0,0125; ц = 0,0001 (1); 0,0005 (2); 0,001 (5); 0,005 (4)
Для начальной функции "солитон" проводились расчеты на равномерных сетках. На рис. 1 в момент времени £ = 1,0 показано решение задачи для уравнения третьего порядка типа Кортевега-де Фриса-Бюргерса в зависимости от значения параметра ¡1 при третьей производной. Увеличение этого параметра происходит постепенно от линии 1 к линии 4- При закрепленном значении параметра
при второй производной с уменьшением коэффициента при третьей производной снижается скорость движения фронта волны и амплитуда осцилляций на фронте. Таким образом, можно предположить, что для полной модифицированной системы уравнений подобный "осциллятор" может существенно изменить поведение решения при увеличении числа Рейнольдса.
Обращение к моделированию течений вязкой несжимаемой жидкости в [15] было связано с простым переходом от системы уравнений Навье-Стокса для вязкого сжимаемого теплопроводного газа к системе уравнений при постоянной плотности. Получающаяся при этом система уравнений является переопределенной, и уравнение состояния дает линейную связь давления и температуры. Многочисленные работы по динамике вязких несжимаемых жидкостей в течение всего предыдущего времени, как правило, были основаны на классической системе уравнений Навье-Стокса, состоящей из динамических уравнений и уравнения неразрывности. В работах [15, 16] была предложена и реализована методика численного исследования задач динамики вязких несжимаемых жидкостей и газов в рамках классической модели Навье-Стокса на основе системы уравнений, сохраняющих прямую и обратную связь динамических уравнений и уравнения энергии и при численных расчетах, обеспечивающих высокую точность выполнения закона сохранения массы. Проведены численные исследования нестационарных течений в плоских каналах при различных граничных и начальных условиях и установлено существование различных неклассических структур решения. Так, например, при малых числах Рейнольдса, когда жидкость инжектируется в канал, заполненный такой же жидкостью, течение характеризуется плавным переходом от постоянного профиля скорости на входе к профилям параболического вида с постепенным нелинейным уменьшением максимального значения скорости в центре канала вплоть до полного затухания, при этом течение становится установившимся.
и
О 0,2 0,4 0,6 0,8 1
а б
Рис. 2. Нестационарное течение в плоском канале
На рис. 2,6 для канала длиной в 5 раз больше его высоты и для числа Рейнольдса Ее = 10 при £ = 10 изображена поверхность продольной составляющей вектора скорости, на которой показано несколько секущих плоскостей, а на рис. 2,а — соответствующие им профили продольной составляющей вектора скорости. Существенное изменение структуры течения в канале, связанное с возникновением отрывных вихревых течений при входе в канал, приводит к полной его блокировке при Ее = 140 даже для короткого канала, длина которого равна его высоте Ьс = 1, а также для каналов большей протяженности. На рис. 3 приведены поля направлений векторов скорости в набегающем равномерном потоке, в канале (0,5 < х < 1,5) и в следе за каналом для различных моментов времени.
Взаимодействие вихрей приводит к возникновению стоячей волны клиновидной формы, которая поддерживается подсосом жидкости в канал из бесконечного затопленного пространства, при этом сохраняется симметрия решения.
Развитие течения в канале длиной Ьс = 5 для 11е = 700, когда входящий поток не параллелен стенкам, а имеет угол атаки а = 45°, приведено на рис. 4, где показаны фрагменты полного поля направлений векторов скорости в начале канала для различных моментов времени.
О 0,25 0,5 0,75 1 1,25 1,5 1,75 2
Рис. 3. Поля направлений вектора скорости при течении в канале с числом Рейнольдса Но = 140
в моменты времени < = 20 (а); 40 (б); 80 (в)
Расчеты при числах Рейнольдса от 1400 до 14 000 показали, что при симметричном начальном профиле скорости характер изменения решения вдоль канала не изменяется и с течением времени устанавливается, что еще раз подтвердило невозможность моделирования переходных и турбулентных течений на основе уравнений Навье-Стокса.
Среди работ по моделированию процессов массо- и теплопереноса отметим также направление, связанное с изучением волновых процессов в газовых средах, как неподвижных, так и движущихся. Эти исследования находятся в русле бурно развивающейся в последние годы отрасли науки — вычислительной аэроакустики, ориентирующейся на применение неупрощенных моделей явлений, основанных на гидродинамических уравнениях Эйлера и Навье-Стокса.
Рис. 4. Поля направлений вектора скорости при течении в канале с числом Рейнольдса Но = 700 и угле атаки 45° в моменты времени < = 10 (а); 20 (б); 30 (в); 40 (г); 50 (с?); 100 (е)
Одной из исследованных в лаборатории в последние годы проблем была проблема резонатора Гарт-мана. Еще в начале XX в. было экспериментально обнаружено, что при помещении препятствия в виде полого цилиндра в сверхзвуковую недорасширенную струю (давление на ее начальном участке превышает атмосферное) могут наблюдаться периодические пульсации, причем лишь при определенных положениях преграды в струе. Изучению этого автоколебательного процесса, нашедшего применение в ряде технических устройств, было посвящено большое количество работ, как экспериментальных, так и теоретических. Однако упрощенные теоретические схемы не могли описать данное явление во всей его полноте, и стало очевидным, что ее решение следует искать лишь на путях численного моделирования. Численное решение задачи о взаимодействии сверхзвуковой струи с введенным в нее препятствием решалось на основе уравнений Эйлера с использованием разностных схем первого (Годунова), второго и третьего порядков [17-19]. Анализ полученных решений показал, что они воспроизводят автоколебательный процесс, причем именно в тех режимах взаимодействия струи с преградой, в которых он наблюдался в экспериментах. Количественное сопоставление расчетных и экспериментальных данных показало также их хорошее согласование по амплитуде и частоте пульсаций. Анализ полей течений позволил провести качественное описание процесса автоколебаний, в котором важную роль играют периодическое затекание определенной массы газа в полую преграду с последующим ее выплескиванием. Проведено обширное параметрическое исследование процесса в зависимости от динамических параметров струи и геометрии преграды. Пример, демонстрирующий характер полученных численных решений, представлен на рис. 5.
50 100 150 200 250 300 t
Рис. 5. Периодический характер давления по дну полого цилиндра длины AD, обтекаемого сверхзвуковой недорасширенной струей (М = 2, ре/ра = 1,2), при расстоянии между срезом сопла и устьем цилиндра 2,25D]
D — диаметр сопла
Близкой по характеру и также весьма сложной проблемой аэроакустики струй является проблема дискретного тона, излучаемого сверхзвуковой струей. Еще в 50-х годах XX в. было экспериментально обнаружено, что на фоне сплошного шума, производимого струей при ее истечении в атмосферу, при определенных условиях могут появиться дискретные составляющие большой интенсивности. Теоретические исследования тонального шума струй ранее проводились на основе линейных моделей [20]. Показано, что линейные решения способны описать ряд особенностей изучаемого явления, находящих подтверждение в эксперименте, включая, например, такой интересный эффект, как образование "зон молчания" в акустическом поле струи. Это демонстрирует рис. 6, на котором построены расчетное и экспериментальное (ИТПМ РАН) распределения фаз пульсаций давления во внешнем акустическом поле сверхзвуковой недорасширенной струи (М = 1,5, ре/ра = 1,5) по длине струи. Характер расчетной и экспериментальной кривых одинаков, хотя имеет место некоторый их сдвиг относительно друг друга вдали от среза сопла. Точки А в расчете и В в эксперименте соответствуют скачкообразному изменению фазы на 180°. В таких точках амплитуда пульсаций давления обращается в нуль, а геометрическое место таких точек и образует "зону молчания" (заштрихованная полоса в нижней части рисунка).
В то же время линейные модели в силу своей природы ограничены необходимостью привлекать априорные предположения для замыкания петли обратной связи в автоколебательном процессе. Моделирование истечения звуковых недорасширенных струй в атмосферу на основе нелинейных уравнений Эйлера позволило воспроизвести автоколебательный процесс без внесения каких-либо дополнительных гипотез, хотя генерирующие процесс возмущения имели не физическую, а вычислительную природу [21].
Также проводилась работа по моделированию распространения ударных и взрывных волн в атмосфере и их взаимодействия с различными препятствиями, в частности фокусировки волн в поло-
0; 360°
0; 360°
У о 2 -
-1-1_I_ |_|_
0 2 4 6 5 10
1
Рис. 6. Сравнение расчетных и экспериментальных данных по тональному излучению сверхзвуковой струи
стях различной геометрии. Эта проблема связана с необходимостью оценки масштабов и характера разрушений в замкнутом пространстве. Очевидно, вычислительный эксперимент был бы наиболее результативным при прямом моделировании взрыва в интересующем исследователей объеме и задании начальных и граничных условий, соответствующих конкретной ситуации. Однако такое моделирование очень трудоемко, а в условиях сложной пространственной геометрии практически невыполнимо. Сложность задачи заключается, прежде всего, в необходимости включить в расчетную модель очень большой объем информации в качестве определяющих параметров, а также в неустановленности иерархической значимости этих параметров. Значимость каждого параметра зависит от конкретных условий моделируемого режима, и в ее определении может помочь расчет достаточно большого количества режимов с набором типичных начальных и граничных условий.
Поэтому целесообразно проведение численного исследования процессов вхождения ударных волн в клиновидную и коническую полости. Полученные в расчетах распределения параметров в объеме и на ограничивающих объем поверхностях дают возможность определять, задание каких параметров является наиболее существенным при оценке возможных последствий взрыва и анализе степени достоверности прогнозов.
В рамках модели идеального совершенного газа (уравнения Эйлера) методом Годунова были проведены расчеты кумуляции взрывных и плоских ударных волн в конических и клиновидных полостях при изменении параметров волн от близких к звуковым (число Маха падающей ударной волны 1,05) до волны с числом Маха 2,4, за которой имеет место сверхзвуковой поток [22-26]. Для использования результатов расчетов при оценках развития процессов в физических экспериментах и в натурных условиях рассматривалось несколько способов задания параметров взрывных волн, входящих в полость, и варьировались внутренняя и внешняя геометрии объемов, в которых происходит их распространение. Рассчитывалось вхождение в полость плоской или взрывной волны. Взрывная волна рассматривалась как результат мгновенного разлета (расширения) области газа (относительно малого размера) с повышенными параметрами.
Максимум давления на волне при ее входе в полость соответствовал давлению, которое в физическом эксперименте имеет возникшая при мгновенном выделении заданной энергии волна на расстоянии, соответствующем входу в полость. В качестве характеристики кумуляции в полости принималось значение степени фокусировки, определяемое как отношение максимального значения параметра, полученного внутри объема в процессе расчета, к значению, которое получилось бы при нормальном отражении от плоской стенки ударной волны с параметрами, равными локальным значениям их на входе в полость. Так, например, для волны, у которой при входе в полость давление на фронте равно 1,8, соответствующая волна с плоским фронтом имеет давление отражения, равное 3,1.
Анализ результатов расчетов показывает, что для слабой волны имеет место заметное понижение степени фокусировки в вершине, а максимальные значения давления в привершинной области конуса
примерно в три раза больше, чем для клина. Расчеты вхождения плоской волны в полость сопоставлялись с экспериментальными результатами. Характер изменения результатов кумуляции в зависимости от геометрии полости и числа Маха входящей плоской волны сходны. Обнаружены особенности распределения параметров в зоне фокусировки в зависимости от конкретных условий. Максимальные значения параметров не всегда реализуются на оси или на стенке полости. Сопоставление кумуляции в клиновидной и конической полостях показывает, что характер изменения процесса кумуляции с изменением геометрии фокусирующей преграды при разной интенсивности волн сохраняется.
Получена картина взаимодействия волн с поверхностями, ограничивающими полости. Получены поля давления внутри и вне полостей на разных стадиях течения и изменение давления, плотности и температуры в фиксированных точках поверхностей (датчиках) (рис. 7, 8).
Рис. 7. Отражение взрывной волны от внутренней поверхности конуса в моменты времени £ = 0,8 (а); 1,3 (б). Картина изобар (угол раствора конуса 90°; энергия инициирования сферической волны 19 кДж)
Рис. 8. Изменение параметров течения в вершине конуса по времени (угол раствора конуса 60°; энергия инициирования сферической волны 8,5 кДж)
При сопоставлении и анализе результатов выявлены особенности кумуляции взрывных и плоских волн в конических и клиновидных полостях для разной формы внешнего расчетного поля. Картина
фокусировки изменяется не только количественно, но и качественно. Характер изменения зависит от конкретной формы преграды и интенсивности волны.
Также были исследованы процессы, возникающие при набегании волнового фронта на природное препятствие в виде массива растительности [27]. Моделирование проводилось в рамках единого методического подхода, использующего уравнения Эйлера для описания движения масс воздуха как на открытой местности, так и в пределах массива. В последнем случае в уравнения вводятся массовые силы, действующие со стороны растительности. Численное решение задачи получено методом Годунова. Исследованы два вида набегающего волнового фронта: плоская ударная волна и нелинейный акустический импульс, моделирующий фронт сферической взрывной волны на больших расстояниях от источника. Для различных типов растительности (хвойный, лиственный лес) исследованы особенности процесса взаимодействия, включающего проникновение фронта в массив, частичное отражение от передней границы массива и дифракцию над верхней границей. Результаты численных решений указывают на образование пары "восходящий поток-нисходящий поток" в верхней части массива (в кронах деревьев). Существование этой структуры находит подтверждение в результатах экспериментальных исследований (рис. 9).
Рис. 9. Распределение давления при набегании (справа) ударной волны (М = 3) на массив хвойного леса. Граница массива при у = 10, высота леса Н = 1, фронт волны находится в положении у = 19,5
Была продолжена работа по численному исследованию течений газа в каналах и соплах [28, 29]. Был разработан модифицированный, полностью неявный метод (МПНМ) решения системы алгебраических уравнений, возникающих при симметричной аппроксимации на 19-точечном пространственном шаблоне краевой задачи для линейного дифференциального уравнения в частных производных эллиптического типа. Метод был применен к расчету пространственных безвихревых течений газа в соплах. Такие течения описываются квазилинейным уравнением потенциала скорости смешанного типа, поэтому перенос численных методов, разработанных для линейных эллиптических уравнений, на уравнение потенциала не является тривиальным. К основным достоинствам нового алгоритма относятся высокая скорость сходимости и приемлемая точность. Для расчета смешанных течений, содержащих большие сверхзвуковые скорости, предложен алгоритм, представляющий собой комбинирование МПНМ и метода приближенной факторизации. С помощью реализованных алгоритмов показана возможность построения сопел с прямоугольным выходом, по тяговым характеристикам близких к осесимметричным соплам.
Перечисленные выше задачи так или иначе связаны с работой различных технических устройств. В то же время численное газодинамическое моделирование с успехом используется для решения задач астрофизики. В лаборатории в разные годы исследовались проблемы взаимодействия потоков вещества в космическом пространстве: соударение солнечного ветра с потоком межзвездной среды, обтекание солнечным ветром ионосфер комет, взаимодействие звездных ветров в тесных двойных звездах.
Применительно к астрофизическим проблемам газодинамические уравнения приходится усложнять за счет учета различных физических и химических процессов, характерных для условий космо-
са. Так в проблеме обтекания динамической ионосферы кометы потоком солнечного ветра уравнения движения сплошной среды могут быть применены лишь к заряженной компоненте вещества (плазме), тогда как течение нейтральной компоненты рассматривается как свободно-молекулярное. При этом вследствие ионизации нейтральной компоненты различными механизмами в уравнениях появляются новые члены, описывающие источники массы, импульса и энергии.
Полученные на основе такого подхода численные решения хорошо согласуются — качественно и количественно — с результатами измерений, проведенных на борту космических аппаратов. Это касается как результатов, полученных при пролете аппаратов у кометы Галлея еще в 1986 г., так и более поздних данных, полученных в первой половине 90-х годов для кометы Григга-Шеллерупа ([30], рис. 10).
V, км 400 Н
300 -
200 I-11
400 -
300 L 16
Рис. 10. Распределение скорости солнечного ветра по траектории космического аппарата Giotto при его прохождении вблизи кометы Григга-Шеллерупа (сплошные линии — расчет, точки — эксперимент на борту Giotto (A.Johnstone et al.), SCET — время в системе отсчета, связанной с Giotto, 10 июля 1992 г.)
К числу таких результатов принадлежит также обнаруженное сначала в расчетах и подтвержденное в экспериментах наличие протяженной, почти застойной области, в которой плотность вещества меняется на несколько порядков ("мантия тяжелых ионов", по терминологии астрофизиков). Учитывая, что две названные кометы являются весьма разными по критерию газопроизводительности, можно сказать, что разработанные модели пригодны для широкого класса комет.
СПИСОК ЛИТЕРАТУРЫ
1. Кокошинская Н. С., Павлов Б. М., П ас ко н о в В. М. Численное исследование сверхзвукового обтекания тел вязким газом. М.: Изд-во Моск. ун-та, 1980.
2. Белова О.Н., Кокошинская Н.С., Лузянина Т.Б., Пасконов В.М. Численное моделирование ламинарного обтекания тел вязким газом. М.: Изд-во Моск. ун-та, 1986.
3. Пасконов В.М., Петухова Т.П., Русаков C.B. Численное моделирование нестационарных ламинарных течений вязкого газа. М.: Изд-во Моск. ун-та, 1986.
4. Русаков C.B. Разностные сплайн-схемы для задач тепло- и массопереноса. Иркутск: Изд-во Иркутского ун-та, 1990.
5. Пасконов В. М., Русаков C.B., Чудов И. И. Об одной разностной сплайн-схеме для решения уравнений Навье-Стокса в естественных переменных // Численные методы в математической физике. М.: Изд-во Моск. ун-та, 1996. С. 33-38.
6. Paskonov V.M., Rusakov S.V., Chudov I.I. Difference spline scheme for Navier-Stokes equations in natural variables // Comp. Math, and Modeling. 1997. 8. N 2. P. 112-117.
7. Пасконов В.М., Фортова C.B. Течение вблизи кормы продольно обтекаемого тела при поперечном вдуве с его боковой поверхности // Изв. РАН. Механика жидкости и газа. 1996. № 2. С. 157-163.
8. Paskonov V.M., Fortova S.E. Flow around the afterbody of a plate in a longitudinal gas stream with transverse injection from the lateral surface // Fluid dynamic. 1996. 31. N. 2. P. 296-301.
9. Пасконов В. M., Фортова C.B., Мурашкина К. Б. Численное исследование течения в следе за пластиной конечной толщины. Проблемы математической физики. М.: Диалог МГУ, 1998. С. 184-192.
10. Пасконов В. М., Фортова С. В., Мурашкина К. Б. Численное исследование течения в следе за пластиной конечной толщины при сверхзвуковых вдувах с ее боковой поверхности // Вестн. Моск. ун-та. Сер. 15. Вычисл. матем. и киберн. 1998. № 3. С. 42-46.
11. Пасконов В.М. Модификация уравнений Навье-Стокса для расчета вязких течений разностными методами // Обратные задачи естествознания. М.: Изд-во МГУ, 1997. С. 189-198.
12. Paskonov V.M. Modified Navier-Stockes equation for finite-difference computation of viscous flow / / Comp. Math, and Modeling. 1997. 8. N 4. P. 400-407.
13. Мурашкина К.Б., Пасконов В.М. Численное исследование краевой задачи для нелинейного уравнения типа Кортевега-де Фриса-Бюргерса // Прикладная математика и информатика. № 10. М.: МАКС Прес, 2002. С. 5-14.
14. Murashkina К.В., Paskonov V.M. Numerical analysis of the boundary-value problem for a nonlinear Korteweg-de Vries-Burgers equation // Computation mathematics and Modeling. 2003. 4(2). April-June. P. 85-92.
15. Бе рези h С. Б., Пасконов В.М. Численное исследование вдува вязкого несжимаемого газа в плоский канал на основе уравнений Навье-Стокса // Вычислительные методы и программирование. НИВЦ МГУ. 2003. 4(1). С. 5-17.
16. Бе резин С. Б., П ас ко н о в В. М. Неклассические решения классической задачи о течении вязкой несжимаемой жидкости в плоском канале // Прикладная математика и информатика. № 17. М.: МАКС Пресс, 2004. С. 13-30.
17. Базаров С.Б., Росляков Г.С., Садков Ю.Н., Шустова М.В. Сверхзвуковые нерасчетные струи и их взаимодействие с преградами // Мат. моделирование. 1996. 8. № 6. С. 19-23.
18. Росляков Г. С., Садков Ю.Н. Автоколебания при взаимодействии сверхзвуковых течений с преградами // Прикл. математика и информатика. № 3. М.: МАКС Пресс, 1999. С. 71-79.
19. Росляков Г. С., Садков Ю.Н. Осесимметрические течения идеального газа и их взаимодействие с преградами // Прикл. математика и информатика. № 3. М.: МАКС Пресс, 2001. С. 57-62.
20. Лебедев М.Г. Линейные решения в задаче об излучении дискретного тона сверхзвуковыми струями // Прикл. математика и информатика. № 10. М.: МАКС Пресс, 2002. С. 70-89.
21. Садков Ю.Н. Численное моделирование сверхзвуковых осесимметричных газовых струй // Прикл. математика и информатика. № 15. М.: МАКС Пресс, 2003. С. 39-45.
22. Андрианов А.Н., Базаров С.Б., Ефимкин К. Н. Решение двумерных задач газовой динамики методом Годунова на параллельных ЭВМ с помощью языка НОРМА. Препринт ИПМ. № 9. М., 1997.
23. Andrianov A.N., Bazarov S.B.,Bugerya A.B., Efimkin K.N. Solution of three-dimensional problems of gas dynamics on multiprocessors computers // Computational Mathematics and Modeling. 1999. 10. N 2. P. 140-150.
24. Базаров С.Б., Набоко И. M. Моделирование кумуляции волн в полостях // Прикл. математика и информатика. № 8. М.: МАКС Пресс, 2001. С. 55-70.
25. Базаров С. Б., Набоко И.М. Моделирование взаимодействия волн с полостями // Мат. моделирование. 2002. 14. № 9. С. 34-40.
26. Базаров С.Б., Гусев П.А., Набоко И.М., Петухов В.А., Солнцев О. И. Кумуляция взрывных волн в конусе в физическом и вычислительном экспериментах // Научные труды Института теплофизики экстремальных состояний ОИВТ РАН. Вып. 5. М.: Изд-во ОИВТ, 2003. С. 333-338.
27. Лебедев М.Г., Ситник В. В. Моделирование взаимодействия волнового фронта с массивом растительности // Прикл. математика и информатика. М.: МАКС Пресс (в печати).
28. Росляков Г. С., Ф е д о р е н к о В. В. Исследование пространственного течения газа в сопле полностью неявным методом // Методы математического моделирования. М.: Изд-во МГУ, 1998. С. 7686.
29. Федоренко В. В. Полностью неявные методы расчета пространственных течений газа в соплах // Мат. моделирование. 2000. 12. № 4. С. 83-104.
30. Lebedev M.G. Comet Grigg-Skjellerup atmosphere interaction with the oncoming solar wind // Astrophys. Space Sei. 2000. 274. N 1-2. P. 221-300.
Поступила в редакцию 10.08.04
УДК 519.632+517.518 Е. Г. Дьяконов
АППРОКСИМАЦИИ В УСИЛЕННЫХ ПРОСТРАНСТВАХ СОБОЛЕВА ДЛЯ НЕКОТОРЫХ ТРЕХМЕРНЫХ ОБЛАСТЕЙ С НЕРЕГУЛЯРНОЙ ГРАНИЦЕЙ
(кафедра общей математики факультета ВМиК)
1. Введение. Настоящая работа тесно связана со статьей [1], в которой рассматривались аналогичные вопросы для более узких энергетических пространств и более простых трехмерных областей. Рассматриваемые усиленные пространства Соболева С1'1 строятся на базе классического пространства Ил21(<5) = Н1(С}) для случая модельной ограниченной области (5 С И.3 с, вообще говоря, нелип-шицевой кусочно-гладкой границей Г. Основное внимание уделяется, например, случаю, когда два гладких двумерных блока, входящих в состав Г, могут касаться и образовывать внутренний (относительно (5) двугранный нулевой угол. В работе, следуя [1-7], анализируются некоторые неоднородные модификации (усиления по отношению к пространству Н1{Ц)) С1'1 = С1'1 (Х'3'; X'2'), использующие наряду с основной трехмерной структурой (блоком) X'3' = (5 еще и двумерные каркасы X'2' (многообразия, подструктуры), на которых требуется дополнительная по отношению к ^1(Х'3') гладкость следов (отметим, что конструкции этих относительно новых энергетических пространств имеют родство с известными инженерными приемами усиления жесткости исходной упругой структуры, аналоги которых легко могут быть найдены и в объектах, создаваемых самой природой). В зависимости от геометрии х<3) получены прямые и обратные теоремы о следах элементов из С1'1 на X'2' (нормальная обратимость оператора следа) даже без привлечения соответствующих теорем о следах для элементов из ^1(Х'3'). На их основе доказана и центральная теорема об аппроксимациях элементов этих необычных энергетических пространств при помощи гладких функций на Х'3', имеющая фундаментальное значение как для построения аппроксимаций типа сплайнов, так и для теории проекционно-сеточных методов решения стационарных и спектральных задач в пространствах типа С1'1. Можно сделать вывод, что наличие подходящих усилений в сингулярном блоке в некотором смысле компенсирует его плохую геометрию.
2. Усиленные пространства Соболева С в случае модельного трехмерного
сингулярного блока. Ниже основной трехмерный блок х<3) = (5 берется таким, чтобы, с одной стороны, можно было продемонстрировать специфику учета наличия сингулярных точек границы (точек границы, в окрестностях которых не выполняется условие Липшица), а с другой стороны, чтобы и геометрические описания не стали чересчур громоздкими. С этой целью мы берем блок (5 гомеоморфным треугольной призме, изображенной на рисунке; его граница разбивается на пять гладких блоков — образов граней призмы, обозначенных на рисунке. Эти блоки называем квазигранями; их пересечения
