Télécharger la plaquette Cast3M


fr en

Télécharger pbefe1.eso

Retour à la liste

Numérotation des lignes : 

pbefe1

  1. C PBEFE1    SOURCE    PV        22/04/22    21:15:11     11344          
  2.       SUBROUTINE PBEFE1 (XMAT,VAR0,VARF,SIG0,SIGF,DEPST,XCAR,IB,XEL,
  3.      &  ISTEP)
  4. C
  5. C====&===1=========2=========3=========4=========5=========6=========7==
  6. C COMMENTAIRES :   SUBROUTINE PERMETTANT DE METTRE EN OEUVRE LE 
  7. C                           MODELE EFEM 
  8. C
  9. C TRAITS       : 
  10. C
  11. C
  12. C
  13. C AUTEUR       : F. RICCARDI - CEA/DEN/DANS/DM2S/SEMT/EMSI
  14. C                B. RICHARD  - CEA/DEN/DANS/DM2S/SEMT/EMSI
  15. C====&===1=========2=========3=========4=========5=========6=========7==
  16. C
  17. C-------------------------------------------------------------------------------------------
  18. C-----DECLARATION GENERALES
  19. C-------------------------------------------------------------------------------------------
  20.          IMPLICIT INTEGER(I-N)
  21.          IMPLICIT REAL*8(A-H,O-Z)
  22. C-------------------------------------------------------------------------------------------
  23. C-----DECLARATION SEGMENTS
  24. C-------------------------------------------------------------------------------------------
  25.  
  26. C-----REAL 8
  27.          REAL*8 EPSIM(2,2),UNIT(2,2),SIGMAM(2,2),XEL(3,3)
  28.          REAL*8 XMAT(*),SIG0(4),SIGF(4),VAR0(*),VARF(*),DEPST(4),XCAR(*)
  29.          REAL*8 AD(4),XRHO(4),XECR(4),TD(2),XOUV(2),XNIE(2),XIE(2)
  30.          REAL*8 XN(2),XT(2),P1(2),P2(2),P3(2)
  31.          REAL*8 BNOD1(1,2),BNOD2(1,2),BNOD3(1,2),BGRAD(3,2)
  32.          REAL*8 XANG(3),G(3,2),D(3,3),N(3,2)
  33.          REAL*8 EPAUX(2,2),EPSVOIGT(3,1),XG(2)
  34.          REAL*8 SIGVOI(3,1),M(2,3),T(2),SIGMA(2,2),VEP(2,2)
  35.          REAL*8 VAP(2),DP_MAX(2),VP(1,2),FF(3)
  36.          REAL*8 VPT(2),SIGAUX(2,2),T_CN0(1,1)
  37.          REAL*8 XNT(1,2),XTT(1,2),T_CN(1),T_CT(1)
  38.          REAL*8 EPSD(3),EPSENR_V(3),EPSENR(2,2),SIG(2,2),T_CN1(2)
  39.          REAL*8 UN1(2,2),UN2(2,2),EPSI_MAT(2,2),XOUVT(1,2)
  40.          REAL*8 T_C(2),T_U(2),EPS1(2,2),EPS2(2,2),EPS4(2,2),EPS7(2,2)
  41.          REAL*8 EPS1P(2,2),EPS2P(2,2),EPS4P(2,2),EPS7P(2,2),TEND1(2,2)
  42.          REAL*8 TEND2(2,2),TEND4(2,2),TEND7(2,2),TEND(2,2),TENDP(2)
  43.          REAL*8 K(2,2),XDG(3,2),F(2,2),RES(2),KINV(2,2),XDGN(3)
  44.          REAL*8 XNP(2),XNPT(1,2),DUN(2),EPNL(2,2)
  45. C-------------------------------------------------------------------------------------------
  46. C-----DECLARATION PARAMETRES NUMERIQUES
  47. C-------------------------------------------------------------------------------------------
  48. C
  49. C-----NOMBRE MAXIMUM D ITERATIONS POUR LA CONVERGENCE SUR LES CONTRAINTES PLANES
  50.          ITMAX1 = 1000
  51. C-----NOMBRE MAXIMUM D ITERATIONS POUR LA CONVERGENCE SUR LA CONTINUITE
  52.          ITMAX2 = 1000
  53. C
  54. C-----CHOIX DE L ECROUISSAGE
  55. C       NECR = 1 -- ECROUISSAGE DE BARGELLINI
  56. C       NECR = 2 -- ECROUISSAGE LINEAIRE
  57. C       NECR = 3 -- ECROUISSAGE EXPONENTIEL
  58.          NECR = 3
  59.          XKAPPA = 70.0D0
  60. C-----TOLERANCE SUR LA CONDITION DE CONTRAINTES PLANES
  61.          TOL1 = 1.0D-6
  62. C
  63. C-----TOLERANCE SUR LA CONDITION DE CONTINUITE
  64.          TOL2 = 1.0D-6
  65. C
  66. C-------------------------------------------------------------------------------------------
  67. C-----CONSTRUCTION DE LA MATRICE IDENTITE D ORDRE DEUX
  68. C-------------------------------------------------------------------------------------------
  69. C
  70.          DO I=1,2
  71.              DO J=1,2
  72.                  IF (I.EQ.J) THEN
  73.                      UNIT(I,J) = 1.0D0
  74.                  ELSE
  75.                      UNIT(I,J) = 0.0D0
  76.                  ENDIF
  77.              ENDDO
  78.          ENDDO
  79. C-------------------------------------------------------------------------------------------
  80. C-----RECUPERATION DES PROPRIETES MATERIAUX
  81. C-------------------------------------------------------------------------------------------
  82. C
  83. C-----MODULE D YOUNG
  84.          YOUN = XMAT(1)
  85. C
  86. C-----COEFFICIENT DE POISSON
  87.          XNU = XMAT(2)  
  88. C
  89. C-----COEFFICIENTS DE LAME
  90.          XLAM = XNU*YOUN/((1.0D0+XNU)*(1.0D0-(2.0D0*XNU)))
  91.          XMU = YOUN/(2.0D0*(1.0D0+XNU))
  92. C
  93. C-----CALCUL DE LA MATRICE DE HOOKE EN CONTRAINTES PLANES
  94.          AM = (YOUN)/(1.0D0-(XNU**2.0D0))
  95. C
  96.          D(1,1) = AM
  97.          D(1,2) = AM*XNU  
  98.          D(1,3) = 0.0D0
  99.          D(2,1) = AM*XNU
  100.          D(2,2) = AM
  101.          D(2,3) = 0.0D0
  102.          D(3,1) = 0.0D0
  103.          D(3,2) = 0.0D0
  104.          D(3,3) = AM*(1-XNU)/2.0D0
  105. C
  106. C-----PARAMETRES NON LINEAIRES
  107.          ALP1 = XMAT(5)
  108.          BET1 = XMAT(6)
  109.          AD(1) = XMAT(7)
  110.          AD(2) = XMAT(8)
  111.          AD(3) = XMAT(9)
  112.          AD(4) = XMAT(10)
  113.          XFT = XMAT(11)
  114.          XECR0 = XMAT(12)
  115.          XNIE(1) = XMAT(13)
  116.          XNIE(2) = XMAT(14)
  117.          XUCR = XMAT(15)
  118.          XINDE = INT(XMAT(16))
  119.          XIE(1) = XMAT(17) 
  120.          XIE(2) = XMAT(18)
  121. C
  122. C-------------------------------------------------------------------------------------------
  123. C-----RECUPERATION DES VARIABLES INTERNES
  124. C-------------------------------------------------------------------------------------------    
  125. C-----DEFORMATION TOTALE
  126.          EPSIM(1,1) = VAR0(5)
  127.          EPSIM(2,2) = VAR0(6)
  128.          EPSIM(1,2) = VAR0(7)
  129.          EPSIM(2,1) = VAR0(7)
  130. C-----DEFORMATION TRANSVERSALE
  131.          EPSZZ = VAR0(19)        
  132. C-----VECTEUR TRACTION DANS LE REPERE GLOBAL
  133.          DO I=1,2
  134.              T_U(I) = VAR0(10 + I)
  135.          ENDDO
  136. C-----VECTEUR DISCONTINUTE DANS LE REPERE GLOBAL
  137.          DO I=1,2
  138.              XOUV(I) = VAR0(12 + I)
  139.          ENDDO           
  140. C-----FLAG ACTIVATION DE LA DISCONTINUITE
  141.          XNFLA = VAR0(31)
  142. C-----NORMALE A LA FISSURE ET TRANSPOSEE DE LA NORMALE
  143.          XN(1) = VAR0(33)
  144.          XN(2) = VAR0(34)
  145.          XNT(1,1) = XN(1)
  146.          XNT(1,2) = XN(2)
  147. C-----ECROUISSAGE POUR LA CONDITION DE CONTINUITE
  148.          Z_U = VAR0(36)  
  149. C-----OUVERTURE MAXIMALE DANS LE REPERE DE LA DISCONTINUITE
  150.          OUVN_MAX = VAR0(37)             
  151. C-----OUVERTURE COURANTE DANS LE REPERE DE LA DISCONTINUITE
  152.          OUVN = VAR0(38)
  153.          OUVT = VAR0(39)         
  154. C-----VECTEUR TRACTION DANS LE REPERE DE LA FISSURE
  155.          T_UN = VAR0(43)
  156.          T_UT = VAR0(44)
  157. C-----NORME VP
  158.          XNORM_VP = VAR0(50)
  159. C-----ENERGIE DISSIPEE
  160.          XPSI = VAR0(47)                         
  161. C
  162.          VPT(1) = VAR0(48)
  163.          VPT(2) = VAR0(49)      
  164. C-------------------------------------------------------------------------------------------
  165. C-----PRECALCUL POUR LA GESTION DE LA NORMALE
  166. C-------------------------------------------------------------------------------------------
  167. C
  168. C-----STOCKAGE DES POINTS
  169.          DO J=1,2
  170.              P1(J) = XEL(J,1)
  171.              P2(J) = XEL(J,2)
  172.              P3(J) = XEL(J,3)
  173.          ENDDO
  174. C-----CALCUL DE LA SURFACE DE L'ELEMENT FINI
  175.          XA0 = 0.5D0*((P2(1)*P3(2))-(P3(1)*P2(2))+(P3(1)*P1(2))-
  176.      &   (P1(1)*P3(2))+(P1(1)*P2(2))-(P1(2)*P2(1)))
  177.          XA = ABS(XA0)
  178.          VARF(59) = (SQRT(XA))/2.0D0
  179. C
  180. C-----CALCUL DU CENTRE DE GRAVITE DE L ELEMENT
  181.          XG(1) = (P1(1)+P2(1)+P3(1))/3.0D0
  182.          XG(2) = (P1(2)+P2(2)+P3(2))/3.0D0
  183. C
  184. C-----CALCUL DES VECTEURS GRADIENT DEFINIS POUR CHAQUE NOEUD DE L'ELEMENT
  185.          BNOD1(1,1) = -1.0D0*(P3(2)-P2(2))
  186.          BNOD1(1,2) = P3(1) - P2(1)
  187.          XNORBN1 = SQRT((BNOD1(1,1)**2) + (BNOD1(1,2)**2)) 
  188. C     
  189.          BNOD2(1,1) = -1.0D0*(P1(2)-P3(2))
  190.          BNOD2(1,2) = P1(1)-P3(1)
  191.          XNORBN2 = SQRT((BNOD2(1,1)**2) + (BNOD2(1,2)**2)) 
  192. C
  193.          BNOD3(1,1) = -1.0D0*(P2(2)-P1(2))
  194.          BNOD3(1,2) = P2(1)-P1(1)
  195.          XNORBN3 = SQRT((BNOD3(1,1)**2) + (BNOD3(1,2)**2))
  196. C
  197. C-----CONSTRUCTIONS DE LA MATRICE GRADIENT
  198.          DO J=1,2
  199.              BGRAD(1,J) = BNOD1(1,J)
  200.              BGRAD(2,J) = BNOD2(1,J)
  201.              BGRAD(3,J) = BNOD3(1,J)
  202.          ENDDO
  203. C
  204. C-------------------------------------------------------------------------------------------
  205. C-----ON REALISE LE PREDICTEUR POUR VERIFIER LA CONDITION
  206. C-----D ACTIVATION DE LA DISCONTINUITE
  207. C-------------------------------------------------------------------------------------------
  208. C
  209. C-----ACTUALISATION DE LA DEFORMATION TOTALE
  210.          EPSIM(1,1) = EPSIM(1,1) + DEPST(1)
  211.          EPSIM(2,2) = EPSIM(2,2) + DEPST(2)
  212.          EPSIM(1,2) = EPSIM(1,2) + DEPST(4)/(2.0D0)
  213.          EPSIM(2,1) = EPSIM(2,1) + DEPST(4)/(2.0D0)
  214. C
  215. C-----REECRITURE SOUS LA FORME DE VOIGT
  216.          EPSVOIGT(1,1) = EPSIM(1,1)
  217.          EPSVOIGT(2,1) = EPSIM(2,2)
  218.          EPSVOIGT(3,1) = 2.0D0*EPSIM(1,2)
  219. C
  220. C-----CAS OU LA DISCONTINUITE N EST PAS ACTIVE
  221.          IF (XNFLA.EQ.0.0D0) THEN
  222. C--------ON CALCUL LES CONTRAINTES ELASTIQUES            
  223.              CALL PRDFE1(D,EPSVOIGT,3,3,1,SIGVOI)
  224. C
  225. C--------ON REECRIT SOUS LA FORME MATRICIELLE
  226.              CALL VECMAT(SIGVOI,SIGMA)            
  227. C--------ON INITIALISE LA DEFORMATION TRANSVERSALE
  228.              EPSZZ = -XNU/YOUN*(SIGVOI(1,1)+SIGVOI(2,1))
  229. C--------ON INITIALISE LA CONTRAINTES TRANSVERSALE - POUR STABILITE NUMERIQUE
  230.              SIG_ZZ = XLAM*((EPSIM(1,1)+EPSIM(2,2))+EPSZZ)+(2*XMU*EPSZZ)
  231. C--------INTEGRATION DU NONLOCAL POUR LE CALCUL DE LA NORMALE
  232. C
  233.              IF (ISTEP.EQ.0) THEN 
  234. C
  235.                 EPNL1 = EPSIM(1,1)
  236.                 EPNL2 = EPSIM(2,2)
  237.                 EPNL3 = EPSIM(1,2)
  238. C
  239.                 VARF(51) = EPNL1
  240.                 VARF(52) = EPNL2 
  241.                 VARF(53) = EPNL3 
  242. C
  243. C
  244.              ELSEIF (ISTEP.EQ.1) THEN 
  245. C
  246.                 VARF(51) = EPSIM(1,1)
  247.                 VARF(52) = EPSIM(2,2) 
  248.                 VARF(53) = EPSIM(1,2)
  249. C
  250.                 DO I=1,50 
  251. C
  252.                   VARF(I)=VAR0(I) 
  253. C
  254.                 ENDDO 
  255. C
  256.              GOTO 2000 
  257. C
  258.              ELSEIF (ISTEP.EQ.2) THEN 
  259. C
  260.                 EPNL1 =VAR0(51)
  261.                 EPNL2 =VAR0(52) 
  262.                 EPNL3 =VAR0(53)  
  263. C
  264.                 VARF(51) = EPNL1
  265.                 VARF(52) = EPNL2 
  266.                 VARF(53) = EPNL3
  267. C
  268.              ENDIF 
  269. C
  270.              EPNL(1,1) = EPNL1
  271.              EPNL(2,2) = EPNL2
  272.              EPNL(1,2) = EPNL3
  273.              EPNL(2,1) = EPNL3
  274. C
  275. C--------CALCUL VALEUR PROPRES ET VECTEURS PROPRES
  276.              CALL DIAFE1(EPNL,VAP,VEP)
  277. C
  278. C--------IDENTIFICATION DE LA CONTRAINTE PRINCIPALE MAJEURE ET DE SON INDICE
  279.              EPSI_MAX= MAX(VAP(1),VAP(2))
  280. C
  281.              DO JJ=1,2
  282.                  IF (VAP(JJ).EQ.EPSI_MAX) THEN
  283.                      INDI = JJ
  284.                  ENDIF
  285.              ENDDO
  286. C
  287. C--------IDENTIFICATION DU VECTEUR PROPRE
  288. C
  289.              DO J=1,2
  290.                 DP_MAX(J) = VEP(J,INDI)
  291.              ENDDO
  292. C
  293.              XNDP = SQRT((DP_MAX(1)**2.0D0)+(DP_MAX(2)**2.0D0))
  294. C 
  295.              XN(1) = XNIE(1)
  296.              XN(2) = XNIE(2)
  297. C
  298.              CALL TRSPOR(XN,1,2,XNT)
  299. C             
  300. C--------MATRICE DE LA NORMALE - SIMPLICITE DE CALCUL
  301.              N(1,1) = XN(1)  
  302.              N(1,2) = 0.0D0
  303.              N(2,1) = 0.0D0   
  304.              N(2,2) = XN(2)
  305.              N(3,1) = XN(2) 
  306.              N(3,2) = XN(1)             
  307. C
  308. C--------TRANSPOSE DE LA MATRICE DE LA NORMALE
  309.              DO II=1,2
  310.                  DO JJ=1,3
  311.                      M(II,JJ) = N(JJ,II)
  312.                  ENDDO
  313.              ENDDO             
  314. C
  315. C--------VECTEUR TRACTION T = SIGMA.N POUR L INITIATION DE LA DISCONTINUITE
  316.              CALL PRDFE1(M,SIGVOI,2,3,1,T)
  317.              CALL PRDFE1(XNT,T,1,2,1,T_CN0)             
  318. C
  319.              VP(1,1) = 0.0D0
  320.              VP(1,2) = 0.0D0
  321. C
  322.              IF (XINDE.EQ.1) THEN
  323. C
  324.                 FP1 = (1.0D0*XNIE(1)*P1(1))+(XNIE(2)*P1(2))-
  325.      &          ((1.0D0*XNIE(1)*XIE(1))+(XNIE(2)*XIE(2)))
  326. C                 
  327.                 FP2 = (1.0D0*XNIE(1)*P2(1))+(XNIE(2)*P2(2))-
  328.      &          ((1.0D0*XNIE(1)*XIE(1))+(XNIE(2)*XIE(2)))
  329. C
  330.                 FP3 = (1.0D0*XNIE(1)*P3(1))+(XNIE(2)*P3(2))-
  331.      &          ((1.0D0*XNIE(1)*XIE(1))+(XNIE(2)*XIE(2)))     
  332. C
  333. C           2 - VERIFICATION DES SIGNES
  334.                 CALL SIGEK(FP1,F1)
  335.                 CALL SIGEK(FP2,F2)
  336.                 CALL SIGEK(FP3,F3)
  337. C                
  338.                 FF(1) = F1
  339.                 FF(2) = F2
  340.                 FF(3) = F3
  341. C                
  342.                 DO II = 1,3
  343.                     IF (FF(II).GT.0.0D0) THEN
  344.                         DO J = 1,2
  345.                             VP(1,J) = VP(1,J) + BGRAD(II,J)
  346.                         ENDDO
  347.                     ENDIF
  348.                 ENDDO
  349. C            
  350. C--------CONDITION D INITIATION
  351.                 IF ((T_CN0(1,1).GT.XFT)) THEN
  352.                     T_UN = XFT
  353.                     T_UT = 0                
  354.                     OUVN = 1.0D-7
  355.                     VAR0(38) = OUVN
  356.                     OUVT = 1.0D-12
  357.                     VAR0(39) = OUVT
  358.                     Z_U = 0.0D0
  359.                     XNFLA = 1.0D0
  360.                ELSE
  361.                     XNFLA = 0.0D0                        
  362.                ENDIF                 
  363.             ENDIF
  364. C
  365.              CALL TRSPOR(VP,2,1,VPT)
  366.              XNORM_VP = SQRT((VPT(1)**2)+(VPT(2)**2))
  367. C             
  368.         ENDIF
  369. C--------FIN DE LA PREDICTION ELASTIQUE NORMALE
  370. C
  371. C--------DEFINITION DE LA MATRICE G
  372.          G(1,1) = -1.D0*VPT(1)/(2.0D0*XA)
  373.          G(1,2) = 0.0D0/(2.0D0*XA)
  374.          G(2,1) = 0.0D0/(2.0D0*XA)
  375.          G(2,2) = -1.D0*VPT(2)/(2.0D0*XA)
  376.          G(3,1) = -1.D0*VPT(2)/(2.0D0*XA)
  377.          G(3,2) = -1.D0*VPT(1)/(2.0D0*XA)
  378. C
  379. C--------CALCUL DE DERIVE T_CN
  380.          CALL PRDFE1(D,G,3,3,2,XDG)
  381.          CALL PRDFE1(XDG,XN,3,2,1,XDGN)
  382.          CALL VECMAT(XDGN,F)
  383. C
  384. C--------CALCUL DU VECTEUR TANGENT A LA FISSURE
  385.          XT(1) = -XN(2)
  386.          XT(2) = XN(1)
  387. C
  388. C         CALL TRSPOR(XT,1,2,XTT)
  389.  
  390.          XTT(1,1) = XT(1)
  391.          XTT(1,2) = XT(2)
  392.  
  393. C         print*,'XTT=',XTT
  394. C         print*,'XT=',XT
  395.  
  396. C         stop
  397.  
  398. C-------------------------------------------------------------------------------------------
  399. C-----ON REALISE DE NOUVEAU UNE ETAPE DE PREDICTION ELASTIQUE EN PRENANT EN COMPTE 
  400. C-----UNE EVENTUELLE DISCONTINUITE
  401. C-------------------------------------------------------------------------------------------
  402. C
  403. C-----CALCUL DU TERME G X [U]
  404.          CALL PRDFE1(G,XOUV,3,2,1,EPSD)
  405. C
  406. C-----CALCUL DU TERME EPSILON + G X [U]
  407.          DO II=1,3
  408.              EPSENR_V(II) = EPSVOIGT(II,1) + EPSD(II)
  409.          ENDDO
  410. C
  411. C-----REECRITURE EN MATRCIEL
  412.          CALL VECDEF(EPSENR_V,EPSENR)
  413. C
  414. C-----CALCUL DES CONTRAINTE DANS LE CONTINUUM
  415.          DO II=1,2
  416.              DO JJ=1,2
  417.                  SIG(II,JJ) = XLAM*(EPSENR(1,1)+EPSENR(2,2)+
  418.      &          EPSZZ)*UNIT(II,JJ)+
  419.      &          2.0D0*XMU*EPSENR(II,JJ)
  420.              ENDDO
  421.          ENDDO
  422. C
  423. C-----CALCUL DU VECTEUR TRACTION
  424.          CALL PRDFE1(SIG,XN,2,2,1,T_CN1)
  425.          CALL PRDFE1(XNT,T_CN1,1,2,1,T_CN)
  426.          CALL PRDFE1(XTT,T_CN1,1,2,1,T_CT)
  427.  
  428. C         CALL PRODMA(SIG,XN,2,2,1,T_CN1)
  429. C         CALL PRODMA(XNT,T_CN1,1,2,1,T_CN)
  430. C         CALL PRODMA(XTT,T_CN1,1,2,1,T_CT)
  431.  
  432. C         print*,'XNT=',XNT
  433. C         print*,'XTT=',XTT
  434.  
  435. C         print*,'T_CN=',T_CN(1)
  436. C         print*,'T_CT=',T_CT(1)
  437. C         STOP
  438. C
  439. C-------------------------------------------------------------------------------------------
  440. C-----BOUCLE SUR LA CONDITION DE CONTRAINTE PLANES - START
  441. C-------------------------------------------------------------------------------------------
  442. C
  443.          DO II=1,ITMAX1
  444. C
  445. C--------SI PAS DE DISCONTINUITE ON RENVOIE LA SOLUTION ELASTIQUE
  446.              IF (XNFLA.EQ.1.0D0) THEN
  447. C-----------INITIALISATION DU RESIDU SUR LA CONTINUITE
  448.                  RES(1) = T_CN(1) - T_UN
  449.                  RES(2) = T_CT(1) - T_UT
  450. C
  451. C-------------------------------------------------------------------------------------------
  452. C-------------BOUCLE SUR LA CONDITION DE CONTINUITE - START
  453. C-------------------------------------------------------------------------------------------
  454.                      DO III=1,ITMAX2                     
  455. C                        
  456. C----------------CALCUL DE LA CORRECTION DE DISCONTINUITE
  457. C                       
  458.                         DO JJ = 1,2
  459.                         DO KK = 1,2
  460.                             K(JJ,KK) = 1.0D0*YOUN*UNIT(JJ,KK)/
  461.      &                                               SQRT(2.0D0*XA)
  462.                         ENDDO
  463.                         ENDDO
  464.                             K(1,2) = 1.0D0*YOUN / 10.0D0
  465.                             K(2,1) = 1.0D0*YOUN / 10.0D0
  466.                         CALL HINV22(K,KINV)
  467.                         CALL PRDFE1(KINV,RES,2,2,1,DUN)                       
  468. C                        
  469.                         DUNN = DUN(1)
  470.                         DUNT = DUN(2)
  471. C
  472. C----------------ACTUALISATION DE LA DISCONTINUITE (RELAXATION NUMERIQUE)
  473.                          OUVN = OUVN + (1.0D-3*DUNN)
  474.                          OUVT = OUVT + (1.0D-3*DUNT)
  475. C
  476. C----------------CHARGE OU DECHARGE (POUR LA CONDITION DE CONTINUITE)
  477.                          PHI_U = 1.0D0
  478. C
  479. C
  480. C----------------CONDITION D UNILATERALITE
  481.                          IF (OUVN.LT.0.0D0) THEN
  482.                              OUVN = 1.0D-12                                  
  483.                              GOTO 9
  484.                          ENDIF      
  485. C
  486. C----------------ACTUALISATION DE L OUVERTURE NORMALE MAXIMALE
  487.                          OUVN_MAX =  MAX(OUVN,VAR0(37))
  488. C
  489. C----------------SWITCH VERS LA LOI DISCRETE CONSIDEREE
  490. C                         
  491. C                       3. ECROUISSAGE EXPONENTIEL
  492.                          IF (NECR.EQ.3) THEN
  493. C
  494.                            IF (PHI_U.GT.0.0D0) THEN
  495. C
  496.                IF ((II.EQ.1).AND.(III.EQ.1)) THEN
  497.                   VMOUV0 = SQRT((VAR0(38)**2.0D0)+(VAR0(39)**2.0D0))
  498.                   VMOUV = VMOUV0
  499.                ELSE
  500. C
  501.                  VMOUV = SQRT((OUVN**2.0D0)+(OUVT**2.0D0))
  502.                ENDIF                  
  503.               T_UN = XFT*(OUVN/VMOUV)*EXP(-1.0D0*VMOUV/XKAPPA*XFT)
  504.               T_UT = XFT*(OUVT/VMOUV)*EXP(-1.0D0*VMOUV/XKAPPA*XFT)
  505.  
  506.               Z_U = T_UN
  507.                              ELSE
  508. C
  509.                                  T_UN = (Z_U/OUVN_MAX)*OUVN
  510. C                       
  511.                              ENDIF
  512.                          ENDIF 
  513. C
  514. C----------------CALCUL DU VECTEUR DISCONTINUITE
  515.                          DO JJ=1,2
  516.                             XOUV(JJ) = (OUVN*XN(JJ))+(OUVT*XT(JJ))
  517.                          ENDDO
  518. C
  519. C----------------CALCUL DU TERME G X [U]
  520.                          CALL PRDFE1(G,XOUV,3,2,1,EPSD)
  521. C      
  522. C----------------CALCUL DU TERME EPSILON + G X [U]
  523.                          DO JJ=1,3
  524.                              EPSENR_V(JJ) = EPSVOIGT(JJ,1) + EPSD(JJ)
  525.                          ENDDO
  526. C
  527. C----------------REECRITURE EN MATRCIEL
  528.                          CALL VECDEF(EPSENR_V,EPSENR)
  529. C
  530. C----------------CALCUL DES CONTRAINTE DANS LE CONTINUUM
  531.                          DO JJ=1,2
  532.                              DO KK=1,2
  533.                                  SIG(JJ,KK) = XLAM*((EPSENR(1,1)+
  534.      &                          EPSENR(2,2))+ EPSZZ)*UNIT(JJ,KK)+
  535.      &                          2.0D0*XMU*EPSENR(JJ,KK)
  536.                              ENDDO
  537.                          ENDDO
  538. C
  539. C----------------CALCUL DU VECTEUR CONTRAINTE
  540.                          CALL PRDFE1(SIG,XN,2,2,1,T_CN1)
  541.                          CALL PRDFE1(XNT,T_CN1,1,2,1,T_CN)
  542.                          CALL PRDFE1(XTT,T_CN1,1,2,1,T_CT)
  543.  
  544. C                         print*,'T_CN=',T_CN(1)
  545. C                         print*,'T_CT=',T_CT(1)
  546. C                         STOP
  547. C
  548. C----------------CALCUL DU RESIDU ACTUALISE
  549.                       RES(1) = T_CN(1) - T_UN
  550.                       RES(2) = T_CT(1) - T_UT
  551. C
  552. C
  553.                       XNORMRES = SQRT((RES(1)*RES(1))+(RES(2)*RES(2)))
  554. C
  555. C----------------CALCUL DE L ENERGIE DISSIPEE
  556.                          XPSI = T_UN*OUVN
  557. C
  558. C----------------TEST DE CONVERGENCE
  559.                          IF ((XNORMRES).LT.TOL2) THEN
  560. C
  561. C-------------------CALCUL DE L OUVERTURE DANS LE REPERE GLOBAL
  562.                              DO JJ=1,2
  563.                                XOUV(JJ) = (OUVN*XN(JJ))+(OUVT*XT(JJ))
  564.                              ENDDO  
  565. C
  566. C-------------------CALCUL DES VECTEURS TRACTIONS DANS LE CONTINUUM ET
  567. C-------------------SUR LA DISCONTINUITE
  568.                              DO JJ=1,2
  569.                        T_C(JJ) = (T_CN(1)*XN(JJ))+(T_CT(1)*XT(JJ))
  570.                        T_U(JJ) = (T_UN*XN(JJ))+(T_UT*XT(JJ))
  571.                              ENDDO 
  572. C-------------------BALISE DE SORTIE
  573.                              GOTO 9
  574. C
  575.                          ENDIF
  576. C
  577.                      ENDDO
  578. C
  579. C-------------------------------------------------------------------------------------------
  580. C-------------BOUCLE SUR LA CONDITION DE CONTINUITE - END
  581. C-------------------------------------------------------------------------------------------
  582.                  ENDIF
  583. C
  584.  9              CONTINUE
  585. C
  586. C
  587. C-----------CALCUL DES CONTRAINTES DE CAUCHY A CONVERGENCE
  588.                  DO JJ=1,2
  589.                     XOUV(JJ) = (OUVN*XN(JJ))+(OUVT*XT(JJ))
  590.                  ENDDO
  591. C
  592.                  CALL PRDFE1(G,XOUV,3,2,1,EPSD)
  593. C      
  594.                  DO JJ=1,3
  595.                      EPSENR_V(JJ) = EPSVOIGT(JJ,1) + EPSD(JJ)
  596.                   ENDDO
  597. C
  598.                  CALL VECDEF(EPSENR_V,EPSENR)
  599. C
  600.                  DO JJ=1,2
  601.                      DO KK=1,2
  602.                          SIG(JJ,KK) = XLAM*((EPSENR(1,1)+
  603.      &                  EPSENR(2,2))+ EPSZZ)*UNIT(JJ,KK)+
  604.      &                  2.0D0*XMU*EPSENR(JJ,KK)
  605.                      ENDDO
  606.                  ENDDO
  607. C
  608. C-----------CALCUL CONTRAINTES HORS PLAN
  609.                  SIG_ZZ = XLAM*((EPSENR(1,1)+EPSENR(2,2))+
  610.      &          EPSZZ)+(2*XMU*EPSZZ)
  611. C
  612. C-----------TEST DE CONVERGENCE SUR LA CONDITION DE CONTRAINTES PLANES
  613.                  IF (ABS(SIG_ZZ).LT.TOL1) THEN
  614. C
  615.                      GOTO 99
  616. C  
  617.                  ELSE
  618. C
  619.                      EPSZZ = EPSZZ - (SIG_ZZ/(XLAM+(2*XMU)))
  620.                  ENDIF
  621. C
  622. C
  623.          ENDDO
  624. C-------------------------------------------------------------------------------------------
  625. C-----BOUCLE SUR LA CONDITION DE CONTRAINTE PLANES - END
  626. C-------------------------------------------------------------------------------------------
  627. C
  628.  99     CONTINUE
  629. C
  630. C-------------------------------------------------------------------------------------------
  631. C-----CONTRAINTES DE CAUCHY EN SORTIE POUR CALCUL DES EFFORTS INTERNES
  632. C-------------------------------------------------------------------------------------------
  633.          IF (XNFLA.EQ.0.0D0) THEN
  634.              CALL VECMAT(SIGVOI,SIGMAM)
  635.          ELSE                    
  636.              DO JJ=1,2
  637.              DO KK=1,2   
  638.              SIGMAM(JJ,KK) = SIG(JJ,KK)
  639.              ENDDO
  640.              ENDDO 
  641.          ENDIF
  642. C
  643.          SIGF(1) = SIGMAM(1,1)
  644.          SIGF(2) = SIGMAM(2,2)
  645.          SIGF(3) = 0.0D0
  646.          SIGF(4) = SIGMAM(1,2)
  647. C
  648. C
  649. C-------------------------------------------------------------------------------------------
  650. C-----VARIABLES INTERNES FINALES
  651. C-------------------------------------------------------------------------------------------
  652. C-----VECTEUR TRACTION CALCULES PAR LA LOI DISCRETE DANS LE REPERE GLOBAL
  653.          DO I=1,2
  654.              VARF(10+I) = T_U(I)
  655.          ENDDO
  656. C
  657. C-----VECTEUR DISCONTINUITE DANS LE REPERE GLOBAL
  658.          DO I=1,2
  659.              VARF(12+I) = XOUV(I)
  660.          ENDDO
  661. C
  662. C-----DEFORMATIONS TOTALES
  663.          VARF(5) = EPSIM(1,1)
  664.          VARF(6) = EPSIM(2,2)     
  665.          VARF(7) = EPSIM(1,2)
  666. C-----DEFORMATION TRANSVERSALE
  667.          VARF(19) = EPSZZ
  668. C-----FLAG ACTIVATION DISCONTINUITE
  669.          VARF(31) = XNFLA
  670. C-----VECTEUR DE LA NORMALE DANS LE REPERE GLOBAL
  671.          VARF(33) = XN(1)
  672.          VARF(34) = XN(2)
  673. C-----ECROUISSAGE LIE A LA LOI DISCRETE
  674.          VARF(36) = Z_U                
  675. C-----OUVERTURE NORMALE MAXIMALE
  676.          VARF(37) = OUVN_MAX               
  677. C-----OUVERTURE NORMALE COURANTE
  678.          VARF(38) = OUVN
  679.          VARF(39) = OUVT         
  680. C-----VECTEUR TRACTION DISCRET
  681.          VARF(43) = T_UN
  682.          VARF(44) = T_UT                 
  683. C-----NORME DE VP
  684.          VARF(50) = XNORM_VP
  685.          IF (XNFLA.EQ.1.0D0) THEN
  686.          VARF(45) = XN(1)
  687.          VARF(46) = XN(2)
  688.          ENDIF
  689. C
  690. C-----ENERGIE DISSIPEE 
  691. C
  692.          VARF(47) = XPSI         
  693.          VARF(48) = VPT(1)
  694.          VARF(49) = VPT(2)       
  695. C-----CONTRAINTES                
  696.          VARF(54) = SIGF(1)
  697.          VARF(55) = SIGF(2)
  698.          VARF(56) = SIGF(3)
  699.          VARF(57) = SIGF(4)              
  700.  2000    CONTINUE
  701. C
  702.          RETURN
  703.          END
  704.  
  705.  
  706.  
  707.  
  708.  

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