Содержание к диссертации
Введение
Глава 1. Вычислительные методы в механике однофазной жидкости 15
1.1 Введение 16
1.2 Основные уравнения и метод численного решения 19
1.3 Моделирование турбулентности 38
1.4 Ускорение гидродинамических расчетов 70
1.5 Базис валидации для задач механики однофазной жидкости 81
1.6 Пакет программ ЛОГОС 105
1.7 Выводы к главе 1 109
Глава 2. Неявный метод численного решения уравнений Навье-Стокса для расчета многофазных течений со свободной поверхностью 111
2.1 Введение 112
2.2 Основные уравнения модели и описание метода 118
2.3 Параллельная реализация метода в пакете программ ЛОГОС 131
2.4 Верификация метода 151
2.5 Учет сил гравитации в задачах со свободной поверхностью 171
2.6 Исследование свойств схем дискретизации уравнения переноса объемной доли 186
2.7 Оценка параметров метода 199
2.8 Выводы к главе 2 213
Глава 3. Моделирование цунами космогенного происхождения 215
3.1 Введение 216
3.2 Астероидно-кометная опасность и методы расчета соударения с гидросферой 223
3.3 Моделирование космогенных цунами источниками различных типов 251
3.4 Моделирование падения метеорита в воду при вхождении под углом 261
3.5 Моделирование падения метеорита в озеро Чебаркуль в 2013 году 274
3.6 Выводы к главе 3 289
Глава 4. Моделирование цунами оползневого происхождения 291
4.1 Введение 292
4.2 Сеточные модели для моделирования цунами на основе уравнений Навье-Стокса 295
4.3 Верификация метода для расчета цунами оползневого типа 311
4.4 Исторические и экспедиционные данные о цунами на острове Монтсеррат 319
4.5 Моделирование Монтсерратского цунами 2003 года 336
4.6 Выводы к главе 4 364
Заключение 365
Работы автора по теме диссертации 367
Список литературы
- Моделирование турбулентности
- Параллельная реализация метода в пакете программ ЛОГОС
- Моделирование космогенных цунами источниками различных типов
- Исторические и экспедиционные данные о цунами на острове Монтсеррат
Моделирование турбулентности
Механика сплошных сред в ходе своего развития и приложения к различным областям науки и техники, а также разрастающихся связей с другими частями физики, химии, биологии, океанологии и другими областями науки, начала подразделяться на самостоятельные направления своего развития. В область самостоятельных исследований можно выделились такие направления, как гидродинамика, аэродинамика, газовая динамика, теплогидравлика, геофизика и другие, различие между которыми обусловлено предпосылками и допущениями при модификации основных уравнений динамики жидкости применительно к расчетам данного класса задач. В свою очередь, модификация базовой системы уравнений и физико-математической модели численного решения на ее основе не могли не сказаться на разработке алгоритмов математического моделирования.
Первые численные задачи динамики жидкостей (одномерные) исторически решались в переменных Лагранжа. Это объясняется тем, что численные схемы в переменных Эйлера являются не столь точными в задачах с множеством контактных границ, а при моделировании многокомпонентных течений может наблюдаться свойство «нефизичного» увеличения диффузии в смешанных ячейках счетной области. При моделировании сжимаемых течений в переменных Эйлера ухудшается математическое определение физической системы вследствие того, что координаты представляют фиксированные точки пространства, а не среды. Численные схемы, построенные на основе подхода Эйлера, обладают свойством устойчивости счета, хотя в некоторых задачах могут уступать по точности, особенно в области границы разделов двух сред.
Однако все эти преимущества переменных Лагранжа оказываются малосущественными при описании неустойчивых и сдвиговых течений (турбулентность, струйные течения, пульсации, опрокидывание волн и др.). При моделировании таких течений в переменных Лагранжа элементарные ячейки могут в очень сильной степени вырождаться и принимать сколь угодно сложную структуру, до такой степени, что счет на такой расчетной сетке становится невозможным. В некоторых практических задачах газодинамики, например в физике высокотемпературной плазмы, «лагранжев» подход используется до того момента времени, до которого позволяет расчетная сетка, после чего переходят на «эйлерову» методику. Для борьбы с искажением сеток применяются методы переинтерполяции, использование подвижных сеток, нерегулярных сеток и метод частиц. Для трехмерных задач эта проблема встает еще более остро.
По мере увеличения вычислительных мощностей и усовершенствования «эйлеровых» методик, особенно в части аккуратного счета контактных границ, им стало отдаваться все большее предпочтение. Что касается гидродинамики поверхностных волн, то «лагранжев» подход здесь применять, вообще говоря, очень проблематично - главной причиной является опрокидывание волн, которое сразу приведет расчетную сетку в негодность. За основу трехмерного численного моделирования гидродинамики поверхностных волн целесообразно принять именно «эйлеров» подход, на базе которого и строится все дальнейшее изложение.
Подходы к построению численной метода расчета динамики жидкости можно условно разделить на два направления, широко применяемые в настоящее время, а именно решение полной системы уравнений Навье-Стокса и ее расщепление.
Решение полной системы применяют для расчета задач аэродинамики - описания транс-, сверхзвуковых и гиперзвуковых течений вязкого сжимаемого газа. Численные схемы, применяемые для расчета этого класса задач, отличаются от схем, используемых для моделирования течений несжимаемой жидкости. Точность определения характеристик сжимаемых течений напрямую зависит от типов выбранных схем расчета потоков, которые и порождают различные версии методов дискретизации для решения уравнений Навье-Стокса. На практике здесь обычно применяют схемы, имеющие ясную физическую интерпретацию и использующие решения задачи Римана [Годунов, 1976; Флетчер, 1991; Куликовский и др., 2001].
При дискретизации полной системы уравнений Навье-Стокса для расчета слабосжимаемых течений, характеризующихся низкими числами Маха, возникают хорошо известные проблемы, связанные с существенно различающимися временами конвективного переноса и распространения акустических возмущений, а также с малыми изменениями относительного давления. В этом случае при использовании классических методов решения полной системы уравнений Навье-Стокса, исторически разрабатываемых для решения задач сверхзвуковой аэродинамики, резко ухудшаются свойства их сходимости и устойчивости, а некоторые методы вообще становятся непригодными для расчета течений с числами Маха, стремящимися к нулю. Более того, здесь практически становится невозможным использование явных схем расчета в силу неприемлемо жесткого ограничения на шаг по времени. Использование же неявных схем, которые в линейном приближении являются безусловно устойчивыми, снимает проблему лишь отчасти по причине существенного снижения области устойчивости по мере снижения числа Маха и замедления сходимости алгоритмов даже при оптимальном выборе шага по времени. Проблему усугубляет тот факт, что для существенно дозвуковых течений характерны чрезвычайно малые относительные изменения давления, вследствие чего расчет градиентов давления содержит операцию вычисления малых разностей близких между собой абсолютных значений статического давления Для решения обозначенных проблем было создано значительное количество различных вычислительных алгоритмов [Стрелец&Шур, 1988; Choi&Merkle, 1985, 1993; Weiss&Smith, 1995; Turkel, 1987, 1993, 1999, 2002; Флетчер, 1991; Merkle, 1995; Ferziger&Peric, 2002], и, судя по накопленному к настоящему времени опыту расчетов самых разнообразных течений [Флетчер, 1991; Ferziger&Peric, 2002; Быстров и др., 2005], одним из наиболее эффективных методов рассматриваемого класса является алгоритм SIMPLE, предложенный в [Patankar&Spalding, 1972]. Наряду с возможностью обобщения этого метода на случай сильносжимаемых течений его важным достоинством является эффективный способ неявной аппроксимации уравнений, позволяющий применять его при сколь угодно больших временах расчета. Именно этот метод и был положен в основу технологии расчета гидродинамики волн цунами космогенного и оползневого происхождения на базе полных уравнений Навье-Стокса (в случае турбулентных течений - уравнений RANS/DES/LES моделей), разработанной диссертантом и реализованной на произвольных неструктурированных сетках в рамках отечественного пакета программ ЛОГОС.
В параграфе 1.2 данной главы диссертации изложена дискретизация системы уравнений Навье-Стокса на основе метода SIMPLE, применяемые численные схемы для аппроксимации конвективных и диффузионных потоков. Представлены особенности численной реализации на произвольных неструктурированных сетках. В параграфе 1.3 данной главы диссертации изложены особенности решения систем линейных алгебраических уравнений, возникающих в результате дискретизации системы уравнений Навье-Стокса. Представлены способы ускорения расчетов. В параграфе 1.4 детально проанализированы методы расчета турбулентности, приводятся результаты калибровки и исследование применимости различных подходов к моделированию турбулентности, исследуется поведение различных схем аппроксимации порядка точности не менее второго на неструктурированных сетках, состоящих из многогранников произвольной формы. Даны рекомендации по выбору численных схем и входящих в модели констант. В параграфе 1.5 представлен базис задач валидации для турбулентных течений вязкой несжимаемой жидкости, который представляет самостоятельный интерес и может оказаться весьма полезным для исследователей и пользователей различных пакетов программ математического моделирования. Описание пакета программ ЛОГОС, в рамках которого реализована разработанная технология, представлено в параграфе 1.6. Полученные результаты суммированы в параграфе 1.7.
Параллельная реализация метода в пакете программ ЛОГОС
Результаты показывают, что использование в основной области неструктурированных сеток, требует уменьшения характерного размера ячеек в основной области, что заметно повышает необходимое общее число ячеек. Для получения приемлемых результатов на полиэдральной сетке число ячеек в основной области должно быть, как минимум, в 2-3 раза больше, чем на гексагональной сетке, а характерный размер полиэдральных ячеек должен быть в 1.3 раза меньше, чем характерный размер шестигранников. Для сетки, составленной из треугольных призм, приемлемый результат был получен лишь на сетке, содержащей в 4 раз больше ячеек с характерным размером 0.025 м, что в 1.4 раза меньше, чем характерный размер ячеек блочно-структурированной сетки. Использование тетраэдральной сетки требует уменьшения характерного размера в 1.6 раза, до 0.021 м. Чтобы обеспечить такой характерный размер тетраэдров в основной области LES требуется около 6 млн расчетных ячеек, что в 6 раз больше, чем число ячеек на блочно-структурированной сетке.
Необходимо сказать, что в данной задаче рассматривались лишь осредненные характеристики течения, а именно осредненный коэффициент трения. В случае расчета нестационарных характеристик течения использование неструктурированных сеток может потребовать дальнейшего уменьшения характерного размера ячеек.
Также необходимо отметить, что расчетная область данной задачи не содержала регионов с «плохой» геометрией (острые кромки, плохообтекаемые элементы) и регионов с сильно изменяющимися градиентами скоростей, поэтому схема 0.9CD+0.1UD была вполне приемлема для устойчивого счета. В других случаях ее применение совсем не гарантирует устойчивости с такими коэффициентами смешения. Увеличение в ее формулировке доли «противопоточности» (до 20-30%) не решает проблему, так как схема становится сильно диссипативной. Отыскание схемы, которая обеспечивает приемлемую диссипацию и монотонность на произвольных неструктурированных сетках, является темой для дальнейших исследований. Альтернативная низкодиссипативная численная схема По сравнению с турбулентной вязкостью, определяемой при использовании RANS-моделей, подсеточная вязкость, как правило, имеет на порядок меньшие значения, сопоставимые со значениями схемной вязкости, возникающей из-за решения дифференциальных уравнений на разностных сетках с конечным уровнем разрешения. В этом состоит главная проблема метода LES – величина подсеточной вязкости, которую необходимо привнести в уравнение движения напрямую зависит от схемной вязкости используемой численной схемы – чем больше схемная вязкость, тем меньше должен быть уровень дополнительной подсеточной вязкости. Стоит отметить, что далеко не каждая численная схема вообще способна обеспечить приемлемое значение схемной вязкости (или иначе – численной диссипации), не превышающее реальной диффузии, вызванной наличием вихрей малых масштабов.
Использование центрально-разностных схем для моделирования турбулентных течений приводит к численным осцилляциям [Jasak et al., 1999], которые, как правило, развиваются в пограничном слое и вблизи границ вход/выход. Уменьшение осцилляций и увеличение устойчивости центрально-разностной схемы можно добиться ее смешением с противопоточной [Leonard, 1979]. Данный способ широко применим на практике, однако, как показано выше, даже внесение 10% противопоточности увеличивает численную диссипацию схемы примерно в 2 раза.
Другой способ решения проблемы представлен в работе [Weinman& Valentino, 2010]. В ней предложена центрально-разностная схема с добавлением искусственной вязкости. Показано, что данная схема приводит к достаточно хорошим результатам для свободных течений при использовании гексагональных сеток, проверка ее диссипативности в случае пристеночных течений не приводится. В [Strelets, 2001] описывается центрально-разностная схема с автоматическим смешением с противопоточной схемой, причем фактор смешения строится таким образом, что центрально-разностная схема включается в тех областях, где она наиболее устойчива и необходима для точного разрешения крупномасштабных вихревых структур. В пристеночной области, а также в области свободной от отрывных вихрей используется противопоточная схема первого порядка. Хорошим способом монотонизации схем высокого порядка является применение ограничителей. В работе [Jasak et al., 1999] описывается схема GAMMA, основанная на центрально-разностной схеме, ограниченной по критерию локальной ограниченности (Convection Boundedness Criterion, CBС [Leonard, 1991]) с применением диаграмм нормализованной переменной NVD (Normalized Variable Diagram). Для данной схемы гарантируется монотонность, однако ограничение по CBС вносит дополнительную диффузию. В проведенном в рамках диссертационной работы исследовании, представленном выше, показано, что схема GAMMA слишком диссипативна для областей LES-моделирования.
В рамках цикла исследований по моделированию турбулентности при руководстве и непосредственном определяющем участии диссертанта предложена альтернативная низкодиссипативная схема, основанная на центральных разностях [СВ24]. Ее формулировка выполнена с использованием диаграмм NVD, однако по сравнению с GAMMA в ней применяется более «мягкое» ограничение, основанное на увеличении доли противопоточности схемы лишь при увеличении отклонения нормализованной переменной от ее допустимого интервала. Также в схеме используется повышение уровня ограниченности в зоне пограничного слоя, позволяющее устранять численные осцилляции в зонах больших градиентов скоростей.
Для предложенной схемы показано, что на задаче конвективного переноса резкого фронта пассивного скаляра она не приводит к возникновению численных осцилляций по сравнению с центрально-разностной, а также не «размывает» фронт, как центрально-разностная схема с добавлением противопоточности. Cравнение численной диссипации предложенной схемы и калибровка для нее констант модели LES для свободных течений проведено на задаче о вырождении изотропной турбулентности. Калибровка для пристенных течений проведена на задаче о турбулентном течении в плоском канале, решение которой выполнено на последовательно сгущающихся расчетных сетках.
Моделирование космогенных цунами источниками различных типов
По-другому дело обстоит с многофазной гидродинамикой, являющейся частым явлением во многих приложениях, таких как течения в ядерных реакторах, теплообменниках, двигателях, сопровождаемые конденсацией, испарением, кипением и другими процессами. Численные проблемы, возникающие при моделировании таких явлений гораздо сложнее, а высокие скорости, наличие турбулентности и фазовые превращения лишь усугубляют проблему.
Волны цунами, тоже не являются исключением, особенно это касается цунами несейсмического происхождения. При моделировании космогенного цунами, когда в воду на огромной скорости падает разогретый до нескольких сотен градусов метеорит, необходимо учитывать процессы испарения воды, тепло-, массообмен с атмосферой, влияние воды на процесс остывания тела, а также возможное взаимодействие нескольких отколовшихся фрагментов. Вулканическое цунами, порождаемое взрывом подводных и надводных вулканов, выбрасывает в холодную воду большие массы горячего газа и лавы, и построить модель для описания процесса генерации волны (в случае подводного вулкана добавляется еще эволюция и схлопывание газового пузыря), а также возможного обвала фрагмента дна - отнюдь не тривиальная задача.
При описании цунами обычно выделяют три стадии: генерация волны, распространение и накат. До настоящего времени моделирование этих стадий осуществляется по разным моделям в силу различия протекающих в них физических процессов, с одной стороны, и невозможности применения полных трехмерных уравнений гидродинамики из-за ограниченности вычислительных ресурсов, с другой. Однако повсеместный рост вычислительных мощностей и их планомерное увеличение до экзафлопсного уровня [сайт проекта TOP500] позволяет (хотя бы частично) пересмотреть данную концепцию. Это позволит продвинуться в понимании физики цунами и перейти на качественно новый уровень их изучения, что особенно важно при изучении цунами несейсмического происхождения.
Описание распространения волн цунами в открытом океане проработано наиболее полно и в настоящее время является уже делом техники. Это связано с тем, что высота цунами в открытом океане намного меньше глубины. В этом случае фазовой и амплитудной дисперсией пренебрегают и применяют для описания распространения линейную или нелинейную теорию длинных волн [Пелиновский, 1996; Левин&Носов, 2005]. Нелинейные эффекты начнут проявляться при достижении волной шельфовой зоны, а также при распространении вдоль береговой линии.
Физическая постановка наката цунами на берег не представляет сложности, однако известные трудности представляют учет обрушения, турбулентности и описание взаимодействия воды и воздуха. Для расчета затопления имеются численные модели, качество применения которых ограничивается степенью детальности топографии района и доступностью вычислительных ресурсов [Левин&Носов, 2005].
В противоположность стадиям распространения и наката совсем иная ситуация обстоит с источником цунами, описание которого до сих пор является нерешенной задачей. Даже для цунами сейсмического происхождения, которые моделируют больше всего, напрямую процесс генерации не считают, а используют механизм генерации возмущения на поверхности воды, форма которого аналогична остаточным деформациям дна. Сами же остаточные деформации рассчитываются на основании параметров очага землетрясения по формулам Окады [Okada, 1985]. Эта модель, естественно, не учитывает возможность дополнительного вклада в источник, например от оползней и обвалов морского дна, спровоцированных землетрясением, а такой вклад может быть весьма существенным [Gusiakov, 2001].
Что касается источников цунами космогенного и оползневого типа, то современные представления о них, а тем более модели, сложно назвать однозначными. Для генерации цунами космогенного происхождения используются в основном параметрические модели [Kharif&Pelinovsky, 2005; Ward&Asphaug, 2000; Левин&Носов, 2005; Badescu&Isvoranu, 2011].
Возможности современных программных кодов для расчета падения тела в воду при больших скоростях можно оценить по результатам работы [Pierazzo et al., 2008]. В работе представлено описание ряда многофункциональных программных комплексов и заложенных в них моделей. Здесь при падении стеклянного шарика со скоростью 4,64 км/с моделируется только первоначальная стадия образования каверны на поверхности до начала её схлопывания. Схлопывание каверны и эволюция газового пузыря не рассматривались. При моделировании по всем кодам, кроме одного, атмосфера не учитывалась, а сжимаемость воды вычислялась по табличному уравнению состояния. В такой постановке средняя погрешность вычислений этой стадии по сравнению с экспериментом составила 15%. Погрешность расчета программой SOVA, используемой для моделирования космогенных цунами в [Shuvalov&Trubetskaya, 2002, 2007; Shuvalov, 2012] составила около 25%. Учет атмосферы при моделировании по одной из представленных программ сразу же привел к погрешности в 50%. Представленные факты говорят сами за себя - необходимо дальнейшее развитие моделей и методов численного моделирования для данного класса задач.
С источниками цунами оползневого типа ситуация гораздо лучше. Цунамигенные оползни условно можно поделить на три категории: надводные, частично погруженные в воду и подводные. Начальное положение оползня является основой для выбора физико-математической модели, приемлемой для описания генерации и распространения цунами. Надводные и частично подводные оползни целесообразно описывать трехфазной системой «жидкость-воздух-оползень», тогда как для подводного оползня достаточно ограничится двухфазной или двухслойной моделью со слоями различных плотностей. При этом сам оползень моделируется недеформируемым твердым телом или системой таких тел, а также несжимаемой жидкостью либо отдельным слоем со своими значениями коэффициентов плотности и вязкости [Imamura&Imteaz, 1995; Heinrich et al., 1998; Watts&Grilli, 2003; Федотова и др., 2004; Chubarov et al., 2011].
Поверхностные волны, порождаемые сходом оползня, также имеют свою специфику. Зарождение волны в прибрежной зоне с малой глубиной осуществляется в достаточно длительное время, сравнимое с временем перемещения оползня, а характерный размер оползня зачастую сравним с глубиной. В отличие от цунами сейсмического происхождения цунами оползневого типа являются более короткими [Dutykh&Dias, 2009], что обуславливает необходимость учета дисперсии волн. Для моделирования таких волн используют нелинейно-дисперсионные уравнения теории мелкой воды, способные воспроизводить дисперсию. Эти уравнения решаются конечно-разностными методами, построенными на базе схем второго порядка точности [Федотова, 2004]. Однако по причине того, что системы такого типа содержат смешанные производные высокого порядка, построение эффективных численных алгоритмов для их решения является нетривиальной задачей. Проблематика применения нелинейно-дисперсионных уравнений для моделирования цунами оползневого типа обсуждается в [Федотова и др., 2004; Chubarov et al., 2011; Гусев и др., 2013].
Исторические и экспедиционные данные о цунами на острове Монтсеррат
Определяющей характеристикой для течений со свободной поверхностью является сила гравитации. Сила гравитации в области свободной поверхности терпит резкие изменения вследствие существенного изменения плотности среды, что приводит к разрыву градиента давления, который в случае покоя среды полностью уравновешивает действие гравитационных сил [Ландау&Лифшиц, 1988]. Построение численного алгоритма, обеспечивающего корректный учет силы гравитации и расчет значений градиента давления, является нетривиальной задачей. Особенно это касается расчетных сеток с «коллокированным»FF5 расположением неизвестных величин, которое в основном используется на практике, однако приводит к слабой связи поля скорости и давления [Ferziger&Peric, 2002; Jasak, 1996; Флетчер, 1991]. Использование «коллокированного» расположения неизвестных величин подразумевает определение давления и скорости в одном и том же месте (как правило, центр ячейки), что приводит к возникновению четно-нечетных осцилляций, избавиться от которых можно путем использования метода типа Рхи-Чоу [Rhie&Chow, 1983].
Вопросу построения численного алгоритма, обеспечивающего отсутствие численных осцилляций в случае наличия неоднородного поля силы тяжести вследствие неоднородного поля плотности, посвящен целый ряд работ [Gu et al, 1991; Mencinger, 2012; Храбрый, 2014]. В [Gu et al, 1991] рассматривается вопрос учета объемной силы в уравнении сохранения импульса, в качестве которой может выступать и сила гравитации. С целью исключения осцилляций в поле скорости и давления предлагается использовать алгоритм, схожий с поправкой Рхи-Чоу [Rhie&Chow, 1983], что решает проблему. Однако в статье не рассматриваются задачи с сильным разрывом в поле объемной силы, в случае наличия которых возникает вопрос корректного расчета градиента давления, значение которого должно полностью уравновешивать объемную силу, когда среда находится в покоящемся состоянии.
Теоретическое рассмотрение проблемы учета силы гравитации представлено в [Mencinger, 2012]. Вкупе с идеями, заимствованными из [Gu et al, 1991], предлагается выражение для интерполяции объемной силы и давления на внутренние грани расчетной сетки, которое обеспечивает условие равновесия. Однако в рассмотренных задачах
Термин «коллокированный» в численных методах означает, что все неизвестные величины хранятся в центре расчетных ячеек. Проблема совмещенного хранения величин скорости и давления при численном решении уравнений Навье-Стокса хорошо известна [Флетчер, 1991]. Данную трудность можно преодолеть введением так называемой «разнесенной» сетки, когда давление хранится в центре ячейки, а скорости в ее вершинах. Однако данный подход невозможно применить на произвольных неструктурированных сетках, состоящих из многогранников произвольной формы, хотя бы по причине того, что заранее неизвестно, сколько вершин будет у ячейки, поэтому здесь хранят величины только в центрах ячеек. продемонстрировано лишь отсутствие осцилляций в поле скорости и не представлено анализа полученных форм свободной поверхности. Также не рассмотрена эффективность применения алгоритма при использовании произвольных неструктурированных сеток. К тому же в [Mencinger, 2012] используется поправка типа Рхи-Чоу для нестационарного слагаемого уравнения сохранения импульса, предложенная в [Majumdar, 1988], однако в [Яцевич и др., 2015] показано, что ее использование может приводить к искажению формы свободной поверхности.
В [Храбрый, 2014] представлена эффективная схема вычисления градиента давления в случае наличия сил гравитации, которая основана на интерполяции значения градиента давления с учетом плотностей среды в смежных ячейках. Данный алгоритм позволяет устранить нефизичные осцилляции в поле скорости вблизи положения свободной поверхности, однако работоспособность алгоритма продемонстрирована лишь на ортогональных расчетных сетках, линии которой параллельны направлению силы тяжести.
В настоящем параграфе формулируется алгоритм построения уравнения для давления, который основан на методе Рхи-Чоу. Алгоритм строится посредством замены в уравнении давления интерполяции силы гравитации её прямым дискретным аналогом. Выражения для прямого дискретного аналога формулируется на основе гидростатического приближения, что обеспечивает корректное поле давления на произвольной неструктурированной сетке. Для обеспечения равновесия силы гравитации и градиента давления в случае покоя среды предложен алгоритм, основанный на замене градиента давления в уравнении движения его модификацией, содержащей учет действия гравитационной силы.
Эффективность предложенных решений исследуется на примере численного решения задачи о равновесии двухфазной среды в поле силы тяжести, имеющей аналитическое решение. Показано, что предложенные численные схемы, в отличие от существующих, позволяют получить корректное гидростатическое давление и обеспечить равновесие силы гравитации и градиента давления на неструктурированных расчетных сетках. Точность предсказания форм свободной поверхности с использованием предложенных численных алгоритмов оценивается путем решения задачи о колебании жидкости в поле силы гравитации и задачи об обрушении столба жидкости.