8
$\begingroup$

I want to solve the plane stress problem of the following stress boundary:

enter image description here

Note that $\sigma_x=10$, etc. in the figure above represents the stress boundary condition rather than the displacement condition.

I built a system of equations based on the stress balance equation and the deformation compatibility equation.

Needs["NDSolve`FEM`"]
Ω = 
  RegionDifference[Rectangle[{-1, -1}, {1, 1}], Disk[{0, 0}, 0.1]];
Region[Ω]


op = { D[σx[x, y], x] + D[τxy[x, y], y], 
  D[σy[x, y], y] + D[τxy[x, y], x], 
  2 D[τxy[x, y], x, y] + D[σx[x, y], x, x] + 
   D[σy[x, y], y, 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;*)

Γ = {DirichletCondition[{σx[x, y] == 
      0., σy[x, y] == 0., τxy[x, y] == 0.}, 
    x^2 + y^2 == 0.1^2], 
   DirichletCondition[{σx[x, y] == 10., σy[x, y] == 
      0., τxy[x, y] == 0.}, x == 1 && -1 <= y <= 1], 
   DirichletCondition[{σx[x, y] == -10., σy[x, y] == 
      0., τxy[x, y] == 0.}, x == -1 && -1 <= y <= 1], 
   DirichletCondition[{σx[x, y] == 0., σy[x, y] == 
      10., τxy[x, y] == 0.}, y == 1 && -1 <= x <= 1], 
   DirichletCondition[{σx[x, y] == 
      0., σy[x, y] == -10., τxy[x, y] == 0.}, 
    y == -1 && -1 <= x <= 1]};


{ufun, vfun, wfun} = 
 NDSolveValue[{op == {0, 0, 
     0}, Γ}, {σx, σy, τxy}, {x, 
    y} ∈ Ω,  StartingStepSize -> 0.1, 
  MaxStepSize -> 0.01, WorkingPrecision -> 20]

ContourPlot[ufun[x, y], {x, y} ∈ Ω, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic, 
 PlotPoints -> 30, WorkingPrecision -> 20, Contours -> Range[0, 5, 1],
  PlotRange -> Full]

But the solution result is obviously wrong: enter image description here

How can I use Mathematica to solve this kind of plane stress problem?

There is a similar post here, but I would like to know if there is a general method to solve this type of problem that does not require additional processing skills. Thank you.

$\endgroup$
  • 3
    $\begingroup$ Replace x^2 + y^2 == 0.1 by x^2 + y^2 == 0.1^2. $\endgroup$ – Henrik Schumacher Jun 20 at 7:15
  • $\begingroup$ @HenrikSchumacher Thank you very much for your prompt, I have updated the question. $\endgroup$ – A little mouse on the pampas Jun 20 at 8:33
  • 4
    $\begingroup$ @PleaseCorrectGrammarMistakes Why don't you use plane stress operator? $\endgroup$ – Alex Trounev Jun 20 at 10:50
  • 2
    $\begingroup$ A question about your operator. Is it dimensionally correct? Where did the operator come from? Perhaps search for Beltrami-Mitchell to find relevant material online. $\endgroup$ – PaulCommentary Jun 21 at 2:45
  • 1
    $\begingroup$ @PleaseCorrectGrammarMistakes op[[3]]=op[[1]]+op[[2]], so we have 3 linear dependent equations op=={0,0,0}. It is the reason to take 2 functions u[x,y], v[x,y] to solve the system op[[1]]==1, op[[2]]==0. $\endgroup$ – Alex Trounev Jul 3 at 10:06
10
$\begingroup$

A more suitable solution can be found in the Mathematica documentation for solving plane stress in the section structural mechanics for NDEigensystem.

Specify a plane stress PDE:

plane stress pde

{vals, funs} = 

NDEigensystem[{ps, DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0]}, {u[x, y], v[x, y]}, {x, y} ∈ Ω, 9];

vals

{351.293, 369.64, 495.516, 1479.33, 2021.45, 2113.61, 2171.36,
2451.13, 3434.16}

Show[{Graphics3D[
     {Gray, 
      GraphicsComplex[{{-1, -1, 0}, {1, -1, 0}, {1, 1, 0}, {1, -1, 
         0}}, Line[{{1, 2, 3, 4, 1}}]]}],
    Plot3D[Sqrt[Total[#^2]], {x, y} ∈ Ω, 
     ColorFunction -> "TemperatureMap", Axes -> False, 
     Mesh -> False]}, Boxed -> False] & /@ funs

Plot3D of the Eigenfunctions of the problem

The realized solution is now an arbitrary linear combination of the Eigenfunctions combined to solve the boundary conditions.

Mind I have selected material with material properties from the Mathematica example.

From FiniteElementProgramming section Coupled PDEs comes right in the example Deformation of a Beam under Load

Clear[u, v, x, y]
op = {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}]};

mesh["Wireframe"]

plot of the mesh

The following are all step from the example that are already abstracted for use in varied cases:

Subscript[Γ, 
  u] = {NeumannValue[{u[x, y] == 0.}, x^2 + y^2 == 0.1^2], 
   NeumannValue[{u[x, y] == 10.}, x == 1 && -1 <= y <= 1], 
   NeumannValue[{u[x, y] == -10.}, x == -1 && -1 <= y <= 1], 
   NeumannValue[{u[x, y] == 0.}, y == 1 && -1 <= x <= 1], 
   NeumannValue[{u[x, y] == 0.}, y == -1 && -1 <= x <= 1]};
Subscript[Γ, 
  v] = {NeumannValue[{v[x, y] == 0.}, x^2 + y^2 == 0.1^2], 
   NeumannValue[{v[x, y] == 0.}, x == 1 && -1 <= y <= 1], 
   NeumannValue[{v[x, y] == 0.}, x == -1 && -1 <= y <= 1], 
   NeumannValue[{v[x, y] == 10.}, y == 1 && -1 <= x <= 1], 
   NeumannValue[{v[x, y] == -10.}, y == -1 && -1 <= x <= 1]};

vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u, v}, {x, y}}];

sd = NDSolve`SolutionData["Space" -> ToNumericalRegion[mesh]];

methodData = InitializePDEMethodData[vd, sd]

Length[mesh["Coordinates"]]*
 Length[NDSolve`SolutionDataComponent[vd, "DependentVariables"]]
methodData["DegreesOfFreedom"]

720

formula 12

diffusionCoefficients = 
  "DiffusionCoefficients" -> {{{{-(Y/(1 - ν^2)), 
        0}, {0, -((Y (1 - ν))/(2 (1 - ν^2)))}}, {{0, -((
         Y ν)/(1 - ν^2))}, {-((Y (1 - ν))/(
         2 (1 - ν^2))), 
        0}}}, {{{0, -((Y (1 - ν))/(2 (1 - ν^2)))}, {-((
         Y ν)/(1 - ν^2)), 
        0}}, {{-((Y (1 - ν))/(2 (1 - ν^2))), 
        0}, {0, -(Y/(1 - ν^2))}}}} /. {Y -> 10^3, ν -> 
     33/100};

initCoeffs = 
 InitializePDECoefficients[vd, sd, {diffusionCoefficients}]

initBCs = 
 InitializeBoundaryConditions[vd, 
  sd, {Subscript[Γ, u], Subscript[Γ, v]}]

discretePDE = DiscretizePDE[initCoeffs, methodData, sd]

{load, stiffness, damping, mass} = discretePDE["SystemMatrices"]

MatrixPlot[stiffness]

stiffness matrix matrixplot

split = Span @@@ 
  Transpose[{Most[# + 1], Rest[#]} &[methodData["IncidentOffsets"]]]

{1 ;; 360, 361 ;; 720}

discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd]

DeployBoundaryConditions[{load, stiffness}, discreteBCs]

And now the time-consuming step. I do not have enough time to verify the boundary conditions in depths. May by my transfer from the given ones are not too suitable.

Short[solution = LinearSolve[stiffness, load]]

ufun = ElementMeshInterpolation[{mesh}, solution[[split[[1]]]]]
vfun = ElementMeshInterpolation[{mesh}, solution[[split[[2]]]]]

ContourPlot[ufun[x, y], {x, y} ∈ mesh, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic]
ContourPlot[vfun[x, y], {x, y} ∈ mesh, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic]

dmesh = ElementMeshDeformation[mesh, {ufun, vfun}]

Show[{
  mesh["Wireframe"],
  dmesh["Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

Since after the material selection only the region, the boundary conditions have to be formulated properly it is not very much effort left after understanding what is done in the given abstracted steps from Wolfram Inc.. Vary the MaxCellMeasure value.

Excuse for the inconvenience. It seems to be an error in Mathematica 12.0, corrected in 12.1.

A workaround is presented in how-do-i-use-low-level-fem.

\[CapitalOmega] = 
 ImplicitRegion[-1 <= x <= 1 && -1 <= y <= 1 && 
   Sqrt[x^2 + y^2] >= 0.1, {x, y}]

RegionPlot[\[CapitalOmega], PlotRange -> {{-1.1, 1.1}, {-1.1, 1.1}}]

Regionplot

Needs["NDSolve`FEM`"]
{state} = 
  NDSolve`ProcessEquations[{Laplacian[u[x, y], {x, y}] == 1, 
    DirichletCondition[u[x, y] == 0, True]}, 
   u, {x, y} \[Element] \[CapitalOmega], Method -> {"FiniteElement"}];

femdata = state["FiniteElementData"]
femdata["Properties"]
methodData = femdata["FEMMethodData"];
bcData = femdata["BoundaryConditionData"];
pdeData = femdata["PDECoefficientData"];
variableData = state["VariableData"];
solutionData = state["SolutionData"][[1]];

(FiniteElementData["<" 1387 ">"]) ({"BoundaryConditionData", "FEMMethodData", "PDECoefficientData",
"Properties", "Solution"}
)

pdeData["All"]

({{{{1}}, {{{{0}, {0}}}}}, {{{{{-1, 0}, {0, -1}}}}, {{{{0}, {0}}}}, {{{{0, 0}}}}, {{0}}}, {{{0}}}, {{{0}}}})

discretePDE = DiscretizePDE[pdeData, methodData, solutionData]
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"]

(DiscretizedPDEData["<" !(* TagBox[ TooltipBox["1387", ""Total degrees of freedom"", TooltipStyle->"TextStyling"], Annotation[#, "Total degrees of freedom", "Tooltip"]& ]) ">"])

SparseArray

MatrixPlot[stiffness]

Matrixplot of the stiffness matrix

discreteBCs = 
  DiscretizeBoundaryConditions[bcData, methodData, solutionData];

DeployBoundaryConditions[{load, stiffness}, discreteBCs]

solution = LinearSolve[stiffness, load];

[![mesh = methodData\["ElementMesh"\];
ifun = ElementMeshInterpolation\[{mesh}, solution\]][12]][12]

Plot3D of the distorted plate:

Plot3D[ifun[x, y], {x, y} \[Element] mesh]

Plot3D of the distorted plate

This looks pretty much like the solution without the hole in the middle superposed with the distortion caused by the fixed whole.

Another solution is

r = ImplicitRegion[-1 <= x <= 1 && -1 <= y <= 1 && 
   Sqrt[x^2 + y^2] >= 0.1, {{x, -2, 2}, {y, -2, 2}}]

op = {Inactive[
       Div][{{0, -((nu*Y)/(1 - nu^2))}, {-((1 - nu)*Y)/(2*(1 - nu^2)),
          0}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] + 
     Inactive[
       Div][{{-(Y/(1 - nu^2)), 
         0}, {0, -((1 - nu)*Y)/(2*(1 - nu^2))}}.Inactive[Grad][
        u[x, y], {x, y}], {x, y}], 
    Inactive[
       Div][{{0, -((1 - nu)*Y)/(2*(1 - nu^2))}, {-((nu*Y)/(1 - nu^2)),
          0}}.Inactive[Grad][u[x, y], {x, y}], {x, y}] + 
     Inactive[
       Div][{{-((1 - nu)*Y)/(2*(1 - nu^2)), 
         0}, {0, -(Y/(1 - nu^2))}}.Inactive[Grad][
        v[x, y], {x, y}], {x, y}]} /. {Y -> 10^3, nu -> 33/100};

Subscript[\[CapitalGamma], D] = 
  DirichletCondition[{u[x, y] == 0., 
    v[x, y] == 
     0.}, (x == -1 && y == -1) || (x == -1 && y == 1) || (x == 1 && 
      y == -1) || (x == 1 && y == 1)];

force = -40; (*stress is 20, surface area is 2*)

{ufun, vfun} = 
  NDSolveValue[{op == {NeumannValue[force, x == 1 || x == -1], 
      NeumannValue[-force, y == -1 || y == 1]}, 
    Subscript[\[CapitalGamma], D]}, {u, v}, {x, y} \[Element] r];

Deformation in the plane:

mesh = ufun["ElementMesh"];
Show[{mesh["Wireframe"["MeshElement" -> "BoundaryElements"]], 
  NDSolve`FEM`ElementMeshDeformation[mesh, {ufun, vfun}][
   "Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

mesh deformation plot

The first example solves with the NeumannValues set and the DirichletValues implicit. This one uses both explicit. This shows both stresses in the same direction and therefore inward and outward combined. This time the center hole moves with the deformed plate and the force is somehow appearing not so super uniform but incremental and therefore maximal in the middle of the sides. All four corners remain fixed in the response. The hole too is not deformed.

This collects the necessary questions that have to be answered to give great solution. This kind of problem belongs most often to the class of complete problems. Despite DirichletValue and NeumannValue given are other stiffnesses needed to be exact in an overall defined problem.

I presented several examples from the Mathematica documentation. Not each is great and matches the question or performs straight through.

Subscript[\[CapitalGamma], D] = 
  DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 
   Sqrt[x^2 + y^2] <= 0.1];

{ufun, vfun} = 
  NDSolveValue[{op == {NeumannValue[force, x == 1 || x == -1], 
      NeumannValue[-force, y == -1 || y == 1]}, 
    Subscript[\[CapitalGamma], D]}, {u, v}, {x, y} \[Element] r];

mesh = ufun["ElementMesh"];
Show[{mesh["Wireframe"["MeshElement" -> "BoundaryElements"]], 
  NDSolve`FEM`ElementMeshDeformation[mesh, {ufun, vfun}][
   "Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

deformed mesh plot

Subscript[\[CapitalGamma], D] = 
  DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 
   Sqrt[x^2 + y^2] <= 
     0.1 || (x == -1 && x == 1 && y == -1 && y == 1)];

gives no difference to the former definition of the DirichletValue.

Subscript[\[CapitalGamma], D] = 
  DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, 
   Sqrt[x^2 + y^2] <= 0.1];

{ufun, vfun} = 
  NDSolveValue[{op == {NeumannValue[Sign[x]*force, x == 1 || x == -1],
       NeumannValue[-Sign[y]*force, y == -1 || y == 1]}, 
    Subscript[\[CapitalGamma], D]}, {u, v}, {x, y} \[Element] r];

mesh = ufun["ElementMesh"];
Show[{mesh["Wireframe"["MeshElement" -> "BoundaryElements"]], 
  NDSolve`FEM`ElementMeshDeformation[mesh, {ufun, vfun}][
   "Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

Meshplot of the deformed and undeformed mesh

or turned by 90 degree to matched the sketch given. Or the minus of the force exerted changed in x and y.

As the path of my presentation went there is much more calculated that displacement by NDSolveValue in there runs and can be displayed.

Show[{Graphics[{Blue, Arrow[{{0, 0}, {0.1, 0}}]}], gr, 
  Graphics[Table[{Red, Arrow[{{k/6, -1.3}, {k/6, -1.1}}]}, {k, -6, 6, 
     2}]], Graphics[
   Table[{Red, Arrow[{{k/6, 1.3}, {k/6, 1.1}}]}, {k, -6, 6, 2}]], 
  Graphics[Table[{Red, Arrow[{{-1.1, k/6}, {-1.3, k/6}}]}, {k, -6, 6, 
     2}]], Graphics[
   Table[{Red, Arrow[{{1.1, k/6}, {1.3, k/6}}]}, {k, -6, 6, 2}]], 
  Graphics[{Red, Inset[Subscript[\[Sigma], x] == 10, {1.3, 0.1}], 
    Inset[Subscript[\[Sigma], x] == 10, {-1.3, 0.1}], 
    Inset[Subscript[\[Sigma], y] == -10, {0.15, 1.35}], 
    Inset[Subscript[\[Sigma], y] == -10, {0.15, -1.4}], Blue, 
    Inset[Text[Style["b", FontSize -> 24]], {0.15, 0.9}], 
    Inset[Text[Style["b", FontSize -> 24]], {0.15, -0.9}], 
    Inset[Text[Style["a", FontSize -> 24]], {-0.9, 0.15}], 
    Inset[Text[Style["a", FontSize -> 24]], {0.9, 0.15}], 
    Inset[Text[Style["r", FontSize -> 12]], {0., -0.0625}]}]}]

result mesh distorted and undistorted

-sigma_y must be up or down or change sign, see my change for force. Same for sigma_x.

Corrected version:

Show[{Graphics[{Blue, Arrow[{{0, 0}, {0.1, 0}}]}], gr, 
  Graphics[Table[{Red, Arrow[{{k/6, -1.3}, {k/6, -1.1}}]}, {k, -6, 6, 
     2}]], Graphics[
   Table[{Red, Arrow[{{k/6, 1.3}, {k/6, 1.1}}]}, {k, -6, 6, 2}]], 
  Graphics[Table[{Red, Arrow[{{-1.1, k/6}, {-1.3, k/6}}]}, {k, -6, 6, 
     2}]], Graphics[
   Table[{Red, Arrow[{{1.1, k/6}, {1.3, k/6}}]}, {k, -6, 6, 2}]], 
  Graphics[{Red, Inset[Subscript[\[Sigma], x] == 10, {1.3, 0.1}], 
    Inset[Subscript[\[Sigma], x] == -10, {-1.3, 0.1}], 
    Inset[Subscript[\[Sigma], y] == -10, {0.15, 1.35}], 
    Inset[Subscript[\[Sigma], y] == 10, {0.15, -1.4}], Blue, 
    Inset[Text[Style["b", FontSize -> 24]], {0.15, 0.9}], 
    Inset[Text[Style["b", FontSize -> 24]], {0.15, -0.9}], 
    Inset[Text[Style["a", FontSize -> 24]], {-0.9, 0.15}], 
    Inset[Text[Style["a", FontSize -> 24]], {0.9, 0.15}], 
    Inset[Text[Style["r", FontSize -> 12]], {0., -0.0625}]}]}]

results with load stresses undeformed and deformed mesh

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you very much, but please note that I used stress boundary conditions instead of displacement boundary conditions. $\endgroup$ – A little mouse on the pampas Jul 3 at 7:53
  • 2
    $\begingroup$ Stress is nothing else than force exerted per surface. Nothing is lost in the information. I allow for buckling other allow for translation. Mind at least the hole in the middle is fixed. $\endgroup$ – Steffen Jaeschke Jul 3 at 10:07
  • $\begingroup$ Thank you very much for your answers. I understand that you use the NeumannValue function to apply stress boundary conditions. I did not know this knowledge before. $\endgroup$ – A little mouse on the pampas Jul 6 at 1:34
7
+100
$\begingroup$

Your model appears to have quarter symmetry. If one can take advantage of symmetry, it will be a smaller model and may even be easier to setup. A good place to start to find a good setup is PDEModels Overview. Clicking on the Plane Stress will bring you to a verified operator.

It could be useful to use FEMAddOns to difference two boundary meshes so that it will be easy to refine the mesh at the hole.

Build a Quarter Symmetry Mesh

The following workflow will build a quarter symmetry mesh with refinement near the hole.

ResourceFunction["FEMAddOnsInstall"][];
Needs["FEMAddOns`"];
bmesh1 = ToBoundaryMesh[Rectangle[{0, 0}, {1, 1}]];
bmesh2 = ToBoundaryMesh[Disk[{0, 0}, 0.1], 
   MaxCellMeasure -> {"Length" -> .005}];
bmesh = BoundaryElementMeshDifference[bmesh1, bmesh2];
bmesh["Wireframe"];
mesh = ToElementMesh[bmesh];
mesh["Wireframe"]

Mesh

Use Plane Stress Operator From Documentation

The Mathematica documentation provides a plane stress and a plane strain form of the operator. Since the OP diagram shows stress boundary conditions versus displacement boundary conditions, we choose the planes stress operator. I will assume a Young's modulus of 100 and a Poisson ratio of 1/3.

ClearAll[ν, Y]
op = {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}]} /. {Y -> 100, ν -> 1/3};

Set up and solve the PDE system

The OP diagram is very similar to the Kirsch Plate Verification Benchmark. You may find a description in the PDF and PPT files here. The modified Kirsch boundary conditions diagram are shown below (note values are not the same as the OP).

Modified Kirsch Problem

On the x and y symmetry planes, we use Dirichlet Conditions to create the "roller type boundary condition" and fix the u and v displacement, respectively. We then can apply stress NeumannValues on the top (negative for compression) and right boundary (positive for tension) as shown in the following workflow:

dcx = DirichletCondition[u[x, y] == 0., x == 0];
dcy = DirichletCondition[v[x, y] == 0., y == 0];
{ufun, vfun} = 
  NDSolveValue[{op == {NeumannValue[10, x == 1], 
      NeumannValue[-10, y == 1]}, dcx, dcy}, {u, 
    v}, {x, y} \[Element] mesh];
Show[{
  mesh["Wireframe"[ "MeshElement" -> "BoundaryElements"]],
  ElementMeshDeformation[mesh, {ufun, vfun}][
   "Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]
ContourPlot[ufun[x, y], {x, 0, 1}, {y, 0, 1}, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic]
ContourPlot[vfun[x, y], {x, 0, 1}, {y, 0, 1}, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic]

Simulation Results

With the assumed parameters, we are near the limit of deforming the mesh.

Verification

To show that this method gives reasonable results, I will verify the solution verus Kirsch plate benchmark. Since the Kirsch plate benchmark assumes an infinitely long plate, we will expect some end effects. Some useful references will be the previously mentioned COMSOL benchmark and this fracturemechanics.org website. Additionally, it will be useful to download @user21's VonMisesStress funtion located at this answer.

The system we will model is a finite plate in uniaxial tension as shown below:

Kirsch Verification Problem

We will use to @user21's suggestion to create a more accurate mesh using the numerical region.

(*ResourceFunction["FEMAddOnsInstall"][]*) (* Uncomment if you need \
to update version *)
Needs["FEMAddOns`"];
r1 = Rectangle[{0, 0}, {1, 1}];
r2 = Disk[{0, 0}, 0.1];
bmesh1 = ToBoundaryMesh[r1];
bmesh2 = ToBoundaryMesh[r2, MaxCellMeasure -> {"Length" -> .005}];
bmesh = BoundaryElementMeshDifference[bmesh1, bmesh2];
bmesh["Wireframe"];
(* Incorporating user21 suggestion for better accuracy *)
rdiff = RegionDifference[r1, r2];
nr = ToNumericalRegion[rdiff];
SetNumericalRegionElementMesh[nr, bmesh];
mesh = ToElementMesh[nr, MaxCellMeasure -> {"Length" -> .04}];
mesh["Wireframe"]

Now, setup and solve the PDE system.

(* set material parameters *)
materialParameters = {Y -> 2.1*^11, ν -> 0.3};
(* set up factor matrix to be used in subsequent stress calcs *)
pfac = Y/(1 - ν^2)*{{1, ν, 0}, {ν, 1, 0}, {0, 
     0, (1 - ν)/2}};
fac = pfac /. materialParameters;
ClearAll[ν, Y]
op = {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}]} /. materialParameters;
dcx = DirichletCondition[u[x, y] == 0., x == 0];
dcy = DirichletCondition[v[x, y] == 0., y == 0];
{ufun, vfun} = 
  NDSolveValue[{op == {NeumannValue[1000, x == 1], 0}, dcx, dcy}, {u, 
    v}, {x, y} ∈ mesh];
Show[{
  mesh["Wireframe"[ "MeshElement" -> "BoundaryElements"]],
  ElementMeshDeformation[mesh, {ufun, vfun}][
   "Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]
ContourPlot[ufun[x, y], {x, 0, 1}, {y, 0, 1}, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic]
ContourPlot[vfun[x, y], {x, 0, 1}, {y, 0, 1}, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic]

Here, we modify @user21's answer slightly to get additional stress outputs.

ClearAll[VonMisesStress]
VonMisesStress[{uif_InterpolatingFunction, vif_InterpolatingFunction},
   fac_] := 
 Block[{dd, df, mesh, coords, dv, ux, uy, vx, vy, ex, ey, gxy, sxx, 
   syy, sxy}, dd = Outer[(D[#1[x, y], #2]) &, {uif, vif}, {x, y}];
  df = Table[Function[{x, y}, Evaluate[dd[[i, j]]]], {i, 2}, {j, 2}];
  (*the coordinates from the ElementMesh*)
  mesh = uif["Coordinates"][[1]];
  coords = mesh["Coordinates"];
  dv = Table[df[[i, j]] @@@ coords, {i, 2}, {j, 2}];
  ux = dv[[1, 1]];
  uy = dv[[1, 2]];
  vx = dv[[2, 1]];
  vy = dv[[2, 2]];
  ex = ux;
  ey = vy;
  gxy = (uy + vx);
  sxx = fac[[1, 1]]*ex + fac[[1, 2]]*ey;
  syy = fac[[2, 1]]*ex + fac[[2, 2]]*ey;
  sxy = fac[[3, 3]]*gxy;
  {ElementMeshInterpolation[{mesh}, sxx], 
   ElementMeshInterpolation[{mesh}, syy], 
   ElementMeshInterpolation[{mesh}, sxy], 
   ElementMeshInterpolation[{mesh}, 
    Sqrt[(sxy^2) + (syy^2) + (sxx^2)]]}]
{sxxfn, syyfn, sxyfn, vmsfn} = VonMisesStress[{ufun, vfun}, fac];
ContourPlot[vmsfn[x, y], {x, y} \[Element] mesh, 
 RegionFunction -> Function[{x, y, z}, (1/10)^2 < x^2 + y^2], 
 Contours -> 40, ColorFunction -> "TemperatureMap", 
 AspectRatio -> Automatic, PlotPoints -> All, PlotRange -> {0, 3000}, 
 PlotLegends -> Automatic]

Stress Plots

Finally, we can verify the simulation results versus the analytical solution for an infinite plate.

Plot[{1000/2*(2 + (0.1/y)^2 + 3*(0.1/y)^4), sxxfn[0, y]}, {y, 0.1, 1},
  PlotRange -> {900, 3000}, Frame -> True, 
 FrameLabel -> {{"\!\(\*SubscriptBox[\(σ\), \(xx\)]\)", 
    None}, {"y", 
    Style["Analytical Normal Stress vs Simulated @ x=0", Larger]}}, 
 LabelStyle -> Directive[Bold], 
 PlotStyle -> {Directive[Opacity[0.5], Thickness[0.014], Orange], 
   Directive[Dashed, Brown]}, 
 PlotLegends -> 
  Placed[SwatchLegend[{"Analytical", "Simulated"}, 
    LegendMarkers -> "Line", LegendFunction -> "Frame", 
    LegendLayout -> "Column"], {{0.7, 0.75}, {0.5, 1}}]]

Verification Result

Aside from the deviation at the end, the analytical and simulated results match quite closely.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you very much, I want to know how you use these stress boundary conditions, which I did not understand very clearly. $\endgroup$ – A little mouse on the pampas Jul 3 at 7:55
  • 1
    $\begingroup$ @PleaseCorrectGrammarMistakes I put some additional explanation about the boundary conditions. Let me know if that helps. $\endgroup$ – Tim Laska Jul 3 at 16:39
  • $\begingroup$ Thank you for your answers. I probably understood that you applied stress boundary conditions using the NeumannValue function. I did not have this knowledge before. $\endgroup$ – A little mouse on the pampas Jul 6 at 1:30
6
$\begingroup$

This is not an answer but a comment on Tim's answer. Tim's answer is just fine as it is. However, I'd like to take the opportunity to show how to create a mesh that is an even more accurate representation of the geometry; the additional accuracy is most likely not needed in this case but it makes a nice example to show the functionality.

Create a boundary ElementMesh with a refined cut out:

ResourceFunction["FEMAddOnsInstall"][];
Needs["FEMAddOns`"];
r1 = Rectangle[{0, 0}, {1, 1}];
r2 = Disk[{0, 0}, 0.1];
bmesh1 = ToBoundaryMesh[r1];
bmesh2 = ToBoundaryMesh[r2, MaxCellMeasure -> {"Length" -> .005}];
bmesh = BoundaryElementMeshDifference[bmesh1, bmesh2];
bmesh["Wireframe"];

Create a NumericalRegion from the symbolic region difference and the corresponding boundary ElementMesh:

rdiff = RegionDifference[r1, r2];
nr = ToNumericalRegion[rdiff];
SetNumericalRegionElementMesh[nr, bmesh]

Construct a full ElementMesh:

mesh = ToElementMesh[nr];
mesh["Wireframe"]

enter image description here

Compute the difference of the numerical region area and the exact symbolic area:

NIntegrate[1, {x, y} \[Element] mesh] - Area[rdiff]
(* 2.3297*10^-8 *)

Compare to the difference in area between the numeric discretization of the boundary ElementMesh and the exact symbolic area:

NIntegrate[1, {x, y} \[Element] ToElementMesh[bmesh]] - Area[rdiff]
(* 2.65977*10^-6 *)

So, we can squeeze out two orders of magnitude additional accuracy. Consult the documentation for more information on Numerical Regions and Region Approximation Quality or the reference page to ToNumericalRegion.

I have updated the FEMAddOns documentation to include this example.

| improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ Thanks for the enhancement! I will make sure that I incorporate it in the future. A quick question. Does your VonMiseStress function posted here still work? I was having trouble and was wondering if you had updated it. Do I need a structure mesh? Thanks. $\endgroup$ – Tim Laska Jul 3 at 17:59
  • 2
    $\begingroup$ I think that I missed the factor matrix specification, fac in the VonMisesStress function. I think that it is working now. $\endgroup$ – Tim Laska Jul 3 at 18:57

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.