Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
§ 23] РЕШЕНИЕ ДВУМЕРНЫХ ЗАДАЧ ГАЗОВОЙ ДИНАМИКИ 363
Уравнения. Исходной для аппроксимации выбирается эйлерова дивергентная форма уравнений газовой динамики (3). Система координат, естественно, связана с обтекаемым телом. На краевых условиях не останавливаемся (это — тема отдельного разговора).
Сетка. Область расчета (обычно, прямоугольник) покрывается равномерной (для простоты) сеткой, ячейки которой нумеруются парами индексов (г, j).
Счетные величины. Состояние среды описывается сеточными функциями и1} j, и? j, р 1j, VZ11 j. Эти величины относятся к центрам ячеек и представляют собой приближенные значения компонент скорости и, V, плотности вещества р и плотности полной энергии
e + (u2 + v2)/ 2. Уравнение состояния используется в виде
р = Р(е, р), где е = w — (и2 + v2)/2. В дальнейшем мы будем ис-
пользовать величины типа pf\ j, понимая под ними вспомогательные числа, полученные из уравнения состояния очевидным образом.
Шаги по х и у считаем, для простоты, равными и обозначаем h, шаг по времени — т, хотя он, конечно, не фиксирован, а выбирается на каждом слое в зависимости от реализовавшихся значений и, у, р, w (из условий устойчивости и прочих). Схема метода крупных частиц явная. Как и в PIC-методе, стандартный шаг численного интегрирования состоит из двух этапов:
I этап: (и, v, р, w)?y-* (м, v, р, w)iJt
II этап: (м, v, р, w)t (и, v, р, н')'^1.
На первом этапе учитываем силы давления, пренебрегая переносом. Используем простую аппроксимацию уравнений (10):
7 (Pu-P"/) =0’ '
7 Plj(H) - </) +1 (Рї+uxj - tf-i/2,/) = °>
\ Plj(vIj-vIj) +І (Plj+ii2- Plj-m) = °>
где р"+1/2,у. Plj+т — полусуммы значений в центрах ячеек. Уравнение энергии имеет вид
7 Pl j^ij - Kj) + J [Срм)7+і/2,у - (PU)I- 1/2,у] +
+ Jl(Py)Ij^n - (Pv)Ij-1/2)=0-
364
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
-1/л і
Здесь возможны варианты: можно р, и, v порознь интерполировать с центров ячеек на их стороны, а можно интерполировать произведения (ри), (pv). Таким образом, первый этап очень прост и не содержит каких-либо нестандартных приемов аппроксимации.
Несколько сложнее и своеобразнее реализация второго этапа. На этом этапе учитываются процессы переноса: схема зависит от направления потока в данной точке и приобретает явно несимметричных характер. Второй этап начинается вычислением скоростей на сторонах ячейки: u*i+ll2j = Ui j + ui+1 j, v)j+ll2. Они используются только для определения направления потока. Затем вычисляются скорости в серединах сторон ячейки на основе отрезка ряда Тейлора: й(х± h/2) = й(х) ± (h/2)ux.
Однако при вычислении ui+112 у, например, можно использовать разложение как в точке (г, j), так и в точке (г + 1, /') — это определяется направлением потока. Предпочтение отдается тому направлению, откуда «приносится информация», т.е. откуда течет газ, попадающий в точку (г + 1/2, /'). В результате мы получаем
Н) + (“,+W ” “i-i.y)/4> uUuitj > О, “i+i.y - («і+2,7 - «i.yV4» uUiIlj < °-
Таким образом, для аппроксимации их используется центральная разность (второй порядок точности) в точке разложения. Вышеизложенный принцип вычисления величин в точках (г + 1/2, /') и (г, j + 1/2) (в последнем случае, очевидно, играет роль знак v] j+l/2) используется для вычисления всех остальных величин, фигурирующих в формулах для потоков через соответствующую сторону ячейки. Вводится, однако, дополнительная корректировка: если знаки и]+ij2j и мі+1ду противоположны, потоки всех величин (массы, импульса, энергии) через сторону ячейки (і + 1/2, /') считаются равными нулю (аналогично для потоков П( у+1/'2).
Теперь уравнение для р аппроксимируется следующим образом:
A2Cp^y1 ~ Pij) 4- (П/ + 1/2 у — П;_1/2^.) + (П; у + 1/2 — П; у._1/2) = О, где
П/ + 1/2,у = At М/ + 1/2J Р/ + 1/2,у’ у + 1/2= hi Vij + i/2 Pi1 j + 1/2*
Здесь уравнение записано в форме, подчеркивающей связь с законом сохранения массы ячейки h2p; величины П имеют смысл потоков массы через границу ячейки за время шага т. Подчеркнем, что каждый поток П;+1/2_ у, например, вычисляется для разделяющей ячейки (г, у) и (г + I, j) стороны независимо от того, является ли
РЕШЕНИЕ ДВУМЕРНЫХ ЗАДАЧ ГАЗОВОЙ ДИНАМИКИ
365
она правой для одной ячейки или левой для другой. Это свойство обеспечивает дивергентность схемы.
Остальные уравнения (законов переноса импульса и полной энергии) имеют общую форму:
(Р Q)t + (PuQ)x + (PvQ)y = °>
где Q принимает значения и, v, w соответственно. Эти уравнения аппроксимируются по одной и той же схеме. Мы уже имеем значения потоков массы ри (на правой и левой границах ячейки) и pv (на верхней и нижней границах ячейки) — это величины П, вычисление которых продемонстрировано выше.
Определим теперь правила вычисления величин Q на этих же границах. При этом используем величины Q в центрах ячеек, где они вычисляются естественным образом: Qij = Ui j, и т.п. Итак,
Qi+1i2,j = {Qt,j при ui+U2J> 0; Qi+1,/ при й1+1/2/<0},
Qi,/+1/2= (Qi,/ при viJ+ll2>0; QiJ+i при Vi j+m<0}.
После этого уравнение аппроксимируется просто (напомним, что PjTy1 уже известно):
