3
$\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$
  • 1
    $\begingroup$ So, where is the question? $\endgroup$ – Henrik Schumacher Jan 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$ – Please Correct GrammarMistakes Jan 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 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 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$ – Please Correct GrammarMistakes Jan 22 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

| improve this answer | |
$\endgroup$
  • $\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$ – Please Correct GrammarMistakes Jan 25 at 4:15
  • $\begingroup$ There is another disadvantage for this method, which can't deal with the boundary condition with shear stress. $\endgroup$ – Please Correct GrammarMistakes Jan 25 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$ – Alex Trounev Jan 25 at 13:13
  • 1
    $\begingroup$ @PleaseCorrectGrammarMistakes See update to my answer. $\endgroup$ – Alex Trounev Jan 25 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$ – Alex Trounev Jan 26 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.