### How to parameterize the linear elasticity example from tutorial?

348
views
0
5 months ago by
Hello all,

I am new to FeNiCS and beginner at coding in general. I was able to follow the Linear Elasticity example from the tutorial (3.3), and now I am trying to modify the code such that I can model the displacement according to variable values of young's modulus/poisson ratio. I am new to coding, so I am not sure how I can begin to change my constants into parameterized variables. Could someone show me an example program where this is done, or give me a push in the right direction?

Thank you,
Joshua
Community: FEniCS Project
Do you want the material parameters to vary in space or dependent on the strain/displacement or do you want something completely different?
written 5 months ago by klunkean
The second one, I would like the material parameters to be dependent on the displacement in the z direction.
written 5 months ago by Joshua Lee

0
5 months ago by
If you can express your material parameter as a closed mathematical expression, you can just use python functions. Note, however, that the problem will become non-linear. Consider the following MWE:
from fenics import *

mesh = UnitCubeMesh(4,4,4)
FE = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
VFE = VectorElement(FE)
V = FunctionSpace(mesh, VFE)

u = Function(V)
delu = TrialFunction(V)
v = TestFunction(V)

def bottom(x, on_boundary):
return near(x[2],0.)
def top(x, on_boundary):
return near(x[2],1.)

dbc = [DirichletBC(V, Constant((0.,0.,0.)), bottom),
DirichletBC(V.sub(2), Constant(.1), top)]

eye = Identity(3)
eyOey = outer(eye,eye)
eye4 = as_tensor(.5*(eye[i,k]*eye[j,l]+eye[i,l]*eye[j,k]), (i,j,k,l))
D4 = eye4 - 1./3.*eyOey

# bulk and shear modulus, shear modulus as function of u
K = 160e+9

def G(u):
# just an example
G0 = 80e+9
uzMax = .1
return G0*(1.-u[2]/uzMax)

#stiffness tensor
def C4(u):
return K*eyOey + G(u)*D4

#strain
#stress
def sigma(u):
return as_tensor(C4(u)[i,j,k,l]*eps[k,l], (i,j))

#nonlinear form

#gateaux derivative
J = derivative(F, u, delu)

solve(F==0, u, dbc, J=J)

File('disp.pvd') << u​
Thank you! This is very helpful. Could you please explain what the

"eye = Identity(3)
eyOey = outer(eye,eye)
eye4 = as_tensor(.5*(eye[i,k]*eye[j,l]+eye[i,l]*eye[j,k]), (i,j,k,l))
D4 = eye4 - 1./3.*eyOey"

part of the code is doing?

I also tried running this code through Docker, but I got the error message:

"SyntaxError: Non-ASCII character '\xe2' in file ft06_elasticity_parameterized.py on line 52, but no encoding declared; see http://python.org/dev/peps/pep-0263/ for details"

I tried following the link, but I just got sent to a 404. Do you know why I got this error?

Thank you!

written 5 months ago by Joshua Lee
The four lines define the following quantities:
eye - Second order unit tensor  $\mathbf{1}$1
eyOey - Dyadic product of second order unit tensor with itself $\mathbf{1}\otimes\mathbf{1}$11  (Trace projector)
eye4 -  Fourth order symmetric unit tensor $\mathbf{1}^{4,\mathrm{s}}=\frac{1}{2}\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)\mathbf{e}_i\otimes\mathbf{e}_j\otimes\mathbf{e}_k\otimes\mathbf{e}_l$14,s=12 (δikδjl+δilδjk)eiejekel
D4 - Fourth order deviatoric projector  $\mathbf{D}^{4,\mathrm{s}}=\mathbf{1}^{4,\mathrm{s}}-\frac{1}{3}\mathbf{1}\otimes\mathbf{1}$D4,s=14,s13 11

With these quantities we can write Hooke's tangent tensor as
$\mathbf{C}^4=K\mathbf{1}\otimes\mathbf{1}+G\mathbf{D}^{4,\mathrm{s}}$C4=K11+GD4,s

Where  $K$K  and  $G$G  are bulk and shear modulus respectively.

As for your error: This is a simple copy-paste error. I don't know what happens in detail, but try to go to the very last line of the pasted code and delete the very last character. There is some invisible character which gets pasted and does not belong there.
written 5 months ago by klunkean
Hi klunkean,

Could you help me with another issue? I am having trouble extracting the sigma function, similar as to how the original tutorial plotted the stress intensity. I would like to plot the stress intensity as a function of displacement.

written 4 months ago by Joshua Lee
1
Hi again.

If your stress is declared as above you could project it onto a finite element space after computing the displacements. Like this:
solve(F==0, u, dbc, J=J)
tensorFE = TensorElement(FiniteElement('DG', mesh.ufl_cell(), 0))
tensorFunSpace = FunctionSpace(mesh,tensorFE)

stress = project(sigma(u), tensorFunSpace)
stressFile = File('stress.pvd')
stressFile << stress​

Then you can view stress.pvd in paraview. Hope that helps.
written 4 months ago by klunkean