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

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