Télécharger d2vlfa.eso

Retour à la liste

Numérotation des lignes :

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

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