Problems with mixed elements at Multiphasic flow in porous media

variational formulation
1

Hi,

I´m studying multiphase flow in porous media using Brinkman equation:

image

image

image

I have a problem with solving this system of equations using mixed elements, more specific in the saturation equation (third equation).

The saturation must be between 0 <= Sw <=1 in all domain and due to numerical approximations it’s blowing up for higher values. This is the formulation used:

V = VectorElement(“P”, mesh.ufl_cell(), 2) # velocity element
Q = FiniteElement(“P”, mesh.ufl_cell(), 1) # pressure element
DG = FiniteElement(“DG”, mesh.ufl_cell(), order - 1) # saturation element
Element = MixedElement([V, Q, DG])
mixed_space = FunctionSpace(mesh, Element)

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)

L = L1 + L2 + L3 # Sum of all variational forms

J = derivative(L, U, dU) # Jacobian

problem = NonlinearVariationalProblem(L, U, bcs=bcs, J=J)
solver = NonlinearVariationalSolver(problem)

while t < T:
t += float(dt)
U0.assign(U)
solver.solve()

My attempt to solve this issue:

After solving the non-linear problem, enforce the saturation between 0 <= Sw <=1. To do so, I need to access de saturation space and verify the saturation values. However, I don’t know how to do this with mixed elements. If this task weren’t in mixed elements, I would do this:

S1.vector()[S.vector() < 0] = 0

S1.vector()[S.vector() > 1] = 1

Where S1 would be a saturation function of DG space

So, how can I do something similar with mixed elements or there is a better solution to enforce 0 <= Sw <=1?

Thanks