Weak form of the Navier-Cauchy equation

[luc]

Hello! I am trying to obtain the weak form of the Navier-Cauchy equation, which is

- \rho \omega ^2 \textbf{U} - \mu \nabla ^2 \textbf{U} - (\mu + \lambda) \nabla (\nabla \cdot \textbf{U}) = \textbf{F}

and can be written in the component form

-(2 \mu +\lambda) \frac{\partial ^2 U_1}{\partial x_1 ^2} - \mu \frac{\partial ^2 U_1}{\partial x_2 ^2} - (\mu + \lambda) \frac{\partial ^2 U_2}{\partial x_1 \partial x_2} - \rho \omega ^2 U_1 = F_1

-(2 \mu +\lambda) \frac{\partial ^2 U_2}{\partial x_2 ^2} - \mu \frac{\partial ^2 U_2}{\partial x_1 ^2} - (\mu + \lambda) \frac{\partial ^2 U_1}{\partial x_1 \partial x_2} - \rho \omega ^2 U_2 = F_2

The general procedure is to multiply the PDE by a test function \textbf{v} in the space \textbf{V}, or v in the space V, and integrate it over the domain \Omega. I will proceed with the component form, for I believe it is easier for me to understand. Setting \textbf{F} = 0 and rearranging the terms

-(2 \mu +\lambda) \int_\Omega v \left[ \frac{\partial ^2 U_1}{\partial x_1 ^2} + \frac{\partial ^2 U_2}{\partial x_2 ^2} \right]dxdy - \mu \int_\Omega v \left[ \frac{\partial ^2 U_1}{\partial x_2 ^2} + \frac{\partial ^2 U_2}{\partial x_1 ^2} \right]dxdy -(\mu + \lambda)\int_\Omega v \left[ \frac{\partial ^2 U_2}{\partial x_1 \partial x_2} + \frac{\partial ^2 U_1}{\partial x_1 \partial x_2} \right]dxdy - \rho \omega ^2 \int_\Omega v \left[ U_1+U_2 \right]dxdy = 0

From Green’s theorem I know that
\int_{\Omega} \left(v \frac{\partial ^2 u}{\partial x ^2} \right)dxdy = \int_\Gamma \left(vu\hat{n}_x \right)ds - \int_{\Omega} \left( \frac{\partial v}{\partial x} \frac{\partial u}{\partial x} \right)dxdy

Which is sufficient to deal with the first and second integrals. However, I do not know how to proceed with the cross derivatives \partial ^2 / \partial x_1 \partial x_2 of the third integral. Can someone help me with this?

[kamensky]

I find it easier to think about this by recognizing that

-\mu\nabla^2\mathbf{U} - (\mu + \lambda)\nabla(\nabla\cdot\mathbf{U}) = -\nabla\cdot\pmb{\sigma}\text{ ,}

where

\pmb{\sigma} = \lambda(\nabla\cdot\mathbf{U})\mathbf{I} + 2\mu\operatorname{sym}(\nabla\mathbf{U})\text{ .}

Then you can use the following integration-by-parts identity:

-\int_\Omega(\nabla\cdot\pmb{\sigma})\cdot\mathbf{v}\,d\mathbf{x} = \int_\Omega \pmb{\sigma}:\nabla\mathbf{v}\,d\mathbf{x} - \int_{\partial\Omega}(\pmb{\sigma}\mathbf{n})\cdot\mathbf{v}\,d\mathbf{s}\text{ ,}

where the boundary term is often left out, to weakly enforce the traction-free boundary condition, \pmb{\sigma}\mathbf{n} = \mathbf{0}.

If that integration-by-parts identity is confusing, I wrote up a series of exercises to guide students through its derivation, which can be found in an appendix of the “mathematical background” notes from my graduate course (although, for those stumbling on this link in the far future, I will likely rearrange that page over subsequent iterations of the course, starting in a month or two).

[luc]

Thank you very much for showing this identity @kamensky! I’ll be sure to read through the notes you wrote too. In order for me to write this on FEniCS, I’d like to hear from you if I’m guessing it right. I believe I should thus solve the following problem

- \int_\Omega (\rho \omega ^2 \mathbf{v} \cdot \mathbf{U})d \mathbf{x} -\int_\Omega(\nabla\cdot\pmb{\sigma})\cdot\mathbf{v}\,d\mathbf{x} = 0

which becomes

\int_{\Omega} (-\rho \omega ^2 \mathbf{v} \cdot \mathbf{U} + \pmb{\sigma}:\nabla\mathbf{v} )\,d\mathbf{x} = \int_{\partial\Omega}(\pmb{\sigma}\mathbf{n})\cdot\mathbf{v}\,d\mathbf{s}\text{ .}

Is this right?

[kamensky]

Yes, that looks right. A helpful reference might be @bleyerj’s tutorial here. (If you’re wondering about the discrepancy between the gradient of \mathbf{v} in my post above vs. symmetric gradient of \mathbf{v} in k_form from the tutorial, note that, because \pmb{\sigma} is symmetric, \pmb{\sigma}:\operatorname{sym}(\nabla\mathbf{v}) = \pmb{\sigma}:\nabla\mathbf{v}.)

[nicoguaro]

Did you ask the exact same question in SciComp.Exchange? If so, please provide a link next time so people know and efforts are not doubled.

[luc]

Yes I did @nicoguaro, and I was planning to write about it here, but I haven’t had the time to test both answers in FEniCS yet. I’ll share the link anyway.
I don’t see why ‘doubling efforts’ is a problem though. Both @kamensky and the user who answered me in Sci. Comp. Exchange provided me with valuable answers, which might be equivalent, but I could not figure it out by myself at first sight.

[nicoguaro]

From your perspective is not a problem. But from the perspective of the person that invest some time to answer you it might. I am just saying that letting the community know is good so they can check what have been answered and decide after that.

[luc]

Hello again. So, now I had the time to put it all together in a code. As pointed out by @nicoguaro, I have posted the same question in another forum, which led me to a different variational formulation than the one from this thread.

In order to test both of them, I set up a simulation where there is a Dirac delta function in the middle of the domain (it may not be best approach though), and see if the results are identical.

From the answer provided in this link, I derived another variational formulation:
\int_{\Omega} \left[ (2 \mu + \lambda) (\nabla \cdot \mathbf{U}) \nabla \cdot \mathbf{v} - \rho \omega^2 \mathbf{U} \cdot \mathbf{v} \right] d \Omega = \int_{\partial \Omega} \left[ (2 \mu + \lambda) (\nabla \cdot \mathbf{U}) \mathbf{v} \cdot \mathbf{n} \right]dS

I wrote, where g = (\nabla \cdot \mathbf{U}) \mathbf{n} = 0

from fenics import *
import matplotlib.pyplot as plt
import numpy as np

#Mesh parameters
nx = 100
ny = 100
x_1, y_1 = [-1, -1]
x_2, y_2 = [1, 1]

#Material properties
omega = 1
E, nu = Constant(1e5), Constant(0.)
rho = Constant(1e-3)
mu = E/2./(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

#mesh and function space
mesh = RectangleMesh(Point(x_1, y_1), Point(x_2, y_2), nx, ny, "crossed")
V = VectorFunctionSpace(mesh, 'CG', 1)

#boundary conditions
u0 = Constant((0.0, 0.0))
bc = DirichletBC(V, u0, "on_boundary")

u = TrialFunction(V)
v = TestFunction(V)
#==================================================================
g = Constant((0.0, 0.0))

a = div(u)*div(v)*(2*mu + lmbda)*dx - rho*omega*omega*inner(u,v)*dx
L = (2*mu + lmbda)*inner(g,v)*ds
A, b = assemble_system(a, L, bc)

u = Function(V)
A, b = assemble_system(a, L, bc)
delta = PointSource(V, Point(0., 0.),1)
delta.apply(b)
solve(A, u.vector(), b)

plot(u)
plt.show()

which gives

Although it resembles what I was expecting, I’m not 100% convinced by this result.

Now, from @kamensky’s answer and following some instructions from the above-mentioned tutorial I wrote, where T = \mathbf{\sigma \cdot n} = 0

#=================================================================
#Delta Dirac Function
class Delta(UserExpression):
    def __init__(self, eps, x0, **kwargs):
        self.eps = eps
        self.x0 = x0
        UserExpression.__init__(self, **kwargs) 
    def eval(self, values, x):
        eps = self.eps
        values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
        values[1] = 0
        
    def value_shape(self): return (2, )

T = Constant((0.0, 0.0))
sigma = 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.geometric_dimension())
delta = Delta(eps=1E-4, x0=np.array([0., 0.]), degree=5)

a = (inner(sigma, grad(v)) - rho*omega*omega*inner(u,v))*dx
L = inner(T,v)*ds + inner(delta, v)*dx

u = Function(V)
solve(a == L, u, bc)
plot(u)
plt.show()

Unfortunately, plugging this bit to the main code gives the error message TypeError: cannot unpack non-iterable bool object
what could it be?

Any comments on my approach are much appreciated!

[bhaveshshrimali]

As for your error, it’s here

@luc

sigma = 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.geometric_dimension())
delta = Delta(eps=1E-4, x0=np.array([0., 0.]), degree=5)

a = (inner(sigma, grad(v)) - rho*omega*omega*inner(u,v))*dx

Your bilinear form isn’t correct due to sigma being evaluated on the test function, and not the trial function. Consider changing the v in your expression of sigma to a u or equivalently make it a lambda function, and then call on u

sigma = lambda v: 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.geometric_dimension())
a = (inner(sigma(u), grad(v)) - rho*omega*omega*inner(u,v))*dx

As for the variational formulation, I may need to take a look more carefully.

[luc]

Brilliant @bhaveshshrimali, your suggestion did the trick. Please, take a look at the right square in this figure:

What intrigues me the most is the amplitude of the delta function. In the left it’s 3.6e+6, whereas in the right it’s 3.4e-7. I was expecting it to be 1 in both cases instead.

Anyway, if you would take a look at the variational formulation, that would be great. If not, I believe I should stick to the formulation with the stress tensor \sigma.