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

I have the following equations $$\frac{\partial{N_x}}{\partial{x}}+\frac{\partial{N_{xy}}}{\partial{y}}=0$$ $$\frac{\partial{N_y}}{\partial{y}}+\frac{\partial{N_{xy}}}{\partial{x}}=0$$ where $$$$\begin{bmatrix}N_{x}\\N_{y}\\N_{xy}\end{bmatrix}=\begin{bmatrix}A_{11}(x)&A_{12}(x)&0\\A_{12}(x)&A_{22}(x)&0\\0&0&A_{66}(x)\end{bmatrix}\begin{bmatrix}\frac{\partial{u}}{\partial{x}}\\\frac{\partial{v}}{\partial{y}}\\\frac{\partial{u}}{\partial{y}}+\frac{\partial{v}}{\partial{x}}\end{bmatrix}$$$$ $$A_{ij}$$ are varies only on $$x$$ direction. The boundary conditions are shown bellow. I'm trying to apply the complete panel BCs.

Can someone help me to model this problem. I looked the FEM package, but I still can't solved. Thank you in advance. I'm looking for solution to u(x,y) e 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 = -1;

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) Abs[x] + 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["NDSolveFEM"]
omega = Rectangle[{-a/2, -b/2}, {a, b}];
mesh = ToElementMesh[omega];

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], {x, 1}], D[Ny[x, y], {y, 1}] +
Nxy[x, y], {y, 1}]};
gammaD = {DirichletCondition[u[x, y] == u0, x == a/2],
DirichletCondition[u[x, y] == -u0, x == -a/2]};
NDSolve[{PDEs == {0, 0}, gammaD}, {u,
v}, {x, y} \[Element] mesh];


I'm getting the following error:

NDSolve::femcnmd: The PDE coefficient {{0.000508 (-8.22842*10^11 Cos[Times[<<2>>]]^3 Sin[1.5708 Abs[<<1>>]] (Abs^\[Prime])[x]-4.04788*10^10 Cos[1.5708 Abs[<<1>>]] Sin[Times[<<2>>]]^3 (Abs^\[Prime])[x])+0.000254 (<<1>>)+<<1>>+1/2 (-0.001016 (-8.22842*10^11 Power[<<2>>] Sin[<<1>>] (<<1>>^(<<1>>))[<<1>>]-4.04788*10^10 Cos[<<1>>] Power[<<2>>] (<<1>>^(<<1>>))[<<1>>])-0.000508 (2.32093*10^11 Power[<<2>>] (<<1>>^(<<1>>))[<<1>>]-7.45065*10^11 Power[<<2>>] Power[<<2>>] (<<1>>^(<<1>>))[<<1>>]+1.62623*10^10 Power[<<2>>] (<<1>>^(<<1>>))[<<1>>]+1.67953*10^10 Plus[<<3>>])-0.000508 (-2.32093*10^11 Power[<<2>>] (<<1>>^(<<1>>))[<<1>>]+7.45065*10^11 Power[<<2>>] Power[<<2>>] (<<1>>^(<<1>>))[<<1>>]-1.62623*10^10 Power[<<2>>] (<<1>>^(<<1>>))[<<1>>]+1.67953*10^10 Plus[<<3>>])),<<1>>}} does not evaluate to a numeric matrix of dimensions {1,2} at the coordinate {-0.4625,-0.4625}; it evaluated to {{0.,5.79879*10^6 (Abs^\[Prime])[-0.4625]}} instead.


Can someone help me? Thank you all in advance.

• Obvious typos: use Pi instead of pi, use Sqrt[x^2] instead of Abs[x]. Then put omega = Rectangle[{-a/2, -b/2}, {a/2, b/2}]; and mesh = ToElementMesh[omega, MaxCellMeasure -> 0.01]. u0=-1 is too much, take u0=-0.01. – Alex Trounev Nov 13 '19 at 11:31
• How do you know that Ny=0 at y=b/2? – Alex Trounev Nov 13 '19 at 13:57
• Your equations do not have the standard form $\partial \sigma _{ik}/\partial x_k=0$. Or is there a typo in PDEs? – Alex Trounev Nov 13 '19 at 14:25
• $N_{ij}$ are $\int^{-h/2}_{h/2}\sigma_{ij}dz$ – Professor P. Cosmo Klunk Nov 13 '19 at 14:36
• I fixed all the typos. See my answer. – Alex Trounev Nov 13 '19 at 14:45

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["NDSolveFEM"]
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} =
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]}


• BCondition DirichletCondition[v[x, y] == 0, y == 0] means that function v[x,y] is antisymmetrical around axes y=0. Also Your demand Nxy=0 looks strange, because with such a deformation we have what we have. – Alex Trounev Nov 13 '19 at 16:30
• Thank you, Alex! You rock! – Professor P. Cosmo Klunk Nov 14 '19 at 12:18
• @DiegoMagela You're welcome! – Alex Trounev Nov 14 '19 at 13:28

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.

• It worked. Thank you so much! – Professor P. Cosmo Klunk Nov 13 '19 at 1:55
• Glad to hear that! You're welcome. – Henrik Schumacher Nov 13 '19 at 2:18