Télécharger la plaquette Cast3M


fr en

Télécharger metall.eso

Retour à la liste

Numérotation des lignes : 

metall

  1. C METALL    SOURCE    CB215821  24/04/12    21:16:42     11897          
  2. C METALL    SOURCE    CB215821  19/06/17    21:15:16     10228          
  3.       SUBROUTINE METALL(iwrk52, iwrmet)
  4. C
  5. C     Modele METALLURGIE : LEBLOND / KOISTINEN
  6. C     Appelle par COML6
  7. C     En entrees : Le segment de travail WRK52 (voir include DECHE)
  8. C                  Le segment IMODEL du modele elementaire (metallurgie)
  9. C     En sortie : ecrit les proportions de phases finales
  10. C     dans le tableau WRK52.VARF
  11. C
  12.       IMPLICIT INTEGER(I-N)
  13.       IMPLICIT REAL*8(A-H,O-Z)
  14.  
  15.       PARAMETER (NPHMAX=20)
  16.       REAL*8     delta(NPHMAX)
  17.       REAL*8     Prop(NPHMAX)
  18.       REAL*8     VAR0L(NPHMAX)
  19.       REAL*8     TYPREA(NPHMAX)
  20.       INTEGER    NPROD(NPHMAX)
  21.       INTEGER    NREAC(NPHMAX)
  22.       REAL*8     CRITER(NPHMAX)
  23.  
  24.       CHARACTER*4 MOT4a,MOT4b,MOT4c
  25.  
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC SMCHAML
  29. -INC CCHAMP
  30. -INC DECHE
  31. -INC CCREEL
  32. -INC SMMODEL
  33. -INC SMLENTI
  34. -INC SMLMOTS
  35.  
  36.       wrk52  = iwrk52
  37.       wrkmet = iwrmet
  38.  
  39.       XPREC  = XZPREC*1000
  40.       ITMAXI = 1000
  41.  
  42. C     Convergence si la norme infini de la matrice des coeffs < 1
  43. C     On se donne de la securite => 0.5 < 1
  44.       borne = 0.5D0
  45.  
  46. C     Theta-methode
  47.       theta = 0.5D0
  48.       utheta= MIN(MAX(1.D0 - theta,0.D0),1.D0)
  49.  
  50.       NPHASE=VAR0(/1)
  51.       NTYPES=TYPES(/2)
  52.       if( NPHASE .gt. NPHMAX ) then
  53.         CALL ERREUR(5)
  54.       endif
  55.  
  56.       do i = 1, NTYPES
  57.         IF     ( TYPES(i) .eq. 'KOIS' ) THEN
  58.           TYPREA(i) = 1
  59.         ELSE IF( TYPES(i) .eq. 'LEBL' ) THEN
  60.           TYPREA(i) = 2
  61.         ELSE
  62.           CALL ERREUR(5)
  63.         ENDIF
  64.  
  65. C       NREAC(i) = Valeur de la position du reactif de la reaction
  66. C                   numero i dans le tableau PHASES
  67. C       Optimisation pour eviter COMPARE_STRING (LENTISSIME)
  68.  
  69.         MOT4a=REACTI(i)
  70.         MOT4c=PRODUI(i)
  71.         DO III=1,NPHASE
  72.           MOT4b=PHASES(III)
  73.           IF(MOT4a .EQ. MOT4b)THEN
  74.             NREAC(i) = III
  75.           ENDIF
  76.           IF(MOT4c .EQ. MOT4b)THEN
  77.             NPROD(i) = III
  78.           ENDIF
  79.         ENDDO
  80.       enddo
  81.  
  82.       SomPha= 1.D0
  83.       do ii = 1, NPHASE
  84.         Xflot   = VAR0(ii)
  85.         SomPha  = SomPha - Xflot
  86.         varf(ii)= Xflot
  87.       enddo
  88.       if(SomPha .GT. XPREC) then
  89.         REAERR(1) = SomPha - 1.D0
  90.         REAERR(2) = SomPha
  91.         CALL ERREUR(1079)
  92.         RETURN
  93.       endif
  94.  
  95.       DT   = tempf -temp0
  96.       ndec = 1
  97.  
  98.       T0   = ture0(1)
  99.       Tf   = turef(1)
  100.  
  101. C     IINTEG=1
  102.  10   continue
  103.       do ii = 1, NPHASE
  104.         VAR0L(ii) = VAR0(ii)
  105.       enddo
  106.       DTL = DT / FLOAT(ndec)
  107.  
  108. C     Boucle sur les sous-pas de temps
  109.       do idec = 1, ndec
  110. C       Determination du temps debut et fin local
  111.         if    ( idec .eq. 1    ) then
  112.           temp0L = temp0
  113.           tempfL = temp0L + DTL
  114.  
  115.         elseif( idec .eq. ndec ) then
  116.           temp0L = tempfL
  117.           tempfL = tempf
  118.  
  119.         else
  120.           temp0L = tempfL
  121.           tempfL = temp0L + DTL
  122.         endif
  123.  
  124.         alpha0 = (temp0L - temp0)/DT
  125.         alphaf = (tempfL - temp0)/DT
  126.         ulpha0 = 1.D0 - alpha0
  127.         ulphaf = 1.D0 - alphaf
  128.  
  129.         ture0L = T0 + alpha0*(Tf - T0)
  130.         turefL = T0 + alphaf*(Tf - T0)
  131.  
  132.         do ini = 1, NPHASE
  133.           Prop(ini)   = VAR0L(ini)
  134.           delta(ini)  = 0.D0
  135.           CRITER(ini) = 0.D0
  136.         enddo
  137.  
  138. C       Boucle de convergence
  139.         iconv = 1
  140.  30     continue
  141.           icont = 1
  142.           rmax  = 0.D0
  143.  
  144. C         Boucle sur les differentes reactions
  145.           do i = 1, NTYPES
  146.             ireac = NREAC(i)
  147.             iprod = NPROD(i)
  148.  
  149.             if( TYPREA(i) .eq. 1 ) then
  150. C             Cas KOISTINEN
  151.               aMS1 = valma0(icont  )
  152.               aMS2 = valmat(icont  )
  153.               aKM1 = valma0(icont+1)
  154.               aKM2 = valmat(icont+1)
  155.               icont= icont + 2
  156.  
  157.               aMS1L = aMS1 + alpha0*(aMS2 - aMS1)
  158.               aKM1L = aKM1 + alpha0*(aKM2 - aKM1)
  159.               aMS2L = aMS1 + alphaf*(aMS2 - aMS1)
  160.               aKM2L = aKM1 + alphaf*(aKM2 - aKM1)
  161.  
  162. C             La reaction a lieu seulement si T°C < MS
  163.               if(ture0L .lt. aMS1L  .OR. turefL .lt. aMS2L) then
  164. C               Tpoint corrige localement en cas de chevauchement avec MS
  165.                 TpoiLo = ( MIN(turefL,aMS1L) - MIN(ture0L,aMS2L) ) / DTL
  166.  
  167. C               critere de convergence local
  168.                 crit1 = ABS(DTL*aKM2L*TpoiLo*theta)
  169.                 CRITER(ireac) = CRITER(ireac) + crit1
  170.                 CRITER(iprod) = CRITER(iprod) + crit1
  171.                 critma = MAX(CRITER(ireac), CRITER(iprod))
  172.  
  173.                 IF( critma .ge. borne ) THEN
  174. C                 On sous-decoupe le pas avant de recommencer entierement
  175.                   ndec = (INT(critma / borne) + 1 ) * ndec
  176.                   GOTO 10
  177.                 ENDIF
  178.  
  179.                 AijL =utheta*MIN(VAR0L(ireac)*aKM1L*TpoiLo,0.D0) +
  180.      &                 theta*MIN( Prop(ireac)*aKM2L*TpoiLo,0.D0)
  181.  
  182.               else
  183.                 AijL = 0.D0
  184.               endif
  185.  
  186.             elseif ( TYPREA(i) .eq. 2 ) then
  187. C             Cas LEBLOND
  188.               PEQ1 = valma0(icont  )
  189.               PEQ2 = valmat(icont  )
  190.               TAU1 = valma0(icont+1)
  191.               TAU2 = valmat(icont+1)
  192.               F1   = valma0(icont+2)
  193.               F2   = valmat(icont+2)
  194.               icont= icont + 3
  195.  
  196.               PEQ1L = PEQ1 + alpha0*(PEQ2 - PEQ1)
  197.               TAU1L = TAU1 + alpha0*(TAU2 - TAU1)
  198.               F1L   = F1   + alpha0*(F2   - F1  )
  199.               PEQ2L = PEQ1 + alphaf*(PEQ2 - PEQ1)
  200.               TAU2L = TAU1 + alphaf*(TAU2 - TAU1)
  201.               F2L   = F1   + alphaf*(F2   - F1  )
  202.  
  203.               UPEQ1L =MIN(MAX(1.D0 - PEQ1L,0.D0),1.D0)
  204.               UPEQ2L =MIN(MAX(1.D0 - PEQ2L,0.D0),1.D0)
  205.  
  206. C             critere de convergence local
  207.               IF (A_EGALE_B(TAU2L,0.D0)) THEN
  208.                 CALL ERREUR(26)
  209.                 RETURN
  210.               ENDIF
  211.               crit1 = ABS(DTL*F2L*utheta*        PEQ2L /TAU2L)
  212.               crit2 = ABS(DTL*F2L*utheta*(1.D0 - PEQ2L)/TAU2L)
  213.               CRITER(ireac) = CRITER(ireac) + crit1
  214.               CRITER(iprod) = CRITER(iprod) + crit2
  215.               critma = MAX(CRITER(ireac), CRITER(iprod))
  216.  
  217.               IF( critma .ge. borne ) THEN
  218. C               On sous-decoupe le pas avant de recommencer entierement
  219.                 ndec = (INT(critma / borne) + 1 ) * ndec
  220.                 GOTO 10
  221.               ENDIF
  222.  
  223. C             TAU Local doit etre superieur a XPETIT :
  224.               if( ABS(TAU1L) .lt. XPETIT .OR.
  225.      &          ABS(TAU2L) .lt. XPETIT ) THEN
  226.                 INTERR(1) = i
  227.                 CALL ERREUR(1081)
  228.                 RETURN
  229.               endif
  230.  
  231.               AijL =utheta*MIN(
  232.      &                (UPEQ1L*VAR0L(iprod)-PEQ1L*VAR0L(ireac))*F1L/TAU1L
  233.      &                         ,0.D0) +
  234.      &               theta*MIN(
  235.      &                (UPEQ2L* Prop(iprod)-PEQ2L* Prop(ireac))*F2L/TAU2L
  236.      &                         ,0.D0)
  237.  
  238.             else
  239.               CALL ERREUR(21)
  240.               RETURN
  241.             endif
  242.  
  243. C           Test pour verifier que le REACTIF calcule sera bien
  244. C           compri entre 0. et 1. Correction de AijL si besoin
  245.             PhTEST = VAR0L(ireac) + delta(ireac) + (DTL*AijL)
  246.             if(PhTEST .lt. 0.D0) then
  247.               AijL= (       XPREC - VAR0L(ireac) - delta(ireac)) / DTL
  248.             else if(PhTEST .gt. 1.D0) then
  249.               AijL= (1.D0 - XPREC - VAR0L(ireac) - delta(ireac)) / DTL
  250.             endif
  251.  
  252. C           Test pour verifier que le PRODUIT calcule sera bien
  253. C           compri entre 0. et 1. Correction de AijL si besoin
  254.             PhTEST = VAR0L(iprod) + delta(iprod) - (DTL*AijL)
  255.             if(PhTEST .lt. 0.D0) then
  256.               AijL= (VAR0L(iprod) - XPREC        + delta(iprod)) / DTL
  257.             else if(PhTEST .gt. 1.D0) then
  258.               AijL= (VAR0L(iprod) + XPREC - 1.D0 + delta(iprod)) / DTL
  259.             endif
  260.  
  261.             xincr        = DTL*AijL
  262.             delta(ireac) = delta(ireac) + xincr
  263.             delta(iprod) = delta(iprod) - xincr
  264.           enddo
  265. C         FIN Boucle sur les differentes reactions
  266.  
  267. C         Actualisation des proportions de phase
  268.           do ii = 1, NPHASE
  269. C           Calcul du residu pour chaque phase (Convergence)
  270.             Pactu= VAR0L(ii) + delta(ii)
  271.             resid= ABS(Prop(ii)-Pactu)
  272.             rmax = MAX(resid,rmax)
  273.  
  274.             Prop(ii)  = Pactu
  275.             delta(ii) = 0.D0
  276.             CRITER(ii)= 0.D0
  277.           enddo
  278.  
  279. C         Test de convergence
  280.           if( iconv .LE. ITMAXI ) THEN
  281. C           Convergence = On sort
  282.             if( rmax .lt. XPREC ) GOTO 20
  283.  
  284.           else
  285. C           Pas Convergence = On sous-decoupe le pas avant de recommencer entierement
  286.             ndec = ndec * 2
  287.             GOTO 10
  288.           endif
  289.  
  290.           iconv = iconv + 1
  291.         GOTO 30
  292.  20     continue
  293. C       SORTIE Boucle de convergence
  294.  
  295. C       La prop en debut du pas N = La prop en fin de pas N-1
  296.         do i = 1, NPHASE
  297.           var0L(i) = Prop(i)
  298.         enddo
  299.       enddo
  300. C     FIN Boucle sur les sous-pas de temps
  301.  
  302. C     Convergence de l'integration
  303. C     XMAX = 0.D0
  304. C     do ii = 1, NPHASE
  305. C       XDIFF    = ABS(Prop(ii) - varf(ii))
  306. C       XMAX     = MAX(XDIFF,XMAX)
  307. C     enddo
  308. C     IF(XMAX .GT. 1.D-4)THEN
  309. C       ndec = ndec * 2
  310. C      PRINT *,'NON-CONVERGENCE',XMAX,IINTEG
  311. C       IINTEG = IINTEG + 1
  312. C       GOTO 10
  313. C     ENDIF
  314.  
  315. C     Ecriture du resultat final
  316.       do ii = 1, NPHASE
  317.         varf(ii) = Prop(ii)
  318.       enddo
  319.  
  320.       END
  321.  
  322.  
  323.  
  324.  

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