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.

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

2 Answers 2

Reset to default
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$
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$
    – Alex Trounev
    Commented Nov 13, 2019 at 16:30
  • $\begingroup$ Thank you, Alex! You rock! $\endgroup$
    – Ponyboy Curtis
    Commented Nov 14, 2019 at 12:18
  • $\begingroup$ @DiegoMagela You're welcome! $\endgroup$
    – Alex Trounev
    Commented Nov 14, 2019 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$
2
  • $\begingroup$ It worked. Thank you so much! $\endgroup$
    – Ponyboy Curtis
    Commented Nov 13, 2019 at 1:55
  • $\begingroup$ Glad to hear that! You're welcome. $\endgroup$
    – Henrik Schumacher
    Commented Nov 13, 2019 at 2:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.