This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Incorrect wall shear stress computations compared to analytical solution

0 votes

Dear FEniCS comunity,

I want to compute the wall shear stresses resulting from CFD problems in FEniCS. For this I use:

# Define vectorspace for computing the wall shear stress
vector = VectorFunctionSpace(mesh, 'CG', velocity_degree - 1)
w = TestFunction(vector)
shear_stress = Function(vector)

# Define the stress tensor and vector
sigma = 2 * nu * sym(grad(uv_tmp))
T = sigma * n

# Define the wall shear stress
Tt = T - inner(T, n) * n
Lt = (1 / FacetArea(mesh)) * inner(w, Tt) * ds

assemble(Lt, tensor=shear_stress.vector())

Where shear_stress holds the final result. To my understanding the physics underlying this code should be correct for Newtonian fluids.

For benchmarking I compute the wall shear stress in a straight tube with a radius of 0.0105 [m]. The viscosity of the fluid is 8.4E-7 [Pa*s] (Because the density is scaled to 1 [kg/m3]). I simulate Poiseuille flow by prescribing a parabolic flow profile at the inlet:

vel_x_in = Constant(0.)
vel_y_in = Constant(0.)
vel_z_in = Expression(
    (pow(x[0], 2) + pow(x[1], 2) - ' +
    'pow(diameter / 2, 2)) * 2 * velocity / pow(diameter / 2, 2)',
    velocity=0.1, diameter=D, degree=2)

Here, velocity, is the mean velocity of the Poiseuille profile = 0.1 [m/s] = 0.5 * v_max.

Given the problem characteristics, the analytic value of the wall shear stress should be equal to: nu * du/dr = nu * (4 * v) / R = 8.4E-7 * 0.4/0.0105 = 3.2E-5 Pa (Note: the numerical value of the kinematic and dynamic viscosity are the same as the density = 1). However, my code returns a wall shear stress of around 6.7E-5 Pa.

When estimating the wall shear stress from the output files by hand (by estimating the term du/dr in paraview), I get wall shear stress values that are very close to the analytical solutions, indicating the simulations themselves are correct. Hence, I believe there is something wrong in my wall shear stress computations.

I cannot explain this discrepancy between analytical and numerical wall shear stress. Can anybody help me solve this problem?

Thank you in advance!
Sjeng

P.S. I do not know if it is relevant: I use fenics 2017.1 in combination with the oasis solver-library: https://github.com/mikaem/Oasis

asked May 24, 2017 by SQ FEniCS Novice (320 points)

2 Answers

0 votes
 
Best answer

Based on the answer I got on this question, I believe I am able to correctly compute the wall shear stress. I will put the rationale/ solution below for the sake of completeness.

The problem with the method as stated in the question is that the wall shear stress should be defined something like:

$$\int_{\Gamma}\tau_w\cdot v\,d\Gamma = \int_{\Gamma}(T - (T\cdot n) * n)\cdot v\,d\Gamma$$

where $$\tau_w$$ is the wall shear stress vector, v a test function and $$(T - (T\cdot n) * n)$$ the expression for the wall shear stress.

Trying to solve this with:

assemble(Lt = (1 / FacetArea(mesh)) * inner(w, Tt) * ds

(where Tt is the expression for the wall shear stress), only works for a zeroth-order discontinuous Galerkin method, as the mass matrix of this function space is a diagonal matrix where the non-zero entries are equal to the FacetArea. Hence, assembling this vector can be viewed as multiplying the expression for the wall shear stress (inner(w, Tt) * ds) by the inverse of the mass matrix (1 / FacetArea(mesh)). Clearly this only works when the mass matrix is a diagonal matrix. Hence, for all other function spaces this method is simply wrong. Moreover, using the DG0 function space is not the ideal in my opinion, as the nodes of this space are not located on the wall but in the center of the tetrahedrons, hence, your wall shear stress approximation will be slightly off. This effect might be neglected for very small boundary elements, but why would you accept this inaccuracy if you have all information available to get a better approximation (or maybe get the same inaccuracy with much larger elements).

The way I can now compute the wall shear stress correctly is by using the following code (Note: my TrialFunction is shear_stress, my TestFunction is w. Both are defined in the function space VectorFunctionSpace(mesh, 'CG', velocity_degree - 1))

# Define the stress tensor and vector
sigma = 2 * nu * sym(nabla_grad(uv_tmp))
T = sigma * n

# Define the wall shear stress
Tt = T - inner(T, n) * n
Lt = inner(w, Tt) * ds

# Define a left hand side for projection 
# CONSTANT throughout program if mesh does not change
a = inner(shear_stress, w) * ds
A_wss = assemble(a, keep_diagonal=True)
A_wss.ident_zeros()

# Assemble the left hand side
# Should assembled each time
b = assemble(Lt)

# Compute the wall shear stress
solve(A_wss, shear_stress_.vector(), b)

Note that in assembling the A matrix I specify keep_diagonal=True and A_wss.ident_zeros(). This is necessary because the expression for the wall shear stress is only defined on the boundary domain. Consequently, the problem is ill-posed: only some columns of the mass matrix have entries (boundary nodes), whereas most parts don't. Hence, for the empty columns we put a 1 on the diagonal with ident_zeros(), which makes the matrix invertible (For a better explanation see again the answer to my other question). This procedure shouldn't affect the solution of solve(A_wss, shear_stress_.vector(), b).

answered May 31, 2017 by SQ FEniCS Novice (320 points)
0 votes

If there's such a large discrepancy, then your code is probably not doing what you want it to. However, the general approach of directly evaluating a boundary flux at a Dirichlet boundary is not really best practice. A better way to extract fluxes is discussed here:

https://www.ices.utexas.edu/media/reports/2011/1122.pdf

and also in Exercise 8 of Section 2.12 of the T.J.R. Hughes FEM book. I tried throwing together a FEniCS implementation just now (with version 2016.2) for heat conduction, and it looks like it at least gets the sign right, but, use at your own risk...

from dolfin import *
from mshr import *
domain = Rectangle(Point(0.0,0.0), Point(1.0,1.0))
mesh = generate_mesh(domain,10)
V = FunctionSpace(mesh,"Lagrange",1)
class Bdry(SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary
bcs = [DirichletBC(V,Constant(0.0),Bdry()),]
u = TrialFunction(V)
v = TestFunction(V)
def R(u,v):
    return (inner(grad(u),grad(v)) - Constant(1000.0)*v)*dx
F = R(u,v)
a = lhs(F)
L = rhs(F)
u = Function(V)
solve(a==L,u,bcs)
class AntiBdry(SubDomain):
    def inside(self,x,on_boundary):
        return x[0] > DOLFIN_EPS and x[0] < 1.0-DOLFIN_EPS \
            and x[1] > DOLFIN_EPS and x[1] < 1.0-DOLFIN_EPS
uBdry = TrialFunction(V)
antiBCs = [DirichletBC(V,Constant(0.0),AntiBdry()),]
FBdry = (uBdry-1.0)*v*ds - R(u,v)
aBdry = lhs(FBdry)
LBdry = rhs(FBdry)
ABdry = assemble(aBdry,keep_diagonal=True)
bBdry = assemble(LBdry)
[bc.apply(ABdry,bBdry) for bc in antiBCs]
flux = Function(V)
solve(ABdry,flux.vector(),bBdry)
uFile = File("u.pvd")
uFile << u
fluxFile = File("flux.pvd")
fluxFile << flux
answered May 27, 2017 by kamensky FEniCS Novice (150 points)

Hi kamensky,

I looked at the article, but it seems that it is not really suited for my problem: I need to compute the stress vector resulting from the velocity field which is defined as:
$$\vec{\tau}=\sigma\cdot\vec{n}$$
where \sigma is the stress tensor. Thus I am not really computing a flux of any kind. I found another solution that seems to work for me (which I will post below).

Anyway, thank you for taking the time to answer my question.

Note that the traction is a flux for the Navier--Stokes problem. See Section 5.2 on fluid--structure interaction from the previously-linked report, or equations (8)--(13) from this report

https://www.ices.utexas.edu/media/reports/2004/0431.pdf

for a scalar problem with fluxes specified with a similar notation. (See also (31)--(33) from the above-linked report for flux extraction.) The heuristic to think of is: "If the Dirichlet BC were not specified strongly, what Neumann (flux) BC would I have needed to apply weakly on that part of the boundary to end up with the solution that I got?" An advantage of computing fluxes this way is that they satisfy conservation properties of the form (34). (In the case of a traction flux, this would be conservation of momentum.)

Indeed, you are right. I will give it a second look and see if this would be better for my aims. Thanks!

...