C PORMAO SOURCE CB215821 17/07/21 21:15:26 9513 SUBROUTINE PORMAO(VELA,MATE,IFOU,IDIM,TXR,XLOC, . XGLOB,D1HOOK,ROTHOO,DDHOOK,LHOOK,COBMA,XMOB,KCAS,IRET) C----------------------------------------------------------------------- C C MATRICE DE HOOK DES ELEMENTS MASSIFS ORTHOTROPES OU ANISOTROPES C SPECIAL POUR MILIEU POREUX C C ENTREES C VELA() = materiau dans un tableau de travail C MATE = Nom du materiau C IFOU = num{ro d'harmonique de fourier: IFOUR de CCOPTIO C IDIM = DEFINIT SI ON EST EN 2D OU 3D C TXR = COS-DIRECTEURS DES AXES LOCAUX /REPERE GLOBAL C LHOOK = TAILLE DE LA MATRICE DE HOOKE C KCAS = 1 SI ON VEUT LA MATRICE POUR ELLE MEME C = 2 SI ON VEUT LA MATRICE POUR L'INVERSER ENSUITE C XLOC | C XGLOB |= TABLEAU DE TRAVAIL C D1HOOK| C C SORTIES C DDHOOK(LHOOK,LHOOK) = matrice de hooke C IRET = 1 si option existante 0 SINON C C INSPIRE DE DOHMAO C----------------------------------------------------------------------- C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(XZER=0.0D0) PARAMETER(UN=1.D0,DEUX=2.D0) CHARACTER*8 MATE C DIMENSION VELA(*),DDHOOK(LHOOK,*),TXR(IDIM,*) DIMENSION XLOC(3,3),D1HOOK(LHOOK,*) DIMENSION XGLOB(3,3),ROTHOO(LHOOK,*) DIMENSION COBMA(*),COBAUX(10) C IRET=1 CALL ZERO(D1HOOK,LHOOK,LHOOK) CALL ZERO(XGLOB,3,3) CALL ZERO(XLOC,3,3) CALL ZERO(ROTHOO,LHOOK,LHOOK) CALL ZERO(DDHOOK,LHOOK,LHOOK) CALL ZERO(COBMA,LHOOK,1) CALL ZERO(COBAUX,LHOOK,1) C 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 YG1=VELA(1) YG2=VELA(2) XNU12=VELA(3) G12=VELA(4) COBAUX(1)=VELA(7) COBAUX(2)=VELA(8) XMOB =VELA(9) XNU21=(YG2/YG1)*XNU12 AUX=UN-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)=G12 C= -> Contraintes planes et KCAS=2 ELSEIF(IFOU.EQ.-2.AND.KCAS.EQ.2)THEN YG1=VELA(1) YG2=VELA(2) XNU12=VELA(3) YG3=VELA(10) XNU23=VELA(11) XNU13=VELA(12) G12=VELA(4) COBAUX(1)=VELA(7) COBAUX(2)=VELA(8) XMOB =VELA(9) XNU21=(YG2/YG1)*XNU12 XNU32=(YG3/YG2)*XNU23 XNU31=(YG3/YG1)*XNU13 AUX=(UN-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31 . -DEUX*XNU21*XNU32*XNU13) AUX1=AUX/YG1 AUX2=AUX/YG2 AUX3=AUX/YG3 D1HOOK(1,1)=(UN-XNU23*XNU32)/AUX1 D1HOOK(1,2)=(XNU21+XNU31*XNU23)/AUX1 D1HOOK(2,1)=D1HOOK(1,2) D1HOOK(1,3)=(XNU31+XNU21*XNU32)/AUX1 D1HOOK(3,1)=D1HOOK(1,3) D1HOOK(2,2)=(UN-XNU13*XNU31)/AUX2 D1HOOK(2,3)=(XNU32+XNU12*XNU31)/AUX2 D1HOOK(3,2)=D1HOOK(2,3) D1HOOK(3,3)=(UN-XNU12*XNU21)/AUX3 D1HOOK(4,4)=G12 C= -> Deformations planes, Axisymetrie ELSEIF(IFOU.EQ.-1.OR.IFOU.EQ.0) THEN YG1=VELA(1) YG2=VELA(2) YG3=VELA(3) XNU12=VELA(4) XNU23=VELA(5) XNU13=VELA(6) G12=VELA(7) COBAUX(1)=VELA(10) COBAUX(2)=VELA(11) COBAUX(3)=VELA(12) XMOB =VELA(13) XNU21=(YG2/YG1)*XNU12 XNU32=(YG3/YG2)*XNU23 XNU31=(YG3/YG1)*XNU13 AUX=(UN-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31 . -DEUX*XNU21*XNU32*XNU13) AUX1=AUX/YG1 AUX2=AUX/YG2 AUX3=AUX/YG3 D1HOOK(1,1)=(UN-XNU23*XNU32)/AUX1 D1HOOK(1,2)=(XNU21+XNU31*XNU23)/AUX1 D1HOOK(2,1)=D1HOOK(1,2) D1HOOK(1,3)=(XNU31+XNU21*XNU32)/AUX1 D1HOOK(3,1)=D1HOOK(1,3) D1HOOK(2,2)=(UN-XNU13*XNU31)/AUX2 D1HOOK(2,3)=(XNU32+XNU12*XNU31)/AUX2 D1HOOK(3,2)=D1HOOK(2,3) D1HOOK(3,3)=(UN-XNU12*XNU21)/AUX3 D1HOOK(4,4)=G12 C= -> Serie de Fourier et 3D ELSEIF(IFOU.EQ.1.OR.IFOU.EQ.2)THEN YG1=VELA(1) YG2=VELA(2) YG3=VELA(3) XNU12=VELA(4) XNU23=VELA(5) XNU13=VELA(6) G12=VELA(7) G23=VELA(8) G13=VELA(9) IF(IFOU.EQ.1) THEN COBAUX(1)=VELA(12) COBAUX(2)=VELA(13) COBAUX(3)=VELA(14) XMOB =VELA(15) ELSEIF (IFOU.EQ.2) THEN COBAUX(1)=VELA(16) COBAUX(2)=VELA(17) COBAUX(3)=VELA(18) XMOB =VELA(19) ENDIF XNU21=(YG2/YG1)*XNU12 XNU32=(YG3/YG2)*XNU23 XNU31=(YG3/YG1)*XNU13 AUX=(UN-XNU12*XNU21-XNU23*XNU32-XNU13*XNU31 . -DEUX*XNU21*XNU32*XNU13) AUX1=AUX/YG1 AUX2=AUX/YG2 AUX3=AUX/YG3 D1HOOK(1,1)=(UN-XNU23*XNU32)/AUX1 D1HOOK(1,2)=(XNU21+XNU31*XNU23)/AUX1 D1HOOK(2,1)=D1HOOK(1,2) D1HOOK(1,3)=(XNU31+XNU21*XNU32)/AUX1 D1HOOK(3,1)=D1HOOK(1,3) D1HOOK(2,2)=(UN-XNU13*XNU31)/AUX2 D1HOOK(2,3)=(XNU32+XNU12*XNU31)/AUX2 D1HOOK(3,2)=D1HOOK(2,3) D1HOOK(3,3)=(UN-XNU12*XNU21)/AUX3 D1HOOK(4,4)=G12 D1HOOK(5,5)=G13 D1HOOK(6,6)=G23 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.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 ===== 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(13) D32=VELA(14) D33=VELA(15) IF (D33.EQ.0.D0) D33=D11*1.E0-6 D41=VELA(4) D42=VELA(5) D43=VELA(16) D44=VELA(6) D1HOOK(1,1)=D11 - ((D31**2)/D33) D1HOOK(2,1)=D21 - ((D31*D32)/D33) D1HOOK(1,2)=D1HOOK(2,1) D1HOOK(2,2)=D22 - ((D32**2)/D33) D1HOOK(4,1)=D41 - ((D31*D43)/D33) D1HOOK(1,4)=D1HOOK(4,1) D1HOOK(4,2)=D42 - ((D32*D43)/D33) D1HOOK(2,4)=D1HOOK(4,2) D1HOOK(4,4)=D44 - ((D43**2)/D33) COBAUX(1)=VELA( 9) COBAUX(2)=VELA(10) COBAUX(3)=0.D0 COBAUX(4)=VELA(11) XMOB =VELA(12) 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(1,2)=D1HOOK(2,1) D1HOOK(2,2)=VELA(3) D1HOOK(3,1)=VELA(13) D1HOOK(1,3)=D1HOOK(3,1) D1HOOK(3,2)=VELA(14) D1HOOK(2,3)=D1HOOK(3,2) D1HOOK(3,3)=VELA(15) D1HOOK(4,1)=VELA(4) D1HOOK(1,4)=D1HOOK(4,1) D1HOOK(4,2)=VELA(5) D1HOOK(2,4)=D1HOOK(4,2) D1HOOK(4,3)=VELA(16) D1HOOK(3,4)=D1HOOK(4,3) D1HOOK(4,4)=VELA(6) COBAUX(1)=VELA( 9) COBAUX(2)=VELA(10) COBAUX(3)=0.D0 COBAUX(4)=VELA(11) XMOB =VELA(12) C= -> Deformations planes et Axisymetrie ELSEIF(IFOU.EQ.-1.OR.IFOU.EQ.0) THEN D1HOOK(1,1)=VELA(1) D1HOOK(2,1)=VELA(2) D1HOOK(1,2)=D1HOOK(2,1) D1HOOK(2,2)=VELA(3) D1HOOK(3,1)=VELA(4) D1HOOK(1,3)=D1HOOK(3,1) D1HOOK(3,2)=VELA(5) D1HOOK(2,3)=D1HOOK(3,2) D1HOOK(3,3)=VELA(6) D1HOOK(4,1)=VELA(7) D1HOOK(1,4)=D1HOOK(4,1) D1HOOK(4,2)=VELA(8) D1HOOK(2,4)=D1HOOK(4,2) D1HOOK(4,3)=VELA(9) D1HOOK(3,4)=D1HOOK(4,3) D1HOOK(4,4)=VELA(10) COBAUX(1)=VELA(13) COBAUX(2)=VELA(14) COBAUX(3)=VELA(16) COBAUX(4)=VELA(15) XMOB =VELA(17) 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(1,2)=D1HOOK(2,1) D1HOOK(2,2)=VELA(3) D1HOOK(3,1)=VELA(4) D1HOOK(1,3)=D1HOOK(3,1) D1HOOK(3,2)=VELA(5) D1HOOK(2,3)=D1HOOK(3,2) D1HOOK(3,3)=VELA(6) D1HOOK(4,1)=VELA(7) D1HOOK(1,4)=D1HOOK(4,1) D1HOOK(4,2)=VELA(8) D1HOOK(2,4)=D1HOOK(4,2) D1HOOK(4,3)=VELA(9) 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) COBAUX(1)=VELA(16) COBAUX(2)=VELA(17) COBAUX(3)=VELA(19) COBAUX(4)=VELA(18) XMOB =VELA(20) ELSE D1HOOK(5,1)=VELA(11) D1HOOK(1,5)=D1HOOK(5,1) D1HOOK(5,2)=VELA(12) D1HOOK(2,5)=D1HOOK(5,2) D1HOOK(5,3)=VELA(13) D1HOOK(3,5)=D1HOOK(5,3) D1HOOK(5,4)=VELA(14) D1HOOK(4,5)=D1HOOK(5,4) D1HOOK(5,5)=VELA(15) D1HOOK(6,1)=VELA(16) D1HOOK(1,6)=D1HOOK(6,1) D1HOOK(6,2)=VELA(17) D1HOOK(2,6)=D1HOOK(6,2) D1HOOK(6,3)=VELA(18) D1HOOK(3,6)=D1HOOK(6,3) D1HOOK(6,4)=VELA(19) D1HOOK(4,6)=D1HOOK(6,4) D1HOOK(6,5)=VELA(20) D1HOOK(5,6)=D1HOOK(6,5) D1HOOK(6,6)=VELA(21) COBAUX(1)=VELA(28) COBAUX(2)=VELA(29) COBAUX(3)=VELA(30) COBAUX(4)=VELA(31) COBAUX(5)=VELA(32) COBAUX(6)=VELA(33) XMOB =VELA(34) ENDIF 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.-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 (par rapport aux axes d'armatures) C ===== D1HOOK(1,1)=VELA(1) COBAUX(1)=VELA(4) XMOB =VELA(5) 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 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 C IDIM2=IDIM IF(IFOU.EQ.1) IDIM2=3 DO j=1,IDIM DO k=1,IDIM2 cc=XZER 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 C MATRICE DE TRANSFORMATION cbp IF (IDIM.EQ.2) THEN IF (IDIM.EQ.2.AND.IFOU.NE.1) THEN ROTHOO(1,1)=XGLOB(1,1)*XGLOB(1,1) ROTHOO(1,2)=XGLOB(1,2)*XGLOB(1,2) ROTHOO(1,4)=XGLOB(1,1)*XGLOB(1,2) ROTHOO(2,1)=XGLOB(2,1)*XGLOB(2,1) ROTHOO(2,2)=XGLOB(2,2)*XGLOB(2,2) ROTHOO(2,4)=XGLOB(2,1)*XGLOB(2,2) ROTHOO(3,3)=UN ROTHOO(4,1)=DEUX*XGLOB(1,1)*XGLOB(2,1) ROTHOO(4,2)=DEUX*XGLOB(1,2)*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(5,6)=XGLOB(1,2) c ROTHOO(6,5)=XGLOB(2,1) c ROTHOO(6,6)=XGLOB(2,2) c ENDIF cbp ELSEIF(IDIM.EQ.3)THEN ELSE DO 100 IC=1,3 DO 100 IL=1,3 ROTHOO(IL,IC)=XGLOB(IL,IC)*XGLOB(IL,IC) 100 CONTINUE DO 110 IL=1,3 ROTHOO(IL,4)=XGLOB(IL,1)*XGLOB(IL,2) ROTHOO(IL,5)=XGLOB(IL,2)*XGLOB(IL,3) ROTHOO(IL,6)=XGLOB(IL,1)*XGLOB(IL,3) 110 CONTINUE DO 120 IC=1,3 ROTHOO(4,IC)=DEUX*XGLOB(1,IC)*XGLOB(2,IC) ROTHOO(5,IC)=DEUX*XGLOB(2,IC)*XGLOB(3,IC) ROTHOO(6,IC)=DEUX*XGLOB(1,IC)*XGLOB(3,IC) 120 CONTINUE DO 130 IL=4,6 IL1=IL-3 IL2=IL1+1 IF(IL2.GT.3)IL2=IL2-3 DO 130 IC=4,6 IC1=IC-3 IC2=IC1+1 IF(IC2.GT.3)IC2=IC2-3 ROTHOO(IL,IC)=XGLOB(IL1,IC1)*XGLOB(IL2,IC2)+ . XGLOB(IL1,IC2)*XGLOB(IL2,IC1) 130 CONTINUE DO 135 IC=1,6 AA=ROTHOO(6,IC) ROTHOO(6,IC)=ROTHOO(5,IC) ROTHOO(5,IC)=AA 135 CONTINUE DO 136 IL=1,6 AA=ROTHOO(IL,6) ROTHOO(IL,6)=ROTHOO(IL,5) ROTHOO(IL,5)=AA 136 CONTINUE ENDIF C 6 - CALCUL DE LA MATRICE DE HOOKE DANS LE REPERE GLOBAL C ET DE COBAUX C ========================================================= C CALL PRODT(DDHOOK,D1HOOK,ROTHOO,LHOOK,LHOOK) * DO 140 IL=1,LHOOK DO 140 IC=1,LHOOK COBMA(IL)=COBMA(IL)+ROTHOO(IC,IL)*COBAUX(IC) 140 CONTINUE RETURN END