ВИХРЕВАЯ ГИДРОДИНАМИКА: НОВЫЙ ПОДХОД К МОДЕЛИРОВАНИЮ ГЕОСИСТЕМ Текст научной статьи по специальности «Математика»

ВЕСТНИК ПЕРМСКОГО УНИВЕРСИТЕТА

2018 Математика. Механика. Информатика Вып. 1(40)

МЕХАНИКА МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

УДК 551.24:550.31:556.013:532.50: 531.5:519.63

Вихревая гидродинамика: новый подход к моделированию геосистем

В. И. Гунин

Центр моделирования геосистем "МоГеос"

Россия, 670034, г. Улан-удэ, пр. 50 лет Октября, 38-19

vigunin@list.ru

Для вхождения в эру устойчивого развития необходимо научиться управлять природой, т.е. оптимизировать ее влияние на человека и ее отклик на антропогенные воздействия. Для этого нужно уметь моделировать все природные процессы, а также научиться прогнозировать их развитие во времени и пространстве, постепенно вводя элементы управления. Как известно, все природные процессы в той или иной степени связаны с вихревыми гидродинамическими потоками, вызванными перераспределением вещества в гравитационном поле. Используя математический аппарат, базирующийся на системе уравнений Навье-Стокса, сложно описать вихревые (кольцевые) структуры, формирующиеся во всех геосферах нашей планеты под действием гравитационного поля. Предложен математический аппарат вихревой гидродинамики (ВГ), основанный на новом фундаментальном свойстве гравитационного поля - векторном потенциале, для которого го1А Ф 0, позволяющий адекватно описывать вихревые структуры и глубже понимать их сущность. С помощью анализа системы уравнений (ВГ) и результатов, полученных при решении ряда задач, сделаны обобщения по некоторым физическим процессам и явлениям, распространенным в различных геосферах Земли.

Ключевые слова: гравитация; векторный потенциал; конвекция; вихревые структуры; математическая модель; численный эксперимент; тепломассоперенос; мантия; литосфера; кора; плюм; диа-пир; рифт.

DOI: 10.17072/1993-0550-2018-1-5-18

Введение

В последние десятилетия активно развивается идея (концепция) устойчивого развития. Существует множество определений и трактовок этой концепции, но ни политики, ни экономисты, ни философы, разработавшие эти определения, до сих пор не могут понять и объяснить суть этой концепции, она становится просто модной теорией для управленческой "элиты". Некоторые страны уже считают, что вышли на путь устойчивого развития. В России объявлен даже конкурс "Регионы -устойчивое развитие". Но природа своей непредсказуемостью, стихийные бедствия, терроризм, экономические кризисы, войны и др. расставляют все на свои места.

С математической точки зрения суть концепции устойчивого развития - это оптимальное управление природой т. е. необходи-

© Гунин В. И., 2018

мо предсказать событие (явление), спрогнозировать его развитие и оптимизировать последствия его воздействия. Под термином природа понимается вся окружающая среда, включая человека. Это минимаксная задача в результате решения которой при минимальном воздействии на окружающую среду можно получить от нее максимальную пользу. Для выполнения этого, на первый взгляд простого определения, необходимо проделать громадную работу, которая по времени займет десятилетия или даже столетия. Необходимо разработать математические модели всех процессов, проходящих в окружающей среде и связи между ними, научиться их прогнозировать во времени и пространстве, а затем постепенно, изучая и корректируя, вводить элементы управления этими процессами.

Для моделирования процессов в окружающей среде нужны адекватные математические модели и эффективные численные методы. На наш взгляд, необходимо использо-

вать разные подходы для описания и решения этих задач, что в свою очередь позволит глубже изучить и понять природные явления и более точно определить направление движения к истине. Только тогда можно будет говорить о начале вхождения в эру устойчивого развития.

Ни для кого не секрет, что весь наш мир, все, что нас окружает, и мы сами находимся под воздействием гравитации. Все процессы, протекающие в природе, в той или иной степени связаны с гидродинамическими потоками, вызванными перераспределением вещества в гравитационном поле. Основной силой переноса вещества в природе является конвекция, которая в свою очередь зависит от сил, вызванных какими-то внешними воздей-

ствиями (вынужденная конвекция), и внутренних сил, вызванных плотностными неод-нородностями, возникающими за счет градиентов температуры и концентрации вещества (свободная конвекция).

Вихревые, кольцевые структуры довольно широко распространены во всех геосферах Земли и космическом пространстве. Особенно ясно это стало при дешифровке космических снимков в конце XX в. Размеры этих структур могут достигать сотен километров. Считается, что в земной коре данные структуры контролируют крупные месторождения металлов и других полезных ископаемых, а в гидросфере и атмосфере являются основными переносчиками энергии и вещества.

Рис. 1. Вихревые структуры, образованные при формировании Земной коры

В земной коре вихревыми называют такие структуры, в которых существенная роль принадлежит спиралевидным, дугообразным и эксцентричным кольцевым формам. В этих структурах всегда присутствуют две морфологически различные части: меньшая - внутренняя, и большая - внешняя. Внутренняя часть (часто округлой формы) занимает не более одной пятой, а нередко - меньше одной десятой всей площади структуры, и может

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

Рис. 2. Вихревые (кольцевые) структуры, а - в атмосфере, б - на льду озера Байкал (космический снимок) Генезис вихревых, кольцевых структур ного внедрения магмы до космического удар-пока не ясен. Существует множество гипотез но-метеоритного. Многие связывают образо-их образования от вращательно-поступатель- вание структур с вращательными движениями

блоков земной коры в горизонтальном направлении. В большинстве случаев в центре или в пределах данных структур фиксируются интрузивные тела, которые, скорее всего, и являются причиной образования вихревых, кольцевых и винтовых структур, но механизм формирования условий для их образования пока остается дискуссионным [2] (рис. 1).

С другой стороны, в атмосфере и гидросфере все мы явно наблюдаем различные вихревые структуры - циклоны, антициклоны, торнадо, бора и тому подобное. Это позволяет предположить, что и все вихревые структуры, связанные с глубинными процессами, развиваются по тому же сценарию, что и на поверхности Земли (рис. 2).

Но до сих пор эти процессы, идущие во всех оболочках Земли, недостаточно изучены, так как для описания и моделирования вихревых структур нет адекватного математического аппарата.

Используемые методы, проблемы:

новый подход

Для описания гидродинамики в задачах со свободной и вынужденной конвекцией в сплошных (вязких) и пористых средах широко используется система уравнений Навье-Стокса, в которой эти силы описываются совместно -одним уравнением. Это приводит к трудностям при ее решении, так как она не является системой типа Коши-Ковалевской, т. е., нет эволюционного уравнения для давления. Имеются и другие особенности, возникающие при использовании численных методов решения данной системы уравнений, например наличие малого параметра при второй производной может приводить к существенному изменению гладкости решения, и ряд других.

Уравнения Навье-Стокса, вытекающие из уравнения Больцмана, полученного на основе уравнений механики и статистики, непосредственно описывают только ламинарные течения, поэтому с помощью этой системы уравнений сложно, а в большей части невозможно описать вихревые структуры, которые формируются в гравитационном поле во всех средах при возникновении плотностной стратификации.

Для описания вихревых структур (турбулентности), привлекают дополнительные гипотезы и предположения, которые повышают степень этой системы уравнений, запутывают физический смысл вихревых процес-

сов и усложняют численные методы решения задач.

Для избавления от этих недостатков существует подход, основанный на исключении давления из уравнений за счет применения операции ротор (rot) к уравнению для вектора скорости, и переходу описания скоростей через функцию тока и вихрь. Такой переход повышает порядок дифференциального уравнения по пространственным переменным, но за счет введения функции тока дает автоматическое удовлетворение уравнения неразрывности, при этом система уравнений усложняется и обычно решается раздельным - двухполевым методом. Приравнивание операции ротора скорости к вихрю (которые по определению тождественно равны) не позволяет вычислить функцию тока через другую независимую функцию, что приводит к новым такого же порядка трудностям, связанным с заданием граничных условий для вихря [10].

Поэтому возникает необходимость в построении новой математической модели, позволяющей избежать указанных затруднений. Для вычисления функции тока и затем определения через нее вектора скорости можно использовать независимую векторную функцию - генератор вихрей. Похожий метод использовался некоторыми исследователями для решения задач в пористых средах, но только для двумерных уравнений [9].

В 1999 г. автором был разработан математический аппарат трехмерной Вихревой Гидродинамики (ВГ) для пористых сред, который затем был обобщен и для вязких сред. Предлагаемый новый математический аппарат основан на расщеплении конвекции на две составляющие - свободную, вызванную градиентом плотности и записанную через векторную функцию тока Т - векторный потенциал, и вынужденную, вызванную источниками, записанную через скалярную функцию U - скалярный потенциал, что позволяет раздельно вычислить их значения, а затем получить параметры суммарного конвективного потока.

На базе этого математического аппарата автором разработан пакет программ на языке Fortran, с помощью которого можно решать широкий круг трехмерных задач тепломассо-переноса в пористых и вязких средах на персональных компьютерах.

За последние годы были получены решения для ряда теоретических и практических задач, которые опубликованы в цен-

тральных журналах РАН, доложены на Российских и международных конференциях.

В данной работе сделан вывод системы уравнений вихревой гидродинамики и дано обоснование предлагаемого подхода. На анализе результатов решения задач тепломассо-переноса в пористых и вязких средах показаны возможности информационной технологии, разработанной на основе (ВГ) и сделаны обобщения по некоторым физическим процессам и явлениям, распространенным в различных геосферах Земли.

Обоснование нового подхода

и вывод уравнений (ВГ)

Из теории поля известно, что любое векторное поле M можно разложить на две составляющие. M = Mi + M2. Поле Mi = grad U - источниковое или дивергентное поле, где U - произвольная функция, имеющая первую производную, и выполняется равенство rot grad U = 0. Функцию U называют скалярным потенциалом векторного поля M. Тогда, div Mi = div grad U = M, ^ V2U = AU = M, где k1 - константа пропорциональности, 5 -плотность источников.

Поле M2 = rot A - вихревое или циркуляционное поле, тогда, div rot A = 0. Функцию A называют векторным потенциалом векторного поля M. Существует бесконечное число векторных функций A, описывающих одно и то же поле M, такие функции могут отличаться друг от друга с точностью до градиента некоторой скалярной функции ф, так как rot grad ф = 0. Это приводит к следующему соотношению: rot rot A = k2W ^ grad div A -V2A = k2W. Так как существует бесконечное множество функций A, описывающих одно и то же поле M, то можно выбрать из них такое, которое удовлетворяет равенству div A = 0 (калибровка Кулона) ^ V2A = AA = k2W, где k2 - константа пропорциональности, W - вектор плотности вихревых линий (линий тока), а плотность (густота) линий векторного поля пропорциональна его интенсивности.

Все процессы, протекающие в природе, в той или иной степени связаны с гидродинамическими потоками, вызванными перераспределением вещества в гравитационном поле. Известно, что генератором свободной конвекции является неоднородность плотности среды, а интенсивность конвекции пропорциональна градиенту этой плотности [1]. Значит, векторную функцию W можно

представить в виде градиента скалярной функции, например р (плотность среды), тогда W = gradp. Свободно конвективные потоки вещества, вызванные градиентом плотности, формирующим гидродинамический диполь, удобно записать с помощью векторной функции тока ^ - векторный потенциал, а потоки вынужденной конвекции, вызванные источниками, через функцию U - скалярный потенциал (давление, напор).

Тогда перенос среды (вещества) можно записать в виде векторного уравнения Пуассона AU = k1 5, описывающего источниковое или дивергентное поле, записанное через скалярную функцию U и A^ = k2Vp - описывающее вихревое или циркуляционное поле, записанное через векторную функцию которое в декартовой (ортогональной) системе координат распадается на три скалярных уравнения инвариантных относительно системы координат:

A^x = k2 5p/5x, = k2 dp%,

A^z = k25p/5z. (1)

В общем случае векторное поле скоростей, вызванное источниками и генераторами вихрей, можно записать так:

V = Vi + V2 = grad U + rot (2)

Это равенство представляет собой одну из форм теоремы Гельмгольца.

На этом основании автором предложена новая трехмерная система уравнений, позволяющая разделить конвекцию на две составляющие - свободную и вынужденную - и, вычислив их значения, получить параметры суммарного конвективного потока. Эту систему можно записать так:

= k2 Vp; (3)

AU = k15; (4)

V = grad U + rot (5)

Система уравнений (3-5) представляет собой систему уравнений Вихревой Гидродинамики (ВГ), записанную в векторной форме.

В общем случае плотность может зависеть от потенциальной функции p=F(U), т.е. от внешнего воздействия, но это достаточно узкий класс задач. В поле гравитации в основном формируются свободно-конвективные потоки, и плотность зависит от других скалярных функций - температуры и концентрации вещества p=F(T,C), для описания полей которых вводятся дополнительные уравнения.

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

изводных по сравнению с системой уравнении Навье-Стокса, но позволяет более точно вычислить конвективную составляющую V, которая играет определяющую роль в задачах гидродинамики и вычислении полей температуры и концентрации вещества, в задачах теп-ломассопереноса. Для решения задач со сво-бодноконвективными течениями уравнение (4) можно исключить, и порядок системы уравнений понизится. Кроме того, система (ВГ) является классической системой уравнений в частных производных эллиптического типа и позволяет более эффективно использовать вычислительные методы для ее решения.

При выводе уравнения для векторного потенциала использовалась калибровка Кулона, но существует другой подход, в котором используется калибровка Лоренца (его часто применяют при выводе уравнения для векторного потенциала в электродинамике). Если использовать калибровку Лоренца, то в уравнении появляется дополнительный член, а при нем - коэффициент (волновое число), который для задач гидродинамики сложно определить и, скорее всего, он должен быть разным для каждой задачи, что усложнит их решение. В дальнейшем, возможно, для каких-то задач гидродинамики такой подход может быть использован.

Данную систему уравнений также можно использовать для решения задач электромагнитной динамики, через векторный и скалярный потенциалы определять напряженность электрического или магнитного поля, используя их коэффициенты сопротивления или проводимости. Так как все напряжения возникают за счет неоднородности среды, то через перемещения (мгновенные скорости) можно определить поле деформаций, дилата-цию, девиатор (тензор напряжений) и таким образом решать задачи прочности. В будущем данная система уравнений будет апробирована на этих задачах.

Математическая модель тепломассопереноса в пористых и вязких средах

Добавив к полученной системе уравнений вихревой гидродинамики (ВГ) уравнения теплопроводности, массопереноса, уравнение состояния, граничные и начальные условия получим следующую полную сопряженную систему уравнений тепломассопереноса в пористых и вязких средах.

Эта модель в терминах функции тока, давления, температуры и концентрации выглядит так (рис. 3):

Рис. 3. Система уравнений вихревой гидродинамики ВГ (пояснения в тексте)

ство компонент; ^ - время; х, у, ъ - координаты,

Здесь ¥ 1,2,3 - проекция вектора функции тока на координатные оси; Т - температура; и -

потенциальная функция; - концентрация

/-го компонента в системе, / =1, т, т - количе-

Ух,Уу,У2;№х;№у;№ъ,ах,ау,аъ,Лх,Лу,Лъ,Вх,Ву,Въ,а

х,ау,аъ, - проекции векторов скоростей конвекции (фильтрации) и тепловой волны, коэффициентов текучести (фильтрации), гидро-

динамического сопротивления, гидродисперсии, температуропроводности на координатные оси; Д, Д*, р, р0, Рф - термический и концентрационный коэффициенты объемного расширения, коэффициент сжимаемости, исходная и текущая плотность среды (флюида); Т0,С0,и0,Т,С,и - фоновые и текущие температура, концентрация, и потенциальная функция (давление - Р или напор - Н); О -обобщенный коэффициент дисперсии (диффузии), или коэффициент гидродисперсии. К

- обобщенная функция кинетики, превращения веществ. х, у, I) - функция источников или генераторов, формирующих внешние (сторонние) поля, например транзитный поток или внутренние поля, сформированные в предыдущие моменты времени 1=1

В уравнение (11) для задач с пористыми средами используется обобщенный коэффициент дисперсии, или коэффициент гидродисперсии О = Ом + Ок, где Бм - коэффициент молекулярной диффузии, Ок -коэффициент диффузии или механической дисперсии:

Бк =ЛУ\ У= VI %, (14) V' - вектор действительной скорости фильтрации, Я - параметр механической дисперсии, X - пористость. Для вектора скорости тепловой волны можно записать следующее выражение:

С

W = V , (15)

С

Сп

где V - вектор скорости фильтрации, Сф, Сп

- объемная теплоемкость флюида и флюидо-насыщенных пород.

К - обобщенная функция кинетики превращения веществ.

К =

а1к"

др1

дг

+ 81 + N1 .

Здесь: а1к - стехиометрический коэффициент для /-го элемента в к-ой реакции осаждения-

др]

растворения минералов.

дг

- прибыль или

убыль у-го минерала за счет реакций растворения - осаждения, замещения.

др- = ±Гх( К, П р/ п СЯ - п р/ п сЯ* > (17)

дг ^ к * 1 г ц

где Сч - молярная концентрация элементов в

поровом водном растворе в моль/л, q = Р - молярная плотность минералов в единице объема в моль/дц3, у'=1^ п, Кк - константа равновесия к-й реакции растворения, осаждения г,^-го минералов 8,ге [1,..._)], 3 -стехиометрический коэффициент для s-го или г-го минерала в к-ой реакции. 1 < s,r < п, Я - стехиометрический коэффициент 1-го или q-го элемента раствора в к-й реакции. 1 < 1,q < т 81 - прибыль или убыль /-го элемента в растворе за счет равновесных и обменных реакций, / = 1^т.

(18)

^=±Х[к г п сЯ'г -п с?

Г V 1 ц У

где Кг - константа равновесия г-ой реакции диссоциации или растворения осаждения 1-го вещества. Я г1 - стехиометрический коэффициент для 1-го или q-го элемента в г-ой реакции: 1 < < т,

N1 - прибыль-убыль /-го элемента за счет сорбции-десорбции или радиоактивного распада:

N1 = а1 С1 , (19)

где а1 - коэффициент сорбции-десорбции при бесконечной емкости, или радиоактивного распада >го элемента. С1 - текущая концентрация >го элемента.

Количество твердого вещества так же, как и растворенного, рассчитывается в моль/дц3. Пористость X в точке определяется по формуле

(20)

х=1 -Ер,е

(16)

У

где QJ - мольный объем j-го минерала дм3/моль, р - молярная плотность j-го минерала моль/дм3.

Пористость учитывается при образовании твердых осадков реакций, осаждение минералов прекращается, если пористость достигает заданного минимума, и при расчете коэффициента дисперсии на каждом новом временном шаге. Расчет уравнений химической кинетики, записанный на основе закона действующих масс, проводится итерационным методом с использованием случайных чисел. Для решения выписанной системы уравнений задаются граничные и начальные условия:

y,2) г=0 = р1(х у,2);

^ г =Ф ;

т\ г = а; т(х, у, 2)

г=0

= К2 (х y, 2);

C Г - R;

с(Х y, ZJ,-0 - F3 (х, y, z)

Р| Г -Н; Р(х, y, z) ,-0 - F4 (х, y, z), (21)

где, Ф, G, R, H, F1-4 - известные функции.

Это система уравнений шестого и более порядка, в зависимости от количества компонентов (m). Первые четыре уравнения (6, 7, 8, 9) описывают гидродинамику, причем первые три записанные через функцию тока - свободную конвекцию, вызываемую внутренними силами (градиентом плотности), а четвертое, записанное через скалярный потенциал, -вынужденную конвекцию, вызываемую любыми внешними силами. Пятое уравнение -теплопроводности (10), шестое - массопере-носа (11). Прямая связь между гидродинамическим блоком и уравнениями тепломассопе-реноса осуществляется через скорости конвекции (фильтрации) уравнение (12) - обратная через плотность среды (флюида) уравнение состояния (13). Для решения данной системы задаются начальные и граничные условия (21). Любые начальные состояния и внешние воздействия (транзитные потоки и т.д.) или взаимодействия разных веществ между собой и их фазовые переходы задаются с помощью функций источников, записанных в правой части для каждого уравнения (14-20).

Представленная модель эллиптическо-параболического типа, является классической системой уравнений, для численного решения которой подходят большинство разработанных и используемых расчетных методов и схем. Повышенный порядок системы уравнений практически не влияет на методику и технологию разработки программного кода. Это позволило автору на основе экономичных конечно-разностных схем, основанных на расщеплении уравнений по координатам осям и физическим процессам, используя абсолютно устойчивый локально-одномерный метод решения [11], разработать пакет программ на языке Fortran, с помощью которого можно решать широкий круг трехмерных задач теп-ломассопереноса в пористых и вязких средах на персональных компьютерах.

При решении задач используются три безразмерных комплекса. Первые два это тепловое и концентрационное числа Фурье Fot = ax/L2, Foc = Dt/ L2. В уравнениях, описывающих свободную и вынужденную конвекцию, возникающую за счет градиента плотности среды (флюида) или внешних сил, коэффициент пропорциональности k2 = (g/^)x(p/L) =

g/vL определяет плотность потока линий тока, характеризующий его интенсивность, в общем случае является тензором второго ранга и имеет размерность 1/(м2сек). Физический смысл этого коэффициента можно интерпретировать как текучесть, а его обратное значение к = 1/к2, - как гидродинамическое сопротивление (чем больше сопротивление среды, тем больше плотность источников). Поэтому третий безразмерный комплекс можно назвать критерием плотности распределения линий тока вихревой гидродинамики (ВГ), Плт = (g/ц,)x(p/L)xL2т. В уравнениях этих трех комплексов а, Б, g - коэффициенты температуропроводности, диффузии, динамической и кинематической вязкости, ускорения свободного падения, т, L2, L - характерное время, площадь и размер в решаемой задаче.

Основную нагрузку при интерпретации огромного объема результатов вычислений несет пакет графических программ визуализации, разработанный автором и позволяющий в динамике просматривать значения любой функции, вычисленной при решении поставленной задачи.

Численная реализация

Для решения системы уравнений, описывающей процессы конвективного тепло-массопереноса в трехмерных областях, необходимо использовать численные методы, позволяющие в каждом узле дискретной сетки, покрывающей модельную область, рассчитать значения искомых функций. Для численной реализации краевых задач используется абсолютно устойчивый локально-одномерный метод решения [11]. Суть его заключается в расщеплении уравнения по координатам и последующем решении одномерных уравнений методом прогонки, в сумме являющихся решением многомерного уравнения. Эта неявная схема решения имеет суммарную аппроксимацию порядка о(т + % 2 ), где т - шаг по времени, % - шаг сетки. Для аппроксимации конвективных членов в уравнениях (10, 11) используется схема с направленными разностями [12]. Кроме того, эти уравнения расщепляются по процессам, т. е. в каждый момент времени вначале учитывается вклад только кондуктивного и дисперсионного членов, рассчитанных по неявной схеме, а затем к ним добавляются значения конвективных составляющих, рассчитанных по явной схеме.

При решении уравнений Пуассона (6-9) применяется метод установления с нижней релаксацией [10]. Процесс считается установившимся при условии

8 <еЧт, (22)

где 8 - максимальное по всем узлам сетки значение невязки, - максимальное по модулю значение функции тока или потенциала, £ = 0,01 - 0,001. Количество итераций для установления процесса зависит от шага по времени и параметров среды. При малом шаге по времени Т и/или однородной структуре среды процесс устанавливается за 1-2 итерации, при большем шаге Т и/или неоднородной структуре среды или при относительно больших коэффициентах проницаемости (текучести) в начальные моменты времени итераций требуется больше, но затем количество итераций быстро уменьшается и стабилизируется в пределах 3-5.

Схему решения краевой задачи на каждом шаге по времени можно представить в следующем порядке:

1. По распределению температуры Тп-1,

Сп-1

5 и если есть источники, функции ип-1, согласно уравнению (13), определяется плотность - рп

2. По уравнениям (6-8) рассчитываются

функции тока ¥п при известной правой части.

3. Определяется распределение потенциальной функции ип, при наличии не нулевых граничных условий или источников по уравнению (9).

4. По функциям тока и потенциалу вычисляются скорости конвекции (фильтрации)

V", V", V". Для пористых сред по скоростям фильтрации вычисляются скорости тепловой волны Ж", Ж", Ж", коэффициенты гидродинамической дисперсии О", О", О" и действительные скоро-

сти фильтрации V х V у, V

5. Используя скорости конвекции или тепловой волны (для пористых сред), определяется температура т .

6. Используя скорости конвекции и диффузии или коэффициенты гидродинамической дисперсии и действительные скорости фильтрации (для пористых сред), определяется концентрация растворенных веществ С".

7. Если в системе учитывается кинетика превращения веществ или взаимодействия вода-порода, сорбция-десорбция (для пористых сред), то рассчитывается функция источника для уравнения (11).

Этот цикл повторяется для всех последующих моментов времени с использованием результатов решения уравнений Пуассона как начальных условий с предыдущего временного слоя.

Начальные и граничные условия

Для скалярных функций T,C,U - задаются стандартные начальные условия, значения функций в точке. Для проекций функции тока задаются нулевые начальные условия

I о = 0, или, если существует транзитный поток, вызванный какими-то внешними силами, значения для любой из этих функций - в виде числовой последовательности с заданным шагом для формирования потока с определённой скоростью в рассматриваемой области, например, горизонтальный поток вдоль координаты x, 10 = K*Ax*N, где К - коэффициент пропорциональности, Ax - шаг по координате x, N - количество точек разбиения. Граничные условия первого, второго, третьего или четвертого рода задаются для скалярных функций, хотя при выносе границ за пределы области влияния процессов их можно свести к условиям первого и второго рода. Для скоростей на гранях расчетной области используются условия, записанные через проекции функции тока. Для условий непротекания и проскальзывания г = 0. Для нормальных скоростей условие протекания, а для касательных (тангенциальных) скоростей условие прилипания S^/Sn | г = 0.

Условия применения

Модель реализована в виде пакета программ на языке Fortran для ПК типа IBM. Трансляция проводится с помощью компилятора Fortran PowerStation, позволяющего получать оптимизированные 32-разрядные приложения, работающие в операционных системах DOS и Windows.

Время расчета задачи определяется количеством точек в расчетной области и количеством моментов времени, на которое ведется расчет. Так, для области, разбитой сеткой на 45000 точек (50x30x30), при решении задачи в пористых средах на 7300 моментов времени (что составляет 20 лет, если шаг по

времени - 1 сутки) время счета на AMD-2000 составляло 150 минут (2,5 часа). При использовании кинетики взаимодействия вода-порода время расчета резко возрастает и намного увеличивается объем результатов вычислений.

В задачах с высоковязкими средами, для области разбитой сеткой на 72000 точек (45x40x40) и с корректировкой вязкости в зависимости от температуры и концентрации среды, при расчете на 16000 моментов времени с шагом 100 - 1000 лет требуется 90 минут (1,5 часа) расчетного времени на персональном компьютере Pentium G620-2.6GHz. Это вполне приемлемо для задач такого класса. При увеличении или уменьшении количества точек или моментов времени время счета пропорционально увеличивается или уменьшается.

Результаты расчета по каждой функции и заданного момента времени выводятся в файлы прямого доступа, затем просматриваются с помощью графических программ, разработанных автором, позволяющих наблюдать изменение вычисленных функций в динамике. При каждом расчете в файл прямого доступа выводятся 20-30 срезов по времени для каждой функции, что в сумме составляет 75-150 Mb дисковой памяти. Тоновые рисунки результатов расчетов представлены в виде

блок-диаграмм и (или) горизонтальных срезов по вертикальной координате, вертикальных по горизонтальным для определенного момента времени. Разметка координат - в метрах или километрах. Стрелками показано поле скоростей фильтрации (конвекции) в пяти диапазонах, выделенных цветом. Для каждого диапазона написано максимальное значение скорости, соответствующее максимальной длине стрелки на рисунке.

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

Анализ результатов

С использованием системы уравнений (ВГ), был решен и опубликован в Российских и зарубежных изданиях ряд задач для пористых, высоковязких и низковязких сред, результаты некоторых из них проанализированы в данной работе.

Рис. 4. Слева - динамика подъема плюма на три момента времени с распределением концентрации вещества среды (белым цветом выделен расплав) и поля скоростей конвекции (максимальная скорость стрелки желтого цвета): а) через 24 млн лет, Утах = 0,28 м/год; б) через 32 млн лет, Утшх = 0,42 м/год; с) через 48 млн лет, Утшх = 0,8 м/год, после начала подъема плюма. Верхняя часть блоков срезана по подошве литосферы. Справа - распределение: а) температуры ТтПп = 200°С, Ттах = 3100°С, и поля скоростей конвекции Утах = 0,8 м/год (желтый цвет стрелок); б) концентрации вещества среды (белым цветом выделен расплав), в линзе расплава у подошвы литосферы через 100 млн лет

Для высоковязких сред анализировались задачи по оценке условий формирования и эволюции нижнемантийного плюма (рис. 4) и условий формирования и развития диапиров в системе литосфера - кора (рис. 5) [4, 5, 6].

Для пористой среды задачи о становлении интрузива как одного из механизмов формирования условий для образования кольцевых и вихревых структур [3] (рис. 6).

Рис 5. Распределение температуры при различной вязкости окружающих пород и плотности вещества диапира при относительно высокой его вязкости (1-3 порядка ниже вязкости окружающих пород) после подъема диапира; а) вязкость окружающих пород равна 1023 П-с, Ттах = 1730С; б) вязкость окружающих пород в десять раз меньше, Ттах = 1690 С; с) вещество диапира претерпевало фазовые переходы; д) вещество диапира не претерпевало фазовые переходы

Рис 6. Схематичный блоковый разрез области расчета: а) распределение температуры Ттах = 1000 °С, Тт/п = 5 С, концентрации вещества и поля скоростей фильтрации флюида Утах = 0,003 м/год в начальный момент времени после внедрения интрузива в пористую среду верхней коры; б) распределение вынесенного потоком флюида из интрузива вещества и поле скоростей фильтрации через 200 000 лет, срез сделан на глубине 5 000м, Утах = 0,007 м/год

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

льду и газовых выбросов в летний период из озера Байкал. [7, 8] (рис. 7).

Рис 7. Формирование кольцевых структур на льду озера Байкал: а) распределение температуры и поля скоростей конвекции при образовании линзы у подошвы ледяного покрова; б) распределение плотности воды при формировании кольцевой структуры подо льдом озера Байкал (лед отрисован белым цветом)

В этих задачах моделировалось влияние источника тепла и/или химического вещества, снижающих плотность мантийного (плюм) и литосферо-корового (диапир) материала для сплошной среды, или флюидов для пористой среды. Для гидросферы задавался источник с повышенной температурой или пониженной плотностью (газ).

Во всех случаях шел одинаковый процесс формирования восходящего и нисходящего потоков (струй), который состоял из нескольких этапов развития. Над заданным источником появлялась куполообразная структура с пониженной плотностью и вязкостью.

Для каждого координатного направления формировался гидродинамический диполь и стратификация становилась неустойчивой (неустойчивость Рэлея-Тейлора). За счет горизонтального градиента плотности возникали выталкивающие силы (силы Архимеда) и куполообразная структура всплывала, формируя восходящий конвективный поток торообразного вида с восходящей ветвью в центре структуры (канале) и нисходящей в окружающем массиве (среде).

При больших размерах источника и бо-

лее низкой вязкости среды возможно образование мульти-диполя и появление нескольких симметрично расположенных куполообразных структур (рис. 5,б). Скорость подъема потока зависела от его плавучести и вязкости для сплошной среды и проницаемости среды и вязкости флюида для пористой. За счет вертикального градиента плотности появлялись тангенциальные силы, вызывающие вращение структуры.

В центральной и головной частях поток вращался против часовой стрелки (циклонический вихрь), а на периферии, где шло погружение, и в подошве структуры -по часовой стрелке (антициклонический вихрь) (рис. 8).

При достижении потоком более вязкой среды (преграды) подъем его прекращался, вещество потока растекалось по его подошве, формируя линзу с вращением в центральной ее части против часовой стрелки, а на периферии - по часовой стрелке. Если преграда отсутствовала, то поток мог внедриться в менее плотную среду (гидросферу, атмосферу) и сформировать кольцеобразную структуру (рис. 6,б, 7,б).

Рис. 8. Распределение плотности среды и поля скоростей при подъеме конвективного потока (струи) для трех моментов времени развития процесса. В верхней и нижней областях рисунка показаны горизонтальные срезы в головной (для каждого) и подошвенной (для всех) моментов времени подъема струи. Максимальная скорость потока - желтый малиновый цвет в центре струи, минимальная - белый цвет на периферии

Анализ системы уравнений и результатов, полученных при решении этих задач, позволил сделать некоторые обобщения по причинам и следствиям динамики процессов, развивающихся в различных геосистемах.

Выводы

На основе этого анализа и ранее опубликованных материалов можно сделать следующие выводы. Всеми процессами и явлениями в нашем мире управляет гравитация. Это основной источник движения и преобразования вещества в природе. Гравитационное поле имеет не только скалярный потенциал (нормальную составляющую), но и векторный (вихревую составляющую) rot A Ф 0, а генератором вихрей является Vp - градиент плотности. Вертикальный градиент плотности вызывает тангенциальные силы, приводящие различные среды к вращению в горизонтальной плоскости, а горизонтальный градиент плотности

вызывает их вертикальные движения (подъем, погружение). Это новое фундаментальное свойство гравитационного поля.

Причиной появления (действия) силы Кориолиса является не вращение Земли, а тангенциальная составляющая гравитационного поля. Это можно подтвердить с помощью простого примера. В ясную морозную погоду, когда влажность воздуха минимальна, а значит его плотность максимальная, дым, поднимаясь из печной трубы, буквально вкручивается в воздух. Если на следующий день, при повышении температуры, а значит и влажности воздуха, плотности их становятся близкими, то дым стелется по земле, успевая только чуть подняться над трубой. Это говорит о том, что причиной такого поведения дыма является разность плотности дыма и воздуха, чем она больше, тем с большей скоростью и выше поднимается дым из трубы (рис. 9).

Рис. 9. Подъем дыма из печной трубы: а) - ясная морозная погода, влажность воздуха минимальна, его плотность максимальная рд << рв; б) - при повышение температуры, а значит и влажности воздуха, плотности их близки рд ~рв

Во всех пограничных слоях разных сред (оболочек Земли) существует скачок плотности, который активизирует (генерирует) тангенциальные силы гравитационного поля, а те в свою очередь формируют вихревые потоки. Менее плотные верхние слои плывут по нижним слоям с вращением по часовой стрелке, формируя антициклонический вихрь. Движение этих потоков осложняется рельефом нижнего, более плотного, слоя. Такие потоки - придонные (приповерхностные) течения - существуют в атмосфере, реках, озерах, морях и океанах.

За счет этих же сил верхняя твердая оболочка Земли (литосфера) разбита на плиты, а те в свою очередь - на более мелкие

блоки верхней более твердой коры. Все эти плиты и блоки движутся относительно друг друга, вращаясь на более плотном основании и/или перемещаясь в вертикальном направлении, при подъеме формируя материки, а при погружении - озера, моря и океаны. Так как эти плиты и блоки плотно упакованы, то никаких горизонтальных перемещений у них быть не может, к тому же материки, моря и океаны являются тонкой пленкой на поверхности литосферных плит и составляют 1,5-2 % от ее мощности-толщины. Иллюзию же перемещения создает движение в виде погружения одного и подъема другого с обоюдным вращением двух соседних блоков (рис. 10).

Рис. 10. Современные сейсмичность, вулканизм и границы плит: материк и океан могут располагаться на одной плите. Желтые точки - очаги землетрясений, красные треугольники - активные наземные вулканы

В северном полушарии Земли большая часть коры испытывает вздымание (подъем) с вращением против часовой стрелки (циклонический вихрь) поэтому основная часть материков находится в этом полушарии. В южном полушарии большая часть коры погружа-

ется с вращением по часовой стрелке (антициклонический вихрь) и покрыта океанами. Результирующая тангенциальная сила направлена с запада на восток, она приводит планету к вращению в том же направлении, т.е. с запада на восток (рис. 11).

Рис. 11. Распределение материков и океанов на поверхности планеты Земля. Слева геоид - уровневая поверхность потенциала силы тяжести, имеет форму груши вытянутой к северному полюсу, желтый цвет - максимальная гравитация, синий - минимальная. Справа - глобус

Более того, собственные колебания (их еще называют свободные) всех тел и объектов на нашей планете являются результатом действия вихревой составляющей гравитационного поля, так как между поверхностью опоры, на которую они опираются или к которой крепятся, существует градиент плотности, приводящий к вращающему моменту, возбуждающего колебания. Поэтому собственные колебания тел, объектов и систем - результат действия тангенциальных сил гравитационного поля.

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

Список литературы

1. Гебхарт Б., Джалурия Й., Махаджан Р., Саммакия Б. Свободноконвективные течения, тепло- и массообмен: в 2 кн. / пер. с англ. М.: Мир, 1991. С. 678.

2. Глуховский М.В. Кольцевые структуры юго-востока Сибири и их возможная природа // Геотектоника. 1978. № 4. С. 50-63.

3. Гунин В.И. Становление интрузива как возможный механизм формирования условий для образования кольцевых и вихревых структур (численный эксперимент) // Известия вузов. Геология и разведка. 2006. № 4. С. 40-44.

4. Гунин В.И. Оценка условий формирования и эволюции нижнемантийного плюма на

основе численного эксперимента // Геодинамика формирования подвижных поясов земли: матер. Междунар. науч. конф. Екатеринбург, 24-26 апреля 2007. С 75-78.

5. Гунин В.И. Оценка условий формирования-развития диапиров в литосфере и земной коре на основе численного эксперимента // Граниты и эволюция Земли: I Междунар. конф. Улан-Удэ, 26-29 августа 2008.С. 97-101.

6. Гунин В.И. Новая информационная технология и ее возможности при моделировании геосистем // Геодинамика и тектонофизика. 2011. Т. 2, № 4. С. 356-377. URL: http://gt.crust.irk.ru (дата обращения: 12.11.2017).

7. Гунин В.И. Оценка причин и условий формирования кольцевых структур на льду Байкала по результатам численного эксперимента // Подземная гидросфера: матер. Всерос. совещания по подземным водам востока России. Иркутск: Изд-во ООО "Географ", 2012. С. 513-517.

8. Гунин В.И. Оценка воздействия газогидротермальной деятельности Байкальского рифта на акваторию озера по результатам численного эксперимента // Геодинамика и тектонофизика. 2014. Т. 5, № 3. С. 763-775. URL: http://gt.crust.irk.ru (дата обращения: 20.11.2017).

9. Кирюхин А.В., Сугробов В.М. Модели теп-лопереноса в гидротермальных системах Камчатки. М.: Наука, 1987. 150 с.

10. Тарунин Е.Л. Численный эксперимент в задачах свободной конвекции, Иркутск: Изд-во Иркут. ун-та, 1990. 223 с. (41).

11. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1975. 735 с.

12. Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. 429 с.

Vortex hydrodynamics: a new approach to modeling geosystems

V. I. Gunin

Center for Modeling Geosystems "MoGeos"; 38-19, 50 let Oktyabrya st., Ulan-Ude, 670034, Russia vigunin@list.ru

Using a mathematical apparatus based on the system of Navier-Stokes equations, it is difficult to describe vortex (circular) patterns formed under the action of gravitational field in all geospheres of our planet. The paper presents a new mathematical apparatus of vortex fluid dynamics (CFD) software, which allows us to adequately describe these patterns and get a deeper understanding of their essence. Based on the analysis of the proposed system of equations and the results of a number of model problems, generalizations are made on some physical processes and phenomena prevalent in different geospheres of the Earth.

Keywords: gravity; vector potential; convection; vortex structures; mathematical model; numerical experiment; heat and mass transfer; mantle; lithosphere; crust; diaper; rift.