Hi!
I have a working example of a simple incompressible Navier-Stokes flow and I have a working example of a single phase flow in porous medium. Now I would like to make a couple system, where some part of the domain is fluid flow defined by NS and some part is a porous medium flow.
Is there any FEniCS example out there, on how to couple the two models?
Thank you!
Can you be more specific in the porous medium? Are you just solving -k\nabla p = v and \nabla \cdot v= 0 there? Sounds like you can just split up in two subdomains. And define NS is one and Darcy in the other. You might not get exactly the correct boundary conditions if you use continuous elements though. What did you try so far?
You can read these on how to couple. The code is probably outdated, but the idea is the same:
Thank you for you response.
I am exactly just solving the system you mention. I have tried splitting the domain up into two subdomains, but perhaps my biggest issue is, how to formulate the incompressible Naver-Stokes system in the mixed finite element formulation.
The first thesis you sent - I have already look at it very much, and the FEniCS functions used there are simply too outdated and the code is not complete enough to understand all variables.
I will however take a look at the second thesis - thank you for your recommendation.
You might want to add some of your code to get help as well. For the coupled Navier Stokes system (and not IPCS or other incremental schemes) you might want to look at this 2D Stationary Navier-Stokes equation at @dokken s response. Here you should just add the time dependency to have the NS instead of Stokes. Solving everything at the same time is more convenient than doing projections and corrections both within the fluid and between the fluid and the porous medium.
This might be a good starting point for you:
Here is my code. Under "set settings line 14-30, you can see what works, what “sort of” works and what does not work. I believe my problems are mostly the function space and dimensions needed for the mixed formulation of Navier-Stokes. I tried the stationary Navier-Stokes equation from @dokken, but something is still wrong there… Unde "define weak formulation for Navier-Stokes you can see all my failed attempts. Sigh!
The mesh is a square where the upper half is a domain marked 0 and the lower half is marked 1.
#region: Load libraries----------------------------
from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div, max_value
from progress.bar import Bar
import numpy as np
import mshr
import random
import sys
#import time
#endregion----------------------------------------
#region: set settings------------------------------
#chose one:
#single porous domain (the other is undefined)
run_single_porous_example = False
#two different porous mediums
run_double_porous_example = False
#single Navier-Stokes domain (the other is undefined)
run_single_NS_example = True #sort of work
#two different NS flows
run_double_NS_example = False #do not work
#mixed Navier-Stokes and porous medium
run_mixed_example = False #do not work
#endregion-----------------------------------------
#region: define square mesh ------------------------------
domain = mshr.Rectangle(Point(0.0), Point(1,1))
porous = mshr.Rectangle(Point(0,0),Point(1,0.5))
domain.set_subdomain(1, porous)
mesh = mshr.generate_mesh(domain, 30)
d = mesh.topology().dim()
mf = MeshFunction("size_t",mesh,d,mesh.domains())
dx=Measure('dx', domain=mesh, subdomain_data=mf)
#endregion ----------------------------------------
#region: define boundaries ------------------------
class Top_wall(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 1.0))
class Bottom_wall(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 0.0))
class Mid_Wall(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 0.5))
class Mid_Wall_with_hole(SubDomain):
def inside(self, x, on_boundary):
return (near(x[1], 0.5) and between(x[0], (0, 0.2)) or near(x[1], 0.5) and between(x[0], (0.8, 1)) )
class leftWall_fluid(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], 0.0) and x[1]>0.5-DOLFIN_EPS)
class rightWall_fluid(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], 1.0) and x[1]>0.5-DOLFIN_EPS)
class leftWall_porous(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], 0.0) and x[1]<0.5+DOLFIN_EPS)
class rightWall_porous(SubDomain):
def inside(self, x, on_boundary):
return (near(x[0], 1.0) and x[1]<0.5+DOLFIN_EPS)
#call boundaries
top = Top_wall()
bottom = Bottom_wall()
midwall = Mid_Wall()
midwall_hole = Mid_Wall_with_hole()
leftwall_fluid=leftWall_fluid()
rightwall_fluid=rightWall_fluid()
leftwall_porous=leftWall_porous()
rightwall_porous=rightWall_porous()
#mark boundaries
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1,0)
#names for easyness
Top=1
Bottom=2
middle=3
middle_hole=4
inlet_fluid=5
outlet_fluid=6
inlet_porous=7
outlet_porous=8
top.mark(boundaries, Top)
bottom.mark(boundaries, Bottom)
if run_single_porous_example or run_single_NS_example: #close middle wall
midwall.mark(boundaries,middle)
else: #open middle wall
midwall_hole.mark(boundaries, middle_hole)
leftwall_fluid.mark(boundaries, inlet_fluid)
rightwall_fluid.mark(boundaries, outlet_fluid)
leftwall_porous.mark(boundaries, outlet_porous)
rightwall_porous.mark(boundaries, inlet_porous)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
boundary_file = File("boundaries.pvd")
boundary_file << boundaries
#endregion ----------------------------------------
#region: define parameters ------------------------
t = 0.0; dt = 0.0256; T = 1
force = Constant((0., 0.))
nu = Constant(0.01)
k=0.5
mu=nu
g = Constant(0)
#parameters["ghost_mode"] = "shared_facet"
n = FacetNormal(mesh)
thing = as_vector((-n[1], n[0]))
#endregion----------------------------------------
# #region: define functionspaces and functions-------
# #define functionspaces
V = FunctionSpace(mesh, 'RT', 1) # velocity space, Raviart–Thomas space
Q = FunctionSpace(mesh, 'DG', 0) # pressure space, Discontinuous Lagrange space
U = Function(V)#velocity, current time step
P = Function(Q)#pressure, current solution
# #and mixed functionspace
if run_single_NS_example:
# Ve = VectorElement("P",triangle, 2)
# Qe = FiniteElement("P",triangle, 1)
# X = FunctionSpace(mesh, MixedElement([Ve, Qe]))
P2v = VectorElement("CG",triangle, 2)
P1 = FiniteElement("CG", triangle, 1)
X = FunctionSpace(mesh, MixedElement([P2v, P1]))
else:
RT = FiniteElement("RT", mesh.ufl_cell(), 1)
DG = FiniteElement("DG", mesh.ufl_cell(), 0)
X = FunctionSpace(mesh, RT * DG) #mixed space
# #define function
# xh = Function(X)
# #define trial and test functions
(u, p) = TrialFunctions(X) #for velocity, pressure
(v, q) = TestFunctions(X) #for velocity, pressure
# #endregion-----------------------------------------
#region: define boundary conditions----------------
pd=Constant(1)
bcp1 = DirichletBC(X.sub(1), pd, boundaries, inlet_fluid)
bcp2 = DirichletBC(X.sub(1), pd, boundaries, inlet_porous)
bcu1 = DirichletBC(X.sub(0), (0,0), boundaries, Top)
bcu2 = DirichletBC(X.sub(0), (0,0), boundaries, Bottom)
if run_single_porous_example: #only boundaries for lower domain
bcu3 = DirichletBC(X.sub(0), (0,0), boundaries, middle)
bc = [bcp2, bcu2, bcu3]
if run_double_porous_example: #boundaries for whole domain
bcu4 = DirichletBC(X.sub(0), (0,0), boundaries, middle_hole)
bc = [bcp1, bcp2, bcu1, bcu2, bcu4]
if run_single_NS_example: #only boundaries for upper domain
bcp3 = DirichletBC(X.sub(1), Constant(0), boundaries, outlet_fluid) #zero on outlet for NS
bcu3 = DirichletBC(X.sub(0), (0,0), boundaries, middle)
bc = [bcp1, bcp3, bcu1, bcu3]
if run_double_NS_example:
bcp3 = DirichletBC(X.sub(1), Constant(0), boundaries, outlet_fluid) #zero on outlet for NS
bcp4 = DirichletBC(X.sub(1), Constant(0), boundaries, outlet_porous) #zero on outlet for NS
bcu4 = DirichletBC(X.sub(0), (0,0), boundaries, middle_hole)
bc = [bcp1, bcp2, bcp3, bcp4, bcu1, bcu2, bcu4]
if run_mixed_example:
bcp3 = DirichletBC(X.sub(1), Constant(0), boundaries, outlet_fluid) #zero on outlet for NS
bcu4 = DirichletBC(X.sub(0), (0,0), boundaries, middle_hole)
bc = [bcp1, bcp2, bcp3, bcu1, bcu2, bcu4]
#endregion---------------------------------------
#region: define weak formulation Navier-Stokes------------
# Define symmetric gradient
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*nu*epsilon(u) - p*Identity(len(u))
#attempt1:
# u, p = split(xh)
# F1_ns = dot(dot(u, nabla_grad(u)), v)*dx \
# + inner(sigma(u, p), epsilon(v))*dx \
# + dot(p*n, v)*ds - dot(nu*nabla_grad(u)*n, v)*ds \
# - dot(force, v)*dx
# a1_ns = lhs(F1_ns)
# L1_ns = rhs(F1_ns)
#attempt2:
#u, p = split(xh)
#F1_ns = nu*inner(grad(u), grad(v))*dx(0) + dot(dot(grad(u), u), v)*dx(0) \
# - p*div(v)*dx(0) - q*div(u)*dx(0)
# a1_ns = lhs(F1_ns)
# L1_ns = rhs(F1_ns)
#attempt3:
# a1_ns = nu*inner(grad(u), grad(v))*dx(0) - inner(p, div(v))*dx(0) - inner(q, div(u))*dx(0)
# L1_ns = inner(Constant((0,0)), v)*dx(0)
#attempt4:
dS = Measure("dS", domain=mesh, subdomain_data=boundaries)
F1_ns= 2.0 * mu * inner(epsilon(u), epsilon(v)) * dx(0) \
+ mu*pow(k, Constant(-0.5))*dot(u("+"), thing("+"))*dot(v("+"), thing("+"))*dS(inlet_fluid) \
- dot(force, v) * dx(0) + g*q*dx(0) + p * div(v) * dx(0) + q * div(u) * dx(0)
a1_ns = lhs(F1_ns)
L1_ns = rhs(F1_ns)
#endregion--------------------------------------------
#region: define weak formulation porous media-------------
a = inner(mu/0.7*u,v)*dx(1) #darcy 1
bv = -div(v)*p*dx(1) #darcy 2
bu = -div(u)*q*dx(1) #conservation of mass
F1_porous = a+bv+bu + g*q*dx(1)-dot(v,n)*pd*ds(inlet_porous)
a1_porous = lhs(F1_porous)
L1_porous = rhs(F1_porous)
a2 = inner(mu/0.2*u,v)*dx(0) #darcy 1
bv2 = -div(v)*p*dx(0) #darcy 2
bu2 = -div(u)*q*dx(0) #conservation of mass
F1_porous2 = a2+bv2+bu2 + g*q*dx(0)-dot(v,n)*pd*ds(5)
a1_porous2 = lhs(F1_porous2)
L1_porous2 = rhs(F1_porous2)
#endregion----------------------------------------
#region: define weak formulation to solve--------
if run_single_porous_example:
a1 = a1_porous
L1 = L1_porous
if run_double_porous_example:
a1 = a1_porous + a1_porous2
L1 = L1_porous + L1_porous2
if run_single_NS_example:
a1 = a1_ns
L1 = L1_ns
if run_double_NS_example:
a1 = a1_ns + a1_ns2
L1 = L1_ns + a1_ns2
if run_mixed_example:
a1 = a1_ns + a1_porous
L1 = L1_ns + L1_porous
#endregion-------------------------------------
#region: create save files--------------------------
#get scripts name
name_of_self = sys.argv[0].split('/')[-1].split('.')[0]
#for storing solutions
directory_xdmf='results/'+name_of_self+'/'
xdmffile_velocity = XDMFFile(directory_xdmf+'velocity.xdmf')
xdmffile_pressure = XDMFFile(directory_xdmf+'pressure.xdmf')
xdmffile_temperature = XDMFFile(directory_xdmf+'temperature.xdmf')
xdmffile_velocity.parameters["flush_output"] = True
xdmffile_pressure.parameters["flush_output"] = True
xdmffile_temperature.parameters["flush_output"] = True
#endregion-----------------------------------------
#region: solve problem---------------------------
set_log_active(False)
bar = Bar('Processing', max=T/dt)
xh = Function(X)
u, p = split(xh)
while t +dt < T + DOLFIN_EPS:
dt = min(T-t, dt)
solve(a1 == L1, xh, bc)
U.assign(xh.sub(0, deepcopy=True))
P.assign(xh.sub(1, deepcopy=True))
xdmffile_velocity.write(U, t); xdmffile_pressure.write(P, t)
#update
t += dt
bar.next()
bar.finish()
#endregion-----------------------------------------