Télécharger fibs3d.eso

Retour à la liste

fibs3d
C FIBS3D    SOURCE    GOUNAND   26/08/19    21:15:02     12620          
      subroutine fibs3d(Df,Lf,lech,mw,rtec,rhof,eo,vf1,vf2,Gf,
     #cokf1,cokf2,cowp1,cowp2,cofp,cow03,cowf,mink01,maxk01,mink02,
     #maxk02,minwp1,maxwp1,minwp2,maxwp2,minfp,maxfp,minw03,maxw03,
     #minwf,maxwf,sigf06,vwpl6,wmax,wmoy,etw,ncf,sigfr,wmax6,le,
     #nc0,wmoy0,etw0,vw30,rig0t,wmax0,vw3,rmax1,rmax2,rmax3,phimoy,
     #wferm,sigdec,iloc,wpsigt,wmiloc,sigmaxt,vecw33,sigfx6,wmaxd,
     #wmaxd0,epstf6,Ef,fu,young,xmt,dtiso,epict,rt00,cvar,wplmax)
 
      implicit real*8 (a-h,o-z)
      implicit integer (i-n)
      integer err1,nvari,ncalc,rmax,r,rmaxi,rmax1,rmax2,rmax3
      integer i,j,k,l,iloc(3)
      parameter (PI=3.141592653589793D0)
      real*8 incremwlis,precision1,precision2
      parameter (incremwlis=2.5d-6,precision1=1.0d-8,precision2=1.0d-8)
 
 
      real*8 rhof,Df,Lf,phicrit,Rtec,mw,Rtb,lech,le(3),Gf
      real*8 young,xmt,epict,dt3
 
      parameter (nvari=7,phicrit=55.d0,npos=10)
c attention npos doit être redefini dans disp3d
c Rtec est celui de lech
C       parameter(mu=0.5d0,Td=6.4d0,Tmax=8.d0,rhof=0.02,
C      # Ef=210000.d0,Df=2.0d-4,Hf=8000000.d0,sk=3.0d-2,
C      #  F0=5.0d-6,alphad=30.d0,M0=8.d0,
C      # Rtec=4.d0,fu=2800.d0,L0=5.0d-3,
C      # Lf=0.013d0,lech=1.d0,mw=12.d0)
c       parameter (ncalc=max(int(Lf/(incremwlis)),1000))
      parameter (ncalc=1000)
      parameter (rmaxi=10)
      real*8 vecw(3),wmax0(3)
      real*8 Eo(3),Vwpl6(6),vwpl33(3,3),vw3(3)
      real*8 veo(3,3),vecw33(3,3)
      real*8 stri(npos,2*npos-2),xg(npos,2*npos-2),yg(npos,2*npos-2)
      real*8 zg(npos,2*npos-2),vei(3)
      real*8 phi1,phi1d(npos,2*npos-2),f1(npos,2*npos-2)
      real*8 pscal(npos,2*npos-2),normv
      real*8 si0(npos,2*npos-2),si1,Reiw(npos,2*npos-2)
      real*8 wmaxa(3,3),wmaxac(3,3),wmaxact(3),wmax6(6),wmax33(3,3)
      real*8 rig0(3,rmaxi),sig0(3),vwfi(3,ncalc)
      real*8 sigmax(3,rmaxi),sigf(3,rmaxi,ncalc),sigf06(6)
      real*8 sigf033(3,3),sigf03(3)
      real*8 ratio(npos,2*npos-2),vwfi0(3),wfis(3)
      real*8 wpsig(3,rmaxi),wmi0(1),cvar(3),wplmax(3)
      real*8 wmiloc(3),phimoy(3),nf(npos,2*npos-2)
      real*8 Romeg(npos,2*npos-2),Sr(3),Sigb(3),Sig(3),wfist(3)
      real*8 nc(rmaxi),Rt(rmaxi),nct(rmaxi),Romeg1(npos,2*npos-2)
      real*8 wmi(rmaxi),wt(3,ncalc),sigt(3,ncalc),wmul(3,rmaxi)
 
c variables pour la reduction du modele
      real*8 k0m, wp,fp,f03,w03,wf
      real*8 cokf1(10),cokf2(10),cowp1(10),cowp2(10)
      real*8 cofp(10),cowf(10),cow03(10)
      real*8 maxfp,minfp,maxwp1,minwp1,maxwf,minwf,maxw03,minw03
      real*8 maxwp2,minwp2,B,vf1(3),vf2(3)
      real*8 maxk02,mink02,maxk01,mink01,phi2d(npos,2*npos-2)
      real*8 tri(npos,2*npos-2,3,3),Reiw1,Reiw2,Reiw3,vst1(3),vst2(3)
      real*8 vst3(3),theta1,gamma1,theta2,gamma2,theta3,gamma3
 
c    variable ajoutée pour contrainte totale
      real*8 wpsigt(3),rig0t(3),sigmaxt(3),sigfr(3),rt00
      logical loc(3),dtiso
      real*8 varft(nvari,ncalc),rig0td(3),gfap(3)
      real*8 wmr(rmaxi),wer(rmaxi),wmoy(3),etw(3),wmax(3)
      real*8 sigfp33(3,3),sigfx33(3,3),sigfx6(6),ncf(3)
      real*8 dw3(3),vw30(3),wmoy0(3),etw0(3),nc0(3)
      real*8 wferm(3),sigdec(3),sigfe(3),epse13(3),vecwj(3)
      real*8 wmaxd(3),wmaxd0(3),w1(3),rtmin,epstf6(6),epse133(3,3)
      real*8 vecs33(3,3)
* pv
      s1 = 0.d0
      do i=1,3
       sig(i)=-1.d0
      enddo
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
c CONSTRUCTION DES 3 COURBES CONTRAINTES - OUVERTURE TOTALES DANS LA BASE
c CCCCCCCCCCCCCCCC  DE FISSURE ACTUELLE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
 
c recuperation des ouvertures max
      call x6x33(wmax6,wmax33)
      call x6x33(vwpl6,vwpl33)
c     vw3 : vecteur des valeurs principales de fissures
c     vecw33 matrice des vecteurs d'orientation des fissures actuelles
      call b3_v33(vwpl33,vw3,vecw33)
 
C        print*,vecw33(1,1),vecw33(2,1),vecw33(3,1),"vw1"
C        print*,vecw33(1,2),vecw33(2,2),vecw33(3,2),"vw2"
C        print*,vecw33(1,3),vecw33(2,3),vecw33(3,3),"vw3"
 
 
c on passe le tenseur d'ouverture max dans la base de fissure actuelle
c pour comparer les valeurs max et actuelles
c      vecw33t=transpose(vecw33)
      wmaxa=matmul(transpose(vecw33),wmax33)
      wmaxac=matmul(wmaxa,vecw33)
 
c on récupère les valeurs diag
c c'est cette valeur qu'on va décomposer en multifissure
      do i=1,3
        do j=1,3
        if(i.eq.j) then
        wmaxact(i)=wmaxac(i,j)
        end if
        end do
      end do
 
c recuperation du tenseur des contraintes fibres du pas d avant
c il sert a calculer une eventuelle decharge
      call x6x33(sigf06,sigf033)
      call b3_v33(sigf033,sigf03,vecs33)
c      print*,sigf03(1),sigf03(2),sigf03(3),"contraintes fibres 0"
 
c tenseur d orientation des fibres dans la base d'orientations principales
      if((Eo(1).eq.0.d0).or.(Eo(2).eq.0.d0).or.(Eo(3).eq.0.d0)) then
        print*,"arret du prigramme dans fibs3d.eso"
        print*,"Donnez des valeurs du tenseur d'orientation non nulles"
        stop
        err1=1
      end if
 
c vecteurs principaux du tenseur d'orientation des fibres dans la base du maillage
c on donne 2 vecteurs et on calcule le 3ème
      veo(1,1)=vf1(1)
      veo(2,1)=vf1(2)
      veo(3,1)=vf1(3)
 
      veo(1,2)=vf2(1)
      veo(2,2)=vf2(2)
      veo(3,2)=vf2(3)
 
      veo(1,3)=veo(2,1)*veo(3,2)-veo(3,1)*veo(2,2)
      veo(2,3)=veo(3,1)*veo(1,2)-veo(1,1)*veo(3,2)
      veo(3,3)=veo(1,1)*veo(2,2)-veo(2,1)*veo(1,2)
 
      if(abs(
     #  veo(1,1)*veo(1,2)+veo(2,1)*veo(2,2)+veo(3,1)*veo(3,2)).gt.
     #  precision1) then
        print*,"arret du prigramme dans fibs3d.eso"
        print*,"Donnez vecteurs orthogonaux d'orientation des fibres "
c       read*
        stop
        err1=1
      end if
 
c    normalisation des vecteurs
      do i=1,3
          normv=sqrt(veo(1,i)**2+veo(2,i)**2+veo(3,i)**2)
          veo(1,i)=veo(1,i)/normv
          veo(2,i)=veo(2,i)/normv
          veo(3,i)=veo(3,i)/normv
      end do
c  discrétisation de l'espace en angles solides
      call disp3d(stri,xg,yg,zg,tri,npos)
c traitement des fissures
      do i=1,3
c       initialisation
        loc(i)=.false.
c       wmi=0.d0
        if(i.eq.1) then
            rmax=rmax1
        else if(i.eq.2) then
            rmax=rmax2
        else if(i.eq.3) then
            rmax=rmax3
        end if
 
        if(vw3(i).lt.precision2) then
            sigfr(i)=0.d0
            wmax(i)=0.d0
            rig0t(i)=0.d0
            wmoy(i)=0.d0
            etw(i)=0.d0
            wmi=0.d0
            wmul=0.d0
            wmr=0.d0
            wer=0.d0
            nc=0.d0
            nct=0.d0
 
            if(wmaxact(i).lt.precision2) then
                wmoy(i)=0.d0
                etw(i)=0.d0
                ncf(i)=0.d0
                sigfr(i)=0.d0
                if(nc0(i).gt.0.d0) then
                    ncf(i)=nc0(i)
                else
                    ncf(i)=0.d0
                end if
            else
                ncf(i)=nc0(i)
            end if
 
c           calcul de phimoy
 
            phimoy(i)=0.d0
 
 
            go to 70
        end if
 
        if((wmaxact(i).gt.0.d0).and.(nc0(i).eq.0.d0)) then
            nc0(i)=1.d0
        end if
 
c       vecteurs orientations de fissure
c       vecw(j,i) : orientation de la fissure i (normé? --> oui)
        vecw(1)=vecw33(1,i)
        vecw(2)=vecw33(2,i)
        vecw(3)=vecw33(3,i)
 
 
c       on construit les rayons de l ellipse et les angles une seule fois pour toute
c       c'est les memes pour tous
        do n=1,npos
          do k=1,2*npos-2
c           vecteur sommet des triangles vst(sommet,coordonnée)
c           sommet 1
            vst1(1)=tri(n,k,1,1)
            vst1(2)=tri(n,k,1,2)
            vst1(3)=tri(n,k,1,3)
c             sommet 2
            vst2(1)=tri(n,k,2,1)
            vst2(2)=tri(n,k,2,2)
            vst2(3)=tri(n,k,2,3)
c           sommet 3
            vst3(1)=tri(n,k,3,1)
            vst3(2)=tri(n,k,3,2)
            vst3(3)=tri(n,k,3,3)
c           vecteur dans la direction ei (il est déjà normalisé dans disp3d)
            vei(1)=xg(n,k)
            vei(2)=yg(n,k)
            vei(3)=zg(n,k)
 
c           passage de vei dans la base du maillage puis de fissuration
c           NON vecw est dans la base du maillage, il suffit de passer
c           vei dans la base du maillage
c           veo est la matrice de passage P du maillage a l ellipse
c           ici on veut passer de l'ellipse au maillage donc P A
            vei=matmul(veo,vei)
c           vei=matmul(vecw33,vei)
c           produit scalaire entre l'orientation de la fissure et de ei
            pscal(n,k)=abs(sum(vei*vecw))
c           calcul de l'angle entre la fissure et la direction ei
            phi1=acos(pscal(n,k))
            phi1d(n,k)=180.d0*phi1/pi
c           calculs sont effectues dans la base de l ellipse
c           calcul des angles d'orientation de chaque sommet du triangle n,k
            theta1=acos(vst1(3))
            gamma1=atan2(vst1(2),vst1(1))
            theta2=acos(vst2(3))
            gamma2=atan2(vst2(2),vst2(1))
            theta3=acos(vst3(3))
            gamma3=atan2(vst3(2),vst3(1))
 
c rayon de l'ellipse d'orientation des fibres dans la direction des sommets
      Reiw1 = (Eo(1) ** 2 * Eo(3) ** 2 * sin(theta1)**2*sin(gamma1)** 2
     #- Eo(2) ** 2 * Eo(3) ** 2 * sin(theta1) ** 2*sin(gamma1) ** 2-Eo(1
     #) ** 2 * Eo(2) ** 2 * sin(theta1) ** 2 + Eo(2) ** 2 * Eo(3) ** 2 *
     #sin(theta1) ** 2+ Eo(1) ** 2 * Eo(2) ** 2) ** (-0.1D1 / 0.2D1) * E
     #o(3) * Eo(2) * Eo(1)
 
      Reiw2 = (Eo(1) ** 2 * Eo(3) ** 2 * sin(theta2)**2*sin(gamma2)** 2
     #- Eo(2) ** 2 * Eo(3) ** 2 * sin(theta2) ** 2*sin(gamma2) ** 2-Eo(1
     #) ** 2 * Eo(2) ** 2 * sin(theta2) ** 2 + Eo(2) ** 2 * Eo(3) ** 2 *
     #sin(theta2) ** 2+ Eo(1) ** 2 * Eo(2) ** 2) ** (-0.1D1 / 0.2D1) * E
     #o(3) * Eo(2) * Eo(1)
 
      Reiw3 = (Eo(1) ** 2 * Eo(3) ** 2 * sin(theta3)**2*sin(gamma3)** 2
     #- Eo(2) ** 2 * Eo(3) ** 2 * sin(theta3) ** 2*sin(gamma3) ** 2-Eo(1
     #) ** 2 * Eo(2) ** 2 * sin(theta3) ** 2 + Eo(2) ** 2 * Eo(3) ** 2 *
     #sin(theta3) ** 2+ Eo(1) ** 2 * Eo(2) ** 2) ** (-0.1D1 / 0.2D1) * E
     #o(3) * Eo(2) * Eo(1)
 
            Reiw(n,k)=(Reiw1+Reiw2+Reiw3)/3.d0
 
c           stockage de reiw*omegai
            Romeg(n,k)=Reiw(n,k)*stri(n,k)
            Romeg1(n,k)=Reiw(n,k)*stri(n,k)*pscal(n,k)
c           terme proportionnel à la densité volumique de fibres dans une direction ei
c           car on ne sait pas calculer la somme rei*omegai
c           ratio(n,k)=Reiw(n,k)*stri(n,k)*rhof
c           proportionnel au nombre de fibres du paquet n,k traversant la fissure i
c           nf(n,k)=4*ratio(n,k)*pscal(n,k)/(pi*Df**2)
c           proportionnel a l ancrage moyen des fibres du paquet ei pondéré par le nombre de fibres traversant la fissure
c           Lm(n,k)=(Lf/2-Lf/(4*cos(phi1)))*nf(n,k)
c           l angle moyen reflete les fibres qui traversent effectivement le plan
            phi2d(n,k)=Reiw(n,k)*stri(n,k)*phi1d(n,k)*pscal(n,k)
c           fin de la 1ere boucle sur n,k qui sert a connaitre reiw(n,k) et phi1d(n,k)
          end do
c         fin de boucle sur k
        end do
c        fin de boucle sur n
c     angle moyen
      phimoy(i)=sum(phi2d)/sum(romeg1)
 
      if(vw3(i).gt.wplmax(i)) then
      wplmax(i)=vw3(i)
      end if
 
 
c      if((vw3(i).ge.wmaxact(i)).and.(vw3(i).ge.vw30(i))) then
c      if(vw3(i).ge.wmaxact(i)) then
 
      if(vw3(i).lt.wferm(i)) then
      go to 100
      end if
 
      if((vw3(i).ge.wplmax(i))) then
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c on ouvre la fissure i
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c       modif
c      wpsigt0(i)=wpsigt(i)
 
      rtmin=Rtec*(lech/(0.5d0*sqrt(2*pi)))**(1/mw)
      do r=1,rmax
        if(r.eq.1) then
            Rt(r)=rt00
        else
            Rt(r)=Rtec*(lech*2.d0**(r-1)/le(i))**(1/mw)
        end if
 
        if(rt(r).lt.rtmin) then
        rt(r)=rtmin
        end if
      end do
c   Rtec valeur pour lech quand r=1 on fait les fissures a le/2
 
c   initialisation
      wmi=0.d0
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c ici on test si on etait dans une phase de localisation
c si c est le cas on ne reconstruit pas les courbes on va
c directement dans la phase de localisation
cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
      if(iloc(i).eq.1) then
        go to 80
      end if
 
 
 
      do r=1,rmax
        Rtb=Rt(r)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c construction des courbes contraintes ouvertures de fissures pour chaque rt
c pour 1 fissure fictive incrementale dans la base de fissuration actuelle
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        do l=1,ncalc
 
c           récupération de la contrainte du pas fictif d'avant
            if(l.eq.1) then
                sig0(1)=0.d0
                sig0(2)=0.d0
                sig0(3)=0.d0
                vwfi0(1)=0.d0
                vwfi0(2)=0.d0
                vwfi0(3)=0.d0
            else
                sig0(1)=varft(1,l-1)
                sig0(2)=varft(2,l-1)
                sig0(3)=varft(3,l-1)
                vwfi0(1)=varft(4,l-1)
                vwfi0(2)=varft(5,l-1)
                vwfi0(3)=varft(6,l-1)
            end if
 
c           fissure fictive
            vwfi(i,l)=l*incremwlis
 
c           et ici on recommence donc une boucle ou on va calculer les contraintes
c           jusqu au pic pour tous les Rt
            do n=1,npos
                do k=1,2*npos-2
c               calcul de la force moyenne pour une fibre traversant la fissure i inclinée de phi1d
c               calcul des points caractéristiques de la courbe force ouverture de fissure pour Td et phid
c               rigidité initiale
                if(phi1d(n,k).le.phicrit) then
                    call reconf(cokf1,k0m,Rtb,phi1d(n,k),mink01,maxk01)
                else
                    call reconf(cokf2,k0m,Rtb,phi1d(n,k),mink02,maxk02)
                end if
c               ouverture au pic
                if(phi1d(n,k).le.phicrit) then
                    call reconf(cowp1,wp,Rtb,phi1d(n,k),minwp1,maxwp1)
                else
                    call reconf(cowp2,wp,Rtb,phi1d(n,k),minwp2,maxwp2)
                end if
c               force au pic
                call reconf(cofp,fp,Rtb,phi1d(n,k),minfp,maxfp)
c               force intermédiaire post pic
                f03=0.3d0*fp
c               ouverture pour f=0.3fp
                call reconf(cow03,w03,Rtb,phi1d(n,k),minw03,maxw03)
                if(w03.le.wp) then
                    w03=3.0d0*wp
                end if
c               ouverture finale (f=0)
                call reconf(cowf,wf,Rtb,phi1d(n,k),minwf,maxwf)
                if(wf.le.w03) then
                    wf=3.0d0*w03
                end if
c               calcul de la force dans la fissure i pour les fibres orientées selon ei
                if(vwfi(i,l).le.wp) then
                    B=(k0m*wp-fp)/(fp*wp)
                    f1(n,k)=k0m*vwfi(i,l)/(B*vwfi(i,l)+1.d0)
                else
                 if(vwfi(i,l).lt.wf) then
                    f1(n,k)=fp*exp(log(f03/fp)*(vwfi(i,l)-wp)/(w03-wp))
c                   print*,2,f1(n,k),wp,wf
                 else
                    f1(n,k)=precision1
                 end if
                end if
c               calcul du terme à l'intérieur de la somme de la contrainte
                si0(n,k)=Reiw(n,k)*stri(n,k)*pscal(n,k)*f1(n,k)
            end do
c           fin de boucle sur k
          end do
c        fin de boucle sur n
c        valeurs fictives
         si1=sum(si0)
         sigf(i,r,l)=4.d0*rhof*si1/(pi*Df**2*(sum(Romeg)))
 
 
c        stockage de la rigidité initiale
         if(l.eq.1) then
            rig0(i,r)=sigf(i,r,l)/incremwlis
         end if
            if(rig0(i,r).ne.rig0(i,r).or.(rig0(i,r).eq.0.d0)) then
 
            print*, "Probleme de rigidite initiale",rig0(i,r)
            err1=1
         end if
 
         if((sigf(i,r,l).lt.sig0(i)).or.(vwfi(i,l).gt.(vw3(i)
     # +2.d0*incremwlis))) then
c        on calcule la contrainte jusqu'au pic uniquement ou l ouverture actuelle
         wpsig(i,r)=vwfi0(i)
         sigmax(i,r)=sig0(i)
         go to 10
c       on passe au rt d'apres (on sort de l) pour les rt autres que la premiere fissure
         end if
 
c stockage des variables fictives du pas
         varft(1,l)=sigf(1,r,l)
         varft(2,l)=sigf(2,r,l)
         varft(3,l)=sigf(3,r,l)
         varft(4,l)=vwfi(1,l)
         varft(5,l)=vwfi(2,l)
         varft(6,l)=vwfi(3,l)
         end do
c        fin de boucle sur l
   10    continue
C               if(i.eq.3) then
C         print*,sigf(i,r,l),r,l,rig0(i,r),rt(r),si1
C         end if
        end do
c      fin de boucle sur les Rtb
 
 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c on recommence une boucle de fissures fictives pour avoir la loi totale
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       nct=0.d0
       do l=1,ncalc
c récupération de la contrainte du pas fictif d'avant
       if(l.eq.1) then
       sig0(1)=0.d0
       sig0(2)=0.d0
       sig0(3)=0.d0
       vwfi0(1)=0.d0
       vwfi0(2)=0.d0
       vwfi0(3)=0.d0
       wmi0(1)=0.d0
c       wt(i,l-1)=0.d0
c       sigt(i,l-1)=0.d0
       else
       sig0(1)=varft(1,l-1)
       sig0(2)=varft(2,l-1)
       sig0(3)=varft(3,l-1)
       vwfi0(1)=varft(4,l-1)
       vwfi0(2)=varft(5,l-1)
       vwfi0(3)=varft(6,l-1)
       wmi0(1)=varft(7,l-1)
       end if
c fissure fictive totale cette fois ci
       vwfi(i,l)=l*incremwlis
 
      call endo3f(vwfi(i,l),gf,young,
     # xmt,dtiso,dt3,epict,rt(1),le(i))
 
c       wk1=(Gf-gfap(i))/Rt(1)
c       wk1=(Gf)/Rt(1)
 
      if(l.eq.1) then
      epse13(i)=vwfi(i,l)/le(i)
      else
      epse13(i)=wt(i,l-1)/le(i)
      end if
 
      sigfe(i)=(1.d0-dt3)*
     #rhof*max(min(Ef*epse13(i)*cos(phimoy(i)*
     #pi/180.d0)**2,fu*cos(phimoy(i)*pi/180.d0)**2),-fu)
 
*      sigfe(i)=0.d0
 
       sigb(i)=Rt(1)*(1.d0-dt3)
       sig(i)=sigb(i)+sigf(i,1,l)*dt3+
     #sigfe(i)
 
c      print*,sig(i),sigb(i),sigf(i,1,l),dt3,sigfe(i),"precons"
 
       Sr(i)=lech/exp(log(sig(i)/Rtec)*mw)
 
c on fixe la limite
 
       s1=Rtec*(lech*8.d0/(Lf*cos(pi*phimoy(i)/180.d0)))**(1/mw)
 
 
 
c      if(Sr(i).lt.Lf*cos(pi*phimoy(i)/180.d0)/4) then
       if(sig(i).gt.s1) then
       Sr(i)=Lf*cos(pi*phimoy(i)/180.d0)/4
c       wmiloc(i)=vwfi(i,l)
       sig(i)=s1
       loc(i)=.true.
       end if
 
       if(sig(i).lt.rt(2)) then
 
c      il n y a qu une seule fissure
       wfist(i)=vwfi(i,l)
       nct(1)=1.d0
       wmi(1)=vwfi(i,l)
       wt(i,l)=wfist(i)
       sigt(i,l)=sigf(i,1,l)
       else
c      on multifissure
 
          do r=1,rmax
          if(r.eq.1) then
          nct(r)=1.d0
          wmi(1)=vwfi(i,l)
          else
              if(Sr(i).le.le(i)/(2.d0**(r-1))) then
                if(le(i)/Sr(i).lt.2.d0**(r)) then
                nct(r)=(le(i)/Sr(i)-(2.d0**(r-1)))/2.d0
                else
                nct(r)=2.d0**(r-2)
                end if
              do k=1,ncalc
 
c ici pour trouver les ouvertures de micro fissures on peut aussi equilibrer avec sigf(i,1,l) au lieu de sig(i)
c ca paraitrait peut etre plus logique ?
C               if(sig(i)-sigf(i,r,k).le.0.d0) then
C               if(k.eq.1) then
C        wmi(r)=(sig(i)-(sigf(i,r,k)-((sigf(i,r,k)-0)/
C      #incremwlis)*vwfi(i,k)))/((sigf(i,r,k)-0)/incremwlis)
C               else
C       wmi(r)=(sig(i)-(sigf(i,r,k)-((sigf(i,r,k)-sigf(i,r,k-1))/
C      #incremwlis)*vwfi(i,k)))/((sigf(i,r,k)-sigf(i,r,k-1))/incremwlis)
C               end if
C               go to 20
C               end if
 
              if(sigf(i,1,l)*dt3-sigf(i,r,k).le.0.d0) then
              if(k.eq.1) then
       wmi(r)=(sigf(i,1,l)*dt3-(sigf(i,r,k)-((sigf(i,r,k)-0)/
     #incremwlis)*vwfi(i,k)))/((sigf(i,r,k)-0)/incremwlis)
              else
      wmi(r)=(sigf(i,1,l)*dt3-(sigf(i,r,k)-((sigf(i,r,k)-sigf(i,r,k-1))/
     #incremwlis)*vwfi(i,k)))/((sigf(i,r,k)-sigf(i,r,k-1))/incremwlis)
 
 
              end if
              go to 20
              end if
 
              end do
 
  20          continue
 
C               if(i.eq.3) then
C               print*,wmi(r),r,nct(r),"pre",sr(i),sig(i),rt(r)
C               end if
 
              if(wmi(r).lt.0.d0) then
              print*,wmi(r)
              read*
              end if
              wmi(r)=nct(r)*wmi(r)
              else
c              wmi(r)=0.d0
              go to 30
              end if
        end if
 
 
        end do
c        fin de boucle sur les Rt
   30  continue
c       ici on va calculer des segments de droite pour reconstruire notre courbe
        wfist(i)=sum(wmi)
        wt(i,l)=wfist(i)
        sigt(i,l)=sigf(i,1,l)
 
 
       end if
 
c       fin du si on multifissure
 
c stockage de la rigidité initiale
      if(l.eq.1) then
      rig0t(i)=sigf(i,1,l)/incremwlis
      end if
 
      if(l.gt.1) then
      if(rig0t(i).ne.rig0t(i).or.(rig0t(i).eq.0.d0)) then
      print*, "Probleme de rigidite initiale",rig0t(i)
      err1=1
      end if
      end if
 
 
      if(l.gt.1) then
      if((vwfi(i,l).gt.wpsig(i,1)).or.(loc(i)).or.(wt(i,l-1).gt.vw3(i)
     # )) then
 
c      print*, vwfi(i,l),wpsig(i,1),loc(i),wt(i,l-1),vw3(i),"w",vwfi0(i)
 
c      on calcule la contrainte jusqu'au pic uniquement
      wpsigt(i)=vwfi0(i)
c modif
C       if(wpsigt(i).lt.wpsigt0(i)) then
C       wpsigt(i)=wpsigt0(i)
C       end if
      wmiloc(i)=wmi0(1)
      sigmaxt(i)=sig0(i)
      go to 40
c      on sort de la boucle des deuxiemes fissures fictives l
      end if
      end if
 
c stockage des variables fictives du pas
      varft(1,l)=sigf(1,1,l)
      varft(2,l)=sigf(2,1,l)
      varft(3,l)=sigf(3,1,l)
      varft(4,l)=wfist(1)
      varft(5,l)=wfist(2)
      varft(6,l)=wfist(3)
      varft(7,l)=wmi(1)
 
      end do
c      fin de boucle sur les fissures fictives l
 
   40 continue
c       print*,l,"l1"
cccccccccccccccccccccccccccccccccccccccccccccccc
c on passe au calcul reel ici
cccccccccccccccccccccccccccccccccccccccccccccccc
      if(vw3(i).le.wpsigt(i)) then
ccccccccccccccccccccccccccccccccccccccccccccccccccc
c       on charge
cccccccccccccccccccccccccccccccccccccccccccccccccccc
      wmul=0.d0
      wmr=0.d0
      wer=0.d0
      nc=0.d0
 
c on calcule la contrainte relle dans la fissure i et l'ouverture de la première fissure
c on peut calculer l ouverture w1 car wt a été construit en imposant vwfi dans la premiere fissure
      if(vw3(i).gt.0.d0) then
 
      do l=1,ncalc
        if(wt(i,l)-vw3(i).gt.0.d0) then
           if(l.eq.1) then
            sigfr(i)=(sigt(i,l)-0)/(wt(i,l)-0)*
     #vw3(i)+sigt(i,l)-(sigt(i,l)-0)/(wt(i,l)-0)*
     #wt(i,l)
 
      w1(i)=vwfi(i,l)/wt(i,l)*vw3(i)
           else
 
           sigfr(i)=(sigt(i,l)-sigt(i,l-1))/(wt(i,l)-wt(i,l-1))*
     #vw3(i)+sigt(i,l)-(sigt(i,l)-sigt(i,l-1))/(wt(i,l)-wt(i,l-1))*
     #wt(i,l)
 
      w1(i)=(vwfi(i,l)-vwfi(i,l-1))/(wt(i,l)-wt(i,l-1))*vw3(i)+
     #vwfi(i,l)-(vwfi(i,l)-vwfi(i,l-1))/(wt(i,l)-wt(i,l-1))*wt(i,l)
           end if
           go to 50
        end if
      end do
 
 
 
 
 
      else
      sigfr(i)=0.d0
      sigb(i)=0.d0
      end if
50    continue
c       print*,l,"l2"
 
      call endo3f(w1(i),gf,young,
     # xmt,dtiso,dt3,epict,rt(1),le(i))
 
c calcul de la déformation totale de l element pour prendre en compte lacontrainteelastique des fibres
      call x6x33(epstf6,epse133)
 
      do j=1,3
      vecwj(1)=vecw33(1,j)
      vecwj(2)=vecw33(2,j)
      vecwj(3)=vecw33(3,j)
      epse13(j)=sum(matmul(vecw,epse133)*vecw)
      sigfe(j)=(1.d0-dt3)*
     #rhof*max(min(Ef*epse13(j)*cos(phimoy(j)*
     #pi/180.d0)**2,fu*cos(phimoy(j)*pi/180.d0)**2),-fu)
      end do
 
*       sigfe(i)=0.d0
 
       sigb(i)=Rt(1)*(1.d0-dt3)
       sig(i)=sigb(i)+sigfr(i)*dt3+
     #sigfe(i)
 
c      print*,sig(i),s1,"calc reel"
 
 
c      print*,sig(i),sigb(i),sigfr(i),dt3,sigfe(i)
 
 
       if(sig(i).gt.0.d0) then
         Sr(i)=lech/exp(log(sig(i)/Rtec)*mw)
        s1=Rtec*(lech*8.d0/(Lf*cos(pi*phimoy(i)/180.d0)))**(1/mw)
c        print*,sr(i),Lf*cos(pi*phimoy(i)/180.d0)/4.d0
        if(sig(i).gt.s1) then
         sig(i)=s1
         Sr(i)=Lf*cos(pi*phimoy(i)/180.d0)/4.d0
         wpsigt(i)=vw3(i)
         sigmaxt(i)=min(sigfr(i),sig(i))
c         sigmaxt(i)=sig(i)
c          sigmaxt(i)=sigf03(i)
c         print*,sigfr(3),sig(3),sigfe(i),sigb(i)
 
         end if
       else
       Sr(i)=le(i)
       end if
 
c ici l ouverture de la fissure 1 est calculée de maniuere exacte
c mais pour les autres on suppose que l endo du beton est le meme pour tous
      if(sig(i).lt.rt(2)) then
         nc(1)=1.d0
         wmul(i,1)=vw3(i)
       else
c       on multifissure
        nc(1)=1.d0
        wmul(i,1)=w1(i)
          do r=2,rmax
 
 
              if(Sr(i).le.le(i)/(2.d0**(r-1))) then
                if(le(i)/Sr(i).lt.2.d0**(r)) then
                nc(r)=(le(i)/Sr(i)-(2.d0**(r-1)))/2.d0
                else
                nc(r)=2.d0**(r-2)
                end if
                if(r.eq.1) then
                nc(r)=1.d0
                end if
              do k=1,ncalc
c ici sig(i) est remplacé par sigfr(i)
c ca ne change apparemment rien
                if(sigfr(i)*dt3-sigf(i,r,k).lt.0.d0) then
                 if(k.eq.1) then
       wmul(i,r)=sigfr(i)*dt3*incremwlis/sigf(i,r,k)
                 else
      wmul(i,r)=(sigfr(i)*dt3-(sigf(i,r,k)-((sigf(i,r,k)-sigf(i,r,k-1))/
     #incremwlis)*vwfi(i,k)))/((sigf(i,r,k)-sigf(i,r,k-1))/incremwlis)
c fin remplacement
                 end if
                go to 60
                end if
              end do
  60   continue
 
c              print*,k,"k2",r,"r"
 
 
 
              else
             wmul(i,r)=0.d0
             nc(r)=0.d0
             if(r.eq.1) then
             nc(1)=1.d0
             wmul(i,1)=vw3(i)
             end if
              end if
C         if(i.eq.3) then
C         print*,r,nc(r),sr(i),"calc",sig(i),s1
C         end if
        end do
c        fin de boucle sur les Rt
      end if
c fin du si sig(i).lt.rt(2)
 
 
 
       if(vw3(i).gt.0.d0) then
c calcul de l'ouverture moyenne
       do r=1,rmax
 
       wmr(r)=nc(r)*wmul(i,r)
       end do
c       print*,wmr,nc,"nc",sig(i),rt(2)
c       print*,Sr(i),le(i)/(2.d0**(1-1)),"test1"
c       print*,le(i)/Sr(i),2.d0**(1),"test2"
       wmoy(i)=sum(wmr)/sum(nc)
 
       if(wmoy(i).ne.wmoy(i)) then
       print*,"pb wmr,nc",sum(wmr),sum(nc),rmax,sig(i),rt(2),Sr(i),
     #le(i)
       print*,"elements finis trop petits"
       read*
       end if
c calcul de l'ecart type d ouverture
       do r=1,rmax
        if(nc(r).eq.0.d0) then
        wer(r)=0.d0
        else
        wer(r)=nc(r)*(wmul(i,r)-wmoy(i))**2
        end if
       end do
       Etw(i)=sqrt(sum(wer)/sum(nc))
       cvar(i)=etw(i)/wmoy(i)
       else
       wmoy(i)=0.d0
       etw(i)=0.d0
       wmr=0.d0
       wer=0.d0
       cvar(i)=0.d0
       end if
c ouverture max
       wmax(i)=wmul(i,1)
 
c nombre total de fissures
       ncf(i)=sum(nc)
 
       if(ncf(i).lt.nc0(i)) then
       ncf(i)=nc0(i)
       end if
 
 
 
        if(sig(i).ge.s1) then
         go to 80
         end if
 
 
c sauvegarde de la contrainte de decharge
c      sigdec(i)=sigf03(i)
      sigdec(i)=sigfr(i)
 
 
 
 
 
**  AM
**       else
        endif
ccccccccccccccccccccccccccccccccccccccccc
c       on localise la fissure en ouvrant
ccccccccccccccccccccccccccccccccccccccccc
80     continue
 
**  AM
 
       if( ((VW3(I).LE.WPSIGT(I)).and.(SIG(I).GE.S1)) . OR.
     &     (VW3(I).GT.WPSIGT(I)) ) then
       iloc(i)=1
 
c le nombre fissure ne peut plus augmenter
       ncf(i)=nc0(i)
c seule la fissure 1 continue de s ouvrir les autres ne se dechargent pas...
       wfis(i)=wmiloc(i)+(vw3(i)-wpsigt(i))
C        if(vw3(i).lt.wpsigt(i)) then
C        wfis(i)=vw3(i)
C        end if
 
 
C et donc ici on peut recalculer la contrainte comme provenant d'une mono fissure
c rt a bien ete construit dans tous les cas
       Rtb=Rt(1)
        do n=1,npos
          do k=1,2*npos-2
c calcul de la force moyenne pour une fibre traversant la fissure i inclinée de phi1d
c calcul des points caractéristiques de la courbe force ouverture de fissure pour Td et phid
c rigidité initiale
      if(phi1d(n,k).le.phicrit) then
      call reconf(cokf1,k0m,Rtb,phi1d(n,k),mink01,maxk01)
      else
      call reconf(cokf2,k0m,Rtb,phi1d(n,k),mink02,maxk02)
      end if
c ouverture au pic
      if(phi1d(n,k).le.phicrit) then
      call reconf(cowp1,wp,Rtb,phi1d(n,k),minwp1,maxwp1)
      else
      call reconf(cowp2,wp,Rtb,phi1d(n,k),minwp2,maxwp2)
      end if
c force au pic
      call reconf(cofp,fp,Rtb,phi1d(n,k),minfp,maxfp)
c force intermédiaire post pic
      f03=0.3d0*fp
c ouverture pour f=0.3fp
      call reconf(cow03,w03,Rtb,phi1d(n,k),minw03,maxw03)
        if(w03.le.wp) then
        w03=3.0d0*wp
        end if
c ouverture finale (f=0)
      call reconf(cowf,wf,Rtb,phi1d(n,k),minwf,maxwf)
        if(wf.le.w03) then
        wf=3.0d0*w03
        end if
c calcul de la force dans la fissure i pour les fibres orientées selon ei
       if(wfis(i).le.wp) then
       B=(k0m*wp-fp)/(fp*wp)
       f1(n,k)=k0m*wfis(i)/(B*wfis(i)+1.d0)
       else
         if(wfis(i).lt.wf) then
         f1(n,k)=fp*exp(log(f03/fp)*(wfis(i)-wp)/(w03-wp))
c        print*,2,f1(n,k),wp,wf
         else
         f1(n,k)=precision1
         end if
       end if
c calcul du terme à l'intérieur de la somme de la contrainte
        si0(n,k)=Reiw(n,k)*stri(n,k)*pscal(n,k)*f1(n,k)
 
         end do
c         fin de boucle sur k
        end do
c        fin de boucle sur n
c  calc de la contrainte
      si1=sum(si0)
c      wk1=(Gf-gfap(i))/Rt(1)
c             wk1=(Gf)/Rt(1)
      sigfr(i)=4.d0*rhof*si1/(pi*Df**2*(sum(Romeg)))
 
      if(sigfr(i).gt.sigmaxt(i)) then
      sigfr(i)=sigmaxt(i)
      end if
 
 
       if(wfis(i).gt.0.d0) then
c calcul de l'ouverture moyenne
C        do r=1,rmax
C        nc(r)=ncr0(r)
C          if(r.eq.1) then
C          wmr(1)=wfis(i)
C          else
C          wmr(r)=nc(r)*wmul(i,r)
C          end if
C        end do
C        wmoy(i)=sum(wmr)/(sum(nc))
         dw3(i)= vw3(i)-vw30(i)
 
c calcul de l'ecart type d ouverture
C        do r=1,rmax
C          if(r.eq.1) then
C          wer(1)=(wfis(i)-wmoy(i))**2
C          else
C          wer(r)=nc(r)*(wmul(i,r)-wmoy(i))**2
C          end if
C        end do
C        Etw(i)=sqrt(sum(wer)/sum(nc))
         if(nc0(i).eq.1.d0) then
         etw(i)=0.d0
         wmoy(i)=vw3(i)
         else
c         wmoy(i)=wmoy0(i)+dw3(i)/nc0(i)
         wmoy(i)=wmoy0(i)
C on conserve le coefficient de variation lors de la localisation
         etw(i)=etw0(i)
         end if
 
       else
       wmoy(i)=0.d0
       etw(i)=0.d0
       wmr=0.d0
       wer=0.d0
       end if
c ouverture max
       wmax(i)=wfis(i)
c      print*,wmax(i),i,"loc",vw3(i),vw30(i),wmaxact(i),rmax,le(i)
c      print*, wpsigt(i),sigmaxt(i),sum(Romeg),sigfr(i),wfis(i),nc0(i),
c     # wmiloc(i)
c       print*,"loc",dw3(i),vw30(i),etw(i),etw0(i),wmoy(i),i
 
 
c sauvegarde de la contrainte pour decharge
c      sigdec(i)=sigf03(i)
      sigdec(i)=sigfr(i)
 
c      print*,sigdec(i),i,vw3(i),vw30(i),wmaxact(i)
 
      if(vw3(i).le.wferm(i)) then
      go to 100
      end if
 
C       if(sigfr(i).lt.0.d0) then
C       print*,sigfr(i),vw3(i),vw30(i),wmaxact(i),wfis(i),
C      #wmiloc(i),vw3(i),wpsigt(i)
C       read*
C       end if
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       end if
c       fin calcul avec w qui augmente
ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
 
**  AM
**      else
        endif
 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      on a vw3(i).lt.wmaxact(i) donc on a fermé la fissure on est dans une phase décharge
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
100   continue
**  AM
              IF(    (VW3(I).LT.WFERM(I)) .OR.
     &             (VW3(I).LT.WPLMAX(I))) THEN
 
 
C       if(vw30(i).ge.vw3(i)) then
C       sigfr(i)=sigf03(i)-rig0t(i)*(vw30(i)-vw3(i))
C       else
C       print*,"comme par hasard",vw30(i),vw3(i)
C       read*
C       if(sigf03(i)-rig0t(i)*(vw30(i)-vw3(i)).gt.sigf03(i)) then
C       sigfr(i)=sigf03(i)
C       else
C       sigfr(i)=sigf03(i)-rig0t(i)*(vw30(i)-vw3(i))
C       end if
C       end if
 
c modif
c      sigfr(i)=sigf03(i)-rig0t(i)*(wmaxact(i)-vw3(i))
 
 
      if(vw30(i).ge.vw3(i)) then
      sigfr(i)=sigf03(i)-rig0t(i)*(vw30(i)-vw3(i))
c      else if(wmaxact(i).ge.vw3(i)) then
c      sigfr(i)=sigf03(i)-rig0t(i)*(wmaxact(i)-vw3(i))
      else if(sigf03(i)-rig0t(i)*(vw30(i)-vw3(i)).ge.sigf03(i)) then
      sigfr(i)=sigf03(i)-rig0t(i)*(vw30(i)-vw3(i))
c      else if(sigf03(i)-rig0t(i)*(wmaxact(i)-vw3(i)).ge.sigf03(i)) then
c      sigfr(i)=sigf03(i)-rig0t(i)*(wmaxact(i)-vw3(i))
      end if
 
c      print*,sigfr(i),sigf03(i),sigdec(i)
 
 
      if((sigfr(i).lt.0.d0).and.(sigf03(i).gt.0.d0)) then
      sigfr(i)=0.d0
      if(vw3(i).gt.wferm(i)) then
      wferm(i)=vw3(i)
      end if
      end if
 
      if(sigfr(i).lt.0.d0) then
      sigfr(i)=0.d0
      else if(sigfr(i).gt.sigdec(i)) then
      sigfr(i)=sigdec(i)
      end if
 
      if(vw3(i).lt.wferm(i)) then
      sigfr(i)=0.d0
      end if
 
 
      ncf(i)=nc0(i)
      dw3(i)= vw3(i)-vw30(i)
 
 
 
      if(nc0(i).eq.1.d0) then
      etw(i)=0.d0
      wmoy(i)=vw3(i)
      wmax(i)=vw3(i)
      else
      wmoy(i)=wmoy0(i)+dw3(i)/nc0(i)
c      Etw(i)=etw0(i)/wmoy0(i)*wmoy(i)
      Etw(i)=cvar(i)*wmoy(i)
      end if
      wmax(i)=wmax0(i)+dw3(i)/nc0(i)
      if(wmax(i).lt.0.d0) then
      wmax(i)=0.d0
      end if
 
      end if
 
70    continue
 
       if(ncf(i).lt.nc0(i)) then
       print*,"le nombre de fissures a diminue"
       print*,ncf(i),nc0(i),i,vw3(i),vw30(i),wmaxact(i)
       stop
       end if
 
      if(wmax(i).ge.wmaxd0(i)) then
      wmaxd(i)=wmax(i)
      else
      wmaxd(i)=wmaxd0(i)
      end if
      end do
c      fin de boucle sur les fissures
 
C       print*,sigfr(1),sigfr(2),sigfr(3),"sig"
C       print*,vw3(1),vw3(2),vw3(3),"vw"
C       print*,ncf(1),ncf(2),ncf(3),"nc"
C       print*,etw(1),etw(2),etw(3),"etw"
 
      if(etw(1).ne.etw(1).or.etw(2).ne.etw(2).or.etw(3).ne.etw(3)) then
      print*,"probleme ecart type fissure" ,etw(1),etw(2),etw(3)
      print*,Etw0(1),vw3(1),dw3(1),wmoy(1),nc0(1),wmaxact(1),wmoy0(1),
     #ncf(1),iloc(1)
      print*,Etw0(2),vw3(2),dw3(2),wmoy(2),nc0(2),wmaxact(2),wmoy0(2),
     #ncf(2),iloc(2)
      print*,Etw0(3),vw3(3),dw3(3),wmoy(3),nc0(3),wmaxact(3),wmoy0(3),
     #ncf(3),iloc(3)
      read*
      end if
 
C       if(nc(3).gt.16.d0) then
C       print*,'pb'
C       read*
C       end if
 
c  reconstruction du tenseur des contraintesde fibres avec endo dans la base fixe
      do i=1,3
        do j=1,3
        if(i.eq.j) then
        sigfp33(i,j)=sigfr(i)
        else
        sigfp33(i,j)=0.d0
        end if
        end do
      end do
c      sigfx33=matmul(matmul(transpose(vecw33),sigfp33),vecw33)
c     vecw33 est la matrice de passage de la base maillage a fissuration
c     ici on veut faire l inverse cad : P A tP
      sigfx33=matmul(matmul(vecw33,sigfp33),transpose(vecw33))
      call x33x6(sigfx33,sigfx6)
 
C       print*,"wmax",wmax
C       print*,"wmoy",wmoy
C       print*,"etw",etw
C       print*,"ncf",ncf
C       if(sum(sigfx6*sigfx6).gt.0.d0) then
C        print*,"sigfx6",sigfx6
C        end if
 
C          print*,"sig(3)",sig(3),"sigfr",sigfr(3),s1
C      #,"dt3","wmaxd",wmaxd(3),dt3
      return
 
      end
 
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales