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

348

views

0

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

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

### 1 Answer

0

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
eps = sym(grad(u))
#stress
def sigma(u):
return as_tensor(C4(u)[i,j,k,l]*eps[k,l], (i,j))
#nonlinear form
F = inner(grad(v),sigma(u))*dx
#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!

"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}$

eyOey - Dyadic product of second order unit tensor with itself $\mathbf{1}\otimes\mathbf{1}$

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$

D4 - Fourth order deviatoric projector $\mathbf{D}^{4,\mathrm{s}}=\mathbf{1}^{4,\mathrm{s}}-\frac{1}{3}\mathbf{1}\otimes\mathbf{1}$

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}}$

Where $K$

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.

eye - Second order unit tensor $\mathbf{1}$

**1**eyOey - Dyadic product of second order unit tensor with itself $\mathbf{1}\otimes\mathbf{1}$

**1**⊗**1**(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$

**1**^{4,s}=12 (`δ`_{ik}`δ`_{jl}+`δ`_{il}`δ`_{jk})`e`_{i}⊗`e`_{j}⊗`e`_{k}⊗`e`_{l}D4 - Fourth order deviatoric projector $\mathbf{D}^{4,\mathrm{s}}=\mathbf{1}^{4,\mathrm{s}}-\frac{1}{3}\mathbf{1}\otimes\mathbf{1}$

`D`^{4,s}=**1**^{4,s}−13**1**⊗**1**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}}$

`C`^{4}=`K`**1**⊗**1**+`G``D`^{4,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.

Thank you for your help!

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.

Thank you for your help!

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:

Then you can view stress.pvd in paraview. Hope that helps.

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

you have been very helpful, thank you again for your help!

written
4 months ago by
Joshua Lee

Please login to add an answer/comment or follow this question.