Hi everyone,
I am trying to solve a fluid flow problem using MixedFunctionSpace, however, I am not sure how to create trial function and test function for the problem. I try to follow various methods in the forum but always got different errors. Can anyone suggest where in the code am I wrong?
My problem looks like this,
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
mesh = RectangleMesh(Point(-0.08, -0.02), Point(0.08, 0.02), 160, 80)
#plot(mesh)
# define meshing parameters
x1 = -0.08
y1 = -0.02
x2 = 0.08
y2 = 0.01
x3 = -0.08
y3 = 0.01
x4 = 0.08
y4 = 0.02
regions = MeshFunction('size_t', mesh, mesh.topology().dim())
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
class Fluid(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.01
fluid = Fluid()
class Solid(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.01
solid = Solid()
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.02)
top = Top()
class LeftF(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], -0.08) and between(x[1], (0.01, 0.02))
leftf = LeftF()
class RightF(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.08) and between(x[1], (0.01, 0.02))
rightf = RightF()
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.01)
interface = Interface()
class LeftS(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], -0.08) and between(x[1], (-0.02, 0.01))
lefts = LeftS()
class RightS(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.08) and between(x[1], (-0.02, 0.01))
rights = RightS()
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], -0.02)
bottom = Bottom()
regions.set_all(0)
fluid.mark(regions, 1)
solid.mark(regions, 2)
boundaries.set_all(0)
top.mark(boundaries, 1)
leftf.mark(boundaries, 2)
rightf.mark(boundaries, 3)
interface.mark(boundaries, 4)
lefts.mark(boundaries, 5)
rights.mark(boundaries, 6)
bottom.mark(boundaries, 7)
# define a submesh, composed of region 0 and 1
fluid = MeshFunction('size_t', mesh, mesh.topology().dim())
fluid.set_all(0)
fluid.array()[regions.array() == 1] = 1
submeshfluid = SubMesh(mesh, fluid, 1)
solid = MeshFunction('size_t', mesh, mesh.topology().dim())
solid.set_all(0)
solid.array()[regions.array() == 2] = 2
submeshsolid = SubMesh(mesh, solid, 2)
# define new meshfunctions on this submesh with the same values as their original mesh
fluid_regions = MeshFunction('size_t', submeshfluid, submeshfluid.topology().dim())
fluid_regions.set_all(0)
fluid_boundaries = MeshFunction('size_t', submeshfluid, submeshfluid.topology().dim() - 1)
fluid_boundaries.set_all(0)
solid_regions = MeshFunction('size_t', submeshsolid, submeshsolid.topology().dim())
solid_regions.set_all(0)
solid_boundaries = MeshFunction('size_t', submeshsolid, submeshsolid.topology().dim() - 1)
solid_boundaries.set_all(0)
top.mark(fluid_boundaries, 1)
leftf.mark(fluid_boundaries, 2)
rightf.mark(fluid_boundaries, 3)
#interface.mark(fluid_boundaries, 4)
interface.mark(solid_boundaries, 4)
lefts.mark(solid_boundaries, 5)
rights.mark(solid_boundaries, 6)
bottom.mark(solid_boundaries, 7)
File("2SUBMESHTEST/Results/fluid_boundaries.pvd") << fluid_boundaries
File("2SUBMESHTEST/Results/solid_boundaries.pvd") << solid_boundaries
dsf = Measure("ds")[fluid_boundaries]
dss = Measure("ds")[solid_boundaries]
UF = VectorFunctionSpace(submeshfluid, 'P', 2)
PF = FunctionSpace(submeshfluid, 'P', 1)
TS = FunctionSpace(submeshsolid, 'P', 1)
W = MixedFunctionSpace(UF, PF, TS)
(v,q,z) = TestFunctions(W)
w = Function(W)
u = w.sub(0)
p = w.sub(1)
dxf = Measure("dx")[fluid_regions]
dxs = Measure("dx")[solid_regions]
bcs = [DirichletBC(W.sub_space(0).sub(0), Constant(0), top),
DirichletBC(W.sub_space(0).sub(1), Constant(0), top),
DirichletBC(W.sub_space(0).sub(0), Constant(0), interface),
DirichletBC(W.sub_space(0).sub(1), Constant(0), interface),
DirichletBC(W.sub_space(1), Constant(80), leftf),
DirichletBC(W.sub_space(1), Constant(0), rightf)]
# Surface normal
nf = FacetNormal(submeshfluid)
ns = FacetNormal(submeshsolid)
#############################
## for solving steady flow ##
#############################
# Define parameters used in the weak form
mu = 0.89
rho = 1000
k = 0.6071
cp = 4220
alp = k/rho/cp
muS = 1
rhoS = 3000
kS = 1.5
cpS = 220
alpS = kS/rhoS/cpS
Fl = 0.00005
# Weak form of the steady state fluid problem
F = rho*inner(grad(u)*u, v)*dxf \
+ mu*inner(grad(u), grad(v))*dxf \
+ inner(grad(p),v)*dxf \
- div(u)*q*dxf
# Jacobian matrix to be supplied to the non linear solver
J = derivative(F, w) # Jacobian
# Call the non'linear solver
solve(F == 0, w, bcs, J=J)
#############################
# Plot solution
plt.figure()
a = plot(u)
plt.colorbar(a)
plt.figure()
b = plot(p)
plt.colorbar(b)