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