Télécharger cfl5.eso

Retour à la liste

Numérotation des lignes :

  1. C CFL5 SOURCE BP208322 15/06/22 21:15:55 8543
  2. SUBROUTINE CFL5(ICAS,IPMAIL,MELE,IVAM1,IVAM2,MELV1,MELV2,N2)
  3. *---------------------------------------------------------------------*
  4. *
  5. * calcul du pas de temps CFL
  6. *
  7. * elements poutre,tuyau,barre,cerc,coque2,coque3,dkt
  8. *
  9. *
  10. * entree
  11. * icas : cas à traiter
  12. * = 1 calcul du pas de temps complet ivam1 avec cara
  13. * = 2 calcul du pas de temps lorsque cson est donne ivam2
  14. * = 3 calcul du pas de temps lorsque la taille est donnée ivam1 si cara
  15. * = 4 calcul de la vitesse du son ivam1 donné
  16. * = 5 calcul du parametre de taille ivam1 si cara
  17. *
  18. * ipmail : pointeur vers le maillage a traiter
  19. * mele : numero de l'élément finis dans nomtp
  20. * ivam1 : pointeur vers mptval du cham1 actif
  21. * ivam2 : pointeur vers mptval du cham2 actif
  22. * n2 : nombre de composante en sortie
  23. *
  24. * sortie
  25. * melv1 : melval de la première composante du chamelem resultat
  26. * inactif en sortie
  27. * melv2 : melval de la deuxième composante du chamelem resultat
  28. * inactif en sortie
  29. *
  30. *
  31. *---------------------------------------------------------------------*
  32. IMPLICIT INTEGER(I-N)
  33. IMPLICIT REAL*8(A-H,O-Z)
  34. *
  35. DIMENSION V13(3),V12(3),V23(3)
  36. *
  37. -INC CCOPTIO
  38. -INC CCREEL
  39. -INC CCHAMP
  40. -INC SMCHAML
  41. -INC SMELEME
  42. -INC SMCOORD
  43. *
  44. *
  45. SEGMENT MPTVAL
  46. * ipos pointeur vers la sous zone du mchelm
  47. * nsof ??
  48. INTEGER IPOS(NS) ,NSOF(NS)
  49. * ival pointeur vers le melval de la composante
  50. * =0 si il n'est pas présente
  51. * ncosou = nbrfac + nbrobl
  52. INTEGER IVAL(NCOSOU)
  53. * continent le type de composante
  54. CHARACTER*16 TYVAL(NCOSOU)
  55. ENDSEGMENT
  56. POINTEUR MPTVA1.MPTVAL,MPTVA2.MPTVAL
  57. *
  58. *
  59. *
  60. MPTVA1 = IVAM1
  61. MPTVA2 = IVAM2
  62. *
  63. *
  64. * branchement en fonction de l'élément fini
  65. *
  66. * 0 5 0 5 0
  67. GOTO (99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  68. 2 99,99,99,99,99,99,27,27,29,99,99,99,99,99,99,99,99,99,99,99,
  69. * bar
  70. 4 99,29,99,44,99,46,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  71. 6 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  72. 8 99,99,99,29,99,99,99,99,99,99,99,99,27,99,99,99,99,99,99,99,
  73. 1 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,
  74. 2 99,99,99,99,99,99,99),MELE
  75. *
  76. * cas de la poutre POUT (mince) et de TIMO (épaisse)
  77. *
  78. 29 CONTINUE
  79. * ================ calcul de la vitesse du son
  80. * le resultat est stocké dans melval avec n1ptel =1
  81. IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN
  82. * recherche des paramètre matériau
  83. * module d'young
  84. MELVA3 = MPTVA1.IVAL(1)
  85. * densite
  86. MELVA4 = MPTVA1.IVAL(3)
  87. SEGACT MELVA3,MELVA4
  88. *
  89. N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2))
  90. N1PTEL = 1
  91. N2PTEL = 0
  92. N2EL = 0
  93. SEGINI MELVAL
  94. MCSON = MELVAL
  95. * boucle sur les éléments pour calculer la vitesse du son
  96. DO 2903 I=1,N1EL
  97. * on prend les valeurs moyennes sur les éléments
  98. YOU1 = 0.D0
  99. I3 = MIN(I,MELVA3.VELCHE(/2))
  100. DO 2901 J=1,MELVA3.VELCHE(/1)
  101. YOU1 = YOU1 + MELVA3.VELCHE(J,I3)
  102. 2901 CONTINUE
  103. YOU1 = YOU1 / MELVA3.VELCHE(/1)
  104. *
  105. RO1 = 0.D0
  106. I4 = MIN(I,MELVA4.VELCHE(/2))
  107. DO 2902 J=1,MELVA4.VELCHE(/1)
  108. RO1 = RO1 + MELVA4.VELCHE(J,I4)
  109. 2902 CONTINUE
  110. RO1 = RO1 / MELVA4.VELCHE(/1)
  111. IF (RO1.EQ.0.D0) THEN
  112. SEGDES MELVA4,MELVA3
  113. SEGSUP MELVAL
  114. CALL ERREUR(855)
  115. RETURN
  116. ENDIF
  117. *
  118. IF (YOU1.EQ.0.D0) THEN
  119. SEGDES MELVA4,MELVA3
  120. SEGSUP MELVAL
  121. CALL ERREUR(856)
  122. RETURN
  123. ENDIF
  124. VELCHE(1,I) = SQRT(YOU1/RO1)
  125. * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i)
  126. 2903 CONTINUE
  127. SEGDES MELVA4,MELVA3
  128. IF (ICAS.EQ.4) THEN
  129. * cas ou seule la vitesse du son est demandée
  130. MELV2 = 0
  131. MELV1 = MELVAL
  132. SEGDES MELVAL
  133. RETURN
  134. ENDIF
  135. ELSE IF (ICAS.EQ.2) THEN
  136. * recuperation du champ
  137. * SEGACT MPTVA2
  138. MELVA1 = MPTVA2.IVAL(1)
  139. SEGACT MELVA1
  140. MCSON = MELVA1
  141. SEGDES MPTVA2
  142. ENDIF
  143. * ================ paramètre geometrique
  144. * stocké dans un melval mtaille
  145. IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN
  146. * on verifie que le champ de caractéristiques existe
  147. IF (IVAM1.EQ.0) THEN
  148. MOTERR(1:8)='CARACTER'
  149. CALL ERREUR(145)
  150. RETURN
  151. ENDIF
  152. MELEME = IPMAIL
  153. SEGACT MELEME
  154. N1EL = NUM(/2)
  155. N1PTEL = 1
  156. N2PTEL = 0
  157. N2EL = 0
  158. SEGINI MELVAL
  159. MTAIL1 = MELVAL
  160. CALL CFLTAI(MTAIL1,IPMAIL,MELE)
  161. SEGDES MELEME
  162. * second paramètre geométrique
  163. IOBL = 3
  164. IF (ICAS.EQ.5 .OR. ICAS.EQ.2) IOBL = 0
  165. SEGINI, MELVA2=MELVAL
  166. MTAIL2 = MELVA2
  167. IF (MELE.EQ.29) THEN
  168. * récupération de la section des inerties et de la torsion
  169. * section
  170. MELVA3 = MPTVA1.IVAL(4+IOBL)
  171. * inertie z
  172. MELVA4 = MPTVA1.IVAL(3+IOBL)
  173. * inertie y
  174. MELVA5 = MPTVA1.IVAL(2+IOBL)
  175. * torsion
  176. MELVA6 = MPTVA1.IVAL(1+IOBL)
  177. SEGACT MELVA3,MELVA4,MELVA5,MELVA6
  178. DO 2905 I=1,N1EL
  179. I3 = MIN(I,MELVA3.VELCHE(/2))
  180. I4 = MIN(I,MELVA4.VELCHE(/2))
  181. I5 = MIN(I,MELVA5.VELCHE(/2))
  182. I6 = MIN(I,MELVA6.VELCHE(/2))
  183. * remaarque les valeurs ne sont pas moyennées sur l'élément
  184. DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I)
  185. DUM1 = DUM1 * MELVA3.VELCHE(1,I3)
  186. DUM1=DUM1/MAX(MELVA4.VELCHE(1,I4),MELVA5.VELCHE(1,I5),
  187. & MELVA6.VELCHE(1,I6))
  188. DUM1=SQRT(DUM1) * 0.10132D0
  189. MELVA2.VELCHE(1,I)=MELVA2.VELCHE(1,I) * DUM1
  190. 2905 CONTINUE
  191. SEGDES MELVA3,MELVA4,MELVA5,MELVA6
  192. ELSEIF (MELE.EQ.84) THEN
  193. * récupération de la section des inerties et de la torsion
  194. * section
  195. MELVA3 = MPTVA1.IVAL(4+IOBL)
  196. * inertie z
  197. MELVA4 = MPTVA1.IVAL(3+IOBL)
  198. * inertie y
  199. MELVA5 = MPTVA1.IVAL(2+IOBL)
  200. * torsion
  201. MELVA6 = MPTVA1.IVAL(1+IOBL)
  202. SEGACT MELVA3,MELVA4,MELVA5,MELVA6
  203. DO 8405 I=1,N1EL
  204. I3 = MIN(I,MELVA3.VELCHE(/2))
  205. I4 = MIN(I,MELVA4.VELCHE(/2))
  206. I5 = MIN(I,MELVA5.VELCHE(/2))
  207. I6 = MIN(I,MELVA6.VELCHE(/2))
  208. * remarque les valeurs ne sont pas moyennées sur l'élément
  209. DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I)
  210. * 0.9 = 2**1.5 / XPI 0.7 coef de securite
  211. DUM2 = MELVA2.VELCHE(1,I)*0.9D0*0.7D0
  212. DUM1 = DUM1 * MELVA3.VELCHE(1,I3)
  213. DUM1=DUM1/MAX(MELVA4.VELCHE(1,I4),MELVA5.VELCHE(1,I5),
  214. & MELVA6.VELCHE(1,I6))
  215. DUM1=SQRT(DUM1) * 0.10132D0
  216. * correction pour tenir compte du cisaillement Blevins page 175 eq 8.30
  217. MELVA2.VELCHE(1,I)=MELVA2.VELCHE(1,I) * DUM1 + DUM2
  218. 8405 CONTINUE
  219. SEGDES MELVA3,MELVA4,MELVA5,MELVA6
  220. ELSEIF(MELE.EQ.42) THEN
  221. * cas du tuyau
  222. * récupération des caractéristiques
  223. * epai
  224. MELVA3 = MPTVA1.IVAL(1+IOBL)
  225. * rayo
  226. MELVA4 = MPTVA1.IVAL(2+IOBL)
  227. * pour les autres on ne ient pas compte des modification
  228. * qui assouplissent le tuyau donc omega max diminue
  229. SEGACT MELVA3,MELVA4
  230. DO 4205 I=1,N1EL
  231. I3 = MIN(I,MELVA3.VELCHE(/2))
  232. I4 = MIN(I,MELVA4.VELCHE(/2))
  233. REXT = MELVA4.VELCHE(1,I4)
  234. RINT = REXT - MELVA3.VELCHE(1,I3)
  235. * remarque les valeurs ne sont pas moyennées sur l'élément
  236. DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I)
  237. DUM1 = DUM1 * 2.D0 / (REXT*REXT+RINT*RINT)
  238. DUM1=SQRT(DUM1) * 0.10132D0
  239. MELVA2.VELCHE(1,I)=MELVA2.VELCHE(1,I) * DUM1
  240. 4205 CONTINUE
  241. SEGDES MELVA3,MELVA4
  242. ENDIF
  243. IF (ICAS.EQ.5) THEN
  244. MELV1 = MTAIL1
  245. MELV2 = MTAIL2
  246. SEGDES MELVAL,MELVA2
  247. RETURN
  248. ENDIF
  249. ELSE IF (ICAS.EQ.3) THEN
  250. * recuperation du champ
  251. * SEGACT MPTVA2
  252. MELVA1 = MPTVA2.IVAL(1)
  253. MELVA2 = MPTVA2.IVAL(2)
  254. SEGACT MELVA1,MELVA2
  255. MTAIL1 = MELVA1
  256. MTAIL2 = MELVA2
  257. SEGDES MPTVA2
  258. ENDIF
  259. * ================ pas de temps cfl
  260. IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN
  261. * recuperation de la vitesse du son
  262. * et du paramètre de taille
  263. MELVA1 = MCSON
  264. MELVA2 = MTAIL1
  265. MELVA3 = MTAIL2
  266. * creation du melval résultat
  267. N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2),
  268. & MELVA3.VELCHE(/2))
  269. N1PTEL = 1
  270. N2EL = 0
  271. N2PTEL = 0
  272. SEGINI MELVAL
  273. *
  274. DO 2904 I=1,N1EL
  275. I1 = MIN(I,MELVA1.VELCHE(/2))
  276. I2 = MIN(I,MELVA2.VELCHE(/2))
  277. I3 = MIN(I,MELVA3.VELCHE(/2))
  278. DUM1=MIN(MELVA2.VELCHE(1,I2),MELVA3.VELCHE(1,I3))
  279. VELCHE(1,I)= DUM1 / MELVA1.VELCHE(1,I1)
  280. * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i)
  281. 2904 CONTINUE
  282. MELV1 = MELVAL
  283. MELV2 = 0
  284. SEGDES MELVAL
  285. IF (ICAS.EQ.1) THEN
  286. SEGSUP MELVA1,MELVA2,MELVA3
  287. ELSE IF (ICAS.EQ.2) THEN
  288. SEGSUP MELVA2,MELVA3
  289. SEGDES MELVA1
  290. ELSE
  291. * icas=3
  292. SEGSUP MELVA1
  293. SEGDES MELVA2,MELVA3
  294. ENDIF
  295. RETURN
  296. ENDIF
  297. *
  298. *
  299. * cas du coq2
  300. *
  301. 44 CONTINUE
  302. * ================ calcul de la vitesse du son
  303. * le resultat est stocké dans melval avec n1ptel =1
  304. IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN
  305. * recherche des paramètre matériau
  306. * module d'young
  307. MELVA3 = MPTVA1.IVAL(1)
  308. * densite
  309. MELVA4 = MPTVA1.IVAL(3)
  310. * nu
  311. MELVA5 = MPTVA1.IVAL(2)
  312. * SEGDES MPTVA1
  313. SEGACT MELVA3,MELVA4,MELVA5
  314. *
  315. N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2),
  316. & MELVA5.VELCHE(/2))
  317. N1PTEL = 1
  318. N2PTEL = 0
  319. N2EL = 0
  320. SEGINI MELVAL
  321. MCSON = MELVAL
  322. * boucle sur les éléments pour calculer la vitesse du son
  323. DO 4404 I=1,N1EL
  324. * on prend les valeurs moyennes sur les éléments
  325. YOU1 = 0.D0
  326. I3 = MIN(I,MELVA3.VELCHE(/2))
  327. DO 4401 J=1,MELVA3.VELCHE(/1)
  328. YOU1 = YOU1 + MELVA3.VELCHE(J,I3)
  329. 4401 CONTINUE
  330. YOU1 = YOU1 / MELVA3.VELCHE(/1)
  331. *
  332. RO1 = 0.D0
  333. I4 = MIN(I,MELVA4.VELCHE(/2))
  334. DO 4402 J=1,MELVA4.VELCHE(/1)
  335. RO1 = RO1 + MELVA4.VELCHE(J,I4)
  336. 4402 CONTINUE
  337. RO1 = RO1 / MELVA4.VELCHE(/1)
  338. *
  339. DNU1 = 0.D0
  340. I5 = MIN(I,MELVA5.VELCHE(/2))
  341. DO 4403 J=1,MELVA5.VELCHE(/1)
  342. DNU1 = DNU1 + MELVA5.VELCHE(J,I5)
  343. 4403 CONTINUE
  344. DNU1 = DNU1 / MELVA5.VELCHE(/1)
  345.  
  346. IF (RO1.EQ.0.D0) THEN
  347. SEGDES MELVA4,MELVA3,MELVA5
  348. SEGSUP MELVAL
  349. CALL ERREUR(855)
  350. RETURN
  351. ENDIF
  352. *
  353. IF (YOU1.EQ.0.D0) THEN
  354. SEGDES MELVA4,MELVA3,MELVA5
  355. SEGSUP MELVAL
  356. CALL ERREUR(856)
  357. RETURN
  358. ENDIF
  359. IF (DNU1.GE.1.D0) THEN
  360. SEGDES MELVA4,MELVA3,MELVA5
  361. SEGSUP MELVAL
  362. CALL ERREUR(36)
  363. RETURN
  364. ENDIF
  365. *
  366. VELCHE(1,I) = SQRT(YOU1/(RO1*(1-DNU1*DNU1)))
  367. * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i)
  368. 4404 CONTINUE
  369. SEGDES MELVA4,MELVA3,MELVA5
  370. IF (ICAS.EQ.4) THEN
  371. * cas ou seule la vitesse du son est demandée
  372. MELVA2 = 0
  373. MELV1 = MELVAL
  374. SEGDES MELVAL
  375. RETURN
  376. ENDIF
  377. ELSE IF (ICAS.EQ.2) THEN
  378. * recuperation du champ
  379. MELVA1 = MPTVA2.IVAL(1)
  380. SEGACT MELVA1
  381. MCSON = MELVA1
  382. SEGDES MPTVA2
  383. ENDIF
  384. * ================ paramètre geometrique
  385. * stocké dans un melval mtaille
  386. IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN
  387. MELEME = IPMAIL
  388. SEGACT MELEME
  389. N1EL = NUM(/2)
  390. N1PTEL = 1
  391. N2PTEL = 0
  392. N2EL = 0
  393. SEGINI MELVAL
  394. MTAIL1 = MELVAL
  395. CALL CFLTAI(MTAIL1,IPMAIL,MELE)
  396. SEGDES MELEME
  397. * second paramètre geométrique
  398. IOBL = 3
  399. IF (ICAS.EQ.5 .OR. ICAS.EQ.2) IOBL = 0
  400. SEGINI, MELVA2=MELVAL
  401. MTAIL2 = MELVA2
  402. * récupération de l' epaisseur
  403. MELVA3 = MPTVA1.IVAL(1+IOBL)
  404. SEGACT MELVA3
  405. DO 4406 I=1,N1EL
  406. I3 = MIN(I,MELVA3.VELCHE(/2))
  407. * remarque les valeurs ne sont pas moyennées sur l'élément
  408. DUM1 = MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I)
  409. DUM1 = DUM1 / (3.23D0 * MELVA3.VELCHE(1,I3))
  410. MELVA2.VELCHE(1,I)=DUM1
  411. 4406 CONTINUE
  412. SEGDES MELVA3
  413. SEGDES MPTVA1
  414. *
  415. IF (ICAS.EQ.5) THEN
  416. MELV1 = MTAIL1
  417. MELV2 = MTAIL2
  418. SEGDES MELVAL
  419. RETURN
  420. ENDIF
  421. ELSE IF (ICAS.EQ.3) THEN
  422. * recuperation du champ
  423. MELVA1 = MPTVA2.IVAL(1)
  424. MELVA2 = MPTVA2.IVAL(2)
  425. SEGACT MELVA1,MELVA2
  426. MTAIL1 = MELVA1
  427. MTAIL2 = MELVA2
  428. SEGDES MPTVA2
  429. ENDIF
  430. * ================ pas de temps cfl
  431. IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN
  432. * recuperation de la vitesse du son
  433. * et du paramètre de taille
  434. MELVA1 = MCSON
  435. MELVA2 = MTAIL1
  436. MELVA3 = MTAIL2
  437. * creation du melval résultat
  438. N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2))
  439. N1PTEL = 1
  440. N2EL = 0
  441. N2PTEL = 0
  442. SEGINI MELVAL
  443. *
  444. DO 4405 I=1,N1EL
  445. I1 = MIN(I,MELVA1.VELCHE(/2))
  446. I2 = MIN(I,MELVA2.VELCHE(/2))
  447. I3 = MIN(I,MELVA3.VELCHE(/2))
  448. DUM1=MIN(MELVA2.VELCHE(1,I2),MELVA3.VELCHE(1,I3))
  449. VELCHE(1,I)=DUM1/MELVA1.VELCHE(1,I1)
  450. * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i)
  451. 4405 CONTINUE
  452. MELV1 = MELVAL
  453. MELV2 = 0
  454. SEGDES MELVAL
  455. IF (ICAS.EQ.1) THEN
  456. SEGSUP MELVA1,MELVA2,MELVA3
  457. ELSE IF (ICAS.EQ.2) THEN
  458. SEGSUP MELVA2,MELVA3
  459. SEGDES MELVA1
  460. ELSE
  461. SEGSUP MELVA1
  462. SEGDES MELVA2,MELVA3
  463. ENDIF
  464. RETURN
  465. ENDIF
  466. *
  467. *
  468. *
  469. * cas de la barre
  470. *
  471. 46 CONTINUE
  472. * ================ calcul de la vitesse du son
  473. * le resultat est stocké dans melval avec n1ptel =1
  474. IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN
  475. * recherche des paramètre matériau
  476. * module d'young
  477. MELVA3 = MPTVA1.IVAL(1)
  478. * densite
  479. MELVA4 = MPTVA1.IVAL(3)
  480. SEGDES MPTVA1
  481. SEGACT MELVA3,MELVA4
  482. *
  483. N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2))
  484. N1PTEL = 1
  485. N2PTEL = 0
  486. N2EL = 0
  487. SEGINI MELVAL
  488. MCSON = MELVAL
  489. * boucle sur les éléments pour calculer la vitesse du son
  490. DO 4603 I=1,N1EL
  491. * on prend les valeurs moyennes sur les éléments
  492. YOU1 = 0.D0
  493. I3 = MIN(I,MELVA3.VELCHE(/2))
  494. DO 4601 J=1,MELVA3.VELCHE(/1)
  495. YOU1 = YOU1 + MELVA3.VELCHE(J,I3)
  496. 4601 CONTINUE
  497. YOU1 = YOU1 / MELVA3.VELCHE(/1)
  498. *
  499. RO1 = 0.D0
  500. I4 = MIN(I,MELVA4.VELCHE(/2))
  501. DO 4602 J=1,MELVA4.VELCHE(/1)
  502. RO1 = RO1 + MELVA4.VELCHE(J,I4)
  503. 4602 CONTINUE
  504. RO1 = RO1 / MELVA4.VELCHE(/1)
  505. IF (RO1.EQ.0.D0) THEN
  506. SEGDES MELVA4,MELVA3
  507. SEGSUP MELVAL
  508. CALL ERREUR(855)
  509. RETURN
  510. ENDIF
  511. *
  512. IF (YOU1.EQ.0.D0) THEN
  513. SEGDES MELVA4,MELVA3
  514. SEGSUP MELVAL
  515. CALL ERREUR(856)
  516. RETURN
  517. ENDIF
  518. VELCHE(1,I) = SQRT(YOU1/RO1)
  519. * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i)
  520. 4603 CONTINUE
  521. SEGDES MELVA4,MELVA3
  522. IF (ICAS.EQ.4) THEN
  523. * cas ou seule la vitesse du son est demandée
  524. MELVA2 = 0
  525. MELV1 = MELVAL
  526. SEGDES MELVAL
  527. RETURN
  528. ENDIF
  529. ELSE IF (ICAS.EQ.2) THEN
  530. * recuperation du champ
  531. MELVA1 = MPTVA2.IVAL(1)
  532. SEGACT MELVA1
  533. MCSON = MELVA1
  534. SEGDES MPTVA2
  535. ENDIF
  536. * ================ paramètre geometrique
  537. * stocké dans un melval mtaille
  538. IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN
  539. MELEME = IPMAIL
  540. SEGACT MELEME
  541. N1EL = NUM(/2)
  542. N1PTEL = 1
  543. N2PTEL = 0
  544. N2EL = 0
  545. SEGINI MELVAL
  546. MTAIL1 = MELVAL
  547. MTAIL2 = 0
  548. CALL CFLTAI(MTAIL1,IPMAIL,MELE)
  549. SEGDES MELEME
  550. IF (ICAS.EQ.5) THEN
  551. MELV1 = MTAIL1
  552. MELV2 = 0
  553. SEGDES MELVAL
  554. RETURN
  555. ENDIF
  556. ELSE IF (ICAS.EQ.3) THEN
  557. * recuperation du champ
  558. MELVA1 = MPTVA2.IVAL(1)
  559. SEGACT MELVA1
  560. MTAIL1 = MELVA1
  561. SEGDES MPTVA2
  562. ENDIF
  563. * ================ pas de temps cfl
  564. IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN
  565. * recuperation de la vitesse du son
  566. * et du paramètre de taille
  567. MELVA1 = MCSON
  568. MELVA2 = MTAIL1
  569. * creation du melval résultat
  570. N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2))
  571. N1PTEL = 1
  572. N2EL = 0
  573. N2PTEL = 0
  574. SEGINI MELVAL
  575. *
  576. DO 4604 I=1,N1EL
  577. I1 = MIN(I,MELVA1.VELCHE(/2))
  578. I2 = MIN(I,MELVA2.VELCHE(/2))
  579. VELCHE(1,I)=MELVA2.VELCHE(1,I2)/MELVA1.VELCHE(1,I1)
  580. * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i)
  581. 4604 CONTINUE
  582. MELV1 = MELVAL
  583. MELV2 = 0
  584. SEGDES MELVAL
  585. IF (ICAS.EQ.1) THEN
  586. SEGSUP MELVA1,MELVA2
  587. ELSE IF (ICAS.EQ.2) THEN
  588. SEGSUP MELVA2
  589. SEGDES MELVA1
  590. ELSE
  591. SEGSUP MELVA1
  592. SEGDES MELVA2
  593. ENDIF
  594. RETURN
  595. ENDIF
  596. *
  597. *
  598. * cas du coq3 et du dkt (meme avec excentrement)
  599. * Rq: dst aussi dans les cas simples
  600. *
  601. 27 CONTINUE
  602. * ================ calcul de la vitesse du son
  603. * le resultat est stocké dans melval avec n1ptel =1
  604. IF (ICAS.EQ.1.OR.ICAS.EQ.3.OR.ICAS.EQ.4) THEN
  605. * recherche des paramètre matériau
  606. * module d'young
  607. MELVA3 = MPTVA1.IVAL(1)
  608. * densite
  609. MELVA4 = MPTVA1.IVAL(3)
  610. * nu
  611. MELVA5 = MPTVA1.IVAL(2)
  612. * SEGDES MPTVA1
  613. SEGACT MELVA3,MELVA4,MELVA5
  614. *
  615. N1EL = MIN(MELVA3.VELCHE(/2),MELVA4.VELCHE(/2),
  616. & MELVA5.VELCHE(/2))
  617. N1PTEL = 1
  618. N2PTEL = 0
  619. N2EL = 0
  620. SEGINI MELVAL
  621. MCSON = MELVAL
  622. * boucle sur les éléments pour calculer la vitesse du son
  623. DO 2704 I=1,N1EL
  624. * on prend les valeurs moyennes sur les éléments
  625. YOU1 = 0.D0
  626. I3 = MIN(I,MELVA3.VELCHE(/2))
  627. DO 2701 J=1,MELVA3.VELCHE(/1)
  628. YOU1 = YOU1 + MELVA3.VELCHE(J,I3)
  629. 2701 CONTINUE
  630. YOU1 = YOU1 / MELVA3.VELCHE(/1)
  631. *
  632. RO1 = 0.D0
  633. I4 = MIN(I,MELVA4.VELCHE(/2))
  634. DO 2702 J=1,MELVA4.VELCHE(/1)
  635. RO1 = RO1 + MELVA4.VELCHE(J,I4)
  636. 2702 CONTINUE
  637. RO1 = RO1 / MELVA4.VELCHE(/1)
  638. *
  639. DNU1 = 0.D0
  640. I5 = MIN(I,MELVA5.VELCHE(/2))
  641. DO 2703 J=1,MELVA5.VELCHE(/1)
  642. DNU1 = DNU1 + MELVA5.VELCHE(J,I5)
  643. 2703 CONTINUE
  644. DNU1 = DNU1 / MELVA5.VELCHE(/1)
  645.  
  646. IF (RO1.EQ.0.D0) THEN
  647. SEGDES MELVA4,MELVA3,MELVA5
  648. SEGSUP MELVAL
  649. CALL ERREUR(855)
  650. RETURN
  651. ENDIF
  652. *
  653. IF (YOU1.EQ.0.D0) THEN
  654. SEGDES MELVA4,MELVA3,MELVA5
  655. SEGSUP MELVAL
  656. CALL ERREUR(856)
  657. RETURN
  658. ENDIF
  659. IF (DNU1.GE.1.D0) THEN
  660. SEGDES MELVA4,MELVA3,MELVA5
  661. SEGSUP MELVAL
  662. CALL ERREUR(36)
  663. RETURN
  664. ENDIF
  665. *
  666. VELCHE(1,I) = SQRT(YOU1/(RO1*(1-DNU1*DNU1)))
  667. * write(6,*) 'Element', i , 'Cson' , VELCHE(1,i)
  668. 2704 CONTINUE
  669. SEGDES MELVA4,MELVA3,MELVA5
  670. IF (ICAS.EQ.4) THEN
  671. * cas ou seule la vitesse du son est demandée
  672. MELVA2 = 0
  673. MELV1 = MELVAL
  674. SEGDES MELVAL
  675. RETURN
  676. ENDIF
  677. ELSE IF (ICAS.EQ.2) THEN
  678. * recuperation du champ
  679. MELVA1 = MPTVA2.IVAL(1)
  680. SEGACT MELVA1
  681. MCSON = MELVA1
  682. SEGDES MPTVA2
  683. ENDIF
  684. * ================ paramètre geometrique
  685. * stocké dans un melval mtaille
  686. IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.5) THEN
  687. MELEME = IPMAIL
  688. SEGACT MELEME
  689. N1EL = NUM(/2)
  690. N1PTEL = 1
  691. N2PTEL = 0
  692. N2EL = 0
  693. SEGINI MELVAL
  694. MTAIL1 = MELVAL
  695. CALL CFLTAI(MTAIL1,IPMAIL,MELE)
  696. SEGDES MELEME
  697. * second paramètre geométrique
  698. IOBL = 3
  699. IF (ICAS.EQ.5 .OR. ICAS.EQ.2) IOBL = 0
  700. SEGINI, MELVA2=MELVAL
  701. MTAIL2 = MELVA2
  702. * récupération de l' epaisseur
  703. MELVA3 = MPTVA1.IVAL(1+IOBL)
  704. SEGACT MELVA3
  705. SEGACT MELEME
  706. DO 2706 I=1,N1EL
  707. I3 = MIN(I,MELVA3.VELCHE(/2))
  708. * remarque les valeurs ne sont pas moyennées sur l'élément
  709. V12L = 0.D0
  710. V13L = 0.D0
  711. V23L = 0.D0
  712. DO 2761 J=1,IDIM
  713. V13(J)=XCOOR((IDIM+1)*(NUM(3,I3)-1)+J)-
  714. & XCOOR((IDIM+1)*(NUM(1,I3)-1)+J)
  715. V12(J)=XCOOR((IDIM+1)*(NUM(2,I)-1)+J)-
  716. & XCOOR((IDIM+1)*(NUM(1,I)-1)+J)
  717. V23(J)=XCOOR((IDIM+1)*(NUM(3,I)-1)+J)-
  718. & XCOOR((IDIM+1)*(NUM(2,I)-1)+J)
  719. V12L = V12L + V12(J)*V12(J)
  720. V13L = V13L + V13(J)*V13(J)
  721. V23L = V23L + V23(J)*V23(J)
  722. 2761 CONTINUE
  723. D2min = 0.D0
  724. D2max = 0.D0
  725. Dmax = 0.D0
  726. D2min = MIN(V12L,V13L,V23L)
  727. D2max = MAX(V12L,V13L,V23L)
  728. Dmax = SQRT(D2max)
  729. * write(6,*) 'D2min', D2min , 'D2max' , D2max, 'Dmax', Dmax
  730. * write(6,*) 'Element', i , 'L' , MELVA2.VELCHE(1,i)
  731. DUM1 = SQRT(D2min-MELVA2.VELCHE(1,I)*MELVA2.VELCHE(1,I))
  732. DUM1 = DUM1 / Dmax - 0.5D0
  733. DUM1 = 5.3706D0 * DUM1 * DUM1
  734. DUM1 = 3.9402D0 - 0.8687D0*MELVA2.VELCHE(1,I)/Dmax - DUM1
  735. DUM1 = Dmax / (XPI * MELVA3.VELCHE(1,I3) * DUM1)
  736. DUM1 = DUM1 * MELVA2.VELCHE(1,I)
  737. IF (MELE.EQ.28 .AND. MPTVA1.IVAL(2+IOBL).NE.0) THEN
  738. * récupération de l'excentricité
  739. MELVA4 = MPTVA1.IVAL(2+IOBL)
  740. SEGACT MELVA4
  741. I4 = MIN(I,MELVA4.VELCHE(/2))
  742. DUM1=DUM1/SQRT(1D0+12D0*MELVA4.VELCHE(1,I4)*MELVA4.VELCHE(1,I4)
  743. & /(MELVA3.VELCHE(1,I3)*MELVA3.VELCHE(1,I3)))
  744. * On applique un coefficient empirique de 1.5 qui tient compte des
  745. * différences de discrétisation par rapport au modèle non excentré
  746. DUM1 = DUM1 * 1.5D0
  747. SEGDES MELVA4
  748. ENDIF
  749. MELVA2.VELCHE(1,I)=DUM1
  750. 2706 CONTINUE
  751. SEGDES MELEME
  752. SEGDES MELVA3
  753. SEGDES MPTVA1
  754. *
  755. IF (ICAS.EQ.5) THEN
  756. MELV1 = MTAIL1
  757. MELV2 = MTAIL2
  758. SEGDES MELVAL
  759. RETURN
  760. ENDIF
  761. ELSE IF (ICAS.EQ.3) THEN
  762. * recuperation du champ
  763. MELVA1 = MPTVA2.IVAL(1)
  764. MELVA2 = MPTVA2.IVAL(2)
  765. SEGACT MELVA1,MELVA2
  766. MTAIL1 = MELVA1
  767. MTAIL2 = MELVA2
  768. SEGDES MPTVA2
  769. ENDIF
  770. * ================ pas de temps cfl
  771. IF (ICAS.EQ.1.OR.ICAS.EQ.2.OR.ICAS.EQ.3) THEN
  772. * recuperation de la vitesse du son
  773. * et du paramètre de taille
  774. MELVA1 = MCSON
  775. MELVA2 = MTAIL1
  776. MELVA3 = MTAIL2
  777. * creation du melval résultat
  778. N1EL = MAX(MELVA1.VELCHE(/2),MELVA2.VELCHE(/2))
  779. N1PTEL = 1
  780. N2EL = 0
  781. N2PTEL = 0
  782. SEGINI MELVAL
  783. *
  784. DO 2705 I=1,N1EL
  785. I1 = MIN(I,MELVA1.VELCHE(/2))
  786. I2 = MIN(I,MELVA2.VELCHE(/2))
  787. I3 = MIN(I,MELVA3.VELCHE(/2))
  788. DUM1=MIN(MELVA2.VELCHE(1,I2),MELVA3.VELCHE(1,I3))
  789. VELCHE(1,I)=DUM1/MELVA1.VELCHE(1,I1)
  790. * write(6,*) 'Element', i , 'Dtcfl' , VELCHE(1,i)
  791. 2705 CONTINUE
  792. MELV1 = MELVAL
  793. MELV2 = 0
  794. SEGDES MELVAL
  795. IF (ICAS.EQ.1) THEN
  796. SEGSUP MELVA1,MELVA2,MELVA3
  797. ELSE IF (ICAS.EQ.2) THEN
  798. SEGSUP MELVA2,MELVA3
  799. SEGDES MELVA1
  800. ELSE
  801. SEGSUP MELVA1
  802. SEGDES MELVA2,MELVA3
  803. ENDIF
  804. RETURN
  805. ENDIF
  806. *
  807. 99 CONTINUE
  808. MOTERR(1:4)=NOMTP(MELE)
  809. MOTERR(9:12)='CFL5'
  810. CALL ERREUR(86)
  811. *
  812. RETURN
  813. END
  814.  
  815.  
  816.  
  817.  
  818.  
  819.  
  820.  
  821.  
  822.  
  823.  
  824.  

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