Télécharger la plaquette Cast3M


fr en

Télécharger trihm31.eso

Retour à la liste

Numérotation des lignes : 

trihm31

  1. C TRIHM31   SOURCE    CHAT      05/01/13    03:47:13     5004
  2.       SUBROUTINE TRIHM31(IGAU,ITEL,NBNO,XEL,SHPTOT,SHP,RHOS,RHOF,B11
  3.      #,B22,B12,F11,F12,SFLU,SCEL,POIGAU,VKL12,VKL22,VKL23,VKL33,VKL42,
  4.      #VKL43,VKL44,E111,E112,E121,E122,E221,E222,LRE,
  5.      #REL,ISDJC,IRET)
  6. C=======================================================================
  7. C
  8. C    CALCULE LA MATRICE  DE  MASSE     DANS  LE  CAS  PLAN  POUR
  9. C        LA FORMULATION (37) HOMOGENE POUR ELEMENT TRH6
  10. C=======================================================================
  11. C     NBDL = NOMBRE DE DDL PAR NOEUD
  12. C  INPUT
  13. C     IGAU=NUMERO DU POINT DE GAUSS
  14. C     ITEL=NUMERO DE L ELEMENT DANS NOMTP
  15. C     NBNO=NOMBRE DE NOEUDS
  16. C     XEL =COORDONNEES  DE L ELEMENT
  17. C     RHOS = MASSE DU  SOLIDE
  18. C     RHOF = MASSE VOLUMIQUE DU  FLUIDE
  19. C     B11,B22,B12 =  PERMEABILITE ACOUSTIQUE DU MILIEU
  20. C     F11,F12  = NORME
  21. C     SFLU = SURFACE  FLUIDE  DANS LA CELLULE ELEMENTAIRE
  22. C     SCEL = SURFACE  DE LA CELLULE ELEMENTAIRE
  23. C     POIGAU=MINTE1.POIGAU(IGAU)
  24. C     VKL12=-((COEFPR*COEFPI)/(RHOF*C**2))*SFLU/SCEL
  25. C     VKL22=-(COEFPI**2)/(RHOF*SCEL)
  26. C     VKL23= COEFPI/SCEL
  27. C     VKL33= 1/SCEL
  28. C     LRE =NOMBRE DE D.D.L DE LA MATRICE  DE  RIGIDITE
  29. C     SHPTOT(6,NBNO,NBGAU)=FONCTIONS DE FORMES ET DERIVEES
  30. C     ISDJC= INDICATEUR  DU SIGNE DU  JACOBIEN
  31. C  ZONE DE TRAVAIL
  32. C     SHP(6,NBNO)=TABLEAU DE TRAVAIL
  33. C OUTPUT
  34. C     REL=MATRICE DE RIGIDITE
  35. C     ISDJC= INDICATEUR  DU SIGNE DU  JACOBIEN
  36. C     IRET=INDICATEUR = 1 : SUCCES
  37. C                     = 0 : ECHEC (ITEL INCOMPATIBLE AVEC LA FORMULATION
  38. C                     = 2 : ECHEC ( JACOBIEN NUL )
  39. C=======================================================================
  40.       IMPLICIT INTEGER(I-N)
  41.       IMPLICIT REAL*8 (A-H,O-Z)
  42.       DIMENSION XEL(3,*),SHP(6,*),SHPTOT(6,NBNO,*),REL(LRE,*)
  43.       IF (ITEL.EQ.157)              GOTO 10
  44. C
  45. C     ERREUR : TYPE D' ELEMENT  INCOMPATIBLE AVEC LA FORMULATION
  46. C
  47.       IRET = 0
  48.       GOTO 666
  49.   10  CONTINUE
  50. C
  51. C     SHP(1,I) : FONCTION DE FORME
  52. C     SHP(2,I) : DERIVEE % X DE LA FONCTION DE FORME
  53. C     SHP(3,I) : DERIVEE % Y DE LA FONCTION DE FORME
  54. C
  55. C     d²/dX²
  56. C
  57.       SHP(4,1)=4D0
  58.       SHP(4,2)=-8D0
  59.       SHP(4,3)=4D0
  60.       SHP(4,4)=0D0
  61.       SHP(4,5)=0D0
  62.       SHP(4,6)=0D0
  63. C
  64. C     d²/dXdY
  65. C
  66.       SHP(5,1)=4D0
  67.       SHP(5,2)=-4D0
  68.       SHP(5,3)=0D0
  69.       SHP(5,4)=4D0
  70.       SHP(5,5)=0D0
  71.       SHP(5,6)=-4D0
  72. C
  73. C     d²/dY²
  74. C
  75.       SHP(6,1)=4D0
  76.       SHP(6,2)=0D0
  77.       SHP(6,3)=0D0
  78.       SHP(6,4)=0D0
  79.       SHP(6,5)=4D0
  80.       SHP(6,6)=-8D0
  81. C
  82.       E211=E121
  83.       E212=E122
  84. C
  85.       DO 101 NP=1,NBNO
  86.       SHP(1,NP)=SHPTOT(1,NP,IGAU)
  87.       SHP(2,NP)=SHPTOT(2,NP,IGAU)
  88.       SHP(3,NP)=SHPTOT(3,NP,IGAU)
  89.   101 CONTINUE
  90.       IDIM=2
  91.       CALL JACOBI(XEL,SHP,IDIM,NBNO,DJAC)
  92.       IF (DJAC.EQ.0.) GOTO 667
  93.       IRET = 1
  94.       IF (DJAC.LT.0.) ISDJC = ISDJC + 1
  95.       DJAC=ABS(DJAC)*POIGAU
  96. C
  97. C     NBDL : NOMBRE DE D.D.L PAR NOEUD  ( = 4 )
  98.       NBDL = 4
  99. C
  100. C     TERMES EN P*PI
  101. C
  102.       IX1=0
  103.       IY1=0
  104.       DO   102 IX=2,LRE ,NBDL
  105.       IX1=IX1 + 1
  106.       DO   103 IY=1,IX  ,NBDL
  107.       IY1=IY1 + 1
  108.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL12*SHP(1,IX1)*SHP(1,IY1)
  109.       REL(IX,IY) = REL(IY,IX)
  110.   103 CONTINUE
  111.       IY1=0
  112.   102 CONTINUE
  113.       DO   104 IX=2,LRE-NBDL,NBDL
  114.       IX1=IX +NBDL-1
  115.       DO   105 IY=IX1,LRE ,NBDL
  116.       REL(IX,IY) = REL(IX-1,IY+1)
  117.       REL(IY,IX) = REL(IX,IY)
  118.   105 CONTINUE
  119.   104 CONTINUE
  120. C
  121. C     TERMES EN PI*PI
  122. C
  123.       IX1=0
  124.       IY1=0
  125.       DO   106 IX=2,LRE ,NBDL
  126.       IX1=IX1 + 1
  127.       DO   107 IY=2,IX  ,NBDL
  128.       IY1=IY1 + 1
  129. c      REL(IY,IX) = REL(IY,IX) + DJAC*VKL22*((B11*SHP(2,IX1)*SHP(2,IY1)+
  130. c     #B22*SHP(3,IX1)*SHP(3,IY1)+B12*(SHP(2,IX1)*SHP(3,IY1)+SHP(3,IX1)*
  131. c     #SHP(2,IY1)))) -DJAC*VKL42/RHOF*((E111*SHP(4,IX1)+E121*SHP(5,IX1))
  132. c     #*SHP(4,IY1)+(E112*SHP(5,IX1)+E122*SHP(6,IX1))*SHP(5,IY1)+(E211*
  133. c     #SHP(4,IX1)+E221*SHP(5,IX1))*SHP(5,IY1)+(E212*SHP(5,IX1)+E222*
  134. c     #SHP(6,IX1))*SHP(6,IY1))
  135.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL22*((B11*SHP(2,IX1)*SHP(2,IY1)+
  136.      #B22*SHP(3,IX1)*SHP(3,IY1)+B12*(SHP(2,IX1)*SHP(3,IY1)+SHP(3,IX1)*
  137.      #SHP(2,IY1)))) -DJAC*VKL42/RHOF*((E111*SHP(4,IX1)+E121*SHP(5,IX1))
  138.      #*SHP(4,IY1)+(E112*SHP(5,IX1)+E122*SHP(6,IX1)+E221*
  139.      #SHP(5,IX1)+E211*SHP(4,IX1))*SHP(5,IY1)+(E212*SHP(5,IX1)+E222*
  140.      #SHP(6,IX1))*SHP(6,IY1))
  141.       REL(IX,IY) = REL(IY,IX)
  142.   107 CONTINUE
  143.       IY1=0
  144.   106 CONTINUE
  145. C
  146. C     TERMES EN PI*(UX,UY)
  147. C
  148.       D11=(SCEL-B11)*F12
  149.       D22=(SCEL-B22)*F12
  150.       D12=-B12*F12
  151.       IX1=0
  152.       IY1=0
  153.       DO   108 IX=3,LRE ,NBDL
  154.       IX1=IX1 + 1
  155.       DO   109 IY=2,IX  ,NBDL
  156.       IY1=IY1 + 1
  157.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL23*(D11*SHP(2,IY1)+
  158.      #D12*SHP(3,IY1))*SHP(1,IX1) -DJAC*VKL43*(SHP(4,IY1)*E111*SHP(2,IX1)
  159.      #+SHP(5,IY1)*(E112*SHP(3,IX1)+E211*SHP(2,IX1))+SHP(6,IY1)*E212*
  160.      #SHP(3,IX1))
  161.       REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL23*(D12*SHP(2,IY1)+
  162.      #D22*SHP(3,IY1))*SHP(1,IX1)-DJAC*VKL43*(SHP(4,IY1)*E121*SHP(2,IX1)
  163.      #+SHP(5,IY1)*(E122*SHP(3,IX1)+E221*SHP(2,IX1))+SHP(6,IY1)*E222*
  164.      #SHP(3,IX1))
  165.       REL(IX,IY) = REL(IY,IX)
  166.       REL(IX+1,IY) = REL(IY,IX+1)
  167.   109 CONTINUE
  168.       IY1=0
  169.   108 CONTINUE
  170.       IX1=1
  171.       IY1=0
  172.       DO   110 IX=2+NBDL,LRE ,NBDL
  173.       IX1=IX1 + 1
  174.       DO   111 IY=3,IX  ,NBDL
  175.       IY1=IY1 + 1
  176.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * ( D11*SHP(2,IX1) +
  177.      #D12*SHP(3,IX1) ) * SHP(1,IY1)-DJAC*VKL43*(SHP(4,IX1)*E111
  178.      #*SHP(2,IY1)+SHP(5,IX1)*(E112*SHP(3,IY1)+E211*SHP(2,IY1))
  179.      #+SHP(6,IX1)*E212*SHP(3,IY1))
  180.       REL(IY+1,IX) = REL(IY+1,IX) + DJAC*VKL23 * (D12*SHP(2,IX1)+
  181.      #D22*SHP(3,IX1) ) * SHP(1,IY1)-DJAC*VKL43*(SHP(4,IX1)*E121
  182.      #*SHP(2,IY1)+SHP(5,IX1)*(E122*SHP(3,IY1)+E221*SHP(2,IY1))
  183.      #+SHP(6,IX1)*E222*SHP(3,IY1))
  184.       REL(IX,IY) = REL(IY,IX)
  185.       REL(IX,IY+1) = REL(IY+1,IX)
  186.   111 CONTINUE
  187.       IY1=0
  188.   110 CONTINUE
  189. C
  190. C     TERMES EN (UX,UY)*(UX,UY)
  191. C
  192.       H11=RHOS + RHOF*(SFLU-B11)*F11
  193.       H22=RHOS + RHOF*(SFLU-B22)*F11
  194.       H12=  (-1.)*RHOF*B12*F11
  195.       IX1=0
  196.       IY1=0
  197.       DO   112 IX=3,LRE ,NBDL
  198.       IX1=IX1 + 1
  199.       DO   113 IY=3,IX  ,NBDL
  200.       IY1=IY1 + 1
  201.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL33*H11*SHP(1,IY1)*SHP(1,IX1)
  202.      #-DJAC*VKL44*RHOF*(SHP(2,IX1)*E111*SHP(2,IY1)+
  203.      #SHP(3,IX1)*E112*SHP(3,IY1))
  204.       REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL33*H12*SHP(1,IY1)*SHP(1,IX1)
  205.      #-DJAC*VKL44*RHOF*(SHP(2,IX1)*E121*SHP(2,IY1)+
  206.      #SHP(3,IX1)*E122*SHP(3,IY1))
  207.       REL(IY+1,IX+1)=REL(IY+1,IX+1)+DJAC*VKL33*H22*SHP(1,IY1)*SHP(1,IX1)
  208.      #-DJAC*VKL44*RHOF*(SHP(2,IX1)*E221*SHP(2,IY1)+
  209.      #SHP(3,IX1)*E222*SHP(3,IY1))
  210.       REL(IX,IY) = REL(IY,IX)
  211.       REL(IY+1,IX) = REL(IY,IX+1)
  212.       REL(IX+1,IY) = REL(IY,IX+1)
  213.       REL(IX,IY+1) = REL(IY+1,IX)
  214.       REL(IX+1,IY+1) = REL(IY+1,IX+1)
  215.   113 CONTINUE
  216.       IY1=0
  217.   112 CONTINUE
  218.       GOTO 666
  219. C
  220. C     MESSAGE D ERREUR : ELEMENT A SURFACE NULLE
  221. C
  222.   667 CONTINUE
  223.       IRET = 2
  224.       GOTO 666
  225. C
  226.   666 CONTINUE
  227.       RETURN
  228.       END
  229.  
  230.  
  231.  

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