invar.f

关于网格自动生成delaunay算法的

*********************************************************************
*     Evaluates the current value of the yeild stress function      
*
*********************************************************************
      subroutine invar(ncrit,props,stres,
     $     devia,sint3,steff,theta,varj2,yield)
      implicit real*8(a-h,o-z)
      dimension devia(4),props(7),stres(4),ncrit(3)
      parameter(ragian= 1.745329251994327d-02)
      parameter(root3= 1.732050807569d0)
      smean= (stres(1)+stres(2)+stres(4))/3.0d0
      devia(1)= stres(1)-smean
      devia(2)= stres(2)-smean
      devia(3)= stres(3)
      devia(4)= stres(4)-smean
      varj2= devia(3)*devia(3) + 0.5d0*
     $     (devia(1)*devia(1)+devia(2)*devia(2)+devia(4)*devia(4))
      varj3= devia(4)*(devia(4)*devia(4)-varj2)
      steff= sqrt(varj2)
      if(ncrit(1).ne.2) then
         if(steff.lt.1.0e-30) then
            sint3= 0.0d0
         else
            sint3= -3.0d0*root3*varj3/(2.0d0*varj2*steff)
            if(sint3.gt.1.0d0) then
               sint3=  1.0d0
            elseif(sint3.lt.-1.0d0) then
               sint3= -1.0d0
            endif
         endif
         theta= asin(sint3)/3.0d0
      endif
      if(ncrit(1).eq.2) then
         yield= root3*steff
      elseif(ncrit(1).eq.1) then
         yield= 2.0d0*cos(theta)*steff
      elseif(ncrit(1).eq.3) then
         phira= props(6)*ragian
         snphi= sin(phira)
         yield= smean*snphi+steff*(cos(theta)-sin(theta)*snphi/root3)
      elseif(ncrit(1).eq.4) then
         phira= props(6)*ragian
         snphi= sin(phira)
         yield= 6.0d0*smean*snphi/(root3*(3.0d0-snphi))+steff
      endif
      return
      end
*
**********************************************************************
*     evaluates the flow vector                                      
**********************************************************************
      subroutine yieldf(ncrit,props,nstr1,
     $     devia,steff,theta,varj2,avect)
      implicit real*8(a-h,o-z)
      parameter(pi    = 3.14159265358979d+00)
      parameter(ragian= 1.74532925199433d-02)
      dimension avect(4),devia(4),props(7),ncrit(3)
      dimension veca1(4),veca2(4),veca3(4)
      if(steff.lt.1.0d-30) then
         return
      endif
      root3= sqrt(3.0d0)
      frict= props(6)
      tanth= tan(theta)
      tant3= tan(3.0d0*theta)
      sinth= sin(theta)
      costh= cos(theta)
      cost3= cos(3.0d0*theta)
      veca1(1)= 1.0d0
      veca1(2)= 1.0d0
      veca1(3)= 0.0d0
      veca1(4)= 1.0d0
      do 10 istr1= 1,nstr1
         veca2(istr1)= devia(istr1)/(2.0d0*steff)
 10   continue
      veca2(3)= devia(3)/steff
      veca3(1)= devia(2)*devia(4) + varj2/3.0d0
      veca3(2)= devia(1)*devia(4) + varj2/3.0d0
      veca3(3)= -2.0d0*devia(3)*devia(4)
      veca3(4)= devia(1)*devia(2)-devia(3)*devia(3)+varj2/3.0d0
      if(ncrit(1).eq.1) then
         cons1= 0.0d0
         abthe= abs(theta/pi*180.0d0)
         if(abthe.gt.29.0d0) then
            cons2= root3
            cons3= 0.0d0
         else
            cons2= 2.0d0*(costh+sinth*tant3)
            cons3= root3*sinth/(varj2*cost3)
         endif
      elseif(ncrit(1).eq.2) then
         cons1= 0.0d0
         cons2= root3
         cons3= 0.0d0
      elseif(ncrit(1).eq.3) then
         cons1= sin(frict*ragian)/3.0d0
         abthe= abs(theta/pi*180.0d0)
         if(abthe.gt.29.0d0) then
            cons3= 0.0d0
            if(theta.gt.0.0d0) then
               plumi=-1.0d0
            else
               plumi= 1.0d0
            endif
            cons2= 0.5d0*(root3+plumi*cons1*root3)
         else
            cons2= costh*((1.0d0+tanth*tant3)+
     $           cons1*(tant3-tanth)*root3)
            cons3= (root3*sinth+3.0d0*cons1*costh)/(2.0d0*varj2*cost3)
         endif
      elseif(ncrit(1).eq.4) then
         snphi= sin(frict*ragian)
         cons1= 2.0d0*snphi/(root3*(3.0d0-snphi))
         cons2= 1.0d0
         cons3= 0.0d0
      endif
      do 50 istr1= 1,nstr1
         avect(istr1)= cons1*veca1(istr1) +
     $                 cons2*veca2(istr1) +
     $                 cons3*veca3(istr1)
 50   continue
      return
      end
*
*********************************************************************
*     plastic d vector    dvector= d(elastic)*avect                 
*********************************************************************
      subroutine flowpl(ncrit,props,ntype,hdpara,nstr1,gepst,
     $     avect,abeta,dvect)
      implicit real*8(a-h,o-z)
      dimension avect(4),dvect(4),props(7)
      dimension ncrit(3),hdpara(20)
      eel= props(1)
      por= props(2)
      para= eel/(1.0d0+por)
      if(ntype.eq.1) then
         coeff= eel*por*(avect(1)+avect(2))/(1.0d0-por*por)
      else
         coeff= eel*por*(avect(1)+avect(2)+avect(4))
     $        /((1.0d0+por)*(1.0d0-2.0d0*por))
      endif
      dvect(1)= para*avect(1) + coeff
      dvect(2)= para*avect(2) + coeff
      dvect(3)= 0.5d0*avect(3)*eel/(1.0d0+por)
      if(ntype.eq.2) then
         dvect(4)= para*avect(4) + coeff
      else
         dvect(4)= 0.0d0
      endif
      denom= hard(ncrit,hdpara,gepst)
      do istr1= 1,nstr1
         denom= denom + avect(istr1)*dvect(istr1)
      enddo
      abeta= 1.0d0/denom
      return
      end
*
*********************************************************************
*     estimate the rate of strain-hardening                         
*********************************************************************
      real*8 function hard(ncrit,hdpara,gepst)
      implicit real*8(a-h,o-z)
      integer ncrit(3)
      real*8 hdpara(20)
      hardid= ncrit(2)
      if(hardid.eq.1) then
         hard= hdpara(1)
      elseif(hardid.eq.2) then
         nline= ncrit(3)
         do 10 i=1,nline
            iline= 2*(i-1)+1
            if(gepst.lt.hdpara(iline)) then
               hard= hdpara(iline+1)
               goto 20
            endif
 10      continue
         hard= hdpara(nline*2)
 20      continue
      elseif(hardid.eq.3) then
         sigmy= hdpara(1)
         strny= hdpara(2)
         alpha= hdpara(3)
         ann  = hdpara(4)
         if(gepst.lt.1.0d-6) then
            hard= ((1.0d-6/(alpha*strny))**(1.0d0/ann))*
     $           sigmy/(ann*1.0d-6)
         else
            hard= ((gepst/(alpha*strny))**(1.0d0/ann))*
     $           sigmy/(ann*gepst)
         endif
      elseif(hardid.eq.4) then
         sigmy= hdpara(1)
         strny= hdpara(2)
         alpha= hdpara(3)
         ann  = hdpara(4)
         hard= ((gepst/(alpha*strny)+1.0d0)**(1.0d0/ann))*
     $           sigmy/(ann*(gepst+sigmy*alpha))
      endif
      return
      end