%%Program for Calculation of normal strain and transverse deflection of a compsoite beam in three-point bending configuration. clc; close all; clear all; digits(32); format short e; %---------------------------------------------------------------------- %enter layer stacking here! l=[0 0]; %enter the mechanical properties here in SI format t = 0.00136; E1=41.8E+9; G12=1.63E+09; E2=14E+9; E3=E2; n21=0.28; n31=n21; n23=.428/.254*n21; G13=G12; G23=3.28/4.7*G12; tf=0.35e-3; Ef=26e9; nf=0.3; Gf=1.2e9; G12f=Gf; tc=3.2e-3; Gc=.022e9; b1=(.03); bw=48.2e-3; y1=.05/2; l1=.5334; tb=0.5e-3; Eb=2.5e9; G12b=24e6; %enter the force vector here! forcev=[778.5 667.23 444.82 222.41]; %--------------------------------------------------------------------- fss=size(forcev); for R=1:fss(1,2) q1(R)= forcev(1,R); all_ang = l; k = {}; N=length(all_ang); h=zeros(N+1); h(N+1)=N*t/2; lstf=zeros(3*N,3); for k = 1 : N %%Compliance Matrix elements S = [1/E1 -n21/E1 -n31/E1 0 0 0; -n21/E1 1/E2 -n23/E2 0 0 0; -n31/E1 -n23/E2 1/E3 0 0 0; 0 0 0 1/G23 0 0; 0 0 0 0 1/G13 0; 0 0 0 0 0 1/G12]; %%Stiffness Matrix elements C = inv(S); %Reduced Stiffness Q = zeros(3,3); Q(1,1) = C(1,1) - (C(1,3))^2 / C(3,3); Q(2,2) = C(2,2) - (C(2,3))^2 / C(3,3); Q(1,2) = C(1,2) - C(1,3)*C(2,3) / C(3,3); Q(2,1) = Q(1,2); Q(3,3) = C(6,6); %Transformatioin Matrix m = cosd (all_ang(k)); n = sind (all_ang(k)); T = zeros(3,3); T(1,1)=m^2; T(1,2)=n^2; T(1,3)=2*m*n; T(2,1)=n^2; T(2,2)=m^2; T(2,3)=-2*m*n; T(3,1)=-m*n; T(3,2)=m*n; T(3,3)=m^2-n^2; %Transformed Stiffness q=zeros(3,3); q([1:2],[1:2])=Q([1:2],[1:2]); q(3,3)=2*Q(3,3); q=inv(T)*q*T; for i = 1:3 q(i,3)=q(i,3)/2; end lstf([3*k-2:3*k],[1:3])=q([1:3],[1:3]); h(k)=(k-N/2-1)*t; end % To Calculate A B and D Matrices for Uniform Laminate A=zeros(3,3); B=zeros(3,3); D=zeros(3,3); for i=1:3 for j=1:3 q([1:3],[1:3])=lstf([1:3],[1:3]); q(1,1)=q(1,1)-q(1,2)^2./q(2,2); q(1,3)=q(1,3)-q(1,2)*q(2,3)./q(2,2); q(3,1)=q(1,3); q(3,3)=q(3,3)-q(2,3)^2./q(2,2); A(i,j) = q(i,j) * (h(2) - h(1)); B(i,j) = 1/2*(q(i,j) * (h(2)^2 - h(1)^2)); D(i,j) = 1/3*(q(i,j) * (h(2)^3 - h(1)^3)); for k = 2 : N q([1:3],[1:3]) = lstf( [3*k-2:3*k] , [1:3] ); q(1,1)=q(1,1)-q(1,2)^2./q(2,2); q(1,3)=q(1,3)-q(1,2)*q(2,3)./q(2,2); q(3,1)=q(1,3); q(3,3)=q(3,3)-q(2,3)^2./q(2,2); A(i,j) = q(i,j) * (h(k+1) - h(k)) + A(i,j); B(i,j) = 1/2*(q(i,j) * (h(k+1)^2 - h(k)^2)) + B(i,j); D(i,j) = 1/3*(q(i,j) * (h(k+1)^3 - h(k)^3)) + D(i,j); end end end LamnStf=zeros(6,6); a=zeros(3,3); b=zeros(3,3); d=zeros(3,3); LamnStf([1:3],[1:3])=A([1:3],[1:3]); LamnStf([4:6],[4:6])=D([1:3],[1:3]); LamnStf([1:3],[4:6])=B([1:3],[1:3]); LamnStf([4:6],[1:3])=B([1:3],[1:3]); LamnCmp=inv(LamnStf); a([1:3],[1:3])=LamnCmp([1:3],[1:3]); b([1:3],[1:3])=LamnCmp([1:3],[4:6]); d([1:3],[1:3])=LamnCmp([4:6],[4:6]); if rem(length(all_ang),2)==0 A= A.*[1 1 0 ; 1 1 0 ; 0 0 1]; B= B.*[0 0 0 ; 0 0 0 ; 0 0 0]; D= D.*[1 1 1 ; 1 1 1 ; 1 1 1]; a= a.*[1 1 0 ; 1 1 0 ; 0 0 1]; b= b.*[0 0 0 ; 0 0 0 ; 0 0 0]; d= d.*[1 1 1 ; 1 1 1 ; 1 1 1]; end A; B; D; q55=G13; q44=G23; A55=0; A45=0; A44=0; for i=1:N q44b(i)=q44*(cosd(l(i)))^2 + q55*(sind(l(i)))^2; q45b(i)=(q55-q44)*(cosd(l(i)))*(sind(l(i))); q55b(i)=q55*(cosd(l(i)))^2 + q44*(sind(l(i)))^2; A44=A44+q44b(i)*t; A45=A45+q45b(i)*t; A55=A55+q55b(i)*t; end A55s=A44/(A44*A55-A45*A45); Gxzf=1/(A55s*t*length(l)); A55n=1/A55s; A11w=2*Ef*tf/(1-nf^2); A66w=2*G12f*tf; A55w=Gc*tc+2*tf*Gf; A55f=G12*t*length(l)+G12b*tb*length(l); E33w(R)=bw.^3.*A11w./12; E33(R)=(2.*A(1,1).*y1.^2 -4.*B(1,1).*y1 +3.*D(1,1)).*b1 + (2.*Eb*tb.*y1.^2).*b1./3+ bw.^3.*A11w./12; E77(R)= 3*b1.*A55 + A(3,3).*b1; E77s(R)=2*b1./A55s+bw.*A55w+A66w.*bw; %A55=1/A55s E77s1(R)=bw.*A55w+A66w.*bw; %full field strain np=20; epsxxl(R)=((q1(R)/2).*y1)./E33(R); xx=linspace(0,l1./2,np); epsxx(R,:)=epsxxl(1,R).*xx; xxo=[xx l1./2+xx]; epsxxo(R,:)=[epsxx(R,:) fliplr(epsxx(R,:))]; %full field deflection for r=1:np vmzzb(R,r)=((q1(R)*l1.^3)./(48.*E33(R))).*(3.*(xx(r)./l1)-4.*(xx(r)./l1).^3); vmzzs(R,r)=((q1(R)*l1)./(2.*E77s(R))).*(xx(r)./l1); vmzzcom(R,r)=(vmzzb(R,r)./vmzzs(R,r)); vmzz(R,r)=((q1(R)*l1.^3)./(48.*E33(R))).*(3.*(xx(r)./l1)-4.*(xx(r)./l1).^3)+ ((q1(R)*l1)./(2.*E77s(R))).*(xx(r)./l1); end vmzzo(R,:)=[vmzz(R,:) fliplr(vmzz(R,:))]; end figure (1) fii=plot(xxo,1e6*epsxxo); xlabel('Length(m)'); ylabel('{\epsilon_{11}}(\mum/m)'); axis([0 .58 0 750]) set(gca,'YLim',[0 750]) set(gca,'YTick',(0:100:750)) figure (2) plot(xxo,-1e3*vmzzo) fii=plot(xxo,-1e3*vmzzo); xlabel('Length(m)'); ylabel('{\nu} (mm)'); axis([0 .58 -1 0]) set(gca,'YLim',[-1 0]) set(gca,'YTick',(-1:.1:0))