Computing stress, strain integrals over a domain

Question

Hello

I have been successful at creating a linear elasticity solver in FEniCS. I compute the displacement field in a Function object u_sol.

I am now interested in calculating volume averaged stress and strain over the domain. In essence I need to calculate (typed in Latex notation) $\langle \sigma \rangle = \frac{1}{\Omega} \int_{\Omega}{\sigma d \Omega}$.

I have the stress defined as a form as follows:

def sigma(u):
    eps = epsilon(u)
    stress = 2.0*Mu*eps + Lambda*tr(eps)*Identity(u.geometric_dimension())
    return stress

def epsilon(u):
    strain = 0.5*(nabla_grad(u) + nabla_grad(u).T)
    return strain

Is it okay for me to write avgStress_11 = assemble(sigma(u_sol)[0,0]*dx) to compute the average stress (component by component)? I tried this, but the results I obtain are numerically wrong even though the displacement field seems to be okay.

I can provide more details of the code as, and if, required.

Any help will be much appreciated.

Thanks

Answer 1

Maybe try to do a projection first, i.e.

sigmaTensor = project(sigma(u_sol),TensorFunctionSpace(mesh,'DG',0))

assuming you are using CG1 for u_sol. I also think that

sym(grad(u))

is equivalent to your epsilon function.

Comments

Thanks KristianE.

Could you kindly also tell how to set up the integral?
I tried doing the following steps, but I get an error message saying Can only integrate scalar expressions. The integrand is a tensor expression with value rank 2 and free indices ()

T = TensorFunctionSpace(mesh,'DG',0)
sigmaTensor = project(sigma(u_sol),T)
sigmaInt = assemble(sigmaTensor*dx)

Also, mathematically, what is the project command doing?

---

The projection does Galerkin projection. If you are using CG1 for u, you should not lose any information by projecting to DG0.

I think you can use brackets as well as the sub function, where
sigmaTensor[0,0] is the same as sigmaTensor.sub(0),
sigmaTensor[0,1] is the same as sigmaTensor.sub(1),
sigmaTensor[1,0] is the same as sigmaTensor.sub(2) and
sigmaTensor[1,1] is the same as sigmaTensor.sub(3).

---

KristianE, I don't get the point of your answer. Mathematically, the stress sigma(u) should already be in DG0, according to your assumption u \in CG1. So the projection should only change the type of sigma(u) from a ufl.something to a dolfin.Function, but not the actual value of the assembled integral, or am I missing something?

So, debmukh, my answer would be that your original approach should have worked and that maybe something else is wrong with your code, which we can only help you with if you post it!

---

Gregor, you are absolutely right. I thought,

from dolfin import *
mesh = RectangleMesh(0,0,1,1,3,3)
u = Function(VectorFunctionSpace(mesh,'CG',1))
assemble(sym(grad(u))[0,0]*dx)

would produce an error message, but it does not. That is what you get for not providing a minimal example with your question :)
I do not see the point of calculating the average stress component by component though. Normally one is interested in the maximum value of the von mises stress.