C FLDO3D SOURCE FD218221 24/02/07 21:15:11 11834 subroutine fldo3d(XMAT,NMAT,sig0,sigf,deps,nstrs,VAR0,VARF, # NVARI,teta1,teta2,dt,ierr1,iso,mfr,ifou,istep,epst0,epstf, # tetaref,NBELAS3D,NB_MAT_BASE,NB_MAT_SUPP,NB_VAR_BASE,NB_VAR_SUPP, # NB_DECAL,NBNMAX3D,NBNB3D,IDIMB3D,XE3D,NBVIA3D,INLVIA3D, # NMAT0,NMAT1,NMAT2,NMAT3,NVAR0,NVAR1,NVAR2,NVAR3,NB_PARA_FIBRE_U) c A Sellier oct 2022 c NMAT0 fin des parametres elastiques c NMAT1 fin du modele de la matrice c NMAT2 fin des parametres pour les modules complementaires de fluendo (exple rag Multon ou fibre Gontero) c NMAT3 fin des parametres pour les renforts repartis c NMAT fin des parametres pour les champs de phase via Helmholtz c NVAR0 0 c NVAR1 fin des variables internes pour la matrice c NVAR2 fin des variables internes pour les modules complementaires c NVAR3 fin des variables internes pour les renforts repartis c NVARI fin des variables internes pour les champs de phases via Helmholtz c calcul de l ecoulement visco-elasto-plastique: Sellier mars 2014 c istep pour la phase non locale : juillet 2015 c actualisation DTT : aout 2018 c 08/05/2020 : hypothese non locale c istep=0 : calcul local c istep=1 : preparation du 1er passge non local c istep=2 : dernier passage non local effectue c istep=3 : passage non local intermediaire en non local iteratif c tables de dimension fixe pour resolution des sytemes lineaires implicit real*8 (a-h,o-z) implicit integer (i-n) c******** declaration des variables externes *************************** integer nmat,nstrs,NVARI,nbelas3d,ierr1,mfr real*8 xmat(nmat),sig0(nstrs),sigf(nstrs),deps(nstrs) c integer IVALMAT3D(nmat) real*8 epst0(nstrs),epstf(nstrs) real*8 var0(NVARI),varf(NVARI) real*8 dt,teta1,teta2,tetaref c variable logique pour matrice init isotrope logical iso integer ifou,istep c *** nombre de parametres materiaux ******************************* integer NB_MAT_BASE,NB_MAT_SUPP c *** nombre de variables internes ********************************* c nombre de variables internes integer NB_VAR_BASE,NB_VAR_SUPP c *** nombre de variables facultative entre les obligatoires et les helmholtz integer NB_DECAL c**************Coordonnes des noeuds de l element fini ***************** c ligne a inclure dans les modeles utilisant les noeuds integer NBNMAX3D,NBNB3D,IDIMB3D real*8 XE3D(3,NBNMAX3D) c*********************************************************************** c************* parametres pour les renforts **************************** c lignes a inclure dans les modeles utilisants les renforts c include './nombre_renforts.h' -INC HNBRREN c include './declare_materiaux_renforts.h' -INC HDECREN c*********************************************************************** c************** parametres pour Helmholtz ****************************** c lignes a inclure dans les modeles utilisant Helmholtz c include './nombre_helmholtz.h' -INC HNBRHEL c include './declare_materiaux_helmholtz.h' -INC HDECHEL c*********************************************************************** c*********************************************************************** c parametres obligatoires pour pointer correctement les renforts et c les formulations d helmholtz integer NMAT0,NMAT1,NMAT2,NVAR0,NVAR1,NVAR2,NVAR3 c*********************************************************************** c *** tableau pour la resolution des systemes lineaires ************ integer ngf parameter (ngf=65) real*8 XN(ngf),BN(ngf),ANN(ngf,(ngf+1)) integer ipzero(ngf) c *** parametres materiaux ***************************************** c parametres materiaux donnes real*8 hyd0,hyds real*8 young00,nu00,rt00,reft00,delta00,beta00,ept00 real*8 gft00,epsm00,psik,xflu,tauk00,taum00,taum0,nrjm,dt80,nrjw real*8 biotw00,xnsat00,poro2,vrag00,ekdg,gfr00,kshr real*8 vw2,mvgn,epc0,ekdc real*8 tetar,tetas,dfmx real*8 taar,nrjg,srsrag,alatr,epssr real*8 mdtt00,wdtt,pdtt real*8 vref,vmax,cvrt00 real*8 tdef,nrjp,srsdef,vdef00,cna,ttdd,ttrp,nrjd,tfid,tdid real*8 exmd,exnd,cnab,cnak,ssad,ttdf,nrjf,nsul real*8 ekds,skdw,tkvg real*8 ttrg,tdyn real*8 cshr,kdef,cdef,vvdef real*8 hplt00,hplw00,hplw,hplg00,hpls00 real*8 tetarw c variables pour la deformation thermique transitoire real*8 tdtt c parametres materiaux deduits des donnes c deformation elastique de reference pour le fluage real*8 epser,epc1,epc00,eweb c **** variables locales a fldo3d ******************************* c tenseur unitaires de reference real*8 vref33(3,3),vref33t(3,3) c variables generales integer i,j,k,l integer ncal,nmax,nmax1 c nmax nombre maxi de sous iteration en retour radial c nmax1 nombre maxi de sous iteration sans couplage des criteres diffus localises parameter(nmax=200,nmax1=190) integer nc c nombre de criteres parameter (nc=13) real*8 alphadp(3),betadp(3),factif(nc) logical actif(nc),complet(nc) real*8 vnorm logical ppas real*8 deps6(6),epstf6(6),epst06(6),epstf33(3,3),epst033(3,3) real*8 theta logical plast_seule,fl3d,end3d real*8 xwl2,tau_trac_wl2,xwl2max,cwrt,tau_trac_wl2_max real*8 sig13(3) real*8 young,nu,rt,rc,reft,rc1,delta,beta,gft,gfr real*8 ept,epsm,xnsat,biotw,Vvgel,kgel,cgel real*8 lambda,Mub real*8 raideur66(6,6),souplesse66(6,6) real*8 xmt logical dtiso real*8 nc33(3,3) real*8 dth00,dth0,poro1,vw1,dth1,CTHp1,CTHv1 real*8 epsk006(6),epsm006(6),sigef06(6) real*8 sigke06(6),sigm06(6) real*8 dtr00,dtr0,dfl00,dfl0 real*8 bw0,pw0,bg0,pg0,bs0,ps0 real*8 epspc600(6),epspc60(6),epspc6p(6),depspc6(6) real*8 epspt600(6),epspt60(6),epspt6p(6),depspt6(6),depspt6p(6) real*8 epsvt600(6),epsvt60(6),depsvt6(6) real*8 epspg600(6),epspg60(6),epspg6p(6),depspg6(6) real*8 epsps600(6),epsps60(6),epsps6p(6),depsps6(6) real*8 epspc6(6),epspt6(6),epsvt6(6),epspt33p(3,3) real*8 epspt33(3,3),epspt3(3),vepspt33(3,3),vepspt33t(3,3) real*8 epspg6(6),epsps6(6) real*8 pw,bw,srw2,CWv2 real*8 epse06(6) real*8 tauk2,taum2 real*8 dfin00,CMp00 real*8 cc03(3) real*8 cke66(6,6),ckk66(6,6),ckm66(6,6) real*8 xkm real*8 depse6(6),epse16(6) real*8 depsk6(6),epsk16(6) real*8 depsm6(6),epsm16(6) real*8 sig16(6),sigke16(6) real*8 At,St,M1,E1,vdef0,E2,AtF,StF,M1F,E1F,vdefF,E2F integer err1 real*8 hydr real*8 sig133(3,3),vsig133(3,3),vsig133t(3,3),sig16p(6) real*8 rts00,rts,rtw00,rtw,refw00,refw,rtg00,rtg real*8 hplt,hplg,hpls,mgel,mdef,bdef,dphidef real*8 dth2,CTHp2,CTHv2,Cwp2,CTHp,CTHv real*8 srw1,CWp1,CWv1,CWv,CWp real*8 tauk,taum real*8 aar0,aar1,def0,def1,bg1,pg1,bs1,ps1 real*8 epsk06(6),epsm06(6) integer na real*8 depleqc,epleqc0,epleqc1 real*8 dflu1 real*8 umdt,sigp,dti,dtr real*8 dt3(3),dgt3(3),dgc3(3),dwt3(3),dwc3(3),dst3(3),dsc3(3) real*8 sigf6(6),sigmf6(6) real*8 sigrm33(3,3),sigrf33(3,3),sigrf33p(3,3),sigrf6p(6) real*8 sigrm33p(3,3),sigrm6p(6),sigrm6(6) real*8 sigrh6p(6),sigrh33p(3,3),sigrh33(3,3),sigrh6(6),sigf6d(6) real*8 mdtt logical errr real*8 zero6(6),zero3(3) real*8 wplt06(6),wpltx06(6),wplt6(6),wpltx6(6) real*8 wpl3(3),vwpl33(3,3),vwpl33t(3,3) real*8 wplx3(3),vwplx33(3,3),vwplx33t(3,3) real*8 dr3(3),nc3(3) real*8 errgf1,errgf0 real*8 alphag,DCDW,alphas,rt_app,rc_app,tau_trac real*8 dtw0,dtg0,dts0,dcw0,dcg0,dcs0,kwrt,kwrc,href,gft_app real*8 dtl00,dtl0,dtl1 real*8 epspmf33(3,3),epspmf6(6) real*8 bw1,pw1,mwat real*8 Kb c eau dans les CSH et coeff de DTT real*8 wcsh0,wcshf,CDTT c directions principales des increments de deformations real*8 deps33(3,3),deps3(3),vdeps33(3,3),vdeps33t(3,3) real*8 taumdtt real*8 coeff_grand,coeff_petit parameter (coeff_grand=1.0d20,coeff_petit=1.0d-20) c poro meca / pressions real*8 bgel,dphigel,dphiw c increment des contraintes effectives real*8 dsigef6(6) c pression w csh real*8 pwcsh0,pwcshf,dpg,dps,dpw c activite de la refermeture localisee logical refermtl c endo pre pic compression real*8 dcpp,epc,xmc,dcpp00,dcpp0 logical dciso c traiement anti battement logical actif0(nc),actif1(nc) c resistances effectives et endos aux pics real*8 rceff,epleqc00,rteff,dpict,dpicc,rteff_app,rceff_app c taux de franchissement des seuils par les contraintes capillaires real*8 fshr6(6),fshr33(3,3),fshr33p(3,3),fshr006(6),fshr06(6) c caracteristique ecoulement initial (liquide) real*8 taum_ja,youn_ja,nu_ja,ss_ja c contrainte d endo hydrique real*8 skdw00,epse_tild6(6) c exposant de couplage pour cwp real*8 T1,T2,Tref c coeff de dilatation thermique des renforts real*8 dalpr c coeff drucker prager jeune age real*8 dlja c evolution jeune age Kelvin real*8 tauk_ja,tauk0,psi_ja,psi00 c resistances jeune age real*8 rcja,rtja c couplage des endommagements traction compression real*8 altc,dct3(3),dct0,dcc3(3),dcc0 c tau de cisaillement de Drucker Pragger real*8 taudp,taudpmax c indicateur de battement sur la resolution des criteres logical battcr3d c autorisation de resolution couplee en multi criteres logical CPTCR(nc,nc) c numeros de la varible non locale a traiter integer nl,ivari c controle de l hypothese d ecoulement logical garder,travail_fissure,endo_interface_explicit parameter (endo_interface_explicit=.true.) c calcul des deformations non locales real*8 eps_nl06(6),w_nl06(6),w_nlx06(6) real*8 eps_nlf6(6),w_nlf6(6),w_nlxf6(6) real*8 deps_nl3(3),deps_nl6(6),deps_l6(6) real*8 deps_nl33(3,3),source_nl_pt3(3),source_nl_pt6(6) real*8 v_source_nl_pt33(3,3) real*8 source_nl_pt33(3,3),Depspl_nl6(6) real*8 eps_nl6(6),eps_nl33(3,3),eps_nl3(3) real*8 veps_nl33t(3,3),veps_nl33(3,3) real*8 eptmasq033(3,3),eptmasq33(3,3),eptmasq33p(3,3) real*8 eptmasq06(6),eptmasq6(6) real*8 alpha_nl real*8 xx33(3,3) c logique mise au point plastc logical log_plastc,log_plastt c romain declarations pour les fibres c include './declare_fibres.h' -INC HDECFIB c variables supplementaires pour le traitement non local des increment de ept00 real*8 depspt33(3,3),depspt3(3),vdepspt33(3,3),vdepspt33t(3,3) c variables pour les traitements non locaux specifiques a fldo3d logical log_H_TWL2 integer Num_H_TWL2 logical log_H_DNL(6) integer Num_H_DNL(6) c autres variables pour non local logical ENDO_NL integer IVARNL c longueur caracteristique deduite de la dissusion et capaci te d Helm holtz (cf. charge_materiaux_helmholtz.h) real*8 LcH logical Method_H,Method_N c*********************************************************************** c initialisation des parametres c nombre de renforts actifs NRENF00=int(max(xmat(NMAT1),0.d0)) c ** indicateur de premier passage if (var0(1).ne.1.) then ppas=.true. else ppas=.false. end if c ****************************************************************** c initialisation d une matrice identité do i=1,3 do j=1,3 if(i.eq.j) then vref33(i,j)=1.d0 vref33t(i,j)=1.d0 else vref33(i,j)=0.d0 vref33t(i,j)=0.d0 end if end do end do c initialisation de la variable de controle d erreur ierr1=0 c coeff pour theta methode theta=0.5d0 c pour ne pas avoir de fluage mettre plast_seule à vrai plast_seule=.false. c pour fluage fl3d=.true. c pour endo end3d=.true. c detecteur de battement lors du retour plastique battcr3d=.false. c *********** chargement des parametres materiaux depuis xmat ****** c l hydratation est considéree en fin de pas pour ne pas avoir c a recuperer celle du pas precedent les exposants de De Shutter c sont dans hydramat, il faur recompiler pour les modifier C do i=1,nmat C print*,'xmat(',i,')=',xmat(i) C end do C read* c Module d young young00=xmat(nbelas3d+88) if( young00.eq.0.) then c .or.(IVALMAT3D(nbelas3d+88).eq.0)) c Print*,'Donner Young YORF pour l hydratation de reference HREF' young00=xmat(1) c print*,'E(href)=',young00 c err1=1 end if c coefficient de Poisson nu00=xmat(nbelas3d+89) if( nu00.eq.0.) then c .or.(IVALMAT3D(nbelas3d+89).eq.0)) then c Print*,'Donner coefficient de Poisson NURF pour l hydratation' c print*,'de reference HREF' c err1=1 nu00=xmat(2) c print*,'nu(href)=',nu00 end if c print*,'alpha fluendo3d=',alpha c hydratation hydr=xmat(nbelas3d+1) c hydratation seuil hyds=xmat(nbelas3d+2) c resistances et pression limite rt00=xmat(nbelas3d+3) if(rt00.eq.0.) then print*,'Il manque la resistance a la traction RT' err1=1 end if c seuil pour la refermeture des fissures reft00=xmat(nbelas3d+4) if(reft00.lt.rt00/100.d0) then print*,'La contrainte de refermeture REF est trop faible' err1=1 end if c resistance en compression-par cisaillement rc00=xmat(nbelas3d+5) if(rc00.eq.0.) then print*,'Il manque la resistance en compression RC ' err1=1 end if c coeff drucker prager delta00=xmat(nbelas3d+6) if(delta00.eq.0.) then print*,'Il manque DELT pour Drucker Prager ' err1=1 end if c coeff dilatance de cisaillement beta00=xmat(nbelas3d+7) if (beta00.eq.0.) then print*,'Il manque BETA pour Drucker Prager ' err1=1 end if c deformation au pic de traction ept00=xmat(nbelas3d+8) if (ept00.eq.0.) then ept00=rt00/young00 end if c taux d ecrouissage des phases effectives / gel hplg00=xmat(nbelas3d+9) if (hplg00.eq.0.) then hplg00=0.03*young00 end if c coeff de concentration de contrainte capillaire pour le critere hydrique cshr=xmat(nbelas3d+10) c compressibilite du gel kgel=xmat(nbelas3d+11) c energie de fissuration en traction directe gft00=xmat(nbelas3d+12) if (gft00.eq.0.) then print*,'Il manque l energie de fissuration GFT ' err1=1 end if c deformation caracteristique du potentiel de fluage epsm00=xmat(nbelas3d+13) if (epsm00.eq.0.) then print*,'Il manque le potentiel de fluage EPSM ' err1=1 end if c raideur relative Kelvin / Young psi00=xmat(nbelas3d+14) if (psi00.eq.0.) then print*,'Il manque le multiplicateur PSIK pour Kelvin ' err1=1 end if c coeff amplicateur à 2/3 rc pour le fluage xflu=xmat(nbelas3d+15) if (xflu.eq.0.) then xflu=1.0 end if c temps caracteristique pour Maxwell taum00=xmat(nbelas3d+17) if(taum00.eq.0.) then print*,'Il manque le temps caracterisqtique de fluage TAUM' err1=1 end if c temps caracteristique pour Kelvin tauk00=xmat(nbelas3d+16) if(tauk00.eq.0.) then print*,'Il manque le temps caracterisqtique de fluage TAUK' err1=1 end if c nrj activation du potentiel de fluage nrjm=xmat(nbelas3d+18) c endommagement thermique a 80°C dt80=xmat(nbelas3d+19) c porosite poro2=xmat(nbelas3d+22) c coeff de biot pour l eau biotw00=xmat(nbelas3d+20) c module de biot pour le non sature xnsat00=xmat(nbelas3d+21) c volume maximal de rgi vrag00=xmat(nbelas3d+23) c volume d eau dans les pores vw2=xmat(nbelas3d+24) if(vw2.eq.0.) then vw2=poro2 end if c coeff de couplage endo hydrique endo de compression DCDW=xmat(nbelas3d+25) c exposant de Van Genuchten pour la pression capillaire mvgn=xmat(nbelas3d+26) if(mvgn.eq.0.) then mvgn=0.5d0 end if c deformation au pic de compression epc0=xmat(nbelas3d+27) c deformation cracateristique endo de compression ekdc=xmat(nbelas3d+28) if(ekdc.eq.0.) then ekdc=1.d0 end if c deformation caracteristique pour l endo de rgi ekdg=xmat(nbelas3d+29) if(ekdg.eq.0.) then ekdg=0.003d0 end if c energie de fissuration pour l endo de traction gfr00=xmat(nbelas3d+30) if(gfr00.eq.0.) then print*,'Il manque l energie de refermeture de fissure GFR' err1=1 end if c volume des vides accessibles pour les RGI vvgel=xmat(nbelas3d+31) c coeff de concentration de contrainte des RGI cgel=xmat(nbelas3d+32) if(cgel.eq.0.) then c par défaut on prend la solution de la sphère élastique cgel=2.0d0 end if c temperature de reference des parametres materiaux (celsius) tetar=xmat(nbelas3d+33) c temperature seuil pour l endo thermique (celsius) tetas=xmat(nbelas3d+34) c endommagement maximum par fluage dfmx=xmat(nbelas3d+35) c temps caracterisqtique de la rag à tref taar=xmat(nbelas3d+36) c nrj d'activation de la rgi nrjg=xmat(nbelas3d+37) c seuil de saturation minimal pour avoir la rgi srsrag=xmat(nbelas3d+38) c constante de calage de la defom therm transitoire mdtt00=xmat(nbelas3d+39) c volume de reference pour Weibull vref=xmat(nbelas3d+40) c volume d integration maxi pour Weibull vmax=xmat(nbelas3d+41) if(vmax.eq.0.) then c on met VMAX = VREF pour emepcher l effet d echelle probabiliste VMAX=VREF end if c coefficient de variation de la resistance a la traction cvrt00=xmat(nbelas3d+42) c temps cracteristique pour la def tdef=xmat(nbelas3d+43) c energie d activation de precipitation de la def nrjp=xmat(nbelas3d+44) c seuil de saturation pour declancher la def srsdef=xmat(nbelas3d+45) c quantite maximale de def pouvant etre realise vdef00=xmat(nbelas3d+46) c teneur en alcalin pour la def cna=xmat(nbelas3d+47) c temperature de reference pour la precipitation de la def (Celsius) ttrp=xmat(nbelas3d+57) c energie d activation des processus de dissolution des phases primaires nrjd=xmat(nbelas3d+56) c temps caracteristique pour la fixation des aluminiums en temperature tfid=xmat(nbelas3d+55) c temps caracteristique pour la dissolution de l ettringite primaire tdid=xmat(nbelas3d+54) c temperature de reference pour la dissolution de la def (Celsius) ttdd=xmat(nbelas3d+53) c exposant de la loi de couplage precipitation def alcalins exmd=xmat(nbelas3d+52) c exposant de la loi de couplage temperature de dissolution alcalins exnd=xmat(nbelas3d+51) c concentration en alcalin de blocage de la def cnab=xmat(nbelas3d+50) c concentration caracteristique pour les lois de couplage de la def cnak=xmat(nbelas3d+49) c rapport molaire S03 / Al2O3 du ciment ssad=xmat(nbelas3d+48) c defcaracteristique pour l endo du a la def ekds=xmat(nbelas3d+58) c contrainte caracteristique pour l endo hydrique skdw00=xmat(nbelas3d+59) * pourquoi skdw ne pourrait pas etre nul? ** if(skdw00.eq.0.) then ** print*,'Il manque la contrainte SKDW pour l endo capillaire ' ** err1=1 ** end if c tetak pour exposant de van genuchten en fonction de teta tkvg=xmat(nbelas3d+60) c temperature minimale pour la fixation des alu en cas de def ttdf=xmat(nbelas3d+61) c energie activation pour la vitesse de fixation des alu en cas de def nrjf=xmat(nbelas3d+62) c nombre de moles de sulfates par unite de volume nsul=xmat(nbelas3d+63) c temps caracteristique pour l amplification dynamique tdyn=xmat(nbelas3d+64) c temperature caracteristique pour la rag (le gel) ttrg=xmat(nbelas3d+65) c module d ecrouissage pour la resistance a la pression de DEF hpls00=xmat(nbelas3d+66) c module de compressibilite pour la def kdef=xmat(nbelas3d+67) c volume des vides pour la def Vvdef=xmat(nbelas3d+68) c coeff de concentration de contrainte pour la def cdef=xmat(nbelas3d+69) c module d crouissage pour la resistance effective a la pression d eau hplw00=xmat(nbelas3d+70) c temperature de reference pour les coeffs de l isotherme tetarw=xmat(nbelas3d+71) c coeff de couplage endo de rag endo de compression alphag=xmat(nbelas3d+72) c coeff de couplage enddo de def endo de compression alphas=xmat(nbelas3d+73) c effet capillaire sur traction apparente kwrt=xmat(nbelas3d+74) c effet capillaire sur compression apparente kwrc=xmat(nbelas3d+75) c hydratation de reference pour la donnee des parametres mecaiques href=xmat(nbelas3d+76) if(href.eq.0.) then print*,'le degre d hydratation de reference HREF est nul ' err1=1 end if c temps caracteristique pour la deformation thermique transitoire tdtt=xmat(nbelas3d+77) c eau pour la DTT wdtt=xmat(nbelas3d+78) c pression d eau de reference dans les csh pdtt=xmat(nbelas3d+79) c contrainte seuil pour initier le fluage de maxwell ss_ja=xmat(nbelas3d+80) c temps caracteristique etat liquide taum_ja=xmat(nbelas3d+81) c young jeune age youn_ja=xmat(nbelas3d+82) c poisson jeune age nu_ja=xmat(nbelas3d+83) c delta drucker prager jeune age dlja=xmat(nbelas3d+84) c rc jeune age rcja=xmat(nbelas3d+85) c rt jeune age rtja=xmat(nbelas3d+86) c supplement de dilatation thermique des acier (alphar-alpha) dalpr=xmat(nbelas3d+87) c couplage traction cisaillement altc =xmat(nbelas3d+90) c avancement latent rag alatr=xmat(nbelas3d+91) c seuil de deformation pour endo de rag epssr=xmat(nbelas3d+92) c module d ecrouissage traction directe hplt00=0.d0 c nrj activation pour l eau nrjw=17000.d0 c amplitude reversible du fluage jeune age psi_ja=psi00*coeff_grand c temps de kelvin jeune age tauk_ja=taum_ja*coeff_grand c recuperation des parametres pour les fibres c include './charge_materiaux_fibres.h' -INC HMATFIB c ***********test de conformite sur les donnees pour la matrice **** if(err1.eq.1) then print*,'Incoherence dans les donnees de fluendo3d' ierr1=1 return end if c ** verif validite de la deformation au pic de traction ept00=max(rt00/young00,ept00) c ** verif coherence deformation pic de compression epc1=rc00/young00 c eps pic comp ne peut pas etre inferieur à rc/E epc00=max(epc0,epc1) c ** verif validite de la dilatance if(beta00.gt.sqrt(3.d0)) then print*,'Dilatance BETA excessive dans fluendo3d' ierr1=1 return end if c ** verif presence d un ecrouissage pour la rag if(hplg00.le.0.d0) then if(vrag00.gt.0.) then print*,'Taux d ecrouissage HRAG doit etre > 0' ierr1=1 return end if end if c ** verif presence d un ecrouissage pour la def if(hpls00.le.0.d0) then if((vdef00.gt.0.).or.(nsul.gt.0.)) then print*,'Taux d ecrouissage HDEF doit etre > 0' ierr1=1 return end if end if c ** test nbr maxi de criteres actifs pour le couplage if(nc.ne.13) then err1=1 print*,'pb nbr max de criteres ds fldo3d .ne. 13' return end if c *** initialisation des parametres calculees a partir des donnees * c ** stockage hydratation if (ppas) then c au 1er passage on suppose hyd0=hydr hyd0=hydr c on initialise var03d en cas de sous incrementation par fluage var0(68)=hyd0 else c on recupere l hydratation du pas precedent hyd0=var0(68) end if varf(68)=hydr c ** initialisation des variables internes associee au volume d eau c eau et porosite (initialise debut=fin 1er pas) if (ppas) then var0(69)=vw2 var0(70)=poro2 end if c stockage pour le prochain pas varf(69)=vw2 varf(70)=poro2 c tables a zero do j=1,6 zero6(j)=0.d0 end do do j=1,3 zero3(j)=0.d0 end do c****************** parametres pour Helmholtz ************************** c lignes a inclure dans les modeles utilisant des formulations d helmholtz c include './charge_materiaux_helmholtz.h' -INC HMATHEL do NL=1,NB_HELM c lecture des variable internes et rangements dans Helm0(NL,ii) c et Helm1(NL,ii) avec ii=1,NB_VARI_PAR_HELM c include './lire_vari_helmholtz.h' -INC HLVIHEL end do c*********************************************************************** c*************** parametres pour les renforts ************************** c lignes a inclure dans les modeles utilisant des renforts c include './charge_materiaux_renforts.h' -INC HMATREN c*********************************************************************** c********** traitement des autorisations d activation des criteres ***** do i=1,nc do j=1,nc c initialisation de la compatibilite des criteres a faux CPTCR(i,j)=.false. if(i.ne.j) then if ((i.le.9).and.(j.le.9)) then c tous les criteres de pression sont inter compatibles CPTCR(i,j)=.true. else if ((i.le.12).and.(j.le.12).and. # (i.ge.10).and.(j.ge.10)) then c les macro fissures sont compatibles entre elles CPTCR(i,j)=.true. else if((i.eq.13).and.(j.ge.10).and.(j.le.12)) then c DP et macros fissures incompatibles CPTCR(i,j)=.false. else c autres couplages interdits CPTCR(i,j)=.false. end if else c cas des criteres seuls CPTCR(i,i)=.true. end if c on desactive tous les criteres c CPTCR(i,j)=.false. end do end do c ****************************************************************** c affichage des options du calcul c print*,'fluage',fl3d c print*,'endo',end3d c*********************************************************************** c*********************************************************************** c chargement increment de deformation imposee et deformation finale c remarque 4-5-6 sont des gama jusqu ici, ensuite des epsilons if(nstrs.lt.6) then do i=1,nstrs deps6(i)=deps(i) epstf6(i)=epstf(i) epst06(i)=epst0(i) end do do i=nstrs+1,6 deps6(i)=0.d0 epstf6(i)=0.d0 epst06(i)=0.d0 end do else do i=1,6 deps6(i)=deps(i) epstf6(i)=epstf(i) epst06(i)=epst0(i) end do end if c passage gama en epsilon do i=4,6 deps6(i)=0.5d0*deps6(i) epstf6(i)=0.5d0*epstf6(i) epst06(i)=0.5d0*epst06(i) end do c*************** traitement non local de TWL2(VARI numero 72) ********** ivari_twl2=72 c on teste l existence d un champ de phase pour cette variable et c on recupere son numero Num_H_TWL2 call exhmtz(istep,NBVIA3D,INLVIA3D,NB_HELM, #ivari_twl2,log_H_TWL2,Num_H_twl2) c si le champ de phase existe if(log_H_TWL2) then c calcul approché de l'exposant de Weibull (m) if(cvrt00.ne.0.) then if (cvrt00.gt.1.) then print*,'Dans fluendo3d, erreur de donnees car le ' print*,'coeff de variation de Rt pour Weibull est >1 ' ierr1=1 return end if eweb=1.2d0/cvrt00-0.2d0 c print*,'eweb=',eweb else print*,'Donner le coeff de variation de Rt pour Weibull' print*,'Mettre Vmax a zero pour annuler l effet de WL2' print*,'ou bien Vmax a Vref pour utiliser WL2' ierr1=1 return end if c print*,'Etape 2-1 ds fldo3d Istep ',istep c istep=1 : 1er passage pour preparer le calcul non local if(istep.eq.1) then c avec cwrt du pas precedent et contraintes pas actuel if (ppas.or.(hydr.le.hyds)) then if ((vmax.ne.0.).and.(eweb.ne.0.).and.(vref.ne.0.)) then xweb=1.d0/eweb cwrt=(vref/vmax)**xweb else cwrt=1.d0 end if c 1er passage il faut initialiser cwrt var0(76)=cwrt else c on prend le dernier cwrt estime cwrt=var0(76) end if c on stocke cwrt dans varf pour les sous increments locaux varf(76)=cwrt c variable pour le stockage du max de xwl2 non local (mxwl ds idvar4) c (a mettre absolument a un car utilise avec cette hypothese dans unpas) varf(73)=1.d0 c variable pour le stockage du taux de charge maxi sur la structure (mtwl ds idvar4) c (a mettre absolument a un car utilise avec cette hypothese dans wl2d3d) varf(74)=1.d0 else c nouveau passage non local c *** recup variable non locales WL2 ************************* c recup de xwl2=//(s/r)**m .psi .dv=Veq.(max(s)/rt)**m xwl2=var0(ivari_twl2) c lecture de vari supplementaires associees a helmholtz ierr1=1 return if(ierr1.eq.1) then print*,'Pb calcul Helmholtz twl2 ds fldo3d' return else c lecture des variables internes complementaires pour Helmholtz NL=Num_H_TWL2 c IVARI=IVARI c include './lire_vari_helmholtz.h' -INC HLVIHEL c les var0 de cette formulation sont alors dans HELM0(Nl,i) end if c recup du maxi sur le modele de sigma/rt_app tau_trac_wl2=var0(72) c recup du maxi sur le modele xwl2 non local xwl2max=var0(73) c tau de charge maxi tau_trac_wl2_max=var0(74) c indicateur de localisation de la fissure c inutile de l affecter par l hydratation cela a ete fait lorsque istep=1 dtl0=var0(75) c endommagement localise probabble fin de pas (estime pour istep=1) dtl1=var0(119) c coeff de Weibull rt pas precedent cwrt=var0(76) c pression capillaire initiale pw0=var0(81) c saturation en eau * biot bw0=var0(80) c ** calcul de la nouvelle resistance a la traction ********** c contrainte et resistance ne sont pas utilisees a cette etape c on les met a zero do i=1,3 sig13(i)=0.d0 end do call wl2d3d(2,vmax,vref,sig13,eweb,0.d0,cwrt, # xwl2,tau_trac_wl2,dtl0,dtl1,xwl2max,tau_trac_wl2_max, # bw0,pw0,kwrt) c on stocke les variables non locales pour verif eventuelle varf(71)=xwl2 c variable non locale : sigm_maxi /rt varf(72)=tau_trac_wl2 c stockage de la valeur maxi de smax /rt varf(74)=tau_trac_wl2_max c endo fin de pas c on copie endommagement fin de pas utilisee pour cwrt pour le c debut de pas suivant car dtl calcule en fin de programme c n est pas celui utilise pour calculer cwrt c DTRA debut de pas varf(75)=dtl1 c DTL0 varf(119)=varf(75) c coeff d effet d echelle sur rt c print*,'cwrt ds fldo3d apres istep=2',cwrt varf(76)=cwrt c ecriture des variables internes HELM1 associees a la formulation NL de Helmholtz c include './ecrire_vari_helmholtz.h' -INC HEVIHEL end if else c pas de calcul non local pour wl2 c ** rt00 n est pas affectee par la methode WL2 ************** c istep=0 :le calcul sans effet d echelle Weibull cwrt=1.d0 c liste des variables non locales : XWL2 varf(71)=0.d0 c variable non locale : sigma maxi /rt TWL2 varf(72)=0.d0 c max de xwl2 non utilisse on y laisse 1. MXWL varf(73)=1.d0 c contrainte relative maxi non locale non utiliseee MTWL varf(74)=0.d0 c endommagement maxi de traction DTMX varf(75)=var0(75) c modif de rt par effet d echelle inactif CWRT varf(76)=cwrt end if c*********************************************************************** c influence de l hydratation sur les parametres materiau c les resistances locales sont prise egales à la resistance c macroscopique idem pour les contraintes de refermeture rts00=rt00 rtg00=rt00 rtw00=rt00 refw00=reft00 c pas d effet d hydratation sur coeff de def therm transitoire mdtt=mdtt00 c controle calcul raideur iso if(.not.iso) then print*,' cas imprevu dans fluendo3d ' print*,' matrice de raideur initiale non isotrope' ierr1=1 return end if c modification des parametres materiau en fonction de l hydratation call hydm3d(hyd0,hydr,hyds,young00,young,nu00,nu,rt00,rt, #rteff,reft00,reft,rc00,rc,rceff,delta00,delta,beta00,beta,gft00, #gft,ept00,ept,epsm00,epsm,xnsat00,xnsat,biotw00,biotw,iso,lambda, #Mub,Kb,rtg00,rtg,raideur66,souplesse66,xmt,dtiso,cwrt,err1,teta2, #href,tetar,tkvg,gfr00,gfr,rts00,rts,rtw00,rtw,hplt00,hplt,hplg00, #hplg,hpls00,hpls,epc00,epc,xmc,dciso,dpict,dpicc,taum_ja,youn_ja, #nu_ja,ss_ja,ssfl,taum00,taum0,skdw00,skdw,dlja,tauk_ja, #tauk00,tauk0,psi_ja,psi00,psik,rcja,rtja) c traitement de l erreur eventuelle if(err1.eq.1) then print*,'Pb lors de la mise a jour avec hydramat ds as3d' ierr1=1 return end if c print*,'istep',istep,'rt',rt c print*,'Etape 4 ds fldo3d apres hydramat' c*********************************************************************** c Nombre maxi de fissures localisees en fonction de la direction if(.not.ppas) then c chargement des longueurs des elements dans les directions des c renforts do i=1,NRENF00 lcr(i)=var0(NVAR2+(i-1)*NB_VARI_PAR_RENF+13) end do end if c calcul du nb max de fissures sur la zone d integration call nfin3d(ppas,lcr,nc33,NRENF00,vecr,NB_RENF,longr,lsr, # lfr,deqr,syr,taur,rhor,rt00,XE3D,NBNMAX3D,NBNB3D,IDIMB3D, # NB_HELM,TAILH,log_H_RENF,Num_H_RENF,err1) if(err1.eq.1) then print*,'Erreur lors de l appel a nfin3d dans fldo3d' print*,'pour istep=',istep ierr1=1 return end if c sauvegarde des longueurs des elements ds les directions des renforts do i=1,NRENF00 varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+13)=longr(i) if(longr(i).eq.0.) then print*,'Ds fldo3d longr',i,'=', longr(i) print*,' istep',istep,'ppas',ppas ierr1=1 return end if end do c print*,'Etape 5 ds fldo3d apres nfin3d' c*********************************************************************** c chargement des principales variables internes c (etat du squelette solide) do i=1,6 epsk006(i)=var0(i+7) epsm006(i)=var0(i+13) epsvt600(i)=var0(31+i) sigef06(i)=var0(i+49) sigke06(i)=var0(i+55) sigm06(i)=var0(i+61) end do c endommagement effectif par fluage dfl00=var0(78) c endommagement thermique dth00=var0(77) c endommagement iso prepic de traction (DTPP) dtr00=var0(79) c endo iso pre pic de cisaillement (DCPP) dcpp00=var0(122) c pression capillaire bw0=var0(80) pw0=var0(81) c pression gel bg0=var0(82) pg0=var0(83) c pression def bs0=var0(84) ps0=var0(85) c tenseurs des deformations plastiques et criteres hydrique do j=1,6 c deformations plastiques epspc600(j)=var0(19+j) epspt600(j)=var0(25+j) epspg600(j)=var0(37+j) epsps600(j)=var0(43+j) c critere pression capillaire fshr006(j)=var0(129+j) end do c def plastique equivalente Drucker Prager epleqc00=var0(94) c*********************************************************************** c influence du degre d hydratation sur les variables internes call hydv3d(ierr1,hyd0,hydr,hyds,dth00,dth0,dtr00,dtr0, # dfl00,dfl0,dcpp00,dcpp0,epspt600,epspt60,epspc600,epspc60, # epsvt600,epsvt60,epspg600,epspg60,epsps600,epsps60,epsk006, # epsk06,epsm006,epsm06,epleqc00,epleqc0,fshr006,fshr06) c print*,'Etape 6 ds fldo3d apres hydv3d' c print*,'deformation totale initiale' c call afic6(epspt60) c*********************************************************************** c increments non locaux des deformations plastiques de traction c (DPTi=VARI numero 140+i) ENDO_NL=.true. do i=1,6 c numero de la variable interne a traiter par Helmholtz (a partie de idvar4) ivarinl=140+i c test du traitement non local + assigantion logique log_H_DNL(i) c et recuperation du numero de la formulation non locale associee c dans Num_H_DNL(i) call exhmtz(istep,NBVIA3D,INLVIA3D,NB_HELM, # ivarinl,log_H_DNL(i),Num_H_DNL(i)) c on verifie que les 3 deformations plastiques principales de c traction soient traitees par la formulation d Helmholtz ENDO_NL=(ENDO_NL.and.log_H_DNL(i)) end do c print*,'fluendo ENDO_NL etape 1',ENDO_NL if(ENDO_NL) then c longueur non locale moyenne pour l endommagement LcH=1.d0 do i=1,6 LcH=LcH*max(LcarH(Num_H_DNL(i),1),LcarH(Num_H_DNL(i),2), # LcarH(Num_H_DNL(i),3)) end do LcH=LcH**(1./6.) c indicateur de la methode de calcul des dimensions Method_H=.true. Method_N=.false. C print*,'LcH',LcH,'Method_H',Method_H,'Method_N',Method_N C read* c on est dans la procedure incrementale non locale do i=1,6 c position de la variable non locale ivarinl=140+i c recuperation de la valeur non locale eps_nl6(i)=var0(ivarinl) end do else LcH=1.d0 Method_H=.false. Method_N=.true. end if c*********************************************************************** c influence de la temperature sur les parametres materiau et c actualisation de l endo thermique debut de pas poro1=var0(70) vw1=var0(69) call thrm3d(teta1,nrjm,tetas,tetar,DT80,dth0, # dth1,CTHp1,CTHv1,poro1,vw1,nrjw) c endommagement thermique en fin de pas call thrm3d(teta2,nrjm,tetas,tetar,DT80,dth1, # dth2,CTHp2,CTHv2,poro2,vw2,nrjw) varf(77)=dth1 c calcul des coeffs moyens pour le potentiel de fluage CTHp=0.5d0*(CTHp1+CTHp2) CTHv=0.5d0*(CTHv1+CTHv2) c*********************************************************************** c calcul des avancements physico-chimiques c*********************************************************************** c **** calcul de la pression capillaire due a l eau en fin de pas ** call watr3d(mfr,biotw,poro2,vw2,xnsat,mvgn,pw,bw,srw2,teta2, # tetarw,tkvg,pw0,Mwat,dphiw,bw0,bw1) varf(80)=bw1 c modif eventuelle des viscosites de Kelvin en fonction de srw2 if(poro1.ne.0.) then srw1=vw1/poro1 else srw1=1.d0 end if CWv1=1.d0/srw1 CWv2=1.d0/srw2 CWv=0.5d0*(CWv1+CWv2) c modif potentiel de fluage de Maxwell en fonction de srw Cwp2=srw2 CWp1=srw1 CWp=0.5d0*(CWp1+Cwp2) C c pas d effet de l eau sur le potentiel C Cwp2=1.d0 C CWp1=1.d0 C Cwp=1.d0 c *** cas de l eau des CSH **************************************** if(ppas) then c initialisation de la quantite d'eau dans les CSH a la valeur c de reference wcsh0=wdtt pwcsh0=0.d0 else wcsh0=var0(120) pwcsh0=var0(121) end if c prise en compte de la pression dans les csh pour la DTT c cf. F. Manzoni, T. Vidal, A. Sellier, X. Bourbon, and G. Camps, c “On the origins of transient thermal deformation of concrete,” c Cem. Concr. Compos., vol. 107, 2020, c https://doi.org/10.1016/j.cemconcomp.2019.103508 call wdtt3d(wcsh0,wcshf,tdtt,mdtt,wdtt,pdtt,CTHv, # teta1,teta2,pw,dt,CDTT,pwcsh0,pwcshf,tetarw,err1) if(err1.eq.1) then print*,'Pb lor du calcul de la DTT dans fldo3d' ierr1=1 return end if varf(120)=wcshf varf(121)=pwcshf varf(129)=cdtt c ** Modification eventuelle de la viscosite de Kelvin (T,Sr) ****** tauk=tauk0*CWv/CTHv c print*,'fldo3d tauk0',tauk0,cwv,cthv c ** Modification de la viscosite de Maxwell (T,Sr) **************** c pour Sr, deja integre dans CC via Cwp taum=taum0/CTHv c ** Modification de la viscosite de Maxwell (T,Sr) pour la DTT **** if((CDTT.gt.coeff_petit).and.(hydr.gt.hyds)) then c dtt active c print*,'Fluendo3d CDTT',CDTT taumdtt=taum/CDTT else c pas de dtt avent percolation taumdtt=taum*coeff_grand end if c print*,'ds fluando3d taum=',taum c print*,'ds fluando3d taumdtt=',taumdtt c **** avancement de la rag **************************************** aar0=var0(86) bg0=var0(82) call rag3d(taar,nrjg,ttrg,aar0,srw2,srsrag,teta1, # teta2,dt,aar1,vrag00,Kb,Kgel,bgel,mgel,dphigel,bg0,bg1, # vvgel,cgel,Rt,alatr,err1) if(err1.eq.1) then print*,'Pb lors du calcul de rag dans fldo3d' ierr1=1 return end if c print*,'fldo3d',dphigel,bgel,mgel c read* varf(86)=aar1 varf(82)=bg1 c **** avancement de la def *************************************** bs0=var0(84) E1=var0(87) M1=var0(88) At=var0(89) St=var0(90) E2=var0(91) def0=var0(92) vdef0=var0(93) if((vdef00.gt.0.).or.(nsul.gt.0.)) then call def3d(ppas,tdef,nrjd,def0,srw2,srsdef,teta1,teta2,dt, # vdef00,def1,CNa,nrjp,ttrp,tfid,ttdd,tdid,exmd,exnd, # cnab,cnak,ssad,At,St,M1,E1,vdef0,E2,AtF,StF,M1F,E1F,vdefF,E2F, # ttdf,nrjf,nsul,ierr1,Kb,Kdef,Vvdef,Cdef,Rt,mdef,bdef,dphidef, # bs0,bs1) if(ierr1.eq.1) then print*,'pb lors du calcul de la def dans fldo3d' return end if else bs1=bs0 E1F=E1 M1F=M1 AtF=At StF=St E2F=E2 def1=def0 vdefF=vdef0 bdef=0.d0 dphidef=0.d0 mdef=0.d0 end if varf(84)=bs1 varf(87)=E1F varf(88)=M1F varf(89)=AtF varf(90)=StF varf(91)=E2F varf(92)=def1 varf(93)=vdefF c ****************************************************************** c prise en compte de l effet de pression (coeffs kwrt et kwrc) c et de l effet d echelle probabiliste (coeff cwrt) sur la c resistance apparente de traction et en compression rt_app=(rteff-bw*pw*kwrt)*cwrt*(1.d0-dpict) rc_app=(rceff-bw*pw*kwrc)*(1.d0-dpicc) c deformation elastique de reference pour le fluage prise a 1/3 c de rc sat effective epser=(rc00/young00)/3.d0 c ****************************************************************** c reevaluation de la deformation compatible avec l etat c de contrainte debut de pas call hydc3d(souplesse66,sigke06,epsk06,psik, # epse06,sigef06,fl3d,plast_seule) c print*,'Etape 7 ds fldo3d apres hydc3d' c*********************************************************************** c tir visco elastique c*********************************************************************** c ** initialisation des deformations ******************************* do i=1,6 c increment deformation elastique epse16(i)=epse06(i) c increment deformation kelvin epsk16(i)=epsk06(i) c deformation maxwell epsm16(i)=epsm06(i) c deformation visqueuse transitoire epsvt6(i)=epsvt60(i) c deformations plastiques cisaillement epspc6(i)=epspc60(i) c deformations plastiques traction epspt6(i)=epspt60(i) c deformations plastiques gel epspg6(i)=epspg60(i) c deformations plastiques def epsps6(i)=epsps60(i) end do c ** directions principales des deformations imposees et matrice de passage associee call x6x33(deps6,deps33) call b3_v33(deps33,deps3,vdeps33) c do i=1,3 c print*,'fldo3d deps3:',i,deps3(i) c end do c construction matrice de passage inverse call traps1(vdeps33t,vdeps33,3) c ** effet du chargement initial sur le potentiel de fluage ******* call dffn3d(sigm06,delta,rc_app,rt_app,xflu, # ssfl,dfin00,CMp00,dfmx,hydr,hyds,err1) if(err1.eq.1)then print*,'Erreur lors du calcul de CMP dans dffn3d fldo3d' ierr1=1 return end if c print*,'fldo3d cmp=',cmp00 c ** effet de l endommagement capillaire sur le potentiel de fluage call ampl3d(souplesse66,sigm06,fshr06,skdw,DCDW,epse06,sigef06, # epse_tild6,refw) c ** actualisation des coeffs de consolidation avec chargement **** C cf. Alain Sellier, Stéphane Multon, Laurie Buffo-Lacarrière, Thierry Vidal, Xavier Bourbon, Guillaume Camps, C Concrete creep modelling for structural applications: non-linearity, multi-axiality, hydration, temperature and drying effects, C Cement and Concrete Research, Volume 79, 2016,Pages 301-315,ISSN 0008-8846, C https://doi.org/10.1016/j.cemconres.2015.10.001. call conso3d(xkm,epsm,epser,epsm06,epse_tild6,cc03, # vdeps33,CMp00,CTHp,CTHv,Cwp) c do i=1,3 c print*,'fldo3d cc3:',i,cc03(i) c end do c ** endommagement par fluage ************************************ call dflu3d(xkm,cc03,dfl0,dflu1,dfin00) c stockage endo de fluage varf(78)=dflu1 c ** initialisation des variables pour le tir visco elastique **** c pour le fluage do i=1,3 alphadp(i)=0.d0 betadp(i)=0.d0 end do c pour les criteres do i=1,nc actif(i)=.false. complet(i)=.true. factif(i)=0.d0 end do c ** initialisation du nombre de criteres actifs (pas de criteres actifs lors du tir) NA=0 c ** assemblage des matrices de couplage de fluage ds base fixe call coupl3d(ann,xn,bn,ngf,err1,deps3,vdeps33,vdeps33t,dt, # cc03,taum,taumdtt,tauk,psik,epse06,epsk06,raideur66,Mgel,bgel, # dphigel,Cgel,Hplg,Mdef,bdef,dphidef,Cdef,Hpls,Mwat,bw,dphiw, # cshr,alphadp,betadp,actif,factif,nc,dsigef6,depse6,depsk6, # depsm6,depspc6,depspt6,depsvt6,depspg6,depsps6,.true.,ipzero, # dpg,dps,dpw,zero6,.false.,complet,NA,CMp00) c ** traitement de l erreur eventuelle lors du tir if(err1.eq.1) then print*,'pb assemblage couplagf3d ds fldo3d' ierr1=1 return end if c print*,'Etape 8 ds fldo3d apres coupl3d' c ** recuperation des increments et actualisation des deformations do i=1,6 c increment deformation elastique epse16(i)=epse16(i)+depse6(i) c increment deformation kelvin epsk16(i)=epsk16(i)+depsk6(i) c deformation maxwell epsm16(i)=epsm16(i)+depsm6(i) c deformation visqueuse transitoire (dtt) epsvt6(i)=epsvt6(i)+depsvt6(i) end do c ******* recuperation deformation equivalente de compression ****** epleqc1=epleqc0 c ***** actualisation des pressions ******************************** c gel rag pg1=pg0+dpg c def ps1=ps0+dps c depression capillaire pw1=pw0+dpw c ***** etat des contraintes effectives apres tir visco elastique ** do i=1,6 sig16(i)=sigef06(i)+dsigef6(i) sigke16(i)=sigke06(i) do j=1,6 sigke16(i)=sigke16(i)+raideur66(i,j)*depsk6(j)*psik end do C print*,'Fluendo3d etape 8-2 de tir visco elastique' C print*,'sig16(',i,')=',sig16(i) C print*,'sigke16(',i,')=',sigke16(i) C print*,'deps6(',i,')=',deps6(i) C print*,'depse6(',i,')=',depse6(i) C print*,'depsk6(',i,')=',depsk6(i) C print*,'depsm6(',i,')=',depsm6(i) end do c*********************************************************************** c test des criteres de plasticité c*********************************************************************** c initialisation du compteur de sousiterations locales NCAL=0 c initialisation indicateur de plasticite en cisaillement log_plastc=.false. c ** increment des iteration locales 20 NCAL=NCAL+1 c print*,'******numero sous iter retour radial ds Fluendo3d:',ncal if(NCAL.eq.1) then c diagonalisation du vecteur des contraintes call x6x33(sig16,sig133) c print*,'fldo3d sig133' c call afic33(sig133) call b3_v33(sig133,sig13,vsig133) c construction matrice de passage inverse call traps1(vsig133t,vsig133,3) c print*,'etape 2-0-1 matrice de passage ds fldo3d' c call afic33(vsig133) else c les directions principales ne changent pas lors du retour call chrep6(sig16,vsig133,.false.,sig16p) do j=1,3 sig13(j)=sig16p(j) end do end if c ------------------------------------------------------------------ c passage de la loi de comportement dans la base principale c des contraintes if(.not.iso) then print*,'programmer changement base loi de comp ds fldo3d' ierr1=1 return end if c passage des resistances et contraintes de refermeture c dans la base principale des contraintes if(.not.iso) then print*,'programmer changement base resistances ds fldo3d' read* ierr1=1 end if c passage dans la base principale du tir visco elastique call chrep6(epspc6,vsig133,.false.,epspc6p) call chrep6(epspt6,vsig133,.false.,epspt6p) call chrep6(epspg6,vsig133,.false.,epspg6p) call chrep6(epsps6,vsig133,.false.,epsps6p) c conditions de positivite des valeurs normales pour les fissures do j=1,3 epspg6p(j)=max(epspg6p(j),0.d0) epsps6p(j)=max(epsps6p(j),0.d0) epspt6p(j)=max(epspt6p(j),0.d0) end do c ****** traitement anti-battement ********************************* if (ncal.ge.3) then c stockage des criteres actifs precedents do i=1,nc actif0(i)=actif1(i) end do end if if (ncal.ge.2) then c stockage des criteres actifs precedents do i=1,nc actif1(i)=actif(i) end do end if c ****** evaluation des criteres de plasticite ********************* call crir3d(bgel,pg1,bdef,ps1,bw1,pw1,rteff,cwrt,kwrt,reft, # rceff,kwrc,delta,beta,rteff_app,rceff_app,actif,factif,nc,err1, # na,sig13,cgel,hplg,cdef,hpls,cshr,hplw,epspg6p,epsps6p,epspt6p, # epspc6p,alphadp,betadp,refermtl,taudp,taudpmax,complet,CPTCR) c traitement de l erreur eventuelle if(err1.eq.1) then print*,'Pb crir3d dans Fluendo3d' ierr1=1 return end if c print*,'Etape 9 ds fldo3d apres crir3d' c ****** traitement anti-battement ********************************* battcr3d=.false. if(ncal.ge.3) then do i=1,nc if(actif(i).and.actif0(i)) then c print*,'Fluendo3d Battement critere',i ,'iter',ncal battcr3d=.true. c on active la compatibilité de ce critere avec tous c les autres pour la suite des calculs do j=1,nc CPTCR(i,j)=.true. CPTCR(j,i)=.true. end do end if end do end if if(na.gt.0) then C print*,'Etat des criteres dans fldo3d: na=',na C do i=1,nc C print*,'actif(',i,'):',actif(i),'factif',factif(i) C end do if(actif(13)) then log_plastc=.true. end if end if c read* c ****** ecoulement plastique le cas échéant *********************** if(na.gt.0) then c *calcul des increments dans la base principale des contraintes c *re-assemblage des matrices de fluage avec increment c deformations totales nulles et initialisation de depsp c on résoud le retour radial sans reconsiderer les termes c sources car ils ont deja été consideres lors du tir dphigel=0.d0 dphidef=0.d0 dphiw=0.d0 c reevaluation des coeff de consolidation dans la direction c principale des contraintes effectives avec les deformations c initiales (pour conserver les CC du tir) cf. call conso3d(xkm,epsm,epser,epsm06,epse_tild6,cc03, # vsig133,CMp00,CTHp,CTHv,Cwp) c print*,'fluendo3 etape2 cmp=',cmp00 c assemblage du systeme et resolution sans increment de deformations (zero3) c et sans la source de fluage due aux def elastiques initiales (zero6) call coupl3d(ann,xn,bn,ngf,err1,zero3,vsig133,vsig133t,dt, # cc03,taum,taumdtt,tauk,psik,zero6,zero6,raideur66,Mgel,bgel, # dphigel,Cgel,Hplg,Mdef,bdef,dphidef,Cdef,Hpls,Mwat,bw,dphiw, # cshr,alphadp,betadp,actif,factif,nc,dsigef6,depse6,depsk6, # depsm6,depspc6,depspt6,depsvt6,depspg6,depsps6,.false., # ipzero,dpg,dps,dpw,epspt6p,refermtl,complet,na,CMp00) c affichage def plastique de Drucker Prager C do i=1,6 C print*,'depspc6(',i,')=',depspc6(i) C end do c traitement de l erreur eventuelle if(err1.eq.1) then ierr1=1 print*,'Pb lors du retour radial ds Fluendo3d' print*,'Etat des criteres lors du pb:' do i=1,nc print*,'critere:',actif(i),'valeur',factif(i) end do print*,'contrainte effective matrice' do i=1,3 print*,'sig(',i,')=',sig13(i) end do return end if c print*,'Etape 10-1 ds fldo3d apres coupl3d iter',ncal else c pas d ecoulement plastique do i=1,6 depse6(i)=0.d0 depsk6(i)=0.d0 depsm6(i)=0.d0 depsvt6(i)=0.d0 depspg6(i)=0.d0 depsps6(i)=0.d0 depspc6(i)=0.d0 depspt6(i)=0.d0 dsigef6(i)=0.d0 end do dpg=0.d0 dps=0.d0 dpw=0.d0 c print*,'Etape 10-2 ds fldo3d apres coupl3d iter',ncal end if c ****************************************************************** c ***** actualisation des deformations (base fixe) ***************** do i=1,6 c elastique epse16(i)=epse16(i)+depse6(i) varf(i+1)=epse16(i) c kelvin epsk16(i)=epsk16(i)+depsk6(i) varf(i+7)=epsk16(i) c stockage depsk pour mise au point varf(122+i)=varf(i+7)-var0(i+7) if((dt.eq.0.).and.(abs(varf(122+i)).gt.1.0d-10)) then print*,'Pb dans fldo3d Etape 10-3 avec dt=0' print*,'depsk(',i,')=',varf(122+i) print*,'devrait etre nul !' ierr1=1 return end if c maxwell epsm16(i)=epsm16(i)+depsm6(i) varf(i+13)=epsm16(i) c cas des deformations plastiques c cisaillement et dilatance epspc6(i)=epspc6(i)+depspc6(i) varf(19+i)=epspc6(i) c print*,'actualisation varf epspc',i,varf(19+i) c traction base fixe epspt6(i)=epspt6(i)+depspt6(i) varf(25+i)=epspt6(i) c dtt epsvt6(i)=epsvt6(i)+depsvt6(i) varf(31+i)=epsvt6(i) c rag epspg6(i)=epspg6(i)+depspg6(i) varf(37+i)=epspg6(i) c rsi epsps6(i)=epsps6(i)+depsps6(i) varf(43+i)=epsps6(i) end do c ****** actualisation deformation equivalente de compression ****** depleqc=0.d0 do i=1,6 depleqc=depleqc+(epspc6(i)-epspc60(i))**2 end do depleqc=sqrt(depleqc*2.d0/3.d0) c print*,'depleqc',depleqc c actualisation et stockage epleqc1=max(epleqc0+depleqc,epleqc0) varf(94)=epleqc1 c ****** actualisation des pressions ******************************* c gel rag pg1=pg1+dpg varf(83)=pg1 c stockage pression maxi du gel varf(140)=max(var0(140),pg1) c def ps1=ps1+dps varf(85)=ps1 c depression capillaire pw1=pw1+dpw varf(81)=pw1 c **** etat des contraintes effectives apres tir visco elastique do i=1,6 sig16(i)=sig16(i)+dsigef6(i) do j=1,6 sigke16(i)=sigke16(i)+raideur66(i,j)*depsk6(j)*psik end do c stockage dans variables internes varf(49+i)=sig16(i) varf(55+i)=sigke16(i) C print*,'apres test criter iter',ncal C print*,'sig16(',i,')=',sig16(i) C print*,'sigke16(',i,')=',sigke16(i) end do c ****************************************************************** c sous iteration en cas de plasticite if(na.ne.0) then c on re-teste les criteres apres cette iteration if(ncal.le.nmax) then if (ncal.ge.nmax1) then c on commence a afficher les pb eventuels print*,'Iter dans Fluendo3d:',ncal c print*,'Avant ecoulement :', ncal+1 do i=1,nc if(actif(i)) then print*,'critere',i,actif(i),factif(i) if((i.ge.10).and.(i.le.12)) then do j=1,3 print*,'defpt',j,epspt6p(j) end do end if end if end do goto 30 do i=1,6 print*,'sig(',i,')=',sig16(i) end do print*, 'Apres ecoulement:' , ncal do i=1,6 print*,'epspt6(',i,')=',epspt6(i) end do do i=1,6 print*,'epspc6(',i,')=',epspc6(i) end do do i=1,6 print*,'epsvt6(',i,')=',epsvt6(i) end do do i=1,6 print*,'epspg6(',i,')=',epspg6(i) end do do i=1,6 print*,'epsps6(',i,')=',epsps6(i) end do print*,'----------------------------' c read* 30 continue end if c on va faire une sous iteration plastique goto 20 else print*,'Approximation plasticite ds Fluendo3d' print*,'apres',ncal,' iterations' do i=1,nc print*,'critere',actif(i),factif(i) end do print*,'Nbr maxi de sous iter atteinte dans Fluendo3d' print*,'Augmenter NMAX ds Fluendo3d ou modifier le Pb' ierr1=1 return end if end if c*********************************************************************** c fin des procedures visco-elasto-plastiques c print*,'Etape 11 ds fldo3d apres convergence',ncal c*********************************************************************** C print*,'dt ds fldo3d=',dt C print*,'Fluendo3d etape 10-2 a la convergence locale' C do i=1,6 C print*,'sig16(',i,')=',sig16(i) C print*,'sigke16(',i,')=',sigke16(i) C print*,'deps6(',i,')=',deps6(i) C print*,'depse6(',i,')=',depse6(i) C print*,'depsk6(',i,')=',depsk6(i) C print*,'depsm6(',i,')=',depsm6(i) C print*,'depsvt6(',i,')=',depsvt6(i) C print*,'depspc6(',i,')=',depspc6(i) C print*,'depspt6(',i,')=',depspt6(i) C print*,'depspg6(',i,')=',depspg6(i) C print*,'depsps6(',i,')=',depsps6(i) C end do C read* c*********************************************************************** c actualisation des taux de franchissement par depression capillaire c passage en matrice 3x3 call x6x33(fshr06,fshr33) c passage en base principale des contraintes call chre3(fshr33p,fshr33,vsig133) c actualisation des taux de passage des seuils en contrainte do i=1,3 fshr33p(i,i)=max(fshr33p(i,i),factif(6+i),0.d0) end do c retour en base fixe et stockage call chre3(fshr33,fshr33p,vsig133t) c retour en vecteur 6 call x33x6(fshr33,fshr6) c stockage do i=1,6 varf(129+i)=fshr6(i) end do c*********************************************************************** c actualisation des ouvertures de fissures actuelles et maxi de traction c recupereration des ouvertures initiales do i=1,6 c ouverture actuelle wplt06(i)=var0(96+i) c print*,'av ma:',wplt06(i),var0(96+i),96+i c ouverture maxi sur le temps wpltx06(i)=var0(102+i) c print*,'av ma:',wpltx06(i),var0(102+i),102+i end do c mise a jour des ouvertures des 3 fissures localisees if(ENDO_NL) then c ************ deformation plastique non locale ***************** alpha_nl=1.d0 c rangement des déformation plastiques non locales en tableau 3*3 call x6x33(eps_nl6,eps_nl33) c diagonalisation tenseur def non locales call b3_v33(eps_nl33,eps_nl3,veps_nl33) c creation de la matrice de passage inverse call traps1(veps_nl33t,veps_nl33,3) c Passage de la déformation locale dans la base principale de la déformation non locale c Rangement des déformation plastiques locales en tableau 3*3 call x6x33(epspt6,epspt33) call chre3(epspt33p,epspt33,veps_nl33) c Indicateur et deformation non-locale equivalente c Chargement du epsilon non-local equivalent initial do i=1,6 eptmasq06(i)=var0(146+i) end do call x6x33(eptmasq06,eptmasq033) do i=1,3 do j=1,3 eptmasq33p(i,j)=0.d0 end do if((epspt33p(i,i).gt.alpha_nl*eps_nl3(i)).and.(epspt33p(i,i) # .gt.0.d0).and.(eps_nl3(i).gt.0.d0))then eptmasq33p(i,i)=eps_nl3(i) end if end do c retour en base fixe call chre3(eptmasq33,eptmasq33p,veps_nl33t) call x33x6(eptmasq33,eptmasq6) do i=1,6 varf(146+i)=eptmasq6(i) end do c ************ fin deformation plastique non locale ************* c print*,'fluendo endo_nl av majw3d', Method_N,LcH,Method_H c l endommagement et les ouvertures sont calculee avec la c deformation plastique non locales call majw3d(eptmasq06,eptmasq6,wplt06,wplt6,wpltx06,wpltx6, # wpl3,vwpl33,vwpl33t,wplx3,vwplx33,vwplx33t,XE3D,NBNMAX3D, # NBNB3D,IDIMB3D,Method_N,LcH,Method_H,err1) c print*,'fluendo endo_nl ap majw3d', Method_N,LcH,Method_H else call majw3d(epspt60,epspt6,wplt06,wplt6,wpltx06,wpltx6, # wpl3,vwpl33,vwpl33t,wplx3,vwplx33,vwplx33t,XE3D,NBNMAX3D, # NBNB3D,IDIMB3D,Method_N,LcH,Method_H,err1) end if if(err1.eq.1)then print*,'Pb lors du calcul ds majw3d dans fldo3d' ierr1=1 return end if do i=1,6 varf(96+i)=wplt6(i) varf(102+i)=wpltx6(i) c print*,'fldo3d ap maj finale',wpltx6(i) end do c ouverture maxi actuelle varf(109)=max(wpl3(1),wpl3(2),wpl3(3)) c indicateur de progression de la fissure if(abs(varf(109)-var0(109)).gt.coeff_petit) then travail_fissure=.true. else travail_fissure=.false. end if C print*,'Etape 12 ds fldo3d apres majw3d',ncal C read* c*********************************************************************** c calcul des contraintes dans les fibres ( Romain Gontero) c include './contraintes_dans_les_fibres.h' -INC HSIGFIB c*********************************************************************** c contraintes dans solide et rgi en fin de pas avec c prise en compte de l endo thermique et de fluage umdt=(1.d0-dth1)*(1.d0-dflu1) c resultante des pressions intraporeuses RGI et c Capillaire (depression) sigp=-bw1*pw1-bg1*pg1-bs1*ps1 c print*,bw1,pw1,bg1,pg1,bs1,ps1 c print*,'fluendo sigp',sigp c effet sur la contrainte apparente en non sature do i=1,6 if(i.le.3) then c prise en compte de la pression rgi sigf6(i)=(sig16(i)+sigp)*umdt else sigf6(i)=sig16(i)*umdt end if c print*,'sigf6',sigf6(i),' sigp',sigp,'umdt',umdt end do c*********************************************************************** c prise en compte de l'endommagement mécanique if (end3d.and.(.not.plast_seule)) then c chargement de la variable de controle d erreur de gf errgf0=var0(110) c endo pre pic traction dtr=dtr0 c endo pre pic de compression dcpp=dcpp0 c actualisation de gft en fonction de la depression capillaire gft_app=gft*rt_app/rt c print*,'fldo3d Method_H', Method_H c read* c calcul des contraintes endommagee et des endommagements call kend3d(wpl3,vwpl33,vwpl33t,wplx3,vwplx33,vwplx33t, # gft_app,gfr,iso,sigf6,sigf6d,err1,reft,young,nu,souplesse66, # xmt,dtiso,dtr,dt3,dr3,dct3,0.,dcc3,nc3,nc33,ept,errgf1, # errgf0,epspc6,0.d0,ekdc,epspg6,ekdg,alphag,dgt3,dgc3,fshr6, # dwt3,dwc3,skdw,0.d0,epsps6,dst3,dsc3,ekds,alphas,rteff_app, # ann,bn,xn,ipzero,ngf,tau_trac,xmc,dciso,dcpp,taudp,taudpmax, # XE3D,NBNMAX3D,NBNB3D,IDIMB3D,altc,epssr,Method_H,Method_N,LcH) c stockage variable de controle de l erreur sur Gft varf(110)=errgf1 c endo pre pic varf(79)=dtr c endo pre pic de compression - cisaillement varf(122)=dcpp c traitement des erreurs if(err1.eq.1) then print*,'pb calcul des endommagements dans fluendo3d' ierr1=1 return end if c contraintes totale dans la matrice apres endommagement do i=1,6 sigmf6(i)=sigf6d(i) varf(61+i)=sigmf6(i) end do else c pas d endommagement dcpp=0.d0 dtr=0.d0 varf(79)=dtr varf(122)=dcpp do i=1,3 dt3(i)=0.d0 nc3(i)=1.d0 dr3(i)=0.d0 dwt3(i)=0.d0 dwc3(i)=0.d0 dgt3(i)=0.d0 dgc3(i)=0.d0 dst3(i)=0.d0 dsc3(i)=0.d0 dcc3(i)=0.d0 dct3(i)=0.d0 end do do i=1,6 sigmf6(i)=sigf6(i) varf(61+i)=sigmf6(i) end do end if c*********************************************************************** c combinaison fibres matrice endommagee (Romain Gontero 2022) c include './combinaison_matrice_fibres.h' -INC HCOMFIB c*********************************************************************** c contrainte dans la matrice (eventuellement fibree) do i=1,6 c debut romain if(ifibre.eq.1) then sigmf6(i)=(1.d0-rhof)*sigmf6(i)+sfib(i) end if c print*,"combinaison contrainte fibres" c fin romain varf(61+i)=sigmf6(i) end do c*********** variables initiales pour les renforts ********************* c lignes a copier si le modele utilise des renforts if (NRENF00.ne.0) then do i=1,NRENF00 c include './lire_vari_renforts.h' -INC HLVIREN c recuperation des non locales a utiliser pour les renforts c pour les renforts le non local est fait sur SER(i) de numero c NVAR2+(i-1)*NB_VARI_PAR_RENF+26 c on teste si ce renfort est traite par helmholtz if(log_H_RENF(i)) then c recuperation des variables pour une traitement non local du renfort ivarnl=NVAR2+(i-1)*NB_VARI_PAR_RENF+26 c numero de la formulation non locale nl=Num_H_RENF(i) c position des variables de Helmholtz associees a cette formulation if(istep.ne.0) then c recup des non locales associees aux renforts if(istep.ne.1) then c recuperation de la contrainte non locale intermediaire sigr_nl(i)=var0(ivarnl) c recuperation du gradient de la contrainte sur Hi2 grad_sigr(i)=Helm0(nl,2) c survie du renfort coeff_nlr(i)=Helm0(nl,3) c survie de l interface coeff_nli(i)=Helm0(nl,4) c coeff de taux dans l erreur Helm1(nl,6)=Helm0(nl,6) c contrainte non locale de l iteration non locale precedente sigr_nl0(i)=Helm0(nl,7) c recuperation de la variable de controle des hypotheses Pb_nl(i)=Helm0(nl,8) c hypothese ecoulement local hypo_nl(i)=Helm0(nl,9) c copie du coeff de normalisation de taux helm1(nl,10)=helm0(nl,10) else c cas istep=1, on a pas encore fait d etape non locale c contrainte non locale de l iteration non locale precedente sigr_nl0(i)=Helm0(nl,1) c recuperation de la contrainte non locale convergee precedente sigr_nl(i)=sigr_nl0(i) c recuperation du gradient de la contrainte non locale convergee precedente sur YU2i grad_sigr(i)=Helm0(nl,2) c on recharge les hypotheses sur endo et refermeture if(ppas) then c survie du renfort coeff_nlr(i)=1.0d0 c survie de l interface coeff_nli(i)=1.0d0 else c survie du renfort coeff_nlr(i)=Helm0(nl,3) c survie de l interface coeff_nli(i)=Helm0(nl,4) end if c hypothese pas d ecoulement plastique Pb_nl(i)=0.d0 hypo_nl(i)=0.d0 end if else c traitement local du renfort grad_sigr(i)=0.d0 coeff_nlr(i)=1.d0 coeff_nli(i)=1.d0 Pb_nl(i)=0.d0 hypo_nl(i)=0.d0 end if else c traitement local du renfort grad_sigr(i)=0.d0 coeff_nlr(i)=1.d0 coeff_nli(i)=1.d0 Pb_nl(i)=0.d0 hypo_nl(i)=0.d0 end if end do end if c ****** variables initiales pour WL2 ****************************** if (log_H_TWL2) then c resultat non local xwl2=var0(71) c max de alpha xwl2max=var0(73) c max du taux de charge tau_trac_wl2_max=var0(74) c endo localise dtl00=var0(75) c coeff d effet d echelle precedent cwrt=var0(76) end if c************* Etapes non locales intermediaires *********************** if((istep.ne.0).and.(istep.ne.2)) then c *** reinitialisation des variables locales pour le calcul reel c car les sous iterations non locales ne doivent pas conserver c les variables locales non convergees lors de la copie de varf c sur var0 lors du passage dans UNPAS en non local sauf si c istep est egal a 2 (il s agit alors du dernier passage non local) do i=1,NVARI varf(i)=var0(i) end do c **** etape non locale intermediaire pour WL2 ***************** c 'TWL2' est la variable a moyenner pour la methode wl2 ivarnl=ivari_twl2 call exhmtz(istep,NBVIA3D,INLVIA3D,NB_HELM, # ivarnl,log_H_TWL2,Num_H_TWL2) c numero de la formulation associee NL=Num_H_TWL2 if (log_H_TWL2) then c ** le calcul de xwl2 est a faire ************************** c print*,'fldo3d, istep=1 avec cwrt=',cwrt c recuperation des contraintes dans la matrice call x6x33(sigmf6,sig133) call b3_v33(sig133,sig13,vsig133) c prise en compte de la dilution d endommagement lors de l hydratation call hydrvari3d(dtl00,dtl0,hyd0,hydr,hyds) c recuperation de l indicateur de localisation fin du nouveau pas dtl1=max(dt3(1),dt3(2),dt3(3)) c taux de charge en traction actualise pour la nouvelle deformation tau_trac_wl2=tau_trac c pression capillaire pw0=pw1 c saturation en eau * biot bw0=bw1 c 1er passage dans wl2 call wl2d3d(1,vmax,vref,sig13,eweb,rt_app,cwrt,xwl2, # tau_trac_wl2,dtl0,dtl1,xwl2max,tau_trac_wl2_max,bw0, # pw0,kwrt) c *** on actualise des vari a passer dans non local pour Weibull c les variables de Weibull (alpha(xwl2),smax(tau_trac_wl2)) c alpha varf(ivarinl)=xwl2 c sigma/rt maxi pour max(structure) print*,'a refaire avec les helm1...' ierr1=1 return varf(72)=tau_trac_wl2 c max xwl2 un pour multiplier par xwl2 non local maxi c ds la procedure non locale varf(73)=1.d0 c max twl2 un pour multiplier par twl2 non local maxi c ds la procedure non locale varf(74)=1.d0 c reconduction de l endommagement localise debut de ne c pas affecte par hydratation varf(75)=dtl0 c reconduction de cwrt en cas de localisation varf(76)=cwrt c pression capillaire fin de pas pour rt istep=2 varf(81)=pw1 c saturation en eau * biot fin de pas pour istep=2 varf(80)=bw1 c endommagement fin de pas varf(119)=dtl1 else c ** rt00 n est pas affectee par la methode WL2 ************* c pas de non local pour xwl2 cwrt=1.d0 c liste des variables non locales : alpha varf(71)=0.d0 c variable non locale : sigma maxi /rt varf(72)=0.d0 c max de xwl2 non utilisse on y laisse 1. varf(73)=1.d0 c contrainte relative maxi non locale non utiliseee varf(74)=0.d0 c pas d endo de traction utilisable pour xwl2 varf(75)=0.d0 c modif de rt par effet d echelle inactif varf(76)=cwrt end if c *** etape non locale intermediaire pour les renforts ********* if (NRENF00.ne.0) then c intialisation des tenseurs pour les renforts call prrf3d(dalpr,teta1,teta2,tetaref,epst06,epst033, # epstf6,epstf33,epspmf6,epspt6,epspmf33) c calcul de la source non locale pour chaque renfort do i=1,NRENF00 c cas des renforts traites par helmholtz if (log_H_RENF(i)) then c numero de la formulation d Helmholtz pour le renfort i nl=Num_H_RENF(i) c recuperation des non locales a utiliser pour les renforts c pour les renforts le non local est fait sur SERi de numero ivarnl=NVAR2+(i-1)*NB_VARI_PAR_RENF+26 if(rhor(i).gt.0.) then c etape non locale de rnfr3d on calcule la source c non locale pour le pas suivant call rnfr3d(istep,NB_RENF,i,epstf33,vecr,epsr0(i), # epsrf(i),eplr0(i),eplrf(i),yor(i),syr(i),sigre0(i), # sigref(i),hplr(i),tor(i),ekr(i),skr(i),ATRR(i),khir(i), # gamr(i),sprec(i),teta1,teta2,dt,ppas,theta,eprm0(i), # eprmf(i),ttaref(i),rhor(i),mu_r0(i),fl3d,errr,xnr(i), # xmuthr(i),eprk0(i),eprkf(i),tokr(i),yksyr(i), # plast_seule,ann,xn,bn,ngf,ipzero,epspmf33,spre0(i), # spref(i),epst033,depsa(i),epleq0(i),epleqf(i), # epsmf(i),epsm0(i),sigr_nl(i),coeff_nlr(i),eper0(i), # epart(i),eparv(i),hypo_nl(i),plr_iso,log_H_RENF(i), # garder,Pb_nl(i),travail_fissure) c traitement de l erreur if(errr) then print*,'erreur dans rnfr3d' ierr1=1 return end if c mise a jour de la def anelastique pour la source de c l iteration suivante if(garder) then c deformation visqueuse eparv(i)=eprkf(i)+eprmf(i) c deformation anelastique totale epart(i)=eparv(i)+eplrf(i) else c deformation visqueuse initiale eparv(i)=eprk0(i)+eprm0(i) c deformation anelastique totale initiale epart(i)=eparv(i)+eplr0(i) end if c stockage de la longueur de l element fini varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+13)=longr(i) c calcul des endommagements et de la correction de source call endorf(sigref(i),eplrf(i),syr(i),sur(i), # epsrf(i),epu(i),hplr(i),yor(i),wpr(i),longr(i), # damrt0(i),damrtf(i),damrt1(i),damrt2(i),damrc0(i), # damrcf(i),refrf(i),eplrt(i),eplrc(i),plr_iso, # diffr0(i),diffrf(i),yor(i),deqr(i),Hintr(i), # grad_sigr(i),tauie(i),taur(i),dami0(i),damif(i), # epsmf(i),eplr0(i),sigre0(i),log_H_RENF(i),epsm0(i), # epsr0(i),depsa(i),Et0(i),Etf(i),Hs0(i),Hsf(i),sigrf(i), # spref(i),epart(i),eparv(i),source_nl(i),coeff_nlr(i), # taui0(i),tauif(i),coeff_nli(i),eprkf(i), # eprmf(i),hypo_nl(i),istep,ierr_renfort,garder, # endo_interface_explicit) c traitement erreur endo renfort c include './traitement_erreur_renforts.h' -INC HERRREN c sauvegarde de la contrainte a moyenner varf(ivarnl)=source_nl(i) c survie du renfort : coeff_nlr(i) Helm1(nl,3)=coeff_nlr(i) c survie de l interface : coeff_nli(i) Helm1(nl,4)=coeff_nli(i) c hypotheses d ecoulements realises : hypo_nl(i) Helm1(nl,9)= hypo_nl(i) c endommagement de l interface pour le calcul implicite du gradient endommagement d interface : damif(i) varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+20)=damif(i) c contrainte de cisaillement implicite varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+19)=tauiF(i) c sauvegarde de la contrainte if(garder) then c les hypotheses sont coherentes on garde la solution c valeur de la contrainte non locale avant ecoulement (SERi) Helm1(nl,1)=sigr_nl(i) c valeur de la contrainte endommagee (SNRi) varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+14)= sigrf(i) c actualisation de la diffusion non locale (DFRi) varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+18)=diffrf(i) c stockage du module tangent du renfort ETRi varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+21)=Etf(i) c stockage module secant interface HSRi varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+22)=Hsf(i) c deformation anelastique EARi varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+24)=epart(i) c stockage de la deformation visqueuse du renfort EVRi varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+25)=eparv(i) else c les hypotheses etaient incoherentes, on reprend la c valeur de la contrainte locale avant l etape non locale if(istep.ne.1) then c les hypotheses etaient incoherente on reprend la c valeur de la contrainte locale avant l etape non locale Helm1(nl,1)=sigr_nl0(i) else c comme on a pas encore fait d etape non locale c on prend la contrainte calculée localement Helm1(nl,1)=sigrf(i) endif end if if(endo_interface_explicit.and.(istep.eq.1)) then varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+20)=dami0(i) varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+19)=taui0(i) end if else c pas de renfort num i c rho(i)=0 do j=1,9 c mise a zero des vari du renfort varf(ivarnl)=0.d0 c valeur de la contrainte apres ecoulement sigrf(i)=0.d0 varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+14)=sigrf(i) c actualisation de la diffusion non locale diffrf(i)=0.d0 varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+18)=diffrf(i) c stockage du module du renfort Etf(i)=0.d0 varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+21)=Etf(i) c stockage module interface Hsf(i)=0.d0 varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+22)=Hsf(i) c stockage de l hypothese de deformation anelastique epart(i)=0.d0 varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+24)=epart(i) c stockage de l autre hypothese eparv(i)=0.d0 varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+25)=eparv(i) end do end if end if end do end if end if c******** fin l etape non locale intermidiaire ************************* c******** Derniere etape non locale ou calcul local ******************** c ***** calcul des contraintes dans les renforts ****************** if ((NRENF00.ne.0).and.((istep.eq.0).or.(istep.eq.2))) then c calcul du taux volumique d armatures rhov=0.d0 do j=1,NRENF00 rhov=rhov+rhor(j) end do if (rhov.gt.1.d0) then print*,'Pb car taux d armature > 1 dans fluendo3d ! :(' ierr1=1 return end if c initialisation des tenseurs de contraintes homogeneisee do j=1,3 do k=1,3 sigrm33(j,k)=0.d0 sigrf33(j,k)=0.d0 end do end do c intialisation des tenseurs pour les renforts call prrf3d(dalpr,teta1,teta2,tetaref,epst06,epst033, # epstf6,epstf33,epspmf6,epspt6,epspmf33) c calcul des contraintes axiales dans chaque renfort do i=1,NRENF00 c recuperation des non locales a utiliser pour les renforts c pour les renforts le non local est fait sur SER de nuemero ivarnl=NVAR2+(i-1)*NB_VARI_PAR_RENF+26 if(rhor(i).gt.0.) then c recuperation du tir elastique non local if((istep.eq.2).and.log_H_RENF(i)) then c derniere estimation pour les renforts non locaux cc sauvegarde du tir non local precedent car istep=2 cc avec la precedente contrainte non locale c numero de la formulation non locale nl=Num_H_RENF(i) c contrainte non locale sigr_nl(i)=var0(ivarnl) c on utilse helm1(nl,1) pour stocker le resultat non local Helm1(nl,1)=sigr_nl(i) c recuperation du gradient grad_sigr(i)=Helm0(nl,2) c sauvegarde du gradient de sigr pour le pas suivant helm1(nl,2)=grad_sigr(i) c copie de la valeur initiale en cas de reprise helm1(nl,7)=sigr_nl(i) c coeff du gradient pour l erreur Helm1(nl,10)=0.25d0*(deqr(i)/sur(i)) c recuperation des predicteurs pour le pas suivant Etf(i)=Et0(i) Hsf(i)=Hs0(i) else if (istep.eq.0) then continue else if ((istep.eq.2).and.(.not.log_H_RENF(i))) then c certains aciers sont calcules en local continue else c recuperation des predicteurs pour le pas suivant print*,'erreur d aiguillage dans fldo3d' print*,'istep',istep,'Helmholtz(',nl,')=', # log_H_RENF(i) ierr1=1 return end if c contrainte effective axiale dans les renforts call rnfr3d(istep,NB_RENF,i,epstf33,vecr,epsr0(i), # epsrf(i),eplr0(i),eplrf(i),yor(i),syr(i),sigre0(i), # sigref(i),hplr(i),tor(i),ekr(i),skr(i),ATRR(i),khir(i), # gamr(i),sprec(i),teta1,teta2,dt,ppas,theta,eprm0(i), # eprmf(i),ttaref(i),rhor(i),mu_r0(i),fl3d,errr,xnr(i), # xmuthr(i),eprk0(i),eprkf(i),tokr(i),yksyr(i), # plast_seule,ann,xn,bn,ngf,ipzero,epspmf33,spre0(i), # spref(i),epst033,depsa(i),epleq0(i),epleqf(i), # epsmf(i),epsm0(i),sigr_nl(i),coeff_nlr(i),eper0(i), # epart(i),eparv(i),hypo_nl(i),plr_iso,log_H_RENF(i), # garder,Pb_nl(i),travail_fissure) c traitement erreur rnfr3d if(errr) then print*,'erreur dans rnfr3d' ierr1=1 return end if c action suivant conformite hypothese if(.not.garder) then print*,'Pb dans fldo3d' print*,'pour istep=2 on devrait garder' print*,'la solution du sous pas non local' ierr1=1 return end if c deformation visqueuse eparv(i)=eprkf(i)+eprmf(i) c deformation anelastique totale epart(i)=eparv(i)+eplrf(i) c nouvelle deformation elastique eperf(i)=epsrf(i)-epart(i) c print*,'fldo3d:sigref(',i,')=',sigref(i) c calcul des endommagements et actualisation lnl, dsource call endorf(sigref(i),eplrf(i),syr(i),sur(i), # epsrf(i),epu(i),hplr(i),yor(i),wpr(i),longr(i), # damrt0(i),damrtf(i),damrt1(i),damrt2(i),damrc0(i), # damrcf(i),refrf(i),eplrt(i),eplrc(i),plr_iso, # diffr0(i),diffrf(i),yor(i),deqr(i),Hintr(i), # grad_sigr(i),tauie(i),taur(i),dami0(i),damif(i), # epsmf(i),eplr0(i),sigre0(i),log_H_RENF(i),epsm0(i), # epsr0(i),depsa(i),Et0(i),Etf(i),Hs0(i),Hsf(i),sigrf(i), # spref(i),epart(i),eparv(i),source_nl(i),coeff_nlr(i), # taui0(i),tauif(i),coeff_nli(i),eprkf(i), # eprmf(i),hypo_nl(i),istep,err1,garder, # endo_interface_explicit) C if (sigrf(i).ne. sigref(i)) then C print*,'Ds fluendo3d sigrf(i).ne. sigref(i)' C print*,i,sigrf(i),sigref(i) C ierr1=1 C return C end if c print*,'Dans fldo3d,sigrf(',i,')=',sigrf(i) c traitement erreur endommagement if(err1.eq.1) then ierr1=1 return end if if((istep.eq.2).and.log_H_RENF(i)) then c rappel numero de la formulation non locale nl=Num_H_RENF(i) c contrainte non locale derniere iteration Helm1(nl,1)=sigr_nl(i) c gradient de la contrainte Helm1(nl,2)=grad_sigr(i) c survie du renfort Helm1(nl,3)=coeff_nlr(i) c survie de l interface Helm1(nl,4)=coeff_nli(i) c coeff de taux dans l erreur Helm1(nl,6)=1.d0/sur(i) c valeur de la contrainte non locale utilisee Helm1(nl,7)=sigr_nl(i) c numero sous iteration non locale Helm1(nl,8)=Pb_nl(i) c moyenne des hypotheses d ecoulements realises Helm1(nl,9)=hypo_nl(i) c print*,'fluendo istep',istep,'errsigr',Helm1(nl,4) end if else epsrf(i)=0.d0 eperf(i)=0.d0 eplrf(i)=0.d0 sigref(i)=0.d0 eprmf(i)=0.d0 mu_r0(i)=0.d0 eprkf(i)=0.d0 spref(i) =0.d0 eplrc(i) =0.d0 damrtf(i)=0.d0 damrt1(i)=0.d0 damrt2(i)=0.d0 sigrf(i)=0.d0 epleqf(i) =0.d0 damrcf(i) =0.d0 eplrt(i)=0.d0 diffrf(i)=0.d0 tauif(i)=0.d0 damif(i)=0.d0 Etf(i)=0.d0 Hsf(i)=0.d0 refrf(i)=0.d0 epart(i)=0.d0 eparv(i)=0.d0 end if c print*,'fldo3d Etape 14 istep2 eper:',eperf(i) c stockage des variables internes des renforts c include './ecrire_vari_renforts.h' -INC HEVIREN c print*,'Etape 13-3 istep 2 ds fluendp3d' c print*,'vari sigrf',i,varf(nvar2+(i-1)*NB_VARI_RENF+14),sigrf(i) c **** tenseur des contraintes dues aux renforts ************ do j=1,3 do k=1,3 if(rhor(i).gt.0.) then sigrm33(j,k)=sigrm33(j,k)+rhor(i)*sigrf(i)* # vecr(i,j)*vecr(i,k) end if end do end do c fin de la boucle sur le renfort i end do c mise a zero des vari des renforts non actifs if (NRENF00.lt.NB_RENF) then do i=(NRENF00+1),NB_RENF do j=1,NB_VARI_PAR_RENF varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+j)=0.d0 end do end do end if else c pas de renforts ou etape non locale intermediaire rhov=0.d0 c **** tenseur des contraintes dues aux renforts ************ do j=1,3 do k=1,3 sigrm33(j,k)=0.d0 end do end do c variables internes des renforts if(NRENF00.eq.0) then c mise a zero des vari des renforts do i=1,NB_RENF do j=1,NB_VARI_PAR_RENF varf(NVAR2+(i-1)*NB_VARI_PAR_RENF+j)=0.d0 end do end do end if end if c combinaison des contraintes dans les renforts et la matrice c print*,'contrainte dans les aciers' call x33x6(sigrm33,sigrm6) c call afic6(sigrm6) c print*,'contraintes dans la matrice ' c call afic6(sigmf6) do i=1,6 sigf6d(i)=(1.d0-rhov)*sigmf6(i)+sigrm6(i) end do c print*,'contrainte moyenne BA' c call afic6(sigf6d) c*********************************************************************** c calcul des endommagements equivalents pour affichage c*********************************************************************** if(end3d.and.(.not.plast_seule)) then c endommagement equivalent de traction dtl1=max(dt3(1),dt3(2),dt3(3)) if(istep.eq.0) then varf(75)=var0(119) varf(119)=dtl1 end if dtw0=max(dwt3(1),dwt3(2),dwt3(3)) dtg0=max(dgt3(1),dgt3(2),dgt3(3)) dts0=max(dst3(1),dst3(2),dst3(3)) dtra=1.d0-(1.d0-dtr)*(1.d0-dflu1)*(1.d0-dth1)* # (1.d0-dtl1)* # (1.d0-dtw0)* # (1.d0-dtg0)* # (1.d0-dts0) c endommagement equivalent de compression dcw0=max(dwc3(1),dwc3(2),dwc3(3)) dcg0=max(dgc3(1),dgc3(2),dgc3(3)) dcs0=max(dsc3(1),dsc3(2),dsc3(3)) dct0=max(dct3(1),dct3(2),dct3(3)) dcc0=max(dcc3(1),dcc3(2),dcc3(3)) dcom=1.d0-(1.d0-dflu1)*(1.d0-dth1)*(1.d0-dcpp)* # (1.d0-dcw0)* # (1.d0-dcg0)* # (1.d0-dcs0)* # (1.d0-dct0)* # (1.d0-dcc0) else c cas sans endommagement dtra=0.d0 dcom=0.d0 end if varf(95)=dtra varf(96)=dcom c endo meca traction varf(111)=1.d0-(1.d0-dtl1)*(1.d0-dtr) c endo hydrique traction varf(112)=dtw0 c endo rag traction varf(113)=dtg0 c end do def traction varf(114)=dts0 c endo meca compression varf(115)=1.d0-(1.d0-dct0)*(1.d0-dcc0)*(1-dcpp) c endo hydrique compression varf(116)=dcw0 c endo rag compression varf(117)=dcg0 c endo def compression varf(118)=dcs0 c endo meca prepic varf(122)=dcpp c endo meca compression par couplage avec traction (buckling) varf(136)=dct0 c contrainte equivalente de Drucker Prager varf(137)=taudp c endommagement de compression par cisaillement varf(138)=dcc0 c*********************************************************************** c affectation dans le tableau de sortie des contraintes c et prise encompte des contribution des renforts do i=1,nstrs sigf(i)=sigf6d(i) c print*,'sigf(',i,')=',sigf(i) end do c indicateur de fin de 1er pas if(ppas.and.((istep.eq.2).or.(istep.eq.0))) then varf(1)=1.d0 else if(ppas) then varf(1)=0.d0 else varf(1)=var0(1) end if end if c*********************************************************************** c VARIABLES INTERNES A ACTUALISER POUR LES FORMULATIONS DE HELMHOLTZ C*********************************************************************** c sauvegarde des vari de helmholtz HELM1(NL,ii=1..NBR_VARI_HELM) do NL=1,NB_HELM c include './ecrire_vari_helmholtz.h' -INC HEVIHEL end do c autre variables internes a sauver en cas d etape non locale if(ENDO_NL) then c residu de deformation plastique locale a diffuser if((istep.eq.1).or.(istep.eq.3)) then do i=1,6 c source locale varf(140+i)=epspt6(i) c -epspt60(i) end do else do i=1,6 c stockage valeur non locale varf(140+i)=eps_nl6(i) end do end if end if C*********************************************************************** c stockage du pas de temps pour verif varf(139)=dt c*********************************************************************** c traitement erreur residuelle if(err1.eq.1) then ierr1=1 return end if return c fin c*********************************************************************** c*********************************************************************** c affichages pour etude des cas ou dt = 0 (dt est stocke dans vari c (139) si dt = 0, dans castem les vari non convergees de la c dernieresi dt = 0, dans castem les vari non convergees de c la derniere sous iteration avant la convergence forcee conservee do i=1,6 if((varf(139).eq.0.).and.(.not.ppas)) then c le pas de temps actuel est nul on affiche les var0 du pas c precedent converge print*,'Fluendo3d Etape 14-1 avec dt=0 au pas actuel' print*,'depsk0(',i,')=',var0(122+i) print*,'depskf(',i,')=',varf(122+i) print*,'var0(',i,')=',var0(7+i) print*,'varf(',i,')=',varf(7+i) c read* c err1=1 end if end do do i=1,6 if((var0(139).eq.0.).and.(.not.ppas).and. # (abs(varf(7+i)-var0(7+i)).gt.1.0d-10)) then c on teste l evolution des vari visqueuses sur le pas c suivant la non covergence forcee, si on constate une c evolution des vari au pas actuel on l affiche print*,'Fluendo3d Etape 14-2 avec dt=0 au pas precedent' print*,'depsk0(',i,')=',var0(122+i) print*,'depskf(',i,')=',varf(122+i) print*,'var0(',i,')=',var0(7+i) print*,'varf(',i,')=',varf(7+i) c err1=1 end if end do c dernier traitement des erreurs if(err1.eq.1) then ierr1=1 return end if c cas ou on souhaite stopper le calcul si convergence forcee if((varf(139).eq.0.).and.(.not.ppas)) then print*,'on sort de fldo3d avec le pas de temps nul' print*,'---------------------------------------------' c ierr1=1 end if c variables internes fin de passage C do i=1,6 C write(*,'(I4,1Xe10.3,1xe10.3)') 19+i,var0(19+i),varf(19+i) C end do C if(log_plastc) then C print*,'---on vient de plastifier en compression------' C c read* C end if return end c***********************************************************************