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    FD218221  24/02/07    21:15:09     11834          
  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.          if(vdef.eq.0.) then       
  100.             print*,'Pb dans les donnees de DEF dans def3d'
  101.             print*,'nsul=0 et vdef=0 ds vdef3d' 
  102.             err1=1
  103.             return
  104.          end if
  105.        end if   
  106.  
  107. c       print*, 'coeff efficacite def dans def3d', coeff1       
  108.  
  109. c       print*,'Def3d Sc,Ac mol/m3',Sc, AC
  110.  
  111. c      initialisation de At,St,M1,E1,vdef0,E2,HG si premier pas
  112.        if(ppas) then
  113. c        on surcharge les variables internes pour initialiser
  114. c        le pr chimique de la def en fonction des donnees materiaux       
  115.          if(Sc.ge.3.d0*Ac) then
  116.             E1=Ac
  117.             M1=0.d0
  118.             E2=0.d0
  119.             vdef0=0.d0
  120.             At=0.d0
  121.             St=SC-3.d0*E1
  122.          else if (Sc.le.Ac) then
  123.             E1=0.d0
  124.             M1=Sc
  125.             E2=0.d0
  126.             vdef0=0.d0
  127.             At=Ac-M1
  128.             St=0.d0
  129.          else
  130.             E1=(Sc-Ac)/2.d0
  131.             M1=(3.d0*Ac-Sc)/2.d0
  132.             E2=0.d0
  133.             vdef0=0.d0
  134.             At=0.d0
  135.             St=0.d0
  136.          end if
  137.       end if
  138. c      print*,'Def3d:At,St,M1,E1,vdef0,E2'
  139. c      print*,At,St,M1,E1,vdef0,E2
  140.  
  141. c     sulfo-aluminates destabilisable en temperature       
  142.       A1=E1+M1+E2
  143.       S1=3.d0*(E1+E2)+M1
  144.  
  145. c     calcul des coeff cinetiques pour la dissolution/ precipitation       
  146. c     passage des temperatures en Kelvin
  147.       temp1=0.5d0*(teta1+teta2)+273.15d0
  148.  
  149. c     temperature seuil de destabilisation fonction de la concentration en Na
  150.       if(Cna.ge.Cnak) then
  151.          tempd=273.15d0+ttdd*((Cnak/CNa)**exnd)
  152.       else
  153.          tempd=273.15d0+ttdd
  154.       end if
  155. c     temperature minimale pour la fixation des alu
  156.       tempf=273.15d0+ttdf     
  157.  
  158.       if(temp1.gt.tempd) then
  159. c      dissolution des phases sulfo alumineuse
  160. c      coeff cinetique de dissolution precipitation
  161. c      influence de la temperature
  162.        CTD=max((exp(-(nrjd/8.31d0)*(1.d0/temp1-1.d0/tempd)))-1.d0,0.d0)
  163. c      influence de la temperature sur la fixation des alu       
  164.        CTF=max((exp(-(nrjf/8.31d0)*(1.d0/temp1-1.d0/tempf)))-1.d0,0.d0)
  165. c      influence des alcalins sur la dissolution
  166.        CCD=CNa/Cnak 
  167. c      influence des alcalins sur la fixation des alu
  168.        if(CNa.gt.0.) then       
  169.           CCF=(Cnak/CNa)**exmd
  170.        else
  171. c         fixation instantanee
  172.           CCF=1.d6
  173.        end if          
  174. c      influence de l'humidite sur la dissolution
  175.        CHD=1.d0
  176. c      influence de l humidite sur la fixation       
  177.        CHF=1.d0  
  178. c      influence THC sur dissolution
  179.        CD=CCD*CTD*CHD
  180. c      influence THC sur la fixation
  181.        CF=CCF*CTF*CHF         
  182. c      actualisation de l aluminium ds les phases sulfates
  183.        A1F=A1*exp(-dt*CD/tdid)
  184. c      actualisation de l aluminium libre
  185.        if(CF.ne.0.) then
  186. c          la fixation de l alu est possible 
  187. c          rapport des temps caractéristiques        
  188.            kfd=(tfid/tdid)*(CD/CF) 
  189.            if(kfd.eq.1.) then
  190.                  kfd=1.0001d0
  191.                  CD=kfd*CD 
  192.            end if            
  193.            AtF=At*exp(-dt*CF/tfid)+(A1*kfd/(kfd-1.d0))*
  194.      #    (exp(-dt*CF/tfid)-exp(-dt*CD/tdid))
  195.        else
  196. c          la fixation des alu est impossible, tout l alu dissous reste
  197. c          disponible pour la reaction secondaire
  198.            AtF=At+(A1-A1F)
  199.        end if
  200. c      actualisation des sulfates libres
  201. c      toutes les phases disparaissent avec la meme cinetique
  202.        StF=St+(3.d0*(E1+E2)+M1)*(1.d0-exp(-dt*CD/tdid)) 
  203. c      actualisation des variables internes
  204.        if(A1.ne.0.) then
  205.          cA1=A1F/A1
  206.        else
  207.          cA1=1.d0
  208.        end if         
  209.        E1F=E1*cA1
  210.        E2F=E2*cA1
  211.        M1F=M1*cA1        
  212.        dvdef=(E2F-E2)*VMaft*coeff1
  213.       else
  214. c       precipitation des phases sulfo-alumineuses
  215. c       temprerature de reference pour la precipitation      
  216.         tempr=ttrp+273.15d0 
  217. c       coeff cinetique pour la precipitation
  218. c       influence de le temperature
  219.         CTP=((1.d0-exp(-(nrjd/8.31d0)*(1.d0/temp1-1.d0/tempd)))/
  220.      #  (1.d0-exp(-(nrjd/8.31d0)*(1.d0/tempr-1.d0/tempd))))*
  221.      #  exp(-(nrjp/8.31d0)*(1.d0/temp1-1.d0/tempr))
  222. c       influence de l humidite 
  223.         if(srsdef.lt.1.d0) then    
  224.            CHP=exp(-(1.d0-Sr1)/(1.d0-srsdef))  
  225.         else
  226.            CHP=0.d0
  227.         end if
  228. c       influence des alcalins
  229.         CCP=(1.d0-min(CNa,Cnab)/Cnab)**exmd
  230.         CP=CTP*CCP*CHP
  231.         if((CP.gt.0.).and.(St.gt.0.)) then
  232. c          temps caracteristique de precipitation         
  233.            taupre=tdef/CP
  234.            if(M1.gt.0.) then
  235. c              il faut d abord convertir m1 
  236.                if (St.gt.(2.d0*M1)) then
  237. c                  il ya assez de st pour consommer tout m1
  238.                    dt1=-(taupre/2.d0)*log(1.d0-(2.d0*M1/St))
  239. c                   print*,'def3d: st>2m1: dt1=',dt1
  240.                    if (dt1.ge.dt) then
  241. c                      pas le temps de tout concertir
  242.                        E2F=E2+(St/2.d0)*(1.d0-exp(-2.d0*dt/taupre))
  243.                        StF=St-2.d0*(E2F-E2)
  244.                        AtF=At
  245.                        M1F=M1-(E2F-E2)
  246.                        E1F=E1 
  247.                    else
  248. c                      on a le temps de consommer tout m1 et une partie de At
  249. c                      on va jusque a t0+dt1 en consommant m1
  250.                        E21=E2+(St/2.d0)*(1.d0-exp(-2.d0*dt1/taupre))
  251.                        St1=St-2.d0*(E21-E2)
  252.                        At1=At
  253.                        M1F=M1-(E21-E2)
  254.                        E1F=E1
  255. c                      ensuite on crée du E2 a partir de At et St uniquement
  256.                        if(St1.le.(3.d0*At1)) then
  257. c                        pas assez de St pour concertir tout At
  258.                          dt2=dt-dt1
  259.                          E2F=E21+(St1/3.d0)*(1.d0-exp(-3.d0*dt2/taupre))
  260.                          StF=St1-3.d0*(E2F-E21)
  261.                          AtF=At1-(E2F-E21)
  262.                        else
  263. c                        pas assez de At pour consommer tout St
  264.                          dt2=-(taupre/3.d0)*log(1.d0-(3.d0*At1/St1))
  265. c                         print*,'def3d: temps pour at=0',dt2
  266. c                         print*,taupre,at1,st1
  267.                          dt2=min(dt-dt1,dt2)
  268.                          E2F=E21+(St1/3.d0)*(1.d0-exp(-3.d0*dt2/taupre))
  269.                          AtF=At1-(E2F-E21) 
  270.                          StF=St1-3.d0*(E2F-E21)                         
  271.                        end if                       
  272.                    end if                       
  273.                else
  274. c                  il n y a pas assez de st pour consommer tout m1 
  275. c                  il faut un temps infini pour consommer m1, les autres cas
  276. c                  ne sont donc pas à considerer
  277.                    E2F=E2+(St/2.d0)*(1.d0-exp(-2.d0*dt/taupre))
  278.                    StF=St-2.d0*(E2F-E2)
  279.                    AtF=At
  280.                    M1F=M1-(E2F-E2)
  281.                    E1F=E1                                
  282.                end if               
  283.            else
  284. c              il n y a plus de m1 a convertir  
  285.                if(St.le.(3.d0*At)) then
  286. c                pas assez de St pour concertir tout At
  287.                  E2F=E2+(St/3.d0)*(1.d0-exp(-3.d0*dt/taupre))
  288.                  StF=St-3.d0*(E2F-E2)
  289.                  AtF=At-(E2F-E2)
  290.                  M1F=0.d0
  291.                  E1F=E1
  292.                else
  293. c                pas assez de At pour consommer tout St
  294.                  dt2=-(taupre/3.d0)*log(1.d0-(3.d0*At/St))
  295. c                 print*,'def3d: st<3At et m1=0: dt2=',dt2                 
  296.                  dt2=min(dt,dt2)
  297.                  E2F=E2+(St/3.d0)*(1.d0-exp(-3.d0*dt2/taupre))
  298.                  Stf=St-3.d0*(E2F-E2)
  299.                  AtF=At-(E2F-E2)
  300.                  M1F=0.d0
  301.                  E1F=E1                       
  302.                end if           
  303.            end if 
  304. c          evolution de vdef par precipiatation de phases neoformees
  305. c          a l endroit d afm1 pour la conversion directe
  306. c          a un nouvel endroit pour aft2
  307.            dvdef=(E2F-E2)*VMaft*coeff1+(M1F-M1)*VMafm
  308.            if (coeff1.eq.0.) then
  309.              dvdef=0.d0
  310.            end if
  311.         else
  312.            E1F=E1
  313.            E2F=E2
  314.            M1F=M1        
  315.            AtF=At
  316.            StF=St
  317.            dvdef=0.d0           
  318.         end if
  319.        end if  
  320.  
  321. C      calcul du volume des phases neoformees
  322.        VDEFF=VDEF0+DVDEF 
  323. c       VDEF1=0.d0       
  324.        def1=VDEFF/vdef00
  325.        dphidef=dvdef    
  326.  
  327. c     **** calcul des coeffs poromécaniques ****************************       
  328.  
  329.       if(vdeff.ne.0.) then
  330.          vdef1=vdeff
  331. c        volume de gel a considerer
  332.          bs1=2.d0*vdef1/(1.d0+vdef1)
  333.          bdef=0.5d0*(bs0+bs1)
  334.          Ks1=Kb/(1.d0-bdef)         
  335.          if (vdef1.gt.0.d0) then
  336.             mdef=((bdef-vdef1)/Ks1+vdef1/Kdef)**(-1)
  337.          else
  338. c           pas de def
  339.             bs1=0.d0
  340.             bdef=0.5d0*(bs0+bs1)
  341.             mdef=0.d0
  342.          end if  
  343. c        prise en compte du volume des vides connectés         
  344.          mdef=mdef/(1.0d0+mdef*Vvdef*Cdef/Rt)       
  345.       else
  346. c       pas de def donc il ne faur pas generer de pression
  347.         bdef=0.d0
  348.         bs1=0.d0
  349.         mdef=0.d0
  350.       end if  
  351.  
  352.  
  353. c       print*,'e1', E1,E1F
  354. c       print*,'e2',E2,E2F
  355. c       print*,'m1',M1,M1F       
  356. c       print*,'at',At,AtF
  357. c       print*,'st',St,StF
  358. c       print*,'vdef',vdef0,vdeff
  359. c       print*,'adef',def1
  360. c       read*
  361.  
  362. c      print*, "vdef00,def1,vdef1"
  363. c      print*,vdef00,def1,vdef1 
  364. c      read*      
  365.  
  366.       return
  367.       end
  368. c***********************************************************************      
  369.  
  370.  
  371.  
  372.  

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