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    MB234859  26/01/26    21:15:08     12457          
  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,BETA,KSEUIL,EPS_C_CHAPE
  66.       INTEGER ISTRS,JSTRS,KCAS,IRTD
  67.       REAL*8 XZERO,UN,DEUX,XPETIT,PT_DROITE,EPS_C_BARRE
  68.       REAL*8 DINI,D,DT,DC,EPSTIL,EPSTIM,ALFAT,ALFAC,GAMMA
  69.       REAL*8 FTRA,FCOM,GFT,GFC,LCAR
  70.       REAL*8 EPC2,XK,EPCU,K1,K2,KMAX,KPIC
  71.       REAL*8 EPSD0,ACOM,BCOM,ATRA,BTRA
  72.       PARAMETER (XZERO=0.D0 , UN=1.D0 , DEUX=2.D0, XPETIT=1.D-12)
  73.       PARAMETER (YPETIT=1.D0, RA2=SQRT(DEUX))
  74. C
  75. C     recuperation des variables initiales dans les tableaux
  76. C
  77.       YOUN = XMAT(1)
  78.       XNU  = XMAT(2)
  79.       BETA = XMAT(5)
  80.       DINI = VAR0(2)
  81.       KSEUIL = XZERO
  82.       EPS_C_CHAPE = -1.D0*UN
  83. * Numéro de la loi 
  84.       jnplas = wrk53.INPLAS
  85. * Mazars origine appel : 'MAZARS'
  86.       IF (jnplas .EQ. 30) THEN
  87.          ACOM = XMAT(6)
  88.          BCOM = XMAT(7)
  89.          ATRA = XMAT(8)
  90.          BTRA = XMAT(9)
  91.          EPSD0= XMAT(10)
  92.       ELSE 
  93.          FTRA   = XMAT(6)
  94.          FCOM  = XMAT(7)
  95.          KPIC   = XMAT(8)
  96.          EPSD0 = FTRA / YOUN
  97. * Mazars RTC appel : 'MAZARS_RTC' -> 193
  98.          IF (jnplas .EQ. 193) THEN
  99.             GFC   = XMAT(9)
  100.             LCAR  = XMAT(10)
  101.             GFT  = XMAT(11)
  102.             EPC2  = XMAT(12)
  103.             XK = (1.05D0*YOUN*KPIC)/FCOM
  104.             EPCU = (DEUX*GFC)/(LCAR*FCOM)
  105.      1       -(EPC2-KPIC)
  106.             K1 = FCOM/(EPCU-EPC2)
  107.             K2 = FCOM + (K1*EPC2)
  108.             BTRA = LCAR * FTRA / 
  109.      1             (GFT - (LCAR * YOUN * (EPSD0**DEUX) / DEUX))
  110. * Mazars initial appel : 'MAZARS_INI' -> 194
  111.          ELSE IF (jnplas .EQ. 194) THEN
  112.             SRES = XMAT(9)
  113.             TAU = XMAT(10)
  114.             BCOM = UN/(KPIC * XNU * RA2)
  115.             ACOM = (FCOM/(KPIC*YOUN)-(BCOM*EPSD0)) / 
  116.      1     ((EXP (BCOM*EPSD0 - UN))- (BCOM*EPSD0))
  117.             ATRA = UN - (SRES/(YOUN*EPSD0))
  118.             BTRA = (UN + TAU ) / EPSD0
  119. * Mazars linéaire appel : 'MAZARS_LIN' -> 195
  120.          ELSE IF (jnplas .EQ. 195) THEN
  121.             KMAX = XMAT(9)
  122.             BCOM = UN/(KPIC * XNU * RA2)
  123.             ACOM = (FCOM/(KPIC*YOUN)-(BCOM*EPSD0)) / 
  124.      1     ((EXP (BCOM*EPSD0 - UN))- (BCOM*EPSD0))
  125. * Mazars RT appel : 'MAZARS_RT' -> 196
  126.          ELSE IF (jnplas .EQ. 196) THEN
  127.             GFT = XMAT(9)
  128.             LCAR = XMAT(10)
  129.             BCOM = UN / (KPIC * XNU * RA2)
  130.             ACOM = (FCOM/(KPIC*YOUN)-(BCOM*EPSD0)) / 
  131.      1            ((EXP (BCOM*EPSD0 - UN))- (BCOM*EPSD0))
  132.             BTRA = LCAR * FTRA / 
  133.      1             (GFT - (LCAR * YOUN * (EPSD0**DEUX) / DEUX))
  134.          END IF
  135.       END IF 
  136. C
  137. C     calcul des seuils
  138. C
  139. C        calcul de la deformation totale
  140. C
  141.       DO 100 ISTRS=1,NSTRS1
  142.          EPSILO(ISTRS)=EPSILI(ISTRS)+DEPST(ISTRS)
  143. 100   CONTINUE
  144. C
  145. C        calcul des deformations principales
  146. C
  147. C        on reecrit les deformations sous forme matricielle
  148. C
  149. C Modif L.Bode - 14/10/92
  150. C Rajout de IFOUR2 en argument de ENDOCA
  151. *      print*,'on appelle ENDOCB'
  152.       CALL ENDOCB (EPSILO,EPS33,2,IFOUR2)
  153. *      print*,'apres endocb'
  154. C Fin modif L.Bode
  155. C
  156. C        et on diagonalise
  157. C
  158. C Modif L.Bode - 14/10/92
  159. C Pour les elts Coques, on travaille en contraintes planes => JDIM =2
  160. C Pour les elts Massifs JDIM =IDIM
  161. *      print*,'avant JACOB3'
  162.       CALL JACOB3 (EPS33,JDIM,EPSIPP,VALP33)
  163. *      print*,'apres JACOB3'
  164. C Fin modif L.Bode
  165.       IF (ISTEP .EQ. 0 .OR. ISTEP.EQ.2) THEN
  166. C
  167. C       on calcule la matrice de hooke et les contraintes ppales
  168. C
  169.          CMATE = 'ISOTROPE'
  170.          KCAS=1
  171. C Modif L.Bode - 14/10/92
  172. C IFOUR --> IFOUR2 dans appel DOHMAS
  173. *         print*,'avant dohmas'
  174.          CALL DOHMAS(XMAT,CMATE,IFOUR2,NSTRS1,KCAS,DDHOOK,IRTD)
  175. *         print*,'apres dohmas'
  176. C Fin modif L.Bode
  177.          DO 200 ISTRS=1,3
  178.             r_z = XZERO
  179.             DO 210 JSTRS=1,3
  180.               r_z = r_z + DDHOOK(ISTRS,JSTRS)*EPSIPP(JSTRS)
  181. 210         CONTINUE
  182.             SIGP(ISTRS)= r_z
  183. 200      CONTINUE
  184.       END IF
  185. C
  186. C        on complete la deformation dans le cas des contraintes planes
  187. C
  188. C Modif L.Bode - 14/10/92
  189. C IFOUR remplace par IFOUR2
  190.       IF (IFOUR2.EQ. -2) THEN
  191.          EPSIPP(3)= -(EPSIPP(1) + EPSIPP(2))*XNU / (UN-XNU)
  192.       END IF
  193. C Fin modif L.Bode
  194. C
  195. C        on calcule le epsilontild
  196. C
  197.       EPSTIL=MAX( XZERO , EPSIPP(1) )**DEUX +
  198.      1       MAX( XZERO , EPSIPP(2) )**DEUX +
  199.      2       MAX( XZERO , EPSIPP(3) )**DEUX
  200.       EPSTIL=SQRT (EPSTIL)
  201.       epstil=max(xpetit,epstil)
  202.       IF (ISTEP .EQ. 0) THEN
  203.          EPSTIM=EPSTIL
  204.          VARF(1)=EPSTIL
  205.       ELSE IF (ISTEP .EQ. 1) THEN
  206.          VARF(2)=DINI
  207.          VARF(1)=EPSTIL
  208.       ELSE IF (ISTEP .EQ. 2) THEN
  209.          EPSTIM=VAR0(1)
  210.          VARF(1)=EPSTIM
  211.       ELSE
  212.          PRINT*,'DANS MAZARS ISTEP = 0,1,2 ET PAS ',ISTEP
  213.       END IF
  214.       IF ( (ISTEP .EQ. 0) .OR. (ISTEP .EQ. 2)) THEN
  215. C
  216. C      on calcule l'endommagement et les contraintes
  217. C
  218. C        on calcule ALFAT ALFAC DT et DC puis D
  219. C        dans le cas ou le seuil initial est depasse
  220. C
  221.          IF (jnplas .EQ. 193) THEN
  222.          EPS_C_CHAPE = EPSTIM/(RA2*XNU)
  223.          KSEUIL = (0.05D0*XK*KPIC/(1.05D0+XK*(XK - DEUX)))
  224.          ENDIF
  225.          IF (( EPSTIM .GT. EPSD0).OR.(EPS_C_CHAPE .GT. KSEUIL )) THEN
  226. C
  227. C      calcul de l'endommagement
  228. C
  229. C        on calcule le signe des contraintes elastiques
  230. C
  231.             DO 300 ISTRS=1,3
  232.                IF (SIGP(ISTRS).LT. XZERO) THEN
  233.                   SIGPC(ISTRS)=SIGP(ISTRS)
  234.                   SIGPT(ISTRS)=XZERO
  235.                ELSE
  236.                   SIGPT(ISTRS)=SIGP(ISTRS)
  237.                   SIGPC(ISTRS)=XZERO
  238.                END IF
  239. 300         CONTINUE
  240.             TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
  241.             TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
  242. C
  243. C        on calcule les deformations dues aux contraintes positives
  244. C
  245.             DO 400 ISTRS=1,3
  246.                EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
  247. 400         CONTINUE
  248. C
  249. C        on en deduit ALFAT et ALFAC
  250. C
  251.             ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
  252.      1              MAX(XZERO,EPSIPP(2))*EPSILT(2) +
  253.      2              MAX(XZERO,EPSIPP(3))*EPSILT(3)
  254.             ALFAT = ALFAT/(EPSTIL*EPSTIL)
  255.             ALFAC = UN - ALFAT
  256. C
  257. C       modification pour la bi ou tricompression
  258. C
  259.             IF (TRSIGC.LT. -YPETIT .AND. TRSIGT.LT.YPETIT) THEN
  260.                GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+
  261.      1                                 SIGPC(3)*SIGPC(3)
  262.                GAMMA=-SQRT(GAMMA)/TRSIGC
  263.                EPSTIM=EPSTIM*GAMMA
  264.             END IF
  265. C
  266. C        amelioration de la reponse en cisaillement pour beta > 1.
  267. C
  268.             IF (BETA .GT. UN) THEN
  269.                IF ( ALFAT .GT. XPETIT ) THEN
  270.                   ALFAT=ALFAT**BETA
  271.                END IF
  272.                IF ( ALFAC .GT. XPETIT ) THEN
  273.                   ALFAC=ALFAC**BETA
  274.                END IF
  275.             END IF
  276. C
  277. C   trois lois d'évolution:
  278. C      ATRA > 0 : loi de mazars classique
  279. C      -10<ATRA<0 : loi d'évolution exponentielle modifiée pour le GF
  280. C      ATRA < -10 : loi d'évolution linéaire  --> BTRA est alors
  281. C               la déformation pour laquelle la contrainte s'annule C
  282. C
  283. C Loi originellement implémentée
  284.             IF (jnplas .EQ. 30) THEN
  285.                IF (EPSTIM .GT. EPSD0) THEN
  286.                   IF (ATRA .GT. XZERO) THEN
  287.                   DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
  288.      1              ATRA*EXP(BTRA*(EPSD0-EPSTIM))
  289.                   ELSE IF (ATRA . GT. -10.D0) THEN
  290.                   DT=UN - epsd0/epstim*EXP(BTRA*(EPSD0-EPSTIM))
  291.                   ELSE
  292.                      IF(EPSTIM .LT. BTRA) THEN
  293.                         DT=UN - EPSD0*(BTRA - EPSTIM)/EPSTIM/
  294.      1                   (BTRA - EPSD0)
  295.                      ELSE
  296.                         DT= UN
  297.                      ENDIF
  298.                   END IF
  299.                   DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
  300.      1            ACOM*EXP(-BCOM*(EPSTIM-EPSD0))
  301.                ELSE
  302.                   DT=XZERO
  303.                   DC=XZERO
  304.                END IF
  305.             ELSE
  306. c Nouvelle implémentation : M de Carlan 
  307. c Partie traction
  308.                IF (jnplas .EQ. 193.OR.jnplas.EQ.196) THEN
  309.                   IF (EPSTIM .GT. EPSD0) THEN
  310.                      DT=UN - epsd0/epstim*EXP(BTRA*(EPSD0-EPSTIM))
  311.                   ELSE
  312.                      DT=XZERO
  313.                   END IF
  314.                ELSE IF (jnplas .EQ. 194) THEN
  315.                   IF (EPSTIM .GT. EPSD0) THEN
  316.                      DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
  317.      1              ATRA*EXP(BTRA*(EPSD0-EPSTIM))
  318.                   ELSE
  319.                      DT=XZERO
  320.                   END IF
  321.                ELSE 
  322.                   IF (EPSTIM .GT. EPSD0) THEN
  323.                      IF(EPSTIM .LT. KMAX) THEN
  324.                         DT=UN - EPSD0*(KMAX - EPSTIM)/EPSTIM/
  325.      1                   (KMAX - EPSD0)
  326.                      ELSE
  327.                         DT= UN
  328.                      ENDIF
  329.                   ELSE
  330.                      DT= XZERO
  331.                   END IF
  332.                END IF
  333. c Partie compression
  334.                IF (jnplas .EQ. 194.OR.jnplas.EQ.195
  335.      1              .OR.jnplas.EQ.196) THEN
  336.                   IF (EPSTIM .GT. EPSD0) THEN
  337.                      DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
  338.      1                  ACOM*EXP(-BCOM*(EPSTIM-EPSD0))
  339.                   ELSE
  340.                      DC= XZERO
  341.                   END IF
  342.                ELSE
  343.                   IF (EPS_C_CHAPE .GT. KSEUIL ) THEN
  344.                      EPS_C_BARRE = EPS_C_CHAPE / KPIC
  345.                      PT_DROITE = XZERO
  346.                      IF (EPS_C_CHAPE .LT. KPIC) THEN
  347.                         FACT = (XK*EPS_C_BARRE-EPS_C_BARRE*EPS_C_BARRE)
  348.      1                       *FCOM
  349.                         DENOM=((UN+(XK-DEUX)*EPS_C_BARRE)*YOUN
  350.      1                       *EPS_C_CHAPE)
  351.                         PT_DROITE = FACT / DENOM
  352.                      ELSE IF (EPS_C_CHAPE .LT. EPC2) THEN
  353.                         PT_DROITE = FCOM / (YOUN * EPS_C_CHAPE)
  354.                      ELSE IF (EPS_C_CHAPE .LT. EPCU) THEN
  355.                         PT_DROITE = K2 / (YOUN * EPS_C_CHAPE) - 
  356.      1                         K1 / YOUN
  357.                      ELSE
  358.                         PT_DROITE = XZERO
  359.                      END IF
  360.                      DC = UN - PT_DROITE
  361.                   ELSE
  362.                      DC = XZERO
  363.                   END IF
  364.                END IF
  365.             END IF
  366.             DC = MAX(DC,XZERO)
  367.             DT = MAX(DT,XZERO) 
  368.             ALFAT = MAX(ALFAT, XZERO)
  369.             ALFAC = MAX(ALFAC, XZERO)
  370.             D = ALFAT*DT + ALFAC*DC
  371. C
  372. C           on borne la valeur de D a 0.99999999 pour avoir une raideur residuelle evitant des points flottants
  373. C
  374. *pv         D=MIN ( D , UN-1.D-20 )
  375.             D=MIN ( D , UN-1.D-8 )
  376.          ELSE
  377.             D=XZERO
  378.          END IF
  379. C
  380. C           on teste la croissance de D
  381. C
  382.          D=MAX ( D , DINI )
  383. C
  384. C           on le stocke dans les variables internes finales
  385. C
  386.          VARF(2)= D
  387. C
  388. C           on calcule les contraintes finales
  389. C
  390.          CALL MATVE1 (DDHOOK,EPSILO,NSTRS1,NSTRS1,SIGF,2)
  391.          DO 500 ISTRS=1,NSTRS1
  392.             SIGF(ISTRS)=SIGF(ISTRS)*(UN-D)
  393. 500      CONTINUE
  394. C
  395. C         et les deformations inelastiques finales
  396. C
  397.          DO 600 ISTRS=1,NSTRS1
  398.             EPINF(ISTRS)=EPSILO(ISTRS)*D
  399. 600      CONTINUE
  400.  
  401.       ENDIF
  402.  
  403.       RETURN
  404.       END
  405.  
  406.  
  407.  
  408.  
  409.  

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