Методы контрольных объемов

8.4.3. Метод С.К. Годунова

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

Метод С.К. Годунова свободен от этого недостатка. Параметры газа на границах контрольных объемов в этом методе устанавливаются автомодельным решением линеаризованной задачи газовой динамики (см. главу 7) [Численные]. Для различных приложений используются различные автомодельные решения. Ниже рассмотрим вариант метода С.К. Годунова в применении к задачам о течении идеального газа. В этих задачах для определения параметров на границах контрольных объемов используется автомодельное решение задачи о распаде произвольного разрыва в инертном газе [Рождественский].

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

- в обе стороны (из области 1 в область 2 и из области 2 в область 1) распространяются ударные волны , разделенные контактным разрывом k (рисунок 8.20а);

- справа налево (в область 1) распространяется ударная волна d, за ней следует контактный разрыв k, а в область повышенного давления (область 2) распространяется волна разрежения r (рисунок 8.20б);

- в обе стороны распространяются волны разрежения , разделенные контактным разрывом k (рисунок 8.20в);

- в обе стороны распространяются волны разрежения , разделенные областью нулевого давления (вакуумом)(рисунок 8.20г).

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

1. Вычисляем значения скорости распространения по расчетной области ударной волны , волны разрежения , значение скорости , соответствующей возможности возникновения области с нулевым значением давления (вакуума):

; (8.127)

; (8.128)

. (8.129)

Здесь обозначено ; . Кроме того, для идеального газа можно принять ;

2.Определяем вариант развития распада разрыва параметров после разрушения мембраны:

- в обе стороны от границы раздела двух сред распространяются ударные волны, если выполняется условие

, (8.130)

- вправо от границы раздела двух сред распространяется ударная волна, а влево – волна разрежения, если выполняются условия

, (8.131)

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

, (8.132)

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

; (8.133)

3.Определяем значение давления на контактном разрыве (внутри области распада разрыва) решением нелинейного уравнения

. (8.134)

В (8.134) обозначено

.

Здесь и далее индекс l принимает значения l=1, 2. Уравнение (8.134) может решаться любым численным методом. Однако при использовании этого уравнения для численного решения двух- и трехмерных уравнений газовой динамики важным является условие экономичности метода. В качестве начального приближения для Р при решении уравнения (8.134) можно принять

; (8.135)

4.Определяем значения массовых скоростей для левой (распространяющейся справа налево) и для правой (распространяющейся слева направо) волн

- для ударной волны

, (8.136)

- для волны разрежения

; (8.137)

5. Определяем скорость движения контактного разрыва U

; (8.138)

6. Определяем скорости D левой и правой волн (для волны разрежения – скорости для крайних характеристик) и плотности газа R на этих волнах

- на ударной волне

, (8.139)

, (8.140)

- на волне разрежения

, (8.141)

, (8.142)

, (8.143)

. (8.144)

Рассмотренный алгоритм расчета значений параметров газа на контактном разрыве (уравнения (8.127)…(8.144)) является основой метода С.К. Годунова. Рассмотрим решение уравнений газовой динамики, записанной в нестационарной одномерной постановке

,

,

, (8.145)

.

Решение системы уравнений (8.145) осуществляется на сетке, представленной на рисунке 8.17а. Рассмотрим объемы с номерами I-1 и I. Если рассматривать процессы на границе i объемов I-1 и I I-1 и I в течение периода времени, соответствующее условию Куранта-Фридрихса-Леви ( ), то взаимодействие потоков газа, размещенных по обе стороны от границы i, можно рассматривать как взаимодействие двух потоков, размещенных по обе стороны от мембраны М (рисунок 8.19) в задаче о распаде произвольного разрыва газодинамических параметров в ударной трубе. Параметры на границе i объемов в течение времени от до определяются как параметры на контактном разрыве и могут быть вычислены по алгоритму (8.127)…(8.144)).

С учетом выше записанного явная конечно-разностная схема расчета параметров в объеме с номером I по методу С.К. Годунова может быть записана следующим образом

,

, (8.146)

.

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

8.4.4. Метод крупных частиц

Основой метода крупных частиц является метод частиц в ячейках, изложенный в работах Ф. Харлоу [ВычМетодыГидро]. Основы метода крупных частиц разработаны академиком О.М. Белоцерковским и его научной школой [Белоц].

Метод крупных частиц следует отнести к методам физического расщепления. Рассмотрим основные этапы алгоритма для системы уравнений течения газа в канале переменной площади F, записанной в виде

,

,

, (8.147)

.

Здесь - источниковые члены, F – площадь поперечного сечения канала ().

В соответствии с теорией метода крупных частиц система уравнений (8.147) может быть расщеплена следующим образом:

Эйлеров этап

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

,

,

, (8.148)

;

Лагранжев этап

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

Заключительный этап

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

,

,

, (8.149)

.

Основная конечно-разностная реализация этапов метода крупных частиц для системы уравнений (8.147) в явном виде записывается:

Эйлеров этап

,

; (8.150)

Лагранжев этап

,

,

,

, (8.151)

;

Заключительный этап

,

,

, (8.152)

.

Недостатком основной схемы метода крупных частиц [ДавыдовБелоцерковский] является ее относительно невысокая устойчивость, особенно при расчете областей с малыми дозвуковыми скоростями (числа Маха М<0,3). Ниже излагаются алгоритмы метода крупных частиц, изложенные в [ЧислЭксперимРДТТ] и позволяющие повысить устойчивость численного счета до чисел Куранта Ku 1. Следует отметить, что схемы [ЧислЭксперимРДТТ] успешно используются при решении задач внутренней баллистики в одномерной и многомерной нестационарной постановке, в том числе, и для многокомпонентного газа.

Основные модификации метода крупных частиц, используемые в схемах [ЧислЭксперимРДТТ], состоят в следующем.

1. Анализ показывает, что относительно невысокая устойчивость схемы (8.150)…(8.152) обусловлена неудачной аппроксимацией дифференциальных уравнений, решаемых на Эйлеровом этапе метода. Применение полностью неявных схем на Эйлеровом этапе является решением проблемы повышения устойчивости. Однако такой подход представляется громоздким. В [ЧислЭкспер] решение уравнений (8.150) осуществляется после решения уравнения

. (8.153)

Уравнение (8.153) может быть получено из уравнений (8.147). Его аппроксимацию можно выполнить следующими способами

, (8.154)

. (8.155)

Схемы (8.154), (8.155) явные, и вычисляемые по этим уравнениям значения должны быть использованы вместо значений в уравнениях Эйлерового этапа (8.150). Записанные модификации позволяют увеличить устойчивость явной схемы метода крупных частиц до значений Ku 0,5.

Допустимое значение числа Куранта Ku 0,5 можно обеспечить и в случае, если переписать уравнение (8.153) в виде

. (8.156)

Конечно-разностная аппроксимация уравнения (8.156) позволяет записать уравнение для следующим образом

. (8.157)

2. Применение неявных схем аппроксимации уравнения (8.153) для давления позволяет увеличить допустимое число Куранта до значений Ku 1,0. Предварительно необходимо выполнить следующие преобразования уравнения.

В соответствии с теорией аппроксимации конечно-разностных схем (см. раздел 8.3.2) можно записать

.

В то же время частную производную можно получить из уравнения (8.153)

.

С учетом последних двух уравнений (8.153) может быть переписано в виде

.

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

. (8.158)

Здесь обозначено:

,

,

,

.

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

.

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

3. Существенным резервом повышения точности и устойчивости счета на Эйлеровом этапе метода крупных частиц является конечно-разностная аппроксимация с вычислением промежуточных значений давления на границах крупных частиц (контрольных объемов). Объяснение этого факта имеется в [Патанкар] и состоит в следующем. Если поле давлений имеет распределение параметров с некоторой пульсацией давления между соседними частицами, то схема Эйлерова этапа (8.150) нисколько не отреагирует на наличие пульсаций, поскольку .

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

,

; (8.159)

В уравнениях (8.159) значения давлений и скоростей на границах крупных частиц (контрольных объемов) устанавливаются по результатам решения задачи о распаде слабого разрыва газодинамических параметров. Для границы с номером i (граница между крупными частицами с номерами I-1 и I) расчетные соотношения для и записываются в виде

,

, (8.160)

.

Вычислительные алгоритмы Эйлерового этапа (8.159), (8.160) используются при решении пространственных задач внутренней баллистики, в том числе, и на неортогональных сетках [ЧисЭксп,ВнуБал,ПроблемыТермо].

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

Отметим требования к моделям, вычислительным алгоритмам и к программному комплексу, реализованным при создании комплекса программ [АлиевКаталог]:

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

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

- элементарная вычислительная ячейка (крупная частица, контрольный объем) может иметь произвольное число граней, конкретное количество которых задается в исходных данных. Грани ячейки могут соединять соседние вычислительные ячейки, быть непроницаемой границей, быть перфорированной границей, быть поверхностью прогревающегося или горящего топлива. Следует обеспечить возможность вычисления потоков по граням ячеек с использованием всех алгоритмов, имеющихся в эксплуатируемом пакете программ (например, граница области с окружающей средой – критическое сечение соплового блока и т.п.);

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

При проведении внутрибаллистических расчетов в произвольной области (система воспламенения, объем, канал и т.д.) могут использоваться уравнения газодинамики в виде [АлиевАндреевВестник??]:

,

, i=1, N,

,

, (8.161)

,

.

В системе уравнений (8.161) последовательно записаны уравнения неразрывности для суммарной плотности продуктов сгорания, уравнения неразрывности для i – ой компоненты плотности (в продуктах сгорания всего N компонентов), уравнения сохранения количества движения (3 уравнения – для составляющих скорости ), уравнение сохранения энергии E (сумма внутренней и кинетической энергий). Здесь коэффициенты Ламе, имеющие разные значения в разных координатных системах. Коэффициенты , , , , , также зависят от выбранной системы координат. В частности, для цилиндрической системы координат справедливы соотношения:

;

.

Записанная выше система уравнений справедлива для произвольной координатной системы, элементарный объем в которой представлен на рис. 7.3 (на рисунке грани элементарного объема обозначены соответственно ).

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

Сформулируем модель элементарного объема, входящего в расчетную область внутри камеры двигателя:

- местоположение элементарной ячейки устанавливается в цилиндрической системе координат, привязанной к центральной оси двигателя (ось x);

- элементарная ячейка может иметь несколько граней (границ), часть из которых соединяет ячейку со смежными с ней ячейками (эти грани полностью открыты), часть граней (границ) – имеют частичную перфорацию, часть граней (границ) – непроницаемые твердые границы, на части граней (границах) – поступают продукты сгорания топлива;

- в элементарной ячейке могут размещаться источники массы и энергии, например, гранулированный воспламенительный состав (это, в частности, необходимо при моделировании элементарной ячейки как корпуса воспламенителя);

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

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

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

На первом (Эйлеровом) этапе дивергентные члены уравнений сохранения не учитываются. Вычисляются значения давлений на границах ячеек (рисунок 8.21, границы АВ, CD). На втором этапе (Лагранжевом) вычисляются потоки на границах смежных ячеек. На заключительном этапе вычисляется баланс масс, энергий и количества движения в каждой расчетной ячейке. Отметим, что при такой трактовке вычислительного метода логика алгоритмов при расчете пространственных течений остается той же, что и при расчете одномерных течений.

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

- на первом этапе для всех границ рассчитываются значения давлений на них;

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

- на третьем этапе для каждой расчетной ячейки рассчитываются потоки массы, количества движения и энергии, а также вычисляются внутренние источники массы, количества движения и энергии;

- на четвертом шаге выбирается шаг по времени и производится расчет масс, энергии и количества движения в ячейке;

- на последнем (заключительном) шаге вычисляются значения газодинамических величин в каждой ячейке.

В вычислительном алгоритме можно объединить первые два (или три) этапа метода.

Рассмотрим подробнее реализацию записанных выше этапов.

1. Значения давлений на границах расчетных ячеек (для ячейки с центром в О2 это границы АВ и CD) определяются, если известны значения газодинамических величин ( и др.) в ячейках с центрами в О1, О2, О3. Как и в методе С.К.Годунова [Год] будем предполагать, что на границах ячеек взаимодействуют два «полубесконечных» потока. Если время взаимодействия этих потоков ограничено (шаг вычислений по времени согласован с шагом по пространственной координате условием Куранта-Фридрихса-Леви), то погрешность вычисления давления на границе смежных ячеек будет допустимой. Как и в задаче о распаде произвольного разрыва параметров можно записать два элементарных решения для значения давления на границе. Для определенности при решении элементарной задачи о распаде разрыва принимается, что давление, например, по левую сторону от границы больше давления по правую сторону от границы. Кроме того, полагается, что разрыв параметров «слабый» [Рождеств]. Далее записывается зависимость давления на границе с учетом варианта развития процесса. В первом варианте в оба направления от границы распространяются ударные волны. Во втором варианте в одном направлении от границы распространяется ударная волна, а в противоположном направлении волна разрежения.

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

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

,

,

, (8.162)

.

Для записи уравнений в конечно-разностном виде систему уравнений полезно переписать в виде

,

,

, (8.163)

.

В таком виде при записи уравнения (8.163) для элементарного объема (при конечно-разностном решении уравнений) с ребрами вместо выражений , , можно использовать значения площадей граней элементарного объема, перпендикулярных соответственно осям , а вместо произведения использовать значение элементарного объема. Вычисление частных производных (или ) можно выполнить, используя значение градиента давления по нормали :

.

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

- для каждой грани вычисляются значения градиентов , , , (расчет производится по известным значениям давлений и скоростей в центре ячейки и на рассматриваемой грани, а также по известному значению расстояния от центра ячейки до рассматриваемой грани);

- по известным значениям , , , и направляющих косинусов , , , определяются значения производных , ;

- при выбранном шаге по времени из уравнений (8.163) находим промежуточные значения скоростей и энергии, просуммировав результаты по каждой грани;

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

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

3. На следующих этапах вычислительного алгоритма решаются уравнения

,