Télécharger la plaquette Cast3M


fr en

Télécharger fatig2.eso

Retour à la liste

Numérotation des lignes : 

fatig2

  1. C FATIG2    SOURCE    PV090527  25/01/07    14:42:38     12115          
  2.       SUBROUTINE FATIG2(ITCONT,ITTEMP,IPMODE,IPMSTA,ICF1,xre1,xre2,
  3.      &ICLE,NCLE,CLE,ZECRIT,ICHOUT)
  4.       IMPLICIT INTEGER(I-N)
  5.       IMPLICIT REAL*8 (A-H,O-Z)
  6.  
  7. -INC PPARAM
  8. -INC CCOPTIO
  9. -INC SMCOORD
  10. -INC SMMODEL
  11. -INC SMCHAML
  12. -INC SMELEME
  13. -INC SMEVOLL
  14. -INC SMLREEL
  15.       PARAMETER (MCRIT=5)
  16.       SEGMENT,MLCARF
  17.         integer lcarfa(2*MCRIT,N1)
  18.       ENDSEGMENT
  19.       SEGMENT,MCYSIG
  20.         integer lcysig(nbrobl,ncycl)
  21.         real*8 sigcyc(nbrobl,ncycl),pcyc(ncycl)
  22.       ENDSEGMENT
  23.       SEGMENT,MRECYC
  24.          real*8 ycyc(ncycl)
  25.       ENDSEGMENT
  26.       SEGMENT,MDEVCY
  27.          real*8 sdcyc(nbrobl-1,ncycl)
  28.       ENDSEGMENT
  29.       SEGMENT,MDEVSI
  30.          real*8 sd(nbrobl)
  31.       ENDSEGMENT
  32.       LOGICAL LOG0,LOG1,dcarf1,dcarf2,d_cle,lcas1
  33.       CHARACTER*4 COFA(2*MCRIT),CLE(NCLE)
  34.       real*8 cofa1(NCLE-1),cofa2(NCLE-1)
  35.       DATA COFA/'ADVK','BDVK','APAP','BPAP','ASIN','BSIN','ACRO','BCRO',
  36.      &'A_DC','B_DC'/
  37.  
  38.  
  39.  
  40.       SQ2 = dsqrt(2.d0)
  41.       SQ3S2 = dsqrt(1.5D0)
  42.       SQ3 = dsqrt(3.d0)
  43.  
  44.       if (ipmsta.gt.0) then
  45.        lcas1 = .false.
  46.       else
  47.        lcas1 = .true.
  48.       endif
  49.  
  50.       mmodel = ipmode
  51.       segact mmodel
  52.       n1 = kmodel(/1)
  53.       n2 = 0
  54.       l1 = 16
  55.       n3 = 6
  56.       segini mchel2
  57.       mchel2.titche(1:8) = 'FATIGUE '
  58.  
  59.       if (ICF1.gt.0) then
  60.         mchel1= ICF1
  61.         segact mchel1
  62. *       mchel2.titche(9:16) = mchel1.titche(9:16)
  63.       endif
  64.       segini mlcarf
  65.  
  66.       if(icle.ge.2.and.icle.le.6) then 
  67.          n=0
  68.          segini mevoll
  69.          IEVTEX = 'EVOLUTION VIDE'
  70.          mevnul = mevoll
  71.          segdes mevoll
  72.       endif
  73.  
  74.  
  75.  
  76.       do ik = 1,n1
  77.         imodel = kmodel(ik)
  78.         segact imodel
  79.         mchel2.conche(ik) = conmod
  80.         mchel2.imache(ik) = imamod
  81.         mchel2.ifoche = ifour
  82.         mchel2.infche(ik,4) = infmod(7)
  83.         mchel2.infche(ik,6) = 5
  84.  
  85.         if(ICF1.gt.0) then
  86.         do ic = 1,mchel1.imache(/1)
  87.          if (mchel1.imache(ic).eq.imamod) then
  88.            mchaml = mchel1.ichaml(ic) 
  89.            segact mchaml
  90.            n2 = nomche(/2)     
  91.            do inom = 1,n2
  92. * controler les noms des caracteristiques du critere             
  93.              melval = ielval(inom)
  94.              if(icle.le.6) then
  95.                do jcr = 1,mcrit
  96.               if(nomche(inom)(1:4).eq.cofa(2*jcr-1)(1:4)) then
  97.                              segact melval
  98.                       lcarfa(2*jcr-1,ik) = melval
  99.               endif
  100.               if(nomche(inom)(1:4).eq.cofa(2*jcr)(1:4)) then
  101.                              segact melval
  102.                       lcarfa(2*jcr,ik) = melval
  103.               endif
  104.                enddo
  105.              endif             
  106.  
  107.            enddo
  108.            segdes mchaml
  109.          endif 
  110.         enddo
  111. * verification        
  112.           d_cle = .true.
  113.           do jcr = 1,mcrit
  114.             dcarf1 = .false.
  115.             dcarf2 = .false.
  116.             if (lcarfa(2*jcr-1,ik).gt.0) dcarf1 = .true. 
  117.             if (lcarfa(2*jcr,ik).gt.0)   dcarf2 = .true.
  118.             if (icle.eq.1) then
  119.              d_cle = dcarf1 .and. dcarf2 .and. d_cle
  120.             else if (icle.ge.2.and.icle.le.6.and.jcr.eq.icle-1) then
  121.               d_cle = dcarf1 .and. dcarf2
  122.             endif
  123.           enddo        
  124.  
  125.           if (.not.d_cle) then
  126.             call erreur(472)
  127.             return
  128.           endif
  129.         endif
  130.  
  131. * sorties
  132.          if(icle.eq.1) then
  133.            n2 = ncle - 1
  134.          elseif(icle.ge.2.and.icle.le.6) then
  135.            n2 = 2
  136.          else 
  137.            n2 = 1
  138.          endif
  139.            segini mchaml
  140.            mchel2.ichaml(ik) = mchaml
  141.          meleme = imamod
  142.          segact meleme
  143.          nbelem = num(/2)
  144.          nbgs = infele(4)
  145.            n1ptel= nbgs
  146.            n1el = nbelem
  147.            n2ptel = 0
  148.            n2el = 0
  149.          if(icle.eq.1) then
  150.            do je = 1,n2
  151.             segini melval
  152.             ielval(je) = melval
  153.             nomche(je) = cle(je+1)
  154.             typche(je) = 'REAL*8'
  155.            enddo
  156.          elseif(icle.gt.1) then
  157.           segini melval
  158.           ielval(1) = melval
  159.           nomche(1) = cle(icle)
  160.            typche(1) = 'REAL*8'
  161.          endif
  162.  
  163.          if(icle.ge.2.and.icle.le.6) then 
  164.            n2ptel = nbgs
  165.            n2el = nbelem
  166.            n1ptel = 0
  167.            n1el = 0
  168.            segini melval
  169.          ielval(2) = melval
  170.          nomche(2) = 'PTAU'
  171.          typche(2) = 'POINTEUREVOLUTIO'
  172.               segini kevoll
  173.               ielche(1,1) = kevoll
  174.               kevdvk = kevoll
  175.               jg = 2
  176.               segini mlreel
  177.               iprogx = mlreel
  178.               segini mlreel
  179.               iprogy = mlreel
  180.               NUMEVX = 4
  181.               NUMEVY='REEL'
  182.               NOMEVX = 'P'
  183.               NOMEVY = 'TAU'
  184.               TYPX = 'LISTREEL'
  185.               TYPY = 'LISTREEL'
  186.          endif
  187.  
  188.       enddo
  189.  
  190.       if (lcas1) then
  191.         mtable = ITCONT
  192.         mtab1 = ittemp
  193.       else
  194.         mmode1 = ipmsta
  195.         segact mmode1 
  196.         n1sta = mmode1.kmodel(/1)
  197.         mchelm = itcont
  198.         segact mchelm
  199.         if (imache(/1).ne.n1sta) then
  200. *          write(6,*) 'correspondance MCHELM et MMODEL stationnaires ?'
  201.           call erreur(21)
  202.           return
  203.         endif
  204.       endif 
  205.  
  206.  
  207.       i0 = 0 
  208.       X0 = 0.D0
  209.       LOG0 = .TRUE.
  210.       ip0 = 0
  211.       I1 = 0
  212.       x1 = 0.d0
  213.       LOG1 = .TRUE.
  214.       IP1 = 0
  215.  
  216.       if (lcas1) then
  217.       CALL DIMEN7 (ittemp,ntemps)
  218.       do jr =1,2
  219.         if(jr.eq.1) xreu = xre1
  220.         if(jr.eq.2) xreu = xre2
  221. * presuppose indice 0 t=0.
  222.       if(xreu.eq.0.d0) then
  223.        intc = 0
  224.       elseif(xreu.gt.0.d0) then
  225.          xd1 = xreu
  226.        do ind1 = 1,ntemps
  227.            I0 = ind1 - 1
  228.            CALL ACCTAB(ITTEMP,'ENTIER',I0,X0,' ',LOG0,IP0,
  229.      &                       'FLOTTANT',I1,X1,' ',LOG1,IP1)
  230.           if (ierr.ne.0) return
  231.          xu = xreu - x1
  232.          if (xu.gt.0.d0.and.xu.le.xd1) then
  233.            xd1 = xu
  234.          else if (dabs(xu).le.xd1) then
  235.            goto 14
  236.          else
  237.            I0 = I0 - 1
  238.            goto 14
  239.          endif
  240.        enddo        
  241.  14    continue 
  242.        intc = I0
  243.       endif
  244.  
  245.         if(jr.eq.1) i0temd = intc
  246.         if(jr.eq.2) then
  247.          if(xre2.gt.0) then
  248.            i0temf = intc
  249.          else
  250.            i0temf = ntemps -1
  251.          endif
  252.          ncycl = i0temf - i0temd + 1
  253.         endif
  254.       enddo
  255.  
  256.       else
  257.         ncycl = int(n1sta/n1)
  258.       endif
  259.  
  260.       DO ik = 1,n1
  261.          isk = 0
  262.          imodel = kmodel(ik)
  263.          meleme = imamod
  264.  
  265. * sorties 
  266.          mcham2 = mchel2.ichaml(ik)
  267.          nomid = lnomid(4)
  268.          segact nomid
  269.          nbrobl = lesobl(/2)
  270.          segini mdevsi
  271.  
  272.          segini mcysig
  273.          segini MDEVCY
  274.          segini mrecyc
  275.  
  276.         if (lcas1) then
  277.           I0 = i0temd - 1
  278.         else
  279.           IS0 = 1 
  280.         endif
  281.  
  282.         ktem = 0
  283.         DO jcyc = 1,ncycl
  284.  
  285.           ktem = ktem + 1
  286.  
  287.          if (lcas1) then
  288.            I0 = I0 + 1
  289.            if (I0.gt.i0temf) then
  290. *             write(6,*) 'hepepep',I0,I0temf,i0temd,ncycl
  291.              call erreur(21)
  292.              return
  293.            endif
  294.            CALL ACCTAB(ITCONT,'ENTIER',I0,X0,' ',LOG0,IP0,
  295.      &                       'MCHAML ',I1,X1,' ',LOG1,IP1)
  296.            MCHELM = ip1
  297.            SEGACT MCHELM 
  298.  
  299.           do is = 1,imache(/1)
  300.            if(imamod.eq.imache(is)) isk = is
  301.           enddo
  302.           if (isk.eq.0) then
  303.            call erreur(472)
  304.            return
  305.           endif
  306.           mchaml = ichaml(isk) 
  307.  
  308.          else
  309.            IS1 = 0
  310.            do 24 isou = IS0,n1sta
  311.              imode1 = mmode1.kmodel(isou)
  312. c* test rustique, on peut utiliser objmod
  313.              if (imode1.nefmod.ne.nefmod.or.imode1.imatee.ne.imatee
  314.      &.or.imode1.cmatee.ne.cmatee) goto 24
  315.              ipt1 = imode1.imamod
  316.              if (ipt1.itypel.ne.itypel) goto 24
  317.              IS1 = isou
  318.              goto 25
  319.   24       continue
  320.   25       continue
  321.            if (IS1.eq.0.and.ktem.ne.ncycl) then
  322. *             write(6,*) 'pas de tranche ',ktem,' stationnaire'
  323.              call erreur(21)
  324.              return
  325.            endif
  326.            IS0 = IS1 + 1
  327.            do im = IS1,n1sta
  328.              if (imache(im).eq.imode1.imamod) goto 27
  329.            enddo
  330. *           write(6,*) 'pas de mchaml stationnaire'
  331.            call erreur(21)
  332.            return
  333.   27       continue
  334.            mchaml = ichaml(im)
  335.            if (conche(im).ne.imode1.conmod) then
  336. *             write(6,*) 'perplexe ??' 
  337.              call erreur(21)
  338.              return
  339.            endif
  340.  
  341.          endif
  342.  
  343.          segact mchaml
  344. * controle des composantes de contraintes
  345.          n2 = nomche(/2)
  346. * on travaille a priori avec les composantes obligatoires
  347.          do iobl = 1, nbrobl
  348.           do imch = 1,nomche(/2)
  349.             if(lesobl(iobl).eq.nomche(imch)) then
  350.               lcysig(iobl,ktem) = ielval(imch)
  351.               melval = ielval(imch)
  352.               segact melval
  353.             endif
  354.           enddo
  355.          enddo
  356.  
  357.          segdes mchaml
  358.         ENDDO
  359.  
  360.          meleme = imamod
  361.          nbelem = num(/2)
  362.          nbgs = infele(4)
  363.          nstrs = infele(16)
  364.          mfr = infele(13)
  365.          DO ib = 1,nbelem
  366.           do igau = 1,nbgs
  367.  
  368. * caracteristiques critere
  369.             IF(ib.eq.1.and.igau.eq.1) THEN
  370. * kich : d un point de vue pratique on attend des constantes
  371.             if(icle.eq.1) then
  372.               do jcr = 1,mcrit
  373.                melva1 = lcarfa(2*jcr-1,ik)
  374.                IGMN=MIN(IGAU,melva1.VELCHE(/1))
  375.                IBMN=MIN(IB  ,melva1.VELCHE(/2))
  376.                cofa1(jcr) = melva1.velche(igmn,ibmn)*(-1)
  377.  
  378.                melva2 = lcarfa(2*jcr,ik)
  379.                IGMN=MIN(IGAU,melva2.VELCHE(/1))
  380.                IBMN=MIN(IB  ,melva2.VELCHE(/2))
  381.                cofa2(jcr) = melva2.velche(igmn,ibmn)
  382.               enddo
  383.  
  384.             elseif(icle.ge.2.and.icle.le.6) then
  385.                melva1 = lcarfa(2*icle-3,ik)
  386.                IGMN=MIN(IGAU,melva1.VELCHE(/1))
  387.                IBMN=MIN(IB  ,melva1.VELCHE(/2))
  388.                cofa1(icle-1) = melva1.velche(igmn,ibmn)*(-1)
  389.  
  390.                melva2 = lcarfa(2*icle-2,ik)
  391.                IGMN=MIN(IGAU,melva2.VELCHE(/1))
  392.                IBMN=MIN(IB  ,melva2.VELCHE(/2))
  393.                cofa2(icle-1) = melva2.velche(igmn,ibmn)
  394.             endif
  395.              cofa1(6) = 0.d0
  396.              cofa2(6) = 0.d0
  397.             ENDIF
  398.  
  399. * trajet de chargement
  400.             DO icyc = 1,ncycl 
  401.  
  402.               do iobl = 1,nbrobl
  403.             MELVAL = lcysig(iobl,icyc)
  404.           if (melval.gt.0) then
  405.             IGMN=MIN(IGAU,VELCHE(/1))
  406.             IBMN=MIN(IB  ,VELCHE(/2))
  407.             sd(iobl)=VELCHE(IGMN,IBMN)
  408.             sigcyc(iobl,icyc) = VELCHE(IGMN,IBMN)
  409.           else 
  410.             sd(iobl) = 0.d0
  411.             sigcyc(iobl,icyc) = 0.d0
  412.           endif
  413.               enddo
  414.  
  415.       PCYC(ICYC) = (SIGCYC(1,ICYC)+SIGCYC(2,ICYC)+SIGCYC(3,ICYC))/3
  416.          if(ib.eq.1.and.igau.eq.1.and.icyc.eq.1) then
  417.            pcymax = pcyc(1)
  418.            pcymin = pcyc(1)
  419.          else
  420.            pcymax = max(pcymax,pcyc(icyc))
  421.            pcymin = min(pcymin,pcyc(icyc))
  422.          endif
  423.        SDCYC(1,icyc)   = (SIGCYC(1,icyc)-SIGCYC(2,icyc))/SQ2
  424.        SDCYC(2,icyc)=(SIGCYC(1,icyc)+SIGCYC(2,icyc)-2.*PCYC(icyc))*SQ3S2
  425.        SDCYC(3,icyc)   =  SIGCYC(4,icyc)
  426.        if(nbrobl.gt.4) then
  427.          SDCYC(4,icyc)   =  SIGCYC(5,icyc)
  428.          SDCYC(5,icyc)   =  SIGCYC(6,icyc)
  429.        endif
  430.  
  431.       if (igau.eq.6) then
  432. *       write(6,*) 'f2-6-icyc',icyc,(sigcyc(iu,icyc),iu = 1,5)
  433.       endif
  434.  
  435. * boucle icyc
  436.             ENDDO
  437.  
  438.       if(nbrobl.le.1) then
  439. *       write(6,*) 'DVPA2, manquent composantes contraintes'
  440.        interr(1) = imodel
  441.        interr(2) = imamod
  442.        call erreur(973)
  443.        return
  444.       endif
  445.  
  446.  
  447. * calculs criteres
  448.              if(icle.eq.1) then
  449.                do jl = 2,ncle
  450.       call FATIG3(sigcyc,pcyc,ycyc,nbrobl,ncycl,jl,
  451.      & cofa1(jl-1),cofa2(jl-1),ycri,SDCYC,ib,igau)
  452.           melval = mcham2.ielval(jl-1)
  453.           velche(igau,ib) = ycri
  454.                enddo
  455.              else
  456.       call FATIG3(sigcyc,pcyc,ycyc,nbrobl,ncycl,icle,
  457.      &cofa1(icle-1),cofa2(icle-1),ycri,SDCYC,ib,igau)
  458.           melval = mcham2.ielval(1)
  459.           velche(igau,ib) = ycri
  460.              endif
  461.  
  462. * sorties
  463.        IF (ICLE.GT.1.and.ICLE.LT.7) THEN
  464.  
  465.         if (ib.eq.1.and.igau.eq.1) then
  466.           melval = mcham2.ielval(2)
  467.           kevdvk = ielche(1,1)
  468.         endif
  469.  
  470.         if (ycri.gt.zecrit) then
  471.             melval = mcham2.ielval(2)
  472.               n = 2
  473.               segini mevoll
  474.               ielche(igau,ib) = mevoll
  475.               ITYEVO='REEL'
  476.               IEVTEX='CYCLE P/TAU '
  477.               IEVTEX(13:20) = mchel1.titche(9:16)
  478.               segini kevoll
  479.               ievoll(1) = kevoll
  480.               ievoll(2) = kevdvk
  481.           if(icle.eq.2.or.icle.eq.3) then
  482.             jg = ncycl
  483.           else
  484.             jg = 1
  485.           endif
  486.               segini mlreel
  487.               iprogx = mlreel
  488.               segini mlree1
  489.               iprogy = mlree1
  490.               TYPX = 'LISTREEL'
  491.               TYPY = 'LISTREEL'
  492.  
  493.           if(icle.eq.2.or.icle.eq.3) then
  494. c       Critère de Dang Van et Papadopoulos
  495.               NUMEVY='REEL'
  496.               NOMEVX = 'P'
  497.               NOMEVY = 'TAU'
  498.             do jcyc = 1,ncycl
  499.              mlree1.prog(jcyc) = ycyc(jcyc)
  500.              prog(jcyc) = pcyc(jcyc)
  501.             enddo
  502.           else
  503.             mlree1.prog(1) = ycyc(2)
  504.             prog(1) = ycyc(1)
  505.             if(icle.eq.4) then
  506. c       Critère de Sines
  507.               NUMEVY='REEL'
  508.               NOMEVX = 'P moyenne'
  509.               NOMEVY = 'sqrt(J2),a'
  510.             elseif(icle.eq.5) then
  511. c       Critère de Crossland
  512.               NUMEVY='REEL'
  513.               NOMEVX = 'P max'
  514.               NOMEVY = 'sqrt(J2),a'
  515.             elseif(icle.eq.6) then
  516. c       Critère de Deperrois
  517.               NUMEVY='REEL'
  518.               NOMEVX = 'P max'
  519.               NOMEVY = 'A(psi)'
  520.             endif
  521.           endif
  522.  
  523.           segdes mlreel,mlree1
  524.           segdes kevoll
  525.           segdes mevoll
  526.  
  527.          else 
  528.            melval = mcham2.ielval(2)
  529.            ielche(igau,ib) = mevnul
  530.          endif
  531.         ENDIF
  532.  
  533. * boucle igau
  534.           enddo
  535. * boucle ib
  536.          ENDDO
  537.  
  538.  
  539.          if(icle.ge.2.and.icle.le.6) then
  540.           kevoll = kevdvk
  541.           mlreel = iprogx
  542.           mlree1 = iprogy
  543. *           WRITE(6,*) 'MINMAX',PCYMIN,PCYMAX
  544.            if (abs(pcymin).le.1.e-6.and.abs(pcymax).le.1.e-6) then
  545.              pcymin = -1.D0
  546.              pcymax = 1.D0
  547.            elseif (PCYMAX.ne.0.) then
  548.              if(abs((pcymax - pcymin)/pcymax).le.0.1) then
  549.                pcymin = pcymin - 0.1*abs(pcymax)
  550.                pcymax = pcymax + 0.1*abs(pcymax)
  551.              endif
  552.            endif
  553.            prog(1) = pcymin
  554.            prog(2) = pcymax
  555.            mlree1.prog(1) = cofa1(icle-1)*(-1)*pcymin + cofa2(icle-1)
  556.            mlree1.prog(2) = cofa1(icle-1)*(-1)*pcymax + cofa2(icle-1)
  557.           segdes mlreel,mlree1
  558.           segdes kevoll
  559.          endif
  560.  
  561.          do lt = 1,ktem
  562.           do iobl = 1, nbrobl
  563.            melval = lcysig(iobl,lt)
  564.            if(melval.gt.0) segdes melval
  565.           enddo
  566.          enddo
  567.          do jcr = 1,mcrit
  568.            melva1 = lcarfa(2*jcr-1,ik)
  569.            if (melva1.gt.0) segdes melva1
  570.            melva2 = lcarfa(2*jcr,ik)
  571.            if (melva2.gt.0) segdes melva2
  572.          enddo
  573.          segsup mdevsi
  574.          segsup mcysig
  575.          segsup mrecyc
  576.         segdes imodel
  577.         mchaml = mchel2.ichaml(ik)
  578.         if(icle.eq.1) then
  579.          do jf = 1,ncle-1
  580.            melval = ielval(jf)
  581.            segdes melval
  582.          enddo
  583.         elseif(icle.ge.2.and.icle.le.6) then
  584.          do jf = 1,2
  585.            melval = ielval(jf)
  586.            segdes melval
  587.          enddo
  588.         endif
  589.         segdes mchaml
  590. * boucle ik
  591.       ENDDO
  592.  
  593.        if(ICF1.gt.0) segdes mchel1
  594.       segdes mmodel,mchel2
  595.       if (lcas1) then
  596.         I0 = I0temd - 1
  597.         do ic = 1,ncycl
  598.           I0 = I0 + 1
  599.            CALL ACCTAB(ITCONT,'ENTIER',I0,X0,' ',LOG0,IP0,
  600.      &                       'MCHAML ',I1,X1,' ',LOG1,IP1)
  601.            MCHELM = ip1
  602.            SEGDES MCHELM
  603.         enddo 
  604.       else
  605.         do jt = 1, mmode1.kmodel(/1)
  606.           imode1 = mmode1.kmodel(jt)
  607.           segdes imode1
  608.         enddo
  609.         segdes mmode1,mchelm
  610.       endif
  611.       ICHOUT = MCHEL2
  612.       RETURN
  613.       END
  614.  
  615.  
  616.  
  617.  
  618.  
  619.  
  620.  
  621.  

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