Check boundary condition of 2D-linear elasticity problem

General
,
1

Hi all!

I just implemented the exact “Orthotropic linear elasticity” example from Orthotropic linear elasticity — Numerical tours of continuum mechanics using FEniCS master documentation. After solving the equation, I want to check whether my traction boundary condition on the Top boundary.

With the following extension (where T0 and T1 are the components of the traction force T which is 0, 1e-03 per default in this example):

Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
sig = Function(Vsig, name="Stress")
sig.assign(project(sigma(u), Vsig))
import numpy as np
for right_boundary_point in [Point(x[0], x[1]) for x in mesh.coordinates() if np.isclose(x[0], L)]:
    sigma_tensor = np.array(sig(right_boundary_point)).reshape(2, 2)
    print(f"sigma_tensor[:, 1]: {sigma_tensor[:, 1]}")
    print(f"T: [{T0}, {T1}]")

I checked the values of the stress tensor on the boundary, multiplied by (0, 1), which is the outer normal. The results are not the exact values of T, but something like [-3.83081901e-07, 1.00894364e-03] for example.

Now in my application I have a lot larger values and T is around [0, 1e6]. With this value I get the following results (just examplewise for one point on the boundary): [-8.29929781e+02, 1.02096271e+06], which in the first component doesn’t match T at all.

Can someone tell me the reason? Is it just a numerical rounding error (maybe related to the mesh or something) or am I understanding the boundary condition wrong?

Best regards and thank you in advance!

2

Are you using a coarse mesh?

In FEM, we don’t identically satisfy PDEs nor boundary conditions (in most cases). Rather, they are satisfied in a weighted sense, tested against all test functions. So a missmatch is to be expected. You would expect this missmatch to converge to zero with mesh refinement.

Additionally, extracting boundary stresses from the FEM output is actually a bit of a sensitive issue. If you’re interested in high-fidelity results, or need to couple the stresses to a different simulation, that is something to be aware of. Keyword to look for here is ‘flux extraction’.

3

You’re right, I increased the mesh resolution from N = 50 (image top row) to N = 200 (image bottom row), and the error decreased. (I plotted sigma(u)*n(x) for the boundary points and I used Vsig = TensorFunctionSpace(mesh, “P”, degree=1) now)

image
image1580×1656 120 KB

It’s good to know, that the boundary condition is just approximated and doesn’t have to be satisfied exactly (as I assumed beforehand).

4

@Elias_Reutelsterz Did you have a look at this in the FEniCS tours, by any chance ? Computing consistent reaction forces — Numerical tours of continuum mechanics using FEniCS master documentation

Reading your interrogations, I think it might hold the answers you are looking for.