Télécharger la plaquette Cast3M


fr en

Télécharger lcls3d.eso

Retour à la liste

Numérotation des lignes : 

lcls3d

  1. C LCLS3D    SOURCE    FD218221  24/02/07    21:15:17     11834          
  2.        subroutine lcls3d(IETAP,NGF,ERR1,AA,XX,A2,BB,AAI,BBI,AAP,IPZERO,
  3.      # NBRINC,NINC,NDIMG,NDIMA,LOG_FLU,PLAST,ITER,DETH,DEPH,DEPST3,
  4.      # JEA,JSE,AK,BK,AM,EE0,EK0,DEPSE1,DEPSM1,DEPSK1,FRAC,DEPST1,
  5.      # NTMAX,FT,NBRACT,LNUMCRT,ACTIFT,PLASTT,DPFTDS,DGFTDS,REST,
  6.      # PFT_GRAD_PFT,PFT_GRAD_GFT,DEPSC1,NCMAX,NBRACC,FC,LNUMCRC,
  7.      # ACTIFC,PLASTC,DPFCDS,DGFCDS,RESC,PFC_GRAD_PFC,PFC_GRAD_GFC,
  8.      # DEPSL1,NLMAX,FL,NBRACL,LNUMCRL,ACTIFL,PLASTL,DPFLDS,DGFLDS,
  9.      # RESL,PFL_GRAD_PFL,PFL_GRAD_GFL,Gch,bg,dbgVg,dbgVd,dbwPw,
  10.      # reduc_resg,JSP,Je1epg,Mbg,NDIMPL,RgR,Jsgep,
  11.      # PFG_GRAD_PFG,PRECISION3D,AFFICHE,CONV_FORCEE)
  12.  
  13.  
  14.  
  15. c     calcul des deformation elastiques (A.Sellier 2021/02/26)
  16. c     modif AS 2022/02/10
  17. c     adaptation variables generalisees 30/03/2022
  18.  
  19. c     calcul des variables internes locales dans  inclusions et matrice
  20. c     construction de la matrice de localisation integrant le fluage
  21. c     resolution dans la base principale de l increment
  22.  
  23.       implicit real*8 (a-h,o-z)
  24.       implicit integer (i-n)
  25.  
  26. c     --- variables externes -------------------------------------------
  27.  
  28.  
  29. c     nbre maxi d inclusions      
  30.       integer NBRINC,NINC,NDIMG,NDIMA,NTMAX,NCMAX,NLMAX
  31. c     nbre de type de criteres de plasticite
  32.       integer NBTYPL      
  33.       logical AFFICHE
  34.       integer IETAP,ngf,err1
  35.       real*8 AA(NGF,NGF+1),BB(NGF),XX(NGF),A2(NGF*NGF)
  36.       real*8 AAI(NGF,NGF+1),BBI(NGF)
  37.       integer ipzero(ngf)
  38.  
  39. c     dimension de l espace de minimisation des residus plastiques
  40.       integer NDIMPL
  41. c     matrice d bgpg / d ea ds l espace de minimisation
  42.       real*8 Jsgep(NDIMG,NDIMPL)
  43. c     Vecteur R.grad(R) dans l espace de minimisation
  44.       real*8 RgR(NDIMPL)
  45. c     Module de Biot des gels Mbg(numero phase , debut ou fin de pas)
  46.       real*8 Mbg(0:NBRINC,0:1)  
  47.  
  48. c     preparation de la matrice inv(A).Jea
  49.       real*8 AAP(NDIMG,NDIMPL)
  50. c     matrice jsp=dseff g/d eps p (base ecoulement plastique)
  51.       real*8 jsp(NDIMG,NDIMPL)
  52.  
  53.       real*8 FRAC(0:NBRINC)
  54. c     coeff de biot des gels
  55.       real*8 bg(0:NBRINC,0:1)      
  56. c     fluage
  57.       real*8 ak(NDIMG),bk(NDIMG),am(NDIMG)
  58.       real*8 ek0(NDIMG),ee0(NDIMG)   
  59.       logical LOG_FLU,PLAST
  60.       integer ITER
  61. c     thermo hydro chimie      
  62.       real*8 DETH(0:NBRINC),DEPH(0:NBRINC)    
  63. c     indicateur de diffusion produits chimiques dans fissure
  64.       real*8 GCH(0:NBRINC)
  65. c     increment de deformation homogeneisee
  66.       real*8 DEPST3(3) 
  67. c     calcul visco elastique forcee
  68.       logical CONV_FORCEE   
  69. c     donnes poro-mechaniques (eau, gel, dechrage)      
  70.       real*8 dbwPw(0:NBRINC),dbgVg(0:NBRINC),dbgVd(0:NBRINC)
  71.  
  72.  
  73. c     derivee du critere par rapport à la contrainte
  74.       real*8 DPFTDS(NTMAX,NDIMG),DPFCDS(NCMAX,NDIMG)  
  75. c     derivee du critere par rapport à la contrainte
  76.       real*8 DGFTDS(NTMAX,NDIMG),DGFCDS(NCMAX,NDIMG)      
  77. c     activite des criteres      
  78.       logical ACTIFT(NTMAX),ACTIFC(NCMAX)
  79. c     table du critere de rankine et dp  
  80.       real*8 FT(NTMAX),FC(NCMAX)
  81. c     numeror critre retenu pour les tensions refermeture      
  82.       integer LNUMCRT(NTMAX),LNUMCRC(NCMAX)
  83. c     gradient pondéré des fonctions seuil et pseudo potentiels pour les critères de  Traction diffuse
  84.       real*8 PFT_GRAD_PFT(NDIMPL),PFT_GRAD_GFT(NDIMPL) 
  85. c     gradient pondéré des fonctions seuil et pseudo potentiels pour les critères de Cisaillement diffus
  86.       real*8 PFC_GRAD_PFC(NDIMPL),PFC_GRAD_GFC(NDIMPL) 
  87. c     gradiant pondéré des fonctions pour les criteres de traction localisee aubords extérieurs des phases
  88.       real*8 PFL_GRAD_PFL(NDIMPL),PFL_GRAD_GFL(NDIMPL)
  89. c     par de fgradf due a la pression de gel pour les criteres de decollement
  90.       real*8 PFG_GRAD_PFG(NDIMPL)
  91. c     coeff de relaxation du residu
  92.       real*8 reduc_resg 
  93. c     jacobienne Je1epg dans la base des ecoulements plastiques
  94.       real*8 Je1epg(NDIMG,NDIMPL)      
  95.  
  96.  
  97. c     variables dans espace des variables generalise
  98. C     EPS(1-3)=INCLUSION
  99. C     EPS(4-6)=INTERFACE RADIALE
  100. C     EPS(7-9)=INTERFACE ORTHORADIALE
  101. C     EPS(10-12)=INFINI MATRICE    
  102.  
  103.  
  104. c     Res vecteur des residus des criteres
  105.       real*8 Rest,resc,resl
  106.  
  107. c     vecteur deformations solutions       
  108.       real*8 DEPSE1(NDIMG),DEPSM1(NDIMG),DEPSK1(NDIMG)
  109. c     vecteur solution def anelastiques      
  110.       real*8 DEPST1(NDIMG),DEPSC1(NDIMG)
  111. c     precision 
  112.       real*8 PRECISION3D
  113. c     jacobienne depse/depsimp si un seul type dinclusion
  114.       real*8 JEA(NDIMG,NDIMA)
  115. c     jacobienne d sigma / d epse      
  116.       real*8 JSE(NDIMG,NDIMG)
  117. c     indicateur de plasticite elementaire      
  118.       logical PLASTT,PLASTC,PLASTL 
  119. c     nombre de criteres actif, numero d ordre du dernier traitement 
  120.       integer NBRACT
  121. c     nombre de criteres actif cisaillement, numero d ordre du dernier traitement 
  122.       integer NBRACC
  123. c     nombre de citeres localises actifs
  124.       integer NBRACL
  125. c     criteres de fissures localisees
  126.       real*8 FL(NLMAX),DEPSL1(NDIMG)
  127.       integer LNUMCRL(NLMAX)
  128.       logical ACTIFL(NLMAX)
  129.       real*8 DPFLDS(NLMAX,NDIMG),DGFLDS(NLMAX,NDIMG)
  130.  
  131.  
  132. c     tableau de classement des criteres 
  133.       integer NORDRE(NDIMG)      
  134.  
  135. c   -- variables locales -----------------------------------------------
  136.       real*8 detc0,detc1,determinant
  137.       integer i,j,NUMA,iordre,idir
  138.       logical affiche_local
  139.       integer NDIM_SOLVE   
  140.  
  141. c     -- intialisations ------------------------------------------------
  142.  
  143.       ERR1=0      
  144.       affiche_local=affiche
  145. c      affiche_local=.true.
  146.       if(affiche_local) then
  147.         print*,'On arrive dans lcls3d'
  148.       end if    
  149.  
  150. c     -- localisation de la deformation imposee ------------------------  
  151.  
  152.       if((IETAP.eq.0).and.(.not.LOG_FLU)) then
  153.  
  154. c     ----------TIR ELASTIQUE (IETAP=0) --------------------------------
  155.  
  156. c        tir elastique calcul direct des depse sans fluage
  157.          if(affiche_local) then
  158.              write(*,'(/,A50,1X,A6,I2)') 
  159.      #       'Calcul direct des DEPSE dans lcls3d', 
  160.      #       'IETAP=',ietap
  161.          end if
  162. c        prise en compte des conditions imposees sur epse generalise et MK
  163.          call elin3d(NBRINC,NDIMG,DETH(1),DEPH(1),DETH(0),DEPH(0),
  164.      #   DEPST3,dbwpw,dbgVg,JEA,DEPSE1,affiche_local,err1)
  165.          do i=1,NDIMG
  166. c           vecteur des inconnues en fluage si une seule inclusion
  167.             DEPSM1(i)=0.d0
  168. c           vecteur des inconnues en fluage si une seule inclusion
  169.             DEPSK1(i)=0.d0                   
  170. c           vecteur des inconnues plastique traction si une seule inclusion           
  171.             DEPST1(i)=0.d0    
  172. c           vecteur des inconnues plastique compression si une seule inclusion           
  173.             DEPSC1(i)=0.d0
  174. c           vecteur des inconnues plastique localisee si une seule inclusion           
  175.             DEPSL1(i)=0.d0             
  176.          end do
  177. c        pas d ecoulement des gels car les fissures n ont pas evoluees        
  178.          dbgVd(0)=0.d0
  179.          dbgVd(1)=0.d0
  180.  
  181.       else if ((IETAP.eq.0).and.LOG_FLU) then
  182.  
  183. c     ----------TIR VISCO-ELASTIQUE (IETAP=0) --------------------------
  184.  
  185.          if(affiche_local) then
  186.              write(*,'(/,A50,1X,A6,I2)') 
  187.      #       'Calcul visco-elastique des DEPSE dans lcls3d', 
  188.      #       'IETAP=',ietap
  189.          end if      
  190. c        construction du systeme lineaire NDIMG*NDIMG pour tir visco elastique
  191.          call a1ti3d(NBRINC,NDIMG,AK,BK,AM,JEA,NGF,AA,ERR1,
  192.      #    AFFICHE_local)
  193. c        construction du second membre
  194.          call bbti3d(NBRINC,NDIMG,ak,bk,am,DETH(1),DEPH(1),DETH(0),
  195.      #   DEPH(0),ek0,ee0,DEPST3,dbwPw,dbgVg,JEA,ngf,BB,ERR1,
  196.      #   affiche_local)     
  197. c        resolution du systeme lineaire             
  198.          call gaus3d(NDIMG,aa,xx,bb,ngf,err1,ipzero)
  199.          if(err1.eq.1) then
  200.              NDIM_SOLVE=NDIMG
  201.              goto 20 
  202.          end if 
  203. c        recuperation des increments elastiques et fluage
  204.          do i=1,NDIMG
  205.             DEPSE1(i)=xx(i)
  206. c           vecteur des inconnues en fluage si une seule inclusion
  207.             DEPSM1(i)=am(i)*(ee0(i)+depse1(i)/2.d0)
  208. c           vecteur des inconnues en fluage si une seule inclusion
  209.             DEPSK1(i)=ak(i)*(ee0(i)+depse1(i))-bk(i)*ek0(i)                    
  210. c           vecteur des inconnues plastique traction si une seule inclusion           
  211.             DEPST1(i)=0.d0    
  212. c           vecteur des inconnues plastique compression si une seule inclusion           
  213.             DEPSC1(i)=0.d0
  214. c           vecteur des inconnues plastique localisee si une seule inclusion           
  215.             DEPSL1(i)=0.d0            
  216.          end do
  217. c        pas d ecoulement  des gels car les fissures n ont pas evoluees     
  218.          dbgVd(0)=0.d0
  219.          dbgVd(1)=0.d0       
  220.       else 
  221.  
  222. c     -------------- ELASTO VISCO PLASTICITE (IETAP=4) ----------------- 
  223.  
  224.          if(affiche_local) then
  225.              write(*,'(A50,1X,A6,I2)') 
  226.      #       'Calcul de l ecoulement plastique dans lcls3d', 
  227.      #       'IETAP=',ietap
  228.          end if
  229. C        construction du systeme lineaire pour maintenir les criteres
  230. c        actifs a zero lors du tir visco plastique
  231.  
  232. c        dimension du sous systeme a resoudre         
  233.          NDIM_SOLVE=NDIMG
  234.  
  235. c        matrice de couplage  AA=1-Je1ee dans le cas avec ecoulement visco
  236. c        elasto plastique en base generalisee des ecoulements plastiques       
  237.          call a1re3d(NBRINC,NDIMG,NDIMA,NGF,AA,AK,BK,AM,JEA,
  238.      #   AFFICHE_local,ERR1)
  239.  
  240. c        inversion de la matrice de couplage AAI*AA=I -> AAI=1/(1-Je1ee)    
  241.          call minv3d(NGF,AA,AAI,BBI,XX,IPZERO,NDIM_SOLVE,
  242.      #   ERR1,AFFICHE_local)
  243.  
  244. c        jacobienne de e1 par rapport à ep Je1ep ep en base  d ecoulement (3*12=32)
  245.          call jepg3d(NBRINC,NINC,NDIMG,NDIMA,NDIMPL,FRAC,GCH,Mbg,bg,
  246.      #   JEA,Je1epg,AFFICHE,ERR1)
  247.  
  248. c        jacobienne des pressions de gel par rapport aux deformatins plastique
  249.          call jbgp3d(NBRINC,NINC,NDIMG,NDIMPL,FRAC,GCH,Mbg,bg,
  250.      #   Jsgep,AFFICHE,ERR1)
  251.  
  252. c        --- calcul des ecoulement plastiques par methode du gradient -- 
  253.  
  254.          call denc3d(ngf,NDIMG,NDIMA,NTMAX,NCMAX,NLMAX,
  255.      #   nbract,ft,lnumcrt,dpftds,dgftds,REST,PFT_GRAD_PFT,PFT_GRAD_GFT,
  256.      #   nbracc,fc,lnumcrc,dpfcds,dgfcds,RESC,PFC_GRAD_PFC,PFC_GRAD_GFC,
  257.      #   nbracl,fl,lnumcrl,dpflds,dgflds,RESL,PFL_GRAD_PFL,PFL_GRAD_GFL,
  258.      #   jse,jea,aai,aap,jsp,Je1epg,DEPST1,DEPSL1,DEPSC1,reduc_resg,
  259.      #   NDIMPL,RgR,ninc,NBRINC,Jsgep,ACTIFT,ACTIFC,ACTIFL,PFG_GRAD_PFG,
  260.      #   precision3d,err1,affiche)
  261.  
  262.  
  263. c        --- actualisation du second membre  ---------------------------
  264.  
  265.          call bbre3d(NBRINC,NDIMG,NDIMA,ngf,BB,bg,DEPST1,
  266.      #   DEPSC1,DEPSL1,JEA,frac,Gch,affiche_local)    
  267.  
  268. c         --- resolution des deformations elastiques -------------------
  269.  
  270.          do i=1,NDIMG
  271.             xx(i)=0.d0
  272.             do j=1,NDIMG
  273.                 xx(i)=xx(i)+aai(i,j)*bb(j)
  274.             end do
  275.          end do
  276.  
  277.          if(err1.eq.1) then             
  278.              goto 20
  279.          else
  280.  
  281. c           --recuperation des increments elastiques et visco elastiques
  282.  
  283.             do i=1,NDIMG
  284. c               elastique
  285.                 DEPSE1(i)=xx(i)
  286. c               Maxwell
  287.                 DEPSM1(i)=am(i)*(depse1(i)/2.d0)
  288. c               Kelvin
  289.                 DEPSK1(i)=ak(i)*(depse1(i))                                
  290.             end do
  291.  
  292. c           ---- modification des volumes de decharge du gel -----------        
  293.  
  294.             call dbgv3d(NBRINC,NINC,NDIMG,GCH,FRAC,dbgVd,bg,
  295.      #      DEPST1,DEPSL1,ERR1,AFFICHE)
  296.  
  297.          end if                
  298.  
  299.       end if
  300.  
  301.  
  302. c     **** affichage fin de lcls3d *********************************       
  303.       if (affiche_local) then
  304.             do i=1,NDIMG
  305.                 if(DEPST1(i).lt.0.) then
  306.                    write(*,'(1X,A39,I2,1X,A9,I2,A2,E10.3)') 
  307.      #             'Dans lcls3d, on referme la fissure ',i,
  308.      #             ' et epst(',i,')=',depst1(i)
  309. c                  read*
  310.                 end if
  311.             end do
  312.       end if
  313.  
  314. c     **** traitement de l erreur enventuelle **************************       
  315.  
  316. 20    if(err1.eq.1) then
  317.             write(*,'(A36,/,3(A7,I2,1X,A6,I2,1X,A3,E10.3,/))')
  318.      #     'Erreur Gauss3d dans lcls3d avec:',
  319.      #     'NBRACT',NBRACT,'NUMCR',LNUMCRT(1),'FT',FT(LNUMCRT(1)),
  320.      #     'NBRACC',NBRACC,'NUMCR',LNUMCRC(1),'FC',FC(LNUMCRC(1)),
  321.      #     'NBRACL',NBRACL,'NUMCR',LNUMCRL(1),'FL',FL(LNUMCRL(1))
  322.             write(*,'(3(A10))') 'Numero','direction','FT'
  323.             do i=1,NTMAX
  324.                 write(*,'(2(I10),E10.3,L2)') i ,LNUMCRT(i),
  325.      #          FT(LNUMCRT(i)),ACTIFT(LNUMCRT(i))
  326.             end do
  327.             write(*,'(3(A10))') 'Numero','direction','FC'  
  328.             do i=1,NCMAX
  329.                 write(*,'(2(I10),E10.3,L2)') i ,LNUMCRC(i),
  330.      #          FC(LNUMCRC(i)),ACTIFC(LNUMCRC(i))
  331.             end do 
  332.             write(*,'(3(A10))') 'Numero','direction','FL'
  333.             do i=1,NLMAX
  334.                 write(*,'(2(I10),E10.3,L2)') i ,LNUMCRL(i),
  335.      #          FL(LNUMCRL(i)),ACTIFL(LNUMCRL(i))
  336.             end do        
  337.             write(*,'(3(A10))') 'AK','BK','AM'
  338.             do i=1,NDIMG
  339.                 write(*,'(3E10.3)') AK(i),BK(i),AM(i)
  340.             end do
  341.             print*, 'Matrice de couplage'
  342.             write(*,'(A3)',advance='no') 'AA'
  343.             do i=1,NDIM_SOLVE
  344.                 write(*,'(9X,I2)',advance='no') i
  345.             end do
  346.             write(*,'(A11)') 'BB'
  347.             do i=1,NDIM_SOLVE
  348.                 write(*,'(1X,I2)',advance='no') i
  349.                 do j=1,NDIM_SOLVE
  350.                     write (*,'(1X,E10.3)',advance='no') AA(i,j)
  351.                 end do
  352.                 write(*,'(1X,A1,1X,E10.3)') '=',BB(i)
  353.             end do                  
  354.             return
  355.       else 
  356. c       print*,'Dans lcls3d'
  357. c       do i=1,nt
  358. c           write(*,10) i,xx(i)
  359. c       end do
  360. c       read*
  361.         continue
  362.       end if 
  363.  
  364. c     *** detection des inconsistances plastiques **********************
  365. c      test determinant du systeme
  366. c      call detq3d(aa,a2,ngf,nt,determinant,ifault)
  367. c      print*,'Dans lcls3d Determinant du systeme', determinant
  368. c      read*    
  369.  
  370.       return
  371.       end 
  372.  
  373. c***********************************************************************
  374.  
  375.  
  376.  
  377.  
  378.  
  379.  

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