### Wall Shear Stress calculation error

306
views
0
6 months ago by
Hi,

I'm trying to compute wall shear stress after applying Navier-Stokes to a 3D cylinder. I've followed the suggestion in https://fenicsproject.org/qa/13629/incorrect-stress-computations-compared-analytical-solution but I'm getting the error:

Traceback (most recent call last):
File "wss_error.py", line 137, in <module>
A_wss.ident_zeros()
AttributeError: 'Vector' object has no attribute 'ident_zeros

Fenics version - 1.6.0

Can anyone explain why this is happening?

Many thanks!

from dolfin import *

# Define function spaces
V = VectorFunctionSpace(mesh, "CG", 1)
Q = FunctionSpace(mesh, "CG", 1)

# Test and trial functions
v = TestFunction(V)
q = TestFunction(Q)
u = TrialFunction(V)
p = TrialFunction(Q)

# Set parameter values
dt = 0.1                              # time step
mu = 0.00345                           # dynamic viscosity (Ns/m2)
rho = (1.056*10**(-3))                 # blood density (g/l)
nu = mu/rho                            # kinematic viscosity
u0 = Constant((0., 0., 0.))            # initial velocity
p0 = Constant(0.)                      # initial pressure
r = 1.0                                 # radius
D = 2.0                                  # Diameter
l = 10.8                                  # Lenght of the tube
n = FacetNormal(mesh)

# Functions
u0 = interpolate(u0, V)
u1 = Function(V)
p0 = interpolate(p0, Q)
p1 = interpolate(p0, Q)
k  = Constant(dt)
n  = FacetNormal(mesh)
w= TestFunction(V)
shear_stress= Function(V)
Tt= Function(V)

parabolicprof ='600*(1-((pow(x[0],2))+(pow(x[1],2)))/pow('+str(r)+',2))'
finlet = Expression ((('0'),('0'),(str(parabolicprof))))

# Define boundary conditions
noslip  = DirichletBC(V, Constant ((0,0,0)),boundaries,3)
inflow  = DirichletBC(V, finlet, boundaries,2)
bcu = [noslip, inflow]
f = Constant((0, 0, 0))   # force
def epsilon(v):return 0.5*(grad(v)+transpose(grad(v)))
a11 = (1/k)*inner(v,u)*dx
a12 = 2.0*inner(epsilon(v), nu*epsilon(u))*dx-inner(v, nu*grad(u).T*n)*ds
a13 = inner(v, grad(u)*u0)*dx
L1 = (1/k)*inner(v,u0)*dx+1/rho*inner(epsilon(v), p0*Identity(u.cell().d))*dx - 1/rho*inner(v, p0*n)*ds

# Pressure correction
a2 = inner(grad(q), grad(p))*dx
L2 = inner(grad(q), grad(p0))*dx - rho*(1/k)*q*div(u1)*dx

# Velocity correction
a3 = inner(v, u)*dx
L3 = inner(v, u1)*dx - k/rho*inner(v, grad(p1 - p0))*dx

#Assemble matrices
A11 = assemble(a11)
A12 = assemble(a12)
A13 = assemble(a13)
A2 = assemble(a2)
A3 = assemble(a3)
A1 = A11+A12+A13
b = assemble(L1)
parameters['allow_extrapolation'] = True

# Compute tentative velocity step
begin("Computing tentative velocity")
b1 = assemble(L1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1, "gmres", "ilu")
end()

# Velocity correction
begin("Computing velocity correction")
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "gmres", "ilu")
end()

# Define the stress tensor and vector
sigma = 2 * mu * sym(nabla_grad(u1))
Tensor = sigma * n

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

# Define a left hand side for projection
a = inner(shear_stress, w) * ds
A_wss = assemble(a, keep_diagonal=True)
A_wss.ident_zeros()

# Assemble the left hand side
b = assemble(Lt)

# Compute the wall shear stress
solve(A_wss, shear_stress_.vector(), b)​
Community: FEniCS Project
1
shear_stressis a Function (not a TrialFunction) and then A_wss is a Vector and a Vector object doesn't have the ident_zeros() method
written 6 months ago by Hernán Mella
Hi,

Probably I didn't explain myself in the best way because I've started to work with fenics quite recently. I did understand that Vector doesn't have ident_zeros() method, despite it was mentioned has working in the link I've posted.

So, the question is, if wall shear stress is a Vector, how can I compute the solution?
written 6 months ago by IVelho
I have answered your question below
written 6 months ago by Hernán Mella

### 1 Answer

0
6 months ago by
Try something like this:

# Define a left hand side for projection
u = TrialFunction(V)
a = inner(u, w) * ds
A_wss = assemble(a, keep_diagonal=True)
A_wss.ident_zeros()

# Assemble the left hand side
b = assemble(Lt)

# Compute the wall shear stress
solve(A_wss, shear_stress.vector(), b)​​
plot(shear_stress, interactive=True)
Please login to add an answer/comment or follow this question.

Similar posts:
Search »