I am having difficulties implementing a Neumann value when numerically solving the Navier equation using NDEigensystem.

The Navier equation is given by

$\nabla^2 \vec u + (p^2 - 1) \nabla(\nabla\cdot \vec u) =- \omega^2 \vec u$

where $\vec u(x_i,x_j,x_k):\mathbb R^3\rightarrow\mathbb R^3$.

This can be re-expressed as the divergence of a flux term

$ \nabla\cdot (\nabla u_{i} + \partial_{x_i} \vec u + (p^2-2) \nabla.\vec u \hat i) =-\omega^2\vec u_i$ for each component $i$. Assuming zero flux over the boundary implies the following purely Neumann boundary condition.

$ \hat n\cdot (\nabla u_x + \partial_x \vec u + \frac{\lambda}{\mu} \nabla.\vec u)= 0 \\ \hat n\cdot (\nabla u_y + \partial_y \vec u + \frac{\lambda}{\mu} \nabla.\vec u)= 0 \\ \hat n\cdot (\nabla u_z + \partial_z \vec u + \frac{\lambda}{\mu} \nabla.\vec u)= 0 $

I am having difficulty understanding the NeumannValue specification in Mathematica. As I understand the documentation, for a one dimensional case it interprets a BC of the form: $n\cdot(c\nabla u + \alpha u -\gamma)$. However, it does not specify in full for the higher dimensional cases.

I am currently using the following code

(*Region = ImplicitRegion[(x | y | z) \[Element] Reals && x^2 + y^2 + z^2 <= 1, {x, y, z}];*)
Region = Prism[{{0, Sqrt[3]/4, -(1/2)}, {-(1/2), -(Sqrt[3]/4), -(1/2)}, {1/2, -(Sqrt[3]/4), -(1/2)}, {0, Sqrt[3]/4, 1/2}, {-(1/2), -(Sqrt[3]/4), 1/2}, {1/2, -(Sqrt[3]/4), 1/2}}];
u[x_, y_, z_] := {ux[x, y, z], uy[x, y, z], uz[x, y, z]};

CouplingTensor = ( {{q^2 D[ux[x, y, z], x, x], q D[uy[x, y, z], x, y], D[uz[x, y, z], x, z]},
{q D[ux[x, y, z], x, y], q^2 D[uy[x, y, z], y, y], D[uz[x, y, z], y, z]},
{D[ux[x, y, z], x, z], D[uy[x, y, z], y, z], D[uz[x, y, z], z, z]}} );

NavierEquation = Laplacian[u[x, y, z], {x, y, z}] + (p^2 - 1) CouplingTensor.{1, 1, 1};

SPDE = NavierEquation /. {p -> 18146/12354, q -> 0.5};
{vals, funcs} = NDEigensystem[SPDE + Table[NeumannValue[0, True], {i, 1, 3}],
    {ux[x, y, z], uy[x, y, z], uz[x, y, z]},
    {x, y, z} \[Element] Region, EigenNO];

where SPDE is the PDE $\nabla^2 \vec u + (p^2 - 1) \nabla(\nabla\cdot \vec u)$, and Region is an arbitrary volume. I am not sure if this is employing the correct boundary conditions. How do I enforce the full boundary condition?

  • 1
    $\begingroup$ Can you edit your question to contain the complete code? $\endgroup$ – user21 Nov 25 '15 at 21:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.