C AMADE2 SOURCE CHAT 05/01/12 21:21:02 5004 C AMADE2.eso SOURCE bald MAR 28/03/95 00:00:00 SUBROUTINE AMADE2(IB,IGAU,NSTRS,SIG0,EPIN0,VAR0,NVARI, & DEPST,IFOURB,XMAT,NMATT,IVAL,SIGF,DEFP,VARF,KERRE) C C----------------------------------------------------------------------- C C PLASTICITE MODELE 2D AMADEI-SAEB ELEMENTS JOINTS C C ENTREES C C IB = NUMERO DE L'ELEMENT C IGAU = NUMERO DU POINT DE GAUSS C NSTRS = NOMBRE DE COMPOSANTES DE CONTRAINTES C SIG0(NSTRS) = CONTRAINTES INITIALES (AU PAS PRECEDENT) C EPIN0(NSTRS) = DEFORMATIONS INITIALES INEL. (AU PAS PRECEDENT) C VAR0(NVARI) = VARIABLES INTERNES DEBUT C VAR0(1) = DEFORMATION PLASTIQUE EQUIVALENTE: EPSE C VAR0(2) = INCR. DEF. PLAST. DUE A L'OUVERTURE SEULE: EPOU C VAR0(3) = ETAT DU JOINT: STAT C VAR0(4) = U EQUIVALENTE: UEQU C VAR0(5) = DEFORMATION TOTALE: EPS1 C VAR0(6) = DEFORMATION TOTALE: EPS2 C NVARI = NOMBRE DE VARIABLES INTERNES C DEPST(NSTRS) = INCREMENTS DE DEFORMATION TOTALE C XMAT(NCOMAT) = COMPOSANTES DE MATERIAU C NMATT = NOMBRE DE COMPOSANTES DE MATERIAU C IVAL(NCOMAT) = INDICE DES COMPOSANTES DE MATERIAU C C SORTIES C C SIGF(NSTRS) = CONTRAINTES FINALES C DEFP(NSTRS) = INCREMENTS DE DEFORMATIONS INELASTIQUES FINALES C VARF(NVARI) = VARIABLES INTERNES FINALES C KERRE = INDICE D'ERREUR C C C----------------------------------------------------------------------- c c déclaration des variables c IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) C Include contenant quelques constantes dont XPI : -INC CCREEL REAL*8 I0, K1, K2, KNI REAL*8 KNN, KNT, KTN, KTT c DIMENSION SIG0(*), EPIN0(*), VAR0(*), DEPST(*), XMAT(*) DIMENSION SIGF(*), DEFP(*), VARF(*) DIMENSION RI0E(2,2), RI0P(2,2), DCON(2) c c paramètres pour le calcul c FATT=XPI/180.D0 c NN=100 c K1=1.5D0 K2=4.0D0 c c tolérances pour définir la condition de plasticité selon sigma c et la condition de nullité des contraintes c EPSS=1.D-6 EPSI=1.D-6 TOL1=1.D-20 TOL2=1.D-15 c c extraction des caractéristiques du materiau c FIMU=XMAT(5)*FATT SGMT=XMAT(6) I0=XMAT(7)*FATT S0=XMAT(8) B0=XMAT(9) UR=XMAT(10) UP=XMAT(11) KNI=XMAT(12) FI0=XMAT(13)*FATT VM=XMAT(14) c c extraction de la variable EPOU c EPOU=VAR0(2) c c extraction de la variable UEQU c UEQU=VAR0(4) c c extraction des déformations totales c EPS1=VAR0(5) EPS2=VAR0(6) c c calcul des petits pas c DU=DEPST(1)/NN DV=DEPST(2)/NN c c incréments de déformation plastique finale c DEFP(1)=0.D0 DEFP(2)=0.D0 c c contraintes finales c SIGF(1)=0.D0 SIGF(2)=0.D0 c c variables internes finales c VARF(1)=0.D0 VARF(2)=0.D0 VARF(3)=0.D0 VARF(4)=0.D0 VARF(5)=0.D0 VARF(6)=0.D0 c c ********************************************************************** c ********** test pour verifier la condition initiale du joint ********* c ********************************************************************** c IF(EPIN0(2).GT.EPSS)THEN c c le joint est ouvert c IF(DEPST(2).LT.(-EPIN0(2)-EPSI))THEN c c il y a rapprochement des lèvres avec compression (1 cas) c R0=ABS(EPIN0(2)/DEPST(2)) c NN=INT(ABS((EPIN0(2)+DEPST(2))/DV)) c EPOU=-EPIN0(2) c EPIN0(1)=EPIN0(1)+R0*DEPST(1) EPIN0(2)=0.D0 c DEFP(1)=R0*DEPST(1) DEFP(2)=EPOU c EPS1=EPS1+R0*DEPST(1) EPS2=EPS2+R0*DEPST(2) c ELSE IF(DEPST(2).GT.(-EPIN0(2)+EPSS))THEN c c il y a rapprochement des lèvres sans fermeture, ou elles restent c bloquées, ou elles s' éloignent (3 cas) c SIGF(1)=0.D0 SIGF(2)=0.D0 c DEFP(1)=DEPST(1) DEFP(2)=DEPST(2) c DEDE1=EPIN0(1)+DEPST(1) DEDE2=EPIN0(2)+DEPST(2) c VARF(1)=SQRT(2.D0/3.D0*(DEDE1**2+DEDE2**2)) c VARF(2)=DEPST(2) c VARF(3)=1.D0 c VARF(4)=VAR0(4) c VARF(5)=EPS1+DEPST(1) VARF(6)=EPS2+DEPST(2) c RETURN c ELSE c c il y a fermeture sans compression (1 cas) c SIGF(1)=0.D0 SIGF(2)=0.D0 c DEFP(1)=DEPST(1) DEFP(2)=DEPST(2) c DEDE1=EPIN0(1)+DEPST(1) DEDE2=EPIN0(2)+DEPST(2) c VARF(1)=SQRT(2.D0/3.D0*(DEDE1**2+DEDE2**2)) c VARF(2)=DEPST(2) c VARF(3)=0.D0 c VARF(4)=VAR0(4) c VARF(5)=EPS1+DEPST(1) VARF(6)=EPS2+DEPST(2) c RETURN c END IF c ELSE IF((EPIN0(2).LE.EPSS).AND.(SIG0(2).GE.-TOL2))THEN c c le joint est fermé mais sans compression c IF(DEPST(2).GT.EPSS)THEN c c les lèvres s' éloignent (1 cas) c SIGF(1)=0.D0 SIGF(2)=0.D0 c DEFP(1)=DEPST(1) DEFP(2)=DEPST(2) c DEDE1=EPIN0(1)+DEPST(1) DEDE2=EPIN0(2)+DEPST(2) c VARF(1)=SQRT(2.D0/3.D0*(DEDE1**2+DEDE2**2)) c VARF(2)=DEPST(2) c VARF(3)=1.D0 c VARF(4)=VAR0(4) c VARF(5)=EPS1+DEPST(1) VARF(6)=EPS2+DEPST(2) c RETURN c ELSE c c les lèvres restent bloquées ou il y a compression (2 cas) c END IF c ELSE c c le joint est fermé avec compression c END IF c c ********************************************************************** c ******************* fin du test et debut du calcul ******************* c ********************************************************************** c DO 300 KAPPA=1,NN c c definition des deux cas (sigma=SGMT) c IF(ABS(SIG0(2)).LT.ABS(SGMT))THEN c c cas avec sigma < SGMT c c prédiction élastique c IF(UEQU.LT.UR)THEN c c c' est le cas élastique si U < UR c c definition des quantités necessaires au calcul c VUPU=((1.D0-SIG0(2)/SGMT)**K2)*TAN(I0) AS=1.D0-((1.D0-SIG0(2)/SGMT)**K1) SR=S0-SIG0(2)*TAN(FI0) TP=-SIG0(2)*TAN(FIMU+ATAN(VUPU))*(1.D0-AS)+AS*SR TR=TP*(B0+(1.D0-B0)*SIG0(2)/SGMT) c AA1=(SIG0(2)*(1.D0-AS)*K2)/(SGMT*(COS(FIMU+ATAN(VUPU))**2 $ )) AA2=((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0) AA3=1.D0+((1.D0-SIG0(2)/SGMT)**(2.D0*K2))*(TAN(I0)**2) AA4=K1*SIG0(2)/SGMT AA5=((1.D0-SIG0(2)/SGMT)**(K1-1.D0))*TAN(FIMU+ATAN(VUPU)) $ *AA4 AA6=SR/SGMT*K1*((1.D0-SIG0(2)/SGMT)**(K1-1.D0)) AA7=AS*TAN(FI0) dTPS1=-(1.D0-AS)*TAN(FIMU+ATAN(VUPU))+AA1*AA2/AA3+AA5+AA6 $ -AA7 dTRS1=dTPS1*(B0+(1.D0-B0)*SIG0(2)/SGMT)+TP*(1.D0-B0)/SGMT AA8=-UEQU*K2/SGMT*((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0 $ ) AA8=AA8+((VM**2)*KNI)/((KNI*VM-SIG0(2))**2) c c calcul de la matrice Dt equivalente (2D) c KNN=1.D0/AA8 KNT=-((1.D0-SIG0(2)/SGMT)**K2)*TAN(I0)*KNN KTN=UEQU/UP*KNN*dTPS1 KTT=UEQU/UP*KNT*dTPS1+TP/UP c ELSE c c c' est le cas élastique si U >= UR c c definition des quantités necessaires au calcul c VUPU=0.D0 AS=1.D0-((1.D0-SIG0(2)/SGMT)**K1) SR=S0-SIG0(2)*TAN(FI0) TP=-SIG0(2)*TAN(FIMU)*(1.D0-AS)+AS*SR TR=TP*(B0+(1.D0-B0)*SIG0(2)/SGMT) c AA1=((1.D0-SIG0(2)/SGMT)**(K1-1.D0))*TAN(FIMU)*K1*SIG0(2) $ /SGMT AA2=SR/SGMT*K1*((1.D0-SIG0(2)/SGMT)**(K1-1.D0)) AA3=AS*TAN(FI0) dTPS1=-(1.D0-AS)*TAN(FIMU)+AA1+AA2-AA3 dTRS1=dTPS1*(B0+(1.D0-B0)*SIG0(2)/SGMT)+TP*(1.D0-B0)/SGMT AA4=-UR*K2/SGMT*((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0) AA4=AA4+((VM**2)*KNI)/((KNI*VM-SIG0(2))**2) c c calcul de la matrice Dt equivalente (2D) c KNN=1.D0/AA4 KNT=0.D0 KTN=UEQU/UP*KNN*dTPS1 KTT=TP/UP c END IF c c appel à CALC2 pour calculer les incréments de contraintes c CALL CALC2(KNN,KNT,KTN,KTT,SIG0,DU,DV,NSTRS, $ DELTA,RI0E,DCON) c c calcul de l' incrément de deformation plastique c DEP1=0.D0 DEP2=0.D0 c IF(UEQU.LT.UP)THEN c c il n' y a pas de plasticité c GO TO 200 c ELSE IF((UEQU.GE.UP).AND.(UEQU.LT.UR))THEN c c il y a de la plasticité: trait decroissant de la courbe tau c c test pour vérifier si continuer en plasticité ou reprendre le trait c lineaire c EFFE=ABS(SIG0(1))-(TR-TP)/(UR-UP)*(UEQU-UP)-TP c IF(ABS(SIG0(1)).LT.TOL1)THEN c DTEF=ABS(DCON(1)) c ELSE c DTEF=SIG0(1)*DCON(1)/ABS(SIG0(1)) c END IF c DTEF=DTEF-(TR-TP)/(UR-UP)*DELTA DTEF=DTEF-((UEQU-UP)/(UR-UP)*(dTRS1-dTPS1)+dTPS1)*DCON(2) c IF((EFFE.LT.0.D0).OR.(DTEF.LT.0.D0))GO TO 200 c c definition des quantités necessaires au calcul c VUPU=((1.D0-SIG0(2)/SGMT)**K2)*TAN(I0) AS=1.D0-((1.D0-SIG0(2)/SGMT)**K1) SR=S0-SIG0(2)*TAN(FI0) TP=-SIG0(2)*TAN(FIMU+ATAN(VUPU))*(1.D0-AS)+AS*SR TR=TP*(B0+(1.D0-B0)*SIG0(2)/SGMT) c AA1=(SIG0(2)*(1.D0-AS)*K2)/(SGMT*(COS(FIMU+ATAN(VUPU))**2 $ )) AA2=((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0) AA3=1.D0+((1.D0-SIG0(2)/SGMT)**(2.D0*K2))*(TAN(I0)**2) AA4=K1*SIG0(2)/SGMT AA5=((1.D0-SIG0(2)/SGMT)**(K1-1.D0))*TAN(FIMU+ATAN(VUPU)) $ *AA4 AA6=SR/SGMT*K1*((1.D0-SIG0(2)/SGMT)**(K1-1.D0)) AA7=AS*TAN(FI0) dTPS1=-(1.D0-AS)*TAN(FIMU+ATAN(VUPU))+AA1*AA2/AA3+AA5+AA6 $ -AA7 dTRS1=dTPS1*(B0+(1.D0-B0)*SIG0(2)/SGMT)+TP*(1.D0-B0)/SGMT AA8=-UEQU*K2/SGMT*((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0 $ ) AA8=AA8+((VM**2)*KNI)/((KNI*VM-SIG0(2))**2) c c calcul de la matrice Dt equivalente (2D) c KNN=1.D0/AA8 KNT=-((1.D0-SIG0(2)/SGMT)**K2)*TAN(I0)*KNN KTN=KNN/(UP-UR)*(dTPS1*(UEQU-UR)+(UP-UEQU)*dTRS1) KTT=(TP-TR)/(UP-UR) KTT=KTT+KNT/(UP-UR)*(dTPS1*(UEQU-UR)+(UP-UEQU)*dTRS1) c c appel à CALC2 pour calculer les incréments de contraintes c CALL CALC2(KNN,KNT,KTN,KTT,SIG0,DU,DV,NSTRS, $ DELTA,RI0P,DCON) c c appel à RISPL2 pour calculer les incréments de déformation plastique c CALL RISPL2(RI0E,SIG0,DCON,DU,NSTRS,SGMT,DUp) c DEP1=DUp DEP2=0.D0 c ELSE IF(UEQU.GE.UR)THEN c c il y a de la plasticité parfaite: trait horizontal de la courbe tau c c test pour vérifier si continuer en plasticité ou reprendre le trait c lineaire c EFFE=ABS(SIG0(1))-TR c IF(ABS(SIG0(1)).LT.TOL1)THEN c DTEF=ABS(DCON(1)) c ELSE c DTEF=SIG0(1)*DCON(1)/ABS(SIG0(1)) c END IF c DTEF=DTEF-dTRS1*DCON(2) c IF((EFFE.LT.0.D0).OR.(DTEF.LT.0.D0))GO TO 200 c c definition des quantités necessaires au calcul c VUPU=0.D0 AS=1.D0-((1.D0-SIG0(2)/SGMT)**K1) SR=S0-SIG0(2)*TAN(FI0) TP=-SIG0(2)*TAN(FIMU)*(1.D0-AS)+AS*SR TR=TP*(B0+(1.D0-B0)*SIG0(2)/SGMT) c AA1=((1.D0-SIG0(2)/SGMT)**(K1-1.D0))*TAN(FIMU)*K1*SIG0(2) $ /SGMT AA2=SR/SGMT*K1*((1.D0-SIG0(2)/SGMT)**(K1-1.D0)) AA3=AS*TAN(FI0) dTPS1=-(1.D0-AS)*TAN(FIMU)+AA1+AA2-AA3 dTRS1=dTPS1*(B0+(1.D0-B0)*SIG0(2)/SGMT)+TP*(1.D0-B0)/SGMT AA4=-UR*K2/SGMT*((1.D0-SIG0(2)/SGMT)**(K2-1.D0))*TAN(I0) AA4=AA4+((VM**2)*KNI)/((KNI*VM-SIG0(2))**2) c c calcul de la matrice Dt equivalente (2D) c KNN=1.D0/AA4 KNT=0.D0 KTN=KNN*dTRS1 KTT=0.D0 c c appel à CALC2 pour calculer les incréments de contraintes c CALL CALC2(KNN,KNT,KTN,KTT,SIG0,DU,DV,NSTRS, $ DELTA,RI0P,DCON) c c calcul de l' incrément de deformation plastique c DEP1=DU DEP2=0.D0 c END IF c ELSE IF(ABS(SIG0(2)).GE.ABS(SGMT))THEN c c cas avec sigma >= SGMT c c il n' y a pas de plasticité c c prédiction élastique c c definition des quantités necessaires au calcul c VUPU=0.D0 AS=1.D0 SR=S0-SIG0(2)*TAN(FI0) TP=SR TR=TP dTPS1=-TAN(FI0) AA1=((VM**2)*KNI)/((KNI*VM-SIG0(2))**2) c c calcul de la matrice Dt equivalente (2D) c KNN=1.D0/AA1 KNT=0.D0 KTN=UEQU/UP*KNN*dTPS1 KTT=TP/UP c c appel à CALC2 pour calculer les incréments de contraintes c CALL CALC2(KNN,KNT,KTN,KTT,SIG0,DU,DV,NSTRS, $ DELTA,RI0E,DCON) c c calcul de l' incrément de deformation plastique c DEP1=0.D0 DEP2=0.D0 c IF(UEQU.LT.UP)GO TO 200 c c il y a de la plasticité parfaite: trait horizontal de la courbe tau c c test pour vérifier si continuer en plasticité ou reprendre le trait c lineaire c EFFE=ABS(SIG0(1))-TP c IF(ABS(SIG0(1)).LT.TOL1)THEN DTEF=ABS(DCON(1)) c ELSE c DTEF=SIG0(1)*DCON(1)/ABS(SIG0(1)) c END IF c DTEF=DTEF-dTPS1*DCON(2) c IF((EFFE.LT.0.D0).OR.(DTEF.LT.0.D0))GO TO 200 c c definition des quantités necessaires au calcul c VUPU=0.D0 AS=1.D0 SR=S0-SIG0(2)*TAN(FI0) TP=SR TR=TP dTPS1=-TAN(FI0) AA1=((VM**2)*KNI)/((KNI*VM-SIG0(2))**2) c c calcul de la matrice Dt equivalente (2D) c KNN=1.D0/AA1 KNT=0.D0 KTN=KNN*dTPS1 KTT=0.D0 c c appel à CALC2 pour calculer les incréments de contraintes c CALL CALC2(KNN,KNT,KTN,KTT,SIG0,DU,DV,NSTRS, $ DELTA,RI0P,DCON) c c calcul de l' incrément de deformation plastique c DEP1=DU DEP2=0.D0 c END IF c 200 CONTINUE c c chargement de la matrice de rigidité au debut du pas c IF(KAPPA.EQ.1)THEN c RTT=KTT RTN=KTN c END IF c c chargement des valeurs à la fin des sub-pas c SIG0(1)=SIG0(1)+DCON(1) SIG0(2)=SIG0(2)+DCON(2) c EPS1=EPS1+DU EPS2=EPS2+DV c IF(SIG0(2).GE.-TOL2) THEN c c on a plasticité selon la composante normale: le joint se ouvre c ou il reste bloqué c SIG0(1)=0.D0 SIG0(2)=0.D0 c UEQU=UEQU c DEFP(1)=DEFP(1)+DU DEFP(2)=DEFP(2)+DV c EPIN0(1)=EPIN0(1)+DU EPIN0(2)=EPIN0(2)+DV c ELSE c c on n' a pas de plasticité selon la composante normale: le joint c reste fermé c UEQU=ABS(UEQU+DELTA) c DEFP(1)=DEFP(1)+DEP1 DEFP(2)=DEFP(2)+DEP2 c EPIN0(1)=EPIN0(1)+DEP1 EPIN0(2)=EPIN0(2)+DEP2 c DPP1=DEFP(1) DPP2=DEFP(2) c EPP1=EPIN0(1) EPP2=EPIN0(2) c END IF c 300 CONTINUE c c ********************************************************************** c ******** fin de la boucle du calcul et preparation des sorties ******* c ********************************************************************** c c test pour vérifier si le deux tau sont presque nulles c TOL3=1.D-1*ABS(RTT*(UEQU-VAR0(4))+RTN*DEPST(2)) c IF(ABS(SIG0(1)).LT.TOL3)THEN c SIG0(1)=0.D0 c END IF c c chargement des valeurs des contraintes finales c SIGF(1)=SIG0(1) SIGF(2)=SIG0(2) c c état du joint c IF(EPIN0(2).GT.EPSS)THEN c c le joint est ouvert c VARF(3)=1.D0 c ELSE IF((EPIN0(2).LE.EPSS).AND.(SIG0(2).GE.-TOL2))THEN c c le joint est fermé mais sans compression c DEFP(1)=DPP1 DEFP(2)=DPP2 c EPIN0(1)=EPP1 EPIN0(2)=EPP2 c VARF(3)=0.D0 c ELSE c c le joint est fermé avec compression c VARF(3)=2.D0 c END IF c c chargement de la variable EPSE c VARF(1)=SQRT(2.D0/3.D0*(EPIN0(1)**2+EPIN0(2)**2)) c c chargement de la variable EPOU c VARF(2)=DEFP(2) c c chargement de la variable UEQU c VARF(4)=UEQU c c chargement des déformations totales c VARF(5)=EPS1 VARF(6)=EPS2 c RETURN c END