Simply supported beam with uniform load


65
views
0
19 days ago by
Aby  
Hello All,
I am trying to create a simply supported beam as shown in the following figure:


I want to apply load to the top facet of the beam. For this I have so far written the following code but it is giving JIT failure. I am unable to find my mistake and thus will be very grateful if I could please get some help.

Thank you very much.

Code so far:
from __future__ import print_function
#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
Created on Thu Jul 26 13:49:52 2018
@author: Aby
"""
from fenics import *

#select top facet of the beam
class Top(SubDomain):
    def inside(self,x,on_boundary):
        return near(x[1],0.2)

top = Top()  

L = 1
W = 0.2
mu = 1
lambda_ = 1.25

#creating mesh and function space
mesh = BoxMesh(Point(0,0,0),Point(L,W,W),10,2,2)
V = VectorFunctionSpace(mesh,"P",1)

boundries = FacetFunction("size_t", mesh, mesh.topology().dim() - 1)
boundries.set_all(0)
top.mark(boundries,1)
ds = Measure('ds', domain=mesh, subdomain_data=boundries)

#defining boundary conditions:
tol = 1E-14  

def fixed_boundary_left(x,on_boundary):
    return on_boundary and x[0]<tol and x[1]<tol

def roller_boundary_bottom(x,on_boundary):
    return on_boundary and x[0]==1 and x[1]<tol

bc_fixed=DirichletBC(V,Constant((0,0,0)),fixed_boundary_left)    
bc_roller=DirichletBC(V.sub(1),Constant(0),roller_boundary_bottom)

bc = [bc_fixed,bc_roller]

def epsilon(u):
    strain_u=0.5*(nabla_grad(u)+nabla_grad(u).T)
    return strain_u

def sigma(u):
    stress_u=lambda_*nabla_div(u)*Identity(d)+2*mu*epsilon(u)
    return stress_u

#defining the variational formulation
u = TrialFunction(V)
v = TestFunction(V)
d = u.geometric_dimension()  
T = Constant((0,0,0))  
f = Constant((0.0,-0.2,0.0))   
a = inner(sigma(u),epsilon(v))*dx
L = dot(T,v)*dx+dot(f,v)*ds

#computing the solution
u = Function(V)
solve(a==L,u,bc)

plot(u,title='Displacement',mode='displacement')
vtkfile_u = File('Displacement.pvd')
vtkfile_u << u​


Community: FEniCS Project

1 Answer


3
19 days ago by
Running your code, I don't get a JIT error.
I only get the warning that no domains for the BCs
bc_fixed=DirichletBC(V,Constant((0,0,0)),fixed_boundary_left)    
bc_roller=DirichletBC(V.sub(1),Constant(0),roller_boundary_bottom)

are found.

This is because your domain definitions only include one node. Add the pointwise kwarg to the definition

bc_fixed=DirichletBC(V,Constant((0,0,0)),fixed_boundary_left, method='pointwise')
bc_roller=DirichletBC(V.sub(1),Constant(0),roller_boundary_bottom, method='pointwise')

*and* drop the on_boundary since it defaults to false for pointwise BCs:

def fixed_boundary_left(x,on_boundary):
    return x[0]<tol and x[1]<tol

def roller_boundary_bottom(x,on_boundary):
    return x[0]==1. and x[1]<tol
Hello Klunkean,

Thank you very much for your guidance. It did solve my problem.

I have one more quick question:
I did not understand your comment: "This is because your domain definitions only include one node. "
Will it not include all nodes on the left bottom edge? can you please explain a bit more about this?

Thank you very much.
written 19 days ago by Aby  
1
Yeah I didn't think that through enough. It's obviously all nodes on the edge. But the reasoning still applies (from the docs):

"Available methods are: topological approach (default), geometric approach, and pointwise approach. The topological approach is faster, but will only identify degrees of freedom that are located on a facet that is entirely on the boundary."
written 19 days ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »