Error trying to solve a 2D elasticity problem using FEM package

Use RealAbs instead of Abs. Abs is supposed to be the absolute value function for complex functions and thus Mathematica refused to take its derivative.

PDEs do not have the standard form $\frac {\partial \sigma _{ik}}{\partial x_k}=0$. After fixing all the typos, I also added the boundary condition for v[x,y].

E1 = 147 10^9; E2 = 10.3 10^9; G12 = 
 7 10^9; nu12 = 0.27; nu21 = (E2*nu12)/E1; t = 
 0.127 10^-3; a = 1; b = 1; u0 = -.01;

Son = {{1/E1, -nu12/E1, 0}, {-nu21/E2, 1/E2, 0}, {0, 0, 1/G12}};
Qon = Inverse[Son];
angles = {{0, 45}, {0, -45}, {0, -45}, {0, 45}};
num = Dimensions[angles][[1]];
h = num*t;
pos = Table[0, num + 1];
pos[[1]] = -h/2;
For[i = 2, i <= num + 1, i++, pos[[i]] = pos[[i - 1]] + t;]
mA = ConstantArray[0, {3, 3}];
mB = ConstantArray[0, {3, 3}];
mD = ConstantArray[0, {3, 3}];

For[i = 1, i <= num, i++, T0 = angles[[i, 1]] (Pi/180);
  T1 = angles[[i, 2]] (Pi/180);
  theta[x_] := (2/a) (T1 - T0) Sqrt[x^2] + T0;
  m = Cos[theta[x]];
  n = Sin[theta[x]];
  Q11 = Qon[[1, 1]];
  Q12 = Qon[[1, 2]];
  Q22 = Qon[[2, 2]];
  Q66 = Qon[[3, 3]];
  Qxx = m^4*Q11 + n^4*Q22 + 2*m^2*n^2*Q12 + 4*m^2*n^2*Q66;
  Qyy = n^4*Q11 + m^4*Q22 + 2*m^2*n^2*Q12 + 4*m^2*n^2*Q66;
  Qxy = m^2*n^2*Q11 + m^2*n^2*Q22 + (m^4 + n^4)*Q12 + -4*m^2*n^2*Q66;
  Qss = m^2*n^2*Q11 + m^2*n^2*Q22 - 2*m^2*n^2*Q12 + (m^2 - n^2)^2*Q66;
  Qxs = m^3*n*Q11 - m*n^3*Q22 + (m*n^3 - m^3*n)*Q12 + 
    2*(m*n^3 - m^3*n)*Q66;
  Qys = m*n^3*Q11 - m^3*n*Q22 + (m^3*n - m*n^3)*Q12 + 
    2*(m^3*n - m*n^3)*Q66;
  Qoff = {{Qxx, Qxy, Qxs}, {Qxy, Qyy, Qys}, {Qxs, Qys, Qss}};
  mA = mA + Qoff*(pos[[i + 1]] - pos[[i]]);
  mB = mB + Qoff*(pos[[i + 1]]^2 - pos[[i]]^2);
  mD = mD + Qoff*(pos[[i + 1]]^3 - pos[[i]]^3);];
mB = mB/2;
mD = mD/3;

Needs["NDSolve`FEM`"]
omega = Rectangle[{-a/2, -b/2}, {a/2, b/2}];
A11[x_] = mA[[1, 1]]; A12[x_] = mA[[1, 2]]; A16[x_] = mA[[1, 3]];
A22[x_] = mA[[2, 2]]; A26[x_] = mA[[2, 3]]; A66[x_] = mA[[3, 3]];
Nx[x_, y_] = 
  A11[x] D[u[x, y], {x, 1}] + A12[x] D[v[x, y], {y, 1}] + 
   A16[x] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);
Ny[x_, y_] = 
  A12[x] D[u[x, y], {x, 1}] + A22[x] D[v[x, y], {y, 1}] + 
   A26[x] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);
Nxy[x_, y_] = 
  A16[x] D[u[x, y], {x, 1}] + A26[x] D[v[x, y], {y, 1}] + 
   A66[x] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

PDEs = {D[Nx[x, y], {x, 1}] + D[Nxy[x, y], {y, 1}], 
   D[Ny[x, y], {y, 1}] + D[Nxy[x, y], {x, 1}]};
gammaD = {DirichletCondition[u[x, y] == u0, x == a/2], 
   DirichletCondition[u[x, y] == -u0, x == -a/2]};


mesh = ToElementMesh[omega, MaxCellMeasure -> 0.01]; {U, V} = 
 NDSolveValue[{PDEs == {0, 0}, gammaD, 
   DirichletCondition[v[x, y] == 0, y == 0]}, {u, 
   v}, {x, y} \[Element] mesh];

{ContourPlot[U[x, y], {x, y} \[Element] mesh, Contours -> 20, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic], 
 ContourPlot[V[x, y], {x, y} \[Element] mesh, Contours -> 20, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic]}