Derivative of stress field with respect to the displacement field

1

Hello everyone, I want to compute the derivative of the stress field with respect to the displacement field. But I am getting the zero value of the derivative everywhere in the mesh.
I am using the following governing differential equation of linear elasticity as:

The equations governing small elastic deformations of a body \Omega can be written as

\begin{aligned} -\nabla \cdot \sigma &=f \text { in } \Omega, \\ \sigma &=\lambda \operatorname{tr}(\varepsilon) I+2 \mu \varepsilon, \\ \varepsilon &=\frac{1}{2}\left(\nabla u+(\nabla u)^{\top}\right), \end{aligned}

where \sigma is the stress tensor, f is the body force per unit volume, \lambda and \mu are Laml’e’s elasticity parameters for the material in \Omega, I is the identity tensor, tr is the trace operator on a tensor, \varepsilon is the symmetric strain-rate tensor (symmetric gradient), and u is the displacement vector field. We have here assumed isotropic elastic conditions.
In matrix notation, I have the following relation-
\sigma = CBu

\frac{\partial \sigma}{\partial u}=CB
where C=constitutive matrix

B= strain displacement matrix

u=displacement vector
The minimal working code is:

from dolfin import *
L, H = 3, 1
Nx, Ny = 30, 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, "crossed")

eps = lambda _u:sym(grad(_u))
sigma = lambda _u: 2.0 * mu * sym(grad(_u)) + lmbda * tr(sym(grad(_u))) * Identity(len(_u))

E, nu, rho_g = Constant(1e5), Constant(0.3), 1e-3

mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

f = Constant((0, -rho_g))

U = VectorFunctionSpace(mesh, 'Lagrange', degree=1)
u, v = TrialFunction(U), TestFunction(U)
a = inner(sigma(u), eps(v))*dx
l = inner(f, v)*dx

left = CompiledSubDomain("near(x[0], 0.) && on_boundary")
bc = DirichletBC(U, Constant((0.,0.)), left)
u_sol = Function(U, name="Displacement")
solve(a == l, u_sol, bc)

# Find the derivative of stress with respect to displacement
u_var = variable(u_sol)
stress = sigma(u_var)
dstress_du = diff(stress, u_var)
dstress_du[0,0,0](1,1)

Can anyone help me how to compute the derivative \frac{\partial \sigma}{\partial u}? Thanks in advance.

1 Like
2

Hi @dokken and Dr.@kamensky,
Can you please have a look at it? I will really appreciate it.

3

The stress depends on spatial derivatives of \mathbf{u}, not its point values, so a zero partial derivative is what you expect. If you see something like "\partial\sigma / \partial u = CB" in another reference, that is likely referring to partial derivatives of the stress with respect to displacement degrees of freedom (i.e., the entries of u_sol.vector()). The strain–displacement matrix often appears in an intermediate step of the algorithm for assembling the system of algebraic equations for linear elasticity, but the code to perform assembly is automatically generated in the standard FEniCS workflow, so it would be uncommon to need direct access to the strain–displacement matrix.

4

Thanks for the response Dr. @kamensky. I know that in the normal FEniCS workflow it is quite un-natural to get the strain displacement matrix. I am trying to solve the problem of stress constrained topology optimisation and would like to get the derivative of von-mises stress with respect to the design variable. With the application of chain rule, I get the following:
\frac{\partial \sigma_{vm}}{\partial \phi} = \frac{\partial \sigma_{vm}}{\partial \sigma}\frac{\partial \sigma}{\partial \boldsymbol{u}}\frac{\partial \boldsymbol{u}}{\partial \phi}
I have the solution for the first and last term, but I am unable to figure out a way to find the middle term \frac{\partial \sigma}{\partial \boldsymbol{u}}.
Can you please suggest a suitable path forward. I am already looking into the documentation of dolfin-adjoint, but would like to have my custom implementation in standard FEniCS for better understanding of the adjoint method.