Математические заметки СВФУ Июль—сентябрь, 2017. Том 24, № 3
УДК 519.63
ТЕЧЕНИЕ И ПЕРЕНОС В ПЕРФОРИРОВАННЫХ И ТРЕЩИНОВАТЫХ ОБЛАСТЯХ С НЕОДНОРОДНЫМИ ГРАНИЧНЫМИ УСЛОВИЯМИ РОБИНА
У. С. Гаврильева, В. Н. Алексеев, М. В. Васильева
Аннотация. Рассматриваются задачи переноса и течения в перфорированных и трещиноватых областях. Система уравнений описывается уравнением Стокса для моделирования течения флюида в области и уравнением переноса концентрации некоторого вещества. Перенос концентрации дополняется неоднородными граничными условиями третьего рода, которые моделируют происходящую реакцию на гранях моделируемого объекта. Для численного решения задачи строится конечно-элементная аппроксимация уравнения. Для получения устойчивого решения задачи переноса используется метод SUPG (streamline upwind Petrov—Galerkin) для стабилизации классического метода Галеркина. Вычислительная реализация основана на вычислительной библиотеке FEniCS. Представлены результаты решения модельной задачи в перфорированных и трещиноватых областях. Численно исследованы различные режимы неоднородных граничных условий.
DOI 10.25587/SVFU.2018.3.10890 Ключевые слова: уравнение переноса, задача Стокса, перфорированная область, трещиноватая область, численное моделирование, граничное условие Робина, численная стабилизация, SUPG, метод конечных элементов.
Введение
В работе рассматривается численное решение задачи конвективного и диффузионного переноса, возникающего при моделировании течений жидкости в пористых средах, а также в трещинах при рассмотрении их на микроуровне [1-5]. Основной особенностью рассматриваемого процесса является учет реакции на стенках среды (перфорациях или стенках трещин), которая описывается неоднородным граничным условием третьего рода (граничное условие Робина) [6, 7]. Течение жидкости в области описывается уравнением Стокса [8]. Для аппроксимации уравнения Стокса используется смешанный метод конечных элементов, позволяющий строить аппроксимации, удовлетворяющие LBB-условию (inf-sup) [9,10]. Система уравнений дополняется уравнением переноса концентрации некоторого вещества. При решении таких задач возникают трудности для расчета уравнения переноса в случае, когда конвективный перенос
Работа выполнена за счет гранта РНФ 17—71—20055.
© 2017 Гаврильева У. С., Алексеев В. Н., Васильева М. В.
существенно преобладает над диффузионным [11—14]. Стандартная аппроксимация с использованием классического метода Галеркина может привести к возникновению осцилляций в решении задачи. Для борьбы с осцилляция-ми используются методики стабилизации [11,12]. Среди классических методов стабилизации можно выделить такие методы, как противопотоковая аппроксимация конвективного слагаемого посредством введения искусственной диффузии, метод SUPG (streamline upwind Petrov-Galerkin) или разрывный метод Галеркина с противопотоковыми аппроксимациями конвективного слагаемого [15-17]. В данной работе будем использовать метод SUPG. Для простоты будем предполагать что свойства жидкости не зависят от концентрации, следовательно, полученную систему уравнений можем решать последовательно: (1) расчет течения жидкости посредством решения уравнения Стокса и (2) расчет нестационарной задачи переноса (уравнения конвекции-диффузии) при заданных скоростях течения [18,19].
Работа состоит из четырех разделов. В разд. 1 представлена математическая модель рассматриваемого процесса. В разд. 2 представлена конечно-элементная аппроксимация. В разд. 3, 4 описаны результаты решения модельной задачи в перфорированных и трещиноватых областях, а также численно исследованы различные режимы неоднородных граничных условий.
Рассматривается реактивный перенос и течение в области О (рис. 1). Однофазный поток несжимаемой жидкости описывается системой уравнений Стокса, которая состоит из уравнения движения и уравнения неразрывности
где и(ж) — скорость течения жидкости, р(х) — динамическое давление и ^ — динамическая вязкость.
Уравнения потока (1.1) дополняются граничными условиями непротекания на границе твердого тела и жидкости Г,; = (Г2 и Г4 и Г5)
Поток обусловлен граничными условиями, наложенными на дО/Г, — внешние границы О:
Жидкость содержит растворенные частицы с некоторой молярной концентрацией с(х, £), которая подвергается конвекции обусловленного перемещения флюида, молекулярной диффузии и линейной гетерогенной реакции на границе твердой и жидкой сред Г,;. Уравнение для концентрации с(х, £) описывается нестационарным уравнением конвекции-диффузии
1. Математическая модель
-^Au + Vp = 0, x е о, V • u = 0, x е о,
(1.1)
u = 0, x е rs;.
u = g, x е дО/Г8;.
— + V • (uc) - W2c = 0
dt
(1.2)
Ч,' Г5
/ V
гчЛ ,
л
/
Рис. 1. Область с перфорациями (сверху) и трещина (снизу)
с учетом граничных условий на границе твердой и жидкой сред Г^
-^Ус • п = Ж(с - со), х е Г^, (1.3)
и граничными условиями на дП/Г81
с(х,£) = Н(х), х е 0П/Г31, (1.4)
где @ — коэффициент молекулярной диффузии, ^(х) — коэффициент поверхностной реакции и п(х) — внешняя нормаль к .
2. Аппроксимация методом конечных элементов
Для численного решения задачи проведем аппроксимацию уравнений переноса и течения с использованием метода конечных элементов. Для неизвестных скорости и давления введем соболевские пространства V(П), ((П), состоящие из векторных функций и и скалярной величины р. Для системы уравнений (1.1) запишем вариационную постановку в следующем виде: найти (и,р) е V(П) х ((П) такие, что
J ^У и • Уvdx - J р vdx = 0, V е V5 (П), J uqdx = 0, д е ( (П), (2.1) п п п
где V(П) и ((П) — пространства тестовых функций. В качестве базисных функций будем использовать элементы Тэйлора — Худа (ТН), удовлетворяющие условию устойчивости.
Проведем дискретизацию уравнения конвективного и диффузионного переноса (1.2) с помощью классического метода Галеркина, в котором тестовые и пробные функции совпадают. Пусть с, г е ((П), вариационную постановку с учетом граничного условия Робина запишем следующим образом:
+ /и.ус"«г„х + /+ /Х(с-+' -со)гds = 0,
п г8,
(2.2)
-zdx + J и п
сп+1 - с
п
для ^ € Здесь для аппроксимации по времени использована устойчивая неявная разностная схема.
Систему уравнений запишем в матричной форме
сп+1 _ _п
М-+ Асп+1 = (2.3)
т
где
М = \rnij] = J фiфj ¿ж, ^ = [/] = J Жс0ф^ ¿в,
О Г3г
А = ] = J (и, Уф^ф,- ¿ж + У (^ Vфi, Уф,-) ¿ж + J Жфiфj ¿в
О О Гзг
и с = ^ ciфi. В качестве базисных функций фi будем использовать стандартные
i= 1
линейные базисные функции.
Представленная аппроксимация может приводить к появлению осцилляций в решении, поэтому необходимо использовать специальные схемы аппроксимации со стабилизацией. В отличие от метода с искусственной диффузией, метод 8иРС обеспечивает более высокий порядок аппроксимации. Основной идеей 8иРС является модификация тестовых функций с учетом направления течения. В методе 8ИРС используется вариационная постановка: найти сп+1 € ^ такое, что
I С + ~ С г йх+ I и-Усп+1гг1х + 1'(9Чсп+1Уг)в,х ^ Ж(сп+1 — с0)£ гЬ = О, О О О Г3г
(2.4)
где г £ г = (г + ' Уг) и ^ — диаметр ячейки сетки. Таким образом, поскольку пробное пространство (3 и тестовое пространство (3 не совпадают, этот метод является методом Петрова — Галеркина (РС, Ре1гоу-Са1егкт). Полученная аппроксимация уже не может быть интерпретирована как схема с искусственной диффузией, поскольку само уравнение остается прежним, а меняется только тестовая функция 2.
3. Численные результаты для перфорированной области
Рассмотрим моделирование течения и переноса в перфорированной области О = [0,1] х [0, 5]. Расчет проведем для Ттаж = 0.7 с шагом по времени т = 0.02. Зададим следующие параметры математической модели: вязкость ^ = 1.0, коэффициент молекулярной диффузии 'З = 0.005. Моделирование проведем для различных коэффициентов реакции Ж = 10-3, Ж = 10-4, Ж = 10-5. Для численного решения строим треугольную расчетную сетку, которая содержит 18766 ячеек и 9993 узлов (рис. 2).
Рассмотрим моделирование двух модельных задач в перфорированной расчетной области (см. рис. 2). На рис. 3 представлены распределения давления и скоростей, которые будем использовать для задания конвективного переноса.
Рис. 2. Треугольная расчетная сетка для области О = [0,1] х [0, 5] для перфорированной области. Задача 1
Рис. 3. Распределение давления (р) и скоростей (их , иу) для задачи Стокса в перфорированной области
Задача 1а. Втекает жидкость с некоторой концентрацией загрязнителя (рис 4)).
Задача 1Ь. Втекает чистая жидкость (нулевая концентрация загрязнителя) , на перфорациях происходит реактивный массообмен за счет неоднородных граничных условий Робина (рис. 5)).
На рис. 4 представлено распределение концентрации на конечный момент времени в перфорированной области для задачи 1а. Рисунок иллюстрирует перенос концентрации при задании втока некоторой концентрации, где на перфорациях с0 = 0, Ж = 10-3 и Ж = 10-4. Для данной задачи в конечный момент моделирования наивысшая концентрация распределяется до четверти рассматриваемой геометрии, при Ж = 10-5 доходит до половины геометрии.
На рис. 5 рассматривается аналогичный случай для задачи 1Ь, где втекает
Рис. 4. Распределение концентрации с различными коэффициентами реакции = 10-3, = 10-4, = 10-5 на конечный момент времени £ = 0.7 в перфорированной области. Задача 1а
Рис. 5. Распределение концентрации с различными коэффициентами реакции = 10-3, = 10-4, = 10-5 на конечный момент времени £ = 0.7 в перфорированной области. Задача 1Ь
чистая жидкость с нулевой концентрацией. Изолинии концентрации для значений с = 0.3 (задача 1а) и с = 0.5 (задача 1Ь) при разных коэффициентах реакции на конечный момент времени представлены на рис. 6. На рис. 7 изображен вертикальный срез по середине области для х = 5. Из полученных результатов можно сделать вывод, что при наименьшем значении коэффициента реакции на порах концентрация огибает их и распространяется дальше.
Рис. 6. Изолинии концентрации для значении с = 0.3 (рисунок сверху для задачи 1а) и с = 0.5 (рисунок снизу для задачи 1Ь) при различных коэффициентах реакции Ж = 10-3(красная линия), Ж = 10-4(зеленая линия), Ж = 10-5(синяя линия) на конечный момент времени £ = 0.7 в перфорированной области
Рис. 7. Вертикальный срез для концентрации при разных коэффициентах реакции Ж = 10-3, Ж = 10-4, Ж = 10-5 на конечный момент времени £ = 0.7 в перфорированной области. Сверху задача 1а, снизу задача 1Ь
4. Численные результаты для области с трещиной
Рассмотрим случай, когда геометрия является областью трещины. Для моделирования будем использовать значения параметров, аналогичные представленным выше для перфорированной области. Моделирование проведем для различных коэффициентов реакции Ж = 10-4, Ж = 10-5 и Ж = 10-6. Для численного решения строим треугольную расчетную сетку, которая содержит 57849 ячеек и 29692 узлов (рис. 8).
Рассмотрим моделирование двух модельных задач в расчетной области (см. рис. 8).
Задача 2а. Втекает жидкость с некоторой концентрацией загрязнителя (см. рис. 4).
Задача 2Ъ. Втекает чистая жидкость (нулевая концентрация загрязнителя) , на верхней и нижней границе трещины происходит реактивный массообмен за счет неоднородных граничных условий Робина (см. рис. 5).
Рис. 8. Треугольная расчетная сетка П = [0,1] х [0,10] для области трещины. Задача 2
ш АГ| А» I
Рис. 9. Распределение давления (р) и скоростей (их, иу) для задачи Стокса для области трещины
В задаче 2 рассматривается область трещины. На рис. 9 представлены распределения давления и скоростей, которые будем использовать для задания конвективного переноса. На рис. 10 показано распределение концентрации с коэффициентами реакции Ж = 10-3, Ж = 10-4 и Ж = 10-5 при с0 = 0 на верхней и нижней границах (задача 2а). Отметим, что коэффициент реакции не сильно влияет на концентрацию. На рис. 11 рассматривается аналогичный случай (задача 2Ь) для коэффициентов реакции Ж = 10-4, Ж = 10-5, Ж = 10-6 и с0 = 1 на верхней и нижней границах. В данном случае коэффициенты реакции уже сильно влияют на распределение концентрации.
На рис. 12 представлены изолинии концентрации при значении с = 0.5 (задача 2а) и с = 0.3 (задача 2Ь) на конечный момент времени. На рис. 13 показан вертикальный срез по середине области х = 5. Из результатов для задачи 2 видно, что при наименьшем значении коэффициента реакции верхние и нижние границы не препятствуют распространению концентрации.
Вычислительная реализация построена с использованием библиотеки
Рис. 10. Распределение концентрации для задачи Стокса с разными коэффициентами реакции Ж = 10-3, Ж = 10-4, Ж = 10-5 на конечный момент времени £ = 0.7 секунд. Задача 2а
Рис. 11. Распределение концентрации для задачи Стокса с разными коэффициентами реакции Ж = 10-4, Ж = 10-5, Ж = 10-6 на конечный момент времени £ = 0.7 секунд. Задача 2Ь
РЕшСЯ. Для построения геометрической области используется программа СМЯЫ. Визуализация происходит с использованием программы Paraview.
5. Заключение
В данной работе рассмотрены задачи конвективного и диффузионного переноса с учетом неоднородных граничных условий Робина. Для численного решения построены схемы аппроксимации со стабилизацией для конечно-
Рис. 12. Изолинии концентрации на значении с = 0.5 для задачи 2а (рисунок сверху) и с = 0.3 для радачи 2Ь (рисунок снизу) при разных коэффициентах реакции на конечный момент времени £ = 0.7. Задача 2а: Ж = 10-3 (синяя линия), Ж = 10-4 (зеленая линия), Ж = 10-5 (красная линия). Задача 2Ь: Ж = 10-4 (красная линия), Ж = 10-5 (зеленая линия), Ж = 10-6 (синяя линия)
Рис. 13. Вертикальный срез для концентрации при разных коэффициентах реакции на конечный момент времени t = 0.7. Задача 2a: Ж = 10-3, Ж = 10-4 и Ж = 10-5. Задача 2b: Ж = 10-4, Ж = 10-5 и Ж = 10-6
элементной аппроксимации задачи переноса. Проведено математическое моделирование и представлены результаты численного решения модельных задач, возникающего при моделировании течений жидкости в пористых средах, а также в области трещин, при рассмотрении их на микроуровне для различных коэффициентах реакции.
ЛИТЕРАТУРА
1. Chung E. T. et al. Multiscale model reduction for transport and flow problems in perforated domains //J. Comput. Appl. Math. 2018. V. 330. P. 519-535.
2. Chung E. T., Vasilyeva M., Wang Y. A conservative local multiscale model reduction technique for Stokes flows in heterogeneous perforated domains //J. Comput. Appl. Math. 2017. V. 321. P. 389-405.
3. Васильева М. В., Прокопьев Г. А. Численное решение задачи двухфазной фильтрации с
неоднородными коэффициентами методом конечных элементов // Мат. заметки СВФУ. 2017. Т. 24. № 2. С. 46-62.
4. Chung E. T., Iliev O., Vasilyeva M. V. Generalized multiscale finite element method for non-Newtonian fluid flow in perforated domain // AIP Conference Proceedings. AIP Publishing, 2016. V. 1773. N 1. P. 100001.
5. Chung E. T. et al. Generalized multiscale finite element methods for problems in perforated heterogeneous domains // Appl. Anal. 2016. V. 95. N 10. P. 2254-2279.
6. Battiato I. et al. Hybrid models of reactive transport in porous and fractured media // Adv. Water Resources. 2011. V. 34. N 9. P. 140-1150.
7. Battiato I., Tartakovsky D. M. Applicability regimes for macroscopic models of reactive transport in porous media //J. contaminant hydrology. 2011. V. 120. P. 18-26.
8. Tomin P., Lunati I. Hybrid Multiscale Finite Volume method for two-phase flow in porous media // J. Comput. Physics. 2013. V. 250. P.'293-307.
9. Brezzi F., Fortin M. Mixed and hybrid finite element methods. New York: Springer Science & Business Media, 2012. (Springer Ser. Comput. Math.; V. 15).
10 Raviart P. A., Thomas J. M. A mixed finite element method for 2-nd order elliptic problems // Mathematical aspects of finite element methods. Berlin; Heidelberg: Springer-Verl, 1977. P. 292-315.
11. Donea J., Huerta A. Finite element methods for flow problems. Chichester: John Wiley & Sons, 2003.
12. Brooks A. N. A Petrov-Galerkin finite element formulation for convection dominated flows: diss. California Institute of Technology, 1981.
13. Vabishchevich P. N., Vasil'eva M. V. Explicit-implicit schemes for convection-diffusion-reaction problems // Numer. Anal. Appl. 2012. V. 5, N 4. P. 297-306.
14. Afanas'eva N. M., Vabishchevich P. N., Vasil'eva M. V. Unconditionally stable schemes for convection-diffusion problems // Russian Mathematics. 2013. V. 57, N 3. P. 1-11.
15. Brooks A. N., Hughes T. J. R. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations // Comp. Methods Appl. Mechanics and Engineering. 1982. V. 32, N 1-3. P. 199-259.
16. Rivire B., Wheeler M. F. Discontinuous Galerkin methods for flow and transport problems in porous media // Intern. J. Numer. Methods in Biomedical Engineering. 2002. V. 18, N 1. P. 63-68.
17. Cockburn B., Shu C. W. The local discontinuous Galerkin method for time-dependent convection-diffusion systems // SIAM J. Numer. Anal. 1998. V. 35, N 6. P. 2440-2463.
18. Алексеев В. Н., Васильева М. В., Степанов С. П. Итерационные методы решения для задачи течения и переноса в перфорированных областях // Вестн. СВФУ. 2016. № 5. С. 67-79.
19. Васильева М. В., Васильев В. И., Тимофеева Т. С. Численное решение методом конечных элементов задач диффузионного и конвективного переноса в сильно гетерогенных пористых средах // Уч. зап. Казан. ун-та. Сер. Физико-математические науки. 2016. Т. 158, № 2, С. 243-261.
Статья поступила 6 июля 2017 г.
Гаврильева Уйгулаана Семеновна, Алексеев Валентин Николаевич, Васильева Мария Васильевна
Северо-Восточный федеральный университет имени М. К. Аммосова, Институт математики и информатики, ул. Кулаковского, 42, Якутск 677000
[email protected], [email protected], [email protected]
Математические заметки СВФУ Июль—сентябрь, 2017. Том 24, № 3
UDC 622.691.4
FLOW AND TRANSPORT IN PERFORATED
AND FRACTURED DOMAINS
WITH ROBIN BOUNDARY CONDITIONS
U. S. Gavrilieva, V. N. Alekseev, and M. V. Vasilyeva
Abstract: We consider transport and flow problems in perforated and fractured domains. The system of equations is described by the Stokes equation for modeling fluid flow and the equation for the concentration transfer of a certain substance. Concentration is supplemented by inhomogeneous boundary conditions of the third type which simulate the occurring reaction on the faces of the modeled object. For the numerical solution of the problem, a finite-element approximation of the equation is constructed. To obtain a sustainable solution to the transport problem, the SUPG (streamline upwind Petrov—Galerkin) method is used to stabilize the classical Galerkin method. The computational implementation is based on the Fenics computational library. The results of solving the model problem in perforated and fractured domains are presented. Numerical studies of various regimes of heterogeneous boundary conditions were carried out.
DOI 10.25587/SVFU.2018.3.10890 Keywords: transport equation, Stokes problem, perforated domain, fractured domain, numerical modeling, Robin boundary condition, numerical stabilization, SUPG, finite element method.
REFERENCES
1. Chung E. T., Leung W. T., Vasilyeva M., and Wang Y., "Multiscale model reduction for transport and flow problems in perforated domains," J. Comput. Appl. Math., 330, 519—535 (2018).
2. Chung E. T., Vasilyeva M., and Wang Y., A conservative local multiscale model reduction technique for Stokes flows in heterogeneous perforated domains," J. Comput. Appl. Math., 321, 389-405 (2017).
3. Vasilyeva M. V. and Prokopyev G. A., "Numerical solution of the two-phase filtration problem with nonhomogenuous coefficients by the finite element method," Mat. Zamet. SVFU, 24, No. 2, 46-62 (2017).
4. Chung E. T., Iliev O., and Vasilyeva M. V., "Generalized multiscale finite element method for non-Newtonian fluid flow in perforated domain," in: AIP Conf. Proc., 1773, No. 1, pp. 100001. AIP Publ. (2016).
5. Chung E. T., Efendiev Y., Li G., and Vasilyeva M. V., "Generalized multiscale finite element methods for problems in perforated heterogeneous domains," Appl. Anal., 95, No. 10, 22542279 (2016).
6. Battiato I., Tartakovsky D. M., and Scheibe T. D., "Hybrid models of reactive transport in porous and fractured media," Adv. Water Resources, 34, No. 9, 140-1150 (2011).
The authors were supported by the Russian Science Foundation (Grant 17-71-20055).
© 2017 U. S. Gavrilieva, V. N. Alekseev, M. V. Vasilyeva
7. Battiato I. and Tartakovsky D. M., "Applicability regimes for macroscopic models of reactive transport in porous media," J. Contaminant Hydrology, 120, 18—26 (2011).
8. Tomin P. and Lunati I., "Hybrid Multiscale Finite Volume method for two-phase flow in porous media," J. Comput. Phys., 250, 293-307 (2013).
9. Brezzi F. and Fortin M., Mixed and Hybrid Finite Element Methods, Springer, New York (1991) (Springer Ser. Comput. Math.; V. 15).
10. Raviart P. A. and Thomas J. M., "A mixed finite element method for 2nd order elliptic problems," in: Mathematical Aspects of Finite Element Methods, pp. 292-315, Springer, Berlin; Heidelberg (1977).
11. Donea J. and Huerta A., Finite Element Methods for Flow Problems, John Wiley & Sons, Chichester (2003).
12. Brooks A. N., A Petrov-Galerkin finite element formulation for convection dominated flows: Thes. ... doct. phylosophy (civil engineering), California Inst. Technology (1981).
13. Vabishchevich P. N. and Vasil'eva M. V., "Explicit-implicit schemes for convection-diffusion-reaction problems," Numer. Anal. Appl., 5, No. 4, 297-306 (2012).
14. Afanas'eva N. M., Vabishchevich P. N., and Vasil'eva M. V., Unconditionally stable schemes for convection-diffusion problems," Russ. Math., 57, No. 3, 1-11 (2013).
15. Brooks A. N. and Hughes T. J. R., "Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations," Comput. Methods Appl. Mech. Eng., 32, No. 1-3, 199-259 (1982).
16. Riviere B. and Wheeler M. F., "Discontinuous Galerkin methods for flow and transport problems in porous media," Int. J. Numer. Methods Biomed. Eng., 18, No. 1, 63-68 (2002).
17. Cockburn B. and Shu C. W., "The local discontinuous Galerkin method for time-dependent convection-diffusion systems," SIAM J. Numer. Anal., 35, No. 6, 2440-2463 (1998).
18. Alekseev V. N., Vasilyeva M. V., and Stepanov S. P., "Iterative methods of solving the flow and transfer problem in perforated domains [in Russian]," Vestn. SVFU, No. 5, 67-79 (2016).
19. Vasilyeva M. V., Vasilyev V. I., and Timofeeva T. S., "Numerical solution by the finite element method of the diffusion and convection problems in strongly heterogenuous media [in Russian]," Uch. Zap. Kazansk. Univ., Ser. Fiz.-Mat. Nauki, 158, No. 2, 243-261 (2016).
Submitted July 6, 2017
Uygulaana S. Gavrilieva, Valentin N. Alekseev, Maria V. Vasilyeva M. K. Ammosov North-Eastern Federal University, Institute of Mathematics and Informatics, 42 Kulakovsky Street, Yakutsk 677891, Russia