Modelització numèrica de problemes acoblats de flux en medi porós i flux de Stokes

PROJECTE O TESINA D’ESPECIALITAT Títol Modelització numèrica de problemes acoblats de flux en medi porós i flux de Stokes 727-TES-CA-6431 Autor/a Pedro Postigo Laguna Tutor/a Sonia Fern á ndez M é ndez Departament Departament de M atemàtica Aplicada Intensificació Enginyeria Computacional Data Juny de 2014 Numerical modeling of porous medium flow and Stokes flow coupled problems Pedro Postigo Laguna, Tutor: Sonia Fern´andez-M´endez Numerical modeling of porous medium flow and Stokes flow coupled problems Pedro Postigo Laguna Tutor: Sonia Fernández-Méndez 2 Abstract Porous medium flow and Stokes flow coupled problems are recurrent in civil engineering applications. Such is the case of the study of the groundwater flow, or the erosion of a material and the consequences that it can have on a structure. To solve these problems, the standard procedure is to model them using the Finite Element Method, which is implemented in the majority of computer software used in civil engineering. In this kind of applications, it is also usual the need to consider the transient effects, which can make the interface move and change with time, increasing the complexity of the problem. When it is the case, the classical Finite Element Method approach can present several disadvantages: an increase in the number of elements needed to model correctly the interface, and the need to redo the mesh when the interface moves. These disadvantages are translated into higher computation costs and larger computation time required to compute the solution. In this work, a levelset function based strategy is proposed, that avoids this kind of problems, computing the solution faster and using less computer resources. First, the equations governing the Darcy and Stokes coupled flow problem are deduced, applying the compatibility conditions between the two domains. The compatibility conditions are what make the problem coupled. The usual compatibility conditions consist in making the velocity normal to the interface and the tangent velocity equal for the porous domain and the free flow domain. These equations are independent from the numerical strategy to solve them. Second, the element integration and the system matrices assembling process are modified, to adapt it to the possibility that the element is cut by the interface. Some implementation details are also commented, for example, the way to obtain the integration points in a cut element, or how to compute the integrals over the interface line, which come from the compatibility conditions. Finally, the usual error convergence tests are applied to the model, to ensure that it really converges to the real solution when the element size is reduced. To compute the error of the model, an analytical solution for a simple problem is deduced, and then it is compared with the computed solution. Other tests are performed to the model to ensure that the model is robust when the complexity of the problem increases. These tests also help to debug programming errors. In this work it is also presented a practical application of the model, computing numerically the solution to a problem that has not an analytical solution. To conclude, the possibilities of this kind of models are discussed, highlighting the advantages that it has in front of the classical FEM used to solve this kind of problems in the context of civil engineering. In all the work only the two dimensional problem is studied, so it is explained here the steps that would be necessary to implement a 3D levelset based model. Keywords: Coupled flow problem, Stokes equation, Darcy equation, Levelset function, Extended Finite Element Method. Modelización numérica de problemas acoplados de flujo en medio poroso y flujo de Stokes Pedro Postigo Laguna Tutor: Sonia Fernández-Méndez 3 Resumen Los problemas acoplados de flujo en medio poroso y flujo de Stokes son recurrentes en aplicaciones de ingeniería civil, por ejemplo, para estudiar el flujo de agua subterránea, o determinar la erosión de un material y las consecuencias que puede tener ésta sobre una estructura. Para resolver estos problemas, el procedimiento estándar es modelizarlos usando el Método de los Elementos Finitos, que está implementado en la mayoría de los paquetes informáticos usados en ingeniería civil. En este tipo de aplicaciones, suele ser habitual tener en cuenta los efectos transitorios, lo que provoca que la interfaz pueda moverse con el tiempo, incrementando la complejidad del problema. Cuando es este el caso, el Método de los Elementos Finitos clásico puede presentar diversos inconvenientes: un aumento del número de elementos necesarios para modelizar correctamente la interfaz, y la necesidad de rehacer la malla en cada iteración del movimiento de la misma. Estas desventajas se traducen en un aumento del coste computacional y un aumento del tiempo necesario para calcular la solución. En este trabajo, se propone una estrategia basada en funciones de nivel, que evita este tipo de desventajas, y calcula la solución más rápidamente y consumiendo menos recursos computacionales. Primero, se deducen las ecuaciones que gobiernan el problema acoplado de Darcy y Stokes, aplicando condiciones de compatibilidad entre los dos dominios. Estas condiciones de compatibilidad son las que provocan el acoplamiento del problema. Las condiciones de compatibilidad normales consisten en imponer que la velocidad normal y tangencial a la interfaz sea igual para el dominio poroso y el de flujo libre. Estas ecuaciones son independientes de la estrategia usada para resolverla. En segundo lugar, se modifica el sistema de integración en el elemento y el proceso de ensamblaje de las matrices, para adaptarlos a la posibilidad de que el elemento esté cortado por la interfaz. Algunos detalles sobre la implementación se comentan en este capítulo, por ejemplo, cómo se obtienen los puntos de integración para un elemento cortado, o cómo se calculan las integrales sobre el segmento de interfaz, que provienen de las condiciones de compatibilidad. Por último, se realizan sobre el modelo los test habituales de convergencia del error, para asegurar que este converge hacia la solución exacta cuando se reduce el tamaño del elemento. Para calcular el error del modelo, se deduce una solución analítica para un problema simple, y luego se compara esta con la solución del modelo. Se han realizado también otros test, para asegurar que el modelo es robusto a medida que aumenta la complejidad de la solución. Estos test también ayudan a depurar errores de programación. También se presenta en este trabajo una aplicación práctica del modelo, calculando una solución numérica para un problema sin solución analítica. Para concluir, se comentan las posibilidades que tiene este modelo, destacando las ventajas que tiene frente a la manera clásica de modelar estos problemas usando el Método de los Elementos Finitos en el contexto de la ingeniería civil. En todo el trabajo se considera el problema en 2 dimensiones, pero también se explica aquí las consideraciones necesarias para implementar el modelo en 3 dimensiones. Palabras clave: Problemas de flujo acoplado, Ecuación de Stokes, Ecuación de Darcy, Función de nivel, Método de los Elementos Finitos Extendido. Contents 1 Introduction 5 1.1 Aim of the work . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.1 Example 1: Groundwater and surface water flow . . . . 5 1.1.2 Example 2: Erosion . . . . . . . . . . . . . . . . . . . . 7 1.2 Standard FEM and Extended Finite Element Method . . . . . 7 1.2.1 Levelset functions . . . . . . . . . . . . . . . . . . . . . 8 1.3 Outline, tasks and contributions . . . . . . . . . . . . . . . . . 11 2 Darcy-Stokes coupled problem 13 2.1 Weak form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3 Discretization 16 3.1 Element choice . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2 Linear system of equations . . . . . . . . . . . . . . . . . . . . 18 3.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3.1 Element integration . . . . . . . . . . . . . . . . . . . . 20 3.3.2 Interface integration . . . . . . . . . . . . . . . . . . . 24 4 3.4 Assembling, boundary conditions and degrees of freedom reduction 25 4 Model validation and convergence test 26 4.1 Error convergence . . . . . . . . . . . . . . . . . . . . . . . . . 27 5 Numerical application 29 6 Possibilities of the model 31 6.1 3D Extension . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 6.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 7 Conclusions 33 A Code 35 B Tests to the program 69 B.1 Analytical solutions . . . . . . . . . . . . . . . . . . . . . . . . 69 B.2 ”Decoupled” system computation . . . . . . . . . . . . . . . . 72 B.3 Residuals computation . . . . . . . . . . . . . . . . . . . . . . 73 5 Notation Vectors, functions and vector functions: uwill be used for scalar functions. uwill be used for vector functions. uwill be used for vector constants. Matrices: Kwill be used to indicate a matrix. [K]ij will be used to express the coefficient of a matrix. 6 1. Introduction 1.1. Aim of the work The aim of this work is to set the methodology to implement a levelset based numerical scheme to solve fluid flow and porous medium flow coupled problems. Free flow and porous medium flow coupled problems are recurrent in engineering, even more in civil engineering. In the following subsections two types of such problems will be presented. At the end of this work, these problems will be recovered, explaining the advantages that a levelset based numerical scheme has in solving these problems. 1.1.1. Example 1: Groundwater and surface water flow Groundwater and surface water problems are usual in civil engineering, for example, to calculate the recharge of an aquifer because of the rain, or the effect of the construction of a well on the aquifer. The same set of PDEs can model the transport of contaminant particles through groundwater or surface water. The levelset method that will be proposed in this work can be effective to solve these problems, specially when the interface between soil and water changes, or if such interface has a complex shape. Figure 1shows an scheme of a groundwater flow problem. 7 Figure 1: Example of a groundwater flow problem 8 CreMatStokes: creates the matrices that come from the Stokes equation taking into account the levelset function. CreMatDarcy: creates the matrices that come from the Darcy’s equation taking into account the levelset function. CreMatInterface: creates the matrices that come from the coupling equations taking into account the levelset function. BoundaryConditions: imposes the boundary conditions, and eliminates the unnecessary degrees of freedom ErrorL2Stokes and ErrorL2Darcy : computes the errors in the two domains. Only the major functions are described in this list. Other functions used to create the mesh or plotting the results have also been adapted to introduce the levelset description of the interface. The code of all the functions can be found in A. 2. Darcy-Stokes coupled problem In this section, the numerical model that leads to the solution of the coupled problem is derived. Let us consider a domain Ω like the one that is shown in figure 6. Such domain is divided into two sub-domains, Ωfand Ωp, separated by an interface Γ. This interface has at each point a normal vector n(always pointing outside of the domain Ωf) and a tangent vector t. In the domain Ωf, which is the 15 one that corresponds to the fluid, the Stokes equations are applied: −ν∆u+∇p=bfin Ωf, (2a) −∇·u= 0 in Ωf, (2b) (−ν∇u+pI)·n=ϕn+ (αu·t)ton Γ, (2c) u=uDon ∂Ωf\Γ (2d) where uis the free flow velocity function, pis the pressure function, ϕis the potential function involved in the Darcy equation (see equation 3), and nand tare the normal and tangent vectors to the interface, as explained before. Iis the identity matrix, and αis a material parameter. In the domain Ωp, which is the one that corresponds to the porous medium, the Darcy’s equations are applied: ∇·(k∇ϕ) = 0 in Ωp, (3a) ∂ϕ ∂np =−u·non Γ, (3b) ϕ=ϕDon ∂Ωp\Γ (3c) where ϕis a potential function, kis the permeability parameter and np=−n is the normal vector pointing outside of the porous domain Ωp. Equations (2c) and (3b) are compatibility conditions between the velocities and pressures in Ωfand Ωpon the interface Γ. Equations (2d) and (3c) are Dirichlet boundary conditions on the boundary ∂Ω\Γ. The implementation of any other kind of boundary conditions on ∂Ω\Γ does not add any difficulty as in standard FEM. 16 Figure 6: Domain Ω under consideration 2.1. Weak form The weak form of the problem, derived by applying integration by parts and the proper boundary conditions, is the following: Find u,pand ϕsuch that u=uDon ∂Ωf\Γand:                  ZΩf ν∇u:∇vdΩ−ZΩf p∇·vdΩ +ZΓ ϕv·ndΓ + ZΓ α(u·t)(v·t)dΓ = (f,v)Ωf∀v −(∇·u,q)Ωf= 0 ∀q (4) The weak form of the Darcy equation is the following: Find ϕ,usuch that ϕ=ϕDon ∂Ωp\Γand: ZΩp k∇v·∇ϕdΩ + ZΓ vu·ndΓ = 0 ∀v(5) 17 Because we are interested in solving the coupled problem, u,pand ϕ must be a solution of both equations (4) and (5) at the same time. 3. Discretization 3.1. Element choice For the discretization of the problem, two different types of elements have been chosen: elements type ”mini”, which consists on a triangular element with 4 nodes (one node in each vertex and another node in the center of the element) and standard triangular elements. The ”mini” element will be used to discretize the Stokes’ velocity in the domain Ωf, while the standard triangular element will be used for all the other variables: pressure (p) in Ωf and potential (ϕ) in Ωp. The nodes for the reference element are shown in figure 7. The ”mini” element satisfies the LBB condition (or inf-sup condition), which must hold in order to guarantee the existence of a stable finite element approximate solution (uh,ph) to the steady stokes problem, see for instance Donea and Huerta[2]. For these element types the shape functions for the standard element are N(ξ, η) = ï1−(ξ+η), ξ, η ò(6) The linear approximation with a bubble can be organized to interpolate 18 (0,0) (1,0) (0,1) (1/3,1/3) (0,0) (1,0) (0,1) Figure 7: Elements in the reference element space of coordinates. the two components of the velocity as Nu(ξ, η) =    1−(ξ+η) 0 ξ0η0 27ξη(1 −ξ−η) 0 0 1 −(ξ+η) 0 ξ0η0 27ξη(1 −ξ−η)    (7) The shape functions (6) and (7) are expressed in the reference element space of coordinates. All the integration process takes place in the reference element, and then the results are transformed to the physical space of coordinates using the isoparametric transformation. This is the standard process in FEM. Now, at each element in the Stokes domain the velocity is approximated as 19 u(x, y)≈uh(ξ, η) = ˜ Nu(ξ, η)u(8) with u=ïu1 x, u1 y, u2 x, u2 y, u3 x, u3 y, u4 x, u4 yòT where (ui x, ui y), i = 1..3 are the nodal values at the 3 vertexes and (u4 x, u4 y) are the bubble function coefficients. Other variables are approximated with standard linear basis functions, that is, p(x, y)≈ph(ξ, η) = N(ξ, η)p=N(ξ, η)        p1 p2 p3        (9) ϕ(x, y)≈ϕh(ξ, η) = N(ξ, η)ϕ=N(ξ, η)        ϕ1 ϕ2 ϕ3        (10) where piand ϕi,i= 1..3 are the nodal values of pand ϕrespectively. 3.2. Linear system of equations Replacing (8), (9),(10) and v=Nuvin (4) the following system of equations is obtained: îK+αMΓ tóu+Gp+BΓ fϕ=f(11) 20 where K,MΓ t,Gand BΓ fare obtained assembling the following elemental matrices, Ke=ZΩf∩Ωeñ∂Nu ∂x ôT∂Nu ∂x +ñ∂Nu ∂y ôT∂Nu ∂y dxdy MΓe=ZΓ∩Ωe [Nu]TtNudxdy BΓe f, =ZΓ∩Ωe [Nu]TnNdxdy (12) for each element Ωeintersecting Ωf In the same way, replacing (8), (10) and v=Nuvin (5), the following system of equations is obtained: Dϕ+BI pu= 0 (13) where matrices Dand BΓ pare obtained by assembly of De=ZΩp∩Ωeñ∂N ∂x ôT∂N ∂x +ñ∂N ∂y ôT∂N ∂y dxdy BΓe p=ZΓ∩Ωe NTnNudxdy (14) for each element Ωeintersecting Ωp. The systems of equations (11) and (13) can be expressed together and solved at the same time, if the matrices are rearranged in the following way:        K+αMΓ tG BΓ f GT0 0 BΓ p0D               u p ϕ        =        ff 0 0        (15) Now this is the system that has to be solved in order to find the solution to the coupled problem. 21 3.3. Implementation 3.3.1. Element integration In the levelset method, the element integration must be changed slightly when it comes to a cut element. The function that computes the integrals over the element has two different steps. At first, the function determines if the element is cut by the interface or if it is in the Stokes domain or in the Darcy’s domain . To do so the nodal values of the levelset function (see section 1.2.1) at the element is used: If all the nodes have a positive value of the levelset function, then the element is strictly in the Stokes domain, and the matrices corresponding to the Stokes equations are computed. If all the nodes have a negative value of the levelset function, then the element is strictly in the Darcy domain, and the matrices corresponding to the Darcy equation are computed. In this case, for simplicity, an auxiliary levelset function is defined, that has the negative value of the original levelset function. This way, if all the nodes have a positive value of this auxiliary levelset function, they will have a negative value of the original levelset function, and therefore, the element will be in the Darcy domain. If there is a node in the element that has sign different to the other nodes, this means that the element is cut by the interface, and the matrices that come from conditions on the interface have to be computed. 22 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 8: Elements marked with a blue star are in the Stokes domain, while elements with a red star are in the Darcy domain. Elements not marked must be considered as Stokes and Darcy elements at the same time. The elements cut by the interface have to be considered as elements in the Stokes domain and in the Darcy domain at the same time. Because of that, the matrices corresponding to the Stokes equations and the Darcy equations have to be computed too. Figure 8illustrates the process described before, for the levelset funtion in figure 4. Elements that are marked with a blue star are elements strictly in the Stokes domain, and elements marked with a red star are strictly in the Darcy domain. Elements that have not any mark must be considered as Stokes and Darcy domain elements at the same time. For the standard elements (that is, the ones that are not cut by the interface) the integration is computed as usual. In this case, a Gauss integration 23 0 0.2 0.4 0.6 0.8 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 9: Gauss quadrature used to integrate the elements. with 7 points is used, the coordinates of which are shown in figure 9. For the integration on the cut elements, the process has several intermediate steps. The first step consist on obtaining the coordinates of the points where the interface cuts the element boundary. For linear triangular elements, there will always be two cut points. The segment joining these two points divides the element in two regions, with two possible situations, as shown in figure 10: the ”positive” area (that is, the area that has a positive value of the levelset function) is similar to a triangle or a quadrilateral. The second step consists on obtaining the integration points. To do so, the element is divided into three triangles as it is shown in figure 11 and a quadrature is applied to each triangle, what gives the integration points for the whole element. The integration points are ordered in a way that the first points of the list are the ones corresponding to the ”positive” area and the ones corresponding to the ”negative” area are in the bottom of the list. Only the integration points for the positive area are considered. Recall that −Φ is considered for Darcy. The third and last step is computing the shape functions on the integra24 5. Numerical application In this section, a more realistic example. A rectangular domain is considered, that can model a horizontal section of a river, with two islands of a porous material in the middle of it. Dirichlet boundary conditions are applied: zero velocity boundary conditions are applied on the bank of the river, to model the no-slip boundary condition that is usually applied in these cases for a viscous fluid. On the inlet section a prescribed velocity will be applied, so the discharge will be determined. A parabolic velocity profile has been applied, because it is the profile that would have in a channel without perturbations. In this case, because of continuity, the discharge on each section of the river would the same, so the discharge on the outlet section would also be determined. This is true if the boundaries are far enough from the porous domain. Figure 13 shows an scheme of the problem, and its boundary conditions. In figures 14 and 15 the horizontal and vertical components of the velocity in the Stokes domain are shown. Figure 15 shows that, because of the no-slip boundary condition, the vertical velocity is higher in the middle of the river, and it gets lower further of the middle, as expected. Figure 14 shows that the horizontal velocity is almost null except near the porous medium. That is because the velocity in the porous medium is lower than in the free flow domain, and the flow field has to adapt to this fact. 31 V=(0,0) V=(0,0) V=(0,Vd) V=(0,Vd) Figure 13: Scheme of the numerical application, with its boundary conditions −0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 vx −0.15 −0.1 −0.05 0 0.05 0.1 Figure 14: Velocity field for the problem described in section 5 32 −0.5 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 vy computed −0.25 −0.2 −0.15 −0.1 −0.05 0 Figure 15: Velocity field for the problem described in section 5 6. Possibilities of the model 6.1. 3D Extension The extension of the model to three dimensions is conceptually easy, following the same principles that have been used to the discretization and integration in 2D, and have been explained in sections 1.2.1,3.3.1 and 3.3.2. Section 1.2.1 describes the concept of a levelset function, that for the case of 3D, would have the following form: Φ(x) = f(x, y, z)(21) and would divide the 3D space in two regions, having different sign of such function. 33 For the element integration described in section 3.3.1, all the concepts can be extended to 3D, taking into account that the elements would be tetrahedrons instead of triangles. The two cases for the cut elements in 2D are the same in 3D, as well the way to divide the cut element in sub-regions to obtain the integration points. For the interface integrals, it will be necessary to take into account that the interface would be a surface instead of a line. Because of that, the parametrization that have been described in 3.3.2 would be a 2D parametrization, and the Gauss integration points would be the corresponding to a 2D quadrature. Apart from these considerations, all the other concepts and processes would be the same for 3D that for 2D. 6.2. Applications The levelset method that has been proposed in this work is suitable to solve flow interaction between surface and ground water. In Discacciati et al. [3] a groundwater flow model is proposed, with the same equations used in this work. In this case, we propose to solve the interaction between surface and groundwater flow with a single iteration using a levelset method, which would be translated into a save in time and computational effort. Erosion problems are a good example of an application of the advantages of the levelset method. Because of the erosion phenomenon, the interface between the porous medium and the fluid will change with time, so the levelset method is a good choice against other methods because it avoids 34 the need to redo the mesh (the advantages of the levelset method in moving interface problems are explained in 1.2.1). In Cottereau and D´ıez [4] the authors use a levelset method to solve the erosion of a bridge pier. To properly solve the erosion problem, it in not enough with the model that we have proposed, because a model for the erosion is needed. Apart from that, all the numerical considerations regarding to the numerical discretization and integration would be the same. 7. Conclusions Along this work, the author has proposed a levelset based solving strategy to compute the solution of Darcy and Stokes coupled problems, which is more efficient than the standard FEM when these problems present a complex geometry or in the case of moving interfaces. With the standard convergence test, it has been proved that the model tends to the solution after a few refinements of the mesh, as it is expected. This guarantees that the model will compute a good solution when applied to a problem that cannot be solved analytically. After that, a problem with no analytical solution has been solved. After all the process, it has been proved in this work that the levelset based strategy to solve Darcy and Stokes coupled problems avoids several disadvantages that are presented with the standard FEM strategy, giving the same results without a significant increase in the computation cost nor the error of the solution. 35 Bibliography [1] Fries TP, Belytschko T. The extended/generalized finite element method: An overview of the method and its applications. Int. J. Numer. Meth. Engng 2010; 84:253304 [2] Donea J. and Huerta A., Finite Element Methods for Flow Problems, JOHN WILEY & SONS, INC., Englewood Cliffs, NJ. Fluid mechanics. [3] Discacciati M. et al. Mathematical and numerical models for coupling surface and groundwater flows Applied Numerical Mathematics 2002; 43:5774 [4] Cottereau R. and D´ıez P. Numerical modeling of erosion using an improvement of the extended finite element method. EJECE 15/2011. Erosion in geomaterials, pages 1187 to 1206. 36