C CFLUN2 SOURCE KICH 09/09/11 21:15:55 6505 SUBROUTINE CFLUN2(wrk52,wrk53,wrk54,wrk2,wrk3, & IB,IGAU,NBPGAU,NBGMAT,NELMAT,IFOURB) * * modele fluage type Norton dep/dt = C sig^n t^m * traite sigf = sig0 + k (deps - dep) * on pourrait separer deviateur et terme spherique * IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC CCOPTIO -INC DECHE -INC SMLREEL -INC SMEVOLL * * * SEGMENT WRK2 REAL*8 TRAC(LTRAC) ENDSEGMENT * SEGMENT WRK3 REAL*8 WORK(LW),WORK2(LW2) ENDSEGMENT * dimension spri0(8),delep0(8),spri1(8),delep1(8), &DIV(8) * cas isotrope ip1 = 4 * youn0 = xmat0(1) sigy0 = xmat0(ip1+1) xc0 = xmat0(ip1+2) xn0 = xmat0(ip1+3) xm0 = xmat0(ip1+4) ips0 = int(xmat0(ip1+5)) ipe0 = int(xmat0(ip1+6)) x2mu0= xmat0(1)/(1.+xmat0(2)) if(ib.eq.1.and.igau.eq.1) then * write(6,*) 't0' , youn0, sigy0, xn0 ,xm0 ,gk0,pk0 endif youn1 = xmat(1) sigy1 = xmat(ip1+1) xc1 = xmat(ip1+2) xn1 = xmat(ip1+3) xm1 = xmat(ip1+4) ips1 = int(xmat(ip1+5)) ipe1 = int(xmat(ip1+6)) x2mu1= xmat(1)/(1.+xmat(2)) if(ib.eq.1.and.igau.eq.1) then * write(6,*) 't1' ,ips1,ipe1 endif * delt = tempf - temp0 if (delt.lt.0..or.tempf.lt.0.) then moterr(1:16) = conm moterr(17:24) = 'CFLUN2-5' call erreur(943) return endif C---------CARACTERISTIQUES GEOMETRIQUES--------------------------------- C C COQUES C ALFAH=1.D0 IF(MFR.EQ.3.OR.MFR.EQ.5.OR.MFR.EQ.9) THEN EP1=work(1) IF(MFR.NE.5) ALFAH=work(2)**2 ENDIF C---------COQUES AVEC CT------------------------------------------------ C ON NE TRAVAILLE QUE SUR LES 6 PREMIERES COMPOSANTES IF(MFR.EQ.9) THEN NCONT=6 ELSE NCONT=NSTRS ENDIF * calcul increments de contrainte * remarque on utilise les caracteristiques elastiques a la date finale CALL CALSIG(depst,DDAUX,NSTAUX,CMATE,VALMAT,VALCAR, & N2EL,N2PTEL,MFR,IFOURB,IB,IGAU,EPAIST, & NBPGAU,MELE,NPINT,NBGMAT,NELMAT,SECT,LHOOK,TXR, & XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI,dsigt,IRTD) IF(IRTD.NE.1) THEN KERRE=69 GOTO 1010 ENDIF * determine direction et sens de sigf et eflu DO I=1,NSTRS DSIGT(I)=SIG0(I) + dsigt(i) ENDDO C---------CAS DES POUTRES----------------------------------------------- IF(MFR.EQ.7) THEN DIV(1)=1.D0/work(4) DIV(2)=1.D0 DIV(3)=1.D0 DIV(4)=work(10)/work(1) DIV(5)=work(11)/work(2) DIV(6)=work(12)/work(3) IF(DIV(4).EQ.0.D0) DIV(4)=1.D-10/SQRT(work(1)*work(4)) IF(DIV(5).EQ.0.D0) DIV(5)=1.D-10/SQRT(work(2)*work(4)) IF(DIV(6).EQ.0.D0) DIV(6)=1.D-10/SQRT(work(3)*work(4)) DO I=1,NCONT DSIGT(I)= DSIGT(I)*DIV(I) ENDDO ENDIF * * raisonne en deviateur trsig0 = (dsigt(1) + dsigt(2) + dsigt(3)) * 0.33333333333d0 spri0(1) = dsigt(1) - trsig0 spri0(2) = dsigt(2) - trsig0 spri0(3) = dsigt(3) - trsig0 do is = 4,nstrs spri0(is) = dsigt(is) enddo C----------------------------------------------------------------------- C CALCUL DE LA CONTRAINTE EQUIVALENTE SEQ C----------------------------------------------------------------------- seqtot=VONMIS0(spri0,NSTRS,MFR,IFOURB,EP1,ALFAH) if (seqtot - sigy1.gt.0.) then seq0 = seqtot else * pas de termes inelastiques do ic = 1,nstrs sigf(ic) = dsigt(ic) enddo varf(1) = var0(1) varf(2) = var0(2) goto 1002 return endif if (xn1.ge.0..and.xc1.ge.0..and.xm1.ge.0..AND.x2mu1.ge.0.) then else moterr(1:16) = conm moterr(17:24) = 'CFLUEN-1' c write(6,*) xn1,xc1,xm1,x2mu1 call erreur(943) return endif * point fixe pour determiner le multiplicateur de sigtot/seqtot icaz = 1 do ipfx = 1,50 if(ib.eq.1.and.igau.eq.1) then c write(6,*) 'pt fixe' , ipfx, seq0,seqtot,icaz endif goto (70,80) icaz 70 continue * fonction delr0 = (seq0 ** xn1) * (tempf ** xm1) xmult = delr0 * xc1 * delt * x2mu1 seq01 = seqtot - xmult goto 90 * fonction associee 80 continue xmult = seqtot - seq0 delr0 = xmult / delt/xc1 / (tempf ** xm1) / x2mu1 seq01 = delr0 ** (1/xn1) goto 90 90 continue if(ib.eq.1.and.igau.eq.1) then * write(6,*) 'delr0' , delr0,xmult,seq01, icaz endif if (ipfx.eq.1) then if (seq01.lt.0) then icaz = 2 seq01 = seqtot/2. else if (seq01.eq.0.) then icaz =1 seq01 = seqtot/2. endif endif varseq = abs((seq01 - seq0) / seq0) if (ib.eq.1.and.igau.eq.1) then c write(6,*) 'variation relative', varseq endif seq0 = seq01 if (seq0 .lt.0.) then c write(6,*) 'erreur point fixe', seqtot, seq0 moterr(1:16) = conm moterr(17:24) = 'CFLUN2-6' call erreur(943) return endif if (varseq.lt.1.e-6) goto 100 enddo 100 if (seqtot - seq0 .lt.0.) then c write(6,*) 'erreur point fixe', seqtot, seq0 moterr(1:16) = conm moterr(17:24) = 'CFLUN2-7' call erreur(943) return endif c au final 1000 continue do ic = 1,nstrs sigf(ic) = dsigt(ic) * seq0 / seqtot enddo varf(1) = var0(1) + xc1*(seq0**xn1)*(tempf**xm1)*delt varf(2) = var0(2) * position vis a vis des abaques if (ips1.gt.0) then temcri =0.d0 mevoll = ips1 segact mevoll kevoll = ievoll(1) segact kevoll mlree1 = iprogx mlree2 = iprogy segact mlree1,mlree2 nds = mlree2.prog(/1) * suppose valeurs de contraintes decroissantes et temps croissants do jds=1,nds-1 if (mlree2.prog(jds).le.mlree2.prog(jds+1).or. &mlree1.prog(jds).ge.mlree1.prog(jds+1)) then moterr(1:16) = conm moterr(17:24) = 'CFLUN2-8' call erreur(943) return endif if(mlree2.prog(jds).ge.seq0.and.seq0.gt.mlree2.prog(jds+1))then tosig = (mlree2.prog(jds) - seq0)/ & (mlree2.prog(jds) - mlree2.prog(jds+1)) * interpole logarithmiquement utemp = tosig * (log(mlree1.prog(jds+1)) & - log(mlree1.prog(jds))) + log(mlree1.prog(jds)) temcri = exp(utemp) * goto 1001 endif enddo 1001 if (temcri.gt.0) then varf(2) = var0(2) + delt/temcri if (varf(2).gt.1) then write(6,*) 'detruire prochaine etape', ipmail, conm,ib,igau endif endif endif 1002 continue if(ib.eq.1.and.igau.eq.1) then c write(6,*) 't0' ,sigf(3),varf(1),varf(2),depst(3),dsigt(3) endif 1010 continue RETURN END