Télécharger la plaquette Cast3M


fr en

Télécharger cmazar.eso

Retour à la liste

Numérotation des lignes : 

cmazar

  1. C CMAZAR    SOURCE    PV090527  24/10/09    21:15:03     12032          
  2.       SUBROUTINE CMAZAR (WRK52,WRK53,WRK54,WRKK2,NSTRS1,NVARI,
  3.      1                   ICARA,JDIM,IFOUR2)
  4. C MAZARS    SOURCE    AM        98/12/23    21:38:30     3409
  5. C
  6. C
  7. C     variables en entree
  8. C
  9. C
  10. C     WRK0     pointeur sur un segment deformation au pas precedent
  11. C
  12. C     WRK1     pointeur sur un segment increment de deformation
  13. C
  14. C     WRKK2     pointeur sur un segment variables internes au pas precedent
  15. C
  16. C     WRK5      pointeur sur un segment de deformations inelastiques
  17. C
  18. C     XMATER     constantes du materiau
  19. C
  20. C     NSTRS1      nombre de composantes dans les vecteurs des contraintes
  21. C                                        et les vecteurs des deformations
  22. C
  23. C     NVARI      nombre de variables internes (doit etre egal a 2)
  24. C
  25. C     NMATT      nombre de constantes du materiau
  26. C
  27. C     ISTEP      flag utilise pour separer les etapes dans un calcul non local
  28. C                ISTEP=0 -----> calcul local
  29. C                ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
  30. C                ISTEP=2 -----> calcul non local etape 2 on continue le calcul
  31. C                               a partir des seuils moyennes
  32. C
  33. C Modif L.Bode - 14/10/92
  34. C Nouveaux parametres en entree
  35. C     JDIM       Dimension de travail
  36. C                ( Coques JDIM =2 , Massifs JDIM = IDIM )
  37. C     IFOUR2     Type de formulation
  38. C                ( Coques IFOUR2 = -2 => contraintes planes ,
  39. C                  Massifs IFOUR2 = IFOUR)
  40. C
  41. C     variables en sortie
  42. C
  43. C     VARINF      variables internes finales
  44. C
  45. C     SIGMAF      contraintes finales
  46. C
  47. C     C. LA BORDERIE MARS 1992
  48. C     declaration des variables
  49. C
  50. C
  51.       IMPLICIT INTEGER(I-N)
  52.       IMPLICIT REAL*8(A-H,O-Z)
  53. -INC PPARAM
  54. -INC CCOPTIO
  55. -INC DECHE
  56. *
  57.       SEGMENT WRKK2
  58.         REAL*8 EPSILI(NSTRS1)
  59.         REAL*8 EPSILO(NSTRS1)
  60.       ENDSEGMENT
  61.  
  62.       INTEGER NSTRS1,NVARI
  63.       REAL*8 EPS33(3,3),EPSIPP(3),EPSILT(3),VALP33(3,3)
  64.       REAL*8 SIGP(3),SIGPT(3),SIGPC(3),TRSIGT,TRSIGC
  65.       REAL*8 YOUN,XNU,EPSD0,ACOM,BCOM,ATRA,BTRA,BETA
  66.       INTEGER ISTRS,JSTRS,KCAS,IRTD
  67.       REAL*8 XZERO,UN,DEUX,XPETIT
  68.       REAL*8 DINI,D,DT,DC,EPSTIL,EPSTIM,ALFAT,ALFAC,GAMMA
  69.       PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)
  70.  
  71. C
  72. C
  73. C     recuperation des variables initiales dans les tableaux
  74. C
  75. C
  76.       N=NSTRS1
  77.       CMATE = 'ISOTROPE'
  78.       YOUN = XMAT(1)
  79.       XNU  = XMAT(2)
  80.       EPSD0= XMAT(5)
  81.       ACOM = XMAT(6)
  82.       BCOM = XMAT(7)
  83.       ATRA = XMAT(8)
  84.       BTRA = XMAT(9)
  85.       BETA = XMAT(10)
  86.       DINI = VAR0(2)
  87.  
  88. C
  89. C     calcul des seuils
  90. C
  91. C
  92. C        calcul de la deformation totale
  93. C
  94. C      SEGINI WRK3
  95.       DO 100 ISTRS=1,NSTRS1
  96.          EPSILO(ISTRS)=EPSILI(ISTRS)+DEPST(ISTRS)
  97. 100   CONTINUE
  98. C
  99. C        calcul des deformations principales
  100. C
  101. C
  102. C        on reecrit les deformations sous forme matricielle
  103. C
  104. C Modif L.Bode - 14/10/92
  105. C Rajout de IFOUR2 en argument de ENDOCA
  106. *      print*,'on appelle ENDOCB'
  107.       CALL ENDOCB (EPSILO,EPS33,2,IFOUR2)
  108. *      print*,'apres endocb'
  109. C Fin modif L.Bode
  110. C
  111. C        et on diagonalise
  112. C
  113. C Modif L.Bode - 14/10/92
  114. C Pour les elts Coques, on travaille en contraintes planes => JDIM =2
  115. C Pour les elts Massifs JDIM =IDIM
  116. *      print*,'avant JACOB3'
  117.       CALL JACOB3 (EPS33,JDIM,EPSIPP,VALP33)
  118. *      print*,'apres JACOB3'
  119. C Fin modif L.Bode
  120.       IF (ISTEP .EQ. 0 .OR. ISTEP.EQ.2) THEN
  121. C
  122. C       on calcule la matrice de hooke et les contraintes ppales
  123. C
  124.          CMATE = 'ISOTROPE'
  125.          KCAS=1
  126. C Modif L.Bode - 14/10/92
  127. C IFOUR --> IFOUR2 dans appel DOHMAS
  128. *         print*,'avant dohmas'
  129.          CALL DOHMAS(XMAT,CMATE,IFOUR2,NSTRS1,KCAS,DDHOOK,IRTD)
  130. *         print*,'apres dohmas'
  131. C Fin modif L.Bode
  132.          DO 200 ISTRS=1,3
  133.             SIGP(ISTRS)= XZERO
  134.             DO 210 JSTRS=1,3
  135.                SIGP(ISTRS)=SIGP(ISTRS)+DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS)
  136. 210         CONTINUE
  137. 200      CONTINUE
  138.       END IF
  139. C
  140. C        on complete la deformation dans le cas des contraintes planes
  141. C
  142. C Modif L.Bode - 14/10/92
  143. C IFOUR remplace par IFOUR2
  144.       IF (IFOUR2.EQ. -2) THEN
  145.          EPSIPP(3)= -(EPSIPP(1) + EPSIPP(2))*XNU / (UN-XNU)
  146.       END IF
  147. C Fin modif L.Bode
  148. C
  149. C        on calcule le epsilontild
  150. C
  151.       EPSTIL=MAX( XZERO , EPSIPP(1) )**2 +
  152.      1       MAX( XZERO , EPSIPP(2) )**2 +
  153.      2       MAX( XZERO , EPSIPP(3) )**2
  154.       EPSTIL=SQRT (EPSTIL)
  155.       epstil=max(xpetit,epstil)
  156.       IF (ISTEP .EQ. 0) THEN
  157.          EPSTIM=EPSTIL
  158.          VARF(1)=EPSTIL
  159.       ELSE IF (ISTEP .EQ. 1) THEN
  160.          VARF(2)=DINI
  161.          VARF(1)=EPSTIL
  162.       ELSE IF (ISTEP .EQ. 2) THEN
  163.          EPSTIM=VAR0(1)
  164.          VARF(1)=EPSTIM
  165.       ELSE
  166.          PRINT*,'DANS MAZARS ISTEP = 0,1,2 ET PAS ',ISTEP
  167.       END IF
  168.       IF ( (ISTEP .EQ. 0) .OR. (ISTEP .EQ. 2)) THEN
  169. C
  170. C      on calcule l'endommagement et les contraintes
  171. C
  172. C
  173. C        on calcule ALFAT ALFAC DT et DC puis D
  174. C        dans le cas ou le seuil initial est depasse
  175. C
  176.          IF ( EPSTIM .GT. EPSD0) THEN
  177. C
  178. C      calcul de l'endommagement
  179. C
  180. C
  181. C        on calcule le signe des contraintes elastiques
  182. C
  183.             DO 300 ISTRS=1,3
  184.                IF (SIGP(ISTRS).LT. XZERO) THEN
  185.                   SIGPC(ISTRS)=SIGP(ISTRS)
  186.                   SIGPT(ISTRS)=XZERO
  187.                ELSE
  188.                   SIGPT(ISTRS)=SIGP(ISTRS)
  189.                   SIGPC(ISTRS)=XZERO
  190.                END IF
  191. 300         CONTINUE
  192.             TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
  193.             TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
  194. C
  195. C        on calcule les deformations dues aux contraintes positives
  196. C
  197.             DO 400 ISTRS=1,3
  198.                EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
  199. 400         CONTINUE
  200. C
  201. C        on en deduit ALFAT et ALFAC
  202. C
  203.             ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
  204.      1              MAX(XZERO,EPSIPP(2))*EPSILT(2) +
  205.      2              MAX(XZERO,EPSIPP(3))*EPSILT(3)
  206.             ALFAT = ALFAT/(EPSTIL*EPSTIL)
  207.             ALFAC = UN - ALFAT
  208. C
  209. C       modification pour la bi ou tricompression
  210. C
  211.             IF (TRSIGC.LT. -XPETIT .AND. TRSIGT.LT.XPETIT) THEN
  212.                GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+
  213.      1                                 SIGPC(3)*SIGPC(3)
  214.                GAMMA=-SQRT(GAMMA)/TRSIGC
  215.                EPSTIM=EPSTIM*GAMMA
  216.             END IF
  217. C
  218. C        amelioration de la reponse en cisaillement pour beta > 1.
  219. C
  220.             IF (BETA .GT. UN) THEN
  221.                IF ( ALFAT .GT. XPETIT ) THEN
  222.                   ALFAT=ALFAT**BETA
  223.                END IF
  224.                IF ( ALFAC .GT. XPETIT ) THEN
  225.                   ALFAC=ALFAC**BETA
  226.                END IF
  227.             END IF
  228.  
  229. C
  230. C        on calcule DT et DC puis D
  231. C        dans le cas ou le seuil initial est depasse
  232. C
  233. C        on est oblige de verifier car on a pu multiplier par gamma
  234. C
  235. C            IF (EPSTIM .GT. EPSD0) THEN
  236. C               DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
  237. C     1            ATRA*EXP(-BTRA*(EPSTIM-EPSD0))
  238. C               DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
  239. C     1            ACOM*EXP(-BCOM*(EPSTIM-EPSD0))
  240. C            ELSE
  241. C               DT=XZERO
  242. C               DC=XZERO
  243. C            END IF
  244. CLB 2010/01/20
  245. C            trois lois d'évolution:
  246. C             ATRA > 0 : loi de mazars classique
  247. C           -10<ATRA<0 : loi d'évolution exponentielle modifiée pour le GF
  248. C             ATRA < -10 : loi d'évolution linéaire  --> BTRA est alors
  249. C                  la déformation pour laquelle la contrainte s'annule C
  250. C
  251.             IF (EPSTIM .GT. EPSD0) THEN
  252.                IF (ATRA .GT. 0.D0) THEN
  253.                  DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
  254.      1              ATRA*EXP(BTRA*(EPSD0-EPSTIM))
  255.                ELSE IF (ATRA . GT. -10.) THEN
  256.                  DT=UN - epsd0/epstim*EXP(BTRA*(EPSD0-EPSTIM))
  257.                ELSE
  258.                  IF(EPSTIM .LT. BTRA) THEN
  259.                    DT=UN - EPSD0*(BTRA - EPSTIM)/EPSTIM/(BTRA - EPSD0)
  260.                  ELSE
  261.                    DT=1.D0
  262.                  ENDIF
  263.                END IF
  264.                DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
  265.      1            ACOM*EXP(-BCOM*(EPSTIM-EPSD0))
  266.             ELSE
  267.  
  268.                DT=XZERO
  269.                DC=XZERO
  270.             END IF
  271.             D = ALFAT*DT + ALFAC*DC
  272. C
  273. C           on borne la valeur de D a 0.99999999 pour avoir une raideur residuelle evitant des points flottants
  274. C
  275. *pv         D=MIN ( D , UN-1.D-20 )
  276.             D=MIN ( D , UN-1.D-8 )
  277.          ELSE
  278.             D=XZERO
  279.          END IF
  280. C
  281. C           on teste la croissance de D
  282. C
  283.          D=MAX ( D , DINI )
  284. C
  285. C           on le stocke dans les variables internes finales
  286. C
  287.          VARF(2)= D
  288. C
  289. C           on calcule les contraintes finales
  290. C
  291.          CALL MATVE1 (DDHOOK,EPSILO,NSTRS1,NSTRS1,SIGF,2)
  292.          DO 500 ISTRS=1,NSTRS1
  293.             SIGF(ISTRS)=SIGF(ISTRS)*(UN-D)
  294. 500      CONTINUE
  295. C
  296. C         et les deformations inelastiques finales
  297. C
  298.          DO 600 ISTRS=1,NSTRS1
  299.             EPINF(ISTRS)=EPSILO(ISTRS)*D
  300. 600      CONTINUE
  301.  
  302.       ENDIF
  303.       RETURN
  304.       END
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  

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