Télécharger devlfa.eso

Retour à la liste

Numérotation des lignes :

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

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