How to explain the meaning of `NeumannValue` function in detail

Question

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).

Additional information:

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.

Answer 1

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))), 
         0}, {0, -(Y/(1 - ν^2))}}.Inactive[Grad][
        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_] := 
 TraditionalForm[
  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.