Télécharger la plaquette Cast3M


fr en

Télécharger sicscal.eso

Retour à la liste

Numérotation des lignes : 

sicscal

  1. C SICSCAL   SOURCE    OF166741  25/11/04    21:16:06     12349          
  2.       SUBROUTINE SICSCAL(WRK52,WRK53,WRK54,WRK22,IB,IGAU,NVARI,NBPGAU
  3.      & ,necou,iecou)
  4. C     Modèle ONERA Scalaire pour le SICf/SiC tissé
  5. C
  6. C     variables en entree
  7. C
  8. C     WRK0,KRK1,WRK5  pointeurs sur des segments de travail
  9. C
  10. C
  11. C     WRK0:
  12. C            XMAT(NMATT): tableau composantes materiaux
  13. C
  14. C     WRK52:
  15. C                 Dans le repère global
  16. C            SIG0(NSTRS) : contraintes au debut du pas
  17. C            SIGF(NSTR) : CONTRAINTES à l afin du pas
  18. C            VAR0(NVARI) : variables internes au debut du pas
  19. C            DEPST(NSTRS): increment de deformation totale
  20. C            EPST0(NSTRS): deformation totale au début
  21. C
  22. C            EPSTF(NSTR): déformations totales finales
  23. C
  24. C     WRK22:
  25. C            XXE: coordonnees de l'element en double precison
  26. C
  27. C     WTRAV:
  28. C            TXR: cosinus directeurs des axes locaux de l'element massif
  29. C
  30. C
  31. C     WRK5:
  32. C
  33. C
  34. C
  35. C     NSTRS      nombre de composantes dans les vecteurs des contraintes
  36. C                                        et les vecteurs des deformations
  37. C
  38. C     NVARI      nombre de variables internes (doit etre egal a 3)
  39. C
  40. C     NMATT      nombre de constantes du materiau
  41. C
  42. C
  43. C
  44. C     variables en sortie
  45. C
  46. C     VARF      variables internes finales dans WRK0
  47. C
  48. C     SIGF      contraintes finales dans WRK0
  49. C
  50.       IMPLICIT INTEGER(I-N)
  51.       IMPLICIT REAL*8(A-H,O-Z)
  52.  
  53. -INC PPARAM
  54. -INC CCOPTIO
  55. -INC CCREEL
  56. -INC DECHE
  57. -INC TECOU
  58.  
  59.       SEGMENT WRK22
  60.         REAL*8 XXE(3,NBNNbi)
  61.       ENDSEGMENT
  62.  
  63.       DIMENSION ZDAUX(6,6)
  64.       DIMENSION H1 (6,6), H2(6,6), H3(6,6)
  65.       DIMENSION H01 (6,6), H02(6,6), H03(6,6)
  66. *     MATRICE DES SOUPLESSES INITIALES
  67.       DIMENSION SEFF(6,6),CEFF(6,6)
  68. *      sortH CONTRAINTES DANS LE REPÈRE ORTHOTROPE PUIS NON
  69. *      DORTH :  DEFO DAN S LE REPÈRE OrTHO
  70. *      sigfv : CONTRAINTE RÉORGANISÉES ORDRE CASTEM
  71. *
  72. *      EPSIN : Défo inel dans le repère ortho ordre onera
  73. *      puis dans le repère glob aordre castem
  74. *      epsel : DEFO ELASTIQUE CALCULÉE POUR L'OBTENTION DES CONTRAINTES
  75. *
  76.       DIMENSION SORTH(6),  SIGFV(6), EPSIN(6),EPSINV(6),
  77.      &        DORTH(6), EPSEL(6),EPSTFV(6),
  78.      &       SORTH1(6),SORTH2(6)
  79. *     defoth : DEFO DANS LE REPERE ORTHO SOus FORME (3X3)
  80. *       POUR LE CALCUL DES DEFO PROPRES PUIS PARTIE >0 DES DÉFO PROPRES
  81. **      xprop : DIRECTIONS PROPRES/REPERE ortho
  82. *       xpinv : MATRICE INVERSE DE CHANGEMENT DE REPERE XPROP
  83. *      xPINVT  : MATRICE TRANSPOSÉE DE XPINV
  84. *      defotht : MATRICE TEMPORAIRE POUR LE CHANGEMENT DE REPÈRE
  85. *      epsclo : DEFO DE FERMETURE (CLOSURE)
  86. *      EPSBAR : DEFO - DEFO CLOSURE
  87.       DIMENSION DEFOTH(3,3),PROP(3),XPROP(3,3),XPINV(3,3),
  88.      &       DEFOTHT(3,3),XPINVT(3,3),EPSCLO(6),
  89.      &       EPSBAR(6),EPSIM(6)
  90.       DIMENSION ZK1(6,6),ZK2(6,6),ZK3(6,6),TERME(6,6)
  91.       DIMENSION EPSS1(6),EPSR1(6)
  92. *     Déformation résiduelles et stockées
  93.       DIMENSION EPSS(6),EPSR(6)
  94. *     DRIVING FORCES NORMALES ET TRANSVERSALES
  95. *     dom : VARIABLES D'ENDOMAGEMENT
  96.       DIMENSION YN(3),YT(3),DOM(3)
  97. *     DEPSIF (DELTA EPSILON I f :
  98. *     LA FERMETURE COMMENCE POUR epsbar = DELTA EPSILON
  99. *     finit POUR  epsbar = DELTA EPSILON
  100. *     ZNUACT  INDICE D'ACTIVATION =1 : ACTIVE
  101. *                                =0 : PASSIVE
  102.       DIMENSION DEPSIF(3),ZNUACT(3)
  103. * Pseudo contraintes pour le calcul des defo residuelles EPSR et stockées EPSS
  104.       DIMENSION SOMS(6,6),SOMR(6,6),
  105.      $ SISIS(6),SISIR(6)
  106. * Tableau de paramètres
  107.       DIMENSION Y0N(3),Y0C(3),Y0T(3),AIF(3),PT(3),PN(3),
  108.      $ YCT(3),DCT(3),DCN(3),DEPSIF0(3),YCN(3)
  109.       DATA ( ( H1(I,J), J=1,6), I=1,6)
  110.      &    /1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  111.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  112.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  113.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  114.      &     0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
  115.      &     0.0, 0.0, 0.0, 0.0, 0.0, 1.0/
  116.       DATA ( ( H2(I,J), J=1,6), I=1,6)
  117.      &    /0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  118.      &     0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
  119.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  120.      &     0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
  121.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  122.      &     0.0, 0.0, 0.0, 0.0, 0.0, 1.0/
  123. *
  124.       DATA ( ( H3(I,J), J=1,6), I=1,6)
  125.      &    /0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  126.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  127.      &     0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
  128.      &     0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
  129.      &     0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
  130.      &     0.0, 0.0, 0.0, 0.0, 0.0, 0.0/
  131.  
  132. * pour le modèle scalaire
  133. * Attention on ne rentre  pas les racines carrée des valeurs
  134. *  Y0 et YC telles que définies  formules (2.21) et (2.24)
  135. *  dans SEMT/LM2S/RT/05-034
  136. *    PHENOMENOLOGICAL COEFFICIENTS
  137. *
  138. * scalar damage law parameters
  139. * G1DC, G1Y0, G1YC, G1P
  140. * G2DC, G2Y0, G2YC, G2P
  141. * G3DC, G3Y0, G3YC, G3P
  142. * paramètre pour le calcul des déformations résiduelles
  143. * modèle scalaire (5 en pseudo tensoriel)
  144. ***** INITIALIZATION VARIABLES
  145.       lhookl = LHOOK
  146. * Proprietes materiau
  147.        YG1   = XMAT(1)
  148.        YG2   = XMAT(2)
  149.        YG3   = XMAT(3)
  150.        XNU12 = XMAT(4)
  151.        XNU23 = XMAT(5)
  152.        XNU13 = XMAT(6)
  153.        G12   = XMAT(7)
  154.        G23   = XMAT(8)
  155.        G13   = XMAT(9)
  156. * Nombre de param obligatoires pour le modèle 15 + ALP1 ALP2 ALP3 et RHO =19
  157.        NOBL = 19
  158. * On extrait les valeurs facultatives du modèle
  159.       H1N = XMAT(NOBL + 1)
  160.       IF(IVALMA(NOBL+1).EQ.0.) H1N = 1.
  161.       H1HP = XMAT(NOBL + 2)
  162.       IF(IVALMA(NOBL+2).EQ.0.) H1HP = 0.7
  163.       H1P = XMAT(NOBL + 3)
  164.       IF(IVALMA(NOBL+3).EQ.0.) H1P = 0.45
  165.       H2N = XMAT(NOBL + 4)
  166.       IF(IVALMA(NOBL+4).EQ.0.) H2N = 1.
  167.       H2HP = XMAT(NOBL + 5)
  168.       IF(IVALMA(NOBL+5).EQ.0.) H2HP = 0.7
  169.       H2P = XMAT(NOBL + 6)
  170.       IF(IVALMA(NOBL+6).EQ.0.) H2P = 0.45
  171.       H3N = XMAT(NOBL+7)
  172.       IF(IVALMA(NOBL+7).EQ.0.) H3N = 1.
  173.       H3P = XMAT(NOBL + 8)
  174.       IF(IVALMA(NOBL+8).EQ.0.) H3P = 0.7
  175.       Y0N(1) = XMAT(NOBL + 9)
  176.       IF(IVALMA(NOBL+9).EQ.0.) Y0N(1) =173.205
  177.       Y0N(2) = XMAT(NOBL + 10)
  178.       IF(IVALMA(NOBL+10).EQ.0.) Y0N(2) =173.205
  179.       Y0N(3) = XMAT(NOBL + 11)
  180.       IF(IVALMA(NOBL+11).EQ.0.) Y0N(3) =173.205
  181.       YCN(1) = XMAT(NOBL + 12)
  182.       IF(IVALMA(NOBL+12).EQ.0.) YCN(1) =1870.83
  183.       YCN(2) = XMAT(NOBL + 13)
  184.       IF(IVALMA(NOBL+13).EQ.0.) YCN(2) =1870.83
  185.       YCN(3) = XMAT(NOBL + 14)
  186.       IF(IVALMA(NOBL+14).EQ.0.) YCN(3) =1870.83
  187.       Y0T(1) = XMAT(NOBL + 15)
  188.       IF(IVALMA(NOBL+15).EQ.0.) Y0T(1) =31.623
  189.       Y0T(2) = XMAT(NOBL + 16)
  190.       IF(IVALMA(NOBL+16).EQ.0.) Y0T(2) =31.623
  191.       Y0T(3) = XMAT(NOBL + 17)
  192.       IF(IVALMA(NOBL+17).EQ.0.) Y0T(3) =173.205
  193.       YCT(1) = XMAT(NOBL + 18)
  194.       IF(IVALMA(NOBL+18).EQ.0.) YCT(1) =1870.83
  195.       YCT(2) = XMAT(NOBL + 19)
  196.       IF(IVALMA(NOBL+19).EQ.0.) YCT(2) =1870.83
  197.       YCT(3) = XMAT(NOBL + 20)
  198.       IF(IVALMA(NOBL+20).EQ.0.) YCT(3) =1870.83
  199.       DCT(1) = XMAT(NOBL + 21)
  200.       IF(IVALMA(NOBL+21).EQ.0.) DCT(1) =4.
  201.       DCT(2) = XMAT(NOBL + 22)
  202.       IF(IVALMA(NOBL+22).EQ.0.) DCT(2) =4.
  203.       DCT(3) = XMAT(NOBL + 23)
  204.       IF(IVALMA(NOBL+23).EQ.0.) DCT(3) =4.
  205.       DCN(1) = XMAT(NOBL + 24)
  206.       IF(IVALMA(NOBL+24).EQ.0.) DCN(1) =4.
  207.       DCN(2) = XMAT(NOBL + 25)
  208.       IF(IVALMA(NOBL+25).EQ.0.) DCN(2) =4.
  209.       DCN(3) = XMAT(NOBL + 26)
  210.       IF(IVALMA(NOBL+26).EQ.0.) DCN(3) =4.
  211.       PN(1) = XMAT(NOBL + 27)
  212.       IF(IVALMA(NOBL+27).EQ.0.) PN(1) =1.
  213.       PN(2) = XMAT(NOBL + 28)
  214.       IF(IVALMA(NOBL+28).EQ.0.) PN(2) =1.
  215.       PN(3) = XMAT(NOBL + 29)
  216.       IF(IVALMA(NOBL+29).EQ.0.) PN(3) =1.
  217.       PT(1) = XMAT(NOBL + 30)
  218.       IF(IVALMA(NOBL+30).EQ.0.) PT(1) =1.2
  219.       PT(2) = XMAT(NOBL + 31)
  220.       IF(IVALMA(NOBL+31).EQ.0.) PT(2) =1.2
  221.       PT(3) = XMAT(NOBL + 32)
  222.       IF(IVALMA(NOBL+32).EQ.0.) PT(3) =1.
  223.       B = XMAT(NOBL + 33)
  224.       IF(IVALMA(NOBL+33).EQ.0.) B =1.
  225.       DELTAL = XMAT(NOBL + 34)
  226.       IF(IVALMA(NOBL+34).EQ.0.) DELTAL =0.
  227.       TZERO = XMAT(NOBL + 35)
  228.       IF(IVALMA(NOBL+35).EQ.0.) TZERO =0.
  229.       DEPSIF0(1) = XMAT(NOBL + 36)
  230.       IF(IVALMA(NOBL+36).EQ.0.) DEPSIF0(1) =0.0003
  231.       DEPSIF0(2) = XMAT(NOBL + 37)
  232.       IF(IVALMA(NOBL+37).EQ.0.) DEPSIF0(2) =0.0003
  233.       DEPSIF0(3) = XMAT(NOBL + 38)
  234.       IF(IVALMA(NOBL+38).EQ.0.) DEPSIF0(3) =0.0003
  235.       AIF(1) = XMAT(NOBL + 39)
  236.       IF(IVALMA(NOBL+39).EQ.0.) AIF(1) =0.5
  237.       AIF(2) = XMAT(NOBL + 40)
  238.       IF(IVALMA(NOBL+40).EQ.0.) AIF(2) =0.5
  239.       AIF(3) = XMAT(NOBL + 41)
  240.       IF(IVALMA(NOBL+41).EQ.0.) AIF(3) =0.5
  241.       ETA1 = XMAT(NOBL + 42)
  242.       IF(IVALMA(NOBL+42).EQ.0.) ETA1  =0.1
  243.       ETA2 = XMAT(NOBL + 43)
  244.       IF(IVALMA(NOBL+43).EQ.0.) ETA2  =0.1
  245.       ETA3 = XMAT(NOBL + 44)
  246.       IF(IVALMA(NOBL+44).EQ.0.) ETA3 =0.0
  247. *
  248. * Pour debuggage
  249.       PARAM=0
  250. * Variables d'endommagement
  251.       DOM1=VAR0(2)
  252.       DOM2=VAR0(3)
  253.       DOM3=VAR0(4)
  254. ****** CALCUL DE LA DEFORMATION INITIALE A PARTIR
  255. ****** DES CONTRAINTES INITIALES
  256. **********Debuggage
  257.       IF (PARAM.EQ.1) THEN
  258.          WRITE(IOIMP,66770) IB,IGAU
  259. 66770    format(////////2x,'element  ',i6,2x,'point  ',i3//)
  260.  
  261.          WRITE (IOIMP,*) 'Increment des deformations '
  262.          WRITE (IOIMP,99999) (DEPST(ILOOP), ILOOP=1,6)
  263.  
  264.          WRITE (IOIMP,*) 'Contraintes au debut de l''iteration'
  265.          WRITE (IOIMP,99999) (SIG0(ILOOP), ILOOP=1,6)
  266.       ENDIF
  267. **********Debuggage
  268. * Controle si il faut calculer la matrice de hook
  269. *      PRINT *,'IB,IGAU,N2PTEL,N2EL',IB,IGAU,N2PTEL,N2EL
  270.       IF (IB.EQ.1.AND.IGAU.EQ.1) THEN
  271.          GOTO 100
  272.       ELSEIF (N2PTEL.EQ.1.AND.N2EL.EQ.1) THEN
  273.          GOTO 200
  274.       ELSEIF (N2PTEL.EQ.1.AND.N2EL.NE.1)  THEN
  275.          IF (IGAU.EQ.1) THEN
  276.             GOTO 100
  277.          ELSE
  278.             GOTO 200
  279.          ENDIF
  280.       ELSE
  281.          GOTO 100
  282.       ENDIF
  283. 100   CONTINUE
  284. * Calcul de la matrice de hook pour un materiau orthotrope: C0
  285.       IPERR=1
  286.       CALL ZERO (DDAUX,LHOOKL,LHOOKL)
  287.       XNU21=(YG2/YG1)*XNU12
  288.       XNU32=(YG3/YG2)*XNU23
  289.       XNU31=(YG3/YG1)*XNU13
  290.       AUX=(1.D0-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31
  291.      &    -2.D0*XNU21*XNU32*XNU13)
  292.       AUX1=AUX/YG1
  293.       AUX2=AUX/YG2
  294.       AUX3=AUX/YG3
  295.       DDAUX(1,1)=(1.D0-XNU23*XNU32)/AUX1
  296.       DDAUX(1,2)=(XNU21+XNU31*XNU23)/AUX1
  297.       DDAUX(2,1)=DDAUX(1,2)
  298.       DDAUX(1,3)=(XNU31+XNU21*XNU32)/AUX1
  299.       DDAUX(3,1)=DDAUX(1,3)
  300.       DDAUX(2,2)=(1.D0-XNU13*XNU31)/AUX2
  301.       DDAUX(2,3)=(XNU32+XNU12*XNU31)/AUX2
  302.       DDAUX(3,2)=DDAUX(2,3)
  303.       DDAUX(3,3)=(1.D0-XNU12*XNU21)/AUX3
  304.       DDAUX(4,4)=G23
  305.       DDAUX(5,5)=G13
  306.       DDAUX(6,6)=G12
  307. *
  308. * On recopie la matrice de hook pour l'inverser
  309. *
  310. 200   CONTINUE
  311.       DO 205 JLOOP=1,LHOOKL
  312.          DO 206 ILOOP=1,LHOOKL
  313.             ZDAUX(ILOOP,JLOOP)= DDAUX(ILOOP,JLOOP)
  314. 206      CONTINUE
  315. 205   CONTINUE
  316. * Inversion de la matrice de hook
  317.       TPREC= 0.D1
  318.       IPERR=0
  319.       CALL INVALM(ZDAUX,LHOOKL,LHOOKL,IPERR,TPREC)
  320.       IF (IPERR.NE.0) THEN
  321.         WRITE (IOIMP,*) 'ERREUR DANS L''INVERSION DE LA MATRICE DE HOOK'
  322.       ENDIF
  323. * Calcul des tenseurs d'endommagement H01, H02, H03
  324.       CALL ZERO (H01,LHOOKL,LHOOKL)
  325.       CALL ZERO (H02,LHOOKL,LHOOKL)
  326.       CALL ZERO (H03,LHOOKL,LHOOKL)
  327.       H01(1,1)=H1(1,1)*H1N*ZDAUX(1,1)
  328.       H01(5,5)=H1(5,5)*H1HP*ZDAUX(5,5)
  329.       H01(6,6)=H1(6,6)*H1P*ZDAUX(6,6)
  330.       H02(2,2)=H2(2,2)*H2N*ZDAUX(2,2)
  331.       H02(4,4)=H2(4,4)*H2HP*ZDAUX(4,4)
  332.       H02(6,6)=H2(6,6)*H2P*ZDAUX(6,6)
  333.       H03(3,3)=H3(3,3)*H3N*ZDAUX(3,3)
  334.       H03(4,4)=H3(4,4)*H3P*ZDAUX(4,4)
  335.       H03(5,5)=H3(5,5)*H3P*ZDAUX(5,5)
  336.  
  337. *****  VI sont rangées dans l'ordre de l'ONERA cad rep orthotrope
  338. ****** 12 ==> 6,    13 ==> 5 et 23 ==> 4
  339.  
  340. * Calcul de la nouvelle défo repère glob ordre cast
  341.       DO 451 ILOOP=1,6
  342.          EPSTFV(ILOOP)=EPST0(ILOOP)+DEPST(ILOOP)
  343.         EPSTF(ILOOP)=EPSTFV(ILOOP)
  344.  451   CONTINUE
  345. * Par contre les depstv cad après les DORTH sont rangés dans l'ordre castem
  346. *  : 12 ==> 4 13 ==> 5 23 ==>6
  347. *
  348.       CALL CICROT (wrk52,wrk53,wrk54,1,EPSTFV,DORTH)
  349. *
  350. * On les réorganise donc pour etre dans l'ordre ONERA
  351.       VARTMP= DORTH(4)
  352.       DORTH(4)= DORTH(6)
  353.       DORTH(6)= VARTMP
  354.       DEFOTH(1,1) = DORTH(1)
  355.       DEFOTH(2,2) = DORTH(2)
  356.       DEFOTH(3,3) = DORTH(3)
  357.       DEFOTH(1,2) = DORTH(6)
  358.       DEFOTH(1,3) =  DORTH(5)
  359.       DEFOTH(2,1) = DEFOTH(1,2)
  360.       DEFOTH(2,3) =  DORTH(4)
  361.       DEFOTH(3,1) = DEFOTH(1,3)
  362.       DEFOTH(3,2) = DEFOTH(2,3)
  363. *
  364.        CALL JACOBA(DEFOTH,3,PROP,XPROP)
  365. C       PRINT *,'VERIF VAL PROP',PROP(1),PROP(2),PROP(3)
  366. C Je ne garde que la partie positive des valeurs propres
  367.        IF(PROP(1).LT.0.) PROP(1)=0.D0
  368.        IF(PROP(2).LT.0.) PROP(2)=0.D0
  369.        IF(PROP(3).LT.0.) PROP(3)=0.D0
  370. C Je remets ma matrice de valeurs propres corrigée axes ortho
  371.       DEFOTH(1,1) = PROP(1)
  372.       DEFOTH(2,2) = PROP(2)
  373.       DEFOTH(3,3) = PROP(3)
  374.       DEFOTH(1,2) = 0.D0
  375.       DEFOTH(1,3) = 0.D0
  376.       DEFOTH(2,1) = 0.D0
  377.       DEFOTH(2,3) = 0.D0
  378.       DEFOTH(3,1) = 0.D0
  379.       DEFOTH(3,2) = 0.D0
  380. C XPROP(3,3) est la matrice donnant les coordonnées des vecteurs propres
  381. C par rapport aux axes d'orthotropie
  382. C 1 er indice direction 2° n° vecteur propre
  383. C on a besoin de la matrice donnant les axes ortho dans le repère
  384. C des vecteurs propres
  385. C
  386.              CALL INVER3(XPROP,XPINV)
  387.              CALL TRSPOD(XPINV,3,3,XPINVT)
  388.              CALL MULMAT(DEFOTHT,DEFOTH,XPINV,3,3,3)
  389.              CALL MULMAT(DEFOTH,XPINVT,DEFOTHT,3,3,3)
  390.  
  391. C en principe on a la matrice de defo efficasse dans les axes ortho defoth
  392. C calcul des driving forces dans le cas scalaire
  393. C mode normal
  394. C
  395.           Y1N =0.5*DDAUX(1,1)*DEFOTH(1,1)*DEFOTH(1,1)
  396.           Y2N =0.5*DDAUX(2,2)*DEFOTH(2,2)*DEFOTH(2,2)
  397.           Y3N =0.5*DDAUX(3,3)*DEFOTH(3,3)*DEFOTH(3,3)
  398. C mode tangentiel
  399. C
  400.           Y1T =0.5*(DDAUX(6,6)*DEFOTH(1,2)*DEFOTH(1,2)+b*
  401.      &         DDAUX(5,5)*DEFOTH(1,3)*DEFOTH(1,3))
  402.           Y2T =0.5*(DDAUX(6,6)*DEFOTH(1,2)*DEFOTH(1,2)+b*
  403.      &         DDAUX(4,4)*DEFOTH(2,3)*DEFOTH(2,3))
  404.           Y3T =0.5*(DDAUX(4,4)*DEFOTH(2,3)*DEFOTH(2,3)+b*
  405.      &         DDAUX(5,5)*DEFOTH(1,3)*DEFOTH(1,3))
  406.  
  407.           IF(Y1N.LT.0.) THEN
  408.            YN(1)=0.
  409.           ELSE
  410.            YN(1)=SQRT(Y1N)
  411.           ENDIF
  412.           IF(Y2N.LT.0.) THEN
  413.            YN(2)=0.
  414.           ELSE
  415.            YN(2)=SQRT(Y2N)
  416.           ENDIF
  417.           IF(Y3N.LT.0.) THEN
  418.            YN(3)=0.
  419.           ELSE
  420.            YN(3)=SQRT(Y3N)
  421.           ENDIF
  422.           IF(Y1T.LT.0.) THEN
  423.            YT(1)=0.
  424.           ELSE
  425.            YT(1)=SQRT(Y1T)
  426.           ENDIF
  427.           IF(Y2T.LT.0.) THEN
  428.            YT(2)=0.
  429.           ELSE
  430.            YT(2)=SQRT(Y2T)
  431.           ENDIF
  432.           IF(Y3T.LT.0.) THEN
  433.            YT(3)=0.
  434.           ELSE
  435.            YT(3)=SQRT(Y3T)
  436.           ENDIF
  437. C Calcul des variables d'endomagemment
  438.         Do 722 IDRF=1,3
  439.          FAFA =YN(IDRF)-Y0N(IDRF)
  440.  
  441.         IF (FAFA.LT.0.) FAFA=0.
  442.         GIN = DCN(IDRF)*(1.D0 - EXP(-1.*
  443.      &        ((FAFA/YCN(IDRF))**PN(IDRF))))
  444. C
  445.         FAFA =YT(IDRF)-Y0T(IDRF)
  446.         IF (FAFA.LT.0.) FAFA=0.
  447.         GIT = DCT(IDRF)*(1.D0 - EXP(-1.*
  448.      &        ((FAFA/YCT(IDRF))**PT(IDRF))))
  449. C
  450.         DOM(IDRF) = GIN + GIT
  451.  722    CONTINUE
  452.       DOM1= MAX(DOM(1),VAR0(2))
  453.       DOM2= MAX(DOM(2),VAR0(3))
  454.       DOM3= MAX(DOM(3),VAR0(4))
  455.  
  456.       IF (PARAM.EQ.1) THEN
  457.         WRITE (IOIMP,*) 'Increment des deformations orth. '
  458.         WRITE (IOIMP,99999) (DORTH(ILOOP), ILOOP=1,6)
  459.       ENDIF
  460. *
  461. *       On récupère les deformations totales dans les var. internes
  462. *       repère ortho ordre ONERA
  463. *
  464.       DO 500 ILOOP =1,6
  465.         VARF(ILOOP+4)=DORTH(ILOOP)
  466.  500   CONTINUE
  467. * On peut calculer EPSILON CLosure
  468. * en fonction de deltal et tzero
  469. * et de la temperature TURE0(1)ou TUREF(1)du segment WRK52 de l'include DECHE
  470. *
  471.       DO 701 IDOU=1,6
  472.        EPSCLO(IDOU) = DELTAL*(TURE0(1)-TZERO)
  473.        EPSBAR(IDOU) = DORTH(IDOU)-EPSCLO(IDOU)
  474.        EPSIM(IDOU) = (DORTH(IDOU)+VAR0(IDOU+4))/2.
  475.  701  CONTINUE
  476. *
  477. * Calcul des valeurs délimitant les  déformations pour
  478. * lesquelles les fissures commencent et finissent à/de se fermer
  479. * en vue du calcul de l'index de désactivation
  480.       DEPSIF(1)= (1.D0 + (AIF(1)*DOM1))*DEPSIF0(1)
  481.       DEPSIF(2)= (1.D0 + (AIF(2)*DOM2))*DEPSIF0(2)
  482.       DEPSIF(3)= (1.D0 + (AIF(3)*DOM3))*DEPSIF0(3)
  483. *
  484. * Calcul de l'index d'activation 1 si actif 0 si passif
  485. *
  486.       DO 703 IDO2=1,3
  487.        IF(EPSBAR(IDO2).GE.DEPSIF(IDO2)) THEN
  488.        ZNUACT(IDO2)= 1.D0
  489.        ELSEIF(EPSBAR(IDO2).GT.(-1.*(DEPSIF(IDO2)))) THEN
  490.        ZNUACT(IDO2)= 0.5D0*(1.D0-(COS(0.5D0*XPI*(EPSBAR(IDO2)
  491.      &              +DEPSIF(IDO2))/DEPSIF(IDO2))))
  492.        ELSE
  493.        ZNUACT(IDO2)= 0.D0
  494.        ENDIF
  495.  703   CONTINUE
  496. *
  497. *  Calcul des déformations résiduelles et stockées
  498. *
  499. * Calcul de la variations du dommage
  500.       DELTDO1 = DOM1 - VAR0(2)
  501.       DELTDO2 = DOM2 - VAR0(3)
  502.       DELTDO3 = DOM3 - VAR0(4)
  503. * Calcul de la variation d'indice d'activation
  504.       DELTNU1 = ZNUACT(1) - VAR0(11)
  505.       DELTNU2 = ZNUACT(2) - VAR0(12)
  506.       DELTNU3 = ZNUACT(3) - VAR0(13)
  507. *
  508. * Calcul de Seff =somme sur i de nui di hi0
  509. *
  510.        DO 744 INDJ3=1,6
  511.           Do 745 INDJ4=1,6
  512.           SEFF(INDJ3,INDJ4) = ZDAUX(INDJ3,INDJ4)
  513.      &       +ZNUACT(1)*DOM1*H01(INDJ3,INDJ4)+
  514.      &      ZNUACT(2)*DOM2*H02(INDJ3,INDJ4)+
  515.      &      ZNUACT(3)*DOM3*H03(INDJ3,INDJ4)
  516.            CEFF(INDJ3,INDJ4) =SEFF(INDJ3,INDJ4)
  517.  745     CONTINUE
  518.  744     CONTINUE
  519. C
  520.       CALL INVALM(CEFF,LHOOKL,LHOOKL,IPERR,TPREC)
  521. *   Déformation stockées (S) résiduelle (R)
  522. *
  523.       CALL MULMAT (TERME,H01,CEFF,6,6,6)
  524.       CALL MULMAT (ZK1,CEFF,TERME,6,6,6)
  525.       CALL MULMAT (TERME,H02,CEFF,6,6,6)
  526.       CALL MULMAT (ZK2,CEFF,TERME,6,6,6)
  527.       CALL MULMAT (TERME,H03,CEFF,6,6,6)
  528.       CALL MULMAT (ZK3,CEFF,TERME,6,6,6)
  529.        DO 741 INDI=1,6
  530.        SISIS(INDI)=0.D0
  531.        SISIR(INDI)=0.D0
  532.         DO 742 INDJ=1,6
  533.         SOMS(INDI,INDJ)= DELTNU1*DOM1*ZK1(INDI,INDJ)+
  534.      $              DELTNU2*DOM2*ZK2(INDI,INDJ)+
  535.      $              DELTNU3*DOM3*ZK3(INDI,INDJ)
  536.         SISIS(INDI)=SISIS(INDI)-SOMS(INDI,INDJ)*EPSIM(INDJ)
  537.         SOMR(INDI,INDJ)= ETA1*ZNUACT(1)*DELTDO1*ZK1(INDI,INDJ)+
  538.      $              ETA2*ZNUACT(2)*DELTDO2*ZK2(INDI,INDJ)+
  539.      $              ETA3*ZNUACT(3)*DELTDO3*ZK3(INDI,INDJ)
  540.         SISIR(INDI)=SISIR(INDI)+SOMR(INDI,INDJ)*EPSIM(INDJ)
  541.   742   CONTINUE
  542.  
  543.  
  544.   741   CONTINUE
  545.       CALL MULMAT (EPSR1,ZDAUX,SISIR,6,1,6)
  546.       CALL MULMAT (EPSS1,ZDAUX,SISIS,6,1,6)
  547.  
  548. *      Je retranche aux defos les defos residuelle et stockées
  549. *      j'ai dejà mis dans DORTH la defo ini
  550.        DO 746 INDI = 1,6
  551.        EPSIN(INDI)=VAR0(13+INDI)+EPSR1(INDI) +EPSS1(INDI)
  552.         EPSEL(INDI) = DORTH(INDI)- EPSIN(INDI)
  553.  746   CONTINUE
  554. *
  555. * Calcul des contraintes a la fin du pas
  556. * SORTH contient les contraintes a la fin du pas
  557. * dans le repere orthotrope
  558.        CALL MULMAT (SORTH1,CEFF,DORTH,6,1,6)
  559.        CALL MULMAT (SORTH2,DDAUX,EPSIN,6,1,6)
  560.        DO 749 INDIA=1,6
  561.        SORTH(INDIA) = SORTH1(INDIA) - SORTH2(INDIA)
  562.  749   CONTINUE
  563.  
  564.       IF (PARAM.EQ.1) THEN
  565.          WRITE (IOIMP,*) 'Contraintes calculees (repere orth.)'
  566.          WRITE (IOIMP,99999) (SORTH(ILOOP), ILOOP=1,6)
  567.       ENDIF
  568. * ATTENTION:
  569. * On reorganise les contraintes en cisaillement pour
  570. * calculer les contraintes dans le repere global
  571. *
  572.       SIGFV(1)=SORTH(1)
  573.       SIGFV(2)=SORTH(2)
  574.       SIGFV(3)=SORTH(3)
  575.       SIGFV(4)=SORTH(6)
  576.       SIGFV(5)=SORTH(5)
  577.       SIGFV(6)=SORTH(4)
  578.       EPSINV(1)=EPSIN(1)
  579.       EPSINV(2)=EPSIN(2)
  580.       EPSINV(3)=EPSIN(3)
  581.       EPSINV(4)=EPSIN(6)
  582.       EPSINV(5)=EPSIN(5)
  583.       EPSINV(6)=EPSIN(4)
  584.       VARF(1)= 0.D0
  585.       DO 186 INDI = 1,6
  586.       VARF(13+INDI)= EPSIN(INDI)
  587.       VARF(1)= VARF(1)+ EPSIN(INDI)
  588.  186  CONTINUE
  589.  
  590.       IF (PARAM.EQ.1) THEN
  591.          WRITE (IOIMP,*) 'Apres SICCNT e reorg'
  592.          WRITE (IOIMP,*) 'Contraintes reorganisees'
  593.          WRITE (IOIMP,99999) (SIGFV(ILOOP), ILOOP=1,6)
  594.       ENDIF
  595. * On appele SICROT pour trouver les contraintes dans le repère global
  596. * On utilise SORTH pour garder le resultat
  597. * (on peut pas passer SIGF à la subroutine car il est dans un segment)
  598. *
  599.       CALL CICROT (wrk52,wrk53,wrk54,0,SIGFV,SORTH)
  600. *
  601.       IF (IERR.EQ.1) THEN
  602.       PRINT *,'KERRE =1 '
  603.           KERRE=1
  604.       ENDIF
  605. * On recopie SORTH in SIGF
  606.       DO 1000 ILOOP=1,6
  607.          SIGF(ILOOP)=SORTH(ILOOP)
  608. 1000  CONTINUE
  609.       IF (PARAM.EQ.1) THEN
  610.          WRITE (IOIMP,*) 'Contraintes calculees'
  611.          WRITE (IOIMP,99999) (SIGF(ILOOP), ILOOP=1,6)
  612.       ENDIF
  613.       VARF(2)=  DOM1
  614.       VARF(3)=  DOM2
  615.       VARF(4)=  DOM3
  616.       VARF(11)=  ZNUACT(1)
  617.       VARF(12)=  ZNUACT(2)
  618.       VARF(13)=  ZNUACT(3)
  619. 99999    format(2x,'  ',(6(1x,3pe12.5),/))
  620. 99998    format(2x,'  ',(4(1x,3pe12.5),/))
  621. 99997    format(2x,'  ',(1x,3pe12.5),/)
  622. 99996    format(2x,'  ',(3(1x,1pe12.3),/))
  623.       RETURN
  624.       END
  625.  
  626.  
  627.  

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