1 : C CUBHM1 SOURCE CHAT 05/01/12 22:32:31 5004 2 : C CUBHM1 SOURCE LOFT 90/01/26 21:08:44 3 : SUBROUTINE CUBHM1(IGAU,ITEL,NBNO,XEL,SHPTOT,SHP,RHOS,RHOF,B11,B22, 4 : #B12,SFLU,SCEL,POIGAU,VKL12,VKL22,VKL23,VKL33,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 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 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.127) 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 SHP(4,I) : DERIVEE % Z DE LA FONCTION DE FORME 54 : C 55 : DO 101 NP=1,NBNO 56 : SHP(1,NP)=SHPTOT(1,NP,IGAU) 57 : SHP(2,NP)=SHPTOT(2,NP,IGAU) 58 : SHP(3,NP)=SHPTOT(3,NP,IGAU) 59 : SHP(4,NP)=SHPTOT(4,NP,IGAU) 60 : 101 CONTINUE 61 : IDIM=3 61 : CALL JACOBI(XEL,SHP,IDIM,NBNO,DJAC) 62 : 63 : C 64 : C 65 : A1 = XEL(3,2) - XEL(3,1) 66 : A1 = ABS(A1) 67 : A2 = XEL(3,3) - XEL(3,1) 68 : A2 = ABS(A2) 69 : A3 = XEL(3,4) - XEL(3,1) 70 : A3 = ABS(A3) 71 : A4 = XEL(3,5) - XEL(3,1) 72 : A4 = ABS(A4) 73 : A5 = XEL(3,6) - XEL(3,1) 74 : A5 = ABS(A5) 75 : A6 = XEL(3,7) - XEL(3,1) 76 : A6 = ABS(A6) 77 : A7 = XEL(3,8) - XEL(3,1) 78 : A7 = ABS(A7) 79 : C 80 : 81 : A3 = (0.25D0)*(A1+A2+A3+A4+A5+A6+A7) 82 : C 83 : IF (DJAC.EQ.0.) GOTO 667 84 : IRET = 1 85 : IF (DJAC.LT.0.) ISDJC = ISDJC + 1 86 : DJAC=ABS(DJAC)*POIGAU 87 : C WRITE(6,*)'DJAC=',DJAC 88 : C 89 : C 90 : C FONCTIONS DE FORME EN Z 91 : C 92 : ZH=SHP(1,5)+SHP(1,6)+SHP(1,7)+SHP(1,8) 93 : C 94 : Z1=1.D0-3.D0*ZH*ZH+2.D0*ZH*ZH*ZH 95 : Z2=3.D0*(ZH**2)-2.D0*ZH*ZH*ZH 96 : Z3=(A3*ZH)*(1.D0-2.D0*ZH+ZH*ZH) 97 : Z4=(A3*ZH)*(ZH*ZH-ZH) 98 : C 99 : C 100 : S1=SHP(1,1)+SHP(1,4)+SHP(1,5)+SHP(1,8) 101 : S2=SHP(1,2)+SHP(1,3)+SHP(1,6)+SHP(1,7) 102 : C 103 : T1=SHP(1,1)+SHP(1,2)+SHP(1,5)+SHP(1,6) 104 : T2=SHP(1,3)+SHP(1,4)+SHP(1,7)+SHP(1,8) 105 : C 106 : C 107 : C FONCTIONS DE FORME POUR LES DEPLACEMENTS 108 : C 109 : SHP(5,1)=S1*T1*Z1 110 : SHP(5,2)=S2*T1*Z1 111 : SHP(5,3)=S2*T2*Z1 112 : SHP(5,4)=S1*T2*Z1 113 : SHP(5,5)=S1*T1*Z2 114 : SHP(5,6)=S2*T1*Z2 115 : SHP(5,7)=S2*T2*Z2 116 : SHP(5,8)=S1*T2*Z2 117 : 118 : SHP(5,9)=S1*T1*Z3 119 : SHP(5,10)=S2*T1*Z3 120 : SHP(5,11)=S2*T2*Z3 121 : SHP(5,12)=S1*T2*Z3 122 : SHP(5,13)=S1*T1*Z4 123 : SHP(5,14)=S2*T1*Z4 124 : SHP(5,15)=S2*T2*Z4 125 : SHP(5,16)=S1*T2*Z4 126 : C 127 : C TERMES EN P*PI 128 : C NBDL : NOMBRE DE D.D.L PAR NOEUD ( = 6 ) 129 : C 130 : NBDL = 6 131 : IX1=0 132 : IY1=0 133 : DO 102 IX=2,LRE ,NBDL 134 : IX1=IX1 + 1 135 : DO 103 IY=1,IX ,NBDL 136 : IY1=IY1 + 1 137 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL12*SHP(1,IX1)*SHP(1,IY1) 138 : REL(IX,IY) = REL(IY,IX) 139 : 103 CONTINUE 140 : IY1=0 141 : 102 CONTINUE 142 : DO 104 IX=2,LRE-NBDL,NBDL 143 : IX1=IX +NBDL-1 144 : DO 105 IY=IX1,LRE ,NBDL 145 : REL(IX,IY) = REL(IX-1,IY+1) 146 : REL(IY,IX) = REL(IX,IY) 147 : 105 CONTINUE 148 : 104 CONTINUE 149 : C 150 : C TERMES EN PI*PI 151 : C 152 : B33=SFLU 153 : IX1=0 154 : IY1=0 155 : DO 106 IX=2,LRE ,NBDL 156 : C WRITE(6,*)'IX= ',IX 157 : IX1=IX1 + 1 158 : C WRITE(6,*)'IX1= ',IX1 159 : DO 107 IY=2,IX ,NBDL 160 : C WRITE(6,*)'IY= ',IY 161 : IY1=IY1 + 1 162 : C WRITE(6,*)'IY1= ',IY1 163 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL22*(B11*SHP(2,IX1)*SHP(2,IY1) 164 : #+B22*SHP(3,IX1)*SHP(3,IY1)+B33*SHP(4,IY1)*SHP(4,IX1)) 165 : REL(IX,IY) = REL(IY,IX) 166 : 107 CONTINUE 167 : IY1=0 168 : 106 CONTINUE 169 : C 170 : C TERMES EN PI*(UX,RY) 171 : C 172 : D11 = (SCEL - B11) 173 : D22 = (SCEL - B22) 174 : D33=SFLU 175 : C 176 : IX1=0 177 : IY1=0 178 : DO 108 IX=3,LRE ,NBDL 179 : IX1=IX1 + 1 180 : DO 109 IY=2,IX ,NBDL 181 : IY1=IY1 + 1 182 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * (D11*SHP(2,IY1)*SHP(5,IX1)) 183 : REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL23 * (D11*SHP(2,IY1) 184 : #*SHP(5,IX1+8)) 185 : REL(IX,IY) = REL(IY,IX) 186 : REL(IX+1,IY) = REL(IY,IX+1) 187 : 109 CONTINUE 188 : IY1=0 189 : 108 CONTINUE 190 : IX1=1 191 : IY1=0 192 : DO 110 IX=2+NBDL,LRE ,NBDL 193 : IX1=IX1 + 1 194 : DO 111 IY=3,IX ,NBDL 195 : IY1=IY1 + 1 196 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * D11*SHP(2,IX1)*SHP(5,IY1) 197 : REL(IY+1,IX) = REL(IY+1,IX) + DJAC*VKL23* 198 : #D11*SHP(2,IX1)*SHP(5,IY1+8) 199 : REL(IX,IY) = REL(IY,IX) 200 : REL(IX,IY+1) = REL(IY+1,IX) 201 : 111 CONTINUE 202 : IY1=0 203 : 110 CONTINUE 204 : C 205 : C TERMES EN PI*(UY,RX) 206 : C 207 : IX1=0 208 : IY1=0 209 : DO 114 IX=5,LRE ,NBDL 210 : IX1=IX1 + 1 211 : DO 115 IY=2,IX ,NBDL 212 : IY1=IY1 + 1 213 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * D22*SHP(3,IY1) 214 : #* SHP(5,IX1) 215 : REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL23 * 216 : #D22*SHP(3,IY1) * SHP(5,IX1+8)*(-1.D0) 217 : REL(IX,IY) = REL(IY,IX) 218 : REL(IX+1,IY) = REL(IY,IX+1) 219 : 115 CONTINUE 220 : IY1=0 221 : 114 CONTINUE 222 : IX1=1 223 : IY1=0 224 : DO 116 IX=2+NBDL,LRE ,NBDL 225 : IX1=IX1 + 1 226 : DO 117 IY=5,IX ,NBDL 227 : IY1=IY1 + 1 228 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * D22*SHP(3,IX1) 229 : #* SHP(5,IY1) 230 : REL(IY+1,IX) = REL(IY+1,IX) + DJAC*VKL23 * 231 : #D22*SHP(3,IX1) * SHP(5,IY1+8)*(-1.D0) 232 : REL(IX,IY) = REL(IY,IX) 233 : REL(IX,IY+1) = REL(IY+1,IX) 234 : 117 CONTINUE 235 : IY1=0 236 : 116 CONTINUE 237 : C 238 : C TERMES EN (UX,RY)*(UX,RY) 239 : C 240 : H11=RHOS + RHOF*(SFLU-B11) 241 : H22=RHOS + RHOF*(SFLU-B22) 242 : IX1=0 243 : IY1=0 244 : DO 112 IX=3,LRE ,NBDL 245 : IX1=IX1 + 1 246 : DO 113 IY=3,IX ,NBDL 247 : IY1=IY1 + 1 248 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1) 249 : REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1+8) 250 : REL(IY+1,IX) = REL(IY+1,IX)+DJAC*VKL33*H11*SHP(5,IY1+8)*SHP(5,IX1) 251 : REL(IY+1,IX+1)=REL(IY+1,IX+1)+DJAC*VKL33*H22*SHP(5,IY1+8) 252 : #*SHP(5,IX1+8) 253 : REL(IX,IY) = REL(IY,IX) 254 : REL(IX+1,IY) = REL(IY,IX+1) 255 : REL(IX,IY+1) = REL(IY+1,IX) 256 : REL(IX+1,IY+1) = REL(IY+1,IX+1) 257 : 113 CONTINUE 258 : IY1=0 259 : 112 CONTINUE 260 : C 261 : C TERMES EN (UY,RX)*(UY,RX) 262 : C 263 : H11=RHOS + RHOF*(SFLU-B11) 264 : H22=RHOS + RHOF*(SFLU-B22) 265 : IX1=0 266 : IY1=0 267 : DO 118 IX=5,LRE ,NBDL 268 : IX1=IX1 + 1 269 : DO 119 IY=5,IX ,NBDL 270 : IY1=IY1 + 1 271 : REL(IY,IX) = REL(IY,IX) + DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1) 272 : REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1+8) 273 : #*(-1.D0) 274 : REL(IY+1,IX) = REL(IY+1,IX)-DJAC*VKL33*H11*SHP(5,IY1+8)*SHP(5,IX1) 275 : REL(IY+1,IX+1)=REL(IY+1,IX+1)+DJAC*VKL33*H22*SHP(5,IY1+8) 276 : #*SHP(5,IX1+8) 277 : REL(IX,IY) = REL(IY,IX) 278 : REL(IX+1,IY) = REL(IY,IX+1) 279 : REL(IX,IY+1) = REL(IY+1,IX) 280 : REL(IX+1,IY+1) = REL(IY+1,IX+1) 281 : 119 CONTINUE 282 : IY1=0 283 : 118 CONTINUE 284 : GOTO 666 285 : C 286 : C MESSAGE D ERREUR : ELEMENT A SURFACE NULLE 287 : C 288 : 667 CONTINUE 289 : IRET = 2 290 : GOTO 666 291 : C 292 : 666 CONTINUE 293 : RETURN 294 : END 295 : 296 :
© Cast3M 2003 - All rights reserved.
Disclaimer