Télécharger la plaquette Cast3M


fr en

Télécharger excpha.eso

Retour à la liste

Numérotation des lignes : 

excpha

  1. C EXCPHA    SOURCE    GOUNAND   25/05/05    21:15:04     12259          
  2.       subroutine excpha
  3.  
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8(A-H,O-Z)
  6.  
  7.  
  8. -INC PPARAM
  9. -INC CCOPTIO
  10. -INC SMCHAML
  11. -INC SMCHPOI
  12. -INC SMMODEL
  13. -INC CCREEL
  14. -INC SMCOORD
  15. -INC SMELEME
  16.  
  17.       segment icpr  (NCPR)
  18.       segment icpr1 (NCPR)
  19.       segment icpr2 (NCPR)
  20.  
  21.       SEGMENT SEXCP
  22.         REAL*8  tdebu (node)
  23.         REAL*8  tprop (node)
  24.         REAL*8  qrest (node)
  25.         REAL*8  qtot  (node)
  26.         REAL*8  qsecon(node)
  27.         INTEGER iact  (node)
  28.       ENDSEGMENT
  29.  
  30.       CHARACTER*(LOCOMP) MOCOMP
  31.       LOGICAL     LOG1
  32.  
  33. C Un peu de marge sur la precision des calculs (pb aix)
  34.       XZPREL=XZPREC*1.D2
  35.  
  36. C lecture du modele
  37.       call lirobj('MMODEL  ',mmodel,1,iretou)
  38.       call actobj('MMODEL  ',mmodel,1)
  39.       if(ierr.ne.0) return
  40.  
  41. C lecture du materiau
  42.       call lirobj('MCHAML  ',IPIN,1,iretou)
  43.       call actobj('MCHAML  ',IPIN,1)
  44.       if(ierr.ne.0) return
  45.       CALL REDUAF(IPIN,mmodel,mchelm,0,IR,KER)
  46.       IF(IR   .NE. 1) CALL ERREUR(KER)
  47.       IF(IERR .NE. 0) RETURN
  48.  
  49. C lecture des proportions de phase
  50.       call lirobj('MCHAML  ',IPIN,1,iretou)
  51.       call actobj('MCHAML  ',IPIN,1)
  52.       if(ierr.ne.0) return
  53.       CALL REDUAF(IPIN,mmodel,mchel1,0,IR,KER)
  54.       IF(IR   .NE. 1) CALL ERREUR(KER)
  55.       IF(IERR .NE. 0) RETURN
  56.  
  57. C lecture des chaleur latentes reduies aux noeuds
  58.       call lirobj('MCHAML  ',IPIN,1,iretou)
  59.       call actobj('MCHAML  ',IPIN,1)
  60.       if(ierr.ne.0) return
  61.       CALL REDUAF(IPIN,mmodel,mchel2,0,IR,KER)
  62.       IF(IR   .NE. 1) CALL ERREUR(KER)
  63.       IF(IERR .NE. 0) RETURN
  64.  
  65. C lecture du chpoint de temperatures initiales
  66.       call lirobj('CHPOINT ',mchpoi,1,iretou)
  67.       CALL ACTOBJ('CHPOINT ',mchpoi,1)
  68.       if(ierr.ne.0) return
  69.  
  70. C lecture du chpoint d'increment de temperatures proposées
  71.       call lirobj('CHPOINT ',mchpo1,1,iretou)
  72.       CALL ACTOBJ('CHPOINT ',mchpo1,1)
  73.       if(ierr.ne.0) return
  74.  
  75. C fin des lectures
  76.  
  77. C Un petit test de coherence mais sans consequence :
  78.       ifopo3=mchpoi.ifopoi
  79.       if (mchpoi.ifopoi .ne. mchpo1.ifopoi) then
  80.         moterr(1:8)='CHPOINT'
  81.         interr(1)=mchpoi.ifopoi
  82.         interr(2)=mchpo1.ifopoi
  83.         interr(3)=ifour
  84. c-dbg      write(ioimp,*) '1132 excpha',mchpoi,mchpo1
  85.         call erreur(1132)
  86.         ifopo3=ifour
  87.       end if
  88.  
  89.       XUN  =1.D0
  90.       NCPR=nbpts
  91.       segini,icpr,icpr1,icpr2
  92.  
  93. C algorithmes on boucle sur les sous-modeles
  94. C puis si la rigidite avait ete prise en compte on teste le sens
  95. C et la valeur de la chaleur latente proposee
  96. C sinon on teste la non traverse de la temperature de changement de
  97. C phase et au besoin on active le blocage (temperature imposee et
  98. C flux de chaleur latente limite)
  99.  
  100. * +--------------------------------------------------------------------+
  101. * |    Preparation de QSECOND :                                        |
  102. * |     CHPOINT de 'LX'                                                |
  103. * |     Remplace les chaleurs latentes dans RESO : UNILATER            |
  104. * +--------------------------------------------------------------------+
  105.       N1    = mmodel.kmodel(/1)
  106.       ib    = 0
  107.       do 200 meri = 1,N1
  108.         imodel = mmodel.kmodel(meri)
  109.  
  110.         nfor=imodel.formod(/2)
  111.         call place(imodel.formod,nfor,iplac,'CHANGEMENT_PHASE')
  112.         if  (iplac .eq. 0) goto 200
  113.         if(tymode(2) .ne. 'MAILLAGE')then
  114.           call erreur(5)
  115.         endif
  116.         ipt2=ivamod(2)
  117.         do 201 nel=1,ipt2.num(/2)
  118. C         Noeud 1 : 'LX'
  119. C         Noeud 2 : 'Inconnues classiques'
  120.           ia = ipt2.num(1,nel)
  121.           if(icpr2(ia).eq.0) then
  122.             ib       = ib+1
  123.             icpr2(ia)= ib
  124.           endif
  125.   201   continue
  126.   200 continue
  127.       node = ib
  128.  
  129. * +--------------------------------------------------------------------+
  130. * |    NOEUDS du CHPOINT debu et proposé                               |
  131. * +--------------------------------------------------------------------+
  132.       segini,icpr
  133.       ib = 0
  134.       do 80 I=1,mchpoi.ipchp(/1)
  135.         msoupo=mchpoi.ipchp(i)
  136. C       Le 'LX' du champ debut ne nous interesse pas
  137.         if (msoupo.nocomp(1) .EQ. 'LX  ') goto 80
  138.         ipt1  =msoupo.igeoc
  139.         mpoval=msoupo.ipoval
  140.         do 81 j=1,ipt1.num(/2)
  141.           ia = ipt1.num(1,j)
  142.           if(icpr(ia).eq.0) then
  143.             ib = ib +1
  144.             icpr(ia)=ib
  145.           endif
  146.    81   continue
  147.    80 continue
  148.       node=max(node,ib)
  149.  
  150.       ib = 0
  151.       segini,icpr1
  152.       do 90 I=1,mchpo1.ipchp(/1)
  153.         msoupo=mchpo1.ipchp(i)
  154.         ipt1  =msoupo.igeoc
  155.         do 91 j=1,ipt1.num(/2)
  156.           ia = ipt1.num(1,j)
  157.           if(icpr1(ia).eq.0) then
  158.             ib       = ib +1
  159.             icpr1(ia)= ib
  160.           endif
  161.    91   continue
  162.    90 continue
  163.       node=max(node,ib)
  164.  
  165.       segini,SEXCP
  166.  
  167. * +--------------------------------------------------------------------+
  168. * |    Boucle sur les SOUS-ZONES du MODELE                             |
  169. * +--------------------------------------------------------------------+
  170.       do 100 mo=1,N1
  171.         imodel = mmodel.kmodel(mo)
  172.         meleme = imodel.imamod
  173.  
  174.         nfor=imodel.formod(/2)
  175.         call place(imodel.formod,nfor,iplac,'CHANGEMENT_PHASE')
  176.         if  (iplac .eq. 0) goto 100
  177.  
  178. *       +------------------------------------------------------------------------+
  179. *       |    Recherche des composantes PRIMALE + 'LX' du CHPOINT debu et propose |
  180. *       +------------------------------------------------------------------------+
  181.         nomid=imodel.lnomid(1)
  182.         MOCOMP=nomid.lesobl(1)
  183.  
  184. C       remise a zero qui va bien
  185.         IF(mo .GT. 1)THEN
  186.           CALL ZERO(tdebu(1),node,1)
  187.           CALL ZERO(tprop(1),node,1)
  188.         ENDIF
  189.  
  190.         do 82 I=1,mchpoi.ipchp(/1)
  191.           msoupo=mchpoi.ipchp(i)
  192. C         Le 'LX' du champ debut ne nous interesse pas
  193.           if (msoupo.nocomp(1) .EQ. 'LX  ') goto 82
  194.           call place(msoupo.nocomp,msoupo.nocomp(/2),iplac,MOCOMP)
  195.           if  (iplac .eq. 0) then
  196.             moterr = MOCOMP
  197.             call erreur(181)
  198.             return
  199.           endif
  200.           ipt1  =msoupo.igeoc
  201.           mpoval=msoupo.ipoval
  202.           do 85 k=1,ipt1.num(/2)
  203.             ia             = ipt1.num(1,k)
  204.             tdebu(icpr(ia))= mpoval.vpocha(k,iplac)
  205.    85     continue
  206.    82   continue
  207.  
  208.         do 92 I=1,mchpo1.ipchp(/1)
  209.           msoup1= mchpo1.ipchp(i)
  210.           iplac = 1
  211.           if (msoup1.nocomp(1) .EQ. 'LX  ') goto 94
  212.           call place(msoup1.nocomp,msoup1.nocomp(/2),iplac,MOCOMP)
  213.           if  (iplac .eq. 0) then
  214.             moterr = MOCOMP
  215.             call erreur(181)
  216.             return
  217.           endif
  218.    94     continue
  219.           ipt1  = msoup1.igeoc
  220.           mpova1= msoup1.ipoval
  221.           do 95 k=1,ipt1.num(/2)
  222.             ia       = ipt1.num(1,k)
  223.             if(icpr(ia) .EQ. 0)then
  224. C             Cas des 'LX'
  225.               tprop(icpr1(ia))= mpova1.vpocha(k,iplac)
  226.             else
  227.               tprop(icpr1(ia))= mpova1.vpocha(k,iplac) + tdebu(icpr(ia))
  228.             endif
  229.    95     continue
  230.    92   continue
  231.  
  232. *       +--------------------------------------------------------------+
  233. *       |    Recherche du MCHAML correspondant a 'PRIM'                |
  234. *       +--------------------------------------------------------------+
  235.         do 51 mchm=1,imache(/1)
  236.           if(conche(mchm) .ne. conmod) goto 51
  237.           if(imache(mchm) .eq. meleme) then
  238.             mchaml=ichaml(mchm)
  239.             goto 52
  240.           endif
  241.    51   continue
  242.         call erreur(472)
  243.         return
  244.    52   continue
  245.  
  246.         do 56 n2=1,nomche(/2)
  247.           if(nomche(n2).eq.'PRIM') then
  248.             if (typche(n2).ne.'REAL*8') then
  249.               moterr(1:16)  = typche(n2)
  250.               moterr(17:20) ='PPHA'
  251.               moterr(21:36) = mchelm.titche
  252.               call erreur(552)
  253.               RETURN
  254.             endif
  255.             melval=ielval(n2)
  256.             goto 57
  257.           endif
  258.    56   continue
  259.         moterr(1:8) = 'PRIM'
  260.         call erreur(677)
  261.         return
  262.  
  263.    57   continue
  264. C       On suppose 'PRIM' constant par sous-modele : verification
  265.         IF(velche(/1) .NE. 1 .AND. velche(/2).NE. 1)THEN
  266.           MOTERR(1:4)='PRIM'
  267.           L1          =mchelm.TITCHE(/1)
  268.           L2          =MIN(20,L1)
  269.           MOTERR(5:L2)=mchelm.TITCHE
  270.           CALL ERREUR(1065)
  271.           RETURN
  272.         ENDIF
  273.  
  274.         T_chph = melval.velche(1,1)
  275.  
  276. *       +--------------------------------------------------------------+
  277. *       |    Recherche du MCHAML correspondant a 'PPHA'                |
  278. *       +--------------------------------------------------------------+
  279.         do 61 mchm=1,mchel1.imache(/1)
  280.           if(mchel1.conche(mchm) .ne. conmod) goto 61
  281.           if(mchel1.imache(mchm) .eq. meleme) then
  282.             mcham1=mchel1.ichaml(mchm)
  283.             goto 62
  284.           endif
  285.    61   continue
  286.         call erreur(472)
  287.         return
  288.  
  289.    62   continue
  290.         do 66 n2=1,mcham1.nomche(/2)
  291.           if(mcham1.nomche(n2).eq.'PPHA') then
  292.             if (typche(n2).ne.'REAL*8') then
  293.               moterr(1:16)  = mcham1.typche(n2)
  294.               moterr(17:20) ='PPHA'
  295.               moterr(21:36) = mchel1.titche
  296.               call erreur(552)
  297.               RETURN
  298.             endif
  299.             melva1=mcham1.ielval(n2)
  300.             goto 67
  301.           endif
  302.    66   continue
  303.         moterr(1:8) = 'PPHA'
  304.         call erreur(677)
  305.         return
  306.    67   continue
  307.  
  308. *       +---------------------------------------------------------------------+
  309. *       |    Recherche du MCHAML correspondant au flux equivalent (aux NOEUDS)|
  310. *       +---------------------------------------------------------------------+
  311.         do 71 mchm=1,mchel2.imache(/1)
  312.           if(mchel2.conche(mchm) .ne. conmod) goto 71
  313.           if(mchel2. imache(mchm) .eq. meleme) then
  314.             mcham2=mchel2.ichaml(mchm)
  315.             goto 72
  316.           endif
  317.    71   continue
  318.         call erreur(472)
  319.         return
  320.    72   continue
  321.  
  322.         nomid = lnomid(2)
  323.         do 76 n2=1,mcham2.nomche(/2)
  324.           if(mcham2.nomche(n2).eq.nomid.lesobl(1)) then
  325.             if (mcham2.typche(n2).ne.'REAL*8') then
  326.               moterr(1:16)  = mcham2.typche(n2)
  327.               moterr(17:20) = nomid.lesobl(1)
  328.               moterr(21:36) = mchel2.titche
  329.               call erreur(552)
  330.               RETURN
  331.             endif
  332.             melva2=mcham2.ielval(n2)
  333.             goto 77
  334.           endif
  335.    76   continue
  336.         moterr(1:8) =nomid.lesobl(1)
  337.         call erreur(677)
  338.         return
  339.    77   continue
  340.  
  341. C       remise a zero qui va bien
  342.         IF(mo .GT. 1)THEN
  343.           CALL ZERO(qrest(1),node,1)
  344.           CALL ZERO(qtot(1) ,node,1)
  345.         ENDIF
  346.  
  347. C       Fabrication de qtot : chaleur latente totale   cumulee par point
  348. C       Fabrication de qrest: chaleur latente restante cumulee par point
  349.         noMAX1=melva1.velche(/1)
  350.         meMAX1=melva1.velche(/2)
  351.         noMAX2=melva2.velche(/1)
  352.         meMAX2=melva2.velche(/2)
  353.         do 110 mel=1,meleme.num(/2)
  354.           do 111 nod = 1,meleme.num(/1)
  355.             ia        = meleme.num(nod,mel)
  356.             ib        = icpr(ia)
  357.             pp        = melva1.velche(MIN(nod,noMAX1),MIN(mel,meMAX1))
  358.             umpp      = MAX(MIN(XUN-pp,XUN),XZERO)
  359.             ql        = melva2.velche(min(nod,noMAX2),min(mel,meMAX2))
  360.             IF(ql .EQ. XZERO)THEN
  361. C             Protection pour les chaleurs latentes strictement egales a 0.D0 !!
  362.               ql      = XPETIT/XZPREL
  363.             ENDIF
  364.             qr        = ql*umpp
  365.  
  366.             qrest(ib) = qrest(ib) + qr
  367.             qtot (ib) = qtot (ib) + ql
  368.   111     continue
  369.   110   continue
  370.  
  371.         ipt2   = imodel.ivamod(2)
  372.         do 120 iel=1,ipt2.num(/2)
  373. C         Noeud 1 : 'LX'
  374. C         Noeud 2 : 'Inconnues classiques'
  375.           No_LX = ipt2.num(1,iel)
  376.           No_INC= ipt2.num(2,iel)
  377.  
  378.           iex   = icpr1(No_LX)
  379.           iri   = icpr2(No_LX)
  380.           ideb  = icpr(No_INC)
  381.           ipro  = icpr1(No_INC)
  382.  
  383.           tdeb  = tdebu(ideb)
  384.           IF(ipro .NE. 0)THEN
  385.             tpro  = tprop(ipro)
  386.           ELSE
  387. C           Cas a la 1ere iteration
  388.             tpro  = tdeb
  389.           ENDIF
  390.  
  391.           q_tot = qtot(ideb)
  392.           q_res = qrest(ideb)
  393.           q_dej = q_tot-q_res
  394.           q_pre = q_tot*XZPREL
  395.  
  396. C         Pour eviter les soucis de denormalisation si T_chph ~ XZERO
  397.           LOG1  =(tdeb.EQ.T_chph) .OR.
  398.      &      (ABS(T_chph-tdeb).LT.MAX(XZPREL*ABS(TDEB+T_chph),XPETIT))
  399.           ibloc = 0
  400.  
  401.           if(iex .ne. 0) then
  402. *           +----------------------------------------------------------+
  403. *           |    Le 'LX' existe dans T proposée sur le noeud No_LX     |
  404. *           +----------------------------------------------------------+
  405.             rea_LX= tprop(iex)
  406.             if(LOG1) then
  407. C             Temperature initiale = T_chph
  408.               if    (rea_LX.GT.XZERO .AND. q_res.GT.q_pre) then
  409. C               la temperature cherche a monter
  410. C               PRINT*,'-E:',No_LX,tdeb,T_chph,tpro
  411.                 ibloc =-1
  412.               elseif(rea_LX.LT.XZERO .AND. q_dej.GT.q_pre) then
  413. C               la temperature cherche a descendre
  414. C               PRINT*,'-F:',No_LX,tdeb,T_chph,tpro
  415.                 ibloc = 1
  416.               endif
  417.  
  418.             elseif(rea_LX.GT.XZERO) then
  419. C             La temperature cherche a monter
  420.               if(tdeb.LT.T_chph) then
  421. C               PRINT*,'-G:',No_LX,tdeb,tpro,q_tot,q_res,q_pre
  422.                 ibloc =-1
  423.               endif
  424.  
  425.             elseif(rea_LX.LT.XZERO) then
  426. C             la temperature cherche a descendre
  427.               if(tdeb .GT. T_chph) then
  428. C               PRINT*,'-G:',No_LX,tdeb,tpro,q_tot,q_dej,q_pre
  429.                 ibloc = 1
  430.               endif
  431.             endif
  432.  
  433.           else
  434. *           +-------------------------------------------------------------+
  435. *           |    Le 'LX' n''existe pas dans T proposée sur le noeud No_LX |
  436. *           +-------------------------------------------------------------+
  437.             if(LOG1) then
  438. C             Temperature initiale = T_chph
  439.               tsens = tpro - T_chph
  440.               if    (tsens.GT.XZERO .AND. q_res.GT.q_pre) then
  441. C               La temperature cherche a monter
  442.                 ibloc =-1
  443. C               PRINT*,'-AA:',No_LX,tdeb,T_chph,tpro,q_tot,q_res,q_res
  444.               elseif(tsens.LT.XZERO .AND. q_dej.GT.q_pre) then
  445. C               la temperature cherche a descendre
  446.                 ibloc = 1
  447. C               PRINT*,'-BB:',No_LX,tdeb,T_chph,tpro,q_tot,q_res,q_res
  448.               endif
  449.  
  450.             elseif((T_chph-tdeb) * (T_chph-tpro) .LT. XZERO) then
  451. C             tdeb et tprop de chaque cote de T_chph : blocage actif
  452.               if    (tdeb .LT. T_chph) then
  453.                 ibloc =-1
  454. C               PRINT*,'-C:',No_LX,(T_chph-tdeb),(T_chph-tpro)
  455.               elseif(tdeb .GT. T_chph) then
  456.                 ibloc = 1
  457. C               PRINT*,'-D:',No_LX,tdeb,T_chph,tpro
  458.               else
  459.                 CALL ERREUR(5)
  460.               endif
  461.  
  462.             elseif(tdeb.LT.T_chph .AND. q_dej.GT.q_pre)then
  463. C             apparemment il reste de la phase
  464.               ibloc = 1
  465. C             PRINT*,'-A:',No_LX,tdeb,T_chph,tpro,q_tot,q_res,q_res
  466.  
  467.             elseif(tdeb.GT.T_chph .AND. q_res.GT.q_pre)then
  468. C             apparemment on n'a pas tout consomme
  469.               ibloc =-1
  470. C             PRINT*,'-B:',No_LX,tdeb,T_chph,tpro,q_tot,q_res,q_res
  471.             endif
  472.           endif
  473.  
  474. C         Les conditions mises sont unilaterales, RESO s'en sort en ne
  475. C         gardant que la premiere 'excitee' dans le sens en question
  476.           IF    (ibloc.EQ. 1)THEN
  477.             iact(iri)     = No_LX
  478.             qsecon(iri)   = q_dej
  479. C           PRINT*,'iact(iri)',No_LX,ibloc,iact(iri),qsecon(iri)
  480.           ELSEIF(ibloc.EQ.-1)THEN
  481.             iact(iri)     = No_LX
  482.             qsecon(iri)   =-q_res
  483. C           PRINT*,'iact(iri)',No_LX,ibloc,iact(iri),qsecon(iri)
  484.           ELSE
  485.             iact(iri)     = 0
  486.             qsecon(iri)   = XZERO
  487.           ENDIF
  488.  
  489. C         fin de la boucle sur les elements du MAILLAGE de la rigidite
  490.   120   continue
  491. C       fin de la boucle sur les sous-modeles
  492.   100 continue
  493.  
  494. *    +--------------------------------------------------------------+
  495. *    | Creation du resultat : CHPOINT de 'LX'                       |
  496. *    |   on dispose maintenant du tableau qsecon qui donne la       |
  497. *    |   chaleur latente a mettre en 'LX' sur les noeuds de iact    |
  498. *    +--------------------------------------------------------------+
  499.       nbelem = 0
  500.       do 300 ia=1,node
  501.         if(iact(ia) .ne. 0) nbelem = nbelem + 1
  502.   300 continue
  503.  
  504.       nbref  = 0
  505.       nbsous = 0
  506.       nat    = 1
  507.       if(nbelem .eq. 0) then
  508.         nbnn   = 0
  509.         nsoupo = 0
  510.         segini,mchpo3
  511.  
  512.       else
  513.         nbnn   = 1
  514.         segini,ipt4
  515.         ipt4.itypel = 1
  516.  
  517.         nsoupo = 1
  518.         nc     = 1
  519.         n      = nbelem
  520.         segini,mchpo3,msoupo,mpoval
  521.         mchpo3.ipchp(1)    = msoupo
  522.         msoupo.nocomp(1)   ='LX'
  523.         msoupo.igeoc       = ipt4
  524.         msoupo.ipoval      = mpoval
  525.  
  526.         ipo=0
  527.         do 301,ia=1,node
  528.           iact1 = iact(ia)
  529.           if (iact1 .eq. 0) goto 301
  530.           ipo             = ipo + 1
  531.           ipt4.num(1,ipo) = iact1
  532.           vpocha(ipo,1)   = qsecon(icpr2(iact1))
  533.  301    continue
  534.       endif
  535.  
  536.       mchpo3.mochde(1:72)='chpoint cree par EXCP'
  537.       mchpo3.mtypoi      ='contmult'
  538.       mchpo3.jattri(1)   = 2
  539.       mchpo3.ifopoi      = ifopo3
  540.  
  541.       segsup icpr,icpr1,icpr2,SEXCP
  542.  
  543.       call actobj('CHPOINT ',mchpo3,1)
  544.       call ecrobj('CHPOINT ',mchpo3)
  545.  
  546. c      return
  547.       end
  548.  
  549.  

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