Télécharger la plaquette Cast3M


fr en

Télécharger def3d.eso

Retour à la liste

Numérotation des lignes : 

def3d

  1. C DEF3D     SOURCE    PV090527  24/06/12    21:15:03     11940          
  2.       subroutine def3d(ppas,tdef,nrjd,def0,sr1,srsdef,teta1,teta2,dt,
  3.      # vdef00,def1,CNa,nrjp,ttrp,tfid,ttdd,tdid,exmd,exnd,cnab,cnak,
  4.      # ssad,At,St,M1,E1,vdef0,E2,AtF,StF,M1F,E1F,vdeff,E2F,ttdf,nrjf,
  5.      # nsul,err1,Kb,Kdef,Vvdef,Cdef,Rt,mdef,bdef,dphidef,bs0,bs1)
  6.  
  7. c      sous programme de calcul de l avancement chimique de def
  8. c      version de mars 2017 sans possibilite de formation de M2 et
  9. c      conversion possible de M1 en E2
  10.  
  11.  
  12. c   ********************************************************************
  13.        implicit real*8 (a-h,o-z)
  14.        implicit integer (i-n)
  15. c   ********************************************************************
  16.       logical ppas
  17.       real*8 tdef,nrjd,def0,sr1,srsdef,teta1,teta2      
  18.       real*8 alphadef,Ear, xk,CNamin,CNa,ARNa,CNamin1
  19.       real*8 nrjp,ttrp,tfid,ttdd
  20.       real*8 tdid,exmd,exnd,cnab,cnak,ssad,vdef0,vdeff
  21.       real*8 At,St,M1,E1,E2,AtF,StF,M1F,E1F,E2F,ttdf,nrjf,nsul
  22.       integer err1
  23.       real*8 Kb,Kdef,Vvdef,Cdef,Rt,Mdef,Bdef,dphidef,bs0,bs1
  24.       real*8 Ks1
  25.  
  26.       real*8 A1,S1
  27.       real*8 VMaft,VMafm,Sc,Ac,E21,AT1,ST1
  28.       real*8 coeff1,taudis,taupre,taufix,kfd,dvdef
  29.       real*8 CCD,CCF,CCP,CTD,CTF,CTP,CHD,CHF,CHP,CD,CF,CP
  30.       real*8 vdef1
  31.  
  32.  
  33. c      print*,"tdef,nrjd,def0,sr1,srsdef,teta1,teta2,dt"
  34. c      print*,tdef,nrjd,def0,sr1,srsdef,teta1,teta2,dt     
  35. c      print*,"nrjp,ttrp,tfid,ttdd"
  36. c      print*,nrjp,ttrp,tfid,ttdd      
  37. c      print*,"tdid,exmd,exnd,cnab,cnak,ssad"
  38. c      print*,tdid,exmd,exnd,cnab,cnak,ssad
  39. c      read*
  40.  
  41. c     **** avancement chimique de la def *******************************
  42.  
  43. c     cf .Chemical modelling of Delayed Ettringite Formation for assessment
  44. c     of affected concrete structures (Sellier / Multon)
  45. c     in  https://doi.org/10.1016/j.cemconres.2018.03.006
  46.  
  47.  
  48. c     teneur en aluminium et sulfates en fonction de VDEF et SSAD
  49.       VMaft=715.d-6
  50.       VMafm=254.6d-6
  51.       err1=0     
  52.       if(ssad.ge.3.d0) then  
  53. c        assez de sulfates pour ne faire que de la def      
  54.          Ac=vdef00/VMaft
  55.          SC=ssad*Ac
  56.       else
  57. c        les sulfates limitent la formation de la def      
  58.          Sc=3.d0*vdef00/VMaft
  59.          if(ssad.ne.0.) then
  60.             Ac=Sc/ssad
  61.          else
  62.            if (vdef00.eq.0.) then
  63.              Sc=0.d0
  64.              Ac=0.d0
  65.            else
  66. c            on suppose qu il y a juste assez d aluminium pour ne faire
  67. c            que de la def           
  68.              Ac=Sc/3.d0
  69.            end if
  70.          end if
  71.        end if
  72.        coeff1=1.d0
  73.  
  74.        if (nsul.ne.0.) then
  75.          if (Sc.ne.0.) then
  76.             if (nsul.le.Sc) then
  77. c             le volume de def calcule sera apparent( il contient des vides)
  78. c             le coeff d efficacite est superieur a un pour les AFT         
  79.               coeff1=Sc/nsul
  80. c             mais la chimie est conforme a nsul et ssad          
  81.               Sc=nsul
  82.               Ac=Sc/ssad
  83.             else
  84.               print*,' Donnees NSUL et VDEF incompatibles dans def3d'
  85.               print*,'NSUL > Sc(VDEF) theoriquement impossible'
  86.               err1=1
  87.               return
  88.             end if
  89.          else
  90. c          comme Sc=0 on prend coeff1=0. et Sc=nsul, la chimie tourne
  91. c          mais la def n est pas efficace
  92.            Sc=nsul
  93.            Ac=Sc/ssad  
  94.            coeff1=0.d0           
  95.          end if         
  96.        else
  97. c        nsul=0 mais vdef ne 0
  98. c        on suppose que c est vdef qui est la donnee correcte
  99. * tentative correction PV: vdef n'existe pas
  100. **       if(vdef.eq.0.) then       
  101.          if(vdef00.eq.0.) then       
  102.             print*,'Pb dans les donnees de DEF dans def3d'
  103.             print*,'nsul=0 et vdef=0 ds vdef3d' 
  104.             err1=1
  105.             return
  106.          end if
  107.        end if   
  108.  
  109. c       print*, 'coeff efficacite def dans def3d', coeff1       
  110.  
  111. c       print*,'Def3d Sc,Ac mol/m3',Sc, AC
  112.  
  113. c      initialisation de At,St,M1,E1,vdef0,E2,HG si premier pas
  114.        if(ppas) then
  115. c        on surcharge les variables internes pour initialiser
  116. c        le pr chimique de la def en fonction des donnees materiaux       
  117.          if(Sc.ge.3.d0*Ac) then
  118.             E1=Ac
  119.             M1=0.d0
  120.             E2=0.d0
  121.             vdef0=0.d0
  122.             At=0.d0
  123.             St=SC-3.d0*E1
  124.          else if (Sc.le.Ac) then
  125.             E1=0.d0
  126.             M1=Sc
  127.             E2=0.d0
  128.             vdef0=0.d0
  129.             At=Ac-M1
  130.             St=0.d0
  131.          else
  132.             E1=(Sc-Ac)/2.d0
  133.             M1=(3.d0*Ac-Sc)/2.d0
  134.             E2=0.d0
  135.             vdef0=0.d0
  136.             At=0.d0
  137.             St=0.d0
  138.          end if
  139.       end if
  140. c      print*,'Def3d:At,St,M1,E1,vdef0,E2'
  141. c      print*,At,St,M1,E1,vdef0,E2
  142.  
  143. c     sulfo-aluminates destabilisable en temperature       
  144.       A1=E1+M1+E2
  145.       S1=3.d0*(E1+E2)+M1
  146.  
  147. c     calcul des coeff cinetiques pour la dissolution/ precipitation       
  148. c     passage des temperatures en Kelvin
  149.       temp1=0.5d0*(teta1+teta2)+273.15d0
  150.  
  151. c     temperature seuil de destabilisation fonction de la concentration en Na
  152.       if(Cna.ge.Cnak) then
  153.          tempd=273.15d0+ttdd*((Cnak/CNa)**exnd)
  154.       else
  155.          tempd=273.15d0+ttdd
  156.       end if
  157. c     temperature minimale pour la fixation des alu
  158.       tempf=273.15d0+ttdf     
  159.  
  160.       if(temp1.gt.tempd) then
  161. c      dissolution des phases sulfo alumineuse
  162. c      coeff cinetique de dissolution precipitation
  163. c      influence de la temperature
  164.        CTD=max((exp(-(nrjd/8.31d0)*(1.d0/temp1-1.d0/tempd)))-1.d0,0.d0)
  165. c      influence de la temperature sur la fixation des alu       
  166.        CTF=max((exp(-(nrjf/8.31d0)*(1.d0/temp1-1.d0/tempf)))-1.d0,0.d0)
  167. c      influence des alcalins sur la dissolution
  168.        CCD=CNa/Cnak 
  169. c      influence des alcalins sur la fixation des alu
  170.        if(CNa.gt.0.) then       
  171.           CCF=(Cnak/CNa)**exmd
  172.        else
  173. c         fixation instantanee
  174.           CCF=1.d6
  175.        end if          
  176. c      influence de l'humidite sur la dissolution
  177.        CHD=1.d0
  178. c      influence de l humidite sur la fixation       
  179.        CHF=1.d0  
  180. c      influence THC sur dissolution
  181.        CD=CCD*CTD*CHD
  182. c      influence THC sur la fixation
  183.        CF=CCF*CTF*CHF         
  184. c      actualisation de l aluminium ds les phases sulfates
  185.        A1F=A1*exp(-dt*CD/tdid)
  186. c      actualisation de l aluminium libre
  187.        if(CF.ne.0.) then
  188. c          la fixation de l alu est possible 
  189. c          rapport des temps caractéristiques        
  190.            kfd=(tfid/tdid)*(CD/CF) 
  191.            if(kfd.eq.1.) then
  192.                  kfd=1.0001d0
  193.                  CD=kfd*CD 
  194.            end if            
  195.            AtF=At*exp(-dt*CF/tfid)+(A1*kfd/(kfd-1.d0))*
  196.      #    (exp(-dt*CF/tfid)-exp(-dt*CD/tdid))
  197.        else
  198. c          la fixation des alu est impossible, tout l alu dissous reste
  199. c          disponible pour la reaction secondaire
  200.            AtF=At+(A1-A1F)
  201.        end if
  202. c      actualisation des sulfates libres
  203. c      toutes les phases disparaissent avec la meme cinetique
  204.        StF=St+(3.d0*(E1+E2)+M1)*(1.d0-exp(-dt*CD/tdid)) 
  205. c      actualisation des variables internes
  206.        if(A1.ne.0.) then
  207.          cA1=A1F/A1
  208.        else
  209.          cA1=1.d0
  210.        end if         
  211.        E1F=E1*cA1
  212.        E2F=E2*cA1
  213.        M1F=M1*cA1        
  214.        dvdef=(E2F-E2)*VMaft*coeff1
  215.       else
  216. c       precipitation des phases sulfo-alumineuses
  217. c       temprerature de reference pour la precipitation      
  218.         tempr=ttrp+273.15d0 
  219. c       coeff cinetique pour la precipitation
  220. c       influence de le temperature
  221.         CTP=((1.d0-exp(-(nrjd/8.31d0)*(1.d0/temp1-1.d0/tempd)))/
  222.      #  (1.d0-exp(-(nrjd/8.31d0)*(1.d0/tempr-1.d0/tempd))))*
  223.      #  exp(-(nrjp/8.31d0)*(1.d0/temp1-1.d0/tempr))
  224. c       influence de l humidite 
  225.         if(srsdef.lt.1.d0) then    
  226.            CHP=exp(-(1.d0-Sr1)/(1.d0-srsdef))  
  227.         else
  228.            CHP=0.d0
  229.         end if
  230. c       influence des alcalins
  231.         CCP=(1.d0-min(CNa,Cnab)/Cnab)**exmd
  232.         CP=CTP*CCP*CHP
  233.         if((CP.gt.0.).and.(St.gt.0.)) then
  234. c          temps caracteristique de precipitation         
  235.            taupre=tdef/CP
  236.            if(M1.gt.0.) then
  237. c              il faut d abord convertir m1 
  238.                if (St.gt.(2.d0*M1)) then
  239. c                  il ya assez de st pour consommer tout m1
  240.                    dt1=-(taupre/2.d0)*log(1.d0-(2.d0*M1/St))
  241. c                   print*,'def3d: st>2m1: dt1=',dt1
  242.                    if (dt1.ge.dt) then
  243. c                      pas le temps de tout concertir
  244.                        E2F=E2+(St/2.d0)*(1.d0-exp(-2.d0*dt/taupre))
  245.                        StF=St-2.d0*(E2F-E2)
  246.                        AtF=At
  247.                        M1F=M1-(E2F-E2)
  248.                        E1F=E1 
  249.                    else
  250. c                      on a le temps de consommer tout m1 et une partie de At
  251. c                      on va jusque a t0+dt1 en consommant m1
  252.                        E21=E2+(St/2.d0)*(1.d0-exp(-2.d0*dt1/taupre))
  253.                        St1=St-2.d0*(E21-E2)
  254.                        At1=At
  255.                        M1F=M1-(E21-E2)
  256.                        E1F=E1
  257. c                      ensuite on crée du E2 a partir de At et St uniquement
  258.                        if(St1.le.(3.d0*At1)) then
  259. c                        pas assez de St pour concertir tout At
  260.                          dt2=dt-dt1
  261.                          E2F=E21+(St1/3.d0)*(1.d0-exp(-3.d0*dt2/taupre))
  262.                          StF=St1-3.d0*(E2F-E21)
  263.                          AtF=At1-(E2F-E21)
  264.                        else
  265. c                        pas assez de At pour consommer tout St
  266.                          dt2=-(taupre/3.d0)*log(1.d0-(3.d0*At1/St1))
  267. c                         print*,'def3d: temps pour at=0',dt2
  268. c                         print*,taupre,at1,st1
  269.                          dt2=min(dt-dt1,dt2)
  270.                          E2F=E21+(St1/3.d0)*(1.d0-exp(-3.d0*dt2/taupre))
  271.                          AtF=At1-(E2F-E21) 
  272.                          StF=St1-3.d0*(E2F-E21)                         
  273.                        end if                       
  274.                    end if                       
  275.                else
  276. c                  il n y a pas assez de st pour consommer tout m1 
  277. c                  il faut un temps infini pour consommer m1, les autres cas
  278. c                  ne sont donc pas à considerer
  279.                    E2F=E2+(St/2.d0)*(1.d0-exp(-2.d0*dt/taupre))
  280.                    StF=St-2.d0*(E2F-E2)
  281.                    AtF=At
  282.                    M1F=M1-(E2F-E2)
  283.                    E1F=E1                                
  284.                end if               
  285.            else
  286. c              il n y a plus de m1 a convertir  
  287.                if(St.le.(3.d0*At)) then
  288. c                pas assez de St pour concertir tout At
  289.                  E2F=E2+(St/3.d0)*(1.d0-exp(-3.d0*dt/taupre))
  290.                  StF=St-3.d0*(E2F-E2)
  291.                  AtF=At-(E2F-E2)
  292.                  M1F=0.d0
  293.                  E1F=E1
  294.                else
  295. c                pas assez de At pour consommer tout St
  296.                  dt2=-(taupre/3.d0)*log(1.d0-(3.d0*At/St))
  297. c                 print*,'def3d: st<3At et m1=0: dt2=',dt2                 
  298.                  dt2=min(dt,dt2)
  299.                  E2F=E2+(St/3.d0)*(1.d0-exp(-3.d0*dt2/taupre))
  300.                  Stf=St-3.d0*(E2F-E2)
  301.                  AtF=At-(E2F-E2)
  302.                  M1F=0.d0
  303.                  E1F=E1                       
  304.                end if           
  305.            end if 
  306. c          evolution de vdef par precipiatation de phases neoformees
  307. c          a l endroit d afm1 pour la conversion directe
  308. c          a un nouvel endroit pour aft2
  309.            dvdef=(E2F-E2)*VMaft*coeff1+(M1F-M1)*VMafm
  310.            if (coeff1.eq.0.) then
  311.              dvdef=0.d0
  312.            end if
  313.         else
  314.            E1F=E1
  315.            E2F=E2
  316.            M1F=M1        
  317.            AtF=At
  318.            StF=St
  319.            dvdef=0.d0           
  320.         end if
  321.        end if  
  322.  
  323. C      calcul du volume des phases neoformees
  324.        VDEFF=VDEF0+DVDEF 
  325. c       VDEF1=0.d0       
  326.        def1=VDEFF/vdef00
  327.        dphidef=dvdef    
  328.  
  329. c     **** calcul des coeffs poromécaniques ****************************       
  330.  
  331.       if(vdeff.ne.0.) then
  332.          vdef1=vdeff
  333. c        volume de gel a considerer
  334.          bs1=2.d0*vdef1/(1.d0+vdef1)
  335.          bdef=0.5d0*(bs0+bs1)
  336.          Ks1=Kb/(1.d0-bdef)         
  337.          if (vdef1.gt.0.d0) then
  338.             mdef=((bdef-vdef1)/Ks1+vdef1/Kdef)**(-1)
  339.          else
  340. c           pas de def
  341.             bs1=0.d0
  342.             bdef=0.5d0*(bs0+bs1)
  343.             mdef=0.d0
  344.          end if  
  345. c        prise en compte du volume des vides connectés         
  346.          mdef=mdef/(1.0d0+mdef*Vvdef*Cdef/Rt)       
  347.       else
  348. c       pas de def donc il ne faur pas generer de pression
  349.         bdef=0.d0
  350.         bs1=0.d0
  351.         mdef=0.d0
  352.       end if  
  353.  
  354.  
  355. c       print*,'e1', E1,E1F
  356. c       print*,'e2',E2,E2F
  357. c       print*,'m1',M1,M1F       
  358. c       print*,'at',At,AtF
  359. c       print*,'st',St,StF
  360. c       print*,'vdef',vdef0,vdeff
  361. c       print*,'adef',def1
  362. c       read*
  363.  
  364. c      print*, "vdef00,def1,vdef1"
  365. c      print*,vdef00,def1,vdef1 
  366. c      read*      
  367.  
  368.       return
  369.       end
  370. c***********************************************************************      
  371.  
  372.  
  373.  
  374.  
  375.  

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