CC BY
24
8. Берштейн Л.С., Боженюк А.В. Определение нечетких внутренне устойчивых, внешне устойчивых множеств и ядер нечетких ориентированных графов // Известия РАН. ТиСу.
Системы уравнений. Рассмотрим систему уравнений вязкой несжимаемой жидкости [1]
где v = iava - вектор скорости среды (парные строчные греческие буквы традиционно для тензорного исчисления означают суммирование от 1 до размерности пространства n) в момент времени t в точке с декартовыми координатами Xі, 1 < j < n ; i j , v1, 1 < j < n - декартовы единичные орты и
компоненты вектора v вдоль соответствующих координатных направлений; v - кинематический коэффициент вязкости; p - плотность среды (p = const); p - давление; ф - гравитационный потенциал.
Сеточными методами система (1), (2) может быть разрешена с помощью MAC-техники [1]
где ~ - промежуточное, не удовлетворяющее уравнению неразрывности (1) поле скоростей; х - шаг по времени; v и p - удовлетворяющее уравнению неразрывности (1) поле скоростей и давление на новом временном слое соответственно.
Но для адекватного описания гидродинамики мелкого моря нельзя не учитывать плотностную неоднородность (p ф const). Так, например, в эстуариях пресные воды впадающих в моря рек, растекаясь поверх соленых морских вод, могут «забрасывать» примеси на расстояния, заметно превышающие расстояния, получающиеся без учета различия в плотностях. В то же время различие в плотностях таково, что в уравнениях вязкой сжимаемой (p ф const) жидкости [1]
1999. № 1.
В.С. Васильев
ТРЕХМЕРНАЯ СЕТОЧНАЯ МОДЕЛЬ ГИДРОДИНАМИКИ МЕЛКОГО МОРЯ. I
div v = 0,
v't + ia div(vav-v grad va)=-p-1grad p - grad ф,
(l)
(2)
(З)
(4)
pt + div(pv) = 0,
(5)
(pv) + ia div(vapv-|a(grad va+dv/ ara))=- grad (p-|a'div v)-p grad ф, (6)
где ц и ц' - первый и второй динамические коэффициенты вязкости, допустимо ограничиться приближением возмущенной плотности [2-5]
P = Po +8 > H/p0 << 1 ’ (7)
где p0 - некоторое среднее значение плотности.
Это позволяет обойтись модификацией MAC-техники вместо использования методов решения уравнений сжимаемой жидкости. Так, в [2] система из урав-
нения неразрывности сжимаемой жидкости (5), в котором, однако, требуется выполнение и уравнения неразрывности несжимаемой жидкости (1), и уравнений движения несжимаемой жидкости (2) в однородном поле тяжести
grad ф = ^,
где g - ускорение свободного падения, разрешается посредством следующей модификации:
~ = v + х(ia divRe-1grad va -vav)+ Fr-1(l-pc/p)g/g), p = p-xdiv(pv),
x div(p 1 grad ф) = div~, (8)
v = ~ - xp-1 grad Ф, (9)
где Re и Fr - числа Рейнольдса и Фруда; g = |g|; grad Ф = grad p -pcg, z < e
- избыточное давление, обусловленное возмущением; pc = p0 (1 + к (e - z)) >
z < e - гидростатическое распределение плотности (стратификация носит линейный характер и при к> 0 устойчива); e = e(t, x, y) - функция возвышения уровня по отношению к невозмущенной эквипотенциали (в гравитационном поле).
В [3, 4] из системы уравнений (5), (6) с учетом (7) и grad p = grad Ф + pog
выводится
div v = 0, (10)
po (v't + ia div(va v - V grad va )) = - grad Ф + gS. (11)
При малых S можно считать зависимости S от концентрации примеси с или от температуры среды T линейными:
8 = кс с или S = -Kr(T - T),
где T0 - температура, которой соответствует невозмущенная плотность p0.
Тогда из уравнения динамики примеси
c't + div(cv - D grad c)= 0, (12)
где D - коэффициент диффузии рассматриваемой примеси, или из уравнения переноса тепла (в правой части опущено слагаемое с тепловыделением из-за трения, так как для рассматриваемой задачи оно пренебрежимо мало)
cv (T/ + div(Tv)) = acp div(grad T), (13)
где c и c p - удельные теплоемкости; a - коэффициент температуропроводности среды; p0cpa - коэффициент теплопроводности, получается уравнение динамики S
S t + div(vS) = G div(grad S), (14)
где G = D для уравнения (12) или G = acpIcv - для уравнения (13).
Система (10), (11), (14) разрешается посредством следующей модификации: ~ = v + х( ia div(vgrad va - vav)+ p01gs), S = S + xdiv(Ggrad S - vS),
xp-1 div(grad ф) = div~, (15)
v = ~ -xp-1grad Ф. (16)
В [6] рассматривается нелинейная зависимость S и от концентрации, и от
температуры одновременно, однако разрешение системы получающихся уравнений ориентировано на использование переменных «функция тока - вихрь», что не лишено проблем в задачах со свободными поверхностями и при использовании криволинейных сеток [7].
Целью данной работы является построение сеточных аппроксимаций трехмерной модели гидродинамики мелкого моря, обеспечивающей получение решений, удовлетворяющей физическим законам сохранения, в условиях отношения протяженности к глубине 104...105, перепада глубин до 102, колебаний уровня от -100 % до нескольких сотен процентов от глубины.
Исходя из сказанного, сеточные аппроксимации строятся для уравнений сжимаемой жидкости (5), (6), но с ориентацией на разрешение последней какой-либо модификацией MAC-техники.
Криволинейная сетка. Пусть введена система криволинейных координат & Ц , ^). Значения сеточных функций и выражений, относящиеся к узлам, ребрам, граням или ячейкам, будем снабжать целочисленными индексами (i, j, к) соответствующих узлов, полуцелочисленными индексами (i + i, j, к), (i, j + i, к) или (i, j, к +1) соответствующих ребер, граней - (i, j + -1, к +1), (i + -1, j, к +1) или ( ± j + ±к) и ячеек - (i + 1 j + к + ,2) .
На трехмерной криволинейной сетке число внутренних узлов P , число внутренних ребер E , число внутренних граней Hj , число ячеек C, число всех граней H, число всех ребер E и число всех узлов P связаны соотношением
P, <1 Ej <1 Hj < C <1H <1E < P.
С другой стороны, на сетке с достаточно большим числом узлов по каждому направлению
Pj * ± Ej * 3 Hj * C * 3 H * 3 E * P . (17)
Поток вектора v через координатную поверхность £,J = const может быть выражен только через одну-единственную контравариантную компоненту V вектора v. Работа вектора v по координатной линии j-го направления может быть выражена только через одну-единственную ковариантную компоненту Vj вектора v. Если задавать контравариантные компоненты сеточного вектора v на гранях сетки, причем компоненты V1 - на гранях сеточной координатной поверхности
\ = const, v2 - ^ = const, v3 - Q = const, то в выражении уравнения неразрывности для конечного объема - ячейки сетки (шестигранник) можно обойтись парой троек компонент. Соответственно, если задавать ковариантные компоненты сеточного вектора v на ребрах сетки, причем компоненты v1 - на ребрах 1-го координатного направления, v2 - 2-го, v3 - 3-го, то в выражении циркуляции сеточного вектора v по контуру - периметру грани ячейки можно обойтись парой соответствующих грани двоек компонент.
Иное дело - перенос не скалярной (масса), а векторной (импульс) величины (уравнения движения вместо уравнения неразрывности), когда через каждую из граней вида £, = const, ^ = const, Q = const переносятся все три компоненты вектора. Это потребует либо задания на каждой из граней всех трех компонент, либо употребления для каких-то из них сеточных проекторов. Но если компоненты вектора должны удовлетворять какому-то балансовому соотношению, то рассчитывать на то, что оно будет обеспечено проекторами, не приходится. То есть все три компоненты должны определяться в результате решения сеточной системы уравнений. Но если скалярные величины (плотность, давление) задаются во внутренних точках ячеек сетки, то возникает несоответствие (17) размерностей массивов сеточных скаляров и векторов. Употребление сеточных проекторов для скалярных величин не столь проблематично, как для векторных (тем более удовлетворяющих какому-то закону сохранения), но создает неоднородности на границе. Задание же скалярных величин не только на ячейках сетки, но и на других ее компонентах (гранях, ребрах, узлах) потребует балансовых соотношений не только для ячеек, но для конечных объемов, охватывающих эти другие компоненты. А это, по сути дела, введение в рассмотрение новой сетки взамен исходной, с новым набором конечных объемов. Более того, возможна такая интерпретация, что конечным объемом для уравнения баланса импульса будет не только конечный объем, отличный от конечного объема для уравнения баланса массы, но и конечные объемы для уравнений баланса различных компонент импульса будут различны между собой. Здесь не утверждается, что сказанное является непреодолимой проблемой. Напротив, сравнение результатов моделирования при разрешении системы сеточных уравнений в контравариантных, ковариантных или декартовых компонентах составляют предмет отдельной публикации. Здесь же положим следующее. Скалярные величины (плотность, давление, потенциал, коэффициенты вязкости) задаются в ячейках, векторные (скорость) - всеми тремя компонентами в узлах. Тем самым снимается проблема несоответствия размерностей массивов: число узлов на достаточно большой сетке примерно равно числу ячеек. Но тогда утрачивается преимущество контравариант-ных компонент перед декартовыми: выражение сеточной дивергенции вектора для ячейки требует по 24 (по 3 от 8 узлов) и декартовых, и контравариантных, и ковариантных компоненты вместо 6 контравариантных (по 1 от 6 граней) при задании последних на гранях. С другой стороны, разрешение системы в декартовых компонентах имеет и свои преимущества. Во-первых, окончательное представление результатов все равно должно быть в декартовых компонентах, то есть не требуется соответствующих преобразований, что на искривленных сетках сопровождается внесением ошибок вплоть до потери аппроксимации. Во-вторых, разложение по глобальному декартовому базису вместо локального контравариантного (или ковариантного) позволяет «удерживать» консервативность на уровне каждого уравнения в отдельности, а не «размазывает» ее по системе уравнений в целом. Поэтому дальнейшее из-
ложение ведется для случая разложения векторных величин по декартовому базису, хотя все нижеследующие результаты воспроизводимы и в случае контравариантных, и ковариантных компонент, но это также предмет отдельной публикации.
Сопряженность дифференциальных операторов. Аргументами дифференциального оператора div являются векторы, значениями - скаляры, для grad наоборот. Теорема Грина
I f div v dV = | f (v, dS)-|(v, grad f) dV (18)
V S V
выражает сопряженность операторов div и -grad на финитных функциях, равных 0 или 0 на границе области S
(f ,div v)l = (- grad f, v)2,
где V - ограниченная пространственно-односвязная область; S - замкнутая регулярная граница этой области; f и v - однозначные и непрерывные вместе со всеми частными производными, встречающимися в объемных или поверхностных интегралах, в V и на S скалярная и векторная функции.
В основе скалярного произведения в пространстве скаляров лежит
обычное арифметическое умножение, а в основе (...,...)2 - скалярное произведение векторов в обычном евклидовом смысле.
Пусть сеточная аппроксимация пространственных производных на ячейке (например, слагаемых div) выражается только через значения дифференцируемого выражения в узлах, смежных ячейке (в угловые скобки будем заключать аппроксимируемые выражения)
=1 f с(g---) (А + с(g--:) (А +
г+1 і+1,k+1 4 f г+1,j+^2,k+|\J/г,j,k г+2 j+1,к+j\f/i,j,k+1
+
+ с(g-+-) (А + C(g-++) (А
+ Ci + 2, j + 2,к+V 'i, j+1,к + Ci+i j+^,к+f /i, j+1,к+1
+с(?1Дк+J+и,к+с^+1к+|(А+и,к+1 +
+с+с+и+j .1,., .1,.+с+с+и j ^+и+1 ) • <19>
где g - одна из координат x, у или z; для производной/’х используются коэффициен-
ты C(х±±±), дляf’y - C(у±±±), для f'z - C(z±±±); А - обозначение конечного объема Jh^h^h^; J = D(x,y,z)/D(^,ц,^) - якобиан преобразования координат (xy,z) ^(#,Ц, ^); для сеточной функции, заданной в узлах, f) . к = f-у-,к .
Тогда следующая сопряженная сеточная аппроксимация пространственных производных в узлах (например, компонент grad)
(fg ^ i, j,k =- 4 ^ f). - |, j - ±,к - 1 + A - 1,j - ±,к+* +
+ c (g+-+) (A + c (g+--) (A +
i - 2 j+i к - i\ J li - jj+± к - 1 i - 2,j+1 к+i\J li -1 j+1 к+1
+ с(g-++) (A + с(g-+-) (A +
+ Ci+1 j - 1,к - ЛJ /i+1, j - 1,к -1 + Ci+1 j -1 к+ЛJ / i+1, j -1 к+i +
+ І.Д.*-+ІІ+2-к-2 + ^.-/+)2.к + 2^А+І/+2-к+1 )
(20)
обеспечивает выполнение на сеточном уровне теоремы Грина (18). На неограниченной сетке или для сеточных функций, в граничных узлах (г, ], к), равных / ■ к = 0 (для наглядности, чтобы не загромождать выражение суммированием по граничным элементам)
Х , ч(^^) і+1, і+і,к+1 ^г'+2-і+2-к+2
к+2)
-X Лм( < д) і
(і-І,к)
І,к '
где суммирование ^ означает суммирование по ячейкам, а ^ - по узлам.
(г+1,]'+1,к+1) (г,]',к)
Учитывая, что = Б(/, у, г)/ £>(^, С), Гу J = б(х /, г)/ Д(^, г|, с),
/[3 = Б(х, у, /)/Б(^, ^, С), а Б(у, 9, %)/Б(^, ^, С) - антисимметричная полилинейная форма производных от у , 9 их, достаточно общим представлением коэффициентов С(я±±±) можно считать
г—) = -V---) - V(г|г---) - V(<Сг—)
.1^1^! +1, /+1,к+1 1 г-4.1 °,-4.1 ,-4.1 г-4-1,
+ІІ+^,к+
С С С С
С(я+--)
і+і І+і,к+і
я—+) _ —+) _ с>(ля—+) , с>(Сг—+)
+іі+і к+і і+2,І+і к+і і+і і+ік+1 і+^ і+і,к+1!
я-+-) — _еЙг-+-) _і_ с(ля-+-) _ с>(Сг-+-)
+1 і+|, к+| і+|, /+1 к+| і+2-і+^,к+1 і+і і+ік+1 !
я-++) — _еЙг-++) _і_ с(ля-++) , с>(Сг-++)
+ і І+|, к+2 1+ІІ+2’к+і + *+2-і+|,к+2 + '+1-і+|,к+і !
_ оЙг+--) _ с>(ля+--) _ о(Сг+--)
г+2,]+1,к+2 г'+!■ 3+|.к+1 г+2,1+|.к+2 г'+Ы+1-к+2 ’
+-+) _ оЙг+-+) _ +-+) , е(сг+-+)
г+1]+1,к+2 г+1,3+1,к+1 г+1^.,-+1,к+2 г+]+^,к+^ ’
С(%++-) - ++-) _|_ ++-) _ е(сг++-)
г++2,к+1 г+1,]+1,к+1 г+1,]+1,к+1 г+1,]+1,к+^ ’
/-'(я+++) - сЙя+++) _|_ с(ля+++) , е(Сг+++)
г+^,]+^,к+ ^ г+^,]+^,к+^ г+^,]+^к+^ г+2,]+^,к+^ '
Если относящиеся к ячейкам величины 5,(^я±±±,), , ,
^ г+2,] + ^,к +1 ’
±±±2^ по сути дела относятся к граням ячеек, то есть
С>(^я+±±) _ С>(^я-±±) _ с>(^я*±±) с>(ля±+±) _ с>(ля±-±) _ с>(ля±*±)
г - 2,]+^,к+2 г+^, ]+^,к+1 г, ]+^,к+^ ’ г+1^.,' - 1,к+ 2 г+^^./-+^,к+2 г+2’ 3*+^ ;
(21)
с>(г|я ±±±) і+2-і+і,к+1 ’
о(с?±±+) _ с’(Ся±±-) _ о(Сг±±ф)
^. і і, і — О. і і, і — О. і її,
1+2-і + і к-1 1++ ік+2 1++|,к’
(22)
то для сеточных функций, в граничных узлах (г, ], к), равных нулю (/ , к = 0), имеет место
2
2
X (Ы = 0. (23)
(+/,к+ 1Г * 1+12,3+^+1
Обозначим через №, №), /(^г), >х), /(лу), /(лг), /(Сх), /(Су), /(Сг) алгебраические дополнения элементов х^, у^, г'^, х^ , у^ , г^ , хС , уС и г’^ соответственно в якобиане / = Б(х, у, г)/ Б(£„ ^ с) .
Если относящиеся к граням величины (22) аппроксимируют соответствующие алгебраические дополнения
j+!,к+1 i+i,},к+i \ ^ 4 i+i,it+1
х:г2±я,к =(j(ewtrj+,к - ^
то, во-первых, (19) и (20) действительно аппроксимируют соответствующие пространственные производные, а, во-вторых, на сеточном уровне выполняется теорема Гаус-са-Остроградского (являющаяся частным случаем теоремы Грина (18) при f = const):
j div v dV = f(v, dS), (25)
V S
в том числе и по отдельным слагаемым оператора div (23).
Таким образом, аппроксимация пространственных производных (19) действительно может быть принята в качестве аппроксимации слагаемых оператора div. При этом сопряженная аппроксимация для компонент grad (20) в корректор-этапах MAC-техники (4), (9), (16)
v j = ~1,м -ТР-1 (grai Р А,м/ (А) j •
где (grad Р Аj = КРХА)i,hk + jР'у А. j,к + Чp'zАi,Лк ;
i • j • k - единичные декартовы орты,
обеспечивает на сеточном уровне положительную определенность операторов этапов расчета поправки к давлению MAC-техники (3), (8), (15). Действительно, в
, I .
(+и+1к+i)
,I ,.(div”А>i+i,i+!,к+4pi+w+!,к+i • <26>
i ixi t-li I
где (divvАУ 1 . , 1 = (u'xA) 1 . . 1 +lv'yA), + (wZAV 1 . , 1 •
\ h+^-,j+|,к+j \ x h+j,j+|,к + 1 \ y li+1,j+1,к+1 ' z h+j,i + ^,к+
может быть выделена (за исключением суммирования по граничным элементам) знакоопределенная сумма
I (W), ДрУ^2,. +<pZa>2Лк1/(4,м = I |(gradра>,,м|2/!A>,.
(i,M)V ,i, 1,М ^, " (i,M) ,Л ' ,i,
БИБЛИОГРАФИЧЕСКИИ СПИСОК
1. Флетчер К. Вычислительные методы в динамике жидкостей: В 2-х т/: Пер. с англ.- М.: Мир, 1991.
2. Белоцерковский С.А., Гущин В.А., Коньшин А.Г. Метод расщепления для исследования течений стратифицированной жидкости со свободной поверхностью // Ж. вычисл. ма-тем. и матем. физики. 1987. Т. 27. № 4. С. 594-609.
3. Лойцянский Л.Г. Механика жидкости и газа. Изд. 5-е, перераб.- М.: Гл. ред. физ.-мат. лит. изд-ва «Наука», 1978.- 736 с.
4. ИевлевВ.М. Численное моделирование турбулентных течений.- М.: Наука, 1990.- 216 с.
5. Габов С.А., Свешников А.Г. Задачи динамики стратифицированных жидкостей.- М.: Наука. Гл. ред. физ.-мат. лит., 1986.- 288 с.
6. Марчук Г.И., Кочергин В.П. Саркисян А.С. и др. Математические модели циркуляции в океане.- Новосибирск, 1980.- 285 с.
7. Васильев В.С. Динамические переменные на криволинейных сетках // Известия ТРТУ. 2004. № 4.- С. 191-201.
В.С. Васильев
ТРЕХМЕРНАЯ СЕТОЧНАЯ МОДЕЛЬ ГИДРОДИНАМИКИ МЕЛКОГО МОРЯ. II
Сеточные источники. Одной из простейших сеточных аппроксимаций, удовлетворяющих изложенному в [1], будет такая, у которой относящиеся к граням величины (22) [1] относятся к граням в целом, не различая углов граней,
V Йг •—) = V *-+) = V *+-) = V (£г *++) = V ) = / Т )Ь Ь \
и.+к+1 и] +1 к+| г J+2к+1 г,/+2к+| г,/+1к+| \ Л С/.,.+х к+х ,
V(лг-•-) = V(лг-•+) = V(лг+*-) = V(лг+,+) = V(лг) = / Т(лг)И к\
г+1 3,к+| г+^3к+1 ;■+ 1;,к+-2 .+1., к+1 .+1 .,к+1 \ £ С/.+х . к+х ’
2 >./! 2
V(сг—'•) = V(сг-+,) = V(Сг+-,) = V(Сг++,) = V(Сг) = / Т(Сг)й ^ (1)
^+1 /+2, к = Ч+2,/+1 к = ^+± /+± к = ^+1,/'+^2,к = ^+2. +1 к = \Т .+х .+X к ‘ (1)
2 2 ’
Если ребра криволинейной сетки [1] - прямолинейные отрезки, то следующие аппроксимации алгебраических дополнений [1] - выражения точные (здесь Ф = г, ст = /, т = к - целые)
(в(/, 1=((«-/сглХд),^ 1 =
=і (і -( ли Иыш ч^,5) )-
2 \\У /ф,у+1,к+1 V /ф,у,кЛ\5/ф,у, к+1 V5/ф,у+1,Ат/
-(л(5) /л(5) >/*\(5) /„\(5) ))
V /ф,у,к+1 V /ф,у+1,кЛ\5/ф,у+1,к+1 \5/ф,у,к//’
(в(/, * > же, фа) ,+ ^ 1 =((/&- /кЫ ,+ ^ 1 =
=1 ((/)(,+).о.* *1 -(/)(",«)(^(,+1,о.* -(* )!",>.* *1 >-
-( А(л) ^ Л>/*\М ^*\))
V /;+1,а,к V /;,а,к+1Л\5А'+1,а,к+1 V5/;,ст,к/Г
{D(f, gVD(5, пМ),+1, у+1,х ={(Т\*л - /л *5)^5йл>'+Ху, 1,т =
= ) (((/)(+)1,.7-+1,т -(45,1 *)(^>+1,т Ч^(?,.,-,х)-
- (/)(5)+1,т ^/)(+1,л,>(^(+)1,.,-+1,т Ч*)(3,т))’ (2)
