program vplate
! F90 version of plate code
! to include voigt effects
! Limits: 10 material sets, 20 plies

! Equation references are to Roylance Module 15

real E_1(10)	! Modulus in fiber direction
real E_2(10)	! Modulus in transverse direction
real G_12(10)	! In-plane shear modulus
real nu_12(10)	! Principal Poisson's ratio
real thk(10)	! Ply thickness
real v_frac(10)	! Creep fraction

integer k		! Ply index
real theta(20)	! Ply orientation (degrees)	

real eps(6)		! vector of strains and curvatures
real eps0(3)	! Laminate in-plane strains
real kappa(3)	! Laminate curvature
real sig(3)		! Ply stresses in 1-2 frame
real sigbar(3)	! Ply stresses in x-y frame
real z(21)		! Vertical coordinate of bottom of ply

real A(3,3)		! Transformation matrix
real A_inv(3,3)	! Inverse transformation matrix
real C(6,6)		! Laminate compliance matrix
real D(3,3)		! Ply stiffness matrix in 1-2 frame
real D_bar(3,3)	! Ply stiffness matrix in x-y frame
real E(6,6)		! Laminate stiffness matrix
real N(6)		! Vector of tractions and moments
real R(3,3)		! Reuter's matrix
real R_inv(3,3)	! Inverse Reuter's matrix
real S(3,3)		! Ply compliance matrix in 1-2 frame
real S_bar(3,3)	! Ply compliance matrix in x-y frame

real eps_l_c(6)		! Laminate creep strain
real eps0_l_c(3)	! Laminate creep strain vector
real kappa_l_c(3)	! Laminate creep curvature vector
real eps_p_e_xy(3)	! Ply elastic strain in x-y frame
real eps_p_e_12(3)	! Ply elastic strain in 1-2 frame
real eps_p_lc_xy(3)	! Ply laminate creep strain in x-y frame
real eps_p_lc_12(3)	! Ply laminate creep strain in 1-2 frame
real eps_pc(3,20)	! Ply creep strain history array
real eps_p_c_12(3)	! Ply creep strain, 1-2 frame
real sig_k_12(3)	! Ply stress, 1-2 frame
real sig_k_xy(3)	! Ply stress, x-y frame
real N_c(3)			! Equivalent laminate creep traction
real M_c(3)			! Equivalent laminate creep bending moment
real N_l_c(6)		! Equivalent laminate creep load

integer m		! material property set index
integer mat(20)	! ply material index
integer Nplies	! number of plies in layup

data R/1.,3*0.,1.,3*0.,2./
data R_inv/1.,3*0.,1.,3*0.,.5/
data E/36*0./

data eps_l_c/6*0./
data eps_pc/60*0./

print*, 'Viscoelastic laminate code, rev 9/11/99'
print*

! ----------------------------------------------------------------------					
! Get material properties

do m = 1,10

  print*, 'Material set ',m,': '
  print*, 'Enter modulus in fiber direction (-1 to stop): '
  read*, E_1(m)
    if (E_1(m) .lt. 0.0 ) exit

  print*, 'Enter modulus in transverse direction: '
  read*, E_2(m)

  print*, 'Enter principal Poisson ratio: '
  read*, nu_12(m)

  check = abs((E_1(m)-E_2(m))/E_1(m) )	   ! check for isotropy
  if (check.lt.0.001) then
    G_12(m) = E_1(m)/(2.*(1.+ nu_12(m) ))  ! calculate G_12 if isotropic
  else
    print*, 'Enter shear modulus: '
	read*, G_12(m)  ! otherwise read G_12
  end if				

  print*, 'Enter ply thickness: '
  read*, thk(m)

  print*, 'Enter creep fraction: '
  read*, v_frac(m)

end do

! ----------------------------------------------------------------------					
! Define layup (ply orientations)

z(1)=0.	  ! vertical coordinate, bottom of first ply

print*, 'Define layup sequence, starting at bottom...'
print*, '   (use negative material set number to stop)'

! loop over plies
do k = 1,20

  write (*,50) k
50  format (/' enter material set number for ply number',i3,': ')

  read (5,*) m
    if (m.lt.0) exit
  mat(k) = m

  print*, 'enter ply angle: '
  read*, theta(k)

  z(k+1)=z(k)+thk(m)
  Nplies = k

end do

! Loop over plies again,
! adjust ply coordinates to measure from centerline

z0 = z(Nplies+1)/2.	! midpoint coordinate

do k = 1, Nplies+1
  z(k) = z(k) - z0
end do


! ----------------------------------------------------------------------					
! Loop over plies,
! compute contributions to stiffness matrix (Eqn. 25, p.9)

do k = 1,Nplies

  z1= z(k+1) - z(k)
  z2= 0.5 * ( z(k+1)**2 - z(k)**2 )
  z3=(1./3.) * ( z(k+1)**3 - z(k)**3)
  call get_D_bar (D_bar)

  do i=1,3
  do j=1,3

    E(i,j)    =  E(i,j) +     D_bar(i,j)*z1
    E(i+3,j)  =  E(i+3,j) +   D_bar(i,j)*z2
    E(i,j+3)  =  E(i,j+3) +   D_bar(i,j)*z2
    E(i+3,j+3)=  E(i+3,j+3) + D_bar(i,j)*z3

  end do
  end do

end do

! ----------------------------------------------------------------------					
! print stiffness matrix E(6,6)

write(*,120)
120 format(/' laminate stiffness matrix:',/)

do i=1,6

  write(*,130) (E(i,j),j=1,6)
130 format (4x,3e12.4,2x,3e12.4)
    if (i.eq.3) print*

end do

! obtain and print laminate compliance matrix C(6,6)

call invert (E,6,6,C)	   ! Numerical recipes inversion routine

write(*,310)
310   format(/' laminate compliance matrix:',/)

do i=1,6

  write(*,130) (C(i,j),j=1,6)
  if (i.eq.3) print*

end do

! ----------------------------------------------------------------------
! obtain traction-moment vector	N(6)

 print*
 print*, 'input tractions and moments...'
 print*

 print*, '   Nx: ' ; read *, N(1)
 print*, '   Ny: ' ; read *, N(2)
 print*, '  Nxy: ' ; read *, N(3)
 print*, '   Mx: ' ; read *, N(4)
 print*, '   My: ' ; read *, N(5)
 print*, '  Mxy: ' ; read *, N(6)

! ----------------------------------------------------------------------
! compute and print resulting strains and rotations

eps = matmul (C,N)

write(6,150) (eps(i),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//)


! ----------------------------------------------------------------------
! compute ply stresses

write(6,160)
160   format (/' stresses:',/2x,'ply',5x,'sigma-1',	 &
        5x,'sigma-2',4x,'sigma-12'/)

do i=1,3  ! get midplane strain and curvature vectors

  eps0(i)  = eps(i)
  kappa(i) = eps(i+3)

end do

do k=1,Nplies	! loop over plies (Eqn. 16, p. 8)

  zctr=(z(k)+z(k+1))/2.		! center coordinate of ply
  call get_D_bar (D_bar)	! this also gets A for current ply

  sigbar = matmul(D_bar,eps0) + zctr * matmul(D_bar,kappa)	! x-y frame
  sig = matmul(A,sigbar)  ! transform to 1-2 frame

  write(6,200) k,sig
  200 format (3x,i2,3e12.4)

end do

! ----------------------------------------------------------------------
! ----------------------------------------------------------------------


! adaptation of Brinson-Tuttle viscoelastic response
! assumes simple voigt behavior in 2- and 12-directions

! get time parameters, begin time loop

! all plies have same relaxation time
print*, 'Enter relaxation time' ; read*, tau

told = 0.
t = tau/100.

do while (t.le.(100.*tau))	! step over +/- 2 decades from tau

dt = t - told

! recover midplane strain and curvature vectors

do i=1,3
   eps0(i) = eps(i)
   kappa(i) = eps(i+3)
   eps0_l_c(i) = eps_l_c(i)
   kappa_l_c(i) = eps_l_c(i+3)
end do

! get creep response of individual plies

! zero-initialize creep load vector

do i=1,3
	N_c(i)=0.
	M_c(i)=0.
end do

do k=1,Nplies

  zctr = (z(k)+z(k+1))/2	! ply center coordinate
  thk_p = z(k+1)-z(k)		! ply thickness

  call get_A(A) ; call get_A_inv(A_inv) ; call get_D(D)

  ! ply elastic strain
  eps_p_e_xy = eps0 + zctr*kappa
  eps_p_e_12 = matmul(R,matmul(A,matmul(R_inv,eps_p_e_xy)))

  ! contribution of laminate creep strain
  eps_p_lc_xy = eps0_l_c + zctr*kappa_l_c
  eps_p_lc_12 = matmul(R,matmul(A,matmul(R_inv,eps_p_lc_xy)))

  ! recover ply creep strain from history
  do i=1,3
    eps_p_c_12(i) = eps_pc(i,k)
  end do

  ! current ply stress
  sig_k_12 = matmul (D,(eps_p_e_12 + eps_p_lc_12 - eps_p_c_12))

  ! independent ply creep strain - voigt model
  ! only 2- and 12-components effected

  m = mat(k)
  C_v_2 =  v_frac(m) / E_2(m)
  C_v_12 = v_frac(m) / G_12(m)

  eps_p_c_12(2) = sig_k_12(2) * C_v_2 * (1. - exp( -dt/tau )) &
				  + eps_p_c_12(2) * exp( -dt/tau )

  eps_p_c_12(3) = sig_k_12(3) * C_v_12 * (1. - exp( -dt/tau )) &
				  + eps_p_c_12(3) * exp( -dt/tau )

  ! store in history for later recovery
  eps_pc(2,k) = eps_p_c_12(2)
  eps_pc(3,k) = eps_p_c_12(3)

  ! stresses corresponding to ply creep strains

  sig_k_12 = matmul (D,eps_p_c_12)
  sig_k_xy = matmul (A_inv,sig_k_12)

  ! contributions to equivalent laminate creep load
  N_c = N_c + thk_p * sig_k_xy
  M_c = M_c + thk_p * zctr * sig_k_xy

end do	! end loop over plies

do i=1,3	! form 6x1 equivalent laminate creep load

  N_l_c(i) = N_c(i)
  N_l_c(i+3) = M_c(i)

end do

! equivalent laminate creep strain
eps_l_c = matmul(C,N_l_c)

write (*,400) t,eps_l_c
!write (*,400) t,eps_p_e_12(2),eps_p_lc_12(2),eps_p_c_12(2), &
!  	 sig_k_12(2),sig_k_xy(1),N_l_c(1),eps_l_c(1)
400 format (7e12.4)

! update time, end time loop

told = t
t = t*(10**(1./4.))	! 4 logarithmic time steps per decade

end do

! ----------------------------------------------------------------------
! ----------------------------------------------------------------------
!
! internal routines
					
contains

   subroutine get_A (A)
   ! Transformation matrix for x-y --> 	1-2	(Eqn. 3.27, p. 99)
   dimension A(3,3)

   thet = theta(k) * 3.14159/180.

    sc = sin(thet)*cos(thet)
    s2 = (sin(thet))**2
    c2 = (cos(thet))**2

    A(1,1) = c2
    A(2,1) = s2
    A(3,1) = -1.*sc

    A(1,2) = s2
    A(2,2) = c2
    A(3,2) = sc

    A(1,3) = 2.*sc
    A(2,3) = -2.*sc
    A(3,3) = c2 - s2
	
	end subroutine get_A
	
	! ----------------------------------------------------------------------					

	subroutine get_A_inv (A_inv)
	! Transformation matrix for 1-2 --> x-y
	dimension A_inv(3,3)

    thet = theta(k) * 3.14159/180.

    sc = sin(thet)*cos(thet)
    s2 = (sin(thet))**2
    c2 = (cos(thet))**2

    A_inv(1,1) = c2
    A_inv(2,1) = s2
    A_inv(3,1) = sc

    A_inv(1,2) = s2
    A_inv(2,2) = c2
    A_inv(3,2) = -1.*sc

    A_inv(1,3) = -2.*sc
    A_inv(2,3) = 2.*sc
    A_inv(3,3) = c2 - s2

	end subroutine get_A_inv
	
! ----------------------------------------------------------------------					

  subroutine get_D(D)
  ! Stiffness matrix for current ply in 1-2 frame
  dimension D(3,3)

  call get_S(S)
  call invert (S,3,3,D)	! Numerical Recipes matrix inversion

  end subroutine get_D

! ----------------------------------------------------------------------					

  subroutine get_D_bar (D_bar)
  ! Stiffness matrix for current ply in x-y frame
  dimension D_bar(3,3)

  ! First, compute compliance matrix S_bar in x-y frame
  ! [S_bar] = [R][A]-1[R]-1[S][A] (Eqn. 3.56, p.114)

  call get_A (A); call get_A_inv (A_inv); call get_S(S)

  S_bar = matmul(R,matmul(A_inv,matmul(R_inv,matmul(S,A))))

  ! invert S_bar to obtain D_bar	

  call invert (S_bar,3,3,D_bar)	! Numerical Recipes matrix inversion

  end subroutine get_D_bar

! ----------------------------------------------------------------------

  subroutine get_S (S)
  ! Compliance matrix for current ply in 1-2 frame	(Eqn. 3.55, p. 119)
  dimension S(3,3)

	m=mat(k)

    S(1,1) = 1./E_1(m)
    S(2,1) = -nu_12(m) / E_1(m)
    S(3,1) = 0.

    S(1,2) = S(2,1)
    S(2,2) = 1./E_2(m)
    S(3,2) = 0.

    S(1,3) = 0.
    S(2,3) = 0.
    S(3,3) = 1./G_12(m)
	
   end subroutine get_S

! ----------------------------------------------------------------------					

end program vplate