dolfin team mailing list archive

[Melanie Jahny <question157341@xxxxxxxxxxxxxxxxxxxxx>]

New question #157341 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/157341

I'd like to compute a flow about an inner boundary.
My testproblem is in a unit square. I solve stokes-equation there and want to
compute the flow in the middle of the mesh.
I start with inflow -1 at the upper boundary(y=1).
The flow at the bottom boundary ( y=0) is also
-1, but when I compute the flow at y = 0.5, the flow value becomes very big
and it seems that it depends much on the discretization of the mesh.
Can anybody help me and find my mistake???
Here is the code:

mesh = Rectangle(0,0,1,1,10, 10)
mesh.order()

V = VectorFunctionSpace(mesh, "CG", 2)
W = FunctionSpace(mesh, "CG", 1) 

VW = V*W

(u,p) = TrialFunctions(VW)
(phi,psi) = TestFunctions(VW)


# --------------
# boundary conditions
def boundary_noslip(x, on_boundary):
    return on_boundary     
    
noslip_const = Constant((0.0, 0.0)) 
noslip = DirichletBC(VW.sub(0), noslip_const, boundary_noslip)

def boundary_in_left(x, on_boundary):
    return x[1] == 1.0 and on_boundary
    
u_in_left = Expression(("0.0", "-1.0")) 
inflow_left  = DirichletBC(VW.sub(0), u_in_left, boundary_in_left)

def boundary_pressure_down(x, on_boundary):    
    return x[1] == 0.0 and on_boundary
    
p_const_down = Constant(0.0)
p_down = DirichletBC(VW.sub(1), p_const_down, boundary_pressure_down)

bcs = [noslip, inflow_left, p_down] 

# stokes equation
f = Constant((0, 0))
a = (inner(grad(phi), grad(u)) - div(phi)*p + psi*div(u))*dx
L = inner(phi, f)*dx

problem = VariationalProblem(a, L, bcs)
(u,p) = problem.solve().split(True) 

# inner boundary flow
normal = FacetNormal(mesh)
boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)

bc_mid = compile_subdomains('x[1] == 0.5')
bc_mid.mark(boundary_parts,0) 

flow = (u('+')[0]*normal('+')[0] - u('-')[0]*normal('-')[0] + (u('+')[1]*normal('+')[1]) - \
                (u('-')[1]*normal('-')[1]))*dS(0)
                
integral_flow = assemble(flow, mesh=mesh, interior_facet_domains = boundary_parts)

-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.