Télécharger la plaquette Cast3M


fr en

Télécharger denc3d.eso

Retour à la liste

Numérotation des lignes : 

denc3d

  1. C DENC3D    SOURCE    FD218221  24/02/07    21:15:09     11834          
  2.         subroutine denc3d(ngf,NDIMG,ndima,NTMAX,NCMAX,NLMAX,
  3.      #  nbract,ft,lnumcrt,dpftds,dgftds,REST,PFT_GRAD_PFT,PFT_GRAD_GFT,
  4.      #  nbracc,fc,lnumcrc,dpfcds,dgfcds,RESC,PFC_GRAD_PFC,PFC_GRAD_GFC,
  5.      #  nbracl,fl,lnumcrl,dpflds,dgflds,RESL,PFL_GRAD_PFL,PFL_GRAD_GFL,
  6.      #  jse,jea,aai,aap,jsp,Je1epg,DEPST1,DEPSL1,DEPSC1,reduc_resg,
  7.      #  NDIMPL,RgR,ninc,NBRINC,Jsgep,ACTIFT,ACTIFC,ACTIFL,PFG_GRAD_PFG,
  8.      #  precision3d,err1,affiche)
  9.  
  10. c       calcul des ecoulements dans l espace des deformations plastiques
  11.  
  12.         implicit real*8 (a-h,o-z)
  13.         implicit integer (i-n)
  14.  
  15. c       dimension espace generalise
  16.         integer NDIMG,ndima,nbract,err1,ngf,NTMAX,NCMAX,NLMAX
  17.         integer nbracc,nbracl,ninc,NBRINC
  18. c       activite des criteres      
  19.         logical ACTIFT(NTMAX),ACTIFC(NCMAX),ACTIFL(NLMAX)
  20.         real*8 aai(ngf,ngf+1)
  21.  
  22. c       criteres de traction 
  23.         real*8 ft(NTMAX),rest
  24. c       criteres de cisaillement
  25.         real*8 fc(NCMAX),resc 
  26. c       traction localisee
  27.         real*8 fl(NCMAX),resl       
  28.         integer lnumcrt(NTMAX),lnumcrc(NCMAX),lnumcrl(NLMAX)
  29. c       inverse de a
  30.         real*8 jea(NDIMG,ndima)
  31. c       regiddite generalisee
  32.         real*8 jse(NDIMG,NDIMG)
  33. c       derivee du critere par rapport à la contrainte
  34.         real*8 DPFTDS(NTMAX,NDIMG),DPFCDS(NCMAX,NDIMG)
  35.         real*8 DPFLDS(NLMAX,NDIMG)  
  36. c       derivee du critere par rapport à la contrainte
  37.         real*8 DGFTDS(NTMAX,NDIMG),DGFCDS(NCMAX,NDIMG) 
  38.         real*8 DGFLDS(NLMAX,NDIMG)
  39.  
  40. c       dimension de l espace de minimisation des residus plastiques
  41.         integer NDIMPL
  42. c       jacobienne Je1epg dans la base des ecoulements plastiques
  43.         real*8 Je1epg(NDIMG,NDIMPL) 
  44. c       matrice jsp=dseff/depsp generalise
  45.         real*8 jsp(NDIMG,NDIMPL)
  46. c       preparation de la matrice inv(A).Jea
  47.         real*8 AAP(NDIMG,NDIMPL)
  48.  
  49. c       gradient des fonctions seuil et pseudo potentiels, direction ecoulement
  50.         real*8 PFT_GRAD_PFT(NDIMPL),PFT_GRAD_GFT(NDIMPL) 
  51. c       gradient des fonctions seuil et pseudo potentiels, direction ecoulement
  52.         real*8 PFC_GRAD_PFC(NDIMPL),PFC_GRAD_GFC(NDIMPL)
  53. c       gradient des fonctions seuil et pseudo potentiels, direction ecoulement
  54.         real*8 PFL_GRAD_PFL(NDIMPL),PFL_GRAD_GFL(NDIMPL)
  55. c       par de fgradf due a la pression de gel pour les criteres de decollement
  56.         real*8 PFG_GRAD_PFG(NDIMPL)
  57. c       precision3
  58.         real*8 precision3d
  59. c       def plastique
  60.         real*8 DEPST1(NDIMG),DEPSC1(NDIMG),DEPSL1(NDIMG)
  61. c       coeff de relaxation du residu
  62.         real*8 reduc_resg
  63. c       matrice d bgpg / d ea ds l espace de minimisation
  64.         real*8 Jsgep(NDIMG,NDIMPL)
  65. c       Vecteur R.grad(R) dans l espace de minimisation
  66.         real*8 RgR(NDIMPL)       
  67.  
  68. c       controle affichage        
  69.         logical affiche,affiche_local
  70.         integer icr,jcr,idir,ityp
  71.  
  72. c       calcul de lambda
  73.         real*8 lambda,numerateur,denom
  74.  
  75.         affiche_local=affiche
  76. c        affiche_local=.true.
  77.         if(affiche_local) then
  78.             print*,'Dans denc3d'
  79.         end if
  80.  
  81. c----------------- preparation des matrice  pour le retour radial ------
  82.  
  83. c       on est en base generali des deformations plastique (3*12=32 ep)
  84. c       01-12: d eps pt
  85. c       13-24: d eps pc
  86. c       25-36: d eps pl
  87. c       avec un seule inclusion NDIMG=12
  88. c       et avec 3 types de plasticite NDIMPL=3*NDIMG=36
  89. c       l indice g est ici relatif a la base de dimansion NDIMPL
  90.  
  91. c       ---------- calcul de la matrice AAp=inv(1-Je1ee)*Je1ep(g) ------
  92.  
  93. c       AAp=d epse_e(epsg) / d eps_p(ecoul)
  94.         do i=1,NDIMG
  95.             do j=1,NDIMPL
  96.                 AAp(i,j)=0.d0
  97.                 do k=1,NDIMG
  98.                      AAp(i,j)=AAp(i,j)+AAi(i,k)*Je1epg(k,j)
  99.                 end do
  100. c               print*,i,j,aap(i,j)               
  101.             end do
  102.         end do  
  103.  
  104. c       ---------matrice Jsp=H*AAp -------------------------------------
  105.  
  106. c       Jsp= d sigma(epsg)/ d eps_p(ecoul)
  107.         do i=1,NDIMG
  108.             do j=1,NDIMPL
  109.                 Jsp(i,j)=0.d0
  110.                 do k=1,NDIMG
  111.                     Jsp(i,j)=Jsp(i,j)+Jse(i,k)*AAp(k,j)
  112.                 end do
  113.             end do
  114.         end do
  115.  
  116. c----------------calcul des gradient ponderes elementaires ------------- 
  117.  
  118. c       les gradients  sont calcules ici dans l espace des variables
  119. c       d ecoulement plastique de dimension = NBTYP_PLAST*NDIMG
  120. 10      do i=1,NDIMPL
  121. c           initialisation du f*gradient de f
  122.             RgR(i)=0.d0
  123.         end do   
  124.  
  125. c       --------- criteres de traction diffuse -------------------------
  126.         rest=0.d0       
  127.         if(nbract.ne.0) then
  128.             do jcr=1,nbract
  129.                 icr=lnumcrt(jcr) 
  130.                 rest=rest+FT(icr)**2
  131.                 call fgrf3d(NDIMG,NDIMPL,NTMAX,DPFTDS,icr,JSP,FT,
  132.      #          PFT_GRAD_PFT,affiche_local)
  133. c               actualisation du gradient complet  
  134.                 do idir=1,NDIMPL
  135.                     RgR(idir)=RgR(idir)+PFT_GRAD_PFT(idir)
  136.                 end do                   
  137.             end do
  138.         end if
  139.  
  140. c       --------- criteres de cisaillement -----------------------------
  141.         resc=0.d0        
  142.         if(nbracc.ne.0) then
  143.             do jcr=1,nbracc
  144.                 icr=lnumcrc(jcr) 
  145.                 resc=resc+FC(icr)**2
  146.                 call fgrf3d(NDIMG,NDIMPL,NCMAX,DPFCDS,icr,JSP,FC,
  147.      #          PFC_GRAD_PFC,affiche_local)
  148. c               actualisation du gradient complet  
  149.                 do idir=1,NDIMPL
  150.                     RgR(idir)=RgR(idir)+PFC_GRAD_PFC(idir)
  151.                 end do                   
  152.             end do
  153.         end if
  154.  
  155. c       --------- criteres de traction decollement localise ------------
  156.         resl=0.d0
  157.         if(nbracl.ne.0) then
  158.             do jcr=1,nbracl
  159.                 icr=lnumcrl(jcr) 
  160.                 resl=resl+FL(icr)**2
  161. c               actualisation de la part de gradient du a sigma'
  162.                 call fgrf3d(NDIMG,NDIMPL,NLMAX,DPFLDS,icr,JSP,FL,
  163.      #          PFL_GRAD_PFL,affiche_local)
  164. c               mise a jour du gradient global dans l espace des 
  165.                 do idir=1,NDIMPL
  166.                     RgR(idir)=RgR(idir)+PFL_GRAD_PFL(idir)
  167.                 end do 
  168. c               actualisation de la part de gradient du a digma gel
  169.                 call fgrf3d(NDIMG,NDIMPL,NLMAX,DPFLDS,icr,JSGEP,FL,
  170.      #          PFG_GRAD_PFG,affiche_local)
  171. c               mise a jour du gradient global dans l espace des 
  172.                 do idir=1,NDIMPL
  173.                     RgR(idir)=RgR(idir)+PFG_GRAD_PFG(idir)
  174.                 end do 
  175.             end do
  176.         end if
  177.  
  178. c------------ multiplicateur plastique ---------------------------------
  179.  
  180. c       ---denominateur-------------------------------------------------
  181.  
  182. c       initialisation du produit scalaire fgradf*fgradF
  183.         denom=0.d0
  184.  
  185. c       ---- prise en compte des ecoulements de traction diffuse -------
  186.         if(nbract.ne.0) then
  187.             do i=1,NBRACT
  188.                 icr=lnumcrt(i) 
  189. c               on limite le produit sur les ndimg du critere de traction 
  190.                 if(ninc.eq.1) then
  191.                     idebut=0
  192.                 else
  193. c                   pour les idebut des autres inclusions il faut
  194. c                   pouvoir calculer idebut
  195.                     print*,'Il manque un traitement  ninc>1 ds denc3d'
  196.                     err1=1
  197.                     return
  198.                 end if
  199. c               contribution de ce critere sur le denominateur
  200.                 do idir=1,NDIMG
  201.                     denom=denom+
  202.      #              RgR(idebut+idir)*DGFTDS(icr,idir)*FT(icr)
  203.                 end do 
  204.             end do
  205.         end if
  206.  
  207. c       ---- prise en compte des ecoulements de cisaillement diffuse ---
  208.         if(nbracc.ne.0) then 
  209.             do i=1,nbracc
  210.                 icr=lnumcrc(i)
  211.                 if(ninc.eq.1) then
  212.                     idebut=NDIMG
  213.                 else
  214. c                   pour les idebut des autres inclusions il faut
  215. c                   pouvoir calculer idebut
  216.                     print*,'Il manque un traitement  ninc>1 ds denc3d'
  217.                     err1=1
  218.                     return
  219.                 end if
  220. c               contribution de ce critere sur le denominateur
  221.                 do idir=1,NDIMG
  222.                     denom=denom+
  223.      #              RgR(idebut+idir)*DGFCDS(icr,idir)*FC(icr)
  224.                 end do             
  225.             end do
  226.         end if
  227.  
  228. c       ----- prise en compte des ecoulements de de decollement --------
  229.         if(nbracl.ne.0) then
  230.             do i=1,NBRACL
  231.                 icr=lnumcrl(i)  
  232.                 if(ninc.eq.1) then
  233.                     idebut=2*NDIMG
  234.                 else
  235. c                   pour les idebut des autres inclusions il faut
  236. c                   pouvoir calculer idebut
  237.                     print*,'Il manque un traitement  ninc>1 ds denc3d'
  238.                     err1=1
  239.                     return
  240.                 end if
  241. c               contribution de ce critere sur le denominateur
  242.                 do idir=1,NDIMG
  243.                     denom=denom+
  244.      #              RgR(idebut+idir)*DGFLDS(icr,idir)*FL(icr)
  245.                 end do
  246.             end do
  247.         end if
  248.  
  249. c       ---numerateur --------------------------------------------------
  250.  
  251.         numerateur=-(rest+resc+resl)
  252.  
  253. c       ---multiplicateur plastique generalise--------------------------
  254.         if(denom.ne.0.) then
  255.           lambda=numerateur/denom
  256.         else
  257.           if(numerateur.lt.0.) then
  258.             print*,'Dans denc3d le denominateur de lambda est nul'
  259.             print*,'NBRACT traction diffuse',NBRACT
  260.             do i=1,NTMAX
  261.                 icr=lnumcrt(i)
  262.                 write(*,'(a4,i3,1x,a4,i3,1x,a3,e10.3,L2)')
  263.      #           'num:',i,'dir:',icr,'Ft:',ft(icr),actift(icr)
  264.             end do
  265.             print*,'NBRACT cisaillement',NBRACC
  266.             do i=1,NCMAX
  267.                 icr=lnumcrc(i)
  268.                 write(*,'(a4,i3,1x,a4,i3,1x,a3,e10.3,L2)')
  269.      #           'num:',i,'dir:',icr,'Fc:',fc(icr),actifc(icr)
  270.             end do
  271.             print*,'NBRACT traction decollement localise',NBRACL
  272.             do i=1,NLMAX
  273.                 icr=lnumcrl(i)
  274.                 write(*,'(a4,i3,1x,a4,i3,1x,a3,e10.3,L2)')
  275.      #           'num:',i,'dir:',icr,'Fl:',fl(icr),actifl(icr)
  276.             end do
  277. c           arret
  278.             err1=1
  279.             return
  280.           else
  281.             print*,'Mise a zero du multiplicateur dans denc3d'
  282.             lambda=0.d0
  283.           end if
  284.         end if
  285.  
  286. c----------- deformations plastiques -----------------------------------
  287.  
  288. c       prise en compte de la relaxation du residu                
  289.         lambda=lambda*reduc_resg
  290.  
  291. c       initialisation des ecoulements
  292.         do idir=1,NDIMG
  293.             DEPST1(idir)=0.d0
  294.             DEPSC1(idir)=0.d0
  295.             DEPSL1(idir)=0.d0
  296.         end do
  297.  
  298. c       ---- prise en compte des ecoulements de traction diffuse -------
  299.         if(nbract.ne.0) then
  300.             do i=1,NBRACT
  301.                 icr=lnumcrt(i) 
  302. c               contribution de ce critere sur le denominateur
  303.                 do idir=1,NDIMG
  304.                     DEPST1(idir)=DEPST1(idir)+
  305.      #              lambda*DGFTDS(icr,idir)*FT(icr)
  306.                 end do 
  307.             end do
  308.         end if
  309.  
  310. c       ---- prise en compte des ecoulements de cisaillement diffuse ---
  311.         if(nbracc.ne.0) then 
  312.             do i=1,nbracc
  313.                 icr=lnumcrc(i)
  314. c               contribution de ce critere sur le denominateur
  315.                 do idir=1,NDIMG
  316.                     DEPSC1(idir)=DEPSC1(idir)+
  317.      #              lambda*DGFCDS(icr,idir)*FC(icr)
  318.                 end do             
  319.             end do
  320.         end if
  321.  
  322. c       ----- prise en compte des ecoulements de decollement -----------
  323.         if(nbracl.ne.0) then
  324.             do i=1,NBRACL
  325.                 icr=lnumcrl(i)  
  326. c               contribution de ce critere sur le denominateur
  327.                 do idir=1,NDIMG
  328.                     DEPSL1(idir)=DEPSL1(idir)+
  329.      #              lambda*DGFLDS(icr,idir)*FL(icr)
  330.                 end do
  331.             end do
  332.         end if
  333.  
  334.  
  335.         if(affiche_local) then
  336.             do i=1,NDIMG
  337.                 write(*,'(A6,I2,A2,E10.3)')
  338.      #          'DEPST(',i,')=',DEPST1(i)
  339.             end do
  340.             do i=1,NDIMG
  341.                 write(*,'(A6,I2,A2,E10.3)')
  342.      #          'DEPSC(',i,')=',DEPSC1(i)
  343.             end do
  344.             do i=1,NDIMG
  345.                 write(*,'(A6,I2,A2,E10.3)')
  346.      #          'DEPSL(',i,')=',DEPSL1(i)
  347.             end do
  348.         end if
  349.  
  350.  
  351.         return
  352.         end
  353.  
  354.  
  355.  
  356.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales