# DirichletCondition in NDSolve causes asymmetry (Finite Element Method)

Edit

Thanks for the responses. I now agree with xzczd and Bill Watts who kindly gave answers and user21 who wrote useful comments. Both my cases below are correct. On reflection I think I made an elementary mistake. I noted the need to give a Dirichlet Condition for the u-direction to provide a reference for displacements. I should have just tied down one point using the condition

 DirichletCondition[u[x, y] == 0, x == 0 && y == 0]


DirichletCondition[u[x, y] == 0, x == 0]


This latter condition fixes the entire edge y == 0. Fixing this edge introduces asymmetry. I apologise if I have caused pointless concerns. I have also added the single point DirichletCondition as Case 3 below. This gives the same result as Case 1 but fixes a starting point. Thanks to all that helped.

Original Question

If a DirichletCondition is not specified in NDSolve then a warning is given. However, if I specify a condition then I seem to get the wrong answer. Below I have Case 1 where I get the warning and then Case 2 where I don't but I think the answer is wrong. Here is a minimum working example involving a stress computation. We need the following code

Needs["NDSolveFEM"]
ClearAll[planeStress];
planeStress[
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))}, {-((Y*(1 - ν))/(2*(1 \
- ν^2))), 0}} . Inactive[Grad][v[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}]}


Case 1. I make a rectangular grid and apply a force to part of the top surface. The bottom is prevented from moving vertically and is free horizontally. There are no other horizontal forces or constraints.

L = 1; (* Length *)
h = 0.2; (* Height *)
Y = 20 10^10;(* Modulus of elasticity *)
ν = 33/100 ;(* Poission ratio *)
stress = 100;
mesh = ToElementMesh[Rectangle[{0, 0}, {L, h}]];
mesh["Wireframe"]


Now the solver

{uif, vif} = NDSolveValue[{
planeStress[Y, ν] == {0,
NeumannValue[-stress, L/4 <=  x <= 3 L/4 && y == h]},

DirichletCondition[v[x, y] == 0, 0 <=  x <=  L && y == 0]
},
{u, v}, Element[{x, y}, mesh]];


I get the message

NDSolveValue::femibcnd: No DirichletCondition or Robin-type NeumannValue was specified for {u}; the result may not be unique.

The results are good if I Plotthe horizontal displacement along the bottom edge I get

Plot[uif[x, 0], {x, 0, L}]


This has an arbitrary starting point but appears good and retains the symmetry of the problem. To test further I check the derivative (which I need to calculate the stress) and see if the left and right side are similar

du = Head@D[uif[x, y], x];
Plot[du[x, 0], {x, 0, L}]
Plot[du[x, 0] - du[L - x, 0], {x, 0, L/2}]


So far so good. I have a symmetric solution which is very good with the left and right sides differing only by numerical noise. However, I do have an arbitrary starting point for the horizontal displacement.

Case2. Now I try and fix the starting point for horizontal displacements by giving an additional DirichletCondition at one end of the mesh.

{uif, vif} = NDSolveValue[{
planeStress[Y, ν] == {0,
NeumannValue[-stress, L/4 <=  x <= 3 L/4 && y == h]},

DirichletCondition[v[x, y] == 0, 0 <=  x <=  L && y == 0],
DirichletCondition[u[x, y] == 0, x == 0]
},
{u, v}, Element[{x, y}, mesh]];


This is where everything begins to go wrong. Plotting the displacement along the bottom edge gives

Plot[uif[x, 0], {x, 0, L}]


This starts from zero. Checking the derivative again I see problems

du = Head@D[uif[x, y], x];
Plot[du[x, 0], {x, 0, L}]
Plot[du[x, 0] - du[L - x, 0], {x, 0, L/2}]


Now there is a clear distortion and the left and right hand sides are not symmetric. What has gone wrong? Is there a workaround so that I can fix the location of the start of the displacements? Thanks

Edit

In the comments user21 suggests that in the second case the Dirichlet condition is applied to just one side so that this is an asymmetric case automatically. However, on physical grounds I think the symmetry should still be there particularly in the derivatives (the stresses).

Only vertical forces are applied. There are horizontal displacements due to the Poisson Ratio effect. This will cause the horizontal length of the block to increase. When considering stresses this effect will be symmetric. All this is confirmed in Case 1. When we come to Case 2 then when the Poisson Ratio expansion occurs the expansion to the left is prevented by the stop there. However, there are no other horizontal constraints so the block should push against the stop and all the expansion should be to the right. Thus in Case 2 the displacements should be the same as in Case 1 except for everything moving to the right. The derivatives on the left and right (the stresses) should just be the same. This is an argument from the physical stress analysis viewpoint. From the numerical viewpoint I think a constant should be added to the displacements. I suspect a gradient is being added to the displacements. I am surprised that in the first case there is a solution at all since there is nothing to tie down the system. I hope that adds some more background.

Edit 2 User xzczd has asked for the boundary condition in math notation. What I want is to define the stress on the upper surface. Thus

$$100 = \sigma_{yy} = \frac{E}{1 - \nu^2} (\frac{\partial v(x, y)}{\partial y} + \nu \frac{\partial u(x, y)}{\partial x})$$

for $1/4 < x < 3/4$ and $y = 0.2$.

Here $\sigma_{yy}$ is the stress in the vertical direction. I have the stress going downwards so this may be -σyy depending on the sign convention for the boundary. Also modulus of elasticity $E$ and Poisson Ratio $ν$ are E = 20 *10^10 and ν =0.33

** Edit and Case 3**

One point at x = 0 and y = 0 fixed.

{uif, vif} = NDSolveValue[{
planeStress[Y, ν] == {0,
NeumannValue[-stress, L/4 <= x <= 3 L/4 && y == h]},
DirichletCondition[v[x, y] == 0, 0 <= x <= L && y == 0],
DirichletCondition[u[x, y] == 0, x == 0 && y == 0]
},
{u, v},
Element[{x, y}, mesh]
];


The same checks as above

Plot[uif[x, 0], {x, 0, L}]


du = Head@D[uif[x, y], x];
Plot[du[x, 0], {x, 0, L}]
Plot[du[x, 0] - du[L - x, 0], {x, 0, L/2}]


No warning message and a symmetric result as expected.

• Why do you think this should be symmetric? You have a single boundary condition orthogonal to the force. The other side x==L is free. – user21 Mar 27 '18 at 13:16
• @user21 I have no horizontal forces only vertical forces. Therefore I think the datum for horizontal displacements is arbitrary. As in the first case. – Hugh Mar 27 '18 at 13:20
• @xzczd, I am not sure which zero valued NeumannValue you refer to. In any case zero valued NeumannValues are eliminated at parser level. The never go into any equation. Also careful with Activate in these equations as sub matrices A12, etc have asymmetric off diagonal components which are hard maintain in active form (see: Formal Partial Differential Equations) I have not checked is your usage of Activate is affecting the equation in that way. – user21 Mar 28 '18 at 5:35
• @xzczd, it is important to realize that NeumannValue[0, whatever] does not contribute anything ever. It is taken out at parser level. Now, let's assume you have NeumannValue[something, whatever] and DirichletCondition[u==someting, whatever] then the DirichletCondition will trump the NeumannValue. For coupled systems the u DirichletCondition will overwrite the first NeumannValue that is part of the the first equation and the a v DirichletCondition will overwrite the second NeumannValue which is part of the second equation. – user21 Mar 28 '18 at 6:40
• @user21 I think we are debating the correct issue here. When is it correct for the NeumannValue to be trumped by the DirichletCondition? Possibly we need a simpler case to investigate this further. I will try and come up with one. It seems @ xzczd is able to make choices which perhaps is the way to go. – Hugh Mar 28 '18 at 6:48

There seems to be nothing wrong with the result given by NDSolve. The following FDM approach gives the same result. I've used pdetoae for the generation of difference equation:

L = 1;(*Length*)
h = 2/10;(*Height*)
Y = 20 10^10;(*Modulus of elasticity*)
ν = 33/100;(*Poission ratio*)
stress = 100 10^10;

lhsy = planeStress[Y, ν] /. Inactive[Div][a_, __] :> {0, 1}.a /. Inactive -> Identity;
lhsx = planeStress[Y, ν] /. Inactive[Div][a_, __] :> {1, 0}.a /. Inactive -> Identity;
bcxR = lhsx == 0 /. x -> L // Thread;
bcxL2 = lhsx == 0 /. x -> 0 // Thread(*//First*)// Last;
bcyL2 = lhsy == 0 /. y -> 0 // Thread // First;
bcyR = lhsy == {0, Piecewise[{{-stress, L/4 <= x <= (3 L)/4}}]} /. y -> h // Thread //
SimplifyPWToUnitStep;
{bcyL1, bcxL1} = With[{u = u[x, y], v = v[x, y]},
{v == 0 /. y -> 0, u == 0 /. x -> 0}];
{domain@x, domain@y} = {{0, L}, {0, h}};
points = 50;
{grid@x, grid@y} = Array[# &, points, domain@#] & /@ {x, y};
difforder = 2;
(* Definition of pdetoae isn't included in this post,
ptoafunc = pdetoae[{u, v}[x, y], grid /@ {x, y}, difforder];
del = #[[2 ;; -2]] &;
ae = del /@ del@# & /@
ptoafunc[planeStress[Y, ν] == 0 /. Inactive -> Identity // Thread];
aebc@x = del /@ ptoafunc[Flatten@{bcxL1, bcxL2, bcxR}];
aebc@y = ptoafunc[{bcyL1, bcyL2, bcyR} // Flatten];
varlst = Outer[#[#2, #3] &, {u, v}, grid@x, grid@y] // Flatten;
{b, m} = CoefficientArrays[Flatten@{ae, aebc /@ {x, y}}, varlst];
sollst = LinearSolve[m, -N@b];
{func@u, func@v} =
ListInterpolation[#, domain /@ {x, y}] & /@ ArrayReshape[sollst, {2, points, points}];
Plot[-du[x, 0], {x, 0, L}]
Plot[-du[x, 0] + du[L - x, 0], {x, 0, L/2}]


As mentioned in the comment above, it's not quite clear to me which Neumann b.c. should be chosen at $$x=0$$ and $$y=0$$, so I just found out the combination that produces the same result as FEM approach by trial and error. This is also the most… Er… good-looking result in my view. You can check other results yourself by modifying First to Last and Last to First in bcxL2 and bcyL2.

• Very many thanks for this. I think you are asking the correct question when you debate which Neumann b.c should be chosen. I will look at this in more detail. Can you also make my Case 1 above (which I think is the correct answer)? – Hugh Mar 28 '18 at 6:40
• @hugh You just need to take away //Last from bcxL2 and remove all the bcxL1 from the code. – xzczd Mar 28 '18 at 6:46
• Thanks for your help. I think I made an elementary mistake and the problem is now resolved. See the edit in my question. Sorry for possibly wasting your time but I have learnt much from your analysis. Thanks – Hugh Apr 10 '18 at 8:56
• @hugh Er… can you explain a bit deeper? – xzczd Apr 10 '18 at 9:03
• I constrained the entire side y==0. If you just constrain one point (which I had meant to do) then it gives a reference to the u-direction which was my objective. If you constrain an entire side then although the forces are applied in the v-direction the Poisson effects cause asymmetry. If you constrain just one point then this does not influence the solution. – Hugh Apr 10 '18 at 9:07

I believe both of your solutions are correct which we can show with additional plots. Your original solution clearly looks correct even though MMA gives you a warning. Let's look at you second soln. I have mostly used your code, but I added MaxCellMeasure->0.00001 to your mesh to smooth the plots and I changed your NDSolveValue function to:

sol = NDSolve[{planeStress[Y, \[Nu]] == {0,
NeumannValue[-stress, L/4 <= x <= (3 L)/4 && y == h]},
DirichletCondition[v[x, y] == 0, 0 <= x <= L && y == 0],
DirichletCondition[u[x, y] == 0, x == 0]}, {u[x, y],
v[x, y]}, {x, y} \[Element] mesh];

u[x_, y_] = u[x, y] /. sol[[1]];
v[x_, y_] = v[x, y] /. sol[[1]];


This is your second (assymetric) solution and I match the plots you made, but let's go further. Check your DirichletCondition u=0 at x=0

Plot[u[0, y], {y, 0, h}]


That's close enough to 0. Now v=0 at y=0

Plot[v[x, 0], {x, 0, L}]


Also looks good. Now calculate the stresses so we can check.

eq1 = D[u[x, y], x] == (1/Y)*(sx - \[Nu]*sy);
eq2 = D[v[x, y], y] == (1/Y)*(sy - \[Nu]*sx);

sols = Solve[{eq1, eq2}, {sx, sy}] // Flatten // Simplify;

\[Sigma]x[x_, y_] = sx /. sols;
\[Sigma]y[x_, y_] = sy /. sols;


Look at the vertical stress at y=h where you have specified the stress.

Plot[\[Sigma]y[x, h], {x, 0, L}, PlotRange -> All]


Which I believe matches the vertical stress at the top you want and even appears to be symmetric. Check the stress value at the middle.

\[Sigma]y[L/2, h]
(*-99.9982*)


You will get closer to the value of 100 with a finer mesh. While I am no theoretical mathematician, one thing I have always learned is that a solution that satisfies the differential equation and also satisfies the boundary conditions is the correct solution and this assymetric solution does that, so evidently the constraint u=0 at x=0 while leaving the other end free causes an assymmetric solution that seems to be correct.

The solution that you are after, and the one that makes the most sense to me is to change your DirichletCondition to u=0 at x = L/2. When you do that, even though x = L/2 is not an edge, MMa solves it with no complaints about unspecified bc's and all your values become symmetric like your first solution, except that the plot of u` is shifted.

• Thanks for your help. Sorry for the slow reply I have been away and thinking about the problem. You are correct. I have edited the post to show where my thinking went wrong. With your last paragraph you were onto my mistake. – Hugh Apr 10 '18 at 8:53