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.

Problems with defining pressure BC for Stokes equation

More
1 year 2 months ago #4643 by Raphael_S
Hy everybody,I try to solve the Stokes equation. If I define a velocity BC (U on surface ‘bot’) everything works fine and the expected result and the calculated result correspond. But, if I define a pressure BC, the results do not correspond to my expectations. On the surface ‘left’ I defined a pressure of 2 MPa and on the right surface I defined a pressure of 0.2 MPa. The velocity on the other surfaces is 0. My problem is, that after solving the Stokes equation, the pressure field on the left and right surface, does not have the value of the BC. So also, the velocity field is quite strange.So, maybe somebody can help me to define the pressure BC correctly.Many thanks to all of you.

Followed you can find my code: 
# V = Periodic(VectorH1(mesh,order=3, dirichlet='bot|top'))  # velocity
Q = H1(mesh,order=2, dirichlet='left|right') # pressure

X = V * Q 

# obtain test and trial functions
u, p = X.TrialFunction()
v, q = X.TestFunction()
# space for solution data
gfu = GridFunction(X)
# boundary condition
#uin = CoefficientFunction((U,0,0))
#gfu.components[0].Set(uin, definedon=mesh.Boundaries("bot"))
cf = CoefficientFunction([(2) if bc=="left" else (0.2) if bc=="right"  else 0 for bc in mesh.GetBoundaries()])
    # Set BC
gfu.components[1].Set(cf, definedon=mesh.Boundaries("left|right"))



nu = (eta/rho)  # kinematic viscosity

SetNumThreads(24) #Anzahl der Kerne zum Parallelisieren
with TaskManager():
    # right hand side of the equation
    f = LinearForm(X)
    f.Assemble()

    # left hand side of the equation
    stokes = nu*InnerProduct(grad(u), grad(v))+div(u)*q+div(v)*p - 1e-10*p*q
    a = BilinearForm(X, symmetric=True)
    a += stokes*dx
    a.Assemble()

    # calculate solution data
    inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")
    res = f.vec.CreateVector()
    res.data = f.vec - a.mat*gfu.vec
    gfu.vec.data += inv_stokes * res

    # calculate drag force
    res.data = a.mat*gfu.vec
    F_drag_j = InnerProduct(res,gfu_help.vec) * rho

# density free pressure --> pressure
gfu.components[1].vec.data.FV().NumPy()[:] = gfu.components[1].vec.data.FV().NumPy() * rho


# pressure solution field
pressure = gfu.components[1]
# Define a coefficient function for val

# Draw velocity
#Draw(Norm(gfu.components[0]), mesh, "velocity",sd=3)   
#
Draw((gfu.components[0]), mesh, "velocity",sd=3)          # Draw velocity
# Draw the pressure
#Draw(pressure,mesh)                                         # Draw pressure
More
1 year 2 months ago #4644 by hvwahl
Hi Raphael_S,

could you please provide a minimal  (non)-working example? I ran the code below and this set the pressure correctly (although I'm not sure the equations make much sense in this form).

Best wishes, Henry

from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
V = VectorH1(mesh, order=3, dirichlet='bottom|top')  # velocity
Q = H1(mesh,order=2, dirichlet='left|right') # pressure
X = V * Q 
(u, p), (v, q) = X.TnT()

gfu = GridFunction(X)
cf = CoefficientFunction([(2) if bc=="left" else (0.2) if bc=="right" else 0 for bc in mesh.GetBoundaries()])
gfu.components[1].Set(cf, definedon=mesh.Boundaries("left|right"))

f = LinearForm(X)
f.Assemble()

stokes = 0.1*InnerProduct(grad(u), grad(v))+div(u)*q+div(v)*p - 1e-10*p*q
a = BilinearForm(X, symmetric=True)
a += stokes*dx
a.Assemble()

inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")
res = f.vec.CreateVector()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res

print(gfu.components[1](mesh(1,0.5)))
print(gfu.components[1](mesh(0,0.5)))
 
More
1 year 2 months ago #4647 by Raphael_S
Hy Henry,  Many thanks for your support. As I see, I did not add the whole code. Sorry for this. Many thanks that you support me. Best wishes, Raphael

from netgen.occ import *
#from netgen.webgui import Draw as DrawGeo
from ngsolve import *
from netgen import*

# Velocity
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

box = Box((0,0,0), (3,1,1))
geo = box

geo.faces.Min(X).name='left'
geo.faces.Max(X).name='right'

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)).Curve(3)
#DrawGeo(geo)
Draw(mesh)

# Define a finite element space
V = VectorH1(mesh,order=3,dirichlet='front|back|bot|top') # velocity
# V = Periodic(VectorH1(mesh,order=3, dirichlet='bot|top'))  # velocity
Q = H1(mesh,order=2, dirichlet='left|right') # pressure

X = V * Q

# obtain test and trial functions
u, p = X.TrialFunction()
v, q = X.TestFunction()
# space for solution data
gfu = GridFunction(X)
# boundary condition
#uin = CoefficientFunction((U,0,0))
#gfu.components[0].Set(uin, definedon=mesh.Boundaries("bot"))
cf = CoefficientFunction([(2) if bc=="left" else (0.2) if bc=="right"  else 0 for bc in mesh.GetBoundaries()])
    # Set BC
gfu.components[1].Set(cf, definedon=mesh.Boundaries("left|right"))

nu = (eta/rho)  # kinematic viscosity

SetNumThreads(24) #Anzahl der Kerne zum Parallelisieren
with TaskManager():
    # right hand side of the equation
    f = LinearForm(X)
    f.Assemble()

    # left hand side of the equation
    stokes = nu*InnerProduct(grad(u), grad(v))+div(u)*q+div(v)*p - 1e-10*p*q
    a = BilinearForm(X, symmetric=True)
    a += stokes*dx
    a.Assemble()

    # calculate solution data
    inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")
    res = f.vec.CreateVector()
    res.data = f.vec - a.mat*gfu.vec
    gfu.vec.data += inv_stokes * res

    # calculate drag force


# density free pressure --> pressure
gfu.components[1].vec.data.FV().NumPy()[:] = gfu.components[1].vec.data.FV().NumPy() * rho

# pressure solution field
pressure = gfu.components[1]
# Define a coefficient function for val

# Draw velocity
#Draw(Norm(gfu.components[0]), mesh, "velocity",sd=3)   
#
#Draw((gfu.components[0]), mesh, "velocity",sd=3)          # Draw velocity
# Draw the pressure
Draw(pressure)                                         # Draw pressure

pInitLeft = gfu.components[1](mesh(0,0.5,0.5))
pInitRight = gfu.components[1](mesh(3,0.5,0.5))
## data visualisation / export

velocity = gfu.components[0]
nameOfFile = 'test'
    #nameOfFile = 'Initial_Conditions'
    # VTKOutput object
vtk = VTKOutput(ma=mesh,
                    coefs=[pressure,velocity],
                    names = ["pressure","velocity"],
                    filename=nameOfFile,
                    subdivision=0,
                    legacy = True)
    # Exporting the results:
vtk.Do()

 
More
1 year 2 months ago #4649 by hvwahl
Dear Raphael,

if you post code that isn't working as expected, please make a minimal working example, i.e., remove all code that isn't relevant to the issue and that runs on a laptop. Something like this:

from netgen.occ import *
from ngsolve 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

geo = Box((0,0,0), (3,1,1))

geo.faces.Min(X).name='left'
geo.faces.Max(X).name='right'
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.5))

V = VectorH1(mesh,order=2,dirichlet='front|back|bot|top') # velocity
Q = H1(mesh,order=1, dirichlet='left|right') # pressure
X = V * Q
u, p = X.TrialFunction()
v, q = X.TestFunction()

gfu = GridFunction(X)

cf = CoefficientFunction([(2) if bc=="left" else (0.2) if bc=="right" else 0 for bc in mesh.GetBoundaries()])
gfu.components[1].Set(cf, definedon=mesh.Boundaries("left|right"))

a = BilinearForm(X, symmetric=True)
a += ((eta/rho)*InnerProduct(grad(u), grad(v))+div(u)*q+div(v)*p - 1e-10*p*q) * dx
a.Assemble()

inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")
res = gfu.vec.CreateVector()
res.data = - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res
print(gfu.components[1](mesh(0, 0.5, 0.5)))
print(gfu.components[1](mesh(3, 0.5, 0.5)))

which again results in the expected output for me.

Best wishes,
Henry
The following user(s) said Thank You: Raphael_S
Time to create page: 0.114 seconds