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

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