annular and longitudinal stresses of cylinder

8 weeks ago by
Hello, I am novice at FEniCS, I want to know how could I compute the annular and longitudinal stresses of cylinder? In the-equations-of-linear-elasticity there is a Von_Mises stress:

from fenics import *

# Scaled variables
L = 1; W = 0.2
mu = 1
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
g = gamma

# Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 10, 3, 3)
V = VectorFunctionSpace(mesh, 'P', 1)

# Define boundary condition
tol = 1E-14

def clamped_boundary(x, on_boundary):
    return on_boundary and x[0] < tol

bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    #return sym(nabla_grad(u))

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution
plot(u, title='Displacement', mode='displacement')

# Plot stress
s = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  # deviatoric stress
von_Mises = sqrt(3./2*inner(s, s))
V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)
plot(von_Mises, title='Stress intensity')

# Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
Community: FEniCS Project

1 Answer

7 weeks ago by
Alright, you have the stress defined as a python function for a displacement in a cartesian basis.

First thing we need are the conversion formulas. If your longitudinal axis is called \(z\) and coincides with your third cartesian axis then also the respective stress components coincide.
For the angular stress (hoop stress) we have:

For both these scalar functions we need a suitable function space to project on:
W = FunctionSpace(mesh, 'DG', 0)​

The longitudinal component is straight-forward:

sigma_zz = project(sigma(u)[2,2], W)

For the hoop stress we first need to define the angle as a function of the cartesian coordinates:

x = SpatialCoordinate(mesh)
def theta(x):
    return tan(x[1]/x[0])

With these formulas at hand we can perform the projection:

sigma_tt = project(sin(theta(x))**2*sigma(u)[0,0] - 
                   cos(theta(x))**2*sigma(u)[1,1] , W)

Hope this helps.
Please login to add an answer/comment or follow this question.

Similar posts:
Search »