Télécharger la plaquette Cast3M


fr en

Télécharger rnfr3d.eso

Retour à la liste

Numérotation des lignes : 

rnfr3d

  1. C RNFR3D    SOURCE    FD218221  24/02/07    21:15:25     11834          
  2.       subroutine rnfr3d(istep,nbrenf,numr,epstf33,vecr,epsr0,epsrf,
  3.      #eprp0,eprpf,your,syr,sigre0,sigref,hplr,tomr,ekr,skr,ATRref,khir,
  4.      #gamr,sprec,teta1,teta2,dt,ppas,theta,eprm0,eprmf,ttaref,rhor,
  5.      #mu00,fl3d,errr,xnr,xmuthr,eprk0,eprkf,tokr,Kr,plast_seule,
  6.      #ann,xn,bn,ngf,ipzero,epspmf33,spre0,spref,epst033,depsa,epleq0,
  7.      #epleqf,epsmf,epsm0,sigr_nl,coeff_nlr,epsre0,epart,eparv,
  8.      #hypo,plr_iso,METHODNL,garder,Pb_nl,travail_fissure)
  9.  
  10. c     calcul de la deformation et de la contrainte axiale dans un renfort
  11.  
  12. c     2 variables de controle doivent être renseignees :
  13. c     - GARDER qui indique si le resultat doit être garder ou non
  14. c     - PB_NL qui indique si une hypothese a ete modifiee pour les 
  15. c     operateurs et la source non locale, si PB_NL=1 alors GARDER=.false.
  16.  
  17. c     si PB_NL=1 comme on repart de la solution precedente on peut ajuster
  18. c     les hypotheses d ecoulement via la variable HYPO, cela en fonction de la
  19. c     solution non locale actuelle. Donc même si GARDER=.false. et PB_NL=1
  20. c     HYPO peut prendre l'une des 3 valeurs (1 ecoulement plastique positif, 
  21. c     0 pas d'ecoulement plastique, -1 ecoulement plastique negatif)
  22.  
  23. c     Pour que GARDER puisse être mis a .true. il faut que toutes les hypothese
  24. c     d ecoulement soient vraies au niveau de la structure, il y a donc une detection
  25. c     des PB_NL lors de l etape non local, si au moins un PB_NL est detecte tout est
  26. c     faut et GARDER=.false., si tous les PB_NL sont à 0 alors on peut
  27. c     GARDER la solution de l ecoulement plastique. Comme l'ecoulement plastique
  28. c     est realise, on poursuit alors les etapes non locales avec une hypothese
  29. c     sans ecoulement (HYPO=0).
  30.  
  31. c     la permutation de PB_NL à 0 n est possible qu'en fin de procedure
  32. c     que si HYPO n a pas changé pour ce point
  33.  
  34. c     si PB_NL=0 au retour de l etape non locale la nouvelle contrainte non locale
  35. c     est gardee definitivement car toutes les hypotheses d ecoulement etaient
  36. c     coherentes
  37.  
  38.  
  39.  
  40.       implicit real*8 (a-h,o-z)
  41.       implicit integer (i-n)
  42.  
  43.       integer istep,nbrenf,numr,err1
  44.       logical plast_seule
  45. c     def finale, orientation      
  46.       real*8 epst033(3,3),epstf33(3,3),vecr(nbrenf,3),epspmf33(3,3)
  47. c     def totale, plastique...      
  48.       real*8 epsr0,epsrf,source_nl,eprp0,eprpf,your,syr,sigr_nl
  49.       real*8 sigre0,sigref,depsa,hplr,tomr,ekr,skr,spre0,spref
  50.       real*8 ATRref,khir,gamr,sprec,coeff_nlr
  51.       real*8 epse0,epsre0,epart,epsref,epsraf,eparv
  52.       real*8 teta1,teta2,dt,theta,epsmf,epsm0,hypo,Pb_nl
  53. c     deformation plastique cumulee pour la plsticite isotrope
  54.       real*8 epleq0,epleqf,Pb_nl0       
  55.       logical ppas
  56.       real*8 eprm0,eprmf,ttaref,rhor,mu00,eprk0,eprkf,tokr,Kr,K,epprr
  57.       logical fl3d,errr,garder
  58.       real*8 xnr,xmuthr
  59.       integer ngf,errgauss
  60.       real*8 ann(ngf,ngf+1),bn(ngf),xn(ngf)
  61.       integer ipzero(ngf) 
  62.       real*8 epspm
  63.       logical METHODNL,travail_fissure
  64.  
  65. c     variables locales      
  66.       real*8 epsrf3(3),gamr3(3)
  67.       real*8 epsr,ynorm,Teta,Tetaref
  68.       real*8 kt,km,kf,mucr,mu0,ATR1,epsmk,cc0,depsm,dsigr,f1,sigeq,f2
  69.       real*8 epsemin,dtk,xkk,ake,akk,bk,xmm,ame,amm,bm
  70.       real*8 depse1,depse2,depspm
  71.       real*8 tol_syr
  72.       real*8 epsre_nl,siger
  73.  
  74.       real*8 epspm3(3),gamp3(3),depspm3(3),dgamp3(3)
  75.       real*8 epsmf3(3),gamf3(3),epsm03(3),gam03(3)
  76.       parameter (epsemin=1.0d-6,tol_syr=1.000001)
  77.       logical plr_iso  
  78. c     prepa resolution par Jacobien
  79.       real*8 Jf(3,3)
  80.  
  81.  
  82.  
  83. c     intialisation      
  84.       errr=.false.
  85.       garder=.true.
  86.       Pb_nl0=Pb_nl
  87.  
  88. c     plasticité isotrope
  89. c      plr_iso=.true.      
  90. c      plr_iso=.false.       
  91. c      if ((istep.eq.0).or.(istep.eq.2))then
  92. c       print*,'rnfr3d theta',theta,'ttaref',ttaref
  93. c       print*,'hplr,tomr,ekr,skr,ATRref,khir,gamr,sprec,teta1,teta2,dt' 
  94. c       print*, hplr,tomr,ekr,skr,ATRref,khir,gamr,sprec,teta1,teta2,dt
  95. c       print*,'DEBUT theta,eprm0,ttaref,rhor'
  96. c       print*, theta,eprm0,ttaref,rhor
  97. c       print*,'eprk0,eprkf,tokr,K'
  98. c       print*, eprk0,eprkf,tokr,K
  99. c       read*
  100. c      end if
  101.  
  102. c********* deformations imposee par la matrice en adherence parfaite ***
  103.  
  104. c     ** deformation finale dans la matrice suivant l axe du renfort ***      
  105.       epsmf=0.d0       
  106.       do i=1,3
  107.           epsmf3(i)=0.d0
  108.           do j=1,3
  109.               epsmf3(i)=epsmf3(i)+epstf33(i,j)*vecr(numr,j)
  110.           end do
  111.           epsmf=epsmf+epsmf3(i)*vecr(numr,i)
  112.       end do
  113. c     deformation orthogonale 
  114.       ynorm=0.d0
  115.       do i=1,3
  116.           gamf3(i)=epsmf3(i)-epsmf*vecr(numr,i)
  117.           ynorm=ynorm+gamf3(i)**2
  118.       end do
  119. c     tension de la matrice dans la direction du renfort
  120.       epsmf=epsmf+0.5d0*ynorm 
  121.  
  122. c     ** deformation initiale dans la matrice suivant l axe du renfort * 
  123.       epsm0=0.d0       
  124.       do i=1,3
  125.           epsm03(i)=0.d0
  126.           do j=1,3
  127.               epsm03(i)=epsm03(i)+epst033(i,j)*vecr(numr,j)
  128.           end do
  129.           epsm0=epsm0+epsm03(i)*vecr(numr,i)
  130.       end do
  131. c     deformation orthogonale 
  132.       ynorm=0.d0
  133.       do i=1,3
  134.           gam03(i)=epsm03(i)-epsm0*vecr(numr,i)
  135.           ynorm=ynorm+gam03(i)**2
  136.       end do
  137. c     tension initiale de la matrice dans la direction du renfort
  138.       epsm0=epsm0+0.5d0*ynorm 
  139.  
  140. c     ** increment de deformation de la matrice suivant l axe du renfort
  141.       depsm=epsmf-epsm0
  142.  
  143. c********   deformation axiale finale et increment *********************
  144.  
  145.       if ((istep.eq.0).or.(.not.METHODNL)) then 
  146. c        deformation finale dans la renfort en adherence parfaite      
  147.          epsrf=epsmf 
  148. c        increment de deformation du renfort = celui de la matrice
  149.          deps=depsm       
  150.       else if (METHODNL) then
  151. c        etape non local initiale (istep=1) ou intermediaire (istep=3)
  152. c        deformation elastique initiale
  153.          if(ppas) then
  154.              sigr0=spre0
  155.              epsre0=sigr0/your             
  156.          end if
  157. c        deformation finale du renfort      
  158.          if(istep.eq.1) then
  159. c            1er passage en calcul local 
  160. c            deformation elastique finale avec tir elastique
  161.              if(spre0.eq.spref) then
  162. c                pas de variation de precontrainte, on continue le chargement 
  163. c                on adopte l increement de la deformation locale comme
  164. c                tir elastique           
  165.                  epsref=epsre0+depsm
  166.              else
  167. c                la precontrainte vient d etre mise a jour             
  168.                  if(spref.ge.0.) then
  169.                      if(coeff_nlr.ne.1.) then
  170.                          print*,'Remise en tension impossible dans'
  171.                          print*,'rnfr3d car cable endommage'
  172.                          print*,'damr=',1.d0-coeff_nlr 
  173.                          errr=.true. 
  174.                          return
  175.                      else                         
  176.                          epsref=spref/your
  177.                      end if
  178.                  else
  179.                     print*,'Pb spref<0 ds rnfr3d'
  180.                     errr=.true.
  181.                     return
  182.                  end if
  183.              end if        
  184. c            deformation resultante de la somme des increments pour le tir            
  185.              epsrf=epsref+eprk0+eprm0+eprp0            
  186.          else          
  187. c            istep>1 
  188. c            contrainte effective
  189.              siger=sigr_nl/coeff_nlr
  190. c            deformation elastique non locale
  191.              epsre_nl=siger/Your
  192. c            deformation finale supposee  (cf endorf pour la
  193. c            deformation anelastique hypothetique eparv)
  194.              if(plr_iso) then
  195.                  print*,'plasticite iso non prevue dans endorf3d'
  196.                  err1=1
  197.                  errr=.true.
  198.                  return
  199.              else
  200. c                ecrouissage cinematique 
  201. c                utilisation des hypotheses pour le predicteur        
  202.                  if(hypo.eq.1.) then   
  203. c                    hypothese ecoulement positif 
  204.                      if(( siger.gt.(syr+hplr*eprp0)).and.(siger.gt.0.))
  205.      #                   then
  206. c                        hypothese ok                     
  207.                          epsrf=epsre_nl+eparv+(siger-syr)/hplr
  208.                          Pb_nl=0.d0                           
  209.                      else  
  210. c                        hypothese erronee
  211.                          hypo=0.d0
  212. c                        la covergence etant assuree par les tirs elastiques
  213. c                        on annonce pas de Pb_nl mais on ne conserve pas la
  214. c                        solution                         
  215.                          garder=.false.
  216.                          Pb_nl=0.d0                           
  217.                          return                      
  218.                      end if
  219.                  else if (hypo.eq.-1.) then
  220. c                    hypothese ecoulement negatif                 
  221.                      if((siger.lt.(-syr+hplr*eprp0)).and.(siger.lt.0.))
  222.      #                   then
  223. c                        hypothese ok                     
  224.                          epsrf=epsre_nl+eparv+(siger+syr)/hplr
  225.                      else  
  226. c                        hypothese erronee                     
  227.                          hypo=0.d0
  228. c                        la covergence etant assuree par les tirs elastiques
  229. c                        on annonce pas de Pb_nl mais on ne conserve pas la
  230. c                        solution                         
  231.                          garder=.false.
  232.                          Pb_nl=0.d0  
  233.                          return
  234.                      end if
  235.                  else
  236. c                    hypothese tir elastique                  
  237.                      epsrf=epsre_nl+epart
  238. c                    on continue                    
  239.                  end if
  240. c                affichage du traitement  en cas d hypothese fausse               
  241.                  if((hypo*siger).lt.0.) then
  242.                    print*,'dans rnfr3d'                 
  243.                    print*,'hypo',hypo 
  244.                    print*,'siger',siger
  245.                    print*,'epsrf',epsrf
  246.                    print*,'siger',siger
  247.                    print*,'eprp0',eprp0
  248.                    print*,'syr+hplr*eprp0',syr+hplr*eprp0
  249.                    print*,'-syr+hplr*eprp0',-syr+hplr*eprp0
  250.                    print*,'hypo',hypo,'hypo',hypo
  251.                    print*,'(siger-syr)/hplr',(siger-syr)/hplr
  252.                    print*,'(siger+syr)/hplr',(siger+syr)/hplr
  253.                    print*,'hplr',hplr
  254.                    print*,'syr/hplr',syr/hplr
  255. c                  errr=.true.
  256. c                  return
  257.                  end if                  
  258.              end if
  259.          end if         
  260.       else
  261.          print*,'istep & methnl ont des valeurs imprevues ds rnfr3d'
  262.          errr=.true.
  263.          return         
  264.       end if
  265. c     increment de deformation supposee
  266.       deps=epsrf-epsr0
  267.  
  268. c     **** calcul de la contrainte effective dans les renforts *********
  269.  
  270. c     istep: 0 calcul local
  271. c     istep: 1 premiere etape non locale
  272. c     istep: 3 etape non locale intermediaire
  273. c     istep: 2 derniere etape non local
  274.  
  275. c     increment de deformation
  276. c     print*,epsr0,epsrf
  277. c     read*
  278.  
  279.       if ((spre0.ne.spref).and.(.not.ppas)) then
  280. c         on est en train de faire varier la precontrainte
  281. c         on ne calcule que la variation de deformation elastique
  282. c         pour qu elle soit compatible avec la precontrainte imposee
  283. c           print*,'rnfr3d:spre0,spref',spre0,spref
  284.           sigr0=spref
  285.           if(sigr0.gt.0.) then
  286.              sigre0=sigr0/coeff_nlr
  287.           else
  288.              print*,'Pb spref<0 ds rnfr3d'
  289.              errr=.true.
  290.              return
  291.           end if
  292. c         taux de chargement maxi        
  293.           mu00=abs(sigre0/syr)
  294.           mu00=max(mu0,mu00)
  295. c         contrainte effective finale        
  296.           sigref=sigre0
  297. c          print*,'rnfr3d:sigrf',sigrf
  298. c         def elastique
  299.           epse0=sigref/your          
  300. c         actualisation de la def plastique         
  301.           eprpf=eprp0
  302. c         actualisation de la def visqueuse  de Maxwell       
  303.           eprmf=eprm0 
  304. c         actualisation de la def visqueuse  de Kelvin       
  305.           eprkf=eprk0    
  306. c         actualisation deformation plastique cumulee
  307.           epleqf=epleq0          
  308. c         fin de traitement si variation de precontrainte
  309.       else
  310. c         pas de variation de la precontrainte      
  311.           if (ppas) then
  312. c            initialisation de la variable interne si premier pas          
  313. c            print*,'initialisation precontrainte renfort:',sprec
  314.              sigr0=sprec
  315.              sigre0=sigr0
  316. c            initialisation des autres variables internes du renfort
  317.              eprp0=0.d0
  318.              eprm0=0.d0
  319.              eprk0=0.d0
  320.              epleq0=0.d0
  321. c            maj taux de chragement maxi             
  322.              mu0=abs(sigre0/syr)
  323.              mu00=max(mu0,mu00)        
  324. c            read*
  325.           end if
  326. c         la precontrainte n a pas variee donc on calcule normalement        
  327. c         on teste si le calcul de relaxation est nécessaire ou pas
  328. c          if (rhor.gt.0.) then       
  329.              if((dt.eq.0.).or.(tomr.eq.0.).or.(ekr.eq.0.)
  330.      #           .or.(.not.fl3d).or.(plast_seule)) then   
  331. c                calcul elastoplastique : tir elastique
  332.                  sigref=sigre0+your*deps
  333. c                test du critere de plasticite
  334.                  if(plr_iso) then
  335. c                    plasticite isotrope                 
  336.                      f1=sigref-(syr+(hplr*epleq0))
  337.                  else
  338. c                    plasticite cinematique                 
  339.                      f1=(sigref-hplr*eprp0)-syr
  340.                  end if
  341.                  if(f1.gt.0.d0) then
  342. c                   franchissement du critere en traction, on plastifie
  343.                     dep=f1/(your+hplr)
  344.                     sigref=sigref-your*dep
  345. c                   print*,'plastification ds rnfr3d',epsr0,epsrf,
  346. c     #             eprp0,eprpf,your,syr,sigre0,sigrf
  347. c                   read*
  348.                  else
  349.                      if(plr_iso) then 
  350.                          f1=-sigref-(syr+(hplr*epleq0))
  351.                      else    
  352.                          f1=-(sigref-hplr*eprp0)-syr
  353.                      end if                         
  354.                      if (f1.gt.0.d0) then
  355. c                      franchissement du critere en compression, on plastifie
  356.                        dep=-f1/(your+hplr) 
  357.                        sigref=sigref-your*dep
  358.                      else
  359. c                      pas d increment plastique
  360.                        dep=0.d0
  361.                      end if                
  362.                  end if
  363. c                actualisation de la def plastique         
  364.                  eprpf=eprp0+dep
  365.                  depsa=dep
  366. c                actualisation de la def visqueuse         
  367.                  eprmf=eprm0
  368.                  eprkf=eprk0
  369. c                coeff de consolidation debut de pas
  370. c                prise en compt du taux de chargement maxi dans le temps
  371.                  mu0=abs(sigre0/syr)
  372.                  mu00=max(mu0,mu00)  
  373. c                actualisation de la deformation plastique cumulee
  374.                  epleqf=epleq0+abs(dep)                 
  375.               else
  376. c                calcul visco elastique : tir visco elastique
  377. c                print*,'rnfr3d calcul visco elastique'
  378. c                deformation elastique debut de pas
  379.                  epse0=sigre0/your
  380. c                coeff de consolidation debut de pas
  381. c                prise en compt du taux de chargement
  382.                  mu0=min(abs((sigre0-hplr*eprp0)/syr),1.d0)
  383.                  mu00=max(mu0,mu00)
  384.                  if(khir.gt.1.) then
  385.                      mucr=0.6667d0*khir/(khir-1.d0)
  386.                      km=mucr/(mucr-mu0)
  387.                      if(km.lt.0.) then
  388.                           print*,'km<0 ds rnfr3d',km
  389.                           print*,mu0,sigr0,syr
  390.                           errr=.true.
  391.                           return
  392.                      end if
  393.                  else
  394.                      km=1.d0
  395.                  end if
  396. c                prise en compte de la temperature
  397. c                energie d activation fonction du chragement
  398.                  ATR1=ATRref*exp(gamr*max(mu00,xmuthr))
  399.                  Teta=0.5d0*(teta1+teta2)+273.15d0
  400.                  Tetaref=ttaref+273.15d0
  401. c                coeff  d amplification de la def differee         
  402. c                kt=exp(-(ATR1/8.31d0)*(1.d0/Teta-1.d0/Tetaref)) 
  403.                  if(Teta.gt.Tetaref) then 
  404.                      kt=exp(ATR1*((Teta-Tetaref)**xnr))
  405.                  else
  406.                      if(Teta.lt.Tetaref) then
  407.                         kt=exp(-ATR1*((Tetaref-Teta)**xnr))
  408.                      else
  409.                         kt=1.d0
  410.                      end if
  411.                  end if
  412.                  if(kt.gt.200.) then
  413.                     print*,'Amplification thermique rnfr3d :',kt
  414.                     print*,'Valeur tres grande'
  415.                     print*,'verifier les donnees ATR et XNR'
  416.                     print*,ATRref,Teta,xnr,Tetaref,mu00
  417.                     errr=.true.
  418.                  end if
  419. c                potentiel de deformation differee
  420.                  epsmk=ekr*km*kt
  421. c                deformation elastique caracteristique
  422.                  epsek=skr/your
  423. c                coefficient de deformation differee
  424.                  kf=epsmk/epsek 
  425. c                print*,'rnfr3d kf:',kf,epsmk,epsek,ekr,km,kt         
  426. c                coeff de consolidation
  427.                  if (abs(epse0).gt.epsemin) then
  428.                      x1=eprm0/epse0
  429.                  else
  430.                      x1=eprm0/sign(epsemin,epse0)
  431. c                    print*,'rnfr3d:',sign(epsemin,epse0)
  432.                  end if
  433.                  if(x1.gt.0.d0) then
  434. c                    deformation elastique dans le sens de la def de consolidation : ok         
  435.                      cc0=(exp(x1/kf))/kf
  436. c                    print*,x1,kf,cc0
  437.                  else
  438. c                    les deformation epse et epsm sont de signe opposée: 
  439. c                    il faut attendre que la def visqueuse revienne du bon signe
  440. c                    en bloquant la consolidation a sa valeur minimale         
  441.                      cc0=1.d0/kf
  442.                  end if
  443. c                print*,'rnfr3d cc0:',cc0
  444. c                resolution du tir visco elastique
  445. c                depsm=dt*(epse0+theta*deps)/(tomr*cc0*(1.d0+theta*dt/(tomr*cc0)))
  446.                  taum=tomr*cc0
  447.                  tauk=tokr/kt
  448.                  K=Kr/kt
  449. c                construction du Jacobien de la solution analytique
  450. c                Jf=d(depse,depsk)/d(v,epse0,epsk0)
  451.                  call jacobflu3d(tauk,K,taum,0.d0,1.d0,dt,Jf,err1) 
  452.                  if(err1.eq.1) then
  453.                      print*,'Erreur prepa Jacobienne ds rnfr3d'
  454.                      errr=.true.
  455.                      return
  456.                  end if
  457. c               *** tir : calcul des increments sans deformation plastique              
  458. c               preparation de la matrice de couplage 
  459.                 do i=1,3
  460.                     do j=1,3           
  461.                         if(i.eq.1) then
  462.                             ann(i,j)=1.d0
  463.                         else
  464.                             ann(i,j)=0.d0
  465.                         end if
  466.                     end do
  467.                     bn(i)=0.d0
  468.                 end do
  469. c               1ere ligne depse
  470.                 depse1=Jf(1,1)*deps/dt+Jf(1,2)*epse0+Jf(1,3)*eprk0
  471. c               2eme ligne depsk
  472.                 depsk1=Jf(2,1)*deps/dt+Jf(2,2)*epse0+Jf(2,3)*eprk0
  473. c               deduction depsm               
  474.                 depsm1=deps-depse1-depsk1
  475. c               actualisation des variables internes
  476.                 eprmf=eprm0+depsm1
  477.                 eprkf=eprk0+depsk1         
  478.                 eprpf=eprp0
  479.                 epleqf=epleq0
  480.                 dsigr=your*depse1
  481.                 sigref=sigre0+dsigr
  482.                 depsa=depsk1+depsm1
  483. c               print*,'rnfr3d---',sigr0,eprmf,eprf,sigrf         
  484. c               test de plasticité
  485.                 if(plr_iso) then
  486. c                  plasticité isotrope              
  487.                    sigeq=sigref
  488.                 else
  489. c                  plasticite cinematique              
  490.                    sigeq=sigref-hplr*eprpf
  491.                 end if
  492.                 if(sigeq.gt.0.d0) then
  493.                   if(plr_iso) then
  494. c                    plasticite isotrope                  
  495.                      f1=sigeq-(syr+hplr*epleq0)
  496.                      Ep=your
  497.                      Hp=-hplr                     
  498.                   else
  499. c                    plasticite cinematique                  
  500.                      f1=sigeq-syr
  501.                      Ep=your
  502.                      Hp=-hplr
  503.                   end if
  504.                 else
  505.                   if (plr_iso) then
  506. c                    plasticite isotrope
  507.                      f1=-sigeq-(syr+hplr*epleq0)
  508.                      Ep=-your
  509.                      Hp=hplr                     
  510.                   else
  511. c                    plasticite cinematique                  
  512.                      f1=-sigeq-syr
  513.                      Ep=-your
  514.                      Hp=hplr
  515.                   end if
  516.                 end if
  517.                 if(f1.gt.0.d0) then
  518. c                print*,'ecoulement plastique dans rnfr3d'
  519. c                print*,'sigre0,eprmf,eprf,sigrf,f1',sigr0,eprmf,eprf,sigref,f1
  520. c                read*           
  521. c                ecoulement plastique couple a la relaxation
  522. c                preparation de la matrice de couplage 
  523.                  do i=1,4
  524.                   do j=1,4           
  525.                     if(i.eq.1) then
  526.                         ann(i,j)=1.d0
  527.                     else
  528.                         ann(i,j)=0.d0
  529.                     end if
  530.                   end do
  531.                   bn(i)=0.d0
  532.                  end do
  533. c                1ere ligne Elastique        
  534. c                couplage avec depsp
  535.                  ann(1,4)=Jf(1,1)/dt               
  536. c                2eme ligne Kelvin
  537. c                affectation dans la matrice de couplage         
  538.                  ann(2,4)=Jf(2,1)/dt
  539. c                3eme ligne Maxwell
  540.                  ann(3,1)=1.d0
  541.                  ann(3,2)=1.d0
  542.                  ann(3,4)=1.d0 
  543. c                4eme ligne plasticite
  544.                  ann(4,1)=Ep
  545.                  ann(4,4)=Hp
  546.                  bn(4)=-f1         
  547. c                resolution du systeme lineaire                   
  548.                  call gaus3d(4,ann,xn,bn,ngf,errgauss,ipzero) 
  549.                  if(errgauss.eq.1) then
  550.                    errr=.true.
  551.                    print*,'Pb avec gaus3d viscoplasticite ds rnfr3d'
  552.                    return
  553.                  end if 
  554. c                recuperation des solution
  555.                  depse2=xn(1)
  556.                  depsk2=xn(2)
  557.                  depsm2=xn(3)
  558.                  depsp =xn(4)
  559. c                actualisation des variables internes 
  560.                  eprkf=eprkf+depsk2    
  561.                  eprmf=eprmf+depsm2
  562.                  eprpf=eprpf+depsp
  563.                  dsigr=your*depse2
  564.                  sigref=sigref+dsigr 
  565.                  depsa=depsk2+depsm2+depsp 
  566. c                actualisation deformation plastique cumulee
  567.                  epleqf=epleq0+abs(depsp)               
  568. c                verif criteres de plasticite
  569.                  if(plr_iso) then
  570. c                  plasticite isotrope               
  571.                    sigeq=sigref
  572.                  else
  573. c                  plasticite cinematique               
  574.                    sigeq=sigref-hplr*eprpf
  575.                  end if
  576.                  if(sigeq.gt.0.d0) then
  577.                   if(plr_iso) then
  578. c                   plasticite isotrope                
  579.                     if(sigeq.gt.((syr+hplr*epleqf)*tol_syr)) then
  580.                         print*,'Pd de plasticite ds rnfr3d'
  581.                         print*,sigeq,'>',syr+hplr*epleqf
  582.                         errr=.true.
  583.                     end if
  584.                   else
  585. c                   plasticite cimenmatique                
  586.                     if(sigeq.gt.syr*tol_syr) then
  587.                         print*,'Pd de plasticite ds rnfr3d'
  588.                         print*,sigeq,'>',syr
  589.                         errr=.true.
  590.                     end if
  591.                   end if
  592.                 else
  593.                  if (plr_iso) then
  594. c                   plasticite isotrope
  595.                     if(sigeq.lt.-((syr+hplr*epleqf)*tol_syr)) then
  596.                         print*,'Pd de plasticite ds rnfr3d'
  597.                         print*,sigeq,'<',-(syr+hplr*epleqf)
  598.                         errr=.true.
  599.                     end if  
  600.                  else
  601. c                   plasticite cinematique                
  602.                     if(sigeq.lt.-syr*tol_syr) then
  603.                         print*,'Pd de plasticite ds rnfr3d'
  604.                         print*,sigeq,'<',-syr
  605.                         errr=.true.
  606.                     end if 
  607.                  end if                    
  608.                 end if
  609.                 if(errr) then
  610.                  do i=1,4
  611.                   do j=1,4
  612.                     print*,'ann(',i,j,')=',ann(i,j)
  613.                   end do
  614.                   print*,'bn(',i,')=',bn(i)
  615.                  end do
  616.                  return
  617.                 end if
  618.              end if   
  619.             end if 
  620. c           else
  621. c            print*,'Renfort avec section nulle ds rnfr3d'
  622. c           si rhor = 0           
  623. c            sigre0=0.d0
  624. c           taux de chargement maxi        
  625. c            mu00=abs(sigre0/syr)
  626. c           contrainte finale        
  627. c            sigref=sigre0
  628. c           actualisation de la def plastique         
  629. c            eprpf=eprp0
  630. c           actualisation de la def visqueuse  de Maxwell       
  631. c           eprmf=eprm0 
  632. c           actualisation de la def visqueuse  de Kelvin       
  633. c            eprkf=eprk0
  634. c           deformation plastique equivalente            
  635. c            epleqf=0.0d0            
  636. c           fin de traitement si rho=0
  637. c            continue     
  638. c           end if       
  639.       end if   
  640.  
  641. c***********************************************************************     
  642.  
  643.  
  644. c      print*,'renfort num:',numr 
  645. c     if(istep.eq.2) then     
  646. c     print*,'FIN theta,eprm0,eprmf,ttaref,rhor'
  647. c     print*, theta,eprm0,eprmf,ttaref,rhor
  648. c     print*,'-----------------------------------------------'
  649. c      print*,'dans rnfr3d'
  650. c*    read*  
  651. c     end if   
  652. C       print*,'fin rnfr3d, renfort :', numr
  653. C       print*,'epsr0',epsr0
  654. C       print*,'epsrf',epsrf
  655. C       print*,'eprp0',eprp0,'eprpf',eprpf
  656. C       print*,'eprm0',eprm0,'eprmf',eprmf
  657. C       print*,'mu00',mu00
  658. C       print*,'sigre0',sigre0,'sigref',sigref 
  659. C       print*,'------------------------------------------------'
  660. C       read*
  661.  
  662. c******** controle de l hypothese d ecoulement *************************
  663.  
  664.       if(istep.ne.1) then   
  665.         if(Pb_nl0.eq.0.) then
  666. c           les hypothese sont stabilisees au niveau GLOBAL: on les garde                 
  667.             Pb_nl=0.d0
  668. c           on peut repartir avec ces valeurs
  669.             garder=.true.
  670.         else
  671. c           on ne garde pas encore la solution car certains EF changent
  672. c           d hypothese        
  673.             garder=.false.
  674.             if((eprpf.gt.eprp0).and.(hypo.eq.1.))then
  675. c               les hypotheses locales sont correctes mais pb niveau global
  676.                 hypo=1.d0
  677. c               mais ok au niveau local, on en informe le niveau global
  678.                 Pb_nl=0.d0 
  679.             else if((eprpf.lt.eprp0).and.(hypo.eq.-1.))then
  680. c               les hypotheses locales sont correctes mais pb niveau global
  681.                 hypo=-1.d0
  682. c               mais ok au niveau local, on en informe le niveau global
  683.                 Pb_nl=0.d0  
  684.             else if((hypo.ne.0.).and.travail_fissure)then
  685. c               on garde l hypothese d ecoulement plastique dans
  686. c               la fissure en cours d 'ouverture quoi qu'il arrive
  687.                 Pb_nl=0.d0                                 
  688.             else if((eprpf.eq.eprp0).and.(hypo.ne.0.).and.
  689.      #          (.not.travail_fissure))then
  690. c               hypothese erronee sur ecoulement hors fissure
  691.                 hypo=0.d0
  692. c               on conserve hypo=0 hors fissure pour la suite                
  693.                 Pb_nl=0.d0
  694.             else if((abs(eprpf-eprp0).lt.epsemin).and.(hypo.eq.0.))
  695.      #          then
  696. c               hypothese quasi ok
  697.                 Pb_nl=0.d0
  698.             else if((abs(eprpf-eprp0).ge.epsemin).and.(hypo.eq.0.))
  699.      #          then     
  700.                 if(travail_fissure) then
  701. c                   on privilegie l hypothese d ecoulement                
  702.                     if (eprpf.gt.eprp0) then
  703. c                        il y a eu plastification avec hypo=0                   
  704.                          hypo=1.d0
  705. c                        on gardera cette hypothese                         
  706.                          Pb_nl=0.d0
  707.                     else 
  708. c                        il y a eu plastification avec hypo=0                   
  709.                          hypo=-1.d0
  710. c                        on gardera cette hypothese                         
  711.                          Pb_nl=0.d0
  712.                     end if
  713.                 else
  714. c                   on gardera hypo=0              
  715.                     hypo=0.d0
  716.                     Pb_nl=0.d0
  717.                 end if                        
  718.             else          
  719.                  print*,'cas imprevu dans endorf'
  720.                  print*,'eprpf',eprpf
  721.                  print*,'eprp0',eprp0                
  722.                  print*,'hypo_actu',hypo
  723.                  print*,'Pb_nl',Pb_nl
  724.                  print*,'travail_fissure',travail_fissure
  725.                  err1=1
  726.                  errr=.true.
  727.                  return         
  728.             end if
  729.         end if
  730.       else
  731. c        istep=1 on impose les hypotheses d ecoulement pour la
  732. c        premiere etape non locale en fonction du tir local
  733. c         if (eprpf.eq.eprp0) then
  734. c             on suppose pas  d ecoulement  plastique        
  735.               hypo=0.d0
  736. c         else if (eprpf.gt.(eprp0+epsemin)) then
  737. c             ecoulement plastique  positif
  738. c              hypo=1.d0
  739. c         else if (eprpf.lt.(eprp0-epsemin)) then
  740. c             ecoulement plastique negatif         
  741. c              hypo=-1.d0
  742. c         end if
  743. c        on repart de la solution precedente pour construire 
  744. c        l operateur tangent de NRM dans UNPAS avec la solution convergee 
  745. c        du pas precedent sauf pour le premier pas ou on garde le tir local
  746.          if (.not.ppas) then         
  747.             garder=.false.
  748.          else
  749.             garder=.true.
  750.          end if
  751. c        on ne sait pas encore si l hypothese est coherente         
  752.          Pb_nl=0.d0
  753.       end if         
  754.  
  755. c***********************************************************************      
  756.  
  757.       return
  758.  
  759.       end     
  760.  
  761.  
  762.  
  763.  
  764.  

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