# Weak form of the Navier-Cauchy equation
**URL:** https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694
**Category:** variational formulation
**Created:** [July 3, 2020, 4:20pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694 "2020-07-03T16:20:17Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@luc](https://fenicsproject.discourse.group/u/luc)
#### Post date: [July 3, 2020, 4:20pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694/1 "2020-07-03T16:20:17Z")
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?
---
### Author:  [@kamensky](https://fenicsproject.discourse.group/u/kamensky)
#### Post date: [July 3, 2020, 4:51pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694/2 "2020-07-03T16:51:04Z")
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](https://david-kamensky.eng.ucsd.edu/teaching/mae-207-fea-for-coupled-problems) (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).
---
### Author:  [@luc](https://fenicsproject.discourse.group/u/luc)
#### Post date: [July 3, 2020, 7:26pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694/3 "2020-07-03T19:26:40Z")
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?
---
### Author:  [@kamensky](https://fenicsproject.discourse.group/u/kamensky)
#### Post date: [July 3, 2020, 10:13pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694/4 "2020-07-03T22:13:01Z")
Yes, that looks right. A helpful reference might be @bleyerj’s tutorial [here](https://comet-fenics.readthedocs.io/en/latest/demo/modal_analysis_dynamics/cantilever_modal.py.html). (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}.)
---
### Author:  [@nicoguaro](https://fenicsproject.discourse.group/u/nicoguaro)
#### Post date: [July 9, 2020, 3:24pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694/5 "2020-07-09T15:24:11Z")
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.
---
### Author:  [@luc](https://fenicsproject.discourse.group/u/luc)
#### Post date: [July 9, 2020, 7:09pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694/6 "2020-07-09T19:09:23Z")
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](https://scicomp.stackexchange.com/questions/35485/weak-form-of-the-navier-cauchy-equation) 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.
---
### Author:  [@nicoguaro](https://fenicsproject.discourse.group/u/nicoguaro)
#### Post date: [July 10, 2020, 2:50pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694/7 "2020-07-10T14:50:33Z")
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.
---
### Author:  [@luc](https://fenicsproject.discourse.group/u/luc)
#### Post date: [July 11, 2020, 8:02pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694/8 "2020-07-11T20:02:27Z")
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](https://scicomp.stackexchange.com/questions/35485/weak-form-of-the-navier-cauchy-equation), 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
```python
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
```python
#=================================================================
#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!
---
### Author:  [@bhaveshshrimali](https://fenicsproject.discourse.group/u/bhaveshshrimali)
#### Post date: [July 11, 2020, 8:30pm UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694/9 "2020-07-11T20:30:24Z")
As for your error, it’s here
> [@luc](#):
>
> ```auto
> 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`
```python
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.
---
### Author:  [@luc](https://fenicsproject.discourse.group/u/luc)
#### Post date: [July 12, 2020, 1:24am UTC](https://fenicsproject.discourse.group/t/weak-form-of-the-navier-cauchy-equation/3694/10 "2020-07-12T01:24:12Z")
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.