Another problem with my incrementing

Question

Thomas on 1 May 2015

0
Link

Direct link to this question

https://nl.mathworks.com/matlabcentral/answers/214517-another-problem-with-my-incrementing

Commented: Geoff Hayes on 1 May 2015

I am having trouble incrementing my dh values all the way through. It's a lot of code, and I've worked on it all, day, but can't find out what I'm not doing right.

%   Fiber
Vf     = 0.6; 
Ef     = 250;               % GPa
rho_f  = 1.8;               % g/(cm^3)
nu_f   = 0.2;
Gf     = Ef / (2*(1 + nu_f)); % GPa
%   Matrix
Vm     = 0.4;
Em     = 3.5;               % GPa
rho_m  = 1.2;               % g/(cm^3)
nu_m   = 0.35;
Gm     = Em / (2*(1 + nu_m)); % GPa
%   Aluminum
E_Al   = 72;   % GPa
rho_Al = 2.78; % g/(cm^3)
nu_Al  = 0.33;
G_Al   = E_Al / (2*(1 + nu_Al)); % GPa
% 3. Please clearly state your assumptions in the development of your
% assessment. Tabulate your material properties in a consistent set of
% units (SI!!).
%
% Part (a)
% Determine the corresponding A, B and D matrix for baseline Al panel and
% each case of composite panels, in which the ply angle in the S2Ep layer
% is assumed to 0, 30, 45, 60, 90 degrees. Each layer thickness =1mm.
% Givens:
%   S2Ep layer
theta = [0 30 45 60 90]; % deg
% GLARE Layout
%   AL
%   S2Ep theta deg.
%   AL
%   S2Ep theta deg.
%   AL
% Calculating Composite (S2Ep) Properties
E_L   = (Vf*Ef) + (Vm*Em);             % GPa
E_T   = (Ef*Em) / ((Ef*Vm) + (Em*Vf)); % GPa
nu_LT = (Vf*nu_f) + (Vm*nu_m);
nu_TL = (E_T*nu_LT) / E_L;
rho_c = (Vf*rho_f) + (Vm*rho_m);
G_LT  = (Gf*Gm) / ((Gf*Vm) + (Gm*Vf)); % GPa
% Calculating Composite (S2Ep) Q Matrix
q11    = E_L / (1-(nu_LT*nu_TL));       % GPa
q22    = E_T / (1-(nu_LT*nu_TL));       % GPa
q12    = nu_LT*E_T / (1-(nu_LT*nu_TL)); % GPa
q66    = G_LT;                          % GPa
Q_S2Ep = [ q11  q12   0;
           q12  q22   0;                % GPa
            0    0   q66];
% Caculating Aluminum Q Matrix
q11_Al = E_Al / (1-(nu_Al^2));          % GPa
q12_Al = nu_Al*q11_Al;                     % GPa
q66_Al = E_Al / (2*(1+nu_Al));          % GPa
Q_Al   = [ q11_Al  q12_Al   0;
           q12_Al  q11_Al   0;                % GPa
Ortho_theta   = [45 -45 -45 45];
  for i = 1:length(Ortho_theta)
     qbar11(i) = (q11.*cosd(Ortho_theta(i)).^4) +...
     (q22.*sind(Ortho_theta(i)).^4) + ...
     (2.*(q12+(2.*q66)).*sind(Ortho_theta(i)).^2.*cosd(Ortho_theta(i)).^2);
     qbar12(i) = (q11 + q22 - ...
     (4.*q66)).*sind(Ortho_theta(i)).^2.*cosd(Ortho_theta(i)).^2 + ...
     (q12.*(cosd(Ortho_theta(i)).^4 + sind(Ortho_theta(i)).^4));
     qbar13(i) = (q11 - q12 - ...
     (2.*q66)).*sind(Ortho_theta(i)).*cosd(Ortho_theta(i)).^3 - ...
     ((q22 - q12 - ...
     (2.*q66)).*(cosd(Ortho_theta(i)).*sind(Ortho_theta(i)).^3));
     qbar21(i) = qbar12(i);
     qbar22(i) = (q11.*sind(Ortho_theta(i)).^4) + ...
     (q22.*cosd(Ortho_theta(i)).^4) + ...
     (2.*(q12 + (2.*q66)).*sind(Ortho_theta(i)).^2.*cosd(Ortho_theta(i)).^2);
      qbar23(i) = (q11 - q12 - ...
     (2.*q66)).*sind(Ortho_theta(i)).^3.*cosd(Ortho_theta(i)) - ...
     ((q22 - q12 - ...
     (2*q66)).*(cosd(Ortho_theta(i)).^3.*sind(Ortho_theta(i))));
      qbar31(i) = qbar13(i);
      qbar32(i) = qbar23(i);
      qbar33(i) = (q11 + q22 - (2*q12) - ...
          (2.*q66)).*sind(Ortho_theta(i)).^2.*cosd(Ortho_theta(i)).^2 + ...
          (q66.*(cosd(Ortho_theta(i)).^4 + sind(Ortho_theta(i)).^4));
      Qbar(1,1,i) = qbar11(i);
      Qbar(1,2,i) = qbar12(i);
      Qbar(1,3,i) = qbar13(i);
      Qbar(2,1,i) = qbar21(i);
      Qbar(2,2,i) = qbar22(i);
      Qbar(2,3,i) = qbar23(i);
      Qbar(3,1,i) = qbar31(i);
      Qbar(3,2,i) = qbar32(i);
      Qbar(3,3,i) = qbar33(i);
  end
  % Calculating Orthotropic A, B, and D Marices
  j = 1;
  dh = [0.01:0.0001:1];
  for i = 1:length(dh)
      Ortho_heights=[-(.7+dh(i)) -.7 0 .7  .7+dh(i)];
      A(:,:,i) = (Ortho_heights(2) - Ortho_heights(1)).*Qbar(:,:,j) + ...
             (Ortho_heights(3) - Ortho_heights(2)).*Qbar(:,:,j+1) + ...
             (Ortho_heights(4) - Ortho_heights(3)).*Qbar(:,:,j+2) + ...
             (Ortho_heights(5) - Ortho_heights(4)).*Qbar(:,:,j+3);
      B(:,:,i) = (1/2)*...
      ((((Ortho_heights(2)).^2 - (Ortho_heights(1)).^2)).*Qbar(:,:,j) + ...
      (((Ortho_heights(3)).^2 - (Ortho_heights(2)).^2)).*Qbar(:,:,j+1) + ...
      (((Ortho_heights(4)).^2 - (Ortho_heights(3)).^2)).*Qbar(:,:,j+2) + ...
      (((Ortho_heights(5)).^2 - (Ortho_heights(4)).^2)).*Qbar(:,:,j+3));
      if B(:,:,i) < 1*10^-6
          B(:,:,i) = 0;
      end
      D(:,:,i) = (1/3)*...
      ((((Ortho_heights(2)).^3 - (Ortho_heights(1)).^3)).*Qbar(:,:,j) + ...
      (((Ortho_heights(3)).^3 - (Ortho_heights(2)).^3)).*Qbar(:,:,j+1) + ...
      (((Ortho_heights(4)).^3 - (Ortho_heights(3)).^3)).*Qbar(:,:,j+2) + ...
      (((Ortho_heights(5)).^3 - (Ortho_heights(4)).^3)).*Qbar(:,:,j+3));
      if ((-(1*10^-6) < D(1,3,i)) && (D(1,3,i) < 1*10^-6))
          flag = 1;
      else
          flag = 0;
      end
      if flag == 1
          breakindex = i;
          break
      end
  end
  i = dh(breakindex)
  fprintf('\nGLARE, [45/-45/-45/45]\n\n');
  fprintf('\nA Matrix\n');
  display(A(:,:,i))
  fprintf('\nB Matrix\n');
  display(B(:,:,i))
  fprintf('\nD Matrix\n');
  display(D(:,:,i))
  ABD_Ortho(:,:,i) = [A(:,:,i) B(:,:,i);
                      B(:,:,i) D(:,:,i)];
  fprintf('------------------------------------------');