Télécharger prpfib.eso

Retour à la liste

prpfib
C PRPFIB    SOURCE    FD218221  24/02/07    21:15:23     11834          
      subroutine prpfib(p,Lf,mu,
     # Td,Tmax,Ef,Hf,Df,sk,F0,alphad,M0,Rt,fu,
     # L0,d,e,fy,mwor,nligne,niter,nangl,nrt,nlongfi,
     # nangredu,phit,phicrit)    
 
 
c      implicit none
        implicit real*8 (a-h,o-z)
        implicit integer (i-n)     
 
 
c   ***************************************************************
c   ************************ MODIFICATIONS ************************
C   ***************************************************************
c        variables pour la reduction de modele
 
c   ***************************************************************
c   ***************** NOMS DES VARIABLES/PARAMETRES ***************
c   ***************************************************************
 
c   i : boucle sur le nombre ligne
c   j : boucle sur les angles
c   k : boucle de sous-itération
c   l : boucle sur les longueurs ancrées
c   nligne : nombre de lignes pour 1 fibre
c   niter : nombre d'itérations max sur k
c   incremf : incrément de force élastique (MN) 
c   incremld : incrément de longueur décollée (m)
c   incremsp : incrément d'extraction de fibre (m)
c   fel : force élastique imposée (MN)
c   ld1 : longueur décollée coté 1 (m)
c   ld2 : longueur décollée coté 2 (m)
c   sp1 : extraction du coté le moins ancré de la fibre (m)
c   fcrit1 : force de début de décollement coté moins ancré (MN)
c   fcrit2 : force de début de décollement coté plus ancré (MN)
c   L1 : longueur ancrée variable par ecaill jusqu'au pic de force coté le moins ancré (m)
c   L2 : longueur ancrée variable par ecaill jusqu'au pic de force coté le plus ancré (m)
c   dls1 : incrément d'écaillage coté le moins ancré  (m)
c   dls2 : incrément d'écaillage coté le plus ancré  (m)
c   ls1 : longueur d'écaillage coté le moins ancré (m)
c   ls2 : longueur d'écaillage coté le plus ancré (m)  
c   F : force pour une fibre droite (MN)
c   w1 : glissement de décollement ou élastique lié au coté 1 (m) 
c   w2 : glissement de décollement ou élastique lié au coté 2 (m) 
c   Ef : module d'young de la fibre (MPa)
c   Df : diamètre de la fibre (m)
c   Af : section de la fibre (m^2)
c   Hf : module de rigidité de l'interface (MPa/m)
c   kc : coefficient de diffusion de l'equation differentielle des fibres (1/m)
c   Lf : longueur totale de la fibre (m)
c   Tmax : contrainte admissible maximale par l'interface fibre-matrice (MPa)
c   Td : contrainte de frottement fibre-matrice après décohésion (MPa)
c   sk : glissement caractéristique, pilote la pente de l'extraction (m)
c   F0 : force d'about de fibre (MN)
c   alpha : angle d'ouverture du tétraèdre d'écaillage et du cone de rupture du béton en rad
c   Rt : résistance à la traction de la matrice (MPa)
c   mu : coefficient de frottement fibre matrice pour fibres orientées
c   phi : angle d'inclinaison des fibres // normale à la fissure en rad
c   Fd1 : force de décollement côté le moins ancré (MN)
c   sd1 : glissement de décollement côté le moins ancré (m)
c   sd2 : glissement de décollement côté le plus ancré (m)
c   La1 : ancrage initial côté le moins ancré (m)
c   La2 : ancrage initial côté le plus ancré (m)
c   lsf1 : longueur d'écaillage sauvegardé côté le moins ancré (m)
c   lsf2 : longueur d'écaillage sauvegardé côté le plus ancré (m)
c   sel1 : glissement élastique côté le moins ancré (m)
c   sel2 : glissement élastique côté le plus ancré (m)
c   ldind1 : sert à indiquer si un décollement s'est déjà produit sur le côté le moins ancré (m)
c   ldind2 : sert à indiquer si un décollement s'est déjà produit sur le côté le plus ancré (m)
c   phid : valeur angle inclinaison fibre // normale à la fissure en °
c   alphad : valeur angle d'ouverture du tétraèdre d'écaillage et du cone de rupture du béton en °
c   sdf1 : glissement de décollement sauvegardé côté le moins ancré (m)
c   sdf2 : glissement de décollement sauvegardé côté le plus ancré (m)
c   Ffo : force de la fibre orientée (MN)
c   Lf1 : ancrage de la fibre sauvegardé côté le moins ancré (m)
c   Lf2 : ancrage de la fibre sauvegardé côté le plus ancré (m)
c   ldf1 : longueur décollée de la fibre sauvegardée côté le moins ancré (m)
c   ldf2 : longueur décollée de la fibre sauvegardée côté le plus ancré (m)
c   wf1 : ouverture de fissure totale sauvegardée (m)
c   Ff : force sauvegardée avant orientation (MN)
c   lsk1,lsk2 : paramètre servant à calculer la pente de la droite d'écaillage 
c   Fco : force admissible sur la surface du cone de rupture final (MN)
c   Fmax : force maximale admisssible par la fibre (MN)
c   fy : limite élastique de la fibre (MPa)
c   rupt : indicateur de rupture de la fibre ou du cone, ou d'extraction  totale passe à true si il y a rupture
c   mi1,mi2 : positions des valeurs de la courbe f(w) d'une fibre encradrant la valeur de wlist lors du calcul de la force moy
c   n : Numéro de ligne de la courbe force moyenne ouverture de fissure
c   nangl +1 : nombre de discrétisation de l'angle
c   nlongfi : nombre de discrétisation de la longueur ancrée des fibres
c   nwlis : nombre d'ouvertures de fissures pour l'intrpoation pour calcul des forces moyennes
c   nombre d'ouvertures de fissures lors du calul de Fmoy
c   vark : valeurs de k lorsque l'on sort de la boucle d'équilibre d'écaillage
c   nlign1 : valeur optimisée du nombre de lignes i pour le calcul des forces
c   incremwlis : incrément d'ouverture de fissure pour le calcul des forces moyennes (m)
c   precision1 : sert à ne pas avoir 0 là ou on veut pas
c   npic : position de wlist pour laquelle on a la force moyenne au pic
c   n0longfi : nombre de longueur de fibre pour calcul de la force moyenne ne prenant pas en compte les longueurs < 1mm
c   M0 : module d'écrouissage liée à l'abrasion (MPa)
c   L0 : longueur caractéristique prenant en compte l'influence de la longueur ancrée sur l'abrasion (m)
c   somFL : somme des forces pour un angle et tous les ancrages pour calcul fmoy 
c   ffo1 : tableau des forces orientées dépendant de i,phi,L (MN)
c   ffo2 : tableau des forces orientées dépendant de n,phi,L (MN)
c   FmoyL : Force moyenne pour une orientation (MN)
c   wlist : liste d'ouverture de fissure de la force moyenne (m)
c   Dwf : différence entre l'ouverture de wlist et de wf1 servant a intrpoer
c   fu : contrainte ultime pour les fibres (MPa)
c   wint0,wint1 : valeurs de wf1 qui encadrent la valeur de wlist pour intrpoation (m)
c   fint0,fint1 : valeurs de ffo1 qui correspondent à wint0 et wint1 pour intrpoation (MN)
c   fres : valeur intrpoée de la force pour calcul de la force moyenne (MN)
c   k0 : rigidité initiale pour une fibre servant a calculer s0 (MN/m)
c   fk0 : valeur de force au pas 2 servant a calculer k0 (MN)
c   sk0 : glissement correspondant à fk0 (m)
c   Fpic : force moyenne au pic pour chaque cas (MN)
c   k0m : rigidité initiale moyenne au pic (MN/m) 
c   fcrit0 : force de début de décollement pour la longueur ancrée initiale pour optimisation du nombre de pas phase élastique (MN)
c   nloc : localisation de la force moyenne maximale
c   Fmoylis : liste de force moyenne 1d servant à localiser le pic (MN)
c   wpic : ouverture de fissure au pic de force moyenne (m)
c   wfin : ouverture de fissure pour laquelle Fmoy vaut 0 pour la première fois après le 1er pas (m)
c   spind : indicateur de mise en extraction
c   rotule : mettre a false pour ne pas prendre en compte la rotule d'extraction des fibres orientées
c   finf : indicateur marquant wfin
c   p : numéro de la ligne des listes (Td,phid,point caractéristique des courbes)
c   lisTd : Liste des Td pour lesquels les forces moyennes ont été calculées (MPa)
c   lisphid : Liste des phid pour lesquels les forces moyennes ont été calculées (MPa)
c   lisfp : Liste des forces moyennes au pic correspondant a Td phid (MN)
c   lisk0m : Liste des rigidités initiales moyennes correspondant a Td phid (MN/m)
c   liswf : Liste des ouvertures de fissures moyennes finales (Fmoy=0) correspondant a Td phid (m)
c   lisw03 : Liste des ouvertures de fissures moyennes pour F=0,3*fpic correspondant a Td phid (m)
c   f03 : 0.3*fpic (MN)
 
 
c        parametres
         integer i,j,k,l,t,mi1,mi2,n,nligne,niter,nangl,nlongfi,nwlis
         integer vark,nlign1,npic,n0longfi,nrt,nangredu,iind
         real*8 incremf,incremld,incremsp,Lf,incremwlis,pi,precision1 
         real*8 phit,phicrit
 
         parameter (incremwlis=5.0d-6)
         parameter (PI=3.141592653589793D0,precision1=1.0d-10)
 
           real*8 M0,L0
           real*8 Ef,Df,Af,Hf,kc,Tmax,Td,sk,F0,alpha,Rt,mu,phi
           real*8 La1,La2
           real*8 phid,alphad
           real*8 lsk1,lsk2
           real*8 Fmax
           real*8 somFL,fy,s0,nligner
           real*8 Fr
           real*8 fu,Fpasav
           real*8 wint0,wint1,fint0,fint1,fres,k0,fk0,sk0,scrit
           real*8 fcrit0,wpic       
 
C cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
           logical rupt,spind,rotule,finf,ldind1,ldind2,ecail1  
           parameter(rotule=.true.,ecail1=.false.)
C c        variables pour la reduction de modele
           real*8 f03,w03,Fpic,k0m,wfin
           integer p,d,e
*
       SEGMENT MWOR
           real*8 fel(mligne),ld1(mligne,miter),ld2(mligne,miter)
           real*8 sp1(mligne),fcrit1(mligne,miter),fcrit2(mligne,miter)        
           real*8 L1(mligne,miter),L2(mligne,miter),dls1(mligne,miter)
           real*8 dls2(mligne,miter),ls1(mligne,miter),ls2(mligne,miter)
           real*8 F(mligne,miter),w1(mligne,miter),w2(mligne,miter)
           real*8 Fd1(mligne,miter),sd1(mligne,miter),lsf1(mligne)
           real*8 sel1(mligne,miter),sdf1(mligne)
           real*8 Fo(mligne,miter),ldf2(mligne),Ffo(mligne)
           real*8 lsf2(mligne),Lf2(mligne),sdf2(mligne),Fco(mligne)
           real*8 sd2(mligne,miter),sel2(mligne,miter)
           real*8 Dwf(mligne),Ff(mligne),Lf1(mligne),ldf1(mligne)   
           real*8 ffo2(mligne,mangl1,mlongfi)
           real*8 Ffo1(mligne,mangl1,mlongfi)
           real*8 FmoyL(mligne,mangl1),wf1(mligne,mangl1,mlongfi)
           real*8 Fmoylis(mligne),wlist(mligne)                                 
*                                       
          real*8 lisRt(mrtang1),lisphid(mrtang1)
          real*8 lisfp(mrtang1),lisk0m(mrtang1)
          real*8 liswf(mrtang1)
          real*8 lisw03(mrtang1)                
*               
          real*8 lisrt1(mangnrt),lisphid1(mangnrt)
          real*8 lisrt2(mangnan)
          real*8 lisphid2(mangnan)  
          real*8 liswp1(mangnrt),liswp2(mangnan)     
          real*8 lisk01(mangnrt),lisk02(mangnan)                
       ENDSEGMENT
*
 
C          print*,Ef,"Ef"
C          print*,Lf,"Lf"
C          print*,Df,"df"
C          print*,Hf,"hf"
C          print*,Tmax,"tmax"
C          print*,td,"td"
C          print*,sk,"sk"
C          print*,F0,"f0"
C          print*,alphad,"alphad"
C          print*,Rt,"Rt"
C          print*,mu,"mu"
C          print*,M0,"M0"
C          print*,fu,"fu"
C          print*,L0,"L0"
C          print*,fy,"fy"
 
         nwlis=int(Lf/(incremwlis))      
 
c        initialisation de variables
         ldf1(1)=0.d0
         lsf1(1)=0.d0 
         sp1(1)=0.d0
         sdf1(1)=0.d0
 
 
         ldf2(1)=0.d0
         lsf2(1)=0.d0 
         sdf2(1)=0.d0                
 
         ldind1=.false.
         ldind2=.false.
 
c parametres increments 
          incremld=Lf/200.d0         
          incremsp=Lf/1000.d0                
 
 
c section de la fibre
        Af=(PI*Df**2)/4
c coefficient de l exponentielle        
        kc=sqrt(4.d0*Hf/(Ef*Df))
c conversion de l'angle du cone en radian
        alpha=PI*alphad/180
 
c     boucle sur les angles
      do j=0,nangl
      if(j.lt.nangredu) then
      phid=phicrit*j/(nangredu-1)
      else
c      phid=phit*j/(nangl)
      if(j.eq.nangredu) then
      phid=phicrit
      else
      phid=phit*j/(nangl)
      end if
      end if
 
 
C       if(j.eq.0) then
C       ffo2=0.d0
C       wlist=0.d0
C       FmoyL=0.d0
C       Fmoylis=0.d0
C       Ffo1=0.d0
C       wf1=0.d0
C       wlist=0.d0
C       end if
 
 
c      print*,j,phid,nangredu
c      read*
 
c      phid=15.d0*j
      phi=PI*phid/180  
      finf=.false.      
c        longueurs des ancrages
c       boucle sur les longueurs ancrées
c on exclue les fibres d ancrage inf a x mm --> non
c int(2.d0*x*nlongfi/Lf)+1
c        n0longfi=int(2.d0*0.0d-3*nlongfi/Lf)+1
        n0longfi=1        
c        print*,n0longfi
        do l=n0longfi,nlongfi
        La1=(Lf/2.d0)*l/nlongfi
        La2=Lf-La1
        Lf1(1)=La1
        Lf2(1)=La2                  
        spind=.false.
 
        lsk1=La1/5.d0        
        lsk2=La2/5.d0        
c     ici on va optimiser le nombre de lignes pour chaque longueur ancrée
c         calcul de la force critique pour L
          fcrit0=Tmax*PI*Df*(exp(2.d0*kc*La1)-1.d0)/
     #    (kc*(exp(2.d0*kc*La1)+1.d0))
          incremf=fcrit0/10.d0
c glissement de debut de decollement
          scrit=Tmax/Hf
c         calcul du nombre ligne nécessaires
          nligner=fcrit0/incremf+la1/incremld+la1/incremsp
          nlign1=int(nligner)
c          print*,nlign1
 
          do i=1,nlign1
 
             if(i.eq.1) then
             rupt=.false.
             s0=0.d0
             end if
 
             if(i.eq.1) then
             fel(i)=0.d0
             else
             fel(i)=fel(i-1)+incremf
             end if
 
             do k=1,niter
c il suffit d'un seul ecaill plus grand que l'ancrage pour casser              
c            longueur du coté le moins ancré             
C                  if(i.eq.1) then
C                      L1(i,k)=La1
C                  else
C                      if(k.eq.1) then
C                      L1(i,k)=Lf1(i-1)
C                      else
C                        if ((Lf1(i-1).lt.dls1(i,k-1)).or.(Lf1(i-1) 
C      #                 .eq.0.d0).or.(rupt)) then
C                        L1(i,k)=0.d0
C                        rupt=.true.
C c                       exit
C                        go to 10
C                        else
C                        L1(i,k)=Lf1(i-1)-dls1(i,k-1)
C                        end if
C                      end if
C                  end if
c            longueur du coté le mois ancré    
 
 
c             print*,tmax ici ca marche          
                 if(i.eq.1) then
                     if(k.eq.1) then
                     L1(i,k)=La1
                     else
                       if((La1.le.ls1(i,k-1)).or.(rupt)) then
                       L1(i,k)=0.d0
                       rupt=.true.
                       vark=k
c                       exit
                       go to 10
                       else
                       L1(i,k)=La1-dls1(i,k-1)
                       end if
                     end if
                 else
                     if(k.eq.1) then
                     L1(i,k)=Lf1(i-1)
                     else
                       if((Lf1(i-1).le.dls1(i,k-1)).or.(rupt)
     #                 .or.(Lf1(i-1).eq.0.d0)) then
                       L1(i,k)=0.d0
                       rupt=.true.
                       vark=k
c                       exit
                       go to 10
                       else
                       L1(i,k)=Lf1(i-1)-dls1(i,k-1)
                       end if
                     end if
                 end if                 
c            longueur du coté le plus ancré             
                 if(i.eq.1) then
                     if(k.eq.1) then
                     L2(i,k)=La2
                     else
                       if((La2.le.ls2(i,k-1)).or.(rupt)) then
                       L2(i,k)=0.d0
                       rupt=.true.
                       vark=k
c                       exit
                       go to 10
                       else
                       L2(i,k)=La2-dls2(i,k-1)
                       end if
                     end if
                 else
                     if(k.eq.1) then
                     L2(i,k)=Lf2(i-1)
                     else
                       if((Lf2(i-1).le.dls2(i,k-1)).or.(rupt)
     #                 .or.(Lf2(i-1).eq.0.d0)) then
                       L2(i,k)=0.d0
                       rupt=.true.
                       vark=k
c                       exit
                       go to 10
                       else
                       L2(i,k)=Lf2(i-1)-dls2(i,k-1)
                       end if
                     end if
                 end if
 
c               print*,tmax
 
c            calcul de la force critique coté moins ancré
             fcrit1(i,k)=Tmax*PI*Df*(exp(2.d0*kc*L1(i,k))-1.d0)/
     #       (kc*(exp(2.d0*kc*L1(i,k))+1.d0))
 
 
 
 
c            calcul de la force critique coté plus ancré
             fcrit2(i,k)=Tmax*PI*Df*(exp(2.d0*kc*L2(i,k))-1.d0)/
     #       (kc*(exp(2.d0*kc*L2(i,k))+1.d0))
 
 
c********************************************                 
c            coté le moins ancré
c********************************************
             if(i.eq.1) then
             ldind1=.false.
             ldind2=.false.
             end if
 
             if((fel(i).ge.fcrit1(i,k)).or.(ldind1)) then  
             ldind1=.true.
c            phase de décollement
c            initialisation ld1
                 if(i.eq.1) then 
                 ld1(i,k)=0.d0
                 else
                     if(k.eq.1) then
                     ld1(i,k)=ldf1(i-1)+incremld
                     else
                       if(ldf1(i-1)+incremld-dls1(i,k-1).le.0.d0) 
     #                 then
C                        print*,"decolllement  a 0",ldf1(i-1),dls1(i,k-1)
C      #                  ,i,k
                       ld1(i,k)=0.d0
                       else
                       ld1(i,k)=ldf1(i-1)+incremld-dls1(i,k-1)
                       end if       
                     end if                       
                 end if
c                 print*,ld1(i,k),i,k
c            calcul de la force de decollement
             Fd1(i,k)=PI*Df*Td*ld1(i,k)+Tmax*PI*Df*
     #       (exp(2.d0*kc*L1(i,k))-exp(2.d0*kc*ld1(i,k)))/
     #       (kc*(exp(2.d0*kc*L1(i,k))+exp(2.d0*kc*ld1(i,k))))
 
c            calcul du glissement de decollement
             sd1(i,k)=Tmax/Hf+((PI*Df*Td*ld1(i,k)**2)/(2.d0*Ef*Af))+
     #       (Tmax*PI*Df*ld1(i,k)/(kc*Af*Ef))*((exp(2.d0*kc*L1(i,k))-
     #       exp(2.d0*kc*ld1(i,k)))/(exp(2.d0*kc*L1(i,k))+
     #       exp(2.d0*kc*ld1(i,k))))
c                      print*,sdf1(i-1),sd1(i,k),l1(i,k),ld1(i,k)     
                 sp1(i)=0.d0  
c                force glissement sauvegardés  
                 if(i.eq.1) then 
                 F(i,k)=0.d0
                 w1(i,k)=0.d0
c                initialisation de variables
                 sdf1(i)=0.d0
                 sd1(1,k)=0.d0 
                 sp1(i)=0.d0    
                 else 
c                    extraction du coté le moins ancré
c                    seulement si le glissement de décollement reduit ET Ld augmente
                     if((Ld1(i,k).ge.L1(i,k)).or.
     #               ((sd1(i,k)-sdf1(i-1).le.0.d0)
     #                 .and.(ld1(i,k).gt.ldf1(i-1))).or.spind) then
                     s0=f0/k0
                     ld1(i,k)=ldf1(i-1)                     
 
                     if(.not.spind) then
                     iind=i
                     end if
 
                     if((spind).and.(i.gt.iind)) then
                     sp1(i)=sp1(i-1)+incremsp
                     else
                     sp1(i)=0.d0
                     end if
                     sd1(i,k)=sdf1(i-1)
                     w1(i,k)=sdf1(i-1) 
 
c                    si on commence a extraire une fois on extrait toujours apres                     
                     spind=.true.
 
                       if(sp1(i).gt.L1(i,k)) then 
                       F(i,k)=0.d0
                       sp1(i)=sp1(i-1)
                       ld1(i,k)=ldf1(i-1)
                       rupt=.true.
                       else
                       F(i,k)=(Td*(sk/(sk+sp1(i)))+M0*(sp1(i))                     
     #                 *L1(i,k)/(L1(i,k)+L0)**2)*PI*Df*(L1(i,k)-
     #                 sp1(i))+F0     
c                       print*,f(i,k),s0,k0,sp1(i),ld1(i,k),sd1(i,k)
                       ld1(i,k)=ldf1(i-1)
                       end if
                     else 
c                    decollement  du coté le moins ancré   
                     sp1(i)=0.d0               
                     F(i,k)=Fd1(i,k)
                     w1(i,k)=sd1(i,k)
 
                     end if
 
                 end if
 
             else 
c            glissement elastique 
               if(L1(i,k).eq.0.d0) then
               rupt=.true.
               else
               sel1(i,k)=fel(i)*kc*(exp(kc*2.d0*L1(i,k))+1.d0)/
     #         (Hf*PI*Df*(exp(2.d0*kc*L1(i,k))-1.d0))
               end if
c            la longueur decollee vaut 0
             Ld1(i,k)=0.d0
 
 
c            force glissement sauvegardés 
                if(i.eq.1) then
                F(i,k)=0.d0
                w1(i,k)=0.d0
                else
                F(i,k)=fel(i)
                w1(i,k)=sel1(i,k)
                end if
c            initialisation de variables     
                sdf1(i)=0.d0
                sd1(i,k)=0.d0 
                sp1(i)=0.d0  
 
 
             end if
 
c**********************************************                 
c            coté le plus ancré
c**********************************************
             if((F(i,k).ge.fcrit2(i,k)).or.(ldind2)) then  
             ldind2=.true.
c            phase de décollement
c            initialisation ld2
                 if(i.eq.1) then 
                 ld2(i,k)=0.d0
                 else
c                si la force décroit la longueur décollée ld2 stagne
c                sauf si les deux longueurs sont égales
c                si la force croit à la fin du décollement le décollement s arrete quand meme
                  if(((F(i,k)-Ff(i-1).gt.0.d0).or.(l.eq.nlongfi))
     #            .and.(sp1(i).eq.0.d0)) then
                     if(k.eq.1) then
                     ld2(i,k)=ldf2(i-1)+incremld
                     else
                       if(ldf2(i-1)+incremld-dls2(i,k-1).le.0.d0) 
     #                 then
                       ld2(i,k)=0.d0
                       else
                       ld2(i,k)=ldf2(i-1)+incremld-dls2(i,k-1)
                       end if       
                     end if  
                  else
                  ld2(i,k)=ldf2(i-1)  
                  end if                  
                 end if
 
c            la force est toujours celle du coté 1
c            calcul du glissement de decollement
             sd2(i,k)=Tmax/Hf+((PI*Df*Td*ld2(i,k)**2.d0)/(2.d0*Ef*Af))+
     #       (Tmax*PI*Df*ld2(i,k)/(kc*Af*Ef))*(exp(2.d0*kc*L2(i,k))-
     #       exp(2.d0*kc*ld2(i,k)))/(exp(2.d0*kc*L2(i,k))+
     #       exp(2.d0*kc*ld2(i,k)))
 
c                glissement sauvegardés  
                 if(i.eq.1) then 
                 w2(i,k)=0.d0
c                initialisation de variables
                 sdf2(i)=0.d0
                 sd2(i,k)=0.d0                 
                 else 
c                decollement  du coté le plus ancré         
                   if(spind) then
                   sd2(i,k)=sdf2(i-1)
                   w2(i,k)=sdf2(i-1)
                   else                  
                   w2(i,k)=sd2(i,k)
                   end if
                 end if
 
             else 
 
c            glissement elastique 
             sel2(i,k)=F(i,k)*kc*(exp(kc*2.d0*L2(i,k))+1.d0)/
     #       (Hf*PI*Df*(exp(2.d0*kc*L2(i,k))-1.d0))
c            la longueur decollee vaut 0
             Ld2(i,k)=0.d0
c            force glissement sauvegardés 
                  if(i.eq.1) then
                  w2(i,k)=0.d0
                  else
                  w2(i,k)=sel2(i,k)
                  end if
c            initialisation de variables     
             sdf2(i)=0.d0
             sd2(i,k)=0.d0
 
             end if 
 
 
ccccccccccccccccccccccccccccccccccccccccccccccc
c rotule plastique pour les fibres orientées
ccccccccccccccccccccccccccccccccccccccccccccc
C             if(rotule.and.spind) then
C               if(phi.gt.0.d0) then
C             F(i,k)=F(i,k)+fy*Df**2/(6.d0*sin(phi))
C               end if
C             end if
 
C        if(rotule.and.spind) then
C          if(i.gt.1) then
C            if(k.gt.1) then
C c           Fp=2.d0*fy*Df**3/(6.d0*(ls1(i,k-1)+w1(i,k)
C c     #     +sp1(i)))
C c       Fo(i,k)=(F(i,k)*(1.d0+mu*sin(phi/2.d0))+Fp)      
C c     #  /(1.d0-mu*sin(phi/2.d0))  
C         Mp=fy*Df**3/6.d0
C         Fo(i,k)=(F(i,k)*L1(i,k)-6*Mp*mu)/((-3*sin(phi)*mu+cos(phi))
C      #  *L1(i,k))
C            else
C c           Fp=2.d0*fy*Df**3/(6.d0*(lsf1(i-1)+w1(i,k)  
C c     #     +sp1(i)))    
C c       Fo(i,k)=(F(i,k)*(1.d0+mu*sin(phi/2.d0))+Fp)      
C c     #  /(1.d0-mu*sin(phi/2.d0))  
C         Mp=fy*Df**3/6.d0
C         Fo(i,k)=(F(i,k)*L1(i,k)-6*Mp*mu)/((-3*sin(phi)*mu+cos(phi))
C      #  *L1(i,k))     
C            end if
C          else 
C          fo(i,k)=0.d0
C          endif
C        else       
C c       Fo(i,k)=F(i,k)*(1.d0+mu*sin(phi/2.d0))/(1.d0-mu*sin(phi/2.d0)) 
C         Mp=fy*Df**3/6.d0
C         Fo(i,k)=(F(i,k)*L1(i,k)-6*Mp*mu)/((-3*sin(phi)*mu+cos(phi))
C      #  *L1(i,k))         
C        end if
c       Fr=sin(phi/2.d0)*(F(i,k)+Fo(i,k))
 
            if(rotule) then
              if(phid.gt.5.d0) then
            if(w1(i,k).le.scrit) then
            D1=w1(i,k)/scrit
            else
            D1=1.d0
            end if
 
            if(l1(i,k).eq.0.d0) then 
            F(i,k)=0.d0
            else
            F(i,k)=F(i,k)+D1*3.d0*fy*Df**3.d0/(6.d0*sin(phi)*l1(i,k))
            end if
              end if
            end if
 
 
         Fo(i,k)=F(i,k)/(1.d0-mu*sin(phi)) 
         Fr=Fo(i,k)*sin(phi)
ccccccccccccccccccccccccccccccccccccccccccccccc
c ecaill ccccccccccccccccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccc             
c dans l ecaill on rentre la force radiale 
c on ne fait pas de calcul d ecaill si ce n est pas necessaire
 
c mise a zero de ff(i-1)       
      if(i.eq.1) then
      Fpasav=0.d0
      else
      Fpasav=Ff(i-1)
      end if
 
      if((Fo(i,k)-Fpasav.gt.0.d0).and.(.not.rupt)) then
 
      if((L1(i,k).eq.0.d0).and.(i.gt.1)) then     
      ls1(i,k)=lsf1(i-1)
      else
        if(L1(i,k).eq.0.d0) then
        ls1(i,k)=0.d0
        else
        call ecaill(i,k,phi,lsk1,alpha,L1(i,k),ls1(i,k),Df,nlign1,
     #  Rt,Fr,j) 
        end if     
      end if
 
      if((L2(i,k).eq.0.d0).and.(i.gt.1)) then 
      ls2(i,k)=lsf2(i-1)
      else
        if((L2(i,k).eq.0.d0).or.(ecail1)) then
        ls2(i,k)=0.d0
        else
        call ecaill(i,k,phi,lsk2,alpha,L2(i,k),ls2(i,k),Df,nlign1,
     #  Rt,Fr,j) 
        end if     
      end if
 
      else   
      if(i.eq.1) then
        ls1(i,k)=0.d0
        ls2(i,k)=0.d0
        else
        ls1(i,k)=lsf1(i-1)
        ls2(i,k)=lsf2(i-1)
        end if
      end if
 
c        la matrice coté 1 ne s'écaille plus si la force décroit
                 if(i.gt.1) then                 
                   if((Fo(i,k)-Fpasav.lt.0.d0).or.(Fo(i,k).eq.0.d0)
     #             .or.(ls1(i,k).lt.lsf1(i-1)).or.rupt)
     #             then
                   ls1(i,k)=lsf1(i-1)
                   end if
                 else
                 ls1(i,k)=0.d0
                 end if
 
c        la matrice coté 2 ne s'écaille plus si la force décroit
                 if(i.gt.1) then                 
                   if((Fo(i,k)-Fpasav.lt.0.d0).or.(Fo(i,k).eq.0.d0)
     #             .or.(ls2(i,k).lt.lsf2(i-1)).or.rupt)
     #             then
                   ls2(i,k)=lsf2(i-1)
                   end if
                 else
                 ls2(i,k)=0.d0
                 end if                 
 
 
c            coté le moins ancré
             if(i.eq.1) then
             dls1(i,k)=ls1(i,k)
             else
                 if(ls1(i,k)-lsf1(i-1).le.0.d0) then
                 dls1(i,k)=0.d0
                 else
                 dls1(i,k)=ls1(i,k)-lsf1(i-1)
                 end if
             end if
 
c            coté le plus ancré 
             if(i.eq.1) then
             dls2(i,k)=ls2(i,k)
             else
                 if(ls2(i,k)-lsf2(i-1).le.0.d0) then
                 dls2(i,k)=0.d0
                 else
                 dls2(i,k)=ls2(i,k)-lsf2(i-1)
                 end if
             end if
 
C                  if(i.gt.90) then
C                  if((j.eq.7).and.(t.eq.5).and.(l.eq.50)) then
C             print*,ls1(i,k),ls2(i,k),i,k,fo(i,k),Fpasav,fo(i,k)-Fpasav,
C      #dls1(i,k),L1(i,k),Fr
C                  read*
C                  end if 
C                  end if                 
 
 
c         if(isnan(ld1(i,k))) then
c         print*,"ld1",i,ldf1(i-1),ldind1,fel(i),L1(i,k)
C          else if(isnan(ld2(i,k))) then
C           print*,"ld2"        
C          else if(isnan(fcrit1(i,k))) then
C           print*,"fcrit1"          
C          else if(isnan(fcrit2(i,k))) then
C           print*,"fcrit2"          
C          else if(isnan(L1(i,k))) then
C           print*,"L1"  
C          else if(isnan(L2(i,k))) then
C           print*,"L2"  
C          else if(isnan(dls1(i,k))) then
C           print*,"dls1"  
C          else if(isnan(dls2(i,k))) then
C           print*,"dls2" 
C          else if(isnan(ls1(i,k))) then
C           print*,"ls1" 
C          else if(isnan(ls2(i,k))) then
C           print*,"ls2" 
C          else if(isnan(F(i,k))) then
C           print*,"F"
C          else if(isnan(w1(i,k))) then
C           print*,"w1"
C          else if(isnan(w2(i,k))) then
C           print*,"w2"          
C          else if(isnan(Fd1(i,k))) then
C           print*,"Fd1"  
C          else if(isnan(sd1(i,k))) then
C           print*,"sd1"  
C          else if(isnan(sel1(i,k))) then
C           print*,"sel1" 
C          else if(isnan(Fo(i,k))) then
C           print*,"Fo" 
C          else if(isnan(sd2(i,k))) then
C           print*,"sd2"          
C          else if(isnan(sel2(i,k))) then
C           print*,"sel2"
c         end if          
 
C           print*,fel(i),ld1(i,k),fcrit1(i,k),fcrit2(i,k),L1(i,k),L2(i,k)
C      #,sp1(i),dls1(i,k),dls2(i,k),ls1(i,k),ls2(i,k),sd1(i,k)
 
C           print*,fcrit1(i,k),L1(i,k),sp1(i),dls1(i,k),sd1(i,k),Ld1(i,k),
C      # F(i,k),ldf1(i-1),i,k,rupt
 
c test iter k
c posiibilité d enlever le if crit1=100?
             if(k.gt.1) then
             crit1=abs((dls1(i,k)-dls1(i,k-1)))/la1
             else
             crit1=100.d0
             end if
c             print*,dls1(i,k),i,k,j,l,crit1
c             read*
            if(((dls1(i,k).eq.0.d0).or.(crit1.lt.1.0d-6))
     #       .and.(k.gt.1)) then   
            vark=k            
c            exit
            go to 10
            else
            vark=niter
            end if
c fin boucle sur k        
            end do 
 
   10   continue
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c           sauvegarde des variables après itération sur k 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c           force finale orientée
            if(rupt) then
             Ff(i)=0.d0
            ffo(i)=0.d0
            else            
            Ff(i)=Fo(i,vark)
            ffo(i)=fo(i,vark)
            end if
 
c ls est recalculé avec la force actuelle donc il faut prendre celui d'avant 
c mais toutes les autres valeurs sont bonnes car elles sont actualisées avec 
c la valeur de dls(i,k-1)     
             if(((l1(i,vark).eq.0.d0).or.(Ffo(i).eq.0.d0)).and.(i.gt.1)) 
     #       then
             lsf1(i)=lsf1(i-1)
             lsf2(i)=lsf2(i-1)
             else
             lsf1(i)=ls1(i,vark-1) 
             lsf2(i)=ls2(i,vark-1)             
             end if
 
             if(lsf1(i).gt.la1) then
             lf1(i)=0.d0
             rupt=.true.
             else 
             Lf1(i)=La1-lsf1(i)             
             end if   
             if(lsf2(i).gt.la2) then
             lf2(i)=0.d0
             rupt=.true.
             else 
             Lf2(i)=La2-lsf2(i)             
             end if                
             Lf2(i)=La2-lsf2(i)
             ldf1(i)=ld1(i,vark)
             ldf2(i)=ld2(i,vark)                                                                        
             sdf1(i)=sd1(i,vark)
             sdf2(i)=sd2(i,vark)
C c            la matrice ne s'écaille plus si la force décroit (lsf est constant)
 
c****************************************************
c*********** Rupture du cone de matrice**************
c****************************************************
 
c    force admissible par le cone
      alphaco=30.d0*pi/180.d0
c on ne compte pas l ecaill dans le cone
       Fco(i)=PI*(Lf1(i)+lsf1(i)-sp1(i))*tan(alphaco)*(Df+(Lf1(i)
     #+lsf1(i)-sp1(i)))*Rt
         if((Ffo(i).ge.Fco(i)).or.(rupt)) then
c         print*,'cone'
         Ffo(i)=0.d0
         rupt=.true.
         end if
 
c****************************************************
c*********** Rupture de la fibre*********************
c****************************************************
 
c    force admissible par la fibre
         Fmax=(PI*Df**2/4.d0)*fu
         if((Ffo(i).ge.Fmax).or.(rupt)) then
         Ffo(i)=0.d0
         rupt=.true.
         end if          
 
c on ajoute la def des fibres sur la longueurs écaillée pour palier le fait que le glissement
c de décollement diminue
             if(i.gt.1) then
               if(rupt) then
               wf1(i,j+1,l)=wf1(i-1,j+1,l)+incremsp
               else               
             wf1(i,j+1,l)=w1(i,vark)+w2(i,vark)+(lsf1(i)+lsf2(i))*
     #        (1.d0-cos(phi))+sp1(i)+s0
     #       +ffo(i)*(lsf1(i)+lsf2(i))/(Af*Ef)
               end if
              if(wf1(i,j+1,l).lt.1.0d-10) then
              wf1(i,j+1,l)=0.d0
              end if
             else 
             wf1(i,j+1,l)=0.d0
             end if
 
c sauvegarde de la rigidité initiale pour une fibre de longueur ancrée connue
              if((i.eq.2).and.(wf1(i,j+1,l).gt.0.d0)) then
              sk0=wf1(i,j+1,l)
              fk0=fel(i)
              k0=fk0/sk0
 
                if(sk0.eq.0.d0) then
                print*,"rigidite initiale infinie",sk0,fk0  
                read*        
                end if              
              end if
 
c           tau(i)= Ffo(i)/(PI*Df*(Lf1(i)-sp1(i)))              
c          ici ffo et tau contiennent les valeurs pour un angle une longueur ancrée 
c          et toutes les i ouvertures de fissures
c          force unitaire pour chacun des parametres            
           Ffo1(i,j+1,l)=ffo(i)
 
 
 
c         affichage des courbes pour 5 longueurs ancrée            
C            if((j.eq.7).and.(t.eq.5)) then
C           if(l.eq.15) then
C           open(unit=1,file="ffol0.res")
C           write(unit=1,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
C           else if (l.eq.20) then
C           open(unit=1,file="ffol1.res")
C           write(unit=1,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
C           else if (l.eq.30) then
C           open(unit=1,file="ffol2.res")
C           write(unit=1,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
C           else if (l.eq.40) then
C           open(unit=1,file="ffol3.res")
C           write(unit=1,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))    
C           else if (l.eq.50) then
C           open(unit=1,file="ffol4.res")
C           write(unit=1,fmt='(11g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l)),
C      #w1(i,vark),w2(i,vark),lsf1(i),lsf2(i),sp1(i),s0,ffo(i),rupt
C           end if
C            endif
 
c plot 'ffol0.res','ffol1.res','ffol2.res','ffol3.res','ffol4.res'
 
c enregistrement des forces ouvertures de fissures pour plusieurs angles
C           if(l.eq.nlongfi) then
C           if(t.eq.nrt) then
C           if((j.eq.0)) then
C           open(unit=11,file="flee0.res")
C           write(unit=11,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
C           else if (j.eq.1) then
C           open(unit=11,file="flee1.res")
C           write(unit=11,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
C           else if (j.eq.3) then
C           open(unit=11,file="flee2.res")
C           write(unit=11,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
C           else if (j.eq.4) then
C           open(unit=11,file="flee3.res")
C           write(unit=11,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
C           else if (j.eq.6) then
C           open(unit=11,file="flee4.res")
C           write(unit=11,fmt='(2g13.5)') (wf1(i,j+1,l)),(ffo1(i,j+1,l))
C           end if
C           end if
C           end if
 
c    plot 'flee0.res','flee1.res','flee2.res','flee3.res','flee4.res'
c plot 'flee2.res','flee3.res'
 
C          if(isnan(ffo1(i,j+1,l))) then
C          print*,"ffo1"
C          else if(isnan(sp1(i))) then
C          print*,"sp1"
C         else if(isnan(lsf1(i))) then
C         print*,"lsf1"
C         else if(isnan(sdf1(i))) then
C         print*,"sdf1"        
C         else if(isnan(ldf2(i))) then
C         print*,"ldf2"            
C         else if(isnan(Ffo(i))) then
C         print*,"Ffo"   
C C         else if(isnan(lsf2(i))) then
C C         print*,"lsf2"           
C C         else if(isnan(Lf2(i))) then
C C         print*,"Lf2"            
C C         else if(isnan(sdf2(i))) then
C C         print*,"sdf2"        
C C         else if(isnan(Fco(i))) then
C C         print*,"Fco"          
C C         else if(isnan(Ff(i))) then
C C         print*,"Ff"
C         else if(isnan(Lf1(i))) then
C         print*,"Lf1"
C C         else if(isnan(ldf1(i))) then
C C         print*,"ldf1"
C C         else if(isnan(ffo2(i,j+1,l))) then
C C         print*,"ffo2"
C C         else if(isnan(Ffo1(i,j+1,l))) then
C C         print*,"Ffo1"
C C         else if(isnan(wf1(i,j+1,l))) then
C C         print*,"wf1" 
C          end if
 
 
 
 
 
 
C c fin boucle sur i          
           end do     
      do n=1,nwlis 
c      construit une liste d'ouverture de fissure     
       if(n.eq.1) then
       wlist(1)=0.d0
       else 
       wlist(n)=wlist(n-1)+incremwlis
       end if   
c calcule le minimum du carré de la diff entre les listes d ouverture de fissure
c on localise la valeur la plus proche de wlist(n)      
 
         do i=1,nlign1
         Dwf(i)=wlist(n)-wf1(i,j+1,l)
         if(dwf(i).lt.0.d0) then
         mi2=i
         mi1=i-1
c         exit
         go to 20
         end if
         end do
   20   continue       
c      les valeurs de wf1 qui encadrent wlis sont
       wint0=wf1(mi1,j+1,l)
       wint1=wf1(mi2,j+1,l)
c      les valeurs de ffo1 qui correspondent a wint0 et wint1 sont 
       fint0=ffo1(mi1,j+1,l) 
       fint1=ffo1(mi2,j+1,l)   
 
       call intrpo(wint0,wint1,fint0,fint1,wlist(n),fres)    
 
 
 
       if((wint0.eq.wint1).or.(fint0.eq.fint1).or.
     #  (fres.lt.0.d0)) then
       ffo2(n,j+1,l)=0.d0
       else
       ffo2(n,j+1,l)=fres
       end if
 
C           if(t.eq.4) then
C           if(j.eq.0) then
C c          print*,n,fres,l,wint0,wint1,fint0,fint1
C            print*,n,ffo2(n,j+1,l)
C           read*
C           end if
C           end if   
 
c affichage de la force orientée pour 10 angles
C          if((l.eq.2).and.(t.eq.1)) then
C           if(j.eq.0) then
C           open(unit=22,file="fphi0.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (j.eq.1) then
C           open(unit=22,file="fphi1.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (j.eq.2) then
C           open(unit=22,file="fphi2.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (j.eq.3) then
C           open(unit=22,file="fphi3.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (j.eq.4) then
C           open(unit=22,file="fphi4.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (j.eq.5) then
C           open(unit=22,file="fphi5.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (j.eq.6) then
C           open(unit=22,file="fphi6.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (j.eq.7) then
C           open(unit=22,file="fphi7.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l)) 
C           else if (j.eq.8) then
C           open(unit=22,file="fphi8.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (j.eq.9) then
C           open(unit=22,file="fphi9.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (j.eq.10) then
C           open(unit=22,file="fphi10.res")
C           write(unit=22,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           end if
C           end if
 
c         plot 'fphi0.res','fphi1.res','fphi2.res','fphi3.res','fphi4.res','fphi5.res','fphi6.res','fphi7.res','fphi8.res','fphi9.res','fphi10.res'
c plot 'fphi6.res'
 
C            if((j.eq.0).and.(t.eq.3)) then
C           if(l.eq.5) then
C           open(unit=17,file="fmoyphi0.res")
C           write(unit=17,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (l.eq.10) then
C           open(unit=17,file="fmoyphi1.res")
C           write(unit=17,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (l.eq.15) then
C           open(unit=17,file="fmoyphi2.res")
C           write(unit=17,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (l.eq.30) then
C           open(unit=17,file="fmoyphi3.res")
C           write(unit=17,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l))
C           else if (l.eq.40) then
C           open(unit=17,file="fmoyphi4.res")
C           write(unit=17,fmt='(2g13.5)') (wlist(n)),(ffo2(n,j+1,l)) 
C           end if        
C            end if
c   plot 'fmoyphi0.res','fmoyphi1.res','fmoyphi2.res','fmoyphi3.res','fmoyphi4.res'
 
       end do 
c       fin boucle n
c va chercher la force correspondant a la valeur de wlist
c  fin boucle sur l
        end do 
 
C         if(j.eq.5) then
C         print*,sum(test)/(nlongfi-n0longfi+1) ,phid,Lf/4*cos(phi)**2
C         read*
C         end if
 
c       calcul de la force moyenne ancrage variable pour chaque angle
          do n=1,nwlis
             do l=n0longfi,nlongfi
             if (l.eq.n0longfi) then
             somFL=ffo2(n,j+1,n0longfi)
             else 
             somFL=ffo2(n,j+1,l)+somFL
             endif
             end do 
c            fin de boucle sur l
 
          if(n.eq.1) then            
          FmoyL(1,j+1)=0.d0 
          else
          FmoyL(n,j+1)=somFL/(nlongfi-n0longfi+1) 
          end if
          Fmoylis(n)=FmoyL(n,j+1)
 
C           print*,Fmoylis(n),ffo2(n,j+1,l),wf1(n,j+1,l),n,l,t,nwlis,
C      #wlist(n)
 
c affichage de la force moyenne pour 10 angles
C           if(t.eq.4) then
C           if(j.eq.0) then
C           open(unit=19,file="fmoyphii0.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))
C           else if (j.eq.1) then
C           open(unit=19,file="fmoyphii1.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))
C           else if (j.eq.2) then
C           open(unit=19,file="fmoyphii2.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
C           else if (j.eq.3) then
C           open(unit=19,file="fmoyphii3.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))
C           else if (j.eq.4) then
C           open(unit=19,file="fmoyphii4.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
C           else if (j.eq.5) then
C           open(unit=19,file="fmoyphii5.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
C           else if (j.eq.6) then
C           open(unit=19,file="fmoyphii6.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
C           else if (j.eq.7) then
C           open(unit=19,file="fmoyphii7.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
C           else if (j.eq.8) then
C           open(unit=19,file="fmoyphii8.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))
C           else if (j.eq.9) then
C           open(unit=19,file="fmoyphii9.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
C           else if (j.eq.10) then
C           open(unit=19,file="fmoyphii10.res")
C           write(unit=19,fmt='(2g13.5)') (wlist(n)),(FmoyL(n,j+1))  
C           end if
C           end if
 
c         plot 'fmoyphii0.res','fmoyphii1.res','fmoyphii2.res','fmoyphii3.res','fmoyphii4.res','fmoyphii5.res','fmoyphii6.res','fmoyphii7.res','fmoyphii8.res','fmoyphii9.res','fmoyphii10.res'
 
 
c  forces moyennes pour 6 angles et 5 contraintes laterales
C        if(t.eq.1) then
C          if(j.eq.0) then 
C          open(unit=19,file="fmoyt1j0.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          elseif (j.eq.1) then
C          open(unit=19,file="fmoyt1j1.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)          
C          elseif (j.eq.2) then
C          open(unit=19,file="fmoyt1j2.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)   
C          elseif (j.eq.3) then
C          open(unit=19,file="fmoyt1j3.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
C          elseif (j.eq.4) then
C          open(unit=19,file="fmoyt1j4.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
C          elseif (j.eq.5) then
C          open(unit=19,file="fmoyt1j5.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          end if
C        elseif(t.eq.2) then
C          if(j.eq.0) then 
C          open(unit=19,file="fmoyt2j0.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          elseif (j.eq.1) then
C          open(unit=19,file="fmoyt2j1.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)       
C          elseif (j.eq.2) then
C          open(unit=19,file="fmoyt2j2.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
C          elseif (j.eq.3) then
C          open(unit=19,file="fmoyt2j3.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)  
C          elseif (j.eq.4) then
C          open(unit=19,file="fmoyt2j4.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          elseif (j.eq.5) then
C          open(unit=19,file="fmoyt2j5.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
C          end if
C        elseif(t.eq.3) then
C          if(j.eq.0) then 
C          open(unit=19,file="fmoyt3j0.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          elseif (j.eq.1) then
C          open(unit=19,file="fmoyt3j1.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)         
C          elseif (j.eq.2) then
C          open(unit=19,file="fmoyt3j2.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)  
C          elseif (j.eq.3) then
C          open(unit=19,file="fmoyt3j3.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)  
C          elseif (j.eq.4) then
C          open(unit=19,file="fmoyt3j4.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
C          elseif (j.eq.5) then
C          open(unit=19,file="fmoyt3j5.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          end if
C        elseif(t.eq.4) then
C          if(j.eq.0) then 
C          open(unit=19,file="fmoyt4j0.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          elseif (j.eq.1) then
C          open(unit=19,file="fmoyt4j1.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)         
C          elseif (j.eq.2) then
C          open(unit=19,file="fmoyt4j2.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)  
C          elseif (j.eq.3) then
C          open(unit=19,file="fmoyt4j3.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
C          elseif (j.eq.4) then
C          open(unit=19,file="fmoyt4j4.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          elseif (j.eq.5) then
C          open(unit=19,file="fmoyt4j5.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          end if
C        elseif(t.eq.5) then
C          if(j.eq.0) then 
C          open(unit=19,file="fmoyt5j0.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
C          elseif (j.eq.1) then
C          open(unit=19,file="fmoyt5j1.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)          
C          elseif (j.eq.2) then
C          open(unit=19,file="fmoyt5j2.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1) 
C          elseif (j.eq.3) then
C          open(unit=19,file="fmoyt5j3.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          elseif (j.eq.4) then
C          open(unit=19,file="fmoyt5j4.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          elseif (j.eq.5) then
C          open(unit=19,file="fmoyt5j5.res")
C          write(unit=19,fmt='(11g13.5)') wlist(n),FmoyL(n,j+1)
C          end if
C          end if   
 
 
 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc              
c sauvergarde des paramètres pour la réduction du modèle
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 
c Rigidité initiale
       if(n.eq.2) then
       k0m=(FmoyL(2,j+1)-FmoyL(1,j+1))/(wlist(2)-wlist(1))
       end if
 
c force et ouverture fissure fin d extraction       
       if((FmoyL(n,j+1).eq.0.d0).and.(n.gt.1).and.(.not.finf)) then
c       Ffin=FmoyL(n-1,j+1)
       wfin=wlist(n-1)
       finf=.true.
       end if
 
C        open(unit=13,file="fmoyphi3d.res")
C        if(n.eq.1) then
C        write(unit=13,fmt='(11g13.5)')
C        end if
C        write(unit=13,fmt='(11g13.5)') wlist(n),phid,FmoyL(n,j+1)
 
 
c fin boucle sur n   
          end do
 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc              
c sauvergarde des paramètres pour la réduction du modèle
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc   
c force et ouverture de fissure au pic
cc      picloc=maxloc(Fmoylis)      
cc       npic=int(picloc(1))
       call amxloc(Fmoylis,Fmoylis(/1),npic)
       Fpic=FmoyL(npic,j+1) 
       wpic=wlist(npic)
 
       F03=0.3d0*Fpic
 
c recherche de l'ouverture de fissure correspondant a f03 et f06       
       do n=npic,nwlis
C         if((fmoylis(n)-F06.le.0.d0).and.
C      #  (fmoylis(n-1)-F06.ge.0.d0)) then
C         w06=wlist(n)
C         end if
 
        if((fmoylis(n)-F03.le.0.d0).and.
     #  (fmoylis(n-1)-F03.ge.0.d0)) then
        w03=wlist(n)
c        exit
        go to 30
        end if
       end do
   30  continue       
 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccc Déplacement total cccccccccccccccccccccc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c on doit reconstruire la liste de déplacement après qu'on connaisse le pic de multifissuration
c pour gérer la localisation qui se fait au pic de force en fait donc est ce qu on a besoin de faire ca ici?
c et puis est ce qu on devrait pas faire la multifissuration qu a la fin?
 
 
c affichage de la force en 3d en fonction de w et angle
c         splot 'fmoyphi3d.res' u 1:2:3
C            open(unit=13,file="fmoyphi3d.res")
C            if(n.eq.0) then
C            write(unit=13,fmt='(11g13.5)')
C            end if
C            write(unit=13,fmt='(11g13.5)') wlist(n),phid,FmoyL(n,j+1)
 
c          flis(n)=FmoyL(n,j+1)
c          print*,flis            
 
c    2eme fin de boucle sur n pour calcul deplacement            
c            end do            
 
C            open(unit=7,file="kmoyphiredu3d.res")
C            if(j.eq.0) then
C            write(unit=7,fmt='(11g13.5)')
C            end if
C            write(unit=7,fmt='(11g13.5)') Rt,phid,k0m 
c liste de rigidité initiale
                lisRt(p)=Rt 
                lisphid(p)=phid 
                lisk0m(p)=k0m               
 
C            open(unit=15,file="fpmoyphiredu3d.res")
C            if(j.eq.0) then
C            write(unit=15,fmt='(11g13.5)')
C            end if
C            write(unit=15,fmt='(11g13.5)') Rt,phid,Fpic 
c liste de force au pic
                lisFp(p)= Fpic 
 
C            open(unit=18,file="wpmoyphiredu3d.res")
C            if(j.eq.0) then
C            write(unit=18,fmt='(11g13.5)')
C            end if
C            write(unit=18,fmt='(11g13.5)') Rt,phid,wpic
c liste d'ouverture au pic
                if(j.lt.nangredu) then
                lisRt1(d)=Rt 
                lisphid1(d)=phid 
                liswp1(d)= wpic 
 
C            open(unit=181,file="wp1moyphiredu3d.res")
C            if(j.eq.0) then
C            write(unit=181,fmt='(11g13.5)')
C            end if
C          write(unit=181,fmt='(11g13.5)') lisrt1(d),lisphid1(d),liswp1(d)    
                else
c                listd2(1)=listd1(d-1)
                lisRt2(e)=Rt 
                lisphid2(e)=phid 
c                lisphid2(1)=lisphid1(d-1)
                liswp2(e)=wpic
c                liswp2(1)=liswp1(d-1)
 
C            open(unit=182,file="wp2moyphiredu3d.res")
C            if(j.eq.nangredu) then
C            write(unit=182,fmt='(11g13.5)')
C            end if
C          write(unit=182,fmt='(11g13.5)') lisrt2(e),lisphid2(e),liswp2(e)                
                end if          
 
               if(j.lt.nangredu) then
               lisk01(d)= k0m
               d=d+1 
               else
               lisk02(e)=k0m 
                e=e+1                         
               end if 
 
C              print*, k0m,e,j,t,phid
C              read*             
 
 
C            open(unit=21,file="wfmoyphiredu3d.res")
C            if(j.eq.0) then
C            write(unit=21,fmt='(11g13.5)')
C            end if
C            write(unit=21,fmt='(11g13.5)') Rt,phid,wfin
c liste d'ouverture finale
                liswf(p)= wfin
 
c liste d'ouverture a 30% de force 
 
C            open(unit=24,file="w03moyphiredu3d.res")
C            if(j.eq.0) then
C            write(unit=24,fmt='(11g13.5)')
C            end if
C            write(unit=24,fmt='(11g13.5)') Rt,phid,w03           
                lisw03(p)=w03    
 
 
c p numéro de ligne pour la réduction   
c            print*,"Force moyenne ",p,"effectuee"       
            p=p+1            
 
c fin boucle sur j
        end do
 
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c fin partie prpfib      
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc            
 
C C         print*,ld2
C C         print*,ld1         
          return
 
 
       end 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

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