How-to implement a support of a simple supported beam?

I think you should try using a 0 DirichletBC for the displacement in all
directions on the lower left edge. Then you make a subdomain for the lower
right edge and impose a 0 DirichleBC for the displacement but only in the y.
The latter makes it possible for the beam to move in the x direction at that
point.

Johan

On Friday November 18 2011 08:40:57 Tim Berger wrote:
> New question #179257 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/179257
>
> I would like to compute with FEniCS a simple supported beam with two
> supports and an distributed load just like you can see on this pictue
> http://en.wikipedia.org/wiki/File:Bending.svg. But I'm not sure how to
> implement the two supports correctly, especially the one which allow only
> rotation (the right support on the picture).
>
> So far I have managed to implement an example with a cantilever (a beam
> with a full moment connection (like a horizontal flag pole bolted to the
> side of a building)). In that example I have defined the complete right
> boundary as a sub domain and used that sub domain for the Dirichlet
> boundary conditions.
>
> class right_boundary(SubDomain):
>
> def inside(self, x, on_boundary):
>
> return on_boundary and (abs(x[0]-xlength) < tol)
>
> right = right_boundary()
>
> bc_right_0 = DirichletBC(V.sub(0), Constant(0.0), right)
>
> bc_right_1 = DirichletBC(V.sub(1), Constant(0.0), right)
>
> As I understand it these conditions do not allow the rotation of the right
> boundary around the center point. But how can I achieve that?
>
> Should I define only a small section of the right boundary as a sub domain
> and use that with the DirichletBC function?

Thank you very much Johan!

I'm going to try this. (So far I have used sections at center of the left and right boundary as a sub domain - and not the lower edges)

Tim

Are you modelling with an Euler (thin) type beam with one unknown, the
transverse displacement, or a Timoshenko (thick) type beam, with
transverse displacement and a rotation? Are you trying to model the
axial extension as well?

I'm trying to get an idea of which variable in your function space V
is which, it's not totally clear at the moment.

bc_right_0 = DirichletBC(V.sub(0), Constant(0.0), right)
bc_right_1 = DirichletBC(V.sub(1), Constant(0.0), right)

This will constrain both degrees of freedom. To constrain only the
transverse displacement just include one of the above statements (the
right one!).

Jack Hale

On 18 November 2011 16:40, Tim Berger
<email address hidden> wrote:
> New question #179257 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/179257
>
> I would like to compute with FEniCS a simple supported beam with two supports and an distributed load just like you can see on this pictue http://en.wikipedia.org/wiki/File:Bending.svg.
> But I'm not sure how to implement the two supports correctly, especially the one which allow only rotation (the right support on the picture).
>
> So far I have managed to implement an example with a cantilever (a beam with a full moment connection (like a horizontal flag pole bolted to the side of a building)).
> In that example I have defined the complete right boundary as a sub domain and used that sub domain for the Dirichlet boundary conditions.
>
>    class right_boundary(SubDomain):
>
>        def inside(self, x, on_boundary):
>
>            return on_boundary and (abs(x[0]-xlength) < tol)
>
>    right = right_boundary()
>
>    bc_right_0 = DirichletBC(V.sub(0), Constant(0.0), right)
>
>    bc_right_1 = DirichletBC(V.sub(1), Constant(0.0), right)
>
> As I understand it these conditions do not allow the rotation of the right boundary around the center point.
> But how can I achieve that?
>
> Should I define only a small section of the right boundary as a sub domain and use that with the DirichletBC function?
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

I'm trying to implement a Euler type beam (http://en.wikipedia.org/wiki/Euler%E2%80%93Bernoulli_beam_equation) and I would like know the transverse displacement (at the half of the length of the beam).

(In example the supports are exchanged compared with the situation on the picture from the Wikipedia.)

Tim,

As Euler theory is fourth order PDE it will include derivatives up to
second order in the weak form. This will require a subspace of H^2 to
be constructed, with C^1 continuous type finite elements. This is one
of the main reasons that Timoshenko Beam theory is used in FE
implementations. I'm pretty sure DOLFIN doesn't have any C^1 elements
(at the moment) implemented? Which elements are you using?

Jack Hale

On 18 November 2011 18:10, Tim Berger
<email address hidden> wrote:
> Question #179257 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/179257
>
> Tim Berger posted a new comment:
> (In example the supports are exchanged compared with the situation on
> the picture from the Wikipedia.)
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

I have tried to implement the left and right lower edge as subdomains.
for the lower left corner:
llc = compile_subdomains('(std::abs(x[0]) < 1e-03) && (std::abs(x[1]) < 1e-03) && on_boundary')

and for the lower right corner:
lrc = compile_subdomains('(std::abs(x[0]-length) < 1e-03) && (std::abs(x[1]) < 1e-03) && on_boundary')

lrc.length = xlength

while the mesh is defined as:
xlength = 100.0 #[mm]

ylength = 2.0 #[mm]

xelements = 500

yelements = 10

mesh = Rectangle(0, 0, xlength, ylength, xelements, yelements)

On the upper boundary is a evenly distributed load.

top = compile_subdomains('(std::abs(x[1]-length) < 1e-03) && on_boundary')

top.length = ylength

V = VectorFunctionSpace(mesh,'CG',1)

D = mesh.topology().dim()

s_domains = MeshFunction("uint", mesh, D-1)

s_domains.set_all(0)

llc.mark(s_domains,1)

lrc.mark(s_domains,2)

top.mark(s_domains,3)

But when do this I get always the error message:
"Warning: Found no facets matching domain for boundary condition."

It seems that I have to catch at least two vertices with the subdomain definition.
This works:
llc = compile_subdomains('(std::abs(x[0]) < 1e-03) && (std::abs(x[1]) < 0.3) && on_boundary')
lrc = compile_subdomains('(std::abs(x[0]-length) < 1e-03) && (std::abs(x[1]) < 0.3) && on_boundary')
but I would say that this (std::abs(x[1]) < 0.3) catches two vertices in x[1] direction. One element length is here 0.2.

Is there a way how I can grab only the last vertex in each corner?

Try this code:

from dolfin import *

xlength = 100.0 #[mm]
ylength = 2.0 #[mm]
xelements = 500
yelements = 10

mesh = Rectangle(0, 0, xlength, ylength, xelements, yelements)

# Note use of near and batch compilation of seberal subdomains
llc, lrc, top = compile_subdomains(['near(x[0], 0.0) && near(x[1], 0.0)', \
                                    'near(x[0], length) && near(x[1], 0.0)',
                                    'near(x[1], length) && on_boundary'])

lrc.length = xlength
top.length = ylength

# Note use of FacetFunction with default value of 0
facet_domains = FacetFunction("uint", mesh, 0)
top.mark(facet_domains, 3)

V = VectorFunctionSpace(mesh,'CG',1)

# Note use of Vector DirichletBC and "pointwise"
bc_right = DirichletBC(V, Constant((0.0, 0.0)), lrc, "pointwise")
bc_left_0 = DirichletBC(V.sub(1), Constant(0.0), llc, "pointwise")

print "Marked edges:", sum(facet_domains.array()==3)
print "Fixed dofs:", bc_right.get_boundary_values()
print "Fixed dofs:", bc_left_0.get_boundary_values()

Johan

On Friday November 18 2011 11:25:51 Tim Berger wrote:
> Question #179257 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/179257
>
> Status: Answered => Open
>
> Tim Berger is still having a problem:
> I have tried to implement the left and right lower edge as subdomains.
> for the lower left corner:
> llc = compile_subdomains('(std::abs(x[0]) < 1e-03) && (std::abs(x[1]) <
> 1e-03) && on_boundary')
>
> and for the lower right corner:
> lrc = compile_subdomains('(std::abs(x[0]-length) < 1e-03) &&
> (std::abs(x[1]) < 1e-03) && on_boundary')
>
> lrc.length = xlength
>
> while the mesh is defined as:
> xlength = 100.0 #[mm]
>
> ylength = 2.0 #[mm]
>
> xelements = 500
>
> yelements = 10
>
> mesh = Rectangle(0, 0, xlength, ylength, xelements, yelements)
>
> On the upper boundary is a evenly distributed load.
>
> top = compile_subdomains('(std::abs(x[1]-length) < 1e-03) &&
> on_boundary')
>
> top.length = ylength
>
>
> V = VectorFunctionSpace(mesh,'CG',1)
>
> D = mesh.topology().dim()
>
> s_domains = MeshFunction("uint", mesh, D-1)
>
> s_domains.set_all(0)
>
> llc.mark(s_domains,1)
>
> lrc.mark(s_domains,2)
>
> top.mark(s_domains,3)
>
> But when do this I get always the error message:
> "Warning: Found no facets matching domain for boundary condition."
>
> It seems that I have to catch at least two vertices with the subdomain
> definition. This works:
> llc = compile_subdomains('(std::abs(x[0]) < 1e-03) && (std::abs(x[1]) <
> 0.3) && on_boundary') lrc = compile_subdomains('(std::abs(x[0]-length) <
> 1e-03) && (std::abs(x[1]) < 0.3) && on_boundary') but I would say that
> this (std::abs(x[1]) < 0.3) catches two vertices in x[1] direction. One
> element length is here 0.2.
>
> Is there a way how I can grab only the last vertex in each corner?

@Jake:
I'm using V = VectorFunctionSpace(mesh,'CG',1). Continuous Galerkin or Langrange elements with the degree 1 (triangle with nodes at the three vertices).

My whole script is now this:
from dolfin import *

xlength = 100.0 #[mm]
ylength = 2.0 #[mm]
xelements = 500
yelements = 10
mesh = Rectangle(0, 0, xlength, ylength, xelements, yelements)

llc, lrc, top = compile_subdomains(['near(x[0], 0.0) && near(x[1], 0.0)', \
                                    'near(x[0], length) && near(x[1], 0.0)',
                                    'near(x[1], length) && on_boundary'])

lrc.length = xlength
top.length = ylength

facet_domains = FacetFunction("uint", mesh, 0)
top.mark(facet_domains, 3)

tr=Expression((('X1','X2')), X1=0.0, X2=0.0) # [N/mm]
tr.X1=0.0 # [MPa]
tr.X2=-0.2 # [MPa]

V = VectorFunctionSpace(mesh,'CG',1)

# Note use of Vector DirichletBC and "pointwise"
bc_right = DirichletBC(V, Constant((0.0, 0.0)), lrc, "pointwise")
bc_left_1 = DirichletBC(V.sub(1), Constant(0.0), llc, "pointwise")
bc = [bc_left_1, bc_right]

#definition for the variational formulation
u = TrialFunction(V)
w = TestFunction(V)

#material coeff.s of a stainless steal
nu = 0.3 # POISSON ratio
E = 210000.0 #MPa Young modulus
G = E / (2.0 * (1.0 + nu)) # shear modulus
l = 2.0 * G * nu / (1.0-2.0 * nu) # lambda (a Lame constant)

#building the variational form F=0
F = w[i].dx(i) * 2.0 * G * nu / (1.0 - 2.0 * nu) * u[k].dx(k) * dx
F += w[i].dx(j) * G * u[i].dx(j) * dx
F += w[i].dx(j) * G * u[j].dx(i) * dx
F += -w[i] * tr[i] * ds(3)
a = lhs(F)
L = rhs(F)
A = assemble(a, exterior_facet_domains=facet_domains, mesh=mesh)
b = assemble(L, exterior_facet_domains=facet_domains, mesh=mesh)

for condition in bc :
    condition.apply(A, b)

u = Function(V)
solve(A, u.vector(), b)

# write out
file_u = File('elastostatic_deformations.pvd')
file_u<< u

# Plot
plot(u, mode='displacement', interactive=True)

@Johan:
Thank very much you for your help!! Your answers solved my questions.

One little correction: The dimension of tr.X1 and tr.X2 should be N/mm (and not MPa or N/mm^2).

@Jack

Tim is not using any beam theory, he is solving a linear momentum balance in two-dimensions (Navier-Lame equations, div Sigma=0 with Hooke).

On Sat, Nov 19, 2011 at 07:50:53PM -0000, B. Emek Abali wrote:
> Question #179257 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/179257
>
> B. Emek Abali posted a new comment:
> @Jack
>
> Tim is not using any beam theory, he is solving a linear momentum
> balance in two-dimensions (Navier-Lame equations, div Sigma=0 with
> Hooke).

I would suggest using the cbc.twist module from CBC.Solve. It gives
you a range of common models for hyperelasticity and allows very
simple specification of boundary conditions:

  https://launchpad.net/cbc.solve

--
Anders