1 : C BCOQ4O SOURCE CHAT 06/03/29 21:15:44 5360 2 : SUBROUTINE BCOQ4O 3 : & (IGAU,XEL,SHPTOT,SHP,BGENE,DJAC,EXCEN,NOPLAN,IARR,IFS) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8(A-H,O-Z) 6 : ************************************************************************ 7 : * 8 : * B C O Q 4 9 : * --------- 10 : * 11 : * FONCTION: 12 : * --------- 13 : * 14 : * Calcul de la matrice "B" pour un COQ4 de comprtement non-isotrope, 15 : * excentré ou non. 16 : * 17 : * PARAMETRES: (E)=ENTREE (S)=SORTIE 18 : * ----------- 19 : * 20 : * IGAU (E) Numéro du point d'intégration. 21 : * XEL (E) Coordonnées locales des noeuds de l'élément. 22 : * SHPTOT (E) Fonctions de forme et dérivées dans l'espace de 23 : * référence. 24 : * SHPTOT(1,... = fonction, 25 : * SHPTOT(2,... = dérivée par rapport à Qsi, 26 : * SHPTOT(3,... = dérivée par rapport à Eta. 27 : * EXCEN (E) Excentrement. 28 : * NOPLAN (E) = 1 si élément non plan, 29 : * = 0 sinon. 30 : * IFS (E) = 1 si demande de Matrice "B" complète 31 : * (pour SIGMA et EPSILON), 32 : * = 0 sinon. 33 : * SHP (S) Fonctions de forme et dérivées dans l'espace 34 : * géométrique. 35 : * SHP(1,... = fonction, 36 : * SHP(2,... = dérivée par rapport à "X" local, 37 : * SHP(3,... = dérivée par rapport à "Y" local. 38 : * - au point d'intégration n.5 si IFS=1, 39 : * - au point d'intégration IGAU courant sinon. 40 : * DJAC (S) Jacobien au point d'intégration IGAU. 41 : * BGENE (S) Matrice "B" au point d'intégration IGAU. 42 : * IARR (S) Différent de 0 si erreur. 43 : * 44 : PARAMETER (LRE=24,NST=8,NBNO=4) 45 : REAL*8 XEL(3,*),BGENE(NST,*),SHP(6,*),SHPTOT(6,NBNO,*) 46 : * 47 : * CONSTANTES: 48 : * ----------- 49 : * 50 : * NBNO = Nombre de noeuds. 51 : * LRE = Nombre de colonnes de la Matrice "B". 52 : * NST = Nombre de composantes de contraintes. 53 : * 54 : * VARIABLES: 55 : * ---------- 56 : * 57 : * FULL = .TRUE. si on veut une Matrice "B" complète. 58 : * GAUSS = .TRUE. si utilisation d'un point d'intégration n.1 à 4. 59 : * SHP7 = Fonctions d'interpolation mixte pour le cisaillement X-Z 60 : * SHP8 = Fonctions d'interpolation mixte pour le cisaillement Y-Z 61 : * 62 : REAL*8 SHP7(4:5,NBNO),SHP8(4:5,NBNO) 63 : LOGICAL FULL,GAUSS 64 : * 65 : * MODE DE FONCTIONNEMENT: 66 : * ----------------------- 67 : * 68 : * POUR LA RIGIDITE: 69 : * 70 : * Pour la membrane, 71 : * - intégration pleine si élément plan, 72 : * - intégration réduite si élément non plan. 73 : * 74 : * Pour le cisaillement transverse, voir: 75 : * DONEA / LAMAIN (1987) - ISPRA 76 : * (élément de type "A.N.S.") 77 : * 78 : * Pour le couplage membrane-flexion (excentrement), intégration 79 : * réduite analogue à membrane, non justifiée au niveau de la 80 : * Rigidité, mais indispensable pour le calcul pratique des forces 81 : * nodales à partir des contraintes, par "BSIGMA". 82 : * 83 : * POUR LES CONTRAINTES: 84 : * 85 : * Cependant, pour que le champ de contraintes produit soit facile 86 : * à lire par l'utilisateur et à traiter par les opérateurs tels 87 : * que: 88 : * - CHANGER CHPOINT, 89 : * - TRACER, ... 90 : * on inscrit aux points de Gauss 1 à 4 les valeurs de contrainte 91 : * au point d'intégration n.5 pour les parties concernées. 92 : * De plus, on remplit complètement le champ de valeurs au point 93 : * d'intégration n.5 . 94 : * 95 : * AUTEUR, DATE DE CREATION: 96 : * ------------------------- 97 : * 98 : * Michel BULIK 26 Novembre 1996 d'après BCOQ4 écrit par 99 : * Pascal MANIGOT 17 Juin 1991 100 : * 101 : ************************************************************************ 102 : * 103 : IARR=0 104 : FULL = IFS .EQ. 1 105 : GAUSS = IGAU .LE. 4 106 : * 107 : * 108 : * Matrice jacobienne: 109 : * ------------------- 110 : * 111 : DO 110 I=1,NBNO 112 : SHP(1,I)=SHPTOT(1,I,IGAU) 113 : SHP(2,I)=SHPTOT(2,I,IGAU) 114 : SHP(3,I)=SHPTOT(3,I,IGAU) 115 : 110 CONTINUE 116 : * END DO 117 : * 118 : DXDQSI = 0.D0 119 : DXDETA = 0.D0 120 : DYDQSI = 0.D0 121 : DYDETA = 0.D0 122 : DO 120 I=1,NBNO 123 : DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I) 124 : DXDETA=DXDETA+SHP(3,I)*XEL(1,I) 125 : DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I) 126 : DYDETA=DYDETA+SHP(3,I)*XEL(2,I) 127 : 120 CONTINUE 128 : * END DO 129 : DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI 130 : IF(DJAC.LE. 0.D0) THEN 131 : IARR=1 132 : RETURN 133 : END IF 134 : * 135 : * Par inversion analytique de la matrice J nous avons les dérivées 136 : * 137 : DQSIDX= DYDETA/DJAC 138 : DQSIDY=-DXDETA/DJAC 139 : DETADX=-DYDQSI/DJAC 140 : DETADY= DXDQSI/DJAC 141 : * 142 : * puis les dérivées des fonctions de forme dans l'espace géométrique: 143 : * 144 : DO 130 I=1,NBNO 145 : TT=SHP(2,I)*DQSIDX+SHP(3,I)*DETADX 146 : SHP(3,I)=SHP(2,I)*DQSIDY+SHP(3,I)*DETADY 147 : SHP(2,I)=TT 148 : 130 CONTINUE 149 : * END DO 150 : * 151 : * 152 : * Matrice "B": 153 : * ------------ 154 : * 154 : CALL ZERO(BGENE,NST,LRE) 155 : * 156 : IF (GAUSS .OR. FULL) THEN 157 : * Préparation pour cisaillement transverse: 158 : * 159 : N1 = 1 160 : N2 = 2 161 : DX1 = SHPTOT(2,N1,IGAU)*XEL(1,N1)+SHPTOT(2,N2,IGAU)*XEL(1,N2) 162 : DY1 = SHPTOT(2,N1,IGAU)*XEL(2,N1)+SHPTOT(2,N2,IGAU)*XEL(2,N2) 163 : N1 = 3 164 : N2 = 4 165 : DX2 = SHPTOT(2,N1,IGAU)*XEL(1,N1)+SHPTOT(2,N2,IGAU)*XEL(1,N2) 166 : DY2 = SHPTOT(2,N1,IGAU)*XEL(2,N1)+SHPTOT(2,N2,IGAU)*XEL(2,N2) 167 : N1 = 1 168 : N2 = 4 169 : DX3 = SHPTOT(3,N1,IGAU)*XEL(1,N1)+SHPTOT(3,N2,IGAU)*XEL(1,N2) 170 : DY3 = SHPTOT(3,N1,IGAU)*XEL(2,N1)+SHPTOT(3,N2,IGAU)*XEL(2,N2) 171 : N1 = 2 172 : N2 = 3 173 : DX4 = SHPTOT(3,N1,IGAU)*XEL(1,N1)+SHPTOT(3,N2,IGAU)*XEL(1,N2) 174 : DY4 = SHPTOT(3,N1,IGAU)*XEL(2,N1)+SHPTOT(3,N2,IGAU)*XEL(2,N2) 175 : * 176 : D2JAC = 2.D0 * DJAC 177 : * 178 : SHP7(4,1) = (- DYDETA*DY1 + DYDQSI*DY3) / D2JAC 179 : SHP7(4,2) = (- DYDETA*DY1 + DYDQSI*DY4) / D2JAC 180 : SHP7(4,3) = (- DYDETA*DY2 + DYDQSI*DY4) / D2JAC 181 : SHP7(4,4) = (- DYDETA*DY2 + DYDQSI*DY3) / D2JAC 182 : SHP7(5,1) = ( DYDETA*DX1 - DYDQSI*DX3) / D2JAC 183 : SHP7(5,2) = ( DYDETA*DX1 - DYDQSI*DX4) / D2JAC 184 : SHP7(5,3) = ( DYDETA*DX2 - DYDQSI*DX4) / D2JAC 185 : SHP7(5,4) = ( DYDETA*DX2 - DYDQSI*DX3) / D2JAC 186 : * 187 : SHP8(4,1) = ( DXDETA*DY1 - DXDQSI*DY3) / D2JAC 188 : SHP8(4,2) = ( DXDETA*DY1 - DXDQSI*DY4) / D2JAC 189 : SHP8(4,3) = ( DXDETA*DY2 - DXDQSI*DY4) / D2JAC 190 : SHP8(4,4) = ( DXDETA*DY2 - DXDQSI*DY3) / D2JAC 191 : SHP8(5,1) = (- DXDETA*DX1 + DXDQSI*DX3) / D2JAC 192 : SHP8(5,2) = (- DXDETA*DX1 + DXDQSI*DX4) / D2JAC 193 : SHP8(5,3) = (- DXDETA*DX2 + DXDQSI*DX4) / D2JAC 194 : SHP8(5,4) = (- DXDETA*DX2 + DXDQSI*DX3) / D2JAC 195 : * 196 : END IF 197 : * 198 : DO 300 I=1,NBNO 199 : K = 6*(I-1) 200 : * 201 : IF (GAUSS .OR. FULL) THEN 202 : IF (NOPLAN .EQ. 0) THEN 203 : * Membrane: 204 : BGENE(1,K+1)=SHP(2,I) 205 : BGENE(2,K+2)=SHP(3,I) 206 : * Couplage membrane-flexion: 207 : BGENE(1,K+5) = SHP(2,I) * EXCEN 208 : BGENE(2,K+4) = -SHP(3,I) * EXCEN 209 : * Cisaillement de membrane et de couplage membrane-flexion: 210 : BGENE(3,K+1)=SHP(3,I) 211 : BGENE(3,K+2)=SHP(2,I) 212 : BGENE(3,K+4) = -SHP(2,I) * EXCEN 213 : BGENE(3,K+5) = SHP(3,I) * EXCEN 214 : ENDIF 215 : * Flexion: 216 : BGENE(4,K+5)=SHP(2,I) 217 : BGENE(5,K+4)=-SHP(3,I) 218 : BGENE(6,K+4)=-SHP(2,I) 219 : BGENE(6,K+5)=SHP(3,I) 220 : * Cisaillement transverse: 221 : BGENE(7,K+3)=SHP(2,I) 222 : BGENE(7,K+4) = SHP7(4,I) 223 : BGENE(7,K+5) = SHP7(5,I) 224 : BGENE(8,K+3)=SHP(3,I) 225 : BGENE(8,K+4) = SHP8(4,I) 226 : BGENE(8,K+5) = SHP8(5,I) 227 : ENDIF 228 : * 229 : IF (.NOT. GAUSS) THEN 230 : IF (NOPLAN.EQ.1 .OR. FULL) THEN 231 : * Membrane: 232 : BGENE(1,K+1)=SHP(2,I) 233 : BGENE(2,K+2)=SHP(3,I) 234 : * Couplage membrane-flexion: 235 : BGENE(1,K+5) = SHP(2,I) * EXCEN 236 : BGENE(2,K+4) = -SHP(3,I) * EXCEN 237 : * Cisaillement de membrane et de couplage membrane-flexion: 238 : BGENE(3,K+1)=SHP(3,I) 239 : BGENE(3,K+2)=SHP(2,I) 240 : BGENE(3,K+4) = -SHP(2,I) * EXCEN 241 : BGENE(3,K+5) = SHP(3,I) * EXCEN 242 : ENDIF 243 : ENDIF 244 : * 245 : 300 CONTINUE 246 : * END DO 247 : * 248 : * 249 : IF (GAUSS .AND. FULL .AND. NOPLAN .EQ. 1) THEN 250 : * On finit de compléter la Matrice "B", pour les contraintes ou 251 : * les déformations, par des valeurs calculées au centre: 252 : * 253 : DO 500 I=1,NBNO 254 : SHP(1,I)=SHPTOT(1,I, 5 ) 255 : SHP(2,I)=SHPTOT(2,I, 5 ) 256 : SHP(3,I)=SHPTOT(3,I, 5 ) 257 : 500 CONTINUE 258 : * END DO 259 : * FONCTIONS DE FORME DANS LA GÉOMÉTRIE RÉELLE: 259 : CALL JACOBI(XEL,SHP,2,NBNO,DJAC5) 260 : IF(DJAC5 .LE. 0.D0) THEN 261 : IARR=1 262 : RETURN 263 : END IF 264 : * 265 : DO 520 I=1,NBNO 266 : K = 6*(I-1) 267 : * Membrane: 268 : BGENE(1,K+1)=SHP(2,I) 269 : BGENE(2,K+2)=SHP(3,I) 270 : * Couplage membrane-flexion: 271 : BGENE(1,K+5) = SHP(2,I) * EXCEN 272 : BGENE(2,K+4) = -SHP(3,I) * EXCEN 273 : * Cisaillement de membrane et de couplage membrane-flexion: 274 : BGENE(3,K+1)=SHP(3,I) 275 : BGENE(3,K+2)=SHP(2,I) 276 : BGENE(3,K+4) = -SHP(2,I) * EXCEN 277 : BGENE(3,K+5) = SHP(3,I) * EXCEN 278 : 520 CONTINUE 279 : * END DO 280 : * 281 : ENDIF 282 : * 283 : ccc write(*,*) 'BCOQ4O: Matrice B =' 284 : ccc write(*,'(8(1X,1P,G10.3))') ((bgene(iii,jjj),jjj=1,lre),iii=1,nst) 285 : END 286 : 287 : 288 :
© Cast3M 2003 - All rights reserved.
Disclaimer