Télécharger la plaquette Cast3M


fr en

Télécharger inosc2.eso

Retour à la liste

Numérotation des lignes : 

inosc2

  1. C INOSC2    SOURCE    PV        15/04/13    21:15:11     8474
  2. CREAD  INOSC1   SOURCE   A1 GIBI   06/12/89 12:32:19
  3. C INOSC1    SOURCE    CHAT      89/06/01    22:06:52
  4.       SUBROUTINE INOSC2(IPT,IPG,DFREQ,XSI,RMAX,TMAX,
  5.      1               DIMALI,DUCTIL,
  6.      2        ALFA,BETA,GAMA,TETA,imod,idema,epse)
  7.       IMPLICIT INTEGER(I-N)
  8.       IMPLICIT REAL*8(A-H,O-Z)
  9. C
  10. C     =====================================================
  11. C                                                         =
  12. C     CALCUL DU MAXIMUM DE LA REPONSE A  UN SIGNAL  D'UN  =
  13. C     OSCILLATEUR NON LINEAIRE AVEC CONDITIONS INITIALES  =
  14. C     EN VITESSE ET EN DEPLACEMENT NULLES.                =
  15. C                                                         =
  16. C     CREATION : 15/06/90                                 =
  17. C     PROGRAMMEUR : A.PINTO AND P.PEGON [ISPRA(I)]        =
  18. C     INSPIREE SUR LA SUBROUTINE 'INOSC1' DE 03/06/87     =
  19. c
  20. c     I.P 05/98 ajout du modele bilineaire, iterations sur
  21. c               la ductilite, changement du modèle Takeda
  22. C                                                         =
  23. c    imod = 1 : modele takeda
  24. c           2          elastoplastique ecrouissage isotrope
  25. c           3          elastoplastique ecrouissage cinematique
  26. c           4          elastique nl
  27. C     =====================================================
  28. C
  29.  
  30. -INC PPARAM
  31. -INC CCOPTIO
  32. -INC SMLREEL
  33. -INC CCREEL
  34.  
  35. * quelques tableaux pour les variables internes du modèle Takeda
  36.       REAL*8 VAR0(21),VARF(21)
  37. * qques tableaux mour utiliser le meme s_p que DYNE (dychec)
  38.       REAL*8 XABSCI(1,4),XORDON(1,4)
  39.  
  40.  
  41.  
  42.       MLREE1=IPT
  43.       SEGACT MLREE1
  44.       MLREE2=IPG
  45.       SEGACT MLREE2
  46.       NN=MLREE1.PROG(/1)
  47.  
  48. C
  49. C     DIFFERENCE ENTRE LES 2 PREMIERS TEMPS SUCCESSIFS : DLTF
  50. C              (PAS DE TEMP CONSTANT)
  51. C
  52.       DLTF=MLREE1.PROG(2)-MLREE1.PROG(1)
  53. **    precision
  54.       epsi = .05d0
  55.       prec = 1d-3
  56.       iprec = 0
  57.       nitmax = 30
  58.  
  59. C
  60. C     CALCUL DE CONSTANTES UTILES
  61. C
  62.       W=2*XPI*DFREQ
  63.       W2=W*W
  64.  
  65. C
  66. C        CALCULO DOS COEFICIENTES DE IWAN
  67. C        ================================
  68.         DW=XSI*W
  69.         D2=XSI*XSI
  70.         A0=EXP(-DW*DLTF)
  71.         A1=W*SQRT(1.0d0-D2)
  72.         AD1=A1*DLTF
  73.         A2=SIN(AD1)
  74.         A3=COS(AD1)
  75.         A4=(2.0d0*D2-1.0d0)/W2
  76.         A5=XSI/W
  77.         A6=2.0d0*A5/W2
  78.         A7=1.0d0/W2
  79.         A8=(A1*A3-DW*A2)*A0
  80.         A9=-(A1*A2+DW*A3)*A0
  81.         A10=A8/A1
  82.         A30=A0/A1
  83.         A12=A30*A2
  84.         A13=A0*A3
  85.         A14=A10*A4
  86.         A15=A12*A4
  87.         A16=A6*A13
  88.         A17=A9*A6
  89.         A11=A0*(DW*A2/A1+A3)
  90.         A21=A10*DW+A9
  91.         A22=A10
  92.         B11=(-A15-A16+A6)/DLTF-A12*A5-A7*A13
  93.         B12=(A15+A16-A6)/DLTF+A7
  94.         B21=(-A14-A17-A7)/DLTF-A10*A5-A9*A7
  95.         B22=(A14+A17+A7)/DLTF
  96.  
  97.         ELAS0=W2
  98.         ELAS=ELAS0
  99.         ELAS1=ELAS0
  100.         if (imod.eq.1 .and. teta.le.0d0) teta = .0001d0
  101.         ELAS2=ELAS1*TETA
  102.  
  103. * qques initialisations pour etre compatible avec dychec
  104.         xabsci(1,1) = 0.d0
  105.         xordon(1,1) = 0.d0
  106.         xjeu = 0d0
  107.         iperm = imod
  108.         if (imod.eq.4) iperm = -2
  109.         nl = 1
  110.         xamo = 0d0
  111.         iannul = 0
  112.         xvit = 0d0
  113.         nip = 4
  114.  
  115. 400     continue
  116.  
  117.         iterm = 0
  118.         iterm2 = 0
  119.         ibisec = 0
  120. 200     continue
  121.         iterm = iterm + 1
  122.         iterm2 = iterm2 + 1
  123. C
  124. C        INICIALIZACAO DE VARIAVEIS
  125. C        ==========================
  126.         F00=0.D0
  127.         F11=0.D0
  128.         D0=0.D0
  129.         V0=0.D0
  130.         FNMDV1=0.D0
  131.         FNMDV=0.D0
  132.         DES0=0.D0
  133.         FORC0=0.D0
  134.         RMAX=0.D0
  135.         epse= 0.d0
  136.         DY=DIMALI/DUCTIL
  137.         fy = dy*elas0
  138.         ITER=0
  139. * qques parametres de Takeda
  140.         dyn = -dy
  141.         fyn = -fy
  142.         fc = .999d0*fy
  143.         fcn = .999d0*fyn
  144.         pinp = beta*fy
  145.         pinn = beta*fyn
  146.         bp = gama/fy
  147. * qques parametres encore pour dychec
  148.       if (imod.ne.1) then
  149.         xabsci(1,2) = dy
  150.         xordon(1,2) = fy
  151.        if (elas2.ge.0d0) then
  152.          xabsci(1,3) = dy + dy
  153.          xordon(1,3)  = fy + elas2*dy
  154.          xabsci(1,4) = 3d0*dy
  155.          xordon(1,4) = fy + elas2*dy*2d0
  156.         else
  157. *        si raideur negative on considere un plateau a .01fy
  158.           xabsci(1,3) = -.9d0*fy/elas2 + dy
  159.           xordon(1,3) = .1d0*fy
  160.           xabsci(1,4) = 2d0*xabsci(1,3)
  161.           xordon(1,4) = .01d0*fy
  162.         endif
  163.         xdplas = 0d0
  164.         xdplac = 0d0
  165.        endif
  166. C
  167. C     CREATION DE LA BOUCLE
  168. C
  169.       DO 100 I=1,NN
  170. C
  171. C        INICIO DAS INTEGRACOES
  172. C        ======================
  173. *  force exterieure a la fin du pas precedent et courant
  174.          F00=F11
  175.          F11=-MLREE2.PROG(I)
  176. C         INICIO DA CONVERGENCIA
  177. C         ======================
  178. 999         ITER=ITER+1
  179. *  on ajoute les forces nl au 2nd membre
  180.          F0=F00-FNMDV
  181.          F1=F11-FNMDV1
  182.          D1=A11*D0+A12*V0+B11*F0+B12*F1
  183.          V1=A21*D0+A22*V0+B21*F0+B22*F1
  184.  
  185. C
  186. C         CALCULO DOS TERMOS NAO-LINEARES NO PASSO
  187. C         ========================================
  188.       DELTAD = D1 - D0
  189.       deltad=max(xpetit,deltad)
  190.  
  191.       if (imod.ne.1) then
  192.  
  193.         call DYCHEC(d1,XDPLAS,XDPLAC,XJEU,IPERM,XABSCI,
  194.      &           XORDON,nl,xfla,nl,NIP,XVIT,XAMO,iannul)
  195.         forca = -xfla
  196.  
  197.        else
  198. * modele Takeda
  199. ** quelques initialisations
  200.        forca = forc0
  201.       if (iter.eq.1) then
  202.        if (i.eq.1) then
  203.         do 1 ii = 1,21
  204.           var0(ii) = 0.d0
  205. 1       continue
  206.        else
  207.         do 2 m = 1,21
  208.           var0(m) = VARF(m)
  209. 2       continue
  210.        endif
  211.       endif
  212.  
  213.       DO 3 k=1,21
  214.         VARF(k)=VAR0(k)
  215. 3     CONTINUE
  216.       IFR1=INT(VARF(1))
  217.  
  218.       IF(IFR1.EQ.0) THEN
  219.         IFR1 = 1
  220.         VARF(2)=elas0
  221.       ENDIF
  222.       DP1 = DELTAD * VARF(2)
  223.  
  224.       CALL DDNSTH(IFR1,IFC1,DP1,FORCA,VARF(21),
  225.      * VARF(2),elas0,elas2,elas2,fc,fcn,fy,fyn,
  226.      * dy,dyn,alfa,alfa,pinp,pinn,bp,bp,
  227.      * VARF(3),VARF(4),VARF(5),
  228.      * VARF(6),VARF(7),VARF(8),VARF(9),
  229.      * VARF(10),VARF(11),VARF(12),VARF(13),
  230.      * VARF(14),VARF(15),VARF(16),VARF(17),
  231.      * VARF(18),VARF(19),VARF(20))
  232.  
  233.       VARF(1)=DBLE(IFR1)
  234.  
  235.       endif
  236.  
  237. C         ANALISE DA CONVERGENCIA EM FORCA
  238. C         ================================
  239. * forces non lineaires à la fin de l'iteration courante
  240.          FICT=-D1*ELAS0+FORCA
  241.          ERRO=ABS((d1-des0)/deltad)
  242.          FNMDV1=FICT
  243.          DES0=D1
  244.  
  245. C
  246. C        ANALISE DE CONVERGENCIA GLOBAL
  247. C        ==============================
  248.         IF(ERRO.gt.prec .and. ITER.lt.200) goto 999
  249.  
  250. * valeurs  a la fin du pas
  251.         ITER=0
  252.         depse = abs(deltad - (forca - forc0)/elas0)
  253.         epse = epse + depse
  254.         FNMDV=FNMDV1
  255.         D0=D1
  256.         V0=V1
  257.         forc0 = forca
  258.  
  259.          IF(ABS(D1).GT.RMAX)THEN
  260.                 RMAX=ABS(D1)
  261.                 TMAX=MLREE1.PROG(I)
  262.          ENDIF
  263.  
  264.               IF(ITER.GT.200)THEN
  265.       print*,'SPON: temps :',MLREE1.PROG(I),'nb d iterations > 200'
  266.               ENDIF
  267. C
  268. 100   CONTINUE
  269.  
  270.       if (idema.eq.0) goto 300
  271.  
  272.       duct2 = rmax/dy
  273. *
  274. **  petite correction de la limite elastique si elle est trop petite
  275. ** (ca arrive pour des frequences elevees)
  276. **  pour une frequence tres grande ca revient a commencer les iterations
  277. **  avec une ductilite de 1./.9
  278.  
  279. c      if (iterm2.eq.1) print*,'dfreq',dfreq
  280. c      print*,'duct2=',duct2,'iterm=',iterm,'iterm2=',iterm2,
  281. c     &           'freq=',dfreq
  282.  
  283. c       print*,dimali,duct2
  284.  
  285.  
  286.       if (duct2.gt. (5d0*ductil).and.iterm.eq.1
  287.      &                .and. (abs(teta)).lt.1d0) then
  288.        if (iterm2.eq.1 .or.
  289.      &            (duct2.gt. (6d0*ductil).and.iterm2.le.3)) then
  290.         iterm = iterm - 1
  291.         dimali = dimali*(.9d0*ductil -
  292.      &           ((.9d0*ductil - 1)* exp(-.02*(duct2-5d0*ductil))))
  293.         goto 200
  294.        endif
  295.       endif
  296.  
  297.         diaux1 = dimali
  298.  
  299.       if ((abs(duct2 - ductil)).gt.
  300.      &            (epsi*ductil).and.iterm2.le.nitmax) then
  301.        if (iterm.eq.1) dimali = rmax
  302.        if (iterm.ge.2) then
  303.          aa = (duct2 - ductil)*(duct1-ductil)
  304.          if (aa.gt.0.d0. and. ibisec.eq.0) then
  305.             dimali = rmax
  306.          else if (ibisec.eq.0) then
  307. * premiere dichotomie
  308.           ibisec = 1
  309.           dimali = (dimali + dimal1)/2d0
  310.           ducaux = duct1
  311.           diaux2 = dimal1
  312.          else
  313. * dichotomies suivantes
  314.            aa2 = (duct2 - ductil)*(duct3-ductil)
  315.              if (aa.lt.0.d0) then
  316.                dimali = (dimali + dimal1)/2d0
  317.                ducaux = duct1
  318.                diaux2 = dimal1
  319.              else if (aa2.lt.0.d0) then
  320.                dimali = (dimali + dimal3)/2d0
  321.                ducaux = duct3
  322.                diaux2 = dimal3
  323.              else
  324. *  dans le cas ou la dichotomie ne marche pas
  325.               dimali = rmax
  326.              endif
  327.          endif
  328.        endif
  329.            duct1 = duct2
  330.            dimal1 = diaux1
  331.         if (ibisec.eq.1) then
  332.           duct3 = ducaux
  333.           dimal3 = diaux2
  334.         endif
  335.             goto200
  336.       else if  (iterm2.gt.nitmax) then
  337. * on essaie 2 fois encore avec un critere plus petit
  338.          if (iprec.le.1) then
  339.            prec = prec/10d0
  340.            iprec = iprec + 1
  341.            goto 400
  342.          endif
  343.            REAERR(1)=dfreq
  344.            REAERR(2)=duct2
  345.            CALL ERREUR(-326)
  346.       endif
  347.       rmax = dy
  348.  
  349. 300   continue
  350.       SEGDES MLREE1
  351.       SEGDES MLREE2
  352.  
  353.  
  354.       RETURN
  355.       END
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  

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