I have been stuck in a problem for weeks and I could not figure out what is wrong with the formulation which I adopted. I am working with Brinkman equation for multi-phase flow porous media. Basically, Brinkman equation is a combination of Stokes e Darcy:
- \mu \nabla ^{2} \textbf u + \frac {\textbf u} {\textbf{K} \lambda (S)} + \nabla p = 0 \: \: \: \:\: \:\: \:eq(1)
\nabla \cdot \mathbf{u} = 0 \: \: \: \:\: \:\: \:\: \: \: \:\: \:\: \:\: \:\: \:\: \:\: \:\: \:\: \:\: \:\: \:\: \:\: \:\: \:\: \:\: \:eq(2)
\frac{\partial S}{\partial t} = \nabla \cdot\left ( f(S)) \mathbf{u} \right ) \:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\: eq(3)
where f(S), \lambda(S) are functions of saturation, \textbf{K} is the permeability tensor and \textbf{u} , p and S are the velocity, pressure and saturation respectively.
I based on my code to modeling Brinkman multiphase flow at this Darcy multiphase example:
https://github.com/hnarayanan/porous-flow/blob/master/src/python/two-phase-flow/two-phase-flow.py
I implemented some modifications in the code above because I need elements that can deal with Stokes and Darcy simultaneously because there is the porous media and a hole with free flow inside this domain, like the figure 1. According with literature, I can use for velocity and pressure the following elements: Taylor-Hood, MINI, Stab P1-P1 and Stab P2-P2. For the Saturation is recommended the DG formulation.
Figure 1
In my code (below) I have modified the reference code (GitHub) elements of velocity and pressure elements from Raviart-Thomas and Discontinuous Galerkin to Taylor-Hood. I kept the same formulation for Saturation, GD. However it is not working properly, with irregularities in the final results on the Saturation as represented in the figure 2. There are elements with S = 1 surrounded by cells with S < 1, which is physically impossible. Furthermore, I tested the other elements mentioned above with no success.
Figure 2
The variational formulation of Brinkman problem is different from the Darcy by adding the Stokes term (first term in eq (5)):
\int_{\Omega}^{}\nabla \mathbf{u} \cdot q dx= 0 \:\: \:\:\:\:\:\:\:\: eq(4)
\int_{\Omega}^{} \mu \nabla \textbf u \cdot \nabla v \: dx + \int_{\Omega}^{}\frac {\mathbf u} {\textbf{K} \lambda (S)} \cdot v \: dx - \int_{\Omega}^{}p \cdot \nabla p\:dx + \int_{\Gamma}^{}n \cdot vp \: ds = 0 \:\: \:\:\:\:\:\:\:\: eq(5)
As mentioned above, the Saturation equation is the same in both codes. The variational formulation is:
\int_{\Omega}^{} \phi S_{}^{n+1}\:r \:dx -\int_{\Omega}^{} \phi S_{}^{n}\:r \:dx + \Delta t\left \{ -\int_{\Omega}^{}\nabla r f(S_{}^{mid}) \mathbf{u} \:dx +\sum_{E \: \epsilon \: \mathcal{E} }^{} \int_{\partial \Gamma}^{} n \cdot \mathbf{u}_{}^{n+1} f(S_{}^{mid})\:dx \right \} + \Delta t \int_{\Gamma in}^{}n \cdot \mathbf{u}_{}^{n+1} S_{}^{mid}\:dx = 0 \:\: \:\:\:\:\:\:\:\: eq(6)
where q, v, r are test functions, S_{}^{mid} = 0.5(S_{}^{n+1} + S_{}^{n}) and n represents the time step.
My questions are: Can I just modify the function elements from RT and DG to Taylor-Hood without any modification in the variational formulation? If it is necessary, which modification should I do? Does anyone have a clue about what is wrong?
I really appreciate any help. Thanks
My code:
from fenics import *
# Relative viscosity of water w.r.t. crude oil
mu_rel = 1
Kinv = Constant(1 / (1e-6))
# Time step
dt = Constant(0.1)
sbar = Constant(1)
pin = 6894.76 # 1.0
pout = 0.0
mu = 1
t = 0.0
T = 400 * float(dt)
Nx = 20
Ny = 20
# Total mobility
def lmbdainv(s):
return 1.0 / ((1.0 / mu_rel) * s ** 2 + (1.0 - s) ** 2)
# Fractional flow function
def F(s):
return s ** 2 / (s ** 2 + mu_rel * (1.0 - s) ** 2)
class vugg(SubDomain):
def inside(self, x, on_boundary):
return between(x[1], (0.4, 0.6)) and between(x[0], (0.2, 0.4))
mesh = UnitSquareMesh(Nx, Ny)
n = FacetNormal(mesh)
V = VectorElement("P", mesh.ufl_cell(), 2) # velocity element
Q = FiniteElement("P", mesh.ufl_cell(), 1) # pressure element
DG = FiniteElement("DG", mesh.ufl_cell(), 0)
Element = MixedElement([V, Q, DG])
mixed_space = FunctionSpace(mesh, Element)
# Function spaces and functions
V = TestFunction(mixed_space)
dU = TrialFunction(mixed_space)
U = Function(mixed_space)
U0 = Function(mixed_space)
v, q, r = split(V)
u, p, s = split(U)
u0, p0, s0 = split(U0)
s_mid = 0.5 * (s0 + s)
U_space = mixed_space.sub(0).collapse()
P_space = mixed_space.sub(1).collapse()
S_space = mixed_space.sub(2).collapse()
obstacle = vugg()
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
obstacle.mark(domains, 1)
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
left = AutoSubDomain(lambda x: near(x[0], 0.0))
right = AutoSubDomain(lambda x: near(x[0], 1.0))
bottom = AutoSubDomain(lambda x: near(x[1], 0.0))
top = AutoSubDomain(lambda x: near(x[1], 1.0))
# Define boundary markers
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
dx = Measure("dx", domain=mesh, subdomain_data=domains)
bc3 = DirichletBC(mixed_space.sub(0).sub(1), Constant(0.0), boundaries, 2)
bc4 = DirichletBC(mixed_space.sub(0).sub(1), Constant(0.0), boundaries, 4)
bcs = [bc3, bc4]
L1 = (
inner(v, lmbdainv(s_mid) * Kinv * u) * dx(0)
- div(v) * p * dx(0)
+ (mu * inner(grad(u), grad(v))) * dx(1)
- div(v) * p * dx(1)
+ pin * dot(v, n) * ds(1)
+ pout * dot(v, n) * ds(3)
)
L2 = q * div(u) * dx(0) + q * div(u) * dx(1)
un = 0.5 * (inner(u0, n) + sqrt(inner(u0, n) * inner(u0, n)))
un_h = 0.5 * (inner(u0, n) - sqrt(inner(u0, n) * inner(u0, n)))
stabilisation = dt("+") * inner(
jump(r), un("+") * F(s_mid)("+") - un("-") * F(s_mid)("-")
) * dS + dt * r * un_h * sbar * ds(1)
L3 = (
r * (s - s0) * dx(0)
+ r * (s - s0) * dx(1)
- dt * inner(grad(r), F(s_mid) * u) * dx(0)
- dt * inner(grad(r), F(s_mid) * u) * dx(1)
+ dt * r * F(s_mid) * un * ds
+ stabilisation
)
L = L1 + L2 + L3
J = derivative(L, U, dU)
snes_solver_parameters = {
"nonlinear_solver": "snes",
"snes_solver": {
"linear_solver": "lu",
"maximum_iterations": 100,
"absolute_tolerance": 1e-7,
"report": True,
"error_on_nonconvergence": True,
},
}
problem = NonlinearVariationalProblem(L, U, bcs=bcs, J=J)
lower = Function(mixed_space)
upper = Function(mixed_space)
U_lower = Function(U_space)
U_lower.vector()[:] = -10000
U_upper = Function(U_space)
U_upper.vector()[:] = +10000
P_lower = Function(P_space)
P_lower.vector()[:] = -100000
P_upper = Function(P_space)
P_upper.vector()[:] = 100000
S_lower = Function(S_space)
S_lower.vector()[:] = 0
S_upper = Function(S_space)
S_upper.vector()[:] = 1
fa = FunctionAssigner(
mixed_space,
[U_space, P_space, S_space],
)
fa.assign(lower, [U_lower, P_lower, S_lower])
fa.assign(upper, [U_upper, P_upper, S_upper])
problem.set_bounds(lower, upper)
solver = NonlinearVariationalSolver(problem)
solver.parameters.update(snes_solver_parameters)
s_file = File("teste/saturation.pvd")
while t < T:
t += float(dt)
U0.assign(U)
solver.solve()
u, p, s = U.split()
s_file << s