Télécharger devfb1.eso

Retour à la liste

Numérotation des lignes :

  1. C DEVFB1 SOURCE BP208322 19/02/25 21:15:31 10120
  2. SUBROUTINE DEVFB1(ITYP,FTOTB,XPTB,IPALB,IPLIB,XPALB,XVALB,NLIAB,
  3. & NPLB,IND,IND1,INDM1,NPA,NPAM1,IND2,PDT,PDTS2,
  4. & FEXPSM,NPC1,XABSCI,XORDON,NIP,I,iannul)
  5. IMPLICIT INTEGER(I-N)
  6. IMPLICIT REAL*8(A-H,O-Z)
  7. *--------------------------------------------------------------------*
  8. * *
  9. * Operateur DYNE : algorithme de Fu - de Vogelaere *
  10. * ________________________________________________ *
  11. * *
  12. * Calcul des forces de choc sur base B pour les liaisons de *
  13. * type POINT_PLAN. *
  14. * *
  15. * Parametres: *
  16. * *
  17. * e ITYP type de la liaison. *
  18. * es FTOTB Forces exterieures totalisees sur la base B. *
  19. * e XPTB Tableau des deplacements des points *
  20. * e IPALB Renseigne sur la liaison. *
  21. * e IPLIB Tableau contenant les numeros "DYNE" de la liaison. *
  22. * e XPALB Tableau contenant les parametres de la liaison. *
  23. * es XVALB Tableau contenant les variables internes de liaisons. *
  24. * e NLIAB Nombre de liaisons sur la base B. *
  25. * e NPLB Nombre total de points intervenant dans les liaisons. *
  26. * e IND Indice du pas. *
  27. * e IND1 Indice du pas (ou demi-pas si De Vogelaere) precedent. *
  28. * e I numero de la liaison. *
  29. * *
  30. * Auteur, date de creation: *
  31. * *
  32. * Lionel VIVAN : le 22 Septembre 1989 : Creation *
  33. * Bertrand BEAUFILS : le 31 Mai 1990 : Ajout frottement sec*
  34. * *
  35. *--------------------------------------------------------------------*
  36. *
  37. INTEGER IPALB(NLIAB,*),IPLIB(NLIAB,*)
  38. REAL*8 XPALB(NLIAB,*),XPTB(NPLB,4,*),FTOTB(NPLB,*)
  39. REAL*8 XVALB(NLIAB,4,*),FEXPSM(NPLB,NPC1,2,*)
  40. REAL*8 XPTP2(3),XPTPP2(3),XFNT(3),XPSMP(3),XPSMPM(3)
  41. REAL*8 XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  42. *
  43. *--------------------------------------------------------------------*
  44. * --- choc elementaire POINT_PLAN
  45. * avec ou sans amortissement
  46. *
  47. IF (ITYP.EQ.1 .OR. ITYP.EQ.100
  48. & .OR. ITYP.EQ.101 .OR. ITYP.EQ.102) THEN
  49. NPOI = IPLIB(I,1)
  50. IDIM = IPALB(I,3)
  51. IPERM = IPALB(I,4)
  52. XRAID = XPALB(I,1)
  53. XJEU = XPALB(I,2)
  54. XAMO = XPALB(I,3)
  55. IF (ITYP.EQ.1 .OR. ITYP.EQ.102) THEN
  56. ID1 = 3
  57. ELSE
  58. ID1 = 4
  59. XSEUIL = XPALB(I,ID1)
  60. XDPLAS = XPALB(I,(ID1+IDIM+1))
  61. ENDIF
  62. *
  63. c XDEP : deplacement normal = x(t_{ind})*n
  64. c XDEPM1 : deplacement normal = x(t_{ind2})*n
  65. c XVIT : vitesse normale = v(t_{ind})*n
  66. XDEP = 0.D0
  67. DO 10 ID = 1,IDIM
  68. IDD1 = 3 + ID
  69. XDE2 = XPTB(NPOI,IND,ID)
  70. XVALB(I,IND,IDD1) = XDE2
  71. XDE2 = XDE2 + FEXPSM(NPOI,NPA,IND1,ID)
  72. XDEP = XDEP + XDE2 * XPALB(I,ID1+ID)
  73. 10 CONTINUE
  74. XDEPM1 = 0.D0
  75. DO 12 ID = 1,IDIM
  76. XDM2 = XPTB(NPOI,IND2,ID) + FEXPSM(NPOI,NPAM1,INDM1,ID)
  77. XDEPM1 = XDEPM1 + XDM2 * XPALB(I,ID1+ID)
  78. 12 CONTINUE
  79. XVIT = (XDEP - XDEPM1) / PDTS2
  80. XVALB(I,IND,3) = XVIT
  81.  
  82. IF (ITYP.EQ.1) THEN
  83. CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,XFL,IPERM,iannul)
  84. ELSE IF (ITYP.EQ.102) THEN
  85. CALL DYCHAM2(XDEP,XVIT,XJEU,XAMO,XABSCI,XORDON,NIP,
  86. & NLIAB,I,XFL,IPERM,iannul)
  87. ELSE IF (ITYP.EQ.101) THEN
  88. CALL DYCHPL2(XDEP,XVIT,XDPLAS,XRAID,XJEU,XAMO,XSEUIL,
  89. & XABSCI,XORDON,NIP,NLIAB,I,XFL,IPERM,iannul)
  90. ELSE IF (ITYP.EQ.100) THEN
  91. CALL DYCHPL(XDEP,XVIT,XDPLAS,XRAID,XJEU,XAMO,XSEUIL,
  92. & XFL,IPERM,iannul)
  93. ENDIF
  94.  
  95. IF (ITYP.EQ.100 .OR. ITYP.EQ.101) THEN
  96. XPALB(I,(ID1+IDIM+1)) = XDPLAS
  97. XVALB(I,IND,13) = XDPLAS
  98. ENDIF
  99. XVALB(I,IND,1) = XFL
  100.  
  101. DO 14 ID = 1,IDIM
  102. FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFL * XPALB(I,ID1+ID)
  103. 14 CONTINUE
  104.  
  105. *
  106. *--------------------------------------------------------------------*
  107. * --- choc elementaire POINT_PLAN_FROTTEMENT
  108. * avec ou sans amortissement
  109. *
  110. ELSE IF (ITYP.EQ.3 .OR. ITYP.EQ.103) THEN
  111. NPOI = IPLIB(I,1)
  112. IGP = IPALB(I,2)
  113. IDIM = IPALB(I,3)
  114. ID1 = 7
  115. ID2 = ID1 + IDIM
  116. ID3 = ID1 + 2*IDIM
  117. ID4 = ID1 + 3*IDIM
  118. ID5 = ID1 + 4*IDIM
  119. ID6 = ID1 + 5*IDIM
  120. ID7 = ID1 + 6*IDIM
  121. * Si glissement au pas precedent, reactualisation de la position
  122. * origine d'adherence
  123. IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN
  124. DO 30 ID=1,IDIM
  125. XPALB(I,ID5+ID) = XPTB(NPOI,IND,ID) +
  126. & FEXPSM(NPOI,NPA,IND1,ID)
  127. 30 CONTINUE
  128. ENDIF
  129. * Calcul de l'enfoncement et de la vitesse normale
  130. c XDEP : deplacement normal = x(t_{ind})*n
  131. c XVITN : vitesse normale = v(t_{ind})*n = (dx*n)/dt
  132. XDEP = 0.D0
  133. PSN = 0.D0
  134. PSN0 = 0.D0
  135. DO 32 ID = 1,IDIM
  136. IDD1 = 3 + ID
  137. XDE2 = XPTB(NPOI,IND,ID)
  138. XDM2 = XPTB(NPOI,IND2,ID)
  139. XVALB(I,IND,IDD1) = XDE2
  140. XDE2 = XDE2 + FEXPSM(NPOI,NPA,IND1,ID)
  141. XDM2 = XDM2 + FEXPSM(NPOI,NPAM1,INDM1,ID)
  142. XPALB(I,ID2+ID) = XDE2 - XDM2
  143. XPALB(I,ID3+ID) = XDE2 - XPALB(I,ID5+ID)
  144. XDEP = XDEP + XDE2 * XPALB(I,ID1+ID)
  145. PSN = PSN + XPALB(I,ID2+ID) * XPALB(I,ID1+ID)
  146. PSN0 = PSN0 + XPALB(I,ID3+ID) * XPALB(I,ID1+ID)
  147. 32 CONTINUE
  148. * Projette sur le plan tangent la vitesse et la variation de deplacement
  149. * par rapport a la position origine d'adherence sur le plan tangent
  150. * XPALB(I,ID2+ID) = (dx - (dx*n)*n)/dt
  151. * avec dx = x(t_ind)-x(t_ind2)
  152. * XPALB(I,ID3+ID) = (dx_adhe - (dx_adhe*n)*n)/dt
  153. * avec dx_adhe = x(t_ind)-x(t_adherence)
  154. DO 34 ID = 1,IDIM
  155. XPALB(I,ID2+ID) = (XPALB(I,ID2+ID) - PSN * XPALB(I,ID1+ID))
  156. & / PDTS2
  157. XPALB(I,ID3+ID) = XPALB(I,ID3+ID) - PSN0 * XPALB(I,ID1+ID)
  158. 34 CONTINUE
  159. XVITN = PSN / PDTS2
  160. XVALB(I,IND,3) = XVITN
  161. IF (ITYP.EQ.3) THEN
  162. CALL DYCHA4(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,XFN,XFT,XPUS,
  163. & iannul)
  164. ELSE
  165. CALL DYCHA41(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,XFN,XFT,XPUS,
  166. & XABSCI,XORDON,NIP,iannul)
  167. ENDIF
  168. XVALB(I,IND,1) = XFN
  169. XVALB(I,IND,10) = ABS(XFT)
  170. XVALB(I,IND,12) = XPUS
  171. IGPM1 = IPALB(I,2)
  172. IPALB(I,2) = IGP
  173.  
  174. IF (IGP.EQ.1) THEN
  175. PS = 0.D0
  176. DO 20 ID = 1,IDIM
  177. PS = PS + (XPALB(I,ID2+ID)*XPALB(I,ID2+ID))
  178. 20 CONTINUE
  179. XVITT = SQRT(PS)
  180. ELSE
  181. XVITT = 0.D0
  182. ENDIF
  183. XVALB(I,IND,11) = XVITT
  184. * Si glissement, memorisation de la vitesse tangentielle et de la force
  185. * tangentielle
  186. IF (IGP.EQ.1.OR.IGP.EQ.-1) THEN
  187. DO 36 ID = 1,IDIM
  188. XPALB(I,ID6+ID) = XPALB(I,ID2+ID)
  189. XPALB(I,ID7+ID) = XPALB(I,ID4+ID)
  190. 36 CONTINUE
  191. ENDIF
  192. DO 38 ID = 1,IDIM
  193. FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFN * XPALB(I,ID1+ID)
  194. & + XPALB(I,ID4+ID)
  195. 38 CONTINUE
  196.  
  197.  
  198. *--------------------------------------------------------------------*
  199. * --- choc elementaire CERCLE_PLAN_FROTTEMENT
  200. * avec ou sans amortissement
  201. *
  202. ELSE IF (ITYP.EQ.5 .OR. ITYP.EQ.6) THEN
  203. NPOI = IPLIB(I,1)
  204. IGP = IPALB(I,2)
  205. IDIM = IPALB(I,3)
  206. IF (ITYP.EQ.5) THEN
  207. ID1 = 6
  208. ELSE
  209. ID1 = 7
  210. ENDIF
  211. ID2 = ID1 + IDIM
  212. ID3 = ID1 + 2*IDIM
  213. ID4 = ID1 + 3*IDIM
  214. ID5 = ID1 + 4*IDIM
  215. ID6 = ID1 + 5*IDIM
  216. ID7 = ID1 + 6*IDIM
  217. ID8 = ID1 + 1 + 7*IDIM
  218. * calcul du deplacement du point de contact au pas precedent
  219. * XMP = n * Rayon
  220. XMP1 = XPALB(I,ID1+1) * XPALB(I,ID8)
  221. XMP2 = XPALB(I,ID1+2) * XPALB(I,ID8)
  222. XMP3 = XPALB(I,ID1+3) * XPALB(I,ID8)
  223. * XPTPP2 = Ucontact(t) = Ucentre(t) + R(t) ∧ (n * Rayon)
  224. XPTPP2(1) = XPTB(NPOI,IND2,1) +
  225. & ( ( XPTB(NPOI,IND2,5) * XMP3 ) -
  226. & ( XPTB(NPOI,IND2,6) * XMP2 ) )
  227. XPTPP2(2) = XPTB(NPOI,IND2,2) +
  228. & ( ( XPTB(NPOI,IND2,6) * XMP1 ) -
  229. & ( XPTB(NPOI,IND2,4) * XMP3 ) )
  230. XPTPP2(3) = XPTB(NPOI,IND2,3) +
  231. & ( ( XPTB(NPOI,IND2,4) * XMP2 ) -
  232. & ( XPTB(NPOI,IND2,5) * XMP1 ) )
  233. XPSMPM(1) = FEXPSM(NPOI,NPAM1,INDM1,1) +
  234. & ( ( FEXPSM(NPOI,NPAM1,INDM1,5) * XMP3 ) -
  235. & ( FEXPSM(NPOI,NPAM1,INDM1,6) * XMP2 ) )
  236. XPSMPM(2) = FEXPSM(NPOI,NPAM1,INDM1,2) +
  237. & ( ( FEXPSM(NPOI,NPAM1,INDM1,6) * XMP1 ) -
  238. & ( FEXPSM(NPOI,NPAM1,INDM1,4) * XMP3 ) )
  239. XPSMPM(3) = FEXPSM(NPOI,NPAM1,INDM1,3) +
  240. & ( ( FEXPSM(NPOI,NPAM1,INDM1,4) * XMP2 ) -
  241. & ( FEXPSM(NPOI,NPAM1,INDM1,5) * XMP1 ) )
  242. * calcul du deplacement du point de contact au pas courant
  243. * XPTP2 = Ucontact(t+dt) = Ucentre(t+dt) + R(t+dt) ∧ (n * Rayon)
  244. XPTP2(1) = XPTB(NPOI,IND,1) +
  245. & ( ( XPTB(NPOI,IND,5) * XMP3 ) -
  246. & ( XPTB(NPOI,IND,6) * XMP2 ) )
  247. XPTP2(2) = XPTB(NPOI,IND,2) +
  248. & ( ( XPTB(NPOI,IND,6) * XMP1 ) -
  249. & ( XPTB(NPOI,IND,4) * XMP3 ) )
  250. XPTP2(3) = XPTB(NPOI,IND,3) +
  251. & ( ( XPTB(NPOI,IND,4) * XMP2 ) -
  252. & ( XPTB(NPOI,IND,5) * XMP1 ) )
  253. XPSMP(1) = FEXPSM(NPOI,NPA,IND1,1) +
  254. & ( ( FEXPSM(NPOI,NPA,IND1,5) * XMP3 ) -
  255. & ( FEXPSM(NPOI,NPA,IND1,6) * XMP2 ) )
  256. XPSMP(2) = FEXPSM(NPOI,NPA,IND1,2) +
  257. & ( ( FEXPSM(NPOI,NPA,IND1,6) * XMP1 ) -
  258. & ( FEXPSM(NPOI,NPA,IND1,4) * XMP3 ) )
  259. XPSMP(3) = FEXPSM(NPOI,NPA,IND1,3) +
  260. & ( ( FEXPSM(NPOI,NPA,IND1,4) * XMP2 ) -
  261. & ( FEXPSM(NPOI,NPA,IND1,5) * XMP1 ) )
  262. *
  263. * Si glissement au pas precedent, reactualisation de la position
  264. * origine d'adherence
  265. IF (IGP.EQ.1 .OR. IGP.EQ.-1) THEN
  266. DO 50 ID=1,IDIM
  267. XPALB(I,ID5+ID) = XPTP2(ID) + XPSMP(ID)
  268. 50 CONTINUE
  269. ENDIF
  270. * Calcul de l'enfoncement et de la vitesse normale
  271. XDEP = 0.D0
  272. PSN = 0.D0
  273. PSN0 = 0.D0
  274. DO 52 ID = 1,IDIM
  275. IDD1 = 3 + ID
  276. IDD2 = 6 + ID
  277. IDD3 = 15 + ID
  278. XDE2 = XPTP2(ID)
  279. XDM2 = XPTPP2(ID)
  280. XVALB(I,IND,IDD1) = XPTB(NPOI,IND,ID)
  281. XVALB(I,IND,IDD2) = (xDE2 - xdm2 ) / pdts2
  282. XVALB(I,IND,IDD3) = XPTB(NPOI,IND,ID+3)
  283. XDE2 = XDE2 + XPSMP(ID)
  284. XDM2 = XDM2 + XPSMPM(ID)
  285. XPALB(I,ID2+ID) = XDE2 - XDM2
  286. XPALB(I,ID3+ID) = XDE2 - XPALB(I,ID5+ID)
  287. XDEP = XDEP + XDE2 * XPALB(I,ID1+ID)
  288. PSN = PSN + XPALB(I,ID2+ID) * XPALB(I,ID1+ID)
  289. PSN0 = PSN0 + XPALB(I,ID3+ID) * XPALB(I,ID1+ID)
  290. 52 CONTINUE
  291. * Projette la vitesse et la variation de deplacement par rapport a
  292. * la position origine d'adherence sur le plan tangent
  293. DO 54 ID = 1,IDIM
  294. XPALB(I,ID2+ID) = (XPALB(I,ID2+ID) - PSN * XPALB(I,ID1+ID))
  295. & / PDTS2
  296. XPALB(I,ID3+ID) = XPALB(I,ID3+ID) - PSN0 * XPALB(I,ID1+ID)
  297. 54 CONTINUE
  298. IF (ITYP.EQ.5) THEN
  299. CALL DYCHE4(XDEP,IDIM,IGP,XPALB,NLIAB,I,XFN,XFT,XPUS,
  300. & iannul)
  301. ELSE
  302. XVITN = PSN / PDTS2
  303. XVALB(I,IND,3) = XVITN
  304. CALL DYCHA4(XDEP,XVITN,IDIM,IGP,XPALB,NLIAB,I,XFN,XFT,XPUS,
  305. & iannul)
  306. ENDIF
  307. XVALB(I,IND,1) = XFN
  308. XVALB(I,IND,10) = ABS(XFT)
  309. XVALB(I,IND,12) = XPUS
  310. IPALB(I,2) = IGP
  311. * Si glissement, memorisation de la vitesse tangentielle et de la force
  312. * tangentielle
  313. IF (IGP.EQ.1.OR.IGP.EQ.-1) THEN
  314. DO 56 ID = 1,IDIM
  315. XPALB(I,ID6+ID) = XPALB(I,ID2+ID)
  316. XPALB(I,ID7+ID) = XPALB(I,ID4+ID)
  317. 56 CONTINUE
  318. ENDIF
  319. * calcul de la force : F = Fn*n + Ft
  320. DO 58 ID = 1,IDIM
  321. XFOR = ( XFN * XPALB(I,ID1+ID) ) + XPALB(I,ID4+ID)
  322. XFNT(ID) = XFOR
  323. FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFOR
  324. 58 CONTINUE
  325. * calcul du moment : M = R ∧ F
  326. XMPFP1 = ( XMP2 * XFNT(3) ) - ( XMP3 * XFNT(2) )
  327. XMPFP2 = ( XMP3 * XFNT(1) ) - ( XMP1 * XFNT(3) )
  328. XMPFP3 = ( XMP1 * XFNT(2) ) - ( XMP2 * XFNT(1) )
  329. XVALB(I,IND,13) = XMPFP1
  330. XVALB(I,IND,14) = XMPFP2
  331. XVALB(I,IND,15) = XMPFP3
  332. FTOTB(NPOI,4) = FTOTB(NPOI,4) + XMPFP1
  333. FTOTB(NPOI,5) = FTOTB(NPOI,5) + XMPFP2
  334. FTOTB(NPOI,6) = FTOTB(NPOI,6) + XMPFP3
  335.  
  336.  
  337. *--------------------------------------------------------------------*
  338. * --- choc elementaire POINT_PLAN_FLUIDE
  339. *
  340. ELSE IF (ITYP.EQ.7) THEN
  341. NPOI = IPLIB(I,1)
  342. IDIM = IPALB(I,3)
  343. XINER = XPALB(I,1)
  344. XCONV = XPALB(I,2)
  345. XVISC = XPALB(I,3)
  346. XPCEL = XPALB(I,4)
  347. XPCRA = XPALB(I,5)
  348. XJEU = XPALB(I,6)
  349. ID1 = 6
  350. ID2 = ID1 + IDIM
  351. XDEP = 0.D0
  352. DO 70 ID = 1,IDIM
  353. XDE2 = XPTB(NPOI,IND,ID)
  354. XVALB(I,IND,ID) = XDE2
  355. XDE2 = XDE2 + FEXPSM(NPOI,NPA,IND1,ID)
  356. XDEP = XDEP + XDE2 * XPALB(I,ID1+ID)
  357. 70 CONTINUE
  358. XDEPM1 = XPALB(I,ID2+1)
  359. XVITM1 = XPALB(I,ID2+2)
  360. XACCM1 = XPALB(I,ID2+3)
  361. XVIT = (XDEP - XDEPM1) / PDTS2
  362. XACC = (XVIT - XVITM1) / PDTS2
  363. XPALB(I,ID2+1) = XDEP
  364. XPALB(I,ID2+2) = XVIT
  365. XPALB(I,ID2+3) = XACC
  366. XDH = XJEU - XDEP
  367. * Calcul de la force d'inertie
  368. CALL DYFINE(XDH,XDEP,XACC,XJEU,XINER,XFI,iannul)
  369. * Calcul de la force de convection
  370. CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFC,iannul)
  371. * Calcul de la force de viscosite
  372. CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFV,iannul)
  373. * Calcul de la force de perte de charge
  374. CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFP,iannul)
  375. *
  376. XFL = XFI + XFC + XFV + XFP
  377. DO 72 ID = 1,IDIM
  378. FTOTB(NPOI,ID) = FTOTB(NPOI,ID) + XFL * XPALB(I,ID1+ID)
  379. 72 CONTINUE
  380. XVALB(I,IND,IDIM+1) = XVIT
  381. XVALB(I,IND,IDIM+2) = XACC
  382. XVALB(I,IND,IDIM+3) = XFI
  383. XVALB(I,IND,IDIM+4) = XFC
  384. XVALB(I,IND,IDIM+5) = XFV
  385. XVALB(I,IND,IDIM+6) = XFP
  386. *
  387. * --- choc ...........
  388. *
  389. * ELSE IF (ITYP.EQ. ) THEN
  390. * .......
  391. * .......
  392. *
  393. ENDIF
  394. *
  395. END
  396.  
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  

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