Spatially varying elastic modulus in linear elastic beam deflection problem

variational formulation
1

Hi. I am trying to implement a model for a clamped beam based on the tutorial https://comet-fenics.readthedocs.io/en/latest/demo/elasticity/2D_elasticity.py.html , but with a spatially varying elastic modulus (having the Lame parameters be defined as a spatial field that varies along the beam length. I am working with the code below:

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

L = 25.
H = 1.
Nx = 250
Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, 'crossed')

def eps(v):
    return sym(grad(v))

E = Constant(1e5)
nu = Constant(0.3)
model = "plane_stress"

mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)

rho_g = 1e-3
f = Constant((0, -rho_g))

V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)
du = TrialFunction(V)
u_ = TestFunction(V)
d2v = dof_to_vertex_map(V)

def lmbda_x(coords):
    return 30000 + 10.0 * coords[0]

lmbda_func = Function(V)
lmbda_vals = np.apply_along_axis(lmbda_x, 1, mesh.coordinates())
# Want Lame parameter to be the same in x and y direction, so copy the values once.
lmbda_vals = np.hstack([lmbda_vals, lmbda_vals]) 
lmbda_func.vector().set_local(lmbda_vals[d2v])

def sigma(v):
    # return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
    return lmbda_func * tr(eps(v)) * Identity(2) + 2.0*mu*eps(v)

a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx

def left(x, on_boundary):
    return near(x[0], 0.)

bc = DirichletBC(V, Constant((0.,0.)), left)

u = Function(V, name="Displacement")
solve(a == l, u, bc)

plt.figure()
plot(1e3*u, mode="displacement")
plt.savefig('test.png')

print("Maximal deflection:", -u(L,H/2.)[1])
print("Beam theory deflection:", float(3*rho_g*L**4/2/E/H**3))

However, the following error occurs:

Invalid ranks 1 and 2 in product.
Traceback (most recent call last):
  File "/usr/local/lib/python3.6/dist-packages/IPython/core/interactiveshell.py", line 3319, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "<ipython-input-120-463cfb80453d>", line 5, in <module>
    a = inner(sigma(du), eps(u_))*dx
  File "<ipython-input-120-463cfb80453d>", line 3, in sigma
    return lmbda_func*Identity(2) * tr(eps(v))* Identity(2) + 2.0*mu*eps(v)
  File "/usr/local/lib/python3.6/dist-packages/ufl/exproperators.py", line 193, in _mul
    return _mult(self, o)
  File "/usr/local/lib/python3.6/dist-packages/ufl/exproperators.py", line 167, in _mult
    error("Invalid ranks {0} and {1} in product.".format(r1, r2))
  File "/usr/local/lib/python3.6/dist-packages/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Invalid ranks 1 and 2 in product.

which suggests that lmbda_func is of rank 1 and tr(eps(v)) * Identity(2) is of rank 2 (although I didn’t manage to check the latter with .rank() ). If I replace the definition of lmbda_x with the commented line (using a constant Lame parameter), it works nicely.

My question is then: how can I convert the rank 1 lmbda_func into an appropriate form to achieve the spatially varying property? Thanks in advance.

2

I’m not so experienced with continuum mechanics… is there a reason you want \lambda to be vector valued?

Changing lambda_func as follows yields a reasonable result

lmbda_func = Expression("3000.0 + 10.0*x[0]", degree=1)

or equivalently in ufl:

x = SpatialCoordinate(mesh)
lmbda_func = 3000.0 + 10.0*x[0]
2 Likes
3

Thanks for the reply. I am also not an expert in continuum mechanics, but as far as I understand, \lambda should be a scalar (https://en.wikipedia.org/wiki/Lamé_parameters).

As an extension, if I have to define \lambda as a Python function (e.g. the form of \lambda(x) depends on an external numpy array of parameters), what should be the way to modify the code?

4

Assemble it into a scalar Function from a scalar FunctionSpace, rather than a VectorFunctionSpace as your code is currently set up.

1 Like