### Definition of Elemental Stiffness Matrix on FeniCs

265

views

-2

Hello all,

I could not find any detailed information about that how FEniCS defines the elemental stiffness matrix and assemblies them. In the elasticity tutorial (ft_06), the elastic moduli is not defined in the code. There were defined gamma and lambda for stress and strain definitions. Without knowing the elastic moduli, how FeniCS calculates the elemental stiffness matrix ? ( $\int\left[B\right]^T\left[C\right]\left[B\right]dV$∫[

Regards,

I could not find any detailed information about that how FEniCS defines the elemental stiffness matrix and assemblies them. In the elasticity tutorial (ft_06), the elastic moduli is not defined in the code. There were defined gamma and lambda for stress and strain definitions. Without knowing the elastic moduli, how FeniCS calculates the elemental stiffness matrix ? ( $\int\left[B\right]^T\left[C\right]\left[B\right]dV$∫[

`B`]^{T}[`C`][`B`]`d``V`)Regards,

Community: FEniCS Project

### 6 Answers

1

The constitutive relation in the tutorial is given by density strain energy function relation; it is not the way you are trying to define it.

$\varphi\left(E\right)=\frac{\lambda}{2}I_1\left(E\right)^2+\mu I_2\left(E\right)$

$\sigma=\frac{\partial\varphi}{\partial\epsilon}$

If you want to get the material constitutive tensor you can get it as follows:

$\varphi\left(E\right)=\frac{\lambda}{2}I_1\left(E\right)^2+\mu I_2\left(E\right)$

`φ`(`E`)=λ2`I`_{1}(`E`)^{2}+`μ``I`_{2}(`E`)$\sigma=\frac{\partial\varphi}{\partial\epsilon}$

`σ`=∂`φ`∂`ϵ`and $C_{ijkl}=\frac{\partial\sigma}{\partial\epsilon}$`C`_{ijkl}=∂`σ`∂`ϵ`If you want to get the material constitutive tensor you can get it as follows:

`Cijkl = diff(sigma(u), variable(epsilon(u))) `

and to have it locally, you have to project it into the proper tensor space function, for elemental space probably you will need 'DG' 0 spaces.

To get the local elemental stiffness matrix, you will need to define your problem as follows:

```
V2 = FiniteElement(mesh, 'CG', 2)
u = Function(V2, name='Displacement')
# internal potentials
PI_int = varphi * dx
# external potentials
PI_ext = inner(Constant(tp_load), u) * ds
PI_p = PI_int - PI_ext
# first variation of the equilibrium
dPI_p = derivative(PI_p, u, TestFunction(V2))
# Second variation of the equilibrium, tangent problem
d2PI_p = derivative(dPI_p, u, TrialFunction(V2))
# Linearization and Interpolation of the tangent problem (global stiffness matrix)
Kt, ft = assemble_system(d2PI_p, dPI_p, bcs)
##### local tangent stiffness matrix (elemental)
# element iterators
cell_iter = cells(mesh)
# local stiffness assembly function
def kte_calc(tp_cell_e):
return assemble_local(d2PI_p, tp_cell_e[1])
# list of local stiffness matrix
list_Kt_e = map(kte_calc, enumerate(cell_iter))
```

1

This may help you .... to write $\int\left[B\right]^T\left[C\right]\left[B\right]dV$∫[

`B`]^{T}[`C`][`B`]`d``V````
from dolfin import *
import numpy as np
def D_Matrix(E, nu, stressState):
if stressState=='PlaneStress':
a = E/(1-nu**2)
D = np.array([[a, a*nu, 0],
[a*nu, a, 0],
[0, 0, a*(1-nu)/2],
])
else:
a = E/((1+nu)*(1-2*nu))
D = np.array([[a*(1-nu), a*nu, 0],
[a*nu, a*(1-nu), 0],
[0, 0, a*(0.5-nu)],
])
return D
# Material properties....
E = 10.0
nu = 0.2
stressState = 'PlaneStrain'
D_mat = D_Matrix(E,nu,stressState)
D_mat = as_matrix(D_mat)
mesh = UnitSquareMesh(10,10)
V = VectorFunctionSpace(mesh,'CG',2)
u,v = TrialFunction(V), TestFunction(V)
def eps(u):
""" Returns a vector of strains of size (3,1) in the Voigt notation
layout {eps_xx, eps_yy, gamma_xy} where gamma_xy = 2*eps_xy"""
return as_vector([u[i].dx(i) for i in range(2)] +[u[i].dx(j) + u[j].dx(i) for (i,j) in [(0,1)]])
# Element stiffness matrix B'DB....
a = inner(eps(v),D_mat*eps(u))*dx
```

0

In elasticity, the stress tensor:

\( \sigma_{ij} = C_{ijkl}\varepsilon_{kl} \)

is defined by the stiffness tensor \( C_{ijkl} \), which is given by

\( C_{ijkl}=\lambda \delta_{ij}\delta_{kl} + \mu (\delta_{ik}\delta_{jl}+ \delta_{il}\delta_{jk} ) \)

for an isotropic material (like steel, aluminum). The Lame coefficients \(\lambda, \mu\) are related to the Young's modulus \(E\) and shear modulus \(G\) as follows:

\( \lambda = \frac{(E-2G)G}{3G-E} , \)

\(\mu = G . \)

\( \sigma_{ij} = C_{ijkl}\varepsilon_{kl} \)

is defined by the stiffness tensor \( C_{ijkl} \), which is given by

\( C_{ijkl}=\lambda \delta_{ij}\delta_{kl} + \mu (\delta_{ik}\delta_{jl}+ \delta_{il}\delta_{jk} ) \)

for an isotropic material (like steel, aluminum). The Lame coefficients \(\lambda, \mu\) are related to the Young's modulus \(E\) and shear modulus \(G\) as follows:

\( \lambda = \frac{(E-2G)G}{3G-E} , \)

\(\mu = G . \)

0

Of course it is the mathematics of the elasticity. However the

*lambda*and*mu*are just varible names in the code. I can call them as constant1 and constant2 also. In the elasticity tutorial there is not any information about the elastic moduli equation that you also mentioned in your comment. Without defining the the elastic moduli FEniCS does not now stress is inner product of elastic modili and strain. So, therefore I could not understand how could it define the elemental stiffness matrix.0

If I have the correct code in front of me, the lines 41 and 42

```
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
```

tell you that it's the Lamé constants.

Since

$\sigma=C:\varepsilon=\left(\lambda1\otimes1+2\mu\mathbf{1}\right):\varepsilon=\lambda\text{tr}\varepsilon1+2\mu\varepsilon$`σ`=`C`:`ε`=(λ1⊗1+2`μ`**1**):`ε`=λtr`ε`1+2`μ``ε`

0

Ok, I think we all agree that $\sigma=C:\varepsilon$

`σ`=`C`:`ε`but this constitutive equation is not defined in the elasticity tutorial, right ? It was directly given the result of that equation which is the stress. There is only stress and strain definitions :```
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
```

but there is not any information about elastic moduli (C). I am asking that if you do not define it then FEniCS does what is the elastic moduli definition. Somehow it must use it to calculate the elemental stiffness...

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

$\int\sigma\left(\varepsilon\right):\text{grad}\left(v\right)\mathrm{d}x=0$∫

σ(ε):grad(v)dx=0Now just inserting the constitutive law $\sigma=\lambda\text{tr}\varepsilon1+2\mu\varepsilon$

σ=λtrε1+2μεgives$\int\lambda\mathrm{div}\left(u\right)\mathrm{div}\left(v\right)\mathrm{d}x+\int2\mu\text{symgrad}\left(u\right):\text{symgrad}\left(v\right)\mathrm{dx}=0$∫λ

div(u)div(v)dx+∫2μsymgrad(u):symgrad(v)dx=0which is a bilinear form in the displacement $u$

uand the test function $v$v. Thus it can be assembled into a linear matrix-vector equation that can directly be solved for the nodal unknowns. No need for linearization here.For more general cases you start from the balance of linear momentum and linearize the variational form like:

$\int\sigma\left(\varepsilon\right):\text{grad}\left(v\right)\mathrm{d}x\approx\int\sigma_k:\text{grad}\left(v\right)\mathrm{d}x+\int\text{grad}\left(v\right):C_k:\text{grad}\left(\Delta u\right)\mathrm{d}x$∫

σ(ε):grad(v)dx≈∫σ_{k}:grad(v)dx+∫grad(v):C_{k}:grad(Δu)dxand increment the solution $u_{k+1}=u_k+\Delta u$

u_{k+1}=u_{k}+Δu.Only now you really need the tangent operator $C_k$

C_{k}.