ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2011, том 51, № 11, с. 2075-2083
УДК 519.633
РАСЧЕТ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ НА НЕСТРУКТУРИРОВАННЫХ КРИВОЛИНЕЙНЫХ СЕТКАХ1)
© 2011 г. В. М. Головизнин*, В. Н. Котеров**, В.М. Кривцов**
(* 113191 Москва, ул. Б. Тульская, 52, ИБРАЭ РАН; ** 119991 Москва, ул. Вавилова, 40, ВЦ РАН) e-mail:gol@ibrae.as.ru; koterov@ccas.ru; krivtsov@ccas.ru Поступила в редакцию 25.03.2011 г.
Переработанный вариант 23.05.2011 г.
Представлен один вычислительный подход к расчету уравнения теплопроводности, отличающийся в случае трехмерных косоугольных (неортогональных) неструктурированных сеток компактностью сеточного шаблона и безусловной устойчивостью численного алгоритма. Особенностью подхода является использование потоковых функций как самостоятельных зависимых переменных. Рассматриваются в основном гексагональные сетки, каждая ячейка которых может быть непрерывным образом отображена на единичный куб. Приводятся примеры расчетов. Библ. 7. Фиг. 3.
Ключевые слова: численные методы, уравнение теплопроводности, потоковые переменные, неортогональные расчетные сетки, неструктурированные сетки.
1. ВВЕДЕНИЕ
Алгоритмы численного интегрирования уравнений математической физики на криволинейных расчетных сетках развиваются давно и достаточно успешно. Среди различных подходов к построению расчетных схем на косоугольных (неортогональных) сетках можно отметить метод конечных элементов (см. [1]), метод конечных объемов (см. [2]) и метод опорных операторов (см. [3]).
Метод конечных элементов хорошо зарекомендовал себя на линейных параболических и эллиптических уравнениях, но его применение затруднено в тех случаях, когда рассматриваемая задача существенно нелинейна. Это, например, является одной из причин, по которым данный метод не нашел широкого применения в задачах подземной гидродинамики.
Метод конечных объемов (интегроинтерполяционный подход) исходит из аппроксимации законов сохранения и является наиболее универсальным. Однако эта универсальность не позволяет, вообще говоря, строить внутренне согласованные (сопряженные) аппроксимации таких дифференциальных операторов, как дивергенция и градиент, что приводит к отсутствию у разностного аналога оператора Лапласа свойства самосопряженности и вызывает дополнительные вычислительные проблемы при его обращении.
Указанного недостатка лишен метод опорных операторов, который в последнее время и используется для разработок самосогласованных дискретных моделей сплошных сред. Однако при использовании данного метода в случае кубических сеток оператор Лапласа вырождается в так называемый "косой крест", что приводит к расщеплению решения на несколько независимых мод, смещение которых относительно друг друга обусловливает немонотонность решения. Корректная аппроксимация потоков на криволинейных неортогональных сетках, получаемая методом опорных операторов, приводит к существенному увеличению вычислительных шаблонов, что оказывает негативное влияние на скорость сходимости итерационных алгоритмов и заметно уменьшает эффективность алгоритмов распараллеливания.
Минимизации вычислительного шаблона можно достичь за счет расширения используемого набора переменных. Так, в работах [4], [5] наряду с неизвестными значениями температур в качестве основных переменных использовались новые переменные — так называемые тепловые смещения [6]. Предложенная в этих работах аппроксимация уравнения теплопроводности на двумер-
1)Работа выполнена при поддержке Программы № 3 фундаментальных исследований ОМН РАН, Федеральной целевой программы "Научные и научно-педагогические кадры инновационной России" (гос. контракт № 02.740.11.0041) и РФФИ (код проекта 11-01-00444-а).
2075
ных косоугольных сетках помимо компактности шаблона приводила к безусловно устойчивому и экономичному алгоритму численного решения, родственному методам переменных направлений.
Последующий анализ показал, что можно отказаться от тепловых смещений и записать разностные уравнения с сохранением всех достигнутых позитивных качеств в терминах температур и тепловых потоков. Настоящая работа посвящена распространению подобного подхода на случай неструктурированных косоугольных расчетных сеток. В статье рассматриваются трехмерные гексагональные сетки, хотя обобщение на случай сеток из произвольных выпуклых многогранников очевидно и, как представляется, не должно вызывать затруднений.
Далее будет рассматриваться следующая эволюционная задача расчета распределения температуры T в трехмерной области Q с границей Г, на которой заданы общие условия III рода:
cV (r)— + divq = F(r), q = -a(r)gradT, r eO, д t
q • n + к(r)[T - Tw (r)] = Qw (r), r e Г, (1.1)
T = T0, t = 0, r e O
Здесь t — время, г — пространственные координаты, n — внутренняя нормаль к границе Г, q — вектор теплового потока, cV — теплоемкость, а — коэффициент теплопроводности, F — функция объемных тепловых источников, Tw — температура поверхности Г, Qw — тепловой поток через Г, к — коэффициент теплообмена с поверхностью. Если на части границы Г1 с Г задано условие I рода, то Qw = 0, к = да при г е Г1. Если на Г2 с Г задано условие II рода, то к = 0 при г е Г2.
2. РАСЧЕТНАЯ СХЕМА С КОМПАКТНОЙ АППРОКСИМАЦИЕИ ПОТОКОВ НА КРИВОЛИНЕЙНЫХ НЕСТРУКТУРИРОВАННЫХ СЕТКАХ
Покрытие трехмерной области неструктурированной расчетной сеткой с гексагональными ячейками порождает три конечных упорядоченных множества, а именно: множество узлов, нумеруемых далее целым числом n, множество ячеек, нумеруемых числом с, и множество граней, нумеруемых числом f С каждой ячейкой могут быть ассоциированы значения температуры T, теплоемкости cV, коэффициента теплопроводности а и функции тепловых источников F. Каждая гексагональная ячейка содержит 12 ребер, 6 граней и 8 вершин, в которых находятся узлы разностной сетки. Грани ячейки в общем случае не являются плоскими поверхностями, но мы будем предполагать, что для каждой грани расстояние А между парой скрещивающихся диагоналей d13, d24 (фиг. 1a) мало по сравнению с ее характерным размером max (| d13|, |d 241) и грани являются "почти плоскими". При этом можно считать, что вектор S = 0.5d13 х d24 перпендикулярен "плоскости грани", единичный вектор s = S/|S| является нормалью к грани, а величина S (f) = |S| — площадь этой грани.
Одна из проблем, возникающих при аппроксимации задачи (1.1) на неструктурированной сетке, заключается в том, что для правильного представления вектора теплового потока q и корректной аппроксимации величины divq для каждой грани необходимо не только фиксировать направление нормали s, но и указывать связь этого единичного вектора с внутренними (или внешними) нормалями n к границам ячеек сетки. Для решения данной проблемы удобно для каждой грани f ввести понятие ячейки-хозяина и ячейки-соседа, назвав так две ячейки, соседствующие с гранью (см. фиг. 1б). Будем полагать, что каждая грань, лежащая внутри области, имеет как ячейку-хозяина, так и ячейку-соседа, а грани, образующие границу области — только ячейку-хозяина. Примем, что для каждой грани f нормаль s, порождаемая введенным выше вектором S, направлена в сторону ячейки-хозяина, и обозначим через Q (f) = q • s плотность теплового потока через эту грань. Далее поступим следующим образом.
Пусть c1f и с2f — номера ячейки-хозяина граниf и ее ячейки-соседа соответственно. Для произвольной пары грань—ячейка (f, с) определим "знаковую функцию" Z:
+ 1, C1 f = с, Z (f, c) = <!-1, C2f = c,
0 в остальных случаях.
Тогда n = s Z(f, с) — внутренняя нормаль к грани f ячейки с и проекция qn вектора теплового потока q на эту внутреннюю нормаль может быть записана в виде
(a)
(б)
Ячейка-хозяин Су грани f
/2
f=fl
\
А
Уз
Ячейка-сосед c2f грани f
n
з
з
2
n
Фиг. 1.
In = Q(f)Z(f, c).
При этом естественная и корректная аппроксимация дивергенции вектора q для ячейки c будет иметь следующий вид:
divq "-¿h 2 Q(f)s(f)Z(f, c). (.i)
V (C) f tIf (c)
Здесь If (с) — множество граней, образующих ячейку c, F(c) — объем рассматриваемой ячейки с, который с использованием тождества
V(c) = J dV = J idivгdV = -1 | г • ndS,
V(c) V (c) S(c)
обеспечивающего дивергентность аппроксимации объемов расчетных ячеек, может быть вычислен по формуле
V(c) -1
3 f ^If (c) 14
X г» (/) • 8 (/) я № (/, с)\.
В двух последних формулах Б(с) — поверхность объема К(с), г — радиус-вектор точки, п — внутренняя нормаль к поверхности S(c), 1п (/) — множество узлов п, принадлежащих грани /, гп (/) — радиус-векторы узлов, образующих грань/, 8 (/) — нормаль к грани/, направленная в сторону ячейки-хозяина.
Для аппроксимации закона Фурье (1.1) воспользуемся тем, что согласно [7] минимальное значение функционала
Ф = |(|42/ст - 2ТёЦ(Ю. + К^/к - 2Глп)(Г, ^ = 4 ■ п, (2.2)
п г
достигается при выполнении следующих соотношений:
4а = ^гаёТ, г еО, дп/к = Т - Тп, г еГ.
Здесь Т — гладкая функция, заданная в области ^ с гладкой границей Г.
Корректная аппроксимация функционала (2.2) будет являться источником расчетных формул, связывающих значения плотности теплового потока О через грани со значениями температуры Т в ячейках. Выполним ее.
Аппроксимация величины ёгу q, входящей в (2.2), дается формулой (2.1). Для построения дискретной аппроксимации квадрата модуля вектора q в (2.2) определим в узле п ячейки с локальный базис в1, б2, б3, составленный из нормалей 8 к граням/1,/2,/3, образующим трехгранный угол и содержащим рассматриваемый узел (см. фиг. 1б; в данном случае с = с1,/=/1). Аналогичный базис может быть определен и для ячейки с = с2. Дуальным к нему будет нормированный базис б1, 82, 83, определяемый следующим образом:
2 3/1 2 3| 3 1/1 3 1 1 2/1 1 2| ->\
8! = 8 X 8 /к х 8 , в2 = 8 X 8 /к X в , 83 = 8 X 8 /к X в . (2.3)
Заметим, что если ребра, выходящие из узла п, образуют ортогональную систему координат, то векторы б1, 82, 83 будут являться ее базисом.
Пусть дк(п,с), к = 1, 2, 3, — компоненты вектора потока q(n, с) в узле п ячейки с в системе координат, порождаемой дуальным базисом. Тогда имеем
3
Ч (п, с) = X (п, с) 8 к .
(=1
Отсюда, во-первых,
(п,с)2 =[У (п,с)]2 + [х2 (п,с
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.