CC BY
19
www.volsu.ru
DOI: https://doi.Org/10.15688/mpcm.jvolsu.2021.3.5
УДК 519.63, 532.5 ББК 26.222
Дата поступления статьи: 30.05.2021 Дата принятия статьи: 16.07.2021
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ САМОСОГЛАСОВАННОЙ ДИНАМИКИ ПОВЕРХНОСТНЫХ И ГРУНТОВЫХ ВОД
1
Сергей Сергеевич Храпов
см о см
и
Ü
m о с я о,
X
©
Кандидат физико-математических наук,
доцент кафедры информационных систем и компьютерного моделирования,
Волгоградский государственный университет
khrapov@volsu.ru
https://orcid.org/0000-0003-2660-2491
просп. Университетский, 100, 400062 г. Волгоград, Российская Федерация
Аннотация. Построена математическая и численная модели совместной динамики поверхностных и грунтовых вод, в которых учитываются нелинейная динамика жидкости, впитывание воды с поверхности в грунт, фильтрационные течения в грунте и высачивание воды из грунта обратно на поверхность. Динамика поверхностных вод описывается уравнениями Сен-Венана с учетом пространственно-неоднородных распределений рельефа местности, коэффициентов придонного трения и инфильтрации, а также нестационарных источников и стоков воды. Для численного интегрирования уравнений Сен-Венана применяется хорошо апробированный CSPH-TVD метод второго порядка точности, параллельный CUDA-алгоритм которого реализован в виде программного комплекса «EcoGIS-Simulation» для высокопроизводительных вычислений на суперкомпьютерах с графическими сопроцессорами (GPU). Динамика грунтовых вод описывается нелинейным уравнением Буссенеска, обобщенным на случай пространственно-неоднородного распределения параметров пористой среды и поверхности водоупора (границы между водопроницаемым и слабопроницаемым грунтами). Численное решение этого уравнения строится на основе конечно-разностной схемы второго порядка точности, CUDA-алгоритм которой интегрирован в расчетный модуль «EcoGIS-Simulation» и согласован с основными этапами CSPH-TVD метода. Относительное отклонение численного решения от точного решения нелинейного уравнения Буссенеска не превышает 10-4-10-5. В работе проводится сравнение результатов численного моделирования динамики грунтовых вод с аналитическими решениями линеаризованного уравнения Буссенеска, используемыми в качестве расчетных формул
в методиках прогноза уровня грунтовых вод в окрестности водных объектов. Показано, что погрешность этих методик составляет несколько процентов даже для простейшего случая плоскопараллельного течения грунтовых вод с постоянным подпором. На основе полученных результатов сделан вывод, что предложенный в работе метод численного моделирования совместной динамики поверхностных и грунтовых вод может являться более универсальным и эффективным (обладает существенно лучшей точностью и производительностью) по сравнению с существующими методиками расчета зон подтопления, особенно для гидродинамических течений со сложной геометрией и нелинейного взаимодействия встречных потоков жидкости, возникающих в период сезонных паводков при затоплении обширных территорий суши.
Ключевые слова: гидродинамика поверхностных и грунтовых вод, модель мелкой воды, пористые среды, фильтрационные течения, численное моделирование, CSPH-TVD-метод, параллельные вычисления, CUDA-алгоритм, суперкомпьютеры с GPU.
Введение
Подземные воды имеют важное значение при решении практических задач, связанных с определением гидрологических режимов пойменных территорий междуречий и речных долин, водохранилищ и озер, а также других природных и гидротехнических водных объектов. К настоящему времени динамика грунтовых вод и фильтрационных течений хорошо изучена как на основе экспериментальных данных, так и с использованием хорошо известных и апробированных математических моделей (с результами можно ознакомиться, например, в работах [1;5;6; 10; 12; 13; 16; 17; 19;22;25]). На основе этих моделей разработаны методики (см., например, [18]) по расчету стационарных и нестационарных уровней грунтовых вод, а также расходов и скорости фильтрационных течений, используемых при проектировании гидротехнических сооружений и определения границ зон подтопления в окрестности водных объектов для задач кадастрового учета. Как правило, аналитические формулы, используемые в этих методиках, основаны на решениях линеаризованных уравнений, являющихся некоторыми приближениями исходной нелинейной модели динамики грунтовых вод, что, в свою очередь, приводит к увеличению погрешности расчетов, ограничению класса решаемых задач по геометрии водных объектов и структуре течения, а также требует больших временных затрат при выборе параметров моделей и нормировочных коэффициентов.
В настоящее время важным и мощным инструментом изучения гидродинамических течений в природных и технических системах наряду с физическим экспериментом и аналитическими методами исследования математических моделей является вычислительный эксперимент, основанный на численном решении уравнений гидродинамики хорошо апробированными методами и применении параллельных технологий для повышения производительности расчетов [3;4; 11; 14; 21; 24; 26].
Целью данной публикации является разработка эффективной методики расчета зон затопления и подтопления на сложном рельефе местности с произвольной геометрией нестационарных водных объектов, основанной на прямом численном интегрировании уравнений нелинейной динамики поверхностных и грунтовых вод. В части 1 описана математическая модель динамики поверхностных и грунтовых вод, а в части 2 —
численный алгоритм и структура параллельного вычислительного модуля программного комплекса «EcoGIS-Simulatюn». Результаты тестовых расчетов в одномерной и двумерной постановках представлены в части 3.
1. Математическая модель 1.1. Динамика поверхностных вод
Моделирование динамики поверхностных вод осуществляется с учетом всех основных физических и метеорологических факторов: источников воды; рельефа суши и дна водоемов; свойств подстилающей поверхности — придонного трения, инфильтрации (впитывание воды в грунт); внутреннего вязкого трения и ветрового воздействия; вращения Земли — силы Кориолиса; испарения воды. Эти факторы позволяет учитывать модель мелкой воды, основанная на системе уравнений Сен-Венана:
дН
— + У± (Ни) = д, (1)
+ У± (Ни 0 и) = —дНУ±ц + Hf , (2)
где Ъ(х, у) — функция рельефа местности; Н(х, у, ¿) = ц(х, у, ¿) — Ъ(х, у) — толщина слоя поверхностной воды; п — уровень свободной водной поверхности; и(х,у,{) = {их,иу} — вектор скорости течения поверхностного слоя жидкости; д — ускорение свободного падения; У± = {д/дх,д/ду} — дифференциальный оператор набла в плоскости (х,у). Величина д = д^ — д(гп^) — д(еи) + д^ определяет скорость притока/оттока жидкости за счет действия источников/стоков в поверхностном слое воды, где д( 3\х,у,1) — скорость притока воды через гидротехнические сооружения и за счет осадков; д(гп^\х,у,Ь) — функция стока из-за просачивания (инфильтрации) в грунт; д— темп потери воды из-за испарения; д^ир) — скорость высачивания (подъема) грунтовых вод на поверхность. В общем случае для суммарной удельной силы, действующей на слой жидкости в модели мелкой воды, имеем:
? = и + Ь + ^ + С + Ъ, , (3)
где fл — сила придонного трения; — сила вязкого турбулентного трения; ^ — сила Кориолиса; ^ — сила воздействия атмосферных потоков (ветра); ^ = дщ/Н — сила воздействия источников/стоков воды д, втекающих/вытекающих в поверхностный слой жидкости со скоростью щ. Более подробное описание удельных сил из (3), действующих на поверхностные воды и используемых в данной работе, дано в [20; 21; 27].
1.2. Метод расчета динамики грунтовых вод
Важной характеристикой, используемой при описании динамики грунтовых вод, является пористость грунта ф, которая определяет объем пор в единице объема грунта, то есть ф = 1 — у8, где = — относительный объем твердых частиц в единице
объема грунта; да — плотность твердых частиц; д^ — плотностью скелета (сухого) грунта, численно равная массе единицы объема грунта без учета воды в его порах. Наряду с величиной ф используют и коэффициент пористости кр = ф/(1 — ф). Уровень грунтовых вод будем характеризовать величиной цд(х,у,Ь) = Нд(х,у,Ь) + Ьд(х,у), где
Нд(х,у, ¿) — толщина слоя грунтовых вод над рельефом слабопроницаемых грунтов
Ъд (х, у).
В общем трехмерном случае для описания динамики грунтовых вод используется математическая модель движения жидкости в пористой среде, которая, в свою очередь, основана на уравнениях механики многофазных сред [5; 13]:
+ Ъ(кгдг-уг) = 0, (4)
д
+ V (а^íVí 0 V;) = -<х^Рг + авig + ^ + аiцiAvi, (5)
где а — объемное содержание 2-й фазы с плотностью вг и давлением р¿; g = {0, 0, —д} —
Ц Vj — Vi
вектор ускорения силы тяжести; = сс^—-----сила межфазного взаимодействия
Рг 1г
между 2-й и ^'-й фазой; ц — динамическая вязкость г-компоненты среды; I — характерный размер области 2-й фазы, который для пористой среды соответствует характерному размеру пор; вi — некоторый коэффициент, смысл которого выяснится ниже.
Если ограничиться рассмотрением модели многофазной среды, состоящей из двух компонент г = 1 — жидкость (вода) и г = 2 — неподвижный недеформируемый пористый скелет грунта, то V! = V, V2 = 0, а = ф, а2 = 1 — а = У«, в1 = в, Р1 = Р, Ц1 = Ц, 11 = I, в1 = в. Тогда для такой двухкомпонентной среды из системы (4)-(5) можно получить уравнения, описывающие динамику жидкости в пористой среде:
д|в1 + V(фвv) = 0, (6)
^ + V (фвv 0 V) = —фVp + фвg — ФЦV + Ф^. (7)
В пределе Ф ^ 1 из (7) получаем стандартное уравнение движения в приближении однокомпонентной гидродинамики, так как в этом случае размер пор I ^ го.
Поскольку для фильтрационных потоков в пористых грунтах характерная скорость течения мала (|V| ~ 0.01-1 м/сут), то в уравнении (7) можно пренебречь левой частью по сравнению с правой. Кроме того, из-за малости характерного размера пор I можно также пренебречь и последним членом в правой части (7) по сравнению с силой межфазного взаимодействия. С учетом этих предположений из уравнения (7) получим выражение для скорости жидкости в пористой среде:
V = — ^р— вg). (8)
Ц
По определению скорость фильтрации ий связана со скоростью жидкости (скоростью частиц жидкости) в пористой среде V следующим соотношением: ий = •V. Учитывая эту связь и выражение (8), для скорости фильтрации будем иметь:
к
и9 =--^р — вg), (9)
Ц
где к = "фр1 = — коэфф ициент проницаемости, зависящий только от свойств
09
мористого скелета, а кф — коэффициент фильтрации, который определяется экспериментально. Соотношение (9) — хорошо известный закон Дарси, связывающий скорость фильтрации жидкости через мористую среду с напором или градиентом давления.
Таким образом, если определить скорость жидкости в уравнении непрерывности (6) с учетом (8) и (9), то получим уравнение дифф узионного (параболического) типа относительно плотности д:
d (4 в) =v
dt V
— (VP - Qg)
(10)
В общем случае коэфф ициент фильтряции1чвляетс.я функцией пространственных координат кф = кф(х,у, z), а для замыкания уравнения (10) используется уравнение состояния р = р( д).
Так как на практике при построении границ зон затопления и подтопления в окрестности населенных пунктов характерные пространственные и временные масштабы задачи в горизонтальном направлении существенно больше соответствующих масштабов вдоль вертикальной координаты, то для описания динамики поверхностных и грунтовых вод можно воспользоваться приближением мелкой воды (тонкого слоя жидкости). Проводя процедуру усреднения уравнения (10) по z-координате, аналогичную случаю поверхностных вод, то есть интегрируя (10) по z в пределах от Ьд до пд = Нд + Ьд в предположении гидростатического равновесия тонкого слоя грунтовых вод р(z) = дд(цд — z), получим следующее уравнение:
-Н
4= v± (кфНдV^g)+ Qg , (11)
которое описывает эволюцию толщины Нд тонкого слоя жидкости в грунте за счет медленных фильтрационных течений (по закону Дарси), определяемых коэффициентом i i
кф, и действующих в слое грунта источников/стоков воды qg = q^^) + — qg— (gv\ где величины q(inf) и qi^ определяют скорость притока жидкости с поверхности и за счет подземных источников, соответственно, а q^^ и qtgap"> — скорость оттока воды из пористого слоя за счет инфильтрации в слябопрониняемый глубинный слой грунта и высачивания на поверхность (подъема грунтовых вод), соответственно.
В пределе однородных грунтов (кф,4 = const) и плоского дна водоупора (Ьд = = const) уравнение (11) сводится к уравнению Буссинеска, описывающему безнапорное движение подземных вод со свободной поверхностью (см., например, [6; 19]), а аналитические решения этого уравнения, получаемые после процедуры его линеаризации, используются в методиках определения нестационарных уровней грунтовых вод в окрестности водных объектов при повышении или спаде уровня поверхностных вод (см., например, [18]). После линеаризации нелинейное уравнение Буссинеска переходит в уравнение, являющееся аналогом уравнения теплопроводности Фурье:
f = а + ^ , (12)
-t 4
где а = кфН3/4, Д± — дифференциальный оператор лапласиан в плоскости, перпендикулярной оси Oz; U = НдН8 при первом способе линеаризации, а при втором —
U = Щ/2; Hs = ß(Hmin + Hmax)/2 — средняя глубина потока грунтовых вод; ß — коэффициент, который зависит от отношения Hm;n/Hmax и способа линеаризации (см., например, [19, с. 45]); Hmin и Hmax — минимальная и максимальная глубина потока, соответственно.
Аналитические формулы, являющиеся решением уравнения (12) и используемые в методиках расчета нестационарных границ зон подтопления местности грунтовыми водами, сильно зависят как от геометрии задачи, так и от начальных и граничных условий. Например, для плоскопараллельного течения с начальным (бытовым) уровнем грунтовых вод Hg(x,t = 0) = Hmin и граничными условиями вида Hg(х = 0, t) = Hmax, dH
lim g = 0 решение (12) имеет вид:
х^ж дх
U =U + (U2 -Ul)I e-e2dd , (13)
Vя l
X
X
где X = —; U\ = HminHs; U2 = HmaxHs для первого способа линеаризации, а для 2\>at
второго — Ui = Hm.J2, U2 = HlaJ2.
Самосогласованная динамика поверхностных и грунтовых вод определяется параметрами впитывания воды в грунт q(m?) и высачивания на поверхность q^, а также гидравлической связью поверхностного и подземного потоков. Различают две основные стадии фильтрационных течений [6; 19]: 1) свободная фильтрация, при которой отсутствует гидравлическая связь между поверхностными и грунтовыми водами; 2) фильтрация с подпором, при которой фильтрационный и грунтовый потоки объединяются и начинают взаимодействовать друг с другом. Свободная фильтрация наблюдается на начальной (первой) стадии промачивания грунтов под водными объектами, а глубина промачивания 4 ограничена снизу фронтом впитывания, который находится выше уровня грунтовых вод ng. Как только фронт впитывания достигает уровня грунтовых вод, начинается вторая стадия — фильтрация с подпором. Для оценки скорости движения фронта впитывания w(4) в водопроницаемом слое грунта, подстилающем дно водоема с глубиной H, применяют следующее уравнение [19]:
- = f = | + 1) ■ (14)
где Hk — высота капиллярного вакуума. Уравнение (14) следует из уравнения (8) на
др Ар H + Hk
z-компоненту скорости в предположении, что —- ~ —— =-. Интегрируя уравне-
dz Az 4
ние (14) при H + Hk = const, можно оценить время Т4, за которое фронт промачивания достигнет уровня грунтовых вод, залегающих на глубине 4max = Ь — ng:
и
Т4 = (H + Hk)[А — ln(1 + A)] , (15)
где А = 4max/(H + Hk).
Как видно из (14), при 4 ^ 0 скорость впитывания w ^ то, что является следствием линейного приближения, используемого при выводе уравнения (8). В реальных пористых средах за счет нелинейных эффектов и вязкости на глубинах промачивания
£ ~ Нк скорость впитывания ограничена сверху величиной птах, которую в рамках простейшей модели можно оценить исходя из того же уравнения (14) следующим образом:
Птах ~ ф + Н~к.
Свободная фильтрация, в свою очередь, также делится на два типа фильтрационных течений [6; 19] — это так называемые напорная фильтрация и инфильтрация. Если водоем расположен в водопроницаемом грунте, то реализуется случай напорной фильтрации, а при наличии на дне водоема так называемого «экрана», представляющего собой тонкий (покровной) слабопроницаемый слой грунта с коэффициентом фильтрации кс ^ кф и толщиной т, возможны оба типа фильтрационных течений в зависимости от соотношения скоростей фильтрации в этих слоях. С учетом (14) и соотношения между скоростью фильтрации и скоростью жидкости в пористой среде можно получить условия возникновения различных типов фильтрационных течений [19]. В случае кф/кс < Ш происходит напорная фильтрация с полным заполнением пор грунта, а при кф/кс > Ш реализуется инфильтрация с частичным заполнением пор грунта, где Ш = (Н + Нк + т)/т.
Учитывая тот факт, что направление фильтрационных течений имеет наряду с основной вертикальной составляющей и горизонтальную, обусловленную как процессом рассеяния жидкости на порах грунта, так и горизонтальным градиентом при на-
порной фильтрации, а также рассматривая характерные времена процессов Ь > т£, можно ограничиться моделью грунтовых вод со свободной поверхностью (11), в которой величина определяется в зависимости от типа фильтрации.
Для случая напорной фильтрации можно записать:
(гпЛ Г фп(£) Н> 0 и Цд <b, (16)
1 0, иначе,
где значение £ определяется в каждый момент времени посредством численного интегрирования уравнения (14).
При инфильтрации будем иметь:
М) Г кФ, Н> 0 и <b, (17)
\ 0, иначе.
Кроме того, при фильтрации с подпором, то есть на второй стадии, когда уровень грунтовых вод достигает дна водоема (цд = Ь), в уравнении (11) вместо величины цд необходимо использовать у\д = Пд + Н = п.
Величину определим следующим образом:
^ { 0, Нд < 0,
Jinf) I к*, Нд > 0, Уо 1 с\ и ^ с\
где к* — коэффициент фильтрации в слабопроницаемом глубинном слое грунта.
Высачивание воды из грунта обратно на поверхность можно учесть посредством вычисления в каждый момент времени величины А Нд = цд — Ь и корректировки толщины
!д. Тогда для величины д^
слоя жидкости на поверхности Н ив грунте Нд. Тогда для величины получим:
4А на
q(uv) = J —, А н > (ig)
А Нд < 0,
где т — временной шаг в численной модели.
2. Численная схема и программная реализация
Для численного моделирования самосогласованной динамики поверхностных и грунтовых вод использовались СБРИ-ТУО-метод [15; 21] и параллельный код («СиОА-DSW»), подробно описанный в работах [24; 26] и являющийся базовой версией вычислительного модуля программного комплекса «Есо61Б-Б1ти1а1юп», который был дополнен модулем расчета фильтрационных течений в слое грунта («CUDA-DGW»).
2.1. Численный алгоритм
Для численной аппроксимации уравнения (11) проведем стандартную процедуру дискретизации на пространственно-временной сетке (х,у, t) ^ (Xi, yj, tn) (xi+i = Xi + + h, yj+\ = Xj + h, tn+\ = tn + тп, где h — размер ячеек пространственной сетки, тп — временной шаг), заменив непрерывные функции их значениями в узлах этой сетки f(x,y, t) ^ ffj, а их производные — конечно-разностными аналогами. В результате для уравнения (11) с учетом (17) и (19) получим явную консервативную численную схему второго порядка точности по пространственной координате ~ 0(h2):
ЦП+1 _ Tjn | Хп
gi,3 gi,3 +
кф i+l/2,jHgi+i/2,j СЛзг+lJ )
— кф i-l/2,jHgi-i/2,j Wgi,j — ngi-l,j) + кф ij+V2^flij+l/2Cnflij+l — ) — (20)
1 Qn
— h Нп (-глп — -глп \ ' Q
кф i,j-\/2пд i,j-\/2 (ngi,j пд i,j-1)
I !дг,з
4
г,з
где значения функций с дробными пространственными индексами (г ± 1/2, j ± 1/2)
, Mj)±i + Mj)
определяются на границах ячеек как j)±i/2 = — —^- .
Отметим, что представленная численная схема (20) имеет первый порядок точности по времени ~ 0(т), но допускают повышение порядка либо за счет применения явно-неявного алгоритма Кранка — Никольсона [23], либо при использовании процедуры предиктор - корректор (так называемый алгоритм leapfrog) как в CSPH-TVD-методе [21; 27]. Кроме того, для повышения порядка аппроксимации как по времени, так и по пространственным координатам можно использовать численные схемы высокого порядка точности (см., например, [8; 14]).
В общем случае условие устойчивости численного алгоритма для моделирования самосогласованной динамики поверхностных и грунтовых вод можно представить в виде:
К min (T<ydro), x<yrunt, (21)
где 0 < К < 1 — число Куранта; T<ydr0 = min ( —-,--,—^ ) — величи-
г,3 у 2|иП1 К,! +
на временного шага в CSPH-TVD-методе при моделировании динамики поверхностных
вод (см. [21; 27]); x<n'runt) = min ( ^ — | — характерное диффузионное время в
%,з \2кфi,jHnitjJ
численной модели динамики грунтовых вод (20).
2.2. Параллельный вычислительный модуль
При разработке расчетного модуля с рабочим названием «EcoGIS-Simulation», предназначенного для моделирования самосогласованной динамики поверхностных и грунтовых вод («CUDA-DSW> + «CUDA-DGW»), использовалась технология параллельных вычислений CUDA для графических процессоров (GPU). Кроме того, при проведении расчетов на компьютерных платформах с несколькими GPU применялась технология GPUdirect (обмен данными между GPU напрямую по шине PCI-E или NVLINK) совместно с технологией параллельных вычислений OpenMP для систем с общей памятью на центральных процессорах (CPU). Реализация расчетного модуля «EcoGIS-Simulation» на основе гибридных параллельных технологий CUDA-GPUdirect-OpenMP позволяет эффективно использовать компьютерные платформы с несколькими CPU и GPU (пх xCPU+mxGPU) и повысить производительность вычислений в тысячи раз по сравнению с последовательной версией на CPU.
В качестве языка программирования был выбран CUDA C [2], являющийся расширением языка C для программирования графических процессоров. На рисунке 1 показана структура вычислительного модуля «EcoGIS-Simulation» в виде потоковой диаграммы, демонстрирующей последовательность выполнения следующих CUDA-модулей:
• «DSW» — вычисление временного шага (21) и расчет динамики поверхностных вод (1)—(3) на основе CSPH-TVD-метода, подробное описание численного алгоритма представлено в работах [24;26];
• «DGW» — расчет динамики грунтовых вод (11)—(19) на основе численного алгоритма (20).
Рис. 1. Потоковая диаграмма для расчетного модуля «EcoGIS-Simulation». Показаны последовательность выполнения блоков CUDA-ядер на GPU и потоки данных между CPU и GPU
Как видно из рисунка 1, для повышения порядка аппроксимации используется метод предиктор - корректор («leapfrog»), подробно описанный для численной схемы CSPH-TVD в работах [21; 27].
3. Результаты численного моделирования
Величины Ь, Ьд, Н, Нд, п, Ла и пространственные координаты (х,у) будем измерять в метрах (м), а время t в сутках (сут). Для водопроницаемого слоя грунта (песка) примем следующие значения коэффициента фильтрации и пористости: кф = 5 м/сут,
ф = 0, 4. Будем считать, что слабопроницаемый глубинный слой грунта (водоупор) не пропускает воду, то есть к* = 0.
3.1. Одномерные расчеты плоскопараллельных течений грунтовых вод
Рассмотрим сначала задачу о динамике грунтовых вод в одномерном приближении с расчетной областью х € [0,£] (Ь = 100 м) и количеством расчетных ячеек Мх = 1000. Построим численные решения уравнения (11) на основе численной схемы (20) с применением алгоритма предиктор - корректор для двух тестовых задач с различными начальными и граничными условиями:
Тест 1
Нд(х, 0) = 0, Нд(0,1) = Ятах = 2, Нд(Ь, 0) = Ят1п = 0. (22)
Тест 2
Нд (х, 0) = Нтш = 1, Нд (0,1) = Ятах = 2, Нд (Ь, 0) = Ятт = 1. (23)
На рисунке 2 представлены пространственные распределения величины пд(х,{) = = Нд(х,£)+Ъд (Ьд = -1) в различные моменты времени. Полученные численные решения соответствуют двум видам автомодельных решений нелинейного уравнения в частных производных параболического типа (нелинейного уравнения теплопроводности или диффузии) [7; 9], так как в уравнении (11) коэффициент кфНд перед градиентом уровня грунтовых вод приводит к появлению квадратичной нелинейности (~ Щ) и является аналогом коэффициента теплопроводности ае ~ Т (где Т — температура). Для начальных и краевых условий (22) получаем решение в виде нелинейной «тепловой» волны (волны подтопления), распространяющейся в пустоте (по сухому водоупору) и имеющей характерный приблизительно параболический профиль с отличной от нуля производной дцд/дС на фронте волны при Нд(Со) = 0, где С = х — Ь — автомодельная переменная, а Vf — скорость фронта волны. В случае (23) имеем существенно более гладкий приблизительно экспоненциальный профиль пд (х,Ь) в окрестности фронта волны (£0) с
Иш —д = 0 и Иш Нд(С) = Ят;п > 0. Как видно из рисунка 2, скорость распространяю о С с^со
нения волны подтопления (повышения уровня грунтовых вод) в задаче Тест 2 почти в два раза превышает аналогичную скорость для задачи Тест 1.
Отметим, что сравнение численных решений Тест 1 и Тест 2 с точными автомодельными решениями уравнения (11) (см., например, [7; 9]) показывает, что разработанный численный алгоритм моделирования динамики грунтовых вод имеет ожидаемый второй порядок точности ~ 0(х2,к2), а для представленных в работе расчетов относительная погрешность составляет ~ 10-4-10-5.
На рисунке 3 показано сравнение численного решения нелинейного уравнения (11) с аналитическими решениями (13) линеаризованного уравнения (12), используемыми в методиках по расчету границ зон подтопления в окрестности водоемов и населенных пунктов в период паводка. В соответствии с [19] для первого способа линеаризации при определении величины Н3 выбрано значение в ~ 1,075, а для второго — в ~ 0, 85. Как видно из рисунка 3, второй тип аналитических решений лучше соответствует численному, особенно в окрестности подпора поверхностными водами (х < 25 при £ =10 и
х < 50 при £ = 40), где относительная погрешность составляет менее 0,5 %, а в дальней зоне (х > 25 при £ = 10 их > 50 при £ = 40) относительное отклонение не превышает 1,6 %. Для первого типа решения (13) относительная погрешность оказывается больше и составляет ~ 2-3 % как в окрестности подпора, так и в дальней зоне.
Рис. 2. Динамика грунтовых вод в одномерной (плоскопараллельной) модели (11) при заданных начальных и граничных условиях: (22) — Тест 1; (23) — Тест 2
Рис. 3. Пространственное распределение глубины грунтовых вод (Тест 2) в различные моменты времени. Сплошной линией 1 показаны численные решения уравнения (11), штриховой линией 2 — аналитическое решение (13) при первом способе линеаризации (см. уравнение (12)), а пунктирной линией 3 — аналитическое решение (13) при втором способе линеаризации
3.2. Двумерные расчеты совместной динамики поверхностных и грунтовых вод
Теперь рассмотрим задачу о самосогласованной динамике поверхностных и грунтовых вод на модельном рельефе Тест 3, решая совместно уравнения (1)-(2) и (11) на двумерной расчетной сетке с х,у £ [-L/2,L/2] и количеством расчетных ячеек Nx х Ny = 1024 х 1024. Модельный рельеф поверхности суши и водоупорного слоя грунта зададим в виде функций:
3пт v2
b(r) = -2 cos— + 2.5 + 0.5 exp(-r2/80), bg (г) = —, (24)
где г = у/х2 + у2 — расстояние от оси симметрии, измеряемое в метрах (м).
Результаты численного моделирования самосогласованной динамики поверхностных и грунтовых вод представлены на рисунке 4. В начальный момент времени (t = 0) параболический колодец в центре расчетной области наполнен водой, а в грунте (песке) заданы два типа начальных условий: 1) сухой водоупор Hg(x,t = 0) = 0; 2) постоянный уровень грунтовых вод ng(x,t = 0) = 0.4 м. Уровень воды в колодце ц(х, 0) = 4 м, а максимальная глубина 3 м. Параметры грунта такие же, как и в тестовых расчетах (Тест 1 и Тест 2). В качестве модели впитывания поверхностных вод в грунт выбран режим инфильтрации (17).
На рисунке 4 хорошо видны процессы инфильтрации воды с поверхности в грунт, фильтрационные течения в грунте и высачивание воды из грунта обратно на поверхность. В процессе начальной стадии (0 < t < 1) инфильтрации воды происходит понижение уровня поверхностных вод и повышение уровня грунтовых вод в центре расчетной области (под колодцем), а также формируется нелинейная волна подтопления с профилем, подобным случаю Тест 1 (рис. 4a) и Тест 2 (рис. 4b), распространяющаяся в радиальном направлении от оси симметрии. На втором этапе (1 < t < 10) уровни поверхностных и грунтовых вод сравниваются, в дальнейшем происходит совместное понижение уровней воды за счет оттока грунтовых вод из центральных областей в процессе распространения волны подтопления. Третья стадия (10 < t < 100) характеризуется полным осушением колодца и процессом высачивания воды из грунта обратно на поверхность (подъема грунтовых вод) в окрестности осесимметричного рва (г ~ 33), окружающего центральный колодец. На заключительной стадии (t > 100) численное решение стремится к стационарному профилю уровня воды n = ng = const.
Отметим, что из-за кривизны водоупора, цилиндрической геометрии и нестационарности уровня поверхностных вод скорость распространения волны подтопления в задаче Тест 3 оказывается меньше, чем в плоскопараллельной постановке Тест 1 и Тест 2 с постоянным значением подпора (Hmax).
Выводы
В заключение сформулируем основные результаты работы и обсудим области их возможного практического применения.
Построенная математическая модель самосогласованной динамики поверхностных и грунтовых вод позволяет рассматривать гидродинамические течения в окрестности водоемов произвольной геометрии с учетом пространственной неоднородности параметров водопроницаемых грунтов, различные стадии и типы фильтрационных течений (свободная, с подпором, напорная и инфильтрация), а также процессы высачивания воды из
грунта на поверхность. Модель основана на уравнениях гиперболического типа (уравнениях Сен-Венана), описывающих динамику поверхностных вод в приближении теории мелкой воды, и уравнении параболического типа (уравнении Буссинеска), описывающего нелинейную динамику грунтовых вод со свободной поверхностью. На основе математической модели разработаны численный метод второго порядка точности и параллельный CUDA-алгоритм, реализованный в виде программного комплекса для суперкомпьютеров с GPU.
Рис. 4. Самосогласованная динамика поверхностных и грунтовых вод на модельном осесимметричном рельефе для различных начальных состояний: — сухой водоупор Нд(ж,£ = 0) = 0; (Ь) — уровень грунтовых вод пд(ж, £ = 0) = 0.4. На врезке в левом верхнем углу показана 3D модель рельефа
Показано, что аналитические решения линеаризованного уравнения Буссинеска, которые используются в рекомендациях и методиках по расчету границ зон подтопления в окрестности водных объектов в период паводков [6; 18; 19], имеют погрешность в несколько процентов по сравнению с точными решениями даже в самом тривиальном случае плоскопараллельного течения в однородном грунте и требуют подбора коэффициентов модели в зависимости как от структуры начального течения, так и от
параметров подпора поверхностными водами. С учетом неоднородности реальных грунтов, существенно более сложной геометрии водных объектов, нестационарности исходных начальных состояний и нелинейных эффектов, возникающих при взаимодействии гидродинамических потоков поверхностных и грунтовых вод, погрешность такого упрощенного линейного подхода может существенно возрасти и привести к некорректным (неадекватным) оценкам уровней грунтовых вод и границ зон подтопления.
Таким образом, методики, основанные на аналитических решениях линеаризованного уравнения Буссинеска в общем случае, то есть для реальной местности с более сложной геометрией водных объектов и пространственной неоднородности параметров подстилающего слоя грунтов, являются существенно менее точными и более затратными по сравнению с предложенной в работе универсальной методикой численного моделировании самосогласованной динамики поверхностных и грунтовых вод. А применение в наших численных моделях высокопроизводительных параллельных вычислений на суперкомпьютерах с GPU позволяет эффективно рассчитывать динамику зон затопления и подтопления на обширных территориях с высоким пространственным разрешением.
Отметим, что для определения точности и пределов применимости рассматриваемой в работе модели самосогласованной динамики поверхностных и грунтовых вод требуется проведение как натурных экспериментов с реальными грунтами, так и модельных численных расчетов на основе полной трехмерной гидродинамики пористых сред.
Благодарности. Автор признателен А.В. Хоперскову за обсуждение результатов работы и полезные замечания.
ПРИМЕЧАНИЕ
1 Работа выполнена при финансовой поддержке Министерства образования и науки РФ (госзадание № 0633-2020-0003).
СПИСОК ЛИТЕРАТУРЫ
1. Барренблатт, Г. И. Теория нестационарной фильтрации жидкости и газа / Г. И. Бар-ренблатт, В. М. Ентов, В. М. Рыжик. — М. : Недра, 1972. — 288 c.
2. Боресков, А. В. Основы работы с технологией CUDA / А. В. Боресков, А. А. Харламов. — М. : ДМК Пресс, 2010. — 232 c.
3. Булатов, О. В. Регуляризованные уравнения мелкой воды и эффективный метод численного моделирования течений в неглубоких водоемах / О. В. Булатов, Т. Г. Елизарова // Журн. вычисл. матем. и матем. физ. — 2011. — Т. 51, № 1. — C. 170-184.
4. Булатов, О. В. Регуляризованные уравнения мелкой воды для численного моделирования течений с подвижной береговой линией / О. В. Булатов, Т. Г. Елизарова // Журн. вычисл. матем. и матем. физ. — 2016. — Т. 56, № 4. — C. 665-684.
5. Великанова, Ю. В. Гидромеханика многофазных сред / Ю. В. Великанова. — Самара : Изд-во Самар. гос. техн. ун-та, 2009. — 166 c.
6. Веригин, Н. Н. Расчет фильтрационных потерь из рыбохозяйственных водоемов / Н. Н. Веригин, С. В. Васильев, В. С. Саркисян. — М. : Пищевая промышленность, 1977. — 143 c.
7. Галактионов, В. А. Методы построения приближенных автомодельных решений нелинейных уравнений теплопроводности / В. А. Галактионов, А. А. Самарский // Матем. сб. — 1982. — Т. 118, № 3. — C. 291-322.
8. Григорян, Л. А. Параллельная реализация задачи транспорта веществ на основе схем повышенного порядка точности для уравнений диффузии-конвекции / Л. А. Григорян, А. А. Семенякина // Известия Южного федерального университета. Технические науки. —
2014. - № 12 (161). - C. 183-192.
9. Казаков, А. Л. О некоторых точных решениях нелинейного уравнения теплопроводности / А. Л. Казаков, Св. С. Орлов // Тр. ИММ УрО РАН. - 2016. - Т. 22, № 1. -C. 112-123.
10. Лейбензон, Л. С. Движение природных жидкостей и газов в пористой среде / Л. С. Лейбензон. - М. ; Л. : Гос. изд-во технико-теорет. лит., 1947. - 244 с.
11. Локально-двумерные схемы расщепления для параллельного решения трехмерной задачи транспорта взвешенного вещества / А. И. Сухинов, А. Е. Чистяков, В. В. Сидоря-кина, С. В. Проценко, А. М. Атаян // Математ. физика и компьютер. моделирование. -2021. - Т. 24, № 2. - C. 38-53. - DOI: https://doi.Org/10.15688/mpcm.jvolsu.2021.2.4.
12. Маскет, М. Течение однородных жидкостей в пористой среде / М. Маскет. - М. : Институт компьютерных исследований, 2004. - 640 с.
13. Механика насыщенных пористых сред / В. Н. Николаевский, К. С. Басниев, А. Т. Горбунов, Г. А. Зотов. - М. : Недра, 1970. - 339 с.
14. Параллельная реализация задач транспорта веществ и восстановления донной поверхности на основе схем повышенного порядка точности / А. И. Сухинов, А. Е. Чистяков, А. А. Семенякина, А. В. Никитина // Вычислительные методы и программирование: новые вычислительные технологии. - 2015. - Т. 16, № 2. - C. 256-267.
15. Писарев, А. В. Численная схема на основе комбинированного подхода SPH-TVD: проблема моделирования сдвиговых течений / А. В. Писарев, С. С. Храпов, А. В. Хопер-сков // Вестник Волгоградского государственного университета. Серия 1, Математика. Физика. - 2011. - № 2 (15). - C. 138-141.
16. Подземная гидравлика / К. С. Басниев, А. М. Власов, И. Н. Кочина, В. М. Максимов. - М. : Недра, 1986. - 303 с.
17. Полубаринова-Кочина, П. Я. Теория движения грунтовых вод / П. Я. Полубаринова-Кочина. - М. : Недра, 1977. - 664 с.
18. Пособие к СНиП 2.06.15-85 «Прогнозы подтопления и расчет дренажных систем на застраиваемых и застроенных территориях» / А. Ж. Муфтахов, И. В. Коринченко, Н. М. Григорьева, В. И. Сологаев, А. П. Шевчик. - М. : Стройиздат, 1989. - 407 с.
19. Фильтрация из водохранилищ и прудов / С. В. Васильев, Н. Н. Веригин, Г. А. Разумов, Б. С. Шержуков. - М. : Колос, 1975. - 304 с.
20. Численная модель динамики поверхностных вод в русле Волги: оценка коэффициента шероховатости / А. В. Писарев, С. С. Храпов, Е. О. Агафонникова, А. В. Хоперсков // Вестник Удмуртского университета. Математика. Механика. Компьютерные науки. -2013. - № 1. - C. 114-130.
21. Численная схема для моделирования динамики поверхностных вод на основе комбинированного SPH-TVD-подхода / С. С. Храпов, А. В. Хоперсков, Н. М. Кузьмин, А. В. Писарев, И. А. Кобелев // Вычислительные методы и программирование: новые вычислительные технологии. - 2011. - Т. 12, № 1. - C. 282-297.
22. Barrenblatt, G. E. Bask Co^epts in the Theory of Seepage of Homogeneous Liquids in Fissured ^ск^ / G. E. Barrenblatt, I. P. Zheltov, I. N. Ko^ina // Journal of Applied Mathematfcs. - 1960. - Vol. 24, № 5. - P. 852-864.
23. Crank, J. A pra^^al method for numerkal evaluation of solutions of partial differential equations of the heat ^ndu^^ type / J. Crank, P. Nkolson // Pro^ Camb. Phil. So^ -1947. - Vol. 43, № 1. - P. 50-67.
24. Dyakonova, T. Numeral Model of Shallow Water: The Use of NVIDIA CUDA Graphks Processors / T. Dyakonova, A. Khoperskov, S. Khrapov // Commun^^ns in Computer and Information Sderce. - 2016. - Vol. 687. - P. 132-145. - DOI: http://dx.doi.org/10.1007/978-3-319-55669-7_11.
25. Flow and Transport in Fradured Porous Media / P. Dietrkh, R. Helmig, M. Sauter, H. Hotzl, J. Kongeter, G. Teuts^. - Berlin : Springer-Verlag, 2005. - 488 p.
26. Khrapov, S. S. Applkation of Graphks Processing Units for Self-Consistent Modelling of Shallow Water Dynamks and Sediment Transport / S. S. Khrapov, A. V. Khoperskov
// Lobachevskii Journal of Mathematics. - 2020. - Vol. 41. - P. 1475-1484. - DOI: https://doi.org/10.1134/S1995080220080089.
27. The Numerical Simulation of Shallow Water: Estimation of the Roughness Coefficient on the Flood Stage / S. S. Khrapov, A. V. Pisarev, I. A. Kobelev, A. G. Zhumaliev, E. O. Agafonnikova, A. G. Losev, A. V. Khoperskov // Advances in Mechanical Engineering. — 2013. - Vol. 2013. - Article ID: 787016. - DOI: http://dx.doi.org/10.1155/2013/787016.
REFERENCES
1. Barrenblatt G.I., Entov V.M., Ryzhik V.M. Teoriya nestatsionarnoy filtratsii zhidkosti i gaza [Theory of Unsteady Filtration of Liquid and Gas]. Moscow, Nedra Publ., 1972. 288 p.
2. Boreskov A.V., Kharlamov A.A. Osnovy raboty s tekhnologiey CUDA [Basics of Working with CUDA Technology]. Moscow, DMK Press Publ., 2010. 232 p.
3. Bulatov O.V., Elizarova T.G. Regulyarizovannye uravneniya melkoy vody i effektivnyy metod chislennogo modelirovaniya techeniy v neglubokikh vodoemakh [Regularized Shallow Water Equations and an Efficient Method for Numerical Simulation of Shallow Water Flows]. Zhurn. vychisl. matem. i matem. fiz. [Computational Mathematics and Mathematical Physics], 2011, vol. 51, no. 1, pp. 170-184.
4. Bulatov O.V., Elizarova T.G. Regulyarizovannye uravneniya melkoy vody dlya chislennogo modelirovaniya techeniy s podvizhnoy beregovoy liniey [Regularized Shallow Water Equations for Numerical Simulation of Flows with a Moving Shoreline]. Zhurn. vychisl. matem. i matem. fiz. [Computational Mathematics and Mathematical Physics], 2016, vol. 56, no. 4, pp. 665-684.
5. Velikanova Yu.V. Gidromekhanika mnogofaznykh sred [Hydromechanics of Multiphase Media]. Samara, Izd-vo Samar. gos. tekhn. un-ta, 2009. 166 p.
6. Verigin N.N., Vasilyev S.V., Sarkisyan V.S. Raschet filtratsionnykh poter iz rybokhozyaystvennykh vodoemov [Calculation of Filtration Losses From Fishery Reservoirs]. Moscow, Pishchevaya promyshlennost Publ., 1977. 143 p.
7. Galaktionov V.A., Samarskiy A.A. Metody postroeniya priblizhennykh avtomodelnykh resheniy nelineynykh uravneniy teploprovodnosti [Methods for Constructing Approximate Self-Similar Solutions of Nonlinear Heat Conduction Equations]. Matem. sb. [Sbornik: Mathematics], 1982, vol. 118, no. 3, pp. 291-322.
8. Grigoryan L.A., Semenyakina A.A. Parallelnaya realizatsiya zadachi transporta veshchestv na osnove skhem povyshennogo poryadka tochnosti dlya uravneniy diffuzii-konvektsii [Parallel Implementation of the Problem of Transport of Substances Based on Schemes of Increased Order of Accuracy for the Equations of Diffusion-Convection]. Izvestiya Yuzhnogo federalnogo universiteta. Tekhnicheskie nauki, 2014, no. 12 (161), pp. 183-192.
9. Kazakov A.L., Orlov Sv.S. O nekotorykh tochnykh resheniyakh nelineynogo uravneniya teploprovodnosti [On Some Exact Solutions of the Nonlinear Heat Equation]. Tr. IMM UrO RAN, 2016, vol. 22, no. 1, pp. 112-123.
10. Leybenzon L.S. Dvizhenie prirodnykh zhidkostey i gazov v poristoy srede [The Movement of Natural Liquids and Gases in a Porous Medium]. Moscow, Leningrad, Gos. izd-vo tekhniko-teoret. lit., 1947. 244 p.
11. Sukhinov A.I., Chistyakov A.E., Sidoryakina V.V., Protsenko S.V., Atayan A.M. Lokalno-dvumernye skhemy rasshchepleniya dlya parallelnogo resheniya trekhmernoy zadachi transporta vzveshennogo veshchestva [Local Two-Dimensional Splitting Schemes for 3d Suspended Matter Transport Problem Parallel Solution]. Matemat. fizika i kompyuter. modelirovanie [Mathematical Physics and Computer Simulation], 2021, vol. 24, no. 2, pp. 38-53. DOI: https://doi.org/10.15688/mpcm.jvolsu.2021.2A
12. Masket M. Techenie odnorodnykh zhidkostey v poristoy srede [The Movement of Natural Liquids and Gases in a Porous Medium]. Moscow, Institut kompyuternykh issledovaniy Publ., 2004. 640 p.
13. Nikolaevskiy V.N., Basniev K.S., Gorbunov A.T., Zotov G.A. Mekhanika nasyshchennykh poristykh sred [Mechanics of Saturated Porous Media]. Moscow, Nedra Publ., 1970. 339 p.
14. Sukhinov A.I., Chistyakov A.E., Semenyakina A.A., Nikitina A.V. Parallelnaya realizatsiya zadach transporta veshchestv i vosstanovleniya donnoy poverkhnosti na osnove skhem povyshennogo poryadka tochnosti [Parallel Implementation of the Problems of Transport of Substances and Reconstruction of the Bottom Surface Based on Schemes of Higher Order of Accuracy]. Vychislitelnye metody i programmirovanie: novye vychislitelnye tekhnologii, 2015, vol. 16, no. 2, pp. 256-267.
15. Pisarev A.V., Khrapov S.S., Khoperskov A.V. Chislennaya skhema na osnove kombinirovannogo podkhoda SPH-TVD: problema modelirovaniya sdvigovykh techeniy [Numerical Scheme Based on the Combined SPH-TVD Approach: the Problem of Shear Flow Modeling]. Vestnik Volgogradskogo gosudarstvennogo universiteta. Seriya 1, Matematika. Fizika [Science Journal of Volgograd State University. Mathematics. Physics], 2011, no. 2 (15), pp. 138-141.
16. Basniev K.S., Vlasov A.M., Kochina I.N., Maksimov V.M. Podzemnaya gidravlika [Underground Hydraulics]. Moscow, Nedra Publ., 1986. 303 p.
17. Polubarinova-Kochina P.Ya. Teoriya dvizheniya gruntovykh vod [The Theory of Groundwater Movement]. Moscow, Nedra Publ., 1977. 664 p.
18. Muftakhov A.Zh., Korinchenko I.V., Grigoryeva N.M., Sologaev V.I., Shevchik A.P. Posobie k SNiP 2.06.15-85 «Prognozy podtopleniya i raschet drenazhnykh sistem na zastraivaemykh i zastroennykh territoriyakh» [''Flooding Forecasts and Calculation of Drainage Systems in Built-Up and Built-Up Areas"]. Moscow, Stroyizdat Publ., 1989. 407 p.
19. Vasilyev S.V., Verigin N.N., Razumov G.A., Sherzhukov B.S. Filtratsiya iz vodokhranilishch i prudov [Filtration from Reservoirs and Ponds]. Moscow, Kolos Publ., 1975. 304 p.
20. Pisarev A.V., Khrapov S.S., Agafonnikova E.O., Khoperskov A.V. Chislennaya model dinamiki poverkhnostnykh vod v rusle Volgi: otsenka koeffitsienta sherokhovatosti [Numerical Model of the Dynamics of Surface Waters in the Volga Channel: Estimation of the Roughness Coefficient]. Vestnik Udmurtskogo universiteta. Matematika. Mekhanika. Kompyuternye nauki, 2013, no. 1, pp. 114-130.
21. Khrapov S.S., Khoperskov A.V., Kuzmin N.M., Pisarev A.V., Kobelev I.A. Chislennaya skhema dlya modelirovaniya dinamiki poverkhnostnykh vod na osnove kombinirovannogo SPH-TVD-podkhoda [Numerical Scheme for Modeling the Dynamics of Surface Waters Based on the Combined SPH-TVD Approach]. Vychislitelnye metody i programmirovanie: novye vychislitelnye tekhnologii, 2011, vol. 12, no. 1, pp. 282-297.
22. Barrenblatt G.E., Zheltov I.P., Kochina I.N. Basic Concepts in the Theory of Seepage of Homogeneous Liquids in Fissured Rocks. Journal of Applied Mathematics, 1960, vol. 24, no. 5, pp. 852-864.
23. Crank J., Nicolson P. A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of the Heat Conduction Type. Proc. Camb. Phil. Soc., 1947, vol. 43, no. 1, pp. 50-67.
24. Dyakonova T., Khoperskov A., Khrapov S. Numerical Model of Shallow Water: The Use of NVIDIA CUDA Graphics Processors. Communications in Computer and Information Science, 2016, vol. 687, pp. 132-145. DOI: http://dx.doi.org/10.1007/978-3-319-55669-7_11.
25. Dietrich P., Helmig R., Sauter M., Hotzl H., Kongeter J., Teutsch G. Flow and Transport in Fractured Porous Media. Berlin, Springer-Verlag, 2005. 488 p.
26. Khrapov S.S., Khoperskov A.V. Application of Graphics Processing Units for Self-Consistent Modelling of Shallow Water Dynamics and Sediment Transport. Lobachevskii Journal of Mathematics, 2020, vol. 41, pp. 1475-1484. DOI: https://doi.org/10.1134/S1995080220080089.
27. Khrapov S.S., Pisarev A.V., Kobelev I.A., Zhumaliev A.G., Agafonnikova E.O., Losev A.G., Khoperskov A.V. The Numerical Simulation of Shallow Water: Estimation of the Roughness Coefficient on the Flood Stage. Advances in Mechanical Engineering, 2013, vol. 2013, article ID: 787016. DOI: http://dx.doi.org/10.1155/2013/787016.
NUMERICAL MODELING OF SELF-CONSISTENT DYNAMICS OF SHALLOW AND GROUND WATERS
Sergey S. Khrapov
Candidate of Physical and Mathematical Sciences, Associate Professor, Department of Information Systems and Computer Modeling, Volgograd State University khrapov@volsu.ru
https://orcid.org/0000-0003-2660-2491
Prosp. Universitetsky, 100, 400062 Volgograd, Russian Federation
Abstract. A mathematical and numerical model of the joint dynamics of shallow and ground waters has been built, which takes into account the nonlinear dynamics of a liquid, water absorption from the surface into the ground, filtration currents in the ground, and water seepage from the ground back to the surface. The dynamics of shallow waters is described by the Saint-Venant equations, taking into account the spatially inhomogeneous distributions of the terrain, the coefficients of bottom friction and infiltration, as well as non-stationary sources and flows of water. For the numerical integration of Saint-Venant's equations, the well-tested CSPH-TVD method of the second order of accuracy is used, the parallel CUDA algorithm of which is implemented as a software package "EcoGIS-Simulation" for high-performance computing on supercomputers with graphic coprocessors (GPU). The dynamics of groundwater is described by the nonlinear Bussensk equation, generalized to the case of a spatially inhomogeneous distribution of the parameters of the porous medium and the surface of the aquiclude (the boundary between water-permeable and low-permeable soils). The numerical solution of this equation is built on the basis of a finite-difference scheme of the second order of accuracy, the CUDA algorithm of which is integrated into the calculation module of the "EcoGIS-Simulation" software package and is consistent with the main stages of the CSPH-TVD method. The relative deviation of the numerical solution from the exact solution of the nonlinear Boussinesk equation does not exceed 10-4-10-5. The paper compares the results of numerical modeling of the dynamics of groundwater with analytical solutions of the linearized Bussensk equation used as calculation formulas in the methods for predicting the level of groundwater in the vicinity of water bodies. It is shown that the error of these methods is several percent even for the simplest case of a plane-parallel flow of groundwater with a constant backwater. Based on the results obtained, it was concluded that the proposed method for numerical modeling of the joint dynamics of surface and ground waters can be more versatile and efficient (it has significantly better accuracy and productivity) in comparison with the existing methods for calculating flooding zones, especially for hydrodynamic flows with complex geometry and nonlinear interaction of counter fluid flows arising during seasonal floods during flooding of vast land areas.
Key words: hydrodynamics of surface and ground waters, shallow water model, porous media, filtration flows, numerical modeling, CSPH-TVD method, parallel computing, CUDA algorithm, supercomputers with GPU.
