Télécharger la plaquette Cast3M


fr en

Télécharger coupl3d.eso

Retour à la liste

Numérotation des lignes : 

coupl3d

  1. C COUPL3D   SOURCE    FD218221  24/02/07    21:15:07     11834          
  2.         subroutine coupl3d(ann,xn,bn,ngf,err1,deps3,vp33,vp33t,dt,cc3,
  3.      #  taum,taumdtt,tauk,psi,epse06,epsk06,raideur66,Mg,bg,dphig,Cg,Hg,
  4.      #  Md,bd,dphid,Cd,Hd,Mw,bw,dphiw,Cw,alphadp,betadp,actif,factif,
  5.      #  nc,dsigef6,depse6,depsk6,depsm6,depspc6,depspt6,depsvt6,depspg6,
  6.      #  depsps6,tirve,ipzero,dpg,dps,dpw,epspt6p,refermtl,complet,na,CM)
  7.  
  8. c       Calcul des increments de deformations et de contraintes
  9. c       construction de la mattrice de couplage
  10. c       entre les fluages de Kelvin et de Maxwell, les deformations
  11. c       elastiques, totales, plastiques et les contraintes, les contraintes
  12. c       equivalentes et les multiplicateurs plastiques
  13.  
  14. c   ********************************************************************
  15.         implicit real*8 (a-h,o-z)
  16.         implicit integer (i-n)
  17.  
  18.         integer ngf,err1
  19.         real*8 ann(ngf,ngf+1),bn(ngf),xn(ngf)
  20.         integer ipzero(ngf)   
  21.  
  22.         real*8 deps3(3),dt
  23.         real*8 taukref,CM
  24.  
  25.         real*8 tauk,psi,taum,taumdtt,Mg,bg,Cg,Hg,dphig,Md,bd,dphid
  26.         real*8 Cd,Hd,Mw,bw,dphiw,Cw,dpg,dps,dpw
  27.         real*8 epse06(6),epsk06(6),raideur66(6,6),alphadp(3),betadp(3)
  28.         real*8 epspt6p(6)
  29.         integer nc
  30.         logical actif(nc),refermtl,complet(nc)
  31.         real*8 factif(nc)
  32.         real*8 vp33(3,3),vp33t(3,3),cc3(3)
  33.         logical tirve
  34.  
  35.  
  36.  
  37.         real*8 Jf(3,3),epse16(6),epsk16(6)
  38.         real*8 denomm,denomd,coeffm,coeffd,coeffg,coeffw  
  39.         integer errgauss
  40.         real*8 depse6p(6),depse6(6)
  41.         real*8 depsk6p(6),depsk6(6) 
  42.         real*8 depsm6p(6),depsm6(6)
  43.         real*8 depspc6p(6),depspc6(6)           
  44.         real*8 depspt6p(6),depspt6(6) 
  45.         real*8 depsvt6p(6),depsvt6(6)
  46.         real*8 depspg6p(6),depspg6(6)
  47.         real*8 depsps6p(6),depsps6(6) 
  48.         real*8 dsigef6p(6),dsigef6(6)  
  49.         integer i,j,k,l
  50.         integer ni,ncal,nmax,na       
  51.         logical reduc,recal
  52.         real*8 coeff
  53. c       compteur et numero des criteres tronquets en raison du seuil en contrainte        
  54.         integer nbrincg,numincg(3)
  55.         integer nbrincd,numincd(3)
  56.         integer nbrincw,numincw(3)
  57.         real*8 dx6(6),dx16(6)
  58.  
  59. c        print*,'verif coupl3d'
  60. c        if(actif(13)) then
  61. c         do i=1,3
  62. c          print*,alphadp(i),betadp(i)
  63. c         end do
  64. c        end if        
  65.         ncal=1
  66.         nmax=nc
  67. c       initialisation des tableaux pour la resolution
  68. 10      recal=.false.
  69. c       remise a zero de la matrice de couplage globale
  70.         do i=1,ngf
  71.             do j=1,(ngf+1)
  72.                 ann(i,j)=0.d0
  73.             end do
  74.             bn(i)=0.d0
  75.             xn(i)=0.d0
  76.         end do  
  77.  
  78. c       mise a zero des compteurs de criteres tronquets
  79.         nbrincg=0
  80.         nbrincd=0
  81.         nbrincw=0        
  82.  
  83. c       *** ordre des inconnues *(a partir de 25r suppression depsw)****
  84. c       v*=(deps-depspl)/dt 1-3  (si dt=0 on suppose v*=depsimp-depspl)
  85. c       et on imposse depsk=0,depse=v*dt, depsm=0)
  86. c       depse               4-6
  87. c       depsk               7-9
  88. c       depsm*              10-12
  89. c       depsm               13-15
  90. c       depstt              16-18
  91. c       dsigef              19-21
  92. c       dpg                 22
  93. c       dps                 23
  94. c       dpw                 24
  95. c       dsigm               25-27       criteres: 10-12
  96. c       dscrg               28-30       criteres: 1-3
  97. c       dscrs               31-33       criteres: 4-6
  98. c       dscrw               34-36       criteres: 7-9
  99. c       dsDP                37          critere : 13
  100. c       depsg               38-40
  101. c       depss               41-43
  102. c       depst               44-46
  103. c       depsc               47-49
  104. c       dlambdag            50-52       criteres: 1-3
  105. c       dlambdas            53-55       criteres: 4-6
  106. c       dlambdat            56-58       criteres: 10-12
  107. c       dlambdac            59          critere : 13
  108. c       ****************************************************************
  109.  
  110. c       passage des deformation initiales dans la base d increment
  111.         call chrep6(epse06,vp33,.false.,epse16)
  112.         call chrep6(epsk06,vp33,.false.,epsk16)
  113. c        print*,'coupl3d'
  114. c        do i=1,6
  115. c           print*,'epse16',i,epse16(i)
  116. c        end do
  117. c        do i=1,6
  118. c           print*,'epsk16',i,epsk16(i)
  119. c        end do
  120.  
  121. c       remplissage de la matrice de couplage dans une base principale
  122. c       bouclage sur les 3 directions principales de l increment
  123.         do i=1,3
  124.  
  125. c           *** vitesses virtuelles de deformations visco elastique ****
  126.              l=i   
  127. c            inconnues v*
  128.              k=l 
  129.              if(dt.gt.0.) then             
  130.                  ann(l,k)=dt
  131.              else
  132.                  ann(l,k)=1.d0
  133.              end if
  134. c            deformations imposees
  135.              bn(l)=deps3(i)
  136. c            couplage avec depsg
  137.              k=37+i
  138.              ann(l,k)=1.d0
  139. c            couplage avec depss
  140.              k=40+i
  141.              ann(l,k)=1.d0 
  142. c            couplage avec depst
  143.              k=43+i
  144.              ann(l,k)=1.d0 
  145. c            couplage avec depsc
  146.              k=46+i
  147.              ann(l,k)=1.d0 
  148.  
  149. c           *** deformations elastiques ********************************
  150.             l=3+i
  151. c           inconnues depse
  152.             k=l
  153.             ann(l,k)=1.d0 
  154.             if(dt.gt.0.) then
  155. c               couplage avec v*
  156. c               preparation du couplage visco elastique
  157. c               print*,tauk,psi,taum*cc3(i),taumdtt*cc3(i),dt,Jf
  158.                 call jacobflu3d(tauk,psi,taum*cc3(i),taumdtt*cc3(i),
  159.      #          CM,dt,Jf,err1)
  160.                 if(err1.eq.1) then
  161.                     print*,'Pb lors du calcul de Jf dans coupl3d'
  162.                     return
  163.                 end if
  164. c               couplage avec inconnues v*
  165.                 k=i            
  166.                 ann(l,k)=-Jf(1,1)
  167. c               couplage avec connues epse0 epsk0
  168.                 bn(l)=Jf(1,2)*epse16(i)+Jf(1,3)*epsk16(i)
  169.             else
  170. c               le pas de temps est nul la def virtuelle=elastique
  171.                 bn(l)=0.d0  
  172. c               couplage avec la deformation virtuelle
  173.                 k=i
  174.                 ann(l,k)=-1.d0
  175.             end if
  176.  
  177. c           *** deformations de Kelvin *********************************
  178.             l=6+i
  179. c           inconnues depsk
  180.             k=l
  181.             ann(l,k)=1.d0
  182.             if(dt.gt.0.) then
  183. c             couplages avec v*            
  184.               k=i            
  185.               ann(l,k)=-Jf(2,1)
  186. c             couplage avec  epse0 epsk0
  187.               bn(l)=Jf(2,2)*epse16(i)+Jf(2,3)*epsk16(i) 
  188.             else
  189. c             l increment de Kelvin est nul            
  190.               bn(l)=0.d0
  191.             end if              
  192.  
  193. c           *** deformation virtuelle de maxwell (em+edtt) *************
  194.             l=9+i
  195. c           inconnue depsm*
  196.             k=l
  197.             ann(l,k)=1.d0
  198.             if(dt.gt.0.) then
  199. c             couplages avec v*            
  200.               k=i            
  201.               ann(l,k)=-Jf(3,1)
  202. c             couplage avec  epse0 epsk0
  203.               bn(l)=Jf(3,2)*epse16(i)+Jf(3,3)*epsk16(i) 
  204.             else
  205. c             l increment de Kelvin est nul            
  206.               bn(l)=0.d0
  207.             end if 
  208. c            if(dt.gt.0.) then
  209. cc                couplage avec v*
  210. c                 k=i            
  211. c                 ann(l,k)=-dt
  212. cc                couplage avec depse
  213. c                 k=3+i
  214. c                 ann(l,k)=1.d0
  215. cc                couplage avec depsk
  216. c                 k=6+i
  217. c                 ann(l,k)=1.d0
  218. cc                second membre
  219. c                 bn(l)=0.d0                 
  220. cc           else
  221. cc                second membre
  222. c c                bn(l)=0.d0                  
  223. c            end if
  224.  
  225. c           *** deformation consolidante de Maxwell ********************
  226.             l=12+i
  227.             k=l
  228.             ann(l,k)=1.d0
  229.             bn(l)=0.d0            
  230.             if(dt.gt.0.) then
  231.                 denomm=(taum*cc3(i))**(-1)
  232.                 denomd=(taumdtt*cc3(i))**(-1)
  233.                 coeffm=denomm/(denomm+denomd)
  234. c               couplage avec depsm*
  235.                 k=9+i
  236.                 ann(l,k)=-coeffm
  237.             end if
  238.  
  239. c           *** deformation thermique transitoire **********************
  240.             l=15+i
  241. c           inconnue depsdtt            
  242.             k=l
  243.             ann(l,k)=1.d0
  244.             bn(l)=0.d0
  245.             if(dt.gt.0.) then
  246.                 coeffd=1.d0-coeffm
  247. c               couplage avec depsm*
  248.                 k=9+i
  249.                 ann(l,k)=-coeffd
  250.             end if              
  251.  
  252. c           *** contrainte effective dans la matrice *******************
  253. c           inconnue dsigp
  254.             l=18+i
  255.             k=l
  256.             ann(l,k)=1.d0
  257. c           couplage avec les def elastiques
  258.             do j=1,3
  259.                  k=3+j
  260.                  ann(l,k)=-raideur66(i,j) 
  261.             end do
  262.  
  263. c           *** PRESSIONS **********************************************
  264.  
  265. c           elles n ont pas d indice de direction donc on passe une 
  266. c           seule fois par ces lignes
  267.             if(i.eq.1) then            
  268.  
  269. c            *** pression de gel ***************************************
  270. c            inconnue dPg
  271.              l=22
  272.              k=l
  273.              ann(l,k)=1.d0
  274. c            couplage avec production gel
  275.              bn(l)=Mg*dphig
  276.              coeffg=Mg*bg
  277. c             print*,'coupl3d rag',Mg,bg,dphig,coeffg
  278. c             read*
  279. c            couplage avec les deformations
  280.              do j=1,3
  281. c                couplage avec depse            
  282.                  k=3+j                 
  283.                  ann(l,k)=coeffg
  284. c                couplage avec depsk
  285.                  k=6+j                 
  286.                  ann(l,k)=coeffg
  287. c                couplage avec depsm
  288.                  k=12+j
  289.                  ann(l,k)=coeffg
  290. c                couplage avec la depsdtt
  291.                  k=15+j
  292.                  ann(l,k)=coeffg
  293. c                couplage avec depsg
  294.                  k=37+j
  295.                  ann(l,k)=Mg
  296. c                couplage avec depss
  297.                  k=40+j                 
  298.                  ann(l,k)=coeffg
  299. c                couplage avec depst
  300.                  k=43+j                 
  301.                  ann(l,k)=0.d0
  302. c                couplage avec depsc
  303.                  k=46+j                 
  304.                  ann(l,k)=0.d0
  305.              end do 
  306.  
  307. c            *** pression de def ***************************************
  308. c            inconnue
  309.              l=23
  310.              k=l
  311.              ann(l,k)=1.d0
  312. c            couplage avec production gel
  313.              bn(l)=Md*dphid
  314. c            couplage avec les deformations
  315.              coeffd=Md*bd
  316. c             print*,'coupl3d def',Md,bd,dphid,coeffd             
  317.              do j=1,3
  318. c                couplage avec depse            
  319.                  k=3+j               
  320.                  ann(l,k)=coeffd
  321. c                couplage avec depsk
  322.                  k=6+j                 
  323.                  ann(l,k)=coeffd
  324. c                couplage avec depsm
  325.                  k=12+j
  326.                  ann(l,k)=coeffd
  327. c                couplage avec la dtt
  328.                  k=15+j
  329.                  ann(l,k)=coeffd
  330. c                couplage avec depsg
  331.                  k=37+j
  332.                  ann(l,k)=0.d0
  333. c                couplage avec depss
  334.                  k=40+j                 
  335.                  ann(l,k)=Md
  336. c                couplage avec depst
  337.                  k=43+j                 
  338.                  ann(l,k)=0.d0
  339. c                couplage avec depsc
  340.                  k=46+j                 
  341.                  ann(l,k)=0.d0
  342.              end do  
  343.  
  344. c            *** pression d'eau  ***************************************
  345. c            inconnue
  346.              l=24
  347.              k=l
  348.              ann(l,k)=1.d0
  349. c            evolution su le pas
  350.              bn(l)=Mw*dphiw
  351. c            pas de couplage avec les deformations
  352.  
  353. c            *** fin de calcul des PRESSIONS ***************************   
  354.  
  355.             end if
  356.  
  357.             if(.not.tirve) then
  358.  
  359. c            *** contraintes totales dans la matrice *******************
  360.              l=24+i
  361. c            inconnue dsigm
  362.              k=l
  363.              ann(l,k)=1.d0
  364. c            couplage avec contrainte effective
  365.              k=18+i
  366.              ann(l,k)=-1.d0
  367. c            couplage avec pression de gel
  368.              k=22
  369.              ann(l,k)=bg
  370. c            couplage avec pression de def
  371.              k=23
  372.              ann(l,k)=bd
  373. c            couplage avec la pression d eau
  374.              k=24
  375.              ann(l,k)=bw
  376.  
  377. c            *** contrainte critique pour le gel ***********************
  378.              l=27+i 
  379.              k=l
  380.              ann(l,k)=1.d0
  381. c            couplage avec la pression de gel
  382.              k=22
  383.              ann(l,k)=-Cg
  384. c            couplage avec la contrainte totale 
  385.              if(complet(i)) then
  386. c                si le critere est complet
  387.                  k=24+i            
  388.                  ann(l,k)=-1.d0
  389.              end if                           
  390.  
  391. c            *** contrainte critique pour la def ***********************
  392.              l=30+i 
  393.              k=l
  394.              ann(l,k)=1.d0
  395. c            couplage avec la pression de def
  396.              k=23
  397.              ann(l,k)=-Cd
  398. c            couplage avec la contrainte totale
  399.              if(complet(3+i)) then
  400. c                si le critere est complet              
  401.                  k=24+i
  402.                  ann(l,k)=-1.d0
  403.             end if            
  404.  
  405. c           *** contrainte critique pour l eau  ************************
  406.             l=33+i 
  407.             k=l
  408.             ann(l,k)=1.d0
  409. c           couplage avec la pression d eau
  410.             k=24
  411.             ann(l,k)=Cw
  412. c           couplage avec la contrainte totale
  413.             if(complet(6+i)) then
  414. c                si le critere est complet
  415.                  k=24+i
  416.                  ann(l,k)=-1.d0
  417.             end if 
  418.  
  419. c           *** contrainte critique de Drucker Prager ******************
  420.             if(i.eq.1) then
  421.              l=37
  422.              k=l
  423.              ann(l,k)=1.d0
  424. c            couplage avec les contraintes dans la matrice (totale)
  425.              do j=1,3
  426.                  k=24+j      
  427. c                si le critere est complet                  
  428.                  ann(l,k)=-alphadp(j)
  429.              end do
  430.             end if            
  431.  
  432. c           *** deformations plastiques pour le gel ********************
  433.             l=37+i
  434.             k=l
  435.             ann(l,k)=1.0d0
  436. c           couplage avec les multiplicateurs plastiques
  437.             k=49+i
  438.             ann(l,k)=-1.0d0
  439.  
  440.  
  441. c           *** deformations plastique pour la def *********************
  442.             l=40+i
  443.             k=l
  444.             ann(l,k)=1.0d0
  445. c           couplage avec les multiplicateurs plastiques
  446.             k=52+i
  447.             ann(l,k)=-1.0d0
  448.  
  449. c           *** deformation plastique de traction **********************
  450.             l=43+i
  451.             k=l
  452.             ann(l,k)=1.0d0
  453. c           couplage avec les multiplicateurs plastiques
  454.             k=55+i
  455.             ann(l,k)=-1.0d0
  456.  
  457. c           *** deformation plastique de compression *******************
  458.             l=46+i
  459.             k=l
  460.             ann(l,k)=1.0d0
  461. c           couplage avec le multiplicateur plastique de DP
  462.             k=59
  463.             ann(l,k)=-betadp(i)
  464.  
  465. c           ci-dessous le nbr de multiplicateurs doit être égal à nc            
  466.  
  467. c           *** multiplicateurs plastique pour le gel ******************
  468.             l=49+i
  469.             if(.not.actif(i)) then
  470. c                le multiplicateur plastique est mis a zero            
  471.                  k=l
  472.                  ann(l,k)=1.0d0
  473.                  bn(l)=0.d0
  474.             else
  475. c                la ligne est utilisee pour annuler le critere,
  476. c                le multiplicateur n intervient plus sur la ligne 
  477. c                prise en compte de la contrainte critique du gel                
  478.                  k=27+i
  479.                  ann(l,k)=1.0d0
  480. c                valeur du critere a annuler                     
  481.                  bn(l)=-factif(i)
  482. c                prise en compte de l ecrouissage
  483.                  k=37+i
  484.                  ann(l,k)=-Hg                                  
  485.             end if            
  486.  
  487. c           *** multiplicateurs plastiques pour la def *****************
  488.             l=52+i
  489.             if(.not.actif(3+i)) then
  490. c                le multiplicateur plastique est mis a zero            
  491.                  k=l
  492.                  ann(l,k)=1.0d0
  493.                  bn(l)=0.d0
  494.             else
  495. c                la ligne est utilisee pour annuler le critere
  496. c                le multiplicateur n intervient plus sur la ligne
  497. c                prise en compte de la contrainte critique de def                
  498.                  k=30+i
  499.                  ann(l,k)=1.0d0
  500. c                valeur du critere a annuler                     
  501.                  bn(l)=-factif(3+i)
  502. c                prise en compte de l ecrouissage
  503.                  k=40+i
  504.                  ann(l,k)=-Hd                                    
  505.             end if                            
  506.  
  507. c           *** multiplicateurs plastiques pour la traction ************
  508.             l=55+i
  509.             if(.not.actif(9+i)) then
  510. c               le multiplicateur plastique est mis a zero            
  511.                 k=l
  512.                 ann(l,k)=1.0d0
  513.                 bn(l)=0.d0
  514.             else
  515. c               la ligne est utilisee pour annuler le critere
  516. c               le multiplicateur n intervient plus sur la ligne            
  517.                 k=24+i
  518.                 ann(l,k)=1.0d0
  519.                 bn(l)=-factif(9+i)                
  520.             end if
  521.  
  522. c           *** multiplicateur plastique pour drucker Prager ***********
  523.             if(i.eq.1) then
  524.               l=59
  525.               if(.not.actif(13)) then
  526. c                le multiplicateur plastique est mis a zero                
  527.                  k=l
  528.                  ann(l,k)=1.0d0
  529.                  bn(l)=0.d0
  530.               else
  531. c                la ligne est utilisee pour annuler le critere
  532. c                le multiplicateur n intervient plus sur la ligne
  533.                  k=37              
  534.                  ann(l,k)=1.0d0
  535.                  bn(l)=-factif(13)
  536.               end if
  537.             end if
  538.  
  539. c           ************************************************************
  540.  
  541.           end if
  542. c         fin de l ecoulement plastique          
  543.  
  544.         end do
  545. c       fin du tir visco elastique        
  546.  
  547. c***********************************************************************
  548. c       resolution 
  549. c***********************************************************************
  550.  
  551.         if(tirve) then
  552.             ni=24
  553.         else
  554.             ni=59
  555.         end if        
  556.         call gaus3d(ni,ANN,XN,BN,ngf,errgauss,ipzero)
  557.  
  558.         goto 20
  559. c       ********* affichages eventuels *********************************        
  560.         if(.not.tirve) then
  561.           if(actif(13))then        
  562.            do i=1,ni
  563.              print*,'coupl3d :xn(',i,')=',xn(i)
  564.              print*,'coupl3d :bn(',i,')=',bn(i)
  565.              do j=1,ni
  566.                print*,'ann(',i,j,')=',ann(i,j)
  567.              end do
  568.              read*
  569.            end do
  570.            read*
  571.           end if
  572.         end if
  573.  
  574. c       ******** traitement des erreurs eventuelles ********************        
  575. 20      if(errgauss.eq.1) then
  576.              err1=1
  577.              print*,'Pb avec gaus3d dans couplagf3d'
  578.              do i=1,3
  579.               print*,'coupl3d direction:',i,'cc',cc3(i)
  580.               call jacobflu3d(tauk,psi,taum*cc3(i),taumdtt*cc3(i),
  581.      #        CM,dt,Jf,err1)     
  582.               print*,tauk,psi,taum*cc3(i),taumdtt*cc3(i),CM,dt,Jf
  583.              end do
  584.              return
  585.         end if 
  586.  
  587. c       *** reduction eventuelle de l amplitude du retour radial *******
  588.         if(.not.tirve) then
  589.              reduc=.false.
  590.              coeff=1.d0
  591.              do j=1,3
  592.                if(actif(9+j)) then
  593.                 if( (epspt6p(j)+xn(43+j)).lt.0.d0) then
  594.                    if(factif(9+j).lt.0.) then
  595. c                    on est bien dans un cas de refermeture  
  596. c                    on reduit le pas de retour radial                 
  597.                      reduc=.true.
  598.                      coeff=min(epspt6p(j)/abs(xn(43+j)),coeff)
  599.                    else
  600. c                    comme on est dans un cas d ouverture
  601. c                    on supprime ce critere et on recalcule depuis t0
  602.                      actif(9+j)=.false.
  603.                      recal=.true.
  604.                    end if                     
  605.                 end if
  606.               end if
  607.               if((coeff.le.0.d0).or.(coeff.gt.1.)) then
  608.                  err1=1
  609.                  print*,'Pb avec gaus3d dans couplagf3d'
  610.                  print*,'Coeff inconsistant dans coupl3d'
  611.                  print*,'lors des refermetures de fissures'
  612.                  print*,actif(9+j),j,epspt6p(j),xn(43+j),coeff
  613.                  return
  614.               end if
  615.             end do
  616.             if(reduc) then 
  617. c                 print*,'Ds coupl3d reduction du retour:',coeff
  618.                  do i=1,ni
  619.                      xn(i)=coeff*xn(i)
  620.                  end do
  621.             end if
  622. c           positivite des multiplicateurs plastiques de traction diffuse
  623.             do i=1,6
  624.                 if(xn(49+i).lt.0.d0) then
  625.                     actif(i)=.false.
  626.                     recal=.true.
  627.                 end if
  628.             end do
  629.         end if 
  630.         if (recal) then
  631.            if (ncal.le.nmax) then
  632. c              print*,'Reduction nbre de criteres actifs dans coupl3d'
  633. c              nombre de criteres actifs residuels
  634.                na=0               
  635.                do i=1,nc
  636.                  if(actif(i)) then
  637.                      na=na+1
  638.                  end if
  639.                end do
  640.                if(na.eq.0) then
  641.                  print*,'Pb tous les criteres ont ete deactives'
  642.                  err1=1
  643.                  return
  644.                else
  645. c                il reste des criteres a verifier              
  646.                  ncal=ncal+1
  647.                  goto 10
  648.                end if
  649.            else
  650.                print*,'Pb avec gaus3d dans couplagf3d'
  651.                print*,'Nmax atteint par annulation de criteres'
  652.                err1=1
  653.                return  
  654.            end if               
  655.         end if           
  656.  
  657. c***********************************************************************
  658. c       recuperation des increments
  659. c***********************************************************************
  660.  
  661.         dpg=xn(22)
  662.         dps=xn(23)
  663.         dpw=xn(24)
  664. c        print*,'coupl3d:',dpg,dps,dpw
  665.         do i=1,3
  666.              depse6p(i)=xn(3+i)
  667.              depse6p(3+i)=0.d0
  668.              depsk6p(i)=xn(6+i)
  669.              depsk6p(3+i)=0.d0
  670.              depsm6p(i)=xn(12+i)
  671.              depsm6p(3+i)=0.d0
  672.              depsvt6p(i)=xn(15+i)
  673.              depsvt6p(3+i)=0.d0
  674.              dsigef6p(i)=xn(18+i)
  675.              dsigef6p(3+i)=0.d0
  676.              if(.not.tirve) then   
  677.                  depspg6p(i)=xn(37+i)
  678.                  depspg6p(3+i)=0.d0  
  679.                  depsps6p(i)=xn(40+i)
  680.                  depsps6p(3+i)=0.d0 
  681.                  depspt6p(i)=xn(43+i)
  682.                  depspt6p(3+i)=0.d0                 
  683.                  depspc6p(i)=xn(46+i)
  684.                  depspc6p(3+i)=0.d0
  685.             else
  686.                  depspg6p(i)=0.d0
  687.                  depspg6p(3+i)=0.d0  
  688.                  depsps6p(i)=0.d0
  689.                  depsps6p(3+i)=0.d0    
  690.                  depspt6p(i)=0.d0
  691.                  depspt6p(3+i)=0.d0                 
  692.                  depspc6p(i)=0.d0
  693.                  depspc6p(3+i)=0.d0      
  694.             end if
  695.         end do
  696.  
  697. c       *** retour des increments de deformation dans la base fixe *****
  698.  
  699.         call chrep6(depse6p,vp33t,.false.,depse6) 
  700.         call chrep6(dsigef6p,vp33t,.false.,dsigef6) 
  701.         if(dt.gt.0.) then        
  702.             call chrep6(depsk6p,vp33t,.false.,depsk6) 
  703.             call chrep6(depsm6p,vp33t,.false.,depsm6)
  704.             call chrep6(depsvt6p,vp33t,.false.,depsvt6)
  705.         else
  706. c           en cas d increment a pas nul on ne fait pas
  707. c           le changement de base sur la visco elasticite
  708. c           car les increments sont nuls        
  709.             do i=1,6
  710.                 depsk6(i)=0.d0
  711.                 depsm6(i)=0.d0
  712.                 depsvt6(i)=0.d0
  713.             end do
  714.         end if        
  715.  
  716.         if(.not.tirve) then
  717.              call chrep6(depspg6p,vp33t,.false.,depspg6) 
  718.              call chrep6(depsps6p,vp33t,.false.,depsps6)  
  719.              call chrep6(depspt6p,vp33t,.false.,depspt6)             
  720.              call chrep6(depspc6p,vp33t,.false.,depspc6)           
  721.         else
  722. c            cas du tir VE
  723.              do i=1,6
  724.                  depspg6(i)=0.d0             
  725.                  depsps6(i)=0.d0  
  726.                  depspt6(i)=0.d0                 
  727.                  depspc6(i)=0.d0
  728.              end do
  729.         end if
  730. c       fin        
  731.         goto 100
  732.  
  733.         if(dt.eq.0.) then 
  734.          do i=1,3
  735.              dx6(i)=deps3(i)
  736.              dx6(i+3)=0.d0
  737.          end do
  738.          call chrep6(dx6,vp33t,.false.,dx16)
  739.          print*,'dans coupl3d dt=0'
  740.          print*,'mettre la convergence forcee a faux pour pasapas'
  741.          do i=1,6
  742.            print*,'deps6(',i,')=',dx16(i)
  743.          end do         
  744.          do i=1,6
  745.            print*,'depse6(',i,')=',depse6(i)
  746.          end do
  747.          do i=1,6
  748.            print*,'depsk6(',i,')=',depsk6(i)
  749.          end do
  750.          do i=1,6
  751.            print*,'depsm6(',i,')=',depsm6(i)
  752.          end do  
  753.          do i=1,6
  754.            print*,'depsvt6(',i,')=',depsvt6(i)
  755.          end do 
  756.          do i=1,6
  757.            print*,'dsigef6(',i,')=',dsigef6(i)
  758.          end do  
  759.          read*
  760.         end if         
  761. 100     return
  762.         end
  763.  
  764.  
  765.  
  766.  
  767.  

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