How to explain the meaning of NeumannValue function in detail

I saw here that NeumannValue is used to represent the stress boundary condition when solving the plane stress problem.

In the help information of the NeumannValue function, we can see some mathematical explanations:

Locations where Neumann values might be specified are shown in green. They appear on the boundary [PartialD][CapitalOmega] of the region [CapitalOmega] and specify a flux across those edges in the direction of the outward normal.

I want to know how NeumannValue[1000, x == 1] specifically represents the stress boundary $$\sigma_x=1000$$. I want to get a simple explanation of the mathematical principles.

And I want to know the specific calculation details of \[Del].(-c1 \[Del]u[x, y] - \[Alpha]1 u + \[Gamma]1 - c2 \[Del]v[x, y] - \[Alpha]2 v + \[Gamma]2), how is it equal to $$\sigma_{x}=\frac{\mathrm{Y}}{1-v^{2}}\left(\frac{\partial \mathrm{u}(\mathrm{x}, \mathrm{y})}{\partial \mathrm{x}}+v \frac{\partial \mathrm{v}(\mathrm{x}, \mathrm{y})}{\partial \mathrm{y}}\right)$$.

So I want to know how the formula $$-c \nabla u - \alpha u + \gamma$$ is equivalent to the stress $$\sigma$$. After all, this formula only has displacement functions u and v, but it does not include Poisson's ratio and elastic modulus (This is the core point of this question).

Using displacement functions to express stress:

$$\begin{array}{l} \sigma_{x}=\frac{\mathrm{Y}}{1-v^{2}}\left(\frac{\partial \mathrm{u}(\mathrm{x}, \mathrm{y})}{\partial \mathrm{x}}+v \frac{\partial \mathrm{v}(\mathrm{x}, \mathrm{y})}{\partial \mathrm{y}}\right) \\ \sigma_{\mathrm{y}}=\frac{\mathrm{Y}}{1-v^{2}}\left(\frac{\partial \mathrm{v}(\mathrm{x}, \mathrm{y})}{\partial \mathrm{y}}+v \frac{\partial \mathrm{u}(\mathrm{x}, \mathrm{y})}{\partial \mathrm{x}}\right) \\ \sigma_{\mathrm{xy}}=\frac{(1-v) \mathrm{Y}}{2\left(1-v^{2}\right)} \quad\left(\frac{\partial \mathrm{u}(\mathrm{x}, \mathrm{y})}{\partial \mathrm{y}}+\frac{\partial \mathrm{v}(\mathrm{x}, \mathrm{y})}{\partial \mathrm{x}}\right) \end{array}$$

In the above formula, $$Y$$ represents Young's modulus and $$v$$ represents Poisson's ratio.

$$\begin{array}{l} \frac{\partial \sigma_{x}}{\partial x}+\frac{\partial \sigma_{x y}}{\partial y}=0 \\ \frac{\partial \sigma_{y}}{\partial y}+\frac{\partial \sigma_{x y}}{\partial x}=0 \end{array}$$

$$\begin{array}{l} \frac{Y}{2(1+v)}\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}\right)+\frac{Y}{2(1-v)}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0 \\ \frac{Y}{2(1+v)}\left(\frac{\partial^{2} v}{\partial x^{2}}+\frac{\partial^{2} v}{\partial y^{2}}\right)+\frac{Y}{2(1-v)}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)=0 \end{array}$$

On page 99 of this textbook, there is a formula for the stress tensor expressed by the displacement vector:

But what puzzles me is that no matter how to simplify the formula 3.1, it cannot be consistent with the explanation in the help of MMA's NeumannValue function.

• Did you see the bit in the documentation about $\vec{n}\cdot(c\nabla u+\alpha u-\gamma)=g$ ? For me, the confusing thing about NeumannValue is that it's added to the equation and that it's not in a separate boundary condition like DirichletValue. Commented Jul 6, 2020 at 11:09
• If you look below the table, you'll see that Mathematica automatically sums over the different functions. See the bit that starts with For systems $\nabla\cdot\sum_j\dots$ has a $u_j$ in it. Commented Jul 6, 2020 at 11:14
• I cannot reproduce the mathematical details of the cuff, but the only thing I can really recommend is to try and read up on the fundamentals of FEM. The Neumann conditions kinda drop out naturally once you understand how FEM is used to solve PDEs. Honestly, I don't think there's a real shortcut if you want a proper understanding of what's going on. The mathematics isn't all that daunting if you're not afraid of a few integrals. It's also worth reading up about the "weak formulation of a PDE", because that's closely related. Commented Jul 6, 2020 at 11:48
• @flinty, please have a look at the documentation: See detials of NeumannValue, SolvingPDEwithFEM tutorial and the Finite Element Method Usage Tips tutorial Commented Jul 6, 2020 at 12:09
• @flinty I agree 1000 times about your statement "the confusing thing about NeumannValue is that it's added to the equation". I am sure there must be some logic somewhere for this, but as a philistine its quite frustrating. Commented Jul 9, 2020 at 14:43

I think the best way to think about the NeumannValue is to consider to fundamental property balance equation over the domain at equilibrium. In the case of the plane stress operator from Mathematica's Documentation, I will show that we can derive it from a balance of the traction vector over the boundary of the domain. Therefore, the NeumannValue is simply the traction vector on the boundary.

Note on Coefficient Form

The power of the Finite Element Method is its ability to model wide variety of physical phenomena. The system of Partial Differential Equations (PDE) that describe these phenomena come from balance equations of $$fluxes[ = ]\frac{{property}}{{Area \cdot time}}$$ across surfaces of fundamental properties, such as Mass, Momentum, and Energy, over a differential region. NeumannValues are fluxes. When possible, it best to express your PDE in coefficient form as described in the documentation. The Left Hand Side (LHS) contains the "operator" and the Right Hand Side (RHS) is always 0.

$$m\frac{{{\partial ^2}}}{{\partial {t^2}}}u + d\frac{\partial }{{\partial t}}u + \nabla \cdot\left( { - c\nabla u - \alpha u + \gamma } \right) + \beta \cdot\nabla u + au - f = 0$$

By maintaining the discipline of expressing your PDE system in coefficient form, you will be less likely to make errors in defining your NeumannValues.

Note on Neumann Values

I have used many PDE solvers in my work and one always needs to learn the solver's conventions. In particular, are surface normals, by convention, point into or out of the domain or region. With Mathematica, by convention, a NeumannValue is positive if the flux is into the domain. The other convention is to place the NeumannValues on the RHS of the "equation". I put equation in quotes because it is not really an equation but a convention to bring Neumann conditions into the solver.

Why would one want to do this? Since NeumannValues are fluxes, there can be parallel modes of transport. A classic example is combined convective and radiative heat transfer found in the Heat Transfer Tutorial as shown below.

These parallel modes of heat transfer, can independently, concisely, and clearly be expressed as shown in the documentation as:

pde = {HeatTransferModel[T[x, y], {x, y}, k, ρ, Cp, "NoFlow",
"NoSource"] == Γconvective + Γradiation, Γtemp} /. parameters;
Tfun = NDSolveValue[pde, T, {x, y} ∈ Ω2D]


Once you get used to it, it is a neat and transparent way of expressing NeumannValues. Most other solvers would require that you open up and inspect model elements to deduce the intention.

Derivation of the Plane Stress Operator

First, let's reproduce the plane stress operator from the documentation here:

parmop = {Inactive[
Div][({{0, -((Y ν)/(1 - ν^2))}, {-((Y (1 - ν))/(
2 (1 - ν^2))), 0}}.Inactive[Grad][v[x, y], {x, y}]), {x,
y}] + Inactive[
Div][({{-(Y/(1 - ν^2)),
0}, {0, -((Y (1 - ν))/(2 (1 - ν^2)))}}.Inactive[
Grad][u[x, y], {x, y}]), {x, y}],
Inactive[
Div][({{0, -((Y (1 - ν))/(2 (1 - ν^2)))}, {-((Y ν)/(
1 - ν^2)), 0}}.Inactive[Grad][u[x, y], {x, y}]), {x,
y}] + Inactive[
Div][({{-((Y (1 - ν))/(2 (1 - ν^2))),
v[x, y], {x, y}]), {x, y}]};


At equilibrium and in the absence of body forces, the integral of the traction vector over the boundary should be zero as illustrated in the diagram below. This is the fundamental balance equation.

As shown in the Wiki article for Cauchy stress tensor, we can define the traction vector, $${{\mathbf{T}}^{(\hat n)}}$$, in terms of the unit surface normal, $$\hat {\mathbf{n}}$$, and the stress tensor, $$\mathbf{\sigma}$$:

$${{\mathbf{T}}^{(\hat {\mathbf{n}})}} = \hat {\mathbf{n}} \cdot {\mathbf{\sigma }}$$

In equilibrium and in the absence of body forces, the integral of the traction should be {0,0}.

$$\mathop \smallint \limits_{\partial \Omega } {{\mathbf{T}}^{(\hat {\mathbf{n}})}} \cdot dA = \mathop \smallint \limits_{\partial \Omega } \hat {\mathbf{n}} \cdot {\mathbf{\sigma }}dA = \left( {\begin{array}{*{20}{c}} 0 \\ 0 \end{array}} \right)$$

The Gauss Divergence Theorem also applies to tensors:

$$\mathop \smallint \limits_{\partial \Omega } \hat {\mathbf{n}} \cdot {\mathbf{\sigma }}dA = \mathop \smallint \limits_\Omega ( - \nabla \cdot {\mathbf{\sigma }})dV = \left( {\begin{array}{*{20}{c}} 0 \\ 0 \end{array}} \right) \Rightarrow - \nabla \cdot {\mathbf{\sigma }} = \left( {\begin{array}{*{20}{c}} 0 \\ 0 \end{array}} \right)$$

We will show that $$- \nabla \cdot {\mathbf{\sigma }}$$ is the same as Mathematica's plane stress operator. Since the RHS is zero, we will have expressed our PDE system in coefficient form.

Now, we can grab the definition of strain and stress from Hooke's Law Wiki Article. The infinitessimal strain is defined by:

$${\mathbf{\varepsilon }} = \frac{1}{2}[\nabla {\mathbf{u}} + {(\nabla {\mathbf{u}})^T}]$$

We can relate stress to strain by:

$$\left[ {\begin{array}{*{20}{c}} {{\sigma _{11}}}&{{\sigma _{12}}} \\ {{\sigma _{12}}}&{{\sigma _{22}}} \end{array}} \right]{\mkern 1mu} = {\mkern 1mu} \frac{E}{{1 - {\nu ^2}}}\left( {(1 - \nu )\left[ {\begin{array}{*{20}{c}} {{\varepsilon _{11}}}&{{\varepsilon _{12}}} \\ {{\varepsilon _{12}}}&{{\varepsilon _{22}}} \end{array}} \right] + \nu {\mathbf{I}}\left( {{\varepsilon _{11}} + {\varepsilon _{22}}} \right)} \right)$$

Or

$${\mathbf{\sigma }} = \frac{E}{{1 - {\nu ^2}}}\left( {\left( {1 - \nu } \right){\mathbf{\varepsilon }} + \nu {\mathbf{I}}\operatorname{tr} \left( {\mathbf{\varepsilon }} \right)} \right)$$

In Mathematica code:

ϵ =
1/2 (Grad[{u[x, y], v[x, y]}, {x, y}] +
Transpose@Grad[{u[x, y], v[x, y]}, {x, y}]);
σ = Y/(
1 - ν^2) ((1 - ν) ϵ + ν IdentityMatrix[
2] Tr[ϵ]);
hookeop = -Div[σ, {x, y}];


We can show that our stress, $$\mathbf{\sigma}$$, is equivalent to what the OP expressed (note that $${\nu ^2} - 1 = \left( {\nu + 1} \right)\left( {\nu - 1} \right)$$).

pdConv[f_] :=
f /. Derivative[inds__][g_][vars__] :>
Apply[Defer[D[g[vars], ##]] &,
Transpose[{{vars}, {inds}}] /. {{var_, 0} :>
Sequence[], {var_, 1} :> {var}}]]
σ [[1, 1]] // Simplify // pdConv
σ [[2, 2]] // Simplify // pdConv
σ [[1, 2]] // Simplify // pdConv


Now, let's verify that Mathematica's plane stress operator and our Hooke operator are equal.

hookeop == Activate[parmop] // Simplify
(* True *)


I think this is pretty compelling evidence that we derived Mathematica's plane stress operator correctly.

What is the NeumannValue?

To understand the NeumannValue, we go back to our initial balance equation:

$$\mathop \smallint \limits_{\partial \Omega } {{\mathbf{T}}^{(\hat {\mathbf{n}})}} \cdot dA = \mathop \smallint \limits_{\partial \Omega } \hat {\mathbf{n}} \cdot {\mathbf{\sigma }}dA = \left( {\begin{array}{*{20}{c}} 0 \\ 0 \end{array}} \right)$$

We can either think of the NeumannValue as the traction, $${{\mathbf{T}}^{(\hat {\mathbf{n}})}}$$ , on a boundary or as the surface normal dotted with stress tensor, $$\hat {\mathbf{n}} \cdot {\mathbf{\sigma }}$$. In the OP case of NeumannValue[1000, x == 1], we need to look at both the $$x$$ and $$y$$ components. In terms of stress, to represent tensile stress in the $$x$$-direction, we could write the equation as:

$$\left[ {\begin{array}{*{20}{c}} 1&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\sigma _{11}}}&{{\sigma _{12}}} \\ {{\sigma _{12}}}&{{\sigma _{22}}} \end{array}} \right]{\mkern 1mu} = \left[ {\begin{array}{*{20}{c}} {{\sigma _{11}}}&{{\sigma _{12}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\sigma _{11}}}&0 \end{array}} \right]$$

So, {NeumannValue[1000, x==1], 0} represents a tensile stress of magnitude 1000 in the $$x$$ direction.

One generalize the approach of "flux balance" to other areas, such as heat transfer, to obtain a similar understanding of the NeumannValue.

• Thanks Tim, very informative. Can you also consider shear stresses. Their direction is not normal to the boundary. I wrote this question when trying to understand the direction convention for them. You could answer that old question or deal with it here.
– Hugh
Commented Jul 10, 2020 at 11:39
• @Hugh You are welcome. For a pure shear case in 2d , it should look something like $\left[ {\begin{array}{*{20}{c}} 1&0 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 0&{{\sigma _{12}}} \\ {{\sigma _{12}}}&0 \end{array}} \right]{\mkern 1mu} = \left[ {\begin{array}{*{20}{c}} 0&{{\sigma _{12}}} \end{array}} \right]$I will try to take a look your question over the next couple of days and see if I can help. Commented Jul 10, 2020 at 16:17
• @Hugh I took a stab at answering your other question. Let me know if you have any questions. Commented Jul 17, 2020 at 3:02
• I have been looking at your answer over the past few days. The idea that we should consider traction rather than stress is a point that I need to fully understand. It does sound correct. Further, it is expressed in the global coordinate system rather than in a local coordinate system based on outward normal. The whole outward normal business seems a false trail when looking at stress. I will modify my question when I find the fault in my thinking. Thanks for your answer.
– Hugh
Commented Jul 17, 2020 at 11:31