Esope source files

   1 : C CMP2DE    SOURCE    PV        07/11/23    21:15:57     5978           
   2 : CCC
   3 : C **********************************************************************
   4 : CCC
   5 :       SUBROUTINE CARAC_MATE_POWDER2_DENSI(XMAT,XDENSI)
   6 :       IMPLICIT INTEGER(I-N)
   7 :       implicit real*8(a-h,o-z)
   8 :       real*8    xmat(*),xdensi,eta,nnc,epref,aux1,aux2,aux3
   9 :       real*8                  kap,xmu,xmu1,alp1,xmu2,alp2,xmu3,alp3
  10 :       logical                 ogden
  11 :       common /elasdata_ogden/ kap,xmu,xmu1,alp1,xmu2,alp2,xmu3,alp3,
  12 :      .                        ogden
  13 :       real*8               nn1,nn2,eta0,aa1,aa2,bb1,bb2,
  14 :      .                     pia,pib,sigy0,pcc,rrr
  15 :       common /poder2_data/ nn1,nn2,eta0,aa1,aa2,bb1,bb2,
  16 :      .                     pia,pib,sigy0,pcc,rrr
  17 :       real*8            you,xnu
  18 :       common /elasdata/ you,xnu
  19 :       real*8    cohmax,phimax,phimin,phi
  20 : 
  21 :       ogden=.false.
  22 :       you    = xmat( 1)
  23 :       xnu    = xmat( 2)
  24 :       cohmax = xmat( 5)
  25 :       phimax = xmat( 6)
  26 :       eta0   = xmat( 7)
  27 :       phimin = xmat( 8)
  28 :       nnc    = xmat( 9)
  29 :       sigy0  = xmat(10)
  30 :       nn1    = xmat(11)
  31 :       nn2    = xmat(12)
  32 :       eta    = xdensi*eta0
  33 :       aux1   = (phimin - phimax)/(1+eta0**2-2*eta0)
  34 :       aux2   = -2.D0 * aux1
  35 :       aux3   = phimax + aux1
  36 :       phi    = aux1 * eta**2 + aux2 * eta + aux3
  37 :       if (eta.lt.eta0) phi = phimin
  38 :       if (eta.gt.1.D0) phi = phimax
  39 : 
  40 :       aa1    = (1.D0 - eta**2)/(2.D0 + eta**2)
  41 :       if (aa1.lt.1.D-30) aa1=1.D-30
  42 :       aa1    = aa1**nn1
  43 :       if (aa1.lt.1.D-5) aa1=1.D-5
  44 :       aa2    = (eta - 0.98D0*eta0)/(1.D0 - 0.98D0*eta0)
  45 : 
  46 :       if (eta.lt.eta0) then 
  47 :          aa2min  = 0.02D0*eta0/(1.D0 - 0.98D0*eta0)
  48 :          daa2min = 1.D0 /(1.D0 - 0.98D0*eta0)
  49 :          aa2     = aa2min
  50 :       endif
  51 : 
  52 :       coh    = cohmax * aa2**nnc
  53 :       aa2    = aa2**nn2
  54 : 
  55 : c cohesion .le. extremo de la elipse a cortante
  56 :       if (coh.gt.(sqrt(aa2)*sigy0)) then 
  57 :          write(*,*)' Cambio de COH',coh,sqrt(aa2)*sigy0
  58 :          coh = sqrt(aa2)*sigy0
  59 :       endif
  60 : 
  61 : ****** D-P Meschke examples
  62 : * Cauchy
  63 : *      coh = cohmax / xdensi
  64 : * Kirk
  65 : *      coh = cohmax 
  66 : *      phi = phimax
  67 : ****** D-P Meschke examples
  68 : 
  69 :       bb1    = sqrt(2.D0/3.D0)*dtan(phi)
  70 :       bb2    = sqrt(2.D0/3.D0)*coh
  71 :       pib    = -sqrt(6.D0*aa2/aa1)*sigy0
  72 :       pia    = (3.D0*bb2-sqrt(6.D0*aa2)*sigy0)/bb1
  73 :       pcc    = bb2*3.D0/bb1
  74 :       rrr    = abs(pia-pib)/(sqrt(aa2/3.D0)*sigy0)
  75 : 
  76 :       return
  77 :       end
  78 :  
  79 :  
  80 :  
  81 :