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

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