Finite Element Analysis Boundary Condition problem

Question

I'm trying to learn how to use Mathematica to perform FEA. In this case i'm trying to perform a compressive test on a disk. I've started with meshing the component.

Reg1 = Annulus[{0, 0}, {14, 25}];
Mesh1 = ToElementMesh[Reg1, MeshQualityGoal -> 0];

The compressive force is 6.41 Mpa

q = 6.41;

The operator for the plane stress analysis is provided by Mathematica and it's:

op = {Inactive[
   Div][({{0, -((Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(
       2 (1 - \[Nu]^2))), 0}}.Inactive[Grad][
     v[x, y], {x, y}]), {x, y}] + 
 Inactive[
   Div][({{-(Y/(1 - \[Nu]^2)), 
      0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}.Inactive[
      Grad][u[x, y], {x, y}]), {x, y}], 
Inactive[
   Div][({{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((
       Y \[Nu])/(1 - \[Nu]^2)), 0}}.Inactive[Grad][
     u[x, y], {x, y}]), {x, y}] + 
 Inactive[
   Div][({{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))), 
      0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
     v[x, y], {x, y}]), {x, y}]} /. {Y -> 3416, \[Nu] -> 0.33};

And now I have problems with the definition of the boundary conditions:

bc = {0, NeumannValue[-q, x == 0 && y==25],NeumannValue[q, x == 0 && y==-25]};


{uif, vif} = 
NDSolveValue[{op == {0, NeumannValue[-q, x == 0]}}, {u, v}, {14^2 < x^2 + y^2 < 25^2}]


Mesh1 = uif["ElementMesh"];
ElementMeshDeformation[Mesh1, {uif, vif}]["Wireframe"]

Obviously this is wrong but I can't understand why and how to set my boundary conditions. **The image above is referring to another type of mesh.

Answer 1

Seems like one needs to have an extended area for the boundary conditions to work properly, so I used z<24.5 instead of z==-25 as you had it. Changed a bit the code but I hope you get the idea:

<< NDSolve`FEM`
ClearAll[x, y, v, u, mesh1, uif, vif]
Reg1 = Annulus[{0, 0}, {14, 25}];

mesh = ToElementMesh[Reg1, "MaxCellMeasure" -> 0.5]
planeStress = {Inactive[
       Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 \
- \[Nu]^2)), 0}}.Inactive[Grad][v[y, z], {y, z}], {y, z}] + 
     Inactive[
       Div][{{-(Y/(1 - \[Nu]^2)), 
         0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}}.Inactive[Grad][
        u[y, z], {y, z}], {y, z}], 
    Inactive[
       Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, {-((Y*\[Nu])/(1 \
- \[Nu]^2)), 0}}.Inactive[Grad][u[y, z], {y, z}], {y, z}] + 
     Inactive[
       Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)), 
         0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
        v[y, z], {y, z}], {y, z}]} /. {Y -> 10^3, \[Nu] -> 33/100};
{uif, vif} = 
  NDSolveValue[{planeStress == {0, -NeumannValue[60, z >= 24.5]}, 
    DirichletCondition[u[y, z] == 0, z < -24.5], 
    DirichletCondition[v[y, z] == 0, z < -24.5]}, {u, 
    v}, {y, z} \[Element] mesh];
dmesh = ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 1];
Show[{mesh["Wireframe"], 
  dmesh["Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]