Whatever Young's modulus is, the displacement remain the same in MBB beam problem.

Question

I try to test a MBB beam problem. It is a linear elastic problem in a rectangle domain. And I try to choose different Young's modulus to get the maximum displacement of each problem. However the maximum displacement seems to be the same whatever the young's modulus is. Could you please help me. I really confused. THANKS

""" MBB beam. Impose zero displacement at the centers of the left and the right side, and a point load on the middle of the top."""

from dolfin import *

Create mesh and define function space

mesh = RectangleMesh(0.0,0.0,4.0,1.0,160,40)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

Mark boundary subdomians

boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundary_parts.set_all(0)

class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return abs(x[0]) < tol and abs(x[1]-0.5) < tol
Gamma_L = LeftBoundary()

class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return abs(x[0]-4) < tol and abs(x[1]-0.5) < tol
Gamma_R = RightBoundary()

class PointLoad(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return abs(x[0]-2) < tol and abs(x[1]-1)<tol
Gamma_P = PointLoad()
Gamma_P.mark(boundary_parts, 3)

Define functions

u = TrialFunction(V) # Trial displacement
v = TestFunction(V) # Test function
B = Constant((0.0, 0.0)) # Body force per unit volume
T = Constant((0.0, -0.2)) # Point load on the boundary

Kinematics

d = u.geometric_dimension()
I = Identity(d) # Identity tensor
stress = sym(grad(u)) # Green Lagrange tensor

Elasticity parameters

E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2(1 + nu))), Constant(Enu/((1 + nu)(1 - 2nu)))

stress

def sigma(v):
return 2.0musym(grad(v)) + lmbdatr(sym(grad(v)))I

Define Dirichlet boundary (x = 0 or x = 1)

c = Expression(("0.0", "0.0"))
bcl = DirichletBC(V, c, Gamma_L, method='pointwise')
bcr = DirichletBC(V, c, Gamma_R, method='pointwise')
bcp = DirichletBC(V, T, Gamma_P, method='pointwise')
bcs = [bcl, bcr, bcp]

Total potential energy

dss = ds[boundary_parts]
a = inner(sigma(u), sym(grad(v)))dx
L = inner(B, v)dx + inner(T, v)*dss(3)

Solve variational problem

A, b = assemble_system(a, L, bcs)
u = Function(V)
U = u.vector()
solve(A, U, b)

Save solution in VTK format

file = File("displacement.pvd");

file << u;

Plot and hold solution

print U.array().max()
plot(u, mode = "displacement", interactive = True)

Comments

Hi, Did you solve the problem finally? I am stuck on the same thing.

As per below listed suggestions, when I apply pointwise Drichilet BC on the point, I get a displacement "resulting" from the boundary condition in the viper plot.

So, if the BC is (0,-7,0), the max displacement in the viper plot is near 7. So I feel that I am implying displacement instead of load.

Also, what ChrishZhang says about unchanging displacement w.r.t. elastic modulus but changing displacement w.r.t. point load is correct. Kristian says their is some kind of normalising function running in the background, is that correct?

If your code is verified and working now, pls. post it here to help end this thread. Thanks a lot!

The same code (posted in the question) in proper format for anyone who would like to run it and take a look:

from dolfin import *

#Create mesh and define function space

o = Point(0, 0)
l = Point(4, 1)
mesh = RectangleMesh(o, l, 160, 40)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

#Mark boundary subdomians

boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
boundary_parts.set_all(0)

#return (abs(x[0]) < tol and abs(x[1]-0.5) < tol)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0],0.0) and near(x[1],0.5)  
Gamma_L = LeftBoundary()

tol = 1E-14 # tolerance for coordinate comparisons
#return near(x[0],4.0) and near(x[1],0.5)
class RightBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0]-4) < tol and abs(x[1]-0.5) < tol  
Gamma_R = RightBoundary()

class PointLoad(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0]-2) < tol and abs(x[1]-1)<tol
Gamma_P = PointLoad()
Gamma_P.mark(boundary_parts, 3)

#Define functions

u = TrialFunction(V) # Trial displacement
v = TestFunction(V) # Test function
B = Constant((0.0, 0.0)) # Body force per unit volume
T = Constant((0.0, -0.2)) # Point load on the boundary

#Kinematics

d = u.geometric_dimension()
I = Identity(d) # Identity tensor
stress = sym(grad(u)) # Green Lagrange tensor

#Elasticity parameters

E, nu = 10.0, 0.3
mu, lmbda = E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 - 2.0*nu))

#stress

sigma = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(v.cell().d)     

#Define Dirichlet boundary (x = 0 or x = 1)

c = Expression(("0.0", "0.0"))
bcl = DirichletBC(V, c, Gamma_L, method='pointwise')
bcr = DirichletBC(V, c, Gamma_R, method='pointwise')
bcp = DirichletBC(V, T, Gamma_P, method='pointwise')
bcs = [bcl, bcr, bcp]

#Total potential energy

dss = ds[boundary_parts]   
a = inner(sigma, grad(v))*dx        
L = inner(B, v)*dx + inner(T, v)*dss(3)

#Solve variational problem

A, b = assemble_system(a, L, bcs)
u = Function(V)
U = u.vector()
solve(A, U, b)

#Save solution in VTK format

file = File("displacement.pvd");

file << u;

#Plot and hold solution

print U.array().max()
plot(u, mode = "displacement", interactive = True)

Answer 1

Hi,
when you write

bcp = DirichletBC(V, T, Gamma_P, method='pointwise')

and then later

A, b = assemble_system(a, L, bcs)

you are imposing a displacement at the point PointLoad, since it is Dirichlet boundary conditions, so it is not what you want. What you want is Neumann conditions. The problem is, I don't know how to apply Neumann conditions at only one point. It seems one could use the function PointSource but I don't know how to use it for vector-valued functions, so I have the same problem as you.

One solution I found which is not very satisfying is to define a small boundary (but larger than just a point) with Neumann conditions, for instance in this way:

class PointLoad(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0]-2) < 0.1 and abs(x[1]-1)<0.1

then you apply your Neumann conditions as you did

L = inner(B, v)*dx + inner(T, v)*dss(3)

However in your case there was a problem because you had both Dirichlet and Neumann conditions at the same point, which means you need to completely remove the following condition

DirichletBC(V, T, Gamma_P, method='pointwise')

In fact your Neumann condition was inactive because the Neumann boundary was reduced to a point.

Still, I would be glad if someone could explain how to apply PointSource to a vector function because so far my attempts have failed.
For instance if I define

V = VectorFunctionSpace(mesh, 'CG', 1)
delta = PointSource(V.sub(1), Point(0.0,-1.0), -1)
delta.apply(b) 

I get:
TypeError: in method 'PointSource_apply', argument 2 of type 'dolfin::GenericVector &'

Comments

This is what you want.