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

1. Обобщение рассмотренных в разд. 1 данной главы разностных схем типа (1.2) на случай одномерных систем уравнений гиперболического типа.

в случае ац = const, f = 0 производится с сохранением всех их свойств достаточно формально, если систему (1) умножением слева на матрицу П привести к характеристической форме v_f + Av_x = 0, v = Ли, или, в скалярной записи,.

затем для каждого из уравнений (2) выбрать ту или иную разностную схему вида (1.2) из числа рассмотренных выше (этот выбор может быть разным для различных характеристических направлений dx = Л/Л) и сделать обратный переход от v₍ к u = В уравнениях (1), (2) Л ;

диагональная матрица из собственных чисел матрицы А, определяемых характеристическим уравнением Det (А — Е) = О, Г2 — неособенная матрица, строками которой являются линейно независимые левые собственные векторы оЗ/ матрицы А, определяемые с точностью до их длины из совокупности линейных однородных систем уравнений сЗ/ (А — Х/?) = О, Л" ¹ — обратная в П матрица, и/ - инварианты Римана. В этом смысле рассматриваемые ниже разностные схемы можно отнести к сеточно-характеристическим.

Такие разностные схемы для системы (1) в случае ац = const, f = const имеют вид.

где Лд — диагональные матрицы, составленные из анализировавшихся в разд. 1 коэффициентов a?. В частности, аналогом разностных схем (1.2), (1.5), (1.10) с неотрицательными коэффициентами a? будут впервые предложенные К .О. Фридрихсом в работе [55] схемы с неотрицательно-определенными матрицами В д, для которых скалярные произведения.

при произвольных ненулевых векторах z. Как известно, в случае симметричных матриц А условия (4) и (1.10), т. е. оi^v_ui > 0 для матриц Вд из (3), эквивалентны. Однако и для систем уравнений (1) с несимметричными матрицами А, допускающими симметризацию, т. е. приведение (1) к виду A_xV1_t + A₂W_x ⁼ с симметричными матрицами A_iy А₂ и положительноопределенной матрицей А (t — гиперболичные по Фридрихсу, системы [55]), вместо условий (4) также вполне допустимо использовать (1.10). Примерами таких систем являются нестационарные уравнения газовой динамики и магнитогазодинамики, нестационарные уравнения теории упругости и др. (см., например, [20]). Заметим также, что в случае a_{j = const все разностные схемы (3) обладают свойством консервативности, как и исходная система (1), если положить F -Аи.

2. Для обобщения рассматриваемых разностных схем на случай квазилинейной системы (1) с А = A (t, х, u), f =f (/, х, и), в том числе для дивергентной формы таких уравнений.

наиболее универсальными и эффективными способами построения разностных схем типа (3) (в том числе обеспечивающих выполнение разностных аналогов законов сохранения, т. е. консервативных схем) являются методы типа Рунге-Кутты и интегроинтерполяционный (см., например, [68, 79]). Если с использованием этих методов на том или ином сеточном шаблоне построить соответствующее множество разностных схем с тем же, что ив (3), числом свободных параметров, а затем на примере системы с dfj = const перейти к записи (3) (полагая F = Ли), то можно установить связь этих свободных параметров с В? и тем самым в общем квазилинейном случае получить аналоги тех или иных разностных схем из числа рассмотренных в разд. 1 данной главы.

В качестве первого примера такого обобщения рассмотрим широко распространенный неявный шеститочечный шаблон (1.14) и воспользуемся интегроинтерполяционным методом [68].

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

при выборе подходящей интерполяционной формулы первого или второго порядков для вычисления F" *_t^l(²₂ ⁼ F (/^{n + I}/², х_{т ±} ^₂, УД) ^ап' проксимирует с первым или соответственно со вторым порядком точности исходную систему (5). Для ^" т¹^² также необходима некоторая связь со значениями этого вектора в узлах разностной сетки, например, простейшая = <�Рт ^в случае схем первого порядка точности.

Полагая.

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

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

Если потребовать аппроксимации с первым порядком точности разностным выражением (7) следствие из уравнений (5).

получим на неопределенные матричные коэффициенты D" дополнительно к (8) еще одно условие:

Используя это условие в соотношении (9), можно исключить еще какойлибо матричный коэффициент, например, Df = (Е —А)/2 + ADj — (E—A)D, и получить двухпараметрическое семейство консервативных разностных схем второго порядка точности (при надлежащем выборе ip^n+ll²) с произвольными матричными коэффициентами Щ, D}. Этот же результат можно получить, если в соотношениях (10) коэффициенты (а®)/ выразить через (aLi)/, (<*})/ с использованием для рассматриваемого здесь шаблона (1.14) условий второго порядка аппроксимации (1.7) при/: = 2. Выбирая те или иные значения коэффициентов (а^)*, (а!) / из числа рассмотренных в п. 3,5 из разд. 1, можно получить консервативные варианты этих схем для квазилинейных дивергентных систем уравнений (5). Например, при выборе (а})/, (а}) / в соответствии с (1.39), (1.40) при o_t >0 и аналогичными выражениями при a_f < 0 можно получить аналог разностной схемы, ''наиболее близкой" на сеточном шаблоне (1.14) к схеме с положительной аппроксимацией. Такая схема приведена в работе [63]. Наконец, если воспользоваться результатами п. 6 из разд. 1 и выбрать в соотношениях (9), (10) коэффициенты (а*)/, (а?)/, (а|)/ в соответствии с (1.46), получим гибридные варианты этих консервативных разностных схем, имеющие первый порядок точности при 7=1 (и монотонные, если.

а_А выбирать в соответствии с условиями (1.11) — (1.13)) и второй порядок точности при 7 = 0. В частности, если в качестве а_А в (1.46) выбрать схему Б. Карлсона (1.18), а в качестве а_в — схему (1.39), (1.40) при О/ > 0 и аналогичную при а/ < 0, т. е. в (10) положить:

то соотношения (9), (10), (13) при 7=1 переходят в консервативный вариант схемы Б. Карлсона для системы (5), имеющей минимальную аппроксимационную вязкость среда монотонных схем первого порядка точности на сеточном шаблоне (1.14), а при 7 = 0 — в консервативный аналог схемы (1.39), (1.40), которая среда всех разностных схем второго порядка точности для шаблона (1.14) дает наименее осциллирующее решение. В линейном скалярном случае этой схеме при 0<7<1 в пространстве коэффициентов соответствуют отрезки А_ХВ_Ь на рис. 4.2 и отрезки А₂2?₆ на рис. 4.3. Отметим также, что при выборе в качестве а_А разностной схемы (1.18), а в качестве а_в «левого прямоугольника» (1.28) при о_{ > 0 и аналогичной схемы — «правого прямоугольника» при <7/ < 0 получим обобщение на случай системы (5) разностной схемы из работы [93] (по сравнению с (13) она менее удачна при 7<^ 1 из-за плохого поведения на разрывных решениях разностных схем типа (1.28)). В линейном скалярном случае схеме [93] при о₍ >0 соответствуют отрезки А_хВ* на рис. 4.2 и отрезки А₂В_А на рис. 4.3.

В качестве примера на рис. 4.12 приведены результаты расчетов по схеме (9), (10), (13) задачи о распаде разрыва для одномерных нестационарных уравнений газовой динамики:

где р — плотность, р — давление, v — скорость, е — удельная внутренняя энергия, к = 1.4 — показатель адиабаты.

Область интегрирования 0 </<7 *1 выбиралась таким образом, чтобы за рассматриваемое время 0 </ возмущения от точки х = 0 не достигали ее границ. Приведены профили плотности в момент t =2,5, полученные в расчетах с числом Куранта (12) а, = 1, 2 и 4. Кривым 1 и 2 соответствуют расчеты при 7=1 и 7 = 0, точками отмечены расчеты для случая, когда значение параметра у связывалось с поведением численного решения так, чтобы вблизи разрывов у «1, а в области гладкого решения -7⁰. Штриховыми линиями обозначено точное решение. Первоначальный разрыв распадается на волну разрежения (слева), контактный разрыв (в центре) и ударную волну (справа). Нелинейная разностная система уравнений (9), (10), (13) решалась итерациями (с линеаризацией по Ньютону), а на каждой итерации использовалась матричная прогонка для решения соответствующей линейной системы. При числах Куранта а, < 1 до сходимости с е* ~ 10'⁶ требовалось 2−3 итерации, при а, > 1 необходимое число итераций увеличивалось, однако сходимость наблюда;

Рис. 4.12.

лась во всех рассчитывавшихся вариантах (проведено до о. =8). Видно, что при о,<1 в такой гибридной схеме удается полностью избавиться от осцилляций, как и в монотонной схеме Б. Карлсона (7=1, кривая /). При этом ширина зоны размазывания (наклон кривых на разрывах) остается такой же, как и в схеме второго порядка точности (7 = 0, кривая 2). При о, > 1 вначале за фронтом ударной волны, где число Куранта максимальное, а затем и вблизи волны разрежения в схеме второго порядка точности (и как следствие — в гибридной схеме, хотя и в меньшей степени) амплитуда осцилляций и глубина их проникновения увеличиваются, как и в линейном случае (рис. 4.11). Это не означает, конечно, что неявные схемы непригодны для расчетов разрывных решений. Например, довольно часто число Куранта в области интегрирования может изменяться очень сильно (на порядки) и большие его значения (ограничивающие временной шаг в явных схемах) достигаются в области гладкого решения, а у разрывов имеют существенно меньшие значения. Для таких и других аналогичных задач преимущества неявных схем бесспорны. Интересно отметить, что на контактном разрыве свойства всех трех вариантов этой схемы (7 = 0, 7 = 1 и гибридная) практически не зависят от значения числа Куранта.

3. Другим, столь же эффективным подходом к обобщению разностных схем, развитых для линейных уравнений и систем на нелинейный случай (используемым, как правило, в явных схемах), является метод РунгеКутты, на основе которого строятся различные многошаговые схемы с пересчетом (предиктор—корректор). Для схем третьего порядка точности такой подход использовался, например, в работах [79, 103] и ряде последующих работ ([83, 104−108] и др.). Сюда же относятся схемы второго порядка точности Лакса-Вендроффа [58], Маккормака [59] и др. Общим для всех таких схем является использование сеточного шаблона типа (1.19), для которого ниже проводятся обобщения различных разностных схем из п. 5, 6 разд. 1 на случай системы (5).

В качестве первого примера рассмотрим однопараметрическое семейство явных консервативных схем для (5) (для простоты в дальнейшем рассматривается случай $ = 0)

где В = Ь{ | — некоторая диагональная матрица, й вычисляется, например, по известной схеме второго порядка точности Лакса-Вендроффа [58] или аналогичной схеме Маккормака [59], в линейном случае соответствующим точке Вs на рис. 4.2, 4.4, й = и^л или и в зависимости от того, по данным какого временного слоя производится сглаживание с помощью оператора в квадратных скобках из правой части (15). В вычислительной практике широко используются соотношения типа (15), однако обычно со скалярным коэффициентом сглаживания (Ь, — = Ь), подбираемым из различных (часто чисто эвристических) соображений. Использование для всех характеристических направлений системы (5) dx = _fdt коэффициента сглаживания b разных его значений b_t для каждого из условий совместимости, к которым может быть приведена система (5), позволяет заметно расширить возможности для оптимального в том или ином смысле выбора значений этих коэффициентов. Эти оптимальные значения зависят от Of = /т/А и являются разными для разных характеристических направлений.

При сглаживании по данным на слое /=/^я, т. е. при й=и^пв (15), и выборе.

в линейном случае для каждого из условий совместности типа (1.1), к которым приводится (5), увеличение у ^в (16), т. е. увеличение А,-, приводит к тому, что точка на рис. 4.4, соответствующая этой схеме, движется от Вs (при 7 = 0) вниз по вертикальной оси (так как (а°₂)/ ⁼ = (сг5)/ =0) и при 7= 1 совпадает с точкой А_х. Дальнейшее увеличение у приводит лишь к излишнему сглаживанию численного решения. В частности, при bf = (1 — о})/2 получим разностную схему, аналогичную схеме П. Лакса [53] (точка О на рис. 4.4), а при еше больших значениях у в пространстве коэффициентов разностной схеме (15) будут соответствовать точки, расположенные ниже многогранника схем с положительной аппроксимацией (неустойчивые схемы). Отсюда следует ограничение сверху на величину b при использовании скалярного коэффициента сглаживания в (15) (при й = и^п):

При использовании в (15) сглаживания по данным временного слоя t=tⁿ*^l₉ т. е. й=й, линейному скалярному случаю (1.1), >0 в пространстве коэффициентов соответствует луч с началом в точке В₅ (при.

b_t = b = 0), направленный в сторону точки G₂ (при b > 0), ни при каких значениях не пересекающий многогранника разностных схем с положительной аппроксимацией О, А _х А₃А₄А_$ А ₆, хотя и проходящий довольно близко от него (рис. 4.4). Ближайшей к этому многограннику точкой луча B₅G₂ является точка С|, для которой 6 = а/(2 + 5о). Тогда при выборе в (15).

такой разностной схеме будут соответствовать точки отрезка B_SGj на рис. 4.4.

Другими распространенными примерами явных схем с пересчетом, использующих сглаживание на верхнем временном слое, являются некоторые варианты метода коррекции потоков (FCT), предложенные для одного скалярного уравнения (1.1) с переменными, в том числе знакопеременными значениями X [95]. Основная идея этих схем заключается в предварительном вычислении v_m по схеме с положительной аппроксимацией, например,.

и в последующем восстановлении в области гладкою решения второго порядка точности с помощью пересчета.

В пространстве коэффициентов разностных схем на шаблоне (1.19) (рис. 4.4) первому этапу соответствует некоторая точка на отрезке А (ее положение зависит от выбора i>), а второму — переход из нее в некоторую точку на плоскости Bj, j = 1, 2 … схем второго порядка точности. Направление этого перехода и положение точки перехода на плоскости Bj, / = 1, 2, … однозначно определяется видом «анти сглаживающего» оператора при расчете и зависит от величины а= Xr/Л. Такой подход, очевидно, предпочтительнее (15), (17), поскольку в нем заведомо осуществляется (там, где это необходимо) выход на схему с положительной аппроксимацией. Кроме того, схема (15) с вычислением и по двухшаговым схемам типа Лакса-Вендроффа, но существу, является трехшаговой схемой, а с помощью двух пересчетов можно было бы уже построить схему с третьим порядком точности. В качестве монотонной схемы первого порядка точности в явных схемах коррекции потоков предпочтительнее выбирать разностную схему, соответствующую ближайшей к плоскости Bj, / =1,2,… точке А, тогда и заключительному этапу на этой плоскости будет соответствовать точка (разностная схема), не слишком уклоняющаяся от ближайшей к многограннику ОАА₃А₄А_ъА_ь точки В₆. Для случая системы (5) такой разностной схемой будет (15) при и = и, вычислении и по схеме (11) и выборе.

Тогда при у = 1 имеем монотонную схему (11), при 7 = 0 (11), (15), (18) будет иметь второй порядок точности на решениях (5), а при промежуточных значениях у в линейном случае 0 < o_t <1 этой схеме будет соответствовать отрезок АС на рис. 4.4 (точка С! попадает на прямую разностных схем третьего порядка точности только при o_t — 0,5, как это изображено на рис. 4.4).

Рассмотрим еше один пример гибридной явной двухшаговой схемы для (5), соответствующей в линейном скалярном случае однопараметрическому семейству разностных схем, расположенных на отрезке А_ХВ_Ь% концы которого, как показано в п. 4,5 из разд. 1, соответствуют наиболее точной монотонной схеме первого порядка и наименее осциллирующей схеме второго. порядка точности (точки А у и В₆ на рис. 4.4). Такая разностная схема может быть реализована по-разному, например.

Как видно, в данной реализации в качестве предиктора используется неустойчивая схема для вычисления и_т, однако в целом схема (19), (20) устойчива при числах Куранта < 2 и любых значениях у из диапазона [О, 1]. При у = 1 она переходит в некоторую монотонную схему (в частности, в (11) при о. < 1), а при 7=0 имеет второй порядок точности на решениях (5) и аналогична явной, «наименее осциллирующей» схеме второго порядка точности из (63] (в линейном случае (1.37), (1.38)). Этот пример еше раз подчеркивает известный из вычислительной практики факт, что о свойствах многошаговых разностных схем (а также схем с искусственной вязкостью, сглаживанием и т. п.) необходимо судить в целом. а не по отдельным их этапам, т. е. хороший элемент в какой-либо схеме может приводить к совсем нежелательным последствиям в другой схеме (и наоборот, естественно). Это связано с тем обстоятельством, что каждая комбинация значений сеточной функции в узлах сеточного шаблона со свободным параметром при ней, добавляемая к некоторой опорной схеме (точке в пространстве коэффициентов), определяет в этом пространстве некоторую прямую, проходящую через эту точку, — опорную схему. Изменяя опорную схему, тем самым переносим в пространстве коэффициентов эту прямую параллельно самой себе и теперь уже различным значениям свободного параметра могут соответствовать совсем другие по своим свойствам разностные схемы, отличные как от опорной схемы, так и от первоначального семейства (например, эта прямая теперь может почти всюду лежать в области неустойчивых схем, не пересекать множества схем с положительной аппроксимацией и т. п.).

Приведем еше один пример гибридной (трехшаговой) разностной схемы для (5):

которая при 7 = 1 переходит в монотонную схему (11). При у = 0 со скалярным #, а не матричным n^_lGH, как в (24) коэффициентом в выражении для 6, т. е. при.

эта схема является предложенным в работе [79] однопараметрическим семейством разностных схем третьего порядка точности на решениях.

(5), устойчивым (как показано в этой же работе, см. также (1.23)) при о* =т шах |Х/ |/Л <1 и — 1/8.

(aj — 4)/24.

i, т

На рис. 4.13 показана область допускаемых этим условием устойчивости значений параметра g (область между штриховыми линиями D и С₃). Так же штриховыми линиями D и С показаны рекомендуемые в работе [83] на основе анализа диссипативных и дисперсионных свойств разностной схемы типа (21) — (23), (25) (в которой (21), (22) заменены более удобной для обобщения на многомерный случай схемой Маккормака [59]) значения скалярного параметра g в (25): схема с минимальной диссипацией, для которой g = а^_вх (oj^ - 4)/24, точка D на рис. 4.4 и схема с минимальной дисперсией, для которой g = (4 о max + 1) Х X (ff?mx — 4)/120, точка С на рис. 4.4, o_max = г max|(Xf)? 1/Л. Сопоставляя на сеточном шаблоне (1.19) соотношения (21) — (23), (25) в линейном скалярном случае с (1.2), (1.5), (1.22), можно установить связь между свободными параметрами g и а°_₂: g = а® ₂ — а (а + 2) (2а — 1)/24, откуда для варианта схемы (21) — (23), (24), у = 0, ''наиболее близкого" к схемам с положительной аппроксимацией, из (1.45) при, а >0 и аналогичного выражения при о < 0 получены оптимальные значения этого параметра:



Зависимость от |а/1 параметра g_t из выражения (26) показана на рис. 4.13 сплошной линией (кривая С₂). При промежуточных значениях 0 < у < 1 схеме (21) — (24) в линейном скалярном случае и о > 0 соответствует отрезок A i С₂ на рис. 4.4.

На рис. 4.14 в сечении х = 2,77 (через 120 шагов) в зависимости от т = y/h представлены профили плотности р из расчетов по различным схемам автомодельной стационарной задачи о взаимодействии прямолинейной ударной волны АС с контактным разрывом ВС (рис. 4.15), в которой после взаимодействия в точке С образуются преломленные ударная волна l контактный разрыв 2 и веер волн разрежения 3−4. Для данной задачи в (5) имеем.



Роль гиперболической переменной t играет пространственная координата Ху а за х в (5) принята другая пространственная координата у. Начальные значения компонент и, и вектора скорости по осям х иу> внутренней энергии е и давления р выбирались следующими:



Данные на рис. 4.14, а-д соответствуют расчетам при о, = 1, а на рис. 4.14, е — при о* = 2. Штриховыми линиями показано точное решение.

Кривые 1 на рис. 4.14, a-в, д соответствуют расчетам по монотонной схеме первого порядка точности с минимальной аппроксимационной вязкостью (11), кривая 0 на рис. 4.14, а — схема П. Лакса (неоптимальная в смысле величины аппроксимационной вязкости разностная схема — точка 0 на рис. 4.2−4.5). Как и в линейном скалярном случае (рис. 4.7), преимущества схем с положительной аппроксимацией ''наиболее близких" в пространстве коэффициентов к схемам более высокого порядка точности по сравнению с неоптимальными в этом отношении схемами здесь также явно заметны.

Если схема (11) дает вполне удовлетворительное решение и на разрывах типа ударной волны вполне конкурентоспособна со схемами более высокого порядка точности, то схема П. Лакса и другие неоптимальные схемы с положительной аппроксимацией практически непригодны для такого рода задач.

О влиянии близости в пространстве коэффициентов схем второго и более высокого порядка точности к множеству схем с положительной аппроксимацией можно судить из сравнения результатов расчетов по оптимальным в этом отношении схемам второго порядка точности (19), (20), 7 = 0 (кривая 2 на рис. 4.14, а) и третьего порядка точности (21) — (24), 7 = 0 (кривая 3 на рис. 4.14, в) и по неоптимальным схемам второго порядка точности Лакса-Вен дроф фа [58] и аналогичной схеме Маккормака [59] (кривые 2 на рис. 4.14, г, б). Видно, что использование таких оптимальных схем позволяет заметно уменьшить амплитуду осцилляций и глубину их распространения, хотя полностью от них избавиться не удается, как это следует из результатов работы [63] и п. 3 разд. 1.

Результаты расчетов задачи (27) по гибридным сеточно-характеристическим схемам 0 < 7 < 1 ((15), (16) — рис. 4.14,6; (15), (17) — рис. 4.14, г; (И), (15), (18) — рис. 4.14, д; (19), (20)-рис.4.14А1 на рис. 4.14, е — данные расчетов при а. = 2 по монотонной схеме первого порядка точности (19), (20), 7=1. Видно, что для всех гибридных сеточно-характеристических схем удается практически полностью избавиться от осцилляций, при этом ширина зоны размазывания разрывов в этих схемах остается такой же, как и в опорных схемах второго и третьего порядка точности. Несколько лучшие результаты дает схема (19), (20) (рис. 4.14, а), что вполне естественно, поскольку и без регуляризации (кривая 2 на рис. 4.14, а) она дает малоосциллирующее решение. Отметим, яго ее положительные свойства сохраняются вплоть до а, = 2 (рис. 4.14, ё). Кривая 1 на рис. 4.14, г — данные расчетов по схеме (15), (17) при выборе скалярного коэффициента сглаживания b_{ = b *= 0, 3. Видно, что подобрать хороший единый для всех условий совместности коэффициент сглаживания довольно сложно. Это же замечание относится и к схемам (21) — (23), (25).

В рассмотренных выше гибридных схемах был оставлен открытым вопрос о выборе конкретных значений скалярного параметра 0 <7 < 1, поскольку его выбор, обычно связываемый с поведением решения, является независимой задачей. Общие требования к его выбору достаточно очевидны — он должен принимать значения 0 в области гладкого решения и значения 1 вблизи разрывов, для чего прежде всего требуются эффективные алгоритмы локализации положения разрывов. Конкретная реализация выбора 7 может быть разной и от удачного подбора алгоритма для его определения в известной степени зависит окончательный результат. Некоторые способы выбора у (или его аналогов) приводятся, например, в работах [92−100, 109]. Следует подчеркнуть, что при построении всех рассматриваемых здесь гибридных схем существенно использовались характеристические свойства систем уравнений гиперболического типа и единый для всех уравнений параметр обеспечивает при у = 1 выход на ту или иную опорную монотонную разностную схему с минимальной аппроксимационной вязкостью, что отличает рассматриваемый здесь подход от других известных способов регуляризации разрывных численных решений. В представленных на рис. 4.12, 4.14 расчетах алгоритм выбора у был аналогичен [94], например,

и некоторые другие его разновидности. Типичное поведение у в численных расчетах показано на рис. 4.14 штрих-пунктирной линией.

4. Рассмотренные в п. 2, 3 разностные схемы строились, исходя из дивергентных уравнений (5). Они обладают свойством консервативности, как известно [110] играющим важную роль при сквозном расчете разрывных решений. Получение не дивергентных вариантов всех этих схем не представляет труда, например, заменой разностей типа F_m4.j — F_m на.

^mf|/2(«m + l «Um) И Т.Д.

На рис. 4.16 для задачи (14) с точным решение (штриховые линии) сравниваются результаты расчетов (профили плотности р) по дивергентным (рис. 4.16, j) и недивергентным (рис. 4.16,6) вариантам схемы (19), (20). Кривые / соответствуют расчетам по схеме первого порядка точности (у = 1), кривые 2 — варианту схемы второго порядка точности (у =* 0), точки — гибридной схеме (у из (28)). Видно, что при малой интенсивности ударной волны (по отношению давлений до полутора) недивергентные схемы вполне удовлетворительно воспроизводят такое разрывное решение, в частности, по параметрам за ударной волной, хотя скорость движения ударной волны в них воспроизведена несколько хуже, чем в дивергентных вариантах этих схем (т.е. в недивергентных схемах соотношения на разрывах выполняются с заметной ошибкой). Однако при увеличении интенсивности разрывов нсдивергентные схемы становятся практически непригодными к сквозному расчету разрывных решений и требуют при их использовании явного выделения разрывов большой интенсивности. Как видно из представленных на рис. 4.17 профилей скорости и, полученных по дивергентному (кривая 1) и недивергентому (кривая 2) вариантам гибридной схемы (19), (20), (28) для той же задачи (14), но с более интенсивными разрывами (р (0, х) = р (0, х) = 10 при х <0), в этом случае недивергентная схема уже весьма неудовлетворительно воспроизводит разрывы, в том числе и контактный, где скорость и должна быть непрерывной. Точное решение на рис. 4.17 показано штриховой линией.

С другой стороны, йспользование консервативных разностных схем для строго дивергентных уравнений (5) (с нулевой правой частью у = 0) в некоторых случаях также может приводить к непрятным последствиям. Например, если хотя бы в одной из систем координат уравнения могут быть записаны в дивергентной форме (5) с = 0, то при произвольном взаимно-однозначном преобразовании независимых переменных f = f (r, т?), х — х (т, г?), как известно (см. разд. 2 гл. 1), заменой w = u/т?*, Ф =.



Рис. 4.16.



Рис. 4.17.

= (r?, u + rj_x F)/rj_x система (5) приводится к эквивалентной (5) дивергентной в новых переменных г, rj форме (см. разд. 2 гл. 1):

w_T + = 0. (29).

Рассмотрим простейшую задачу Коши для (5), у = 0, а именно и (0,х) = = u° = 1 и®,…, Uj |, и* = const, 1 = 1,…,/. Тогда и для всех t > О и (Г, х) = u°. Однако если сделать некоторую нелинейную замену независимых переменных, например, t = т, х = х (т, т?), Vx ^ 0,°°, то w (0, т?) и Ф (0,т?) становятся функциями rj и любая консервативная схема для (29) из-за ошибок аппроксимации не будет обеспечивать выполнение условия u (f, х) = и⁰, в то время как недивергентные схемы свободны от этого недостатка. Этот эффект особенно заметен вблизи оси цилиндрической, центра сферической и т. д. систем координат, где происходит наиболее резкое изменение преобразованной фукнции Ф и имеют место наибольшие значения ошибок аппроксимации, пропорциональных тем или иным пространственным производным от Ф. Обычно в таких ситуациях используют различные промежуточные между (1) и (5) формы уравнений (например, [18, 20] и др.). Однако возможны и другие подходы. Для дивергентных уравнений (5), помимо задачи отыскания обобщенного решения (аналогом которой является использование консервативных схем сквозного счета), распространен подход, связанный с разбиением области интегрирования на подобласти гладкого решения, разделенные выделяемыми явно поверхностями разрыва. Дивергентная форма уравнений (5) в таком подходе используется только на поверхностях разрывов. Этот подход широко и успешно используется на основе не дивергентных численных методов (особенно для схем первого порядка точности) и позволяет на сравнительно грубых сетках получить достаточно точные результаты, когда число разрывов большой интенсивности сравнительно невелико, поскольку слабые разрывы удовлетворительно воспроизводятся недивергентными схемами. Когда такой подход становится слишком громоздким (из-за большого числа интенсивных разрывов, сложной их формы и т. п.), можно воспользоваться также схемами сквозного счета, однако консервативные варианты схем в этом случае применять только там, где они действительно необходимы, т. е. вблизи разрывов и других поверхностей негладкого решения. В области гладкого решения — использовать недивергентные варианты этих схем, т. е. поступать так, как это делается в гибридных схемах с порядком аппроксимации. Такие локально-дивергентные схемы выгодны еще и тем, что расчет по консервативным вариантам (обычно требующий заметно большего, чем недивергентные варианты этих же схем, числа арифметических операций) проводится лишь в относительно небольшом числе узлов от их общего количества в области интегрирования. Пример расчета задачи о распаде разрыва с р (0, х) = р (0, х) = 10 при х <0 не основе гибридной схемы из работы [64] показан на рис. 4.17 точками. Вблизи разрыва применялась непосредственно оптимальная явная схема из [64], а в области гладкого решения использовался недивергентный вариант этой же схемы (таких точек на слое t = const было в среднем 87—90%). Как видно, такая локально-дивергентная схема по своим свойствам практически не уступает чисто консервативным схемам (19), (20).