CC BY
53
УДК 551.462.8
М.М. Чумаков
Численное моделирование процесса размыва донного грунта в окрестности киля тороса
Постановка задачи и структура модели
В качестве базового инструментария для разработки моделей гидродинамических течений, возникающих при обтекании ледяных килей, в настоящей статье используется вычислительный конечно-элементный комплекс (ВК) ANSYS. В нем для расчетов скорости течения, изменения уровня моря и температуры воды используется система дифференциальных уравнений в частных производных (уравнения Навье-Стокса), описывающая движение жидкости, дополненная уравнениями теплопроводности, гидростатики и состояния. Важным достоинством использования этой системы уравнений является возможность воспроизведения реалистичной трехмерной структуры морских течений, что позволяет решать целый ряд инженерных задач, связанных с моделированием внешних воздействий на гидротехнические сооружения и донный грунт.
Для целей изучения процессов донной эрозии в области расположения килей ледяного образования (тороса) необходимо выполнить моделирование обтекания водой протяженного твердого тела, форма которого аппроксимирует подводную часть тороса, размещенного в достаточно широком плоском канале таким образом, что имеется зазор между его нижней поверхностью и дном канала. Ширина и длина канала зависят от геометрических размеров обтекаемого водой твердого тела, которые, по сути, задают линейный масштаб турбулентного движения жидкости. Ширина канала определяется максимальной шириной тороса и должна быть выбрана такой, чтобы турбулентное движение, возникающее у его стенок, не искажалось из-за воздействия граничных условий, задаваемых на стенках канала.
Трехмерное твердое тело, аппроксимирующее подводную часть тороса, было построено с помощью программы DesignModeler, которая является составной частью ВК ANSYS, на основании натурных данных о морфометрических свойствах торосов, полученных специалистами ААНИИ в ходе проведения экспедиционных работ на акватории Байдарацкой губы в 2010 г. Максимальная ширина тороса составила около 60 м, а длина - около 80 м. Максимальная глубина киля - 10,5 м.
При таких геометрических размерах подводной части тороса ширина канала была принята равной 300 м при длине 350 м. Глубина канала была выбрана равной 12 м, при этом необходимо отметить, что значение глубины является одним из тех параметров модели, которые определяют специфику моделирования динамики вод в том или ином географическом пункте.
На рис. 1 приведен общий изометрический вид снизу расчетной области модели, на котором «каркасными» блоками показаны регионы, имеющие различное пространственное разрешение конечно-элементной сетки, генерируемой в среде вычислительного комплекса ANSYS, а также изображена форма подводной части тороса.
Такое разбиение расчетной области модели на отдельные регионы определяется необходимостью обеспечить приемлемую точность расчетов при ограниченных компьютерных ресурсах. Учитывая, что целью моделирования является исследование воздействия течения на донный грунт в районе килей тороса, т.е. с точки зрения гидродинамики изучение распределения величины касательного напряжения на донной поверхности в области расположения тороса разрешение конечно-элементной сетки должно быть как можно более высоким. На некотором удалении от него допустимо более низкое пространственное разрешение, поскольку это не повлечет потери точности в расчете донного касательного напряжения под его килем.
Ключевые слова:
программный
комплекс ANSYS,
модель обтекания,
ледяное
образование,
тип донного грунта,
неньютоновская
жидкость,
эрозийный
потенциал,
локальный размыв.
Keywords:
ANSYS program complex, streamline model, ice formation, bottom ground type, non-Newtonian liquid,
erosive potential, local wash-out.
№ 3 (14) / 2013
134
Научно-технический сборник • ВЕСТИ ГАЗОВОЙ НАУКИ
Рис. 1. Изометрический вид расчетной области модели обтекания киля
ледяного образования
Для адекватного воспроизведения моделью турбулентных процессов в пристеночном слое расчетная сетка модели модифицируется с помощью применения инструмента среды вычислительного пакета ANSYS «сеточная инфляция», что обеспечивает плавное уменьшение размеров конечных элементов при приближении к стенке цилиндра, при этом сохраняя их гексаэдральную форму.
На рис. 2 приведен фрагмент расчетной области с нанесенной конечно-элементной сеткой, иллюстрирующий «сеточную инфляцию» у стенки киля тороса. В регионе расположения киля тороса среднее значение длины ребра гек-саэдральных элементов составляет 0,5 м.
При выбранном таким образом пространственном разрешении конечно-элементная сетка модели содержит 4324360 элементов.
Рис. 2. Фрагмент расчетной области с нанесенной конечно-элементной сеткой, иллюстрирующий «сеточную инфляцию» у стенки киля тороса
№ 3 (14) / 2013
Современные подходы и перспективные технологии в проектах освоения нефтегазовых месторождений российского шельфа
135
Необходимо отметить, что возможность применения той или иной расчетной сетки в среде ВК ANSYS определяется геометрической формой конечных элементов, которые ее формируют [1]. Для повышения точности численного решения методом конечных элементов дифференциальных уравнений, определяющих динамику жидкости, необходимо чтобы большинство конечных элементов расчетной сетки имело гексаэдральную форму. В принципе допускается использование элементов, имеющих форму тетраэдра, но это ведет к потере точно -сти, так как в этом случае используются аппроксимирующие функции формы более низкого порядка.
Однако независимо от того, какую форму имеют конечные элементы, для получения корректного с математической точки зрения решения необходимо, чтобы их форма была близка к равноугольной (90° для гексаэдра и 60° для тетраэдра). Значение коэффициента асимметрии, автоматически рассчитываемое в среде ВК ANSYS, определяет степень отклонения элементов модельной сетки от равноугольной формы и является критерием ее применимости для моделирования гидродинамических процессов.
В руководстве по использованию ВК ANSYS [1] указывается, что элемент, имеющий значение коэффициента асимметрии большее, чем 0,75, применяться не должен, использование элементов со значениями коэффициента асимметрии в диапазоне 0,75-0,5 дает удовлетворительный результат, в диапазоне 0,5-0,25 -гарантирует хороший результат, а диапазон от 0,25 до 0 соответствует высокому качеству элементов модельной сетки. Допускается наличие незначительного числа элементов с высокими значениями коэффициента асимметрии, но среднее значение для всей расчетной области должно быть меньше 0,5.
Для разработанной расчетной сетки модели обтекания киля ледяного образования среднее значение коэффициента асимметрии равно 0,2 при стандартном отклонении 0,12.
Приведенные данные позволяют сделать вывод о высоком качестве разработанной расчетной сетки, которое гарантирует получение в процессе интегрирования уравнений Навье-Стокса корректного с математической точки зрения описания поля течения.
Получение корректного с математической точки зрения результата моделирования поля
течения не гарантирует его адекватности с физической точки зрения. Реалистичность численных расчетов и, соответственно, возможность их применения для оценки воздействия течения на донный грунт зависит от правильного учета физических процессов, определяющих локальную динамику вод, главным из которых является турбулентность.
В настоящее время в практике инженерных расчетов используются различные модели турбулентности, причем ни одна из них не свободна от определенных недостатков [2]. В связи с этим при разработке модели, описывающей локальную динамику водной среды при ее взаимодействии с конкретным физическим объектом (например, килем торосистого образования), необходимо выбрать для применения ту модель турбулентности, которая при использовании на выбранной расчетной сетке и граничных условиях, отражающих специфику географической локализации рассматриваемого объекта, позволит получить реалистичные оценки распределения скорости течения при меньших затратах вычислительных ресурсов.
Учитывая результаты численных экспериментов, приведенных в [3], в разработанной модели была применена SST модель турбулентности.
Необходимо отметить, что вышеописанная конфигурация расчетной области модели обтекания киля ледяного образования была использована при проведении численных экспериментов по определению оценок эрозийного потенциала. Анализ полученных в результате выполненных расчетов полей скорости течения показал возможность модификации расчетной области модели в целях увеличения пространственного разрешения в области подводной части тороса при сохранении прежних требований к затрачиваемым компьютерным ресурсам. На рис. 3 изображена модифицированная расчетная область модели, примененной в численных экспериментах по определению зон размыва связного грунта под килем торосистого образования.
Ширина канала была принята равной 220 м при длине 200 м. В цилиндрическом регионе, содержащем торос и имеющем радиус 60 м, был применен автоматический метод построения расчетной сетки с максимальной длиной ребра элемента, равной 0,5 м. В других регионах сетка строилась методом Sweep, что позволило сгенерировать элементы в форме
№ 3 (14) / 2013
136
Научно-технический сборник • ВЕСТИ ГАЗОВОЙ НАУКИ
Рис. 3. Изометрический вид снизу модифицированной расчетной области модели обтекания киля ледяного образования
параллелепипедов с высотой 0,5 м и горизонтальными размерами 1 х 1 м.
При таком пространственном разрешении конечно-элементная сетка модели содержит 3280580 элементов, среднее значение коэффициента асимметрии равно 0,19, а стандартное отклонение составляет 0,12.
Такая конфигурация расчетной области модели была использована при моделировании процесса размыва связного грунта, который рассматривался как неньютоновская жидкость.
Во всех проведенных численных экспериментах использовались следующие граничные условия:
• условие «прилипания» на всех твердых границах (боковые, а также верхняя и нижняя грани расчетной области);
• условие равенства нулю статического давления воды на границе водовыпуска;
• скорость воды на границе водовпуска задавалась в виде профиля, характерного для подледного сизигийного приливного течения в Байдарацкой губе.
Результаты численных экспериментов Размыв несвязного грунта
Основная цель использования разработанной в среде ВК ANSYS модели обтекания киля ледяного образования (ОКЛО) состоит в определении расположения зон возможного размыва дна вблизи киля вследствие возникающей пространственной неоднородности скорости течения в придонном слое.
Линейные размеры киля тороса позволяют полагать гранулометрические свойства грунта одинаковыми для всей расчетной области, следовательно, локализация зон с различной интенсивностью размыва дна в области ледяного образования определяется вертикальной изменчивостью вектора скорости течения в придонном слое и имеет пространственную корреляцию с зонами локальных максимумов поля касательного напряжения трения (эрозийного потенциала) на поверхности дна.
Для оценки эрозийного потенциала было выполнено моделирование динамики жидкости с помощью модели ОКЛО, использующей расчетную область, изображенную на рис. 1. Принимая во внимание существование стационарных зон возможной эрозии, вызванных воздействием на донный грунт как подкилевых подковообразных вихрей, так и следных вихрей в областях их генерации и отрыва у боковых стенок тороса, численное моделирование скорости течения было выполнено в стационарной постановке, что позволило существенно сократить время расчета и требуемые вычислительные ресурсы.
Максимальное значение скорости течения воды на границе водовпуска равнялось 1,2 м/с и было приурочено к глубине 4 м, при этом среднее по глубине потока значение скорости составляло 0,7 м/с.
На рис. 4 приведено распределение величины касательного напряжения при минимальном расстоянии ото дна до киля 1,5 м, получен-
№ 3 (14) / 2013
Современные подходы и перспективные технологии в проектах освоения нефтегазовых месторождений российского шельфа
137
Рис. 4. Распределение величины касательного напряжения при минимальном расстоянии ото дна до киля 1,5 м. Красной стрелкой обозначено направление приливного течения
ное на основании стационарного расчета скорости течения.
Из данных рис. 4 следует, что на донной поверхности сформировались зоны локальных максимумов касательного напряжения, приуроченные в основном к наиболее глубоким частям киля тороса (рис. 5). Анализ векторов скоро-
сти течения на продольном (вдоль направления приливного течения) сечении расчетной сетки модели, проходящем через эти области, показал, что увеличение касательного напряжения обусловлено повышением скорости воды, вызванным уменьшением живого сечения канала (рис. 6). Однако наиболее протяженная область
Рис. 5. Киль тороса с нанесенными цветными изолиниями поля касательного напряжения - изометрический вид снизу
№ 3 (14) / 2013
138
Научно-технический сборник • ВЕСТИ ГАЗОВОЙ НАУКИ
Рис. 6. Векторы скорости течения на продольном (вдоль направления приливного течения) сечении расчетной сетки модели. На донной поверхности изображены изолинии поля
касательного напряжения
локального максимума поля касательного напряжения длиной около 20 м при ширине около 3 м связана с образованием подковообразного вихря позади области, обозначенной символом Р1 на рис. 4.
На основании полученного в ходе численного эксперимента распределения величины касательного напряжения с помощью встроенных функций блока постпроцессинга ВК ANSYS были определены площади зон эрозии для разных типов несвязного донного грунта. Под килем тороса мелкий гравий может размываться на площади 82 м2 в областях, отмеченных красным или оранжевым цветом на рис. 4 или оконтуренных изолиниями аналогичного цвета на рис. 5. Размыв грунта, состоящего из очень мелкого гравия, возможен на площади 1482 м2 в областях, расположенных внутри изолиний зеленого цвета, а размыв песчаного грунта (песков) - на всей площади дна под килем тороса.
Важной особенностью полученного в ходе выполненного численного эксперимента поля касательного напряжения является наличие у него области локального максимума, расположенной на расстоянии около 50 м вниз по потоку от тороса (см. рис. 4). Из теоретической гидродинамики, численных и натурных экспериментов [4-7] известно, что при достаточно больших числах Рейнольдса в процессе обтекания
жидкостью цилиндрической преграды вблизи ее боковой поверхности генерируется система вихрей с вертикальной осью (следные вихри), возникновение которых связано с отрывом потока по бокам обтекаемого тела вследствие местного сжатия и увеличения скорости в 1,5—2,0 раза по сравнению со скоростью невозмущенного потока. Морские течения, особенно в области килей ледяных образований, имеют турбулентный характер, который определяет значительную величину числа Рейнольдса (порядка 105-106). В этом случае наблюдаются срыв след-ных вихрей сначала с одной стороны цилиндра, затем с другой, и их миграция вниз по потоку с последующим их дроблением, слиянием и диссипацией. Каждый концентрированный следный вихрь со своим центром низкого давления действует как «вакуумный очиститель» [8] и при определенной интенсивности может вызвать эрозию донного грунта.
Выполненный анализ локализации регионов расчетной области модели ОКЛО, в которых происходила эволюция турбулентных вихрей, позволил выявить наличие пространственной корреляции между вышеописанной областью локального максимума поля касательного напряжения и соответствующим регионом турбулентных движений.
Таким образом, в результате численного моделирования показано, что при заданном
№ 3 (14) / 2013
Современные подходы и перспективные технологии в проектах освоения нефтегазовых месторождений российского шельфа
139
профиле приливного течения вниз по потоку от киля тороса под воздействием следных вихрей формируется область донной эрозии, имеющая ширину около 15 м, в которой возможен размыв грунта, состоящего из очень крупного песка.
Размыв связного грунта
При моделировании размыва связного грунта в области киля ледяного образования процесс донной эрозии рассматривался как изменение формы свободной поверхности более плотной неньютоновской жидкости под воздействием касательных напряжений, вызываемых турбулентным движением над ней менее плотной ньютоновской жидкости (водой).
При такой постановке задачи начальное положение донной поверхности определяется высотой слоя неньютоновской жидкости относительно нижней границы расчетной области модели (см. рис. 3).
В проведенных численных экспериментах начальная высота слоя неньютоновской жидкости была выбрана таким образом, чтобы наименьшее расстояние между наиболее глубокой частью тороса и дном равнялось 1 м.
С учетом реологических свойств глинистых илов неньютоновская жидкость была выбрана удовлетворяющей соотношению Кессона [9].
На рис. 7 изображена изоповерхность раздела двух сред - воды и неньютоновской жидкости, которая отражает форму поверхности дна под килем тороса на 6000-й с расчета. Цвет изоповерхности соответствует величине придонной скорости приливного течения под килем ледяного образования.
Из данных рис. 7 следует, что зоны эрозии связного грунта образовались под наиболее глубокими частями киля тороса, причем их локализация совпадает с областями возможного размыва несвязного грунта, которые отмечены на рис. 4 символами Р1 и Р2.
Однако необходимо отметить отсутствие зоны усиления придонной скорости течения и, соответственно, зоны эрозии связного донного грунта в области, расположенной вниз по потоку от области Р1, где в случае численного эксперимента по моделированию размыва несвязного грунта наблюдались локальные максимумы касательного напряжения (см. рис. 4). Для выяснения причины этого явления были проанализированы векторы скорости течения на продольных (вдоль направления приливного течения) сечениях расчетной сетки модели, пролегавших через область Р1. В результате было установлено, что в отличие от предыдущего численного эксперимента при моделировании совместной динамики течений воды
Рис. 7. Форма донной поверхности на 6000-й с расчета.
Цвет поверхности соответствует величине придонной скорости приливного течения
№ 3 (14) / 2013
140
Научно-технический сборник • ВЕСТИ ГАЗОВОЙ НАУКИ
и неньютоновской жидкости не произошло формирования подковообразного вихря (см. рис. 6), который определял высокие значения касательного напряжения в рассматриваемой области донной поверхности. Возможно, это связано с меньшим расстоянием от киля тороса до донной поверхности.
Важным результатом проведенного численного эксперимента является установление того факта, что в отличие от случая обтекания твердого тела в однофазной среде при моделировании динамики двухфазной среды наибольшие величины завихренности трехмерного поля течения ньютоновской жидкости возникают
Список литературы
1. ANSYS-CFD. Environment: User Guide.
2. Menter F.R. Two-equation eddy-viscosity turbulence models for engineering applications / F.R. Menter // AIAA-Journal. - 1994. -
№ 32 (8). - Р. 1598-1605.
3. Menter F.R. Review of the shear-stress transport turbulence model experience from an industrial perspective / F.R. Menter // International Journal of Computational Fluid Dynamics. - 2009. -
V. 23. - № 4. - Р. 305-316.
4. Aghaee Y. Three Dimensional Numerical Modeling of Flow around Bridge PiersUsing LES and RANS / Y. Aghaee, H. Hakimzadeh // River Flow, Dittrich, Koll, Aberle & Geisenhainer (eds). - 2010. - Р 211-218.
5. Dargahi B. The turbulent flow field around a circular cylinder / B. Dargahi // Experiments in Fluids. - 1989. - № 8. - Р 1-12.
непосредственно над поверхностью раздела фаз, что отражает существенную роль мелкомасштабных турбулентных вихрей в формировании рельефа, а в рассматриваемой интерпретации - в формировании зон эрозии связного грунта.
Это обстоятельство определяет необходимость применения существенно более мелкомасштабной расчетной сетки для адекватного воспроизведения турбулентности и, соответственно, получения реалистичных оценок скорости размыва связного грунта и определения возможной глубины воронки размыва в зависимости от гидродинамического режима акватории.
6. Dargahi B. Controlling mechanism of local scouring / B. Dargahi // Journal of Hydraulic Engineering. ASCE. - 1990. - № 116 (10). -Р 1197-1214.
7. Roulund A. Numerical and experimental investigation of flow and scour around a circular pier / A. Roulund, B. Sumer, J. Fredsoe et. al //
J. Fluid Mech. - 2005. - V 534. - Р 351-401.
8. Breusers H.N.C. Local scour around cylindrical piers / H.N.C. Breusers, G. Nicollet, H.W. Shen // J. Hydraulic Res. - 1977. - V 15. - Р 211-252.
9. Casson W. A flow equation for pigment-oil suspensions of printing of the printing ink type. Rheology of Dispersed Systems / W. Casson. -London: Pergamon,1959.
№ 3 (14) / 2013
