Явный безусловно устойчивый разностный метод расчета течений жидкости Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 1, № 1, 1996

ЯВНЫЙ БЕЗУСЛОВНО УСТОЙЧИВЫЙ РАЗНОСТНЫЙ МЕТОД РАСЧЕТА ТЕЧЕНИЙ ЖИДКОСТИ*

В. Ф. КУРОПАТЕНКО Российский федеральный ядерный центр — Всесоюзный научно-исследовательский институт технической физики

Снежинск, Россия

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

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

Пусть движение жидкости описывается системой законов сохранения в виде

— — = 0 (1)

дЬ дт '

дЬ + дт ' ( )

1(Е+2 и2)+дт(ри '=0' (3)

где V — удельный объем, и — массовая скорость, Р — давление, Е — удельная внутренняя энергия, Ь — время, т — массовая лагранжева координата.

Дифференциальные уравнения (1)-(3) аппроксимируем разностными уравнениями. В литературе известно много разностных схем, применяемых для численного интегрирования системы (1)-(3) и обладающих различными достоинствами и недостатками. Ограничимся обсуждением только двуслойных разностных схем (РС), связывающих решение в точках сетки для двух соседних моментов времени Ьп и Ьп+1 = Ьп + т. Все такие РС делятся на явные и неявные. Разностную схему

ип+1 = сп+1ип,

*© В. Ф. Куропатенко, 1996.

где оператор шага Cn+1(r,h,T) = вTe, T — оператор сдвига, в — матрица в 4-мерном пространстве компонент вектор-функции

U(m,t) = {V(m,t),U(m, t), P(m,t),E(m,t)},

в = -qi, —qi + 1, ■ ■■, q2 — 1,92, qi > 0, q2 > 0, qi + > 1,

будем называть явной, если для любого наперед заданного Q из промежутка [1, N), где N — число точек сетки в области интегрирования, найдется такое т > 0, что q(r) = q1 + q2 < Q. В противном случае разностную схему будем называть неявной.

Как правило, явные РС условно устойчивы [1-3]. Условие устойчивости формулируется в виде ограничений на соотношение шагов сетки по времени и по лагранжевой координате. В случае, когда различные фрагменты рассчитываемой системы имеют разные размеры, шаг по времени выбирается по малому элементу системы и оказывается неоправданно малым для больших элементов той же системы. Таким образом, применение явных условно устойчивых РС для расчета течений с разномасштабными элементами оказывается экономически невыгодным. В качестве примера рассмотрим, как влияет тонкая прослойка на увеличение времени решения на ЭВМ задачи о распространении ударной волны в одномерной системе размером H. Ударная волна, распространяющаяся со скоростью W, проходит всю область H за время T. Требуется определить состояние системы в момент времени 5T. Будем считать ударную волну слабой и положим W = const. В однородных явных условно устойчивых РС ударная волна проходит одну точку сетки за ~ 5 шагов по времени. Если в области интегрирования взять N точек сетки по пространству и на расчет всех величин в одной точке потратить ^ операций, то для расчета задачи потребуется 25^N2 операций. Изменим систему, поместив внутри области H тонкую прослойку толщиной h. Возьмем в прослойке n точек по пространству. Значения h и n та-Hh

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

Nn

этом случае на решение задачи потребуется затратить + и) ^ операций. Таким

п(М + и)Н ^ ^^

образом, количество операций возрастает в <р &-- раз. Если применяемая ЭВМ

такова, что на решение задачи без прослойки при N = 100 требуется ~ 20 мин, то для решения задачи с прослойкой при и = 10, к = 10-3Н потребуется более 30 часов. Если же решение в прослойке требуется найти с более высокой точностью, то и = 10 недостаточно. Но уже при и = 50 время решения задачи с прослойкой возрастет до ~ 250 часов.

Безусловно, при разных — и — количественные результаты будут разными, но каче-

N Н

ственный вывод сохранится: для решения задач с тонкими прослойками требуется затратить заметно больше ресурсов ЭВМ.

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

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

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

В данной работе построена такая явная разностная схема для системы законов сохранения (1)-(3), в которой условие устойчивости трансформируется в условие для выбора числа точек в момент £п, по которым находится решение в заданной точке в момент £га+1. Иными словами, сеточный шаблон, на котором определяются разностные уравнения, зависит от шага т. В целом разностная схема с переменным шаблоном оказывается безусловно устойчивой. Идея построения такой разностной схемы впервые изложена в [4], где построена явная разностная схема, в которой ограничение на соотношение шагов т и к заменяется условием для определения в момент £п основания характеристического треугольника, по которому проводится суммирование значений давлений и скоростей для определения вспомогательных значений этих величин. Область применимости разностной схемы [4] ограничена только гладкими течениями без ударных волн.

Рассмотрим комбинированную разностную схему, в которой для расчета ударных волн используется метод из [3], а адиабатические гидродинамические течения рассчитываются с помощью явной безусловно устойчивой разностной схемы, идейно близкой к разностной схеме из [4].

Аппроксимируем уравнения (1), (2) разностными уравнениями

т/га+1 — Vй ТТ * — ТТ *

Vi+0.5 4+0.5 Тг+1 иг _ д (4)

тк

ТТП+1 ттп Р* _ Р*

Т_+ рг+0.5 рг—0.5 _ д (5)

тк

Рассмотрим два типа решений: адиабатические волны разрежения (Д-волна) и неадиабатические волны сжатия (¿"-волна). В случае ¿"-волны (ди/дх < 0) вспомогательные значения Р * и и * определим следующим образом:

и* _ ТП+1; Р+0.5 _ Рп+0.5- (6)

Далее найдем вспомогательное значение Р¿+05. Предположим, что в сеточном интервале Xi < х < Хг+1 распространяется ударная волна, на фронте которой справедливы уравнения

— V—)W + (Т+ — Т—) _ 0, (7)

(Т+ — ТЦ—^ — (Р+ — Р—) _ 0, (8)

Е+ — Е— + 0.5(Р+ + Р— — V—)_0, (9)

где величины с индексом "-"характеризуют состояние перед фронтом разрыва, а с индексом "+—за фронтом. Будем считать, что

Р _ Рп • Р _ Рп • Е _ Е п Р— _ ^¿+0.5; И— _ ^¿+0.5; Е— _ Ei+0.5,

ди _ и+ — и— _ ип++11 — ип+\

Поскольку Рп+0 5, Рп+05, Еп+05, и;++11, и™+1 известны, то уравнения (7)-(9) вместе с уравнением состояния

Р _ Р(р, Е) (10)

образуют систему четырех уравнений с четырьмя неизвестными Р+, р+, Е+, Ш. Решив эту систему, найдем Р+ и положим

ркб = (11)

Внутренняя энергия Е*++15 определяется из разностного уравнения

Ег+0.5 = Ег+0.5 — -¿+-0.5 + Р г+0.5)( ^¿+0.5 — ^¿+0.5); (12)

совпадающего с уравнением ударной адиабаты. Таким образом, обеспечивается правильное изменение энтропии при переходе от En к En+1 в "размазанной ударной волне". В работе [2] показано, что рассмотренная РС для S-волны условно устойчива и шаги т и h должны удовлетворять условию

Ta < <!3)

где а — массовая скорость звука, определяемая уравнением

В случае R-волны вспомогательные значения Р*, U* будем вычислять по формулам

k2

р+0.5 = 0-5(в0Рг+о.5 + Pi+0.5+v ), (15)

v=-ki

k2

U* = 0.5(eoUr+1 + Y, &Un+v1), (16)

v=-kl

где

ev = hi+0.5+v < 1.

Ta+0.5+v

Значение v в суммах (15), (16) увеличивается до тех пор, пока выполняются условия

kl - 1 k2 - 1

Е в-v < 1; Ё ev < 1. (17)

v=0 v=0

При нарушении условий (17)

k1-i k2-i

^в-v + в-ki > 1, Ë^v + ek2 > 1

v=0 v=0

значения e-k1 и ^k2 пересчитываются по формулам

kl -1 k2 - 1 в-kl = 1 -Y, в-v , ek2 = 1 ^ ev.

v=0 v=0

При решении задачи о распространении ударной волны по однородному веществу шаг по времени т будет выбираться в сеточных интервалах "¿"с S-волной по условию

та=1. (18)

При этом в интервалах с Д-волной может оказаться

Но тогда из (15)—(17) следует, что

р*

р г+0.5

в = ^ > 1. та

туп тт* _ т тп+1

р г+0.5, Ui — Ui ■

(19)

(20)

РС с такими вспомогательными значениями Р * и и * исследована в [3].

Рассмотрим погрешность аппроксимации уравнений (1), (2) разностными уравнениями (4), (5) со вспомогательными величинами (20). Условия точности получим, ограничив главные члены этих выражений

1 д

2 дг2 1 д 2и

т —

-т +

1 д3и

24 дт3 1 д3 Р

К

К2

2 дг2 24 дт3 С помощью (1), (2) преобразуем (21) к виду

' 1 ди 1 д2 и

< £г,

< £2.

(21)

_д_

дт

дт

2 дг Т 24 дт2

К2

1 дР 1 д2Р, 2 ---т +---К2

2 дг +24 дт2

< £1,

< £2.

Умножим каждое из этих неравенств на дт, проинтегрируем по т от 0 до К и запишем полученный результат в виде

1 ди 2 дгт — 1 д2 иК2 24 дт2 < £1К,

1 дР 2 ~дГт — 1 д2 РК2 24 дт2 < £2К,

(22)

ди дР дР

дг дт дг

условия точности примут вид

ди

Заменим на — —— и на —а ——. После этой замены и несложных преобразований

дт

т1 <

2£1К + -2 К2 д 2и дт2

дР дт

т2 <

2£2К +^ к2 д 2Р дт2

2 ди а2 — дт

(23)

Заменим производные в (23) разностями, получим условия точности в виде

т1 < К

2£1К + — |ит — 2иг + Ui-l|

1^+0.5 — Р^о5

(24)

2^2^ + 12 |Рг+1.5 — 2Рг+0.5 + Д-О.б! а2+0.5 |иг+1 - иг1

В случае движений без ускорения и деформации Рг+0.5 = Рг-0.5, иг+1 = и и из (24), (25) следует, что т1 < то, т2 < то, т.е. условия точности не дают ограничений на т.

В адиабатическом течении внутренняя энергия и давление вдоль траектории частицы зависят только от V. Изменения Р и Е при изменении V определяется системой уравнений

1 + Р1 = 0- Р = Р>• <26>

где V = —. Первое из этих уравнений является следствием трех законов сохранения (1)-Р

(3). Запишем его в виде

уп + 1

Еп+1 = Еп - ! Р(V, Е)dV. (27)

V п

Заменим интеграл интегральной суммой

г

Еп+1 = Еп - ^ РкV, Ек)(У^+1 - ^), (28)

к=0

где число слагаемых г выбирается так, чтобы обеспечить нужную точность интегрирования вдоль изэнтропы. Значения 'Ук+1 , Ек+1 определяются из уравнений

к + 1

Vk+l = Vn + - Vn), (29)

Х + 1

Ек+1 = Ек - —^Рк(Vк,Ек)(Vn+1 - Vn), (30)

х + —

а давление Рк+1 — из уравнения состояния

Рк+1 = Р (^к+1,Ек+1)-

При использовании изложенной разностной схемы для расчета характеристик движения жидкости с ударными волнами шаг по времени т выбирается лишь в интервалах сетки, где решение является ¿"-волной:

т = ш1п(тг(5)}, (3—)

здесь тДб") = —. При этом в некоторых интервалах сетки с Д-волной может оказаться, что а

т^г

-—— > —. В этих интервалах вспомогательные значения и* и Рг+0 5 выбираются с помощью уравнений (15), (16).

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

1 ЭТАП. Ударная волна распространяется по веществу А и не дошла до прослойки Б. В этом случае шаг по времени выбирается из условий (31) в тех интервалах сетки,

где течение есть б'-волна. В остальных интервалах сетки, где решение является Л-волной, ограничение на шаг т снимается применением изложенной выше схемы.

2 ЭТАП. Ударная волна проходит через вещество Б. Пусть для определенности шаг ЛБ пространственной сетки в этом тонком слое в к раз меньше Ла. Поскольку в интервалах сетки в веществе Б решение есть б'-волна, то шаг по времени будет выбираться из условия

т < Лв = Ла

ав кав

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

3 ЭТАП. Ударная волна ушла от прослойки Б. В этом случае шаг по времени выбирается, как на этапе 1, то есть снова возрастает в к раз.

Рассмотрим две задачи с тонкими прослойками.

ЗАДАЧА 1. Жесткость вещества Б (прослойка) выше, чем жесткость вещества А (рБсБ > РаСа).

Расчеты проводились для системы 0 < х < 100, причем в областях 0 < х < 50; 50.1 < х < 100 находится вещество А с параметрами Р0 = 0, р0 = 1, Е0 = 0, и0 = 0, а в области 50 < х < 50.1 находится вещество Б с начальными параметрами Р0 = 0, р0 = 10, Е0 = 0, и0 = 0. Уравнения состояния были взяты в виде

А: Р = (7 - 1)рЕ, 7 = 5/3;

Б: Р = (7 - 1)рЕ + Со2к(р - рок), рок = 10, Сок = 5, 7 = 3. В качестве краевого условия при х = 0 бралось и = 0.5.

\ ^^ 6 7^—--—

\ '— -\ <С5 5^^^^ у

^ 4 \ 4

ч N 3 N 3

N. 2 2 /

1 А Б А

т

Рис. 1.

Таким образом, в точке £ = 0, х = 0 имеется сильный разрыв — ударная волна с параметрами Р1 = 0.333, и = 0.5. В момент прихода этой волны на прослойку образуются отраженная и прошедшая ударные волны. После выхода ударной волны на вторую границу прослойки происходит распад разрыва с образованием прямой ударной волны и отраженной волны разрежения. Амплитуда этой ударной волны меньше, чем амплитуда прошедшей ударной волны. Далее начинается взаимодействие волны разрежения с контактной

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

и

Рис. 2.

Результаты расчетов приведены на рис. 1, 2 и в табл. 1. На рис. 1 схематически изображена волновая картина в переменных т, а на рис. 2 — состояния в прослойке и на прошедшей и отраженной волнах в переменных Р, и. Пунктиром (см. рис. 1) обозначены траектории волн разрежения.

Таблица 1

№пп. X Р и х - 50.1 * = 0.1 / Р - Р ^ = Р1

1 50.00 0.333 0.500 -1 0

2 50.10 1.807 0.036 0 4.426

3 50.41 0.249 0.433 3.1 -0.25

4 50.88 0.316 0.487 7.8 -0.05

5 51.40 0.332 0.500 13 -0.003

6 51.93 0.336 0.502 18.3 0.009

7 52.46 0.336 0.502 23.6 0.009

8 52.98 0.335 0.501 28.8 0.006

9 55.09 0.335 0.501 49.9 0.006

10 57.20 0.335 0.501 71 0.006

11 59.31 0.334 0.501 92.1 0.003

12 67.24 0.333 0.500 171.4 0

В таблице приведены зависимости давления и массовой скорости на фронте прошедшей ударной волны в зависимости от координаты фронта волны. Расчеты показывают, что на расстоянии порядка 15 толщин прослойки давление на прошедшей волне меньше, чем на падающей. Максимальное отклонение в этой области достигает 25 %. Затем давление

становится больше. Область избыточного давления находится расстоянии примерно от 15 до 100 толщин прослойки. При ^ > —00 амплитуда прошедшей ударной волны становится равна амплитуде падающей. Максимальное превышение давления на фронте прошедшей волны давлением на фронте падающей волны для данной пары веществ составило 1 %.

Задача решалась по РС 1 из работы [2], по изложенной выше РС и по схеме из [5] с выделением разрывов до выхода прошедшей ударной волны на правую границу системы с координатой х = —00. В изложенной выше схеме затраты времени ЭВМ для решения задачи оказались в 17 раз меньше.

ЗАДАЧА 2. Жесткость вещества Б (прослойки) меньше, чем жесткость вещества А (РбСБ < РАСА).

Расчеты проводились для системы 0 < х < —00, однако в областях 0 < х < 50, 50.1 < х < —00 находилось вещество А с начальными параметрами Р0 = 0, р0 = —0, Е0 = 0, и0 = 0. В области 50 < х < 50.— находилось вещество Б с начальными параметрами Р0 = 0, р0 = —, Е0 = 0, и0 = 0. Для вещества А было взято уравнение состояния

Р = (7 - —)рЕ + С0к(р - Р0к), где р0к = —0, С0к = 5, 7 = 3, для вещества Б — уравнение состояния

Р =(7 - —)рЕ, где 7 = 5/3.

В качестве краевого условия при х = 0 было взято значение и = 5.

В точке £ = 0, х = 0 образуется ударная волна с параметрами Р = 603.5, и = 5.0 В момент прихода этой ударной волны на прослойку образуются отраженная волна разрежения и прошедшая ударная волна. После ее выхода на вторую границу прослойки происхо-

т

Рис. 3.

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

Результаты расчетов показаны на рис. 3, 4 и в табл. 2. На рис. 3 схематически изображена получаемая в задаче 2 волновая картина в переменных те, а на рис. 4 — состояния

и

Рис. 4.

в прослойке и на прошедшей и отраженной волнах в переменных Р, и.

Таблица 2

№пп. X Р и х - 50.1 ^ = 0.1 , Р - Рг * = Рг

1 50.00 603.5 5.000 -1 0

2 50.10 104.6 8.860 0 -1.17

3 50.18 339.9 3.525 0.8 -0.44

4 50.21 554.1 4.754 1.1 -0.082

5 50.41 576.1 4.865 3.1 -0.045

6 50.71 713.1 5.508 6.1 0.18

7 51.00 669.9 5.313 9.0 0.11

8 51.33 637.3 5.161 12.3 0.056

9 51.91 605.1 5.007 18.1 0.0026

10 52.22 600.4 4.984 21.2 -0.0051

11 56.99 597.8 4.972 68.9 -0.0094

12 58.83 607.8 5.020 87.3 0.0071

13 60.36 605.7 5.010 102.6 0.0036

14 69.51 604.4 5.004 194.1 0.0015

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

Расчеты показывают, что в веществе А за прослойкой есть область, где давление на фронте прошедшей ударной волны ниже, чем на фронте падающей (около пяти толщин прослойки), и это различие достигает 44 %. Затем следует область шириной примерно 15 толщин прослойки, где давление на фронте ударной волны выше, чем на фронте падающей. Это превышение достигает 18 %. Далее вплоть до 200 толщин прослойки имеется область отрицательного и положительного отклонения от амплитуды падающей волны, где отличия не превышают 1 %.

Задача решалась по РС 1 из [2] и по изложенной выше схеме до выхода ударной волны на правую границу системы с координатой х = 100. В этой схеме затраты времени ЭВМ для решения задачи оказались в ~ 50 раз меньше.

Возможности изложенной РС проверялись также при расчете волны разрежения в задаче с начальными данными 0 < х < 1, Р0 = 2,р0 = 1.5, Е0 = 0.666, и0 = 0 при £ = 0, краевыми условиями и(0,£) = 0, и(1,£) = 1 и уравнением состояния Р = (7 — 1)рЕ с

7 = 3.

На рис. 5. приводятся профили давления Р и скорости и на один из моментов времени (--точное решение, о — результаты расчетов с та/к = 2 при к = 0.02).

Рис. 5.

Устойчивость этой системы проверялась методом гармоник в акустическом случае и с помощью расчетов в нелинейном случае. Рассчитывалось распространение малых возмущений в покоящейся среде с параметрами: и0 = 0, р0 = 1, Ро = 1, Е0 = 0.5, Р =

(7 — 1)рЕ, 7 = 3, 0 < т < 1. В трех точках (т = 0.25, 0.5, 0.75) значения Е0 и Р0 отлича-

та

лись от указанных на 1 %. Счет велся с — = 10, при к = 0.02. Через 20 шагов по времени

к

возмущения "размазались"и амплитуда возмущений уменьшилась до 0.4%.

Список литературы

[1] РихтмАйЕр Р., Мортон К. Разностные методы решения краевых задач. Мир, М., 1972.

[2] Рождественский Б. Л., Яненко Н. Н. Системы квазилинейных уравнений. Наука, М., 1978.

[3] Куропатенко В. Ф. О разностных методах для уравнений гидродинамики. Труды МИ АН СССР им. В.А. Стеклова, 74, ч. 1, 1966.

[4] Гаджиева В. В., Куропатенко В. Ф. О некоторых явных разностных схемах для уравнений гидродинамики. Журн. вычисл. матем. и матем. физ., 11, №6, 1971, 1603.

[5] Куропатенко В. Ф., ЕськовА Т. Е., Коваленко Г. В., Кузнецова В. И., Михайлова Г. И., ПотАпкин Б. К., САпожниковА Г. Н. Комплекс программ "Вол-на"и неоднородный разностный метод для расчета неустановившихся движений сжимаемых сплошных сред. Вопросы атомной науки и техники, Сер. Математическое моделирование физических процессов, вып. 2, 1989, 9-25.

Поступила в редакцию 8 августа 1996 г.