Télécharger la plaquette Cast3M


fr en

Télécharger cubhm1.eso

Retour à la liste

Numérotation des lignes : 

cubhm1

  1. C CUBHM1    SOURCE    CHAT      05/01/12    22:32:31     5004
  2. C CUBHM1    SOURCE    LOFT      90/01/26    21:08:44
  3.       SUBROUTINE CUBHM1(IGAU,ITEL,NBNO,XEL,SHPTOT,SHP,RHOS,RHOF,B11,B22,
  4.      #B12,SFLU,SCEL,POIGAU,VKL12,VKL22,VKL23,VKL33,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
  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     SFLU = SURFACE  FLUIDE  DANS LA CELLULE ELEMENTAIRE
  21. C     SCEL = SURFACE  DE LA CELLULE ELEMENTAIRE
  22. C     POIGAU=MINTE1.POIGAU(IGAU)
  23. C     VKL12=-((COEFPR*COEFPI)/(RHOF*C**2))*SFLU/SCEL
  24. C     VKL22=-(COEFPI**2)/(RHOF*SCEL)
  25. C     VKL23= COEFPI/SCEL
  26. C     VKL33= 1/SCEL
  27. C     LRE =NOMBRE DE D.D.L DE LA MATRICE  DE  RIGIDITE
  28. C     SHPTOT(6,NBNO,NBGAU)=FONCTIONS DE FORMES ET DERIVEES
  29. C     ISDJC= INDICATEUR  DU SIGNE DU  JACOBIEN
  30. C  ZONE DE TRAVAIL
  31. C     SHP(6,NBNO)=TABLEAU DE TRAVAIL
  32. C OUTPUT
  33. C     REL=MATRICE DE RIGIDITE
  34. C     ISDJC= INDICATEUR  DU SIGNE DU  JACOBIEN
  35. C     IRET=INDICATEUR = 1 : SUCCES
  36. C                     = 0 : ECHEC (ITEL INCOMPATIBLE AVEC LA FORMULATION
  37. C                     = 2 : ECHEC ( JACOBIEN NUL )
  38. C=======================================================================
  39.       IMPLICIT INTEGER(I-N)
  40.       IMPLICIT REAL*8 (A-H,O-Z)
  41.       DIMENSION XEL(3,*),SHP(6,*),SHPTOT(6,NBNO,*),REL(LRE,*)
  42.       IF (ITEL.EQ.127)              GOTO 10
  43. C
  44. C     ERREUR : TYPE D' ELEMENT  INCOMPATIBLE AVEC LA FORMULATION
  45. C
  46.       IRET = 0
  47.       GOTO 666
  48.   10  CONTINUE
  49. C
  50. C     SHP(1,I) : FONCTION DE FORME
  51. C     SHP(2,I) : DERIVEE % X DE LA FONCTION DE FORME
  52. C     SHP(3,I) : DERIVEE % Y DE LA FONCTION DE FORME
  53. C     SHP(4,I) : DERIVEE % Z DE LA FONCTION DE FORME
  54. C
  55.       DO 101 NP=1,NBNO
  56.       SHP(1,NP)=SHPTOT(1,NP,IGAU)
  57.       SHP(2,NP)=SHPTOT(2,NP,IGAU)
  58.       SHP(3,NP)=SHPTOT(3,NP,IGAU)
  59.       SHP(4,NP)=SHPTOT(4,NP,IGAU)
  60.   101 CONTINUE
  61.       IDIM=3
  62.       CALL JACOBI(XEL,SHP,IDIM,NBNO,DJAC)
  63.  
  64. C
  65. C
  66.       A1 = XEL(3,2) - XEL(3,1)
  67.       A1 = ABS(A1)
  68.       A2 = XEL(3,3) - XEL(3,1)
  69.       A2 = ABS(A2)
  70.       A3 = XEL(3,4) - XEL(3,1)
  71.       A3 = ABS(A3)
  72.       A4 = XEL(3,5) - XEL(3,1)
  73.       A4 = ABS(A4)
  74.       A5 = XEL(3,6) - XEL(3,1)
  75.       A5 = ABS(A5)
  76.       A6 = XEL(3,7) - XEL(3,1)
  77.       A6 = ABS(A6)
  78.       A7 = XEL(3,8) - XEL(3,1)
  79.       A7 = ABS(A7)
  80. C
  81.  
  82.       A3 = (0.25D0)*(A1+A2+A3+A4+A5+A6+A7)
  83. C
  84.       IF (DJAC.EQ.0.) GOTO 667
  85.       IRET = 1
  86.       IF (DJAC.LT.0.) ISDJC = ISDJC + 1
  87.       DJAC=ABS(DJAC)*POIGAU
  88. C      WRITE(6,*)'DJAC=',DJAC
  89. C
  90. C
  91. C     FONCTIONS DE FORME EN Z
  92. C
  93.       ZH=SHP(1,5)+SHP(1,6)+SHP(1,7)+SHP(1,8)
  94. C
  95.       Z1=1.D0-3.D0*ZH*ZH+2.D0*ZH*ZH*ZH
  96.       Z2=3.D0*(ZH**2)-2.D0*ZH*ZH*ZH
  97.       Z3=(A3*ZH)*(1.D0-2.D0*ZH+ZH*ZH)
  98.       Z4=(A3*ZH)*(ZH*ZH-ZH)
  99. C
  100. C
  101.       S1=SHP(1,1)+SHP(1,4)+SHP(1,5)+SHP(1,8)
  102.       S2=SHP(1,2)+SHP(1,3)+SHP(1,6)+SHP(1,7)
  103. C
  104.       T1=SHP(1,1)+SHP(1,2)+SHP(1,5)+SHP(1,6)
  105.       T2=SHP(1,3)+SHP(1,4)+SHP(1,7)+SHP(1,8)
  106. C
  107. C
  108. C     FONCTIONS DE FORME POUR LES DEPLACEMENTS
  109. C
  110.       SHP(5,1)=S1*T1*Z1
  111.       SHP(5,2)=S2*T1*Z1
  112.       SHP(5,3)=S2*T2*Z1
  113.       SHP(5,4)=S1*T2*Z1
  114.       SHP(5,5)=S1*T1*Z2
  115.       SHP(5,6)=S2*T1*Z2
  116.       SHP(5,7)=S2*T2*Z2
  117.       SHP(5,8)=S1*T2*Z2
  118.  
  119.       SHP(5,9)=S1*T1*Z3
  120.       SHP(5,10)=S2*T1*Z3
  121.       SHP(5,11)=S2*T2*Z3
  122.       SHP(5,12)=S1*T2*Z3
  123.       SHP(5,13)=S1*T1*Z4
  124.       SHP(5,14)=S2*T1*Z4
  125.       SHP(5,15)=S2*T2*Z4
  126.       SHP(5,16)=S1*T2*Z4
  127. C
  128. C     TERMES EN P*PI
  129. C     NBDL : NOMBRE DE D.D.L PAR NOEUD  ( = 6 )
  130. C
  131.       NBDL = 6
  132.       IX1=0
  133.       IY1=0
  134.       DO   102 IX=2,LRE ,NBDL
  135.       IX1=IX1 + 1
  136.       DO   103 IY=1,IX  ,NBDL
  137.       IY1=IY1 + 1
  138.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL12*SHP(1,IX1)*SHP(1,IY1)
  139.       REL(IX,IY) = REL(IY,IX)
  140.   103 CONTINUE
  141.       IY1=0
  142.   102 CONTINUE
  143.       DO   104 IX=2,LRE-NBDL,NBDL
  144.       IX1=IX +NBDL-1
  145.       DO   105 IY=IX1,LRE ,NBDL
  146.       REL(IX,IY) = REL(IX-1,IY+1)
  147.       REL(IY,IX) = REL(IX,IY)
  148.   105 CONTINUE
  149.   104 CONTINUE
  150. C
  151. C     TERMES EN PI*PI
  152. C
  153.       B33=SFLU
  154.       IX1=0
  155.       IY1=0
  156.       DO   106 IX=2,LRE ,NBDL
  157. C      WRITE(6,*)'IX= ',IX
  158.       IX1=IX1 + 1
  159. C      WRITE(6,*)'IX1= ',IX1
  160.       DO   107 IY=2,IX  ,NBDL
  161. C      WRITE(6,*)'IY= ',IY
  162.       IY1=IY1 + 1
  163. C      WRITE(6,*)'IY1= ',IY1
  164.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL22*(B11*SHP(2,IX1)*SHP(2,IY1)
  165.      #+B22*SHP(3,IX1)*SHP(3,IY1)+B33*SHP(4,IY1)*SHP(4,IX1))
  166.       REL(IX,IY) = REL(IY,IX)
  167.   107 CONTINUE
  168.       IY1=0
  169.   106 CONTINUE
  170. C
  171. C     TERMES EN PI*(UX,RY)
  172. C
  173.       D11 = (SCEL - B11)
  174.       D22 = (SCEL - B22)
  175.       D33=SFLU
  176. C
  177.       IX1=0
  178.       IY1=0
  179.       DO   108 IX=3,LRE ,NBDL
  180.       IX1=IX1 + 1
  181.       DO   109 IY=2,IX  ,NBDL
  182.       IY1=IY1 + 1
  183.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * (D11*SHP(2,IY1)*SHP(5,IX1))
  184.       REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL23 * (D11*SHP(2,IY1)
  185.      #*SHP(5,IX1+8))
  186.       REL(IX,IY) = REL(IY,IX)
  187.       REL(IX+1,IY) = REL(IY,IX+1)
  188.   109 CONTINUE
  189.       IY1=0
  190.   108 CONTINUE
  191.       IX1=1
  192.       IY1=0
  193.       DO   110 IX=2+NBDL,LRE ,NBDL
  194.       IX1=IX1 + 1
  195.       DO   111 IY=3,IX  ,NBDL
  196.       IY1=IY1 + 1
  197.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * D11*SHP(2,IX1)*SHP(5,IY1)
  198.       REL(IY+1,IX) = REL(IY+1,IX) + DJAC*VKL23*
  199.      #D11*SHP(2,IX1)*SHP(5,IY1+8)
  200.       REL(IX,IY) = REL(IY,IX)
  201.       REL(IX,IY+1) = REL(IY+1,IX)
  202.   111 CONTINUE
  203.       IY1=0
  204.   110 CONTINUE
  205. C
  206. C     TERMES EN PI*(UY,RX)
  207. C
  208.       IX1=0
  209.       IY1=0
  210.       DO   114 IX=5,LRE ,NBDL
  211.       IX1=IX1 + 1
  212.       DO   115 IY=2,IX  ,NBDL
  213.       IY1=IY1 + 1
  214.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * D22*SHP(3,IY1)
  215.      #* SHP(5,IX1)
  216.       REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL23 *
  217.      #D22*SHP(3,IY1) * SHP(5,IX1+8)*(-1.D0)
  218.       REL(IX,IY) = REL(IY,IX)
  219.       REL(IX+1,IY) = REL(IY,IX+1)
  220.   115 CONTINUE
  221.       IY1=0
  222.   114 CONTINUE
  223.       IX1=1
  224.       IY1=0
  225.       DO   116 IX=2+NBDL,LRE ,NBDL
  226.       IX1=IX1 + 1
  227.       DO   117 IY=5,IX  ,NBDL
  228.       IY1=IY1 + 1
  229.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL23 * D22*SHP(3,IX1)
  230.      #* SHP(5,IY1)
  231.       REL(IY+1,IX) = REL(IY+1,IX) + DJAC*VKL23 *
  232.      #D22*SHP(3,IX1) * SHP(5,IY1+8)*(-1.D0)
  233.       REL(IX,IY) = REL(IY,IX)
  234.       REL(IX,IY+1) = REL(IY+1,IX)
  235.   117 CONTINUE
  236.       IY1=0
  237.   116 CONTINUE
  238. C
  239. C     TERMES EN (UX,RY)*(UX,RY)
  240. C
  241.       H11=RHOS + RHOF*(SFLU-B11)
  242.       H22=RHOS + RHOF*(SFLU-B22)
  243.       IX1=0
  244.       IY1=0
  245.       DO   112 IX=3,LRE ,NBDL
  246.       IX1=IX1 + 1
  247.       DO   113 IY=3,IX  ,NBDL
  248.       IY1=IY1 + 1
  249.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1)
  250.       REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1+8)
  251.       REL(IY+1,IX) = REL(IY+1,IX)+DJAC*VKL33*H11*SHP(5,IY1+8)*SHP(5,IX1)
  252.       REL(IY+1,IX+1)=REL(IY+1,IX+1)+DJAC*VKL33*H22*SHP(5,IY1+8)
  253.      #*SHP(5,IX1+8)
  254.       REL(IX,IY) = REL(IY,IX)
  255.       REL(IX+1,IY) = REL(IY,IX+1)
  256.       REL(IX,IY+1) = REL(IY+1,IX)
  257.       REL(IX+1,IY+1) = REL(IY+1,IX+1)
  258.   113 CONTINUE
  259.       IY1=0
  260.   112 CONTINUE
  261. C
  262. C     TERMES EN (UY,RX)*(UY,RX)
  263. C
  264.       H11=RHOS + RHOF*(SFLU-B11)
  265.       H22=RHOS + RHOF*(SFLU-B22)
  266.       IX1=0
  267.       IY1=0
  268.       DO   118 IX=5,LRE ,NBDL
  269.       IX1=IX1 + 1
  270.       DO   119 IY=5,IX  ,NBDL
  271.       IY1=IY1 + 1
  272.       REL(IY,IX) = REL(IY,IX) + DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1)
  273.       REL(IY,IX+1) = REL(IY,IX+1)+DJAC*VKL33*H11*SHP(5,IY1)*SHP(5,IX1+8)
  274.      #*(-1.D0)
  275.       REL(IY+1,IX) = REL(IY+1,IX)-DJAC*VKL33*H11*SHP(5,IY1+8)*SHP(5,IX1)
  276.       REL(IY+1,IX+1)=REL(IY+1,IX+1)+DJAC*VKL33*H22*SHP(5,IY1+8)
  277.      #*SHP(5,IX1+8)
  278.       REL(IX,IY) = REL(IY,IX)
  279.       REL(IX+1,IY) = REL(IY,IX+1)
  280.       REL(IX,IY+1) = REL(IY+1,IX)
  281.       REL(IX+1,IY+1) = REL(IY+1,IX+1)
  282.   119 CONTINUE
  283.       IY1=0
  284.   118 CONTINUE
  285.       GOTO 666
  286. C
  287. C     MESSAGE D ERREUR : ELEMENT A SURFACE NULLE
  288. C
  289.   667 CONTINUE
  290.       IRET = 2
  291.       GOTO 666
  292. C
  293.   666 CONTINUE
  294.       RETURN
  295.       END
  296.  
  297.  
  298.  

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