9. Гидродинамические задачи

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

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

В задачах с числом измерений, большим 1, возникают серьезные трудности в планировании вычислительной работы. Дифференциальные уравнения Эйлера и Лагранжа могут быть, конечно, заменены аппроксимирующей системой разностных уравнений. Непрерывно меняющееся время также может быть заменено дискретной последовательностью моментов. Практические трудности выполнения таких вычислений весьма значительны: если в задаче возникают большие смещения жидкости или газа (большие по сравнению с начальными размерами системы) и искажения формы значительны, то начальная сетка лагранжевых точек со временем изменяет свою структуру и весьма часто перестает служить адекватным базисом для приближенного представления вычисляемых дифференциальных выражений (например, для оценки градиентов плотности или давления).

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

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

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

Рассмотрим процесс, имеющий менее «ортодоксальный» характер. Можно предложить по меньшей мере два различных

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

Ниже дается краткое описание работы, проведенной в этом направлении на счетных машинах Паста и автором. За этой работой следовали дополнительные вычислительные эксперименты такого же рода. Развитие метода в основном принадлежит Харлоу [1] из Научной лаборатории Лос Аламоса.

Можно пытаться обосновать гидродинамические вычисления с помощью кинетической модели газа или жидкости. Принимая за модель физической реальности интегро-дифференциальное уравнение Больцмана в гидродинамике, можно вычислить свойства движения «точек» — атомов или молекул, которые имеют данное распределение скоростей, в целом. Эти точки подвержены флуктуациям, и согласуются с макроскопическим движением только их статистические средние.

Хорошо известно, например, что из этой модели можно вывести уравнения Навье — Стокса. Возникает вопрос: реально ли на этой основе вычислить макроскопические величины, такие, как плотность, давление и скорость, в зависимости от времени. В настоящее время ответ на этот вопрос представляется нам отрицательным. Если подходить к уравнениям Больцмана «буквально», считая, что отдельные точки представляют атомы, то для получения средней скорости движения газа в некоторой точке пространства нужно потребовать, скажем, частиц на клетку в интересующем нас пространстве решений. Пусть число клеток при вычислении равно Если на линейной шкале измерений число интервалов равно то где число измерений в данной задаче (для задач в одно-, двух- или трехмерном пространстве, не обладающих специальной симметрией, , а общее число точек в вычислении должно быть равно

При статистическом исследовании, если требовать, чтобы средняя ошибка была порядка 1%, нужно считать Если линейная «разрешимость» должна при этом быть, скажем, 5%, то Даже в одномерной задаче придется рассматривать 200 000 точек. Это намного превосходит возможности современных машин, и поэтому ясно, что вычисления должны вестись с «точками», представляющими не индивидуальные атомы, а довольно большие группы атомов.

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

В этом общем описании мы можем предположить процесс, например, изотермическим или адиабатическим, т. е. задать Градиенты давления будут тогда определяться через градиенты плотности, а эти последние будут вычисляться через полученные одновременно положения системы точек, являющихся, по нашему предположению, центрами масс частиц жидкости. Задача состоит тогда в нахождении правила или рецепта для оценки плотности в точке пространства, если дана только конечная система точек (на практике довольно далеких друг от друга). В предельном случае очень большого числа точек можно просто сосчитать число точек в каком-либо квадрате или кубе решетки и это число будет пропорционально плотности. На практике у нас вообще имеется не очень большое число частиц во всей жидкости, и вопрос состоит в том, как при этом наиболее надежно оценить плотность. Рассмотрим случай двух измерений. Можно считать, что вначале точки расположены каким-нибудь регулярным способом, например в вершинах прямоугольного разбиения или, что лучше, в вершинах разбиения пространства на треугольники. Конечно, после нескольких вычислительных циклов геометрия пространства изменится и полученную плотность можно оценить, например, следующим образом: окружим каждую точку треугольниками, вершинами которых являются ближайшие к ней точки. Отношение площади наименьшего треугольника, внутри которого расположена данная точка, к площади исходного треугольника даст по крайней мере представление об изменений плотности в этой точке.

У этого процесса есть несколько недостатков. Они происходят из-за некоторой «неустойчивости» наших вычислений. Выбор ближайших точек приводит к разрывности

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

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

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

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

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

Силы (отталкивающие) между двумя произвольными точками

где расстояние между двумя точками.

Среднее значение в точке жидкости является в пределе при функцией от локальной плотности: для трех измерений.

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

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

Чтобы испытать некоторые из этих предложений на конкретных задачах, были проведены вычисления на электронной машине МАНИАК в Лос Аламосе и на ее прототипе в Institute for Advanced Study в Приистоне (ср. Паста и Улам [1]). Главной целью этих вычислений было скорее испытать осуществимость таких вычислительных схем, чем получить количественную оценку результатов. Одна из задач относится к изучению движения тяжелой жидкости, расположенной над слоем более легкой, причем обе жидкости расположены в кубе и заданы начальные нерегулярности на границе. Задача состоит в прослеживании этого неустойчивого положения (известного под названием тейлоровской

неустойчивости) для движения в целом (а не только в бесконечно малом), ведущего к постепенному перемешиванию обеих жидкостей.

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

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

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

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

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

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

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

Со временем происходит перемешивание и наша средняя мера перемешивания возрастает. Снова можно ожидать, что после достаточного количества вычислительных экспериментов станет более ясным поведение этой «функции перемешивания» в зависимости от времени или по крайней мере определится время нужное для того, чтобы общее перемешивание (которое первоначально близко к нулю) стало порядка или Это время из соображений размерности должно зависеть inter alia от где глубина сосуда и ускорение более легкой жидкости в более тяжелую.

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

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

Харлоу и Эванс [1] произвели ряд вычислений, используя метод, несколько отличающийся от нашего, но также

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