Télécharger d2vlfa.eso

Retour à la liste

Numérotation des lignes :

  1. C D2VLFA SOURCE BP208322 18/12/20 21:15:14 10048
  2. c
  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. PDTS2 = PDT
  56.  
  57. *--------------------------------------------------------------------*
  58. * BOUCLE SUR LES LIAISONS
  59. *--------------------------------------------------------------------*
  60.  
  61. DO 10 I = 1,NLIAA
  62.  
  63. ITYP = IPALA(I,1)
  64. icond= IPALA(I,2)
  65. iannul= 0
  66.  
  67. IF (ICOND .NE. 1 ) GOTO 199
  68.  
  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. if(jfonct.eq.100) then
  207. cbp calcul eventuel d'un produit de convolution :
  208. c XDEP=\int_0^T h(\tau)*Qj(t-\tau) d\tau
  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,PDT,XCONV)
  220. XFL=XCPLGE*XCONV
  221. elseif(jfonct.eq.101) then
  222. cbp calcul d'un produit de convolution GRANGER_PAIDOUSSIS
  223. IP1=IPALA(jpliai,4)
  224. IP2=IPALA(jpliai,5)
  225. IP3=IPALA(jpliai,6)
  226. VSD=XPALA(jpliai,2)
  227. XA0=XPALA(jpliai,3)
  228. CALL DYCPL2(IP1,IP2,IP3,VSD,XA0,XDEP,NPAS,PDT,XCONV)
  229. XFL=XCPLGE*XCONV
  230. c WRITE(*,*) '>DYCPL2:t=',T,' x=',XDEP,' F1=',XFL
  231. else
  232. c produit XFL=XCPLGE*XDEP
  233. CALL DYCPLD(XDEP,XCPLGE,XFL,iannul)
  234. cbp on mutliplie eventuellement par une fonction temporelle
  235. if(jfonct.ge.1.and.jfonct.le.2) then
  236. XFREQ = XPALA(jpliai,2)
  237. XTIME = DBLE(NPAS-1)*PDT
  238. if(IND.EQ.2) XTIME=XTIME-PDTS2
  239. if(jfonct.eq.1) XFONCT = COS(XFREQ*XTIME)
  240. if(jfonct.eq.2) XFONCT = SIN(XFREQ*XTIME)
  241. XFL = XFL * XFONCT
  242. endif
  243. endif
  244. XVALA(jpliai,IND,1) = XFL
  245. XVALA(jpliai,IND,4) = XDEP
  246. FTest(INA1,IND) =Ftest(INA1,IND) + XFL
  247. *
  248. * ------ force elementaire de type POLYNOMIALE
  249. *
  250. ELSE IF (JTYP.EQ.6) THEN
  251. * nombre de modes "origine"
  252. NMOD = IPALA(I,2)
  253. CALL D2POL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,PDT,
  254. & NPAS,jpliai,NMOD,FTest,IVINIT,iannul)
  255. *
  256. * ------ choc elementaire ...
  257. *
  258. * ELSE IF (JTYP.EQ. ) THEN
  259. * .......
  260. * .......
  261. *
  262. ENDIF
  263.  
  264. * 2> TESTS POUR VOIR SI ON ANNULERA
  265. xff = 0d0
  266. do 104 ik = 1,na1
  267. do 105 jk = 1,4
  268. xff = xff + ( ftest(ik,jk) ** 2)
  269. 105 continue
  270. 104 continue
  271. xff = xff ** .5
  272. if ( ((xff .le. 1e-20 ) .and. ( jliai .gt. 0) )
  273. & .OR. ((xff .gt. 1e-20 ) .and. ( jliai .lt. 0) ) )
  274. & then
  275. iannul = 1
  276. endif
  277.  
  278. 101 CONTINUE
  279. * --> fin de boucle sur j
  280.  
  281. 199 CONTINUE
  282. * >>> FIN CAS DES LIAISONS CONDITIONNELLES <<<
  283. *
  284.  
  285. * ------ choc elementaire POINT_PLAN sans amortissement
  286. *
  287. IF (ITYP.EQ.1) THEN
  288. XRAID = XPALA(I,1)
  289. XJEU = XPALA(I,2)
  290. XNORM = 1.
  291. IF (XJEU.LT.0D0) THEN
  292. XNORM = -1.
  293. XJEU = -1*XJEU
  294. ENDIF
  295. INA1 = IPLIA(I,1)
  296. XDEP = XNORM*Q1(INA1,IND)
  297. CALL DYCHEL(XDEP,XRAID,XJEU, XFL,iannul)
  298. XVALA(I,IND,1) = XNORM*XFL
  299. XVALA(I,IND,4) = XNORM*XDEP
  300. FTOTA(INA1,IND) = FTOTA(INA1,IND) + XNORM*XFL
  301. *
  302. * ------ choc elementaire POINT_PLAN avec amortissement
  303. *
  304. ELSE IF (ITYP.EQ.2) THEN
  305. XRAID = XPALA(I,1)
  306. XJEU = XPALA(I,2)
  307. XAMO = XPALA(I,3)
  308. XNORM = 1.
  309. IPERM = 0
  310. IF (XJEU.LT.XZERO) THEN
  311. XNORM = -1.
  312. XJEU = -1*XJEU
  313. ENDIF
  314. INA1 = IPLIA(I,1)
  315. XDEP = XNORM * Q1(INA1,IND)
  316. IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  317. XVIT = XNORM * Q2(INA1,IND)
  318. ELSE
  319. IND2 = IND + 1
  320. XDEPM1 = XNORM * Q1(INA1,IND2)
  321. XVIT = (XDEP - XDEPM1) / PDTS2
  322. ENDIF
  323. XVALA(I,IND,3) = XNORM*XVIT
  324. CALL DYCHAM(XDEP,XVIT,XRAID,XJEU,XAMO,XFL,IPERM,iannul)
  325. *+*
  326. IF (XDEP.GT.0.D0 .AND. XFL.GT.0.D0) XFL = 0.D0
  327. IF (XDEP.LT.0.D0 .AND. XFL.LT.0.D0) XFL = 0.D0
  328. *+*
  329. XVALA(I,IND,1) = XNORM*XFL
  330. XVALA(I,IND,4) = XNORM*XDEP
  331. FTOTA(INA1,IND) = FTOTA(INA1,IND) + XNORM*XFL
  332. *
  333. * ------ choc elementaire POINT_PLAN_FLUIDE
  334. *
  335. ELSE IF (ITYP.EQ.3) THEN
  336. XINER = XPALA(I,1)
  337. XCONV = XPALA(I,2)
  338. XVISC = XPALA(I,3)
  339. XPCEL = XPALA(I,4)
  340. XPCRA = XPALA(I,5)
  341. XJEU = XPALA(I,6)
  342. INA1 = IPLIA(I,1)
  343. XDEP = Q1(INA1,IND)
  344. IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  345. XVIT = Q2(INA1,IND)
  346. ELSE
  347. IND2 = IND + 1
  348. XDEPM1 = Q1(INA1,IND2)
  349. XVIT = (XDEP - XDEPM1) / PDTS2
  350. ENDIF
  351. IF (XJEU.GT.0D0) THEN
  352. XDH= XJEU - XDEP
  353. XNORM = 1.
  354. ELSE
  355. XDH= XDEP - XJEU
  356. XNORM = -1.
  357. ENDIF
  358. * Calcul de la masse ajoutee
  359. XXIN = -XINER / XDH
  360. FINERT(INA1,IND) = FINERT(INA1,IND) + XXIN
  361. * Calcul de la force de convection
  362. CALL DYFCON(XDH,XDEP,XVIT,XJEU,XCONV,XFCO,iannul)
  363. * Calcul de la force de viscosite
  364. CALL DYFVIS(XDH,XDEP,XVIT,XJEU,XVISC,XFVI,iannul)
  365. * Calcul de la force de perte de charge
  366. CALL DYFPDC(XDH,XDEP,XVIT,XJEU,XPCEL,XPCRA,XFPE,
  367. & iannul)
  368. XFL = (XFCO* XNORM ) + XFVI + XFPE
  369. XVALA(I,IND,1) = XDEP
  370. XVALA(I,IND,2) = XVIT
  371. XVALA(I,IND,3) = XXIN
  372. XVALA(I,IND,4) = XFCO*XNORM
  373. XVALA(I,IND,5) = XFVI
  374. XVALA(I,IND,6) = XFPE
  375. FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
  376. *
  377. * ------ force elementaire de COUPLAGE EN VITESSE
  378. *
  379. ELSE IF (ITYP.EQ.4) THEN
  380. INA1 = IPLIA(I,1)
  381. INA2 = IPLIA(I,2)
  382. XDEP = Q1(INA2,IND)
  383. IF ((NPAS.EQ.1).AND.(IND.EQ.3)) THEN
  384. XVIT = Q2(INA2,IND)
  385. ELSE
  386. IND2 = IND + 1
  387. XDEPM1 = Q1(INA2,IND2)
  388. XVIT = (XDEP - XDEPM1) / PDTS2
  389. ENDIF
  390. XCPLGE = XPALA(I,1)
  391. CALL DYCPLV(XVIT,XCPLGE,XFL,iannul)
  392. c WRITE(*,*) 't=',(NPAS*PDT),'dx/dt=',XVIT,'F2=',XFL
  393. XVALA(I,IND,3) = XVIT
  394. XVALA(I,IND,1) = XFL
  395. XVALA(I,IND,4) = XDEP
  396. FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
  397. *
  398. * ------ force elementaire de COUPLAGE EN DEPLACEMENT
  399. *
  400. ELSE IF (ITYP.EQ.5) THEN
  401. INA1 = IPLIA(I,1)
  402. INA2 = IPLIA(I,2)
  403. XDEP = Q1(INA2,IND)
  404. XCPLGE = XPALA(I,1)
  405. jfonct = ipala(I,3)
  406. cbp calcul eventuel d'un produit de convolution :
  407. c XDEP=\int_0^T h(\tau)*Qj(t-\tau) d\tau
  408. if(jfonct.eq.100) then
  409. IP1=IPALA(I,4)
  410. c t_{n+1} ou t_{n}
  411. c if (IND.eq.1.or.IND.eq.3) then
  412. c t_{n+1} ou t_{0} pour diff centreee ?
  413. c if (IND.eq.1.or.IND.eq.2) then
  414. IP2=IPALA(I,5)
  415. c c t_{n+1/2} ou t_{n-1/2} (seulement pour devogelaere)
  416. c else
  417. c IP2=IPALA(I,6)
  418. c endif
  419. CALL DYCPL1(IP1,IP2,XDEP,NPAS,PDT,XCONV)
  420. XFL=XCPLGE*XCONV
  421. c WRITE(*,*) 't=',(NPAS*PDT),'x=',XDEP,'F1=',XFL
  422. elseif(jfonct.eq.101) then
  423. cbp calcul d'un produit de convolution GRANGER_PAIDOUSSIS
  424. IP1=IPALA(I,4)
  425. IP2=IPALA(I,5)
  426. IP3=IPALA(I,6)
  427. VSD=XPALA(I,2)
  428. XA0=XPALA(I,3)
  429. CALL DYCPL2(IP1,IP2,IP3,VSD,XA0,XDEP,NPAS,PDT,XCONV)
  430. XFL=XCPLGE*XCONV
  431. c WRITE(*,*) '>DYCPL2:t=',T,' x=',XDEP,' F1=',XFL
  432. else
  433. c produit XFL=XCPLGE*XDEP
  434. CALL DYCPLD(XDEP,XCPLGE,XFL,iannul)
  435. cbp on mutliplie eventuellement par une fonction temporelle
  436. if(jfonct.ge.1.and.jfonct.le.2) then
  437. XFREQ = XPALA(I,2)
  438. XTIME = DBLE(NPAS-1)*PDT
  439. if(IND.EQ.2) XTIME=XTIME-PDTS2
  440. if(jfonct.eq.1) XFONCT = COS(XFREQ*XTIME)
  441. if(jfonct.eq.2) XFONCT = SIN(XFREQ*XTIME)
  442. XFL = XFL * XFONCT
  443. endif
  444. endif
  445. XVALA(I,IND,1) = XFL
  446. XVALA(I,IND,4) = XDEP
  447. FTOTA(INA1,IND) = FTOTA(INA1,IND) + XFL
  448. *
  449. * ------ force elementaire de type POLYNOMIALE
  450. *
  451. ELSE IF (ITYP.EQ.6) THEN
  452. * nombre de modes "origine"
  453. NMOD = IPALA(I,2)
  454. CALL D2POL1(Q1,Q2,NA1,IPLIA,XPALA,XVALA,NLIAA,IND,PDT,
  455. & NPAS,I,NMOD,FTOTA,IVINIT,iannul)
  456. *
  457. * ------ choc elementaire ...
  458. *
  459. * ELSE IF (ITYP.EQ. ) THEN
  460. * .......
  461. * .......
  462. *
  463. ENDIF
  464.  
  465. *
  466. * la suite n'est plus utile car on passe iannul aux s_p de calcul des
  467. * forces de liaisons.
  468.  
  469. * si la liaison etait annulee on l'annule
  470. * if ( ( icond.eq. 1 ) .and. ( iannul.eq.1)) then
  471. * on annulle l'increment de ftotb
  472. * do 112 ik = 1,na1
  473. * do 113 jk = 1,4
  474. * ftota (ik,jk) = ftota0(ik,jk)
  475. * 113 continue
  476. * 112 continue
  477.  
  478. * end if
  479.  
  480. 10 CONTINUE
  481. *--------------------------------------------------------------------*
  482. * Fin de Boucle sur le nombre de liaisons
  483. *--------------------------------------------------------------------*
  484. *
  485. END
  486.  
  487.  
  488.  
  489.  
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  

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