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

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