1 : C KPCOQ3 SOURCE CHAT 05/01/13 01:03:59 5004 2 : SUBROUTINE KPCOQ3(XX, XP, XKP, IANT) 3 : C 4 : C Procedure de calcul de la matrice Kppour un element COQ4 5 : C Entrees : XX(3, 3) : REAL*8 : Coordonnees des noeuds 6 : C XP : REAL : Pression 7 : C IANT : INTEGER : 1 si calcul asymétrique, 0 sinon 8 : C Sortie : XKP(18, 18) : REAL*8 : Matrice Kp elementaire 9 : C 10 : C (D'apres "Design variations of nonlinear elastic structures 11 : C subjected to follower forces" 12 : C M.J. Poldneff, I.S. Rai, J.S. Arora) 13 : C 14 : IMPLICIT INTEGER(I-N) 15 : IMPLICIT REAL*8(A-H,O-Z) 16 : DIMENSION XX(3, 3), XKP(18, 18), SKRO(3, 3, 3), 17 : 1 XN(3, 9), XNB1(9), XNB2(9), XNB3(9), XNB(9), 18 : 2 DNB1(9), DNB2(9), DN(2, 3, 9), DX(2, 3, 9), 19 : 3 XTIN(9), A(9), B(9), C(9), D(9), E(9), F(9), XDUM(18,18) 20 : DATA XNB1/1.D0, -1.D0, -1.D0, 0.D0, 0.D0, 21 : 1 0.D0, 0.D0, 0.D0, 0.D0/ 22 : DATA XNB2/0.D0, 0.D0, 1.D0, 0.D0, 0.D0, 23 : 1 0.D0, 0.D0, 0.D0, 0.D0/ 24 : DATA XNB3/0.D0, 1.D0, 0.D0, 0.D0, 0.D0, 25 : 1 0.D0, 0.D0, 0.D0, 0.D0/ 26 : C 27 : C Initialisation du symbole epsilon 28 : C 29 : DO 10 I = 1, 3 30 : DO 10 J = 1, 3 31 : DO 10 K = 1, 3 32 : 10 SKRO(I, J, K) = 0.D0 33 : SKRO(1, 2, 3) = 1.D0 34 : SKRO(1, 3, 2) = -1.D0 35 : SKRO(2, 3, 1) = 1.D0 36 : SKRO(2, 1, 3) = -1.D0 37 : SKRO(3, 1, 2) = 1.D0 38 : SKRO(3, 2, 1) = -1.D0 39 : C 40 : C Calcul de la pression 41 : C 42 : P = XP 43 : * DO 20 I = 1, 3 44 : * 20 P = P + XP(I) 45 : * P = P/4. 46 : * WRITE (*,*) 'Pression : ', P 47 : C 48 : C Fonctions de forme et derivees 49 : C 50 : C Les coefficients sont ranges comme suit : 51 : C indice : 1 2 3 4 5 6 7 8 9 52 : C terme : 1 T2 T1 T1*T2 T2^2 T1^2 T1*T2^2 T1^2*T2 T1^2*T2^2 53 : C 53 : CALL DERIP2(XNB1, 1, DNB1) 53 : CALL DERIP2(XNB1, 2, DNB2) 54 : DO 31 I = 1, 9 55 : XN(1, I) = XNB1(I) 56 : DN(1, 1, I) = DNB1(I) 57 : DN(2, 1, I) = DNB2(I) 58 : 31 CONTINUE 58 : CALL DERIP2(XNB2, 1, DNB1) 58 : CALL DERIP2(XNB2, 2, DNB2) 59 : DO 32 I = 1, 9 60 : XN(2, I) = XNB2(I) 61 : DN(1, 2, I) = DNB1(I) 62 : DN(2, 2, I) = DNB2(I) 63 : 32 CONTINUE 63 : CALL DERIP2(XNB3, 1, DNB1) 63 : CALL DERIP2(XNB3, 2, DNB2) 64 : DO 33 I = 1, 9 65 : XN(3, I) = XNB3(I) 66 : DN(1, 3, I) = DNB1(I) 67 : DN(2, 3, I) = DNB2(I) 68 : 33 CONTINUE 69 : C 70 : C Vecteurs tangents a la coque dans la configuration initiale 71 : C 72 : C Initialisation 73 : DO 40 I = 1, 2 74 : C Boucle sur les parametres 75 : DO 40 J = 1, 3 76 : C Boucle sur les composantes 77 : DO 40 K = 1, 9 78 : C Boucle sur les coefficients des polynomes 79 : DX(I, J, K) = 0.D0 80 : 40 CONTINUE 81 : C Calcul 82 : DO 50 I = 1, 2 83 : C Boucle sur les parametres 84 : DO 50 J = 1, 3 85 : C Boucle sur les composantes 86 : DO 50 K = 1, 9 87 : C Boucle sur les coefficients des polynomes 88 : DO 50 L = 1, 3 89 : C Boucle sur les noeuds 90 : DX(I, J, K) = DX(I, J, K) + XX(J, L)*DN(I, L, K) 91 : 50 CONTINUE 92 : C 93 : C Calcul des termes de la matrice Kp 94 : C 95 : C Initialisation 96 : DO 55 I = 1, 18 97 : DO 55 J = 1, 18 98 : XDUM(I,J) = 0.D0 99 : 55 XKP(I, J) = 0.D0 100 : C Calcul 101 : DO 60 II = 1, 3 102 : DO 60 IL = 1, 3 103 : DO 60 IS = 1, 3 104 : DO 60 IT = 1, 3 105 : XRES = 0.D0 106 : DO 70 J = 1, 3 107 : IFLAG = 0 108 : DO 71 K = 1, 9 109 : A(K) = XN(IL, K) 110 : D(K) = 0.D0 111 : 71 CONTINUE 112 : IF (SKRO(II, J, IS) .NE. 0.) THEN 113 : DO 81 K = 1, 9 114 : B(K) = DN(2, IT, K) 115 : C(K) = DX(1, J, K) 116 : 81 CONTINUE 116 : CALL MULTP2(B, C, E) 117 : DO 72 K = 1, 9 118 : 72 D(K) = SKRO(II, J, IS)*E(K) 119 : IFLAG = 1 120 : ENDIF 121 : IF (SKRO(II, IS, J) .NE. 0.) THEN 122 : DO 82 K = 1, 9 123 : B(K) = DN(1, IT, K) 124 : C(K) = DX(2, J, K) 125 : 82 CONTINUE 125 : CALL MULTP2(B, C, E) 126 : IF (IFLAG .EQ. 1) THEN 127 : DO 73 K = 1, 9 128 : 73 D(K) = D(K) + SKRO(II, IS, J)*E(K) 129 : ELSE 130 : DO 74 K = 1, 9 131 : 74 D(K) = SKRO(II, IS, J)*E(K) 132 : ENDIF 133 : IFLAG = 1 134 : ENDIF 135 : IF (IFLAG .NE. 0) THEN 135 : CALL MULTP2(A, D, XTIN) 135 : CALL INT3P2(XTIN, TING) 136 : XRES = XRES + TING 137 : ENDIF 138 : 70 CONTINUE 139 : XRES = XRES*P 140 : 60 XDUM(6*(IL-1) + II, 6*(IT-1) + IS) = XRES 141 : IF (IANT .EQ. 0) THEN 142 : DO 90 II = 1, 18 143 : DO 90 IJ = II, 18 144 : XKP(II, IJ) = (XDUM(II, IJ) + XDUM(IJ, II))*0.5D0 XRESA = 0.D0 145 : 146 : 90 XKP(IJ, II) = XKP(II, IJ) 147 : ELSE 148 : DO 91 II = 1, 18 149 : DO 91 IJ = 1, 18 150 : 91 XKP(II, IJ) = XDUM(II, IJ) 151 : ENDIF 152 : RETURN 153 : END 154 : 155 : 156 : 157 :
© Cast3M 2003 - All rights reserved.
Disclaimer