Sources Esope | Cast3M

dohmao
C DOHMAO    SOURCE    BP208322  17/03/01    21:17:02     9325           
 
C=======================================================================
C=                            D O H M A O                              =
C=                            -----------                              =
C=                                                                     =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=  Calcul de la matrice de Hooke (elasticite lineaire) dans le cas    =
C=  des materiaux orthotropes et anisotropes                           =
C=                                                                     =
C=  Parametres :  (E)=Entree  (S)=Sortie                               =
C=  ------------                                                       =
C=   VELA     (E)   Valeurs des parametres elastiques du materiau      =
C=   MATE     (E)   Chaine contenant le type du materiau               =
C=   IFOU     (E)   Mode de calcul (cf. IFOUR de CCOPTIO)              =
C=   IDIM     (E)   Dimension du probleme (1,2,3) (cf. CCOPTIO)        =
C=   TXR      (E)   Cosinus directeurs des axes LOCAUX / repere GLOBAL =
C=   XLOC     (S)   Tableau de travail contenant le repere             =
C=                  d'orthotropie  / repere LOCAL                      =                                  
C=   XGLOB    (S)   Tableau de travail contenant le repere             =
C=                  d'orthotropie  / repere GLOBAL                     =                                  
C=   D1HOOK   (S)   Matrice de Hooke dans le repere LOCAL d'orthotropie=
C=   DDHOOK   (S)   Matrice de Hooke dans le repere GLOBAL             =
C=   LHOOK    (S)   Dimension de la matrice de Hooke                   =
C=   KCAS     (E)   =1 si on veut la matrice de Hooke associee au mode =
C=                  de calcul (non necessairement inversible)          =
C=                  =2 si on veut la matrice de Hooke complete pour    =
C=                  l'inverser par la suite (matrice de souplesse)     =
C=   IRET     (S)   =1 si option existante , =0 sinon                  =
C=                                                                     =
C=  Remarque : En DIMEnsion 1, on calcule directement la matrice DDHOOK=
C=  ---------- car les directions d'ORTHOTROPIE correspondent a celles =
C=             du repere GLOBAL.                                       =
C=======================================================================
 
      SUBROUTINE DOHMAO (VELA,MATE,IFOU,IDIM,TXR,XLOC,XGLOB,
     .                   D1HOOK,ROTHOO,DDHOOK,LHOOK,KCAS,IRET)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
      PARAMETER (XZero=0.,XUn=1.,XDeux=2.)
 
      DIMENSION VELA(*),DDHOOK(LHOOK,*),D1HOOK(LHOOK,*),ROTHOO(LHOOK,*),
     .          TXR(IDIM,*),XLOC(3,3),XGLOB(3,3)
      CHARACTER*(*) MATE
 
      IRET=1
      CALL ZERO(DDHOOK,LHOOK,LHOOK)
      CALL ZERO(XLOC,3,3)
 
C  1 - MATERIAU ORTHOTROPE
C =========================
      IF (MATE.EQ.'ORTHOTRO') THEN
C =====
C  1.1 - Definition de la matrice de HOOKE dans les axes d'ORTHOTROPIE
C =====
C= -> Contraintes planes et KCAS=1
        IF (IFOU.EQ.-2.AND.KCAS.EQ.1) THEN
          CALL ZERO(D1HOOK,LHOOK,LHOOK)
          YG1=VELA(1)
          YG2=VELA(2)
          XNU12=VELA(3)
          XNU21=(YG2/YG1)*XNU12
          AUX=XUn/(XUn-XNU12*XNU21)
          D1HOOK(1,1)=YG1*AUX
          D1HOOK(2,1)=XNU21*YG1*AUX
          D1HOOK(1,2)=D1HOOK(2,1)
          D1HOOK(2,2)=YG2*AUX
          D1HOOK(4,4)=VELA(4)
C= -> Contraintes planes et KCAS=2
        ELSEIF (IFOU.EQ.-2.AND.KCAS.EQ.2) THEN
          CALL ZERO(D1HOOK,LHOOK,LHOOK)
          YG1=VELA(1)
          YG2=VELA(2)
          YG3=VELA(7)
          XNU12=VELA(3)
          XNU23=VELA(8)
          XNU13=VELA(9)
          XNU21=(YG2/YG1)*XNU12
          XNU32=(YG3/YG2)*XNU23
          XNU31=(YG3/YG1)*XNU13
          AUX=XUn/( XUn-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31
     .             -XDeux*XNU12*XNU23*XNU31)
          AUX1=AUX*YG1
          AUX2=AUX*YG2
          D1HOOK(1,1)=(XUn-XNU23*XNU32)*AUX1
          D1HOOK(2,1)=(XNU21+XNU23*XNU31)*AUX1
          D1HOOK(3,1)=(XNU31+XNU32*XNU21)*AUX1
          D1HOOK(1,2)=D1HOOK(2,1)
          D1HOOK(2,2)=(XUn-XNU13*XNU31)*AUX2
          D1HOOK(3,2)=(XNU32+XNU31*XNU12)*AUX2
          D1HOOK(1,3)=D1HOOK(3,1)
          D1HOOK(2,3)=D1HOOK(3,2)
          D1HOOK(3,3)=(XUn-XNU12*XNU21)*AUX*YG3
          D1HOOK(4,4)=VELA(4)
C= -> Deformations planes, Axisymetrie, Serie de Fourier et 3D
        ELSEIF (IFOU.EQ.-1.OR.IFOU.EQ.-3.OR.
     .           IFOU.EQ.0.OR.IFOU.EQ.1.OR.IFOU.EQ.2) THEN
          CALL ZERO(D1HOOK,LHOOK,LHOOK)
          YG1=VELA(1)
          YG2=VELA(2)
          YG3=VELA(3)
          XNU12=VELA(4)
          XNU23=VELA(5)
          XNU13=VELA(6)
          XNU21=(YG2/YG1)*XNU12
          XNU32=(YG3/YG2)*XNU23
          XNU31=(YG3/YG1)*XNU13
          AUX=XUn/( XUn-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31
     .             -XDeux*XNU12*XNU23*XNU31)
          AUX1=AUX*YG1
          AUX2=AUX*YG2
          D1HOOK(1,1)=(XUn-XNU23*XNU32)*AUX1
          D1HOOK(2,1)=(XNU21+XNU23*XNU31)*AUX1
          D1HOOK(3,1)=(XNU31+XNU32*XNU21)*AUX1
          D1HOOK(1,2)=D1HOOK(2,1)
          D1HOOK(2,2)=(XUn-XNU13*XNU31)*AUX2
          D1HOOK(3,2)=(XNU32+XNU31*XNU12)*AUX2
          D1HOOK(1,3)=D1HOOK(3,1)
          D1HOOK(2,3)=D1HOOK(3,2)
          D1HOOK(3,3)=(XUn-XNU12*XNU21)*AUX*YG3
          D1HOOK(4,4)=VELA(7)
          IF (IFOU.EQ.1.OR.IFOU.EQ.2) THEN
            D1HOOK(5,5)=VELA(9)
            D1HOOK(6,6)=VELA(8)
          ENDIF
C= -> Modes de calcul 1D
C=    On calcule directement DDHOOK car les directions d'ORTHOTROPIE
C=    correspondent aux direction du repere GLOBAL
        ELSEIF (IFOU.GE.3.AND.IFOU.LE.15) THEN
          IF (IFOU.EQ.6.AND.KCAS.EQ.1) THEN
            DDHOOK(1,1)=VELA(1)
          ELSEIF ((IFOU.EQ.5.OR.IFOU.EQ.10.OR.IFOU.EQ.13).AND.
     .             KCAS.EQ.1) THEN
            YG1=VELA(1)
            YG3=VELA(3)
            XNU13=VELA(6)
            XNU31=(YG3/YG1)*XNU13
            AUX=XUn/(XUn-XNU13*XNU31)
            DDHOOK(1,1)=YG1*AUX
            DDHOOK(3,1)=XNU13*YG3*AUX
            DDHOOK(1,3)=DDHOOK(3,1)
            DDHOOK(3,3)=YG3*AUX
          ELSEIF ((IFOU.EQ.4.OR.IFOU.EQ.8).AND.KCAS.EQ.1) THEN
            YG1=VELA(1)
            YG2=VELA(2)
            XNU12=VELA(4)
            XNU21=(YG2/YG1)*XNU12
            AUX=XUn/(XUn-XNU12*XNU21)
            DDHOOK(1,1)=YG1*AUX
            DDHOOK(2,1)=XNU21*YG1*AUX
            DDHOOK(1,2)=DDHOOK(2,1)
            DDHOOK(2,2)=YG2*AUX
          ELSE
            YG1=VELA(1)
            YG2=VELA(2)
            YG3=VELA(3)
            XNU12=VELA(4)
            XNU23=VELA(5)
            XNU13=VELA(6)
            XNU21=(YG2/YG1)*XNU12
            XNU32=(YG3/YG2)*XNU23
            XNU31=(YG3/YG1)*XNU13
            AUX=XUn/(XUn-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31
     .              -XDeux*XNU12*XNU23*XNU31)
            AUX1=AUX*YG1
            AUX2=AUX*YG2
            DDHOOK(1,1)=(XUn-XNU23*XNU32)*AUX1
            DDHOOK(2,1)=(XNU21+XNU23*XNU31)*AUX1
            DDHOOK(3,1)=(XNU31+XNU32*XNU21)*AUX1
            DDHOOK(1,2)=DDHOOK(2,1)
            DDHOOK(2,2)=(XUn-XNU13*XNU31)*AUX2
            DDHOOK(3,2)=(XNU32+XNU31*XNU12)*AUX2
            DDHOOK(1,3)=DDHOOK(3,1)
            DDHOOK(2,3)=DDHOOK(3,2)
            DDHOOK(3,3)=(XUn-XNU12*XNU21)*AUX*YG3
          ENDIF
C= -> Cas non prevus actuellement
        ELSE
          IRET=0
          RETURN
        ENDIF
C =====
C  1.2 - Definition de la matrice de passage des axes d'ORTHOTROPIE
C        aux axes LOCAUX au point considere (dimensions 2 et 3)
C =====
cbp     IF (IDIM.EQ.2) THEN
        IF (IDIM.EQ.2.AND.IFOU.NE.1) THEN
          IF (IFOU.EQ.-2) THEN
            XLOC(1,1)=VELA(5)
            XLOC(2,1)=VELA(6)
          ELSEIF (IFOU.EQ.-1.OR.IFOU.EQ.-3.OR.IFOU.EQ.0) THEN
            XLOC(1,1)=VELA(8)
            XLOC(2,1)=VELA(9)
c           ELSEIF (IFOU.EQ.1) THEN
c             XLOC(1,1)=VELA(10)
c             XLOC(2,1)=VELA(11)
          ENDIF
          XLOC(1,2)=-XLOC(2,1)
          XLOC(2,2)=XLOC(1,1)
cbp        ELSEIF (IDIM.EQ.3) THEN
        ELSE
          XLOC(1,1)=VELA(10)
          XLOC(2,1)=VELA(11)
          XLOC(3,1)=VELA(12)
          XLOC(1,2)=VELA(13)
          XLOC(2,2)=VELA(14)
          XLOC(3,2)=VELA(15)
          CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR)
        ENDIF
 
C  2 - MATERIAU ANISOTROPE
C =========================
      ELSEIF (MATE.EQ.'ANISOTRO') THEN
C =====
C  2.1 - Definition de la matrice de HOOKE dans les axes d'ANISOTROPIE
C =====
        CALL ZERO(D1HOOK,LHOOK,LHOOK)
C= -> Contraintes planes et KCAS=1
        IF (IFOU.EQ.-2.AND.KCAS.EQ.1) THEN
          D11=VELA(1)
          D21=VELA(2)
          D22=VELA(3)
          D31=VELA(9)
          D32=VELA(10)
          D33=VELA(11)
          IF (D33.EQ.XZero) D33=D11*1.E-6
          D41=VELA(4)
          D42=VELA(5)
          D43=VELA(12)
          D44=VELA(6)
          D1HOOK(1,1)=D11-(D31*D31/D33)
          D1HOOK(2,1)=D21-(D31*D32/D33)
          D1HOOK(4,1)=D41-(D31*D43/D33)
          D1HOOK(1,2)=D1HOOK(2,1)
          D1HOOK(2,2)=D22-(D32*D32/D33)
          D1HOOK(4,2)=D42-(D32*D43/D33)
          D1HOOK(1,4)=D1HOOK(4,1)
          D1HOOK(2,4)=D1HOOK(4,2)
          D1HOOK(4,4)=D44-(D43*D43/D33)
C= -> Contraintes planes et KCAS=2
        ELSEIF (IFOU.EQ.-2.AND.KCAS.EQ.2) THEN
          D1HOOK(1,1)=VELA(1)
          D1HOOK(2,1)=VELA(2)
          D1HOOK(3,1)=VELA(9)
          D1HOOK(4,1)=VELA(4)
          D1HOOK(1,2)=D1HOOK(2,1)
          D1HOOK(2,2)=VELA(3)
          D1HOOK(3,2)=VELA(10)
          D1HOOK(4,2)=VELA(5)
          D1HOOK(3,3)=VELA(11)
          D1HOOK(1,3)=D1HOOK(3,1)
          D1HOOK(2,3)=D1HOOK(3,2)
          D1HOOK(4,3)=VELA(12)
          D1HOOK(1,4)=D1HOOK(4,1)
          D1HOOK(2,4)=D1HOOK(4,2)
          D1HOOK(3,4)=D1HOOK(4,3)
          D1HOOK(4,4)=VELA(6)
C= -> Deformations planes et Axisymetrie
        ELSEIF (IFOU.EQ.-1.OR.IFOU.EQ.-3.OR.IFOU.EQ.0) THEN
          D1HOOK(1,1)=VELA(1)
          D1HOOK(2,1)=VELA(2)
          D1HOOK(3,1)=VELA(4)
          D1HOOK(4,1)=VELA(7)
          D1HOOK(1,2)=D1HOOK(2,1)
          D1HOOK(2,2)=VELA(3)
          D1HOOK(3,2)=VELA(5)
          D1HOOK(4,2)=VELA(8)
          D1HOOK(1,3)=D1HOOK(3,1)
          D1HOOK(2,3)=D1HOOK(3,2)
          D1HOOK(3,3)=VELA(6)
          D1HOOK(4,3)=VELA(9)
          D1HOOK(1,4)=D1HOOK(4,1)
          D1HOOK(2,4)=D1HOOK(4,2)
          D1HOOK(3,4)=D1HOOK(4,3)
          D1HOOK(4,4)=VELA(10)
C= -> Serie de Fourier et 3D
        ELSEIF (IFOU.EQ.1.OR.IFOU.EQ.2) THEN
          D1HOOK(1,1)=VELA(1)
          D1HOOK(2,1)=VELA(2)
          D1HOOK(3,1)=VELA(4)
          D1HOOK(4,1)=VELA(7)
          D1HOOK(1,2)=D1HOOK(2,1)
          D1HOOK(2,2)=VELA(3)
          D1HOOK(3,2)=VELA(5)
          D1HOOK(4,2)=VELA(8)
          D1HOOK(1,3)=D1HOOK(3,1)
          D1HOOK(2,3)=D1HOOK(3,2)
          D1HOOK(3,3)=VELA(6)
          D1HOOK(4,3)=VELA(9)
          D1HOOK(1,4)=D1HOOK(4,1)
          D1HOOK(2,4)=D1HOOK(4,2)
          D1HOOK(3,4)=D1HOOK(4,3)
          D1HOOK(4,4)=VELA(10)
          IF (IFOU.EQ.1) THEN
            D1HOOK(5,5)=VELA(11)
            D1HOOK(6,5)=VELA(12)
            D1HOOK(5,6)=D1HOOK(6,5)
            D1HOOK(6,6)=VELA(13)
          ELSE
            D1HOOK(5,1)=VELA(11)
            D1HOOK(6,1)=VELA(16)
            D1HOOK(5,2)=VELA(12)
            D1HOOK(6,2)=VELA(17)
            D1HOOK(5,3)=VELA(13)
            D1HOOK(6,3)=VELA(18)
            D1HOOK(5,4)=VELA(14)
            D1HOOK(6,4)=VELA(19)
            D1HOOK(1,5)=D1HOOK(5,1)
            D1HOOK(2,5)=D1HOOK(5,2)
            D1HOOK(3,5)=D1HOOK(5,3)
            D1HOOK(4,5)=D1HOOK(5,4)
            D1HOOK(5,5)=VELA(15)
            D1HOOK(6,5)=VELA(20)
            D1HOOK(1,6)=D1HOOK(6,1)
            D1HOOK(2,6)=D1HOOK(6,2)
            D1HOOK(3,6)=D1HOOK(6,3)
            D1HOOK(4,6)=D1HOOK(6,4)
            D1HOOK(5,6)=D1HOOK(6,5)
            D1HOOK(6,6)=VELA(21)
          ENDIF
C= -> Cas non prevus actuellement
        ELSE
          IRET=0
          RETURN
        ENDIF
C =====
C  2.2 - Definition de la matrice de passage des axes d'ANSIOTROPIE
C        aux axes LOCAUX au point considere (dimensions 2 et 3)
C =====
        IF (IDIM.EQ.2) THEN
          IF (IFOU.EQ.-2) THEN
            XLOC(1,1)=VELA(7)
            XLOC(2,1)=VELA(8)
          ELSEIF (IFOU.EQ.-3.OR.IFOU.EQ.-1.OR.IFOU.EQ.0) THEN
            XLOC(1,1)=VELA(11)
            XLOC(2,1)=VELA(12)
          ELSEIF (IFOU.EQ.1) THEN
            XLOC(1,1)=VELA(14)
            XLOC(2,1)=VELA(15)
          ENDIF
          XLOC(1,2)=-XLOC(2,1)
          XLOC(2,2)=XLOC(1,1)
        ELSEIF (IDIM.EQ.3) THEN
          XLOC(1,1)=VELA(22)
          XLOC(2,1)=VELA(23)
          XLOC(3,1)=VELA(24)
          XLOC(1,2)=VELA(25)
          XLOC(2,2)=VELA(26)
          XLOC(3,2)=VELA(27)
          CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR)
        ENDIF
 
C  3 - MATERIAU UNIDIRECTIONNEL
C ==============================
      ELSEIF (MATE.EQ.'UNIDIREC') THEN
C =====
C  3.1 - Definition de la matrice de HOOKE dans les axes d'ANISOTROPIE
C =====
        CALL ZERO(D1HOOK,LHOOK,LHOOK)
        D1HOOK(1,1)=VELA(1)
C =====
C  3.2 - Definition de la matrice de passage des axes d'ANISOTROPIE
C        aux axes LOCAUX au point considere (dimensions 2 et 3)
C =====
        IF (IDIM.EQ.2) THEN
          XLOC(1,1)=VELA(2)
          XLOC(2,1)=VELA(3)
          XLOC(1,2)=-XLOC(2,1)
          XLOC(2,2)=XLOC(1,1)
        ELSEIF (IDIM.EQ.3) THEN
          XLOC(1,1)=VELA(2)
          XLOC(2,1)=VELA(3)
          XLOC(3,1)=VELA(4)
          XLOC(1,2)=VELA(5)
          XLOC(2,2)=VELA(6)
          XLOC(3,2)=VELA(7)
          CALL CROSS2(XLOC(1,1),XLOC(1,2),XLOC(1,3),IRR)
        ENDIF
 
C  4 - MATERIAUX NON PREVUS
C ==========================
      ELSE
        IRET=0
        RETURN
      ENDIF
 
C  5 - DEFINITION DE LA MATRICE DE PASSAGE DES AXES
C      D'ORTHO/ANISOTROPIE AUX AXES DU REPERE GLOBAL
C ===================================================
      IF (IDIM.EQ.1) RETURN
 
      CALL ZERO(XGLOB,3,3)
      CALL ZERO(ROTHOO,LHOOK,LHOOK)
      IDIM2=IDIM
      IF(IFOU.EQ.1) IDIM2=3
      DO j=1,IDIM
        DO k=1,IDIM2
          cc=XZero
          DO i=1,IDIM
            cc=cc+TXR(j,i)*XLOC(i,k)
          ENDDO
          XGLOB(k,j)=cc
        ENDDO
      ENDDO
cbp   en 2D Fourier, vrai TXR = [TXR(2x2) [0] ; [0] 1]   
      IF (IFOU.EQ.1) THEN
        XGLOB(1,3)=XLOC(3,1)
        XGLOB(2,3)=XLOC(3,2)
        XGLOB(3,3)=XLOC(3,3)        
      ENDIF
      IF (IDIM.EQ.2.AND.IFOU.NE.1) THEN
        ROTHOO(1,1)=XGLOB(1,1)*XGLOB(1,1)
        ROTHOO(2,1)=XGLOB(2,1)*XGLOB(2,1)
        ROTHOO(4,1)=XDeux*XGLOB(1,1)*XGLOB(2,1)
        ROTHOO(1,2)=XGLOB(1,2)*XGLOB(1,2)
        ROTHOO(2,2)=XGLOB(2,2)*XGLOB(2,2)
        ROTHOO(4,2)=XDeux*XGLOB(1,2)*XGLOB(2,2)
        ROTHOO(3,3)=XUn
        ROTHOO(1,4)=XGLOB(1,1)*XGLOB(1,2)
        ROTHOO(2,4)=XGLOB(2,1)*XGLOB(2,2)
        ROTHOO(4,4)=XGLOB(1,2)*XGLOB(2,1)+XGLOB(1,1)*XGLOB(2,2)
c         IF (IFOU.EQ.1) THEN
c           ROTHOO(5,5)=XGLOB(1,1)
c           ROTHOO(6,5)=XGLOB(2,1)
c           ROTHOO(5,6)=XGLOB(1,2)
c           ROTHOO(6,6)=XGLOB(2,2)
c         ENDIF
cbp      ELSEIF (IDIM.EQ.3) THEN
      ELSE
        IF (LHOOK.LT.6) THEN
          IRET=0
          RETURN
        ENDIF
        DO i=1,3
          ROTHOO(1,i)=XGLOB(1,i)*XGLOB(1,i)
          ROTHOO(2,i)=XGLOB(2,i)*XGLOB(2,i)
          ROTHOO(3,i)=XGLOB(3,i)*XGLOB(3,i)
          ROTHOO(i,4)=XGLOB(i,1)*XGLOB(i,2)
          ROTHOO(i,5)=XGLOB(i,2)*XGLOB(i,3)
          ROTHOO(i,6)=XGLOB(i,1)*XGLOB(i,3)
          ROTHOO(4,i)=XDeux*XGLOB(1,i)*XGLOB(2,i)
          ROTHOO(5,i)=XDeux*XGLOB(2,i)*XGLOB(3,i)
          ROTHOO(6,i)=XDeux*XGLOB(1,i)*XGLOB(3,i)
        ENDDO
        DO j=4,6
          j1=j-3
          j2=j1+1
          IF (j2.GT.3) j2=j2-3
          DO i=4,6
            i1=i-3
            i2=i1+1
            IF (i2.GT.3) i2=i2-3
            ROTHOO(i,j)= XGLOB(i1,j1)*XGLOB(i2,j2)
     .                  +XGLOB(i1,j2)*XGLOB(i2,j1)
          ENDDO
        ENDDO
        DO i=1,6
          AA=ROTHOO(6,i)
          ROTHOO(6,i)=ROTHOO(5,i)
          ROTHOO(5,i)=AA
        ENDDO
        DO i=1,6
          AA=ROTHOO(i,6)
          ROTHOO(i,6)=ROTHOO(i,5)
          ROTHOO(i,5)=AA
        ENDDO
      ENDIF
 
C  6 - CALCUL DE LA MATRICE DE HOOKE DANS LE REPERE GLOBAL
C =========================================================
      CALL PRODT(DDHOOK,D1HOOK,ROTHOO,LHOOK,LHOOK)
 
      RETURN
      END