amade3
C AMADE3 SOURCE CHAT 05/01/12 21:21:07 5004 C AMADE3.eso SOURCE bald MAR 28/03/95 00:00:00 & DEPST,IFOURB,XMAT,NMATT,IVAL,SIGF,DEFP,VARF,KERRE) C C----------------------------------------------------------------------- C C PLASTICITE MODELE 3D 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 VAR0(7) = DEFORMATION TOTALE: EPS3 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(3,3), RI0P(3,3), DCON(3) 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 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) EPS3=VAR0(7) c c calcul des petits pas c DU1=DEPST(1)/NN DU2=DEPST(2)/NN DV=DEPST(3)/NN c c incréments de déformation plastique finale c DEFP(1)=0.D0 DEFP(2)=0.D0 DEFP(3)=0.D0 c c contraintes finales c SIGF(1)=0.D0 SIGF(2)=0.D0 SIGF(3)=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 VARF(7)=0.D0 c c ********************************************************************** c ********** test pour verifier la condition initiale du joint ********* c ********************************************************************** c IF(EPIN0(3).GT.EPSS)THEN c c le joint est ouvert c c c il y a rapprochement des lèvres avec compression (1 cas) c R0=ABS(EPIN0(3)/DEPST(3)) c NN=INT(ABS((EPIN0(3)+DEPST(3))/DV)) c EPOU=-EPIN0(3) c EPIN0(1)=EPIN0(1)+R0*DEPST(1) EPIN0(2)=EPIN0(2)+R0*DEPST(2) EPIN0(3)=0.D0 c DEFP(1)=R0*DEPST(1) DEFP(2)=R0*DEPST(2) DEFP(3)=EPOU c EPS1=EPS1+R0*DEPST(1) EPS2=EPS2+R0*DEPST(2) EPS3=EPS3+R0*DEPST(3) c ELSE IF(DEPST(3).GT.(-EPIN0(3)+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 SIGF(3)=0.D0 c DEFP(1)=DEPST(1) DEFP(2)=DEPST(2) DEFP(3)=DEPST(3) c DEDE1=EPIN0(1)+DEPST(1) DEDE2=EPIN0(2)+DEPST(2) DEDE3=EPIN0(3)+DEPST(3) c VARF(1)=SQRT(2.D0/3.D0*(DEDE1**2+DEDE2**2+DEDE3**2)) c VARF(2)=DEPST(3) c VARF(3)=1.D0 c VARF(4)=VAR0(4) c VARF(5)=EPS1+DEPST(1) VARF(6)=EPS2+DEPST(2) VARF(7)=EPS3+DEPST(3) c RETURN c ELSE c c il y a fermeture sans compression (1 cas) c SIGF(1)=0.D0 SIGF(2)=0.D0 SIGF(3)=0.D0 c DEFP(1)=DEPST(1) DEFP(2)=DEPST(2) DEFP(3)=DEPST(3) c DEDE1=EPIN0(1)+DEPST(1) DEDE2=EPIN0(2)+DEPST(2) DEDE3=EPIN0(3)+DEPST(3) c VARF(1)=SQRT(2.D0/3.D0*(DEDE1**2+DEDE2**2+DEDE3**2)) c VARF(2)=DEPST(3) c VARF(3)=0.D0 c VARF(4)=VAR0(4) c VARF(5)=EPS1+DEPST(1) VARF(6)=EPS2+DEPST(2) VARF(7)=EPS3+DEPST(3) c RETURN c END IF c ELSE IF((EPIN0(3).LE.EPSS).AND.(SIG0(3).GE.-TOL2))THEN c c le joint est fermé mais sans compression c IF(DEPST(3).GT.EPSS)THEN c c les lèvres s' éloignent (1 cas) c SIGF(1)=0.D0 SIGF(2)=0.D0 SIGF(3)=0.D0 c DEFP(1)=DEPST(1) DEFP(2)=DEPST(2) DEFP(3)=DEPST(3) c DEDE1=EPIN0(1)+DEPST(1) DEDE2=EPIN0(2)+DEPST(2) DEDE3=EPIN0(3)+DEPST(3) c VARF(1)=SQRT(2.D0/3.D0*(DEDE1**2+DEDE2**2+DEDE3**2)) c VARF(2)=DEPST(3) c VARF(3)=1.D0 c VARF(4)=VAR0(4) c VARF(5)=EPS1+DEPST(1) VARF(6)=EPS2+DEPST(2) VARF(7)=EPS3+DEPST(3) 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 et sigma>=SGMT) c IF(ABS(SIG0(3)).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(3)/SGMT)**K2)*TAN(I0) AS=1.D0-((1.D0-SIG0(3)/SGMT)**K1) SR=S0-SIG0(3)*TAN(FI0) TP=-SIG0(3)*TAN(FIMU+ATAN(VUPU))*(1.D0-AS)+AS*SR TR=TP*(B0+(1.D0-B0)*SIG0(3)/SGMT) c AA1=(SIG0(3)*(1.D0-AS)*K2)/(SGMT*(COS(FIMU+ATAN(VUPU))**2 $ )) AA2=((1.D0-SIG0(3)/SGMT)**(K2-1.D0))*TAN(I0) AA3=1.D0+((1.D0-SIG0(3)/SGMT)**(2.D0*K2))*(TAN(I0)**2) AA4=K1*SIG0(3)/SGMT AA5=((1.D0-SIG0(3)/SGMT)**(K1-1.D0))*TAN(FIMU+ATAN(VUPU)) $ *AA4 AA6=SR/SGMT*K1*((1.D0-SIG0(3)/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(3)/SGMT)+TP*(1.D0-B0)/SGMT AA8=-UEQU*K2/SGMT*((1.D0-SIG0(3)/SGMT)**(K2-1.D0))*TAN(I0 $ ) AA8=AA8+((VM**2)*KNI)/((KNI*VM-SIG0(3))**2) c c calcul de la matrice Dt equivalente (2D) c KNN=1.D0/AA8 KNT=-((1.D0-SIG0(3)/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(3)/SGMT)**K1) SR=S0-SIG0(3)*TAN(FI0) TP=-SIG0(3)*TAN(FIMU)*(1.D0-AS)+AS*SR TR=TP*(B0+(1.D0-B0)*SIG0(3)/SGMT) c AA1=((1.D0-SIG0(3)/SGMT)**(K1-1.D0))*TAN(FIMU)*K1*SIG0(3) $ /SGMT AA2=SR/SGMT*K1*((1.D0-SIG0(3)/SGMT)**(K1-1.D0)) AA3=AS*TAN(FI0) dTPS1=-(1.D0-AS)*TAN(FIMU)+AA1+AA2-AA3 dTRS1=dTPS1*(B0+(1.D0-B0)*SIG0(3)/SGMT)+TP*(1.D0-B0)/SGMT AA4=-UR*K2/SGMT*((1.D0-SIG0(3)/SGMT)**(K2-1.D0))*TAN(I0) AA4=AA4+((VM**2)*KNI)/((KNI*VM-SIG0(3))**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 à CALC3 pour calculer les incréments de contraintes c $ DELTA,RI0E,DCON) c c calcul de l' incrément de deformation plastique c DEP1=0.D0 DEP2=0.D0 DEP3=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=SQRT(SIG0(1)**2+SIG0(2)**2)-(TR-TP)/(UR-UP)*(UEQU-UP $ )-TP c IF((ABS(SIG0(1)).LT.TOL1).AND.(ABS(SIG0(2)).LT.TOL1))THEN c DTEF=SQRT((DCON(1)**2)+(DCON(2)**2)) c ELSE c DTEF=(SIG0(1)*DCON(1)+SIG0(2)*DCON(2)) DTEF=DTEF/SQRT(SIG0(1)**2+SIG0(2)**2) c END IF c DTEF=DTEF-(TR-TP)/(UR-UP)*DELTA DTEF=DTEF-((UEQU-UP)/(UR-UP)*(dTRS1-dTPS1)+dTPS1)*DCON(3) 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(3)/SGMT)**K2)*TAN(I0) AS=1.D0-((1.D0-SIG0(3)/SGMT)**K1) SR=S0-SIG0(3)*TAN(FI0) TP=-SIG0(3)*TAN(FIMU+ATAN(VUPU))*(1.D0-AS)+AS*SR TR=TP*(B0+(1.D0-B0)*SIG0(3)/SGMT) c AA1=(SIG0(3)*(1.D0-AS)*K2)/(SGMT*(COS(FIMU+ATAN(VUPU))**2 $ )) AA2=((1.D0-SIG0(3)/SGMT)**(K2-1.D0))*TAN(I0) AA3=1.D0+((1.D0-SIG0(3)/SGMT)**(2.D0*K2))*(TAN(I0)**2) AA4=K1*SIG0(3)/SGMT AA5=((1.D0-SIG0(3)/SGMT)**(K1-1.D0))*TAN(FIMU+ATAN(VUPU)) $ *AA4 AA6=SR/SGMT*K1*((1.D0-SIG0(3)/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(3)/SGMT)+TP*(1.D0-B0)/SGMT AA8=-UEQU*K2/SGMT*((1.D0-SIG0(3)/SGMT)**(K2-1.D0))*TAN(I0 $ ) AA8=AA8+((VM**2)*KNI)/((KNI*VM-SIG0(3))**2) c c calcul de la matrice Dt equivalente (2D) c KNN=1.D0/AA8 KNT=-((1.D0-SIG0(3)/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 à CALC3 pour calculer les incréments de contraintes c $ DELTA,RI0P,DCON) c c appel à RISPL3 pour calculer les incréments de déformation plastique c c DEP1=DU1p DEP2=DU2p DEP3=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=SQRT(SIG0(1)**2+SIG0(2)**2)-TR c IF((ABS(SIG0(1)).LT.TOL1).AND.(ABS(SIG0(2)).LT.TOL1))THEN c DTEF=SQRT((DCON(1)**2)+(DCON(2)**2)) c ELSE c DTEF=(SIG0(1)*DCON(1)+SIG0(2)*DCON(2)) DTEF=DTEF/SQRT(SIG0(1)**2+SIG0(2)**2) c END IF c DTEF=DTEF-dTRS1*DCON(3) 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(3)/SGMT)**K1) SR=S0-SIG0(3)*TAN(FI0) TP=-SIG0(3)*TAN(FIMU)*(1.D0-AS)+AS*SR TR=TP*(B0+(1.D0-B0)*SIG0(3)/SGMT) c AA1=((1.D0-SIG0(3)/SGMT)**(K1-1.D0))*TAN(FIMU)*K1*SIG0(3) $ /SGMT AA2=SR/SGMT*K1*((1.D0-SIG0(3)/SGMT)**(K1-1.D0)) AA3=AS*TAN(FI0) dTPS1=-(1.D0-AS)*TAN(FIMU)+AA1+AA2-AA3 dTRS1=dTPS1*(B0+(1.D0-B0)*SIG0(3)/SGMT)+TP*(1.D0-B0)/SGMT AA4=-UR*K2/SGMT*((1.D0-SIG0(3)/SGMT)**(K2-1.D0))*TAN(I0) AA4=AA4+((VM**2)*KNI)/((KNI*VM-SIG0(3))**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 à CALC3 pour calculer les incréments de contraintes c $ DELTA,RI0P,DCON) c c calcul de l' incrément de deformation plastique c DEP1=DU1 DEP2=DU2 DEP3=0.D0 c END IF c ELSE IF(ABS(SIG0(3)).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(3)*TAN(FI0) TP=SR TR=TP dTPS1=-TAN(FI0) AA1=((VM**2)*KNI)/((KNI*VM-SIG0(3))**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 à CALC3 pour calculer les incréments de contraintes c $ DELTA,RI0E,DCON) c c calcul de l' incrément de deformation plastique c DEP1=0.D0 DEP2=0.D0 DEP3=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=SQRT(SIG0(1)**2+SIG0(2)**2)-TP c IF((ABS(SIG0(1)).LT.TOL1).AND.(ABS(SIG0(2)).LT.TOL1))THEN DTEF=SQRT((DCON(1)**2)+(DCON(2)**2)) c ELSE c DTEF=(SIG0(1)*DCON(1)+SIG0(2)*DCON(2)) DTEF=DTEF/SQRT(SIG0(1)**2+SIG0(2)**2) c END IF c DTEF=DTEF-dTPS1*DCON(3) 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(3)*TAN(FI0) TP=SR TR=TP dTPS1=-TAN(FI0) AA1=((VM**2)*KNI)/((KNI*VM-SIG0(3))**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 à CALC3 pour calculer les incréments de contraintes c $ DELTA,RI0P,DCON) c c calcul de l' incrément de deformation plastique c DEP1=DU1 DEP2=DU2 DEP3=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) SIG0(3)=SIG0(3)+DCON(3) c EPS1=EPS1+DU1 EPS2=EPS2+DU2 EPS3=EPS3+DV c IF(SIG0(3).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 SIG0(3)=0.D0 c UEQU=UEQU c DEFP(1)=DEFP(1)+DU1 DEFP(2)=DEFP(2)+DU2 DEFP(3)=DEFP(3)+DV c EPIN0(1)=EPIN0(1)+DU1 EPIN0(2)=EPIN0(2)+DU2 EPIN0(3)=EPIN0(3)+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 DEFP(3)=DEFP(3)+DEP3 c EPIN0(1)=EPIN0(1)+DEP1 EPIN0(2)=EPIN0(2)+DEP2 EPIN0(3)=EPIN0(3)+DEP3 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(3)) c IF((ABS(SIG0(1)).LT.TOL3).AND.(ABS(SIG0(2)).LT.TOL3))THEN c SIG0(1)=0.D0 SIG0(2)=0.D0 c END IF c c chargement des valeurs des contraintes finales c SIGF(1)=SIG0(1) SIGF(2)=SIG0(2) SIGF(3)=SIG0(3) c c état du joint c IF(EPIN0(3).GT.EPSS)THEN c c le joint est ouvert c VARF(3)=1.D0 c ELSE IF((EPIN0(3).LE.EPSS).AND.(SIG0(3).GE.-TOL2))THEN c c le joint est fermé mais sans compression 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+EPIN0(3)**2)) c c chargement de la variable EPOU c VARF(2)=DEFP(3) 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 VARF(7)=EPS3 c RETURN c END
© Cast3M 2003 - Tous droits réservés.
Mentions légales