1 : C KPCOQ6 SOURCE CHAT 05/01/13 01:04:08 5004 2 : SUBROUTINE KPCOQ6(XX, P, XKP, IANT) 3 : C 4 : C Procedure de calcul de la matrice Kppour un element COQ4 5 : C Entrees : XX(3, 6) : 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(36, 36) : 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, 6), XKP(36, 36), SKRO(3, 3, 3), 17 : 1 XN(8, 30), XNB1(30), XNB2(30), XNB3(30), XNB4(30), XNB(30), 18 : 2 DNB1(30), DNB2(30), DN(2, 8, 30), DX(2, 3, 30), 19 : 3 XTIN(30), A(30), B(30), C(30), D(30), E(30), F(30), 20 : 4 XNB5(30), XNB6(30), XNB7(30), XDUM(36, 36) 21 : DATA XNB1/0.D0, 4.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 22 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 23 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 24 : 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 25 : 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ 26 : DATA XNB2/0.D0, 0.D0, 4.D0, 0.D0, 0.D0, 0.D0, 0.D0, 27 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 28 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 29 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 30 : 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ 31 : DATA XNB3/0.D0, 0.D0, 0.D0, 4.D0, 0.D0, 0.D0, 0.D0, 32 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 33 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 34 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 35 : 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ 36 : DATA XNB4/1.D0, -1.D0, -1.D0, 0.D0, 0.D0, 0.D0, 0.D0, 37 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 38 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 39 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 40 : 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ 41 : DATA XNB5/1.D0, -2.D0, -2.D0, 0.D0, 0.D0, 0.D0, 42 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 43 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 44 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 45 : 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ 46 : DATA XNB6/-0.25, 0.5, 0.D0, 0.D0, 0.D0, 0.D0, 47 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 48 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 49 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 50 : 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ 51 : DATA XNB7/-0.25, 0.D0, 0.5, 0.D0, 0.D0, 0.D0, 52 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 53 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 54 : 1 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 55 : 2 0.D0, 0.D0, 0.D0, 0.D0, 0.D0/ 56 : C 57 : C Initialisation du symbole epsilon 58 : C 59 : DO 10 I = 1, 3 60 : DO 10 J = 1, 3 61 : DO 10 K = 1, 3 62 : 10 SKRO(I, J, K) = 0.D0 63 : SKRO(1, 2, 3) = 1.D0 64 : SKRO(1, 3, 2) = -1.D0 65 : SKRO(2, 3, 1) = 1.D0 66 : SKRO(2, 1, 3) = -1.D0 67 : SKRO(3, 1, 2) = 1.D0 68 : SKRO(3, 2, 1) = -1.D0 69 : C 70 : C Calcul de la pression 71 : C 72 : * P = 0. 73 : * DO 20 I = 1, 3 74 : * 20 P = P + XP(I) 75 : * P = P/3. 76 : C 77 : C Fonctions de forme et derivees 78 : C 79 : C Les coefficients sont ranges comme suit : 80 : C indice : 1 2 3 4 5 6 7 8 9 81 : C terme : 1 T1 T2 T1*T2 T1^2 T2^2 T1*T2^2 T1^2*T2 T1^3 82 : C indice : 10 11 12 13 14 15 83 : C terme : T2^3 T1*T2^3 T1^2*T2^2 T1^3*T2^3 T1^4 T2^4 84 : C indice : 16 17 18 19 20 21 85 : C terme : T1*T2^4 T1^2*T2^3 T1^3*T2^2 T1^4*T2 T1^5 T2^5 86 : C indice : 22 23 24 25 26 87 : C terme : T1*T2^5 T1^2*T2^4 T1^3*T2^3 T1^4*T2^2 T1^5*T2 88 : C indice : 27 28 29 30 89 : C terme : T1^2*T2^5 T1^3*T2^4 T1^4*T2^3 T1^5*T2^2 90 : C 90 : CALL MULQP2(XNB4, XNB5, XNB) 90 : CALL DERQP2(XNB, 1, DNB1) 90 : CALL DERQP2(XNB, 2, DNB2) 91 : DO 31 I = 1, 30 92 : XN(1, I) = XNB(I) 93 : DN(1, 1, I) = DNB1(I) 94 : DN(2, 1, I) = DNB2(I) 95 : 31 CONTINUE 95 : CALL MULQP2(XNB1, XNB4, XNB) 95 : CALL DERQP2(XNB, 1, DNB1) 95 : CALL DERQP2(XNB, 2, DNB2) 96 : DO 32 I = 1, 30 97 : XN(2, I) = XNB(I) 98 : DN(1, 2, I) = DNB1(I) 99 : DN(2, 2, I) = DNB2(I) 100 : 32 CONTINUE 100 : CALL MULQP2(XNB1, XNB6, XNB) 100 : CALL DERQP2(XNB, 1, DNB1) 100 : CALL DERQP2(XNB, 2, DNB2) 101 : DO 33 I = 1, 30 102 : XN(3, I) = XNB(I) 103 : DN(1, 3, I) = DNB1(I) 104 : DN(2, 3, I) = DNB2(I) 105 : 33 CONTINUE 105 : CALL DERQP2(XNB3, 1, DNB1) 105 : CALL DERQP2(XNB3, 2, DNB2) 106 : DO 34 I = 1, 30 107 : XN(4, I) = XNB3(I) 108 : DN(1, 4, I) = DNB1(I) 109 : DN(2, 4, I) = DNB2(I) 110 : 34 CONTINUE 110 : CALL MULQP2(XNB2, XNB7, XNB) 110 : CALL DERQP2(XNB, 1, DNB1) 110 : CALL DERQP2(XNB, 2, DNB2) 111 : DO 35 I = 1, 30 112 : XN(5, I) = XNB(I) 113 : DN(1, 5, I) = DNB1(I) 114 : DN(2, 5, I) = DNB2(I) 115 : 35 CONTINUE 115 : CALL MULQP2(XNB2, XNB4, XNB) 115 : CALL DERQP2(XNB, 1, DNB1) 115 : CALL DERQP2(XNB, 2, DNB2) 116 : DO 36 I = 1, 30 117 : XN(6, I) = XNB(I) 118 : DN(1, 6, I) = DNB1(I) 119 : DN(2, 6, I) = DNB2(I) 120 : 36 CONTINUE 121 : C 122 : C Vecteurs tangents a la coque dans la configuration initiale 123 : C 124 : C Initialisation 125 : DO 40 I = 1, 2 126 : C Boucle sur les parametres 127 : DO 40 J = 1, 3 128 : C Boucle sur les composantes 129 : DO 40 K = 1, 30 130 : C Boucle sur les coefficients des polynomes 131 : DX(I, J, K) = 0.D0 132 : 40 CONTINUE 133 : C Calcul 134 : DO 50 I = 1, 2 135 : C Boucle sur les parametres 136 : DO 50 J = 1, 3 137 : C Boucle sur les composantes 138 : DO 50 K = 1, 30 139 : C Boucle sur les coefficients des polynomes 140 : DO 50 L = 1, 6 141 : C Boucle sur les noeuds 142 : DX(I, J, K) = DX(I, J, K) + XX(J, L)*DN(I, L, K) 143 : 50 CONTINUE 144 : C 145 : C Calcul des termes de la matrice Kp 146 : C 147 : C Initialisation 148 : DO 55 I = 1, 36 149 : DO 55 J = 1, 36 150 : XDUM(I,J) = 0.D0 151 : 55 XKP(I, J) = 0.D0 152 : C Calcul 153 : DO 60 II = 1, 3 154 : DO 60 IL = 1, 6 155 : DO 60 IS = 1, 3 156 : DO 60 IT = 1, 6 157 : XRES = 0.D0 158 : DO 70 J = 1, 3 159 : IFLAG = 0 160 : DO 71 K = 1, 30 161 : A(K) = XN(IL, K) 162 : D(K) = 0.D0 163 : 71 CONTINUE 164 : IF (SKRO(II, J, IS) .NE. 0.D0) THEN 165 : DO 81 K = 1, 30 166 : B(K) = DN(2, IT, K) 167 : C(K) = DX(1, J, K) 168 : 81 CONTINUE 168 : CALL MULQP2(B, C, E) 169 : DO 72 K = 1, 30 170 : 72 D(K) = SKRO(II, J, IS)*E(K) 171 : IFLAG = 1 172 : ENDIF 173 : IF (SKRO(II, IS, J) .NE. 0.D0) THEN 174 : DO 82 K = 1, 30 175 : B(K) = DN(1, IT, K) 176 : C(K) = DX(2, J, K) 177 : 82 CONTINUE 177 : CALL MULQP2(B, C, E) 178 : IF (IFLAG .EQ. 1) THEN 179 : DO 73 K = 1, 30 180 : 73 D(K) = D(K) + SKRO(II, IS, J)*E(K) 181 : ELSE 182 : DO 74 K = 1, 30 183 : 74 D(K) = SKRO(II, IS, J)*E(K) 184 : ENDIF 185 : IFLAG = 1 186 : ENDIF 187 : IF (IFLAG .NE. 0) THEN 187 : CALL MULQP2(A, D, XTIN) 187 : CALL INT6P2(XTIN, TING) 188 : XRES = XRES + TING 189 : ENDIF 190 : 70 CONTINUE 191 : XRES = XRES*P 192 : 60 XDUM(6*(IL-1) + II, 6*(IT-1) + IS) = XRES 193 : IF (IANT .EQ. 0) THEN 194 : DO 90 II = 1, 36 195 : DO 90 IJ = II, 36 196 : XKP(II, IJ) = (XDUM(II, IJ) + XDUM(IJ, II))*0.5D0 197 : 90 XKP(IJ, II) = XKP(II, IJ) 198 : ELSE 199 : DO 91 II = 1, 36 200 : DO 91 IJ = 1, 36 201 : 91 XKP(II, IJ) = XDUM(II, IJ) 202 : ENDIF 203 : RETURN 204 : END 205 : 206 :
© Cast3M 2003 - All rights reserved.
Disclaimer