Télécharger devlfa.eso

Retour à la liste

Numérotation des lignes :

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

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