6
$\begingroup$

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{equation} \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}\end{equation} $A_{ij}$ are varies only on $x$ direction. The boundary conditions are shown bellow. I'm trying to apply the complete panel BCs.

99

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).

Follow my code:

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["NDSolve`FEM`"]
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.

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ Nov 13 '19 at 11:31
  • $\begingroup$ How do you know that Ny=0 at y=b/2? $\endgroup$ Nov 13 '19 at 13:57
  • $\begingroup$ Your equations do not have the standard form $\partial \sigma _{ik}/\partial x_k=0$. Or is there a typo in PDEs? $\endgroup$ Nov 13 '19 at 14:25
  • $\begingroup$ $N_{ij}$ are $\int^{-h/2}_{h/2}\sigma_{ij}dz$ $\endgroup$ Nov 13 '19 at 14:36
  • $\begingroup$ I fixed all the typos. See my answer. $\endgroup$ Nov 13 '19 at 14:45
4
$\begingroup$

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]}

Figure 1

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$ Nov 13 '19 at 16:30
  • $\begingroup$ Thank you, Alex! You rock! $\endgroup$ Nov 14 '19 at 12:18
  • $\begingroup$ @DiegoMagela You're welcome! $\endgroup$ Nov 14 '19 at 13:28
4
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ It worked. Thank you so much! $\endgroup$ Nov 13 '19 at 1:55
  • $\begingroup$ Glad to hear that! You're welcome. $\endgroup$ Nov 13 '19 at 2:18

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.