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

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