1 : C TRIHM3 SOURCE CHAT 05/01/13 03:47:17 5004 2 : SUBROUTINE TRIHM3(IGAU,ITEL,NBNO,XEL,SHPTOT,SHP,RHOS,RHOF,B11,B22, 3 : #B12,F11,F12,SFLU,SCEL,POIGAU,VKL12,VKL22,VKL23,VKL33,LRE, 4 : #REL,ISDJC,IRET) 5 : C======================================================================= 6 : C 7 : C CALCULE LA MATRICE DE MASSE DANS LE CAS PLAN POUR 8 : C LA FORMULATION (37) HOMOGENE 9 : C======================================================================= 10 : C NBDL = NOMBRE DE DDL PAR NOEUD 11 : C INPUT 12 : C IGAU=NUMERO DU POINT DE GAUSS 13 : C ITEL=NUMERO DE L ELEMENT DANS NOMTP 14 : C NBNO=NOMBRE DE NOEUDS 15 : C XEL =COORDONNEES DE L ELEMENT 16 : C RHOS = MASSE DU SOLIDE 17 : C RHOF = MASSE VOLUMIQUE DU FLUIDE 18 : C B11,B22,B12 = PERMEABILITE ACOUSTIQUE DU MILIEU 19 : C F11,F12 = NORME 20 : C SFLU = SURFACE FLUIDE DANS LA CELLULE ELEMENTAIRE 21 : C SCEL = SURFACE DE LA CELLULE ELEMENTAIRE 22 : C POIGAU=MINTE1.POIGAU(IGAU) 23 : C VKL12=-((COEFPR*COEFPI)/(RHOF*C**2))*SFLU/SCEL 24 : C VKL22=-(COEFPI**2)/(RHOF*SCEL) 25 : C VKL23= COEFPI/SCEL 26 : C VKL33= 1/SCEL 27 : C LRE =NOMBRE DE D.D.L DE LA MATRICE DE RIGIDITE 28 : C SHPTOT(6,NBNO,NBGAU)=FONCTIONS DE FORMES ET DERIVEES 29 : C ISDJC= INDICATEUR DU SIGNE DU JACOBIEN 30 : C ZONE DE TRAVAIL 31 : C SHP(6,NBNO)=TABLEAU DE TRAVAIL 32 : C OUTPUT 33 : C REL=MATRICE DE RIGIDITE 34 : C ISDJC= INDICATEUR DU SIGNE DU JACOBIEN 35 : C IRET=INDICATEUR = 1 : SUCCES 36 : C = 0 : ECHEC (ITEL INCOMPATIBLE AVEC LA FORMULATION 37 : C = 2 : ECHEC ( JACOBIEN NUL ) 38 : C======================================================================= 39 : IMPLICIT INTEGER(I-N) 40 : IMPLICIT REAL*8 (A-H,O-Z) 41 : DIMENSION XEL(3,*),SHP(6,*),SHPTOT(6,NBNO,*),REL(LRE,*) 42 : IF (ITEL.EQ.92) GOTO 10 43 : C 44 : C ERREUR : TYPE D' ELEMENT INCOMPATIBLE AVEC LA FORMULATION 45 : C 46 : IRET = 0 47 : GOTO 666 48 : 10 CONTINUE 49 : C 50 : C SHP(1,I) : FONCTION DE FORME 51 : C SHP(2,I) : DERIVEE % X DE LA FONCTION DE FORME 52 : C SHP(3,I) : DERIVEE % Y DE LA FONCTION DE FORME 53 : C 54 : DO 101 NP=1,NBNO 55 : SHP(1,NP)=SHPTOT(1,NP,IGAU) 56 : SHP(2,NP)=SHPTOT(2,NP,IGAU) 57 : SHP(3,NP)=SHPTOT(3,NP,IGAU) 58 : 101 CONTINUE 59 : IDIM=2 59 : CALL JACOBI(XEL,SHP,IDIM,NBNO,DJAC) 60 : IF (DJAC.EQ.0.) GOTO 667 61 : IRET = 1 62 : IF (DJAC.LT.0.) ISDJC = ISDJC + 1 63 : DJAC=ABS(DJAC)*POIGAU 64 : C 65 : C TERMES EN P*PI 66 : C NBDL : NOMBRE DE D.D.L PAR NOEUD ( = 4 ) 67 : C 68 : NBDL = 4 69 : IX1=0 70 : IY1=0 71 : DO 102 IX=2,LRE ,NBDL 72 : IX1=IX1 + 1 73 : DO 103 IY=1,IX ,NBDL 74 : IY1=IY1 + 1 75 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL12*SHP(1,IX1)*SHP(1,IY1) 76 : REL(IX,IY) = REL(IY,IX) 77 : 103 CONTINUE 78 : IY1=0 79 : 102 CONTINUE 80 : DO 104 IX=2,LRE-NBDL,NBDL 81 : IX1=IX +NBDL-1 82 : DO 105 IY=IX1,LRE ,NBDL 83 : REL(IX,IY) = REL(IX-1,IY+1) 84 : REL(IY,IX) = REL(IX,IY) 85 : 105 CONTINUE 86 : 104 CONTINUE 87 : C 88 : C TERMES EN PI*PI 89 : C 90 : IX1=0 91 : IY1=0 92 : DO 106 IX=2,LRE ,NBDL 93 : IX1=IX1 + 1 94 : DO 107 IY=2,IX ,NBDL 95 : IY1=IY1 + 1 96 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL22*((B11*SHP(2,IX1)*SHP(2,IY1)+ 97 : #B22*SHP(3,IX1)*SHP(3,IY1)+B12*(SHP(2,IX1)*SHP(3,IY1)+SHP(3,IX1)* 98 : #SHP(2,IY1)))) 99 : REL(IX,IY) = REL(IY,IX) 100 : 107 CONTINUE 101 : IY1=0 102 : 106 CONTINUE 103 : C 104 : C TERMES EN PI*(UX,UY) 105 : C 106 : D11=(SCEL-B11)*F12 107 : D22=(SCEL-B22)*F12 108 : D12=-B12*F12 109 : IX1=0 110 : IY1=0 111 : DO 108 IX=3,LRE ,NBDL 112 : IX1=IX1 + 1 113 : DO 109 IY=2,IX ,NBDL 114 : IY1=IY1 + 1 115 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * ( D11*SHP(2,IY1) + 116 : #D12*SHP(3,IY1) ) * SHP(1,IX1) 117 : REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL23 * ( D12*SHP(2,IY1) + 118 : #D22*SHP(3,IY1) ) * SHP(1,IX1) 119 : REL(IX,IY) = REL(IY,IX) 120 : REL(IX+1,IY) = REL(IY,IX+1) 121 : 109 CONTINUE 122 : IY1=0 123 : 108 CONTINUE 124 : IX1=1 125 : IY1=0 126 : DO 110 IX=2+NBDL,LRE ,NBDL 127 : IX1=IX1 + 1 128 : DO 111 IY=3,IX ,NBDL 129 : IY1=IY1 + 1 130 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * ( D11*SHP(2,IX1) + 131 : #D12*SHP(3,IX1) ) * SHP(1,IY1) 132 : REL(IY+1,IX) = REL(IY+1,IX) + DJAC*VKL23 * (D12*SHP(2,IX1)+ 133 : #D22*SHP(3,IX1) ) * SHP(1,IY1) 134 : REL(IX,IY) = REL(IY,IX) 135 : REL(IX,IY+1) = REL(IY+1,IX) 136 : 111 CONTINUE 137 : IY1=0 138 : 110 CONTINUE 139 : C 140 : C TERMES EN (UX,UY)*(UX,UY) 141 : C 142 : H11=RHOS + RHOF*(SFLU-B11)*F11 143 : H22=RHOS + RHOF*(SFLU-B22)*F11 144 : H12= (-1.)*RHOF*B12*F11 145 : IX1=0 146 : IY1=0 147 : DO 112 IX=3,LRE ,NBDL 148 : IX1=IX1 + 1 149 : DO 113 IY=3,IX ,NBDL 150 : IY1=IY1 + 1 151 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL33*H11*SHP(1,IY1)*SHP(1,IX1) 152 : REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL33*H12*SHP(1,IY1)*SHP(1,IX1) 153 : REL(IY+1,IX+1)=REL(IY+1,IX+1)+DJAC*VKL33*H22*SHP(1,IY1)*SHP(1,IX1) 154 : REL(IX,IY) = REL(IY,IX) 155 : REL(IY+1,IX) = REL(IY,IX+1) 156 : REL(IX+1,IY) = REL(IY,IX+1) 157 : REL(IX,IY+1) = REL(IY+1,IX) 158 : REL(IX+1,IY+1) = REL(IY+1,IX+1) 159 : 113 CONTINUE 160 : IY1=0 161 : 112 CONTINUE 162 : GOTO 666 163 : C 164 : C MESSAGE D ERREUR : ELEMENT A SURFACE NULLE 165 : C 166 : 667 CONTINUE 167 : IRET = 2 168 : GOTO 666 169 : C 170 : 666 CONTINUE 171 : RETURN 172 : END 173 :
© Cast3M 2003 - All rights reserved.
Disclaimer