Télécharger hbm26.eso

Retour à la liste

Numérotation des lignes :

hbm26
  1. C HBM26 SOURCE CB215821 23/11/02 21:15:05 11779
  2. *
  3. C DYNE26 SOURCE BP208322 19/04/03 21:15:10 10174
  4. SUBROUTINE HBM26(IBAS,KTKAM,KTLIAB,KTPHI,IA1,IB,ICOMP,RIGIDE,ITKM)
  5. IMPLICIT INTEGER(I-N)
  6. IMPLICIT REAL*8(A-H,O-Z)
  7. *--------------------------------------------------------------------*
  8. * *
  9. * Operateur DYNE : algorithme de Fu - de Vogelaere *
  10. * ________________________________________________ *
  11. * *
  12. * Transpose l'information des objets de Castem2000 dans des *
  13. * tableaux de travail. *
  14. * *
  15. * Parametres: *
  16. * *
  17. * e IBAS Table representant une base modale *
  18. * es KTKAM Segment contenant les matrices XK, XASM et XM *
  19. * es KTPHI Segment des deformees modales *
  20. * e KTLIAB Segment des liaisons sur base B *
  21. * es IA1 Compteur *
  22. * e IB Compteur de la sous base *
  23. * es RIGIDE Vrai si l'on a un corps rigide, faux sinon *
  24. * e ITKM >0 si table RAIDEUR_ET_MASSE fournie *
  25. * *
  26. * Auteur, date de creation: *
  27. * *
  28. * Lionel VIVAN, le 24 octobre 1989. *
  29. * *
  30. *--------------------------------------------------------------------*
  31. -INC PPARAM
  32. -INC CCOPTIO
  33. -INC CCREEL
  34. -INC SMCHAML
  35. -INC SMELEME
  36. -INC SMMODEL
  37. *
  38. *-INC TMDYNC.INC
  39. ************************** debut TMDYNC.INC ****************************
  40.  
  41. * TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
  42. * TODO : a extraire dans un include des que stabilise
  43. *
  44. * Segment des variables generalisees:
  45. * -----------------------------------
  46. SEGMENT MTQ
  47. REAL*8 Q1(NT1)
  48. REAL*8 OMEG,XPARA
  49. REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
  50. REAL*8 dX(NT1), dw, dv
  51. ENDSEGMENT
  52. * Q1 : vecteur des inconnues frequentielles de dimension (2h+1)*n
  53. * Q1 = {q_0 q_c1 q_s1 ... q_sh}
  54. * avec q_i vecteur de dimension n ou n=nombre de modes
  55. * OMEG : frequence fondamentale de l'approximation
  56. * XPARA: parametre de continuation (par defaut la frequence)
  57. * \in [PARINI,PARFIN]
  58. * RX : matrice jacobienne = ZZ + dFnl/dX
  59. * JAC : jacobienne des efforts non-lineaires = dFnl/dX
  60. * ZZ : matrice dynamique associee aux matrices modales K, M et C
  61. * lineaires et constantes
  62. * {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
  63. *
  64. *
  65. * Segment contenant les matrices XK, XASM et XM:
  66. * ---------------------------------------------
  67. SEGMENT MTKAM
  68. REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  69. REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
  70. * REAL*8 GAMFIN(NPC2,nl1)
  71. ENDSEGMENT
  72. * XK,XASM et XM : matrices de raideur, amortissement et masse
  73. * GAM et IGAM : matrices pour la FFT et son inverse
  74. * GAMFIN :
  75. *
  76. * Segment des deformees modales:
  77. * ------------------------------
  78. * (idem DYNE)
  79. SEGMENT MTPHI
  80. INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  81. INTEGER IAROTA(NSB)
  82. REAL*8 XPHILB(NSB,NPLSB,NA2,IDIMB)
  83. ENDSEGMENT
  84. *
  85. * Segment descriptif des liaisons en base A:
  86. * ------------------------------------------
  87. * (idem DYNE)
  88. SEGMENT MTLIAA
  89. INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  90. REAL*8 XPALA(NLIAA,NXPALA)
  91. ENDSEGMENT
  92. *
  93. * Segment descriptif des liaisons en base B:
  94. * ------------------------------------------
  95. * (idem DYNE)
  96. SEGMENT MTLIAB
  97. INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  98. REAL*8 XPALB(NLIAB,NXPALB)
  99. REAL*8 XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  100. ENDSEGMENT
  101. *
  102. * Segment representant les chargements exterieurs:
  103. * -----------------------------------------------
  104. SEGMENT MTFEX
  105. REAL*8 FEXA(NT1)
  106. REAL*8 FEXPSM(NPLB,NPC1,2,IDIMB)
  107. INTEGER BAL
  108. ENDSEGMENT
  109. * FEXA : Vecteur des efforts ext. sous la forme de coefficients de
  110. * Fourier et exprimes en base A
  111. * FEXPSM: chargement/deplacement statique lie aux modes negliges
  112. * (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
  113. * compatibilite avec calcul des Fnl.
  114. * BAL : indique s'il s'agit d'un chargement de type balourd
  115. * (cad proportionnel a OMEG**2)
  116. *
  117. * Segment "local" pour DEVLFA:
  118. * ----------------------------
  119. SEGMENT LOCLFA
  120. REAL*8 FTEST(NA1,4)
  121. ENDSEGMENT
  122. *
  123. * Segment "local" pour DEVLB1:
  124. * ----------------------------
  125. SEGMENT LOCLB1
  126. REAL*8 FTEST2(NPLB,6)
  127. ENDSEGMENT
  128. *
  129. * Segment contenant les variables au cours d un pas de temps:
  130. * ----------------------------------------------------------
  131. SEGMENT MTPAS
  132. REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  133. REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  134. REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  135. REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  136. REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
  137. REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB)
  138. ENDSEGMENT
  139. * FTOTA/B/BA : forces sur base A, B et B projetees sur A
  140. * XPTB : deplacement du point d'une liaison en base B
  141. * XVALA/B : grandeurs de la liaison en base A/B a stocker
  142. * FEXB : forces exterieures en base B (a priori uniquement
  143. * pour les moments appliques aux rotations rigides ?)
  144. * XCHPFB : forces de contact en base B (lorsqu'on considere un
  145. * maillage de contact dans certaines liaisons)
  146. * KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
  147. * base A/B (= contributions a dFnl/dX)
  148. *
  149. *
  150. * Segment des points de reference des modes (base A):
  151. * --------------------------------------------------
  152. SEGMENT MPREF
  153. INTEGER IPOREF(NPREF)
  154. ENDSEGMENT
  155. *
  156. * Segment des points en base B:
  157. * -----------------------------
  158. SEGMENT NCPR(XCOOR(/1)/(IDIM+1))
  159. * NCRP(#global) = #local dans XPTB (1er indice)
  160. *
  161. * Segment des parametres numeriques pour la continuation:
  162. * ------------------------------------------------------
  163. SEGMENT PARNUM
  164. CHARACTER*4 TYPS
  165. REAL*8 DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
  166. REAL*8 PARINI,PARFIN
  167. INTEGER ITERMAX,NBPAS
  168. LOGICAL JANAL
  169. ENDSEGMENT
  170. *
  171. * Segment des resultats:
  172. * ---------------------
  173. SEGMENT PSORT
  174. REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  175. REAL*8 VSAVE(NPAS)
  176. LOGICAL ZSAVE(NPAS)
  177. CHARACTER*2 TYPBIF(NBIFU)
  178. REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  179. REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  180. INTEGER CBIF
  181. ENDSEGMENT
  182. * QSAVE(i,j) = Q harmonique i au pas j
  183. * VSAVE(j) = parametre de continuation (si non w) au j-eme pas
  184. * ZSAVE(j) = stabilite au j-eme pas
  185. * LSAVE(1,j) : partie reelle de l'exposant de Floquet
  186. * LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  187. * TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  188. * QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  189. * WBIF2 : partie imaginaire de l'exposant de Floquet
  190. * QPSIR,QPSII : vecteur propre au point de bifurcation
  191.  
  192. * Segment des tableaux de travail:
  193. * -------------------------------
  194. SEGMENT MTEMP
  195. REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  196. REAL*8 T02(NT1+2), TP2(NT1+2)
  197. INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  198. REAL*8 res
  199. REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  200. REAL*8 QOLD(NT1),OMEGOLD
  201. REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  202. REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  203. ENDSEGMENT
  204. * Jacobiennes augmentees
  205. * Ja : [ RX Rw ; dX dw]
  206. * Jaa: [ RX Rw Ra; gx 0 0; dX dw da]
  207.  
  208. * SEGMENT NNNN
  209. * REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
  210. * ENDSEGMENT
  211.  
  212. *************************** fin TMDYNC.INC *****************************
  213.  
  214. *
  215. LOGICAL L0,L1,RIGIDE
  216. CHARACTER*4 NOMTRI(6),NOMAXI(6),NOMPLA(3)
  217. CHARACTER*8 CMOT,TYPRET,MORIGI,CHARRE
  218. REAL*8 XAXROT(3),XROTA(2,3)
  219. *
  220. * si IFOMOD = -1 : modele PLAN
  221. * si IFOMOD = 0 : modele AXIS
  222. * si IFOMOD = 1 : modele FOUR
  223. * si IFOMOD = 2 : modele TRID
  224. *
  225. * Les noms de composante sont
  226. * - en modele PLAN : UX, UY, RT
  227. * - en modele AXIS : UX, UY, RZ
  228. * - en modele FOUR 1 : UR, UZ, UT, RT
  229. * - en modele TRID : UX, UY, UZ, RX, RY, RZ
  230. *
  231. DATA NOMTRI/'UX ','UY ','UZ ','RX ','RY ','RZ '/
  232. DATA NOMAXI/'UR ','UT ','UZ ','RR ','RT ','RZ '/
  233. DATA NOMPLA/'UX ','UY ','RZ '/
  234. *
  235. MTKAM = KTKAM
  236. MTPHI = KTPHI
  237. MTLIAB = KTLIAB
  238. *
  239. NLIAB = IPALB(/1)
  240. NPLB = JPLIB(/1)
  241. NSB = XPHILB(/1)
  242. NPLSB = XPHILB(/2)
  243. NA2 = XPHILB(/3)
  244. IDIMB = XPHILB(/4)
  245. DEUXPI = 2.D0 * XPI
  246. *
  247. IORSB(IB) = IA1 + 1
  248. IAROTA(IB) = 0
  249. IROT = 0
  250. IN = 0
  251.  
  252. ************************************************************************
  253. * table BASE_MODALE
  254. ************************************************************************
  255.  
  256. 10 CONTINUE
  257. IN = IN + 1
  258. TYPRET = ' '
  259. CALL ACCTAB(IBAS,'ENTIER',IN,X0,' ',L0,IP0,
  260. & TYPRET,I1,X1,CHARRE,L1,IBAMOD)
  261. IF (IERR.NE.0) RETURN
  262. * -on a bien un objet de type table
  263. IF (IBAMOD.NE.0) THEN
  264. IF (TYPRET.EQ.'TABLE ') THEN
  265.  
  266. IA1 = IA1 + 1
  267.  
  268. * remplissage de XM et XK diagonale depuis la table BASE_MODALE
  269. * sauf si deja fait car on a une table RAIDEUR_ET_MASSE !
  270. IF (ITKM.LE.0) THEN
  271. CALL ACCTAB(IBAMOD,'MOT',I0,X0,'MASSE_GENERALISEE',L0,IP0,
  272. & 'FLOTTANT',I1,XMASSE,' ',L1,IP1)
  273. IF (IERR.NE.0) RETURN
  274. XM(IA1,1) = XMASSE
  275. CALL ACCTAB(IBAMOD,'MOT',I0,X0,'FREQUENCE',L0,IP0,
  276. & 'FLOTTANT',I1,XFREQ,' ',L1,IP1)
  277. IF (IERR.NE.0) RETURN
  278. OMEGA = XFREQ * DEUXPI
  279. XK(IA1,1) = XMASSE * OMEGA * OMEGA
  280. IF (IIMPI.EQ.333) THEN
  281. WRITE(IOIMP,*)'HBM26 : XM(',IA1,') =',XMASSE
  282. WRITE(IOIMP,*)'HBM26 : XK(',IA1,') =',XK(IA1,1)
  283. ENDIF
  284. ENDIF
  285.  
  286. * si liaison_B existe, remplissage de IPLSB, XPHILB, IAROTA, INMSB...
  287. IF (NLIAB.NE.0) THEN
  288. CALL ACCTAB(IBAMOD,'MOT',I0,X0,'DEFORMEE_MODALE',L0,IP0,
  289. & 'CHPOINT',I1,X1,' ',L1,ICDM)
  290. IF (IERR.NE.0) RETURN
  291. CALL ACTOBJ('CHPOINT',ICDM,1)
  292.  
  293. DO 12 ID = 1,IDIMB
  294. IF (IFOUR.EQ.0 .OR. IFOUR.EQ.1) THEN
  295. CMOT = NOMAXI(ID)
  296. ELSE
  297. IF (IFOMOD.EQ.-1) THEN
  298. CMOT = NOMPLA(ID)
  299. ELSE
  300. CMOT = NOMTRI(ID)
  301. ENDIF
  302. ENDIF
  303. IF (IIMPI.EQ.333)
  304. & WRITE(IOIMP,*)'HBM26 : composante a extraire :',CMOT
  305. ICOMP = 0
  306. DO 14 IP = 1,NPLB
  307. IPOINT = JPLIB(IP)
  308. * On extrait du chpoint ICDM au point IPOINT de composante CMOT
  309. CALL EXTRA9(ICDM,IPOINT,CMOT,0,.FALSE.,XVAL,IRET)
  310. ICOMP = ICOMP + 1
  311. * on ajuste la taille si necessaire
  312. IF(ICOMP.GT.NPLSB) THEN
  313. NPLSB=ICOMP
  314. SEGADJ MTPHI
  315. ENDIF
  316. IPLSB(IP) = ICOMP
  317. * suite a la modif dans extra9, car on attribue une valeur meme
  318. * si le point n'existe pas dans le chpoint
  319. IF (XVAL.NE.0.) THEN
  320. IF ((IBASB(IP).NE.0).AND.(IBASB(IP).NE.IB)) THEN
  321. call erreur (783)
  322. RETURN
  323. ENDIF
  324. IBASB(IP) = IB
  325. ELSEIF ((IB.EQ.NSB).AND.(IBASB(IP).EQ.0)) THEN
  326. IBASB(IP) = IB
  327. ENDIF
  328. XPHILB(IB,ICOMP,IN,ID) = XVAL
  329. IF (IIMPI.EQ.333) THEN
  330. WRITE(IOIMP,*)'HBM26 : IPLSB(',IP,') =',IPLSB(IP)
  331. WRITE(IOIMP,*)'HBM26 : IBASB(',IP,') =',IBASB(IP)
  332. XVA2 = XPHILB(IB,ICOMP,IN,ID)
  333. WRITE(IOIMP,*)'HBM26 : XPHILB(',IB,ICOMP,IN,ID,') =',XVA2
  334. ENDIF
  335. 14 CONTINUE
  336. 12 CONTINUE
  337. ENDIF
  338.  
  339. c * Prise en compte d'un mode de rotation de corps rigide
  340. MORIGI = ' '
  341. CALL ACCTAB(IBAMOD,'MOT',I0,X0,'CORPS_RIGIDE',L0,IP0,
  342. & MORIGI,I1,X1,CMOT,L1,IP1)
  343. IF (IERR.NE.0) RETURN
  344. IF (MORIGI.EQ.'MOT') THEN
  345. IF (CMOT(1:4).EQ.'VRAI') THEN
  346. CALL ACCTAB(IBAMOD,'MOT',I0,X0,'CENTRE_DE_GRAVITE',
  347. & L0,IP0,'POINT',I1,X1,' ',L1,ICDG)
  348. IF (IERR.NE.0) RETURN
  349. IAROTA(IB)=IA1
  350. IROT = IN
  351. ENDIF
  352. ENDIF
  353. GOTO 10
  354. ELSE
  355. CALL ERREUR(491)
  356. RETURN
  357. ENDIF
  358. ENDIF
  359. * -fin du cas ou on a bien un objet de type table
  360. INMSB(IB) = IN - 1
  361. *
  362. ************************************************************************
  363. * Remplissage des fausses deformees modales de rotations
  364. ************************************************************************
  365. *
  366. *50 continue
  367. IF (IAROTA(IB).NE.0) THEN
  368. RIGIDE = .TRUE.
  369. MERR = 0
  370. NPLUS = IN + 1
  371. IF (NPLUS.GT.NA2) THEN
  372. * On reajuste le dimension NA2 de XPHILB
  373. NA2 = NPLUS
  374. SEGADJ MTPHI
  375. ENDIF
  376. DO 18 IP=1,NPLB
  377. IPOINT=JPLIB(IP)
  378. IPOS=IPLSB(IP)
  379. IBBAS= IBASB(IP)
  380. IF (IBBAS.EQ.IB) THEN
  381. DO 20 ID=(IDIM+1),IDIMB
  382. XAXROT(ID-IDIM) = XPHILB(IB,IPOS,IROT,ID)
  383. 20 CONTINUE
  384. * En tridimensionnel l'axe de rotation est le vecteur propre de rotation
  385. * On norme l axe du plan de rotation
  386. CALL DYNE41(XAXROT,MERR,IDIM)
  387. * En bidimensionnel l'axe de rotation est fixe
  388. * Calcul des fausses deformees modales de rotation
  389. CALL DYNE42(XROTA,XAXROT,IPOINT,ICDG,IDIMB,MERR)
  390. DO 22 ID =1,IDIMB
  391. XPHILB(IB,IPOS,IN,ID) = XROTA(1,ID)
  392. XPHILB(IB,IPOS,IN+1,ID)= XROTA(2,ID)
  393. 22 CONTINUE
  394. ENDIF
  395. 18 CONTINUE
  396. ENDIF
  397.  
  398. IF (IIMPI.EQ.333) THEN
  399. WRITE(IOIMP,*)'HBM26 : INMSB(',IB,') =',INMSB(IB)
  400. WRITE(IOIMP,*)'HBM26 : IORSB(',IB,') =',IORSB(IB)
  401. WRITE(IOIMP,*)'HBM26 : IAROTA(',IB,') =',IAROTA(IB)
  402. ENDIF
  403. END
  404.  
  405.  
  406.  
  407.  
  408.  
  409.  
  410.  

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