Télécharger la plaquette Cast3M


fr en

Télécharger crir3d.eso

Retour à la liste

Numérotation des lignes : 

crir3d

  1. C CRIR3D    SOURCE    FD218221  24/02/07    21:15:08     11834          
  2. c***********************************************************************
  3.       subroutine crir3d(bg,pg,bs,ps,bw,pw,rteff,cwrt,kwrt,reft,
  4.      #rceff,kwrc,delta,beta,rteff_app,rceff_app,actif,factif,nc,err1,
  5.      #na,sig13,cg,hplg,cs,hpls,cw,hplw,epspg6p,epsps6p,epspt6p,
  6.      #epspc6p,alphadp3,betadp3,refermtl,taudp,taudpmax,complet,CPTCR)
  7.  
  8. c     criteres de plasticite pour fldo3d,
  9. c     cf. coupl3d pour l ordre des criteres :
  10.  
  11. c------variables -----------N° inconnue ---------N° critere ------------
  12. c       dsigm (rankine)          25-27       criteres: 10-12
  13. c       dscrg (rag)              28-30       criteres: 1-3
  14. c       dscrs (def)              31-33       criteres: 4-6
  15. c       dscrw (capillaire)       34-36       criteres: 7-9
  16. c       dsDP  (Drucker Prager)   37          critere : 13 
  17. c-----------------------------------------------------------------------
  18.  
  19.       implicit real*8 (a-h,o-z)
  20.       implicit integer (i-n)     
  21.  
  22. c     variables externes
  23.       real*8 bg,pg,bs,ps,bw,pw,rt,cwrt,kwrt,reft,rc,kwrc
  24.       real*8 delta,beta,rteff_app,rceff_app
  25.       integer nc
  26.       logical actif(nc),complet(nc)
  27.       logical CPTCR(nc,nc)
  28.       real*8 factif(nc)
  29.       integer err1,na
  30.       real*8 sig13(3)
  31.       real*8 epspg6p(6),epsps6p(6),epspt6p(6),epspc6p(6)
  32.       real*8 alphadp3(3),betadp3(3)
  33.  
  34. c     variables locales      
  35.       real*8 bgpg,bwpw,bsps,sqrt3,sig_th,rt_dp
  36.       integer i
  37.       real*8 sigma,sigeq,siglim,taur,sigms,sigmd3(3),sigc3(3),taueq
  38.       logical confdp,refermtl
  39.       real*8 coeffdp,taudpmax,limit,fdp,xdp
  40.  
  41. c     gestion de la precision
  42.       real*8 precision1,precision2
  43.       parameter(precision1=1.0d-5,precision2=1.0d-6)
  44.  
  45. c     priorisation des criteres en pression
  46.       logical crpr  
  47. c     couplage refermeture avec DP
  48.       logical casr(3),cast(3) 
  49. c     nombre de rankine localises actifs
  50.       integer nr      
  51.  
  52.  
  53. c***********************************************************************
  54. c     calcul des contraintes intraporeuses
  55.       bgpg=bg*pg
  56.       bwpw=bw*pw
  57.       bsps=bs*ps
  58.  
  59. c***********************************************************************
  60. c     calcul des resistances effectives totales (apparentes avant endo)
  61.       rteff_app=rteff-bwpw*kwrt
  62.       rceff_app=rceff-bwpw*kwrc     
  63.  
  64. c***********************************************************************      
  65. c     test de la coherence des pressions intra poreuses
  66.       if((bgpg.lt.0.d0).or.(bwpw.gt.precision2).or.(bsps.lt.0.d0)) then
  67.          print*,'cas de pression intraporeuse imprevue dans crir3d'
  68.          print*,'bgpg:',bgpg,'bwpw:',bwpw,'bsps:',bsps
  69. c         err1=1
  70. c         na=0
  71. c         return
  72.       end if
  73.  
  74. c***********************************************************************      
  75. c     verif debut de l ecrouissage de compression pour eviter les
  76. c     interactions avec la traction directe
  77.       sqrt3=sqrt(3.d0) 
  78.       if(delta.ge.sqrt3) then
  79.          print*,'delta > sqrt(3) impossible dans fldo3d et crir3d'
  80.          err1=1
  81.          na=0
  82.          return
  83.       else 
  84. c        valeur min du criter de dp pour eviter les interactions avec rankine         
  85.          sig_th=rteff_app*((sqrt3+delta)/(sqrt3-delta))
  86. c        non interaction avec rteff_app triaxial
  87.          sig_th=max(sig_th,(rteff_app*3.d0*delta)/(sqrt3-delta))
  88. c        marge pour eviter la confusion numerique         
  89.          sig_th=sig_th*(1.d0+precision2)
  90. c        test de compatibilites des criteres traction compression
  91.          if(sig_th.gt.rceff_app) then
  92.              print*,'Dans crir3d'
  93.              print*,'DP atteignable par traction triaxiale ds crir3d'
  94.              print*,'augmenter Rc ou diminuer ept  (ou ept=rteff/E)'
  95.              print*,'rt:',rteff_app
  96.              print*,'rc:',rceff_app
  97.              print*,'sig_th minimun necessaire pour DP:',sig_th
  98. c             rt max pour eviter l interaction             
  99.               rt_dp=rceff_app*min((((sqrt3+delta)/(sqrt3-delta))**(-1)),
  100.      #        (((3.d0*delta)/(sqrt3-delta))**(-1)))
  101.      #        *(1.d0-2.d0*precision2)    
  102. c              print*,'On adopte rt_dp<rteff_app ds crir3d',rt_dp
  103.              err1=1
  104.              return
  105.          else
  106.               rt_dp=rteff_app
  107.          end if         
  108.       end if    
  109.  
  110.  
  111. c***********************************************************************      
  112. c     test des criteres 
  113. c***********************************************************************
  114.  
  115. c     initialisation des indicateurs par defaut    
  116.       do i=1,nc
  117.         actif(i)=.false.
  118.         complet(i)=.true.
  119.         factif(i)=0.d0
  120.       end do
  121. c     indicateur de refermeture de fissure de traction localisee      
  122.       refermtl=.false.
  123. c     indicateur d activite d au moins un critere de pression
  124.       crpr=.false.      
  125. c     boucle sur les criteres
  126.       na=0 
  127.  
  128.       do i=1,nc
  129. c        print*,'crir3d, test du critere :',i 
  130.  
  131.         if(i.le.3) then
  132.  
  133. c----------critere pression de gel sur scrg-----------------------------
  134.  
  135.             j=i
  136. c           verifier CPTCR avec couplage3d            
  137.             sigma=sig13(j)-bwpw-bgpg-bsps
  138. c           verifier CPTCR avec couplage3d             
  139.             if(sigma.le.0.d0) then  
  140.                  if(epspt6p(j).gt.precision2) then
  141. c                   refermeture possible  sigma est limite par reft
  142.                     if(sigma.lt.-reft) then            
  143.                          sigeq=cg*pg-reft
  144. c                        on indique au solveur que le critere est tronqué de la contrainte totale                        
  145. c                         complet(i)=.false.
  146.                          complet(i)=.true.
  147.                     else
  148.                          sigeq=cg*pg+sigma
  149.                          complet(i)=.true.
  150.                     end if
  151.                  else
  152. c                   la fissure est deja refermee
  153.                     sigeq=cg*pg+sigma
  154.                     complet(i)=.true.
  155.                  end if                            
  156.             else
  157.                  if(sigma.lt.rteff_app) then
  158.                      sigeq=cg*pg+sigma
  159.                      complet(i)=.true.
  160.                  else
  161.                      sigeq=cg*pg+rteff_app
  162. c                     complet(i)=.false.
  163.                          complet(i)=.true.
  164.                  end if
  165.             end if
  166. c           la resistance a la traction dépend de la capillarite
  167. c           et de l ecrouissage a considerer dans coupl3d si ce critere actif            
  168.             siglim=rteff_app+hplg*max(epspg6p(j),0.d0) 
  169. c           precision sur le critere            
  170.             fglim=siglim*precision1            
  171.             factif(i)=sigeq-siglim
  172.             if(factif(i).gt.fglim) then
  173. c               le critere est actif
  174.                 actif(i)=.true.
  175.                 do k=1,nc
  176.                      if((.not.CPTCR(i,k)).and.actif(k)) then
  177. c                        incompatibilite detectee, on desactive le critere                 
  178.                          actif(i)=.false.
  179.                      end if
  180.                 end do
  181.                 if (actif(i)) then
  182.                      crpr=.true.                
  183. c                    actualisation nbr de criteres actifs                
  184.                      na=na+1
  185.                 end if
  186.             end if 
  187. C        if((pg.gt.0.).and.(.not.(actif(i)))) then
  188. C        write (*,'(a16,e10.3,1x,a3,e10.3,1x,a6,i2,a2,e10.3,1x,a7,e10.3)')          
  189. C      #     'Dans crir3d: pg=',pg,'cg=',cg,
  190. C      #     'sigma(',i,')=',sigma,'siglim=',siglim
  191. C        write (*,'(a5,e10.3,1x,a8,e10.3,1x,a10,e10.3)')
  192. C      #    'hplg=',hplg,'epspg6p=',epspg6p(j),'rteff_app=',rteff_app
  193. C        write (*,'(a2,e10.3,a2,l2)')
  194. C      #    'f=',factif(i),'A=',actif(i)  
  195. C        read*
  196. C       end if             
  197.  
  198.         else if ((i.ge.4).and.(i.le.6)) then
  199.  
  200. c----------critere pression de def sur scrs-----------------------------        
  201.  
  202.             j=i-3
  203. c           verifier CPTCR avec couplage3d            
  204.             sigma=sig13(j)-bwpw-bgpg-bsps
  205. c           verifier CPTCR avec couplage3d
  206.             if(sigma.lt.0.d0) then  
  207.                  if(epspt6p(j).gt.precision2) then
  208. c                   refermeture possible  sigma est limite par reft
  209.                     if(sigma.lt.-reft) then            
  210.                          sigeq=cs*ps-reft
  211. c                        on indique au solveur que le critere est tronqué                         
  212. c                         complet(i)=.false.
  213.                          complet(i)=.true.
  214.                     else
  215.                          sigeq=cs*ps+sigma
  216.                          complet(i)=.true.
  217.                     end if
  218.                  else
  219. c                   la fissure est deja refermee
  220.                     sigeq=cs*ps+sigma
  221.                  end if
  222.             else
  223.                  if(sigma.lt.rteff_app) then
  224.                      sigeq=cs*ps+sigma
  225.                      complet(i)=.true.
  226.                  else
  227.                      sigeq=cs*ps+rteff_app
  228. c                     complet(i)=.false.
  229.                          complet(i)=.true.
  230.                  end if
  231.             end if
  232. c           la resistance a la traction dépend de la capillarite
  233. c           et de l ecrouissage a considerer dans coupl3d si critere actif            
  234.             siglim=rteff_app+hpls*max(epsps6p(j),0.d0) 
  235. c           precision sur le critere            
  236.             fglim=siglim*precision1            
  237.             factif(i)=sigeq-siglim
  238.             if(factif(i).gt.fglim) then
  239. c               le critere est actif
  240.                 actif(i)=.true.
  241.                 do k=1,nc
  242.                      if((.not.CPTCR(i,k)).and.actif(k)) then
  243. c                        incompatibilite detectee, on desactive le critere                 
  244.                          actif(i)=.false.
  245.                      end if
  246.                 end do
  247.                 if (actif(i)) then
  248.                      crpr=.true.                
  249. c                    actualisation nbr de criteres actifs                
  250.                      na=na+1
  251.                 end if
  252.             end if
  253.  
  254.        else if ((i.ge.7).and.(i.le.9)) then
  255.  
  256. c----------critere pression d eau sur scrw------------------------------        
  257.  
  258.             j=i-6
  259. c           verifier CPTCR avec couplage3d            
  260.             sigma=sig13(j)-bwpw-bgpg-bsps
  261. c           verifier CPTCR avec couplage3d 
  262.             if(sigma.lt.0.d0) then            
  263.                  if(epspt6p(j).gt.precision1) then
  264. c                   refermeture possible  sigma est limite par reft
  265.                     if(sigma.lt.-reft) then            
  266.                          sigeq=-cw*bwpw-reft
  267. c                        on indique au solveur que le critere est tronqué                         
  268.                          complet(i)=.false.
  269.                     else
  270.                          sigeq=-cw*bwpw+sigma
  271.                     end if
  272.                  else
  273.                      sigeq=-cw*bwpw+sigma
  274.                  end if
  275.             else
  276.                  sigeq=-cw*bwpw
  277.                  complet(i)=.false.
  278.             end if
  279. c           la resistance a la traction dépend de la capillarite           
  280.             siglim=rteff_app
  281. c           precision sur le critere            
  282. c           fglim=siglim*precision1            
  283.             factif(i)=sigeq-siglim 
  284.             factif(i)=max(factif(i),0.d0)
  285. c           pas d ecoulement pour ce critere, car seulement endo  
  286. c           hydrique contraintes effectives            
  287.             actif(i)=.false.    
  288.  
  289.  
  290.        else if ((i.ge.10).and.(i.le.12)) then
  291.  
  292. c----------critere traction localisee ----------------------------------        
  293.  
  294.             j=i-9
  295. c           verifier CPTCR avec couplage3d            
  296.             sigma=sig13(j)-bwpw-bgpg-bsps
  297. c           verifier CPTCR avec couplage3d             
  298.             sigeq=sigma
  299. c           la resistance a la traction dépend de la capillarite
  300. c           et de l ecrouissage a considerer dans coupl3d si critere actif
  301.             if(sigeq.gt.0.) then            
  302.                  siglim=rteff_app
  303.             else
  304.                 if(epspt6p(j).ge.precision2) then
  305. c                   refermeture possible 
  306.                     siglim=reft
  307.                 else
  308.                     siglim=rteff_app
  309.                 end if
  310.             end if 
  311. C            print*,'direction', j,'sigma',sigma,'siglim',siglim               
  312. c           precision sur le critere            
  313.             fglim=siglim*precision1 
  314.             if(sigeq.gt.0.d0) then
  315. c                le critere a annuler est positif            
  316.                  factif(i)=sigeq-siglim
  317.                  if(factif(i).gt.fglim) then
  318. c                    le critere est actif
  319.                      actif(i)=.true.
  320.                      do k=1,nc
  321.                          if((.not.CPTCR(i,k)).and.actif(k)) then
  322. c                             incompatibilite detectee, on desactive le critere  
  323. C                               print*,'Incompatibilite',i ,'avec', k  
  324. C                               read*             
  325.                               actif(i)=.false.
  326.                          end if
  327.                      end do
  328.                      if (actif(i)) then              
  329. c                         actualisation nbr de criteres actifs                
  330.                           na=na+1
  331.                      end if
  332.                 end if               
  333.             else
  334.                 if(epspt6p(j).ge.precision2) then
  335. c                 le critere a annuler est negatif pour ne pas modifier coupl3d            
  336.                   factif(i)=sigeq+siglim
  337.                   if(factif(i).lt.-fglim) then                 
  338.                      actif(i)=.true.
  339. c                    le critere est actif
  340.                      actif(i)=.true.
  341.                      do k=1,nc
  342.                          if(.not.CPTCR(i,k).and.actif(k)) then
  343. c                             incompatibilite detectee, on desactive le critere                 
  344.                               actif(i)=.false.
  345.                          end if
  346.                      end do
  347.                      if (actif(i)) then              
  348. c                         actualisation nbr de criteres actifs                
  349.                           na=na+1
  350. c                         indicateur de refermeture active                          
  351.                           refermtl=.true.                         
  352.                      end if                     
  353.                   end if
  354.                 end if
  355.             end if
  356.  
  357.  
  358.        else if (i.eq.13) then
  359.  
  360. c----------critere de cisaillement -------------------------------------
  361.  
  362. c       ce critere n est traitee qu'apres les refermetures de fissures
  363. c       localisees de traction 
  364. c       contrainte hydrostatique (totale)      
  365.         sigms=0.d0         
  366.         do j=1,3
  367.              sigc3(j)=min(sig13(j)-bwpw-bgpg-bsps,rt_dp)
  368. c            modif en fonction des couplages avec les refermetures             
  369.              if(actif(9+j)) then
  370.                  if(factif(9+j).lt.0.) then
  371. c                    on referme la fissure elle est donc > (-reft) et 
  372. c                    ne va pas varier car doit verifier l autre critere                    
  373.                      sigc3(j)=max(sigc3(j),-reft)
  374.                      casr(j)=.true.
  375.                      cast(j)=.false.
  376.                  else
  377. c                    on ouvre la fissure elle est donc <rteff_app et 
  378. c                    ne va pas varier car doit verifier l autres critere
  379.                      sigc3(j)=min(sigc3(j),rteff_app) 
  380.                      cast(j)=.true.
  381.                      casr(j)=.false.
  382.                  end if
  383.              else  
  384.                  casr(j)=.false.
  385.                  cast(j)=.false.   
  386.              end if                 
  387.              sigms=sigms+sigc3(j)
  388.         end do
  389.         sigms=sigms/3.d0
  390. c        if(sigms.gt.0.) then
  391. c            print*,'Dans crir3d sigms=', sigms
  392. c        end if            
  393. c       deviateur et contrainte de Von Mises         
  394.         taueq=0.d0
  395.         do j=1,3
  396.             sigmd3(j)=sigc3(j)-sigms             
  397.             taueq=taueq+sigmd3(j)**2
  398.         end do
  399.         taueq=sqrt(taueq/2.d0)
  400. c       contrainte de Drucker Prager 
  401.         if(sigms.lt.0.d0) then
  402. c           effet favorable du confinement
  403.             taudp=taueq+delta*sigms
  404.             confdp=.true.
  405.         else
  406. c           effet du deconfinement non pris en compte
  407.             taudp=taueq
  408.             confdp=.false.
  409. c           print*,'ds crir3d pas de confinement sur dp'
  410.         end if
  411. c       Critere de Drucker Prager limite 
  412.         if (confdp) then
  413.              coeffdp=(1.d0/sqrt(3.d0)-delta/3.d0)
  414.         else
  415.              coeffdp=1.d0/sqrt(3.d0)
  416.         end if
  417.         taudpmax=coeffdp*rceff_app
  418.         limit=taudpmax*precision1
  419.         fdp=taudp-taudpmax
  420. c       calcul des coeffs pour l ecoulement de Drucker Prager        
  421.         if((taueq.gt.0.).and.(fdp.ge.limit)) then
  422.              do j=1,3
  423.                  xdp=0.5d0*sigmd3(j)/taueq
  424.                  if (confdp) then
  425.                      alphadp3(j)=xdp+delta/3.d0
  426.                  else
  427.                      alphadp3(j)=xdp
  428.                  end if
  429.                  if(casr(j).or.cast(j)) then
  430. c                    refermeture ou ouverture simultane direction j                 
  431.                      betadp3(j)=0.d0
  432.                  else
  433. c                    cas normal                 
  434.                      betadp3(j)=xdp+beta/3.d0
  435.                  end if                     
  436. c                 print*,'crir3d:',j,alphadp3(j),betadp3(j)
  437.              end do
  438.         else
  439.              do j=1,3
  440.                 alphadp3(j)=0.d0
  441.                 betadp3(j)=0.d0
  442.              end do
  443.         end if
  444. c       test du critere  
  445.         factif(i)=fdp      
  446.         if(fdp.ge.limit) then
  447. c            le critere est actif
  448.              actif(i)=.true.
  449.              do k=1,nc
  450.                  if((.not.CPTCR(i,k)).and.actif(k)) then
  451. c                     incompatibilite detectee, on desactive le critere                 
  452.                       actif(i)=.false.
  453.                  end if
  454.              end do
  455.              if (actif(i)) then              
  456. c                actualisation nbr de criteres actifs                
  457.                  na=na+1
  458. c                on compte les criteres de traction / refermeture actif 
  459.                  nr=0               
  460.                  do k=1,3
  461.                      if(actif(9+k)) then
  462.                        nr=nr+1
  463.                        if(factif(9+k).lt.0.) then
  464. c                           on a une refermeture dans la direction j et 
  465. c                           un ecoulement DP dans cette meme direction
  466. c                           on garde le couplage 
  467.                             continue
  468.                        else
  469. c                           on a un ecoulement dp et une traction dans
  470. c                           la meme direction, on ne garde que DP
  471. c                           on repousse la verification de traction a
  472. c                           une autre sous iteration
  473.                             actif(9+k)=.false.
  474.                             nr=nr-1
  475.                             na=na-1                            
  476.                        end if
  477.                      end if                       
  478.                  end do
  479. c                on comptabilise DP dans les actifs
  480.                  if(nr.eq.3) then
  481. c                    on annule le couplage avec DP car on a trois 
  482. c                    couplages avec les refermetures, on fait dabord les
  483. c                    3 refermetures
  484.                      actif(i)=.false.
  485.                      na=na-1
  486.                      if(na.eq.0) then
  487.                          print*,'Desactivation DP ds crir3d'
  488.                          print*,'Pb ds crir3d na=0 et refermtl=T !'
  489.                          do k=1,13
  490.                               print*,factif(k)
  491.                          end do
  492.                          err1=1
  493.                          return
  494.                     end if                       
  495.                  end if
  496.              end if                
  497.         end if
  498. c        do j=1,3
  499. c          if((sig13(j).lt. -rceff_app).and.(.not.actif(13)).and.
  500. c     #    (.not.refermtl)) then
  501. c               print*,'DP faux et sigeff',j,'=',sig13(j)
  502. c               do k=1,3
  503. c                 print*,'sig',k,'=',sig13(k)
  504. c               end do
  505. c               print*,'rceff_app=',rceff_app
  506. c               err1=1
  507. c               return
  508. c          end if
  509. c        end do
  510. c       fin de test des criteres          
  511.        end if 
  512. c      fin de boucle sur les criteres       
  513.       end do  
  514. C       print*,'Dans crir3d'
  515. C       do i=1,13
  516. C          write(*,'(A7,I3,L1,E10.3)') 'critere', i, actif(i),factif(i)
  517. C       end do 
  518. C       read*         
  519.       return       
  520.       end
  521.  
  522.  
  523.  
  524.  
  525.  
  526.  

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