Hi, I have a simple problem. Im trying to solve a standard incomp. fluid flow governed by Navier-Stokes equations. The DG method I have implemented gives me the right results on several examples. But in the following domain, just the zero solution is found.
Here is the definition of domain and mesh, together with boundary description
from dolfin import *
import mshr
dim = 2
# Define Domain
###############
SquareSize = 5.0
LeftDownCorner = Point(2.48,1.0)
RightTopCorner = Point(2.52,4.0)
def closeToMid(x):
dx = RightTopCorner[0] - x[0]
dy = RightTopCorner[1] - x[1]
divX = RightTopCorner[0] - LeftDownCorner[0]
divY = RightTopCorner[1] - LeftDownCorner[1]
if ((dx >= -DOLFIN_EPS and dx <= divX+DOLFIN_EPS) and (dy >= -DOLFIN_EPS and dy <= divY+DOLFIN_EPS)):
return True
else:
return False
def inlet(x):
if ( near(x[1],0.0) and ( (x[0]>=2.0 and x[0]<=3.0) ) ):
return True
else:
return False
def outlet(x):
if ( near(x[1],SquareSize) and (x[0]>=2.0 and x[0]<=3.0) ):
return True
else:
return False
def noslip(x):
if (not outlet(x) and not inlet(x)) and ( near(x[0],0.0) or near(x[0],SquareSize) or near(x[1],0.0) or near(x[1],SquareSize) or closeToMid(x) ):
return True
else:
return False
geometry = mshr.Rectangle(Point(0.0,0.0), Point(SquareSize, SquareSize)) - mshr.Rectangle(LeftDownCorner,RightTopCorner)
# Build mesh
############
mesh = mshr.generate_mesh(geometry, 20)
def onMid(x,bndry):
return bndry and closeToMid(x)
def inletV(x,bndry):
return bndry and inlet(x)
def outletV(x,bndry):
return bndry and outlet(x)
def noslipV(x,bndry):
return bndry and noslip(x)
# Boundary description
######################
bndry = FacetFunction("size_t", mesh)
for f in facets(mesh):
mp = f.midpoint()
if inlet(mp): bndry[f] = 1 # inflow
elif outlet(mp): bndry[f] = 2 # outflow
elif noslip(mp): bndry[f] = 3 # walls
plot(bndry,interactive = True)
You can run the previous code to see how the domain, inlet, outlet and noslip look like. My hope is that just there is some problem. The result looks like the inlet condition is not fulfiled.
The following is the rest of the code, where are among others defined DirichletBC's and the whole DG method.
# Function spaces
V = VectorFunctionSpace(mesh, 'DG', 2)
Q = FunctionSpace(mesh, 'DG', 1)
W = MixedFunctionSpace([V,Q])
zeroD = Expression(("0.0", "0.0"), t = 0.0)
u_IN = Expression(("0.0","100*(x[0]-2.0)*(3.0-x[0])"))
bc_inV = DirichletBC(W.sub(0), u_IN, inletV,'pointwise')
bc_noslipV = DirichletBC(W.sub(0), zeroD, noslipV, 'pointwise')
bcs = [bc_inV, bc_noslipV]
# Unknown and test functions
w = Function(W)
(u, p) = split(w)
(v, q) = TestFunctions(W)
u_1 = Function(V) # velocity from previous time step
u0 = Expression(("0.0", "0.0"))
u_1.assign(interpolate(u0,V))
# Problem parameters
T_end = 3.14
dt = 0.2
t = dt
sigma = 10.0
mu = 1.36*1.e-4
# Variational formulation
n = FacetNormal(mesh)
ds = Measure("ds", subdomain_data = bndry)
I = Identity(u.geometric_dimension())
D = mu*(grad(u) + grad(u).T)
F = D
def a(u,v):
M = inner(F,grad(v))*dx - inner(avg(F)*n('+'),jump(v))*dS
return M
def J(u,v):
M = sigma*inner(jump(u),jump(v))*dS
return M
def b(p,v):
M = -p*div(v)*dx + avg(p)*dot(jump(v),n('+'))*dS
return M
def c(u,w,v):
P = avg(dot(w,n))
H = conditional(P < 0.0, dot(u('+'),w('+')), dot(u('-'),w('-')))
M = -0.5*inner(grad(v)*u,w)*dx + inner(0.5*H*n('+'),jump(v))*dS
return M
T = (1/dt)*inner(u - u_1,v)*dx + a(u,v) + J(u,v) + b(p,v) + b(q,u) + c(u,u_1,v)
# Solution
t = dt
plt = plot(u, mesh = mesh, mode="color", interactive = True, wireframe = False)
pltP = plot(p, mesh = mesh, mode="color", interactive = True, wireframe = False)
while t<=T_end:
solve(T==0,w,bcs)
(u,p) = w.split(True)
u_1.assign(u)
t += float(dt)
plt.plot(u)
pltP.plot(p)
I'm sorry for such a long code, but in order to provide the whole information I kind of didn't want to cut it more.
So, if you can please help me find the reason why the zero solution is computed, even though there is a nonzero inlet flow, I will be very thanhkfull.
