Definition of Elemental Stiffness Matrix on FeniCs

4 months ago by
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$[B]T[C][B]dV )

Community: FEniCS Project

6 Answers

4 months ago by
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)$φ(E)=λ2 I1(E)2+μI2(E)
$\sigma=\frac{\partial\varphi}{\partial\epsilon}$σ=φϵ and   $C_{ijkl}=\frac{\partial\sigma}{\partial\epsilon}$Cijkl=σϵ

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))

4 months ago by
This may help you .... to write  $\int\left[B\right]^T\left[C\right]\left[B\right]dV$[B]T[C][B]dV 
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],
        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
4 months ago by
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 . \)
4 months ago by
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.
4 months ago by
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.

4 months ago by
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...

I just want to add to Ricardo D. Lahuerta's answer, that in a fully linear problem you don't need to consider the tangent problem, i.e. the linearized version of the weak form of the balance of linear momentum. You directly get a system of linear equations. Let's take a look at this. Consider the simplest weak form of equilibrium:


Now just inserting the constitutive law  $\sigma=\lambda\text{tr}\varepsilon1+2\mu\varepsilon$σ=λtrε1+2με  gives


which is a bilinear form in the displacement  $u$u and 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):Ck:grad(Δu)dx

and increment the solution  $u_{k+1}=u_k+\Delta u$uk+1=uk+Δu  .
Only now you really need the tangent operator  $C_k$Ck .
written 4 months ago by klunkean  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »