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