How do I project stresses from displacements?

General
1

Hello,
I want to solve a linear elastic problem and I have problems to calculate the resulting stresses.

As far as I understand it, first the variational problem is solved to obtain the displacements. Then these displacements are used to calculate the stresses afterwards.

However, the calculation of the displacement is based on the function for the stresses, so one can assume that if the displacements look correct, the stresses must also be correct. Unfortunately, this is not the case, so I suspect that the error is in the subsequent calculation of the stresses.

To calculate the stresses, ‘project’ is used. Here I have already tried to change the degree (1,2,3) as well as the type of functions in the function space (P, CG, DG), but unfortunately this does not provide the correct stresses.

I would be very happy about a helpful comment.

Best regards

Felix

This is my code:


from fenics import *
import matplotlib.pyplot as plt

# Geometry and mesh
mesh = RectangleMesh(Point(0, 0), Point(1.0, 1.0), 10, 10)

# Define function spaces
V = VectorFunctionSpace(mesh, 'P', 1)
T = TensorFunctionSpace(mesh, 'P', 1)  # Space for stresses
u = TrialFunction(V)
v = TestFunction(V)

# Define material parameters
mu = Constant(1.666666667)
lambda_ = Constant(0.576923077)

# Helper function for deformation gradient
def epsilon(u):
    return 0.5 * (grad(u) + grad(u).T)

# Function to compute the stresses
def stress(u):
    return lambda_ * div(u) * Identity(2) + 2 * mu * epsilon(u)

# Define boundary conditions as functions
def left_boundary_condition():
    return DirichletBC(V, Constant((0, 0)), "near(x[0], 0)")

def right_boundary_condition(disp_value):
    return DirichletBC(V, Constant((disp_value, 0)), "near(x[0], 1.0)")

# Define weak form of the problem
F = inner(stress(u), epsilon(v)) * dx
a, L = lhs(F), rhs(F)

u_solution = Function(V)  # Define the solution function

# Apply boundary conditions
bcs = [left_boundary_condition(), right_boundary_condition(0.1)]

# Solve the problem
solve(a == L, u_solution, bcs)

s = stress(u_solution)
sigma = project(s, T)

# Visualize the stresses
c = plot(sigma[0, 0])  # XX-component of the stress
plt.colorbar(c, orientation='vertical', label='Stress XX')
plt.title("Distribution of Stress XX-component")
plot(mesh, color='k', linewidth=0.5)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

# Visualize the displacements
c = plot(u_solution, mode='displacement')
plt.title("Final Displacement Distribution")
plot(mesh, color='k', linewidth=0.5)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()
2

Please format your code using the three backtick syntax (also available in the GUI editor as the </> symbol) like so:

``` your code here ```

The way your code is currently formatted is time intensive to try to reformat and run

What are you expecting to see in the stress distribution? Can you write out in a bit more detail what you mean by:

… unfortunately this does not provide the correct stresses.

Additionally, one should note that in finite element analysis, convergence rates are typically different for displacements and stresses (with the former being higher order convergence). I would also try comparing what you expect to see with a more resolved mesh as a coarse mesh might not show you much.

Another thing you might also note is the difference in the ufl function div() vs nabla_div() and nabla_grad() vs grad(). dokken’s excellent FEniCSx tutorial covers this difference here

A final footnote: If you plan are working with the new version of FEniCS, FEniCSx, you will notice that the “project” function no longer exits. See this link for more information on the reason why that was done. At this point, I would highly suggest all new work you do to be in FEniCSx/Dolfinx. This link will be the closest to the legacy fenics code you have been working with.

1 Like
3

I expect the stresses to be high on the left side and decrease towards the right. But the way the code works now, the stresses are in the corners.

I have also tried other mesh resolutions, but that does not make the stresses concentrate in the corners.

I have not been working with fenics for so long and will have a closer look at your links. Thank you very much.