Stokes flow with two viscosities

mesh
1

Hi, I am a student trying to solve a Stokes system that has a liquid phase bubble within a solid phase, which I represent with two separate viscosities. I am following the example shown here at the bottom example with varying permeabilities.

My code is here:

# Set domain geometries
L = 400
H = 200
resolution = 32

domain = Rectangle(Point(0.0,0.0), Point(L,H))
liquid = Circle(Point(100, 0.5*H), 50)
solid = domain - liquid

domain.set_subdomain(1, liquid) 
domain.set_subdomain(2, solid) 

# Export mesh as xml file, read the xml file
domain_mesh = generate_mesh(domain, resolution)
File("meshes/twomaterials/mesh.xml") << domain_mesh
mesh = Mesh("meshes/twomaterials/mesh.xml")


# Define subdomains ----------------------------
# Sub domain for no-slip (mark whole boundary, inflow and outflow will overwrite)
class Noslip(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary 

# Sub domain for inflow (left)
class Inflow(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < DOLFIN_EPS and on_boundary

# Sub domain for outflow (right)
class Outflow(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > L - 1e3*DOLFIN_EPS and on_boundary

sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

sub_domains.set_all(3)

noslip = Noslip()
noslip.mark(sub_domains, 0)

inflow = Inflow()
inflow.mark(sub_domains, 1)

outflow = Outflow()
outflow.mark(sub_domains, 2)

file = File("meshes/twomaterials/subdomains.xml")
file << sub_domains

# defining problem space -----------------

# Taylor-Hood function space

VE = VectorElement("DG", mesh.ufl_cell(), 2) # previously was CG, not DG
QE = FiniteElement("DG", mesh.ufl_cell(), 1)
WE = VE * QE

W = FunctionSpace(mesh, WE)
V = FunctionSpace(mesh, VE)
Q = FunctionSpace(mesh, QE)

# boundary conditions-----------------------

# No-slip boundary condition for the liquid pocket
noslip = Constant((0, 0))
bc0 = DirichletBC(W.sub(0), noslip, sub_domains, 0) # subdomain marked 0 (noslip)

# Inflow boundary condition for velocity
YMAX = H
inflow = Expression(("4*(x[1]*(YMAX-x[1]))/(YMAX*YMAX)", "0."), YMAX=YMAX, element = V.ufl_element())
bc1 = DirichletBC(W.sub(0), inflow, sub_domains, 1) # subdomain marked 1 (inflow)

# Boundary condition for pressure at outflow
# x0 = 0
zero = Constant(0)
bc2 = DirichletBC(W.sub(1), zero, sub_domains, 2) # subdomain marked 2 (outflow)

# Collect boundary conditions
bcs = [bc0, bc1, bc2]

# "Define subdomain markers and integration measure"
markers = MeshFunction('size_t', mesh, 2, mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=markers)

# Define an Expression for viscosity 
class Viscosity(UserExpression):
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs)
        self.markers = markers
    def eval_cell(self, values, x, cell): 
        if self.markers[cell.index]==1: # liquid
            values[0]==1 # viscosity in liquid is 1
        else:
            values[0]==3e15 # solid viscosity much higher
    
    def value_shape(self): # from https://fenicsproject.discourse.group/t/problems-interpolating-a-userexpression-and-plotting-it/1303/2
        #return (1,)
        return ()
    
eta = Viscosity(markers)

# Define variational problem

I = Identity(2)
def epsilon(u): # strain rate tensor
  return 0.5*(grad(u) + grad(u).T)

def sigma(u,p): # Cauchy stress tensor
  return 2*eta*epsilon(u) - p*I # from https://fenicsproject.org/pub/tutorial/html/._ftut1008.html 

# Define trial and test functions
u, p = TrialFunctions(W)
v, q = TestFunctions(W)
f = Constant((0, 0)) # set f, body force, to 0
a = inner(sigma(u,p), grad(v))*dx + q*div(u)*dx
L = inner(f, v)*dx

# Compute solution
w = Function(W)
solve(a == L, w, bcs)

(u, p) = w.split(True)

plot(u)
plt.show()

plot(p)
plt.show()

The error I keep getting is shown below. Since there are no values being computed, my plots also come out blank. I have spent many hours trying to figure out where it is coming from but I am not even sure what part of the code is causing this malfunction :frowning_face: Any help would be greatly appreciated!

/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:715: RuntimeWarning: divide by zero encountered in double_scalars
  length = a * (widthu_per_lenu / (self.scale * self.width))
/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:715: RuntimeWarning: invalid value encountered in multiply
  length = a * (widthu_per_lenu / (self.scale * self.width))
/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:767: RuntimeWarning: invalid value encountered in less
  short = np.repeat(length < minsh, 8, axis=1)
/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/quiver.py:780: RuntimeWarning: invalid value encountered in less
  tooshort = length < self.minlength
2

I haven’t tried your code. But you’re testing for equality rather than assigning values in your Viscosity class:

    def eval_cell(self, values, x, cell): 
        if self.markers[cell.index]==1: # liquid
            values[0]=1 # viscosity in liquid is 1
        else:
            values[0]=3e15 # solid viscosity much higher

Note the difference between = and ==.

Also if you want to use a DG finite element, the facet terms are missing from your FE formulation. It’s probably easier to use a Taylor-Hood element.