Télécharger la plaquette Cast3M


fr en

Télécharger mater2tg.eso

Retour à la liste

Numérotation des lignes : 

mater2tg

  1. C MATER2TG  SOURCE    CB215821  16/04/21    21:17:44     8920
  2.       SUBROUTINE mater2tg(cmat,ntot,lerror,idimpara3,
  3.      *                  parahot3,npttot3,H66t,sigma6v,deps6v,lcp,
  4.      *                  NOEL,NPT,Kinc)
  5.  
  6. c cmat is the name (characters) of the material
  7. c ntot is the integration point (global number for the whole structure): counter
  8. c idimpara3 is the number of parameters for parahot3
  9. c parahot3 is the array that contains the material parameters and state variables
  10. c npttot3 is the number of integration points (number for the considered material 3D)
  11. c         and is thus the number of columns in parahot3
  12. c H66t is the tangent matrix
  13. c sigma6v is the nominal stress vector
  14. c deps6v is the incremental strains vector
  15.  
  16. c===========================================================================
  17.  
  18. c input
  19. c     CMAT     : name of the material
  20. c     NTOT     : # of the integration point ( global # for the entire structure )
  21. c                 can be NTOT1 or NTOT2 or NTOT3
  22.  
  23. c output
  24. c     H66t     : tangent stiffness matrix for 3D materials (THG add)
  25.  
  26. c input and output
  27. c     PARAHOT3 : the first columns contain the material properties at elevated temperature,
  28. c                the last columns contain strains, stresses, ....
  29.  
  30. c     SIGMA6v  : stress vector for 3D materials (THG add)
  31. c     DEPS6V   : incremental strain vector for 3D materials (THG add)
  32.  
  33. c=======================================================================
  34.  
  35.  
  36. c     This routine is run for every iteration,
  37. c                         for every point of integration.
  38. c
  39. C     At the integration point  NTOT  is the material described by CMAT.
  40. c     This material has, at the point NTOT, the caracteristics PARAHOT
  41. c           at elevated temperature.
  42.  
  43. c     In multiaxial situation
  44. c       We are going to calculate the stress vector sigm6v and the stiffness
  45. c       matrix H66t as a function of the strain incr. vector deps6v,
  46. c       calculated from the strain vector at the beginning of the time step.
  47.  
  48. c i o
  49. c     a       : df/dsigma
  50. c               df/dsigmax, df/dsigmay, df/dtauxy, df/dsigmaz
  51. c   * Ai      : Hardening parameter
  52.  
  53. c THG add concrete model 3D
  54.  
  55.       IMPLICIT REAL*8 (A-B,D-H,O-Z)
  56.       implicit integer (I-K,M,N)
  57.       implicit logical (L)
  58.       implicit character*10 (C)
  59.  
  60.       dimension parahot3(idimpara3,npttot3)
  61.       dimension epstr6v(6), esigmac(6)
  62.       dimension esigma6v(6),H66t(6,6),deps6v(6),sigma6v(6)
  63.       dimension rt1(1,6),rc1(1,6),H66tsave(6,6),H66(6,6)
  64.  
  65.       i1 = 1
  66.       i2 = 2
  67.       i3 = 3
  68.       i4 = 4
  69.       i5 = 5
  70.       i6 = 6
  71.       i7 = 7
  72.       i8 = 8
  73.       i9 = 9
  74.       i10=10
  75.       i11=11
  76.       i12=12
  77.       i13=13
  78.       i14=14
  79.       i15=15
  80.       i16=16
  81.       i17=17
  82.       i18=18
  83.       i19=19
  84.       i20=20
  85.       i21=21
  86.       i22=22
  87.       i23=23
  88.       i24=24
  89.       i25=25
  90.       i26=26
  91.       i27=27
  92.       i28=28
  93.       i29=29
  94.       i30=30
  95.       i31=31
  96.       i32=32
  97.       i33=33
  98.       i34=34
  99.       i35=35
  100.       i36=36
  101.       i37=37
  102.       i38=38
  103.       i39=39
  104.  
  105.       r2 = 2.
  106.       r3 = 3.
  107.  
  108. c     TEMPERATURE
  109.       T = parahot3(idimpara3-39,ntot)
  110.       TMAX = parahot3(idimpara3-38,ntot)
  111.  
  112.  
  113. c     Compressive strength
  114.       fc20 = parahot3(i3,ntot)
  115.       rloc = rKfcsi(t,tmax)
  116.       parahot3(i3,ntot) = rloc * fc20
  117.       fc = parahot3(i3,ntot)
  118.  
  119. c     Tensile strength
  120.       ft20 = parahot3(i4,ntot)
  121.       rloc = rkftco(t,tmax)
  122.       parahot3(i4,ntot) = rloc * ft20
  123.       if (parahot3(i4,ntot).lt.(1.d5)) then
  124. c       minimal value for ft = 0.10 MPa implemented for numerical stability
  125. c       so if ft entered by the user < 0.10 MPa, the concrete is considered as
  126. c       previously cracked with dt>0.96 and ft = 1.00 MPa
  127.         parahot3(i4,ntot) = (1.d6)
  128.         parahot3(idimpara3-41,ntot) =
  129.      .             max(0.96d0,parahot3(idimpara3-41,ntot))
  130.       endif
  131.       ft = parahot3(i4,ntot)
  132.  
  133. c     Strain at peak stress
  134. c     Strain at peak stress according to EC2: epsc1,ec2
  135.       rloc = 1.0d0
  136.       epsc1ec2 = EPSc1co(tmax,rloc)
  137. c     Strain at peak stress according to ENV (minimal value): epsc1,env
  138.       epsc1env = EPSc1co(tmax,-rloc)
  139. c     Strain at peak stress for the instantaneous stress-strain relationship
  140. c     epsc1,ETC = (2 * epsc1,min + epsc1,EC2)/3
  141.       epsc1etc = (2.*epsc1env+epsc1ec2)/3.
  142.       parahot3(i5,ntot) = epsc1etc
  143.  
  144.       epsc1ec220 = EPSc1co(20.d0,rloc)
  145.       epsc1env20 = EPSc1co(20.d0,-rloc)
  146.       epsc1etc20 = (2.*epsc1env20+epsc1ec220)/3.
  147.  
  148. c     Modulus at the origin
  149.       parahot3(i1,ntot) = 2.0d0*parahot3(i3,ntot)/parahot3(i5,ntot)
  150.  
  151. c     Compressive limit of elasticity
  152.       parahot3(i6,ntot) = (0.3d0)*parahot3(i3,ntot)
  153.       fc0 = parahot3(i6,ntot)
  154.  
  155. c     Biaxial compression strength: considered as 1.20 x compressive strength
  156.       parahot3(7,ntot) = (1.2d0)*parahot3(3,ntot)
  157.  
  158. c     thermal strain
  159. c     first, save thermal strain at the precedent step
  160.       parahot3(idimpara3-31,ntot) = parahot3(idimpara3,ntot)
  161. c     then, calculate new thermal strain
  162.       parahot3(idimpara3,ntot) = epsthsi(t,tmax)
  163.  
  164. c     Modulus used in tension
  165.       parahot3(i12,ntot) = parahot3(i1,ntot)
  166.  
  167. c     Material parameter dc: damage in compression at peak stress
  168.       dco = parahot3(i9,ntot)
  169.  
  170. c     Material parameter x,c: adimensional densitive crack energy in compression
  171. c     This parameter is assumed to be constant with temperature
  172.       rxc = parahot3(i10,ntot)
  173.  
  174. c     Model parameter kappa,c1: accumulated compressive plastic strain at peak stress
  175.       denom = (2.-2.*dco)*((parahot3(i8,ntot))-1.)
  176.       parahot3(i13,ntot) = -parahot3(i5,ntot)*(1.-2.*dco)/(denom)
  177.       rkappa1 = parahot3(i13,ntot)
  178.  
  179. c     Model Parameter ac
  180.       parahot3(i14,ntot) = ( -LOG(1-dco)/(rkappa1) )
  181.  
  182. c     Model Parameter bc
  183.       parahot3(i15,ntot) = 2*fc*rxc/((1-rxc)*((fc0*rkappa1)+
  184.      .                     ((fc-fc0)*rkappa1*LOG(r2))))
  185.  
  186. c     Parameter at
  187.       if (ft20.ge.1.d6) then
  188.         parahot3(i16,ntot) = (7.0d0*ft20)/(12.0d0*
  189.      .                    parahot3(i11,ntot))
  190.       else
  191. c     for numerical reason at cannot be equal to 0 even if ft0 = 0
  192.         parahot3(i16,ntot) = (7.0d0*1.0d6)/(12.0d0*
  193.      .                    parahot3(i11,ntot))
  194.       endif
  195.  
  196. c     TRANSIENT CREEP STRAIN
  197.       do iloc=1,6
  198.         epstr6v(iloc) = 0.d0
  199.       enddo
  200. c     Phi function
  201. c     First, phi(ti) is saved as phiprev
  202.       phiprev = parahot3(17,ntot)
  203. c     Then, phi(ti+1) is calculated and saved as phi
  204. c     phi = 2*(epsc1,EC2 - epsc1,min) / (3*(fc/fck))
  205.       phi = 2.*(epsc1ec2-epsc1env)/(3.*(fc/fc20))
  206.       parahot3(17,ntot) = phi
  207. c     Calculation of transient creep strain
  208. c     If the temperature is bellow Tmax (T is decreasing or it is
  209. c     not a first heating) the transient creep strain is unchanged
  210.       if (T.lt.(tmax-1.0d-3)) then
  211.         continue
  212. c     Else, it is a first-time heating
  213. c     so there is an increment of transient creep strain under compressive stress
  214.       else
  215. c       unless the concrete is in the descending branch of the material law
  216. c       if (kappa,c < kappa,c1) (i.e. material in the ascending branch of the
  217. c                                  compressive part of the material law)
  218.         rkappac = parahot3(idimpara3-36,ntot)
  219.         if (rkappac.le.rkappa1) then
  220. c         epstr,i+1 = epstr,i + [phi(Ti+1) - phi(Ti)] * (1 - dc,i) * (sigma,eff-,i / fck)
  221.           dc = parahot3(idimpara3-40,ntot)
  222.           do iloc=i1,i6
  223.             esigmac(iloc) = parahot3(29+iloc,ntot)
  224.             epstr6v(iloc) = parahot3(17+iloc,ntot) + ((phi-phiprev)/
  225.      .                      fc20)*(1.0d0 - dc) * esigmac(iloc)
  226.           end do
  227. c         epstr cannot decrease in absolute value
  228.           do iloc=i1,i6
  229.             if (epstr6v(iloc).lt.parahot3(17+iloc,ntot)) then
  230.               parahot3(17+iloc,ntot) = epstr6v(iloc)
  231.             endif
  232.           end do
  233.         endif
  234.       endif
  235.  
  236. c     CALCULATE THE INCREMENT IN MECHANICAL STRAIN DEPS
  237. c     (the increment in thermal strain is subtracted)
  238. c     Deps = Deps - (Epsth(n)-Epsth(n-1))
  239.       deps6v(1) = deps6v(1) - parahot3(idimpara3,ntot)
  240.      .                      + parahot3(idimpara3-31,ntot)
  241.       deps6v(2) = deps6v(2) - parahot3(idimpara3,ntot)
  242.      .                      + parahot3(idimpara3-31,ntot)
  243.       if (.not.lcp) deps6v(3) = deps6v(3) - parahot3(idimpara3,ntot)
  244.      .                      + parahot3(idimpara3-31,ntot)
  245.  
  246.       lg = .false.
  247.       lbroyden = .false.
  248.       lfulldamage = .false.
  249.  
  250.       n = 1
  251.    23 continue
  252.       do istep = 1,n
  253. c     To improve convergence, the step can be subdivided
  254.  
  255.  
  256. c       ***********************************************************
  257. c           *********************
  258. c           ! 3D concrete model !
  259. c           !    PLASTIC BOX    !
  260. c           *********************
  261.  
  262.         call DPRAN3D(parahot3(i1,ntot),idimpara3,H66t,rt1,rc1,
  263. c            *******
  264.      .             esigma6v,deps6v,lg,ntot,lerror,i1,i2,i3,i4,i5,
  265.      .             i6,n,lnoconv,istep,lcp,lbroyden,H66,lfulldamage)
  266. c       Plastic Box of the 3D plastic-damage concrete model by T. Gernay
  267. c       Yield functions of Drucker Prager in compression and Rankine in tension
  268. c       Calculate the effective stresses and the effective elastoplastic tangent
  269. c       modulus Dt (implicit process)
  270. c       ***********************************************************
  271.  
  272.         if (lnoconv) then
  273.           n = n*2
  274. c         the step is subdivided by factor 2
  275. c         we come back to the state at the previous (converged) step: rkappac = rkappaic
  276.           parahot3(idimpara3-35,ntot) = parahot3(idimpara3-37,ntot)
  277.           parahot3(idimpara3-34,ntot) = parahot3(idimpara3-36,ntot)
  278.           if (n.lt.32) go to 23
  279. c         the step can be divided in maximum 16 "sub-steps"
  280.           if (.not.lbroyden) then
  281. c           if still no convergence when dividing the size of the step => try another algorithm
  282.             lbroyden = .true.
  283.             n = 1
  284.             go to 23
  285.           endif
  286.           if (.not.lfulldamage) then
  287. c           still no convergence => if completely damaged, "remove" the point (sigma = 0)
  288.             dc = parahot3(idimpara3-40,ntot)
  289.             dt = parahot3(idimpara3-41,ntot)
  290.             if ((dt.gt.0.97).or.(dc.gt.0.97)) then
  291.               lfulldamage = .true.
  292.               n = 1
  293.               go to 23
  294.             endif
  295.           endif
  296. c         if we arrive here, all the strategies to improve convergence have been tested with
  297. c         no success => lerror is true (send back an error from the MATERIAL subroutines)
  298.           lerror = .true.
  299.           return
  300.         endif
  301.       enddo
  302.  
  303. c     if ((NOEL.eq.18).and.(NPT.eq.23)) then
  304. *       print *,'SORTIE DPRAN3D contraintes = ',esigma6v
  305. c     endif
  306.  
  307. c     *************************************************************
  308. c         *********************
  309. c         ! 3D concrete model !
  310. c         !     DAMAGE BOX    !
  311. c         *********************
  312.       call DAMTG3D(parahot3(i1,ntot),idimpara3,H66t,rt1,rc1,
  313. c          *******
  314.      .    esigma6v,lg,ntot,lerror,i1,i2,i3,i4,i5,i6,H66,
  315.      .    lfulldamage,lcp)
  316. c     Damage Box of the 3D plastic-damage concrete model by T. Gernay
  317. c     Calculate the nominal stresses and the nominal consistent algorithmic
  318. c     tangent modulus Ct (explicit process)
  319. c     *************************************************************
  320.  
  321. c     if ((NOEL.eq.18).and.(NPT.eq.23)) then
  322. *       print *,'SORTIE DAMTG3D contraintes = ', sigma6v
  323. c     endif
  324.  
  325.  
  326. c     The tangent matrix is symmetrized
  327. c     because SAFIR uses a symmetric matrix for the resolution of the equilibrium
  328. c                T
  329. c     B = ( A + A ) / 2
  330.       do iloc=2,6
  331.         rloc = H66t(1,iloc) + H66t(iloc,1)
  332.         H66t(1,iloc) = rloc/(2.d0)
  333.         H66t(iloc,1) = H66t(1,iloc)
  334.       enddo
  335.       do iloc=3,6
  336.         rloc = H66t(2,iloc) + H66t(iloc,2)
  337.         H66t(2,iloc) = rloc/(2.d0)
  338.         H66t(iloc,2) = H66t(2,iloc)
  339.       enddo
  340.       do iloc=4,6
  341.         rloc = H66t(3,iloc) + H66t(iloc,3)
  342.         H66t(3,iloc) = rloc/(2.d0)
  343.         H66t(iloc,3) = H66t(3,iloc)
  344.       enddo
  345.       do iloc=5,6
  346.         rloc = H66t(4,iloc) + H66t(iloc,4)
  347.         H66t(4,iloc) = rloc/(2.d0)
  348.         H66t(iloc,4) = H66t(4,iloc)
  349.       enddo
  350.       rloc = H66t(5,6) + H66t(6,5)
  351.       H66t(5,6) = rloc/(2.d0)
  352.       H66t(6,5) = H66t(5,6)
  353.  
  354. c     Save the nominal stresses
  355.       do iloc=i1,i6
  356.         sigma6v(iloc) = PARAHOT3(idimpara3-13+iloc,ntot)
  357.       end do
  358.  
  359. c      if (lcp) then
  360. cc     Plane stress model
  361. c        parahot2(idimpara2-i3,ntot) = sigma6v(1)
  362. c        parahot2(idimpara2-i2,ntot) = sigma6v(2)
  363. c        parahot2(idimpara2-i1,ntot) = sigma6v(4)
  364. c        PARAHOT2(idimpara2-i7,ntot)=parahot3(idimpara3-18,ntot)
  365. c        PARAHOT2(idimpara2-i6,ntot)=parahot3(idimpara3-17,ntot)
  366. c        PARAHOT2(idimpara2-i5,ntot)=parahot3(idimpara3-15,ntot)
  367. c      endif
  368.       END
  369.  
  370.  
  371.  
  372.  
  373.  

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