### Wall Shear Stress calculation error

89

views

0

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:

Fenics version - 1.6.0

Can anyone explain why this is happening?

Many thanks!

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 Answer

0

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.

`shear_stress`

is a Function (not a TrialFunction) and then`A_wss`

is a Vector and a Vector object doesn't have the`ident_zeros()`

methodProbably 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?