Heterogeneous elastic properties

1

Hello everyone,

I’m new to FENICS and I haven’t managed to find a solution to my problem yet, thus I would like to request your help. I would like to extend the example ft06 about linear elasticity (https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft06_elasticity.py) to consider the case in which the elastic constants mu and lambda are not spatially constant.

More specifically, my final goal is to have elastic properties which are homogeneous within each finite element, but different from one element to another. I would like to know your suggestions about which is the best way to do this with FENICS.

Thank you

2

Hi, there are several ways of doing this. You could for instance define spatially varying subdomains, or use a DG-0 function.

3

Hi dokken, thank you for your help, this is exactly what I needed. But this raises a new question, how can I get the full stress tensor in a multiple-subdomain problem? (for e.g., plotting a certain component as a 2D plot). With a single domain, I was doing:

T = TensorFunctionSpace(mesh, “Lagrange”, 1)
stress = project(sigma(u), T)
[s11, s12, s21, s22] = stress.split(True)

and then I can access a certain stress component as s12.vector(). But now, the function to compute the stress is sigma(v,i), where i is the subdomain index. Could you suggest me a way around this? Thanks.

4

I dont follow, project(sigma(u), T) is the projection to the full domain. Thus you Get the component for every cell.

5

Now I’m using subdomains, and for that I’m doing the following:

def sigma(v, i):
    lmbda = 2*G[i]*nu/(1-2*nu)
    return lmbda*tr(eps(v))*Identity(2) + 2.0*G[i]*eps(v)

V = VectorFunctionSpace(mesh, 'Lagrange', 1)
du = TrialFunction(V)
u_ = TestFunction(V)
dx = Measure("dx")(subdomain_data=subdomains)
a = inner(sigma(du,0), eps(u_))*dx(0) + inner(sigma(du,1), eps(u_))*dx(1)

The function sigma(v, i) computes the stress in the domain i, characterized by elastic properties defined by G[i]. With this approach, I cannot use this code, which I used for a single-domain problem:

T = TensorFunctionSpace(mesh, “Lagrange”, 1)
stress = project(sigma(u), T)

because project(sigma(u), T) does not consider the subdomains. I’m looking for some similar way of computing the stress, which takes into account my partitioning into subdomains.

Summary

This text will be hidden

6

What is G[i] in this case? Please produce a Minimal work example as described here.

7

G[i] is the shear modulus of the subdomain i.

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


Nx = 10
Ny = 10
Lx = Nx
Ly = Ny
mesh = RectangleMesh(Point(0., 0.), Point(Lx, Ly), Nx, Ny, "crossed")

class Inclusion(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (2, 6)) and between(x[0], (2, 6)))

inclusion = Inclusion()    
subdomains = MeshFunction("size_t", mesh, 2)
subdomains.set_all(0)
inclusion.mark(subdomains, 1)



nu = Constant(0.3)
G = [Constant(1), Constant(10)]

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

def sigma(v,i):
    lmbda = 2*G[i]*nu/(1-2*nu)
    return lmbda*tr(eps(v))*Identity(2) + 2.0*G[i]*eps(v)

f = Constant((0,0))
V = VectorFunctionSpace(mesh, 'Lagrange', 1)
du = TrialFunction(V)
u_ = TestFunction(V)
dx = Measure("dx")(subdomain_data=subdomains)
a = inner(sigma(du,0), eps(u_))*dx(0) + inner(sigma(du,1), eps(u_))*dx(1)
l = inner(f, u_)*dx(0) + inner(f, u_)*dx(1)



def bottom(x, on_boundary): return near(x[1],0.)
def top(x, on_boundary): return near(x[1],Ly)
def left(x, on_boundary): return near(x[0],0.)
def right(x, on_boundary): return near(x[0],Lx)

dirichelt_BCs  = [DirichletBC(V, Constant((0.,0.)), bottom), 
                  DirichletBC(V, Constant((1.,0.)), top), 
                  DirichletBC(V.sub(1), Constant(0.), left),
                  DirichletBC(V.sub(1), Constant(0.), right)]
     

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

T = TensorFunctionSpace(mesh, "Lagrange", 1)
stress = project(sigma(u), T)
[s11, s12, s21, s22] = stress.split(True)
8

If you instead of using subdomains, use a piecewise constant function to define G, it is easy to obtain the stresses. See for instance How to define different materials / importet 3D geometry (gmsh).

9

Thank you, this is exactly what I need.