# The NDSolvevalue of MMA can not be used to solve the finite element problems according to the stress boundary conditions

The NDSolvevalue of MMA can well solve the finite element problems according to the displacement boundary conditions

(*FEMDocumentation/tutorial/SolvingPDEwithFEM*)

Ω=RegionDifference[Rectangle[{-1,-1},{1,1}],Rectangle[{-1/2,-1/2},{1/2,1/2}]];
op={-Derivative[0, 2][u][x, y] - Derivative[2, 0][u][x,
y] - Derivative[1, 1][v][x, y],
-Derivative[1, 1][u][x, y] - Derivative[0, 2][v][x,
y] - Derivative[2, 0][v][x, y]}

Subscript[Γ, D]={DirichletCondition[{u[x,y]==1.,v[x,y]==0.},x==1/2&&-1/2<=y<=1/2],DirichletCondition[{u[x,y]==-1.,v[x,y]==0.},x==-1/2&&-1/2<=y<=1/2],DirichletCondition[{u[x,y]==0.,v[x,y]==-1.},y==-1/2&&-1/2<=x<=1/2],DirichletCondition[{u[x,y]==0.,v[x,y]==1.},y==1/2&&-1/2<=x<=1/2],DirichletCondition[{u[x,y]==0.,v[x,y]==0.},Abs[x]==1||Abs[y]==1]}
{ufun,vfun}=NDSolveValue[{op=={0,0},Subscript[Γ, D]},{u,v},{x,y}∈Ω,  StartingStepSize->0.1,MaxStepSize->0.01, WorkingPrecision->30,InterpolationOrder->All, NormFunction->(Norm[#, 1]&)]
ContourPlot[ufun[x,y],{x,y}∈Ω,ColorFunction->"Temperature",AspectRatio->Automatic,PlotPoints->30,WorkingPrecision->20,Contours->30]


But the ndsolvevalue of MMA can not be used to solve the finite element problems according to the stress boundary conditions

      Clear["Gloabal*"]
Ω =
RegionDifference[Rectangle[{-1, -1}, {1, 1}],
Rectangle[{-1/2, -1/2}, {1/2, 1/2}]];

op = {D[σx[x, y], x] + D[τxy[x, y], y],
D[τxy[x, y], x] + D[σy[x, y], y],
Laplacian[σx[x, y] + σy[x, y], {x, y}]};
(*∂Subscript[σ,xx](x,y)/∂x+∂\
Subscript[τ,xy](x,y)/∂y\[Equal]0 ∂Subscript[\
σ,yy](x,y)/∂y+∂Subscript[τ,xy](x,y)/\
∂x\[Equal]0*)
Subscript[Γ,
D] = {DirichletCondition[{σx[x, y] == 10., σy[x, y] ==
0., τxy[x, y] == 0.},
Abs[x] == 1/2 && -1/2 <= y <= 1/2 || -1/2 <= x <= 1/2 &&
Abs[y] == 1/2],
DirichletCondition[{σx[x, y] == -10., σy[x, y] ==
0., τxy[x, y] == 0.}, Abs[x] == 1 || Abs[y] == 1]}

(*{ufun,vfun,wfun}=NDSolveValue[{op\[Equal]{0,0,0},Subscript[\
Γ,D]},{σx,σy,τxy},{x,0,5},{y,0,1},\
Method\[Rule]{"PDEDiscretization"\[Rule]{"MethodOfLines",{\
"SpatialDiscretization"\[Rule]"FiniteElement"}}}]*)
{ufun, vfun, wfun} =
NDSolveValue[{op == {0, 0, 0},
Subscript[Γ,
D]}, {σx, σy, τxy}, {x,
y} ∈ Ω, StartingStepSize -> 0.1,
MaxStepSize -> 0.01, WorkingPrecision -> 20]
ContourPlot[ufun[x, y], {x, y} ∈ Ω,
ColorFunction -> "Temperature", AspectRatio -> Automatic]


The result of this image is obviously incorrect.

Supplementary information:

Equilibrium differential equation: $$\frac {\partial \sigma _ {\text {x}}} {\partial x} +\frac {\partial \tau _ {\text {xy}}} {\partial y} =0$$

$$\frac {\partial \tau _ {\text {xy}}} {\partial x}+\frac {\partial \sigma _ {\text {y}}} {\partial y} =0$$ Deformation compatibility equation expressed by stress: $$\left( \frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2} \right) \left( \sigma _{\text{x}}+\sigma _{\text{y}} \right) =0$$

Because $$\frac {\partial \tau _ {\text {xy}}} {\partial y}=-\frac {\partial \sigma _ {\text {x}}} {\partial x}$$ and $$\frac {\partial \tau _ {\text {xy}}} {\partial x}=-\frac {\partial \sigma _ {\text {y}}} {\partial y}$$, we can get $$2\frac{\partial ^2\tau _{\text{xy}}}{\partial x\partial x}=-2\left( \frac{\partial ^2\sigma _{\text{x}}}{\partial x^2}+\frac{\partial ^2\sigma _{\text{y}}}{\partial y^2} \right)$$ Therefore, the deformation compatibility equation expressed by stress ( $$\frac {\partial^{2} (\sigma _ {\text {x}} - \mu \sigma _ {\text {y}})} {\partial y^{2}} + \frac {\partial^{2} (\sigma _ {\text {y}} - \mu \sigma _ {\text {x}})} {\partial x^{2}}=2(1+\mu)\frac {\partial^{2} \tau _ {\text {xy}}} {\partial x \partial y}$$) can be simplified as $$\frac {\partial^{2} \sigma _ {\text {x}}} {\partial x^{2}}+\frac {\partial^{2} \sigma _ {\text {x}}} {\partial y^{2}} +\frac {\partial^{2} \sigma _ {\text {y}}} {\partial x^{2}}+\frac {\partial^{2} \sigma _ {\text {y}}} {\partial y^{2}}=0$$.

It can be abbreviated as $$\left( \frac{\partial ^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2} \right) \left( \sigma _{\text{x}}+\sigma _{\text{y}} \right) =0$$

This is also the expression of op[[3]] before the modification of my code: 2 ∂τxy(x,y)/(∂x∂y)+∂σx(x,y)/∂x^2+∂σy(x,y)/∂y^2

It's a mistake because I'm dizzy.

• So, where is the question? – Henrik Schumacher Jan 20 at 9:08
• Because the stress boundary condition is symmetric, the result of solving by the stress boundary condition should also be symmetric, which is generally similar to the result of solving by the symmetric displacement boundary condition, but the result of solving by the stress boundary condition is obviously asymmetric – Montevideo Jan 20 at 9:15
• Have you seen for example this this from the NeumannValue ref. page: NeumannValue #1751664584.? And as Henrik, pointed out, you should add a question? – user21 Jan 20 at 9:25
• You'd need to give a definition of what a symmetric stress boundary condition is. I suspect that this amortizes to a NeumannValue which you use in conjunction with an Inactive version of your stress operator. Unfortunately, I can not help you unless you explain better what you need. – user21 Jan 21 at 6:13
• This is a plane stress problem in elasticity. Because stress boundary conditions and geometric characteristics are symmetrical(op in code) , the final stress distribution result should be symmetric, but the results obtained by using NDSolve are disordered. – Montevideo Jan 22 at 3:35

Problem number 2. I do not understand why the author changed the system of equations, but for the new system there is also a symmetric solution. Differentiating op[[1]] with respect to x and op[[2]] with respect to y, we solve the resulting system using FEM, we find a solution

Ω =
RegionDifference[Rectangle[{-1, -1}, {1, 1}],
Rectangle[{-1/2, -1/2}, {1/2, 1/2}]];

op1 = {D[σx[x, y], x, x] + D[τxy[x, y], y, x],
D[τxy[x, y], x, y] + D[σy[x, y], y, y],
Laplacian[σx[x, y] + σy[x, y], {x, y}]};

Subscript[Γ,
D] = {DirichletCondition[{σx[x, y] == 10., σy[x, y] ==
0., τxy[x, y] == 0.},
Abs[x] == 1/2 && -1/2 <= y <= 1/2 || -1/2 <= x <= 1/2 &&
Abs[y] == 1/2],
DirichletCondition[{σx[x, y] == -10., σy[x, y] ==
0., τxy[x, y] == 0.}, Abs[x] == 1 || Abs[y] == 1]}

{ufun, vfun, wfun} =
NDSolveValue[{op1 == {0, 0, 0},
Subscript[Γ,
D]}, {σx, σy, τxy}, {x,
y} ∈ Ω,
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}]
DensityPlot[ufun[x, y], {x, y} ∈ Ω,
ColorFunction -> "Rainbow", AspectRatio -> Automatic,
PlotRange -> All, PlotPoints -> 100, PlotLegends -> Automatic]


There is another system that can be deduced from the original:

op2 = {D[σx[x, y], x, x] - D[σy[x, y], y, y],
Laplacian[σx[x, y] + σy[x, y], {x, y}],
Laplacian[τxy[x, y], {x, y}] +
D[D[σx[x, y] + σy[x, y], x], y]};


With boundary conditions

bc2={DirichletCondition[{σx[x, y] == 10., σy[x, y] ==
1., τxy[x, y] == 1.},
Abs[x] == 1/2 && -1/2 <= y <= 1/2 || -1/2 <= x <= 1/2 &&
Abs[y] == 1/2],
DirichletCondition[{σx[x, y] == -10., σy[x, y] ==
0., τxy[x, y] == 0.}, Abs[x] == 1 || Abs[y] == 1]};


We have

{ufun, vfun, wfun} =
NDSolveValue[{op2 == {0, 0, 0},
bc2}, {σx, σy, τxy}, {x,
y} ∈ Ω,
Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}}]
{DensityPlot[ufun[x, y], {x, y} ∈ Ω,
ColorFunction -> "Rainbow", AspectRatio -> Automatic,
PlotRange -> All, PlotPoints -> 100, PlotLegends -> Automatic],
DensityPlot[vfun[x, y], {x, y} ∈ Ω,
ColorFunction -> "Rainbow", AspectRatio -> Automatic,
PlotRange -> All, PlotPoints -> 100, PlotLegends -> Automatic],
DensityPlot[wfun[x, y], {x, y} ∈ Ω,
ColorFunction -> "Rainbow", AspectRatio -> Automatic,
PlotRange -> Automatic, PlotPoints -> 100,
PlotLegends -> Automatic]}
`

• Thank you very much, but your differential equation is not the same as the differential equilibrium equation in elasticity (and lack of shear stress function). Are they equivalent? – Montevideo Jan 25 at 4:15
• There is another disadvantage for this method, which can't deal with the boundary condition with shear stress. – Montevideo Jan 25 at 7:23
• @PleaseCorrectGrammarMistakes I see that you have changed the equations in your question. This does not match your previous one now. Therefore, my answer looks strange. You did it on purpose to puzzle me? – Alex Trounev Jan 25 at 13:13
• @PleaseCorrectGrammarMistakes See update to my answer. – Alex Trounev Jan 25 at 19:55
• @PleaseCorrectGrammarMistakes I answered your question about a symmetric solution for $\sigma_x$. If there is another question, then you need to open a new topic. – Alex Trounev Jan 26 at 2:53