C KPCOQ6 SOURCE CHAT 05/01/13 01:04:08 5004 SUBROUTINE KPCOQ6(XX, P, XKP, IANT) C C Procedure de calcul de la matrice Kppour un element COQ4 C Entrees : XX(3, 6) : REAL*8 : Coordonnees des noeuds C XP : REAL : Pression C IANT : INTEGER : 1 si calcul asymétrique, 0 sinon C Sortie : XKP(36, 36) : REAL*8 : Matrice Kp elementaire C C (D'apres "Design variations of nonlinear elastic structures C subjected to follower forces" C M.J. Poldneff, I.S. Rai, J.S. Arora) C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION XX(3, 6), XKP(36, 36), SKRO(3, 3, 3), 1 XN(8, 30), XNB1(30), XNB2(30), XNB3(30), XNB4(30), XNB(30), 2 DNB1(30), DNB2(30), DN(2, 8, 30), DX(2, 3, 30), 3 XTIN(30), A(30), B(30), C(30), D(30), E(30), F(30), 4 XNB5(30), XNB6(30), XNB7(30), XDUM(36, 36) DATA XNB1/0.D0, 4.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ DATA XNB2/0.D0, 0.D0, 4.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ DATA XNB3/0.D0, 0.D0, 0.D0, 4.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ DATA XNB4/1.D0, -1.D0, -1.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ DATA XNB5/1.D0, -2.D0, -2.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ DATA XNB6/-0.25, 0.5, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ DATA XNB7/-0.25, 0.D0, 0.5, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ C C Initialisation du symbole epsilon C DO 10 I = 1, 3 DO 10 J = 1, 3 DO 10 K = 1, 3 10 SKRO(I, J, K) = 0.D0 SKRO(1, 2, 3) = 1.D0 SKRO(1, 3, 2) = -1.D0 SKRO(2, 3, 1) = 1.D0 SKRO(2, 1, 3) = -1.D0 SKRO(3, 1, 2) = 1.D0 SKRO(3, 2, 1) = -1.D0 C C Calcul de la pression C * P = 0. * DO 20 I = 1, 3 * 20 P = P + XP(I) * P = P/3. C C Fonctions de forme et derivees C C Les coefficients sont ranges comme suit : C indice : 1 2 3 4 5 6 7 8 9 C terme : 1 T1 T2 T1*T2 T1^2 T2^2 T1*T2^2 T1^2*T2 T1^3 C indice : 10 11 12 13 14 15 C terme : T2^3 T1*T2^3 T1^2*T2^2 T1^3*T2^3 T1^4 T2^4 C indice : 16 17 18 19 20 21 C terme : T1*T2^4 T1^2*T2^3 T1^3*T2^2 T1^4*T2 T1^5 T2^5 C indice : 22 23 24 25 26 C terme : T1*T2^5 T1^2*T2^4 T1^3*T2^3 T1^4*T2^2 T1^5*T2 C indice : 27 28 29 30 C terme : T1^2*T2^5 T1^3*T2^4 T1^4*T2^3 T1^5*T2^2 C CALL MULQP2(XNB4, XNB5, XNB) CALL DERQP2(XNB, 1, DNB1) CALL DERQP2(XNB, 2, DNB2) DO 31 I = 1, 30 XN(1, I) = XNB(I) DN(1, 1, I) = DNB1(I) DN(2, 1, I) = DNB2(I) 31 CONTINUE CALL MULQP2(XNB1, XNB4, XNB) CALL DERQP2(XNB, 1, DNB1) CALL DERQP2(XNB, 2, DNB2) DO 32 I = 1, 30 XN(2, I) = XNB(I) DN(1, 2, I) = DNB1(I) DN(2, 2, I) = DNB2(I) 32 CONTINUE CALL MULQP2(XNB1, XNB6, XNB) CALL DERQP2(XNB, 1, DNB1) CALL DERQP2(XNB, 2, DNB2) DO 33 I = 1, 30 XN(3, I) = XNB(I) DN(1, 3, I) = DNB1(I) DN(2, 3, I) = DNB2(I) 33 CONTINUE CALL DERQP2(XNB3, 1, DNB1) CALL DERQP2(XNB3, 2, DNB2) DO 34 I = 1, 30 XN(4, I) = XNB3(I) DN(1, 4, I) = DNB1(I) DN(2, 4, I) = DNB2(I) 34 CONTINUE CALL MULQP2(XNB2, XNB7, XNB) CALL DERQP2(XNB, 1, DNB1) CALL DERQP2(XNB, 2, DNB2) DO 35 I = 1, 30 XN(5, I) = XNB(I) DN(1, 5, I) = DNB1(I) DN(2, 5, I) = DNB2(I) 35 CONTINUE CALL MULQP2(XNB2, XNB4, XNB) CALL DERQP2(XNB, 1, DNB1) CALL DERQP2(XNB, 2, DNB2) DO 36 I = 1, 30 XN(6, I) = XNB(I) DN(1, 6, I) = DNB1(I) DN(2, 6, I) = DNB2(I) 36 CONTINUE C C Vecteurs tangents a la coque dans la configuration initiale C C Initialisation DO 40 I = 1, 2 C Boucle sur les parametres DO 40 J = 1, 3 C Boucle sur les composantes DO 40 K = 1, 30 C Boucle sur les coefficients des polynomes DX(I, J, K) = 0.D0 40 CONTINUE C Calcul DO 50 I = 1, 2 C Boucle sur les parametres DO 50 J = 1, 3 C Boucle sur les composantes DO 50 K = 1, 30 C Boucle sur les coefficients des polynomes DO 50 L = 1, 6 C Boucle sur les noeuds DX(I, J, K) = DX(I, J, K) + XX(J, L)*DN(I, L, K) 50 CONTINUE C C Calcul des termes de la matrice Kp C C Initialisation DO 55 I = 1, 36 DO 55 J = 1, 36 XDUM(I,J) = 0.D0 55 XKP(I, J) = 0.D0 C Calcul DO 60 II = 1, 3 DO 60 IL = 1, 6 DO 60 IS = 1, 3 DO 60 IT = 1, 6 XRES = 0.D0 DO 70 J = 1, 3 IFLAG = 0 DO 71 K = 1, 30 A(K) = XN(IL, K) D(K) = 0.D0 71 CONTINUE IF (SKRO(II, J, IS) .NE. 0.D0) THEN DO 81 K = 1, 30 B(K) = DN(2, IT, K) C(K) = DX(1, J, K) 81 CONTINUE CALL MULQP2(B, C, E) DO 72 K = 1, 30 72 D(K) = SKRO(II, J, IS)*E(K) IFLAG = 1 ENDIF IF (SKRO(II, IS, J) .NE. 0.D0) THEN DO 82 K = 1, 30 B(K) = DN(1, IT, K) C(K) = DX(2, J, K) 82 CONTINUE CALL MULQP2(B, C, E) IF (IFLAG .EQ. 1) THEN DO 73 K = 1, 30 73 D(K) = D(K) + SKRO(II, IS, J)*E(K) ELSE DO 74 K = 1, 30 74 D(K) = SKRO(II, IS, J)*E(K) ENDIF IFLAG = 1 ENDIF IF (IFLAG .NE. 0) THEN CALL MULQP2(A, D, XTIN) CALL INT6P2(XTIN, TING) XRES = XRES + TING ENDIF 70 CONTINUE XRES = XRES*P 60 XDUM(6*(IL-1) + II, 6*(IT-1) + IS) = XRES IF (IANT .EQ. 0) THEN DO 90 II = 1, 36 DO 90 IJ = II, 36 XKP(II, IJ) = (XDUM(II, IJ) + XDUM(IJ, II))*0.5D0 90 XKP(IJ, II) = XKP(II, IJ) ELSE DO 91 II = 1, 36 DO 91 IJ = 1, 36 91 XKP(II, IJ) = XDUM(II, IJ) ENDIF RETURN END