The .exe file found in this section is an executable program and can be run from your computer after download is complete.
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.
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 i=1 10 write(6,20) i 20 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 write (6,50) iply 50 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 do 70 i=1,np1 z(i)=z(i)-z0 70 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 transform to x-y axes c 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 matmul (3,3,3,3,3,3,R,Tinv,temp1) call matmul (3,3,3,3,3,3,temp1,Rinv,temp2) call matmul (3,3,3,3,3,3,temp2,S,temp3) call matmul (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 do 80 i=1,3 do 80 j=1,3 Qbar(i,j)=Sbar(i,j) 80 continue call matinv(isol,idsol,3,3,Qbar,3,kwa,det) c save Qbar and Tinv matrices do 90 i=1,3 do 90 j=1,3 Qsave(i,j,iply)=Qbar(i,j) Tsave(i,j,iply)=Tinv(i,j) 90 continue c add to laminate stiffness matrix 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 100 continue c end loop over plies; stiffness matrix now formed 110 continue c———————————————————————- c———————————————————————- c output stiffness matrix write(6,120) 120 format(/’ laminate stiffness matrix:’,/) do 140 i=1,6 write(6,130) (e(i,j),j=1,6) 130 format (4x,3e12.4,2x,3d12.4) if (i.eq.3) write(6,*) 140 continue c———————————————————————- c obtain and print laminate compliance matrix c do 300 i=1,6 c do 300 j=1,6 c Et(i,j)=E(i,j) c300 continue c call matinv(isol,idsol,6,6,Et,6,kwa,det) c write(6,310) c310 format(/’ laminate compliance matrix:’,/) c do 320 i=1,6 c write(6,130) (Et(i,j),j=1,6) c if (i.eq.3) write(6,*) c320 continue c———————————————————————- c obtain traction-moment vector write(6,*) write(6,*) ‘input tractions and moments…’ write(6,*) write(6,*) ’ Nx: ' read (5,*) e(1,7) write(6,*) ’ Ny: ' read (5,*) e(2,7) write(6,*) ’ Nxy: ' read (5,*) e(3,7) write(6,*) ’ Mx: ' read (5,*) e(4,7) write(6,*) ’ My: ' read (5,*) e(5,7) write(6,*) ’ Mxy: ' read (5,*) e(6,7) c———————————————————————- c solve resulting system; print strains and rotations call matinv(isol,idsol,6,7,e,6,kwa,det) write(6,150) (e(i,7),i=1,6) 150 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 write(6,160) 160 format (/’ stresses:’,/2x,‘ply’,5x,‘sigma-1’, * 5x,‘sigma-2’,4x,‘sigma-12’/) do 210 iply=1,N 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) 170 continue 180 continue call matmul (3,3,3,3,3,1,Qbar,eps0,vtemp1) call matmul (3,3,3,3,3,1,Qbar,xkappa,vtemp2) zctr=(z(iply)+z(iply+1))/2. do 190 i=1,3 sigbar(i) = vtemp1(i) + zctr*vtemp2(i) 190 continue call matmul (3,3,3,3,3,1,Tinv,sigbar,sig) write(6,200) iply,sig 200 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 this subroutine performs the multiplication of two c two-dimensional matrices (a(i,j)*b(j,k) = c(i,k)). c c lra - row dimension of “a” (multiplier) matrix c lrb - row dimension of “b” (multiplicand) matrix c lrc - row dimension of “c” (product) matrix c i - actual number of rows in “a” c j - actual number of columns in “a” c k - actual number of columns in “b” c a - multiplier c b - multiplicand c c - product 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) 10 ln = ln + lra nm1 = nm1 + lrb 20 lm = lm + lrc return end subroutine matinv(isol,idsol,nr,nc,a,mra,kwa,det) c c this subroutine finds the inverse and/or solves c simultaneous equations, or neither, and c calculates a determinant of a real matrix. c c isol - communications flag (output) c 1 - inverse found or equations solved c 2 - unable to solve c 3 - input error c idsol - determinant calculation flag (output) c 1 - did not overflow c 2 - overflow c nr - number of rows in input matrix “a” c nc - number of columns in “a” c a - input matrix, first “nr” columns will be inverted c on output, “a” is converted to a-inverse c mra - row dimension of “a” matrix c kwa - work array c 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 10 mdiv = jbmp + 1 iric = ir - ic go to 30 20 mdiv = 1 30 mad = mdiv mser = 1 kser = ir mz = 1 det = 1.0 40 piv = 0. i = mser 50 if (( i - kser).gt.0) go to 70 if((abs(a(i))-piv).le.0.) go to 60 piv = abs(a(i)) ip = i 60 i = i + ibmp go to 50 70 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 80 i = ip j = mser 90 if (ip - mser) 330,120,100 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 det = - det 120 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 if(j - k) 330,320,310 310 psto = a(l) a(l) = a(j) a(j) = psto j = j - ibmp l = l - ibmp go to 300 320 jr = jr - 1 go to 270 330 isol = 3 return 340 det = 0. isol = 2 idsol = 1 return 350 isol = 2 idsol = 2 360 return end
