Off[General::spell];
Off[General::spell1];

E1 = 17.3;
E2 = 1.27;
G23 =  .4;
ν12 =  .43;
α1 = -3.04*(10^(-7));
α2 = 4.08*(10^(-5));
Tf = 155;
Tref = 27;
ΔT = Tref - Tf;
t0 =  .008;
t90 =  .008;
d0 = 1;

σxo[x_, y_, d90_] = 
        A[d90]*Cosh[R*λx*x] + B[d90]*Cosh[R*λy*y] + 
          Kx;
σx90[x_, y_, 
          d90_] = (((t0 + t90)*σxbar - t0*σxo[x, y, d90])/t90);

σyo[x_, y_, d90_] = 
        A[d90]*((ν12*E2*(t0 + t90))/(t0*E2 + t90*E1))*
            Cosh[R*λx*x] + 
          B[d90]*((t90*E2 + t0*E1)/(ν12*E2*(t0 + t90)))*
            Cosh[R*λy*y] + Ky;
σy90[x_, y_, 
          d90_] = (((t0 + t90)*σybar - t0*σyo[x, y, d90])/t90);

FEMσx90 = {{1, 13.76}, {10, 11.98}, {20, 9.16}, {30, 
            8.25}, {40, 6.73}, {50, 5.57}};

λx = (Sqrt[((3*G23)/(t0*(t90^2)))*(((t0*E1 + t90*E2))/(E1*E2))]);
λy = (Sqrt[((3*G23)/((t0^2)*t90))*(((t90*E1 + t0*E2))/(E1*E2))]);
μx = (Sqrt[((3*G23)/(t0*(t90^2)))*(t0 + t90)*ν(12/E1)]);
μy = (Sqrt[((3*G23)/((t0^2)*t90))*(t0 + t90)*ν(12/E1)]);
R = (Sqrt[1 - ((μ(x^2)*μ(y^2))/(λ(x^2)*λ(y^2)))]);
A[
          d90_] = (-((((t0 + t90))/t0) σxbar - 
                      Kx)*(((t0*E1 + t90*E2))/((t0 + t90)*ν12*E2))*Cosh[((R*λy)/(2*d0))] - 
              Ky*((2*d0)/(R*λy))*
                Sinh[((R*λy)/(2*d0))])/(((4*ν12*E2*d0*d90*(t0 + t90))/((t90*E1 + t0*E2)*(R^2)*λx*λy))*
                Sinh[((R*λx)/(2*d90))]*
                Sinh[((R*λy)/(2*d0))] - (((t0*E1 + t90*E2))/((t0 + t90)*ν12*E2))*Cosh[((R*λx)/(2*d90))]*
                Cosh[((R*λy)/(2*d0))]);
B[
          d90_] = (((((t0 + t90))/t0) σxbar - 
                    Kx)*(((t0 + t90)*ν12*E2)/((t90*E1 + t0*E2)))*((2*d90)/(R*λx))*Sinh[((R*λx)/(2*d90))] + 
              Ky*Cosh[((R*λx)/(2*d90))])/(((4*ν12*E2*d0*d90*(t0 + t90))/((t90*E1 + t0*E2)*(R^2)*λx*λy))*
                Sinh[((R*λx)/(2*d90))]*
                Sinh[((R*λy)/(2*d0))] - (((t0*E1 + t90*E2))/((t0 + t90)*ν12*E2))*Cosh[((R*λx)/(2*d90))]*
                Cosh[((R*λy)/(2*d0))]);
Kx = ((μ(x^2)*μ(y^2))/(μ(x^2)*μ(y^2) - λ(x^2)*λ(y^2)))*((1 - λ((y^2)/μ)(y^2)*
                      (E1/(ν12*E2)))*σxbar + (λ((y^2)/μ)(y^2) - (1/ν)12)*σybar);
Ky = ((μ(x^2)*μ(y^2))/(μ(x^2)*μ(y^2) - λ(x^2)*λ(y^2)))*((λ((x^2)/μ)(x^2) - 
                    (E1/(ν12*E2)))*σxbar + (1 - λ((x^2)/μ)(x^2)*(1/ν)12)*σybar);

σxbar =  .102;
σybar = 0;

$TextStyle = {FontFamily -> "Times", FontWeight -> Bold, 
          FontSize -> 18};
Plot[σx90[0, 0, d90]*1000, {d90, 1, 100}, PlotRange -> All, 
      Frame -> True, 
      FrameLabel -> {"Cracks/inch", "Axial Stress"}, 
      PlotLabel -> "90 layer"]
Plot[σxo[0, 0, d90]*1000, {d90, 1, 100}, PlotRange -> All, 
      Frame -> True, 
      FrameLabel -> {"Cracks/inch", "Axial Stress"}, 
      PlotLabel -> "0 layer"]

Cracked Elastic Constants

Clear[σxbar, σybar]

E10 = E1;
E20 = E2;
ν120 = ν12;
ν210 = (E20/E10)*ν120;
E190[1] = E1;
E290[1] = E2;
ν1290[1] = ν12;
ν2190[1] = ν210;
cdsat = 100;

\(Do[A11[
          d90] = \(2*t0*E10\)\/\((1 - ν120*ν210)\) + \(2*t90*E290[d90]\)\/\((1 - ν1290[d90]*ν2190[d90])\); 
A22[
          d90] = \(2*t0*E20\)\/\((1 - ν120*ν210)\) + \(2*t90*E190[d90]\)\/\((1 - ν1290[d90]*ν2190[d90])\); 
A12[
          d90] = \(2*t0*ν210*E10\)\/\((1 - ν120*ν210)\) + \(2*t90*ν2190[d90]*E190[d90]\)\/\((1 - ν1290[d90]*ν2190[d90])\); 
Ex[d90] = \(A11[d90]*A22[d90] - A12[d90]\^2\)\/\(2*\((t0 + t90)\)*A22[d90]\); 

      Ey[d90] = \(A11[d90]*A22[d90] - A12[d90]\^2\)\/\(2*\((t0 + t90)\)*A11[d90]\); 
νxy[d90] = A12[d90]\/A22[d90]; 
νyx[d90] = 
        A12[d90]\/A11[d90];, \ {d90, 1, cdsat, 1}]\)

ϵxbar0[d90_] = 
        (1/E10)*d0*
            d90*(∫_(((-1)/(2*d0)))^((1/(2*d0))))(∫_(((-1)/(2*d90)))^((1/(2*d90))))σxo[x, y, 
                    d90] \[DifferentialD]x \[DifferentialD]y - ν(120/E10)*d0*
            d90*(∫_(((-1)/(2*d0)))^((1/(2*d0))))(∫_(((-1)/(2*d90)))^((1/(2*d90))))σyo[x, y, 
                    d90] \[DifferentialD]x \[DifferentialD]y;
ϵxbar90[d90_] = 
        (1/E290)[d90]*d0*
            d90*(∫_(((-1)/(2*d0)))^((1/(2*d0))))(∫_(((-1)/(2*d90)))^((1/(2*d90))))σx90[x, y, 
                    d90] \[DifferentialD]x \[DifferentialD]y - ν1290[d90]/E190[d90]*d0*
            d90*(∫_(((-1)/(2*d0)))^((1/(2*d0))))(∫_(((-1)/(2*d90)))^((1/(2*d90))))σy90[x, y, 
                    d90] \[DifferentialD]x \[DifferentialD]y;
ϵybar0[d90_] = 
        (1/E20)*d0*
            d90*(∫_(((-1)/(2*d0)))^((1/(2*d0))))(∫_(((-1)/(2*d90)))^((1/(2*d90))))σyo[x, y, 
                    d90] \[DifferentialD]x \[DifferentialD]y - ν(210/E20)*d0*
            d90*(∫_(((-1)/(2*d0)))^((1/(2*d0))))(∫_(((-1)/(2*d90)))^((1/(2*d90))))σxo[x, y, 
                    d90] \[DifferentialD]x \[DifferentialD]y;
ϵybar90[d90_] = 
        (1/E190)[d90]*d0*
            d90*(∫_(((-1)/(2*d0)))^((1/(2*d0))))(∫_(((-1)/(2*d90)))^((1/(2*d90))))σy90[x, y, 
                    d90] \[DifferentialD]x \[DifferentialD]y - ν2190[d90]/E290[d90]*d0*
            d90*(∫_(((-1)/(2*d0)))^((1/(2*d0))))(∫_(((-1)/(2*d90)))^((1/(2*d90))))σx90[x, y, 
                    d90] \[DifferentialD]x \[DifferentialD]y;

\(Do[eqn1[d90] = 
        Coefficient[ϵxbar0[d90], σxbar] - 
            Coefficient[ϵxbar90[d90], σxbar] == 0; 

      eqn2[d90] = 
        Coefficient[ϵxbar0[d90], σybar] - 
            Coefficient[ϵxbar90[d90], σybar] == 0; 

      eqn3[d90] = 
        Coefficient[ϵybar0[d90], σxbar] - 
            Coefficient[ϵybar90[d90], σxbar] == 0; 

      eqn4[d90] = 
        Coefficient[ϵybar0[d90], σybar] - 
            Coefficient[ϵybar90[d90], σybar] == 
          0;, {d90, 2, cdsat, 1}]\)

\(Do[soln[d90] = 
        Solve[{eqn1[d90], eqn2[d90], eqn3[d90], eqn4[d90]}, {E190[d90], 
            E290[d90], ν1290[d90], ν2190[d90]}]; 

      nu1290[d90] = \(soln[d90]\)[\([1, 1, 2]\)]; 

      nu2190[d90] = \(soln[d90]\)[\([1, 2, 2]\)]; 

      Em190[d90] = \(soln[d90]\)[\([1, 3, 2]\)]; 

      Em290[d90] = \(soln[d90]\)[\([1, 4, 2]\)]; 
E190[d90] = Em190[d90]; 

      E290[d90] = Em290[d90]; 
ν1290[d90] = 
        nu1290[d90]; 
ν2190[d90] = nu2190[d90];, {d90, 2, cdsat, 
        1}]\)

ListPlot[Table[{d90, E290[d90]}, {d90, 2, cdsat}], 
      PlotRange -> All, Frame -> True]

Do[Ex[d90] = ((E10 + E290[d90])/2), {d90, 1, cdsat}]

Exi = ((((((((-E1)/E2) + ν(12^2))*((t0/t90) + 1)*((E1/E2) + (t0/t90)))/(E1*(-((E1/E2)*(t0/t90) + 1)*((E1/E2) + (t0/t90)) + ν(12^2)*((((t0/t90) + 1))^2))))))^(-1));

Do[Gmaxm[
          d90] = ((t0/t90) + 
              1)*σ((xbar^2)/(R*λx))*(1/Exi)*(1/(t0/t90))*((((ν1290[d90]^2*(1 + (t0/t90)) - (E190[d90]/E290[d90] + (t0/t90))))^2)/((E190[d90]/E290[d90] - ν1290[d90]^2)*(((E190[d90]/E290[d90] + (t0/t90)))^2))); 
Gmaxt[d90] = 
        (1/(R*λx))*
          (t0/t90)*((((α1 - α2))^2) Δ(T^2) E190[d90]^2*(((-(E190[d90]/E290[d90] + (t0/t90)) + ν1290[d90]*((t0/t90) + 1)))^2))/(Exi*((t0/t90) + 
                    1)*(((-E190[d90])/E290)[d90] + ν(12^2))*(((E190[d90]/E290[d90] + (t0/t90)))^2)); 

Gmax[d90] = Gmaxm[d90] + Gmaxt[d90];, {d90, 1, cdsat, 
        1}]

σfpf = 2500/ .032/1000000;
d90i = 1/4;
ηi = ((2*d90i)/(R*λx));

Gc = 
      Gmax[1]*Tanh[(1/η)i] - 
        (1/η)i*(1 - (((Tanh[(1/η)i]))^2))

Gcrit = Gc /. σxbar -> σfpf

σxbar =  .102;

$TextStyle = {FontFamily -> "Times", FontWeight -> Bold, 
          FontSize -> 14};

ListPlot[Table[{d90, Gcrit/Gmax[d90]}, {d90, 1, 20}], 
        Frame -> True, 
        FrameLabel -> {"Cracks/inch", "Gc/Gmax"}, 
        PlotLabel -> Room Temperature, 
        PlotJoined -> True];

Gcrit/Gmax[27]