Where’s the force sitting?

1

Hello,
It’s me again. I’d like to simulate a cantilever under tension or compression. I modify the file called: demo_hyperelasticity.py from
https://fenicsproject.org/docs/dolfin/latest/python/demos/hyperelasticity/demo_hyperelasticity.py.html`

And I have some questions before and after the simulation.
My first question:

  1. Does that demo suit for my simulation?
    Assume the answer is “Yes”, and then here the modifications:
# a. Mesh
#=======
...
X = 20
Y = 5
Z = 5
canti = Box( Point(0,0,0), Point(X,Y,Z) )
mesh = generate_mesh (canti, 40)

# b. Boundary
# ===========
# ...
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)

c = Constant((0, 0, 0))

bcl = DirichletBC(V, c, left)
#bcr = DirichletBC(V, r, right)
bcs = [bcl] #,bcr]

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], X)

right = Right()

# exterior facets MeshFunction
facets = MeshFunction("size_t", mesh, 1)
facets.set_all(0)
right.mark(facets, 1)
ds = Measure('ds', subdomain_data=facets)

# c. Define functions
# ===================
# ...
B  = Constant((0.0, 0.0, 0.0))  # Body force per unit volume
T  = Constant((1.0, 0.0, 0.0))  # Traction force on the boundary

# d. Total potential energy
# =========================
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds(1)

# e. Solve variational problem
# ============================
solve(F == 0, u, bcs, J=J,
      form_compiler_parameters=ffc_options)

# f. Save solution in VTK format
# ==============================
file = File("displacement.pvd");
file << u;

With this arrangement, I believe I have a cantilever, no movement on the left and it’s pulled to the right, and the force sitting on all of facets “1”.

My second question:
2. Does this arrangement suit my idea?

Assume the answer is “Yes”, then I send the result to Paraview, and I find the figure is not what I expect. I’m pulling the cantilever to the right, and the force sitting on all of the nodes of the right plane (I marked “R”), but the displacement of X on the same plane (“R”) has different sign (it doesn’t make sense), and seems the force is sitting on the corner only.

dispX
dispX882×490 27.9 KB

  1. Why?
    That’s my final question.

Thanks for the enlightenment.

Suess

2

Please encapsulate your code with ```, to make it readable. Please also supply a minimal working example and not just parts of your code. Remove all code that you are not using, and make sure that what you post can be run by anyone else.

1 Like
3

Note that this is the source of your error, as you are using a 3 dimensional mesh, the dimension should be 2, and not 1, to define facets. By changing this to 2, you obtain facet and not edge integrals of the fource, and the expected result.

1 Like
4

I see, thank you Dokken. I get the correct result now.

dispX-OK
dispX-OK882×490 16.2 KB