Needs["NDSolve`FEM`"];
G = 6.894745 10^9;
E1 = 26.25 G;
E2 = 1.49 G;
G12 = 1.04 G;
nu12 = 0.28;
nu21 = (E2*nu12)/E1;
t = 0.0050 .0254;
a = 1; b = 1;
u0 = .01;
Son = {{1/E1, -nu12/E1, 0}, {-nu21/E2, 1/E2, 0}, {0, 0, 1/G12}};
Qon = Inverse[Son];
Do[
Do[
angles = {{angle0, angle1}, {-angle0, -angle1}, {angle0,
angle1}, {-angle0, -angle1}, {angle0,
angle1}, {-angle0, -angle1}, {-angle0, -angle1}, {angle0,
angle1}, {-angle0, -angle1}, {angle0,
angle1}, {-angle0, -angle1}, {angle0, angle1}};
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]] ;
T1 = angles[[i, 2]] ;
theta[x_] := ((2/a) (T1 - T0) Sqrt[x^2] + T0) (Pi/180);
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;
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]];
D11[x] = mD[[1, 1]]; D12[x] = mD[[1, 2]]; D16[x] = mD[[1, 3]];
D22[x] = mD[[2, 2]]; D26[x] = mD[[2, 3]]; D66[x] = mD[[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/2, x == a/2],
DirichletCondition[u[x, y] == u0/2, x == -a/2]
};
omega = Rectangle[{-a/2, -b/2}, {a/2, b/2}];
mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001];
{U, V} =
NDSolveValue[
{
PDEs == {0, 0},
gammaD,
DirichletCondition[v[x, y] == 0, y == b/2],
DirichletCondition[v[x, y] == 0, y == -b/2]
}
,
{u, v}, {x, y} \[Element] mesh
];
Nx[x_, y_] =
A11 D[U[x, y], {x, 1}] + A12 D[V[x, y], {y, 1}] +
A16 (D[U[x, y], {y, 1}] + D[V[x, y], {x, 1}]);
Ny[x_, y_] =
A12 D[U[x, y], {x, 1}] + A22 D[V[x, y], {y, 1}] +
A26 (D[U[x, y], {y, 1}] + D[V[x, y], {x, 1}]);
Nxy[x_, y_] =
A16 D[U[x, y], {x, 1}] + A26 D[V[x, y], {y, 1}] +
A66 (D[U[x, y], {y, 1}] + D[V[x, y], {x, 1}]);
nL = 20; nC = 20;
K = ConstantArray[0, {nL*nC, nL*nC}];
M = ConstantArray[0, {nL*nC, nL*nC}];
Do[
Do[
Do[
Do[
k = (m - 1) nC + n;
l = (p - 1) nC + q;
If[n == q,
K[[k, l]] = (\[Pi]^4 b m^2 p^2)/(2 a^4)
NIntegrate[
D11[x] Sin[(m \[Pi] (x + 0.5 a))/a] Sin[(
p \[Pi] (x + 0.5 a))/a], {x, -a/2, a/2},
Method -> "Trapezoidal"]
+ (\[Pi]^4 n^2 q^2)/(2 b^3)
NIntegrate[
D22[x] Sin[(m \[Pi] (x + 0.5 a))/a] Sin[(
p \[Pi] (x + 0.5 a))/a], {x, -a/2, a/2},
Method -> "Trapezoidal"], K[[k, l]] = 0]
, {q, 1, nC}],
{p, 1, nL}],
{n, 1, nC}],
{m, 1, nL}];
, {angle0, 0, 90, 1}]
, {angle1, 0, 90, 1}]