Télécharger cubhm1.eso

Retour à la liste

Numérotation des lignes :

  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