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

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