Buckling of an hyperelastic beam

dolfinx
1

Hi, I’m trying to solve a buckling problem for an hyperelastic beam.
Here I show my code:

from dolfin import *
import fenics as fe
import time


mesh = RectangleMesh(Point(0.0, 0.0), Point(1.0, 0.1), 200, 50, "crossed")

V = VectorFunctionSpace(mesh, "Lagrange", 2)

left = CompiledSubDomain("near(x[0], side) && on_boundary", side=0.0, tol=10e-15) 
right = CompiledSubDomain("near(x[0], side) && on_boundary", side=1.0, tol=10e-15)

c = Expression(("0.0", "0.0"), degree=2)
d = Expression(("-0.2", "0.0"), degree=2)
t = Expression(("0.0", "0.0"), degree=2)


du = TrialFunction(V)            
v = TestFunction(V)            
u = Function(V)                 

D = mesh.topology().dim()
print(D)
neumann_domain = MeshFunction("size_t", mesh, D-1)
neumann_domain.set_all(0)
CompiledSubDomain("near(x[0], side) && on_boundary", side=1.0, tol=10e-15).mark(neumann_domain, 1)
ds = Measure("ds", subdomain_data=neumann_domain)


bc1 = DirichletBC(V, c, left)
bc2 = DirichletBC(V, d, right)
bcs=[bc1,bc2]

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             
F = I + grad(u)             
C = F.T*F                  

# Invariants of deformation tensors
Ic = tr(C)
J = det(F)

E, nu = 1e6, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

#Neo-Hookean
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

bb = Expression(("0.0", "-1000.0"), degree=2)
# Total potential energy
Pi = psi*dx - dot(t, u)*ds(1) #+ dot(bb,u)*dx

F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-9
prm['newton_solver']['relative_tolerance'] = 1E-9
prm['newton_solver']['maximum_iterations'] = 200

solver.solve()

L2 = inner(u, u) * dx
H10 = inner(grad(u), grad(u)) * dx
energynorm = sqrt(assemble(psi*dx))
# H1 = inner(F, P) * dx
L2norm = sqrt(assemble(L2))
H10norm = sqrt(assemble(H10))

print("L2 norm = %.10f" % norm(u, norm_type="L2"))
print("H1 norm = %.10f" % norm(u, norm_type="H1"))
print("H10 norm = %.10f" % norm(u, norm_type="H10"))


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

#NEO-HOOKEAN

With an imposed displacement of -0.2 on the free right end, I should obtain a solution with buckling, but I only get a solution with 0 vertical displacement.
How can I solve this issue?
Thank you

2

To me, this reads that you are imposing a zero force on boundary, hence the zero solution.

3

Thank you for your prompt reply. I have removed it, however the solution I get is again the one without buckling.

4

What do you mean by removed it? I guess it must be there, just t should be non zero (e.g. with first component equal to -0.2).

5

I think it is not required to have t different from 0. We work with controlled displacement.
bucklbeam
I get this result that is acceptable, however it is not unique. I don’t know how to get the buckling solution.

6

@HumbertHumbert7
I think you can implement a DirichletBC at the right end specifying the displacement you want. The buckling also depends on the length to width ratio. If width>>length, you will see thickening of the beam instead of buckling. If width<<length, you will see buckling. The optimal ratio of length to width in order to get buckling depends on the problem.

7

Thank you, in my code I have used this function
d = Expression((“-0.2”, “0.0”), degree=2)
and the boundary condtion on the right side is bc2 = DirichletBC(V, d, right).

8

@HumbertHumbert7
Yes, that’s right and that is the reason the length has shortened from 1 to 0.8. What I don’t understand is the role of dot(t,u)*ds(1). If t=Expression(("0.0", "0.0"), degree=2), then dot(t,u)*ds(1) must be 0.
You can get buckling if you make the height small and the length large.

9

It is present because I started with another case of buckling in presence of a compression force on that boundary. Then, I passed to a case of controlled displacement, so I put t equal to 0. I couldn’t get buckling in either case.
Regarding the geometry and the displacement, they should produce a buckling solution not only the one with shortening of the beam, but I am not able to understand how can I get that solution.

10

@HumbertHumbert7
You can try reducing the value of E without changing the geometry.

11

I have tried E=1e3 but it’s the same.

12

It’s possible that you need to tweak the initial guess to the nonlinear solver. Seems to me that you are starting from a zero initial guess, and this may favor convergence to the solution in your third post. Try and start from a non-zero initial guess, possibly not symmetric in y.

My understanding is that for this kind of buckling/bifurcation problems, people typically use continuation methods (or similar techniques) to build up a sequence of initial guesses. I don’t have much experience in this field, but maybe this gives you a keyword to look up.

13

Thank you very much. I think this could solve the problem, do you know how to change the initial guess?

14

I found an example where the continuation method is implemeted for an hyperelastic beam: pefarrell / defcon — Bitbucket

I have installed firedrake from Fem on Colab and defcon using !pip install.
Now I don’t know how to run the examples, could you help me?

15

Simply interpolate or project an expression into u, it will get automatically used as initial guess.

16

There is an example folder in the repository. Pick one you are interested in, and copy its code in the notebook. If the example is split into a couple of files, just copy the content of all of them into the same notebook (dropping any import of files you’ve copied already: for instance, drop Bitbucket once you have copied the content of Bitbucket ).
It’s not very clean or very smart, but assuming that you’ll run the example once in your life to see how it goes and move on afterwards, it’s the simplest solution.

17

Thank you very much.
The problem is that if I just copy and paste this: pefarrell / defcon / examples / firedrake / hyperelasticity / hyperelasticity.py — Bitbucket

Then I get this error:

NameError Traceback (most recent call last)
in <cell line: 10>()
8 from slepc4py import SLEPc
9
—> 10 class HyperelasticityProblem(BifurcationProblem):
11 def mesh(self, comm):
12 return RectangleMesh(40, 40, 1, 0.1, comm=comm)

NameError: name ‘BifurcationProblem’ is not defined

18

BifurcationProblem is part of defcon itself (see Bitbucket), so you should be getting in when doing from defcon import * at Bitbucket

Considering that I am the developer of FEM on Colab, I’d be happy to take a look at the notebook to verify that indeed defcon is compatible with FEM on Colab. However, doing this through a series of replies here is inefficient. Please share your notebook with my email address (see website linked in my profile) through Google Colab.

Apart from debugging this minor technical issue related to Colab, keep posting in this topic for further questions, since (1) for the community, it’s important that the discussion remains public, (2) for you, in this way you will get several people answering, and not only me.

19

Thank you. I shared the code with you.
I have already copied that line, but the error remains.

20

The main issue was that you were installing defcon from pypi, i.e. defcon · PyPI , which is another library with the same name.
Installation must be done cloning the (right) defcon from bitbucket.

Unfortunately, running the demo fails with

Defcon started with only 1 process.
At least 2 processes are required (one master, one worker).

which means that you won’t be able to use it on Colab, where you only have 1 process. If you need further able specific to defcon the best course of action would be to get in touch with its authors, because I (and most of the people in this community) don’t know much about it, and will probably not be able to help much further.