residu.f

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

**********************************************************************
*     This subroutine reduces the stresses to the yield surface and
*     evaluates the equivalent nodal forces (for domain region)
**********************************************************************
      subroutine domresidu(mpoin,mnod,mgpf,mavnd,ng,nszf,
     $     ntype,ngpf,x,tu,asdis,bsdis,eload,effstf,disp,iflnod,
     $     xxgpf,wwgpf,dmaxf,ndgpf,idgpf,dldx,bl,bmat,eforc,
     $     alcm,kk,xi,kfx,epstnf,strsgf,epsgf,props,ncrit,hdpara,
     $     epstnfpre)
      implicit real*8(a-h,o-z)
      parameter(radian= 1.74532925199433d-02)
      dimension  x(ng,mpoin+3)
      dimension  xxgpf(ng,mgpf), wwgpf(mgpf), dmaxf(mgpf)
      dimension  ndgpf(0:mgpf),idgpf(mgpf*mavnd),iflnod(mnod)
      dimension  asdis(mpoin*ng),bsdis(mpoin*ng)
      dimension  eload(mpoin*ng),tu(mpoin*ng)
      dimension  effstf(mgpf),epstnf(mgpf),strsgf(4,mgpf),epsgf(4,mgpf)
      dimension  props(7),ncrit(3),hdpara(20)
      dimension  aa(6,6),ai(6,6), bl(6,mnod)
      dimension  dadx(6,6,2), dldx(6,mnod,ng)
      dimension  dmat(4,4),xx(2),xi(ng,mnod)
      dimension  bmat(3,mnod*ng)
      dimension  kfx(mpoin)
      dimension  avect(4),devia(4),dvect(4)
      dimension  stres(4),desig(4),sigma(4),sgtot(4)
      dimension  eps(4)
      dimension  disp(mnod*ng),eforc(mnod*ng)
      dimension  epstnfpre(mgpf)
      nstre= 3
      nstr1= nstre+1
      md = 3
      do iszf=1,nszf
         asdis(iszf)= eload(iszf)
         bsdis(iszf)= bsdis(iszf) + eload(iszf)
         tu(iszf)= tu(iszf)+eload(iszf)
         eload(iszf)= 0.0d0
      enddo
      do 20 igpf= 1,ngpf
         epstnfpre(igpf) = epstnf(igpf)
         eel= props(1)
         por= props(2)
         thick= props(3)
         if(ncrit(1).eq.1.or.ncrit(1).eq.2) then
            uniax= props(5)
         elseif(ncrit(1).eq.3) then
            frict= props(6)
            uniax= props(5)*cos(frict*radian)
         elseif(ncrit(1).eq.4) then
            frict= props(6)
            uniax= 6.0d0*props(5)*cos(frict*radian)/
     $           (sqrt(3.0d0)*(3.0d0-sin(frict*radian)))
         endif
         do ig= 1,ng
            xx(ig)= xxgpf(ig,igpf)
         enddo
         www= wwgpf(igpf)
         ddmax= dmaxf(igpf)
         nnod= ndgpf(igpf)-ndgpf(igpf-1)
         do inod= 1,nnod
            iflnod(inod)= idgpf(ndgpf(igpf-1)+inod)
            ipoin= iflnod(inod)
            do ig= 1,ng
               xi(ig,inod)= x(ig,ipoin) - xx(ig)
            enddo
         enddo
         dm= ddmax
         cm= alcm*dm
         call makeaa(ng,mnod,md,nnod,xi,dm,cm,
     $        kk,aa,bl,ai,dadx,dldx)
         call bmatrx(ng,mnod,md,nnod,bl,ai,dadx,dldx,bmat)
         call dmatrx(ntype,props,nstrs,dmat)
         do i = 1,mnod*2
            disp(i) = 0.0d0
         enddo
         do i = 1,nnod
            disp(i*2-1) = asdis(iflnod(i)*2-1)
            disp(i*2) = asdis(iflnod(i)*2)
         enddo
         do i=1,4
            eps(i) =0.0d0
            stres(i) =0.0d0
         enddo
         do i = 1,nstre
            do inod= 1,nnod*ng
               eps(i) = eps(i) + bmat(i,inod)*disp(inod)
            enddo
         enddo
         do j = 1,nstre
            do i = 1,nstre
               stres(j) = stres(j) + dmat(j,i)*eps(i)
            enddo
         enddo
         if(ntype.eq.1) then
            eps(4) = - por/ eel * (stres(1)+stres(2))
            stres(4) = 0.0d0
         else
            eps(4) = 0.0d0
            stres(4) = por*(stres(1)+stres(2))
         endif
         do i= 1,4
            epsgf(i,igpf)= epsgf(i,igpf) + eps(i)
         enddo
         gepst= epstnf(igpf)
         preys= yslmt(ncrit,hdpara,gepst,uniax,eel)
         do istr1= 1,nstr1
            desig(istr1)= stres(istr1)
            sigma(istr1)= strsgf(istr1,igpf)+stres(istr1)
         enddo
         call invar(ncrit,props,sigma,
     $        devia,sint3,steff,theta,varj2,yield)
         espre= effstf(igpf) - preys
         if(espre.lt.0.0d0) then
            escur= yield-preys
            if(escur.gt.0.0d0) then
               rfact= escur/(yield-effstf(igpf))
               iredu= 1
            else
               iredu= 0
            endif
         else
            escur= yield-effstf(igpf)
            if(escur.gt.0.0d0) then
               rfact= 1.0d0
               iredu= 1
            else
               iredu= 0
            endif
         endif
         if(iredu.eq.1) then
            mstep= escur*8.0d0/uniax+1.0d0
            if(mstep.gt.50) then
               mstep=50
            endif
            astep= dble(mstep)
            reduc= 1.0d0-rfact
            do istr1= 1,nstr1
               sgtot(istr1)= strsgf(istr1,igpf) + 
     $              reduc*stres(istr1)
               stres(istr1)= rfact*stres(istr1)/astep
            enddo
            do 90 istep= 1,mstep
               call invar(ncrit,props,sgtot,
     $              devia,sint3,steff,theta,varj2,yield)
               call yieldf(ncrit,props,nstr1,
     $              devia,steff,theta,varj2,avect)
               call flowpl(ncrit,props,ntype,hdpara,nstr1,
     $              gepst,avect,abeta,dvect)
               agash= 0.0d0
               do istr1= 1,nstr1
                  agash= agash+avect(istr1)*stres(istr1)
               enddo
               dlamd= agash*abeta
               if(dlamd.lt.0.0d0) then
                  dlamd= 0.0d0
               endif
               bgash= 0.0d0
               do istr1= 1,nstr1
                  bgash= bgash+avect(istr1)*sgtot(istr1)
                  sgtot(istr1)= sgtot(istr1) + 
     $                 stres(istr1)-dlamd*dvect(istr1)
               enddo
               epstnf(igpf)= epstnf(igpf)+dlamd*bgash/yield
               gepst= epstnf(igpf)
 90         continue 
            call invar(ncrit,props,sgtot,
     $           devia,sint3,steff,theta,varj2,yield)
            curys= yslmt(ncrit,hdpara,gepst,uniax,eel)
            bring= 1.0d0
            if(yield.gt.curys) then
               bring= curys/yield
            endif
            do istr1= 1,nstr1
               strsgf(istr1,igpf)= bring*sgtot(istr1)
            enddo
            effstf(igpf)= bring*yield
         else
            do istr1= 1,nstr1
               strsgf(istr1,igpf)= strsgf(istr1,igpf)+desig(istr1)
            enddo
            effstf(igpf)= yield
         endif
         do i=1,nnod*ng
            eforc(i) =0.0d0
         enddo
         do inod=1,nnod*ng
            do istre= 1,nstre
               eforc(inod)= eforc(inod) + www*
     $              bmat(istre,inod)*strsgf(istre,igpf)
            enddo
         enddo
         do i = 1,nnod
            iflnod(i) = kfx(iflnod(i))
         enddo
         call formtf(ng,mpoin,mnod,nnod,iflnod,
     $        eforc,eload)
 20   continue
      return
      end
*
**********************************************************************
*     This subroutine reduces the stresses to the yield surface and
*     evaluates the equivalent nodal forces (for boundary region)
**********************************************************************
      subroutine bunresidu(mpoin,mnod,mgpb,mavnd,ng,
     $     ntype,ngpb,x,effstb,disp,iflnod,xxgpb,
     $     asdis,dmaxb,ndgpb,idgpb,dldx,bl,bmat,
     $     alcm,kk,xi,epstnb,strsgb,props,ncrit,hdpara,
     $     epstnbpre)
      implicit real*8(a-h,o-z)
      parameter(radian= 1.74532925199433d-02)
      dimension  x(ng,mpoin+3)
      dimension  xxgpb(ng,mgpb), dmaxb(mgpb)
      dimension  ndgpb(0:mgpb),idgpb(mgpb*mavnd),iflnod(mnod)
      dimension  effstb(mgpb),epstnb(mgpb),strsgb(4,mgpb)
      dimension  props(7),ncrit(3),hdpara(20)
      dimension  aa(6,6),ai(6,6), bl(6,mnod)
      dimension  dadx(6,6,2), dldx(6,mnod,ng)
      dimension  dmat(4,4),xx(2),xi(ng,mnod)
      dimension  bmat(3,mnod*ng)
      dimension  disp(mnod*ng),asdis(mpoin*ng)
      dimension  avect(4),devia(4),dvect(4)
      dimension  stres(4),desig(4),sigma(4),sgtot(4)
      dimension  epstnbpre(mgpb)
      nstre= 3
      nstr1= nstre+1
      md = 3
      do 20 igpb= 1,ngpb
         epstnbpre(igpb) = epstnb(igpb)
         eel= props(1)
         por= props(2)
         thick= props(3)
         if(ncrit(1).eq.1.or.ncrit(1).eq.2) then
            uniax= props(5)
         elseif(ncrit(1).eq.3) then
            frict= props(6)
            uniax= props(5)*cos(frict*radian)
         elseif(ncrit(1).eq.4) then
            frict= props(6)
            uniax= 6.0d0*props(5)*cos(frict*radian)/
     $           (sqrt(3.0d0)*(3.0d0-sin(frict*radian)))
         endif
*     
*---- calculate incremental elastic stress 
*     
         do ig= 1,ng
            xx(ig)= xxgpb(ig,igpb)
         enddo
         ddmax= dmaxb(igpb)
         nnod= ndgpb(igpb)-ndgpb(igpb-1)
         do inod= 1,nnod
            iflnod(inod)= idgpb(ndgpb(igpb-1)+inod)
            ipoin= iflnod(inod)
            do ig= 1,ng
               xi(ig,inod)= x(ig,ipoin) - xx(ig)
            enddo
         enddo
         dm= ddmax
         cm= alcm*dm
         call makeaa(ng,mnod,md,nnod,xi,dm,cm,
     $        kk,aa,bl,ai,dadx,dldx)
         call bmatrx(ng,mnod,md,nnod,bl,ai,dadx,dldx,bmat)
         call dmatrx(ntype,props,nstrs,dmat)
         do i = 1,mnod*2
            disp(i) = 0.0d0
         enddo
         do i = 1,nnod
            disp(i*2-1) = asdis(iflnod(i)*2-1)
            disp(i*2) = asdis(iflnod(i)*2)
         enddo
         do i=1,4
            stres(i) =0.0d0
         enddo
         do j = 1,nstre
            do i = 1,nstre
               do inod= 1,nnod*ng
                  stres(j) = stres(j) + 
     $                 dmat(j,i)*bmat(i,inod)*disp(inod)
               enddo
            enddo
         enddo
         if(ntype.eq.1) then
            stres(4) = 0.0d0
         else
            stres(4) = por*(stres(1)+stres(2))
         endif
         gepst= epstnb(igpb)
         preys= yslmt(ncrit,hdpara,gepst,uniax,eel)
         do istr1= 1,nstr1
            desig(istr1)= stres(istr1)
            sigma(istr1)= strsgb(istr1,igpb)+stres(istr1)
         enddo
         call invar(ncrit,props,sigma,
     $        devia,sint3,steff,theta,varj2,yield)
         espre= effstb(igpb)-preys
         if(espre.le.0.0d0) then
*---- elastic behavior at pre-step
            escur= yield-preys
            if(escur.gt.0.0d0) then
*----- change to plastic behavior at this step
               rfact= escur/(yield-effstb(igpb))
               iredu= 1
            else
*---- continue elastic behavior
               iredu= 0
            endif
         else
*---- plastic behavior at pre-step
            escur= yield-effstb(igpb)
            if(escur.gt.0.0d0) then
*---- continue plastic behavior
               rfact= 1.0d0
               iredu= 1
            else
*---- change to elastic behavior at this step
               iredu= 0
            endif
         endif
*---- reduce the stress to the yield surface
         if(iredu.eq.1) then
            mstep= escur*8.0d0/uniax+1.0d0
            if(mstep.gt.50) then
               mstep=50
            endif
            mstep= 1
            astep= dble(mstep)
            reduc= 1.0d0-rfact
            do istr1= 1,nstr1
               sgtot(istr1)= strsgb(istr1,igpb) + 
     $              reduc*stres(istr1)
               stres(istr1)= rfact*stres(istr1)/astep
            enddo
            do 90 istep= 1,mstep
               call invar(ncrit,props,sgtot,
     $              devia,sint3,steff,theta,varj2,yield)
               call yieldf(ncrit,props,nstr1,
     $              devia,steff,theta,varj2,avect)
               call flowpl(ncrit,props,ntype,hdpara,nstr1,gepst,
     $              avect,abeta,dvect)
               agash= 0.0d0
               do istr1= 1,nstr1
                  agash= agash+avect(istr1)*stres(istr1)
               enddo
               dlamd= agash*abeta
               if(dlamd.lt.0.0d0) then
                  dlamd= 0.0d0
               endif
               bgash= 0.0d0
               do istr1= 1,nstr1
                  bgash= bgash+avect(istr1)*sgtot(istr1)
                  sgtot(istr1)= sgtot(istr1) + 
     $                 stres(istr1)-dlamd*dvect(istr1)
               enddo
               epstnb(igpb)= epstnb(igpb)+dlamd*bgash/yield
               gepst= epstnb(igpb)
 90         continue 
            call invar(ncrit,props,sgtot,
     $           devia,sint3,steff,theta,varj2,yield)
            curys= yslmt(ncrit,hdpara,gepst,uniax,eel)
            bring= 1.0d0
            if(yield.gt.curys) then
               bring= curys/yield
            endif
            do istr1= 1,nstr1
               strsgb(istr1,igpb)= bring*sgtot(istr1)
            enddo
            effstb(igpb)= bring*yield
         else
            do istr1= 1,nstr1
               strsgb(istr1,igpb)= strsgb(istr1,igpb)+desig(istr1)
            enddo
            effstb(igpb)= yield
         endif
*     
*------- end of gauss integral loop
*
 20   continue
      return
      end
*
*********************************************************************
*   estimation of yield surface
*********************************************************************
      real*8 function yslmt(ncrit,hdpara,gepst,uniax,eel)
      implicit real*8(a-h,o-z)
      dimension ncrit(3),hdpara(20)
*
*----- bi-linear
*
      nhardid= ncrit(2)
      if(nhardid.eq.1) then
         yslmt= uniax+hdpara(1)*gepst
      elseif(nhardid.eq.2) then
         nline= ncrit(3)
         yslmt= uniax
         ysnlt= uniax/eel
         do 10 i=1,nline
            iline= 2*(i-1)+1
            if(gepst.lt.hdpara(iline)) then
               yslmt= yslmt + hdpara(iline+1)*(gepst-ysnlt)
               goto 20
            else
               yslmt= yslmt + hdpara(iline+1)*
     $              (hdpara(iline)-ysnlt)
               ysnlt= hdpara(iline)
            endif
 10      continue
 20      continue
      elseif(nhardid.eq.3) then
         sigmy= hdpara(1)
         strny= hdpara(2)
         alpha= hdpara(3)
         ann  = hdpara(4)
         gpslmt= alpha*strny*((uniax/sigmy)**ann)
         if(gepst.lt.gpslmt) then
            yslmt= uniax
         else
            yslmt= sigmy*((gepst/alpha/strny)**(1.0D0/ann))
         endif
      elseif(nhardid.eq.4) then
         sigmy= hdpara(1)
         strny= hdpara(2)
         alpha= hdpara(3)
         ann  = hdpara(4)
         yslmt= sigmy*((1.0d0+gepst/alpha/strny)**(1.0D0/ann))
      endif
      return
      end