Télécharger la plaquette Cast3M


fr en

Télécharger kend3d.eso

Retour à la liste

Numérotation des lignes : 

kend3d

  1. C KEND3D    SOURCE    PV090527  23/01/27    21:15:46     11574          
  2.       subroutine kend3d(wpl3,vwpl33,vwpl33t,wplx3,vwplx33,vwplx33t,
  3.      # gft_app,gfr,iso,sigf6,sigf6d,err1,reft,young,nu,souplesse66,
  4.      # xmt,dtiso,dtr,dt3,dr3,dct3,dcdp,dcc3,nc3,nc33,epict,errgf1,
  5.      # errgf0,epspc6,epeqpc,ekdc,epspg6,ekdg,alphag,dgt3,dgc3,fshr6,
  6.      # dwt3,dwc3,skdw,alphaw,epsps6,dst3,dsc3,ekds,alphas,rteff_app,
  7.      # A,B,X,ipzero,ngf,tau_trac,xmc,dciso,dcpp,taudp,taudpmax,
  8.      # XE3D,NBNMAX3D,NBNB3D,IDIMB3D,altc,epssr,Method_H,Method_N,LcH)
  9.  
  10. c      calcul de l endommagement cf. :
  11. C      A. Sellier, G. Casaux-Ginestet, L. Buffo-Lacarrière, X. Bourbon,
  12. C      Orthotropic damage coupled with localized crack reclosure processing. Part I: Constitutive laws,
  13. C      Engineering Fracture Mechanics,
  14. C      Volume 97, 2013, Pages 148-167, ISSN 0013-7944,
  15. C      https://doi.org/10.1016/j.engfracmech.2012.10.012.
  16.  
  17.       implicit real*8 (a-h,o-z)
  18.       implicit integer (i-n)
  19.  
  20. c     calcul des tailles   
  21.       logical Method_H,Method_N 
  22.       real*8 LcH
  23.       integer NBNMAX3D,NBNB3D,IDIMB3D
  24.       real*8 XE3D(3,NBNMAX3D)
  25. c     etat du materiau
  26.       real*8 sigf6(6),sigf6d(6),reft,rteff_app
  27.       real*8 wpl3(3),vwpl33(3,3),vwpl33t(3,3)
  28.       real*8 wplx3(3),vwplx33(3,3),vwplx33t(3,3)
  29.       real*8 dt66p(6,6),dr66p(6,6),gft_app,long3(3)
  30.       real*8 souplesse66(6,6),epspg6(6),gfr
  31.       real*8 epspc6(6)
  32.       real*8 fshr6(6),dwt3(3),dwc3(3)
  33.       real*8 epsps6(6),dst3(3),dsc3(3)
  34.       real*8 ekdc,epeqpc,young,nu,tau_trac
  35.       real*8 xmt,epc,xmc,dcpp,dcpp1,dcdp
  36.       logical dciso
  37.       logical iso,dtiso
  38.       integer err1,i,j,k,l
  39.       real*8 ekdg,dtr,alphag,skdw,alphaw,ekds,alphas,epssr
  40.  
  41.       real*8 dt3(3),dr3(3),dct3(3),dcc3(3)    
  42.       real*8 sigf3(3),vsigf33(3,3),vsigf33t(3,3)
  43.       real*8 sigft6(6),sigfc6(6),sigfc3(3),sigft3(3)
  44.  
  45.       real*8 rt331(3,3),ref331(3,3)
  46.  
  47.       real*8 sigft6p(6),sigft6d(6),sigfc6p(6),sigfc6r(6),sigfc6ct(6)
  48.       real*8 sigft61(6),sigfc61(6)
  49.  
  50.       real*8 epspg33(3,3),epspg3(3),vepspg33(3,3),vepspg33t(3,3)
  51.       real*8 dgt3(3),dgc3(3),sigfc6d(6)
  52.  
  53.       real*8 epspt6d(6)
  54.       real*8 ekt3(3),ekr3(3)
  55.       real*8 gf1,gf2,gfrmin,gfreff
  56.  
  57.       real*8 umdtr,xx,altc
  58.       real*8 dt03(3)
  59.  
  60. c     seuil init d endo par rgi      
  61.       real*8 epsseuil0
  62.       parameter (epsseuil0=0.d0)
  63.  
  64. c     endo de refermeture de fissure      
  65.       logical endor
  66.  
  67. c     matrice des endo de traction      
  68.       real*8 umdt66(6,6),umdr66(6,6),umdct66(6,6)
  69.       integer errg
  70.  
  71. c     matrice pour methode e Gauss      
  72.       integer ngf
  73.       real*8 X(ngf),B(ngf),A(ngf,(ngf+1))
  74.       integer ipzero(ngf)
  75.  
  76. c     deformation equivalente pour endo grand lors de la reouverture 
  77. c     de fissure (lorsque la fissure se reouvre on a 2*epsmt6    
  78.       real*8 epseq3(3)
  79.  
  80. c     endo maxi et precision 
  81.       real*8 dmaxi,drgimaxi,xngf,precision1
  82.       parameter (precision1=0.d0)
  83.       parameter (dmaxi=1.d0-precision1,drgimaxi=1.0d0-precision1)
  84.       parameter (xngf=6.d0) 
  85.  
  86. c     declaration complementaires juillet 2015
  87.       real*8 wkt,sigefft,beta1,dt1,wkr,xx1,trepsdc 
  88.  
  89. c     anisotropie eventuelle de Gf en presence de multifissuration
  90.       real*8 nc33(3,3),nc331(3,3),nc3(3)
  91. c     endommagement isotrope au pic de traction
  92.       real*8 dtisopic 
  93.  
  94. c     direction dans laquelle est calculee la taille de l element
  95.       real*8 dir3(3)
  96.  
  97. c     longueur de l element dans chaque direction
  98.       real*8 longr(3)
  99.  
  100. c     energie de fissuration recuperee dans la decharge elastique
  101.       real*8 gfmin,gfpp,errgf,errgf1,errgf0,gfap 
  102. c     deformation au pic de traction
  103.       real*8 epict  
  104.  
  105. c     volume de fissure de traction
  106.       real*8 trepst,dx  
  107.  
  108. c     contrainte apres endo de traction
  109.       real*8 sigf61(6)
  110.  
  111. c     endo de compression orthotrope
  112.       real*8 epspc33(3,3),epspc3(3),vepspc33(3,3),vepspc33t(3,3)
  113.  
  114. c     correpsondance avec numero non local
  115.       integer NL1,ID1 
  116.  
  117. c     deformation fictive pour imposer un endo initial
  118.       real*8 test6(6)     
  119.  
  120. c    initialisation de l erreur  d energir pour eviter snap back      
  121.       errgf1=errgf0
  122.  
  123. c***********************************************************************
  124. c     prise en compte de l endommagement de traction
  125. c*********************************************************************** 
  126. c     endo iso max au pic de traction
  127.       if(dtiso) then
  128. c         l endommagement isotrope conduit a utiliser rteff_app=rt_app/(1-dtpic)>rt_app
  129.           dtisopic=1.d0-exp(-1.d0/xmt)
  130.       else
  131.           dtisopic=0.d0
  132.       end if     
  133.  
  134. c     passage de Gf dans la base principale des ouvertures maxi
  135.       call chre3(nc331,nc33,vwplx33) 
  136.  
  137.       do i=1,3
  138. c       energie elastique consommee au pic dans cette direction
  139.         do j=1,3
  140.             dir3(j)=vwplx33(j,i)
  141.         end do
  142. c       taille dans cette direction 
  143. C         print*, longr(i),dir3,Method_N,XE3D,NBNMAX3D,
  144. C      #  NBNB3D,IDIMB3D,Method_H,LcH,err1  
  145.         call tail1d(longr(i),dir3,Method_N,XE3D,NBNMAX3D,
  146.      #  NBNB3D,IDIMB3D,Method_H,LcH,err1)
  147. C         write(*,'(A3,I2,1x,A6,e10.3)') 'dir',i,'longr:',longr(i)
  148. C         print*, longr(i),dir3,Method_N,XE3D,NBNMAX3D,
  149. C      #  NBNB3D,IDIMB3D,Method_H,LcH,err1
  150.         if(err1.eq.1)then
  151.           print*,'Pb lors du calcul de taille1 dans kend3d'
  152.           return
  153.         end if
  154. c       print*,'------------------------------------------'
  155. c       print*,'direction',dir3,' taille',longr(i),'epict' ,epict           
  156. c       gfmin=0.5d0*(rteff_app*(1.d0-dtisopic))*epict*longr(i)
  157. c       gf mini pour controler l ecrouissage negatif post pic
  158.         gfmin=longr(i)*((rteff_app*(1.d0-dtisopic))**2)/young*
  159.      #  (4.d0/xngf+1.d0/(1.d0-dtisopic) )/2.d0 
  160. c       gf pre pic     
  161.         gfap= longr(i)*(rteff_app*(1.d0-dtisopic))*epict/2.d0  
  162. c        print*,'gfmin',gfmin
  163.         if(gfmin.lt.(nc331(i,i)*gft_app)) then
  164. c          le controle de l energie de fissuration sans snap back local est possible
  165.            gfpp=nc331(i,i)*gft_app-gfap
  166. c           print*,'kend3d gfpp cas1',gfpp
  167.         else
  168. c          un snap back local est necessaire on calcule l erreur comise
  169. c          en preservant une decroissance post pic correspondant à gft_app/xngf
  170.            gfpp=gfmin-gfap
  171. c           print*,'kend3d gfpp cas2',gfpp
  172. c          erreur max comise           
  173.            errgf=(gfmin-nc331(i,i)*gft_app)/(nc331(i,i)*gft_app)
  174. c           print*,'dans kend3d erreur gf maxi',errgf
  175.         end if                           
  176. c       ouverture caracteristique ds cette direction  
  177.         wkt=gfpp/(rteff_app*(1.d0-dtisopic))
  178. c       endo localise        
  179.         dt3(i)=wplx3(i)*(wplx3(i)+2.d0*wkt)/(wkt+wplx3(i))**2
  180. c       borne de dt a dmaxi
  181.         dt3(i)=min(dmaxi,dt3(i)) 
  182. c        print*,'kend3d dt(',i,')=',dt3(i)  
  183. c       erreur relative effective        
  184.         errgf1=max(errgf1*dt3(i),errgf1)
  185.       end do
  186. c     calcul de la matrice d endommagement de traction
  187. c      print*,'ngf dans endo3d',ngf
  188. c      print*,'endo3d umdt'
  189.       call umdt3d(young,nu,souplesse66,dt3,umdt66,
  190.      # a,b,x,ipzero,ngf,errg,iso) 
  191.  
  192.  
  193. c***********************************************************************
  194. c     appli de l endo de traction aux contraintes positives
  195. c***********************************************************************
  196.  
  197. c     partition du tenseur des contraintes
  198.       call prtt3d(sigf6,sigf3,vsigf33,vsigf33t,
  199.      #sigft6,sigfc6,sigfc3,sigft3) 
  200. c     ******************************************************************
  201. c     traitement de l endo isotrope de traction pre pic si necessaire
  202. c     contrainte effective maxi      
  203.       sigefft=max(sigft3(1),sigft3(2),sigft3(3)) 
  204. c     partie elastique de l ouverture de fissure
  205.       tau_trac=sigefft/rteff_app
  206.       if(dtiso) then
  207. c       a positionner ici si utile        
  208.         beta1=(tau_trac**xmt)/xmt
  209. c       endommagement        
  210.         dt1=1.d0-dexp(-beta1)
  211. c       borne de dtiso
  212.         dt1=min(dmaxi,dt1)       
  213. c       actualisation de l endommagement iso de traction
  214.         dtr=max(dt1,dtr) 
  215. c        print*,'kend3d:',sigefft,rteff_app,xmt,beta1,dtr        
  216.       else
  217.         dtr=0.d0
  218.       end if
  219.       umdtr=1.d0-dtr
  220. c     ******************************************************************             
  221. c     passage  des contraintes effectives dans la base prin des endo
  222.       call chrep6(sigft6,vwplx33,.false.,sigft6p)  
  223. c     application du tenseur d endommagement aux contraintes de tractions
  224.       do i=1,6
  225.         sigft6d(i)=0.d0
  226.         do j=1,6
  227.              sigft6d(i)=sigft6d(i)+umdt66(i,j)*sigft6p(j)
  228.         end do
  229. c       application endo trac iso pre pic le cas echeant        
  230.         if(dtiso) then
  231. c         application de l endo iso prepic
  232.           sigft6d(i)=sigft6d(i)*umdtr
  233.         end if          
  234.       end do
  235. c     retour des contraintes positives en base fixe
  236.       call chrep6(sigft6d,vwplx33t,.false.,sigft61)  
  237.  
  238. c***********************************************************************
  239. c    reconstruction du tenseur des contraintes
  240. c***********************************************************************
  241.       do i=1,6
  242.          sigf61(i)=sigfc6(i)+sigft61(i)
  243.       end do
  244. c     Nouvelle partition du tenseur des contraintes endommagés par les Dt
  245.       call prtt3d(sigf61,sigf3,vsigf33,vsigf33t,
  246.      #sigft61,sigfc6,sigfc3,sigft3) 
  247.  
  248. c***********************************************************************
  249. c     calcul de la fonction de refermeture de fissure
  250. c***********************************************************************
  251.  
  252. c     materiaux et geometrie
  253.       do i=1,3
  254. c       resistance a la refermeture anisotrope
  255.         do j=1,3
  256.             if(i.eq.j) then
  257.                 ref331(i,i)=reft
  258.             else
  259.                 ref331(i,j)=0.d0
  260.             end if
  261.         end do
  262.       end do
  263. c     calcul des indicateur de refermetures 
  264. c     tailles des elements dans les directions de refermeture      
  265.       call tail3d(long3,vwpl33,Method_N,Method_H,
  266.      # XE3D,NBNMAX3D,NBNB3D,IDIMB3D,LcH,err1)
  267.       if(err1.eq.1) then
  268.          print*,'Pb lors du calcul de tail3d dans endo3d'
  269.          return
  270.       end if
  271. c     nombre  de fissures dans cette direction obtenu
  272. c     par intrpoation entre les nombres de fissure localisees
  273.       call chre3(nc331,nc33,vwpl33)      
  274.       do i=1,3
  275.           gfrmin=2.d0*(ref331(i,i)**2/young)*long3(i)
  276.           if(gfr*nc331(i,i).lt.gfrmin) then
  277. c             on adopte le gfr mini necessaire
  278.               gfreff=gfrmin
  279. c             mais on conserve le nombre de fissures             
  280.               nc3(i)=nc331(i,i)
  281.           else
  282.               gfreff=gfr*nc331(i,i)
  283. c             nombre de fissure dans l element dans cette direction              
  284.               nc3(i)=nc331(i,i)
  285.           end if          
  286.           wkr=gfreff/ref331(i,i)     
  287. c         nouvelles valeurs          
  288.           dr3(i)=wpl3(i)*(wpl3(i)+2.d0*wkr)/(wkr+wpl3(i))**2
  289. c         borne de dr3
  290.           dr3(i)=max(min(dmaxi,dr3(i)),0.d0)      
  291. c          print*,'kend3d wpl3(',i,')=',wpl3(i),'dr3(',i,')=',dr3(i)
  292.       end do
  293. c      print*,'endo3d umdr'
  294.       call umdt3d(young,nu,souplesse66,dr3,umdr66,
  295.      # a,b,x,ipzero,ngf,errg,iso)
  296.  
  297. c***********************************************************************
  298. c     application de la fonction de refermetures aux contraintes negatives
  299. c***********************************************************************
  300.  
  301. c     passage  des contraintes effectives dans la base prin des endo actuels
  302.       call chrep6(sigfc6,vwpl33,.false.,sigfc6p)  
  303. c     application du tenseur d endommagement aux contraintes de tractions
  304.       do i=1,6
  305.          sigfc6r(i)=0.d0
  306.          do j=1,6
  307.             sigfc6r(i)=sigfc6r(i)+umdr66(i,j)*sigfc6p(j)
  308.          end do            
  309.       end do
  310. c     retour des contraintes negatives en base fixe
  311.       call chrep6(sigfc6r,vwpl33t,.false.,sigfc61)
  312.  
  313. c***********************************************************************
  314. c     endommagement de couplage traction directe sur compression
  315. c***********************************************************************
  316.       call buck3d(altc,dr3,dct3,vwpl33,vwpl33t,sigfc61)
  317.  
  318. c***********************************************************************
  319. c     prise en compte de l endommagement du a la pression capillaire
  320.       call endort(fshr6,dwt3,dwc3,dmaxi,drgimaxi,skdw,
  321.      #alphaw,0.d0,sigft61,sigfc61)
  322. c***********************************************************************      
  323.  
  324. c***********************************************************************
  325. c     prise en compte de l endommagement du  a la rag
  326. c      do i=1,3
  327. c         test6(i)=0.003d0
  328. c         test6(3+i)=0.d0
  329. c      end do
  330. c      call endort(test6,dgt3,dgc3,dmaxi,drgimaxi,0.003d0,
  331. c     # 0.15d0,0.d0,sigft61,sigfc61)      
  332.       call endort(epspg6,dgt3,dgc3,dmaxi,drgimaxi,ekdg,
  333.      #alphag,epssr,sigft61,sigfc61)
  334. c*********************************************************************** 
  335.  
  336. c***********************************************************************
  337. c     prise en compte de l endommagement du  a la def
  338.       call endort(epsps6,dst3,dsc3,dmaxi,drgimaxi,ekds,
  339.      #alphas,epsseuil0,sigft61,sigfc61)
  340. c***********************************************************************     
  341.  
  342. c***********************************************************************
  343. c     variation de raideur due aux dilatations transverses
  344. c     induites par le cisaillement 
  345. c***********************************************************************
  346.       call damdp3d(precision1,ekdc,0.d0,epspc6,dcc3,sigfc61,sigft61,
  347.      #dcdp,.true.)     
  348.  
  349. c***********************************************************************
  350. c    endommagement pre pic par cisaillement applicable en compression
  351. c***********************************************************************
  352. 100   if(dciso) then
  353.           xcpp=max(taudp,0.d0)/taudpmax
  354.           xcpp=(xcpp**xmc)/xmc
  355.           dcpp1=1.d0-exp(-xcpp)
  356.           dcpp=max(dcpp,dcpp1)
  357.           do i=1,6
  358.              sigf6d(i)=sigfc61(i)*(1.d0-dcpp)+sigft61(i)
  359.           end do          
  360.       else
  361.           dcpp=0.d0
  362.           do i=1,6
  363.             sigf6d(i)=sigfc61(i)+sigft61(i)
  364. c           print*,'endo3d',sigf6(i),sigf6d(i)
  365.           end do 
  366.       end if
  367. c      print*,'ds endo3d perte de rigidite isotrope comp :',dcpp 
  368. c      print*,xmc,taudp,taudpmax,dcpp1 
  369.  
  370. c***********************************************************************   
  371. c     traitement des erreurs et sortie
  372. c***********************************************************************
  373.       goto 1000 
  374.       print*,'dans endo3d long3:',long3      
  375. 999   do i=1,6
  376.          sigf6d(i)=sigf6(i)
  377.       end do
  378.       print*,'dcpp',dcpp
  379.       print*,'dcdp',dcdp
  380. 998   do i=1,3
  381.         print*,'dt3(',i,')=',dt3(i)
  382.         print*,'dr3(',i,')=',dr3(i)
  383.         print*,'dwt3(',i,')=',dwt3(i)
  384.         print*,'dwc3(',i,')=',dwc3(i)
  385.         print*,'dgt3(',i,')=',dgt3(i)
  386.         print*,'dgc3(',i,')=',dgc3(i)
  387.         print*,'dst3(',i,')=',dst3(i)
  388.         print*,'dsc3(',i,')=',dsc3(i)
  389.         print*,'dcc3(',i,')=',dcc3(i)
  390.         print*,'dct3(',i,')=',dct3(i)
  391.       end do
  392.       read*
  393.  
  394. 1000  return
  395.       end      
  396.  
  397.  

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