I'm trying to get the displacements u(x,y) and v(x,y) of a beam simply suported, the stress distribution in the middle section should be a linear function:

but i'm getting this:

My code is:

Needs["NDSolve`FEM`"];
PS = {
  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}]
  }

L = 2;
q = 6000;
Propiedades = {Y -> 205940000000, \[Nu] -> 30/100};

h1 = 1/2;
h2 = 2;
h3 = 3;

Reg1 = Rectangle[{0, 0}, {L, h1}];
Reg2 = Rectangle[{0, 0}, {L, h2}];
Reg3 = Rectangle[{0, 0}, {L, h3}];

Mesh1 = ToElementMesh[Reg1, MeshQualityGoal -> 0];
Mesh2 = ToElementMesh[Reg2];
Mesh3 = ToElementMesh[Reg3];
{u1, v1, \[Sigma]x1, \[Sigma]y1, \[Tau]xy1} = NDSolveValue[{
     PS == {0, NeumannValue[-q, {0 <= x <= L, y == h1}]},
     \[Sigma]x[x, y] == 
      Y/(1 - \[Nu]^2) (D[u[x, y], x] + \[Nu] D[v[x, y], y]),
     \[Sigma]y[x, y] == 
      Y/(1 - \[Nu]^2) (D[v[x, y], y] + \[Nu] D[u[x, y], x]),
     \[Sigma]xy[x, y] == (Y*\[Nu])/(
       1 - \[Nu]^2) (D[u[x, y], x] + D[v[x, y], y]),
     DirichletCondition[v[x, y] == 0, {x == 0, y == 0}],
     DirichletCondition[v[x, y] == 0, {x == L, y == 0}],
     DirichletCondition[u[x, y] == 0, {x == 0, y == 0}]
     } /. Propiedades, {u, 
    v, \[Sigma]x, \[Sigma]y, \[Sigma]xy}, {x, y} \[Element] Mesh1];

DMesh1 = ElementMeshDeformation[Mesh1, {u1, v1}, 
   "ScalingFactor" -> 2000000];

Row[{
  Show[{Mesh1[
     "Wireframe"[
      "ElementMeshDirective" -> 
       Directive[EdgeForm[Gray], FaceForm[]]]], 
    Graphics[{EdgeForm[Thickness[0.001]], RGBColor[0, 0, 0, 0.1], 
      Reg1}]}, ImageSize -> 300, Epilog -> {
     Scale[Translate[Apoyo2, {-0.5, -0.5}], 0.2],
     Scale[Translate[Apoyo1, {-0.5 + L, -0.5}], 0.2]
     }, PlotRange -> {{-0.1, L + 0.1}, {-0.15, h1 + 0.15}}],
  Show[{
    Mesh1[
     "Wireframe"[
      "ElementMeshDirective" -> 
       Directive[EdgeForm[Gray], FaceForm[]]]],
    DMesh1[
     "Wireframe"[
      "ElementMeshDirective" -> 
       Directive[EdgeForm[RGBColor[0, 0.3, 0.8]], FaceForm[]]]]
    }, ImageSize -> 300]
  }]

Plot[\[Sigma]x1[L/2, y]/1000, {y, 0, h1}, Filling -> Axis, 
 AxesLabel -> {"h[m]", 
   "\!\(\*SubscriptBox[\(\[Sigma]\), \(x\)]\)[kPa]"}, 
 ImageSize -> 400]

I don't why is this happening, but when i set the y coordinate in the Dirichlet Condition to y=h1/2 i get the correct stress distribution