Télécharger la plaquette Cast3M


fr en

Télécharger fefp2.eso

Retour à la liste

Numérotation des lignes : 

fefp2

  1. C FEFP2     SOURCE    OF166741  25/02/21    21:16:27     12166          
  2.  
  3. *******************************************************************
  4.       SUBROUTINE FEFP2(MATE,INPLAS,MELE,MELEME,MINTE,xMATRI,
  5.      .     NBELEM,NBPTEL,NBNN,LRE,MFR,
  6.      .     IVADESP,IVADPI,IVARI,IVAMAT,
  7.      .     IVASTF,IVARIF,IVADPF,LHOOK,IRIGE7,
  8.      .     NDEP,NDEF,NSTRS,NVARI,NMATT,PRECIS,NITMAX,NUPDATE,
  9.      .     KERRE)
  10. *******************************************************************
  11. *
  12. * entrees :
  13. *  mate   = numero de material elastico
  14. *  inplas = numero de material inelastico
  15. *  mele   = numero de elemento finito
  16. *  meleme = puntero a la malla
  17. *  minte  = puntero valores integracion
  18. *  xmatri = puntero de matrices elementares
  19. *  nbelem = numero de elementos
  20. *  nbptel = nombre de puntos por elemento
  21. *  nbnn   =
  22. *  lre    =
  23. *  mfr    =
  24. *  ivadesp= puntero a un segmento mptval de desplazamientos
  25. *  ivadpi = puntero a un segmento mptval de deformaciones inelasticas
  26. *  ivari  = puntero a un segmento mptval de variables internas
  27. *  ivamat = puntero a un segmento mptval de material
  28. *  lhook  = dimensiones de la matriz de hook
  29. *  irige7 = simetria de la matriz tangente
  30. *  ndep   = numero de componentes de desplazamientos
  31. *  ndef   = numero de componentes de deformaciones inelasticas
  32. *  nstrs  = numero de componentes de tensiones
  33. *  nvari  = numero de componentes de variables internas
  34. *  nmatt  = numero de componentes de propiedades materiales
  35. *  precis = precision de las iteraciones internas
  36. *  nitmax = numero maximo de iteraciones
  37. *  nupdate = total (0) / update (1) lagrangian formulation
  38. *
  39. * sorties :
  40. *  ivastf = puntero a un segmento mptval de tensiones
  41. *  ivarif = puntero a un segmento mptval de variables internas
  42. *  ivadpf = puntero a un segmento mptval de deformaciones inelasticas
  43. *  kerre  = indicador de error
  44. *
  45. ************************************************************************
  46.  
  47.       IMPLICIT INTEGER(I-N)
  48.       IMPLICIT REAL*8(A-H,O-Z)
  49.  
  50. -INC PPARAM
  51. -INC CCOPTIO
  52.  
  53. -INC SMCHAML
  54. -INC SMELEME
  55. -INC SMCOORD
  56. -INC SMMODEL
  57. -INC SMINTE
  58. -INC SMRIGID
  59.  
  60. -INC TMPTVAL
  61.  
  62. cc    hay +4 en XMAT para compatibilidad con feap
  63.       SEGMENT WRK0
  64.         REAL*8 XMAT(NMATT+4)
  65.         REAL*8 DDHOOK(LHOOK,LHOOK),DEFPINI(NDEF),DEFPFIN(NDEF)
  66.         REAL*8 SIGF(NSTRS),VAR0(NVARI),VARF(NVARI),BE(NDEF)
  67.         REAL*8 REL(LRE,LRE),REL2(LRE,LRE),SHPWRK(6,NBNN)
  68.         REAL*8 BGENE(NSTRS,LRE)
  69.         REAL*8 XENEW(3,NBNN),XEOLD(3,NBNN),XMOVI(3,NBNN)
  70.         REAL*8 FFI(3,3),FF(3,3),FFIv(9),FFv(9)
  71.       ENDSEGMENT
  72.  
  73.       SEGINI WRK0
  74.  
  75. ******************
  76.       nescri   =0
  77.       nues     =199
  78. c numerical differentiation parameters
  79.       nnumer   =0
  80.       deltax   =0.D0
  81. c line-search parameters
  82.       level    =1
  83.       kmax     =3
  84.       iaugla   =0
  85.       caugla   =0.D0
  86. c number of internal variable 'relative dendity'
  87.       nvardet  =4
  88. c model number
  89.       if      (inplas.eq.114) then
  90.          nmodel   =8
  91.       else if (inplas.eq.115) then
  92.          nmodel   =2
  93.       else if (inplas.eq.116) then
  94.          nmodel   =5
  95. c        DO NOT USE LINE-SEARCH
  96.          kmax     =0
  97.       else if (inplas.eq.117) then
  98.          nmodel   =6
  99.       else
  100.         write(nues,*)'  Model not defined',inplas
  101.       endif
  102. c*****************************************************************
  103. c     INICIO bucle elementos
  104. c*****************************************************************
  105.       DO 1000 IB=1,NBELEM
  106. c*****************************************************************
  107.  
  108.       CALL ZERO(REL,LRE,LRE)
  109.       CALL ZERO(REL2,LRE,LRE)
  110.  
  111. c coord del elemento ib (en t_(n+1) por el form)
  112.       CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XENEW)
  113.  
  114. c desplazamientos impuestos y coordenadas en t_n
  115.       MPTVAL=IVADESP
  116.       DO INO=1,NBNN
  117.         DO ID=1,IDIM
  118.            MELVAL=IVAL(ID)
  119.            IBMN=MIN(IB,VELCHE(/2))
  120.            INMN=MIN(INO,VELCHE(/1))
  121.            XMOVI(ID,INO)=VELCHE(INMN,IBMN)
  122.            XEOLD(ID,INO)=XENEW(ID,INO)-XMOVI(ID,INO)
  123.         ENDDO
  124.       ENDDO
  125. cc
  126. c*****************************************************************
  127. c     INICIO bucle puntos de gauss
  128. c*****************************************************************
  129.       DO 1100 IGAU=1,NBPTEL
  130. c*****************************************************************
  131.  
  132. *  var0 = variables internas iniciales
  133.          MPTVAL=IVARI
  134.          DO IC=1,NVARI
  135.             MELVAL=IVAL(IC)
  136.             IBMN=MIN(IB,VELCHE(/2))
  137.             IGMN=MIN(IGAU,VELCHE(/1))
  138.             VAR0(IC)=VELCHE(IGMN,IBMN)
  139.          enddo
  140.  
  141. *  defplini = deformacion inelastica inicial
  142.          MPTVAL=IVADPI
  143.          DO IC=1,NDEF
  144.             MELVAL=IVAL(IC)
  145.             IBMN=MIN(IB,VELCHE(/2))
  146.             IGMN=MIN(IGAU,VELCHE(/1))
  147.             DEFPINI(IC)=VELCHE(IGMN,IBMN)
  148.          enddo
  149.  
  150. *  xmat = caracteristicas materiales
  151.          MPTVAL=IVAMAT
  152.          DO IC=1,NMATT-5
  153.             MELVAL    =IVAL(IC)
  154.             IGMN      =MIN(IGAU,VELCHE(/1))
  155.             IBMN      =MIN(IB  ,VELCHE(/2))
  156.             XMAT(IC+2)=VELCHE(IGMN,IBMN)
  157.          ENDDO
  158.          XMAT(1)=XMAT(3)
  159.          XMAT(2)=XMAT(4)
  160.          XMAT(3)=0.D0
  161.          XMAT(4)=0.D0
  162. c  compatibilidad con feap, ogden model
  163.          if (nmodel.eq.8) then
  164.             XMAT(10)=XMAT(9)
  165.             XMAT(9) =0.D0
  166.             XMAT(13)=-1.D0
  167.          endif
  168.  
  169. ****************************************
  170. * integrar
  171. ****************************************
  172.  
  173.          CALL HPRIME(XENEW,NBNN,IDIM,SHPTOT,IGAU,SHPWRK,DJAC)
  174.  
  175. C Matriz Finv (d x_old / d x_new)
  176.          CALL ZERO(FFI,3,3)
  177.          DO ID=1,NBNN
  178.           DO IE=1,IDIM
  179.           DO IF=1,IDIM
  180.             FFI(IE,IF)=FFI(IE,IF)-SHPWRK(IF+1,ID)*XMOVI(IE,ID)
  181.           ENDDO
  182.           ENDDO
  183.          ENDDO
  184.          FFI(1,1)=1.D0+FFI(1,1)
  185.          FFI(2,2)=1.D0+FFI(2,2)
  186.          FFI(3,3)=1.D0+FFI(3,3)
  187.          IF (IDIM.EQ.2) THEN
  188.            FFI(3,3)=1.D0
  189.            IF (IFOUR.EQ.0) THEN
  190. c          Cas axisim
  191.              R1=0.D0
  192.              R2=0.D0
  193.              DO ID=1,NBNN
  194.                R1=R1+SHPWRK(1,ID)*XENEW(1,ID)
  195.                R2=R2+SHPWRK(1,ID)*XEOLD(1,ID)
  196.              ENDDO
  197.              FFI(3,3)=R2/(R1+1.E-20)
  198.            ENDIF
  199.          ENDIF
  200.  
  201. C Matrix F (d x_new / d x_old)
  202. C*       CALL ZERO(FF ,3,3)
  203.          DO I=1,3
  204.          DO J=1,3
  205.            FF(I,J)=FFI(I,J)
  206.          ENDDO
  207.          ENDDO
  208.          CALL INVALM(FF,3,3,KERRE,1.D-12)
  209.          IF (KERRE.NE.0) THEN
  210.             WRITE(*,*) ' MATRICE SINGULIERE'
  211.             RETURN
  212.          ENDIF
  213.  
  214. c Almacenamiento vectorial
  215.          DO i = 1,3
  216.            FFv(i) =FF(i,i)
  217.            FFIv(i)=FFI(i,i)
  218.          ENDDO
  219.          FFv(4) =FF(1,2)
  220.          FFv(5) =FF(2,1)
  221.          FFIv(4)=FFI(1,2)
  222.          FFIv(5)=FFI(2,1)
  223.          IF (IDIM.EQ.3) THEN
  224.           write(*,*)' 3D FeFp not available'
  225.           return
  226.          ENDIF
  227.  
  228. C Determinante de Finv
  229.          IF (IDIM.EQ.2) THEN
  230.           DETFFI=(FFI(1,1)*FFI(2,2)-FFI(1,2)*FFI(2,1))*FFI(3,3)
  231.          ELSEIF (IDIM.EQ.3) THEN
  232.           DETFFI=       FFI(1,1)*(FFI(2,2)*FFI(3,3)-FFI(3,2)*FFI(2,3))
  233.           DETFFI=DETFFI-FFI(2,1)*(FFI(1,2)*FFI(3,3)-FFI(3,2)*FFI(1,3))
  234.           DETFFI=DETFFI+FFI(3,1)*(FFI(1,2)*FFI(2,3)-FFI(1,3)*FFI(2,2))
  235.          ENDIF
  236.  
  237. ********************************
  238. c Initialize plastic internal variables in first time step
  239.          aux = abs(DEFPINI(1))+abs(DEFPINI(2))+abs(DEFPINI(3))
  240.          if (aux.lt.1.d-14) then
  241.            call zero(DEFPINI,ndef,1)
  242.            do i = 1,3
  243.              DEFPINI(i) = 1.d0
  244.            end do
  245.          endif
  246.  
  247. c Inicializar la variable interna 'det(FF_tot)' si es Update_Lagrangian
  248.          if ((nupdate.eq.1).AND.(VAR0(nvardet).lt.1.D-14))
  249.      .      VAR0(nvardet) = 1.D0
  250.  
  251. c Push-Forward plastic internal variables
  252.          call pushf35(DEFPINI,FFv,BE)
  253.  
  254. c Copy internal variables to t_n+1
  255.          do i = 1,NVARI
  256.             VARF(i) = VAR0(i)
  257.          end do
  258.  
  259. c  density dependent plasticity, POWDER & PODERCAP
  260.         xdensi=-1.D0
  261.         if (inplas.eq.116.or.inplas.eq.117) then
  262.          xdensi = VAR0(nvardet) * DETFFI
  263.         endif
  264.  
  265. c Compute KIRCHHOFF stress + tangent
  266.          call apf_driver_fefp(BE,VARF,SIGF,DDHOOK,
  267.      .                        NDEF,NVARI,NSTRS,LHOOK,
  268.      .                        XMAT,xdensi,PRECIS,NITMAX,KERRE,
  269.      .                        nescri,nues,nmodel,nnumer,deltax,
  270.      .                        level,kmax,iaugla,caugla,ib,igau,izone)
  271.  
  272. c Total Lagrangian
  273.          if (nupdate.eq.0) then
  274. c           Guardar deformaciones elasticas en referencia
  275.             call pushf35(BE,FFIv,DEFPFIN)
  276.  
  277. c Update Lagrangian
  278.          else
  279. c           Guardar deformaciones elasticas en t_(n+1)
  280.             call equv(DEFPFIN,BE,NDEF)
  281. c           Actualizar el determinante de FF^(-1)_(total)
  282.             VARF(nvardet) = VAR0(nvardet) * DETFFI
  283.             DETFFI        = VARF(nvardet)
  284.          endif
  285.  
  286. c Get CAUCHY stresses and tangent moduli
  287.          do i = 1,NSTRS
  288.            SIGF(i) = SIGF(i)*DETFFI
  289.            do j = 1,NSTRS
  290.              DDHOOK(i,j) = DDHOOK(i,j)*DETFFI
  291.            enddo
  292.          enddo
  293.  
  294. c add pressure contribution, log model (only deviatoric models)
  295. c al ser deviatoric models => det(Fe) = det (F)
  296. c por este motivo se puede calcular con este DETFFI
  297. c los casos no isocoricos calculan a partir de BE, que solo tiene
  298. c la parte elastica
  299.         if (nmodel.eq.8) then
  300.           if (abs(DETFFI).lt.1.D-14) DETFFI = 1.D-14
  301.           uppj  = xmat(1)/3.D0/(1.D0-2.D0*xmat(2)) * DETFFI
  302.           press = uppj * LOG(ABS(1.D0/DETFFI))
  303.           do i = 1,3
  304.            SIGF(i) = SIGF(i) + press
  305.            DDHOOK(i,i) = DDHOOK(i,i) - 2.d0*press + uppj
  306.           enddo
  307.           do i = 4,NSTRS
  308.            DDHOOK(i,i) = DDHOOK(i,i) - press
  309.           enddo
  310.           DDHOOK(1,2) = DDHOOK(1,2) + uppj
  311.           DDHOOK(1,3) = DDHOOK(1,3) + uppj
  312.           DDHOOK(2,3) = DDHOOK(2,3) + uppj
  313.         endif
  314.  
  315.         IF (IRIGE7.EQ.0)THEN
  316. c compute lower part by symmetry
  317.           do i = 2,NSTRS
  318.           do j = 1,i-1
  319.              DDHOOK(i,j) = DDHOOK(j,i)
  320.           enddo
  321.           enddo
  322.         ENDIF
  323. ****************************************
  324. * guardar resultados: Remember: Cauchy stresses
  325. ****************************************
  326.  
  327. *  sigf = tensiones finales
  328.             MPTVAL=IVASTF
  329.             DO IC=1,NSTRS
  330.                MELVAL=MPTVAL.IVAL(IC)
  331.                MELVAL.VELCHE(IGAU,IB)=SIGF(IC)
  332.             enddo
  333.  
  334. *  varf = variables internas finales
  335.             MPTVAL=IVARIF
  336.             DO IC=1,NVARI
  337.                MELVAL=MPTVAL.IVAL(IC)
  338.                MELVAL.VELCHE(IGAU,IB)=VARF(IC)
  339.             enddo
  340.  
  341. *  defpfin = deformations plasticas finales
  342.             MPTVAL=IVADPF
  343.             DO IC=1,NDEF
  344.                MELVAL=MPTVAL.IVAL(IC)
  345.                MELVAL.VELCHE(IGAU,IB)=DEFPFIN(IC)
  346.             enddo
  347.  
  348. ****************************************
  349. c calcula la matriz B = BGENE y el jacobiano DJAC
  350. ****************************************
  351.          XDPGE = 0.D0
  352.          YDPGE = 0.D0
  353.          DIM3  = 1.D0
  354.          CALL BMATST(IGAU,NBPTEL,POIGAU,QSIGAU,ETAGAU,DZEGAU,
  355.      .               MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NIFOUR,DIM3,
  356.      .               XENEW,SHPTOT,SHPWRK,BGENE,DJAC,XDPGE,YDPGE)
  357.          IF (abs(DJAC).LT.1.E-17) then
  358.             write(*,*)' Jacobiano cero, en elem', ib,' gauss',igau
  359.          ELSEIF (DJAC.LT.0.D0) then
  360.             write(*,*)' Jacobiano negativo, en elem', ib,' gauss',igau
  361.             return
  362.          endif
  363.          DJAC=ABS(DJAC)*MINTE.POIGAU(IGAU)
  364. c Ktan = BDB
  365.          IF (IRIGE7.EQ.2)THEN
  366.             CALL BDBSTS(BGENE,DJAC,DDHOOK,LRE,NSTRS,REL)
  367.          ELSE
  368.             CALL BDBST (BGENE,DJAC,DDHOOK,LRE,NSTRS,REL)
  369.          ENDIF
  370. c Ksigma
  371.          DO IA=1,NBNN
  372.           DO IO=1,6
  373.             SHPWRK(IO,IA)=SHPTOT(IO,IA,IGAU)
  374.           ENDDO
  375.          ENDDO
  376.          CALL DEVOLU(XENEW,SHPWRK,MFR,NBNN,IFOUR,
  377.      .               NIFOUR,IDIM,DIM3,RR,DJAC)
  378.          DJAC=ABS(DJAC)*MINTE.POIGAU(IGAU)
  379.          IF(IFOUR.EQ.1) THEN
  380.           IF(NIFOUR.EQ.0) THEN
  381.             CALL THSIG1(SHPWRK,DJAC,SIGF,NBNN,LRE,REL2,RR)
  382.           ELSE
  383.             CALL THSIG2(SHPWRK,DJAC,SIGF,NBNN,LRE,REL2,NIFOUR,RR)
  384.           ENDIF
  385.          ELSEIF(IFOUR.EQ.0) THEN
  386.           CALL THSIG3(SHPWRK,DJAC,SIGF,NBNN,LRE,REL2,RR)
  387.          ELSE
  388.           CALL THSIGH(SHPWRK,DJAC,SIGF,NBNN,IDIM,LRE,REL2)
  389.          ENDIF
  390. c*****************************************************************
  391. c     FIN fin bucle puntos de gauss
  392. 1100  continue
  393. c*****************************************************************
  394. c guarda la matriz de rigidez elemental REL en XMATRI.RE
  395. C*       NLIGRP=LRE
  396. C*       NLIGRD=LRE
  397.          IF (IRIGE7.EQ.2)THEN
  398.             DO IA=1,LRE
  399.             DO IIB=1,IA
  400.               RE(IA,IIB,ib)=REL(IA,IIB)+REL2(IA,IIB)
  401.               RE(IIB,IA,ib)=REL(IIB,IA)+REL2(IA,IIB)
  402.             ENDDO
  403.             ENDDO
  404.          ELSE
  405.             DO IA=1,LRE
  406.             DO IIB=1,IA
  407.               RE(IA,IIB,ib)=REL(IA,IIB)+REL2(IA,IIB)
  408.               RE(IIB,IA,ib)=RE(IA,IIB,ib)
  409.             ENDDO
  410.             ENDDO
  411.          ENDIF
  412. c*****************************************************************
  413. c     FIN bucle elementos
  414. 1000  continue
  415.       segdes xmatri
  416. c*****************************************************************
  417. c desactivar segmentos de trabajo
  418.       SEGSUP WRK0
  419.  
  420.       RETURN
  421.       END
  422.  
  423.  
  424.  

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