1 : C SHB8 SOURCE PV 11/03/07 21:18:13 6885 2 : SUBROUTINE SHB8(ICLE,XE,D,PROPEL,WORK,RE,OUT) 3 : IMPLICIT REAL *8(A-H,O-Z) 4 : 5 : * SUBROUTINE SHB8(ICLE,XE,D,PROPEL,WORK,RE,OUT) 6 : 7 : C 8 : C ELEMENT SHB8 9 : C 10 : 11 : INTEGER IOP,IBID,IDPLA(5),LAG,IRDC 12 : DIMENSION EYG(5),ENU(5) 13 : DIMENSION CMATTMP(6,6) 14 : DIMENSION XE(24),D(*),PROPEL(*),WORK(24),RE(24,24),OUT(*) 15 : DIMENSION XXG5(5),PXG5(5),XCOQ(3,4),BKSIP(3,8,5),B(3,8) 16 : DIMENSION XELOC(24),XCENT(3),PPP(3,3),PPPT(3,3) 17 : DIMENSION XL(3,4),XXX(3),YYY(3),BGL_LOC(6,24) 18 : DIMENSION BGL_LOCT(24,6),TMPTAB(6,24),TMPKE(24,24),CMATLOC(6,6) 19 : DIMENSION XELOCTMP(24) 20 : DIMENSION XXVB(3),HIJ(6),XKSTAB(24,24),TMPKE2(24,24) 21 : DIMENSION GB(8,4),GS(8,4),VB(8,3),XXGB(3,4) 22 : DIMENSION XK11(8,8),XK22(8,8),XK33(8,8),RR2(3,3),XK12(8,8) 23 : DIMENSION XK21(8,8),XK13(8,8),XK23(8,8),XK31(8,8),XK32(8,8) 24 : DIMENSION XR(3,3),XRT(3,3) 25 : C 26 : DIMENSION DEPS(6),DUSDX(9),UE(3,8),SIG(6),XE_LOCP(24) 27 : DIMENSION UDEF(24),XE_LOC(24),XE_LOC12(24) 28 : DIMENSION DEPS_LOC(6),SIG_LOC(6) 29 : C 30 : DIMENSION F(3,8),BSIG(24),SIGMAG(6),ULOC(3,8),PQIALPH(3,4) 31 : DIMENSION QIALPHA(3,4),FHG(3,8),RR12(3,3),RR1(3,3),FHG24(24) 32 : C 33 : DIMENSION DIREC(3),V(3),X56(3),X67(3),X78(3),X85(3),X12(3) 34 : DIMENSION X23(3),X34(3),X41(3),V56_78(3),V85_67(3),VSUP(3) 35 : DIMENSION V12_34(3),V41_23(3),VINF(3) 36 : C 37 : DIMENSION XKSIGTMP1(8,8),XKSIGTMP2(8,8) 38 : C 39 : DIMENSION XCOQ14(3),XCOQ23(3),XCOQ34(3),XCOQ12(3),XNORMALE(3) 40 : DIMENSION XCOQ12_34(3),XCOQ14_23(3),XKP(24,24),BKP(3,4) 41 : C 42 : DIMENSION XMASSE(8,8) 43 : C 44 : DIMENSION FQ(24) 45 : C 46 : C AVRIL 2004: PROGRAMMATION DU SHB8 AVEC B_2_CHAPEAU_BAR 47 : C 48 : DIMENSION CB2BAR(6,6),SIGMOY(6) 49 : C 50 : CCCCCCCCCCCCCC ENTREES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 51 : C ICLE=2 ON CALCULE LA MATRICE DE RAIDEUR 52 : C ICLE=3 ON CALCULE LA MATRICE DE MASSE 53 : C ICLE=4 ON CALCULE LA MATRICE TANGENTE 54 : C ICLE=5 ON CALCULE LES FORCES DE PRESSION 55 : C ICLE=7 ON CALCULE LES CONTRAINTES 56 : C ICLE=8 ON CALCULE BSIGMA 57 : C ICLE=9 ON CALCULE KSIGMA 58 : C ICLE=10 ON CALCULE KP 59 : C ICLE=11 On calcule epsi a partir de deplacement 60 : C icle=12 on calcule epsi à partir des contraintes 61 : C ou sigma à partir de epsi suivan propel(3) 62 : C 63 : C D MATRICE D'ELASTICITE 64 : C PROPEL PROPRIETES ELASTIQUES (1)=YUNG (2)=NU 65 : C (3)=RO (4)=ALPHA 66 : C WORK TABLEAU DE TRAVAIL 67 : CCCCCCCCCCCCCC SORTIE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 68 : C 69 : C RE MATRICE 70 : C OUT VECTEUR FORCE OU VECTEUR CONTRAINTE AUX NOEUDS 71 : C 72 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 73 : C 74 : C INITIALISATIONS 75 : C 76 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 77 : C INFOS: 78 : C XE EST RANGE COMME CA: (XNOEUD1 YNOEUD1 ZNOEUD1, XNOEUD2 YNOEUD2 ZNOEUD2,...) 79 : C DANS SHB8_TEST_NUM: ATTENTION A LA NUMEROTATION DES NOEUDS 80 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 81 : DATA GB/1.D0, 1.D0,-1.D0,-1.D0,-1.D0,-1.D0, 1.D0, 1.D0, 82 : & 1.D0,-1.D0,-1.D0, 1.D0,-1.D0, 1.D0, 1.D0,-1.D0, 83 : & 1.D0,-1.D0, 1.D0,-1.D0, 1.D0,-1.D0, 1.D0,-1.D0, 84 : & -1.D0, 1.D0,-1.D0, 1.D0, 1.D0,-1.D0, 1.D0,-1.D0/ 85 : C 86 : C VB: COORD DES NOEUDS DANS REPERE DE REFERENCE 87 : C 88 : DATA VB/-1.D0, 1.D0, 1.D0,-1.D0,-1.D0, 1.D0, 1.D0,-1.D0, 89 : & -1.D0,-1.D0, 1.D0, 1.D0,-1.D0,-1.D0, 1.D0, 1.D0, 90 : & -1.D0,-1.D0,-1.D0,-1.D0, 1.D0, 1.D0, 1.D0, 1.D0/ 91 : C 92 : C ON DEFINI LES POINTS GAUSS ET LES POIDS 93 : C 94 : XXG5(1) = -0.906179845938664D0 95 : XXG5(2) = -0.538469310105683D0 96 : XXG5(3) = 0.D0 97 : XXG5(4) = 0.538469310105683D0 98 : XXG5(5) = 0.906179845938664D0 99 : C 100 : PXG5(1) = 0.236926885056189D0 101 : PXG5(2) = 0.478628670499366D0 102 : PXG5(3) = 0.568888888888889D0 103 : PXG5(4) = 0.478628670499366D0 104 : PXG5(5) = 0.236926885056189D0 105 : C 106 : UNS8 = 1.D0/8.D0 107 : UNS3 = 1.D0/3.D0 108 : UNS5 = 1.D0/5.D0 109 : UNS6 = 1.D0/6.D0 110 : ZERO = 0.D0 111 : C 112 : C TYPE DE LOI DE COMPORTEMENT: 113 : C IRDC = 1 : SHB8 TYPE PLEXUS 114 : C IRDC = 2 : C.P. 115 : C IRDC = 3 : 3D COMPLETE 116 : C 117 : C IRDC = 1 118 : IRDC=nint(OUT(1)) 119 : IF(ICLE.EQ.2)THEN 120 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 121 : C C 122 : C ON CALCULE LA RAIDEUR : SORTIE DANS RE C 123 : C C 124 : C SI IETAN = 1 , ALORS ON CALCULE AUSSI C 125 : C LA MATRICE TANGENTE PLASTIQUE C 126 : C C 127 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 128 : IETAN = nint(PROPEL(13)) 129 : C 130 : C INTIALISATION LONGUEUR DES COTES 131 : C CALCUL DES COEFF D ELANCEMENT A METTRE DANS LA MATRICE DE CPT 132 : C 133 : XXL1 = 0. 134 : XXL2 = 0. 135 : TT1 = 0. 136 : TT2 = 0. 137 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 138 : C 139 : C STABILISATION ADAPTATIVE EN FONCTION DE LA DISTORTION DE L'ELEMENT 140 : C 141 : DO I=1,3 142 : C DISTANCE ENTRE 1 ET 5 (EPAISSEUR) 143 : TT1 = TT1+(XE(I+12)-XE(I))**2 144 : C DISTANCE ENTRE 3 ET 7 (EPAISSEUR) 145 : TT2 = TT2+(XE(I+18)-XE(I+6))**2 146 : C DISTANCE ENTRE 1 ET 2 147 : XXL1 = XXL1+(XE(I+3)-XE(I))**2 148 : C DISTANCE ENTRE 2 ET 3 149 : XXL2 = XXL2+(XE(I+6)-XE(I+3))**2 150 : ENDDO 151 : XXL1 = SQRT(XXL1) 152 : XXL2 = SQRT(XXL2) 153 : TT1 = 0.5*(SQRT(TT1)+SQRT(TT2)) 154 : COELA1 = 5./6. 155 : COELA2 = 5./6. 156 : C ELANCEMENT DANS DIRECTION 2 157 : ELT = 6.*TT1/XXL1 158 : IF(COELA1.GT.ELT) COELA1=ELT 159 : C ELANCEMENT DANS DIRECTION 1 160 : ELT = 6.*TT1/XXL2 161 : IF(COELA2.GT.ELT) COELA2=ELT 162 : C POUR L'INSTANT, ON NE MET PAS EN SERVICE: 163 : COELA1 = 1. 164 : COELA2 = 1. 165 : C 166 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 167 : C 168 : IF (IETAN.EQ.1) THEN 169 : EYG(1) = PROPEL(3) 170 : EYG(2) = PROPEL(4) 171 : EYG(3) = PROPEL(5) 172 : EYG(4) = PROPEL(6) 173 : EYG(5) = PROPEL(7) 174 : ENU(1) = PROPEL(8) 175 : ENU(2) = PROPEL(9) 176 : ENU(3) = PROPEL(10) 177 : ENU(4) = PROPEL(11) 178 : ENU(5) = PROPEL(12) 179 : DO I=1,5 180 : C SI IDPLA(I)=1, POINT GAUSS (I) PLASTIQUE 181 : IDPLA(I)=nint(RE(I,1)) 182 : ENDDO 183 : ENDIF 184 : C 184 : CALL ZDANUL(RE,576) 184 : CALL ZDANUL(CMATLOC,36) 184 : C CALL ZDANUL(CMATTMP,36) 185 : IELEM=nint(WORK(1)) 186 : C 187 : C PROCEDURE DE TEST DE NUMEROTATION DES ELEMENTS: MARCHE BIEN, MAIS DEJA FAIT EN 188 : C AMONT 189 : C 189 : C CALL SHB8_TEST_NUM(XE,IELEM) 190 : C 191 : C ON DEFINI CMATLOC: MATRICE DE COMPORTEMENT 192 : C 193 : XNU = PROPEL(2) 194 : XNU_B = XNU / (1-XNU) 195 : XLAMBDA = PROPEL(1)*PROPEL(2)/(1-PROPEL(2)*PROPEL(2)) 196 : XMU = 0.5 * PROPEL(1) / (1+PROPEL(2)) 197 : CMATLOC(1,1) = XLAMBDA+2*XMU 198 : CMATLOC(2,2) = XLAMBDA+2*XMU 199 : IF (IRDC.EQ.1) THEN 200 : C COMPORTEMENT SHB8 PLEXUS 201 : CMATLOC(3,3) = PROPEL(1) 202 : ENDIF 203 : C 204 : IF (IRDC.EQ.2) THEN 205 : C COMPORTEMENT C.P. 206 : CMATLOC(3,3) = 0 207 : ENDIF 208 : C 209 : CMATLOC(1,2) = XLAMBDA 210 : CMATLOC(2,1) = XLAMBDA 211 : CMATLOC(4,4) = XMU 212 : CMATLOC(5,5) = XMU 213 : CMATLOC(6,6) = XMU 214 : C 215 : IF (IRDC.EQ.3) THEN 216 : C COMPORTEMENT LOI TRIDIM MMC 3D 217 : XNU = PROPEL(2) 218 : XCOOEF = PROPEL(1)/((1+XNU)*(1-2*XNU)) 219 : CMATLOC(1,1) = (1-XNU)*XCOOEF 220 : CMATLOC(2,2) = (1-XNU)*XCOOEF 221 : CMATLOC(3,3) = (1-XNU)*XCOOEF 222 : CMATLOC(1,2) = XNU*XCOOEF 223 : CMATLOC(2,1) = XNU*XCOOEF 224 : CMATLOC(1,3) = XNU*XCOOEF 225 : CMATLOC(3,1) = XNU*XCOOEF 226 : CMATLOC(2,3) = XNU*XCOOEF 227 : CMATLOC(3,2) = XNU*XCOOEF 228 : CMATLOC(4,4) = (1-2*XNU)*0.5*XCOOEF 229 : CMATLOC(5,5) = (1-2*XNU)*0.5*XCOOEF 230 : CMATLOC(6,6) = (1-2*XNU)*0.5*XCOOEF 231 : ENDIF 232 : C 232 : C CALL SHIFTD(CMATLOC,CMATTMP,36) 233 : C 234 : C CALCUL DE BKSIP(3,8,IP) DANS REPERE DE REFERENCE 235 : C BKSIP(1,*,IP) = VECTEUR BX AU POINT GAUSS IP 236 : C BKSIP(2,*,IP) = VECTEUR BY AU POINT GAUSS IP 237 : C BKSIP(3,*,IP) = VECTEUR BZ AU POINT GAUSS IP 238 : C 238 : CALL CUBKSI(5,XXG5,BKSIP) 239 : C 240 : C DEBUT DE LA BOUCLE SUR LES 5 PTS GAUSS 241 : C 242 : DO IP=1,5 243 : C 244 : C SI POINT PLASTIQUE, ON ITERE AVEC LE MODULE TANGENT 245 : C 246 : C IF (IETAN.EQ.1) THEN 247 : C IF (IDPLA(IP).EQ.1) THEN 248 : C CMATLOC(1,1) = CMATTMP(1,1)*EYG(IP)/PROPEL(1) 249 : C CMATLOC(2,2) = CMATTMP(2,2)*EYG(IP)/PROPEL(1) 250 : C CMATLOC(3,3) = CMATTMP(3,3)*EYG(IP)/PROPEL(1) 251 : C CMATLOC(1,2) = CMATTMP(1,2)*EYG(IP)/PROPEL(1) 252 : C CMATLOC(2,1) = CMATTMP(2,1)*EYG(IP)/PROPEL(1) 253 : C CMATLOC(4,4) = CMATTMP(4,4)*EYG(IP)/PROPEL(1) 254 : C CMATLOC(5,5) = CMATTMP(5,5)*EYG(IP)/PROPEL(1) 255 : C CMATLOC(6,6) = CMATTMP(6,6)*EYG(IP)/PROPEL(1) 256 : C ELSE 257 : C CMATLOC(1,1) = CMATTMP(1,1) 258 : C CMATLOC(2,2) = CMATTMP(2,2) 259 : C CMATLOC(3,3) = CMATTMP(3,3) 260 : C CMATLOC(1,2) = CMATTMP(1,2) 261 : C CMATLOC(2,1) = CMATTMP(2,1) 262 : C CMATLOC(4,4) = CMATTMP(4,4) 263 : C CMATLOC(5,5) = CMATTMP(5,5) 264 : C CMATLOC(6,6) = CMATTMP(6,6) 265 : C ENDIF 266 : C ENDIF 267 : C 268 : C DEFINITION DES 4 POINTS COQUES 269 : C 270 : ZETA = XXG5(IP) 271 : ZLAMB = 0.5*(1.-ZETA) 272 : DO I=1,4 273 : DO J=1,3 274 : XCOQ(J,I) = ZLAMB*XE((I-1)*3+J) 275 : ! + (1.-ZLAMB)*XE((I-1+4)*3+J) 276 : END DO 277 : END DO 278 : C 279 : C CALCUL DE PPP 3*3 PASSAGE DE GLOBAL A LOCAL,COQUE 280 : C XCENT : COORD GLOBAL DU CENTRE DE L'ELEMENT 281 : C 281 : CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,XBID) 282 : C 283 : C CALCUL DE B EN GLOBAL 284 : C 285 : C ATTENTION A L'ORDRE DE EPSILON: 286 : C FARID DANS SON PAPIER: 11 22 33 12 13 23 287 : C HARID DANS PLEXUS: 11 22 33 12 23 13 288 : C 289 : C 290 : C ON RAJOUTE LES TERMES H1,X . G1 , H2,X . G2 291 : C ET H1,Y . G1 , H2,Y . G2 292 : C AVEC H1 = Y.Z H2 = X.Z 293 : C DONC H1,X =0 H2,X = Z 294 : C ET H1,Y =Z H2,Y = 0 295 : C 296 : C DONC IL NE RESTE PLUS QU'A CALCULER G1 ET G2, ET A AJOUTER A BKSIP 297 : C 297 : CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC) 298 : C 299 : C IL FAUT: EPS_LOCAL=BGLOB .U_NODAL_GLOBAL 300 : C ON CALCULE BGL_LOC LA MATRICE B(6,24) UGLOBAL ---> EPS LOCAL 301 : C 301 : CALL ASBGLB(BGL_LOC,B,PPP) 302 : C 303 : C IL FAUT TRANSPOSER BGL_LOC 304 : C 304 : CALL AEQBT(BGL_LOCT,BGL_LOC,6,24) 305 : C 306 : C CALCUL DE LA MATRICE TANGENTE PLASTIQUE: VOIR FICHIER SHB8_MAT_TGE_PLAS.FFF 307 : C 308 : C 309 : C IL NE RESTE PLUS QU'A FAIRE: BGL_LOCT * C * BGL_LOC 310 : C 310 : CALL ZDANUL(TMPTAB,144) 310 : CALL ZDANUL(TMPKE,576) 310 : CALL ZDANUL(TMPKE2,576) 310 : * CALL MULMAT(6,6,24,CMATLOC,BGL_LOC,TMPTAB) 310 : CALL MULMAT( TMPTAB,CMATLOC,BGL_LOC,6,24,6) 310 : * CALL MULMAT(24,6,24,BGL_LOCT,TMPTAB,TMPKE2) 310 : CALL MULMAT(TMPKE2,BGL_LOCT,TMPTAB,24,24,6) 311 : C 312 : C ASSEMBLAGE: KE=KE + POIDS*JACOBIAN*TMPKE 313 : C 314 : DO J=1,8 315 : DO I=1,24 316 : TMPKE(I,(J-1)*3+1)=TMPKE2(I,J) 317 : TMPKE(I,(J-1)*3+2)=TMPKE2(I,J+8) 318 : TMPKE(I,(J-1)*3+3)=TMPKE2(I,J+16) 319 : END DO 320 : END DO 320 : CALL ZDANUL(TMPKE2,576) 321 : DO I=1,8 322 : DO J=1,24 323 : TMPKE2((I-1)*3+1,J)=TMPKE(I,J) 324 : TMPKE2((I-1)*3+2,J)=TMPKE(I+8,J) 325 : TMPKE2((I-1)*3+3,J)=TMPKE(I+16,J) 326 : END DO 327 : END DO 328 : DO J=1,24 329 : DO I=1,24 330 : RE(I,J)=RE(I,J) + 4.*AJAC*PXG5(IP)*TMPKE2(I,J) 331 : END DO 332 : END DO 333 : END DO 334 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 335 : C C 336 : C MATRICE DE STABILISATION : PAS DE BOUCLE SUR LES POINTS DE GAUSS C 337 : C C 338 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 339 : C 340 : C ON DEFINIE LES CONSTANTS DE B_2_CHAPEAU_BAR 341 : C 341 : CALL ZDANUL(CB2BAR,36) 342 : 343 : CB2BAR(1,1) = 1.D0 344 : CB2BAR(1,2) = 1.D0 345 : CB2BAR(1,3) = ZERO 346 : CB2BAR(1,4) = ZERO 347 : CB2BAR(1,6) = ZERO 348 : C 349 : CB2BAR(2,1) = ZERO 350 : CB2BAR(2,2) = ZERO 351 : CB2BAR(2,3) = 1.D0 352 : CB2BAR(2,4) = 1.D0 353 : CB2BAR(2,6) = ZERO 354 : C 355 : CB2BAR(3,1) = ZERO 356 : CB2BAR(3,2) = ZERO 357 : CB2BAR(3,3) = ZERO 358 : CB2BAR(3,4) = ZERO 359 : CB2BAR(3,6) = ZERO 360 : C 361 : CB2BAR(4,1) = ZERO 362 : CB2BAR(4,2) = ZERO 363 : CB2BAR(4,3) = ZERO 364 : CB2BAR(4,4) = ZERO 365 : C 366 : CB2BAR(5,4) = ZERO 367 : CB2BAR(5,5) = ZERO 368 : CB2BAR(5,6) = 1.D0 369 : C 370 : CB2BAR(6,2) = ZERO 371 : CB2BAR(6,5) = ZERO 372 : CB2BAR(6,6) = 1.D0 373 : C 374 : C ON A BESOIN DES VECTEURS GAMMA(8, ALPHA=1 A 4) = GS(8,4) 375 : C 376 : DO I=1,24 377 : XELOC(I) = XE(I) 378 : END DO 379 : C 380 : C REMARQUE: ATTENTION, RR, MATRICE DU CORROTATIONNEL EST RANGEE PAR LIGNES! 381 : C 381 : CALL CUBROT(XELOC,RR2) 381 : CALL VECTRO(RR2,XELOC,1) 381 : CALL CUBEBB(XELOC,B,VOL) 382 : C 383 : C CALCUL DES FONCTIONS DE FORME GS 384 : C XXGB = X * GB 385 : C 386 : DO J=1,3 387 : DO IA=1,4 388 : XXGB(J,IA) = HOUXGB(XELOC(J),IA) 389 : END DO 390 : END DO 391 : C 392 : C GS = (BBB) * XXGB 393 : C 394 : DO I=1,8 395 : DO J=1,4 396 : GS(I,J)=0.D0 397 : END DO 398 : END DO 399 : DO J=1,3 400 : DO IA=1,4 401 : DO I=1,8 402 : GS(I,IA)=GS(I,IA) + B(J,I)*XXGB(J,IA) 403 : END DO 404 : END DO 405 : END DO 406 : C 407 : C GS = GB - GS 408 : C 409 : DO I=1,4 410 : DO J=1,8 411 : GS(J,I)= (GB(J,I) - GS(J,I))*UNS8 412 : END DO 413 : END DO 414 : C 415 : C CALCUL DE XXVB = X * VB 416 : C 417 : XXVB(1)= - XELOC(1) + XELOC(4) + XELOC(7) - XELOC(10) 418 : & - XELOC(13) + XELOC(16) + XELOC(19) - XELOC(22) 419 : XXVB(2)= - XELOC(2) - XELOC(5) + XELOC(8) + XELOC(11) 420 : & - XELOC(14) - XELOC(17) + XELOC(20) + XELOC(23) 421 : XXVB(3)= - XELOC(3) - XELOC(6) - XELOC(9) - XELOC(12) 422 : & + XELOC(15) + XELOC(18) + XELOC(21) + XELOC(24) 423 : C 424 : C CALCUL DES RELATIONS CONTRAINTES ET DEFORMATIONS GENERALISEES 425 : C 426 : HIJ(1) = UNS3*XXVB(2)*XXVB(3)/XXVB(1) 427 : HIJ(2) = UNS3*XXVB(1)*XXVB(3)/XXVB(2) 428 : HIJ(3) = UNS3*XXVB(2)*XXVB(1)/XXVB(3) 429 : HIJ(4) = UNS3*XXVB(3) 430 : HIJ(5) = UNS3*XXVB(1) 431 : HIJ(6) = UNS3*XXVB(2) 432 : C 433 : C CALCUL DES COEFS A METTRE DANS KIJ POUR COMPOSER KSTAB 434 : C 435 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 436 : C ANCIENNE PROG A FARID (PLEXUS): 437 : C XK11_1=(XLAMBDA+2*XMU)*HIJ(1)+XMU*HIJ(2) 438 : C XK11_2=UNS3*((XLAMBDA+2*XMU)*HIJ(1)+XMU*(HIJ(2)+HIJ(3))) 439 : C XK22_1=(XLAMBDA+2*XMU)*HIJ(2)+XMU*HIJ(1) 440 : C XK22_2=UNS3*((XLAMBDA+2*XMU)*HIJ(2)+XMU*(HIJ(1)+HIJ(3))) 441 : C XK33_1=XMU*(HIJ(1)+HIJ(2)) 442 : C XK33_2=UNS3*(PROPEL(1)*HIJ(3)+XMU*(HIJ(1)+HIJ(2))) 443 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 444 : C 445 : C ICI IL FAUT PRENDRE LA MOYENNE DES MODULES D'YOUNG TANGENTS 446 : C 447 : C POUR LE SHB8 PLASTIQUE ET NON_LINEAIRE, METTRE IETAN A 1 DANS RIGID.FF 448 : C SINON METTRE IETAN A 0 449 : C 450 : IF (IETAN.EQ.1) THEN 451 : XYOUNGT = (EYG(1)+EYG(2)+EYG(3)+EYG(4)+EYG(5))/5. 452 : ELSE 453 : XYOUNGT = PROPEL(1) 454 : ENDIF 455 : C 456 : C POUR RESOUDRE LE CISAILLEMENT TRANSVERSE: 457 : C 458 : XYOUNGT = XYOUNGT*0.5*(COELA1+COELA2) 459 : XLAMBDA = XYOUNGT*PROPEL(2)/(1-PROPEL(2)*PROPEL(2)) 460 : XMU = 0.5 * XYOUNGT / (1+PROPEL(2)) 461 : XLAMBDA2= XLAMBDA + 2.D0*XMU 462 : C ON MET LES LIGNES SUIVANTES EN COMMENTAIRE 463 : C 464 : C XK11_1 = (XLAMBDA+2.D0*XMU)*HIJ(1) 465 : C XK11_2 = UNS3*((XLAMBDA+2.D0*XMU)*HIJ(1)) 466 : C XK22_1 = (XLAMBDA+2.D0*XMU)*HIJ(2) 467 : C XK22_2 = UNS3*((XLAMBDA+2.D0*XMU)*HIJ(2)) 468 : C XK33_1 = 0. 469 : C 470 : XK11_1 = (CB2BAR(1,1)*(XLAMBDA2*CB2BAR(1,1) + XLAMBDA*CB2BAR(2,1)) 471 : &+ CB2BAR(2,1)*(XLAMBDA*CB2BAR(1,1) + XLAMBDA2*CB2BAR(2,1)) 472 : &+ XYOUNGT*CB2BAR(3,1)**2)*HIJ(1) + XMU*HIJ(2)*CB2BAR(4,1)**2 473 : XK11_2 = (CB2BAR(1,2)*(XLAMBDA2*CB2BAR(1,2) + XLAMBDA*CB2BAR(2,2)) 474 : &+ CB2BAR(2,2)*(XLAMBDA*CB2BAR(1,2) + XLAMBDA2*CB2BAR(2,2)) 475 : C &+ XYOUNGT*CB2BAR(4,2)**2)*HIJ(1)*UNS3 476 : &+ XYOUNGT*CB2BAR(3,2)**2)*HIJ(1)*UNS3 477 : &+ XMU*UNS3*(HIJ(2)*CB2BAR(4,2)**2 + HIJ(3)*CB2BAR(6,2)**2) 478 : C 479 : XK22_1 = (CB2BAR(1,3)*(XLAMBDA2*CB2BAR(1,3) + XLAMBDA*CB2BAR(2,3)) 480 : &+ CB2BAR(2,3)*(XLAMBDA*CB2BAR(1,3) + XLAMBDA2*CB2BAR(2,3)) 481 : &+ XYOUNGT*CB2BAR(3,3)**2)*HIJ(2) + XMU*HIJ(1)*CB2BAR(4,3)**2 482 : XK22_2 = (CB2BAR(1,4)*(XLAMBDA2*CB2BAR(1,4) + XLAMBDA*CB2BAR(2,4)) 483 : &+ CB2BAR(2,4)*(XLAMBDA*CB2BAR(1,4) + XLAMBDA2*CB2BAR(2,4)) 484 : &+ XYOUNGT*CB2BAR(3,4)**2)*HIJ(2)*UNS3 485 : &+ XMU*UNS3*(HIJ(1)*CB2BAR(4,4)**2 + HIJ(3)*CB2BAR(5,4)**2) 486 : C 487 : C XK33_1 = XMU*(HIJ(2)*CB2BAR(5,5)**2 + HIJ(1)*CB2BAR(6,5)) 488 : XK33_1 = XMU*(HIJ(2)*CB2BAR(5,5)**2 + HIJ(1)*CB2BAR(6,5)**2) 489 : XK33_2 = (CB2BAR(1,6)*(XLAMBDA2*CB2BAR(1,6) + XLAMBDA*CB2BAR(2,6)) 490 : &+ CB2BAR(2,6)*(XLAMBDA*CB2BAR(1,6) + XLAMBDA2*CB2BAR(2,6)) 491 : &+ XYOUNGT*CB2BAR(3,6)**2)*UNS3*HIJ(3) 492 : &+ XMU*UNS3*(HIJ(2)*CB2BAR(5,6)**2 + HIJ(1)*CB2BAR(6,6)**2) 493 : C 494 : C ON MET LES SUIVANTES EN COMMENTAIRE POUR ANNULER LES TERMES CROISES 495 : C 496 : C XK13_1 = XYOUNGT*CB2BAR(3,1)*HIJ(6) 497 : C XK13_2 = XMU*CB2BAR(6,5)*HIJ(6) 498 : C XK23_1 = XYOUNGT*CB2BAR(3,3)*HIJ(5) 499 : C XK23_2 = XMU*CB2BAR(5,5)*HIJ(5) 500 : C XK31_1 = XMU*CB2BAR(6,5)*HIJ(6) 501 : C XK31_2 = XYOUNGT*CB2BAR(3,1)*HIJ(6) 502 : C XK32_1 = XMU*CB2BAR(5,5)*HIJ(5) 503 : C XK32_2 = XYOUNGT*CB2BAR(3,3)*HIJ(5) 504 : C 505 : XK13_1 = ZERO 506 : XK13_2 = ZERO 507 : XK23_1 = ZERO 508 : XK23_2 = ZERO 509 : XK31_1 = ZERO 510 : XK31_2 = ZERO 511 : XK32_1 = ZERO 512 : XK32_2 = ZERO 513 : C 514 : C ON MET LES LIGNES SUIVANTES EN COMMENTAIRE 515 : C 516 : C IF (IRDC.EQ.1) THEN 517 : C COMPORTEMENT SHB8 PLEXUS 518 : C XK33_2 = XYOUNGT*HIJ(3)*UNS3 519 : C ENDIF 520 : C 521 : IF (IRDC.EQ.2) THEN 522 : C COMPORTEMENT C.P. 523 : XK33_2 = 0. 524 : ENDIF 525 : C 526 : IF (IRDC.EQ.3) THEN 527 : C COMPORTEMENT LOI TRIDIM MMC 3D 528 : XK11_1 = XCOOEF*(1-XNU)*HIJ(1) 529 : XK11_2 = UNS3*(XCOOEF*(1-XNU)*HIJ(1)) 530 : XK22_1 = XCOOEF*(1-XNU)*HIJ(2) 531 : XK22_2 = UNS3*(XCOOEF*(1-XNU)*HIJ(2)) 532 : XK33_1 = 0. 533 : XK33_2 = XCOOEF*(1-XNU)*HIJ(3)*UNS3 534 : ENDIF 535 : C 535 : CALL ZDANUL(XK11,64) 535 : CALL ZDANUL(XK22,64) 535 : CALL ZDANUL(XK12,64) 535 : CALL ZDANUL(XK21,64) 535 : CALL ZDANUL(XK13,64) 535 : CALL ZDANUL(XK23,64) 535 : CALL ZDANUL(XK31,64) 535 : CALL ZDANUL(XK32,64) 535 : CALL ZDANUL(XK33,64) 536 : C 537 : DO J=1,8 538 : DO I=1,8 539 : C 540 : C IL FAUT CALCULER K11 K22 K33 K13 K23 K31 K32 MATRICES 8*8 541 : C 542 : XK11(I,J)=XK11_1*GS(I,3)*GS(J,3)+XK11_2*GS(I,4)*GS(J,4) 543 : XK22(I,J)=XK22_1*GS(I,3)*GS(J,3)+XK22_2*GS(I,4)*GS(J,4) 544 : XK33(I,J)=XK33_1*GS(I,3)*GS(J,3)+XK33_2*GS(I,4)*GS(J,4) 545 : C 546 : XK13(I,J)=XK13_1*GS(I,3)*GS(J,1)+XK13_2*GS(I,1)*GS(J,3) 547 : XK23(I,J)=XK23_1*GS(I,3)*GS(J,2)+XK23_2*GS(I,2)*GS(J,3) 548 : XK31(I,J)=XK31_1*GS(I,3)*GS(J,1)+XK31_2*GS(I,1)*GS(J,3) 549 : XK32(I,J)=XK32_1*GS(I,3)*GS(J,2)+XK32_2*GS(I,2)*GS(J,3) 550 : 551 : END DO 552 : END DO 553 : C 554 : C ASSEMBLAGE DE KSTAB 555 : C 555 : CALL ASKSTB(XKSTAB,XK11,XK22,XK33,XK12,XK21,XK13,XK23, 556 : ! XK31,XK32) 557 : C 558 : C 559 : C REMISE EN ORDRE DE KSTAB 560 : C 561 : DO J=1,8 562 : DO I=1,24 563 : TMPKE(I,(J-1)*3+1)=XKSTAB(I,J) 564 : TMPKE(I,(J-1)*3+2)=XKSTAB(I,J+8) 565 : TMPKE(I,(J-1)*3+3)=XKSTAB(I,J+16) 566 : END DO 567 : END DO 567 : CALL ZDANUL(XKSTAB,576) 568 : DO I=1,8 569 : DO J=1,24 570 : XKSTAB((I-1)*3+1,J)=TMPKE(I,J) 571 : XKSTAB((I-1)*3+2,J)=TMPKE(I+8,J) 572 : XKSTAB((I-1)*3+3,J)=TMPKE(I+16,J) 573 : END DO 574 : END DO 575 : C 576 : C IL FAUT REPASSER DANS LE REPERE GLOBAL AVEC RR2^T . KSTAB . RR2 577 : C EN FAIT C'EST RR2. KSTAB .RR2^T CAR RR2 RANGEE PAR LIGNES 578 : C 578 : CALL ASKSGL(XKSTAB,RR2) 579 : C 580 : C AJOUT DE KSTAB A KE 581 : C 582 : DO J=1,24 583 : DO I=1,24 584 : RE(I,J)=RE(I,J)+XKSTAB(I,J) 585 : END DO 586 : END DO 587 : ENDIF 588 : C 589 : IF(ICLE.EQ.7.or.icle.eq.11)THEN 590 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 591 : C icle=11 calcul de epsi C 592 : C ON CALCULE LES CONTRAINTES : SORTIE DANS OUT(30) C 593 : C CONTRAINTES LOCALES DANS CHAQUE COUCHE C 594 : C SUR LA CONFIGURATION 1 C 595 : C LE DEPLACEMENT NODAL A L'AIR D'ETRE DANS WORK(1 A 24) C 596 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 597 : C 598 : C UE: INCREMENT DE DEPLACEMENT NODAL, REPERE GLOBAL 599 : C 600 : C XE: DEBUT DU PAS 601 : irdc=3 602 : DO J=1,8 603 : DO I=1,3 604 : UE(I,J)=WORK((J-1)*3+I) 605 : END DO 606 : END DO 607 : IGRDDEP = nint(D(1)) 608 : LAG = nint(PROPEL(3)) 609 : C 610 : C ON DEFINIT CMATLOC LOI MODIFIEE SHB8 611 : C 612 : if(icle.eq.7) then 612 : CALL ZDANUL(CMATLOC,36) 613 : XLAMBDA = PROPEL(1)*PROPEL(2)/(1-PROPEL(2)*PROPEL(2)) 614 : XMU = 0.5 * PROPEL(1) / (1+PROPEL(2)) 615 : CMATLOC(1,1) = XLAMBDA+2*XMU 616 : CMATLOC(2,2) = XLAMBDA+2*XMU 617 : IF (IRDC.EQ.1) THEN 618 : C COMPORTEMENT SHB8 PLEXUS 619 : CMATLOC(3,3) = PROPEL(1) 620 : ENDIF 621 : C 622 : IF (IRDC.EQ.2) THEN 623 : C COMPORTEMENT C.P. 624 : CMATLOC(3,3) = 0 625 : ENDIF 626 : C 627 : CMATLOC(1,2) = XLAMBDA 628 : CMATLOC(2,1) = XLAMBDA 629 : CMATLOC(4,4) = XMU 630 : CMATLOC(5,5) = XMU 631 : CMATLOC(6,6) = XMU 632 : C 633 : IF (IRDC.EQ.3) THEN 634 : C COMPORTEMENT LOI TRIDIM MMC 3D 635 : XNU = PROPEL(2) 636 : XCOOEF = PROPEL(1)/((1+XNU)*(1-2*XNU)) 637 : CMATLOC(1,1) = (1-XNU)*XCOOEF 638 : CMATLOC(2,2) = (1-XNU)*XCOOEF 639 : CMATLOC(3,3) = (1-XNU)*XCOOEF 640 : CMATLOC(1,2) = XNU*XCOOEF 641 : CMATLOC(2,1) = XNU*XCOOEF 642 : CMATLOC(1,3) = XNU*XCOOEF 643 : CMATLOC(3,1) = XNU*XCOOEF 644 : CMATLOC(2,3) = XNU*XCOOEF 645 : CMATLOC(3,2) = XNU*XCOOEF 646 : CMATLOC(4,4) = (1-2*XNU)*0.5*XCOOEF 647 : CMATLOC(5,5) = (1-2*XNU)*0.5*XCOOEF 648 : CMATLOC(6,6) = (1-2*XNU)*0.5*XCOOEF 649 : ENDIF 650 : endif 651 : C 652 : C CALCUL DE BKSIP(3,8,IP) DANS REPERE DE REFERENCE 653 : C BKSIP(1,*,IP) = VECTEUR BX AU POINT GAUSS IP 654 : C BKSIP(2,*,IP) = VECTEUR BY AU POINT GAUSS IP 655 : C BKSIP(3,*,IP) = VECTEUR BZ AU POINT GAUSS IP 656 : C 656 : CALL CUBKSI(5,XXG5,BKSIP) 657 : C 658 : DO IP=1,5 659 : C 660 : C DEFINITION DES 4 POINTS COQUES 661 : C 662 : ZETA = XXG5(IP) 663 : ZLAMB = 0.5*(1.-ZETA) 664 : DO J=1,3 665 : DO I=1,4 666 : XCOQ(J,I) = ZLAMB*XE((I-1)*3+J) 667 : ! + (1.-ZLAMB)*XE((I-1+4)*3+J) 668 : END DO 669 : END DO 670 : C 671 : C CALCUL DE PPP 3*3 PASSAGE DE GLOBAL A LOCAL,COQUE 672 : C XCENT : COORD GLOBAL DU CENTRE DE L'ELEMENT 673 : C 673 : CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,XBID) 674 : C 675 : C CALCUL DE B : U_GLOBAL ---> EPS_GLOBAL 676 : C 676 : CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC) 677 : C 678 : C CALCUL DE EPS DANS LE REPERE GLOBAL: 1 POUR DEFORMATIONS LINEAIRES 679 : C 2 POUR TERMES CARRES EN PLUS 680 : DO I=1,6 681 : DEPS(I)=0. 682 : ENDDO 683 : IF (LAG.EQ.1) THEN 684 : C ON AJOUTE LA PARTIE NON-LINEAIRE DE EPS 684 : CALL DSDX3D(2,B,UE,DEPS,DUSDX,8) 685 : ELSE 685 : CALL DSDX3D(1,B,UE,DEPS,DUSDX,8) 686 : ENDIF 687 : C 688 : C SORTIE DE DUSDX DANS PROPEL(1 A 9 * 5 PT DE GAUSS) 689 : C POUR UTILISATION ULTERIEURE DANS Q8PKCN_SHB8 690 : C 690 : CALL AEQBT(PPPT,PPP,3,3) 691 : RR12(1,1) = DUSDX(1) 692 : RR12(1,2) = DUSDX(2) 693 : RR12(1,3) = DUSDX(3) 694 : RR12(2,1) = DUSDX(4) 695 : RR12(2,2) = DUSDX(5) 696 : RR12(2,3) = DUSDX(6) 697 : RR12(3,1) = DUSDX(7) 698 : RR12(3,2) = DUSDX(8) 699 : RR12(3,3) = DUSDX(9) 699 : * CALL MULMAT(3,3,3,PPPT,RR12,RR2) 699 : CALL MULMAT(RR2,PPPT,RR12,3,3,3) 699 : * CALL MULMAT(3,3,3,RR2,PPP,RR12) 699 : CALL MULMAT(RR12,RR2,PPP,3,3,3) 700 : DUSDX(1) = RR12(1,1) 701 : DUSDX(2) = RR12(1,2) 702 : DUSDX(3) = RR12(1,3) 703 : DUSDX(4) = RR12(2,1) 704 : DUSDX(5) = RR12(2,2) 705 : DUSDX(6) = RR12(2,3) 706 : DUSDX(7) = RR12(3,1) 707 : DUSDX(8) = RR12(3,2) 708 : DUSDX(9) = RR12(3,3) 709 : DO I=1,9 710 : PROPEL(I+(IP-1)*9)=DUSDX(I) 711 : ENDDO 712 : DO I=1,6 713 : DEPS_LOC(I) = 0. 714 : SIG_LOC(I) = 0. 715 : ENDDO 715 : CALL CHRP3D(PPP,DEPS,DEPS_LOC,2) 716 : C 717 : C CALCUL DE SIGMA DANS LE REPERE LOCAL 718 : C 718 : * CALL MULMAT(6,6,1,CMATLOC,DEPS_LOC,SIG_LOC) 719 : if(icle.eq.7) then 719 : CALL MULMAT(SIG_LOC,CMATLOC,DEPS_LOC,6,1,6) 720 : endif 721 : C 722 : C CONTRAINTES ECRITES SOUS LA FORME: 723 : C [SIG] = [S_11, S_22, S_33, S_12, S_23, S_13] 724 : if(icle.eq.7) then 725 : DO I=1,6 726 : C ON LAISSE LES CONTRAINTES DANS LE REPERE LOCAL POUR LA PLASTICITE 727 : OUT((IP-1)*6+I)=SIG_LOC(I) 728 : END DO 729 : elseif(icle.eq.11)then 730 : DO I=1,6 731 : C ON LAISSE LES CONTRAINTES DANS LE REPERE LOCAL POUR LA PLASTICITE 732 : OUT((IP-1)*6+I)=DEPS_LOC(I) 733 : END DO 734 : endif 735 : END DO 736 : ENDIF 737 : IF(icle.eq.12) then 738 : C propel(3)=1 on veut epsi à partir de sigma 739 : C propel(3)=2 on veut sigma à partir de epsi 740 : C en entrée work contient soit epsilon soit sigma 741 : C work(1,2 ... contrainte 1 du point 1, contrainte 2 sdu point 1 .. 742 : C propel(1)=young 743 : C propel(2)=nu 744 : C propel(3) =1 si on veut des contraintes 745 : C =2 si on veut des epsilons 746 : C 747 : C out en sortie out(1,2 .. epsi d1 du point 1 , epsi 2 du point 1 .. 748 : irdc=3 749 : icas=nint(propel(3)) 750 : C 751 : C ON DEFINIT CMATLOC LOI MODIFIEE SHB8 752 : C 752 : CALL ZDANUL(CMATLOC,36) 753 : XLAMBDA = PROPEL(1)*PROPEL(2)/(1-PROPEL(2)*PROPEL(2)) 754 : XMU = 0.5 * PROPEL(1) / (1+PROPEL(2)) 755 : CMATLOC(1,1) = XLAMBDA+2*XMU 756 : CMATLOC(2,2) = XLAMBDA+2*XMU 757 : IF (IRDC.EQ.1) THEN 758 : C COMPORTEMENT SHB8 PLEXUS 759 : CMATLOC(3,3) = PROPEL(1) 760 : ENDIF 761 : C 762 : IF (IRDC.EQ.2) THEN 763 : C COMPORTEMENT C.P. 764 : CMATLOC(3,3) = 0 765 : ENDIF 766 : C 767 : CMATLOC(1,2) = XLAMBDA 768 : CMATLOC(2,1) = XLAMBDA 769 : CMATLOC(4,4) = XMU 770 : CMATLOC(5,5) = XMU 771 : CMATLOC(6,6) = XMU 772 : C 773 : IF (IRDC.EQ.3) THEN 774 : C COMPORTEMENT LOI TRIDIM MMC 3D 775 : XNU = PROPEL(2) 776 : XCOOEF = PROPEL(1)/((1+XNU)*(1-2*XNU)) 777 : CMATLOC(1,1) = (1-XNU)*XCOOEF 778 : CMATLOC(2,2) = (1-XNU)*XCOOEF 779 : CMATLOC(3,3) = (1-XNU)*XCOOEF 780 : CMATLOC(1,2) = XNU*XCOOEF 781 : CMATLOC(2,1) = XNU*XCOOEF 782 : CMATLOC(1,3) = XNU*XCOOEF 783 : CMATLOC(3,1) = XNU*XCOOEF 784 : CMATLOC(2,3) = XNU*XCOOEF 785 : CMATLOC(3,2) = XNU*XCOOEF 786 : CMATLOC(4,4) = (1-2*XNU)*0.5*XCOOEF 787 : CMATLOC(5,5) = (1-2*XNU)*0.5*XCOOEF 788 : CMATLOC(6,6) = (1-2*XNU)*0.5*XCOOEF 789 : ENDIF 790 : if(icas.eq.2) then 790 : call invalm(cmatloc,6,3 , ier,1.d-12) 791 : CMATLOC(4,4) =1.d0/CMATLOC(4,4) 792 : CMATLOC(5,5) =1.d0/CMATLOC(5,5) 793 : CMATLOC(6,6) =1.d0/CMATLOC(6,6) 794 : endif 795 : write(6,*) ' cmatloc ' 796 : write(6,*) (cmatloc(1,iop),iop=1,6) 797 : write(6,*) (cmatloc(2,iop),iop=1,6) 798 : write(6,*) (cmatloc(3,iop),iop=1,6) 799 : write(6,*) (cmatloc(4,iop),iop=1,6) 800 : write(6,*) (cmatloc(5,iop),iop=1,6) 801 : write(6,*) (cmatloc(6,iop),iop=1,6) 802 : write(6,*) 'epsi ' 803 : C 804 : DO IP=1,5 805 : C 806 : C DEFINITION DES 4 POINTS COQUES 807 : C 808 : write(6,*) (work((ip-1)*6 +iop),iop=1,6) 809 : C 810 : C CALCUL DE SIGMA DANS LE REPERE LOCAL 811 : C 811 : CALL MULMAT(out((ip-1)*6+1),CMATLOC,work((ip-1)*6+1),6,1,6) 812 : END DO 813 : ENDIF 814 : 815 : IF(ICLE.EQ.8)THEN 816 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 817 : C C 818 : C ON CALCULE BSIGMA: SORTIE DANS OUT(24) C 819 : C ENTREE DE SIGMA DANS WORK(DIM=30) C 820 : C ENTREE DU MATERIAU DANS D(1) ET D(2) C 821 : C ENTREE DE QIALPHA PAS PRECEDENT C 822 : C DANS RE(1 A 12) C 823 : C QIALPHA(J,I)=RE(J+(I-1)*3) C 824 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 825 : LAG = nint(D(3)) 826 : IRDC=3 826 : CALL ZDANUL(XKSIGTMP2,64) 827 : DO J=1,8 828 : DO I=1,3 829 : F(I,J) = 0.D0 830 : END DO 831 : END DO 832 : C 833 : C CALCUL DE BKSIP(3,8,IP) DANS REPERE DE REFERENCE 834 : C BKSIP(1,*,IP) = VECTEUR BX AU POINT GAUSS IP 835 : C BKSIP(2,*,IP) = VECTEUR BY AU POINT GAUSS IP 836 : C BKSIP(3,*,IP) = VECTEUR BZ AU POINT GAUSS IP 837 : C 837 : CALL CUBKSI(5,XXG5,BKSIP) 838 : DO IP=1,5 839 : C 840 : C RECHERCHE DE SIGMA DU POINT DE GAUSS GLOBAL 841 : C 842 : DO I=1,6 843 : C 844 : C C'EST LES CONTRAINTES LOCALES POUR POUVOIR TRAITER LA PLASTICITE AVANT 845 : C LES CONTRAINTES SONT PASSEES DANS LA CONFI. A LA FIN DU PAS DANS Q8PKCN 846 : C 847 : SIG_LOC(I)=WORK((IP-1)*6+I) 848 : END DO 849 : ZETA = XXG5(IP) 850 : ZLAMB = 0.5*(1.-ZETA) 851 : DO J=1,3 852 : DO I=1,4 853 : XCOQ(J,I) = ZLAMB*XE((I-1)*3+J) 854 : ! + (1.-ZLAMB)*XE((I-1+4)*3+J) 855 : END DO 856 : END DO 856 : CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,xBID) 857 : C 858 : C PASSAGE DES CONTRAINTES AU REPERE GLOBAL 859 : C 859 : CALL CHRP3D(PPP,SIG_LOC,SIGMAG,1) 860 : C 860 : CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC) 861 : C 862 : C CALCUL DE BQ.SIGMA SI LAGRANGIEN TOTAL 863 : C 864 : IF (LAG.EQ.1) THEN 864 : CALL ASKSIG(SIGMAG,B,XKSIGTMP1) 865 : DO J=1,8 866 : DO I=1,8 867 : XKSIGTMP2(I,J) = XKSIGTMP2(I,J) 868 : ! + 4D0*AJAC*PXG5(IP)*XKSIGTMP1(I,J) 869 : ENDDO 870 : ENDDO 871 : ENDIF 872 : C 873 : C CALCUL DE B.SIGMA EN GLOBAL 874 : C 875 : POIDS = 4D0 * AJAC * PXG5(IP) 876 : DO K=1,8 877 : F(1,K) = F(1,K) + POIDS * 878 : & (B(1,K)*SIGMAG(1)+B(2,K)*SIGMAG(4)+B(3,K)*SIGMAG(6)) 879 : F(2,K) = F(2,K) + POIDS * 880 : & (B(1,K)*SIGMAG(4)+B(2,K)*SIGMAG(2)+B(3,K)*SIGMAG(5)) 881 : F(3,K) = F(3,K) + POIDS * 882 : & (B(1,K)*SIGMAG(6)+B(2,K)*SIGMAG(5)+B(3,K)*SIGMAG(3)) 883 : END DO 884 : C 885 : C SI LAGRANGIEN TOTAL: AJOUT DE FQ A F 886 : C 887 : END DO 888 : IF (LAG.EQ.1) THEN 888 : CALL ZDANUL(TMPKE,576) 889 : DO KK=1,3 890 : DO I=1,8 891 : DO J=1,8 892 : TMPKE(I+(KK-1)*8,J+(KK-1)*8) = XKSIGTMP2(I,J) 893 : ENDDO 894 : ENDDO 895 : ENDDO 895 : CALL ZDANUL(TMPKE2,576) 896 : DO J=1,8 897 : DO I=1,24 898 : TMPKE2(I,(J-1)*3+1)=TMPKE(I,J) 899 : TMPKE2(I,(J-1)*3+2)=TMPKE(I,J+8) 900 : TMPKE2(I,(J-1)*3+3)=TMPKE(I,J+16) 901 : END DO 902 : END DO 902 : CALL ZDANUL(TMPKE,576) 903 : DO I=1,8 904 : DO J=1,24 905 : TMPKE((I-1)*3+1,J)=TMPKE2(I,J) 906 : TMPKE((I-1)*3+2,J)=TMPKE2(I+8,J) 907 : TMPKE((I-1)*3+3,J)=TMPKE2(I+16,J) 908 : END DO 909 : END DO 909 : * CALL MULMAT(24,24,1,TMPKE,PROPEL,FQ) 909 : CALL MULMAT(FQ,TMPKE,PROPEL,24,1,24) 910 : DO K=1,8 911 : F(1,K) = F(1,K) + FQ((K-1)*3+1) 912 : F(2,K) = F(2,K) + FQ((K-1)*3+2) 913 : F(3,K) = F(3,K) + FQ((K-1)*3+3) 914 : END DO 915 : ENDIF 916 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 917 : C 918 : C ON CALCULE LA STABILISATION POUR LE CALCUL DE B.SIGMA 919 : C UE EST DANS PROPEL! 920 : C ON A BESOIN DU MATERIAU AUSSI! 921 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 922 : C 923 : XYOUNG = D(1) 924 : XNU = D(2) 925 : XLAMBDA = XYOUNG*XNU/(1-XNU*XNU) 926 : XMU = 0.5 * XYOUNG / (1+XNU) 927 : C 928 : C DEPLACEMENT NODAL, REPERE GLOBAL : DANS PROPEL(1 A 24) 929 : C 930 : DO I=1,24 931 : XE_LOC(I) = XE(I) 932 : XE_LOCP(I) = XE(I)-PROPEL(I) 933 : XE_LOC12(I) = XE(I)-0.5*PROPEL(I) 934 : ENDDO 935 : C 936 : C REMARQUE: ATTENTION, RR, MATRICE DU CORROTATIONNEL EST RANGEE PAR LIGNES! 937 : C 937 : CALL CUBROT(XE_LOC,RR2) 937 : CALL CUBROT(XE_LOCP,RR1) 937 : CALL CUBROT(XE_LOC12,RR12) 937 : CALL VECTRO(RR2,XE_LOC,1) 937 : CALL VECTRO(RR1,XE_LOCP,1) 937 : CALL VECTRO(RR12,XE_LOC12,1) 938 : C 939 : C ON A BESOIN DES VECTEURS GAMMA(8, ALPHA=1 A 4) = GS(8,4) 940 : C ORIGINE DE XE_LOC12 PAS BONNE, MAIS PAS GRAVE : DIFFERENTIELLES DANS B 941 : C 941 : CALL CUBEBB(XE_LOC12,B,VOL) 942 : C 943 : C 944 : C CALCUL DEPLACEMENT DEFORMANT DANS REPERE 1/2 PAS: 945 : C 946 : DO I=1,24 947 : UDEF(I)=XE_LOC(I)-XE_LOCP(I) 948 : ENDDO 949 : C 950 : C CALCUL DES FONCTIONS DE FORME GS 951 : C XXGB = X * GB 952 : C 953 : DO IA=1,4 954 : DO J=1,3 955 : XXGB(J,IA) = HOUXGB(XE_LOC12(J),IA) 956 : END DO 957 : END DO 958 : C 959 : C GS = (BBB) * XXGB 960 : C 961 : DO J=1,4 962 : DO I=1,8 963 : GS(I,J)=0.D0 964 : END DO 965 : END DO 966 : DO J=1,3 967 : DO IA=1,4 968 : DO I=1,8 969 : GS(I,IA)=GS(I,IA) + B(J,I)*XXGB(J,IA) 970 : END DO 971 : END DO 972 : END DO 973 : C 974 : C GS = GB - GS 975 : C 976 : DO I=1,4 977 : DO J=1,8 978 : GS(J,I)= (GB(J,I) - GS(J,I))*UNS8 979 : END DO 980 : END DO 981 : C 982 : C CALCUL DE XXVB = X * VB 983 : C 984 : XXVB(1)= - XE_LOC12(1) + XE_LOC12(4) + XE_LOC12(7) - XE_LOC12(10) 985 : & - XE_LOC12(13) + XE_LOC12(16) + XE_LOC12(19) - XE_LOC12(22) 986 : XXVB(2)= - XE_LOC12(2) - XE_LOC12(5) + XE_LOC12(8) + XE_LOC12(11) 987 : & - XE_LOC12(14) - XE_LOC12(17) + XE_LOC12(20) + XE_LOC12(23) 988 : XXVB(3)= - XE_LOC12(3) - XE_LOC12(6) - XE_LOC12(9) - XE_LOC12(12) 989 : & + XE_LOC12(15) + XE_LOC12(18) + XE_LOC12(21) + XE_LOC12(24) 990 : C 991 : C CALCUL DES RELATIONS CONTRAINTES ET DEFORMATIONS GENERALISEES 992 : C 993 : HIJ(1) = UNS3*XXVB(2)*XXVB(3) / XXVB(1) 994 : HIJ(2) = UNS3*XXVB(1)*XXVB(3) / XXVB(2) 995 : HIJ(3) = UNS3*XXVB(2)*XXVB(1) / XXVB(3) 996 : HIJ(4) = UNS3*XXVB(3) 997 : HIJ(5) = UNS3*XXVB(1) 998 : HIJ(6) = UNS3*XXVB(2) 999 : C 1000 : C CALCUL DES DEFORMATIONS GENERALISEES 1001 : C 1002 : DO IA=1,4 1003 : DO J=1,3 1004 : AUX=0.D0 1005 : DO I=1,8 1006 : AUX = AUX + GS(I,IA)*UDEF((I-1)*3+J) 1007 : END DO 1008 : PQIALPH(J,IA) = AUX 1009 : END DO 1010 : END DO 1011 : C 1012 : C CALCUL DES CONTRAINTES GENERALISEES 1013 : C 1014 : DO I=1,4 1015 : DO J=1,3 1016 : QIALPHA(J,I)=RE(J+(I-1)*3,1) 1017 : END DO 1018 : END DO 1019 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1020 : C PROGRAMMATION FARID (PLEXUS) C 1021 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1022 : C QIALPHA(1,1) = QIALPHA(1,1) 1023 : C QIALPHA(2,2) = QIALPHA(2,2) 1024 : C QIALPHA(3,3) = QIALPHA(3,3) + XMU*((HIJ(2)+HIJ(1))*PQIALPH(3,3)) 1025 : C QIALPHA(1,2) = QIALPHA(1,2) 1026 : C QIALPHA(2,3) = QIALPHA(2,3) + ((XLAMBDA+2.D0*XMU)*HIJ(2) 1027 : C ! + XMU*HIJ(1))*PQIALPH(2,3) 1028 : C QIALPHA(1,3) = QIALPHA(1,3) + ((XLAMBDA+2.D0*XMU)*HIJ(1) 1029 : C ! + XMU*HIJ(2))*PQIALPH(1,3) 1030 : C QIALPHA(2,1) = QIALPHA(2,1) 1031 : C QIALPHA(3,2) = QIALPHA(3,2) 1032 : C QIALPHA(3,1) = QIALPHA(3,1) 1033 : C QIALPHA(1,4) = QIALPHA(1,4) + UNS3*((XLAMBDA+2.D0*XMU)*HIJ(1) 1034 : C ! + XMU*(HIJ(2)+HIJ(3)))*PQIALPH(1,4) 1035 : C QIALPHA(2,4) = QIALPHA(2,4) + UNS3*((XLAMBDA+2.D0*XMU)*HIJ(2) 1036 : C ! + XMU*(HIJ(1)+HIJ(3)))*PQIALPH(2,4) 1037 : C QIALPHA(3,4) = QIALPHA(3,4) + UNS3*(XYOUNG*HIJ(3) 1038 : C ! + XMU*(HIJ(1)+HIJ(2)))*PQIALPH(3,4) 1039 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1040 : IF (IRDC.EQ.1) THEN 1041 : QIALPHA(1,1) = QIALPHA(1,1) 1042 : QIALPHA(2,2) = QIALPHA(2,2) 1043 : QIALPHA(3,3) = QIALPHA(3,3) 1044 : QIALPHA(1,2) = QIALPHA(1,2) 1045 : QIALPHA(2,3) = QIALPHA(2,3) 1046 : ! + ((XLAMBDA+2.D0*XMU)*HIJ(2))*PQIALPH(2,3) 1047 : QIALPHA(1,3) = QIALPHA(1,3) 1048 : ! + ((XLAMBDA+2.D0*XMU)*HIJ(1))*PQIALPH(1,3) 1049 : QIALPHA(2,1) = QIALPHA(2,1) 1050 : QIALPHA(3,2) = QIALPHA(3,2) 1051 : QIALPHA(3,1) = QIALPHA(3,1) 1052 : QIALPHA(1,4) = QIALPHA(1,4) 1053 : ! + UNS3*((XLAMBDA+2*XMU)*HIJ(1))*PQIALPH(1,4) 1054 : QIALPHA(2,4) = QIALPHA(2,4) 1055 : ! + UNS3*((XLAMBDA+2*XMU)*HIJ(2))*PQIALPH(2,4) 1056 : C PROJECTION LEGAY POUR LE TERME QIALPHA(3,4) 1057 : C COMPORTEMENT SHB8 PLEXUS 1058 : CC QIALPHA(3,4) = QIALPHA(3,4) + XYOUNG*HIJ(3)*UNS3*PQIALPH(3,4) 1059 : C PROJECTION 28 LE TERME QIALPHA(3,4) EST MODIFIE 1060 : QIALPHA(3,4)=QIALPHA(3,4)+UNS3*(HIJ(1)+HIJ(2))*XMU*PQIALPH(3,4) 1061 : ENDIF 1062 : C A MODIFIER POUR LA PROJECTION 28 (CAS DES C.P.) 1063 : IF (IRDC.EQ.2) THEN 1064 : QIALPHA(1,1) = QIALPHA(1,1) 1065 : QIALPHA(2,2) = QIALPHA(2,2) 1066 : QIALPHA(3,3) = QIALPHA(3,3) 1067 : QIALPHA(1,2) = QIALPHA(1,2) 1068 : QIALPHA(2,3) = QIALPHA(2,3) 1069 : ! + ((XLAMBDA+2.D0*XMU)*HIJ(2))*PQIALPH(2,3) 1070 : QIALPHA(1,3) = QIALPHA(1,3) 1071 : ! + ((XLAMBDA+2.D0*XMU)*HIJ(1))*PQIALPH(1,3) 1072 : QIALPHA(2,1) = QIALPHA(2,1) 1073 : QIALPHA(3,2) = QIALPHA(3,2) 1074 : QIALPHA(3,1) = QIALPHA(3,1) 1075 : QIALPHA(1,4) = QIALPHA(1,4) 1076 : ! + UNS3*((XLAMBDA+2*XMU)*HIJ(1))*PQIALPH(1,4) 1077 : QIALPHA(2,4) = QIALPHA(2,4) 1078 : ! + UNS3*((XLAMBDA+2*XMU)*HIJ(2))*PQIALPH(2,4) 1079 : C COMPORTEMENT C.P. 1080 : QIALPHA(3,4) = QIALPHA(3,4) 1081 : ENDIF 1082 : IF (IRDC.EQ.3) THEN 1083 : C COMPORTEMENT LOI TRIDIM MMC 3D 1084 : C A MODIFIER POUR LA PROJECTION 28 (CAS DE LA LOI 3D) 1085 : XCOOEF = XYOUNG/((1+XNU)*(1-2*XNU)) 1086 : QIALPHA(1,1) = QIALPHA(1,1) 1087 : QIALPHA(2,2) = QIALPHA(2,2) 1088 : QIALPHA(3,3) = QIALPHA(3,3) 1089 : QIALPHA(1,2) = QIALPHA(1,2) 1090 : QIALPHA(2,3) = QIALPHA(2,3) 1091 : ! + ((XCOOEF*(1-XNU))*HIJ(2))*PQIALPH(2,3) 1092 : QIALPHA(1,3) = QIALPHA(1,3) 1093 : ! + ((XCOOEF*(1-XNU))*HIJ(1))*PQIALPH(1,3) 1094 : QIALPHA(2,1) = QIALPHA(2,1) 1095 : QIALPHA(3,2) = QIALPHA(3,2) 1096 : QIALPHA(3,1) = QIALPHA(3,1) 1097 : QIALPHA(1,4) = QIALPHA(1,4) 1098 : ! + UNS3*((XCOOEF*(1-XNU))*HIJ(1))*PQIALPH(1,4) 1099 : QIALPHA(2,4) = QIALPHA(2,4) 1100 : ! + UNS3*((XCOOEF*(1-XNU))*HIJ(2))*PQIALPH(2,4) 1101 : QIALPHA(3,4) = QIALPHA(3,4) 1102 : ! + XCOOEF*(1-XNU)*HIJ(3)*UNS3*PQIALPH(3,4) 1103 : ENDIF 1104 : CCCCCCCCCCCCCCCCCCC 1105 : C 1106 : C SAUVEGARDE DES FORCES DE STABILISATION 1107 : C 1108 : DO I=1,4 1109 : DO J=1,3 1110 : RE(J+(I-1)*3,1)=QIALPHA(J,I) 1111 : END DO 1112 : END DO 1113 : C 1114 : C CALCUL DES FORCES DE HOURGLASS FHG 1115 : C 1116 : DO I=1,8 1117 : DO J=1,3 1118 : FHG(J,I)=0.D0 1119 : ENDDO 1120 : END DO 1121 : DO J=1,3 1122 : DO I=1,8 1123 : DO IA=1,4 1124 : FHG(J,I) = FHG(J,I) + QIALPHA(J,IA)*GS(I,IA) 1125 : END DO 1126 : END DO 1127 : END DO 1128 : C ON REPASSE AU REPERE GLOBAL 1129 : DO I=1,3 1130 : DO J=1,8 1131 : FHG24((J-1)*3+I)=FHG(I,J) 1132 : ENDDO 1133 : ENDDO 1133 : CALL VECTRO(RR12,FHG24,2) 1134 : C 1135 : C AJOUT DE LA STABILISATION AU B SIGMA DEJA CALCULE 1136 : C 1137 : DO J=1,3 1138 : DO I=1,8 1139 : F(J,I) = F(J,I) + FHG24((I-1)*3+J) 1140 : ENDDO 1141 : ENDDO 1142 : C 1143 : C ATTENTION A L'ORDRE DE OUT 1144 : C 1145 : DO I=1,3 1146 : DO J=1,8 1147 : OUT((J-1)*3+I) = F(I,J) 1148 : END DO 1149 : END DO 1150 : ENDIF 1151 : IF(ICLE.EQ.5)THEN 1152 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1153 : C C 1154 : C ON CALCULE LES FORCES DE PRESSION C 1155 : C C 1156 : C ENTREES: PROPEL(1) : VALEUR DE LA PRESSION C 1157 : C XE : GEOMETRIE ELEMENT C 1158 : C WORK(1 A 3) : DIRECTION DE PRESSION C 1159 : C PROPEL(2) : 0 RECHERCHE DE LA FACE C 1160 : C 1 OU 2 : FACE INF. OU SUP. (POUR C 1161 : C PRESSION SUIVEUSE) C 1162 : C SORTIE: OUT FORCE DE PRESSION C 1163 : C WORK(4) : NUMERO DE FACE POUR LA PRESSION C 1164 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1165 : DO I=1,8 1166 : DO J=1,3 1167 : F(J,I)=0. 1168 : END DO 1169 : END DO 1170 : C 1171 : C LECTURE: PRESSION A METTRE AUX 4 NOEUDS, FACE ET DIRECTION 1172 : C 1173 : PRESS=PROPEL(1) 1174 : IFACE=nint(PROPEL(2)) 1175 : IF(IFACE.EQ.0)THEN 1176 : DO I=1,3 1177 : DIREC(I)=WORK(I) 1178 : END DO 1179 : ENDIF 1180 : C 1181 : C COORD. DES POINTS MILIEUX 1182 : C 1183 : DO J=1,3 1184 : X56(J) =( XE((5-1)*3+J)+XE((6-1)*3+J) )/2. 1185 : X67(J) =( XE((6-1)*3+J)+XE((7-1)*3+J) )/2. 1186 : X78(J) =( XE((7-1)*3+J)+XE((8-1)*3+J) )/2. 1187 : X85(J) =( XE((8-1)*3+J)+XE((5-1)*3+J) )/2. 1188 : X12(J) =( XE((1-1)*3+J)+XE((2-1)*3+J) )/2. 1189 : X23(J) =( XE((2-1)*3+J)+XE((3-1)*3+J) )/2. 1190 : X34(J) =( XE((3-1)*3+J)+XE((4-1)*3+J) )/2. 1191 : X41(J) =( XE((4-1)*3+J)+XE((1-1)*3+J) )/2. 1192 : END DO 1193 : C 1194 : C NORMALE RENTRANTE FACE SUP: V56_78 ^ V85_67 1195 : C 1196 : IF(IFACE.EQ.0.OR.IFACE.EQ.2)THEN 1197 : DO J=1,3 1198 : V56_78(J)=X78(J)-X56(J) 1199 : V85_67(J)=X67(J)-X85(J) 1200 : END DO 1201 : VSUP(1) = V56_78(2) * V85_67(3) - V56_78(3) * V85_67(2) 1202 : VSUP(2) = V56_78(3) * V85_67(1) - V56_78(1) * V85_67(3) 1203 : VSUP(3) = V56_78(1) * V85_67(2) - V56_78(2) * V85_67(1) 1204 : ENDIF 1205 : C 1206 : C NORMALE RENTRANTE FACE INF: V41_23 ^ V12_34 1207 : C 1208 : IF(IFACE.EQ.0.OR.IFACE.EQ.1)THEN 1209 : DO J=1,3 1210 : V12_34(J)=X34(J)-X12(J) 1211 : V41_23(J)=X23(J)-X41(J) 1212 : END DO 1213 : VINF(1) = V41_23(2) * V12_34(3) - V41_23(3) * V12_34(2) 1214 : VINF(2) = V41_23(3) * V12_34(1) - V41_23(1) * V12_34(3) 1215 : VINF(3) = V41_23(1) * V12_34(2) - V41_23(2) * V12_34(1) 1216 : ENDIF 1217 : IF(IFACE.EQ.0) THEN 1218 : C 1219 : C RECHERCHE DE LA BONNE FACE: PREMIER PASSAGE 1220 : C 1221 : SCALSUP = 0. 1222 : SCALINF = 0. 1223 : DO I=1,3 1224 : SCALSUP = SCALSUP+VSUP(I)*DIREC(I) 1225 : SCALINF = SCALINF+VINF(I)*DIREC(I) 1226 : END DO 1227 : C 1228 : C TEST ERREUR 1229 : C 1230 : IF((SCALSUP*SCALINF).GE.0)THEN 1231 : IF ((SCALSUP*SCALINF).EQ.0)THEN 1232 : WRITE(6,*)'** ERREUR SUR LES SURFACES **' 1233 : WRITE(6,*)'** DE SHB8 POUR PRESSION **' 1234 : WRITE(6,*)'** PRODUIT SCALAIRE ÉGAL A ZERO **' 1235 : WRITE(6,*)'** PROBLÈME DE DIRECTION DE PRESSION **' 1236 : ENDIF 1237 : IF ((SCALSUP*SCALINF).GT.0)THEN 1238 : WRITE(6,*)'** ERREUR SUR LES SURFACES **' 1239 : WRITE(6,*)'** DE SHB8 POUR PRESSION **' 1240 : WRITE(6,*)'** LES 2 FACES SONT POSSIBLES **' 1241 : ENDIF 1242 : STOP 1243 : ENDIF 1244 : C 1245 : C SI SCAL <0 : PRESSION SUR CETTE FACE 1246 : C 1247 : IF(SCALSUP.LT.0)THEN 1248 : WORK(4) = 1 1249 : IFACE = 1 1250 : ENDIF 1251 : IF(SCALINF.LT.0)THEN 1252 : WORK(4) = 2 1253 : IFACE = 2 1254 : ENDIF 1255 : ENDIF 1256 : IF(IFACE.EQ.2)THEN 1257 : C 1258 : C BOUCLE SUR LES 4 NOEUDS DE LA FACE SUPERIEURE 1259 : C 1260 : DO I=5,8 1261 : DO J=1,3 1262 : F(J,I)= PRESS * VSUP(J) / 4. 1263 : END DO 1264 : END DO 1265 : ENDIF 1266 : IF(IFACE.EQ.1)THEN 1267 : C 1268 : C BOUCLE SUR LES 4 NOEUDS DE LA FACE INFERIEURE 1269 : C 1270 : DO I=1,4 1271 : DO J=1,3 1272 : F(J,I)= PRESS * VINF(J) / 4. 1273 : END DO 1274 : END DO 1275 : ENDIF 1276 : C 1277 : C ATTENTION A L'ORDRE DE OUT 1278 : C 1279 : DO J=1,8 1280 : DO I=1,3 1281 : OUT((J-1)*3+I)=F(I,J) 1282 : END DO 1283 : END DO 1284 : ENDIF 1285 : IF(ICLE.EQ.9)THEN 1286 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1287 : C CALCUL DE K.SIGMA C 1288 : C EN ENTREE DANS WORK : SIGMA NORMALEMENT LONGUEUR 30 C 1289 : C PROPEL(1): 1 POUR RE=RE+KSIGMA C 1290 : C PROPEL(1): 0 POUR RE=KSIGMA C 1291 : C SORTIE : RE(24*24)=RE(24*24)+KSIGMA C 1292 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1293 : C 1294 : C CALCUL DE B (1 2 3) AUX 5 POINTS DE GAUSS 1295 : C 1295 : CALL CUBKSI(5,XXG5,BKSIP) 1296 : DO J=1,8 1297 : DO I=1,8 1298 : XKSIGTMP2(I,J)=0. 1299 : ENDDO 1300 : ENDDO 1301 : C 1302 : C 1303 : C CALCUL DE LA CONTRAINTE MOYENNE SUR LES 5 POINTS DE GAUSS POUR LE 1304 : C KSIGMA DE STABILISATION CALCULE PLUS LOIN 1305 : C 1306 : DO J=1,6 1307 : SIGMOY(J)=0.D0 1308 : ENDDO 1309 : C 1310 : DO J=1,6 1311 : DO IP=1,5 1312 : SIGMOY(J)=SIGMOY(J)+WORK((IP-1)*6+J) 1313 : ENDDO 1314 : SIGMOY(J)=UNS5*SIGMOY(J) 1315 : ENDDO 1316 : C 1317 : C DEBUT DE LA BOUCLE SUR LES 5 PTS GAUSS 1318 : C 1319 : DO IP=1,5 1320 : C 1321 : C CALCUL DE B 1322 : C 1322 : CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC) 1323 : C 1324 : C CALCUL DE MATRICE DE PASSAGE POUR POUVOIR CALCULER LES CONTRAINTES DANS 1325 : C LE REPERE GLOBAL 1326 : C 1327 : DO I=1,6 1328 : C C'EST LES CONTRAINTES LOCALES POUR POUVOIR TRAITER LA PLASTICITE AVANT 1329 : SIG_LOC(I)=WORK((IP-1)*6+I) 1330 : END DO 1331 : ZETA = XXG5(IP) 1332 : ZLAMB = 0.5*(1.-ZETA) 1333 : DO J=1,3 1334 : DO I=1,4 1335 : XCOQ(J,I) = ZLAMB*XE((I-1)*3+J) 1336 : ! + (1.-ZLAMB)*XE((I-1+4)*3+J) 1337 : END DO 1338 : END DO 1338 : CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,XBID) 1339 : C 1340 : C PASSAGE DES CONTRAINTES AU REPERE GLOBAL 1341 : C 1341 : CALL CHRP3D(PPP,SIG_LOC,SIGMAG,1) 1341 : CALL ASKSIG(SIGMAG,B,XKSIGTMP1) 1342 : DO J=1,8 1343 : DO I=1,8 1344 : XKSIGTMP2(I,J) = XKSIGTMP2(I,J) 1345 : ! + 4D0*AJAC*PXG5(IP)*XKSIGTMP1(I,J) 1346 : ENDDO 1347 : ENDDO 1348 : ENDDO 1348 : CALL ZDANUL(TMPKE,576) 1349 : DO KK=1,3 1350 : DO I=1,8 1351 : DO J=1,8 1352 : TMPKE(I+(KK-1)*8,J+(KK-1)*8) = XKSIGTMP2(I,J) 1353 : ENDDO 1354 : ENDDO 1355 : ENDDO 1356 : C 1357 : C ON MET DE L'ORDRE: 1358 : C 1358 : CALL ZDANUL(TMPKE2,576) 1359 : DO J=1,8 1360 : DO I=1,24 1361 : TMPKE2(I,(J-1)*3+1)=TMPKE(I,J) 1362 : TMPKE2(I,(J-1)*3+2)=TMPKE(I,J+8) 1363 : TMPKE2(I,(J-1)*3+3)=TMPKE(I,J+16) 1364 : END DO 1365 : END DO 1365 : CALL ZDANUL(TMPKE,576) 1366 : DO I=1,8 1367 : DO J=1,24 1368 : TMPKE((I-1)*3+1,J)=TMPKE2(I,J) 1369 : TMPKE((I-1)*3+2,J)=TMPKE2(I+8,J) 1370 : TMPKE((I-1)*3+3,J)=TMPKE2(I+16,J) 1371 : END DO 1372 : END DO 1373 : IPROPEL=nint(PROPEL(1)) 1373 : IF(IPROPEL.EQ.0) CALL SHIFTD(TMPKE,RE,576) 1373 : IF(IPROPEL.EQ.1) CALL AEQAPB(RE,TMPKE,576) 1374 : C 1375 : C 1376 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1377 : C C 1378 : C MATRICE DE STABILISATION POUR K.SIGMA: PAS DE BOUCLE SUR LES POINTS DE GAUSS C 1379 : C C 1380 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1381 : C 1382 : C 1383 : C 1384 : C ON A BESOIN DES VECTEURS GAMMA(8, ALPHA=1 A 4) = GS(8,4) 1385 : C 1386 : CCC DO I=1,24 1387 : CCC XELOC(I) = XE(I) 1388 : CCC END DO 1389 : C 1390 : C REMARQUE: ATTENTION, RR, MATRICE DU CORROTATIONNEL EST RANGEE PAR LIGNES! 1391 : C 1391 : CCC CALL CUBEROT(XELOC,RR2) 1391 : CCC CALL VECTROT(RR2,XELOC,1) 1391 : CCC CALL CUBE_BB(XELOC,B,VOL) 1392 : C 1393 : C CALCUL DES FONCTIONS DE FORME GS 1394 : C XXGB = X * GB 1395 : C 1396 : CCC DO J=1,3 1397 : CCC DO IA=1,4 1398 : CCC XXGB(J,IA) = HOUXGB(XELOC(J),IA) 1399 : CCC END DO 1400 : CCC END DO 1401 : C 1402 : C GS = (BBB) * XXGB 1403 : C 1404 : CCC DO I=1,8 1405 : CCC DO J=1,4 1406 : CCC GS(I,J)=0.D0 1407 : CCC END DO 1408 : CCC END DO 1409 : CCC DO J=1,3 1410 : CCC DO IA=1,4 1411 : CCC DO I=1,8 1412 : CCC GS(I,IA)=GS(I,IA) + B(J,I)*XXGB(J,IA) 1413 : CCC END DO 1414 : CCC END DO 1415 : CCC END DO 1416 : C 1417 : C GS = GB - GS 1418 : C 1419 : CCC DO I=1,4 1420 : CCC DO J=1,8 1421 : CCC GS(J,I)= (GB(J,I) - GS(J,I))*UNS8 1422 : CCC END DO 1423 : CCC END DO 1424 : C 1425 : C CALCUL DE XXVB = X * VB 1426 : C 1427 : CCC XXVB(1)= - XELOC(1) + XELOC(4) + XELOC(7) - XELOC(10) 1428 : CCC & - XELOC(13) + XELOC(16) + XELOC(19) - XELOC(22) 1429 : CCC XXVB(2)= - XELOC(2) - XELOC(5) + XELOC(8) + XELOC(11) 1430 : CCC & - XELOC(14) - XELOC(17) + XELOC(20) + XELOC(23) 1431 : CCC XXVB(3)= - XELOC(3) - XELOC(6) - XELOC(9) - XELOC(12) 1432 : CCC & + XELOC(15) + XELOC(18) + XELOC(21) + XELOC(24) 1433 : C 1434 : C CALCUL DES RELATIONS CONTRAINTES ET DEFORMATIONS GENERALISEES 1435 : C 1436 : CCC HIJ(1) = UNS3*XXVB(2)*XXVB(3)/XXVB(1) 1437 : CCC HIJ(2) = UNS3*XXVB(1)*XXVB(3)/XXVB(2) 1438 : CCC HIJ(3) = UNS3*XXVB(2)*XXVB(1)/XXVB(3) 1439 : CCC HIJ(4) = UNS3*XXVB(3) 1440 : CCC HIJ(5) = UNS3*XXVB(1) 1441 : CCC HIJ(6) = UNS3*XXVB(2) 1442 : C 1443 : C CALCUL DES COEFS A METTRE DANS KIJ POUR COMPOSER KSTAB 1444 : C 1445 : C 1446 : C ICI IL FAUT PRENDRE LA MOYENNE DES CONTRAINTES SUR LES POINTS DE GAUSS 1447 : C CECI EST DEJA CALCULE PLUS HAUT DANS LA VARIABLE SIGMOY(6) 1448 : C 1449 : C 1450 : C POUR LA PROJECTION 28 ON OBTIENT 1451 : C 1452 : CCC XK11_1 = SIGMOY(1)*HIJ(1) 1453 : CCC XK11_2 = SIGMOY(1)*HIJ(1)*UNS3 1454 : CCC XK11_3 = SIGMOY(6)*HIJ(6) 1455 : CCC XK11_4 = SIGMOY(6)*HIJ(6) 1456 : C 1457 : CCC XK22_1 = SIGMOY(2)*HIJ(2) 1458 : CCC XK22_2 = SIGMOY(2)*HIJ(2)*UNS3 1459 : CCC XK22_3 = SIGMOY(5)*HIJ(5) 1460 : CCC XK22_4 = SIGMOY(5)*HIJ(5) 1461 : C 1462 : CCC XK33_1 = ZERO 1463 : CCC XK33_2 = (SIGMOY(1)*HIJ(1)+SIGMOY(2)*HIJ(2))*UNS3 1464 : C 1465 : CCC XK13_1 = ZERO 1466 : CCC XK13_2 = ZERO 1467 : CCC XK23_1 = ZERO 1468 : CCC XK23_2 = ZERO 1469 : CCC XK31_1 = ZERO 1470 : CCC XK31_2 = ZERO 1471 : CCC XK32_1 = ZERO 1472 : CCC XK32_2 = ZERO 1473 : C 1474 : C 1474 : CCC CALL ZDANUL(XK11,64) 1474 : CCC CALL ZDANUL(XK22,64) 1474 : CCC CALL ZDANUL(XK12,64) 1474 : CCC CALL ZDANUL(XK21,64) 1474 : CCC CALL ZDANUL(XK13,64) 1474 : CCC CALL ZDANUL(XK23,64) 1474 : CCC CALL ZDANUL(XK31,64) 1474 : CCC CALL ZDANUL(XK32,64) 1474 : CCC CALL ZDANUL(XK33,64) 1475 : C 1476 : CCC DO J=1,8 1477 : CCC DO I=1,8 1478 : C 1479 : C IL FAUT CALCULER K11 K22 K33 K13 K23 K31 K32 MATRICES 8*8 1480 : C 1481 : CCC XK11(I,J)=XK11_1*GS(I,3)*GS(J,3)+XK11_2*GS(I,4)*GS(J,4) 1482 : C & +XK11_3*GS(I,3)*GS(J,1)+XK11_4*GS(I,1)*GS(J,3) 1483 : CCC XK22(I,J)=XK22_1*GS(I,3)*GS(J,3)+XK22_2*GS(I,4)*GS(J,4) 1484 : C & +XK22_3*GS(I,3)*GS(J,2)+XK22_4*GS(I,2)*GS(J,3) 1485 : CCC XK33(I,J)=XK33_1*GS(I,3)*GS(J,3)+XK33_2*GS(I,4)*GS(J,4) 1486 : C 1487 : CCC XK13(I,J)=XK13_1*GS(I,3)*GS(J,1)+XK13_2*GS(I,1)*GS(J,3) 1488 : CCC XK23(I,J)=XK23_1*GS(I,3)*GS(J,2)+XK23_2*GS(I,2)*GS(J,3) 1489 : CCC XK31(I,J)=XK31_1*GS(I,3)*GS(J,1)+XK31_2*GS(I,1)*GS(J,3) 1490 : CCC XK32(I,J)=XK32_1*GS(I,3)*GS(J,2)+XK32_2*GS(I,2)*GS(J,3) 1491 : C 1492 : CCC END DO 1493 : CCC END DO 1494 : C 1495 : C ASSEMBLAGE DE KSTAB 1496 : C 1496 : CCC CALL ASSE_KSTAB(XKSTAB,XK11,XK22,XK33,XK12,XK21,XK13,XK23, 1497 : CCC ! XK31,XK32) 1498 : C 1499 : C 1500 : C REMISE EN ORDRE DE KSTAB 1501 : C 1502 : CCC DO J=1,8 1503 : CCC DO I=1,24 1504 : CCC TMPKE(I,(J-1)*3+1)=XKSTAB(I,J) 1505 : CCC TMPKE(I,(J-1)*3+2)=XKSTAB(I,J+8) 1506 : CCC TMPKE(I,(J-1)*3+3)=XKSTAB(I,J+16) 1507 : CCC END DO 1508 : CCC END DO 1508 : CCC CALL ZDANUL(XKSTAB,576) 1509 : CCC DO I=1,8 1510 : CCC DO J=1,24 1511 : CCC XKSTAB((I-1)*3+1,J)=TMPKE(I,J) 1512 : CCC XKSTAB((I-1)*3+2,J)=TMPKE(I+8,J) 1513 : CCC XKSTAB((I-1)*3+3,J)=TMPKE(I+16,J) 1514 : CCC END DO 1515 : CCC END DO 1516 : C 1517 : C IL FAUT REPASSER DANS LE REPERE GLOBAL AVEC RR2^T . KSTAB . RR2 1518 : C EN FAIT C'EST RR2. KSTAB .RR2^T CAR RR2 RANGEE PAR LIGNES 1519 : C 1519 : CCC CALL ASSE_KSTAB_GL(XKSTAB,RR2) 1520 : C 1521 : C AJOUT DE KSTAB A KE 1522 : C 1523 : CCC DO J=1,24 1524 : CCC DO I=1,24 1525 : CCC RE(I,J)=RE(I,J)+XKSTAB(I,J) 1526 : CCC END DO 1527 : CCC END DO 1528 : C 1529 : C 1530 : C 1531 : ENDIF 1532 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1533 : IF(ICLE.EQ.3)THEN 1534 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1535 : C C 1536 : C CALCUL DE LA MATRICE DE MASSE C 1537 : C ENTREE: PROPEL(1) = MASSE VOLUMIQUE C 1538 : C XE(1 A 14)= COORDONNEES DES NOEUDS C 1539 : C SORTIE: RE(24,24) = MATRICE MASSE C 1540 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1541 : XRO = PROPEL(1) 1541 : CALL ZDANUL(RE,576) 1542 : C 1543 : C CALCUL DU PLAN MOYEN: 1544 : C 1545 : DO J=1,3 1546 : DO I=1,4 1547 : XCOQ(J,I) = 0.5*XE((I-1)*3+J)+0.5*XE((I-1+4)*3+J) 1548 : END DO 1549 : END DO 1550 : C CALCUL DE LA SURFACE MOYENNE: 1551 : DO J=1,3 1552 : X12(J) = (XCOQ(J,1)+XCOQ(J,2))*0.5 1553 : X34(J) = (XCOQ(J,3)+XCOQ(J,4))*0.5 1554 : X23(J) = (XCOQ(J,2)+XCOQ(J,3))*0.5 1555 : X41(J) = (XCOQ(J,4)+XCOQ(J,1))*0.5 1556 : END DO 1557 : DO J=1,3 1558 : V12_34(J) = X34(J)-X12(J) 1559 : V41_23(J) = X23(J)-X41(J) 1560 : END DO 1561 : C 1562 : C CALCUL DE V41_23 ^ V12_34: 1563 : C 1564 : VINF(1) = V41_23(2) * V12_34(3) - V41_23(3) * V12_34(2) 1565 : VINF(2) = V41_23(3) * V12_34(1) - V41_23(1) * V12_34(3) 1566 : VINF(3) = V41_23(1) * V12_34(2) - V41_23(2) * V12_34(1) 1567 : C 1568 : C CALCUL DU VOLUME DE L'ELEMENT: 1569 : C 1569 : CALL CUBKSI(5,XXG5,BKSIP) 1569 : CALL CUCALB(BKSIP(1,1,3),XE,B,AJAC) 1570 : XVOLUME = AJAC*8. 1570 : CALL ZDANUL(XMASSE,64) 1571 : XXXMASSE = 0. 1572 : DO I=1,8 1573 : XMASSE(I,I) = (8./27.)*AJAC*XRO 1574 : XXXMASSE = XXXMASSE+XMASSE(I,I) 1575 : ENDDO 1576 : DO I=1,8 1577 : DO J=1,8 1578 : IF(I.NE.J)THEN 1579 : XMASSE(I,J) = (1./64.)*(2.+VB(I,1)*VB(J,1)*2./3.) 1580 : ! *(2.+VB(I,2)*VB(J,2)*2./3.) 1581 : ! *(2.+VB(I,3)*VB(J,3)*2./3.) 1582 : ! *AJAC*XRO 1583 : XXXMASSE = XXXMASSE+XMASSE(I,I) 1584 : ENDIF 1585 : ENDDO 1586 : ENDDO 1586 : CALL ZDANUL(TMPKE,576) 1587 : DO K=1,3 1588 : DO I=1,8 1589 : DO J=1,8 1590 : TMPKE(I+(K-1)*8,J+(K-1)*8) = XMASSE(I,J) 1591 : ENDDO 1592 : ENDDO 1593 : ENDDO 1594 : DO J=1,8 1595 : DO I=1,24 1596 : TMPKE2(I,(J-1)*3+1)=TMPKE(I,J) 1597 : TMPKE2(I,(J-1)*3+2)=TMPKE(I,J+8) 1598 : TMPKE2(I,(J-1)*3+3)=TMPKE(I,J+16) 1599 : END DO 1600 : END DO 1600 : CALL ZDANUL(RE,576) 1601 : DO I=1,8 1602 : DO J=1,24 1603 : RE((I-1)*3+1,J)=TMPKE2(I,J) 1604 : RE((I-1)*3+2,J)=TMPKE2(I+8,J) 1605 : RE((I-1)*3+3,J)=TMPKE2(I+16,J) 1606 : END DO 1607 : END DO 1608 : ENDIF 1609 : IF(ICLE.EQ.10)THEN 1610 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1611 : C C 1612 : C CALCUL DE LA MATRICE KP C 1613 : C ENTREE: PRESSION DANS PROPEL(1) C 1614 : C KTAFL DANS PROPEL(2) C 1615 : C RE=RE+KP C 1616 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 1617 : PRESSI = PROPEL(1) 1618 : IFACE = nint(PROPEL(2)) 1619 : ITEST = 0 1620 : IF (IFACE.EQ.1) THEN 1621 : C 1622 : C DEFINITION DE LA NORMALE N0, FACE INFERIEURE 1623 : C 1624 : DO I=1,4 1625 : DO J=1,3 1626 : XCOQ(J,I) = XE((I-1)*3+J) 1627 : END DO 1628 : END DO 1629 : DO J=1,3 1630 : XCOQ14(J) = (XCOQ(J,1)+XCOQ(J,4))/2. 1631 : XCOQ23(J) = (XCOQ(J,2)+XCOQ(J,3))/2. 1632 : XCOQ34(J) = (XCOQ(J,3)+XCOQ(J,4))/2. 1633 : XCOQ12(J) = (XCOQ(J,1)+XCOQ(J,2))/2. 1634 : ENDDO 1635 : DO J=1,3 1636 : XCOQ12_34(J) = XCOQ34(J)-XCOQ12(J) 1637 : XCOQ14_23(J) = XCOQ23(J)-XCOQ14(J) 1638 : ENDDO 1639 : C 1640 : XNORMALE(1) = XCOQ14_23(2)*XCOQ12_34(3) 1641 : ! -XCOQ14_23(3)*XCOQ12_34(2) 1642 : XNORMALE(2) = XCOQ14_23(3)*XCOQ12_34(1) 1643 : ! -XCOQ14_23(1)*XCOQ12_34(3) 1644 : XNORMALE(3) = XCOQ14_23(1)*XCOQ12_34(2) 1645 : ! -XCOQ14_23(2)*XCOQ12_34(1) 1646 : XNORMALE(1) = -XNORMALE(1) 1647 : XNORMALE(2) = -XNORMALE(2) 1648 : XNORMALE(3) = -XNORMALE(3) 1649 : C 1650 : XXG5(1) = -1. 1651 : XXG5(2) = 0. 1652 : XXG5(3) = 0. 1653 : XXG5(4) = 0. 1654 : XXG5(5) = 0. 1655 : ITEST = 1 1656 : ENDIF 1657 : IF (IFACE.EQ.2) THEN 1658 : C 1659 : C DEFINITION DE LA NORMALE N0, FACE SUPERIEURE 1660 : C 1661 : DO I=1,4 1662 : DO J=1,3 1663 : XCOQ(J,I) = XE(12+(I-1)*3+J) 1664 : END DO 1665 : END DO 1666 : DO J=1,3 1667 : XCOQ14(J)=(XCOQ(J,1)+XCOQ(J,4))/2. 1668 : XCOQ23(J)=(XCOQ(J,2)+XCOQ(J,3))/2. 1669 : XCOQ34(J)=(XCOQ(J,3)+XCOQ(J,4))/2. 1670 : XCOQ12(J)=(XCOQ(J,1)+XCOQ(J,2))/2. 1671 : ENDDO 1672 : DO J=1,3 1673 : XCOQ12_34(J)=XCOQ34(J)-XCOQ12(J) 1674 : XCOQ14_23(J)=XCOQ23(J)-XCOQ14(J) 1675 : ENDDO 1676 : C 1677 : XNORMALE(1)=XCOQ14_23(2)*XCOQ12_34(3) 1678 : ! -XCOQ14_23(3)*XCOQ12_34(2) 1679 : XNORMALE(2)=XCOQ14_23(3)*XCOQ12_34(1) 1680 : ! -XCOQ14_23(1)*XCOQ12_34(3) 1681 : XNORMALE(3)=XCOQ14_23(1)*XCOQ12_34(2) 1682 : ! -XCOQ14_23(2)*XCOQ12_34(1) 1683 : C 1684 : XXG5(1)=+1. 1685 : XXG5(2)=0. 1686 : XXG5(3)=0. 1687 : XXG5(4)=0. 1688 : XXG5(5)=0. 1689 : ITEST=1 1690 : ENDIF 1691 : IF(ITEST.EQ.0)THEN 1692 : WRITE(6,*)'ARRET DANS SHB8, CALCUL DE KP, FACE NON DEFINIE' 1693 : STOP 1694 : ENDIF 1695 : C 1696 : C CALCUL DE B SUR LA FACE INF OU SUP, SUIVANT XXG5 1697 : C 1698 : C DANS UE: ON MET L'EQUIVALENT DE BKSIP, CALCULER 1699 : C POUR LES 4 NOEUDS DE LA FACE IFACE 1700 : C 1700 : CALL CUKSKP(IFACE,UE) 1701 : C 1702 : XNORM = SQRT(XNORMALE(1)*XNORMALE(1)+XNORMALE(2)*XNORMALE(2) 1703 : ! +XNORMALE(3)*XNORMALE(3)) 1704 : DO J=1,3 1705 : XNORMALE(J) = XNORMALE(J)/XNORM 1706 : ENDDO 1706 : CALL CUCBKP(IFACE,UE,XE,BKP,AJAC) 1706 : CALL ASSKP2(IFACE,BKP,XNORMALE,XKP) 1707 : C 1708 : C REMISE EN ORDRE DE KP 1709 : C 1710 : DO I=1,24 1711 : DO J=1,8 1712 : TMPKE(I,(J-1)*3+1)=XKP(I,J) 1713 : TMPKE(I,(J-1)*3+2)=XKP(I,J+8) 1714 : TMPKE(I,(J-1)*3+3)=XKP(I,J+16) 1715 : END DO 1716 : END DO 1716 : CALL ZDANUL(XKP,576) 1717 : DO I=1,8 1718 : DO J=1,24 1719 : XKP((I-1)*3+1,J)=TMPKE(I,J) 1720 : XKP((I-1)*3+2,J)=TMPKE(I+8,J) 1721 : XKP((I-1)*3+3,J)=TMPKE(I+16,J) 1722 : END DO 1723 : END DO 1724 : DO I=1,24 1725 : DO J=1,24 1726 : RE(I,J)=RE(I,J)+XKP(I,J)*XNORM*PRESSI/8. 1727 : ENDDO 1728 : ENDDO 1729 : ENDIF 1730 : RETURN 1731 : END 1732 : 1733 : 1734 : 1735 : 1736 :
© Cast3M 2003 - All rights reserved.
Disclaimer