Télécharger la plaquette Cast3M


fr en

Télécharger dohot1.eso

Retour à la liste

Numérotation des lignes : 

dohot1

  1. C DOHOT1    SOURCE    OF166741  25/02/21    21:15:55     12166          
  2.  
  3.       SUBROUTINE DOHOT1(IVAMAT,NMATT,IVACON,NSTRS,IVARI,NVART,TRAC,
  4.      &             LTRAC,IGAU,IB,MATE,MAPL,XPREC,DTPS,IFOU,LHOOK,
  5.      &             DDHOOK,IRET)
  6. *_______________________________________________________________________
  7. *
  8. *  entrees :
  9. *  --------
  10. *
  11. *     ivamat  pointeur sur un segment mptval contenant les materiaux
  12. *     nmatt   nombre de composantes du materiau
  13. *     ivacon  pointeur sur un segment mptval contenant les contraintes
  14. *     nstrs   nombre de composantes de contraintes
  15. *     ivari   pointeur sur un segment mptval contenant les variables
  16. *             internes
  17. *     nvart   nombre de composantes de variables internes
  18. *     trac    courbe de traction
  19. *     ltrac   taille du tableau trac
  20. *     igau    numero du point de gauss ou l'on veut le mat. de hooke
  21. *     ib      numero de l'element ou l'on veut le mat. de hooke
  22. *     mate    numeror du materiau lineaire
  23. *     mapl    numero du materiau plastique
  24. *     xprec   flottant precision
  25. *     dtps    flottant pas de temps pour les modeles visqueux
  26. *     ifou    option ifour de ccoptio
  27. *     lhook   taille de la matrice de hooke
  28. *
  29. *  sorties :
  30. *  --------
  31. *
  32. *     ddhook  matrice de hooke tangente pour l'element ib
  33. *     iret    = 1 si ok
  34. *             = < 0 : no de l erreur a appeler dans ktanga
  35. *             = 0 : autre erreur
  36. *
  37. *     passage aux nouveaux chamelem par jm campenon le 05/91
  38. *_______________________________________________________________________
  39. *
  40.       IMPLICIT INTEGER(I-N)
  41.       IMPLICIT REAL*8(A-H,O-Z)
  42.  
  43. -INC PPARAM
  44. -INC CCOPTIO
  45.  
  46. -INC SMCHAML
  47.  
  48. -INC TMPTVAL
  49.  
  50.       SEGMENT WRK5
  51.          INTEGER NTRAT,NTRAC
  52.       ENDSEGMENT
  53. *
  54.       SEGMENT WTRAV
  55.          REAL*8 TXR(IDIM,IDIM)
  56.       ENDSEGMENT
  57.  
  58.       SEGMENT WRKMAT
  59.         REAL*8 VALMAT(NMATT)
  60.         REAL*8 D1HO(LHOOK,LHOOK),ROTH(LHOOK,LHOOK)
  61.         REAL*8 XLOC(3,3),XGLOB(3,3)
  62.       ENDSEGMENT
  63. *
  64.       PARAMETER(XZER=0.D0,UN=1.D0,DEUX=2.D0,TROIS=3.D0,QUATRE=4.D0)
  65.       PARAMETER(CINQ=5.D0,SIX=6.D0,HUIT=8.D0,XNEUF=9.D0,DOUZE=12.D0)
  66.       PARAMETER(UNTIER=.33333 33333 33333D0,TRDEMI=1.5D0)
  67.       PARAMETER(DETIER=.66666 66666 66666D0,UNDEMI=0.5D0)
  68.       PARAMETER(UNQU=.25D0)
  69. *
  70.       DIMENSION DDHOOK(LHOOK,*)
  71.       DIMENSION TRAC(*)
  72. *
  73. *     arguments pour l appel a DOHMAS et DOHMAO
  74.       DIMENSION   VELA(2)
  75.       CHARACTER*8 CMATE
  76. *
  77.       DIMENSION XSTRS(8),DEVIAT(8),XMAT(30)
  78.       DIMENSION DDF(8),DDG(8),DHDF(8),DHDG(8)
  79.       DIMENSION XVAR(14),XKTA(3,3)
  80.       DIMENSION PMATRIC(6,6),GMATRIC(6,6),GDMATR(6,6),GINV(6,6),COEF3(6)
  81. *     matrice P pour le modele viscoplastique parfait
  82. *     la matrice sert a passer de la contrainte a la contrainte deviatorique
  83.       DATA PMATRIC /2.D0,-1.D0,-1.D0,0.D0,0.D0,0.D0,
  84.      &              -1.D0,2.D0,-1.D0,0.D0,0.D0,0.D0,
  85.      &              -1.D0,-1.D0,2.D0,0.D0,0.D0,0.D0,
  86.      &               0.D0, 0.D0,0.D0,6.D0,0.D0,0.D0,
  87.      &               0.D0, 0.D0,0.D0,0.D0,6.D0,0.D0,
  88.      &               0.D0, 0.D0,0.D0,0.D0,0.D0,6.D0/
  89. *
  90.       DATA COEF3 /1.D0,1.D0,1.D0,2.D0,2.D0,2.D0/
  91. *
  92. *-----------------------------------------------------------------------
  93. *                formulation massive (MFR=1)
  94. * Pour autres formulations (ex : coques) voir dans dohota.eso
  95. *-----------------------------------------------------------------------
  96.       IF ((MATE.EQ.1).AND.(MAPL.EQ.54)) THEN
  97.          WRK5 = IRET
  98.       ELSE IF ((MATE.EQ.2).OR.(MATE.EQ.3).OR.(MATE.EQ.4)) THEN
  99.          WTRAV = IRET
  100.       ENDIF
  101. *
  102.       IRET=0
  103.       CALL ZERO(DDHOOK,LHOOK,LHOOK)
  104.       CALL ZERO(XMAT,30,1)
  105. *
  106. *-----------------------------------------------------------------------
  107. *                     MATERIAU ISOTROPE
  108. *-----------------------------------------------------------------------
  109.       IF (MATE.EQ.1) THEN
  110. *-----------------------------------------------------------------------
  111. C        Initialisations
  112.          XXHH =0.D0
  113.          SIGY =0.D0
  114.          ALP  =0.D0
  115. C        Initialisation de ILOGEL : sortie de PENT1
  116. C        =1 on doit/peut calculer DDHOOK tangent ; =0 DDHOOK elastique
  117.          ILOGEL = 0
  118. *
  119. C        Contraintes
  120.          MPTVAL=IVACON
  121.          DO 50 ICOMP=1,NSTRS
  122.             MELVAL=IVAL(ICOMP)
  123.             IGMN=MIN(IGAU,VELCHE(/1))
  124.             IBMN=MIN(IB  ,VELCHE(/2))
  125.             XSTRS(ICOMP)=VELCHE(IGMN,IBMN)
  126.   50     CONTINUE
  127. C
  128. C        Deformation plastique
  129. C
  130.          MPTVAL=IVARI
  131.          MELVAL=IVAL(1)
  132.          IGMN=MIN(IGAU,VELCHE(/1))
  133.          IBMN=MIN(IB  ,VELCHE(/2))
  134.          EPS =VELCHE(IGMN,IBMN)
  135. C        si autres variables internes
  136.          IF (NVART.EQ.NSTRS+1) THEN
  137.             DO 60 ICOMP=1,NSTRS
  138.                MELVAL=IVAL(ICOMP+1)
  139.                IGMN=MIN(IGAU,VELCHE(/1))
  140.                IBMN=MIN(IB  ,VELCHE(/2))
  141.                XSTRS(ICOMP)= XSTRS(ICOMP)-VELCHE(IGMN,IBMN)
  142.   60        CONTINUE
  143.          ENDIF
  144. C
  145. C Recuperation des caracteristiques elastiques isotropes
  146. C
  147.          MPTVAL=IVAMAT
  148. C        Module d Young
  149.          MELVAL=IVAL(1)
  150.          IGMN=MIN(IGAU,VELCHE(/1))
  151.          IBMN=MIN(IB  ,VELCHE(/2))
  152.          YOU =VELCHE(IGMN,IBMN)
  153. C        Coefficient de Poisson
  154.          MELVAL=IVAL(2)
  155.          IGMN=MIN(IGAU,VELCHE(/1))
  156.          IBMN=MIN(IB  ,VELCHE(/2))
  157.          XNU =VELCHE(IGMN,IBMN)
  158. *
  159. *-----------------------------------------------------------------------
  160. *               Recuperation des caracteristiques
  161. *          Calcul eventuel de la contrainte equivalente
  162. *-----------------------------------------------------------------------
  163. C  Plastique parfait
  164.          IF (MAPL.EQ.1) THEN
  165. *
  166.             MPTVAL=IVAMAT
  167.             MELVAL=IVAL(3)
  168.             IGMN=MIN(IGAU,VELCHE(/1))
  169.             IBMN=MIN(IB  ,VELCHE(/2))
  170.             SIGY =VELCHE(IGMN,IBMN)
  171. *
  172. C           Calcul de la contrainte equivalente
  173.             XXXX   = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  174.      &               + XSTRS(3)*XSTRS(3)
  175.             YYYY   = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  176.      &               + XSTRS(3)*XSTRS(1)
  177.             SIGET2 = XXXX - YYYY
  178.             IF (NSTRS.GE.4) THEN
  179.               DO 61 I=4,NSTRS
  180.                 SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  181. 61            CONTINUE
  182.             ENDIF
  183.             SIGET = SQRT(SIGET2)
  184. *
  185. C           A t on plastifie ?
  186.             CALL PENT1(EPS,SIGET,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC,
  187.      &                  YOUTA,ILOGEL)
  188.             IF (ILOGEL.EQ.-1) THEN
  189.                IRET=-275
  190.                RETURN
  191.             ENDIF
  192. *
  193. C  Drucker simple : ALFA * Tr(S) +  Seq = K
  194.          ELSE IF (MAPL.EQ.3) THEN
  195. *
  196.             MPTVAL=IVAMAT
  197.             MELVAL=IVAL(3)
  198.             IGMN=MIN(IGAU,VELCHE(/1))
  199.             IBMN=MIN(IB  ,VELCHE(/2))
  200.             XLTR =VELCHE(IGMN,IBMN)
  201.             MELVAL=IVAL(4)
  202.             IGMN=MIN(IGAU,VELCHE(/1))
  203.             IBMN=MIN(IB  ,VELCHE(/2))
  204.             XLCS =VELCHE(IGMN,IBMN)
  205. *
  206.             DEN = ABS(XLCS) + XLTR
  207.             IF(DEN.EQ.0.D0) THEN
  208.               IRET=-366
  209.               RETURN
  210.             ENDIF
  211. C           ALP = ALFA ; XXHH = BETA = 1 ; SIGY = K
  212.             ALP  = (ABS(XLCS) - XLTR)/DEN
  213.             XXHH = UN
  214.             SIGY = DEUX*ABS(XLCS)*XLTR/DEN
  215. *
  216. *            XMAT(1) = 2.0D0*ABS(ALP)*SIGY/DEN
  217. *            XMAT(2) = (ABS(ALP) - SIGY)/DEN
  218. *            XMAT(3) = 1.D0
  219.             XMAT(1) = ALP
  220.             XMAT(2) = XXHH
  221.             XMAT(3) = SIGY
  222.             XMAT(4)=XMAT(1)
  223.             XMAT(5)=XMAT(2)
  224.             XMAT(6)=XMAT(1)
  225.             XMAT(7)=XMAT(2)
  226.             XMAT(8)=XMAT(3)
  227.             XMAT(9)=0.D0
  228. *
  229. C           Calcul de la contrainte equivalente
  230.             XXXX   = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  231.      &               + XSTRS(3)*XSTRS(3)
  232.             YYYY   = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  233.      &               + XSTRS(3)*XSTRS(1)
  234.             SIGET2 = XXXX - YYYY
  235.             IF (NSTRS.GE.4) THEN
  236.               DO 62 I=4,NSTRS
  237.                 SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  238. 62            CONTINUE
  239.             ENDIF
  240.             SIGET = SQRT(SIGET2)
  241. *
  242. C           A t on plastifie ?
  243.             SIDD = ALP * (XSTRS(1)+XSTRS(2)+XSTRS(3)) + SIGET
  244.             CALL PENT1(EPS,SIDD,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC,
  245.      &                  YOUTA,ILOGEL)
  246.             IF (ILOGEL.EQ.-1) THEN
  247.                IRET=-275
  248.                RETURN
  249.             ENDIF
  250. *
  251. C  Cinematique
  252.          ELSE IF (MAPL.EQ.4) THEN
  253. *
  254.             MPTVAL=IVAMAT
  255.             MELVAL=IVAL(3)
  256.             IGMN=MIN(IGAU,VELCHE(/1))
  257.             IBMN=MIN(IB  ,VELCHE(/2))
  258.             SIGY =VELCHE(IGMN,IBMN)
  259. *
  260.             MELVAL=IVAL(4)
  261.             IGMN=MIN(IGAU,VELCHE(/1))
  262.             IBMN=MIN(IB  ,VELCHE(/2))
  263.             XXHH=VELCHE(IGMN,IBMN)
  264. *
  265. C           Calcul de la contrainte equivalente
  266.             XXXX   = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  267.      &               + XSTRS(3)*XSTRS(3)
  268.             YYYY   = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  269.      &               + XSTRS(3)*XSTRS(1)
  270.             SIGET2 = XXXX - YYYY
  271.             IF (NSTRS.GE.4) THEN
  272.               DO 63 I=4,NSTRS
  273.                 SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  274. 63            CONTINUE
  275.             ENDIF
  276.             SIGET = SQRT(SIGET2)
  277. *
  278. C           A t on plastifie ?
  279.             CALL PENT1(EPS,SIGET,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC,
  280.      &                  YOUTA,ILOGEL)
  281.             IF (ILOGEL.EQ.-1) THEN
  282.                IRET=-275
  283.                RETURN
  284.             ENDIF
  285. *
  286. C  Plastique isotrope
  287.          ELSE IF (MAPL.EQ.5) THEN
  288. *
  289. C           Calcul de la contrainte equivalente
  290.             XXXX   = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  291.      &               + XSTRS(3)*XSTRS(3)
  292.             YYYY   = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  293.      &               + XSTRS(3)*XSTRS(1)
  294.             SIGET2 = XXXX - YYYY
  295.             IF (NSTRS.GE.4) THEN
  296.               DO 64 I=4,NSTRS
  297.                 SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  298. 64            CONTINUE
  299.             ENDIF
  300.             SIGET = SQRT(SIGET2)
  301. *
  302. C           A t on plastifie ?
  303.             CALL PENT1(EPS,SIGET,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC,
  304.      &                  YOUTA,ILOGEL)
  305.             IF (ILOGEL.EQ.-1) THEN
  306.                IRET=-275
  307.                RETURN
  308.             ENDIF
  309. *
  310. C  Drucker Prager
  311.          ELSE IF (MAPL.EQ.15) THEN
  312. *
  313.             MPTVAL=IVAMAT
  314.             DO 1106 IC=3,NMATT
  315.               MELVAL=IVAL(IC)
  316.               XMAT(IC-2)=VELCHE(IGMN,IBMN)
  317. 1106        CONTINUE
  318. *
  319. ** ALP = ALFA ou ETA ; XXHH = BETA ou MU ; SIGY = K ou KL+H*EPS
  320.             IF(EPS.EQ.0.D0) THEN
  321.               ALP=XMAT(6)
  322.               XXHH=XMAT(7)
  323.               SIGY=XMAT(8)
  324.             ELSE
  325.               ALP=XMAT(1)
  326.               XXHH=XMAT(2)
  327.               SIGY=XMAT(3)+XMAT(9)*EPS
  328.             ENDIF
  329. *
  330. C           Calcul de la contrainte equivalente
  331.             XXXX   = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  332.      &               + XSTRS(3)*XSTRS(3)
  333.             YYYY   = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  334.      &               + XSTRS(3)*XSTRS(1)
  335.             SIGET2 = XXXX - YYYY
  336.             IF (NSTRS.GE.4) THEN
  337.               DO 66 I=4,NSTRS
  338.                 SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  339. 66            CONTINUE
  340.             ENDIF
  341.             SIGET = SQRT(SIGET2)
  342. *
  343. C           A t on plastifie ?
  344.             SIDD = ALP * (XSTRS(1)+XSTRS(2)+XSTRS(3)) + XXHH * SIGET
  345.             CALL PENT1(EPS,SIDD,YOU,SIGY,XXHH,TRAC,LTRAC,MAPL,XPREC,
  346.      &                  YOUTA,ILOGEL)
  347.             IF (ILOGEL.EQ.-1) THEN
  348.                IRET=-275
  349.                RETURN
  350.             ENDIF
  351. *
  352. C  Viscoplastique parfait
  353.          ELSE IF (MAPL.EQ.43) THEN
  354. *
  355. C           Recuperation des caracteristiques
  356.             MPTVAL=IVAMAT
  357.             MELVAL = IVAL(3)
  358.             IGMN=MIN(IGAU,VELCHE(/1))
  359.             IBMN=MIN(IB  ,VELCHE(/2))
  360.             SIGY = VELCHE(IGMN,IBMN)
  361.             MELVAL = IVAL(5)
  362.             IGMN=MIN(IGAU,VELCHE(/1))
  363.             IBMN=MIN(IB  ,VELCHE(/2))
  364.             CK = VELCHE(IGMN,IBMN)
  365.             MELVAL = IVAL(4)
  366.             IGMN=MIN(IGAU,VELCHE(/1))
  367.             IBMN=MIN(IB  ,VELCHE(/2))
  368.             CN = VELCHE(IGMN,IBMN)
  369.             ILOGEL = 1
  370. *
  371. C  Betocy
  372.          ELSE IF (MAPL.EQ.54) THEN
  373. C
  374.             MPTVAL=IVAMAT
  375.             DO 1107 IC=1,NMATT-2
  376.               MELVAL=IVAL(IC)
  377.               IGMN=MIN(IGAU,VELCHE(/1))
  378.               IBMN=MIN(IB  ,VELCHE(/2))
  379.               XMAT(IC)=VELCHE(IGMN,IBMN)
  380. 1107        CONTINUE
  381. C          Lecture Variables internes
  382.             MPTVAL=IVARI
  383.             DO 1108 IC=1,NVART
  384.               MELVAL=IVAL(IC)
  385.               IGMN=MIN(IGAU,VELCHE(/1))
  386.               IBMN=MIN(IB  ,VELCHE(/2))
  387.               XVAR(IC)=VELCHE(IGMN,IBMN)
  388. 1108        CONTINUE
  389.             ILOGEL = 1
  390. *
  391. C  Rotating Crack
  392.          ELSE IF (MAPL.EQ.55) THEN
  393. C          Lecture materiau
  394.             MPTVAL=IVAMAT
  395.             DO 1109 IC=1,NMATT
  396.               MELVAL=IVAL(IC)
  397.               IGMN=MIN(IGAU,VELCHE(/1))
  398.               IBMN=MIN(IB  ,VELCHE(/2))
  399.               XMAT(IC)=VELCHE(IGMN,IBMN)
  400. 1109        CONTINUE
  401. C          Lecture Variables internes
  402.             MPTVAL=IVARI
  403.             DO 1110 IC=1,NVART
  404.               MELVAL=IVAL(IC)
  405.               IGMN=MIN(IGAU,VELCHE(/1))
  406.               IBMN=MIN(IB  ,VELCHE(/2))
  407.               XVAR(IC)=VELCHE(IGMN,IBMN)
  408. 1110        CONTINUE
  409.             ILOGEL = 1
  410. *
  411. * BCN (start)
  412.          ELSE IF ((MAPL.EQ.111).OR.(MAPL.EQ.112).OR.(MAPL.EQ.113)) THEN
  413.             IF (MAPL.EQ.112.OR.MAPL.EQ.113) THEN
  414. C           J2 (MAPL=112) or Rounded Hyperbolic Mohr-Coulomb (MAPL=113)
  415. C           get auxiliary variable for CTM
  416.               MPTVAL=IVARI
  417.               MELVAL=IVAL(2)
  418.               IGMN=MIN(IGAU,VELCHE(/1))
  419.               IBMN=MIN(IB  ,VELCHE(/2))
  420.               AUX1 =VELCHE(IGMN,IBMN)
  421.             ELSE IF (MAPL.EQ.111) THEN
  422. C           MRS-Lade (MAPL=111)
  423. C           get auxiliary variable for CTM
  424.               MPTVAL=IVARI
  425.               MELVAL=IVAL(1)
  426.               IGMN=MIN(IGAU,VELCHE(/1))
  427.               IBMN=MIN(IB  ,VELCHE(/2))
  428.               EPSCON =VELCHE(IGMN,IBMN)
  429.               MPTVAL=IVARI
  430.               MELVAL=IVAL(2)
  431.               IGMN=MIN(IGAU,VELCHE(/1))
  432.               IBMN=MIN(IB  ,VELCHE(/2))
  433.               EPSCAP =VELCHE(IGMN,IBMN)
  434.               MPTVAL=IVARI
  435.               MELVAL=IVAL(3)
  436.               IGMN=MIN(IGAU,VELCHE(/1))
  437.               IBMN=MIN(IB  ,VELCHE(/2))
  438.               AUX1CON =VELCHE(IGMN,IBMN)
  439.               MPTVAL=IVARI
  440.               MELVAL=IVAL(4)
  441.               IGMN=MIN(IGAU,VELCHE(/1))
  442.               IBMN=MIN(IB  ,VELCHE(/2))
  443.               AUX1CAP =VELCHE(IGMN,IBMN)
  444.             ENDIF
  445. C           The three models: read material constants in XMAT
  446.             XMAT(1)=YOU
  447.             XMAT(2)=XNU
  448.             XMAT(3)=0.D0
  449.             XMAT(4)=0.D0
  450.             MPTVAL=IVAMAT
  451.             DO IC=3,NMATT
  452.               MELVAL=IVAL(IC)
  453.               IGMN=MIN(IGAU,VELCHE(/1))
  454.               IBMN=MIN(IB  ,VELCHE(/2))
  455.               XMAT(IC+2)=VELCHE(IGMN,IBMN)
  456.             ENDDO
  457.             ILOGEL = 1
  458.          ENDIF
  459. * BCN (end)
  460. *
  461. *-----------------------------------------------------------------------
  462. *                 Calcul de la matrice elastique
  463. *-----------------------------------------------------------------------
  464.          VELA(1) = YOU
  465.          VELA(2) = XNU
  466. * MATE = 1 --> CMATE = ISOTROPE
  467.          CMATE = 'ISOTROPE'
  468.          CALL DOHMAS(VELA,CMATE,IFOU,LHOOK,1,DDHOOK,IRETOU)
  469.          IF (IRETOU.EQ.0) THEN
  470.            IRET = -81
  471.            RETURN
  472.          ENDIF
  473.          IRET = 1
  474. *
  475. *        ILOGEL = 0 on a fini, =1 on calcule la matrice tangente
  476.          IF (ILOGEL.EQ.0) RETURN
  477. *        series de Fourier, on rend la matrice elastique
  478.          IF (IFOU.EQ.1) RETURN
  479. *
  480. *-----------------------------------------------------------------------
  481. *                 Calcul de la matrice tangente
  482. *-----------------------------------------------------------------------
  483. c  MODELE MRS_LADE
  484.          IF (MAPL.EQ.111) THEN
  485. *
  486. *        notice de l operateur MODE: "option de calcul : plan defo,axis"
  487. *        sinon, on rend la matrice elastique
  488.            IF ((IFOU.EQ.-1).OR.(IFOU.EQ.0)) THEN
  489.              nescri =0
  490.              nues   =6
  491. C            mrs-lade requiere siempre derivacion numerica
  492.              nnumer=3
  493.              deltax=2.D0**(int(log10(1.D-6)/log10(2.D0)))
  494.              call mtc_MRSMAC(DDHOOK,LHOOK,XSTRS,NSTRS,EPSCON,EPSCAP,
  495.      &                       AUX1CON,AUX1CAP,XMAT,nescri,nues,
  496.      &                       nnumer,deltax,kerre)
  497. *            kerre = 1 : formulacion no disponible
  498.              IF (KERRE.EQ.1) THEN
  499.                IRET = -328
  500.              ENDIF
  501.            ENDIF
  502. *
  503. c  MODELE J2
  504.          ELSE IF (MAPL.EQ.112) THEN
  505. *        notice de l operateur MODE: "option de calcul : plan defo,axis"
  506. *          sinon, on rend la matrice elastique
  507.            IF ((IFOU.EQ.-1).OR.(IFOU.EQ.0)) THEN
  508.              nescri   =0
  509.              nues     =6
  510.              call mtc_j2(DDHOOK,LHOOK,XSTRS,NSTRS,EPS,AUX1,XMAT,
  511.      &                   nescri,nues,kerre)
  512. *          kerre vaut toujours 0
  513.            ENDIF
  514. *
  515. c  MODELE RH_COULOMB
  516.          ELSE IF (MAPL.EQ.113) THEN
  517. *        notice de l operateur MODE: "option de calcul : plan defo,axis"
  518. *          sinon, on rend la matrice elastique
  519.            IF ((IFOU.EQ.-1).OR.(IFOU.EQ.0)) THEN
  520.              nescri   =0
  521.              nues     =6
  522.              call mtc_rhmc(DDHOOK,LHOOK,XSTRS,NSTRS,AUX1,XMAT,
  523.      &                     nescri,nues,kerre)
  524. *          kerre vaut toujours 0
  525.            ENDIF
  526. *
  527. C  Viscoplastique parfait
  528.          ELSE IF (MAPL.EQ.43) THEN
  529. *
  530. *          tout sauf contraintes planes
  531. *          sinon on rend la matrice elastique
  532.            IF ((IFOU.EQ.-1).OR.(IFOU.EQ.-3).OR.(IFOU.EQ.0).OR.
  533.      &      (IFOU.EQ.2).OR.(IFOU.EQ.3).OR.(IFOU.EQ.7).OR.(IFOU.EQ.9)
  534.      &      .OR.(IFOU.EQ.11).OR.(IFOU.EQ.12).OR.(IFOU.EQ.14).OR.
  535.      &      (IFOU.EQ.15)) THEN
  536. *            SIGMA EQUIVALENT: (3/2(s':s'))**0.5D0
  537.              XXXX   = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  538.      &               + XSTRS(3)*XSTRS(3)
  539.              YYYY   = XSTRS(1)*XSTRS(2) + XSTRS(2)*XSTRS(3)
  540.      &               + XSTRS(3)*XSTRS(1)
  541.              SIGET2 = XXXX - YYYY
  542.              IF (NSTRS.GE.4) THEN
  543.                DO 67 I=4,NSTRS
  544.                  SIGET2 = SIGET2 + TROIS*XSTRS(I)*XSTRS(I)
  545. 67             CONTINUE
  546.              ENDIF
  547.              SIGEQ = SQRT(SIGET2)
  548. *
  549.              XMOY=(XSTRS(1)+XSTRS(2)+XSTRS(3))*UNTIER
  550.              DO 200 IA=1,3
  551.                DEVIAT(IA)=XSTRS(IA)- XMOY
  552.   200        CONTINUE
  553.              IF (NSTRS.GE.4) THEN
  554.                DO 300 IA=4,NSTRS
  555.                  DEVIAT(IA)=XSTRS(IA)
  556.   300          CONTINUE
  557.              ENDIF
  558. *
  559.              SYN = ( SIGEQ - SIGY ) / CK
  560.              IF (SYN .GT. XPREC) THEN
  561. *            on ne travaille que si l'evolution est plastique
  562. *
  563.                AUX6 = UNDEMI * YOU / (UN+XNU)
  564. *
  565. *              ON CALCULE LA MATRICE G
  566.                COEF1 = SYN ** CN / SIGEQ*DTPS*2.D0*AUX6/3.D0
  567.                COEF2=SYN**(CN-1.D0)*(SYN/SIGEQ+CN/CK)/SIGET2*DTPS
  568.      &                *2.D0*AUX6
  569.                DO 610 I=1,NSTRS
  570.                  DO 600 J=1,NSTRS
  571.                    GMATRIC(I,J)=COEF1*PMATRIC(I,J)+
  572.      &                       COEF2*DEVIAT(I)*DEVIAT(J)*COEF3(I)
  573.  600             CONTINUE
  574.                  GMATRIC(I,I)=GMATRIC(I,I)+1.D0
  575.  610           CONTINUE
  576.  
  577.                CALL INVHOK(GMATRIC,GINV,NSTRS)
  578. *
  579. *              on calcule InvG . D
  580.                DO 650 I=1,NSTRS
  581.                  DO 640 J=I,NSTRS
  582.                    GDMATR(I,J)=0.D0
  583.                    DO 630 K=1,NSTRS
  584.                      GDMATR(I,J)=GDMATR(I,J)+GINV(I,K)*DDHOOK(K,J)
  585.  630               CONTINUE
  586.                    GDMATR(J,I)=GDMATR(I,J)
  587.  640             CONTINUE
  588.  650           CONTINUE
  589. *              on remet tout dans DHOOK
  590.                DO 670 I=1,NSTRS
  591.                  DO 660 J=1,NSTRS
  592.                    DDHOOK(I,J)=GDMATR(I,J)
  593.  660             CONTINUE
  594.  670           CONTINUE
  595.              ENDIF
  596.            ENDIF
  597. *
  598. C  Betocy
  599.          ELSE IF (MAPL.EQ.54) THEN
  600. *
  601. *          notice de l operateur MODE : "en 2D contraintes planes"
  602. *          sinon on rend la matrice elastique
  603.            IF (IFOU.EQ.-2) THEN
  604.              CALL BETOKT(XMAT,TRAC,XVAR,XSTRS,NTRAT,NTRAC,DDHOOK,KERRE)
  605. *            kerre toujours nul
  606.            ENDIF
  607. *
  608. C  Rotating Crack
  609.          ELSE IF (MAPL.EQ.55) THEN
  610. *
  611. * notice de MATE : "la direction de fissuration change a chaque pas"
  612. * possible en 1D ? dans ROTKTA DDHOOK(4,4), en 1D on va jusque 3
  613. *          sinon on rend la matrice elastique
  614.            IF (IFOU.EQ.-2) THEN
  615.              CALL ROTKTA(XMAT,XVAR,XSTRS,DDHOOK,KERRE)
  616. *            kerre = 2 "impossible d inverser"
  617.              IF (KERRE.EQ.2) IRET = 0
  618.            ENDIF
  619. *
  620. C  Plastique parfait, cinematique, isotrope
  621.          ELSE IF ((MAPL.EQ.1).OR.(MAPL.EQ.4).OR.(MAPL.EQ.5)) THEN
  622. *
  623.            IF (IFOU.EQ.-2) THEN
  624. *          2D contraintes planes
  625.              XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  626.              YYYY = XSTRS(1)*XSTRS(2)
  627.              ZZZZ = XSTRS(4)*XSTRS(4)
  628.              AUX4=UNTIER/(UN-XNU*XNU)
  629.              U=(DEUX*XSTRS(1)-XSTRS(2)+XNU*(DEUX*XSTRS(2)-
  630.      &               XSTRS(1)))*AUX4
  631.              V=(DEUX*XSTRS(2)-XSTRS(1)+XNU*(DEUX*XSTRS(1)-
  632.      &               XSTRS(2)))*AUX4
  633.              W= XSTRS(4)/(UN+XNU)
  634.              A= CINQ*XXXX - HUIT*YYYY
  635.              B= CINQ*YYYY - DEUX*XXXX
  636.              C= DETIER*(XXXX-YYYY)+DEUX*ZZZZ
  637.              D= (A+DEUX*XNU*B)*AUX4  +  SIX*ZZZZ/(UN+XNU)
  638.              XNUM=TROIS*YOU*(YOU-YOUTA)
  639.              DENOM=DEUX*YOUTA*C + (YOU - YOUTA)*D
  640.              BETA=XNUM/DENOM
  641. *
  642.              DDHOOK(1,1)=DDHOOK(1,1)-BETA*U*U
  643.              DDHOOK(2,1)=DDHOOK(2,1)-BETA*V*U
  644.              DDHOOK(4,1)=DDHOOK(4,1)-BETA*W*U
  645.              DDHOOK(1,2)=DDHOOK(1,2)-BETA*U*V
  646.              DDHOOK(2,2)=DDHOOK(2,2)-BETA*V*V
  647.              DDHOOK(4,2)=DDHOOK(4,2)-BETA*W*V
  648.              DDHOOK(1,4)=DDHOOK(1,4)-BETA*U*W
  649.              DDHOOK(2,4)=DDHOOK(2,4)-BETA*V*W
  650.              DDHOOK(4,4)=DDHOOK(4,4)-BETA*W*W
  651. *
  652.            ELSE IF ((IFOU.EQ.4).OR.(IFOU.EQ.8)) THEN
  653. *          DYCZ,GYCZ
  654.              XXXX = XSTRS(1)*XSTRS(1) + XSTRS(2)*XSTRS(2)
  655.              YYYY = XSTRS(1)*XSTRS(2)
  656.              AUX4=UNTIER/(UN-XNU*XNU)
  657.              U=(DEUX*XSTRS(1)-XSTRS(2)+XNU*(DEUX*XSTRS(2)-
  658.      &               XSTRS(1)))*AUX4
  659.              V=(DEUX*XSTRS(2)-XSTRS(1)+XNU*(DEUX*XSTRS(1)-
  660.      &               XSTRS(2)))*AUX4
  661.              A= CINQ*XXXX - HUIT*YYYY
  662.              B= CINQ*YYYY - DEUX*XXXX
  663.              C= DETIER*(XXXX-YYYY)
  664.              D= (A+DEUX*XNU*B)*AUX4
  665.              XNUM=TROIS*YOU*(YOU-YOUTA)
  666.              DENOM=DEUX*YOUTA*C + (YOU - YOUTA)*D
  667.              BETA=XNUM/DENOM
  668. *
  669.              DDHOOK(1,1)=DDHOOK(1,1)-BETA*U*U
  670.              DDHOOK(2,1)=DDHOOK(2,1)-BETA*V*U
  671.              DDHOOK(1,2)=DDHOOK(1,2)-BETA*U*V
  672.              DDHOOK(2,2)=DDHOOK(2,2)-BETA*V*V
  673. *
  674.            ELSE IF ((IFOU.EQ.5).OR.(IFOU.EQ.10).OR.(IFOU.EQ.13)
  675.      &               .OR.(IFOU.EQ.6)) THEN
  676. *          CYDZ,CYGZ,AXCZ,CYCZ
  677. *            comment adapter les formules 2D contraintes planes au 1D CY.. ?
  678. *            pour le moment, on rend la matrice elastique
  679. *
  680.            ELSE
  681. *          tout sauf contraintes planes (1D aussi)
  682. *            deviateur
  683.              XMOY=(XSTRS(1)+XSTRS(2)+XSTRS(3))*UNTIER
  684.              DO 201 IA=1,3
  685.                DEVIAT(IA)=XSTRS(IA)- XMOY
  686.   201        CONTINUE
  687.              IF (NSTRS.GE.4) THEN
  688.                DO 301 IA=4,NSTRS
  689.                  DEVIAT(IA)=XSTRS(IA)
  690.   301          CONTINUE
  691.              ENDIF
  692. *
  693.              XNUM = XNEUF * YOU * ( YOU - YOUTA)
  694.              XXXX = TROIS * ( YOU - YOUTA ) +
  695.      &                DEUX * ( UN + XNU ) * YOUTA
  696.              DENOM= DEUX * SIGET2 * ( UN + XNU ) * XXXX
  697.              BETA = XNUM / DENOM
  698. *
  699.              DO 400 ICOMP1=1,NSTRS
  700.                DO 400 ICOMP2=1,NSTRS
  701.                  DDHOOK(ICOMP1,ICOMP2)= DDHOOK(ICOMP1,ICOMP2) -
  702.      &                         BETA*DEVIAT(ICOMP1)*DEVIAT(ICOMP2)
  703.  400         CONTINUE
  704. *
  705.            ENDIF
  706. *
  707. C  Drucker simple, Drucker Prager
  708.          ELSE IF ((MAPL.EQ.3).OR.(MAPL.EQ.15)) THEN
  709. *
  710. *          si SIGET = 0 (avec tr(S) dans le critere, c est possible)
  711. *          on rend la matrice elastique (sinon on divise par SIGET=0)
  712.            IF (SIGET.EQ.0.) RETURN
  713.            CALL ZERO(DDF,NSTRS,1)
  714.            CALL ZERO(DDG,NSTRS,1)
  715. *
  716.            IF ((IFOU.EQ.-2).OR.(IFOU.EQ.4).OR.(IFOU.EQ.8)) THEN
  717. *          2D contraintes planes et 1D DYCZ GYCZ
  718.              SIGM = (XSTRS(1)-XSTRS(2)*0.5D0)/SIGET
  719.              DDF(1)=XMAT(1) + XMAT(2)*SIGM
  720.              DDG(1)=XMAT(4) + XMAT(5)*SIGM
  721.              SIGM = (XSTRS(2)-XSTRS(1)*0.5D0)/SIGET
  722.              DDF(2)=XMAT(1) + XMAT(2)*SIGM
  723.              DDG(2)=XMAT(4) + XMAT(5)*SIGM
  724.              IF (NSTRS.GE.4) THEN
  725.                DDF(4)=6.D0*XMAT(2)*XSTRS(4)/SIGET
  726.                DDG(4)=6.D0*XMAT(5)*XSTRS(4)/SIGET
  727.              ENDIF
  728. *
  729.            ELSE IF ((IFOU.EQ.5).OR.(IFOU.EQ.10).OR.(IFOU.EQ.13)) THEN
  730. *          1D CYDZ CYGZ et AXCZ
  731.              SIGM = (XSTRS(1)-XSTRS(3)*0.5D0)/SIGET
  732.              DDF(1)=XMAT(1) + XMAT(2)*SIGM
  733.              DDG(1)=XMAT(4) + XMAT(5)*SIGM
  734.              SIGM = (XSTRS(3)-XSTRS(1)*0.5D0)/SIGET
  735.              DDF(3)=XMAT(1) + XMAT(2)*SIGM
  736.              DDG(3)=XMAT(4) + XMAT(5)*SIGM
  737. *
  738.            ELSE IF (IFOU.EQ.6) THEN
  739. *          1D CYCZ
  740.              SIGM = (XSTRS(1)-XSTRS(2)*0.5D0)/SIGET
  741.              DDF(1)=XMAT(1) + XMAT(2)*SIGM
  742.              DDG(1)=XMAT(4) + XMAT(5)*SIGM
  743. *
  744.            ELSE
  745.              SIGM = (XSTRS(1)-XSTRS(2)*0.5D0 -XSTRS(3)*0.5D0)/SIGET
  746.              DDF(1)=XMAT(1) + XMAT(2)*SIGM
  747.              DDG(1)=XMAT(4) + XMAT(5)*SIGM
  748.              SIGM = (XSTRS(2)-XSTRS(1)*0.5D0 -XSTRS(3)*0.5D0)/SIGET
  749.              DDF(2)=XMAT(1) + XMAT(2)*SIGM
  750.              DDG(2)=XMAT(4) + XMAT(5)*SIGM
  751.              SIGM = (XSTRS(3)-XSTRS(1)*0.5D0 -XSTRS(2)*0.5D0)/SIGET
  752.              DDF(3)=XMAT(1) + XMAT(2)*SIGM
  753.              DDG(3)=XMAT(4) + XMAT(5)*SIGM
  754.              IF (NSTRS.GE.4) THEN
  755.                DO 421 I=4,NSTRS
  756.                  DDF(I)=6.D0*XMAT(2)*XSTRS(I)/SIGET
  757.                  DDG(I)=6.D0*XMAT(5)*XSTRS(I)/SIGET
  758.  421           CONTINUE
  759.              ENDIF
  760. *
  761.            ENDIF
  762. *
  763.            DO 422 I=1,NSTRS
  764.              DHDF(I)=0.D0
  765.              DHDG(I)=0.D0
  766.              DO 422 J=1,NSTRS
  767.                DHDF(I)=DHDF(I)+DDHOOK(I,J)*DDF(J)
  768.                DHDG(I)=DHDG(I)+DDHOOK(I,J)*DDG(J)
  769.  422       CONTINUE
  770. *
  771.            DENOM = XMAT(9)* SQRT(2.D0*XMAT(4)*XMAT(4)+
  772.      &                                     XMAT(5)*XMAT(5))
  773.            DO 423 I=1,NSTRS
  774.              DENOM = DENOM + DDF(I)*DHDG(I)
  775.  423       CONTINUE
  776.            DO 320 ICOMP1=1,NSTRS
  777.              DO 320 ICOMP2=1,NSTRS
  778.                DDHOOK(ICOMP1,ICOMP2)=DDHOOK(ICOMP1,ICOMP2)
  779.      &                       -  DHDG(ICOMP1)*DHDF(ICOMP2)/DENOM
  780.  320       CONTINUE
  781. *
  782.          ELSE
  783. *          on rend la matrice elastique pour les autres modeles
  784.          ENDIF
  785. *
  786. *-----------------------------------------------------------------------
  787. *                     MATERIAU UNIDIRECTIONNEL (MATE=4)
  788. *                      Acier unidirectionnel (MAPL=40)
  789. *-----------------------------------------------------------------------
  790.       ELSE IF ((MATE.EQ.4).AND.(MAPL.EQ.40).AND.(IFOU.LE.0)) THEN
  791. *-----------------------------------------------------------------------
  792.          DO 530 IC=1,3
  793.             MPTVAL=IVAMAT
  794.             MELVAL=IVAL(IC)
  795.             IGMN=MIN(IGAU,VELCHE(/1))
  796.             IBMN=MIN(IB  ,VELCHE(/2))
  797.             XMAT(IC)=VELCHE(IGMN,IBMN)
  798.  530     CONTINUE
  799.          V1X=TXR(1,1)*XMAT(2)+TXR(1,2)*XMAT(3)
  800.          V1Y=TXR(2,1)*XMAT(2)+TXR(2,2)*XMAT(3)
  801.          XMAT(2)=V1X
  802.          XMAT(3)=V1Y
  803. C     Lecture Variables internes
  804.          MPTVAL=IVARI
  805.          MELVAL=IVAL(4)
  806.          IGMN=MIN(IGAU,VELCHE(/1))
  807.          IBMN=MIN(IB  ,VELCHE(/2))
  808.          XVAR(1)=VELCHE(IGMN,IBMN)
  809. *
  810.          CALL UNIAKT(XMAT,XVAR,IFOU,DDHOOK,KERRE)
  811. * kerre toujours nul
  812.          IRET = 1
  813. *
  814. *-----------------------------------------------------------------------
  815. *                          AUTRE MATERIAU :
  816. *            UNIDIRECTIONNEL (sauf MAPL=40 et IFOU = -3,-2,-1,0)
  817. *            ORTHOTROPE
  818. *            ANISOTROPE
  819. *-----------------------------------------------------------------------
  820.       ELSE
  821. *-----------------------------------------------------------------------
  822.          IF (MATE.EQ.2) THEN
  823.            CMATE = 'ORTHOTRO'
  824.          ELSE IF (MATE.EQ.3) THEN
  825.            CMATE = 'ANISOTRO'
  826.          ELSE IF (MATE.EQ.4) THEN
  827.            CMATE = 'UNIDIREC'
  828.          ELSE
  829. C        cet operateur ne fonctionne pas pour le materiau %m1:8
  830.            IRET = -834
  831.            RETURN
  832.          ENDIF
  833.          SEGINI,WRKMAT
  834.          MPTVAL = IVAMAT
  835.          SEGACT,MPTVAL
  836.          DO 9004 IM=1,NMATT
  837.             IF (IVAL(IM).NE.0) THEN
  838.                MELVAL=IVAL(IM)
  839.                IBMN=MIN(IB  ,VELCHE(/2))
  840.                IGMN=MIN(IGAU,VELCHE(/1))
  841.                if (ibmn.gt.0.and.igmn.gt.0) then
  842.                  VALMAT(IM)= VELCHE(IGMN,IBMN)
  843. C*             else
  844. C*               VALMAT(IM)=0.D0
  845.                endif
  846. C*          ELSE
  847. C*             VALMAT(IM)=0.D0
  848.             ENDIF
  849.  9004    CONTINUE
  850.          CALL DOHMAO(VALMAT,CMATE,IFOU,IDIM,TXR,XLOC,XGLOB,D1HO,
  851.      &               ROTH,DDHOOK,LHOOK,1,IRETOU)
  852.          SEGSUP,WRKMAT
  853.          IF (IRETOU.EQ.0) THEN
  854.            IRET = -81
  855.          ELSE IF (IRETOU.EQ.1) THEN
  856.            IRET = 1
  857.          ELSE
  858.            IRET = 0
  859.          ENDIF
  860.       ENDIF
  861. *
  862.       RETURN
  863.       END
  864.  
  865.  
  866.  

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