C ELLA23 SOURCE CB215821 25/04/24 21:15:09 12248 SUBROUTINE ELLA23(CARACT,COOR,GAMA,ZXX,XCOR,VALDE1,VALDE2, & VALDE3,VALDE4,ZS,NP,NBELEM,XXPI) C IMPLICIT INTEGER(I-N) COMPLEX*16 ZXX,ZD,ZE,ZI,ZIMOIN,ZS,ZU,ZU1,ZU2,ZCE,ZCP,ZCMU,ZCR, & ZALFA2,ZALFAX,ZALFAR,ZALFA4,ZALFAF,ZZ3,ZALFAA, & ZA56,ZB56,ZA89,ZB89,ZA1112,ZB1112, & ZUX1,ZUY1,ZUZ1,ZUP1,ZX1,ZY1,ZZ1,ZP1, & ZUX2,ZUY2,ZUZ2,ZUP2,ZX2,ZY2,ZZ2,ZP2, & Z1,Z2, & Z10X, Z11X, Z30X, Z31X, Z32X, Z33X, Z34X, Z40X, Z41X, & Z10XL,Z11XL,Z30XL,Z31XL,Z32XL,Z33XL,Z34XL,Z40XL,Z41XL C REAL*8 CARACT,COOR,GAMA,XCOR,VALDE1,VALDE2,VALDE3, & VALDE4,XXPI,R1,R2,XXA,XYA,XZA,XXB,XYB,XZB,XL, & XI1,XI2,XI3,GX,GY,GZ,GG,DELTA,DET,XJ1,XJ2,XJ3, & XK1,XK2,XK3,XX1,XY1,XZ1,XX2,XY2,XZ2, & XXD1,XXD2,XXL1,XXL2,EPSILO REAL*8 CE,CNU,CRO,CRIN,CREX,CKC,CAM,CETA,CROF,CSON, & CSE,CSEF,CC,CI,CIP,COEF, & PRZX1,PIZX1,PRZY1,PIZY1,PRZZ1,PIZZ1,PRZP1,PIZP1, & PRZX2,PIZX2,PRZY2,PIZY2,PRZZ2,PIZZ2,PRZP2,PIZP2, & PPX1,PPY1,PPZ1,PPX2,PPY2,PPZ2,PPP1,PPP2 C -INC CCREEL INTEGER NP,NBELEM,N1,N2,I,IN,IOUI,J,L C C OPERATEUR ELFE LAPLACE ACOU C C CALCUL POUR LA POUTRE Num INP LES VALEURS DE Ua ET DE Ub EN LOCAL C C PARAMETRES : C C CARACT : TABLEAU DES CARACTERISTIQUE DES POUTRES (10,NP) C COOR : TABLEAU DES COORDONNEES DES NOEUDS (3,2*NP) C GAMA : TABLEAU DES VECTEUR DEFINISSANT LE PLAN LOCALE OXY (3,2*NP) C ZXX : VECTEUR (28*NP) DONNANT LE VECTEUR SOLUTION DU TREILLIS C DANS LE REPERE GLOBAL C XCOR : TABLEAU DONNANT LES COORD. DES NOEUDS DE L'OBJET MAILLAGE C DANS LE REPERE GLOBAL C ZS : VALEUR DE S = S0 + I W C NP : NOMBRE DE POUTRES DU TREILLIS C NBELEM : NOMBRE D'ELEMENTS DU SOUS-MAILLAGE C XXPI : VALEUR PRECISE DE PI C C C SORTIES : C C VALDE1 : MODULE DE UX, UY, UZ POUR LES 2 POINTS DE CHAQUE ELEMENT C DU SOUS-MAILLAGE C VALDE2 : PHASE DE UX, UY, UZ POUR LES 2 POINTS DE CHAQUE ELEMENT C DU SOUS-MAILLAGE C VALDE3 : MODULE DE P POUR LES 2 POINTS DE CHAQUE ELEMENT C DU SOUS-MAILLAGE C VALDE4 : PHASE DE P POUR LES 2 POINTS DE CHAQUE ELEMENT C DU SOUS-MAILLAGE C C C AUTEURS : SAINT-DIZIER ET GORCY C DATE : 24 JANVIER 1991 C C DIMENSION CARACT(10,*),COOR(3,*),GAMA(3,*) DIMENSION ZXX(*),XCOR(2,3,NBELEM) DIMENSION VALDE1(2,NBELEM,3),VALDE2(2,NBELEM,3) DIMENSION VALDE3(2,NBELEM,1),VALDE4(2,NBELEM,1) C DIMENSION ZU(28),ZU1(28),ZU2(28) DIMENSION R1(3,3),R2(3,3),ZD(14),ZE(14,14) C ZIMOIN=CMPLX(0.D0,-1.D0) * EPSILO=1.0D-38 EPSILO=xpetit C DO 100 INP = 1 , NP C DO 101 II = 1 , 28 ZU(II) = ZXX(28*(INP-1)+II) 101 CONTINUE C C -------------- ACQUISITION DES CARACTERISTIQUES DE LA POUTRE C --------------------------------------------- C CE = CARACT( 1,INP) CNU = CARACT( 2,INP) CRO = CARACT( 3,INP) CRIN = CARACT( 4,INP) CREX = CARACT( 5,INP) CKC = CARACT( 6,INP) CAM = CARACT( 7,INP) CETA = CARACT( 8,INP) CROF = CARACT( 9,INP) CSON = CARACT(10,INP) C ZI = CMPLX(0.D0 , 1.D0) C CSE = XXPI * ((CREX*CREX)-(CRIN*CRIN)) CSEF = XXPI * (CRIN*CRIN) CC = XXPI * ((CREX*CREX*CREX*CREX)-(CRIN*CRIN*CRIN*CRIN)) / 2.D0 CI = XXPI * ((CREX*CREX*CREX*CREX)-(CRIN*CRIN*CRIN*CRIN)) / 4.D0 CIP = 2.D0 * CI C COEF = 1.D0 +((2.*CROF*CSON*CSON*CRIN) / (CE*(CREX-CRIN))) CSON = CSON / SQRT(COEF) C ZCE = CE*(CMPLX(1.D0,0.D0)+CETA*ZI) ZCP = SQRT(ZCE/CRO) ZCMU = ZCE/(CMPLX(2.D0,0.D0)* > (CMPLX(1.D0,0.D0)+CNU)) ZCR = SQRT(ZCMU/CRO) C ZALFA2= ZS*ZS*CRO*CSE + CAM*ZS ZALFA2= ZALFA2 / (ZCE*CSE) ZALFAX= SQRT(ZALFA2) C ZALFA2= ZS*ZS*CRO*CIP + CAM*ZS ZALFA2= ZALFA2 / (ZCMU*CC) ZALFAR= SQRT(ZALFA2) C ZALFA4= CRO*CRO*CI*ZS*ZS*ZS*ZS/ZCMU/CKC * + CAM*CRO*CI*ZS*ZS*ZS/ZCMU/CKC/CSE * + ((CRO*CSE)+(CROF*CSEF))*ZS*ZS * + CAM*ZS ZALFA4= ZALFA4/(ZCE*CI) ZALFA2= SQRT(ZALFA4) ZALFAF= SQRT(ZALFA2) C ZZ3 = CRO*CI*(CMPLX(1.D0,0.D0)+ZCE/(ZCMU*CKC))*ZS*ZS * + CAM*ZCE*CI*ZS/(ZCMU*CKC*CSE) ZZ3 = ZZ3 / (ZCE*CI*ZALFA2) C ZALFAA = ZS / CSON C N1 = 2*INP-1 N2 = 2*INP XXA = COOR (1,N1) XYA = COOR (2,N1) XZA = COOR (3,N1) XXB = COOR (1,N2) XYB = COOR (2,N2) XZB = COOR (3,N2) C XL = SQRT((XXB-XXA)**2 + (XYB-XYA)**2 + (XZB-XZA)**2) C C ------------------------------ VECTEUR UNITAIRE OX REPERE LOCALE C --------------------------------- XI1 = (XXB-XXA)/XL XI2 = (XYB-XYA)/XL XI3 = (XZB-XZA)/XL C C ------------------------- VECTEUR UNITAIRE OY REPERE LOCALE C --------------------------------- GX = GAMA(1,INP) GY = GAMA(2,INP) GZ = GAMA(3,INP) GG = SQRT(GX*GX + GY*GY + GZ*GZ) GX = GX/GG GY = GY/GG GZ = GZ/GG C DELTA = SQRT (1.D0 - (XI1*GX + XI2*GY + XI3*GZ)**2) C DET=-(GX*XI3-GZ*XI1)**2-(GY*XI1-GX*XI2)**2-(GY*XI3-GZ*XI2)**2 C IF (ABS(DET).LT.1D-12) THEN XJ1 = -XI2 XJ2 = 0.D0 XJ3 = 0.D0 ELSE XJ1 = (XI2*(GY*XI1-GX*XI2)-XI3*(GX*XI3-GZ*XI1))*DELTA/DET XJ2 = (XI3*(GZ*XI2-GY*XI3)-XI1*(GY*XI1-GX*XI2))*DELTA/DET XJ3 = (XI1*(GX*XI3-GZ*XI1)-XI2*(GZ*XI2-GY*XI3))*DELTA/DET END IF C C ---------------------------- VECTEUR UNITAIRE OZ REPERE LOCALE C --------------------------------- XK1 = XI2*XJ3 - XI3*XJ2 XK2 = XI3*XJ1 - XI1*XJ3 XK3 = XI1*XJ2 - XI2*XJ1 C R1(1,1) = XI1 R1(1,2) = XJ1 R1(1,3) = XK1 R1(2,1) = XI2 R1(2,2) = XJ2 R1(2,3) = XK2 R1(3,1) = XI3 R1(3,2) = XJ3 R1(3,3) = XK3 C R2(1,1) = XJ2*XK3 - XJ3*XK2 R2(1,2) = XJ3*XK1 - XJ1*XK3 R2(1,3) = XJ1*XK2 - XJ2*XK1 R2(2,1) = XI3*XK2 - XI2*XK3 R2(2,2) = XI1*XK3 - XI3*XK1 R2(2,3) = XI2*XK1 - XI1*XK2 R2(3,1) = XI2*XJ3 - XI3*XJ2 R2(3,2) = XI3*XJ1 - XI1*XJ3 R2(3,3) = XI1*XJ2 - XI2*XJ1 C C ----------------------------- MATRICE DIMENSIONNEMENT ZD C -------------------------- ZD( 1) = XL**2/ZCP ZD( 2) = XL**2/ZCP ZD( 3) = XL**2/ZCP ZD( 4) = XL/ZCR ZD( 5) = XL**2/ZCP*ZALFAF ZD( 6) = XL**2/ZCP*ZALFAF ZD( 7) = XL**2/ZCP*ZALFAX ZD( 8) = XL**2/ZCP*(ZALFAF**3) ZD( 9) = XL**2/ZCP*(ZALFAF**3) ZD(10) = XL/ZCR*ZALFAR ZD(11) = XL**2/ZCP*(ZALFAF**2) ZD(12) = XL**2/ZCP*(ZALFAF**2) ZD(13) = CMPLX(XL*CROF*CSON,0.D0) ZD(14) = XL*CROF*CSON*ZALFAA C C ------------------------------- MATRICE DES EFFORTS ZE C ---------------------- DO 20 I = 1 , 14 DO 21 J = 1 , 14 ZE(I,J) = CMPLX(0.D0,0.D0) 21 CONTINUE 20 CONTINUE C ZA56 = CMPLX(1.D0,0.D0) ZB56 = CMPLX(1.D0,0.D0) / (ZCMU*CKC*CSE) C ZA89 = (CRO*CSE*ZS*ZS+CAM*ZS)/(ZCMU*ZCMU*CKC*CKC*CSE*CSE) * - CMPLX(1.D0,0.D0)/(ZCE*CI) ZB89 = CRO*ZS*ZS*(CMPLX(1.D0,0.D0)/(ZCMU*CKC)+ > CMPLX(1.D0,0.D0)/ZCE) * + CAM*ZS/(ZCMU*CKC*CSE) C ZA1112 = CMPLX(1.D0,0.D0)/(ZCE*CI) ZB1112 = (CRO*CSE*ZS*ZS+CAM*ZS)/(ZCMU*CKC*CSE) C ZE( 1, 1) = CMPLX(1.D0,0.D0) ZE( 2, 2) = CMPLX(1.D0,0.D0) ZE( 3, 3) = CMPLX(1.D0,0.D0) ZE( 4, 4) = CMPLX(1.D0,0.D0) ZE( 5, 5) = -ZA56 ZE( 5, 9) = ZB56 ZE( 6, 6) = ZA56 ZE( 6, 8) = ZB56 ZE( 7, 7) = CMPLX(1.D0,0.D0)/(ZCE*CSE) ZE( 8, 6) = ZB89 ZE( 8, 8) = ZA89 ZE( 9, 5) = -ZB89 ZE( 9, 9) = ZA89 ZE(10,10) = CMPLX(1.D0,0.D0)/(ZCMU*CC) ZE(11, 3) = ZB1112 ZE(11,11) = -ZA1112 ZE(12, 2) = ZB1112 ZE(12,12) = ZA1112 ZE(13,13) = CMPLX(1.D0,0.D0) ZE(14,14) = -ZS/CSEF C C -- CALCUL DE ZU1 = ZR * ZU C DO 27 I = 1 , 28 ZU1(I) = CMPLX(0.D0,0.D0) ZU2(I) = CMPLX(0.D0,0.D0) 27 CONTINUE C DO 26 I = 1 , 4 DO 25 J = 1 , 3 DO 24 K = 1 , 3 L = (I-1)*3 ZU1(L+J)=ZU1(L+J)+R2(J,K)*ZU(L+K) ZU1(L+14+J)=ZU1(L+14+J)+R2(J,K)*ZU(L+14+K) 24 CONTINUE 25 CONTINUE 26 CONTINUE C ZU1(13) = ZU(13) ZU1(14) = ZU(14) ZU1(27) = ZU(27) ZU1(28) = ZU(28) C C -- CALCUL DE ZU2 = ZE * ZU1 C DO 30 I = 1,12 DO 31 J = 1,12 ZU2(I) = ZU2(I) + ZE(I,J)*ZU1(J) 31 CONTINUE 30 CONTINUE C ZU2(13) = ZE(13,13)*ZU1(13) ZU2(14) = -ZE(14,14)*ZU1(14) C DO 32 I = 1,12 DO 33 J = 1,6 ZU2(I+14) = ZU2(I+14) + ZE(I,J )*ZU1(J+14) ZU2(I+14) = ZU2(I+14) - ZE(I,J+6)*ZU1(J+20) 33 CONTINUE 32 CONTINUE C ZU2(27) = ZE(13,13)*ZU1(27) ZU2(28) = ZE(14,14)*ZU1(28) C C -- CALCUL DE ZU2 = ZD * ZU2 C DO 34 I = 1,14 ZU2(I ) = ZU2(I )/ZD(I) ZU2(I+14) = ZU2(I+14)/ZD(I) 34 CONTINUE C C --------------- NBELEM : NOMBRE D'ELEMENTS DANS LE SOUS MAILLAGE C DO 110 IN = 1 , NBELEM C XX1 = XCOR ( 1,1,IN) XY1 = XCOR ( 1,2,IN) XZ1 = XCOR ( 1,3,IN) C XX2 = XCOR ( 2,1,IN) XY2 = XCOR ( 2,2,IN) XZ2 = XCOR ( 2,3,IN) C CALL ELLA07(XXA,XYA,XZA,XXB,XYB,XZB,XX1,XY1,XZ1,XX2,XY2,XZ2, * IOUI) C IF (IOUI.EQ.1) THEN C XXD1 = SQRT((XXA-XX1)**2+(XYA-XY1)**2+(XZA-XZ1)**2) XXD2 = SQRT((XXA-XX2)**2+(XYA-XY2)**2+(XZA-XZ2)**2) XXL1= XL -XXD1 XXL2= XL -XXD2 C C XXDI : DISTANCE DE A AU PT M DE CALCUL DU DEPLACEMENT C CALL ELLA31(XXD1,ZALFAX,ZALFAR,ZALFAF,ZALFAA,Z10X, * Z11X,Z1,Z2,Z30X,Z31X,Z32X,Z33X,Z34X,ZZ3,Z40X,Z41X) C CALL ELLA31(XXL1,ZALFAX,ZALFAR,ZALFAF,ZALFAA,Z10XL, * Z11XL,Z1,Z2,Z30XL,Z31XL,Z32XL,Z33XL,Z34XL,ZZ3,Z40XL,Z41XL) C ZUX1 = ZU2(21)*Z10XL-ZU2(15)*Z11XL-ZU2(7)*Z10X-ZU2(1)*Z11X C ZUY1 =-(ZZ3*Z31X -Z33X )*ZU2(2) - (ZZ3*Z30X -Z32X )*ZU2(6) * + Z31X *ZU2(12) + Z30X *ZU2(8) * -(ZZ3*Z31XL-Z33XL)*ZU2(16)+(ZZ3*Z30XL-Z32XL)*ZU2(20) * + Z31XL*ZU2(26) - Z30XL*ZU2(22) C ZUZ1 =-(ZZ3*Z31X -Z33X )*ZU2(3) - (ZZ3*Z30X -Z32X )*ZU2(5) * + Z31X *ZU2(11) + Z30X *ZU2(9) * -(ZZ3*Z31XL-Z33XL)*ZU2(17)+(ZZ3*Z30XL-Z32XL)*ZU2(19) * + Z31XL*ZU2(25) - Z30XL*ZU2(23) C ZUP1 = ZU2(28)*Z40XL-ZU2(27)*Z41XL-ZU2(14)*Z40X-ZU2(13)*Z41X C ZUX1 = ZUX1 * ZD(1) ZUY1 = ZUY1 * ZD(2) ZUZ1 = ZUZ1 * ZD(3) ZUP1 = ZUP1 * ZD(13) C ZX1 = R1(1,1)*ZUX1 + R1(1,2)*ZUY1 + R1(1,3)*ZUZ1 ZY1 = R1(2,1)*ZUX1 + R1(2,2)*ZUY1 + R1(2,3)*ZUZ1 ZZ1 = R1(3,1)*ZUX1 + R1(3,2)*ZUY1 + R1(3,3)*ZUZ1 ZP1 = ZUP1 C CALL ELLA31(XXD2,ZALFAX,ZALFAR,ZALFAF,ZALFAA,Z10X, * Z11X,Z1,Z2,Z30X,Z31X,Z32X,Z33X,Z34X,ZZ3,Z40X,Z41X) C CALL ELLA31(XXL2,ZALFAX,ZALFAR,ZALFAF,ZALFAA,Z10XL, * Z11XL,Z1,Z2,Z30XL,Z31XL,Z32XL,Z33XL,Z34XL,ZZ3,Z40XL,Z41XL) C C ZUX2 = ZU2(21)*Z10XL-ZU2(15)*Z11XL-ZU2(7)*Z10X-ZU2(1)*Z11X C ZUY2 =-(ZZ3*Z31X -Z33X )*ZU2(2) - (ZZ3*Z30X -Z32X )*ZU2(6) * + Z31X *ZU2(12) + Z30X *ZU2(8) * -(ZZ3*Z31XL-Z33XL)*ZU2(16)+(ZZ3*Z30XL-Z32XL)*ZU2(20) * + Z31XL*ZU2(26) - Z30XL*ZU2(22) C ZUZ2 =-(ZZ3*Z31X -Z33X )*ZU2(3) - (ZZ3*Z30X -Z32X )*ZU2(5) * + Z31X *ZU2(11) + Z30X *ZU2(9) * -(ZZ3*Z31XL-Z33XL)*ZU2(17)+(ZZ3*Z30XL-Z32XL)*ZU2(19) * + Z31XL*ZU2(25) - Z30XL*ZU2(23) C ZUP2 = ZU2(28)*Z40XL-ZU2(27)*Z41XL-ZU2(14)*Z40X-ZU2(13)*Z41X C ZUX2 = ZUX2 * ZD(1) ZUY2 = ZUY2 * ZD(2) ZUZ2 = ZUZ2 * ZD(3) ZUP2 = ZUP2 * ZD(13) C ZX2 = R1(1,1)*ZUX2 + R1(1,2)*ZUY2 + R1(1,3)*ZUZ2 ZY2 = R1(2,1)*ZUX2 + R1(2,2)*ZUY2 + R1(2,3)*ZUZ2 ZZ2 = R1(3,1)*ZUX2 + R1(3,2)*ZUY2 + R1(3,3)*ZUZ2 ZP2 = ZUP2 C PRZX1 = ZX1 PIZX1 = ZX1*ZIMOIN PRZY1 = ZY1 PIZY1 = ZY1*ZIMOIN PRZZ1 = ZZ1 PIZZ1 = ZZ1*ZIMOIN PRZP1 = ZP1 PIZP1 = ZP1*ZIMOIN PRZX2 = ZX2 PIZX2 = ZX2*ZIMOIN PRZY2 = ZY2 PIZY2 = ZY2*ZIMOIN PRZZ2 = ZZ2 PIZZ2 = ZZ2*ZIMOIN PRZP2 = ZP2 PIZP2 = ZP2*ZIMOIN C IF (ABS(PRZX1).LT.EPSILO.AND.ABS(PIZX1).LT.EPSILO) THEN PPX1 = 0.D0 ELSE PPX1 = ATAN2(PIZX1,PRZX1)*180.D0/XXPI END IF C IF (ABS(PRZY1).LT.EPSILO.AND.ABS(PIZY1).LT.EPSILO) THEN PPY1 = 0.D0 ELSE PPY1 = ATAN2(PIZY1,PRZY1)*180.D0/XXPI END IF C IF (ABS(PRZZ1).LT.EPSILO.AND.ABS(PIZZ1).LT.EPSILO) THEN PPZ1 = 0.D0 ELSE PPZ1 = ATAN2(PIZZ1,PRZZ1)*180.D0/XXPI END IF C IF (ABS(PRZX2).LT.EPSILO.AND.ABS(PIZX2).LT.EPSILO) THEN PPX2 = 0.D0 ELSE PPX2 = ATAN2(PIZX2,PRZX2)*180.D0/XXPI END IF C IF (ABS(PRZY2).LT.EPSILO.AND.ABS(PIZY2).LT.EPSILO) THEN PPY2 = 0.D0 ELSE PPY2 = ATAN2(PIZY2,PRZY2)*180.D0/XXPI END IF C IF (ABS(PRZZ2).LT.EPSILO.AND.ABS(PIZZ2).LT.EPSILO) THEN PPZ2 = 0.D0 ELSE PPZ2 = ATAN2(PIZZ2,PRZZ2)*180.D0/XXPI END IF C IF (ABS(PRZP1).LT.EPSILO.AND.ABS(PIZP1).LT.EPSILO) THEN PPP1 = 0.D0 ELSE PPP1 = ATAN2(PIZP1,PRZP1)*180.D0/XXPI END IF C IF (ABS(PRZP2).LT.EPSILO.AND.ABS(PIZP2).LT.EPSILO) THEN PPP2 = 0.D0 ELSE PPP2 = ATAN2(PIZP2,PRZP2)*180.D0/XXPI END IF C VALDE1(1,IN,1) = ABS(ZX1) VALDE1(1,IN,2) = ABS(ZY1) VALDE1(1,IN,3) = ABS(ZZ1) VALDE1(2,IN,1) = ABS(ZX2) VALDE1(2,IN,2) = ABS(ZY2) VALDE1(2,IN,3) = ABS(ZZ2) C VALDE2(1,IN,1) = PPX1 VALDE2(1,IN,2) = PPY1 VALDE2(1,IN,3) = PPZ1 VALDE2(2,IN,1) = PPX2 VALDE2(2,IN,2) = PPY2 VALDE2(2,IN,3) = PPZ2 C VALDE3(1,IN,1) = ABS(ZP1) VALDE3(2,IN,1) = ABS(ZP2) C VALDE4(1,IN,1) = PPP1 VALDE4(2,IN,1) = PPP2 C END IF C 110 CONTINUE C 100 CONTINUE C END