2
$\begingroup$

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]

enter image description here

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]

enter image description here

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.

$\endgroup$
9
  • 1
    $\begingroup$ So, where is the question? $\endgroup$ Jan 20 '20 at 9:08
  • $\begingroup$ 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 $\endgroup$ Jan 20 '20 at 9:15
  • $\begingroup$ Have you seen for example this this from the NeumannValue ref. page: NeumannValue #1751664584.? And as Henrik, pointed out, you should add a question? $\endgroup$
    – user21
    Jan 20 '20 at 9:25
  • 1
    $\begingroup$ 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. $\endgroup$
    – user21
    Jan 21 '20 at 6:13
  • 1
    $\begingroup$ 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. $\endgroup$ Jan 22 '20 at 3:35
6
+50
$\begingroup$

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]

Figure 2

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]}

Figure 2

$\endgroup$
12
  • $\begingroup$ 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? $\endgroup$ Jan 25 '20 at 4:15
  • $\begingroup$ There is another disadvantage for this method, which can't deal with the boundary condition with shear stress. $\endgroup$ Jan 25 '20 at 7:23
  • 1
    $\begingroup$ @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? $\endgroup$ Jan 25 '20 at 13:13
  • 1
    $\begingroup$ @PleaseCorrectGrammarMistakes See update to my answer. $\endgroup$ Jan 25 '20 at 19:55
  • 1
    $\begingroup$ @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. $\endgroup$ Jan 26 '20 at 2:53

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.