CC BY
93
УДК 519.633.6:519.612.2
А.М. Бубенчиков, Д.К. Фирсов
РАЗНОСТНАЯ СХЕМА ДЛЯ ИНТЕГРИРОВАНИЯ УРАВНЕНИЙ НАВЬЕ-СТОКСА НЕСЖИМАЕМОЙ ЖИДКОСТИ
НА НЕРАЗНЕСЕННОЙ НЕОРТОГОНАЛЬНОЙ СЕТКЕ
В работе представлена конечно-разностная схема для уравнений Навье-Стокса записанных в переменных скорость-давление в криволинейной неортогональной системе координат, согласованной с границей рассматриваемой физической области. Дискретизация уравнений движения выполнена на неразнесенной сетке. Для аппроксимации системы уравнений был применен новый способ, основанный на алгебраическом разложении скоростей на две части. При аппроксимации уравнения неразрывности использовались направленные разности, а центрально-разностный оператор, аппроксимирующий градиент давления, был представлен в виде суммы двух разностных операторов, шаблоны которых направлены в разные стороны. Данный подход позволяет строить систему линейных уравнений, соответствующих разностной аппроксимации уравнений На-вье-Стокса с обратимой матрицей, и исключить появление фиктивных колебаний решения без использования регуляризации или стабилизирующих процедур. В работе приводятся примеры расчетов, демонстрирующие эффективность данного подхода при моделировании течений в трубах сложной геометрии.
Для решения задач гидродинамики широко используются конечно-разностные методы [1]. В задачах, описывающих течение несжимаемой жидкости в областях сложной формы или с подвижной границей, используют неортогональную криволинейную систему координат. Это приводит к необходимости применения неразнесенной сетки для расчета течений в переменных скорость-давление и, как следствие этого, к необходимости использования центрально-разностных аппроксимаций для производных первого порядка. Шаблоны таких аппроксимаций не содержат центральной точки. Поэтому матрицы центрально-разностных аппроксимаций для производных первого порядка, возникающих при дискретизации операторов div и grad на не-разнесенной сетке, имеют ненулевое ядро, содержащее в себе всевозможные периодические по пространственным координатам изменения искомых характеристик, что является одной из причин возникновения фиктивных колебаний расчетных величин и нестабильности решений [2]. Авторы обходят эту проблему, вводя регуляризатор [3-4] или стабилизирующий член в уравнение неразрывности, как, например, это было сделано в методе искусственной сжимаемости [1]. Надежная регуляризация требует априорной информации о распределениях искомых величин, что усложняет поиск решения в конкретных
случаях. В нашей работе предлагается новый способ конечно-разностной аппроксимации уравнений Навье-Стокса для несжимаемой жидкости, позволяющий без ввода регуляризирующих операторов строить матрицу системы линейных уравнений с нулевым ядром.
Конкретной конечно-разностной схеме, перенумеровав узлы сетки, в которых необходимо найти значения, и перенеся в правую часть узлы с известными значениями, можно поставить в соответствие систему алгебраических уравнений [5]. Тогда каждой расчетной величине ставится в соответствие некоторый вектор, а разностному оператору - квадратная матрица. Можно полагать, что конечно-разностная схема, выписанная на заданной сетке, покрывающей область, и соответствующая ей система алгебраических уравнений, есть одно и то же. Удобным способом нумерации узлов, которого мы придерживались при записи систем алгебраических уравнений, соответствующих предложенным разностным схемам, является нумерация элементов в строку, затем переход к другой строке по вертикали и затем в глубину. На сетке
& = = (у;,,,= /рйр,Ив = 1/Мр,/р = р,в = 13}
узлу с /2, /3 можно поставить в соответствие элемент вектора с номером п = /3 + 1)(М2 +1) + /2 +1) + /;.
Система определяющих уравнений
Уравнения Навье-Стокса для случая движения несжимаемой жидкости в декартовой системе координат имеют вид
(ду/ д . ^
р -+ —- (У]У )
д дх1
-цУ V + — р = 0
^ дх/у
д
д . (1)
^у(у) = — у. = 0,
дх.
где у - вектор скорости, заданный в каждой точке области и состоящий из компонент V ; р - давление; ц - вязкость; р - плотность; / -номер координатного направления декартовой системы координат.
На практике исследователям чаще всего приходится изучать пространственные течения в областях сложной формы, для которых не всегда удобно строить ортогональное преобразование исходной физической области Н на каноническую (стандартную) область S, например параллелепипед. В связи с этим возникает необходимость применения неортогональных криволинейных коорди-6
нат, согласованных с границей области течения S. Построим такой переход от декартовой системы координат х1, х2, х3, заданной в Н, к
~ 12 3
неортогональной криволинейной системе координат у , у , у , заданной в S, чтобы якобиан перехода был больше нуля. Формулы
д к д к дук
перехода в этом случае имеют вид: —г = п ,■ —г, где п ,■ = —".
дх1 ' дук ' дх1
Это дает возможность представить (1) в следующей форме:
и1 =п'кУк
{я..' я я Л д
-цУ V +гк^р = 0
^ 1г дук
дУ1 д 1 , 1 ' д -+—- (и V ) + и1 V —- 1п(3)
д ду1 ду1 у
СЗД = 1—{/и1)= 0;3 = 1/с1в1 (пкт) > 0, (2)
3 ду1
уУ = !-дГ (^^Л
3 дук ду1
Здесь и1 - скорости, касательные к криволинейной неортогональной
к ,1 к системе координат, £ - компоненты метрического тензора, пг -
коэффициенты замены координат. При проведении вычислений в неортогональной криволинейной системе координат удобней всего пользоваться неразнесенной равномерной сеткой, т.е. сеткой, в которой все расчетные величины вычисляются в одних и тех же узлах.
Новый способ дискретизации уравнений Навье-Стокса, основанный на алгебраическом разложении скоростей
Для исключения фиктивных колебаний расчетных величин при использовании неразнесенной сетки нами был предложен новый способ аппроксимации уравнений Навье-Стокса, основанный на алгебраическом разложении скоростей на две части. Полагая V = уСу + у2) , перепишем систему (1):
р{М±^( (1+ У2 )])}-цу 2[[ (1+ У2 )]+Ар = 0,(3)
(VI + У2)] = 7) + сИу(у2)] = 0. Граничные условия для у и у2 задаются в точности такие, каким должна удовлетворять физическая скорость V . В системе (3) зна-
чения v1 и v2 определяются неоднозначно, однако суммарное решение V единственно, если (1) имеет единственное решение.
Пусть О - равномерная сетка, покрывающая область 5" с постоянным шагом И1 по координатному направлению у1 . При построении дискретного аналога системы (3) на сетке О будем исходить из положений: (Уу ))О + (Еу(У2 ))О = + БУ2У2 + 01(Ит) + 02 (Ит),
где Бук = (/J ^5^2 -1 ),(к = 1,2) - блочные мат-
рицы, соответствующие аппроксимации операторов сНу в пространственном случае, в которых для расчета производных д/ду применены направленные разности 51 и 5 2 соответственно. Операторы 51 и 5^ выберем так, чтобы шаблоны производных были направлены в разные
стороны. Например, для Оу1 имеем
' д >
ду,
V ^1
1 п
51 = ~г УвТ , а для
О И1 1=0
ВУ2
( д ^
ду1
» 52. =—У-р.7-1, где Т/ - оператор сдвига по расчет; Уо И1 /=о
ной сетке в направлении. на / узлов, в. - коэффициентах направленной разности. Если шаблон направленной разности в приграничной точке выходит за пределы расчетной области, то здесь значения расщепленных скоростей полагаются равными значениям на границе. Данный способ приводит к тому, что матрицы, отвечающие аппроксимациям д/ду, становятся треугольными. Это благоприятно сказывается на обусловленности матрицы всей конечно-разностной схемы. В случае нецентрирован-ных компактных разностей аппроксимацию можно проводить так [6]:
Ак} (з(к)) = (е +1 € -7*(к)€0)-1 X 0,5(50 - s(k)Ц)/Н], (4)
*(к) ^^ € = № -2Е + Т-1) и €0 = "2(Т11 -Т-1),
используя в приграничных узлах конечные разности, направленные в ту же сторону и имеющие тот же порядок точности, что и компактные.
Аппроксимация оператора ^уас1 проводится с использованием центральных разностей. Общий вид такой аппроксимации для любых производных первого порядка может быть задан соотношением:
1 " ( ■ ■ \ " Дг = у Хак Тк - Т-к ), причем Ь = Хак 0 , где а к - коэффициен-
к=1 к=1
ты шаблона (в случае, когда к=1 ак = ). Разложим оператор Дг, аппроксимирующий первую производную вдоль координатного направления 1, на два: Дг = ? (Д? + Д2), где шаблоны Д? и Д2 направлены в разные стороны, таким образом, чтобы Д1х = 0, если х - константа. Поэтому элементы разложения центрально-разностного оператора, переводящие константу в нуль, будут следующими:
Д = -21 ЬЕ - ±акгк 1, Д2 = У| ЬЕ - ^акТ-к
к=1
(5)
В случае компактных центрированных разностных аппроксимаций, например четвертого порядка точности \7], которые имеют вид Дг = е-151, где б = | Е +1(Т? + ТД), а 51 = -^(Т/ - ТД),
2П:
мы можем применить следующее разложение:
Д = б-1 (Т? - Е), Д2 = б-1 1 (Е - Т-1).
у у
(6)
Введем обозначение: Б к =
Ч Дк11
п2 Дк1 п3 Дк,
V
(к = 1,2) и Б = + Б2)
- блочная матрица центрально-разностного оператора аппрок-
симации градиента давления. Под У" =
'(у 1)" 1
(У2) (У3)
будем понимать
приближенное решение для скоростей на сетке О системы дифференциальных уравнений (3) на слое " по времени, а под
'(У ■) Л
(к=1,2) - приближенные решения для V к. Обозна-
У" =
(у 2)
(У3)
(
чим через А" =
О" 0 0 0 О" 0
0 0 Оп
Л
разностный оператор, отвечающим
слою п по времени, по главной диагонали которого стоят О, где О"(р*)п+1 - некоторая разностная аппроксимация следующей конвективно-диффузионной части:
1А V* + -д- ((V 1п( 3) V I
д
\
д г д
р[ — V +— .. . ^
[д^ ду] у ду] ^ [ ду] у ду
Матрица О" имеет размер МхМ, где М - количество точек, в
которых рассчитываются скорости. Число блоков О" на главной
диагонали А" совпадает с размерностью решаемой задачи. Дискретизацию системы (3) в алгебраическом виде, отвечающую неявной схеме, можно записать следующим образом:
(
А
А
1Щ " 0у2 0
А" 0
0
А"
Утл"+Л (Ъ" Л
"+1
"+1
У"+1 =■
, /V А
:(("+1 + V"+1)
VI
^2
У ЪР ,
у р /
(7)
где Ъ
у , Ъуг возникают за счет дискретизации производной по времени, граничных условий для скорости и переноса в правую часть узлов
т. "
с известными значениями. ЪР возникает только за счет граничных условий для скорости и переноса в правую часть узлов с известными
значениями. В системе все расчетные величины V
"+1
V
"+1
Р"
являются векторами в том смысле, что каждой расчетной величине, заданной на некоторой сетке, перенумеровав узлы, можно поставить в соответствие вектор. Таким образом, если скорости вычисляются в М узлах, а давление в N узлах, то размерность матрицы в системе (7) будет равна (6M+N)х(6M+N). Разрешив систему (7), получим: Г А" ± у+1 + V""+1) + ВР"+1 = (ъ" + Ъ"2) ;
[^фуг-1 + ^уг1) = ър. 1 2
Давление находится с порядком точности центральных разностей, аппроксимирующих оператор градиента, а общий порядок схемы - мини-
2
мальныи в аппроксимациях членов уравнении движения и неразрывности, в данной статье первый за счет аппроксимации конвективной части.
Систему линейных уравнений (7) можно решить разными способами: 1) решить непосредственно систему (7); 2) построить итерационный процесс, аналогичный алгоритму SIMPLE [2]; 3) сконструировать альтернативный указанным способам алгоритм.
Первый способ не эффективен для расчета стационарных задач. Второй способ эффективен при решении стационарных задач и менее сложен в реализации, нежели первый. Для упрощения (7) введем обозначе-
' A" 0 1 B" = (B" 0
0 A" J' Я ( 0 B"
л
ния: A =
, где B - некоторая матрица;
D = ( D2
Dv = ((vi,DV2), V" =
fVn 1
V"
V 2
T"
bv =
"
Vi
\V2J
. Тогда система
(7) будет выглядеть следующим образом:
A" Dv
D 0
f—" 1 f-" 1 bV
V p"
V bP J
После проведенных оценок сходимости итерационного процесса, аналогичного SIMPLE [2], для системы (7) было сформулировано Утверждение: если после наложения граничных условий
имеет место следующее: матрицы A и B - положительно определенные; спектр F = D(C")-1 Dv состоит лишь из отрицательных чисел и нуля; матрица Dv(C")-1 D имеет нулевое ядро и
справедливо неравенство
(B" )-1(A" -1F)
< 2, где
= max J (Ax, Ax)
(X, x)=1 V
то итерационный процесс:
7"+1,m+1
V
C" (P
V"
= (E - (B" )-1 A")V"+1,m - (B" )-1 DP" + (B")-1b
"+1, m+1 p"+1,m
= "2 (V
" +1, m +1
) = - ± Dv-V
" +1, m +1 + V2" +1, m +1 )
+ bp
(8)
сходится к решению
f v"+11 p"+1
системы (7). В записи (8) m - номер
внутреннего итерационного слоя.
2
2
Доказательство утверждения. Опустив для краткости часть индексов, запишем итерационный процесс (8) в следующем виде:
( уп+1,т+1 \
Р
п+1,т+1
= (Е - Ж)
( +1 |
Р
п,т+1
+
в
СЛр - С-ОУ • в~1Ь"Р
(9)
(
где Ж =
В- А
О
V 2
ОУ(Е - В А) - { С- БУВ- О
Воспользовавшись тем, что Ц |2 совпадает с максимальным
по модулю собственным значением матрицы Ж, докажем, что все собственные значения матрицы | больше нуля и
ЖI <
(В )-1( лп - -2 ^)
. Для этого покажем, что все собствен-
ные значения матрицы Ж больше либо равны нулю:
¡XV = В - АУ + В-1 ОР,
[АР = С^ ОУ(Е - В-1 А)У - С-1-1 БУВ ОР; ГхУ = В-1 ~АУ + В-1 ОР,
-1—т
[АР = С-1 1 БУУ - АС- ± БУУ;
V = | X-1ВОС- БУ X 2
В - (А -1 ОС-1 ОУ) IV.
Так как А и В - положительно определенные, спектр ^ = ОС - Оу состоит из отрицательных чисел и нуля, то все собственные значения матрицы Ж либо положительны, либо нуль. Поэтому
справедливо следующее неравенство |Ж|I <
(Вп)-1(Ап -1 ^)
Матрица БуВ - О обратима и можно показать обратимость Ж и отсутствие нулевого собственного значения. Из предыдущих фактов следует, что все собственные значения Ж больше нуля. Итерационный процесс (8) будет сходиться, если матрица Е-Ж сжимающая (все собственные значения Е-Ж меньше единицы). Поскольку все собственные значения матрицы Ж больше нуля,
итерационный процесс (8) сходится, если
ВА -1 ^)
< 2.
2
2
Численная реализация
Для того чтобы проиллюстрировать данный подход, были рассмотрены течения жидкости в трубах сложной геометрии. Причиной выбора этих примеров является необходимость интегрировать полные уравнения Навье-Стокса для ньютоновской жидкости на неразнесенной сетке, для которых при обычных подходах требуется вводить регуляризирующий оператор или применять методы установления. Применяемый здесь подход позволяет избежать данных проблем.
При реализации алгоритма для обеспечения сходимости в качестве
СП 7~» П
и В можно использовать как сложные операторы, например ВП = АП, так и упрощенные, например константы. Константы стоит использовать лишь в случае незначительного количества узлов. Для исследования свойств полученной схемы в качестве обеспечивающих сходимость операторов мы пользовались как константами, так и имеющими сложную структуру матрицами. Согласованную с границами области расчетную сетку получали с помощью эллиптического сеточного генератора [8] из имеющего в вычислительном пространстве размер 1x1x1 и покрытого равномерной декартовой сеткой куба:
о = {у = (у,у\,у\)| ув = 'Л,Нр = 1/Щ,/р = ,р=1з}.
Коэффициенты замены системы координат п/ = ду/ / дХ' вычислялись с помощью центральных разностей второго порядка точности во внутренних узлах и трех точечных направленных разностей второго порядка в граничных узлах.
В настоящей работе не ставилась задача решения проблемы схемной вязкости [9] и при аппроксимации конвективной части уравнений Навье-Стокса нами была использована противопоточная схема, имеющая первый порядок точности. Конечно-разностная аппроксимация, полученная так же, как в [2], при этом имела вид:
дУ д , г ( ■ д
--1--Т и'у + и —
дt ду/ I ду
■ ( д ■ д ^ 1п(.Т) V - V —и + и/ —г 1п(.Т)
■
V ) V
ду ■ ду
:Ж(Г )+1 - Ь' =[арЕ -X [ + а/Т-1 ]] - Ь'
V /=° )
о
где
а]Е /2, ¿3] = Мах
аш [/, /2, г3] = Мах
(и1 )" [/1,/2, /3] + Т (и1)"
[/1, /2 , /3 ]
2к,
(и1 )" [/1, /2, /3]+\г-1 (и1У
[/1, /2, У
/ л
2Л,.
а
з 1
р =е (ае + а]ш) + -, ьг =
1=1
И.
т
Под коэффициентами аР, аЕ , аШ мы понимаем диагональные матрицах, а под аР [/1, /2, /3], а1. [/, /2, /3] - элементы, стоящие на главной диагонали матриц аР, а. (п = Е,Ш) и одновременно соответствующие внутренним узлам 0<11<Ы1,0<12<Ы2, 0<3<¥3 расчетной сетки.
Диффузионную часть мы приближали во внутренних узлах конечными разностями второго порядка точности относительно пространственного шага:
1 _д_ / дук
(
д
\
3
33
к= Ш4ЛЛ)
^ Т1 - тк/кь - т$ ](), сЕ = (Е+Т^1,1,
(е+Т- -1
( з
СрЕ (сЕТ+1 + 4Т+) +
I 1=1
ш
а с; =■
2^2
= Е (сЕ + сШ).
=1
Такой способ аппроксимации конвективно-диффузионной части
позволяет строить Ап положительно определенной и обладающей свойством диагонального преобладания.
Для получения давления в граничных узлах применяют экстрапо-ляционные формулы [10]. Такая технология может оказаться очень трудоемкой, поскольку уравнение неразрывности в граничных точках не удовлетворяется и давление в них в пространственном случае при
0
0
т
использовании криволинейной неортогональной системы координат зачастую следует получать за счет матричной прогонки. Поэтому давление на границе мы находим с использованием разностного аналога уравнения неразрывности, полученного с применением разностей, направленных в ту же сторону, что и во внутренних узлах. Значения в точках, лежащих вне области, мы считали равными значениям на границе. В этом случае в граничных точках одна из расщепленных скоростей может удовлетворять уравнению неразрывности, тогда как вторая будет удовлетворять этому уравнению лишь в том случае, когда полная скорость является решением уравнения неразрывности. Все это позволяет надежно определить давление в граничных узлах.
При реализации алгоритма, основанного на разложении скоростей, выбор направленных разностей в операторах Оу, входящих в (7), в общем случае должен быть таким, чтобы их порядок аппроксимации не уменьшал порядка аппроксимации всей системы. Несложное доказательство того, что спектр О(Сп)-1 Оу состоит лишь из отрицательных чисел и нуля, основывается на том, что якобиан замены координат больше нуля; матрицы 5к и Д/ благодаря выбранному способу нумерации узлов являются треугольными. Причем на их диагоналях стоят положительные числа, если матрицы нижнетреугольные, и отрицательные, если матрицы верхнетреугольные. Легко доказать, что спектры треугольных матриц состоят лишь из элементов, находящихся на главной диагонали.
Для завершения доказательства необходимо расписать произведение О(Сп)-1 Оу, дополнив прямоугольные матрицы, входящие в него, нулями до квадратных. Затем оценить знаки собственных значений каждой матрицы в произведении.
При апробации метода использовались разложения центрально-разностных операторов второго порядка точности (5). В качестве направленных разностей при разностной записи оператора &у выбирались трехточечные аппроксимации первой производ-
S ( 3
ной со вторым порядком точности, т.е. 5/ = — V2Е - 27-5 + + ] , а Д/ = Н-( -Е), 5 = (-1)к+1, к = 1,2.
Чтобы показать независимость решения, полученного с помощью (8), от способа аппроксимации производных в операторах Шу и grad, были использованы компактные разности [6 —7], имеющие порядок аппроксимации больше, чем второй, и хорошо зарекомендовавшие себя в расчетах течений сжимаемой среды [11-12]. В качестве центрально-разностного оператора использовались центрированные компактные разности четвертого порядка точности [6] относительно шага по пространству, реализованные на двухточечном шаблоне, которые были разложены по формуле (6). В качестве направленных разностей выбирались нецентрированные компактные разности третьего порядка точности (4). Данные вычислений на неразнесенной сетке сопоставлялись с результатами, полученными с помощью алгоритма, реализованного на разнесенной сетке Аракавы, позволяющей рассчитывать давление в центрах параллелепипедов, а скорости - в их вершинах (рис. 1).
v
v v
Рис. 1. Элемент расчетной сетки
Аппроксимация конвективно-диффузионной части проводилась в последнем случае так же, как для алгоритма, использующего неразнесенную сетку. Кроме того, была использована обычная SIMPLE [2] процедура, как и в двухмерном случае при использовании разнесенной сетки Аракавы [1]. Давление на сетке Аракавы может быть нефизичным, но скорости при этом определяются однозначно.
Сравнение полей скорости показало практически полное совпадение результатов расчетов по алгоритму (8) и по алгоритму на сетке Аракавы для ряда различных гидродинамических ситуаций.
Пример расчетов на существенно неортогональной сетке
По предложенному алгоритму были проведены расчеты стационарного течения ньютоновской жидкости на существенно неортогональной сетке, полученной сеточным генератором из единичного куба декартовой сеткой, размер которой составлял 21x21x41 (рис. 2).
Рис. 2. Расчетная область - труба с несимметрическим вздутием, полученная из единичного куба с сеткой размером 21x21x41
В качестве операторов, обеспечивающих сходимость, использовались константы. Представленные ниже результаты соответ-
ствуют Re =-=400 (где р - плотность, ц - вязкость, ~ -
средняя скорость во входном профиле Пуазеля, а d - диаметр входного отверстия). На стенках ставилось условие прилипания
д
(v=0), а на выходе условие свободного вытекания (^ v = 0, где
дп
п - нормаль к сечению трубы на выходе). Полученные данные приведены в таблице и на рис. 3.
SIMPLE процедура для итерационной схемы (8)
Для расчета течений на сетках со значительным количеством узлов необходимо использовать более сложные предобусловливающие матрицы. Поэтому для численного анализа стационарных течений в протяженных трехмерных областях, требующих значительного количества узлов, был реализован эффективный алгоритм со сложными опера-
Скорости и3, касательные к криволинейной системе координат, д ля направления .у3 и при /=0,5; /=0,902, отнесенные к средней по сечению скорости на входе. Координатная линия расчетного пространства / направлена вдоль канала
Координаты у1 и3 - компактные аппроксимации в (7) и3 - центральные разности в (7) и3 - разнесенная сетка Аракавы
0,95 0,840 0,869 0,817
0,286 1,906 1,912 1,907
0,571 2,073 2,071 2,101
0,762 1,840 1,841 1,884
0,952 1,110 1,116 1,077
1.0 0.8 0.6 0.4 0.2 0.0
о Разнесенная сетка
-Компактные разности в (7)
х Обычные разности в (7)
у =0.02 у=0.12 у =0.24 у =0.46 у =0.68 у =0.90
Рис. 3. Сравнение результатов расчетов, полученных различными способами. Профили скоростей и3 для у2=0,5, а координатная линия у3 направлена вдоль канала
торами, обеспечивающими относительно быструю сходимость и второй порядок аппроксимации для уравнения неразрывности, а также для разностного оператора градиента давления, где, как и в предыдущем случае, в качестве 5к, мы использовали — Е - 27-— + 2Т-2— ), а
И
Ду = — {г— -Е), где — = (-1)к+1, к=1 или 2. За счет противопоточ-
ной дискретизации конвективных членов матрица А п обладает свойством диагонального преобладания и поэтому в качестве Вп в (8) использовалась матрица Ап. 18
Чтобы получить давление во всех точках расчетной области, необходимо сначала по алгоритму (8) рассчитать величину поправки давления. При вычислении этой поправки мы ввели фиктивные узлы, лежащие вне расчетной области. Так, если равномерная сетка по расчетной области в трехмерном случае имела вид Q = {y = (y^, y2, y33 )|
У^р = iphp, hp = 1/ Np, /р = 0, Nß, ß = 1,3} то поправка давления определялась в следующей, расширенной фиктивными узлами области:
^ = = (У1, У2, у33)|ув = /ehe, hß = 1/Np, /р = -1, Nß+1, ß = 1,3}
Во вспомогательных (фиктивных) узлах поправка предполагалась нулевой. В качестве C" мы брали матрицу, близкую к -1 (Dy(d(An))-1 D + Dv2(d(An))-1 D2), где d(An ) - диагональ An :
C" =-L J-11 (((+-1 - E )gJJd ( An )-1 (E - Tj )+
2 j=1
+ (E-Tj)jgj,jd(An)-1(t+, -E)) .
Здесь K— размерность решаемой задачи; L - некоторый эмпирически подбираемый параметр (в большинстве случаев мы полагали L=2),
компенсирующий отсутствие в C аппроксимаций смешанных производных и разницу порядка аппроксимации уравнения неразрывности в Cn и в матрице 2(Dy(d(An))D1 + Dv2(d(An))D2); значения
коэффициентов матрицы Cn в фиктивных узлах брались такими же, как в граничных точках разностной сетки. В результате в трехмерном
Cn ~
получилась семидиагональной, положительно определенной и за счет ввода фиктивных узлов обратимой.
Для решения системы уравнений с участием Bn применялся метод простой итерации [2], а для решения системы линейных уравнений с участием Cn - явный метод Булеева (ЯМБ-3) [13]. Для увеличения скорости сходимости итерационной схемы, так же как в алгоритме SIMPLE [2], на основе поправки давления вычислялись поправки скоростей. В итоге получился следующий алгоритм для расчета стационарных течений несжимаемой жидко-
сти, в котором: n - номер итерации; Д0 j = T+j - T-1)/ 2hj - центрально-разностная аппроксимация производной д/ду1 ; V1 - вектор-столбец i компоненты вектора скорости; K - размерность решаемой задачи; (bVk )j - вектор, получаемый за счет дискретизации по времени и задания граничных условий для VkJ ; Gn -аппроксимация конвективно-диффузионной части.
V> = (Gn)-1 ((bnvk)j - nLДр) k = 1,2; j = 1K;
P' = (Cn)-1 ±/+ /8/nV - bp j;
pV + 1 = Pn - p'; __(10)
Vj = :2(Vj + V/ ), j = 1,K _
(Vj )n+1 = V] + ((Gn )) n™ До,тР' / L, j = 1, K.
Итерационный процесс (10) необходимо продолжать до тех пор, пока норма разницы между значениями на двух итерационных слоях не станет меньше некоторого наперед заданного числа и в то же самое время не выполнится уравнение неразрывности.
В полученном алгоритме (10) для расчета стационарного течения, описанного в «Примере расчетов на существенно неортогональной сетке», потребовалось 930 итераций. При этом
3 N N2 N3 | , . , ,1
SS S S [11, 12, 1з]) - ( [11, 12,13]) < 10-3, а максимум не-
1 =1 11 =0 12 =0 13 =0 '
вязки уравнения неразрывности был меньше 10-4. Разница расхода на входе и выходе составила 0,5 %. На персональном компьютере с одним процессором AMD DURON 700 МГц для расчета потребовалась 31 минута, что примерно в 10 раз меньше, чем при выборе констант в качестве операторов, обеспечивающих сходимость.
Расчеты по алгоритму (10) сравнивались с экспериментальными данными [14]. В эксперименте изучалось движение вязкой несжимаемой жидкости в тороидальной трубе для различных чисел Рейнольдса. Угол разворота ф сегмента тора составлял 180°, а соотношение радиуса кривизны оси тора к радиусу трубы составляло 29,1. Перед входом в трубу имелся прямолинейный участок, который обеспечивал формирование профиля Пуазейля. 20
Расчеты выполнялись для числа Рейнольдса 1060 на сетке размером 41x41x61, покрывающей рассматриваемый сегмент тора. На входе был задан параболический профиль Пуазейля, на стенках ставилось условие прилипания, на выходе - условие свободного вытекания. Результаты оказались достаточно близкими к эксперименту (рис. 4).
■ - Эксперимент
Рис. 4. Сравнение расчетов с экспериментом uz / Uz - соотношение проекции скорости на касательное направление к оси тора по отношению к средней скорости течения на входе; a - радиус трубы
В результате исследований нами был разработан новый подход к дискретизации уравнений Навье-Стокса, который позволяет строить невырожденную матрицу линейной системы уравнений, соответствующую предложенной разностной схеме, что делает алгоритм надежным и исключает фиктивные колебания расчетных величин по расчетной области. Хотя данный подход и приводит к увеличению размеров системы линейных уравнений, но он не требует применения регуляризирующих операторов и методов установления.
ЛИТЕРАТУРА
1. Пейре Роже, Тейлор Томас Д. Вычислительные методы в задачах механики жидкости. Л.: Гидрометеоиздат, 1986.
2. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энергоатомиздат, 1984. 152 с.
3. Rhie C.M., Chow W.L. Numerical Study of the Turbulent Flow Past an Airfoil with Trailing Edge Separation // AIAAA Journal, 1983. Vol. 21, No.11. P. 1525-1532.
4. P1que /., Vasseur X. Multigrid Preconditioned Krylov Subspace Method for Three-dimensional Numerical Solutions of the Incompressible Navier-Stokes Equations // Numerical Algorithms. 1998. Vol. 17. 1, 2. P. 1-32.
5. HagemanL.A. and YoungD.V. Applied Iterative Methods. Academic Press. New York, 1981.
6. Толстых А.И. Об итерационных схемах с нецентрированными компактными аппроксимациями // ДАН. Математика. 1992. Т. 326, № 3. C. 425-430.
7. Толстых А.И. Компактные разностные схемы и их применение в задачах аэрогидродинамики. М.: Наука, 1990. 230 с.
8. Steger /.L., Sorenson R.L. Automatic mesh-point clustering near a boundary in grid generation with elliptic partial differential equations // J. comput. phys. 1979. Vol. 33. P. 405-510.
9. MalUnson Dav1s. False diffusion in numerical fluid mechanics. Univ., of new south woks, School of mechanics and industry. Engineering report. 1972. FMT 1.
10. Захаренков М.Н. Особенности разностных схем решения двумерных уравнений Навье-Стокса, связанные с постановкой граничных условий на твердой поверхности // Журнал вычислительной математики и мат. физики. 1990. Т. 30, № 8. С. 1224-1236.
11. Савельев А.Д. Расчеты течений вязкого газа на основе компактных схем третьего порядка // Журнал вычислительной математики и мат. физики. 1995. Т. 35, № 10. С. 1538-1551.
12. Толстых А.И., Широбоков ДА. О разностных схемах с компактными разностями пятого порядка для пространственных течений вязкого газа // Журнал вычислительной математики и мат. физики. 1996. Т. 36, № 4. С. 71-85.
13. Ильин В.П. Методы неполной факторизации для решения алгебраических систем. М.: Физматлит, 1995. 288 с.
14. Patankar S. V., Pratap V.S. and Spaldng D.B. Prediction of laminar flow and heat transfer in helically coiled pipes // J. Fluid Mech. 1974. Vol. 62, part 3. P. 539-551.
