Code for Laminated Plate Calculations

A Code for Laminated Plate Calculations The computational scheme for laminated plate calculations outlines in Module 15 has been coded in Fortran as listed below. A PC-executable version is also available in the readings section. . Fortran Source c plate - a prompt-driven routine for laminated plate calculations dimension S(3,3),Sbar(3,3),Qbar(3,3),E(6,7),kwa(6), * T(3,3),Tinv(3,3),R(3,3),Rinv(3,3),Et(6,6), * temp1(3,3),temp2(3,3),temp3(3,3), * eps0(3),xkappa(3),sigbar(3),sig(3),vtemp1(3), * vtemp2(3),E1(10),E2(10),Gnu12(10),G12(10),thk(10), * z(21),mat(20),theta(20),Qsave(3,3,20),Tsave(3,3,20) data R/1.,3*0.,1.,3*0.,2./,Rinv/1.,3*0.,1.,3*0.,.5/, * E/42*0./ c---------------------------------------------------------------------c input material set selections 10 20 i=1 write(6,20) i format (' assign properties for lamina type ',i2,'...'/) write(6,*) 'enter modulus in fiber direction...' write(6,*) ' (enter -1 to stop): ' read (5,*) E1(i) if (E1(i) .lt. 0.) go to 30 write(6,*) 'enter modulus in transverse direction: ' read (5,*) E2(i) write(6,*) 'enter principal Poisson ratio: ' read (5,*) Gnu12(i) c check for isotropy check=abs((E1(i)-E2(i))/E1(i)) if (check.lt.0.001) then G12(i)=E1(i)/(2.*(1.+Gnu12(i))) else write(6,*) 'enter shear modulus: ' read (5,*) G12(i) end if write(6,*) 'enter ply thickness: ' read (5,*) thk(i) i=i+1 go to 10 c---------------------------------------------------------------------c define layup 30 iply=1 z(1)=0. write(6,*) 'define layup sequence, starting at bottom...' write(6,*) ' (use negative material set number to stop)' 40 50 write (6,50) iply format (/' enter material set number for ply number',i3,': ') read (5,*) m if (m.lt.0) go to 60 mat(iply)=m write(6,*) 'enter ply angle: ' read (5,*) theta(iply) z(iply+1)=z(iply)+thk(m) iply=iply+1 go to 40 c compute boundary coordinates (measured from centerline) 60 thick=z(iply) N = iply-1 z0 = thick/2. np1=N+1 70 do 70 i=1,np1 z(i)=z(i)-z0 continue c---------------------------------------------------------------------c---------------------------------------------------------------------c loop over plies, form stiffness matrix do 110 iply=1,N m=mat(iply) c form lamina compliance in 1-2 directions (Eqn. 3.55) S(1,1) = 1./E1(m) S(2,1) = -Gnu12(m) / E1(m) S(3,1) = 0. S(1,2) = S(2,1) S(2,2) = 1./E2(m) S(3,2) = 0. S(1,3) = 0. S(2,3) = 0. S(3,3) = 1./G12(m) c---------------------------------------------------------------------- c c transform to x-y axes obtain transformation matrix T (Eqn. 3.27) thet = theta(iply) * 3.14159/180. sc = sin(thet)*cos(thet) s2 = (sin(thet))**2 c2 = (cos(thet))**2 T(1,1) = c2 T(2,1) = s2 T(3,1) = -1.*sc T(1,2) = s2 T(2,2) = c2 T(3,2) = sc T(1,3) = 2.*sc T(2,3) = -2.*sc T(3,3) = c2 - s2 c inverse transformation matrix Tinv(1,1) = c2 Tinv(2,1) = s2 Tinv(3,1) = sc Tinv(1,2) = s2 Tinv(2,2) = c2 Tinv(3,2) = -1.*sc Tinv(1,3) = -2.*sc Tinv(2,3) = 2.*sc Tinv(3,3) = c2 - s2 c transformation [Sbar] = [R][T]-1[R]-1[S][T] (Eqn. 3.56) call call call call matmul matmul matmul matmul (3,3,3,3,3,3,R,Tinv,temp1) (3,3,3,3,3,3,temp1,Rinv,temp2) (3,3,3,3,3,3,temp2,S,temp3) (3,3,3,3,3,3,temp3,T,Sbar) c---------------------------------------------------------------------c invert Sbar (transformed compliance matrix) to obtain c Qbar (transformed stiffness matrix) c start by setting Qbar = Sbar, then call inversion routine 80 do 80 i=1,3 do 80 j=1,3 Qbar(i,j)=Sbar(i,j) continue call matinv(isol,idsol,3,3,Qbar,3,kwa,det) c save Qbar and Tinv matrices 90 do 90 i=1,3 do 90 j=1,3 Qsave(i,j,iply)=Qbar(i,j) Tsave(i,j,iply)=Tinv(i,j) continue c add to laminate stiffness matrix 100 ip1=iply+1 z1= (z(ip1) -z(iply) ) z2= 0.5*(z(ip1)**2-z(iply)**2) z3=(1./3.)*(z(ip1)**3-z(iply)**3) do 100 i=1,3 do 100 j=1,3 E(i,j) = E(i,j) + Qbar(i,j)*z1 xx = Qbar(i,j)*z2 E(i+3,j) = E(i+3,j) + xx E(i,j+3) = E(i,j+3) + xx E(i+3,j+3)= E(i+3,j+3) + Qbar(i,j)*z3 continue c 110 end loop over plies; stiffness matrix now formed continue c---------------------------------------------------------------------c---------------------------------------------------------------------c output stiffness matrix 120 130 140 write(6,120) format(/' laminate stiffness matrix:',/) do 140 i=1,6 write(6,130) (e(i,j),j=1,6) format (4x,3e12.4,2x,3d12.4) if (i.eq.3) write(6,*) continue c---------------------------------------------------------------------c obtain and print laminate compliance matrix c c c c300 do 300 i=1,6 do 300 j=1,6 Et(i,j)=E(i,j) continue c call matinv(isol,idsol,6,6,Et,6,kwa,det) c c310 c c c c320 write(6,310) format(/' laminate compliance matrix:',/) do 320 i=1,6 write(6,130) (Et(i,j),j=1,6) if (i.eq.3) write(6,*) continue c---------------------------------------------------------------------c obtain traction-moment vector write(6,*) write(6,*) write(6,*) write(6,*) read (5,*) write(6,*) read (5,*) write(6,*) read (5,*) write(6,*) read (5,*) write(6,*) read (5,*) write(6,*) read (5,*) 'input tractions and moments...' ' Nx: e(1,7) ' Ny: e(2,7) ' Nxy: e(3,7) ' Mx: e(4,7) ' My: e(5,7) ' Mxy: e(6,7) ' ' ' ' ' ' c---------------------------------------------------------------------c solve resulting system; print strains and rotations 150 call matinv(isol,idsol,6,7,e,6,kwa,det) write(6,150) (e(i,7),i=1,6) format(/' midplane strains:',//3x,'eps-xx =',e12.4, * /3x,'eps-yy =',e12.4,/3x,'eps-xy =',e12.4, * //' rotations:',//3x,'kappa-xx =',e12.4, * /3x,'kappa-yy= ',e12.4,/3x,'kappa-xy =',e12.4//) c---------------------------------------------------------------------c compute ply stresses 160 write(6,160) format (/' stresses:',/2x,'ply',5x,'sigma-1', * 5x,'sigma-2',4x,'sigma-12'/) do 210 iply=1,N 170 180 do 180 i=1,3 eps0(i)=e(i,7) xkappa(i)=e(i+3,7) do 170 j=1,3 Qbar(i,j)=Qsave(i,j,iply) Tinv(i,j)=Tsave(i,j,iply) continue continue call matmul (3,3,3,3,3,1,Qbar,eps0,vtemp1) call matmul (3,3,3,3,3,1,Qbar,xkappa,vtemp2) 190 zctr=(z(iply)+z(iply+1))/2. do 190 i=1,3 sigbar(i) = vtemp1(i) + zctr*vtemp2(i) continue 200 call matmul (3,3,3,3,3,1,Tinv,sigbar,sig) write(6,200) iply,sig format (3x,i2,3e12.4) 210 continue stop end c---------------------------------------------------------------------c---------------------------------------------------------------------c -------- library routines for matrix operations ----------subroutine matmul(lra,lrb,lrc,i,j,k,a,b,c) c c c c this subroutine performs the multiplication of two two-dimensional matrices (a(i,j)*b(j,k) = c(i,k)). c c c c c c c c c lra lrb lrc i j k a b c 10 20 - row dimension row dimension row dimension actual number actual number actual number multiplier multiplicand product of of of of of of "a" (multiplier) matrix "b" (multiplicand) matrix "c" (product) matrix rows in "a" columns in "a" columns in "b" dimension a(1), b(1), c(1) do 20 l = 1,i nm1 = 0 lm = l do 20 m = 1,k c(lm) = 0.0 nm = nm1 ln = l do 10 n = 1,j nm = nm + 1 c(lm) = c(lm) + a(ln)*b(nm) ln = ln + lra nm1 = nm1 + lrb lm = lm + lrc return end subroutine matinv(isol,idsol,nr,nc,a,mra,kwa,det) c c c c c this subroutine finds the inverse and/or solves simultaneous equations, or neither, and calculates a determinant of a real matrix. c c c c c c c isol - communications flag (output) 1 - inverse found or equations solved 2 - unable to solve 3 - input error idsol - determinant calculation flag (output) 1 - did not overflow 2 - overflow c c c c c c c 10 20 30 40 50 60 70 80 90 100 nr - number of rows in input matrix "a" nc - number of columns in "a" a - input matrix, first "nr" columns will be inverted on output, "a" is converted to a-inverse mra - row dimension of "a" matrix kwa - work array det - value of determinant (if idsol = 1) dimension a(1), kwa(1) ir = nr isol = 1 idsol = 1 if(nr.le.0) go to 330 if((ir-mra).gt.0) go to 330 ic = iabs(nc) if ((ic - ir).lt.0) ic = ir ibmp = 1 jbmp = mra kbmp = jbmp + ibmp nes = ir*jbmp net = ic*jbmp if(nc) 10,330,20 mdiv = jbmp + 1 iric = ir - ic go to 30 mdiv = 1 mad = mdiv mser = 1 kser = ir mz = 1 det = 1.0 piv = 0. i = mser if (( i - kser).gt.0) go to 70 if((abs(a(i))-piv).le.0.) go to 60 piv = abs(a(i)) ip = i i = i + ibmp go to 50 if(piv.eq.0.) go to 340 if(nc.lt.0) go to 80 i = ip-((ip - 1)/jbmp)*jbmp j = mser - ((mser - 1)/jbmp)*jbmp jj = mser/kbmp + 1 ii = jj + (ip -mser) kwa(jj) = ii go to 90 i = ip j = mser if (ip - mser) 330,120,100 if ((j - net).gt.0) go to 110 psto = a(i) a(i) = a(j) a(j) = psto i = i + jbmp j = j + jbmp go to 100 110 120 det = - det psto = a(mser) det = det*psto 130 if (det.eq.0.) goto 150 140 psto = 1./psto go to 160 150 idsol = 1 isol = 2 return 160 continue a(mser) = 1.0 i = mdiv 170 if((i - net).gt.0) go to 180 a(i) = a(i)*psto i = i + jbmp go to 170 180 if((mz - kser).gt.0) go to 210 if((mz-mser).eq.0) go to 200 i = mad j = mdiv psto = a(mz) if(psto.eq.0.) go to 200 a(mz) = 0. 190 if((j-net).gt.0) go to 200 a(i) = a(i) - a(j)*psto j = j + jbmp i = i + jbmp go to 190 200 mad = mad + ibmp mz = mz + ibmp go to 180 210 continue c 210 need a test here.....call overfl(ivf) c go to (350,220),ivf ccccccc need at test here, anyhow 220 kser = kser + jbmp if ((kser-nes).gt.0) go to 260 mser = mser + kbmp if(nc.lt.0) go to 230 mdiv = mdiv + ibmp mz = ((mser - 1)/jbmp)*jbmp + 1 mad = 1 go to 40 230 mdiv = mdiv + kbmp if(iric.ne.0) go to 240 mz = mser + ibmp go to 250 240 mz = ((mser - 1)/jbmp)*jbmp + 1 250 mad = mz + jbmp go to 40 260 if(nc.lt.0) return jr = ir 270 if(jr) 330,360,280 280 if(kwa(jr) - jr) 330,320,290 290 k = (jr - 1)*jbmp j = k + ir l = (kwa(jr) - 1)*jbmp + ir 300 310 320 330 340 350 360 if(j - k) 330,320,310 psto = a(l) a(l) = a(j) a(j) = psto j = j - ibmp l = l - ibmp go to 300 jr = jr - 1 go to 270 isol = 3 return det = 0. isol = 2 idsol = 1 return isol = 2 idsol = 2 return end