4
$\begingroup$

I have taken a standard equation from Mathematica help from here for a stress operator. I have also formulated the equivalent equations from engineering texts. When I compare the two equations Using FullSimplify I see they are the same. However, when I use them I get different answers. What is happening? Here are the two equations and the comparison using FullSimplify

ps1 = {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}]};
ps = {Y/(2 (1 + ν)) (D[u[x, y], {x, 2}] + D[u[x, y], {y, 2}]) + 
    Y/(2 (1 - ν)) D[(D[u[x, y], x] + D[v[x, y], y]), x],
   Y/(2 (1 + ν)) (D[v[x, y], {x, 2}] + D[v[x, y], {y, 2}]) + 
    Y/(2 (1 - ν)) D[(D[u[x, y], x] + D[v[x, y], y]), y]};
FullSimplify[Activate[ps1] == -ps]

This gives

*True*

I generate a mesh and then use NDSolveValue and get solutions for both which I plot.

Needs["NDSolve`FEM`"];
Ω = ImplicitRegion[True, {x, y}];
mesh = ToElementMesh[Ω, {{0, 1}, {0, 0.25}}, 
   "MaxCellMeasure" -> 0.001];
{uif1, vif1} = NDSolveValue[{
     ps1 == {0, NeumannValue[1, x == 1]},
     DirichletCondition[u[x, y] == 0, x == 0],
     DirichletCondition[v[x, y] == 0, x == 0]
     } /. {Y -> 10^3, ν -> 33/100},
   {u, v},
   {x, y} ∈ mesh];
{uif, vif} = NDSolveValue[{
     ps == {0, -NeumannValue[1, x == 1]},
     DirichletCondition[u[x, y] == 0, x == 0],
     DirichletCondition[v[x, y] == 0, x == 0]
     } /. {Y -> 10^3, ν -> 33/100},
   {u, v},
   {x, y} ∈ mesh];
Plot3D[{vif1[x, y], vif[x, y]}, {x, y} ∈ mesh, 
 BoxRatios -> {2, 1, 1}, PlotRange -> All]

Mathematica graphics

The two solutions are different why? If I look at the ratio, and ignore diving by zero, we can see that there is a factor of about 1.5 but not a constant value.

   Plot3D[Evaluate[vif1[x, y]/vif[x, y]], {x, y} ∈ mesh, 
     BoxRatios -> {2, 1, 1}, PlotRange -> {All, All, {1, 2}}]

Mathematica graphics

Comparison of the other direction gives a similar ratio. What is happening? Is there some subtle issue over NeumannValues that I am missing? Thanks

EDIT: A bit more

I have also looked at the eigenvalues and vectors. These do not involve the Neumann boundary conditions.

{vals1, vecs1} = 
  NDEigensystem[{ps1, DirichletCondition[u[x, y] == 0, x == 0], 
     DirichletCondition[v[x, y] == 0, x == 0]} /. {Y -> 10^3, \[Nu] ->
       33/100}, {u, v}, {x, y} \[Element] mesh, 10];
{vals, vecs} = 
  NDEigensystem[{-ps, DirichletCondition[u[x, y] == 0, x == 0], 
     DirichletCondition[v[x, y] == 0, x == 0]} /. {Y -> 10^3, \[Nu] ->
       33/100}, {u, v}, {x, y} \[Element] mesh, 10];
TableForm[Transpose[{vals1, vals}], 
 TableHeadings -> {None, {"Help Eqn.", "Textbook Eqn."}}]

The table comparing the eigenvalues gives

Mathematica graphics

where I have called equation ps1 the Help equation and ps the text book equation. Clearly they are very different values. Looking at the first eigenvector I again compare the ratio

{vif1, vif} = {vecs1[[1, 2]], vecs[[1, 2]]};
Plot3D[vif1[x, y]/vif[x, y], {x, y} \[Element] mesh, 
 BoxRatios -> {2, 1, 1}, PlotRange -> {All, All, {0.8, 1.2}}]

Mathematica graphics

These are more similar than the deflection calculation but still significantly different.

I am not sure what to conclude but the stiffness matrices must be different. Is there a good reason for this? Version 10.3 on Windows 7.

Can anyone confirm?

$\endgroup$
2
$\begingroup$

Let's look at a simple example:

pde = Inactive[
   Div][{{0, 1}, {2, 0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}]

the coefficient is

{{0, 1}, {2, 0}}

Now when you activate that you get

Activate[pde]
3*Derivative[1, 1][u][x, y]

This is then what NDSolve parses in the two cases

{state} =
  NDSolve`ProcessEquations[{pde == 0,
    DirichletCondition[u[x, y] == 0, True]},
   u, {x, y} \[Element] Rectangle[]];
state["FiniteElementData"][
  "PDECoefficientData"]["DiffusionCoefficients"]
{{{{0, 1}, {2, 0}}}}

{state} =
  NDSolve`ProcessEquations[{Activate[pde] == 0,
    DirichletCondition[u[x, y] == 0, True]},
   u, {x, y} \[Element] Rectangle[]];
state["FiniteElementData"][
  "PDECoefficientData"]["DiffusionCoefficients"]
{{{{0, 3/2}, {3/2, 0}}}}

So you see this a different PDE model. Inactive is a way to prevent Mathematica to evaluate the PDE too early until there is a chance to parse the PDE before it is evaluated. There is no way to go back from the evaluated coefficients like

{{{{0, 3/2}, {3/2, 0}}}} to {{{{0, 1}, {2, 0}}}}

So, yes, Inactive is needed if either

1) you have unsymmetrical coefficients (which we have for plane stress/strain) That's what we have here or

2) if you need NeumannValue that work with the divergence part of the PDE which is explained here:

FEMDocumentation/tutorial/FiniteElementBestPractice

in the section:

NeumannValue and Formal Partial Differential Equations

More info here and there is a question about this as well. That the inactive stress operator is correct has been discussed here.

I hope this helps.

$\endgroup$
  • $\begingroup$ Thank you for this. I have had one look at your references (i will read in depth) and can see there is a problem in providing the information needed for the Neumann boundary conditions. This is disappointing because I now realise that I can't just "put the equations into Mathematica" and they will work. Setting up equations is going to be deeply difficult. Is there a way we can get around this by working with stress -that would give us Dirichlet conditions? $\endgroup$ – Hugh Dec 10 '15 at 16:33
  • $\begingroup$ I am not sure I understand correctly, could you explain to me how you would want to input the equations? $\endgroup$ – user21 Dec 10 '15 at 17:57
  • $\begingroup$ I would like to enter the relationship between stress and strain and then the equilibrium equations between normal and shear stresses. One could then apply Dirichlet conditions for the stresses. These two sets of equations combined lead to my textbook equation above. If you wish I could put this into another question or as an answer here. $\endgroup$ – Hugh Dec 10 '15 at 19:56
  • $\begingroup$ Maybe it's better to put in a new question. $\endgroup$ – user21 Dec 10 '15 at 21:33
  • $\begingroup$ @Hugh, I just want to make sure that I did not miss a new question? Do you still have a question related to this or did things clear up a bit? $\endgroup$ – user21 Dec 15 '15 at 17:06

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.