Télécharger la plaquette Cast3M


fr en

Télécharger fibs3d.eso

Retour à la liste

Numérotation des lignes : 

fibs3d

  1. C FIBS3D    SOURCE    PV090527  24/06/12    21:15:04     11940          
  2.       subroutine fibs3d(Df,Lf,lech,mw,rtec,rhof,eo,vf1,vf2,Gf,
  3.      #cokf1,cokf2,cowp1,cowp2,cofp,cow03,cowf,mink01,maxk01,mink02,
  4.      #maxk02,minwp1,maxwp1,minwp2,maxwp2,minfp,maxfp,minw03,maxw03,
  5.      #minwf,maxwf,sigf06,vwpl6,wmax,wmoy,etw,ncf,sigfr,wmax6,le,
  6.      #nc0,wmoy0,etw0,vw30,rig0t,wmax0,vw3,rmax1,rmax2,rmax3,phimoy,
  7.      #wferm,sigdec,iloc,wpsigt,wmiloc,sigmaxt,vecw33,sigfx6,wmaxd,
  8.      #wmaxd0,epstf6,Ef,fu,young,xmt,dtiso,epict,rt00,cvar,wplmax)
  9.  
  10.       implicit real*8 (a-h,o-z)
  11.       implicit integer (i-n)   
  12.       integer err1,nvari,ncalc,rmax,r,rmaxi,rmax1,rmax2,rmax3
  13.       integer i,j,k,l,iloc(3)
  14.       parameter (PI=3.141592653589793D0)
  15.       real*8 incremwlis,precision1,precision2
  16.       parameter (incremwlis=2.5d-6,precision1=1.0d-8,precision2=1.0d-8)
  17.  
  18.  
  19.       real*8 rhof,Df,Lf,phicrit,Rtec,mw,Rtb,lech,le(3),Gf
  20.       real*8 young,xmt,epict,dt3
  21.  
  22.       parameter (nvari=7,phicrit=55.d0,npos=10)
  23. c attention npos doit être redefini dans disp3d
  24. c Rtec est celui de lech     
  25. C       parameter(mu=0.5d0,Td=6.4d0,Tmax=8.d0,rhof=0.02,
  26. C      # Ef=210000.d0,Df=2.0d-4,Hf=8000000.d0,sk=3.0d-2,
  27. C      #  F0=5.0d-6,alphad=30.d0,M0=8.d0,
  28. C      # Rtec=4.d0,fu=2800.d0,L0=5.0d-3,
  29. C      # Lf=0.013d0,lech=1.d0,mw=12.d0)
  30. c       parameter (ncalc=max(int(Lf/(incremwlis)),1000))
  31.       parameter (ncalc=1000)
  32.       parameter (rmaxi=10)  
  33.       real*8 vecw(3),wmax0(3)
  34.       real*8 Eo(3),Vwpl6(6),vwpl33(3,3),vw3(3)
  35.       real*8 veo(3,3),vecw33(3,3)
  36.       real*8 stri(npos,2*npos-2),xg(npos,2*npos-2),yg(npos,2*npos-2)
  37.       real*8 zg(npos,2*npos-2),vei(3)
  38.       real*8 phi1,phi1d(npos,2*npos-2),f1(npos,2*npos-2)
  39.       real*8 pscal(npos,2*npos-2),normv
  40.       real*8 si0(npos,2*npos-2),si1,Reiw(npos,2*npos-2)
  41.       real*8 wmaxa(3,3),wmaxac(3,3),wmaxact(3),wmax6(6),wmax33(3,3)
  42.       real*8 rig0(3,rmaxi),sig0(3),vwfi(3,ncalc)
  43.       real*8 sigmax(3,rmaxi),sigf(3,rmaxi,ncalc),sigf06(6)
  44.       real*8 sigf033(3,3),sigf03(3)
  45.       real*8 ratio(npos,2*npos-2),vwfi0(3),wfis(3)
  46.       real*8 wpsig(3,rmaxi),wmi0(1),cvar(3),wplmax(3)
  47.       real*8 wmiloc(3),phimoy(3),nf(npos,2*npos-2)
  48.       real*8 Romeg(npos,2*npos-2),Sr(3),Sigb(3),Sig(3),wfist(3)
  49.       real*8 nc(rmaxi),Rt(rmaxi),nct(rmaxi),Romeg1(npos,2*npos-2)
  50.       real*8 wmi(rmaxi),wt(3,ncalc),sigt(3,ncalc),wmul(3,rmaxi)
  51.  
  52. c variables pour la reduction du modele      
  53.       real*8 k0m, wp,fp,f03,w03,wf
  54.       real*8 cokf1(10),cokf2(10),cowp1(10),cowp2(10)
  55.       real*8 cofp(10),cowf(10),cow03(10)
  56.       real*8 maxfp,minfp,maxwp1,minwp1,maxwf,minwf,maxw03,minw03
  57.       real*8 maxwp2,minwp2,B,vf1(3),vf2(3)
  58.       real*8 maxk02,mink02,maxk01,mink01,phi2d(npos,2*npos-2)     
  59.       real*8 tri(npos,2*npos-2,3,3),Reiw1,Reiw2,Reiw3,vst1(3),vst2(3)
  60.       real*8 vst3(3),theta1,gamma1,theta2,gamma2,theta3,gamma3
  61.  
  62. c    variable ajoutée pour contrainte totale 
  63.       real*8 wpsigt(3),rig0t(3),sigmaxt(3),sigfr(3),rt00    
  64.       logical loc(3),dtiso
  65.       real*8 varft(nvari,ncalc),rig0td(3),gfap(3)
  66.       real*8 wmr(rmaxi),wer(rmaxi),wmoy(3),etw(3),wmax(3)
  67.       real*8 sigfp33(3,3),sigfx33(3,3),sigfx6(6),ncf(3)
  68.       real*8 dw3(3),vw30(3),wmoy0(3),etw0(3),nc0(3)
  69.       real*8 wferm(3),sigdec(3),sigfe(3),epse13(3),vecwj(3)
  70.       real*8 wmaxd(3),wmaxd0(3),w1(3),rtmin,epstf6(6),epse133(3,3)
  71.       real*8 vecs33(3,3)
  72. * pv
  73.       s1 = 0.d0
  74.       do i=1,3     
  75.        sig(i)=-1.d0
  76.       enddo
  77. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  78. c CONSTRUCTION DES 3 COURBES CONTRAINTES - OUVERTURE TOTALES DANS LA BASE 
  79. c CCCCCCCCCCCCCCCC  DE FISSURE ACTUELLE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  80. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  81.  
  82.  
  83. c recuperation des ouvertures max
  84.       call x6x33(wmax6,wmax33)      
  85.       call x6x33(vwpl6,vwpl33)
  86. c     vw3 : vecteur des valeurs principales de fissures
  87. c     vecw33 matrice des vecteurs d'orientation des fissures actuelles
  88.       call b3_v33(vwpl33,vw3,vecw33)
  89.  
  90. C        print*,vecw33(1,1),vecw33(2,1),vecw33(3,1),"vw1"
  91. C        print*,vecw33(1,2),vecw33(2,2),vecw33(3,2),"vw2"
  92. C        print*,vecw33(1,3),vecw33(2,3),vecw33(3,3),"vw3"
  93.  
  94.  
  95. c on passe le tenseur d'ouverture max dans la base de fissure actuelle
  96. c pour comparer les valeurs max et actuelles
  97. c      vecw33t=transpose(vecw33)
  98.       wmaxa=matmul(transpose(vecw33),wmax33)
  99.       wmaxac=matmul(wmaxa,vecw33)
  100.  
  101. c on récupère les valeurs diag
  102. c c'est cette valeur qu'on va décomposer en multifissure
  103.       do i=1,3
  104.         do j=1,3
  105.         if(i.eq.j) then 
  106.         wmaxact(i)=wmaxac(i,j)
  107.         end if
  108.         end do
  109.       end do    
  110.  
  111. c recuperation du tenseur des contraintes fibres du pas d avant  
  112. c il sert a calculer une eventuelle decharge   
  113.       call x6x33(sigf06,sigf033)
  114.       call b3_v33(sigf033,sigf03,vecs33) 
  115. c      print*,sigf03(1),sigf03(2),sigf03(3),"contraintes fibres 0"    
  116.  
  117. c tenseur d orientation des fibres dans la base d'orientations principales     
  118.       if((Eo(1).eq.0.d0).or.(Eo(2).eq.0.d0).or.(Eo(3).eq.0.d0)) then
  119.         print*,"arret du prigramme dans fibs3d.eso"
  120.         print*,"Donnez des valeurs du tenseur d'orientation non nulles"
  121.         stop
  122.         err1=1
  123.       end if      
  124.  
  125. c vecteurs principaux du tenseur d'orientation des fibres dans la base du maillage
  126. c on donne 2 vecteurs et on calcule le 3ème 
  127.       veo(1,1)=vf1(1)
  128.       veo(2,1)=vf1(2)
  129.       veo(3,1)=vf1(3)
  130.  
  131.       veo(1,2)=vf2(1)
  132.       veo(2,2)=vf2(2)
  133.       veo(3,2)=vf2(3)
  134.  
  135.       veo(1,3)=veo(2,1)*veo(3,2)-veo(3,1)*veo(2,2)
  136.       veo(2,3)=veo(3,1)*veo(1,2)-veo(1,1)*veo(3,2)
  137.       veo(3,3)=veo(1,1)*veo(2,2)-veo(2,1)*veo(1,2)
  138.  
  139.       if(abs(
  140.      #  veo(1,1)*veo(1,2)+veo(2,1)*veo(2,2)+veo(3,1)*veo(3,2)).gt.
  141.      #  precision1) then
  142.         print*,"arret du prigramme dans fibs3d.eso"
  143.         print*,"Donnez vecteurs orthogonaux d'orientation des fibres "
  144. c       read*
  145.         stop
  146.         err1=1
  147.       end if
  148.  
  149. c    normalisation des vecteurs
  150.       do i=1,3
  151.           normv=sqrt(veo(1,i)**2+veo(2,i)**2+veo(3,i)**2)
  152.           veo(1,i)=veo(1,i)/normv
  153.           veo(2,i)=veo(2,i)/normv
  154.           veo(3,i)=veo(3,i)/normv
  155.       end do
  156. c  discrétisation de l'espace en angles solides  
  157.       call disp3d(stri,xg,yg,zg,tri,npos)             
  158. c traitement des fissures
  159.       do i=1,3  
  160. c       initialisation
  161.         loc(i)=.false.  
  162. c       wmi=0.d0
  163.         if(i.eq.1) then
  164.             rmax=rmax1
  165.         else if(i.eq.2) then
  166.             rmax=rmax2
  167.         else if(i.eq.3) then
  168.             rmax=rmax3
  169.         end if
  170.  
  171.         if(vw3(i).lt.precision2) then
  172.             sigfr(i)=0.d0
  173.             wmax(i)=0.d0
  174.             rig0t(i)=0.d0
  175.             wmoy(i)=0.d0
  176.             etw(i)=0.d0    
  177.             wmi=0.d0
  178.             wmul=0.d0
  179.             wmr=0.d0 
  180.             wer=0.d0 
  181.             nc=0.d0
  182.             nct=0.d0
  183.  
  184.             if(wmaxact(i).lt.precision2) then
  185.                 wmoy(i)=0.d0
  186.                 etw(i)=0.d0
  187.                 ncf(i)=0.d0
  188.                 sigfr(i)=0.d0
  189.                 if(nc0(i).gt.0.d0) then
  190.                     ncf(i)=nc0(i)
  191.                 else
  192.                     ncf(i)=0.d0
  193.                 end if     
  194.             else  
  195.                 ncf(i)=nc0(i)
  196.             end if
  197.  
  198. c           calcul de phimoy
  199.  
  200.             phimoy(i)=0.d0
  201.  
  202.  
  203.             go to 70
  204.         end if
  205.  
  206.         if((wmaxact(i).gt.0.d0).and.(nc0(i).eq.0.d0)) then
  207.             nc0(i)=1.d0
  208.         end if
  209.  
  210. c       vecteurs orientations de fissure
  211. c       vecw(j,i) : orientation de la fissure i (normé? --> oui)
  212.         vecw(1)=vecw33(1,i)
  213.         vecw(2)=vecw33(2,i)
  214.         vecw(3)=vecw33(3,i)
  215.  
  216.  
  217. c       on construit les rayons de l ellipse et les angles une seule fois pour toute
  218. c       c'est les memes pour tous         
  219.         do n=1,npos
  220.           do k=1,2*npos-2           
  221. c           vecteur sommet des triangles vst(sommet,coordonnée)
  222. c           sommet 1
  223.             vst1(1)=tri(n,k,1,1) 
  224.             vst1(2)=tri(n,k,1,2)
  225.             vst1(3)=tri(n,k,1,3) 
  226. c             sommet 2
  227.             vst2(1)=tri(n,k,2,1) 
  228.             vst2(2)=tri(n,k,2,2)
  229.             vst2(3)=tri(n,k,2,3) 
  230. c           sommet 3
  231.             vst3(1)=tri(n,k,3,1) 
  232.             vst3(2)=tri(n,k,3,2)
  233.             vst3(3)=tri(n,k,3,3)          
  234. c           vecteur dans la direction ei (il est déjà normalisé dans disp3d)
  235.             vei(1)=xg(n,k)
  236.             vei(2)=yg(n,k)
  237.             vei(3)=zg(n,k)
  238.  
  239. c           passage de vei dans la base du maillage puis de fissuration 
  240. c           NON vecw est dans la base du maillage, il suffit de passer 
  241. c           vei dans la base du maillage 
  242. c           veo est la matrice de passage P du maillage a l ellipse
  243. c           ici on veut passer de l'ellipse au maillage donc P A
  244.             vei=matmul(veo,vei)                 
  245. c           vei=matmul(vecw33,vei)          
  246. c           produit scalaire entre l'orientation de la fissure et de ei
  247.             pscal(n,k)=abs(sum(vei*vecw))
  248. c           calcul de l'angle entre la fissure et la direction ei
  249.             phi1=acos(pscal(n,k))  
  250.             phi1d(n,k)=180.d0*phi1/pi  
  251. c           calculs sont effectues dans la base de l ellipse          
  252. c           calcul des angles d'orientation de chaque sommet du triangle n,k
  253.             theta1=acos(vst1(3))
  254.               if(vst1(1).eq.0.d0) then
  255.                 gamma1=pi/2.d0
  256.               else
  257.                 gamma1=atan(vst1(2)/vst1(1))  
  258.               end if
  259.               theta2=acos(vst2(3))
  260.               if(vst2(1).eq.0.d0) then
  261.                 gamma2=pi/2.d0
  262.               else
  263.                 gamma2=atan(vst2(2)/vst2(1))
  264.               end if        
  265.               theta3=acos(vst3(3))
  266.               if(vst3(1).eq.0.d0) then
  267.                 gamma3=pi/2.d0
  268.               else
  269.                 gamma3=atan(vst3(2)/vst3(1))
  270.               end if  
  271.  
  272. c rayon de l'ellipse d'orientation des fibres dans la direction des sommets             
  273.       Reiw1 = (Eo(1) ** 2 * Eo(3) ** 2 * sin(theta1)**2*sin(gamma1)** 2 
  274.      #- Eo(2) ** 2 * Eo(3) ** 2 * sin(theta1) ** 2*sin(gamma1) ** 2-Eo(1
  275.      #) ** 2 * Eo(2) ** 2 * sin(theta1) ** 2 + Eo(2) ** 2 * Eo(3) ** 2 * 
  276.      #sin(theta1) ** 2+ Eo(1) ** 2 * Eo(2) ** 2) ** (-0.1D1 / 0.2D1) * E
  277.      #o(3) * Eo(2) * Eo(1)
  278.  
  279.       Reiw2 = (Eo(1) ** 2 * Eo(3) ** 2 * sin(theta2)**2*sin(gamma2)** 2 
  280.      #- Eo(2) ** 2 * Eo(3) ** 2 * sin(theta2) ** 2*sin(gamma2) ** 2-Eo(1
  281.      #) ** 2 * Eo(2) ** 2 * sin(theta2) ** 2 + Eo(2) ** 2 * Eo(3) ** 2 * 
  282.      #sin(theta2) ** 2+ Eo(1) ** 2 * Eo(2) ** 2) ** (-0.1D1 / 0.2D1) * E
  283.      #o(3) * Eo(2) * Eo(1)     
  284.  
  285.       Reiw3 = (Eo(1) ** 2 * Eo(3) ** 2 * sin(theta3)**2*sin(gamma3)** 2 
  286.      #- Eo(2) ** 2 * Eo(3) ** 2 * sin(theta3) ** 2*sin(gamma3) ** 2-Eo(1
  287.      #) ** 2 * Eo(2) ** 2 * sin(theta3) ** 2 + Eo(2) ** 2 * Eo(3) ** 2 * 
  288.      #sin(theta3) ** 2+ Eo(1) ** 2 * Eo(2) ** 2) ** (-0.1D1 / 0.2D1) * E
  289.      #o(3) * Eo(2) * Eo(1)     
  290.  
  291.             Reiw(n,k)=(Reiw1+Reiw2+Reiw3)/3.d0         
  292.  
  293. c           stockage de reiw*omegai
  294.             Romeg(n,k)=Reiw(n,k)*stri(n,k)
  295.             Romeg1(n,k)=Reiw(n,k)*stri(n,k)*pscal(n,k)  
  296. c           terme proportionnel à la densité volumique de fibres dans une direction ei
  297. c           car on ne sait pas calculer la somme rei*omegai 
  298. c           ratio(n,k)=Reiw(n,k)*stri(n,k)*rhof     
  299. c           proportionnel au nombre de fibres du paquet n,k traversant la fissure i  
  300. c           nf(n,k)=4*ratio(n,k)*pscal(n,k)/(pi*Df**2)        
  301. c           proportionnel a l ancrage moyen des fibres du paquet ei pondéré par le nombre de fibres traversant la fissure        
  302. c           Lm(n,k)=(Lf/2-Lf/(4*cos(phi1)))*nf(n,k)
  303. c           l angle moyen reflete les fibres qui traversent effectivement le plan
  304.             phi2d(n,k)=Reiw(n,k)*stri(n,k)*phi1d(n,k)*pscal(n,k)  
  305. c           fin de la 1ere boucle sur n,k qui sert a connaitre reiw(n,k) et phi1d(n,k)          
  306.           end do 
  307. c         fin de boucle sur k      
  308.         end do 
  309. c        fin de boucle sur n          
  310. c     angle moyen
  311.       phimoy(i)=sum(phi2d)/sum(romeg1) 
  312.  
  313.       if(vw3(i).gt.wplmax(i)) then
  314.       wplmax(i)=vw3(i)
  315.       end if
  316.  
  317.  
  318. c      if((vw3(i).ge.wmaxact(i)).and.(vw3(i).ge.vw30(i))) then  
  319. c      if(vw3(i).ge.wmaxact(i)) then 
  320.  
  321.       if(vw3(i).lt.wferm(i)) then
  322.       go to 100
  323.       end if           
  324.  
  325.       if((vw3(i).ge.wplmax(i))) then        
  326. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  327. c on ouvre la fissure i    
  328. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  329. c       modif
  330. c      wpsigt0(i)=wpsigt(i)
  331.  
  332.       rtmin=Rtec*(lech/(0.5d0*sqrt(2*pi)))**(1/mw)  
  333.       do r=1,rmax
  334.         if(r.eq.1) then
  335.             Rt(r)=rt00
  336.         else
  337.             Rt(r)=Rtec*(lech*2.d0**(r-1)/le(i))**(1/mw)
  338.         end if
  339.  
  340.         if(rt(r).lt.rtmin) then
  341.         rt(r)=rtmin
  342.         end if
  343.       end do   
  344. c   Rtec valeur pour lech quand r=1 on fait les fissures a le/2   
  345.  
  346. c   initialisation      
  347.       wmi=0.d0
  348. cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  349. c ici on test si on etait dans une phase de localisation
  350. c si c est le cas on ne reconstruit pas les courbes on va 
  351. c directement dans la phase de localisation
  352. cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  353.  
  354.       if(iloc(i).eq.1) then
  355.         go to 80
  356.       end if
  357.  
  358.  
  359.  
  360.       do r=1,rmax
  361.         Rtb=Rt(r)
  362. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
  363. c construction des courbes contraintes ouvertures de fissures pour chaque rt
  364. c pour 1 fissure fictive incrementale dans la base de fissuration actuelle   
  365. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc   
  366.         do l=1,ncalc
  367.  
  368. c           récupération de la contrainte du pas fictif d'avant
  369.             if(l.eq.1) then
  370.                 sig0(1)=0.d0
  371.                 sig0(2)=0.d0
  372.                 sig0(3)=0.d0
  373.                 vwfi0(1)=0.d0
  374.                 vwfi0(2)=0.d0
  375.                 vwfi0(3)=0.d0         
  376.             else
  377.                 sig0(1)=varft(1,l-1)
  378.                 sig0(2)=varft(2,l-1) 
  379.                 sig0(3)=varft(3,l-1) 
  380.                 vwfi0(1)=varft(4,l-1)
  381.                 vwfi0(2)=varft(5,l-1)
  382.                 vwfi0(3)=varft(6,l-1)        
  383.             end if 
  384.  
  385. c           fissure fictive
  386.             vwfi(i,l)=l*incremwlis
  387.  
  388. c           et ici on recommence donc une boucle ou on va calculer les contraintes
  389. c           jusqu au pic pour tous les Rt
  390.             do n=1,npos
  391.                 do k=1,2*npos-2            
  392. c               calcul de la force moyenne pour une fibre traversant la fissure i inclinée de phi1d
  393. c               calcul des points caractéristiques de la courbe force ouverture de fissure pour Td et phid
  394. c               rigidité initiale
  395.                 if(phi1d(n,k).le.phicrit) then 
  396.                     call reconf(cokf1,k0m,Rtb,phi1d(n,k),mink01,maxk01)
  397.                 else
  398.                     call reconf(cokf2,k0m,Rtb,phi1d(n,k),mink02,maxk02)
  399.                 end if         
  400. c               ouverture au pic
  401.                 if(phi1d(n,k).le.phicrit) then 
  402.                     call reconf(cowp1,wp,Rtb,phi1d(n,k),minwp1,maxwp1)
  403.                 else
  404.                     call reconf(cowp2,wp,Rtb,phi1d(n,k),minwp2,maxwp2)
  405.                 end if          
  406. c               force au pic
  407.                 call reconf(cofp,fp,Rtb,phi1d(n,k),minfp,maxfp)  
  408. c               force intermédiaire post pic     
  409.                 f03=0.3d0*fp          
  410. c               ouverture pour f=0.3fp  
  411.                 call reconf(cow03,w03,Rtb,phi1d(n,k),minw03,maxw03)
  412.                 if(w03.le.wp) then        
  413.                     w03=3.0d0*wp  
  414.                 end if  
  415. c               ouverture finale (f=0)    
  416.                 call reconf(cowf,wf,Rtb,phi1d(n,k),minwf,maxwf)
  417.                 if(wf.le.w03) then        
  418.                     wf=3.0d0*w03  
  419.                 end if  
  420. c               calcul de la force dans la fissure i pour les fibres orientées selon ei
  421.                 if(vwfi(i,l).le.wp) then
  422.                     B=(k0m*wp-fp)/(fp*wp)
  423.                     f1(n,k)=k0m*vwfi(i,l)/(B*vwfi(i,l)+1.d0)
  424.                 else 
  425.                  if(vwfi(i,l).lt.wf) then
  426.                     f1(n,k)=fp*exp(log(f03/fp)*(vwfi(i,l)-wp)/(w03-wp))
  427. c                   print*,2,f1(n,k),wp,wf
  428.                  else
  429.                     f1(n,k)=precision1 
  430.                  end if
  431.                 end if
  432. c               calcul du terme à l'intérieur de la somme de la contrainte
  433.                 si0(n,k)=Reiw(n,k)*stri(n,k)*pscal(n,k)*f1(n,k)       
  434.             end do 
  435. c           fin de boucle sur k      
  436.           end do 
  437. c        fin de boucle sur n 
  438. c        valeurs fictives 
  439.          si1=sum(si0)   
  440.          sigf(i,r,l)=4.d0*rhof*si1/(pi*Df**2*(sum(Romeg))) 
  441.  
  442.  
  443. c        stockage de la rigidité initiale
  444.          if(l.eq.1) then        
  445.             rig0(i,r)=sigf(i,r,l)/incremwlis
  446.          end if      
  447.             if(rig0(i,r).ne.rig0(i,r).or.(rig0(i,r).eq.0.d0)) then
  448.  
  449.             print*, "Probleme de rigidite initiale",rig0(i,r)
  450.             err1=1
  451.          end if  
  452.  
  453.          if((sigf(i,r,l).lt.sig0(i)).or.(vwfi(i,l).gt.(vw3(i)
  454.      # +2.d0*incremwlis))) then 
  455. c        on calcule la contrainte jusqu'au pic uniquement ou l ouverture actuelle
  456.          wpsig(i,r)=vwfi0(i)
  457.          sigmax(i,r)=sig0(i)
  458.          go to 10
  459. c       on passe au rt d'apres (on sort de l) pour les rt autres que la premiere fissure
  460.          end if         
  461.  
  462. c stockage des variables fictives du pas        
  463.          varft(1,l)=sigf(1,r,l)
  464.          varft(2,l)=sigf(2,r,l)
  465.          varft(3,l)=sigf(3,r,l) 
  466.          varft(4,l)=vwfi(1,l)
  467.          varft(5,l)=vwfi(2,l)
  468.          varft(6,l)=vwfi(3,l)        
  469.          end do 
  470. c        fin de boucle sur l 
  471.    10    continue
  472. C               if(i.eq.3) then              
  473. C         print*,sigf(i,r,l),r,l,rig0(i,r),rt(r),si1
  474. C         end if  
  475.         end do 
  476. c      fin de boucle sur les Rtb  
  477.  
  478.  
  479. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  480. c on recommence une boucle de fissures fictives pour avoir la loi totale
  481. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  482.        nct=0.d0
  483.        do l=1,ncalc    
  484. c récupération de la contrainte du pas fictif d'avant
  485.        if(l.eq.1) then
  486.        sig0(1)=0.d0
  487.        sig0(2)=0.d0
  488.        sig0(3)=0.d0
  489.        vwfi0(1)=0.d0
  490.        vwfi0(2)=0.d0
  491.        vwfi0(3)=0.d0 
  492.        wmi0(1)=0.d0 
  493. c       wt(i,l-1)=0.d0
  494. c       sigt(i,l-1)=0.d0
  495.        else
  496.        sig0(1)=varft(1,l-1)
  497.        sig0(2)=varft(2,l-1) 
  498.        sig0(3)=varft(3,l-1) 
  499.        vwfi0(1)=varft(4,l-1)
  500.        vwfi0(2)=varft(5,l-1)
  501.        vwfi0(3)=varft(6,l-1) 
  502.        wmi0(1)=varft(7,l-1)        
  503.        end if
  504. c fissure fictive totale cette fois ci
  505.        vwfi(i,l)=l*incremwlis   
  506.  
  507.       call endo3f(vwfi(i,l),gf,young,
  508.      # xmt,dtiso,dt3,epict,rt(1),le(i))       
  509.  
  510. c       wk1=(Gf-gfap(i))/Rt(1)
  511. c       wk1=(Gf)/Rt(1)
  512.  
  513.       if(l.eq.1) then
  514.       epse13(i)=vwfi(i,l)/le(i)
  515.       else
  516.       epse13(i)=wt(i,l-1)/le(i)
  517.       end if
  518.  
  519.       sigfe(i)=(1.d0-dt3)*
  520.      #rhof*max(min(Ef*epse13(i)*cos(phimoy(i)*
  521.      #pi/180.d0)**2,fu*cos(phimoy(i)*pi/180.d0)**2),-fu)
  522.  
  523. *      sigfe(i)=0.d0      
  524.  
  525.        sigb(i)=Rt(1)*(1.d0-dt3)    
  526.        sig(i)=sigb(i)+sigf(i,1,l)*dt3+
  527.      #sigfe(i)
  528.  
  529. c      print*,sig(i),sigb(i),sigf(i,1,l),dt3,sigfe(i),"precons"
  530.  
  531.        Sr(i)=lech/exp(log(sig(i)/Rtec)*mw)  
  532.  
  533. c on fixe la limite   
  534.  
  535.        s1=Rtec*(lech*8.d0/(Lf*cos(pi*phimoy(i)/180.d0)))**(1/mw)
  536.  
  537.  
  538.  
  539. c      if(Sr(i).lt.Lf*cos(pi*phimoy(i)/180.d0)/4) then    
  540.        if(sig(i).gt.s1) then   
  541.        Sr(i)=Lf*cos(pi*phimoy(i)/180.d0)/4
  542. c       wmiloc(i)=vwfi(i,l)
  543.        sig(i)=s1
  544.        loc(i)=.true.
  545.        end if
  546.  
  547.        if(sig(i).lt.rt(2)) then 
  548.  
  549. c      il n y a qu une seule fissure
  550.        wfist(i)=vwfi(i,l)
  551.        nct(1)=1.d0
  552.        wmi(1)=vwfi(i,l)
  553.        wt(i,l)=wfist(i)
  554.        sigt(i,l)=sigf(i,1,l)
  555.        else 
  556. c      on multifissure
  557.  
  558.           do r=1,rmax 
  559.           if(r.eq.1) then
  560.           nct(r)=1.d0
  561.           wmi(1)=vwfi(i,l)
  562.           else
  563.               if(Sr(i).le.le(i)/(2.d0**(r-1))) then  
  564.                 if(le(i)/Sr(i).lt.2.d0**(r)) then
  565.                 nct(r)=(le(i)/Sr(i)-(2.d0**(r-1)))/2.d0
  566.                 else
  567.                 nct(r)=2.d0**(r-2)
  568.                 end if 
  569.               do k=1,ncalc   
  570.  
  571. c ici pour trouver les ouvertures de micro fissures on peut aussi equilibrer avec sigf(i,1,l) au lieu de sig(i)
  572. c ca paraitrait peut etre plus logique ?
  573. C               if(sig(i)-sigf(i,r,k).le.0.d0) then
  574. C               if(k.eq.1) then
  575. C        wmi(r)=(sig(i)-(sigf(i,r,k)-((sigf(i,r,k)-0)/
  576. C      #incremwlis)*vwfi(i,k)))/((sigf(i,r,k)-0)/incremwlis)             
  577. C               else 
  578. C       wmi(r)=(sig(i)-(sigf(i,r,k)-((sigf(i,r,k)-sigf(i,r,k-1))/
  579. C      #incremwlis)*vwfi(i,k)))/((sigf(i,r,k)-sigf(i,r,k-1))/incremwlis)
  580. C               end if
  581. C               go to 20               
  582. C               end if
  583.  
  584.               if(sigf(i,1,l)*dt3-sigf(i,r,k).le.0.d0) then
  585.               if(k.eq.1) then
  586.        wmi(r)=(sigf(i,1,l)*dt3-(sigf(i,r,k)-((sigf(i,r,k)-0)/
  587.      #incremwlis)*vwfi(i,k)))/((sigf(i,r,k)-0)/incremwlis)             
  588.               else 
  589.       wmi(r)=(sigf(i,1,l)*dt3-(sigf(i,r,k)-((sigf(i,r,k)-sigf(i,r,k-1))/
  590.      #incremwlis)*vwfi(i,k)))/((sigf(i,r,k)-sigf(i,r,k-1))/incremwlis)
  591.  
  592.  
  593.               end if
  594.               go to 20               
  595.               end if
  596.  
  597.               end do
  598.  
  599.   20          continue  
  600.  
  601. C               if(i.eq.3) then
  602. C               print*,wmi(r),r,nct(r),"pre",sr(i),sig(i),rt(r)
  603. C               end if
  604.  
  605.               if(wmi(r).lt.0.d0) then
  606.               print*,wmi(r)
  607.               read*
  608.               end if
  609.               wmi(r)=nct(r)*wmi(r)
  610.               else       
  611. c              wmi(r)=0.d0
  612.               go to 30
  613.               end if 
  614.         end if
  615.  
  616.  
  617.         end do 
  618. c        fin de boucle sur les Rt
  619.    30  continue        
  620. c       ici on va calculer des segments de droite pour reconstruire notre courbe
  621.         wfist(i)=sum(wmi)
  622.         wt(i,l)=wfist(i)
  623.         sigt(i,l)=sigf(i,1,l)
  624.  
  625.  
  626.        end if  
  627.  
  628. c       fin du si on multifissure      
  629.  
  630. c stockage de la rigidité initiale
  631.       if(l.eq.1) then        
  632.       rig0t(i)=sigf(i,1,l)/incremwlis
  633.       end if      
  634.  
  635.       if(l.gt.1) then
  636.       if(rig0t(i).ne.rig0t(i).or.(rig0t(i).eq.0.d0)) then
  637.       print*, "Probleme de rigidite initiale",rig0t(i)
  638.       err1=1
  639.       end if 
  640.       end if     
  641.  
  642.  
  643.       if(l.gt.1) then       
  644.       if((vwfi(i,l).gt.wpsig(i,1)).or.(loc(i)).or.(wt(i,l-1).gt.vw3(i)
  645.      # )) then 
  646.  
  647. c      print*, vwfi(i,l),wpsig(i,1),loc(i),wt(i,l-1),vw3(i),"w",vwfi0(i)
  648.  
  649. c      on calcule la contrainte jusqu'au pic uniquement
  650.       wpsigt(i)=vwfi0(i)
  651. c modif      
  652. C       if(wpsigt(i).lt.wpsigt0(i)) then
  653. C       wpsigt(i)=wpsigt0(i)
  654. C       end if
  655.       wmiloc(i)=wmi0(1)
  656.       sigmaxt(i)=sig0(i)
  657.       go to 40
  658. c      on sort de la boucle des deuxiemes fissures fictives l
  659.       end if    
  660.       end if      
  661.  
  662. c stockage des variables fictives du pas        
  663.       varft(1,l)=sigf(1,1,l)
  664.       varft(2,l)=sigf(2,1,l)
  665.       varft(3,l)=sigf(3,1,l) 
  666.       varft(4,l)=wfist(1)
  667.       varft(5,l)=wfist(2)
  668.       varft(6,l)=wfist(3)
  669.       varft(7,l)=wmi(1)
  670.  
  671.       end do 
  672. c      fin de boucle sur les fissures fictives l
  673.  
  674.    40 continue     
  675. c       print*,l,"l1"  
  676. cccccccccccccccccccccccccccccccccccccccccccccccc  
  677. c on passe au calcul reel ici 
  678. cccccccccccccccccccccccccccccccccccccccccccccccc
  679.       if(vw3(i).le.wpsigt(i)) then 
  680. ccccccccccccccccccccccccccccccccccccccccccccccccccc      
  681. c       on charge
  682. cccccccccccccccccccccccccccccccccccccccccccccccccccc
  683.       wmul=0.d0
  684.       wmr=0.d0 
  685.       wer=0.d0 
  686.       nc=0.d0
  687.  
  688. c on calcule la contrainte relle dans la fissure i et l'ouverture de la première fissure
  689. c on peut calculer l ouverture w1 car wt a été construit en imposant vwfi dans la premiere fissure
  690.       if(vw3(i).gt.0.d0) then
  691.  
  692.       do l=1,ncalc
  693.         if(wt(i,l)-vw3(i).gt.0.d0) then        
  694.            if(l.eq.1) then
  695.             sigfr(i)=(sigt(i,l)-0)/(wt(i,l)-0)*
  696.      #vw3(i)+sigt(i,l)-(sigt(i,l)-0)/(wt(i,l)-0)*
  697.      #wt(i,l)          
  698.  
  699.       w1(i)=vwfi(i,l)/wt(i,l)*vw3(i)     
  700.            else
  701.  
  702.            sigfr(i)=(sigt(i,l)-sigt(i,l-1))/(wt(i,l)-wt(i,l-1))*
  703.      #vw3(i)+sigt(i,l)-(sigt(i,l)-sigt(i,l-1))/(wt(i,l)-wt(i,l-1))*
  704.      #wt(i,l)
  705.  
  706.       w1(i)=(vwfi(i,l)-vwfi(i,l-1))/(wt(i,l)-wt(i,l-1))*vw3(i)+
  707.      #vwfi(i,l)-(vwfi(i,l)-vwfi(i,l-1))/(wt(i,l)-wt(i,l-1))*wt(i,l)     
  708.            end if
  709.            go to 50
  710.         end if      
  711.       end do
  712.  
  713.  
  714.  
  715.  
  716.  
  717.       else 
  718.       sigfr(i)=0.d0
  719.       sigb(i)=0.d0
  720.       end if
  721. 50    continue
  722. c       print*,l,"l2"
  723.  
  724.       call endo3f(w1(i),gf,young,
  725.      # xmt,dtiso,dt3,epict,rt(1),le(i))    
  726.  
  727. c calcul de la déformation totale de l element pour prendre en compte lacontrainteelastique des fibres       
  728.       call x6x33(epstf6,epse133)
  729.  
  730.       do j=1,3
  731.       vecwj(1)=vecw33(1,j)
  732.       vecwj(2)=vecw33(2,j)
  733.       vecwj(3)=vecw33(3,j)
  734.       epse13(j)=sum(matmul(vecw,epse133)*vecw)
  735.       sigfe(j)=(1.d0-dt3)*
  736.      #rhof*max(min(Ef*epse13(j)*cos(phimoy(j)*
  737.      #pi/180.d0)**2,fu*cos(phimoy(j)*pi/180.d0)**2),-fu)
  738.       end do
  739.  
  740. *       sigfe(i)=0.d0
  741.  
  742.        sigb(i)=Rt(1)*(1.d0-dt3)
  743.        sig(i)=sigb(i)+sigfr(i)*dt3+
  744.      #sigfe(i)
  745.  
  746. c      print*,sig(i),s1,"calc reel"
  747.  
  748.  
  749. c      print*,sig(i),sigb(i),sigfr(i),dt3,sigfe(i)
  750.  
  751.  
  752.        if(sig(i).gt.0.d0) then
  753.          Sr(i)=lech/exp(log(sig(i)/Rtec)*mw)
  754.         s1=Rtec*(lech*8.d0/(Lf*cos(pi*phimoy(i)/180.d0)))**(1/mw)
  755. c        print*,sr(i),Lf*cos(pi*phimoy(i)/180.d0)/4.d0
  756.         if(sig(i).gt.s1) then 
  757.          sig(i)=s1
  758.          Sr(i)=Lf*cos(pi*phimoy(i)/180.d0)/4.d0
  759.          wpsigt(i)=vw3(i)
  760.          sigmaxt(i)=min(sigfr(i),sig(i))
  761. c         sigmaxt(i)=sig(i)
  762. c          sigmaxt(i)=sigf03(i)
  763. c         print*,sigfr(3),sig(3),sigfe(i),sigb(i)
  764.  
  765.          end if
  766.        else
  767.        Sr(i)=le(i)
  768.        end if
  769.  
  770. c ici l ouverture de la fissure 1 est calculée de maniuere exacte
  771. c mais pour les autres on suppose que l endo du beton est le meme pour tous     
  772.       if(sig(i).lt.rt(2)) then 
  773.          nc(1)=1.d0
  774.          wmul(i,1)=vw3(i)
  775.        else 
  776. c       on multifissure
  777.         nc(1)=1.d0
  778.         wmul(i,1)=w1(i)
  779.           do r=2,rmax 
  780.  
  781.  
  782.               if(Sr(i).le.le(i)/(2.d0**(r-1))) then
  783.                 if(le(i)/Sr(i).lt.2.d0**(r)) then
  784.                 nc(r)=(le(i)/Sr(i)-(2.d0**(r-1)))/2.d0
  785.                 else
  786.                 nc(r)=2.d0**(r-2)
  787.                 end if 
  788.                 if(r.eq.1) then
  789.                 nc(r)=1.d0
  790.                 end if               
  791.               do k=1,ncalc
  792. c ici sig(i) est remplacé par sigfr(i)   
  793. c ca ne change apparemment rien           
  794.                 if(sigfr(i)*dt3-sigf(i,r,k).lt.0.d0) then            
  795.                  if(k.eq.1) then
  796.        wmul(i,r)=sigfr(i)*dt3*incremwlis/sigf(i,r,k)                
  797.                  else             
  798.       wmul(i,r)=(sigfr(i)*dt3-(sigf(i,r,k)-((sigf(i,r,k)-sigf(i,r,k-1))/
  799.      #incremwlis)*vwfi(i,k)))/((sigf(i,r,k)-sigf(i,r,k-1))/incremwlis)
  800. c fin remplacement
  801.                  end if
  802.                 go to 60               
  803.                 end if
  804.               end do
  805.   60   continue 
  806.  
  807. c              print*,k,"k2",r,"r" 
  808.  
  809.  
  810.  
  811.               else                
  812.              wmul(i,r)=0.d0
  813.              nc(r)=0.d0
  814.              if(r.eq.1) then
  815.              nc(1)=1.d0
  816.              wmul(i,1)=vw3(i)
  817.              end if
  818.               end if
  819. C         if(i.eq.3) then              
  820. C         print*,r,nc(r),sr(i),"calc",sig(i),s1
  821. C         end if
  822.         end do 
  823. c        fin de boucle sur les Rt            
  824.       end if  
  825. c fin du si sig(i).lt.rt(2)  
  826.  
  827.  
  828.  
  829.        if(vw3(i).gt.0.d0) then      
  830. c calcul de l'ouverture moyenne
  831.        do r=1,rmax
  832.  
  833.        wmr(r)=nc(r)*wmul(i,r)
  834.        end do    
  835. c       print*,wmr,nc,"nc",sig(i),rt(2)
  836. c       print*,Sr(i),le(i)/(2.d0**(1-1)),"test1"
  837. c       print*,le(i)/Sr(i),2.d0**(1),"test2"
  838.        wmoy(i)=sum(wmr)/sum(nc)  
  839.  
  840.        if(wmoy(i).ne.wmoy(i)) then
  841.        print*,"pb wmr,nc",sum(wmr),sum(nc),rmax,sig(i),rt(2),Sr(i),
  842.      #le(i)
  843.        print*,"elements finis trop petits"
  844.        read*
  845.        end if       
  846. c calcul de l'ecart type d ouverture      
  847.        do r=1,rmax
  848.         if(nc(r).eq.0.d0) then
  849.         wer(r)=0.d0
  850.         else
  851.         wer(r)=nc(r)*(wmul(i,r)-wmoy(i))**2
  852.         end if
  853.        end do
  854.        Etw(i)=sqrt(sum(wer)/sum(nc))
  855.        cvar(i)=etw(i)/wmoy(i)
  856.        else
  857.        wmoy(i)=0.d0
  858.        etw(i)=0.d0
  859.        wmr=0.d0
  860.        wer=0.d0
  861.        cvar(i)=0.d0
  862.        end if
  863. c ouverture max 
  864.        wmax(i)=wmul(i,1) 
  865.  
  866. c nombre total de fissures
  867.        ncf(i)=sum(nc)
  868.  
  869.        if(ncf(i).lt.nc0(i)) then
  870.        ncf(i)=nc0(i)
  871.        end if
  872.  
  873.  
  874.  
  875.         if(sig(i).ge.s1) then 
  876.          go to 80 
  877.          end if           
  878.  
  879.  
  880. c sauvegarde de la contrainte de decharge    
  881. c      sigdec(i)=sigf03(i)
  882.       sigdec(i)=sigfr(i)
  883.  
  884.  
  885.  
  886.  
  887.  
  888. **  AM       
  889. **       else 
  890.         endif
  891. ccccccccccccccccccccccccccccccccccccccccc      
  892. c       on localise la fissure en ouvrant
  893. ccccccccccccccccccccccccccccccccccccccccc
  894. 80     continue
  895.  
  896. **  AM
  897.  
  898.        if( ((VW3(I).LE.WPSIGT(I)).and.(SIG(I).GE.S1)) . OR.
  899.      &     (VW3(I).GT.WPSIGT(I)) ) then 
  900.        iloc(i)=1
  901.  
  902. c le nombre fissure ne peut plus augmenter
  903.        ncf(i)=nc0(i)
  904. c seule la fissure 1 continue de s ouvrir les autres ne se dechargent pas...       
  905.        wfis(i)=wmiloc(i)+(vw3(i)-wpsigt(i))
  906. C        if(vw3(i).lt.wpsigt(i)) then
  907. C        wfis(i)=vw3(i)
  908. C        end if
  909.  
  910.  
  911. C et donc ici on peut recalculer la contrainte comme provenant d'une mono fissure
  912. c rt a bien ete construit dans tous les cas       
  913.        Rtb=Rt(1)
  914.         do n=1,npos
  915.           do k=1,2*npos-2               
  916. c calcul de la force moyenne pour une fibre traversant la fissure i inclinée de phi1d
  917. c calcul des points caractéristiques de la courbe force ouverture de fissure pour Td et phid
  918. c rigidité initiale
  919.       if(phi1d(n,k).le.phicrit) then 
  920.       call reconf(cokf1,k0m,Rtb,phi1d(n,k),mink01,maxk01)
  921.       else
  922.       call reconf(cokf2,k0m,Rtb,phi1d(n,k),mink02,maxk02)
  923.       end if       
  924. c ouverture au pic
  925.       if(phi1d(n,k).le.phicrit) then 
  926.       call reconf(cowp1,wp,Rtb,phi1d(n,k),minwp1,maxwp1)
  927.       else
  928.       call reconf(cowp2,wp,Rtb,phi1d(n,k),minwp2,maxwp2)
  929.       end if        
  930. c force au pic
  931.       call reconf(cofp,fp,Rtb,phi1d(n,k),minfp,maxfp)
  932. c force intermédiaire post pic     
  933.       f03=0.3d0*fp        
  934. c ouverture pour f=0.3fp  
  935.       call reconf(cow03,w03,Rtb,phi1d(n,k),minw03,maxw03)
  936.         if(w03.le.wp) then        
  937.         w03=3.0d0*wp  
  938.         end if  
  939. c ouverture finale (f=0)    
  940.       call reconf(cowf,wf,Rtb,phi1d(n,k),minwf,maxwf)
  941.         if(wf.le.w03) then        
  942.         wf=3.0d0*w03  
  943.         end if  
  944. c calcul de la force dans la fissure i pour les fibres orientées selon ei     
  945.        if(wfis(i).le.wp) then
  946.        B=(k0m*wp-fp)/(fp*wp)
  947.        f1(n,k)=k0m*wfis(i)/(B*wfis(i)+1.d0)
  948.        else 
  949.          if(wfis(i).lt.wf) then
  950.          f1(n,k)=fp*exp(log(f03/fp)*(wfis(i)-wp)/(w03-wp))
  951. c        print*,2,f1(n,k),wp,wf
  952.          else
  953.          f1(n,k)=precision1
  954.          end if
  955.        end if
  956. c calcul du terme à l'intérieur de la somme de la contrainte
  957.         si0(n,k)=Reiw(n,k)*stri(n,k)*pscal(n,k)*f1(n,k)
  958.  
  959.          end do 
  960. c         fin de boucle sur k      
  961.         end do 
  962. c        fin de boucle sur n 
  963. c  calc de la contrainte 
  964.       si1=sum(si0)   
  965. c      wk1=(Gf-gfap(i))/Rt(1)
  966. c             wk1=(Gf)/Rt(1)
  967.       sigfr(i)=4.d0*rhof*si1/(pi*Df**2*(sum(Romeg)))
  968.  
  969.       if(sigfr(i).gt.sigmaxt(i)) then
  970.       sigfr(i)=sigmaxt(i)
  971.       end if
  972.  
  973.  
  974.        if(wfis(i).gt.0.d0) then      
  975. c calcul de l'ouverture moyenne
  976. C        do r=1,rmax
  977. C        nc(r)=ncr0(r)
  978. C          if(r.eq.1) then
  979. C          wmr(1)=wfis(i)
  980. C          else
  981. C          wmr(r)=nc(r)*wmul(i,r)
  982. C          end if
  983. C        end do    
  984. C        wmoy(i)=sum(wmr)/(sum(nc)) 
  985.          dw3(i)= vw3(i)-vw30(i)
  986.  
  987. c calcul de l'ecart type d ouverture      
  988. C        do r=1,rmax
  989. C          if(r.eq.1) then
  990. C          wer(1)=(wfis(i)-wmoy(i))**2
  991. C          else
  992. C          wer(r)=nc(r)*(wmul(i,r)-wmoy(i))**2
  993. C          end if
  994. C        end do
  995. C        Etw(i)=sqrt(sum(wer)/sum(nc))
  996.          if(nc0(i).eq.1.d0) then
  997.          etw(i)=0.d0
  998.          wmoy(i)=vw3(i)
  999.          else       
  1000. c         wmoy(i)=wmoy0(i)+dw3(i)/nc0(i)   
  1001.          wmoy(i)=wmoy0(i)      
  1002. C on conserve le coefficient de variation lors de la localisation
  1003.          etw(i)=etw0(i)
  1004.          end if
  1005.  
  1006.        else
  1007.        wmoy(i)=0.d0
  1008.        etw(i)=0.d0
  1009.        wmr=0.d0
  1010.        wer=0.d0
  1011.        end if
  1012. c ouverture max 
  1013.        wmax(i)=wfis(i) 
  1014. c      print*,wmax(i),i,"loc",vw3(i),vw30(i),wmaxact(i),rmax,le(i)
  1015. c      print*, wpsigt(i),sigmaxt(i),sum(Romeg),sigfr(i),wfis(i),nc0(i),
  1016. c     # wmiloc(i)
  1017. c       print*,"loc",dw3(i),vw30(i),etw(i),etw0(i),wmoy(i),i
  1018.  
  1019.  
  1020. c sauvegarde de la contrainte pour decharge
  1021. c      sigdec(i)=sigf03(i)
  1022.       sigdec(i)=sigfr(i)
  1023.  
  1024. c      print*,sigdec(i),i,vw3(i),vw30(i),wmaxact(i)
  1025.  
  1026.       if(vw3(i).le.wferm(i)) then
  1027.       go to 100
  1028.       end if      
  1029.  
  1030. C       if(sigfr(i).lt.0.d0) then
  1031. C       print*,sigfr(i),vw3(i),vw30(i),wmaxact(i),wfis(i),
  1032. C      #wmiloc(i),vw3(i),wpsigt(i)
  1033. C       read*
  1034. C       end if
  1035. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc       
  1036.        end if   
  1037. c       fin calcul avec w qui augmente    
  1038. ccccccccccccccccccccccccccccccccccccccccccccccccccccccc     
  1039.  
  1040.  
  1041. **  AM       
  1042. **      else 
  1043.         endif  
  1044.  
  1045. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
  1046. c      on a vw3(i).lt.wmaxact(i) donc on a fermé la fissure on est dans une phase décharge
  1047. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
  1048.  
  1049. 100   continue  
  1050. **  AM  
  1051.               IF(    (VW3(I).LT.WFERM(I)) .OR. 
  1052.      &             (VW3(I).LT.WPLMAX(I))) THEN
  1053.  
  1054.  
  1055. C       if(vw30(i).ge.vw3(i)) then
  1056. C       sigfr(i)=sigf03(i)-rig0t(i)*(vw30(i)-vw3(i))
  1057. C       else
  1058. C       print*,"comme par hasard",vw30(i),vw3(i)
  1059. C       read*
  1060. C       if(sigf03(i)-rig0t(i)*(vw30(i)-vw3(i)).gt.sigf03(i)) then
  1061. C       sigfr(i)=sigf03(i)
  1062. C       else
  1063. C       sigfr(i)=sigf03(i)-rig0t(i)*(vw30(i)-vw3(i))
  1064. C       end if
  1065. C       end if
  1066.  
  1067. c modif
  1068. c      sigfr(i)=sigf03(i)-rig0t(i)*(wmaxact(i)-vw3(i))
  1069.  
  1070.  
  1071.       if(vw30(i).ge.vw3(i)) then
  1072.       sigfr(i)=sigf03(i)-rig0t(i)*(vw30(i)-vw3(i))
  1073. c      else if(wmaxact(i).ge.vw3(i)) then
  1074. c      sigfr(i)=sigf03(i)-rig0t(i)*(wmaxact(i)-vw3(i))
  1075.       else if(sigf03(i)-rig0t(i)*(vw30(i)-vw3(i)).ge.sigf03(i)) then
  1076.       sigfr(i)=sigf03(i)-rig0t(i)*(vw30(i)-vw3(i))
  1077. c      else if(sigf03(i)-rig0t(i)*(wmaxact(i)-vw3(i)).ge.sigf03(i)) then
  1078. c      sigfr(i)=sigf03(i)-rig0t(i)*(wmaxact(i)-vw3(i))
  1079.       end if
  1080.  
  1081. c      print*,sigfr(i),sigf03(i),sigdec(i)
  1082.  
  1083.  
  1084.       if((sigfr(i).lt.0.d0).and.(sigf03(i).gt.0.d0)) then
  1085.       sigfr(i)=0.d0
  1086.       if(vw3(i).gt.wferm(i)) then
  1087.       wferm(i)=vw3(i)
  1088.       end if
  1089.       end if
  1090.  
  1091.       if(sigfr(i).lt.0.d0) then
  1092.       sigfr(i)=0.d0
  1093.       else if(sigfr(i).gt.sigdec(i)) then
  1094.       sigfr(i)=sigdec(i)      
  1095.       end if
  1096.  
  1097.       if(vw3(i).lt.wferm(i)) then
  1098.       sigfr(i)=0.d0
  1099.       end if 
  1100.  
  1101.  
  1102.       ncf(i)=nc0(i)
  1103.       dw3(i)= vw3(i)-vw30(i)
  1104.  
  1105.  
  1106.  
  1107.       if(nc0(i).eq.1.d0) then
  1108.       etw(i)=0.d0
  1109.       wmoy(i)=vw3(i)
  1110.       wmax(i)=vw3(i)
  1111.       else
  1112.       wmoy(i)=wmoy0(i)+dw3(i)/nc0(i)
  1113. c      Etw(i)=etw0(i)/wmoy0(i)*wmoy(i)
  1114.       Etw(i)=cvar(i)*wmoy(i)
  1115.       end if
  1116.       wmax(i)=wmax0(i)+dw3(i)/nc0(i)
  1117.       if(wmax(i).lt.0.d0) then
  1118.       wmax(i)=0.d0
  1119.       end if
  1120.  
  1121.       end if
  1122.  
  1123. 70    continue 
  1124.  
  1125.        if(ncf(i).lt.nc0(i)) then
  1126.        print*,"le nombre de fissures a diminue"
  1127.        print*,ncf(i),nc0(i),i,vw3(i),vw30(i),wmaxact(i)
  1128.        stop
  1129.        end if
  1130.  
  1131.       if(wmax(i).ge.wmaxd0(i)) then
  1132.       wmaxd(i)=wmax(i)
  1133.       else
  1134.       wmaxd(i)=wmaxd0(i)
  1135.       end if
  1136.       end do 
  1137. c      fin de boucle sur les fissures 
  1138.  
  1139. C       print*,sigfr(1),sigfr(2),sigfr(3),"sig"
  1140. C       print*,vw3(1),vw3(2),vw3(3),"vw"
  1141. C       print*,ncf(1),ncf(2),ncf(3),"nc"
  1142. C       print*,etw(1),etw(2),etw(3),"etw"
  1143.  
  1144.       if(etw(1).ne.etw(1).or.etw(2).ne.etw(2).or.etw(3).ne.etw(3)) then
  1145.       print*,"probleme ecart type fissure" ,etw(1),etw(2),etw(3)   
  1146.       print*,Etw0(1),vw3(1),dw3(1),wmoy(1),nc0(1),wmaxact(1),wmoy0(1),
  1147.      #ncf(1),iloc(1)
  1148.       print*,Etw0(2),vw3(2),dw3(2),wmoy(2),nc0(2),wmaxact(2),wmoy0(2),
  1149.      #ncf(2),iloc(2)
  1150.       print*,Etw0(3),vw3(3),dw3(3),wmoy(3),nc0(3),wmaxact(3),wmoy0(3),
  1151.      #ncf(3),iloc(3)   
  1152.       read*
  1153.       end if
  1154.  
  1155. C       if(nc(3).gt.16.d0) then
  1156. C       print*,'pb'
  1157. C       read*
  1158. C       end if
  1159.  
  1160. c  reconstruction du tenseur des contraintesde fibres avec endo dans la base fixe
  1161.       do i=1,3
  1162.         do j=1,3
  1163.         if(i.eq.j) then 
  1164.         sigfp33(i,j)=sigfr(i)
  1165.         else
  1166.         sigfp33(i,j)=0.d0
  1167.         end if
  1168.         end do
  1169.       end do    
  1170. c      sigfx33=matmul(matmul(transpose(vecw33),sigfp33),vecw33)
  1171. c     vecw33 est la matrice de passage de la base maillage a fissuration
  1172. c     ici on veut faire l inverse cad : P A tP
  1173.       sigfx33=matmul(matmul(vecw33,sigfp33),transpose(vecw33))
  1174.       call x33x6(sigfx33,sigfx6)
  1175.  
  1176. C       print*,"wmax",wmax
  1177. C       print*,"wmoy",wmoy
  1178. C       print*,"etw",etw
  1179. C       print*,"ncf",ncf
  1180. C       if(sum(sigfx6*sigfx6).gt.0.d0) then
  1181. C        print*,"sigfx6",sigfx6
  1182. C        end if
  1183.  
  1184. C          print*,"sig(3)",sig(3),"sigfr",sigfr(3),s1
  1185. C      #,"dt3","wmaxd",wmaxd(3),dt3
  1186.       return
  1187.  
  1188.       end
  1189.  
  1190.  
  1191.  
  1192.  
  1193.  

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