C ECOU21 SOURCE OF166741 25/11/04 21:15:47 12349 SUBROUTINE ECOU21(MATE,INPLAS,MELE,IPMAIL,NBPTEL,IMAT,ICAR, 1 NUMAT,NUCAR,IVASTR,IVARI,IVADEF,IVADET, 1 IVADS,IVAMAT,IVACAR, 2 IPH1,IPH2,IPH3,ITHHER,LHOOK,NSTRS,NVARI,NMATT,NCARR, 3 CMATE,PRECIS,JECHER,IPOTAB,ISTEP,NPINT,JNOID,LOGSUC, 4 N2EL,N2PTEL,NBNO,NBPGAU,LW,IVASTF,IVARIF,IVADEP,KERRE) *********************************************************************** * CAS DES CERAMIQUES *********************************************************************** * entrees : * * mate = numero de materiau elastique * inplas = numero de materiau inelastique * mele = numero element fini * ipmail = pointeur du maillage * nbptel = nombre de points par element * imat = pointeur sur un segment mptval de materiau (utilise par calsig) * icar = pointeur sur un segment mptval de caracteristiques * geometriques (utilise par calsig) * numat = nb de composantes du melval de imat * nucar = nb de composantes du melval de icar * ivastr =pointeur sur un segment mptval de contraintes * ivari =pointeur sur un segment mptval de variables internes * ivadef =pointeur sur un segment mptval de deformations * ivads =pointeur sur un segment mptval de contraintes (increments) * ivamat =pointeur sur un segment mptval de materiau * ivacar =pointeur sur un segment mptval de cacarteristiques geometrique * iph1 = pointeur sur un mchaml de temperatures au debut du pas * iph2 = pointeur sur un mchaml de temperatures a la fin du pas * iph3 = pointeur sur un mchaml de temperatures de reference * ithher = 0 si pas de chargement thermique * = 1 si chargement thermique mais materiau constant * = 2 si chargement thermique et mat. dependant de la temperature * ipch1,ipch2,ipch3,ithher ne servent que pour les materiaux * endommageables de lemaitre quand ils dependent de la temperature * lhook =taille de la matrice de hooke * nstrs =nombre de composantes de contraintes * nvari =nombre de composantes de variables internes * nmatt =nombre de composnates de proprietes de materiau * ncarr =nombre de composnates de caracteristiques geometriques * cmate =nom du materiau * precis =precision dans les iterations internes * jecher =0 ou 1 pour action dans ecoule * jnoid =0 ou 1 pour action dans ecoule * ipotab =pointeur sur segment table * istep =indicateur d'action pour calcul nonlocal * =0 dans le cas d'un calcul local (normal) * =1 ou 2 dans le cas d'un calcul nonlocal * =1 pour calcul des fonctions seuil uniquement * =2 pour calcul des variables dissipatives a partir * des fonctions seuil moyennees prealablement par nloc * * sorties : * ivastf =pointeur sur un segment mptval de contraintes * ivarif =pointeur sur un segment mptval de variables internes * ivadep =pointeur sur un segment mptval de deformations inelastiques * kerre =indicateur d'erreur * * p dowlatyari fev. 1992 * * c. la borderie fev 92 restructuration et reecriture de certains * passages pour une meilleure lisibilite * * avril 92 ajout istep pour le non local * dec 92 modif pour poutres timoschenko * ************************************************************************ IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMCHAML -INC SMELEME -INC SMCOORD -INC SMMODEL -INC SMINTE -INC CCHAMP -INC CECOU c======================================================================= c la variable kerre regit les impressions d erreurs dans plast c toutes erreurs de ecoule gerees dans ce sous programme c kerre=0 tout ok c de 1 a 6 s aligner sur valeurs donnees par ecoinc c = 7 un element tuyau a une epaisseur nulle c = 21 on ne trouve pas d intersection avec la surface de charge c = 22 sig0 a l exterieur de la surface de charge c c anomalies avec la courbe de traction c = 30 limite elastique nulle c = 31 trop de points c = 32 pas assez de points c = 33 pente incorrecte c = 34 module d'young nul c = 35 manque l'origine c = 36 pente a l'origine non egale a e c = 37 manque la courbe de traction c = 38 nu devrait etre nul c c = 48 donnees erronnees pour drucker-prager c = 49 matrice singuliere dans iter internes drucker-prager c = 51 pb dans drucker prager option non disponible c = 52 pb dans drucker prager donnees incompatibles c = 53 pb dans drucker prager solution impossible c = 54 les valeurs admissibles pour istep sont 0 1 ou 2 c = 55 modele non implante en non local c = 56 probleme dans l'integration du modele mazars c = 57 .... c = 58 .... c = 59 .... c = 60 pb donnees du cam-clay c c = 99 cas non encore disponible c======================================================================= -INC TMPTVAL SEGMENT WRK0 REAL*8 XMAT(NCXMAT) ENDSEGMENT * SEGMENT WR00 CHARACTER*16 TYMAT(NCXMAT) REAL*8 XMAT1(NCXMAT),XMAT2(NCXMAT) ENDSEGMENT * SEGMENT WRK1 REAL*8 DDHOOK(LHOOK,LHOOK),SIG0(NSTRS),DEPST(NSTRS) REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI) REAL*8 DEFP(NSTRS),XCAR(ICARA) ENDSEGMENT * SEGMENT WRK2 REAL*8 TRAC(LTRAC) ENDSEGMENT * SEGMENT WRK22 REAL*8 XXE(3,NBNN) ENDSEGMENT * SEGMENT WRK3 REAL*8 WORK(LW),WORK2(LW2) ENDSEGMENT * SEGMENT WRK4 REAL*8 XE(3,NBBB) ENDSEGMENT * SEGMENT WRK5 REAL*8 EPIN0(NSTRS),EPINF(NSTRS),EPST0(NSTRS) ENDSEGMENT * SEGMENT WRK6 REAL*8 BB(NSTRS,NNVARI),R(NSTRS),XMU(NSTRS) REAL*8 S(NNVARI),QSI(NNVARI),DDR(NSTRS),BBS(NSTRS) REAL*8 SIGMA(NSTRS),SIGGD(NSTRS),XMULT(NSTRS),PROD(NSTRS) ENDSEGMENT * SEGMENT WRK7 REAL*8 F(NCOURB,2),W(NCOURB),TRUC(NCOURB) ENDSEGMENT * SEGMENT WRK8 REAL*8 DD(NSTRS,NSTRS),DDV(NSTRS,NSTRS),DDINV(NSTRS,NSTRS) ENDSEGMENT * SEGMENT WRK9 REAL*8 YOG(NYOG),YNU(NYNU),YALFA(NYALFA),YSMAX(NYSMAX) REAL*8 YN(NYN),YM(NYM),YKK(NYKK),YALFA1(NYALF1) REAL*8 YBETA1(NYBET1),YR(NYR),YA(NYA),YKX(NYKX),YRHO(NYRHO) REAL*8 SIGY(NSIGY) INTEGER NKX(NNKX) ENDSEGMENT * SEGMENT WR10 INTEGER IABLO1(NTABO1) REAL*8 TABLO2(NTABO2) ENDSEGMENT * SEGMENT WR11 INTEGER IABLO3(NTABO3) REAL*8 TABLO4(NTABO4) ENDSEGMENT * SEGMENT WR12 REAL*8 EM0(2,NWA(1)),EM1(2,NWA(2)),EM2(2,NWA(3)) REAL*8 EM3(2,NWA(4)),EM4(2,NWA(5)),EM5(2,NWA(6)) REAL*8 EM6(2,NWA(7)),EM7(2,NWA(8)),EM8(2,NWA(9)) REAL*8 SM0(NSTRS),SM1(NSTRS),SM2(NSTRS),SM3(NSTRS) REAL*8 SM4(NSTRS),SM5(NSTRS),SM6(NSTRS),SM7(NSTRS) REAL*8 SM8(NSTRS) ENDSEGMENT * SEGMENT WTRAV REAL*8 DDAUX(LHOOK,LHOOK),VALMAT(NUMAT) REAL*8 VALCAR(NUCAR),DSIGT(NSTRS) REAL*8 TXR(IDIM,IDIM),DDHOMU(LHOOK,LHOOK) REAL*8 XLOC(3,3),XGLOB(3,3) REAL*8 D1HOOK(LHOOK,LHOOK),ROTHOO(LHOOK,LHOOK) ENDSEGMENT * SEGMENT WPOUT REAL*8 X(2),Y(2),Z(2) ENDSEGMENT LOGICAL LOGVIS,LOGIN,LOGRE,LOGSUC LOGICAL LUNI1,LUNI2 DIMENSION BID(6),BID2(6),CRIGI(12),CMASS(12) DIMENSION NWA(9) DIMENSION SIG01(4),VAR01(36) CHARACTER*72 CHARRE CHARACTER*8 CMATE * * mise à disposition des temperatures tini tfin tref * aux points de gauss * TETA1=-1.E35 TETA2=-1.E35 TETREF=-1.E35 TREFA=-1.E35 IF (ITHHER.EQ.1.OR.ITHHER.EQ.2) THEN MCHAM3=IPH1 MCHAM4=IPH2 MCHAM5=IPH3 SEGACT MCHAM3,MCHAM4,MCHAM5 MELVA3=MCHAM3.IELVAL(1) MELVA4=MCHAM4.IELVAL(1) MELVA5=MCHAM5.IELVAL(1) SEGACT MELVA3,MELVA4,MELVA5 ENDIF c c Initialisations de variables c--------------------------------- WRK8 = 0 WRK4 = 0 WRK22 = 0 minte2 = 0 c - mise à zéro des variables du commun NECOU si besoin c - modèles viscoplastiques: c . on récupère le pas de temps c . on récupère le nombre maximal de sous-pas c . on met IND=1 c - initialisation des dimensions des tableaux des segments c Sorties: en plus du commun NECOU, on range les autres données c initialisées dans les COMMON IECOU et XECOU c Sauf pour KERRE,LW,LOGVIS,LUNI1 et LUNI2 qui sont sortis comme c argument de DEFINI c CALL DEFINI(MELE,NCARR,NSTRS,NMATT,CMATE,MATE, . ISTEP,INPLAS,NPINT,IPOTAB,IVADEF, . IPMAIL,IVAMAT, . ITHHER,NUMAT,NUCAR,LOGVIS, . LUNI1,LUNI2,LW,KERRE) IF (KERRE.EQ.999) RETURN c IF ((CMATE.EQ.'ORTHOTRO'.OR.CMATE.EQ.'ANISOTRO'.OR. 1 CMATE.EQ.'UNIDIREC').AND.(MFR.EQ.1.OR.MFR.EQ.31 1 .OR.MFR.EQ.33)) THEN CALL RESHPT(1,NBNO,IELE,MELE,NPINT,IPTR1,IRT1) MINTE2=IPTR1 SEGACT,MINTE2 SEGINI WRK22 ENDIF c IF (LOGVIS) SEGINI WRK8 SEGINI WRK0,WR00,WRK1,WRK2,WRK3,WRK5 IF(MFR.EQ.7.OR.MFR.EQ.13.OR.LUNI1)THEN SEGINI WRK4 ENDIF c SEGINI WTRAV * * boucle sur les elements * DO 1000 IB=1,NBELEM * * Matériaux orthotropes, anisotropes et unidirectionnels * en formulation massive: * - on cherche les coordonnees des noeuds de l element ib * - calcul des axes locaux * Cas particulier de l'ACIER_UNI * CALL DEFROT(CMATE,MFR,NBNN,IB,MELE,LUNI1,IPTR1, . MELEME,WRK4,WRK22,WTRAV) * * boucle sur les points de gauss * DO 1100 IGAU=1,NBPTEL * * -recuperation de valmat et de valcar * -on recupere les contraintes initiales * -on recupere les variables internes * -on recupere les deformations inelastiques initiales si besoin * -on recupere les increments de deformations totales * -on cherche la section de l'element ib * -prise en compte de l'epaisseur et de l'excentrement * dans le cas des coques minces avec ou sans cisaillement * transverse * CALL DEFVAL(NUMAT1,NBPTEL,NDEF, . IMAT,IVACAR,ICAR,IVASTR,IVARI,IVADEF,IVADET, . IVADS,MFR,CMATE,INPLAS,IB,IGAU,IND, . WTRAV,WRK1,WRK5,SECT,EPAIST) * * on recupere les constantes du materiau * en cas de reels, on a directement les valeurs * en cas d'objets, on a les pointeurs eu guise de valeurs * et on calcule les contraintes effectives en milieu poreux * CALL DEFMAT(NMATT,NSTRS,MFR,MELE,INPLAS, . IVAMAT,IB,IGAU,CMATE,MATE,LUNI1,LUNI2, . WRK1,WRK5,WRK0,WR00,WTRAV,CMASS,CRIGI,COB,XMOB, . BID,BID2,KERR0) IF (KERR0.EQ.99) THEN KERRE=99 GOTO 1000 ELSE IF (KERR0.EQ.10) THEN GOTO 1000 ENDIF * * >>>>>>>>>> fin du traitement du materiau * * on recupere les caracteristiques geometriques * CALL DEFCAR(NCARR,IB,IGAU,MFR,MELE,IVACAR, . WRK1) * * quelques impressions si iimpi = 99 * IF(IIMPI.EQ.99) THEN * WRITE(IOIMP,66770) IB,IGAU *66770 format(////////2x,'element ',i6,2x,'point ',i3//) * WRITE(IOIMP,66771) MATE,INPLAS *66771 format('0 mate=',i4,2x,'inplas=',i4/) * WRITE(IOIMP,66772) (SIG0(I),I=1,NSTRS) *66772 format(2x,' sig0 '/(6(1x,1pe12.5))) * WRITE(IOIMP,66773) (VAR0(I),I=1,NVARI) *66773 format(2x,' var0 '/(6(1x,1pe12.5))) * WRITE(IOIMP,66774) (DEPST(I),I=1,NSTRS) *66774 format(2x,' depst '/(6(1x,1pe12.5))) WRITE(IOIMP,66775) (XMAT(I),I=1,NMATT) 66775 format(2x,' xmat '/(6(1x,1pe12.5))) * IF(IVACAR.NE.0)THEN * WRITE(IOIMP,66776) (XCAR(I),I=1,ICARA) *66776 format(2x,' xcar '/(6(1x,1pe12.5))) * ENDIF ENDIF * * mise à disposition des temperatures tini tfin tref * aux points de gauss * IF (ITHHER.EQ.1.OR.ITHHER.EQ.2) THEN IBMN=MIN(IB,MELVA3.VELCHE(/2)) IGMN=MIN(IGAU,MELVA3.VELCHE(/1)) TETA1=MELVA3.VELCHE(IGMN,IBMN) IBMN=MIN(IB,MELVA4.VELCHE(/2)) IGMN=MIN(IGAU,MELVA4.VELCHE(/1)) TETA2=MELVA4.VELCHE(IGMN,IBMN) IBMN=MIN(IB,MELVA5.VELCHE(/2)) IGMN=MIN(IGAU,MELVA5.VELCHE(/1)) TETREF=MELVA5.VELCHE(IGMN,IBMN) ENDIF * *--------------------------------------------------------------------- * * ecoulement * *--------------------------------------------------------------------- IF (INPLAS.EQ.65) THEN * SEGINI WRK7 SEGINI WRK9 IF((MFR.EQ.1).AND.(IFOMOD.EQ.2)) THEN IBIDO = 19 ELSE IBIDO = 14 ENDIF * CAS OU ON NE PREND PAS EN COMPTE LA TEMPERATURE DE TRANSITION * CAD LORSQUE TTRAN = 0 * IF ((XMAT(IBIDO).LE.0.1).AND.(XMAT(IBIDO).GE.-0.1)) THEN * * si le point de gauss est déjà endommagé par endommagement généralisé * on le traite simplement par ceraca IF (VAR0(NVARI-1).EQ.1) THEN CALL CERACA(WRK0,WRK1,WRK5,WTRAV,INPLAS,MFR1,DT,NSTRSS, 1 NVARI,PRECIS,MSOUPA,JECHER,DTT,NSSINC,INV,KERRE, 2 ICARA,IFOURB,CMATE,N2EL,N2PTEL,IB,IGAU,EPAIST, 3 NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK, 4 CRIGI) IND=1 ELSE * si le point de gauss n'a pas un endommagement généralisé * on regarde si il a été fissuré * par ottosen et si non on applique le fluage puis ottosen * si oui on le traite par Ottosen MPTVAL=IVAMAT CALL OTOBO(VAR0,XMAT,IVAL,ITOTO,MFR) IF (ITOTO.EQ.0) THEN CALL CERACA(WRK0,WRK1,WRK5,WTRAV,INPLAS,MFR1,DT,NSTRSS, 1 NVARI,PRECIS,MSOUPA,JECHER,DTT,NSSINC,INV, 2 KERRE,ICARA,IFOURB,CMATE,N2EL,N2PTEL,IB,IGAU, 3 EPAIST,NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT, 4 LHOOK,CRIGI) IND=1 * Ligne suivante à supprimer * IF(IND.EQ.0) THEN * on regarde si on a eu endommagement généralisé * si on n'a pas eu endommagement généralisé on appele ottosen IF (VARF(NVARI-1).NE.1) THEN DO 161 I = 1,NVARI VAR01(I) = VARF(I) 161 CONTINUE DO 535 I=1,NSTRS * PRINT *,'DEPST EPINF-EPIN0 ',I,DEPST(I),(EPINF(I)-EPIN0(I)) DEPST(I) = DEPST(I) -( EPINF(I)-EPIN0(I)) C On remplace SIGF par SIG0 SIG01(I) = SIG0(I) 535 CONTINUE MPTVAL=IVAMAT CALL OTTOSE(INPLAS,SIG01,NSTRSS,DEPST,VAR01,XMAT,IVAL, & NMATT,XCAR,ICARA,NVARI,SIGF,VARF,DEFP,MFR1,KERRE, & IB,IGAU) C on met à jour le variable interne EPSE commune aux deux modèles VARF(1) = VARF(1)+VARF(NVARI) C On calcule l'increment de déformation du pas de temps DO 536 I=1,NSTRS DEFP(I) =DEFP(I)+( EPINF(I)-EPIN0(I)) 536 CONTINUE IND=0 ENDIF * Ligne suivante à supprimer * ENDIF ELSE MPTVAL=IVAMAT CALL OTTOSE(INPLAS,SIG0,NSTRSS,DEPST,VAR0,XMAT,IVAL, & NMATT,XCAR,ICARA,NVARI,SIGF,VARF,DEFP,MFR1,KERRE, & IB,IGAU) VARF(1) = VARF(1)+VARF(NVARI) IND=0 ENDIF ENDIF * ELSE * * CAS OU ON PREND EN COMPTE LA TEMP2RATURE DE TRANSITION * IF(TETA2.GE.XMAT(IBIDO)) THEN MPTVAL=IVAMAT CALL OTOBO(VAR0,XMAT,IVAL,ITOTO,MFR) IF (ITOTO.EQ.0) THEN CALL CERACA(WRK0,WRK1,WRK5,WTRAV,INPLAS,MFR1,DT,NSTRSS, 1 NVARI,PRECIS,MSOUPA,JECHER,DTT,NSSINC,INV, 2 KERRE,ICARA,IFOURB,CMATE,N2EL,N2PTEL,IB,IGAU, 3 EPAIST,NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT, 4 LHOOK,CRIGI) IND=1 ELSE MPTVAL=IVAMAT CALL OTTOSE(INPLAS,SIG0,NSTRSS,DEPST,VAR0,XMAT,IVAL, & NMATT,XCAR,ICARA,NVARI,SIGF,VARF,DEFP,MFR1,KERRE, & IB,IGAU) VARF(1) = VARF(1)+VARF(NVARI) IND=0 ENDIF ELSE IF (VAR0(NVARI-1).EQ.1) THEN CALL CERACA(WRK0,WRK1,WRK5,WTRAV,INPLAS,MFR1,DT,NSTRSS, 1 NVARI,PRECIS,MSOUPA,JECHER,DTT,NSSINC,INV,KERRE, 2 ICARA,IFOURB,CMATE,N2EL,N2PTEL,IB,IGAU,EPAIST, 3 NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK, 4 CRIGI) IND=1 ELSE MPTVAL=IVAMAT CALL OTTOSE(INPLAS,SIG0,NSTRSS,DEPST,VAR0,XMAT,IVAL, & NMATT,XCAR,ICARA,NVARI,SIGF,VARF,DEFP,MFR1,KERRE, & IB,IGAU) VARF(1) = VARF(1)+VARF(NVARI) IND=0 ENDIF ENDIF ENDIF IF (MFR1.EQ.17) THEN IF (KERRE.NE.0.AND.NSSINC.EQ.1) THEN CALL ERREUR(KERRE) ENDIF ENDIF SEGSUP WRK7 SEGSUP WRK9 DTOPTI = MIN(DTOPTI,DTT) NINCMA = MAX(NINCMA,NSSINC) NCOMP = NCOMP + 1 TSOM = TSOM + DTT NSOM = NSOM + NSSINC NINV = NINV + INV TCAR = TCAR + DTT* DTT IF(KERRE.NE.0.AND.KERRE.NE.99) THEN KERR1=1 ENDIF c ELSE IF (INPLAS.EQ.74) THEN * * CHAINE DE MAXWELL * * on commence par recuperer le nombre d'elements dans la chaine * et les proprietes et variables internes associees a des objets CALL MAXTRA(WRK0,WRK1,WRK5,WR12,WTRAV,IB,IGAU, & NBGMAT,NELMAT,NPINT,NWA,NSTRSS,NCHAIN,CMATE,MFR) IF(IERR.NE.0) THEN SEGSUP WR12 GOTO 1789 ENDIF IF (MFR.EQ.3.OR.MFR.EQ.39) THEN CALL MAXGEN(WRK0,WRK1,WRK5,WR12,MFR, 1 IB,IGAU,MELE,NCHAIN,KERRE,DT,CMATE,NWA,TEMP0) ELSE * * MLR 10/08/99 * * ON PASSE LE SEGMENT DE TRAVAIL WTRAV * CALL MAXWEL(WRK0,WRK1,WRK5,WR12,MFR, + IB,IGAU,MELE,NCHAIN,KERRE,DT,CMATE,NWA,TEMP0) ENDIF * * ici gerer les erreurs * CALL MAXTRB(WTRAV,WRK1,WRK5,WR12,NWA,NSTRSS, & NCHAIN,CMATE) SEGSUP WR12 * * FIN DES DIFFERENTS MODELES * ELSE KERRE = 99 ENDIF * * Erreurs * - problèmes de convergence * 1789 CALL DEFER1(JNOID,KERR1,KERRE,LOGSUC) * * - autres problèmes * 1990 CALL DEFER2(INPLAS,MFR,MELE,IB,IGAU, . KERR1,KERRE) IF (KERRE.NE.0) GOTO 99 c c remplissage du segment contenant les contraintes a la fin * ( rearrangement pour milieu poreux ), c les variables internes finales c et les increments de deformations plastiques c CALL DEFSIG(MFR,NDEF, . INPLAS,IND,WRK1,WRK5,WTRAV, . IVASTF,IVARIF,IVADEP,COB,XMOB,IB,IGAU, . CMATE,MATE,MELE,KERRER) IF (KERRER.NE.0) GOTO 1000 c c fin de la boucle sur les points de gauss c 1100 continue c c special poutres et tuyaux sauf timoschenko c CALL DEFPOU(MFR,MELE,MELEME,IB,WRK4,IVASTF) c c fin de la boucle sur les elements c 1000 continue c * FIN: modèles visqueux, on stocke le pas de temps * optimal en indice 'dtopti' * CALL DEFFIN(INPLAS,TSOM,NSOM,NCOMP,NINV,NINCMA, . TCAR,DTOPTI,IPOTAB,KERRE) 99 CONTINUE IF (LOGVIS) SEGSUP WRK8 SEGSUP WRK0,WR00,WRK1,WRK2,WRK3,WRK5,WTRAV IF (WRK4.NE.0) SEGSUP WRK4 IF (WRK22.NE.0) THEN SEGDES MINTE2 SEGSUP WRK22 ENDIF IF (ITHHER.EQ.1.OR.ITHHER.EQ.2) THEN SEGDES MELVA3,MELVA4,MELVA5 SEGDES MCHAM3,MCHAM4,MCHAM5 ENDIF RETURN END