Télécharger la plaquette Cast3M


fr en

Télécharger cinc3d.eso

Retour à la liste

Numérotation des lignes : 

cinc3d

  1. C CINC3D    SOURCE    OF166741  25/11/04    21:15:23     12349          
  2. c A.Sellier avril jeu. 09 sept. 2010 11:55:47 CEST
  3.       SUBROUTINE cinc3d(WRK52,WRK53,WRK54,MWRKXE,
  4.      # nbnb,idimb,teta13d,teta23d,nvarib,nstrsb,ifourb,dtb)
  5.  
  6. C     WRK52,53,54 segments déclarés dans le common DECHE
  7.  
  8. C     XMAT(NCOMAT)  =  COMPOSANTES DE MATERIAU
  9. C     IVAL(NCOMAT)  =  INDICE DES COMPOSANTES DE MATERIAU
  10. C     NCOMAT        = NOMBRE DE COMPOSANTES DE MATERIAU
  11. C     XCAR(ICARA)   =  CARACTERISTIQUES
  12. C     MFR           =  NUMERO DE LA FORMULATION DE L'ELEMENT FINI
  13. C                   = 1       MASSIF
  14. C                   = 3       COQUE MINCE ( COQ2 , COQ3  ET DKT )
  15. C                   = 5       COQUE EPAISSE ( COQ6 , COQ8 )
  16. C                   = 7       POUTRE
  17. C                   = 9       COQUE MINCE AVEC CISAILLEMENT TRANSVERSE  ( COQ4 ET DST )
  18. C                   = 11      LIQUIDE
  19. C                   = 13      TUYAU
  20. C                   = 15      LINESPRING
  21. C                   = 17      TUYAU FISSURE
  22. C                   = 19      RACCORD MASSIF
  23. C                   = 21      RACCORD COQUE
  24. C                   = 23      SURFACE LIBRE
  25. C                   = 25      MEMBRANE
  26. C                   = 27      UNIAXIALE
  27. C                   = 29      THERMIQUE
  28. C                   = 31      INCOMPRESSIBLES
  29. C                   = 33      POREUX
  30. C                   = 35      JOINT
  31. C                   = 37      HOMOGENEISE
  32. C                   = 39      TUYO
  33. C                   = 41      TUYAU ACOUSTIQUE PURE
  34. C                   = 43      RACCORD TUYAU FLUIDE
  35. c*    DDAUX = MATRICE DE HOOKE ELASTIQUE
  36. c*    NSTRS = NBRE DE COMPOSANTES DES DEFORMATIONS
  37. c*    CMATE = NOM DU MATERIAU
  38. c*    VALMAT= TABLEAU DE CARACTERISTIQUES DU MATERIAU
  39. c*    VALCAR= TABLEAU DE CARACTERISTIQUES GEOMETRIQUES
  40. c*    N2EL  = NBRE D ELEMENTS DANS SEGMENT DE HOOKE
  41. c*    N2PTEL= NBRE DE POINTS DANS SEGMENT DE HOOKE
  42. c*    MFR   = NUMERO DE LA FORMULATION
  43. c*    IFOU  = type de formulation
  44. c*    IB    = NUMERO DE L ELEMENT COURANT
  45. c*    IGAU  = NUMERO DU POINT COURANT
  46. c*    EPAIST= EPAISSEUR
  47. c*    NBPGAU= NBRE DE POINTS DE GAUSS
  48. c*    MELE  = NUMERO DE L ELEMENT FINI
  49. c*    NPINT = NBRE DE POINTS D INTEGRATION
  50. c*    NBGMAT= NBRE DE POINTS DANS SEGMENT DE CARACTERISTIQUES
  51. c*    NELMAT= NBRE D ELEMENTS DANS SEGMENT DE CARACTERISTIQUES
  52. c*    SECT  = SECTION
  53. c*    LHOOK = TAILLE DE LA MATRICE DE HOOKE
  54. c*    TXR,XLOC,XGLOB,D1HOOK,ROTHOO,DDHOMU,CRIGI = TABLEAUX UTILISES
  55. c*    POUR LE CALCUL DE LA MATRICE DE HOOKE
  56. C-----------------------------------------------------------------------
  57. C      VARIABLES PASSEES PAR LES COMMONS COPTIO , ECOU  ET NECOU
  58. C
  59. C  IFOUR    INDICE DU TYPE DE PROBLEME
  60. C            -3  DEFORMATIONS PLANES GENERALISEES
  61. C            -2  CONTRAINTES PLANES
  62. C            -1  DEFORMATIONS PLANES
  63. C             0  AXISYMETRIQUE
  64. C             1  SERIE DE FOURIER
  65. C             2  TRIDIMENSIONNEL
  66. C  ITYP      TYPE DE FORMULATION MECANIQUE
  67. C  --------------- ATTENTION ---------------
  68. C     IL EST ACTIF APRES L APPEL DE VISAVI
  69. C  -----------------------------------------
  70. C       ITYP=1   CAS DES ELEMENTS MASSIFS
  71. C       ITYP=2   CAS DES COQUES
  72. C       ITYP=3   CAS DES MEMBRANES
  73. C       ITYP=4   CAS DES CABLES ET DES BARRES
  74. C       ITYP=5   CAS QUELCONQUE
  75. C       ITYP=6   CAS DES CONTRAINTES PLANES
  76. C       ITYP=7   CAS DES COQUES A NU=0. OU CONTRAINTES PLANES
  77. C       ITYP=8   CAS DES MEMBRANES A NU=0. OU CONTRAINTES PLANES
  78. C       ITYP=9   CAS DES COQUES EPAISSES
  79. C       ITYP=10  CAS DES JOINTS
  80. C       ITYP=11  CAS DES POUTRES
  81. C       ITYP=12  CAS DES TUYAUX
  82. C       ITYP=13  CAS DES COQUES AVEC CISAILLEMENT TRANSVERSE
  83. C
  84. C     ISTEP      flag utilise pour separer les etapes dans un calcul non local
  85. C                ISTEP=0 -----> calcul local
  86. C                ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
  87. C                ISTEP=2 -----> calcul non local etape 2 on continue le calcul
  88. C                               a partir des seuils moyennes
  89. C                ISTEP est défini par DEFINI.ESO
  90. C
  91. C     EPIN0      déformations inélastiques initiales
  92. C     EPINF      déformations inélastiques finales
  93. C     epst0      déformation totale initiale
  94. C     EPSTF      déformation totale finale
  95. C
  96. C-----------------------------------------------------------------------
  97. C  SORTIES
  98. C     SIGF(NSTRS)   = CONTRAINTES FINALES
  99. C     VARF(NVARI)   = VARIABLES INTERNES FINALES
  100. C     DEFP          = DEFORMATIONS PLASTIQUES
  101. C     KERRE         = 0    TOUT OK
  102.  
  103.       IMPLICIT INTEGER(I-N)
  104.       IMPLICIT REAL*8(A-H,O-Z)
  105.  
  106. -INC PPARAM
  107. -INC CCOPTIO
  108. -INC DECHE
  109.  
  110. c     segment de coordonnees des noeuds de l element
  111.       integer nbnb,insb,idimb
  112.       SEGMENT MWRKXE
  113.         REAL*8 XE(3,nbnb)
  114.       ENDSEGMENT
  115.  
  116. c     segment des numero des variables a moyenner
  117.       integer NBVIA
  118.       SEGMENT WR14
  119.         INTEGER INLVIA(NBVIA)
  120.       ENDSEGMENT
  121.  
  122. c     temperatures debut et fin de pas , moyenne, pas de temps, volume rgi
  123.       real*8 teta13d,teta23d,temp3d,dt3d
  124. c     variables de transfert de donnees ( a declarer suivant idvar4 et idvisc)
  125.       integer nstrs3d,nvari3d,ierr1,mfr11
  126. c     variable pour passer le numero de l etape non locale
  127.       integer istep3d ,nvarib,nstrsb,ifourb
  128.  
  129. c     ******** dimension des tableaux parametres materiau **************
  130.  
  131. c     nombre de parametres materiau du modele (cf coherence avec idvisc)
  132.       integer NBELAS3D,NBPARC3D,NBRINC3D,NBPPARP3D,NBRTAIL3D
  133. c     nmat1 nombre de parametres materiau sans la taille nmat3d-nbrtail3d
  134.       integer NMAT1,NMAT3D
  135. c     nombre de parametres materiaux
  136.  
  137. c     nombre d inclusions maximum **************************************
  138. c     a regeler aussi dans idvisc et idvar4
  139.       parameter (NBRINC3D=1)
  140. c     ******************************************************************
  141.  
  142. c     nombre de parametres par phases **********************************
  143. c     NBELAS3D nb elastique commun a tout les modeles meca (cf idvisc.eso)
  144. c     NBPARC3D nb de parametres commun
  145. c     NBPPARP3D nb de parametres par phase
  146.       parameter (NBELAS3D=4,NBPARC3D=4,NBPPARP3D=26)
  147. c     nb parametres materiaux maxi pour inclure la taille des elements finis
  148.       parameter (NBRTAIL3D=1)
  149. c     nbre de parametres materiaux sans les tailles des elements
  150.       parameter (NMAT1=NBELAS3D+NBPARC3D+(NBRINC3D+1)*NBPPARP3D)
  151. c     nombre total de parametres
  152.       parameter (NMAT3D=NMAT1+NBRTAIL3D+3)
  153. c     declaration du nombre totale de materiaux
  154.       real*8 xmat3d(NMAT3D)
  155. c     ******************************************************************
  156.  
  157. c     ******** dimension des tableaux variables internes ***************
  158. c     nombre de variables globales isotropes
  159.       parameter (NBVARISOG=1)
  160. c     nombre de variable globales tensorielles
  161.       parameter (NBVARTENG=2)
  162. c     nombre total de variables globales--------------------------------
  163.       parameter (NBVIND3D=NBVARISOG+6*NBVARTENG)
  164. c     ---cas des variables existantes pour toutes les phases -----------
  165. c     nombre de variable isotropes
  166.       parameter (NBVARISOPP=7)
  167. c     nombre de variable tensorielle
  168.       parameter (NBVARTENPP=8)
  169. c     nombre de variables par phase commune a toutes les phases
  170.       parameter (NBVPARP3D=NBVARISOPP+6*NBVARTENPP)
  171. c     --sous total des varis globales et communes-----------------------
  172.       parameter (NVTOT1=NBVIND3D+(NBRINC3D+1)*NBVPARP3D)
  173. c     nombre de variables supplementaires que pour les inclusions-------
  174. c     isotrope
  175.       parameter (NBVARISOPI=0)
  176. c     nombre de variables tensorielles
  177.       parameter (NBVARTENPI=12)
  178. c     nombre total de variables propres aux interfaces
  179.       parameter (NBVPARI3D=NBVARISOPI+6*NBVARTENPI)
  180. c     nombre de variables totales
  181.       parameter (NVARI3D=NVTOT1+NBVPARI3D*NBRINC3D)
  182. c     declaration des tableaux locaux des  variables internes ---------
  183.       real*8 VAR03D(NVARI3D),VARF3D(NVARI3D)
  184. c     ******************************************************************
  185.  
  186. c     ******** dimensionde la table des coordonnees des noeuds *********
  187. c     nbr de noeuds maxi par element
  188.       integer NBNMAX3D
  189.       parameter (NBNMAX3D=20)
  190.       real*8 XE3D(3,NBNMAX3D)
  191.  
  192. c     ********  taille des pseudo vecteurs *****************************
  193. c     tjrs 6 en raison de son utilisation dans incl3d
  194.       parameter (nstrs3d=6)
  195. c     contrainteq
  196.       real*8 sig03d(nstrs3d),sigf3d(nstrs3d)
  197. c     deformations
  198.       real*8 epst03d(nstrs3d),epstf3d(nstrs3d)
  199. c     increments de deformations
  200.       real*8 depst3d(nstrs3d)
  201. c     pas de temps
  202.       real*8 dtb
  203.  
  204. c     variables logiques :
  205. c     isotropie initiale,
  206.       logical iso1
  207. c     traitement de Gf en  local
  208.       logical local11
  209.  
  210. c     ifou type d element cf au debut du fichier
  211. C            -1  DEFORMATIONS PLANES
  212. C             0  AXISYMETRIQUE
  213. C             2  TRIDIMENSIONNEL
  214.       integer ifou11
  215. c     taille du psudo vecteur des tailles
  216.       integer ntail3d
  217.  
  218. c     dimension de l espace physique (1d, 2d ou 3d)
  219.       integer idimb3d
  220.  
  221. c     variable a vrai si il est effectivement prevu de passer
  222.  
  223. c***********************************************************************
  224. c remarque : parametres elastiques en formulation poreux massif (mfr=33)
  225. c    en debut de xmat :
  226. c   'YOUN'  :  module d'Young
  227. c   'NU  '  :  coefficient de poisson
  228. c   'RHO '  :  masse volumique (a la fin en non poreux)
  229. c   'ALPH'  :  coefficient de dilatation thermique
  230. c    à la fin de xmat :
  231. c   'MOB '  :  module de Biot
  232. c   'COB '  :  coefficient de Biot
  233. c   'PERM'  :  perméabilité intrinsèque
  234. c   'VISC'  :  viscosité dynamique du fluide
  235. c   'ALPM'  :  coefficient de couplage pression - température
  236. c   l ordre n est pas le meme si le materiau est orthotrope
  237. c   cf deche.inc pour recuperer les noms exacts des variables
  238. c   transfert des variables dans les tableaux du modele
  239. c***********************************************************************
  240.  
  241. c     ***** valeur de ntail3d
  242.       ntail3d=1
  243.  
  244. c***********************************************************************
  245. c      transfert des parametres materiaux dans le tableau de as3d
  246.        if(mfr.ne.33) then
  247. c       on est pas en poro meca
  248.         if(NMAT3D.lt.nmatt) then
  249.              print*,'pb dimension de xmat dans cinc3d'
  250.              print*,'nmatt',nmatt,'nmat3d',NMAT3D,'nmat1',NMAT1
  251.              do i=1,nmatt
  252.                  print*,'xmat(',i,')=',xmat(i)
  253.              enddo
  254.              read*
  255.              kerre=1
  256.              return
  257.         end if
  258. c       recuperation des parametres materiau
  259.         do i=1,NMAT1
  260.              xmat3d(i)=xmat(i)
  261.         end do
  262. c       recuperation des parametres de taille
  263.         if(ntail3d.ne.0) then
  264.          do i=(NMAT1+1),(NMAT1+ntail3d)
  265.              xmat3d(i)=xmat(i)
  266.          end do
  267.          if (ntail3d.lt.NBRTAIL3D) then
  268.           do i=(NMAT1+ntail3d+1),(NMAT1+NBRTAIL3D)
  269.              xmat3d(i)=0.d0
  270.           end do
  271.          else
  272.           if(ntail3d.gt.NBRTAIL3D) then
  273.             print*,'pb de dimension de ntail3d ds cinc3d'
  274.             kerre=1
  275.             return
  276.           end if
  277.          end if
  278.         else
  279.           do i=(NMAT1+1),(NMAT1+NBRTAIL3D)
  280.              xmat3d(i)=0.d0
  281.           end do
  282.         end if
  283.        else
  284. c       traiter ici la formulation poreux suivant les formulations
  285.         if(NMAT3D.lt.(nmatt-8)) then
  286.          print*,'on est en poreux'
  287.          print*,'Pb affectation des caracteristiques dans cinc3d'
  288.          do i=1,nmatt
  289.              print*,'xmat(',i,')=',xmat(i)
  290.          end do
  291.          kerre=1
  292.          return
  293.         else
  294. c         recuperation des parametres dans le cas poro meca
  295.           do i=1,NMAT3D
  296.              xmat3d(i)=xmat(i)
  297.           end do
  298.         end if
  299.        end if
  300.  
  301. c      **** recuperation des variables internes debut de pas ***********
  302.        if(nvari3d.ne.nvarib) then
  303.              print*,'Pb dimension de var0 dans cinc3d'
  304.              print*,'Verifier dimension table des variables internes'
  305.              print*,'et compatibilite avec idvar4 et incl3d'
  306.              read*
  307.              kerre=1
  308.              return
  309.        end if
  310. c      transfert vers le tableau de fldo3d
  311.        do i=1,nvari3d
  312. c            indicateur de passage de 1er pas
  313.              if(var0(1).ne.1.) then
  314.                  var0(i)=0.d0
  315.              end if
  316. c            autres variables internes
  317.              var03d(i)=var0(i)
  318. c            initialisation de variable internes finales aux initiales
  319.              varf3d(i)=var03d(i)
  320.        end do
  321.  
  322. c      ***** recuperation des contraintes totales debut de pas *********
  323.        if(mfr.ne.33) then
  324. c       on est pas en poreux
  325.         if(nstrsb.lt.nstrs3d) then
  326.             do i=1,nstrsb
  327.                 sig03d(i)=sig0(i)
  328.                 sigf3d(i)=sigf(i)
  329. c               !!!! ATTENTION les depst 3-6 doivent etre des gamas !!!!
  330.                 depst3d(i)=depst(i)
  331.                 epst03d(i)=epst0(i)
  332.                 epstf3d(i)=epstf(i)
  333.             end do
  334.             do i=nstrsb+1,nstrs3d
  335.                 sig03d(i)=0.d0
  336.                 sigf3d(i)=0.d0
  337. c               !!!! ATTENTION les depst 3-6 doivent etre des gamas !!!!
  338.                 depst3d(i)=0.d0
  339.                 epst03d(i)=0.d0
  340.                 epstf3d(i)=0.d0
  341.             end do
  342.         else
  343.             do i=1,nstrs3d
  344.                 sig03d(i)=sig0(i)
  345.                 sigf3d(i)=sigf(i)
  346. c               !!!! ATTENTION les depst 3-6 doivent etre des gamas !!!!
  347.                 depst3d(i)=depst(i)
  348.                 epst03d(i)=epst0(i)
  349.                 epstf3d(i)=epstf(i)
  350.             end do
  351.         end if
  352.        else
  353.          print*,'On est en poreux'
  354.          print*,'Pb affectation des contraintes dans cinc3d'
  355.          print*,'a terminer...'
  356.          do i=1,nstrsb
  357.              print*,'sig(',i,')=',sig0(i)
  358.          end do
  359.          kerre=1
  360.          return
  361.        end if
  362.  
  363.  
  364. c      ********   autres parametres a renseigner ***********************
  365.  
  366. c      initialisation indicateur d erreur
  367.        ierr1=0
  368. c      indicateur isostropie initiale (elastique et de resistance)
  369.        iso1=.true.
  370. c      numero de la formulation (33 pour poreux)
  371.        mfr11=mfr
  372. c      type de formulation
  373.        ifour11=ifourb
  374. c      controle de regularisation en cas d endommagement par methode locale
  375.        local11=.true.
  376. c      pas de temps
  377.        dt3d=dtb
  378. c       if(dt3d.eq.0.) then
  379. c        print*,'ds CINC3D: dt t0 tf', dt3d,tempf,temp0
  380. c       end if
  381. c      numero  pour le traitement local/non local eventuel
  382.        istep3d=istep
  383. c      test de compatibilite
  384.        if(local11.and.(istep3d.gt.0)) then
  385.             print*,'Dans inclusion3d, on ne peut pas avoir local=vrai'
  386.             kerre=1
  387.             return
  388.        end if
  389.  
  390. c      ******* recuperation des coordonnees de noeuds ******************
  391. c      boucle sur les noeuds
  392.        if(nbnb.le.NBNMAX3D) then
  393.            NBNB3D=nbnb
  394.            if(idimb.le.3) then
  395.               idimb3d=idimb
  396.            else
  397.               print*,'pb table de coord des neouds ds cflu3d'
  398.               kerre=1
  399.               return
  400.            end if
  401.        else
  402.           print*,'Element avec + de noeuds que capacite de cflu3d'
  403.           kerre=1
  404.           return
  405.        end if
  406. c      chargement des coordonnees des  noeuds
  407.        do insb=1,NBNB3D
  408.           do j=1,idimb3d
  409. c            print*,'xe(',j,insb,')=',xe(j,insb)
  410.              xe3d(j,insb)=xe(j,insb)
  411.           end do
  412.           if (idimb.lt.3) then
  413.             do j=idimb+1,3
  414.               xe3d(j,insb)=0.d0
  415.             end do
  416.           end if
  417.        end do
  418. c      on double les noeuds si le probleme est 2D
  419.        if ((ifourb.eq.0).or.(ifourb.eq.-1)) then
  420. c        on double les noeuds avec dimension 3 recupere en xmat(nmat1+8)
  421.          dimension3=xmat3d(nmat1+1)
  422. c        verif presence de la 3eme dimension
  423.          if(dimension3.eq.0.) then
  424.            if(ifourb.eq.0) then
  425. c             cas axi
  426.               print*,'Ds Cincl3d il manque la dimension DIMT'
  427.               kerre=1
  428.               return
  429.            else if (ifourb.eq.-1) then
  430. c             cas def plane
  431.               print*,'Ds Cincl3d il manque la dimension DIMZ'
  432.               kerre=1
  433.               return
  434.            else
  435. c             cas inconnu d incl3d
  436.               print*,'Ds Cincl3d EF incompatible'
  437.               kerre=1
  438.               return
  439.            end if
  440.          end if
  441.          NBNB3D=2*nbnb
  442.          do insb=nbnb+1,NBNB3D
  443.             do j=1,2
  444. c               print*,'xe(',j,insb,')=',xe(j,insb)
  445.                 xe3d(j,insb)=xe3d(j,(insb-nbnb))
  446.             end do
  447.             xe3d(3,insb)=dimension3
  448.          end do
  449.        end if
  450. c      ******* appel a incl3d *************************************
  451. c
  452. c      ne pas effacer nom des vraibles originales (pour memoire)
  453. c      call fluag3D(xmat,nmatt,var0,varf,nvari,dt1,DEPST,
  454. c     #nstrs,sigf,mfr,errfluag3D,teta1,teta2)
  455.  
  456.        call incl3d(xmat3d,NMAT3D,sig03d,sigf3d,depst3d,NSTRS3D,
  457.      # var03d,varf3d,NVARI3D,NBELAS3D,teta13d,teta23d,dt3d,ierr1,iso1,
  458.      # MFR11,local11,NMAT1,ifour11,istep3d,epst03d,epstf3d,NBRTAIL3D,
  459.      # NBNMAX3D,NBNB3D,IDIMB3D,XE3D,NBRINC3D,NBPARC3D,NBPPARP3D,
  460.      # NBVARISOG,NBVARTENG,NBVARISOPP,NBVARTENPP,NVTOT1,NBVARISOPI,
  461.      # NBVARTENPI)
  462.  
  463. c      ****** re-affectation des tableaux de variables internes ********
  464.        do i=1,nvarib
  465.             varf(i)=varf3d(i)
  466. c            print*,'cas3d f(',i,')=',varf(i),' i',var0(i)
  467.        end do
  468. c       read*
  469. c      ecriture des contraintes totales fin de pas
  470.        do i=1,nstrsb
  471.           sigf(i)=sigf3d(i)
  472.        end do
  473.  
  474. c***********************************************************************
  475. c       traitement de l erreur d ecoulement
  476.         if(ierr1.eq.0)then
  477.             kerre=0
  478.         else
  479.             kerre=1
  480.         end if
  481. 1000  continue
  482.       return
  483.       end
  484.  
  485.  
  486.  

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