Giving every element a separate Youngs Modulus in linear elasticity problem

1

Hey there fellow friends,
i would like to make a little tool, so that i can calculate the displacements of nodes in a 2D truss system. It is working alright, but now i would like the elements of my system to have different Youngs Moduli. The mesh im using only has 6 elements so it’s not particularly complex. So far this is my code:

import numpy as np
from dolfin import *

E = 1000
nu = 0.25                                       # Poisson's ratio
mu = E/(2*(1+nu))                               # 2. Lame's parameter
lmbda = Constant((E*nu)/((1+nu)*(1-2*nu)))      # 1. Lame's parameter
mesh = "files/mesh/truss.xml"

mesh = Mesh(mesh)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

def cb1(x, on_boundary):
        return near(x[1], -2.0, 0.1)

bc1 = DirichletBC(V, Constant((0, 0)), cb1)

dx = Measure("dx")

def sigma(u, E):
        d = u.geometric_dimension()
        I = Identity(d)
        return lmbda * tr(epsilon(u)) * I + 2 * mu * epsilon(u)
    
def epsilon(z):
        return 0.5 * (nabla_grad(z) + nabla_grad(z).T)

du = TrialFunction(V)
U = Function(V)
u = U.vector()
v = TestFunction(V)
d = U.geometric_dimension()
I = Identity(d)
E = 700
f = Constant((0, 0))
a = inner(sigma(du, E), epsilon(v)) * dx
L = inner(f, v) * dx
K, f_vec = assemble_system(a, L, bc1)
f_vec[4] = 400
f_vec[5] = -400
solve(K, u, f_vec)
u = np.array(u)
print(u)

I think i have to get some kind of loop going somewhere in:

a = inner(sigma(du, E), epsilon(v)) * dx

but maybe some of you guys have some experience with this.
Thanks for the help and best regards!

2

Could you provide the mesh file or an image for me to test?

3

You can use a discontinuous Galerkin FunctionSpace of order 0 to model such things, i.e., you define

DG0 = FunctionSpace(mesh, 'DG', 0)
lmbda = Function(DG0)

and then you have to set the DoF’s of lambda, which correspond to the constant value on each triangle / tetrahedron, to the appropriate value. For this, you might want to take a look here: https://fenicsproject.org/qa/5476/accessing-the-cell-index-of-a-given-dof/
or the dofmap

DG0.dofmap()