Télécharger d2vlfa.eso

Retour à la liste

Numérotation des lignes :

  1. C D2VLFA SOURCE BP208322 18/01/10 21:15:02 9684
  2. C DEVLFA SOURCE IANIS 96/08/08 21:15:06 2249
  3. SUBROUTINE D2VLFA(Q1,Q2,FTOTA,NA1,IPALA,IPLIA,XPALA,XVALA,NLIAA,
  4. & PDT,T,NPAS,IND,FINERT,IVINIT,FTEST,FTOTA0)
  5. IMPLICIT INTEGER(I-N)
  6. IMPLICIT REAL*8(A-H,O-Z)
  7. *--------------------------------------------------------------------*
  8. * *
  9. * Operateur DYNE : algorithme differences centrees *
  10. * ________________________________________________ *
  11. * *
  12. * Calcul des forces de choc base A. *
  13. * *
  14. * Parametres: *
  15. * *
  16. * e Q1(.,.) Vecteur des deplacements generalises. *
  17. * e Q2(.,.) Vecteur des vitesses generalisees. *
  18. * es FTOTA Forces exterieures totalisees sur la base A. *
  19. * e NA1 Nombre total d'inconnues en base A. *
  20. * e IPALA Renseigne sur le type de la liaison. *
  21. * e IPLIA Tableau contenant les numeros "DYNE" de la liaison. *
  22. * e XPALA Tableau contenant les parametres de la liaison. *
  23. * es XVALA Tableau contenant les variables internes des liaisons *
  24. * e NLIAA Nombre de liaisons sur la base A. *
  25. * e PDT pas de temps *
  26. * e T temps *
  27. * e NPAS Numero du pas de temps *
  28. * e IND Indice du demi-pas de temps *
  29. * = 2 si 1er demi-pas = 1 si 2eme demi-pas *
  30. * es FINERT Forces d'inertie *
  31. * e IVINIT =1 si vitesses initiales, =0 sinon *
  32. * *
  33. * Auteur, date de creation: *
  34. * *
  35. * Lionel VIVAN, le 20 aout 1989. *
  36. * E. de LANGRE, 15 fevrier 91 correction jeu negatif *
  37. * Lionel VIVAN, 12 mars 1991 ajout de la liaison FLUIDE *
  38. * Denis ROBERT, 30 avril 1992 ajout liaison POLYNOMIALE *
  39. * *
  40. * remarque: *
  41. * ========= *
  42. * Si jeu negatif (cas particulier de la base A ou il n'y a pas de *
  43. * normale), on renverse les variables avec XNORM. *
  44. * *
  45. *--------------------------------------------------------------------*
  46. *
  47. INTEGER IPALA(NLIAA,*),IPLIA(NLIAA,*)
  48. REAL*8 XPALA(NLIAA,*),Q1(NA1,*),Q2(NA1,*),FTOTA(NA1,*)
  49. REAL*8 XVALA(NLIAA,4,*),FINERT(NA1,*)
  50. PARAMETER (XZERO = 0.D0)
  51. REAL*8 FTest(nA1,4)
  52. REAL*8 ftOTa0(NA1,4)
  53. *
  54. XFIN = 0.D0
  55. c PDT = XDT(NPAS)
  56. c PDTS2 = 0.5D0 * PDT
  57. PDTS2 = PDT
  58.  
  59. *--------------------------------------------------------------------*
  60. * BOUCLE SUR LES LIAISONS
  61. *--------------------------------------------------------------------*
  62.  
  63. DO 10 I = 1,NLIAA
  64.  
  65. ITYP = IPALA(I,1)
  66. icond= IPALA(I,2)
  67. iannul= 0
  68.  
  69. IF (ICOND .NE. 1 ) GOTO 199
  70.  
  71. * >>> CAS DES LIAISONS CONDITIONNELLES <<<
  72.  
  73. * --> boucle sur j
  74. DO 101 j = 4,20
  75.  
  76. jliai = ipala(i,j)
  77. jpliai = abs ( jliai)
  78. if (jliai.EQ.0) goto 101
  79.  
  80. jtyp = ipala(jpliai,1)
  81. do 102 ik = 1,nA1
  82. do 102 jk = 1,4
  83. ftest(ik,jk) = 0d0
  84. ftota0 (ik,jk) = ftota(ik,jk)
  85. 102 continue
  86.  
  87. * 1> CALCUL DES LIAISONS POUR TEST
  88. *
  89. * ------ choc elementaire POINT_PLAN sans amortissement
  90. *
  91. IF (JTYP.EQ.1) THEN
  92. XRAID = XPALA(jpliai,1)
  93. XJEU = XPALA(jpliai,2)
  94. XNORM = 1.
  95. IF (XJEU.LT.0D0) THEN
  96. XNORM = -1.
  97. XJEU = -1*XJEU
  98. ENDIF
  99. INA1 = IPLIA(jpliai,1)
  100. XDEP = XNORM*Q1(INA1,IND)
  101. CALL DYCHEL(XDEP,XRAID,XJEU, XFL,iannul)
  102. XVALA(jpliai,IND,1) = XNORM*XFL
  103. XVALA(jpliai,IND,4) = XNORM*XDEP
  104. FTest(INA1,IND) =Ftest(INA1,IND) + XNORM*XFL
  105. *
  106. * ------ choc elementaire POINT_PLAN avec amortissement
  107. *
  108. ELSE IF (JTYP.EQ.2) THEN
  109. XRAID = XPALA(jpliai,1)
  110. XJEU = XPALA(jpliai,2)
  111. XAMO = XPALA(jpliai,3)
  112. XNORM = 1.
  113. IPERM = 0
  114. IF (XJEU.LT.XZERO) THEN
  115. XNORM = -1.
  116. XJEU = -1*XJEU
  117. ENDIF
  118. INA1 = IPLIA(jpliai,1)
  119. XDEP = XNORM * Q1(INA1,IND)
  120. IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  121. XVIT = XNORM * Q2(INA1,IND)
  122. ELSE
  123. IND2 = IND + 1
  124. XDEPM1 = XNORM * Q1(INA1,IND2)
  125. XVIT = (XDEP - XDEPM1) / PDTS2
  126. ENDIF
  127. XVALA(jpliai,IND,3) = XNORM*XVIT
  128. CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,XFL,IPERM,iannul)
  129. *+*
  130. IF (XDEP.GT.0.D0 .AND. XFL.GT.0.D0) XFL = 0.D0
  131. IF (XDEP.LT.0.D0 .AND. XFL.LT.0.D0) XFL = 0.D0
  132. *+*
  133. XVALA(jpliai,IND,1) = XNORM*XFL
  134. XVALA(jpliai,IND,4) = XNORM*XDEP
  135. FTest(INA1,IND) =Ftest(INA1,IND) + XNORM*XFL
  136. *
  137. * ------ choc elementaire POINT_PLAN_FLUIDE
  138. *
  139. ELSE IF (JTYP.EQ.3) THEN
  140. XINER = XPALA(jpliai,1)
  141. XCONV = XPALA(jpliai,2)
  142. XVISC = XPALA(jpliai,3)
  143. XPCEL = XPALA(jpliai,4)
  144. XPCRA = XPALA(jpliai,5)
  145. XJEU = XPALA(jpliai,6)
  146. INA1 = IPLIA(jpliai,1)
  147. XDEP = Q1(INA1,IND)
  148. IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  149. XVIT = Q2(INA1,IND)
  150. ELSE
  151. IND2 = IND + 1
  152. XDEPM1 = Q1(INA1,IND2)
  153. XVIT = (XDEP - XDEPM1) / PDTS2
  154. ENDIF
  155. IF (XJEU.GT.0D0) THEN
  156. XDH= XJEU - XDEP
  157. XNORM = 1.
  158. ELSE
  159. XDH= XDEP - XJEU
  160. XNORM = -1.
  161. ENDIF
  162. * Calcul de la masse ajoutee
  163. XXIN = -XINER / XDH
  164. FINERT(INA1,IND) = FINERT(INA1,IND) + XXIN
  165. * Calcul de la force de convection
  166. CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFCO,iannul)
  167. * Calcul de la force de viscosite
  168. CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFVI,iannul)
  169. * Calcul de la force de perte de charge
  170. CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFPE,iannul)
  171. XFL = (XFCO* XNORM ) + XFVI + XFPE
  172. XVALA(jpliai,IND,1) = XDEP
  173. XVALA(jpliai,IND,2) = XVIT
  174. XVALA(jpliai,IND,3) = XXIN
  175. XVALA(jpliai,IND,4) = XFCO*XNORM
  176. XVALA(jpliai,IND,5) = XFVI
  177. XVALA(jpliai,IND,6) = XFPE
  178. FTest(INA1,IND) =Ftest(iNA1,IND) + XFL
  179. *
  180. * ------ force elementaire de COUPLAGE EN VITESSE
  181. *
  182. ELSE IF (JTYP.EQ.4) THEN
  183. INA1 = IPLIA(jpliai,1)
  184. INA2 = IPLIA(jpliai,2)
  185. XDEP = Q1(INA2,IND)
  186. IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  187. XVIT = Q2(INA2,IND)
  188. ELSE
  189. IND2 = IND + 1
  190. XDEPM1 = Q1(INA2,IND2)
  191. XVIT = (XDEP - XDEPM1) / PDTS2
  192. ENDIF
  193. XCPLGE = XPALA(jpliai,1)
  194. CALL DYCPLV(XVIT,XCPLGE,XFL,iannul)
  195. XVALA(jpliai,IND,3) = XVIT
  196. XVALA(jpliai,IND,1) = XFL
  197. XVALA(jpliai,IND,4) = XDEP
  198. FTest(INA1,IND) =Ftest(INA1,IND) + XFL
  199. *
  200. * ------ force elementaire de COUPLAGE EN DEPLACEMENT
  201. *
  202. ELSE IF (JTYP.EQ.5) THEN
  203. INA1 = IPLIA(jpliai,1)
  204. INA2 = IPLIA(jpliai,2)
  205. XDEP = Q1(INA2,IND)
  206. XCPLGE = XPALA(jpliai,1)
  207. jfonct = ipala(jpliai,3)
  208. if(jfonct.eq.100) then
  209. cbp calcul eventuel d'un produit de convolution :
  210. c XDEP=\int_0^T h(\tau)*Qj(t-\tau) d\tau
  211. IP1=IPALA(jpliai,4)
  212. c t_{n+1} ou t_{n}
  213. c if (IND.eq.1.or.IND.eq.3) then
  214. c t_{n+1} ou t_{0} pour diff centreee ?
  215. c if (IND.eq.1.or.IND.eq.2) then
  216. IP2=IPALA(jpliai,5)
  217. c c t_{n+1/2} ou t_{n-1/2} (seulement pour devogelaere)
  218. c else
  219. c IP2=IPALA(jpliai,6)
  220. c endif
  221. CALL DYCPL1(IP1,IP2,XDEP,NPAS,PDT,XCONV)
  222. XFL=XCPLGE*XCONV
  223. elseif(jfonct.eq.101) then
  224. cbp calcul d'un produit de convolution GRANGER_PAIDOUSSIS
  225. IP1=IPALA(jpliai,4)
  226. IP2=IPALA(jpliai,5)
  227. IP3=IPALA(jpliai,6)
  228. VSD=XPALA(jpliai,2)
  229. XA0=XPALA(jpliai,3)
  230. CALL DYCPL2(IP1,IP2,IP3,VSD,XA0,XDEP,NPAS,PDT,XCONV)
  231. XFL=XCPLGE*XCONV
  232. c WRITE(*,*) '>DYCPL2:t=',T,' x=',XDEP,' F1=',XFL
  233. else
  234. c produit XFL=XCPLGE*XDEP
  235. CALL DYCPLD(XDEP,XCPLGE,XFL,iannul)
  236. cbp on mutliplie eventuellement par une fonction temporelle
  237. if(jfonct.ge.1.and.jfonct.le.2) then
  238. XFREQ = XPALA(jpliai,2)
  239. XTIME = DBLE(NPAS-1)*PDT
  240. if(IND.EQ.2) XTIME=XTIME-PDTS2
  241. if(jfonct.eq.1) XFONCT = COS(XFREQ*XTIME)
  242. if(jfonct.eq.2) XFONCT = SIN(XFREQ*XTIME)
  243. XFL = XFL * XFONCT
  244. endif
  245. endif
  246. XVALA(jpliai,IND,1) = XFL
  247. XVALA(jpliai,IND,4) = XDEP
  248. FTest(INA1,IND) =Ftest(INA1,IND) + XFL
  249. *
  250. * ------ force elementaire de type POLYNOMIALE
  251. *
  252. ELSE IF (JTYP.EQ.6) THEN
  253. * nombre de modes "origine"
  254. NMOD = IPALA(I,2)
  255. CALL D2POL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,PDT,
  256. & NPAS,jpliai,NMOD,FTest,IVINIT,iannul)
  257. *
  258. * ------ choc elementaire ...
  259. *
  260. * ELSE IF (JTYP.EQ. ) THEN
  261. * .......
  262. * .......
  263. *
  264. ENDIF
  265.  
  266. * 2> TESTS POUR VOIR SI ON ANNULERA
  267. xff = 0d0
  268. do 104 ik = 1,na1
  269. do 105 jk = 1,4
  270. xff = xff + ( ftest(ik,jk) ** 2)
  271. 105 continue
  272. 104 continue
  273. xff = xff ** .5
  274. if ( ((xff .le. 1e-20 ) .and. ( jliai .gt. 0) )
  275. & .OR. ((xff .gt. 1e-20 ) .and. ( jliai .lt. 0) ) )
  276. & then
  277. iannul = 1
  278. endif
  279.  
  280. 101 CONTINUE
  281. * --> fin de boucle sur j
  282.  
  283. 199 CONTINUE
  284. * >>> FIN CAS DES LIAISONS CONDITIONNELLES <<<
  285. *
  286.  
  287. * ------ choc elementaire POINT_PLAN sans amortissement
  288. *
  289. IF (ITYP.EQ.1) THEN
  290. XRAID = XPALA(I,1)
  291. XJEU = XPALA(I,2)
  292. XNORM = 1.
  293. IF (XJEU.LT.0D0) THEN
  294. XNORM = -1.
  295. XJEU = -1*XJEU
  296. ENDIF
  297. INA1 = IPLIA(I,1)
  298. XDEP = XNORM*Q1(INA1,IND)
  299. CALL DYCHEL(XDEP,XRAID,XJEU, XFL,iannul)
  300. XVALA(I,IND,1) = XNORM*XFL
  301. XVALA(I,IND,4) = XNORM*XDEP
  302. FTOTA(INA1,IND) = FTOTA(INA1,IND) + XNORM*XFL
  303. *
  304. * ------ choc elementaire POINT_PLAN avec amortissement
  305. *
  306. ELSE IF (ITYP.EQ.2) THEN
  307. XRAID = XPALA(I,1)
  308. XJEU = XPALA(I,2)
  309. XAMO = XPALA(I,3)
  310. XNORM = 1.
  311. IPERM = 0
  312. IF (XJEU.LT.XZERO) THEN
  313. XNORM = -1.
  314. XJEU = -1*XJEU
  315. ENDIF
  316. INA1 = IPLIA(I,1)
  317. XDEP = XNORM * Q1(INA1,IND)
  318. IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  319. XVIT = XNORM * Q2(INA1,IND)
  320. ELSE
  321. IND2 = IND + 1
  322. XDEPM1 = XNORM * Q1(INA1,IND2)
  323. XVIT = (XDEP - XDEPM1) / PDTS2
  324. ENDIF
  325. XVALA(I,IND,3) = XNORM*XVIT
  326. CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,XFL,IPERM,iannul)
  327. *+*
  328. IF (XDEP.GT.0.D0 .AND. XFL.GT.0.D0) XFL = 0.D0
  329. IF (XDEP.LT.0.D0 .AND. XFL.LT.0.D0) XFL = 0.D0
  330. *+*
  331. XVALA(I,IND,1) = XNORM*XFL
  332. XVALA(I,IND,4) = XNORM*XDEP
  333. FTOTA(INA1,IND) = FTOTA(INA1,IND) + XNORM*XFL
  334. *
  335. * ------ choc elementaire POINT_PLAN_FLUIDE
  336. *
  337. ELSE IF (ITYP.EQ.3) THEN
  338. XINER = XPALA(I,1)
  339. XCONV = XPALA(I,2)
  340. XVISC = XPALA(I,3)
  341. XPCEL = XPALA(I,4)
  342. XPCRA = XPALA(I,5)
  343. XJEU = XPALA(I,6)
  344. INA1 = IPLIA(I,1)
  345. XDEP = Q1(INA1,IND)
  346. IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  347. XVIT = Q2(INA1,IND)
  348. ELSE
  349. IND2 = IND + 1
  350. XDEPM1 = Q1(INA1,IND2)
  351. XVIT = (XDEP - XDEPM1) / PDTS2
  352. ENDIF
  353. IF (XJEU.GT.0D0) THEN
  354. XDH= XJEU - XDEP
  355. XNORM = 1.
  356. ELSE
  357. XDH= XDEP - XJEU
  358. XNORM = -1.
  359. ENDIF
  360. * Calcul de la masse ajoutee
  361. XXIN = -XINER / XDH
  362. FINERT(INA1,IND) = FINERT(INA1,IND) + XXIN
  363. * Calcul de la force de convection
  364. CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFCO,iannul)
  365. * Calcul de la force de viscosite
  366. CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFVI,iannul)
  367. * Calcul de la force de perte de charge
  368. CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFPE,
  369. & iannul)
  370. XFL = (XFCO* XNORM ) + XFVI + XFPE
  371. XVALA(I,IND,1) = XDEP
  372. XVALA(I,IND,2) = XVIT
  373. XVALA(I,IND,3) = XXIN
  374. XVALA(I,IND,4) = XFCO*XNORM
  375. XVALA(I,IND,5) = XFVI
  376. XVALA(I,IND,6) = XFPE
  377. FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
  378. *
  379. * ------ force elementaire de COUPLAGE EN VITESSE
  380. *
  381. ELSE IF (ITYP.EQ.4) THEN
  382. INA1 = IPLIA(I,1)
  383. INA2 = IPLIA(I,2)
  384. XDEP = Q1(INA2,IND)
  385. IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  386. XVIT = Q2(INA2,IND)
  387. ELSE
  388. IND2 = IND + 1
  389. XDEPM1 = Q1(INA2,IND2)
  390. XVIT = (XDEP - XDEPM1) / PDTS2
  391. ENDIF
  392. XCPLGE = XPALA(I,1)
  393. CALL DYCPLV(XVIT,XCPLGE,XFL,iannul)
  394. c WRITE(*,*) 't=',(NPAS*PDT),'dx/dt=',XVIT,'F2=',XFL
  395. XVALA(I,IND,3) = XVIT
  396. XVALA(I,IND,1) = XFL
  397. XVALA(I,IND,4) = XDEP
  398. FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
  399. *
  400. * ------ force elementaire de COUPLAGE EN DEPLACEMENT
  401. *
  402. ELSE IF (ITYP.EQ.5) THEN
  403. INA1 = IPLIA(I,1)
  404. INA2 = IPLIA(I,2)
  405. XDEP = Q1(INA2,IND)
  406. XCPLGE = XPALA(I,1)
  407. jfonct = ipala(I,3)
  408. cbp calcul eventuel d'un produit de convolution :
  409. c XDEP=\int_0^T h(\tau)*Qj(t-\tau) d\tau
  410. if(jfonct.eq.100) then
  411. IP1=IPALA(I,4)
  412. c t_{n+1} ou t_{n}
  413. c if (IND.eq.1.or.IND.eq.3) then
  414. c t_{n+1} ou t_{0} pour diff centreee ?
  415. c if (IND.eq.1.or.IND.eq.2) then
  416. IP2=IPALA(I,5)
  417. c c t_{n+1/2} ou t_{n-1/2} (seulement pour devogelaere)
  418. c else
  419. c IP2=IPALA(I,6)
  420. c endif
  421. CALL DYCPL1(IP1,IP2,XDEP,NPAS,PDT,XCONV)
  422. XFL=XCPLGE*XCONV
  423. c WRITE(*,*) 't=',(NPAS*PDT),'x=',XDEP,'F1=',XFL
  424. elseif(jfonct.eq.101) then
  425. cbp calcul d'un produit de convolution GRANGER_PAIDOUSSIS
  426. IP1=IPALA(I,4)
  427. IP2=IPALA(I,5)
  428. IP3=IPALA(I,6)
  429. VSD=XPALA(I,2)
  430. XA0=XPALA(I,3)
  431. CALL DYCPL2(IP1,IP2,IP3,VSD,XA0,XDEP,NPAS,PDT,XCONV)
  432. XFL=XCPLGE*XCONV
  433. c WRITE(*,*) '>DYCPL2:t=',T,' x=',XDEP,' F1=',XFL
  434. else
  435. c produit XFL=XCPLGE*XDEP
  436. CALL DYCPLD(XDEP,XCPLGE,XFL,iannul)
  437. cbp on mutliplie eventuellement par une fonction temporelle
  438. if(jfonct.ge.1.and.jfonct.le.2) then
  439. XFREQ = XPALA(I,2)
  440. XTIME = DBLE(NPAS-1)*PDT
  441. if(IND.EQ.2) XTIME=XTIME-PDTS2
  442. if(jfonct.eq.1) XFONCT = COS(XFREQ*XTIME)
  443. if(jfonct.eq.2) XFONCT = SIN(XFREQ*XTIME)
  444. XFL = XFL * XFONCT
  445. endif
  446. endif
  447. XVALA(I,IND,1) = XFL
  448. XVALA(I,IND,4) = XDEP
  449. FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
  450. *
  451. * ------ force elementaire de type POLYNOMIALE
  452. *
  453. ELSE IF (ITYP.EQ.6) THEN
  454. * nombre de modes "origine"
  455. NMOD = IPALA(I,2)
  456. CALL D2POL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,PDT,
  457. & NPAS,I,NMOD,FTOTA,IVINIT,iannul)
  458. *
  459. * ------ choc elementaire ...
  460. *
  461. * ELSE IF (ITYP.EQ. ) THEN
  462. * .......
  463. * .......
  464. *
  465. ENDIF
  466.  
  467. *
  468. * la suite n'est plus utile car on passe iannul aux s_p de calcul des
  469. * forces de liaisons.
  470.  
  471. * si la liaison etait annulee on l'annule
  472. * if ( ( icond.eq. 1 ) .and. ( iannul.eq.1)) then
  473. * on annulle l'increment de ftotb
  474. * do 112 ik = 1,na1
  475. * do 113 jk = 1,4
  476. * ftota (ik,jk) = ftota0(ik,jk)
  477. * 113 continue
  478. * 112 continue
  479.  
  480. * end if
  481.  
  482. 10 CONTINUE
  483. *--------------------------------------------------------------------*
  484. * Fin de Boucle sur le nombre de liaisons
  485. *--------------------------------------------------------------------*
  486. *
  487. END
  488.  
  489.  
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  
  497.  

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