c$Id:$ subroutine shell3d(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) c * * F E A P * * A Finite Element Analysis Program c.... Copyright (c) 1984-2005: Robert L. Taylor c All rights reserved c-----[--.----+----.----+----.-----------------------------------------] c Purpose: Quadrilateral shell element for feap c Input parameters set as follows: c Small deformation c ndm = 3 (x,y,z cartesian coordinates at nodes) c ndf = 6 (u-x,u-y,u-z,r-x,r-y,r-z at nodes) c nen = 4 nodes (counterclockwise around element) c Note: 1-direction bisects diagonals between 2-3 element and c 2-direction bisects diagonals between 3-4 element nodes. c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none include 'eldata.h' include 'iofile.h' include 'hdata.h' include 'mdata.h' include 'strnum.h' include 'comblk.h' integer ndm,ndf,nst,isw, tdof, i integer ix(*) real*8 d(*),xl(ndm,*),ul(ndf,*),s(nst,*),p(nst) save c Input material properties if(isw.eq.1) then if(ior.lt.0) write(*,2000) write(iow,2000) call inmate(d,tdof, 24 ,5) c Deactivate dof in element for dof > 6 do i = 7,ndf ix(i) = 0 end do c History for through thickness integrations if(nint(d(102)).gt.1) then nh1 = nh1*nint(d(102)) endif c Construct rotation parameters: u-x = 1; u-y = 2 (same as defaults) ea(1,-iel) = 1 ea(2,-iel) = 2 c Construct rotation parameters: theta-x = 4; theta-y = 5 er(1,-iel) = 4 er(2,-iel) = 5 c Set plot sequence call plqud4(iel) c Set maximum number of stress plots istv = 24 C Remaining options else c Set zero state c Small deformation if(d(18).gt.0.0d0) then if(nint(d(102)).gt.1) then call shl3di(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) else call shl3ds(d,ul,xl,ix,s,p,ndf,ndm,nst,isw) endif c Finite deformation else write(*,*) ' *ERROR* No finite deformation shell' endif endif c Format 2000 format(5x,'T h r e e D i m e n s i o n a l S h e l l', & ' E l e m e n t') end subroutine shl3di(d,ul,xl,ix,s,r,ndf,ndm,nst,isw) c-----[--.----+----.----+----.-----------------------------------------] c Purpose: Quadrilateral shell element for feap c incorporating membrane with normal drilling dof c modified to remove effects of constant strains and c including deep shell curvature corrections to the c discrete kirchhoff quadrilateral plate bending element c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c.... Input parameters set as follows: c ndm = 3 (x,y,z cartesian coordinates at nodes) c ndf = 6 (u-x,u-y,u-z,r-x,r-y,r-z at nodes) c nen = 4 nodes (counterclockwise around element) c Note: 1-direction bisects diagonals between 2-3 element and c 2-direction bisects diagonals between 3-4 element nodes. c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none include 'bdata.h' include 'cdata.h' include 'eldata.h' include 'elplot.h' include 'eltran.h' include 'evdata.h' include 'iofile.h' include 'prld1.h' include 'prstrs.h' include 'shlc16.h' include 'shld16.h' include 'shpf16.h' include 'sstr16.h' include 'comblk.h' integer ndm,ndf,nst,isw integer i,j,k,l,lint,ll,lt,ialph integer i1, j1, ii,jj, nn real*8 a11,a13, a23,a33, xx,yy,zz, hh real*8 dv, dv1,dv2, pen,ggi,ggv,tc, xsj, thk real*8 xyn, x1n,x2n,y1n,y2n, ctan1,ctan3 real*8 shp1i,shp2i,shp3i real*8 shp13, shp23 integer ix(*) real*8 d(*),xl(ndm,*),ul(ndf,nen,*),s(nst,*),r(ndf,*) real*8 dvl(9),eps(6),sig(6) real*8 yl(3,4),vl(6,4,3),tr(3,3),bl(3),bg(3) real*8 gshp1(3,4),gshp2(3,4),gg(6,4),dd(6,6,5) real*8 bmat(4,6,4), dbmat(6,4) save data ialph /2/ c Transfer to correct processor go to (1,2,3,3,3,3,1,3), isw 1 return c Check element for errors 2 call tran3d(xl,yl,tr,ndm) call ckisop(ix,yl,shp,3) return c Compute element tangent array c Compute transformation and midsurface coords 3 call tran3d(xl,yl,tr,ndm) c Compute local coordinates do i = 1,4 ! { do j = 1,3 ! { vl(j ,i,1) = 0.0d0 vl(j+3,i,1) = 0.0d0 do k = 1,3 vl(j ,i,1) = vl(j ,i,1) + tr(j,k)*ul(k ,i,1) vl(j+3,i,1) = vl(j+3,i,1) + tr(j,k)*ul(k+3,i,1) end do ! k } end do ! j } end do ! i } c Get quadratrure data l = d(5) lt = nint(d(102)) call int2d (l ,lint,sg) call int1dl(lt,sgt) c Test for triangular element if( ix(1) .eq. ix(2) .or. ix(2) .eq. ix(3) .or. & ix(3) .eq. ix(4) .or. ix(4) .eq. ix(1) ) then qdflg = .false. call pzero( bm, 162) call pzero(shp1,12*lint) call pzero(shp2,12*lint) shp13 = 0.0d0 shp23 = 0.0d0 a13 = 0.0d0 a23 = 0.0d0 a33 = 0.0d0 x1n = 0.0d0 y1n = 0.0d0 x2n = 0.0d0 y2n = 0.0d0 xyn = 0.0d0 call jtri3d(yl,xsjt) else qdflg = .true. pen = d(60) call jacq3d(yl) endif if(yl(3,1).eq.0.0d0) then ii1 = 3 ii2 = 5 else ii1 = 1 ii2 = 6 endif c Construct integrals of drilling shape functions call pzero(gshp1,12) call pzero(gshp2,12) call pzero(gg ,24) dv = 0.0d0 do l = 1,lint ! { c Form shape functions and their integrals if(qdflg) then call rshp3d(sg(1,l),yl,shp(1,1,l),shp1(1,1,l), & shp2(1,1,l),xsj,3) dvl(l) = xsj*sg(3,l) do j = 1,4 ! { do i = 1,3 ! { gshp1(i,j) = gshp1(i,j) + shp1(i,j,l)*dvl(l) gshp2(i,j) = gshp2(i,j) + shp2(i,j,l)*dvl(l) end do ! i } end do ! j } else call shp2d (sg(1,l),yl,shp(1,1,l),xsj,3,nel,ix,.false.) dvl(l) = xsj*sg(3,l) endif dv = dv + dvl(l) end do ! l } if(isw.eq.5) go to 5 c Compute thickness for element thk = d(14) if(isw.eq.4) go to 4 c Construct modified drilling shape functions if(qdflg) then do j = 1,4 ! { do i = 1,3 ! { dv1 = gshp1(i,j)/dv dv2 = gshp2(i,j)/dv do l = 1,lint ! { shp1(i,j,l) = shp1(i,j,l) - dv1 shp2(i,j,l) = shp2(i,j,l) - dv2 end do ! l } end do ! j } end do ! i } endif if(isw.eq.8) go to 8 c Compute element load vectors bl(1) = 0.0d0 bl(2) = 0.0d0 c Set body loading factors if(int(d(74)).gt.0) then bg(1) = d(11) + prldv(int(d(74)))*d(71) else bg(1) = d(11)*dm endif if(int(d(75)).gt.0) then bg(2) = d(12) + prldv(int(d(75)))*d(72) else bg(2) = d(12)*dm endif if(int(d(76)).gt.0) then bl(3) = d( 8) bg(3) = d(13) + prldv(int(d(76)))*d(73) else bl(3) = d( 8)*dm bg(3) = d(13)*dm endif c Multiply body loads by thickness do i = 1,3 ! { bg(i) = bg(i)*thk end do ! i } c Transform body loads to local system do i = 1,3 ! { do j = 1,3 ! { bl(i) = bl(i) + tr(i,j)*bg(j) end do ! j } end do ! i } tc = 1.0d0 c Transform displacements do i = 1,4 vl(1,i,1) = vl(1,i,1) - yl(3,i)*vl(5,i,1) vl(2,i,1) = vl(2,i,1) + yl(3,i)*vl(4,i,1) end do c Transient factors ctan1 = ctan(1) + d(78)*ctan(2) ctan3 = ctan(3) + d(77)*ctan(2) c Compute local velocity and acceleration do i = 1,4 ! { do j = 1,3 ! { vl(j ,i,2) = 0.0d0 vl(j+3,i,2) = 0.0d0 vl(j ,i,3) = 0.0d0 vl(j+3,i,3) = 0.0d0 do k = 1,3 vl(j ,i,2) = vl(j ,i,2) + tr(j,k)*ul(k ,i,4) vl(j+3,i,2) = vl(j+3,i,2) + tr(j,k)*ul(k+3,i,4) vl(j ,i,3) = vl(j ,i,3) + tr(j,k)*ul(k ,i,5) vl(j+3,i,3) = vl(j+3,i,3) + tr(j,k)*ul(k+3,i,5) end do ! k } end do ! j } end do ! i } c Compute parts nn = 0 do l = 1,lint ! { c Compute loading term do i = 1,4 shp3i = shp(3,i,l) if(qdflg) then shp13 = shp1(3,i,l) shp23 = shp2(3,i,l) end if r(1,i) = r(1,i) + shp3i*bl(1)*dvl(l) r(2,i) = r(2,i) + shp3i*bl(2)*dvl(l) r(3,i) = r(3,i) + shp3i*bl(3)*dvl(l) r(4,i) = r(4,i) + shp13*bl(3)*dvl(l)*tc r(5,i) = r(5,i) + shp23*bl(3)*dvl(l)*tc r(6,i) = r(6,i) -(shp13*bl(1) + shp23*bl(2))*dvl(l)*tc end do ! i c Membrane and bending stiffness part dv1 = 0.5d0*thk*dvl(l)*ctan1 do ll = 1,lt hh = 0.5d0*thk*sgt(1,ll) call str3di(d,yl,vl,3,nel,l, xx,yy,zz,hh,eps,sig,dd,nn,isw) nn = nn + nint(d(15)) c Store time history plot data for element i = 8*(l-1) do j = 1,4 tt(j+i ) = sig(j) tt(j+i+3) = eps(j) end do ! j c Multiply stress and tangent by volume and quadratrue weight dv2 = dv1*sgt(2,ll) do i = 1,4 sig(i) = sig(i)*dv2 do j = 1,4 dd(j,i,1) = dd(j,i,1)*dv2 end do ! j end do ! i c Recover previously computed shape functions do i = 1,4 ! { shp1i = shp(1,i,l) shp2i = shp(2,i,l) bmat(1,1,i) = bm(1,1,i)*hh + shp1i bmat(2,1,i) = bm(2,1,i)*hh bmat(4,1,i) = bm(3,1,i)*hh + shp2i bmat(1,2,i) = bm(1,2,i)*hh bmat(2,2,i) = bm(2,2,i)*hh + shp2i bmat(4,2,i) = bm(3,2,i)*hh + shp1i bmat(1,3,i) = bm(1,3,i)*hh bmat(2,3,i) = bm(2,3,i)*hh bmat(4,3,i) = bm(3,3,i)*hh bmat(1,4,i) = bm(1,4,i)*hh bmat(2,4,i) = bm(2,4,i)*hh bmat(4,4,i) = bm(3,4,i)*hh bmat(1,5,i) = bm(1,5,i)*hh bmat(2,5,i) = bm(2,5,i)*hh bmat(4,5,i) = bm(3,5,i)*hh if(qdflg) then bmat(1,6,i) = bm(1,6,i) - shp1(1,i,l) bmat(2,6,i) = bm(2,6,i) - shp2(2,i,l) bmat(4,6,i) = bm(3,6,i) - shp1(2,i,l) - shp2(1,i,l) endif c Compute stress divergence terms do ii = 1,6 r(ii,i) = r(ii,i) - bmat(1,ii,i)*sig(1) & - bmat(2,ii,i)*sig(2) & - bmat(4,ii,i)*sig(4) end do ! ii end do ! i c Form stiffness i1 = 0 do i = 1,4 ! { c Form stress-displacement matrix (Bi-trans * D) do ii = 1,6 dbmat(ii,1) = bmat(1,ii,i)*dd(1,1,1) & + bmat(2,ii,i)*dd(2,1,1) & + bmat(4,ii,i)*dd(4,1,1) dbmat(ii,2) = bmat(1,ii,i)*dd(1,2,1) & + bmat(2,ii,i)*dd(2,2,1) & + bmat(4,ii,i)*dd(4,2,1) dbmat(ii,4) = bmat(1,ii,i)*dd(1,4,1) & + bmat(2,ii,i)*dd(2,4,1) & + bmat(4,ii,i)*dd(4,4,1) end do ! ii c Loop on columns j1 = i1 do j = i,4 ! { c Compute stiffness contribution do jj = 1,6 ! { do ii = 1,6 ! { s(ii+i1,jj+j1) = s(ii+i1,jj+j1) + + dbmat(ii,1)*bmat(1,jj,j) + + dbmat(ii,2)*bmat(2,jj,j) + + dbmat(ii,4)*bmat(4,jj,j) end do ! ii } end do ! jj } j1 = j1 + ndf end do ! j } i1 = i1 + ndf end do ! i } end do ! ll c Compute Hughes/Brezzi rotation matrix if(pen.gt.0.0d0 .and. qdflg) then do j = 1,4 ! { gg(1,j) = gg(1,j) - shp(2,j,l)*dvl(l) gg(2,j) = gg(2,j) + shp(1,j,l)*dvl(l) gg(6,j) = gg(6,j) - 2.d0*shp(3,j,l)*dvl(l) if(ialph.gt.1) then gg(6,j) = gg(6,j) + (shp1(2,j,l) - shp2(1,j,l))*dvl(l) endif end do ! j } endif end do ! l } c Perform rank one update for Huges/Brezzi term if(pen.gt.0.0d0 .and. qdflg) then c Compute H/B strain xx = 0.0d0 do i = 1,4 ! { do j = 1,6 ! { xx = xx + gg(j,i)*vl(j,i,1) end do ! j } end do ! i } c Compute H/B residual and tangent ggv = pen*d(27)*thk*ctan(1)/dv do i = 1,4 ggi = ggv*xx do j = 1,6 r(j,i) = r(j,i) - ggi*gg(j,i) end do ! j end do ! i do i = 1,nst ! { ggi = ggv*gg(i,1) do j = i,nst ! { s(i,j) = s(i,j) + ggi*gg(j,1) end do ! j } end do ! i } endif c Correct residual and tangent matrix for element warpage i1 = 1 if(yl(3,1).ne.0.0d0) then do i = 1,4 ! { r(4,i) = r(4,i) + yl(3,i)*r(2,i) r(5,i) = r(5,i) - yl(3,i)*r(1,i) j1 = i1 do j = i,4 ! { call proj3d(s(i1,j1),yl(3,i),yl(3,j),nst) j1 = j1 + ndf end do ! j } i1 = i1 + ndf end do ! i } endif c Rotate to global frame call rots3d(s,r,tr,nst,ndf) return c Compute and output element variables 4 nn = 0 do l = 1,lint ! { c Form shape functions call rshp3d(sg(1,l),yl,shp,shp1,shp2,xsj,3) c Modify rotational shape functions do j = 1,4 do i = 1,3 shp1(i,j,1) = shp1(i,j,1) - gshp1(i,j)/dv shp2(i,j,1) = shp2(i,j,1) - gshp2(i,j)/dv end do ! i } end do ! j } do ll = 1,lt hh = 0.5d0*thk*sgt(1,ll) call str3di(d,xl,vl,ndm,nel,1, xx,yy,zz,hh,eps,sig,dd,nn,isw) nn = nn + nint(d(15)) mct = mct - 3 if(mct.le.0) then write(iow,2002) o,head if(ior.lt.0) then write(*,2002) o,head endif mct = 50 end if xx = xx + tr(1,3)*hh yy = yy + tr(2,3)*hh zz = zz + tr(3,3)*hh write(iow,2003) n,xx,sig,ma,yy,eps,zz if(ior.lt.0) then write(*,2003) n,xx,sig,ma,yy,eps,zz endif end do ! ll end do ! l } return c Compute element mass or geometric stifness arrays 5 if(imtyp.eq.1) then c Compute mass do l = 1,lint ! { dv1 = dvl(l)*d(4)*d(14) i1 = 0 do j = 1,4 ! { do i = 1,3 ! { r(i,j) = r(i,j) + shp(3,j,l)*dv1 s(i1+i,i1+i) = r(i,j) end do ! i } i1 = i1 + ndf end do ! j } end do ! l } elseif(imtyp.eq.2) then c Compute geometric stiffness do i = 1,4 ! { do j = 1,3 ! { vl(j ,i,1) = 0.0d0 vl(j+3,i,1) = 0.0d0 do k = 1,3 ! { vl(j ,i,1) = vl(j ,i,1) + tr(j,k)*ul(k ,i,1) vl(j+3,i,1) = vl(j+3,i,1) + tr(j,k)*ul(k+3,i,1) end do ! k } end do ! j } end do ! i } nn = 0 do l = 1,lint ! { c Modify rotational shape functions do j = 1,4 ! { do i = 1,3 ! { shp1(i,j,l) = shp1(i,j,l) - gshp1(i,j)/dv shp2(i,j,l) = shp2(i,j,l) - gshp2(i,j)/dv end do ! i } end do ! j } hh = 0.0d0 call str3di(d,xl,vl,ndm,nel,l, xx,yy,zz,hh,eps,sig,dd,nn,isw) nn = nn + nint(d(15)) do i = 1,4 sig(i) = sig(i)*thk end do ! i i1 = 0 do i = 1,4 !{ j1 = 0 dv1 = (shp(1,i,l)*sig(1) + shp(2,i,l)*sig(4))*dvl(l) dv2 = (shp(1,i,l)*sig(4) + shp(2,i,l)*sig(2))*dvl(l) do j = 1,4 !{ a11 = dv1*shp(1,j,l) + dv2*shp(2,j,l) do k = 1,3 !{ s(i1+k,j1+k) = s(i1+k,j1+k) - a11 end do ! k } j1 = j1 + ndf end do ! j } i1 = i1 + ndf end do ! i } end do ! l } i1 = 1 if(yl(3,1).ne.0.0d0) then do i = 1,4 !{ j1 = i1 do j = i,4 !{ call proj3d(s(i1,j1),yl(3,i),yl(3,j),nst) j1 = j1 + ndf end do ! j } i1 = i1 + ndf end do ! i } endif call rots3d(s,r,tr,nst,ndf) endif return c Compute nodal output quantities 8 call stcn3si(ix,d,yl,ul,tr,hr(nph),hr(nph+numnp),dvl,thk, & ndf,nel,numnp,isw) return c Formats 2002 format(a1,20a4//5x,'S h e l l S t r e s s e s'// & ' elmt x-coord xx-stress yy-stress zz-stress xy-stress', & ' matl y-coord xx-strain yy-strain zz-strain xy-strain', & ' z-coord'/38(' -')) 2003 format(/i5,0p,1f8.3,1p,5e11.3,0p,1f8.2/i5,0p,1f8.3,1p,5e11.3, & 0p,1f8.2/5x,0p,1f8.3,1p,6e11.3/(13x,1p,6e11.3)) end subroutine stcn3si(ix,d,yl,ul,tr,dt,st,dvl,thk, & ndf,nel,numnp,isw) implicit none include 'shpf16.h' include 'shlc16.h' include 'sstr16.h' integer ndf,nel,numnp, i,j,l, ii, ll,lt, nn,isw real*8 xsji, xx,yy,zz, hh, dh, thk integer ix(*) real*8 dt(numnp),st(numnp,*),yl(3,*),tr(3,3),d(*),ul(ndf,*) real*8 vl(6,4),dvl(4) real*8 eps(6),sig(6),momt(6),norm(6),dd(6,6,5),eps0(6) save c Compute membrane and bending stresses for projection. do i = 1,4 ! { do j = 1,3 ! { vl(j ,i) = 0.0d0 vl(j+3,i) = 0.0d0 do l = 1,3 ! { vl(j ,i) = vl(j ,i) + tr(j,l)*ul(l ,i) vl(j+3,i) = vl(j+3,i) + tr(j,l)*ul(l+3,i) end do ! l } end do ! j } end do ! i } lt = nint(d(102)) nn = 0 do l = 1,4 ! { do j = 1,6 norm(j) = 0.0d0 momt(j) = 0.0d0 end do ! j do ll = 1,lt hh = 0.5d0*thk*sgt(1,ll) dh = 0.5d0*thk*sgt(2,ll) call str3di(d,yl,vl,3,nel,l, xx,yy,zz,hh,eps,sig,dd,nn,isw) norm(1) = norm(1) + sig(1)*dh norm(2) = norm(2) + sig(2)*dh norm(3) = norm(3) + sig(4)*dh momt(1) = momt(1) + sig(1)*hh*dh momt(2) = momt(2) + sig(2)*hh*dh momt(3) = momt(3) + sig(4)*hh*dh if(ll.eq.1) then do j = 1,4 ! { xsji = dvl(l)*shp(3,j,l) ii = ix(j) if(ii.gt.0) then st(ii,20) = st(ii,20) + sig(1)*xsji st(ii,21) = st(ii,21) + sig(2)*xsji st(ii,22) = st(ii,22) + sig(4)*xsji endif end do ! j do j = 1,3 eps0(j) = eps(j) end do elseif(ll.eq.lt) then do j = 1,4 ! { xsji = dvl(l)*shp(3,j,l) ii = ix(j) if(ii.gt.0) then st(ii,17) = st(ii,17) + sig(1)*xsji st(ii,18) = st(ii,18) + sig(2)*xsji st(ii,19) = st(ii,19) + sig(4)*xsji endif end do ! j do j = 1,3 eps0(j+3) = (eps(j) - eps0(j))/thk eps0(j ) = (eps(j) + eps0(j))*0.5d0 end do endif nn = nn + nint(d(15)) end do ! ll do j = 1,4 ! { xsji = dvl(l)*shp(3,j,l) ii = ix(j) if(ii.gt.0) then dt(ii) = dt(ii) + xsji st(ii,1) = st(ii,1) + norm(1)*xsji st(ii,2) = st(ii,2) + norm(2)*xsji st(ii,3) = st(ii,3) + norm(3)*xsji st(ii,4) = st(ii,4) + norm(4)*xsji st(ii,5) = st(ii,5) + norm(5)*xsji st(ii,6) = st(ii,6) + momt(1)*xsji st(ii,7) = st(ii,7) + momt(2)*xsji st(ii,8) = st(ii,8) + momt(3)*xsji st(ii,9) = st(ii,9) + momt(4)*xsji st(ii,10) = st(ii,10) + momt(5)*xsji st(ii,11) = st(ii,11) + eps0(1)*xsji st(ii,12) = st(ii,12) + eps0(2)*xsji st(ii,13) = st(ii,13) + eps0(3)*xsji st(ii,14) = st(ii,14) + eps0(4)*xsji st(ii,15) = st(ii,15) + eps0(5)*xsji st(ii,16) = st(ii,16) + eps0(6)*xsji end if end do ! j } end do ! l } end subroutine str3di(d,xl,vl,ndm,nel,l, xx,yy,zz,hh, eps,sig,dd, & nn,isw) implicit none include 'hdata.h' include 'shlc16.h' include 'shld16.h' include 'shpf16.h' include 'sstr16.h' include 'comblk.h' integer ndm,nel,l, i,j, nn,isw, nhv, istrt real*8 xx,yy,zz, hh, ta real*8 d(*), xl(ndm,*),vl(6,*), eps0(6),eps(6),sig(6) real*8 dd(6,6,5) save c Compute membrane and bending strains if(qdflg) then call dktq3d(shp(1,5,l),shp(1,1,l)) else call dktb3d(sg(1,l),xsjt) endif do i = 1,6 ! { eps0(i) = 0.0 eps (i) = 0.0 end do ! i } xx = 0.0d0 yy = 0.0d0 zz = 0.0d0 do j = 1,nel ! { xx = xx + shp(3,j,l)*xl(1,j) yy = yy + shp(3,j,l)*xl(2,j) zz = zz + shp(3,j,l)*xl(3,j) eps0(1) = eps0(1) + shp(1,j,l)*vl(1,j) & - shp1(1,j,l)*vl(6,j) eps0(2) = eps0(2) + shp(2,j,l)*vl(2,j) & - shp2(2,j,l)*vl(6,j) eps0(3) = eps0(3) + shp(1,j,l)*vl(2,j) + shp(2,j,l)*vl(1,j) & - (shp1(2,j,l) + shp2(1,j,l))*vl(6,j) do i = ii1,ii2 ! { eps0(4) = eps0(4) + bm(1,i,j)*vl(i,j) eps0(5) = eps0(5) + bm(2,i,j)*vl(i,j) eps0(6) = eps0(6) + bm(3,i,j)*vl(i,j) end do ! i } end do ! j } c Compute strain at current elevation eps(1) = eps0(1) + hh*eps0(4) eps(2) = eps0(2) + hh*eps0(5) eps(4) = eps0(3) + hh*eps0(6) ta = 0.0d0 c Compute stress at point nhv = nint(d(15)) istrt = nint(d(84)) call modlsd(d,ta,eps,hr(nh1+nn),hr(nh2+nn),nhv,istrt,dd,sig,isw) end subroutine shl3ds(d,ul,xl,ix,s,r,ndf,ndm,nst,isw) c-----[--.----+----.----+----.-----------------------------------------] c Purpose: Quadrilateral shell element for feap c incorporating membrane with normal drilling dof c modified to remove effects of constant strains and c including deep shell curvature corrections to the c discrete kirchhoff quadrilateral plate bending element c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - c.... Input parameters set as follows: c ndm = 3 (x,y,z cartesian coordinates at nodes) c ndf = 6 (u-x,u-y,u-z,r-x,r-y,r-z at nodes) c nen = 4 nodes (counterclockwise around element) c Note: 1-direction bisects diagonals between 2-3 element and c 2-direction bisects diagonals between 3-4 element nodes. c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - implicit none include 'bdata.h' include 'cdata.h' include 'eldata.h' include 'elplot.h' include 'eltran.h' include 'evdata.h' include 'iofile.h' include 'prld1.h' include 'prstrs.h' include 'shlc16.h' include 'shld16.h' include 'shpf16.h' include 'sstr16.h' include 'comblk.h' integer ndm,ndf,nst,isw integer i,j,k,l,lint,ialph, i1, j1, ii,jj,kk real*8 a11,a12,a13, a21,a22,a23, a31,a32,a33, xx,yy,zz real*8 dv, dv1,dv2,dvm, pen,ggi,ggv,tc, xsj, thk,thk3 real*8 xyn, x1n,x2n,y1n,y2n, xn,yn, ctan1,ctan3 real*8 shp1i,shp2i,shp3i real*8 shp11,shp12,shp13, shp21,shp22,shp23 integer ix(*) real*8 d(*),xl(ndm,*),ul(ndf,nen,*),s(nst,*),r(ndf,*) real*8 norm(6),dvl(9),momt(6),eps(6) real*8 yl(3,4),vl(6,4,3),tr(3,3),bl(3),bg(3) real*8 btd(3,6),gshp1(3,4),gshp2(3,4),gg(24),dd(6,6) c real*8 b0(4,2,4),bf(4,3,4), ddb0(4,2), ddbf(4,3) save data ialph /2/ c Transfer to correct processor go to (1,2,3,3,3,3,1,3), isw 1 return c Check element for errors 2 call tran3d(xl,yl,tr,ndm) call ckisop(ix,yl,shp,3) return c Compute element tangent array c Compute transformation and midsurface coords 3 call tran3d(xl,yl,tr,ndm) c Compute local coordinates do i = 1,4 ! { do j = 1,3 ! { vl(j ,i,1) = 0.0d0 vl(j+3,i,1) = 0.0d0 do k = 1,3 vl(j ,i,1) = vl(j ,i,1) + tr(j,k)*ul(k ,i,1) vl(j+3,i,1) = vl(j+3,i,1) + tr(j,k)*ul(k+3,i,1) end do ! k } end do ! j } end do ! i } c Get quadratrure data l = d(5) if(l*l.ne.lint) call int2d(l ,lint,sg) c Test for triangular element if( ix(1) .eq. ix(2) .or. ix(2) .eq. ix(3) .or. & ix(3) .eq. ix(4) .or. ix(4) .eq. ix(1) ) then qdflg = .false. call pzero( bm, 162) shp11 = 0.0d0 shp12 = 0.0d0 shp13 = 0.0d0 shp21 = 0.0d0 shp22 = 0.0d0 shp23 = 0.0d0 a13 = 0.0d0 a23 = 0.0d0 a33 = 0.0d0 x1n = 0.0d0 y1n = 0.0d0 x2n = 0.0d0 y2n = 0.0d0 xyn = 0.0d0 call jtri3d(yl,xsjt) else qdflg = .true. pen = d(60) call jacq3d(yl) endif c if(qdflg .and. yl(3,1).eq.0.0d0) then if(yl(3,1).eq.0.0d0) then ii1 = 3 ii2 = 5 else ii1 = 1 ii2 = 6 endif c Construct integrals of drilling shape functions call pzero(gshp1,12) call pzero(gshp2,12) call pzero(gg ,24) dv = 0.0d0 do l = 1,lint ! { c Form shape functions and their integrals if(qdflg) then call rshp3d(sg(1,l),yl,shp(1,1,l),shp1(1,1,l), & shp2(1,1,l),xsj,3) dvl(l) = xsj*sg(3,l) do j = 1,4 ! { do i = 1,3 ! { gshp1(i,j) = gshp1(i,j) + shp1(i,j,l)*dvl(l) gshp2(i,j) = gshp2(i,j) + shp2(i,j,l)*dvl(l) end do ! i } end do ! j } else call shp2d (sg(1,l),yl,shp(1,1,l),xsj,3,nel,ix,.false.) dvl(l) = xsj*sg(3,l) endif dv = dv + dvl(l) end do ! l } if(isw.eq.5) go to 5 c Compute thickness for element thk = d(14) thk3 = thk**3/12.d0 if(isw.eq.4) go to 4 c Construct modified drilling shape functions if(qdflg) then do j = 1,4 ! { do i = 1,3 ! { dv1 = gshp1(i,j)/dv dv2 = gshp2(i,j)/dv do l = 1,lint ! { shp1(i,j,l) = shp1(i,j,l) - dv1 shp2(i,j,l) = shp2(i,j,l) - dv2 end do ! l } end do ! i } end do ! j } endif if(isw.eq.8) go to 8 c Compute element load vectors bl(1) = 0.0d0 bl(2) = 0.0d0 c Set body loading factors if(int(d(74)).gt.0) then bg(1) = d(11) + prldv(int(d(74)))*d(71) else bg(1) = d(11)*dm endif if(int(d(75)).gt.0) then bg(2) = d(12) + prldv(int(d(75)))*d(72) else bg(2) = d(12)*dm endif if(int(d(76)).gt.0) then bl(3) = d( 8) bg(3) = d(13) + prldv(int(d(76)))*d(73) else bl(3) = d( 8)*dm bg(3) = d(13)*dm endif do i = 1,3 ! { do j = 1,3 ! { bl(i) = bl(i) + tr(i,j)*bg(j) end do ! j } end do ! i } tc = 1.0d0 c Transform displacements do i = 1,4 vl(1,i,1) = vl(1,i,1) - yl(3,i)*vl(5,i,1) vl(2,i,1) = vl(2,i,1) + yl(3,i)*vl(4,i,1) end do c Transient factors ctan1 = ctan(1) + d(78)*ctan(2) ctan3 = ctan(3) + d(77)*ctan(2) c Compute local velocity and acceleration do i = 1,4 ! { do j = 1,3 ! { vl(j ,i,2) = 0.0d0 vl(j+3,i,2) = 0.0d0 vl(j ,i,3) = 0.0d0 vl(j+3,i,3) = 0.0d0 do k = 1,3 vl(j ,i,2) = vl(j ,i,2) + tr(j,k)*ul(k ,i,4) vl(j+3,i,2) = vl(j+3,i,2) + tr(j,k)*ul(k+3,i,4) vl(j ,i,3) = vl(j ,i,3) + tr(j,k)*ul(k ,i,5) vl(j+3,i,3) = vl(j+3,i,3) + tr(j,k)*ul(k+3,i,5) end do ! k } end do ! j } end do ! i } c Compute membrane/load parts do l = 1,lint ! { c Membrane and bending stiffness part dv1 = thk *dvl(l)*ctan1 dv2 = thk3*dvl(l)*ctan1 call stre3d(d,yl,vl,3,nel,l, xx,yy,zz,eps,norm,momt,dd) c Store time history plot data for element i = 6*(l-1) do j = 1,3 tt(j+i ) = norm(j) tt(j+i+3) = momt(j) end do ! j c Add Rayleigh damping stress terms if(d(78).ne.0.0d0) then call rayl3d(d,vl(1,1,2),nel,l, norm,momt,dd) endif c Recover previously computed shape functions i1 = 0 do i = 1,4 ! { shp1i = shp(1,i,l) shp2i = shp(2,i,l) shp3i = shp(3,i,l) if(qdflg) then shp11 = shp1(1,i,l) shp12 = shp1(2,i,l) shp13 = shp1(3,i,l) shp21 = shp2(1,i,l) shp22 = shp2(2,i,l) shp23 = shp2(3,i,l) a13 = -(dd(1,1)*shp11 + dd(1,2)*shp22 & + dd(1,4)*(shp12+shp21))*dv1 a23 = -(dd(2,1)*shp11 + dd(2,2)*shp22 & + dd(2,4)*(shp12+shp21))*dv1 a33 = -(dd(4,1)*shp11 + dd(4,2)*shp22 & + dd(4,4)*(shp12+shp21))*dv1 endif c Compute loading term r(1,i) = r(1,i) + shp3i*bl(1)*dvl(l) r(2,i) = r(2,i) + shp3i*bl(2)*dvl(l) r(3,i) = r(3,i) + shp3i*bl(3)*dvl(l) r(4,i) = r(4,i) + shp13*bl(3)*dvl(l)*tc r(5,i) = r(5,i) + shp23*bl(3)*dvl(l)*tc r(6,i) = r(6,i) -(shp13*bl(1) + shp23*bl(2))*dvl(l)*tc c Compute stress divergence terms r(1,i) = r(1,i) - (shp1i*norm(1) + shp2i*norm(3))*dvl(l) r(2,i) = r(2,i) - (shp2i*norm(2) + shp1i*norm(3))*dvl(l) r(6,i) = r(6,i) + (shp11*norm(1) + shp22*norm(2) & + (shp12 + shp21)*norm(3))*dvl(l) do ii = ii1,ii2 ! { r(ii,i) = r(ii,i) - (bm(1,ii,i)*momt(1) & + bm(2,ii,i)*momt(2) & + bm(3,ii,i)*momt(3))*dvl(l) end do ! ii } c Form stress-displacement matrix (Bi-trans * D) a11 = (dd(1,1)*shp1i + dd(1,4)*shp2i)*dv1 a12 = (dd(1,2)*shp2i + dd(1,4)*shp1i)*dv1 a21 = (dd(2,1)*shp1i + dd(2,4)*shp2i)*dv1 a22 = (dd(2,2)*shp2i + dd(2,4)*shp1i)*dv1 a31 = (dd(4,1)*shp1i + dd(4,4)*shp2i)*dv1 a32 = (dd(4,2)*shp2i + dd(4,4)*shp1i)*dv1 c Form plate stress-displacement matrix do ii = ii1,ii2 ! { btd(1,ii) = (dd(1,1)*bm(1,ii,i) & + dd(1,2)*bm(2,ii,i) & + dd(1,4)*bm(3,ii,i))*dv2 btd(2,ii) = (dd(2,1)*bm(1,ii,i) & + dd(2,2)*bm(2,ii,i) & + dd(2,4)*bm(3,ii,i))*dv2 btd(3,ii) = (dd(4,1)*bm(1,ii,i) & + dd(4,2)*bm(2,ii,i) & + dd(4,4)*bm(3,ii,i))*dv2 end do ! ii } c Inertial parts dvm = dvl(l)*d(4)*d(14)*shp3i do ii = 1,3 ! { r(ii,i) = r(ii,i) - dvm*( vl(ii,i,3) & + d(77)*vl(ii,i,2)) s(i1+ii,i1+ii) = s(i1+ii,i1+ii) + dvm*ctan3 end do ! i } c Loop on columns j1 = i1 do j = i,4 ! { xn = shp(1,j,l) yn = shp(2,j,l) if(qdflg) then x1n = - shp1(1,j,l) y2n = - shp2(2,j,l) x2n = - shp1(2,j,l) y1n = - shp2(1,j,l) xyn = x2n + y1n endif c Compute membrane part s(i1+1,j1+1) = s(i1+1,j1+1) + (a11*xn + a31*yn) s(i1+2,j1+1) = s(i1+2,j1+1) + (a12*xn + a32*yn) s(i1+1,j1+2) = s(i1+1,j1+2) + (a21*yn + a31*xn) s(i1+2,j1+2) = s(i1+2,j1+2) + (a22*yn + a32*xn) s(i1+6,j1+1) = s(i1+6,j1+1) + a13*xn + a33*yn s(i1+6,j1+2) = s(i1+6,j1+2) + a23*yn + a33*xn s(i1+1,j1+6) = s(i1+1,j1+6) + a11*x1n + a21*y2n + a31*xyn s(i1+2,j1+6) = s(i1+2,j1+6) + a12*x1n + a22*y2n + a32*xyn s(i1+6,j1+6) = s(i1+6,j1+6) + a13*x1n + a23*y2n + a33*xyn c Compute bending part do ii = ii1,ii2 ! { do jj = ii1,ii2 ! { do kk = 1,3 ! { s(i1+ii,j1+jj) = s(i1+ii,j1+jj) & + btd(kk,ii)*bm(kk,jj,j) end do ! kk } end do ! jj } end do ! ii } j1 = j1 + ndf end do ! j } i1 = i1 + ndf end do ! i } c Compute Hughes/Brezzi rotation matrix if(ialph.ne.0 .and. qdflg) then j1 = 0 do j = 1,4 ! { gg(j1+1) = gg(j1+1) - shp(2,j,l)*dvl(l) gg(j1+2) = gg(j1+2) + shp(1,j,l)*dvl(l) gg(j1+6) = gg(j1+6) - 2.d0*shp(3,j,l)*dvl(l) if(ialph.gt.1) then gg(j1+6) = gg(j1+6) + (shp1(2,j,l) - shp2(1,j,l))*dvl(l) endif j1 = j1 + ndf end do ! j } endif end do ! l } c Perform rank one update for Huges/Brezzi term if(ialph.gt.0 .and. qdflg) then c Compute H/B strain xx = 0.0d0 j1 = 0 do i = 1,4 ! { do j = 1,6 ! { xx = xx + gg(j1+j)*vl(j,i,1) end do ! j } j1 = j1 + ndf end do ! i } c Compute H/B residual and tangent ggv = pen*d(27)*thk*ctan1/dv do i = 1,nst ! { ggi = ggv*gg(i) r(i,1) = r(i,1) - ggi*xx do j = 1,nst ! { s(i,j) = s(i,j) + ggi*gg(j) end do ! j } end do ! i } endif c Correct residual and tangent matrix for element warpage i1 = 1 if(yl(3,1).ne.0.0d0) then do i = 1,4 ! { r(4,i) = r(4,i) + yl(3,i)*r(2,i) r(5,i) = r(5,i) - yl(3,i)*r(1,i) j1 = i1 do j = i,4 ! { call proj3d(s(i1,j1),yl(3,i),yl(3,j),nst) j1 = j1 + ndf end do ! j } i1 = i1 + ndf end do ! i } endif c Rotate to global frame call rots3d(s,r,tr,nst,ndf) return c Compute and output element variables 4 l = d(6) if(l*l.ne.lint) call int2d(l,lint,sg) do l = 1,lint ! { c Form shape functions if(qdflg) then call rshp3d(sg(1,l),yl,shp,shp1,shp2,xsj,3) c Modify rotational shape functions do j = 1,4 do i = 1,3 shp1(i,j,1) = shp1(i,j,1) - gshp1(i,j)/dv shp2(i,j,1) = shp2(i,j,1) - gshp2(i,j)/dv end do ! i } end do ! j } else call shp2d (sg(1,l),yl,shp,xsj,3,nel,ix,.false.) endif call stre3d(d,xl,vl,ndm,nel,1, xx,yy,zz,eps,norm,momt,dd) mct = mct - 3 if(mct.le.0) then write(iow,2002) o,head if(ior.lt.0) then write(*,2002) o,head endif mct = 50 end if write(iow,2003) n,xx,norm,ma,yy,momt,zz,eps,sigt,sigb if(ior.lt.0) then write(*,2003) n,xx,norm,ma,yy,momt,zz,eps,sigt,sigb endif end do ! l } return c Compute element mass or geometric stifness arrays 5 if(imtyp.eq.1) then c Compute mass do l = 1,lint ! { dv1 = dvl(l)*d(4)*d(14) i1 = 0 do j = 1,4 ! { do i = 1,3 ! { r(i,j) = r(i,j) + shp(3,j,l)*dv1 s(i1+i,i1+i) = r(i,j) end do ! i } i1 = i1 + ndf end do ! j } end do ! l } elseif(imtyp.eq.2) then c Compute geometric stiffness do i = 1,4 ! { do j = 1,3 ! { vl(j ,i,1) = 0.0d0 vl(j+3,i,1) = 0.0d0 do k = 1,3 ! { vl(j ,i,1) = vl(j ,i,1) + tr(j,k)*ul(k ,i,1) vl(j+3,i,1) = vl(j+3,i,1) + tr(j,k)*ul(k+3,i,1) end do ! k } end do ! j } end do ! i } do l = 1,lint ! { c Modify rotational shape functions do j = 1,4 ! { do i = 1,3 ! { shp1(i,j,l) = shp1(i,j,l) - gshp1(i,j)/dv shp2(i,j,l) = shp2(i,j,l) - gshp2(i,j)/dv end do ! i } end do ! j } call stre3d(d,xl,vl,ndm,nel,l, xx,yy,zz,eps,norm,momt,dd) i1 = 0 do i = 1,4 !{ j1 = 0 dv1 = (shp(1,i,l)*norm(1) + shp(2,i,l)*norm(3))*dvl(l) dv2 = (shp(1,i,l)*norm(3) + shp(2,i,l)*norm(2))*dvl(l) do j = 1,4 !{ a11 = dv1*shp(1,j,l) + dv2*shp(2,j,l) do k = 1,3 !{ s(i1+k,j1+k) = s(i1+k,j1+k) - a11 end do ! k } j1 = j1 + ndf end do ! j } i1 = i1 + ndf end do ! i } end do ! l } i1 = 1 if(yl(3,1).ne.0.0d0) then do i = 1,4 !{ j1 = i1 do j = i,4 !{ call proj3d(s(i1,j1),yl(3,i),yl(3,j),nst) j1 = j1 + ndf end do ! j } i1 = i1 + ndf end do ! i } endif call rots3d(s,r,tr,nst,ndf) endif return c Compute nodal output quantities 8 call stcn3sh(ix,d,yl,ul,tr,hr(nph),hr(nph+numnp),dvl, & ndf,nel,numnp) return c Formats 2002 format(a1,20a4//5x,'S h e l l S t r e s s e s'// & ' elmt x-coord xx-stress yy-stress xy-stress 1-stress', & ' 2-stress angle'/' matl y-coord xx-moment yy-moment', & ' xy-moment 1-moment 2-moment angle'/ & ' z-coord xx-strain yy-strain xy-strain xx-curvtr', & ' yy-curvtr xy-curvtr'/13x,' xx-sig(t) yy-sig(t) xy-sig(t)', & ' 1-sig(t) 2-sig(t) angle'/13x,' xx-sig(b) yy-sig(b)', & ' xy-sig(b) 1-sig(b) 2-sig(b) angle'/38(' -')) 2003 format(/i5,0p,1f8.3,1p,5e11.3,0p,1f8.2/i5,0p,1f8.3,1p,5e11.3, & 0p,1f8.2/5x,0p,1f8.3,1p,6e11.3/(13x,1p,6e11.3)) end subroutine dktb3d(sg,xsjt) implicit none include 'shld16.h' include 'shle16.h' integer i,j,k real*8 sg(2),el(3),shn(2,3),bd(3),cd(3),xsjt,shm1,shm2 real*8 a1(3),a2(3),b1(3),b2(3),c1(3),c2(3),d1(3),d2(3) real*8 e1(3),e2(3) save c Compute area coordinates for point el(1) = 0.25d0*(1.d0 - sg(1))*(1.d0 - sg(2)) el(2) = 0.25d0*(1.d0 + sg(1))*(1.d0 - sg(2)) el(3) = 0.50d0*(1.d0 + sg(2)) c Form shape function terms for strains do i = 1,3 ! { bd(i) = b(i)/xsjt cd(i) = c(i)/xsjt end do ! i } do i = 1,3 ! { j = mod(i,3) + 1 k = mod(j,3) + 1 shn(1,i) = (bd(i)*(el(i) - 1.0))*2.d0 shn(2,i) = (cd(i)*(el(i) - 1.0))*2.d0 shm1 = (bd(j)*el(k) + bd(k)*el(j))*2.d0 shm2 = (cd(j)*el(k) + cd(k)*el(j))*2.d0 a1(i) = aa(i)*shm1 a2(i) = aa(i)*shm2 b1(i) = bb(i)*shm1 b2(i) = bb(i)*shm2 c1(i) = cc(i)*shm1 c2(i) = cc(i)*shm2 d1(i) = dd(i)*shm1 d2(i) = dd(i)*shm2 e1(i) = ee(i)*shm1 e2(i) = ee(i)*shm2 end do ! i } c Plate strain displacement matrix do i = 1,3 ! { j = mod(i,3) + 1 k = mod(j,3) + 1 bm(1,3,i) = a1(k) - a1(j) bm(1,4,i) = b1(k) + b1(j) bm(1,5,i) = c1(k) + c1(j) - shn(1,i) bm(2,3,i) = d2(k) - d2(j) bm(2,4,i) =-e2(k) - e2(j) + shn(2,i) bm(2,5,i) =-b2(k) - b2(j) bm(3,3,i) = a2(k) - a2(j) + d1(k) - d1(j) bm(3,4,i) =-e1(k) - e1(j) + shn(1,i) - bm(2,5,i) bm(3,5,i) = c2(k) + c2(j) - shn(2,i) - bm(1,4,i) end do ! i } end subroutine dktq3d(shm,shn) implicit none include 'shld16.h' include 'shle16.h' integer i,j real*8 shm(3,4),shn(3,4) save c Form strain-displacement array for DKQ element do i = 1,4 ! { j = mod(i+2,4) + 1 bm(1,3,i) = aa(i)*shm(1,i) - aa(j)*shm(1,j) bm(1,4,i) = bb(i)*shm(1,i) + bb(j)*shm(1,j) bm(1,5,i) = cc(i)*shm(1,i) + cc(j)*shm(1,j) - shn(1,i) bm(2,3,i) = dd(i)*shm(2,i) - dd(j)*shm(2,j) bm(2,4,i) =-ee(i)*shm(2,i) - ee(j)*shm(2,j) + shn(2,i) bm(2,5,i) =-bb(i)*shm(2,i) - bb(j)*shm(2,j) bm(3,3,i) = aa(i)*shm(2,i) - aa(j)*shm(2,j) & + dd(i)*shm(1,i) - dd(j)*shm(1,j) bm(3,4,i) =-ee(i)*shm(1,i) - ee(j)*shm(1,j) + shn(1,i) & - bm(2,5,i) bm(3,5,i) = cc(i)*shm(2,i) + cc(j)*shm(2,j) - shn(2,i) & - bm(1,4,i) end do ! i } end subroutine hshp3d(shp,shp1,shp2) implicit none include 'shle16.h' integer i,j,l real*8 shp(3,4),shp1(3,4),shp2(3,4),shx(4),shy(4) save c Form shape functions do l = 1,3 ! { do i = 1,4 ! { shx(i) = shp(l,i)*c(i) shy(i) = shp(l,i)*b(i) end do ! i } j = 4 do i = 1,4 ! { shp1(l,i) = shy(i) - shy(j) shp2(l,i) = shx(i) - shx(j) j = i end do ! i } end do ! l } end subroutine jacq3d(xl) implicit none include 'shld16.h' include 'shle16.h' integer i,k real*8 cxx,cyy,cxy, dz, ad, a4, b2,c2, sql, x31,y31, x42,y42 real*8 xl(3,*) save c Form geometric constants for DKQ element cxx = 0. cyy = 0. cxy = 0. do i = 1,4 ! { k = mod(i,4) + 1 b(i) = xl(2,k) - xl(2,i) c(i) = xl(1,i) - xl(1,k) dz = xl(3,k) - xl(3,i) b2 = b(i)*b(i) c2 = c(i)*c(i) cxx = cxx + dz*b2/(b2 + c2) cyy = cyy + dz*c2/(b2 + c2) cxy = cxy + (dz+dz)*b(i)*c(i)/(b2 + c2) sql = (b2 + c2)/0.75 aa(i) = (c(i)+c(i))/sql bb(i) = b(i)*c(i)/sql cc(i) = c2/sql dd(i) =-(b(i)+b(i))/sql ee(i) = b2/sql b(i) = b(i)/8. c(i) = c(i)/8. end do ! i } if(xl(3,1).ne.0.0d0) then x31 = xl(1,3) - xl(1,1) y31 = xl(2,3) - xl(2,1) x42 = xl(1,4) - xl(1,2) y42 = xl(2,4) - xl(2,2) a4 = (x31*y42 - x42*y31)*2.0 ad = a4*a4/4.0 x31 = x31/ad y31 = y31/ad x42 = x42/ad y42 = y42/ad bm(1,1,1) = -cxx*x42 bm(2,1,1) = -cyy*x42 bm(3,1,1) = -cxy*x42 bm(1,2,1) = -cxx*y42 bm(2,2,1) = -cyy*y42 bm(3,2,1) = -cxy*y42 bm(1,6,1) = -cxx/a4 bm(2,6,1) = -cyy/a4 bm(3,6,1) = -cxy/a4 bm(1,1,2) = cxx*x31 bm(2,1,2) = cyy*x31 bm(3,1,2) = cxy*x31 bm(1,2,2) = cxx*y31 bm(2,2,2) = cyy*y31 bm(3,2,2) = cxy*y31 bm(1,6,2) = -cxx/a4 bm(2,6,2) = -cyy/a4 bm(3,6,2) = -cxy/a4 do i = 1,3 ! { bm(i,1,3) = -bm(i,1,1) bm(i,2,3) = -bm(i,2,1) bm(i,6,3) = bm(i,6,1) bm(i,1,4) = -bm(i,1,2) bm(i,2,4) = -bm(i,2,2) bm(i,6,4) = bm(i,6,2) end do ! i } endif end subroutine jtri3d(xl,xsjt) implicit none include 'shld16.h' include 'shle16.h' integer i,j,k real*8 xl(3,*), xsjt,sql,cs,bs save c Form terms for jacobian of a triangle do i = 1,3 ! { j = mod(i,3) + 1 k = mod(j,3) + 1 b(i) = xl(2,k) - xl(2,j) c(i) = xl(1,j) - xl(1,k) sql = (b(i)+b(i))*(b(i)+b(i)) + (c(i)+c(i))*(c(i)+c(i)) cs = c(i)/sql bs = b(i)/sql aa(i) = 6.d0*cs bb(i) = 3.d0*bs*c(i) cc(i) = c(i)*cs - b(i)*(bs+bs) dd(i) =-6.d0*bs ee(i) = b(i)*bs - c(i)*(cs+cs) end do ! i } b(4) = 0.0d0 c(4) = 0.0d0 aa(4) = 0.0d0 bb(4) = 0.0d0 cc(4) = 0.0d0 dd(4) = 0.0d0 ee(4) = 0.0d0 c Store jacobian xsjt = b(1)*c(2) - b(2)*c(1) end subroutine proj3d(s,zi,zj,nst) implicit none integer i, nst real*8 zi, zj real*8 s(nst,*) c Modify stiffness for offset projections c Postmultiply by transformation do i = 1,6 ! { s(i,4) = zj*s(i,2) + s(i,4) s(i,5) = -zj*s(i,1) + s(i,5) end do ! i } c Premultiply using modified terms from postmultiplication do i = 1,6 ! { s(4,i) = zi*s(2,i) + s(4,i) s(5,i) = -zi*s(1,i) + s(5,i) end do ! i } end subroutine rots3d(s,p,t,nst,ndf) c Transform loads and stiffness to global coords. implicit none integer nst,ndf, i,i0,i1, ir,ii, j,j0,j1, jc,jj real*8 s(nst,nst),p(nst),t(3,3),a(3,3),b(6) real*8 dot i0 = 0 do ir = 1,4 ! { do ii = 1,3 ! { b(ii ) = dot(t(1,ii),p(i0+1),3) b(ii+3) = dot(t(1,ii),p(i0+4),3) end do ! ii } do ii = 1,6 ! { p(i0+ii) = b(ii) end do ! ii } j0 = i0 do jc = ir,4 ! { i1 = i0 do i = 1,2 ! { j1 = j0 do j = 1,2 ! { do ii = 1,3 ! { do jj = 1,3 ! { a(jj,ii) = dot(t(1,ii),s(i1+1,jj+j1),3) end do ! jj } end do ! ii } do jj = 1,3 ! { do ii = 1,3 ! { s(ii+i1,jj+j1) = dot(a(1,ii),t(1,jj),3) end do ! ii } end do ! jj } j1 = j1 + 3 end do ! j } i1 = i1 + 3 end do ! i } c Compute symmetric block if(ir.ne.jc) then do i = 1,6 ! { do j = 1,6 ! { s(j0+j,i0+i) = s(i0+i,j0+j) end do ! j } end do ! i } endif j0 = j0 + ndf end do ! jc } i0 = i0 + ndf end do ! ir } end subroutine rshp3d(sg,x,shp,shp1,shp2,xsj,ndm) c Shape function routine for two dimensional elements implicit none integer ndm, i real*8 s2,t2, tp, xsj, xsj1 real*8 sg(2),shp(3,8),x(ndm,*),sx(2,2),shp1(3,4),shp2(3,4) c Form 4-node quadrilateral shape functions shp(1,2) = 0.25d0*(1.d0-sg(2)) shp(1,3) = 0.25d0*(1.d0+sg(2)) shp(1,1) = -shp(1,2) shp(1,4) = -shp(1,3) shp(2,4) = 0.25d0*(1.d0-sg(1)) shp(2,3) = 0.25d0*(1.d0+sg(1)) shp(2,2) = -shp(2,3) shp(2,1) = -shp(2,4) shp(3,1) = shp(1,2)*(1.d0-sg(1)) shp(3,2) = shp(1,2)*(1.d0+sg(1)) shp(3,3) = shp(1,3)*(1.d0+sg(1)) shp(3,4) = shp(1,3)*(1.d0-sg(1)) s2 = (1.d0-sg(1)*sg(1))*0.5d0 t2 = (1.d0-sg(2)*sg(2))*0.5d0 shp(3,5) = s2*(1.d0-sg(2)) shp(3,6) = t2*(1.d0+sg(1)) shp(3,7) = s2*(1.d0+sg(2)) shp(3,8) = t2*(1.d0-sg(1)) shp(1,5) = -sg(1)*(1.d0-sg(2)) shp(1,6) = t2 shp(1,7) = -sg(1)*(1.d0+sg(2)) shp(1,8) = -t2 shp(2,5) = -s2 shp(2,6) = -sg(2)*(1.d0+sg(1)) shp(2,7) = s2 shp(2,8) = -sg(2)*(1.d0-sg(1)) c Construct jacobian and its inverse sx(2,2) = x(1,1)*shp(1,1)+x(1,2)*shp(1,2)+x(1,3)*shp(1,3) & +x(1,4)*shp(1,4) sx(1,2) = x(1,1)*shp(2,1)+x(1,2)*shp(2,2)+x(1,3)*shp(2,3) & +x(1,4)*shp(2,4) sx(2,1) = x(2,1)*shp(1,1)+x(2,2)*shp(1,2)+x(2,3)*shp(1,3) & +x(2,4)*shp(1,4) sx(1,1) = x(2,1)*shp(2,1)+x(2,2)*shp(2,2)+x(2,3)*shp(2,3) & +x(2,4)*shp(2,4) xsj = sx(1,1)*sx(2,2)-sx(1,2)*sx(2,1) xsj1 = xsj if(xsj.eq.0.0d0) xsj1 = 1.0d0 sx(2,2) = sx(2,2)/xsj1 sx(1,1) = sx(1,1)/xsj1 sx(1,2) =-sx(1,2)/xsj1 sx(2,1) =-sx(2,1)/xsj1 c Form global derivatives do i = 1,8 ! { tp = shp(1,i)*sx(1,1)+shp(2,i)*sx(2,1) shp(2,i) = shp(1,i)*sx(1,2)+shp(2,i)*sx(2,2) shp(1,i) = tp end do ! i } c Form rotational and 5-th shape functions call hshp3d(shp(1,5),shp1,shp2) end subroutine stcn3sh(ix,d,yl,ul,tr,dt,st,dvl,ndf,nel,numnp) implicit none include 'shpf16.h' include 'sstr16.h' integer ndf,nel,numnp, i,j,k, ii real*8 xsji, xx,yy,zz integer ix(*) real*8 dt(numnp),st(numnp,*),yl(3,*),tr(3,3),d(*),ul(ndf,*) real*8 vl(6,4),eps(6),momt(6),norm(6),dd(6,6),dvl(4) save c Compute membrane and bending stresses for projection. do i = 1,4 ! { do j = 1,3 ! { vl(j ,i) = 0.0d0 vl(j+3,i) = 0.0d0 do k = 1,3 ! { vl(j ,i) = vl(j ,i) + tr(j,k)*ul(k ,i) vl(j+3,i) = vl(j+3,i) + tr(j,k)*ul(k+3,i) end do ! k } end do ! j } end do ! i } do k = 1,4 ! { call stre3d(d,yl,vl,3,nel,k, xx,yy,zz,eps,norm,momt,dd) do j = 1,4 ! { xsji = dvl(k)*shp(3,j,k) ii = ix(j) if(ii.gt.0) then dt(ii) = dt(ii) + xsji st(ii,1) = st(ii,1) + norm(1)*xsji st(ii,2) = st(ii,2) + norm(2)*xsji st(ii,3) = st(ii,3) + norm(3)*xsji st(ii,4) = st(ii,4) + norm(4)*xsji st(ii,5) = st(ii,5) + norm(5)*xsji st(ii,6) = st(ii,6) + momt(1)*xsji st(ii,7) = st(ii,7) + momt(2)*xsji st(ii,8) = st(ii,8) + momt(3)*xsji st(ii,9) = st(ii,9) + momt(4)*xsji st(ii,10) = st(ii,10) + momt(5)*xsji st(ii,11) = st(ii,11) + eps(1)*xsji st(ii,12) = st(ii,12) + eps(2)*xsji st(ii,13) = st(ii,13) + eps(3)*xsji st(ii,14) = st(ii,14) + eps(4)*xsji st(ii,15) = st(ii,15) + eps(5)*xsji st(ii,16) = st(ii,16) + eps(6)*xsji st(ii,17) = st(ii,17) + sigt(1)*xsji st(ii,18) = st(ii,18) + sigt(2)*xsji st(ii,19) = st(ii,19) + sigt(3)*xsji st(ii,20) = st(ii,20) + sigb(1)*xsji st(ii,21) = st(ii,21) + sigb(2)*xsji st(ii,22) = st(ii,22) + sigb(3)*xsji st(ii,23) = st(ii,23) + max(sigb(4),sigt(4))*xsji st(ii,24) = st(ii,24) + min(sigb(5),sigt(5))*xsji end if end do ! j } end do ! k } end subroutine stre3d(d,xl,vl,ndm,nel,l, xx,yy,zz,eps,norm,momt,dd) implicit none include 'shlc16.h' include 'shld16.h' integer ndm,nel,l, i,j real*8 xx,yy,zz,thk,thk3 real*8 d(*), xl(ndm,*),vl(6,*), eps(6),norm(6),momt(6),temp(4) real*8 dd(6,6),alp(6) include 'shpf16.h' include 'sstr16.h' save c Compute membrane and bending strains if(qdflg) then call dktq3d(shp(1,5,l),shp(1,1,l)) else call dktb3d(sg(1,l),xsjt) endif do i = 1,6 ! { eps(i) = 0.0 end do ! i } xx = 0.0d0 yy = 0.0d0 zz = 0.0d0 do j = 1,nel ! { xx = xx + shp(3,j,l)*xl(1,j) yy = yy + shp(3,j,l)*xl(2,j) zz = zz + shp(3,j,l)*xl(3,j) eps(1) = eps(1) + shp(1,j,l)*vl(1,j) & - shp1(1,j,l)*vl(6,j) eps(2) = eps(2) + shp(2,j,l)*vl(2,j) & - shp2(2,j,l)*vl(6,j) eps(3) = eps(3) + shp(1,j,l)*vl(2,j) + shp(2,j,l)*vl(1,j) & - (shp1(2,j,l) + shp2(1,j,l))*vl(6,j) do i = ii1,ii2 ! { eps(4) = eps(4) + bm(1,i,j)*vl(i,j) eps(5) = eps(5) + bm(2,i,j)*vl(i,j) eps(6) = eps(6) + bm(3,i,j)*vl(i,j) end do ! i } end do ! j } call dmat2d(d,d(31),dd,alp) thk = d(14) thk3 = thk**3/12.d0 norm(1) = (dd(1,1)*eps(1) + dd(1,2)*eps(2) + dd(1,4)*eps(3))*thk norm(2) = (dd(2,1)*eps(1) + dd(2,2)*eps(2) + dd(2,4)*eps(3))*thk norm(3) = (dd(4,1)*eps(1) + dd(4,2)*eps(2) + dd(4,4)*eps(3))*thk temp(1) = norm(1) temp(2) = norm(2) temp(4) = norm(3) call pstr2d(temp,norm(4)) momt(1) = (dd(1,1)*eps(4) + dd(1,2)*eps(5) + dd(1,4)*eps(6))*thk3 momt(2) = (dd(2,1)*eps(4) + dd(2,2)*eps(5) + dd(2,4)*eps(6))*thk3 momt(3) = (dd(4,1)*eps(4) + dd(4,2)*eps(5) + dd(4,4)*eps(6))*thk3 temp(1) = momt(1) temp(2) = momt(2) temp(4) = momt(3) call pstr2d(temp,momt(4)) c Compute surface stresses from bending and in-plane forces thk = 1.d0/thk do i = 1,3 ! { sigt(i) = (norm(i) + 6.d0*momt(i)*thk)*thk sigb(i) = (norm(i) - 6.d0*momt(i)*thk)*thk end do ! i } temp(1) = sigt(1) temp(2) = sigt(2) temp(4) = sigt(3) call pstr2d(temp,sigt(4)) temp(1) = sigb(1) temp(2) = sigb(2) temp(4) = sigb(3) call pstr2d(temp,sigb(4)) end subroutine rayl3d(d,vl,nel,l, norm,momt,dd) c Purpose: Rayleigh damping stresses implicit none include 'shlc16.h' include 'shld16.h' integer nel,l, i,j real*8 thk,thk3 real*8 d(*), vl(6,*), eps(6),norm(6),momt(6), dd(6,6) include 'shpf16.h' include 'sstr16.h' save c Compute membrane and bending strains if(qdflg) then call dktq3d(shp(1,5,l),shp(1,1,l)) else call dktb3d(sg(1,l),xsjt) endif do i = 1,6 ! { eps(i) = 0.0 end do ! i } do j = 1,nel ! { eps(1) = eps(1) + shp(1,j,l)*vl(1,j) & - shp1(1,j,l)*vl(6,j) eps(2) = eps(2) + shp(2,j,l)*vl(2,j) & - shp2(2,j,l)*vl(6,j) eps(3) = eps(3) + shp(1,j,l)*vl(2,j) + shp(2,j,l)*vl(1,j) & - (shp1(2,j,l) + shp2(1,j,l))*vl(6,j) do i = ii1,ii2 ! { eps(4) = eps(4) + bm(1,i,j)*vl(i,j) eps(5) = eps(5) + bm(2,i,j)*vl(i,j) eps(6) = eps(6) + bm(3,i,j)*vl(i,j) end do ! i } end do ! j } thk = d(14)*d(78) thk3 = thk**3/12.d0 norm(1) = norm(1) + (dd(1,i)*eps(i) & + dd(1,2)*eps(2) & + dd(1,4)*eps(3))*thk norm(2) = norm(1) + (dd(2,i)*eps(1) & + dd(2,2)*eps(2) & + dd(2,4)*eps(3))*thk norm(3) = norm(1) + (dd(4,i)*eps(1) & + dd(4,2)*eps(2) & + dd(4,4)*eps(3))*thk momt(1) = momt(1) + (dd(1,1)*eps(4) & + dd(1,2)*eps(5) & + dd(1,4)*eps(6))*thk3 momt(2) = momt(2) + (dd(2,1)*eps(4) & + dd(2,2)*eps(5) & + dd(2,4)*eps(6))*thk3 momt(3) = momt(3) + (dd(4,1)*eps(4) & + dd(4,2)*eps(5) & + dd(4,4)*eps(6))*thk3 end subroutine tran3d(xl,yl,t,ndm) c Compute transformation array and surface coords. implicit none integer ndm, i,j,k real*8 v1, v2, htol, dl1,dl2 real*8 x0(3),xl(ndm,*),yl(3,4),t(3,3) c Compute inplane direction cosines (bisect diagonals) do i = 1,3 ! { t(1,i) = xl(i,3) - xl(i,1) t(2,i) = xl(i,2) - xl(i,4) end do ! i } dl1 = sqrt(t(1,1)**2 + t(1,2)**2 + t(1,3)**2) dl2 = sqrt(t(2,1)**2 + t(2,2)**2 + t(2,3)**2) do i = 1,3 ! { v1 = t(1,i)/dl1 v2 = t(2,i)/dl2 t(1,i) = v1 + v2 t(2,i) = v1 - v2 end do ! i } dl1 = sqrt(t(1,1)**2 + t(1,2)**2 + t(1,3)**2) dl2 = sqrt(t(2,1)**2 + t(2,2)**2 + t(2,3)**2) do i = 1,3 ! { t(1,i) = t(1,i)/dl1 t(2,i) = t(2,i)/dl2 c Compute center (0,0) displacement x0(i) = 0.25*(xl(i,1) + xl(i,2) + xl(i,3) + xl(i,4)) end do ! i } c Compute normal to surface t(3,1) = t(1,2)*t(2,3) - t(2,2)*t(1,3) t(3,2) = t(1,3)*t(2,1) - t(2,3)*t(1,1) t(3,3) = t(1,1)*t(2,2) - t(2,1)*t(1,2) c Compute projected middle surface coordinates do i = 1,4 ! { do j = 1,3 ! { yl(j,i) = 0.0 do k = 1,3 ! { yl(j,i) = yl(j,i) + t(j,k)*(xl(k,i) - x0(k)) end do ! k } end do ! j } end do ! i } c Set offset coordinates to zero if small compared to plan size htol = 0.0 do i = 1,4 ! { htol = max(htol,abs(yl(1,i)),abs(yl(2,i))) end do ! i } htol = htol*1.d-7 do i = 1,4 ! { if(abs(yl(3,i)) .le. htol) yl(3,i) = 0.0d0 end do ! i } end