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$
  • $\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$ – Alex Trounev Nov 13 '19 at 11:31
  • $\begingroup$ How do you know that Ny=0 at y=b/2? $\endgroup$ – Alex Trounev 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$ – Alex Trounev Nov 13 '19 at 14:25
  • $\begingroup$ $N_{ij}$ are $\int^{-h/2}_{h/2}\sigma_{ij}dz$ $\endgroup$ – Magela Nov 13 '19 at 14:36
  • $\begingroup$ I fixed all the typos. See my answer. $\endgroup$ – Alex Trounev 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

| improve this answer | |
$\endgroup$
  • 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$ – Alex Trounev Nov 13 '19 at 16:30
  • $\begingroup$ Thank you, Alex! You rock! $\endgroup$ – Magela Nov 14 '19 at 12:18
  • $\begingroup$ @DiegoMagela You're welcome! $\endgroup$ – Alex Trounev 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.

| improve this answer | |
$\endgroup$
  • $\begingroup$ It worked. Thank you so much! $\endgroup$ – Magela Nov 13 '19 at 1:55
  • $\begingroup$ Glad to hear that! You're welcome. $\endgroup$ – Henrik Schumacher 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.