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

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