Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками

Схема центральных разностей, часто применяемая для численной аппроксимации адвективных слагаемых уравнений движения в моделях геофизической гидродинамики с разнесенными сетками, может является причиной возникновения нефизических осцилляций поля скорости и приводить к неустойчивости. В работе предлож...

Full description

Saved in:

Bibliographic Details
Date:	2008
Main Author:	Нестеров, А.А.
Format:	Article
Language:	Russian
Published:	Інститут гідромеханіки НАН України 2008
Online Access:	https://nasplib.isofts.kiev.ua/handle/123456789/4631
Tags:	Add Tag No Tags, Be the first to tag this record!
Journal Title:	Digital Library of Periodicals of National Academy of Sciences of Ukraine
Cite this:	Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками / А.А. Нестеров // Прикладна гідромеханіка. — 2008. — № 1. — С. 58-68. — Бібліогр.: 10 назв. — рос.

Institution

Digital Library of Periodicals of National Academy of Sciences of Ukraine

_version_	1860214210477162496
author	Нестеров, А.А.
author_facet	Нестеров, А.А.
citation_txt	Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками / А.А. Нестеров // Прикладна гідромеханіка. — 2008. — № 1. — С. 58-68. — Бібліогр.: 10 назв. — рос.
collection	DSpace DC
description	Схема центральных разностей, часто применяемая для численной аппроксимации адвективных слагаемых уравнений движения в моделях геофизической гидродинамики с разнесенными сетками, может является причиной возникновения нефизических осцилляций поля скорости и приводить к неустойчивости. В работе предложен метод расчeта адвективных слагаемых уравнений движения на основе специальной интерполяции поля скорости, с последующим применением какой-либо из схем уравнений переноса, включая TVD схемы. Метод рассматривается в применении к задачам моделирования течений в водных резервуарах со свободной поверхностью. Схема центральних рiзниць, яка часто використовується для чисельної апроксимацiї адвекцiйних доданкiв рiвнянь руху в моделях геофiзичної гiдродинамiки з рознесеними сiтками, може бути причиною виникнення нефiзичних осциляцiй поля швидкостi i призводити до нестiйкостi. В роботi запропоновано метод розрахунку адвекцiйних компонент рiвнянь руху на основi спецiальної iнтерполяцiї поля швидкостi, з наступним використанням будь-якої зi схем рiвнянь переносу, включаючи TVD схеми. Метод розглядається в застосуваннi до задач моделювання течiй у водних резервуарах з вiльною поверхнею. Scheme of central differences, which is often used for numerical approximation of the advection terms in the momentum equations in geophysical hydrodynamics models with staggered grids, can induce non-physical oscillations of the velocity field and cause numerical instabilities. This paper proposes a numerical method for momentum advection approximation, which is based on the specific interpolation of the velocity field with subsequent application of any scheme for a scalar transport equation, including TVD schemes. This method is studied in the application to the problems of flow modeling in the water reservoirs with free surface.
first_indexed	2025-12-07T18:15:43Z
format	Article
fulltext	ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 УДК 532.465 ПРИМЕНЕНИЕ TVD–СХЕМ ДЛЯ АППРОКСИМАЦИИ ПЕРЕНОСА ИМПУЛЬСА В МОДЕЛЯХ ГЕОФИЗИЧЕСКОЙ ГИДРОДИНАМИКИ С РАЗНЕСЕННЫМИ СЕТКАМИ А. Н Е СТЕР ОВ Tropical Marine Science Institute, National University of Singapore Получено15.05.2007 Схема центральных разностей, часто применяемая для численной аппроксимации адвективных слагаемых урав- нений движения в моделях геофизической гидродинамики с разнесенными сетками, может является причиной во- зникновения нефизических осцилляций поля скорости и приводить к неустойчивости. В работе предложен метод расчeта адвективных слагаемых уравнений движения на основе специальной интерполяции поля скорости, с по- следующим применением какой-либо из схем уравнений переноса, включая TVD схемы. Метод рассматривается в применении к задачам моделирования течений в водных резервуарах со свободной поверхностью. Схема центральних рiзниць, яка часто використовується для чисельної апроксимацiї адвекцiйних доданкiв рiвнянь руху в моделях геофiзичної гiдродинамiки з рознесеними сiтками, може бути причиною виникнення нефiзичних осциляцiй поля швидкостi i призводити до нестiйкостi. В роботi запропоновано метод розрахунку адвекцiйних ком- понент рiвнянь руху на основi спецiальної iнтерполяцiї поля швидкостi, з наступним використанням будь-якої зi схем рiвнянь переносу, включаючи TVD схеми. Метод розглядається в застосуваннi до задач моделювання течiй у водних резервуарах з вiльною поверхнею. Scheme of central differences, which is often used for numerical approximation of the advection terms in the momentum equations in geophysical hydrodynamics models with staggered grids, can induce non-physical oscillations of the velocity field and cause numerical instabilities. This paper proposes a numerical method for momentum advection approximation, which is based on the specific interpolation of the velocity field with subsequent application of any scheme for a scalar transport equation, including TVD schemes. This method is studied in the application to the problems of flow modeling in the water reservoirs with free surface. ВВЕДЕНИЕ В задачах геофизической гидродинамики основ- ной вклад в уравнения движения обычно вносят слагаемые, обусловленные гравитационной силой [1]. С целью упрощения их численной аппрокси- мации часто используются, так называемые, ра- знесенные сетки [2]. Средние значения скалярных переменных состояния, таких как температура, со- лeность, плотность, полагаются при этом в центре численной ячейки, а значения компонент вектора скорости – на соответствующих гранях. Дополни- тельным преимуществом применения разнесенных сеток является то, что дискретные консерватив- ные формулировки для уравнений переноса скаля- ров также упрощаются. Проблеме численной ап- проксимации адвективных слагаемых уравнения переноса скаляра уделяется особое внимание, так как к свойствам скалярных переменных, как пра- вило, предъявляются повышенные требования по сравнению со свойствами компонент скорости, в частности, монотонность и позитивность. TVD- схемы [2–4] обладают такими свойствами. В то же время, во многих моделях гидродинамики резер- вуаров со свободной поверхностью (например [5] или [6]) для аппроксимации адвективных слага- емых в уравнениях импульса используется схема центральных разностей, как самая простая консе- рвативная схема 2-го порядка. Непосредственное применение тех же самых схем адвекции, что и для скаляра, не представляется возможным из-за использования разнесенных сеток. Однако ещe в [2] упоминается, что в случае ма- лой вязкости в задачах геофизической гидродина- мики могут возникать проблемы с устойчивостью, такие же, как при решении уравнения Бюргерса. В задачах моделирования течений в резервуарах со свободной поверхностью инерционные слагаемые уравнений движения могут давать существенный вклад по сравнению с баротропными слагаемыми в областях, где скорость или направление пото- ка резко изменяется по тем или иным причинам. Типичными случаями являются: портовые зоны с резкими изменениями береговой линии; каналы и реки с дамбами и порогами; прибрежные регионы морей с большими градиентами глубин; области осушения/затопления. Например, при моделиро- вании поля скорости в канале с поперечной дамбой с помощью модели [5] вблизи кромки дамбы во- зникали нефизические осцилляции поля скорости, выделенные на pис. 1 жирными стрелками. Для изучения влияния инерционных слагаемых гори- зонтальная вязкость в этом примере полагалась равной нулю. Поле скорости моделировалось как 58 c© А. Нестеров, 2008 ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 Рис. 1. Осцилляции поля скорости вблизи острой дамбы: расчeт без горизонтальной вязкости трехмерное; показанное поле – усредненное по глу- бине. Канал длиной 10 км, шириной 2.1 км и глу- биной 10 м покрывался сеткой с горизонтальным разрешением 100 м. У берегов задавались условия скольжения. На двух открытых границах задава- лись уровни 0.1 и −0.1 м. То, что проблема имеет отношение именно к адвективным слагаемым, сле- дует из того, что: пренебрежение адвективными слагаемыми приводит к отсутствию осцилляций (pис. 2); неустойчивость возникает в месте резко- го изменения потока; добавление горизонтальной вязкости уменьшает или предотвращает осцилля- ции; возникновение осцилляций не зависит от ша- га по времени и метода расчeта уровня поверхно- сти (явного [6] или неявного [7]). На pис. 3 пока- зано поле скорости при наличии горизонтальной вязкости (соответственно, с условиями прилипа- ния на берегах) с коэффициентом, рассчитывае- мым по формуле Смагоринского, как в [6]. Одна- ко использование вязкости для подавления осцил- Рис. 2. Поле скорости вблизи острой дамбы: расчeт без инерционных слагаемых ляций является сомнительным, так как не гаран- тирует их отсутствия и может влиять на резуль- тат во всей расчeтной области. В данном примере увеличение вязкости приводит к уменьшению ин- тенсивности вихря позади дамбы, выделенного на рисунках. Очевидно, что схема центральных разностей становится непригодной в вышеупомянутых слу- чаях. Для решения этой проблемы в [8] предла- галась схема ван Леера [9] для аппроксимации одномерных и двумерных уравнений мелкой во- ды. В данной работе предлагается относительно простой метод, позволяющий применить любую из схем уравнений переноса скаляра в разнесенных сетках для аппроксимации адвективных слагае- мых трёхмерных уравнений импульса. При этом возможно использование тех же самыx подпро- грамм расчeта адвекции, что и для скалярных пе- ременных. А. Нестеров 59 ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 Рис. 3. Поле скорости вблизи острой дамбы: расчeт при наличии горизонтальной вязкости 1. ЧИСЛЕННАЯ АППРОКСИМАЦИЯ АДВЕ- КТИВНЫХ СЛАГАЕМЫХ В УРАВНЕНИЯХ ИМПУЛЬСА 1.1. Уравнения модели Уравнения, рассматриваемые в данной работе, описывают гидродинамику трeхмерных водоeмов со свободной поверхностью в гидростатическом приближении. Уравнения записываются в декар- товых координатах x, y по горизонтали и в σ- системе координат по вертикали [6]. Последняя связана с декартовой z-координатой соотношени- ем z − η = σD, где η = η(x, y, t) – уровень поверх- ности; D = H + η – полная глубина; H = H(x, y) – глубина при невозмущённой поверхности. Рассма- триваемые уравнения для горизонтальных компо- нент вектора скорости U и V имеют вид: ∂DU ∂t + ∂DU2 ∂x + ∂DUV ∂y + ∂Uω ∂σ = = −gD ∂η ∂x + fDV + ∂ ∂σ ( ν D ∂U ∂σ ) + Fx, (1) ∂DV ∂t + ∂DUV ∂x + ∂DV 2 ∂y + ∂V ω ∂σ = = −gD ∂η ∂y − fDU + ∂ ∂σ ( ν D ∂V ∂σ ) + Fy. (2) “Вертикальная” (см. [6]) компонента скорости ω в гидростатическом приближении рассчитывается из уравнения неразрывности, которое в σ-системе координат принимает вид: ∂DU ∂x + ∂DV ∂y + ∂ω ∂σ + ∂η ∂t = 0, (3) а уровень поверхности η находится из этого же уравнения, усреднeнного по вертикали. Здесь t – время; g – гравитационное ускорение; f – параметр Кориолиса; ν – коэффициент вертикальной вяз- кости; Fx, Fy – горизонтальные компоненты сил вязкости, которые записываются так же, как в [6]. Далее рассматривается вопрос о применении чис- ленных схем аппроксимации уравнения адвектив- ного переноса скаляра C: ∂DC ∂t + ∂DUC ∂x + ∂DV C ∂y + ∂ωC ∂σ = 0, (4) в разнесенных сетках к адвективным слагаемым уравнений (1), (2). Метод также может быть использован для аппроксимации адвективных сла- гаемых для третьей компоненты скорости в неги- дростатическом случае. 1.2. Cхема центральных разностей Будем рассматривать постоянный в σ-системе раз- мер численной ячейки ∆σ = const, так же как и ∆x = const, ∆y = const. Соответствующий вер- тикальный размер численной ячейки в z-системе зависит от горизонтальных координат и времени: ∆z = ∆z(x, y, t), так как ∆z = D∆σ. Обозначим Un i,j,k компоненту вектора скорости на n-м шаге по времени на грани пространственной ячейки (i, j, k) (pис. 4). Рассматриваемое адвективное слагаемое уравнения (1) состоит из трех компонент, соответ- ствующих переносу в каждом из трех пространс- твенных направлений: ADVU (U)= ∂ ∂x ( DU2 ) , ADVV (U)= ∂ ∂y (DUV ) , ADVω (U) = ∂ ∂σ (Uω) . В явной схеме центральных разностей, которая используется, например, в моделях [5, 6], они ап- проксимируются как: 60 А. Нестеров ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 −ADVU (U)∆x ≈ ≈ ∆zn i− 3 2 ,jU n i−1,j,k + ∆zn i− 1 2 ,jU n i,j,k 2 Ui−1,j,k + Ui,j,k 2 − − ∆zn i−1 2 ,jU n i,j,k + ∆zn i+ 1 2 ,jU n i+1,j,k 2 Ui,j,k + Ui+1,j,k 2 , −ADVV (U)∆y ≈ ≈ ∆zn i−1,j−1 2 V n i−1,j,k + ∆zn i,j−1 2 V n i,j,k 2 × × Ui,j−1,k + Ui,j,k 2 − − ∆zn i−1,j−1 2 V n i−1,j+1,k + ∆zn i,j+ 1 2 V n i,j+1,k 2 × × Ui,j,k + Ui,j+1,k 2 , −ADVω (U)∆σ ≈ ≈ ωi−1,j,k+1 + ωi,j,k+1 2 Ui,j,k + Ui,j,k+1 2 − − ωi−1,j,k + ωi,j,k 2 Ui,j,k−1 + Ui,j,k 2 , где ∆z i− 1 2 , j = ∆zn i−1,j + ∆zn i,j 2 , ∆z i,j− 1 2 = ∆zn i,j−1 + ∆zn i,j 2 . Индекс “k” возрастает здесь в направлении “вниз” (см. pис. 4). Отметим, что каждая из этих формулировок представляет из себя разность вте- кающего и вытекающего потоков и поэтому явля- ется консервативной. 1.3. Краткий обзор схем для адвекции скаляра Для численной аппроксимации уравнения адве- ктивного переноса скаляра (4) в разнесенных се- тках существует множество схем. Следуя работе [10], в одномерном случае они могут быть пред- ставлены как Cn+1 i Dn+1 i = Cn i Dn i + +ULDL ∆t ∆x CL − URDR ∆t ∆x CR, где UL, UR – значения скорости на левой и правой гранях ячейки; DL, CL, DR, CR – глубина и зна- чения скаляра на соответствующих гранях. Зада- ча построения численной схемы адвекции как раз и состоит в определении значений CL, CR. На- пример, в схеме центральных разностей они пред- ставляются как средние арифметические: CR = 1 2 ( Cn i + Cn i+1 ) . Для рассмотренных далее TVD-схем может быть использовано следующее представление: CR = { C+ R (UR ≥ 0), C− R (UR < 0), а C+ R и C− R определяются из уравнений C+ R = Cn i + 1 2 Ψ(r+ i )(1 − c) ( Cn i+1 − Cn i ) , C− R = Cn i+1 − 1 2 Ψ(r−i )(1 + c) ( Cn i+1 − Cn i ) , где c – локальное число Куранта со знаком; ri является отношением градиента скаляра вверх по потоку к локальному градиенту: c = UR ∆t ∆x , r+ i = Cn i − Cn i−1 Cn i+1 − Cn i , r−i = Cn i+2 − Cn i+1 Cn i+1 − Cn i . В данной работе рассмотрены ограничители Ψ из [10], включающие не только TVD-схемы: Ψ(r, c) =                                                                  0 Upwind, 1 Lax-Wendroff max[min(2r, 1), min(r, 2), 0] Superbee, (r + \|r\|)/(1 + r) van Leer, min(r, 1) MINMOD, max [ min ( 2, 2r, 1 + r 2 ) , 0 ] MUSCL, max [ min ( 3(1 + r) + (1 − r)(1 − 2\|c\|) 6 , ULTIMATE-, 2 1 − \|c\| , 2r \|c\| ) , 0 ] QUICKEST,        min (2r/\|c\|, 1) , 0 ≤ r ≤ 1 min ( r, 2 1 − \|c\| ) , r > 1 0, r < 0 Super-C. В описанных выше схемах параметры UL и UR известны. В случае применения к компонентам скорости, когда роль скаляра C играет U или V , эти значения должны быть проинтерполированы. А. Нестеров 61 ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 Рис. 4. Ячейка для U -компоненты скорости 1.4. Схема интерполяции поля скорости Для обеспечения монотонности любой из схем, описанных в предыдущем разделе, дискретное по- ле скорости, которое используется в качестве вхо- дного параметра, должно удовлетворять дискре- тному аналогу уравнения неразрывности (3). В применении к “ячейкам скаляра” Ci−1,j,k и Ci,j,k последнее аппроксимируется как: Un i−1,j,k ∆zn i−2,j + ∆zn i−1,j 2 ∆y− −Un i,j,k ∆zn i−1,j + ∆zn i,j 2 ∆y+ +V n i−1,j,k ∆zn i−1,j + ∆zn i−1,j−1 2 ∆x− −V n i−1,j+1,k ∆zn i−1,j+1 + ∆zn i−1,j 2 ∆x+ +ωn i−1,j,k+1∆x∆y − ωn i−1,j,k∆x∆y+ + ηn+1 i−1,j − ηn i−1,j ∆t = 0, (5) Un i,j,k ∆zn i−1,j + ∆zn i,j 2 ∆y− −Un i+1,j ∆zn i,j + ∆zn i+1,j 2 ∆y+ +V n i,j,k ∆zn i,j + ∆zn i,j−1 2 ∆x− −V n i,j+1,k ∆zn i,j+1 + ∆zn i,j 2 ∆x+ +ωn i,j,k+1∆x∆y − ωn i,j,k∆x∆y+ + ηn+1 i,j − ηn i,j ∆t = 0, (6) Аналогично, для того, чтобы применить какую- либо из приведeнных в предыдущем разделе схем к адвективным слагаемым уравнений импульса, трeхмерное поле скорости и поле уровня поверх- ности, используемые в качестве параметров, необ- ходимо проинтерполировать таким образом, что- бы полученные поля удовлетворяли дискретно- му уравнению неразрывности (3) в применении к соответствующей компоненте скорости. Желае- мым условием при этом является то, чтобы проин- терполированные поля были аппроксимацией 2-го или более высокого порядка исходного поля. Опре- делим геометрическое положение “ячейки компо- ненты скорости” для U так, как показано на pис. 4. В данном случае эта ячейка оказывается такой же призмой, как и ячейка для скаляра, с прямоу- гольным основанием, смещeнной в направлении, соответствующем компоненте скорости. На pис. 4 используются следующие обозначения: Ci,j,k – скаляр (в центре прозрачной ячейки); Ui,j,k – ком- понента вектора скорости, заданного в центре гра- ни между “ячейками скаляра” Ci−1,j,k и Ci,j,k. Компоненты скорости V и ω задаются аналогично. Индекс k возрастает здесь в направлении “вниз”. Компоненты скорости из выражений (5) и (6) входят в качестве параметров в схему дискрети- зации адвекции скаляра, описанных в предыду- щем разделе. Обозначим теперь Un i+ 1 2 ,j,k , Un i− 1 2 ,j,k , V n i− 1 2 ,j,k , V n i+ 1 2 ,j,k , Wn i− 1 2 ,j,k , Wn i+ 1 2 ,j,k – компонен- ты скорости на гранях “Un i,j,k-ячейки” (pис. 5). Известные компоненты скорости, расположенные в центрах граней “скалярных ячеек”, показаны пунктирными стрелками; неизвестные, подлежа- щие определению, – сплошными. Для выполнения дискретного уравнения нера- зрывности, примененного к “Ui,j,k-ячейке”, необхо- димо выполнение следующего условия: 62 А. Нестеров ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 Рис. 5. Ui,j,k ячейка с известными и подлежащими определению компонентами скорости Un i− 1 2 ,j,k ∆zn i−1,j∆y − Un i+ 1 2 ,j,k ∆zn i,j∆y+ +V n i− 1 2 ,j,k × × ∆zn i−1,j + ∆zn i,j + ∆zn i−1,j−1 + ∆zn i,j−1 4 ∆x− −V n i− 1 2 ,j+1,k × × ∆zn i−1,j+1 + ∆zn i,j+1 + ∆zn i−1,j + ∆zn i,j 4 ∆x− −ωn i− 1 2 ,j,k ∆x∆y − ωn i− 1 2 ,j,k+1 ∆x∆y+ + 1 2 ( ηn+1 i−1,j + ηn+1 i,j ) − 1 2 ( ηn i−1,j + ηn i,j ) ∆t = 0. (7) Легко видеть, что если сложить выражения (5) и (6) и затем определить компоненты скорости как Un i− 1 2 ,j,k = Ui−1,j,k∆zn i− 3 2 ,j + Ui,j,k∆zn i−1 2 ,j 2∆zi−1,j , V n i− 1 2 ,j,k = Vi−1,j,k∆zn i−1,j−1 2 + Vi,j,k∆zn i,j−1 2 ∆zn i−1,j−1 2 + ∆zn i,j−1 2 , ωn i− 1 2 ,j,k = ωi−1,j,k + ωi,j,k 2 , (8) где ∆z i− 1 2 , j = ∆zn i−1,j+ ∆zn i,j 2 , ∆zi,j−1 2 = ∆zn i,j−1+ ∆zn i,j 2 , то уравнение (7) выполняется тождественно. Та- кое определение не является единственно возмо- жным, однако далее будет доказано, что такое Рис. 6. Интерполяция поля глубины для U представление имеет 2-й порядок точности аппро- ксимации в случае, если поля уровня и скорости имеют непрерывные производные. Схема дополняется проинтерполированным по- лем уровня поверхности: DUC(i,j) = Di−1,j + Di,j 2 , DUL(i,j) = Di−1,j, DUR(i,j) = Di,j, DUD(i,j) = Di−1,j + Di,j + Di−1,j−1 + Di,j−1 4 , DUU(i,j) = Di−1,j+1 + Di,j+1 + Di−1,j + Di,j 4 , где DUC(i,j) – полная глубина в центре основания “U -ячейки”; DUL(i,j), DUR(i,j), DUD(i,j), DUU(i,j) – глубины в центрах соответственно левого, пра- вого, нижнего и верхнего рeбер (pис. 6). Аналогичным образом, с точностью до пере- становки индексов i и j, определяются дискре- тные поля скорости и уровня для аппроксима- ции адвективных слагаемых уравнения (2) для V - компоненты скорости. При рассмотрении негидро- статических задач выделение “ωi,j,k-ячейки”, как состоящей из верхней половины “Ci,j,k-ячейки” и нижней половины “Ci,j,k−1-ячейки”, приводит к тому, что компоненты вектора скорости на еe гра- нях определяются как средние арифметические: Un i,j,k−1 2 = Un i,j,k−1 + Un i,j,k 2 , V n i,j,k−1 2 = V n i,j,k−1 + V n i,j,k 2 , ωn i,j,k−1 2 = ωn i,j,k−1 + ωn i,j,k 2 . А. Нестеров 63 ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 Рис. 7. Поле скорости вблизи острой дамбы, расчeт по предложенному методу со схемой MUSCL 1.5. Теоретический порядок аппроксимации Покажем, что в случае гладких полей скорости и уровня, предложенная схема имеет 2-й порядок аппроксимации. Перепишем уравнение для U из [8], опуская индексы k и n, в виде Ui− 1 2 ,j = Ui−1,j + Ui,j 4 + + Ui−1,j ∆zi−2,j ∆zi−1,j + Ui,j ∆zi,j ∆zi−1,j 4 . Будем рассматривать высоту численной ячейки ∆z как гладкую функцию ∆z = ∆z(x, ·). Тогда ∆zi−2,j = ∆zi−1,j − ∆z′i−1,j∆x + O ( ∆x2 ) , ∆zi,j = ∆zi−1,j + ∆z′i−1,j∆x + O ( ∆x2 ) . Подставим эти выражения в предыдущее: Рис. 8. Батиметрия для квазиодномерного теста Ui− 1 2 ,j = 1 4 ∆z′i−1,j ∆zi−1,j ∆x (Ui,j − Ui−1,j) + + Ui,j + Ui−1,j 2 + O ( ∆x2 ) . Предположим, что U является гладкой функци- ей U = U(x, ·). Тогда Ui,j − Ui−1,j = ∆xU ′ i− 1 2 ,j + O ( ∆x2 ) . Подставляя в предыдущее, получим Ui− 1 2 ,j = Ui,j + Ui−1,j 2 + O ( ∆x2 ) , что является аппроксимацией 2-го порядка то- чности U -компоненты скорости посередине между двумя соседними узлами, в которых задаются зна- чения Vi−1,j,k и Vi,j,k как среднее арифметическое (см. [2]). Аналогично, уравнение для Vi− 1 2 ,j,k из (8) при- водится к виду Vi− 1 2 ,j = Vi−1,j + Vi,j 2 + + 1 2 ∆zi,j + ∆zi,j−1 − ∆zi−1,j − ∆zi−1,j−1 ∆zi,j + ∆zi,j−1 + ∆zi−1,j + ∆zi−1,j−1 × × (Vi,j−Vi−1,j) = = 1 2 ∆z′i− 1 2 ,j∆x + ∆z′i−1 2 ,j−1∆x + O ( ∆x2 ) ∆zi,j + ∆zi,j−1 + ∆zi−1,j + ∆zi−1,j−1 × ×V ′ i− 1 2 ,j ∆x+ + Vi−1,j + Vi,j 2 + O ( ∆x2 ) = = Vi−1,j + Vi,j 2 + O ( ∆x2 ) , что является аппроксимацией 2-го порядка V - компоненты скорости посередине между двумя со- седними узлами, в которых задаются Vi−1,j,k и Vi,j,k. Отметим, что и третья компонента поля ско- рости, определeнная по (8), также является аппро- ксимацией 2-го порядка. 64 А. Нестеров ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 Рис. 9. Рассчитанные уровни в сравнении с точным решением уравнений (t=1000 с) Рис. 10. Рассчитанные скорости в сравнении с точным решением уравнений (t=1000 с) 2. РЕЗУЛЬТАТЫ ПРИМЕНЕНИЯ МЕТОДА 2.1. Задача о течении в канале с острой дамбой Результат использования предложенной схемы ап- проксимации поля скорости с последующим при- менением схемы адвекции MUSCL в задаче моде- лирования течения в канале с острой дамбой, опи- санной во введении, показан на рис. 7. Коэффици- ент горизонтальной диффузии полагался равным нулю. Предложенная схема не производит осцил- ляций, возникавших при применении схемы цен- тральных разностей (см. в сравнении с pис. 1) и, в отличие от последней, является устойчивой. 2.2. Устойчивость в квазиодномерном случае Целью даного раздела является проверка устой- чивости предложенного метода и его сходимости к точному решению. Исследование этих свойств в трeхмерном случае затруднено вследствие не- линейности схемы. Осцилляции одной из компо- нент скорости, порожденные трeхмерной моделью локально в областях с большими еe градиента- ми, могут проникать во всю расчетную область, а также вызывать осцилляции других компонент скорости и переменных состояния. Устойчивость и сходимость схемы можно исследовать в квазио- дномерном случае течения в канале в резким пе- репадом глубины, описываемым частным случаем уравнений (1)–(3), которые в этом случае сводятся к уравнению для уровня поверхности η = η(x, t) и уравнению для осреднeнной по глубине скоро- сти потока U = U(x, t), при пренебрежении вязко- стными и кориолисовыми силами: ∂η ∂t + ∂UD ∂x = 0, ∂U ∂t + ∂ ∂x ( 1 2 U2 + gη ) = 0, (9) А. Нестеров 65 ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 Рис. 11. Батиметрия для задачи моделирования волнового движения где x – направление вдоль оси канала. В рассматриваемом примере канал с резким изменением батиметрии, показанной на pис. 8, по- крывался сеткой с разрешением ∆x = 10 м. Глу- бина в узлах численной сетки описывалась как Hi =      25.0м, 1 ≤ i ≤ 40, 2.5м, 41 ≤ i ≤ 59, 25.0м, 60 ≤ i ≤ 100. На численное решение влияют только значе- ния в узлах, поэтому глубину не обязательно счи- тать разрывной функцией; можно предполагать еe гладкое (например, линейное) изменение в под- сеточном масштабе. С физической точки зрения следовало бы рассматривать негидростатическую задачу, однако это привело бы к дополнительно- му уравнению, причeм уравнения зависели бы от двух пространственных координат, что существен- но усложнило бы задачу. В стационарном случае уравнения (9) упроща- ются и имеют следующие интегралы: UD = C1 U2 2 + gη = C2    , где C1, C2 – константы. С физической точки зре- ния C1 = Q является удельным расходом (пол- ным расходом, делeнным на ширину канала) че- рез сечение канала. Эти два уравнения сводятся к поиску корня многочлена. Поскольку стационар- ное решение для численной модели заранее неи- звестно, то в качестве начальных условий задава- лись нулевые скорости и уровень, а на открытых границах задавался удельный расход Q и условие излучения для уровня поверхности η. Расход пола- гался равным на втоке (“слева”) и вытоке (“спра- ва”), и с целью согласования начальных и грани- чных условий, а также с целью избежания “удар- ной волны”, описывался как: Q(t) = { 6.25t/400, 0 ≤ t ≤ 400, 6.25, 400 ≤ t. Удельный расход 6.25м2/с соответствует скорости потока 0.25м/с при глубине 25м (батиметрическая глубина на границах). Точное решение, соответ- ствующее стационарному расходу Q = 6.25м2/с при t � 400, имеет вид: (Ui, ηi) =      (0.249м/c, 0.081м), 1 ≤ x ≤ 40, (2.9м/c,−0.345м), 41 ≤ x ≤ 59, (0.249м/c, 0.081м), 60 ≤ x ≤ 100. Сходимость численного решения к этому ста- ционарному решению после окончания “разгона” обуславливается как видом и дискретизацией гра- ничных условий, так и наличием численных дисси- пации и дисперсии. На pис. 9 и 10 показано сравне- ние результатов численного моделирования с при- менением различных схем с точным решением при t=1000 с. В одномерном случае неустойчивость проявляется в возникновении незатухающих коле- баний скорости и уровня, которые могут привести к “развалу” модели, если возможность “осушения” не учитывается. Затухания колебаний, порождае- мых схемами центральных разностей и Super-C, не наблюдалось даже при t=10000 с. Колебания, по- рождаемые схемой Lax-Wendroff, либо затухали, и схема сходилась к точному решению лучше дру- гих схем либо приводила к “развалу” на этапе “ра- згона”, в зависимости от его длительности. Схема ULTIMATE-QUICKEST в процессе “разгона” так- же давала нелинейные колебания, но меньшей ам- плитуды, чем Lax-Wendroff, Super-C или централь- ные разности, после чего сходилась к стационарно- му решению, которое несколько отличалось от то- чного решения. Это отличие было меньше, чем при применении TVD-схем, которые проявили тенден- цию быстрой сходимости к стационарному реше- нию без возникновения колебаний, однако отлича- ющегося от точного. Численная диссипация этих схем привела к аналогу физической диссипации, что выражалось в наличии перепада уровня по- верхности (pис. 9). 2.3. Диссипативные свойства По аналогии с численной диффузией, в приме- нении к переносу импульса “наилучшей” можно считать схему с наименьшей численной диссипа- цией. Для исследования диссипативных свойств 66 А. Нестеров ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 Рис. 12. Диссипация полной относительной энергии системы со временем рассматривалась задача моделирования нелиней- ного волнового движения, описываемого уравени- ями (9) в канале с батиметрией, показанной на pис. 11 и описываемой как H = 25+20 sin(πx/100). Выбор такого профиля обусловлен желанием по- лучить относительно большие значения скорости и еe градиента. В качестве начальных условий за- давалась нулевая скорость и линейное изменение уровня поверхности η = 0.001x− 0.5 \|t=0. Полная энергия системы в предположении, что границы канала непроницаемы, будет E(t) = ∫ 1000 x=0 gη2 + DU2 2 dx = const. (10) В этом можно убедиться, если уравнение ∂ ∂t ( gη2 + DU2 2 ) + ∂ ∂x ( UD ( 1 2 U2 + gη )) = 0, вытекающее из (9), проинтегрировать по длине канала и принять во внимание, что на границах U = 0. Численная диссипация зависит не только от схемы расчета адвективных слагаемых, но и от схемы расчета уровня поверхности, дискрети- зации производных по времени – в данном при- мере использовалась производная 1-го порядка по времени и неявная схема расчeта уровня [7]. Рис. 12 демонстрирует изменение во времени полной энергии системы, отнесенной к полной на- чальной энергии. Точное решение соответствует единице, что следует из уравнения (10). С точ- ки зрения полной диссипации системы, схемы ра- спределились в следующем порядке её возраста- ния: Центральные разности, Lax-Vendroff, Super- C, Superbee, Ultimate-Quickest, MUSCL, Van Leer, Albada, MinMod, 1-го порядка вверх по потоку. ЗАКЛЮЧЕНИЕ Схема центральных разностей для численной ап- проксимации адвективных слагаемых уравнений импульса является непригодной для моделирова- ния течений с резкими изменениями скорости или направления. Работа показала, что даже в задачах геофизической гидродинамики схема дискретиза- ции переноса импульса может существенно влиять на результат, в частности на устойчивость и дис- сипативные свойства. Предложенный метод расче- та адвективных слагаемых уравнений импульса в моделях с регулярными разнесенными сетками оказался эффективным в практическом примене- нии, так как он может использовать те же самые подпрограммы расчета адвекции, что и скаляр- ные переменные. Применение TVD-схем адвекции не производит осцилляций поля скорости и явля- ется устойчивым. Из рассмотренных TVD-схем с точки зрения наименьшей диссипации наилучшим оказалось применение схемы Superbee. В случае, если задача не является критичной к устойчиво- сти, то возможно применение схем Lax-Wendroff или ULTIMATE-QUICKEST. Дальнейший интерес представляет разработка устойчивых схем, обла- дающих меньшей диссипацией, чем представлен- ные в этой работе. Автор выражает искреннюю благодарность А. Нестеров 67 ISSN 1561 -9087 Прикладна гiдромеханiка. 2008. Том 10, N 1. С. 58 – 68 доктору физ.-матем. наук, профессору В.С. Маде- ричу за ценные замечания и оказанную помощь в подготовке этой публикации. 1. Педлоски Дж. Геофизическая гидродинамика. В 2-х томах.– М.: Мир, 1984.– 398 c. (т.1), 416 c. (т.2). 2. Флетчер К. Вычислительные методы в динами- ке жидкостей. Т.1. Основные положения и общие методы.– М.: Мир, 1991.– 502 с. 3. Hirsch C. Numerical computation of internal and external flows, v2. Computational Methods for invi- scid and viscous flows.– New York: Wiley and Sons, 1990.– 714 p. 4. Tannehill J.C., Anderson D.A., Pletcher R.H. Computational Fluid Mechanics and Heat Transfer, 2nd Ed.– Washington: Taylor and Francis, 1997.– 792 p. 5. Кошебуцкий В., Мадерич В., Нестеров A., Хе- линг Р. Моделирование распространения тепла во внутренних водах и прибрежных областях мо- рей // Прикл. гидромеханика.– 2004.– 6(78).– С. 34–44. 6. Blumberg A.F.,Mellor G.L. A description of a three-dimensional coastal ocean circulation // Three- Dimensional Coastal Ocean Models.– 1987, N. Heaps (ed), Washington, D.C., Am. Geoph. Union.– P. 1–16. 7. Casulli V., Cheng R.T. Semi-Implicit Finite Di- fference Methods for Three-Dimensional Shallow Water Flow // Int. J. for Numer. Meth. in Fluids.– 1992.– 15.– P. 629-648. 8. Stelling G.S., Duinmeijer S.P.A. A staggered conservative scheme for every Froude Number in rapi- dly varied shallow water flows // Int. J. for Num. Meth. in Fluids.– 2003.– 43.– P. 1329–1354. 9. Van Leer B. Toward the ultimate conservative difference scheme. V: A second order sequel to Godunov’s method // J. Comput. Phys.– 1979.– 32.– P. 101–136. 10. Fringer O-B., Armfield S.W., Street R.L. Reduci- ng numerical diffusion in interfacial gravity wave si- mulations // Int. J. Numer. Meth. Fluids.– 2005.– 49.– P. 301–329. 68 А. Нестеров
id	nasplib_isofts_kiev_ua-123456789-4631
institution	Digital Library of Periodicals of National Academy of Sciences of Ukraine
issn	1561-9087
language	Russian
last_indexed	2025-12-07T18:15:43Z
publishDate	2008
publisher	Інститут гідромеханіки НАН України
record_format	dspace
spelling	Нестеров, А.А. 2009-12-08T15:23:45Z 2009-12-08T15:23:45Z 2008 Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками / А.А. Нестеров // Прикладна гідромеханіка. — 2008. — № 1. — С. 58-68. — Бібліогр.: 10 назв. — рос. 1561-9087 https://nasplib.isofts.kiev.ua/handle/123456789/4631 Схема центральных разностей, часто применяемая для численной аппроксимации адвективных слагаемых уравнений движения в моделях геофизической гидродинамики с разнесенными сетками, может является причиной возникновения нефизических осцилляций поля скорости и приводить к неустойчивости. В работе предложен метод расчeта адвективных слагаемых уравнений движения на основе специальной интерполяции поля скорости, с последующим применением какой-либо из схем уравнений переноса, включая TVD схемы. Метод рассматривается в применении к задачам моделирования течений в водных резервуарах со свободной поверхностью. Схема центральних рiзниць, яка часто використовується для чисельної апроксимацiї адвекцiйних доданкiв рiвнянь руху в моделях геофiзичної гiдродинамiки з рознесеними сiтками, може бути причиною виникнення нефiзичних осциляцiй поля швидкостi i призводити до нестiйкостi. В роботi запропоновано метод розрахунку адвекцiйних компонент рiвнянь руху на основi спецiальної iнтерполяцiї поля швидкостi, з наступним використанням будь-якої зi схем рiвнянь переносу, включаючи TVD схеми. Метод розглядається в застосуваннi до задач моделювання течiй у водних резервуарах з вiльною поверхнею. Scheme of central differences, which is often used for numerical approximation of the advection terms in the momentum equations in geophysical hydrodynamics models with staggered grids, can induce non-physical oscillations of the velocity field and cause numerical instabilities. This paper proposes a numerical method for momentum advection approximation, which is based on the specific interpolation of the velocity field with subsequent application of any scheme for a scalar transport equation, including TVD schemes. This method is studied in the application to the problems of flow modeling in the water reservoirs with free surface. ru Інститут гідромеханіки НАН України Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками Application of TVD-schemes for approximation of momentum transfer in geophysical models of hydromechanics with spaced grids Article published earlier
spellingShingle	Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками Нестеров, А.А.
title	Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками
title_alt	Application of TVD-schemes for approximation of momentum transfer in geophysical models of hydromechanics with spaced grids
title_full	Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками
title_fullStr	Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками
title_full_unstemmed	Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками
title_short	Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками
title_sort	применение tvd-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками
url	https://nasplib.isofts.kiev.ua/handle/123456789/4631
work_keys_str_mv	AT nesterovaa primenenietvdshemdlâapproksimaciiperenosaimpulʹsavmodelâhgeofizičeskoigidrodinamikisraznesennymisetkami AT nesterovaa applicationoftvdschemesforapproximationofmomentumtransferingeophysicalmodelsofhydromechanicswithspacedgrids

Применение TVD-схем для аппроксимации переноса импульса в моделях геофизической гидродинамики с разнесенными сетками

Institution

Similar Items