книги / Метод крупных частиц в газовой динамике

2) итерации, обеспечивающие введение давления в уравнение импульсаг

ивведение сжимаемости в уравнение неразрывности,

3)расчет конвективных потоков, которые учитываются, если сетка не чисто лагранжева.

Вычисления проводятся в двумерной (плоской или осесимметричной) рас­ четной области, которая покрывается четырехугольной сеткой. Типичная:

счетная ячейка показана на рис. П.4.1. Отметим, что в узлах сетки, как и в методе ALE, определяются координаты х, у , масса М, компоненты скоростиг иу v; в центре яйчеки определяется плотность р, давление р, объем ячейки V,. полная и внутренняя энергии £ , 3. Подробные разностные формулы содержа­ тся в [45]. Основные элементы данного численного подхода были обсуждены

⧧ 2—3 Приложения.

2.Приведем результаты численного расчета, иллюстрирующие возмож­ ности данной методики.

Рассмотрим задачу о развитии взрыва малой интенсивности [45]. Перво­ начальное (при t=0) расположение частиц-маркеров, конфигурация счетной: сетки, поля вектора скорости, изопикнические линии p=const, изотермы и линии равной завихренности изображены на рис. П.4.2. Число расчетных узлов равнялось 26 x 52. На рис. П.4.3 приведены соответствующие данные лагран-

жевых вычислений, полученные методом YAQUI для момента времени t= l .

§ 5. Метод L1NC для расчета неустановившихся несжимаемых течений жидкости

Метод LINC (Lagrangian Method for Incompressible Flow) [46], основанный на использовании лагранжевых координат, позволяет рассчитывать слабо де­ формируемые нестационарные потоки несжимаемой жидкости со свободными границами *). В отличие от методов эйлерова типа, где жидкость движется через фиксированную расчетную сетку, метод LINC базируется на лагранже­ вых координатах (расчетные ячейки здесь перемещаются вместе с жидкостью).

Следует отметить, что лагранжевы уравнения движения решаются здесь аналогичным способом, как и в методе ALE [44], описанном в § 3 Приложения.

Уравнения импульса записываются в тензорной форме

(П .5.1)

где по повторяющимся индексам производится суммирование; DIDt — произ­ водная по времени вдоль траектории данного элемента, р — плотность, и — скорость, р — давление, П — девиатор тензора напряжений, А — ускоре­ ние за счет массовых сил. Систему (П.5.1) дополняют уравнения для лагранжевой переменной

(П.5.2)

где г* — радиус-вектор.

Для решения системы (П.5.1) — (П.5.2) используется лагранжева четы­ рехугольная сетка. Параметры потока определяются в узлах сетки и в центрах ячеек аналогично тому, как это делалось в методе ALE. Давление в каждой ячейке определяют из условия сохранения объема, что приводит к уравнению типа Пуассона, решаемого итерационным методом Гаусса-Зейделя.2*

*) См. также Н i г t С. W., C o o k J. L., В u 1 1 е г Т. D. A Lagrangian Method for Cal­ culating the Dynamics of an Incompressible Fluid with Free Surface.—В сб.: Труды секции по численным методам в газовой динамике второго международного коллоквиума по газодинамике взрыва и реагирующих систем. (Новосибирск, 19—23 августа 1969.) — М.: ВЦ АН СССР, 1971, 2, с. 188—227.

Расчет каждого временного цикла разбивается на следующие этапы:

1) в каждой вершине четырехугольной ячейки рассчитываются величины ускорений (кроме создаваемых за счет сил давления),

2)на основании данных этапа 1 вычисляются члены типа «источника», входящие в уравнение Пуассона для давления,

3)проводится расчет поля давлений,

4)определяются новые значения скоростей и перемещений вершин ячеек,

5)делается шаг по времени.

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

многих схемах лагранжева типа из-за коротковолновых возмущений. Указан­ ные явления приводят к искажению расчетных ячеек, которые принимают лри этом форму песочных часов. В обычной схеме LINC не имеется каких-либо механизмов, тормозящих такие передвижения вершин, что может привести к резким искажениям сетки и потере точности. С этой целью в работе [46] опи­ сывается некоторый механизм регуляризации расчетной сетки, позволяющий замедлить искажение ячеек. Для этого рассматривается демпфирующая сила, которая стремится удержать вершину ячейки в центре масс восьми соседних вершин и определяется по закону Гука. Для этого в уравнении (П.5.1) к вели­ чине A t добавляется член

A'i = k2 (г•— r{) + 2k (и\— и(),

где г' и м'— радиусы-векторы центра массы и скорости восьми соседних вер­ шин, k — некоторый параметр, зависящий от конкретной задачи.

Эффективность такой процедуры показана на рис. П.5.1, где приведен расчет колебаний однородной невязкой жидкости в сосуде (в момент t= 0 на

Рис. П.5.1. Расчетная сетка в методе LINC (задача о расчете колебаний невязкой жидкости в сосуде) [46]. а) Первоначаль­ ная конфигурация; б) расчетная сетка без регуляризации; в)

расчетная сетка при наличии регуляризации.

свободной поверхности жидкости задается возмущение: косинусоидальный им­ пульс давления) [46]. На рис. П.5.1, а показана первоначальная конфигурация разностной сетки, на рис. П.5.1, б и в приводятся конфигурации расчетных сеток соответственно в момент t= 21 без регуляризации и при наличии демпфи­

рующей силы.

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

Здесь Т — коэффициент поверхностного натяжения, R — радиус кривизны поверхности.

Контрольные расчеты показали, в частности, что использование в (П.6.9) Q/fy вместо R if/ несущественно влияет на результаты вычислений.

Таким образом, в каждом вычислительном цикле можно выделить следую­ щие основные этапы:

а) поле скоростей предполагается известным либо как начальное прибли­ жение, либо как результат предыдущего цикла. Поле скоростей соленоидально, т. е. D I /= 0;

б) соответствующее поле давления находится из решения уравнения Пуассона (П.6.9);

в) подставляя полученное поле давления в уравнения (П.6.4) и (П.6.5), находим новое поле скоростей, которое также соленоидально, и т. д.

Если в задаче отсутствуют свободные поверхности, то основная часть вычислительного процесса на этом заканчивается, в противном случае вво­ дятся меченые частицы. Следует еще раз отметить, что, в отличие от PICметода, маркеры здесь не участвуют в вычислениях, а только показывают, например, какие ячейки содержат жидкость (они выделяют свободную поверх­ ность). При этом на первом этапе вычислений предполагается, что координаты набора маркеров известны. Вычислив новое поле скоростей, рассчитывается новое положение свободной поверхности и т. д.

3. Рассмотрим теперь постановку граничных условий.

Ось (плоскость) симметрии: нормальная компонента скорости меняет знак, касательная остается неизменной. Условия для давления получаются из уравнений (П.6.4) и (П.6.5) в зависимости от положения оси (плоскости) симметрии. В случае вертикальной оси из (П.6.4) имеем ф'= Ф 1±£*Дл:, где знак + для случая, когда жидкость находится слева от оси, а —, когда справа. В случае горизонтальной оси ф/ = ф1± ^ Д ^ , где знак + берется, когда жидкость под осью, а —, когда жидкость находится над ней.

Твердая стенка: нормальная компонента скорости остается неизменной, касательная меняет знак. Для давления: вертикальная стенка ф '= ф 1±|&Д*:Е

3 6 8 ПРИЛОЖЕНИЕ

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

4. В качестве примера расчета методом MAC рассмотрим задачу о движении воды под воротами шлюза [47]. Конфигурации потока в последовательные мо­ менты времени приведены на рис. П.6.2. Вначале жидкость покоится, причем ее уровень слева от ворот выше, чем справа. Расчеты проводились на сетке 50x32, длина расчетной области составляла 4,8, высота —-3,0. Коэффициент вязкости v был взят равным 0,01, g = —1,6. На дне резервуара ф =2,5. Рис. П.6.2 отчетливо иллюстрирует динамику'развития течения.

§ 7. SMAC-метод для расчета течений вязкой несжимаемой жидкости

Взаключение остановимся еще на одном численном методе, разработан­ ном под руководством .Харлоу, для расчета неустановивигихся течений вязкой несжимаемой жидкости — SMAC-методё (A Simplified MAC Method) [50].

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

Вметоде SMAC давление вычислять не обязательно. Граничные условия также-упрощены: они ставятся в уравнении импульса только для скорости и для напряжения давления на свободной поверхности (нормального и танген­ циального). В то же время решение уравнения Пуассона (для сохранения массы) требует везде только однородных граничных условий. Основная идея метода заключается в расщеплении каждого временного цикла на две фазы.

1)Вычисляется промежуточное поле скоростей un+1, vn+1 при использо­ вании произвольного поля давления. При этом граничные условия для дав­ ления на свободной поверхности выполнены лишь для нормальных напряже­ ний. Например, если пустая ячейка находится слева или справа от данной поверхностной ячейки, то

Ф /. / = Фг, / (внеш) + (2 V/Д ^ ) (VU /+ , / , — Vi%l - 1/2) .

Соответствующие граничные условия для скорости обеспечивают нужную за­ вихренность поля для промежуточных значений скоростей во всех внутрен­ них точках. Однако эти промежуточные значения не удовлетворяют условию неразрывности, т. е. Ь ,^ Ф 0.

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

Промежуточное поле скоростей ц“+1, vn+l при этом изменяется путем до­ бавления градиента потенциальной функции ф, определяемой из условия вы­

370 ПРИЛОЖЕНИЕ

Как метод MAC, так и SMAC применимы при изучении как плоских, так и осесимметричных течений в различных системах координат и для различных конфигураций областей; возможно распространение этих подходов и на слу­ чай пространственных течений. Подробное описание расчетных формул ме­ тода SMAC приводится в [50].

В качестве примера использования SMAC-метода можно привести эффек­ тивный расчет падающего столба воды. На рис. П.7.1 приводятся положения частиц (левый столбец) и вектора скорости (справа) для последовательных моментов времени (0 ^ ^ 1 3 ,5 ). На стенках задавалось условие свободного проскальзывания.