Code for Laminated Plate Calculations

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): '
if (E1(i) .lt. 0.) go to 30
write(6,*) 'enter modulus in transverse direction: '
write(6,*) 'enter principal Poisson ratio: '
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: '
end if
write(6,*) 'enter ply thickness: '
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,': ')
if (m.lt.0) go to 60
mat(iply)=m
write(6,*) 'enter ply angle: '
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: '
write(6,*) '   Ny: '
write(6,*) '  Nxy: '
write(6,*) '   Mx: '
write(6,*) '   My: '
write(6,*) '  Mxy: '
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
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
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
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
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
```