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