{{{id=6| @interact def Input_Values(t=input_box(0.000125, label="Ply Thickness (m)",type=RDF), E11=input_box(170*10^9, label="E11 (Pa)", type=RDF), E22=input_box(10.0*10^9, label="E22 (Pa)"), G=input_box(13*10^9, label="G (Pa)",type=RDF), v=input_box(0.3, label="Poisson's Ration v",type=RDF), a11=input_box(-0.9*10^-6, label="Coeff Thermal exp, 11 (m/m K)", type=RDF), a22= input_box(27*10^-6, label="Coeff Thermal exp, 22 (m/m K)",type=RDF), b11=input_box(50*10^-6, label= "Coeff. Moisture exp, 11 (m/m %moisture)",type=RDF), b22= input_box(1200*10^-6, label= "Coeff. Moisture exp, 11 (m/m %moisture)",type=RDF), Nxx=input_box(0.0, label="Normal Force xx (N)",type=RDF), Nyy=input_box(0.0,label="Normal Force yy (N)",type=RDF), Nxy=input_box(0.0,label="Normal Shear Force (N)",type=RDF), Mxx=input_box(0.0,label="Moment xx (N-m)",type=RDF), Myy=input_box(0.0,label="Moment yy (N-m)",type=RDF), Mxy=input_box(0.0,label="Shear Moment xy (N-m)",type=RDF), dT=input_box(0.0, label="Change in temperature (C)",type=RDF), dM=input_box(0.0, label="Change in Moisture content (% moisture)",type=RDF), theta=input_box([0.0,45.0,90.0,-45.0,0.0], label="Ply angles in sequence (bottom to top)",type=list), sigma11fT=input_box(170.0*10^9, label="Max stress 11 (Pa)",type=RDF), sigma22yT=input_box(50.0*10^9, label="Max stress 22 (Pa)",type=RDF), tau12y=input_box(75.0*10^6, label="Max shear stress(Pa)",type=RDF)): N=[Nxx,Nyy,Nxy] M=[Mxx,Myy,Mxy] clt(t,E11,E22,v,G,sigma11fT,sigma22yT,tau12y,theta,a11,a22,b11,b22,dT,dM,N,M) def clt(t,E11,E22,v12,G12,sigma11fT,sigma22yT,tau12y,theta,a11,a22,b11,b22,dT,dM,N,M): for i in range(len(theta)): theta[i] = theta[i]*RDF(pi)/180 n = len(theta) z = range(n+1) z[0] = -n/2*t for i in range(1,n+1): z[i] = z[0]+i*t za = range(2*n) za[0] = round(z[0],6) za[2*n-1] = round(z[n],6) for i in range(2,(2*n),2): za[i] = round(z[int(i/2)],6) za[i-1] = round(z[int(i/2)],6) # Calculate the Stiffness matrix Q11 = E11^2/(E11-v12^2*E22) Q12 = v12*E11*E22/(E11-v12^2*E22) Q22 = E11*E22/(E11 - v12^2*E22) Q66 = G12 # Calculate transformed stiffness matrix components and thermal and moisture coefficients for each ply Qt11 = range(n) Qt12 = range(n) Qt16 = range(n) Qt22 = range(n) Qt26 = range(n) Qt66 = range(n) axx = range(n) ayy = range(n) axy = range(n) bxx = range(n) byy = range(n) bxy = range(n) for i in range(n): Qt11[i] = Q11*cos(theta[i])^4 + 2*(Q12 + 2*Q66)*cos(theta[i])^2*sin(theta[i])^2+Q22*sin(theta[i])^4 Qt12[i] = Q12*(cos(theta[i])^4 + sin(theta[i])^4)+(Q11+Q22-4*Q66)*cos(theta[i])^2*sin(theta[i])^2 Qt16[i] = (Q11 - Q12 - 2*Q66)*cos(theta[i])^3*sin(theta[i])-(Q22-Q12-2*Q66)*cos(theta[i])*sin(theta[i])^3 Qt22[i] = Q11*sin(theta[i])^4+2*(Q12+2*Q66)*cos(theta[i])^2*sin(theta[i])^2+Q22*cos(theta[i])^4 Qt26[i] = (Q11 - Q12 - 2*Q66)*cos(theta[i])*sin(theta[i])^3 - (Q22 - Q12 - 2*Q66)*cos(theta[i])^3*sin(theta[i]) Qt66[i] = (Q11 + Q22 -2*Q12 - 2*Q66)*cos(theta[i])^2*sin(theta[i])^2 + Q66*(cos(theta[i])^4 + sin(theta[i])^4) axx[i] = a11*cos(theta[i])^2 + a22*sin(theta[i])^2 ayy[i] = a11*sin(theta[i])^2 + a22*cos(theta[i])^2 axy[i] = 2*cos(theta[i])*sin(theta[i])*(a11 - a22) bxx[i] = b11*cos(theta[i])^2 + b22*sin(theta[i])^2 byy[i] = b11*sin(theta[i])^2 + b22*cos(theta[i])^2 bxy[i] = 2*cos(theta[i])*sin(theta[i])*(b11 - b22) Ab11 = range(n) Bb11 = range(n) Db11 = range(n) Ab12 = range(n) Bb12 = range(n) Db12 = range(n) Ab22 = range(n) Bb22 = range(n) Db22 = range(n) Ab16 = range(n) Bb16 = range(n) Db16 = range(n) Ab26 = range(n) Bb26 = range(n) Db26 = range(n) Ab66 = range(n) Bb66 = range(n) Db66 = range(n) Nt_xx = range(n) Nt_yy = range(n) Nt_xy = range(n) Mt_xx = range(n) Mt_yy = range(n) Mt_xy = range(n) Nm_xx = range(n) Nm_yy = range(n) Nm_xy = range(n) Mm_xx = range(n) Mm_yy = range(n) Mm_xy = range(n) for k in range(n): Ab11[k] = Qt11[k]*(z[k+1]-z[k]) Ab12[k] = Qt12[k]*(z[k+1]-z[k]) Ab22[k] = Qt22[k]*(z[k+1]-z[k]) Ab16[k] = Qt16[k]*(z[k+1]-z[k]) Ab26[k] = Qt26[k]*(z[k+1]-z[k]) Ab66[k] = Qt66[k]*(z[k+1]-z[k]) Bb11[k] = 0.5*Qt11[k]*(z[k+1]^2 - z[k]^2) Bb12[k] = 0.5*Qt12[k]*(z[k+1]^2 - z[k]^2) Bb22[k] = 0.5*Qt22[k]*(z[k+1]^2 - z[k]^2) Bb16[k] = 0.5*Qt16[k]*(z[k+1]^2 - z[k]^2) Bb26[k] = 0.5*Qt26[k]*(z[k+1]^2 - z[k]^2) Bb66[k] = 0.5*Qt66[k]*(z[k+1]^2 - z[k]^2) Db11[k] = 1/3*Qt11[k]*(z[k+1]^3 - z[k]^3) Db12[k] = 1/3*Qt12[k]*(z[k+1]^3 - z[k]^3) Db22[k] = 1/3*Qt22[k]*(z[k+1]^3 - z[k]^3) Db16[k] = 1/3*Qt16[k]*(z[k+1]^3 - z[k]^3) Db26[k] = 1/3*Qt26[k]*(z[k+1]^3 - z[k]^3) Db66[k] = 1/3*Qt66[k]*(z[k+1]^3 - z[k]^3) Nt_xx[k] = dT*(Qt11[k]*axx[k]+Qt12[k]*ayy[k]+Qt16[k]*axy[k])*(z[k+1]-z[k]) Nt_yy[k] = dT*(Qt12[k]*axx[k]+Qt22[k]*ayy[k]+Qt26[k]*axy[k])*(z[k+1]-z[k]) Nt_xy[k] = dT*(Qt16[k]*axx[k]+Qt26[k]*ayy[k]+Qt66[k]*axy[k])*(z[k+1]-z[k]) Mt_xx[k] = dT*0.5*(Qt11[k]*axx[k]+Qt12[k]*ayy[k]+Qt16[k]*axy[k])*(z[k+1]^2-z[k]^2) Mt_yy[k] = dT*0.5*(Qt12[k]*axx[k]+Qt22[k]*ayy[k]+Qt26[k]*axy[k])*(z[k+1]^2-z[k]^2) Mt_xy[k] = dT*0.5*(Qt16[k]*axx[k]+Qt26[k]*ayy[k]+Qt66[k]*axy[k])*(z[k+1]^2-z[k]^2) Nm_xx[k] = dM*(Qt11[k]*bxx[k]+Qt12[k]*byy[k]+Qt16[k]*bxy[k])*(z[k+1]-z[k]) Nm_yy[k] = dM*(Qt12[k]*bxx[k]+Qt22[k]*byy[k]+Qt26[k]*bxy[k])*(z[k+1]-z[k]) Nm_xy[k] = dM*(Qt16[k]*bxx[k]+Qt26[k]*byy[k]+Qt66[k]*bxy[k])*(z[k+1]-z[k]) Mm_xx[k] = dM*0.5*(Qt11[k]*bxx[k]+Qt12[k]*byy[k]+Qt16[k]*bxy[k])*(z[k+1]^2-z[k]^2) Mm_yy[k] = dM*0.5*(Qt12[k]*bxx[k]+Qt22[k]*byy[k]+Qt26[k]*bxy[k])*(z[k+1]^2-z[k]^2) Mm_xy[k] = dM*0.5*(Qt16[k]*bxx[k]+Qt26[k]*byy[k]+Qt66[k]*bxy[k])*(z[k+1]^2-z[k]^2) Ntxx = sum(Nt_xx) Ntyy = sum(Nt_yy) Ntxy = sum(Nt_xy) Mtxx = sum(Mt_xx) Mtyy = sum(Mt_yy) Mtxy = sum(Mt_xy) Nmxx = sum(Nm_xx) Nmyy = sum(Nm_yy) Nmxy = sum(Nm_xy) Mmxx = sum(Mm_xx) Mmyy = sum(Mm_yy) Mmxy = sum(Mm_xy) A11 = sum(Ab11) A12 = sum(Ab12) A22 = sum(Ab22) A16 = sum(Ab16) A26 = sum(Ab26) A66 = sum(Ab66) B11 = sum(Bb11) B12 = sum(Bb12) B22 = sum(Bb22) B16 = sum(Bb16) B26 = sum(Bb26) B66 = sum(Bb66) D11 = sum(Db11) D12 = sum(Db12) D22 = sum(Db22) D16 = sum(Db16) D26 = sum(Db26) D66 = sum(Db66) A = matrix(RDF,3,3,[A11,A12,A16,A12,A22,A26,A16,A26,A66]) B = matrix(RDF,3,3,[B11,B12,B16,B12,B22,B26,B16,B26,B66]) D = matrix(RDF,3,3,[D11,D12,D16,D12,D22,D26,D16,D26,D66]) ABD = block_matrix([A,B,B,D],subdivide=False) abd = (ABD)^-1 L = matrix(RDF,6,1,[(N[0]+Ntxx+Nmxx),(N[1]+Ntyy+Nmyy),(N[2]+Ntxy+Nmxy),(M[0]+Mtxx+Mmxx),(M[1]+Mtyy+Mmyy),(M[2]+Mtxy+Mmxy)]) strain = abd*L exx = range(len(za)) eyy = range(len(za)) exy = range(len(za)) #Calculating Strains for i in range(len(za)): exx[i] = (strain[0,0]+za[i]*strain[3,0]) eyy[i] = (strain[1,0]+za[i]*strain[4,0]) exy[i] = (strain[2,0]+za[i]*strain[5,0]) #Calculating Stresses sxx = range(len(za)) syy = range(len(za)) sxy = range(len(za)) s11 = range(len(za)) s22 = range(len(za)) s12 = range(len(za)) for i in range(0,2*n,2): sxx[i] = (Qt11[i/2]*(exx[i]-dT*axx[i/2]-dM*bxx[i/2])+Qt12[i/2]*(eyy[i]-dT*ayy[i/2]-dM*byy[i/2])+Qt16[i/2]*(exy[i]-dT*axy[i/2]-dM*bxy[i/2])) sxx[i+1] = (Qt11[i/2]*(exx[i+1]-dT*axx[i/2]-dM*bxx[i/2])+Qt12[i/2]*(eyy[i+1]-dT*ayy[i/2]-dM*byy[i/2])+Qt16[i/2]*(exy[i+1]-dT*axy[i/2]-dM*bxy[i/2])) syy[i] = (Qt12[i/2]*(exx[i]-dT*axx[i/2]-dM*bxx[i/2])+Qt22[i/2]*(eyy[i]-dT*ayy[i/2]-dM*byy[i/2])+Qt26[i/2]*(exy[i]-dT*axy[i/2]-dM*bxy[i/2])) syy[i+1] = (Qt12[i/2]*(exx[i+1]-dT*axx[i/2]-dM*bxx[i/2])+Qt22[i/2]*(eyy[i+1]-dT*ayy[i/2]-dM*byy[i/2])+Qt26[i/2]*(exy[i+1]-dT*axy[i/2]-dM*bxy[i/2])) sxy[i] = (Qt16[i/2]*(exx[i]-dT*axx[i/2]-dM*bxx[i/2])+Qt26[i/2]*(eyy[i]-dT*ayy[i/2]-dM*byy[i/2])+Qt66[i/2]*(exy[i]-dT*axy[i/2]-dM*bxy[i/2])) sxy[i+1] = (Qt16[i/2]*(exx[i+1]-dT*axx[i/2]-dM*bxx[i/2])+Qt26[i/2]*(eyy[i+1]-dT*ayy[i/2]-dM*byy[i/2])+Qt66[i/2]*(exy[i+1]-dT*axy[i/2]-dM*bxy[i/2])) s11[i] = cos(theta[i/2])^2*sxx[i] + sin(theta[i/2])^2*syy[i] + 2*cos(theta[i/2])*sin(theta[i/2])*sxy[i] s11[i+1] = cos(theta[i/2])^2*sxx[i+1] + sin(theta[i/2])^2*syy[i+1] + 2*cos(theta[i/2])*sin(theta[i/2])*sxy[i+1] s22[i] = sin(theta[i/2])^2*sxx[i] + cos(theta[i/2])^2*syy[i] - 2*cos(theta[i/2])*sin(theta[i/2])*sxy[i] s22[i+1] = sin(theta[i/2])^2*sxx[i+1] + cos(theta[i/2])^2*syy[i+1] - 2*cos(theta[i/2])*sin(theta[i/2])*sxy[i+1] s12[i] = -cos(theta[i/2])*sin(theta[i/2])*sxx[i]+ cos(theta[i/2])*sin(theta[i/2])*syy[i] + (cos(theta[i/2])^2-sin(theta[i/2])^2)*sxy[i] s12[i+1] = -cos(theta[i/2])*sin(theta[i/2])*sxx[i+1]+ cos(theta[i/2])*sin(theta[i/2])*syy[i+1] + (cos(theta[i/2])^2-sin(theta[i/2])^2)*sxy[i+1] output(sxx,syy,sxy,theta,za) failure(s11,s22,s12,za,sigma11fT,sigma22yT,tau12y) def output(sxx,syy,sxy,theta,za): str_sxx=[] str_syy=[] str_sxy=[] str_za =[] woutx=[] wouty=[] woutxy=[] woutza=[] for i in range(len(za)): str_sxx.append(str(sxx[i])) str_syy.append(str(syy[i])) str_sxy.append(str(syy[i])) str_za.append(str(za[i])) woutx.append(len(str_sxx[i])) wouty.append(len(str_syy[i])) woutxy.append(len(str_sxy[i])) woutza.append(len(str_za[i])) w_outx = max(woutx) w_outy = max(wouty) w_outxy = max(woutxy) w_outza = max(woutza) print "" print "| Ply# | Angle | z (m)"+" "*(w_outza-4)+"| Stress x (Pa)"+" "*(w_outx-12)+"| Stress y (Pa)"+" "*(w_outy-12)+"| Stress xy (Pa)"+" "*(w_outxy-13)+"|" print "+"+"-"*6+ "+"+"-"*7+ "+"+"-"*(w_outza+2)+ "+"+"-"*(w_outx+2)+ "+"+"-"*(w_outy+2)+ "+"+"-"*(w_outxy+2)+"+" for i in range(len(za)): if i%2==0: print "|"+ str((i/2)+1).rjust(6, " ")+"|"+str(round((theta[i/2]*180/RDF(pi)))).rjust(7," ")+"|", else: print "|"+" "*6+"|"+" "*7+"|", print str_za[i].ljust(w_outza+1, " ")+"| "+str_sxx[i].ljust(w_outx," ")+" | "+str_syy[i].ljust(w_outy," ")+" | "+str_sxy[i].ljust(w_outxy," ")+" |" if i%2!=0: print "+"+"-"*6+ "+"+"-"*7+ "+"+"-"*(w_outza+2)+ "+"+"-"*(w_outx+2)+ "+"+"-"*(w_outy+2)+ "+"+"-"*(w_outxy+2)+"+" print "" def failure(s11,s22,s12,za,s11fT,s22yT,s12y): #Applies Tsai-Hill Failure Criterion for i in range(len(za)): f = s11[i]^2/s11fT^2 + s22[i]^2/s22yT^2+s12[i]^2/s12y^2-s11[i]*s22[i]/s11fT^2 if f >= 1: print "Ply # %s failed at z= %s"%(int(i/2)+1,za[i]) /// }}} {{{id=30| /// }}}