Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.

Notice

The forum is in read only mode.

Pressure induced flow

More
1 year 1 month ago - 1 year 1 month ago #4658 by Raphael_S
Hy everybody, I am quiet new to NGSolve and so I need some help. I would like to solve the Stokes equation for a pressure induced flow (cf. flow inside a pipe).  In the tutorials there is a velocity induced flow presented. In my case, I would like to apply a pressure difference on the left and right face of my geometry, which are called inlet and outlet. In addition, I added a Neumann BC on both sides. The major problem now is, that the fluid is not moving in the centre of the geometry (clip through the y axis shows this). Furthermore, if the pressure difference is increased, the max. velocity is not rising.  I am vey thankful for every support.

from ngsolve import *
import netgen.gui
from netgen.occ import *
from netgen.geom2d import SplineGeometry
import numpy as np

U = -5*10**3 # velcoity in mm/s
rho = 857.22*10**-3/10**9   # density of the oil // t/mm^3
eta = 0.022474*10**-6       # dynamic viscosity of the oil // N/mm^2*s


nu = (eta/rho)  # kinematic viscosity
geo = Box((0,0,0), (3,1,1))

geo.faces.Min(X).name='inlet'
geo.faces.Max(X).name='outlet'
geo.faces.Min(Y).name='front'
geo.faces.Max(Y).name='back'
geo.faces.Min(Z).name='bot'
geo.faces.Max(Z).name='top'

mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.1))

Draw (mesh)

V = VectorH1(mesh, order=3, dirichlet="front|back|bot|top")
Q = H1(mesh, order=2,dirichlet="inlet|outlet")

u,v = V.TnT()
p,q = Q.TnT()

a = BilinearForm(V)
a += nu*InnerProduct(grad(u), grad(v))*dx

b = BilinearForm(trialspace=V, testspace=Q)
b += div(u)*p*dx+div(v)*q*dx


a.Assemble()
b.Assemble()

mp = BilinearForm(Q)
mp += SymbolicBFI(1e-10*p*q)
mp.Assemble()

f = LinearForm(V)
f+= CoefficientFunction((1,0,0))* v * ds(definedon='outlet')
f += CoefficientFunction((1,0,0))* v * ds(definedon='inlet')
f.Assemble()


g = LinearForm(Q)
g.Assemble()

gfu = GridFunction(V, name="u")
gfp = GridFunction(Q, name="p")

cf = CoefficientFunction([(2/rho) if bc=="inlet" else (0.2/rho) if bc=="outlet"  else 0 for bc in mesh.GetBoundaries()])
gfp.Set(cf, definedon=mesh.Boundaries("inlet|outlet"))

K = BlockMatrix( [ [a.mat, b.mat.T], [b.mat, None] ] )
C = BlockMatrix( [ [a.mat.Inverse(V.FreeDofs()), None], [None, mp.mat.Inverse()] ] )

rhs = BlockVector ( [f.vec, g.vec] )
sol = BlockVector( [gfu.vec, gfp.vec] )

solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, initialize=False)
# density free pressure --> pressure
gfp.vec.data.FV().NumPy()[:] = gfp.vec.data.FV().NumPy() * rho

pressure = gfp
velocity = gfu

Draw(gfu)
Last edit: 1 year 1 month ago by Raphael_S.
More
1 year 1 month ago - 1 year 1 month ago #4680 by janellmenreich
Dear Raphael,

from a mathematical point of view you can not impose a Dirichlet boundary condition on the pressure, certainly not on the pressure gradient, since $p \in L^2_0$. For Stokes equations we can only set Dirichlet boundary conditions on the velocity. The natural boundary conditions (Neumann) arising from the derivation of the weak formulation, are so called force boundary conditions, namely
$ \int_{\Gamma_N} \left( \mathbf{\tau} - p\mathbf{I} \right) n \cdot \mathbf{v}$.

What you can do, is to think of the pressure gradient as a imposed force acting on the fluid in the whole domain, and not only on the boundary. Therefore it suffices to add the imposed force to the linearform
$\int_{\Omega} \mathbf{f} \cdot \mathbf{v}$.

Best,
Jan

###Code###
from ngsolve import *
from netgen.occ import *

U = -5*10**3 # velcoity in mm/s
rho = 857.22*10**-3/10**9 # density of the oil // t/mm^3
eta = 0.022474*10**-6 # dynamic viscosity of the oil // N/mm^2*s
nu = (eta/rho) # kinematic viscosity
imposed_pressure_gradient = CF((-2, 0, 0)) # pressure decreases in x

geo = Box((0, 0, 0), (3, 1, 1))
geo.faces.Min(X).name = 'inlet'
geo.faces.Max(X).name = 'outlet'
geo.faces.Min(Y).name = 'front'
geo.faces.Max(Y).name = 'back'
geo.faces.Min(Z).name = 'bot'
geo.faces.Max(Z).name = 'top'

mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.25))

V = VectorH1(mesh, order=2, dirichlet="front|back|bot|top")
Q = H1(mesh, order=1)

fes = V*Q
(u, p), (v, q) = fes.TnT()

a = BilinearForm(fes)
a += nu*InnerProduct(grad(u), grad(v))*dx
a += -div(u)*p*dx-div(v)*q*dx
a += 1e-8 * p * q * dx
a.Assemble()

f = LinearForm(fes)
f += -imposed_pressure_gradient * v * dx
f.Assemble()

gfu = GridFunction(fes)
velocity = gfu.components[0]
pressure = gfu.components[1]

inv = a.mat.Inverse(freedofs=fes.FreeDofs(), inverse="umfpack")
gfu.vec.data = inv * f.vec

Draw(velocity, mesh, "velocity")
Draw(pressure, mesh, "pressure")
Last edit: 1 year 1 month ago by janellmenreich.
Time to create page: 0.120 seconds