Научная статья на тему 'Трехмерная сеточная модель гидродинамики мелкого моря. I'

Трехмерная сеточная модель гидродинамики мелкого моря. I Текст научной статьи по специальности «Физика»

CC BY
175
26
i Надоели баннеры? Вы всегда можете отключить рекламу.
i Надоели баннеры? Вы всегда можете отключить рекламу.
iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

Текст научной работы на тему «Трехмерная сеточная модель гидродинамики мелкого моря. I»

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,

+ІІ+^,к+

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

С С С С

С(я+--)

і+і І+і,к+і

я—+) _ —+) _ с>(ля—+) , с>(Сг—+)

+іі+і к+і і+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)

i Надоели баннеры? Вы всегда можете отключить рекламу.