Télécharger la plaquette Cast3M


fr en

Télécharger endorf.eso

Retour à la liste

Numérotation des lignes : 

endorf

  1. C ENDORF    SOURCE    PV090527  23/01/27    21:15:30     11574          
  2.            subroutine endorf(sigref,eplrf,syr,sur,epsrf,epu,hplr,
  3.      #     yor,wpr,lcr,damrt0,damrtf,damrt1,damrt2,damrc0,damrcf,refr1,
  4.      #     epsrt,epsrc,plr_iso,diff0,difff,Er,Dr,Hi,grad_sigr,
  5.      #     tauie,tyr,dami0,damif,epsmf,eplr0,sigr0,METHODNL,epsm0,epsr0,
  6.      #     depsa,Et0,Etf,Hs0,Hsf,sigrf,spref,epart,eparv,source_nl,
  7.      #     coeff_nlr,taui0,tauif,coeff_nli,eprkf,eprmf,
  8.      #     hypo,istep,err1,garder,endo_interface_explicit)  
  9.  
  10. c          tables de dimension fixe pour resolution des sytemes lineaires 
  11.            implicit real*8 (a-h,o-z)
  12.            implicit integer (i-n)     
  13.  
  14. c          variables externes     
  15.            real*8 sigref,eplrf,syr,sur,epu,hplr,yor,wpr,lcr,coeff_nlr
  16.            real*8 damrt1,damrt2,refr1,epsrt,epsrf,diff0,difff
  17.            real*8 Er,Dr,Hi,lcgr,grad_sigr,dami0,damif,taui,tyr,damrtf
  18.            real*8 eplr0,epsmf,damrcf,sigr0,depsa,epsm0,epsr0,eparv
  19.            real*8 Et0,Etf,Hs0,Hsf,damrt0,damrc0,sigrf,source_nl,spref
  20.            real*8 eprkf,eprmf,depspl,hypo,epsrc
  21.            logical plr_iso,METHODNL,garder
  22. c          variables locales           
  23.            real*8 damrt3,damrt4,epr,dmax,xref,epse,epart
  24.            integer irest0
  25.            real*8 precision1,epse_max,damrc3,damrt5,Hid,drf,dr0,Ert
  26.            real*8 tauie,taui0,tauif,coeff_nli,eppr_et,coeffd
  27. c          coeff de reduction d ecrouissage pour la phase 2           
  28.            real*8 deplrdt2,eprk
  29.            integer err1,istep
  30.            parameter (precision1=1.0d-5,dmax=1.d0-precision1)
  31.            real*8 DPI,xx,depsrr
  32.            PARAMETER (DPI=3.141592653589793D0)
  33. c          gestion de la methode de calcul du module du renfort 
  34.            logical tangent,endo_interface_explicit
  35.            parameter (tangent=.false.)
  36. c          coeff correcteur de la longeur pour prendre en compte
  37. c          l elargissement de la zone plasqtique endommagee
  38. c          due a l intrpoation des contraintes non locales lorsque
  39. c          l endo est local
  40.            real*8 coeff_lc
  41.            parameter (coeff_lc=1.0d0)     
  42.  
  43. c*************** test des hypotheses de calculs prevues  ***************
  44. c          test calcul en module tangent          
  45.            if(tangent) then
  46.                 print*,'Module tangent non disponible'
  47.                 print*,'dans endorf'
  48.                 err1=1
  49.                 return
  50.            end if                
  51.  
  52. c          test option de plasticite du renfort           
  53.            if  (plr_iso) then
  54.              print*,'Plasticite isotrope des renforts'
  55.              print*,'incompatible avec endo des des renfort'
  56.              print*,'dans endorf'
  57.              err1=1
  58.              return
  59.            end if 
  60.  
  61. c**************** initialisation *************************************** 
  62.            err1=0            
  63.  
  64. c**************** endommagement du renfort *****************************           
  65.  
  66.            if (sur.ge.syr) then         
  67. c           endommagement de palier d ecrouissage 
  68.             if(garder) then
  69. c              la solution actuelle est issue d'hypotheses coherentes            
  70.                if(sigref.gt.sur) then
  71. c                   pas d endo de compression                    
  72.                     damrc3=0.d0               
  73. c                   cas de la traction               
  74.                     damrt3=1.d0-sur/sigref   
  75.                else 
  76.                     if(sigref.lt.-sur) then
  77. c                       pas d endo de traction                        
  78.                         damrt3=0.d0                    
  79. c                       cas de la compression                    
  80.                         damrc3=1.d0+sur/sigref
  81.                     else
  82. c                       pas d endo                    
  83.                         damrt3=0.d0
  84.                         damrc3=0.d0
  85.                     end if
  86.                end if
  87.                damrt1=max(damrt1,damrt3)
  88.                damrcf=max(damrc0,damrc3)
  89.             else
  90. c              une hypothese incoherente detecteee            
  91. c              on reprend les endos initiaux           
  92.                damrcf=damrc0
  93.             end if
  94.            else
  95.                 print*,'Dans endorf il faur sur > syr'
  96.                 err1=1
  97.                 return
  98.            end if
  99.  
  100. c          increment de deformation plastique           
  101.            depspl=eplrf-eplr0 
  102.  
  103.            if((epu.ge.0.).and.(wpr.gt.0.)) then
  104. c           les donnees pour la rupture sont coherentes  
  105. c           deformation caracteristique de rupture en traction
  106. c           en fonction de la longueur caracteristique d integration lcr
  107. c           - decroissance Gaussienne
  108. c           eprk=epu+(wpr/(lcr*sur))*sqrt(2.d0/DPi)*(yor/(yor+hplr))
  109. c           - decroissance exponentielle
  110.            eprk=epu+(wpr/(coeff_lc*lcr*sur))*(yor/(yor+hplr))
  111. c           - decroissance parabolique
  112. c            eprk=epu+(3.d0*wpr/(2.d0*lcr*sur))*(yor/(yor+hplr))
  113.             if (garder) then 
  114. c              les hypotheses d ecoulement etaient coherentes            
  115.                if(depspl.gt.0.d0) then
  116. c                    cumul des deformation positive pour endo tract               
  117.                      epsrt=epsrt+depspl
  118. c                    remise a zero de la def de refermeture                     
  119.                      epsrc=0.d0
  120.                else
  121. c                    actualisation de la refermeture               
  122.                      epsrc=epsrc+depspl
  123.                end if
  124. c              evolution de l endo de traction
  125.                if(epsrt.gt.epu) then
  126. c                    - endo Gaussien               
  127. c                     xx=((epsrt-epu)/(eprk-epu))**2
  128. c                     damrt5=1.d0-exp(-0.5d0*xx)
  129. c                    - endo exponentielle
  130.                      xx=((epsrt-epu)/(eprk-epu))
  131.                      damrt5=1.d0-exp(-xx)
  132. c                    - endo parabolique
  133. c                     if(epsrt.lt.eprk) then
  134. c                         damrt5=((epsrt-epu)/(eprk-epu))**2
  135. c                     else
  136. c                         damrt5=dmax
  137. c                     end if
  138.                       damrt5=min(damrt5,dmax)
  139.                else
  140.                      damrt5=damrt2
  141.                end if
  142. c              actualisation endo de traction localise               
  143.                damrt2=max(damrt5,damrt2)
  144.             else
  145. c              les hypotheses d ecoulement etaient incoherentes            
  146. c              on reprend le cas converge precedent
  147.                continue
  148.             end if                              
  149.            else
  150.                print*,'Pas de donnees pour la rupture du renfort'
  151.                print*,'dans endorf'
  152.                err1=1
  153.                return
  154.            end if
  155.  
  156. c          endommagement resultant en traction        
  157.            damrt4=1.d0-(1.d0-damrt1)*(1.d0-damrt2)
  158.  
  159. c          condition de croissance pour l endo de traction           
  160.            damrtf=max(damrt0,damrt4)
  161.  
  162. c          endommagement maximum autorisé           
  163.            damrtf=min(damrtf,dmax)
  164.  
  165. c**********  coeff de refermeture de fissure ***************************
  166.  
  167. c       increment total de def plastique necessaire pour refermer
  168. c       - cas de l endo Gaussien                
  169. c        depsrr=sqrt(0.2D1)*sqrt(-log((1.d0-damrt2))*(eprk-epu)**2)
  170. c       - cas de l endo exponentiel  
  171.         depsrr=(epu-eprk)*log(1.d0-damrt2)
  172. c       - cas de l endo parabolique   
  173. c        depsrr=sqrt(damrt2*((eprk-epu)**2))
  174. c       test somme des increments neg - somme des incr positifs               
  175.         if(epsrc.gt.-depsrr) then               
  176. c            le cumul est en faveur de l ouverture  
  177. c            - endo Gaussien          
  178. c             xx=((depsrr+epsrc)/(eprk-epu))**2
  179. c             refr1=exp(-0.5d0*xx)
  180. c            - endo exponentiel 
  181.              xx=((depsrr+epsrc)/(eprk-epu))
  182.              refr1=exp(-xx)
  183. c            - endo parabolique
  184. c             refr1=1.d0-((depsrr+epsrc)/(eprk-epu))**2
  185. c            bornes de la fonction de refermeture
  186.              refr1=min(1.d0,max(precision1,refr1))             
  187.         else
  188. c            le cumul est une refermeture               
  189.              refr1=1.d0
  190.         end if             
  191.  
  192. c**************** endommagement de l  interface ************************ 
  193.  
  194.         if((METHODNL.and.(istep.gt.1).and.(garder).and.
  195.      #     (.not.endo_interface_explicit)).or.
  196.      #     ((METHODNL.and.(istep.eq.1).and.
  197.      #     (endo_interface_explicit))))then
  198. c          les hypotheses etaient coherentes        
  199. c          approximation de la contrainte de cisaillement actuelle
  200.            taui=-grad_sigr*Dr/4.d0
  201. c          test de plasticification de l interface
  202.            if((tyr.eq.0.).and.(taui.ne.0.)) then
  203.                  print*,'Dans endorf'
  204.                  print*,'Resistance de l interface non fournie'
  205.                  print*,'pour endorf appele par fldo3d'
  206.                  print*,'fournir TYRi non nul pour continuer...'
  207.                  err1=1
  208.                  return
  209.            else  
  210. c               contrainte effective d interface           
  211.                 tauie=taui/coeff_nli
  212.                 if(abs(tauie).gt.tyr) then
  213.                     damif=1.d0-tyr/abs(tauie)
  214. c                   condition de croissance de lendo de l interface                 
  215.                     damif=max(damif,dami0)
  216.                     damif=min(damif,dmax)
  217.                 else
  218.                     damif=dami0
  219.                 end if                 
  220. c               raideur de l interface endommagee
  221.                 Hid=Hi*(1.d0-damif) 
  222.            end if
  223.         else if (METHODNL)then             
  224.             if(((.not.endo_interface_explicit).and. 
  225.      #            ((istep.eq.1).or.(.not.garder))).or. 
  226.      #            (endo_interface_explicit.and.(istep.ne.1)))then
  227. c                 on reprend la valeur convergee
  228.                   damif=dami0            
  229.                   Hid=Hi*(1.d0-dami0)
  230.                   tauie=taui0/coeff_nli
  231.             else
  232.                   print*,'endorf: Cas non prevu / interface'
  233.                   print*,METHODNL,endo_interface_explicit
  234.                   print*,istep,garder
  235.                   err1=1
  236.                   return
  237.             end if                        
  238.         else if(.not.METHODNL) then
  239. c           calcul local            
  240.             Hid=Hi
  241.             dami0=0.d0
  242.             damif=dami0 
  243.             taui=0.d0                
  244.             tauie=taui 
  245.         else                
  246. c           on initialise a la valeur convergee precedente 
  247.             tauie=taui0/coeff_nli           
  248.             damif=dami0            
  249.             Hid=Hi*(1.d0-damif)
  250.         end if
  251. c       module de l interface        
  252.         if(METHODNL) then 
  253. c           module secant        
  254.             Hsf=Hid
  255.         else 
  256. c           module initial        
  257.             Hsf=Hi
  258.         end if
  259. c       contrainte dans l interface apres endommagement
  260.         coeff_nli=(1.d0-damif) 
  261.         tauif=tauie*coeff_nli
  262.  
  263.  
  264. c********** coeff de diffusion pour l adherence ************************
  265.  
  266. c       calcul en module secant
  267.         if(.not.garder) then
  268. c           on reprend la contrainte convergee pour calculer les coeffs
  269. c           non lineaires necessaire au calcul de diffusion
  270.             sigref=sigr0
  271.         end if            
  272.         if (sigref.ge.0.d0) then 
  273. c            avec l endo de traction
  274.              coeff_nlr=(1.0d0-damrtf)                         
  275.         else
  276. c            avec l endo de compression et la refermeture
  277.              coeff_nlr=(1.0d0-damrcf)*refr1
  278.         end if
  279.  
  280.  
  281. c       module endommagee actuel pour le renfort 
  282.         Ert=Er*coeff_nlr
  283.  
  284. c       coeff de diffusion
  285.         if(METHODNL) then
  286. c            .and.(istep.ne.1)) then
  287. c            actualisation du coeff pour le pas suivant
  288.              difff=(Ert*Dr)/(4.d0*Hid)    
  289. c            calcul suivant hypothese d ecoulement                    
  290.              difff=difff*(hplr/(hplr+yor*abs(hypo)))
  291. c             difff=4.0d-2             
  292. c         else if (METHODNL) then
  293. c             print*,'Cas imprevu sur lnl dans endorf'
  294. c             err1=1
  295. c             return
  296.          else if ((istep.eq.0).or.(.not.METHODNL)) then
  297. c            calcul local                 
  298.              difff=0.d0 
  299.              Etf=Et0
  300.              Hsf=Hs0 
  301.          else 
  302.              print*,'Cas imprevu sur lnl dans endorf'
  303.              err1=1
  304.              return                 
  305.          end if
  306. c        print*,'endo renfort lnl',lnl,' istep',istep                
  307.  
  308. c        module actuel du renfort       
  309.          if(METHODNL) then       
  310.             Etf=Ert
  311.          else
  312.             Etf=Er
  313.          end if          
  314.  
  315. c*********** contraintes apparentes dans les renforts ******************
  316.  
  317. c        calcul des contraintes apparentes dans les renforts        
  318.          sigrf=sigref*coeff_nlr
  319. c        test coherence de la contrainte        
  320.          if(abs(sigrf).gt.sur*(1.d0+precision1)) then
  321. c            incoherence detectee         
  322.              print*,'Contrainte anormale dans dans endorf'
  323.              print*,'|sirf|>sur',sigrf,sur
  324.              print*,'Coeff de non linearite',coeff_nlr
  325.              print*,'valeur des endommagements'
  326.              print*,'damrt1',damrt1
  327.              print*,'damrt2',damrt2
  328.              print*,'damrtf',damrtf 
  329.              print*,'damrcf',damrcf
  330.              print*,'refermeture',refr1
  331.              print*,'Def plastique actuelle',eplrf
  332.              print*,'Def totale actuelle',epsrf
  333.              if(damrt2.ge.dmax) then
  334.                  print*,'Def de refermeture actuelle',epsrt
  335.              end if
  336.              err1=1
  337.              return
  338.          end if
  339.          if(difff.lt.0.) then
  340.           print*,'diffusion negative dans endorf'
  341.           err1=1
  342.           return
  343.          end if
  344.  
  345. c*********** source non locale *****************************************
  346.  
  347.          if(METHODNL) then 
  348.              if (hypo.ne.0.) then
  349. c               hypothese avec ecoulement plastique
  350.                 eppr_et=syr/hplr
  351.                 if (hypo.eq.1.) then                
  352.                     source_nl=Hid*(epsmf-(eparv -eppr_et))                              
  353.                 else if(hypo.eq.-1.) then
  354.                     source_nl=Hid*(epsmf-(eparv+eppr_et))             
  355.                 else
  356.                     print*,'valeur anormale de hypo ds endorf'
  357.                     err1=1
  358.                     return
  359.                 end if 
  360.              else           
  361. c               hypothese tir visco elastique                
  362.                 source_nl=Hid*(epsmf-epart)            
  363.              end if          
  364. c            coeff multiplicateur des sources
  365.              coeffd=4.d0*difff/Dr    
  366.              source_nl=(epsmf-epart)*Ert
  367. c            print*,'endorf, source_nl=',source_nl
  368. c             print*,'endorf, Difff=',difff  
  369.          else
  370. c            calcul local: source inutile         
  371.              source_nl=0.d0
  372.          end if
  373.  
  374. c***********************************************************************
  375.  
  376.          return
  377.          end
  378.  
  379.  
  380.  

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