Gradient calculation on the boundary with dirichlet bc

Main Forum Old-Forum - do not post here
1

Hi,

I am currently trying to calculate the drag / friction force in a fluid simulation.

In order to obtain the velocity and the pressure field, I have already solved the Stokes / Navier-Stokes equation, as described in the examples on the documentation website.
Now I need the gradient of the velocity normal to a boundary. Unfortunately, the gradient of the velocity is 0 on the boundaries with a dirichlet condition, making it impossible to determine the shear stresses on these surfaces.

I am using the Grad()-function like this:

velocity = gfu.components[0]

c_u = GridFunction(Q)
c_u.Set(velocity[0])
c_w = GridFunction(Q)
c_w.Set(velocity[2])

dudz = Grad(c_u)[2]
dwdx = Grad(c_w)[0]

tau = eta * (dudz + dwdx)

F_fric = Integrate(cf = tau,mesh=mesh,VOL_or_BND=BND,definedon=mesh.Boundaries(‘bottom’))

Is there any way to calculate the gradient of the velocity correctly on a surface with a dirichlet boundary condition?

I also tried to determine the tau tensor via the formula for incompressible fluids:

tau = eta (Grad(velocity)+Grad(velocity).trans)
F_fric = Integrate(cf = taunormal,mesh=mesh,VOL_or_BND=BND,definedon=mesh.Boundaries(‘bottom’))

But the gradient on “bottom” (dirichlet boundary condition) is always calculated to 0, this is also visible in the netgen visualisation.

Using grad() instead of Grad() doesn’t work, either.

I hope, I made my problem clear.
Thank you for your help!

Best regards,
Niklas

2

Hi Niklas,

in principle there shoudn’t be a problem with integrating the Gradient of a GridFunction on a boundary (Dirichlet or not). Could you please provide a MWE?

You might also want to look at this tutorial . Computing forces by testing the bilinear form with an appropriate non-conforming test-function computes forces significantly more accurately (with double order of convergence).

Best wishes,
Henry

3

Hi Henry,

thank you for your fast reply.

Unfortunately I am not able to attach files to this post. Therefore I am posting the source code as text:

[code]from ngsolve import *
from netgen.geom2d import *
from netgen.occ import *

U = 10 # velocity of surface
rho = 890 # density of the oil
eta = 0.0078792 # viscosity of the oil

filename_geometry = ‘wedge.step’

maximum element size of mesh

mesh_max_h = 1e-1

import geometry from given path

geo = OCCGeometry(filename_geometry)

generate and refine mesh

mesh = Mesh(geo.GenerateMesh(maxh=mesh_max_h))

generate finite elemente spaces

V = VectorH1(mesh,order=3, dirichlet=‘bottom|top|NONE’) # velocity
#V = VectorH1(mesh,order=3, dirichlet=‘inlet’) # velocity
Q = H1(mesh,order=2, dirichlet=‘inlet|outlet|right|left’) # 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(“bottom”))

draw the boundary conditions

Draw (Norm(gfu.components[0]), mesh, “velocity”, sd=3)

nu = (eta/rho)

SetNumThreads(8)
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
#stokes = (eta/rho)*InnerProduct(grad(u), grad(v)) + div(u)*q - div(v)*p + InnerProduct(v,div(u)*u) - 1e-12*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

density free pressure → pressure

gfu.components[1].vec.data.FV().NumPy()[:] = gfu.components[1].vec.data.FV().NumPy() * rho

velocity = gfu.components[0]

friction calculation

c_u = GridFunction(Q)
c_u.Set(velocity[0])
c_w = GridFunction(Q)
c_w.Set(velocity[2])

Draw(Norm(Grad(velocity)),mesh,‘velocity gradient’)

dudz = Grad(c_u)[2]
dwdx = Grad(c_w)[0]

tau = eta * (dudz + dwdx)

#Draw(Norm(Grad(velocity)),mesh,‘text’)
#Draw(tau,mesh,‘dudz’)

F_fric = Integrate(cf = tau,mesh=mesh,VOL_or_BND=BND,definedon=mesh.Boundaries(‘bottom’))
print(‘Friction force:’,F_fric)[/code]

Best regards,
Niklas

4

Hi Niklas,

integrating the gradient of an H1-GridFunction over the boundary will give zero (the gradient of H1 is in L2 -more precisely in HCurl- and therefore no -more precisely only the tangential- trace is defined).

Either you first interpolate the gradient of the GridFunction into another H1-GridFunction and then integrate this GridFunction over the boundary, or you use the before mentioned approach with the testfunction. I would recommend the second option, which is more accurate.

Best,
Michael