C BCOQ4O SOURCE CHAT 06/03/29 21:15:44 5360 SUBROUTINE BCOQ4O & (IGAU,XEL,SHPTOT,SHP,BGENE,DJAC,EXCEN,NOPLAN,IARR,IFS) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) ************************************************************************ * * B C O Q 4 * --------- * * FONCTION: * --------- * * Calcul de la matrice "B" pour un COQ4 de comprtement non-isotrope, * excentré ou non. * * PARAMETRES: (E)=ENTREE (S)=SORTIE * ----------- * * IGAU (E) Numéro du point d'intégration. * XEL (E) Coordonnées locales des noeuds de l'élément. * SHPTOT (E) Fonctions de forme et dérivées dans l'espace de * référence. * SHPTOT(1,... = fonction, * SHPTOT(2,... = dérivée par rapport à Qsi, * SHPTOT(3,... = dérivée par rapport à Eta. * EXCEN (E) Excentrement. * NOPLAN (E) = 1 si élément non plan, * = 0 sinon. * IFS (E) = 1 si demande de Matrice "B" complète * (pour SIGMA et EPSILON), * = 0 sinon. * SHP (S) Fonctions de forme et dérivées dans l'espace * géométrique. * SHP(1,... = fonction, * SHP(2,... = dérivée par rapport à "X" local, * SHP(3,... = dérivée par rapport à "Y" local. * - au point d'intégration n.5 si IFS=1, * - au point d'intégration IGAU courant sinon. * DJAC (S) Jacobien au point d'intégration IGAU. * BGENE (S) Matrice "B" au point d'intégration IGAU. * IARR (S) Différent de 0 si erreur. * PARAMETER (LRE=24,NST=8,NBNO=4) REAL*8 XEL(3,*),BGENE(NST,*),SHP(6,*),SHPTOT(6,NBNO,*) * * CONSTANTES: * ----------- * * NBNO = Nombre de noeuds. * LRE = Nombre de colonnes de la Matrice "B". * NST = Nombre de composantes de contraintes. * * VARIABLES: * ---------- * * FULL = .TRUE. si on veut une Matrice "B" complète. * GAUSS = .TRUE. si utilisation d'un point d'intégration n.1 à 4. * SHP7 = Fonctions d'interpolation mixte pour le cisaillement X-Z * SHP8 = Fonctions d'interpolation mixte pour le cisaillement Y-Z * REAL*8 SHP7(4:5,NBNO),SHP8(4:5,NBNO) LOGICAL FULL,GAUSS * * MODE DE FONCTIONNEMENT: * ----------------------- * * POUR LA RIGIDITE: * * Pour la membrane, * - intégration pleine si élément plan, * - intégration réduite si élément non plan. * * Pour le cisaillement transverse, voir: * DONEA / LAMAIN (1987) - ISPRA * (élément de type "A.N.S.") * * Pour le couplage membrane-flexion (excentrement), intégration * réduite analogue à membrane, non justifiée au niveau de la * Rigidité, mais indispensable pour le calcul pratique des forces * nodales à partir des contraintes, par "BSIGMA". * * POUR LES CONTRAINTES: * * Cependant, pour que le champ de contraintes produit soit facile * à lire par l'utilisateur et à traiter par les opérateurs tels * que: * - CHANGER CHPOINT, * - TRACER, ... * on inscrit aux points de Gauss 1 à 4 les valeurs de contrainte * au point d'intégration n.5 pour les parties concernées. * De plus, on remplit complètement le champ de valeurs au point * d'intégration n.5 . * * AUTEUR, DATE DE CREATION: * ------------------------- * * Michel BULIK 26 Novembre 1996 d'après BCOQ4 écrit par * Pascal MANIGOT 17 Juin 1991 * ************************************************************************ * IARR=0 FULL = IFS .EQ. 1 GAUSS = IGAU .LE. 4 * * * Matrice jacobienne: * ------------------- * DO 110 I=1,NBNO SHP(1,I)=SHPTOT(1,I,IGAU) SHP(2,I)=SHPTOT(2,I,IGAU) SHP(3,I)=SHPTOT(3,I,IGAU) 110 CONTINUE * END DO * DXDQSI = 0.D0 DXDETA = 0.D0 DYDQSI = 0.D0 DYDETA = 0.D0 DO 120 I=1,NBNO DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I) DXDETA=DXDETA+SHP(3,I)*XEL(1,I) DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I) DYDETA=DYDETA+SHP(3,I)*XEL(2,I) 120 CONTINUE * END DO DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI IF(DJAC.LE. 0.D0) THEN IARR=1 RETURN END IF * * Par inversion analytique de la matrice J nous avons les dérivées * DQSIDX= DYDETA/DJAC DQSIDY=-DXDETA/DJAC DETADX=-DYDQSI/DJAC DETADY= DXDQSI/DJAC * * puis les dérivées des fonctions de forme dans l'espace géométrique: * DO 130 I=1,NBNO TT=SHP(2,I)*DQSIDX+SHP(3,I)*DETADX SHP(3,I)=SHP(2,I)*DQSIDY+SHP(3,I)*DETADY SHP(2,I)=TT 130 CONTINUE * END DO * * * Matrice "B": * ------------ * CALL ZERO(BGENE,NST,LRE) * IF (GAUSS .OR. FULL) THEN * Préparation pour cisaillement transverse: * N1 = 1 N2 = 2 DX1 = SHPTOT(2,N1,IGAU)*XEL(1,N1)+SHPTOT(2,N2,IGAU)*XEL(1,N2) DY1 = SHPTOT(2,N1,IGAU)*XEL(2,N1)+SHPTOT(2,N2,IGAU)*XEL(2,N2) N1 = 3 N2 = 4 DX2 = SHPTOT(2,N1,IGAU)*XEL(1,N1)+SHPTOT(2,N2,IGAU)*XEL(1,N2) DY2 = SHPTOT(2,N1,IGAU)*XEL(2,N1)+SHPTOT(2,N2,IGAU)*XEL(2,N2) N1 = 1 N2 = 4 DX3 = SHPTOT(3,N1,IGAU)*XEL(1,N1)+SHPTOT(3,N2,IGAU)*XEL(1,N2) DY3 = SHPTOT(3,N1,IGAU)*XEL(2,N1)+SHPTOT(3,N2,IGAU)*XEL(2,N2) N1 = 2 N2 = 3 DX4 = SHPTOT(3,N1,IGAU)*XEL(1,N1)+SHPTOT(3,N2,IGAU)*XEL(1,N2) DY4 = SHPTOT(3,N1,IGAU)*XEL(2,N1)+SHPTOT(3,N2,IGAU)*XEL(2,N2) * D2JAC = 2.D0 * DJAC * SHP7(4,1) = (- DYDETA*DY1 + DYDQSI*DY3) / D2JAC SHP7(4,2) = (- DYDETA*DY1 + DYDQSI*DY4) / D2JAC SHP7(4,3) = (- DYDETA*DY2 + DYDQSI*DY4) / D2JAC SHP7(4,4) = (- DYDETA*DY2 + DYDQSI*DY3) / D2JAC SHP7(5,1) = ( DYDETA*DX1 - DYDQSI*DX3) / D2JAC SHP7(5,2) = ( DYDETA*DX1 - DYDQSI*DX4) / D2JAC SHP7(5,3) = ( DYDETA*DX2 - DYDQSI*DX4) / D2JAC SHP7(5,4) = ( DYDETA*DX2 - DYDQSI*DX3) / D2JAC * SHP8(4,1) = ( DXDETA*DY1 - DXDQSI*DY3) / D2JAC SHP8(4,2) = ( DXDETA*DY1 - DXDQSI*DY4) / D2JAC SHP8(4,3) = ( DXDETA*DY2 - DXDQSI*DY4) / D2JAC SHP8(4,4) = ( DXDETA*DY2 - DXDQSI*DY3) / D2JAC SHP8(5,1) = (- DXDETA*DX1 + DXDQSI*DX3) / D2JAC SHP8(5,2) = (- DXDETA*DX1 + DXDQSI*DX4) / D2JAC SHP8(5,3) = (- DXDETA*DX2 + DXDQSI*DX4) / D2JAC SHP8(5,4) = (- DXDETA*DX2 + DXDQSI*DX3) / D2JAC * END IF * DO 300 I=1,NBNO K = 6*(I-1) * IF (GAUSS .OR. FULL) THEN IF (NOPLAN .EQ. 0) THEN * Membrane: BGENE(1,K+1)=SHP(2,I) BGENE(2,K+2)=SHP(3,I) * Couplage membrane-flexion: BGENE(1,K+5) = SHP(2,I) * EXCEN BGENE(2,K+4) = -SHP(3,I) * EXCEN * Cisaillement de membrane et de couplage membrane-flexion: BGENE(3,K+1)=SHP(3,I) BGENE(3,K+2)=SHP(2,I) BGENE(3,K+4) = -SHP(2,I) * EXCEN BGENE(3,K+5) = SHP(3,I) * EXCEN ENDIF * Flexion: BGENE(4,K+5)=SHP(2,I) BGENE(5,K+4)=-SHP(3,I) BGENE(6,K+4)=-SHP(2,I) BGENE(6,K+5)=SHP(3,I) * Cisaillement transverse: BGENE(7,K+3)=SHP(2,I) BGENE(7,K+4) = SHP7(4,I) BGENE(7,K+5) = SHP7(5,I) BGENE(8,K+3)=SHP(3,I) BGENE(8,K+4) = SHP8(4,I) BGENE(8,K+5) = SHP8(5,I) ENDIF * IF (.NOT. GAUSS) THEN IF (NOPLAN.EQ.1 .OR. FULL) THEN * Membrane: BGENE(1,K+1)=SHP(2,I) BGENE(2,K+2)=SHP(3,I) * Couplage membrane-flexion: BGENE(1,K+5) = SHP(2,I) * EXCEN BGENE(2,K+4) = -SHP(3,I) * EXCEN * Cisaillement de membrane et de couplage membrane-flexion: BGENE(3,K+1)=SHP(3,I) BGENE(3,K+2)=SHP(2,I) BGENE(3,K+4) = -SHP(2,I) * EXCEN BGENE(3,K+5) = SHP(3,I) * EXCEN ENDIF ENDIF * 300 CONTINUE * END DO * * IF (GAUSS .AND. FULL .AND. NOPLAN .EQ. 1) THEN * On finit de compléter la Matrice "B", pour les contraintes ou * les déformations, par des valeurs calculées au centre: * DO 500 I=1,NBNO SHP(1,I)=SHPTOT(1,I, 5 ) SHP(2,I)=SHPTOT(2,I, 5 ) SHP(3,I)=SHPTOT(3,I, 5 ) 500 CONTINUE * END DO * FONCTIONS DE FORME DANS LA GÉOMÉTRIE RÉELLE: CALL JACOBI(XEL,SHP,2,NBNO,DJAC5) IF(DJAC5 .LE. 0.D0) THEN IARR=1 RETURN END IF * DO 520 I=1,NBNO K = 6*(I-1) * Membrane: BGENE(1,K+1)=SHP(2,I) BGENE(2,K+2)=SHP(3,I) * Couplage membrane-flexion: BGENE(1,K+5) = SHP(2,I) * EXCEN BGENE(2,K+4) = -SHP(3,I) * EXCEN * Cisaillement de membrane et de couplage membrane-flexion: BGENE(3,K+1)=SHP(3,I) BGENE(3,K+2)=SHP(2,I) BGENE(3,K+4) = -SHP(2,I) * EXCEN BGENE(3,K+5) = SHP(3,I) * EXCEN 520 CONTINUE * END DO * ENDIF * ccc write(*,*) 'BCOQ4O: Matrice B =' ccc write(*,'(8(1X,1P,G10.3))') ((bgene(iii,jjj),jjj=1,lre),iii=1,nst) END