%Classical Lamination Theory design code modified for use with %tubing %Note: this thing has plots and pauses hit return play with scaling factors clear all close all %set up a diary file diary CLT.dat %units are US customary (lb, in, E in psi) % total laminate definition in matrix below % [ply angles, thicknesses, matl. #] %Set up for two materials % Data in there now is %1-carbon %2-Eglass %Laminate is defined in this matrix l (sorry it looks like a one) % [ angle thick matl #] l=[ 0 1*0.014 1 +45 1*0.014 1 -45 1*0.014 1 0 1*0.014 1 0 1*0.014 1]; % this is the total laminate % cut, paste, edit above to study your laminate of choice %delta temp disp('Temperature Change [Degrees F]') DT = 0.0 % size command to get number of plies n = size(l,1); % Lamina Properties % matrix for engineering constants %E1 E2 v12 G12 a11 a22 E = [17e6 0.98e6 .31 .42e6 -.5e-6 15e-6; %T700 PrePreg 5.84e6 .9e6 .2 .3e6 0.0e-6 0.0e-6]; %E-Glass/Epoxy % a's are CTE's %intiialize the ply distance and ABD matrices NT = zeros(3,1); MT = zeros(3,1); h = zeros(n+1,1); A = zeros(3); B = zeros(3); D = zeros(3); % Form R matrix which relates engineering to tensor strain R = [1 0 0; 0 1 0; 0 0 2]; %find the total thickness disp('Total Laminate Thickness [inches]') total = sum(l,1); thick = total(1,2) % locate the bottom of the first ply h(1) = -thick/2.; imax = n + 1; %loop for rest of the ply distances from midsurf for i = 2 : imax h(i) = h(i-1) + l(i-1,2); end %loop over each ply to integrate the ABD matrices for i = 1:n %ply material ID mi=l(i,3); v21 = E(mi,2)*E(mi,3)/E(mi,1); d = 1 - E(mi,3)*v21; %Q12 matrix Q = [E(mi,1)/d v21*E(mi,1)/d 0; E(mi,3)*E(mi,2)/d E(mi,2)/d 0; 0 0 E(mi,4)]; %ply angle in radians a1=l(i,1)*pi/180; %Form transformation matrices T1 for ply T1 = [(cos(a1))^2 (sin(a1))^2 2*sin(a1)*cos(a1); (sin(a1))^2 (cos(a1))^2 -2*sin(a1)*cos(a1); -sin(a1)*cos(a1) sin(a1)*cos(a1) (cos(a1))^2-(sin(a1))^2 ]; %Form Qxy Qxy = inv(T1)*Q*R*T1*inv(R); % build up the laminate stiffness matrices A = A + Qxy*(h(i+1)-h(i)); B = B + Qxy*(h(i+1)^2 - h(i)^2); D = D + Qxy*(h(i+1)^3 - h(i)^3); %load alphs into and array a=[E(mi,5); E(mi,6); 0.0]; %transform cte's mult by DT to get thermal strain exy exy = (R*inv(T1)*inv(R)*a)*DT; %build up thermal load as well now NT = NT + Qxy*exy*(h(i+1)-h(i)); MT = MT + .5*(Qxy*exy*(h(i+1)^2 - h(i)^2)); %end of stiffness loop end %change the display format for compliance matrix format short e disp('Laminate Matrices A,B & D') A = 1.0*A B = .5*B D = (1/3)*D %Stiffness Matrix disp('Laminate Stiffness Matrix, K') K = [A, B; B, D] %Compliance Matrix disp('Laminate Compliance Matrix, C (Inverse of K)') C = inv(K) %Input Tube Data Here %Length of the span in inches disp('Length of Tube Span [inches]') L=20 %Tube weight per unit length [lbs/in] disp('Tube Weight per unit Length') W=.0094 %Force to deflect the tube disp('Force on Tube [lb]') F=354.1389996 %Bending Moment disp('Max Bending Moment [in-lbs]') Mb=F*(L/2) %Innerdiameter of the tube tm=0.75; %Distance from the center line to the wall => y=Midradius of the tube!!! Rad=(tm+thick)/2; %Resulting Force in x-Direction disp('X-direction Line Load on Tube [lb/in]') Nx=Mb/(Rad^2*pi) %Resulting Shear Force disp('Shear Line Load on Tube [lb/in]') Nxy=F/(pi*Rad) %Effective EI disp('Effective EI of the Tube') EI=pi*Rad^3/(C(1,1)) %Deflection of the tube Simply Supported disp('Max Deflection of the Tube [Inches]') defl=(F*L^3)/(48*EI) %put in mechanical loads here %mech loads Nx=Nx Ny=0.0 Ns=Nxy Mx=0.0 My=0.0 Ms=0.0 % % superimpose mech and thermal loads load = [ NT(1) + Nx; NT(2) + Ny; NT(3) + Ns; MT(1) + Mx; MT(2) + My; MT(3) + Ms]; %compute the strains = compliance times load disp('Overall Laminate Stains [in/in]') e = C*load %compute the natural frequency in Hz: disp('Natural frequency of tube [Hz]') fn = (pi/2)*((32.2*EI)/(W*L^4))^.5 % calc radii of curvature Rx = 1/e(4); Ry = 1/e(5); Rxy= 1/e(6); %____________________________________________________________ % This stuff below is simply to plot the deformed shape % attempt to plot the deformed element % % % plot scaling, magnification factor A=2.0; % define half of unit element u=.5; % ux = [-u; u; u; -u; -u]; % uy = [-u; -u; u; u; -u]; % dx = [-e(1)*u*A; e(1)*u*A; e(1)*u*A; -e(1)*u*A; -e(1)*u*A]; % % dy = [-e(2)*u*A; -e(2)*u*A; e(2)*u*A; e(2)*u*A; -e(2)*u*A]; % dxy= [-e(3)*u*A; -e(3)*u*A; e(3)*u*A; e(3)*u*A; -e(3)*u*A]; % x1= ux+dx+dxy; y1= uy+dy; % hold on plot(ux,uy, '--b') plot(x1,y1, '-r') axis([-1 1 -1 1]) hold off % pause %_______________________________________________________________________ %_______________________________________________________________________ % Given constant curvatures % plot for ME412 % Kx = e(4); Ky = e(5); Kxy = e(6); %plate W L=6; %Constants of integration C2 = (L/2)*Kx; C3 = (L/2)*Ky; % create X and Y over the domain of the function [X,Y] = meshgrid(0:L,0:L); for i = 1:L+1; for j = 1:L+1; w(i,j) = -.5*(Kx*X(i,j)^2+Ky*Y(i,j)^2+.5*Kxy*X(i,j)*Y(i,j))... + C2*X(i,j) + C3*Y(i,j); end end w; %w(1,L) mesh(w) % %______________________________________________________________________ % Now calc stress and strain and failure index using Max strain % % reduction factor for ultimate (pseudo A-basis use .80) RF=.80; % % % allowable strains reduced to account for ultimate strength after impact % row1 is carbon % row2 is E-glass % transverse prperties assumed same % load allowable strains into array % ELU ELUP ETU ETUP ELTU ea = [RF*.014 RF*.012 RF*.007 RF*.031 RF*.0296; RF*.02 RF*.018 RF*.0067 RF*.031 RF*.0296]; % % %zero out results array ERES = zeros(2*n,6); SRES = zeros(2*n,6); % loop over each ply and calculate strain for i=1 : n; %loop over top and bottom of each ply for j=1 : 2; % one is bottom two is top for loc ply = i; loc = j; z = h(i-1+j); %ply strain from midplane strain el= [ e(1)+z*e(4); e(2)+z*e(5); e(3)+z*e(6)]; %ply material ID mi=l(i,3); v21 = E(mi,2)*E(mi,3)/E(mi,1); d = 1 - E(mi,3)*v21; %Q12 matrix Q = [E(mi,1)/d v21*E(mi,1)/d 0; E(mi,3)*E(mi,2)/d E(mi,2)/d 0; 0 0 E(mi,4)]; % %ply angle in radians a1=l(i,1)*pi/180; %Form transformation matrices T1 for ply T1 = [(cos(a1))^2 (sin(a1))^2 2*sin(a1)*cos(a1); (sin(a1))^2 (cos(a1))^2 -2*sin(a1)*cos(a1); -sin(a1)*cos(a1) sin(a1)*cos(a1) (cos(a1))^2-(sin(a1))^2 ]; % load alpha for the ply a=[E(mi,5); E(mi,6); 0.0]; % tranform to 1,2 % subtract off alpha delta T to get mech strain that causes stress ep = R*T1*inv(R)*el - a*DT; %calculate stress in 1,2 coords sp = Q*ep; %failure index now looks at two different materials if ep(1) > 0.0; FI = ep(1)/ea(mi,1); FIF=FI; elseif ep(1) < 0.0; FI = abs( ep(1) )/ea(mi,2); FIF=FI; end if ep(2) > 0.0; F1 = ep(2)/ea(mi,3); elseif ep(2) < 0.0; F1 = abs( ep(2) )/ea(mi,4); end % if F1 > FI; FI = F1; end % % F1 = abs( ep(3) )/ea(mi,5); if F1 > FI ; FIe = F1; elseif F1 < FI; FIe = FI; end %load the results array % note top and botom of every ply! %strain results, FI based on Max Strain %angle,eps1,eps2,gamma12,FI, FIfiber ERES(2*i+j-2,1)=l(i); ERES(2*i+j-2,2)=ep(1); ERES(2*i+j-2,3)=ep(2); ERES(2*i+j-2,4)=ep(3); ERES(2*i+j-2,5)=FIe; ERES(2*i+j-2,6)=FIF; %stress results, FI based on max strain %angle,Sigma1,Sigma2,Tau12, FI, FIfiber SRES(2*i+j-2,1)=l(i); SRES(2*i+j-2,2)=sp(1); SRES(2*i+j-2,3)=sp(2); SRES(2*i+j-2,4)=sp(3); SRES(2*i+j-2,5)=FIe; SRES(2*i+j-2,6)=FIF; end % end disp('Strains at top and bottom of each ply in the laminate') disp('Ply Angle,Strain Longitudinal,Strain Transverse,Shear Strain, Failure Index, Failure Index Fiber') ERES=ERES*1 disp('Stresses at top and bottom of each ply in the laminate') disp('Ply Angle,Stress Longitudinal,Stress Transverse,Shear Stress, Failure Index, Failure Index Fiber') SRES=SRES*1 diary off % %